Biophysical Journal - Torsional Rigidities of Weakly Strained DNAs

Article Outline

Abstract
Introduction
Theory
Results and discussion
Conclusions
Acknowledgments
Appendix
References
Publication Information

Article Information

PubMed

Monte Carlo Simulations of Locally Melted Supercoiled DNAs i...

Dynamic Bending Rigidity of a 200-bp DNA in 4mM Ionic Streng...

Restrained Torsional Dynamics of Nuclear DNA in Living Proli...

…more

Copyright © 2006 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 91, Issue 11, 4166-4179, 1 December 2006

doi:10.1529/biophysj.106.087593

Nucleic Acids

Torsional Rigidities of Weakly Strained DNAs

Bryant S. Fujimoto, Gregory P. Brewood and J. Michael Schurr^,

Department of Chemistry, University of Washington, Seattle, Washington

Address reprint requests to J. Michael Schurr, Dept. of Chemistry, University of Washington, Box 351700, Seattle, WA 98195. Tel.: 206-543-6681; Fax: 206-685-8665.

Abstract

Measurements on unstrained linear and weakly strained large (≥340bp) circular DNAs yield torsional rigidities in the range C=170–230fJ fm. However, larger values, in the range C=270–420fJ fm, are typically obtained from measurements on sufficiently small (≤247bp) circular DNAs, and values in the range C=300–450fJ fm are obtained from experiments on linear DNAs under tension. A new method is proposed to estimate torsional rigidities of weakly supercoiled circular DNAs. Monte Carlo simulations of the supercoiling free energies of solution DNAs, and also of the structures of surface-confined supercoiled plasmids, were performed using different trial values of C. The results are compared with experimental measurements of the twist energy parameter, E_T, that governs the supercoiling free energy, and also with atomic force microscopy images of surface-confined plasmids. The results clearly demonstrate that C-values in the range 170–230fJ fm are compatible with experimental observations, whereas values in the range C≥269fJ fm, are incompatible with those same measurements. These results strongly suggest that the secondary structure of DNA is altered by either sufficient coherent bending strain or sufficient tension so as to enhance its torsional rigidity.

Introduction

The torsional rigidities of many different DNAs have been assessed by a variety of experimental methods. Surprisingly, the reported values span a 4.5-fold range, even though statistical errors in the individual measurements are probably <15% in most cases. Moreover, the variations in torsional rigidity (C) with 1), temperature (T) from 293 to 310K; 2), overall base composition from 34% to 100% GC; 3), NaCl concentration from 0.01 to 1.0M; or 4), Mg²⁺concentration from 0 to 40mM in 0.1M NaCl are all <∼15%, so the large spread in reported values must stem from other causes 1,2,3,4,5,6. All of the larger reported C-values pertain to DNAs that are subject to substantial bending strain in small circles with N≤254bp 7,8,9,10,11,12,13,14,15,16,17,18,19, or substantial tensile stresses in single-molecule pulling experiments 20,21,22,23,24,25,26. There arises now the question of whether the large spread in reported C-values reflects genuinely different torsional rigidities of the differently strained DNAs, or instead simply reflects large systematic errors in one or another of the experimental methods. In this article, we shall present new evidence that the large values of C reported for sufficiently bent and stretched DNAs do not apply to weakly strained, large, circular DNAs. We shall also argue that such large torsional rigidities more likely arise from strain-induced alterations of the DNA secondary structure than from systematic errors in the measurements themselves.

Variations in secondary structure and torsional rigidity

Considerable evidence that duplex DNAs under various conditions exhibit different average secondary structures with different torsional rigidities was presented previously 2,3,5,6,10,27,28,29,30,31,32,33,34,35,36,37,38. Certain local perturbations, such as a particular change in sequence or the binding of a CAP or Sp1 transcriptional activator to its specific site, were found to exert very long-range effects (over several hundred base pairs) on the average secondary structure and torsional rigidity of the flanking DNA 6,36. In addition, the dynamic bending rigidity (A_d) of DNA, which governs Brownian flexure for times t≤1μs, was found to exceed the static equilibrium value (A_eq) by ∼4-fold 37,39,40,41. This important finding implies that two or more differently curved (or superhelical), but slowly interconverting, conformations coexist even in unstrained duplex DNAs. A substantial increase in torsional rigidity upon imposing a coherent bending strain of ∼2°/bp was attributed to the shift of such a prevailing conformational equilibrium toward the more curved structure, which evidently has the greater torsional rigidity 35,37. The very limited available structural information suggests that these different structures are probably conformational substates within the B-family. Additional evidence pertaining to long-range allosteric transitions in DNA secondary structure and their possible role in gene regulation was reviewed previously 36.

Methods to measure the torsional rigidity

Several rather different methods have been employed to determine the torsional rigidities of particular DNAs under various conditions. In the fluorescence polarization anisotropy (FPA) method, the torsional rigidity is obtained by analyzing the time-resolved FPA of intercalated ethidium dye 5. This analysis requires an estimated value of A_d4,5. Experimental values of 1), A_d; 2), the equilibrium persistence length (P_eq); and 3), C for a 200-bp DNA were determined by combining FPA measurements with transient polarization grating (TPG) experiments on the same DNA 37. The TPG method involves diffractive detection of the photo-induced dichroism within the grating, and extends the time window for measurement of the optical anisotropy by 100-fold into the regime where bending and end-over-end tumbling completely dominate the decay 42. For purposes of comparison, A_d is expressed in terms of a dynamic persistence length, , where k is Boltzmann’s constant and T the absolute temperature. For this 200-bp DNA in 4mM ionic strength at 294K, the values C=188fJ fm, =200nm, and =50nm, were obtained. Robinson and co-workers employed electron paramagnetic resonance spin-label methods to assess , and obtained comparably large values in the range =150–170nm for numerous synthetic DNAs in 0.1M NaCl 40,41,43. The values of C obtained from FPA data on long DNAs are very insensitive to the choice of in the range 150–200nm, and all of the C-values quoted below are obtained by using P_d=180nm. Possible systematic errors in C-values obtained by the FPA method are discussed in the Appendix. The maximum possible systematic error is judged to be <20%, and the actual systematic error is likely <10%.

In the topoisomer ratio (TR) method, a short (205≤N≤250bp) linear DNA is circularized by ligation, the resulting topoisomers are separated by gel electrophoresis, and their relative populations are measured. By repeating this experiment in different concentrations of ethidium 8,10, or with DNAs of slightly different length 11, it is possible to extract the difference in free energy between topoisomers, as well as the intrinisic twist per base pair in the absence of ethidium. From such free energy differences, C could be extracted by appropriate modeling 12,13,44.

In the cyclization kinetics (CK) method, one measures the ratio of the rate of formation of circular monomers to the rate of formation of linear plus circular dimers of short (150≤N≤350bp) DNAs 7,14,15,16,17,18,19,45. This ratio provides information about the free energy to form the noncovalently closed circular monomer that is sealed by the ligase. When measured for homologous DNAs of different lengths, this ratio is found to vary with length in an oscillatory manner with an ∼10.4bp period. Again, appropriate modeling of such data provides estimates of both C and the equilibrium bending rigidity () 9,45. It is important to note that the enzyme ligates the two strands of the duplex in successive steps. In the first step, the ligase may tolerate a significant departure from perfect tangential and, especially, azimuthal alignment of the cohesive ends in the noncovalently closed DNA circles upon which it acts. This tolerance may significantly lower the deformational free energy of the noncovalently closed circular DNA that is sampled by the ligase. Because the rate of formation of monomer circles is determined entirely by the initial ligation event, the CK method may provide only a lower bound, rather than the actual, value of C17. Because equilibration of the topoisomers presumably still occurs between the two ligation events, the second ligation step should still yield the desired equilibrium ratio of topoisomers. Hence, the TR results should be valid, even when the CK results significantly underestimate the actual value of C.

The simplest protocols for modeling TR and CK data are valid only in the absence of complicating factors, such as unsuspected directional permanent bends or twisting and bending rigidities that vary with position along the filament. Recent experiments suggest that any directional permanent bends that are present in unstrained native DNA sequences may not significantly affect the probability of forming small circles 17. Numerous attempts to model and characterize such complicating factors have been made 14,15,16,17,18,19. These attempts all invoke the independent dinucleotide step model, wherein the secondary structure and torsional and bending rigidities associated with a given base-pair step are assumed to be completely independent of the sequence of its flanking DNA, or of any altered state of its flanking DNA that is induced by protein or other ligand binding, or by a B-to-Z transition. However, this assumption has been unequivocally contradicted by numerous published NMR structures of small duplexes, which yield different structural parameters for particular steps, e.g., A-A steps, embedded in different flanking sequences 46, as well as by numerous other published experiments 6,36,47,48,49,50,51,52 and unpublished studies of S. A. Winkle (Florida International University, personal communication, 1997) that were reviewed previously 36. Other evidence indicates that dinucleotide step models for directional permanent bends, or for , cannot account satisfactorily for the behavior of all sequences 43,53,54. The most common DNA melting models, namely two-state nearest-neighbor models, can also be regarded as two-state dinucleotide step models. These, too, cannot account for all of the extensive melting data 49. It is important to note that the failure of dinucleotide step models is more likely due to the assumption of a single duplex state than to interactions of appreciably longer range than a single base pair. At the midpoint of a cooperative transition between two duplex conformations, the average domain size is larger, possibly very much larger, than a dinucleotide step, so the spatial range of structural correlations may far exceed the range of the interactions per se 36. For such reasons, the values of C and other parameters of specific subsequences in small circles, which are extracted from CK measurements on multiple DNAs using the dinucleotide step model, are potentially prone to greater systematic error than “overall” values that are obtained without the use of that model.

Analyses of various single-molecule pulling (SMP) experiments provided several new means to assess the torsional rigidities of variously twisted DNAs under tension 20,21,22,23,24,25,26.

Reported values of torsional rigidity

Values of C obtained by different methods are presented in Table 1, where they are ranked in order of increasing coherent bend for circular DNAs, or increasing tension in the case of SMP experiments. Omitted from Table 1 are lower C-values in the range 103–170fJ fm, which were obtained for subsequences of 150–170bp circular DNAs only in the presence of both phased A-tracts and either bound CAP or repeats of certain other sequences 15,16,18,19. As noted above, such values depend heavily on the validity of the independent dinucleotide step model, which has already been contradicted in many cases.

Table 1 C-values obtained by different methods


DNA	Bend (°/bp)	Tension (pN)	Method C(fJ fm)	Reference(s)
Linear viral DNAs	0	0	FPA 150–170	5
Linearized plasmids	0	0	FPA 190–220	2,3,4,5,35,37
Linear 181bp DNAs	0	0	FPA 220	35
Circular plasmids	≤0.40^*	0	FPA 200–230	5,10,34,64,67
Circular (340–350bp)	1.0^†	0	CK 200	45
Circular (237–254bp)	1.45^†	0	CK 240–300	7,9
Circular (247bp)	1.45^†	0	TR 410–420	8,10,35
Circular (205–217bp)	1.7^†	0	TR 320–330	11,12,13
Circular (150–200bp)	1.8–2.4^‡	0	CK 220–340^(b)	14,15,16,17,18,19
Circular (181bp; fresh)	2.0^†	0	FPA 310–330	35
Circular (181bp; 8 mos)	2.0^†	0	FPA 400	35
Linear	None	0.1–2.0	SMP 300	20,21,22,25
Linear	None	0.1–5.0	SMP 350	20,21,22,24
Linear	None	0.3–8.0	SMP 450	20,21,22,23
Linear	None	15–45	SMP 410–440	26

^* RMS coherent bend (°/bp) was estimated for the circular pSA509 plasmid with |σ|=0.047 in 161mM ionic strength, as follows. First, the fluctuation bend energy, MkT, was subtracted from the mean bending energy for λ=0 in Table 2 to determine the energy attributable to coherent bend,

, where M=134 rods,

erg, and

is the mean-squared coherent bend angle between adjacent subunit rods. After solving for

, the root mean-squared coherent bend in °/bp is reckoned as,

, where N=3760bp. This value depends upon superhelix density, but should be independent of DNA length for such large DNAs.
^† For small circles containing N≤350bp, the coherent bend is predominantly planar and is estimated simply as 360/N.
^‡ These values apply when the circle contains no bound CAP or phased repeats of other sequences besides oligo-A tracts.

Upon circularization of the linear 181-bp DNA, its torsional rigidity increased by ∼1.5-fold, as indicated in Table 1. In addition, its intrinsic binding constant for intercalated ethidium increased by ∼4.3-fold, and its circular dichroism spectrum changed significantly 35. These changes imply that the coherent bending strain inherent in these 181-bp circles induced a conformational change to a torsionally stiffer state. However, this conformational change was evidently not complete. After 8 months, the torsional rigidity of these 181-bp circles rose still further from ∼310–330fJ fm to 400fJ fm, whereas the corresponding 181-bp linear DNAs remained unchanged. In addition, the ratio of the intrinsic ethidium binding constant of these “aged” circles to that of the corresponding linear DNAs declined from ∼4.3 to ∼1.6. Within experimental error, the torsional rigidity and ethidium binding constant ratio of these “aged” 181-bp circles matched the corresponding properties of the 247-bp circles, which were determined via the TR method 10,35. Thus, the higher torsional rigidities (410–420 and 400fJ fm) of these 247- and aged 181-bp circles might well represent the true long-term equilibrium values in such small circles with 181≤N≤247 bp.

All of the C-values determined for circular DNAs with N≤247bp via the FPA and TR methods lie in the range 300–430fJ fm. Some of the C-values obtained by the CK method likewise exceed 300fJ fm, although values as low as ∼220–240fJ fm were also reported. All of the C-values for equilibrium unstrained linear DNAs of any length and for circular DNAs with N≥340bp lie in the range 150–230fJ fm. These findings strongly suggest that sufficient coherent bending strain (∼1.5–2.0°/bp) alters the average secondary structure in such a way as to significantly increase the torsional rigidity, whereas significantly smaller bending strains (≤1°/bp) do not.

Analyses of SMP experiments on variously twisted DNAs all yielded C-values in the range 300–450fJ fm (cf. Table 1). The best-fit C-value for each given range of tensions is seen to rise with increasing tension. This latter finding suggests that tension, too, may alter the average secondary structure in such a way as to increase the torsional rigidity.

At present, the reported C-values extend from 103 to 450fJ fm. Till now, only the FPA method was capable of providing C-values for unstrained linear DNAs and weakly strained large circular DNAs (with superhelix density, |σ|≤0.05), as well as for more highly curved small (181-bp) circles. It would be highly desirable to have an independent method to assess the torsional rigidities of unstrained and/or weakly strained DNAs. Comparatively precise measurements of the supercoiling free energy have been made for several different samples of p30δ DNA (4932bp) at low superhelix densities, |σ|≤0.008 33,38,55, where the strain free energy associated with the coherent (superhelical) deformation is ≤0.7kT/1000bp. Simulations of the supercoiling free energies of topoisomers with small linking differences, which correspond to such low superhelix densities, can now be made with sufficient statistical accuracy to allow a substantial narrowing of the range of C-values that are compatible with the experimental data.

Writhe as a probe of the torsional rigidity

In contrast to the torsional rigidity, there is a reasonable consensus concerning the value of the equilibrium bending rigidity () and equilibrium persistence length (=50nm) of DNA. Recent CK data suggest that sequence-dependent directional permanent bends in native DNA sequences contribute negligibly to , in which event for native sequences 17. For a supercoiled DNA, the relative magnitudes of C and govern the partitioning of the superhelical strain into twist and writhe. Atomic force microscopy (AFM) images allow one to estimate the writhe of a surface-confined supercoiled pSA509 DNA (3760bp) under buffer 56. Previous Monte Carlo simulations demonstrated that confinement in a surface flattening potential restricts the available configurations to an extremely small and unrepresentative subset of those exhibited in solution 57. The average writhe of a surface-confined supercoiled DNA is greater than that of the same DNA in solution. Nevertheless, the tertiary structures of the surface-confined supercoiled DNA are influenced by the torsional rigidity, so comparisons of simulated surface-confined supercoiled model DNAs with the observed AFM images may provide some additional information about the range of admissible values of C.

Objectives of this work

The aims of this study are to simulate 1), the supercoiling free energies of model circular DNAs with low superhelix densities in solution; and 2), the mean writhes and typical structures of surface-confined model circular DNAs with native superhelix density by using different trial values of the torsional rigidity. By comparing the results with experimental observations, the range of C-values that is compatible with the experimental data is ascertained.

Theory

Topological and geometrical aspects of supercoiling

Every circular duplex DNA is characterized by its integral linking number ℓ, the number of turns of one single strand around the other. ℓ is a topological invariant that is unaltered by any change in secondary or tertiary structure of the DNA. The extent of deformation of the DNA is characterized by its linking difference, , wherein is the intrinsic twist of the unstrained DNA. For a large circular DNA comprising N bp, , where ((1/10.45) turns/bp) is the intrinsic succession angle between base pairs of the duplex. The superhelix density is defined by , and is typically ∼–0.05 for native plasmid DNAs.

The linking number is partitioned between twist (t) and writhe (w) according to 58,59

(1)

Numerous Monte Carlo studies of supercoiled circular DNAs have been performed in our laboratory 2,44,57,60,61,62,63,64. Detailed descriptions of the mesoscopic models, coordinate systems, and potentials of mean force, and full explications of the algorithms and protocols for uniform sampling of configuration space, ring closure, detection of hard-cylinder overlap, knot detection, and reversible work calculations were provided in one or another of those prior works 44,57,60,63,64. Consequently, only a brief account of our models and methods is presented here, with emphasis on any changes from previously published procedures.

The mesoscopic model

The DNA is here regarded as a chain of M rigid-rod subunits, each containing ν=N/M bp. These subunit rods are labeled consecutively by the index j, j=1, … M. In each subunit (j) is fixed a coordinate frame (x_j,y_j,z_j), the z_j axis of which lies along the bond vector (b_j) that extends from the origin of the jth frame to the origin of the succeeding (j+1)th frame. The Euler rotation that carries a coordinate frame from coincidence with the jth frame to coincidence with the (j+1)th frame is , where the component rotations, , are taken sequentially around the body-fixed z_j, new body-fixed y_j′, and final body-fixed z_j″ axes 65. The net twist of the DNA is given by

(2)

where

(rad) is the net twist of the Euler rotation,

, from the jth to the (j+1)th frame 44,57,60,63.

The writhe is commonly approximated by the discretized Gauss integral,

(3)

wherein r_i and r_j denote the positions of the origins of the ith and jth subunit frames, respectively, in the laboratory frame, and

denotes a unit vector along r_i-r_j66. For the subunit rod length,

=9.54nm, that is used in the present simulations, the discretization implied by Eq. (3) does not provide a sufficiently accurate approximation to the actual writhe of the array of bond vectors. By subdividing each bond vector into 10 collinear subsections and performing the analogous sum over those, the accuracy of the writhe calculation becomes adequate for our purposes 63. Specifically, Eq. (1) was obeyed to within±0.1 turn throughout these simulations, which extend to |Δℓ|=24 turns.

Monte Carlo simulation models and methods

Each rigid-rod subunit of our model DNA is connected to its neighbors at either end by Hookean torsion and bending springs. The subunit rod length, b=9.54nm, was chosen to be slightly <1/5 of a persistence length, P_eq=50nm. Various properties, including the mean-squared writhe of a nicked circular DNA and the supercoiling free energy of a closed circular DNA, were found to become independent of rod length, when b≤(P_eq/5), as shown by simulations of nicked circles 67 and by (unpublished) simulations of closed circles in our laboratory. This choice of b corresponds to ν=28.06bp. Two different model DNAs are simulated. The first model comprises M=170 subunits, corresponding to 4770bp. This is comparable in size to p30δ (4932bp), whose supercoiling free energy was measured in numerous experiments 33,38,55. The second model comprises M=134 subunits, corresponding to 3760bp. This is identical in size to pSA509, which was imaged via AFM 56. The trial twisting torque constants are chosen to be uniform along the filament and to produce one or another torsional rigidity in the range C=170–410fJ·fm. The bending torque constant is always chosen to yield a total bending persistence length P_eq=50nm.

Potential functions

The potential of mean force of a particular configuration is given by

(4)

where

is the torsion potential energy,

is the bending potential energy, and

is the intersubunit potential energy and includes both the hard-cylinder

and electrostatic

interactions between different rigid-rod subunits of the model DNA 57,60,63,64.

is the energy due to an external potential, which confines the DNA near a planar surface, and is used only in simulated transfers of a supercoiled DNA from solution to a surface for comparison with AFM images.

The torsion potential is given by

(5)

where α≡C/|b| is the trial torsion elastic constant between subunits, and θ is a parameter that is used to vary the linking difference, while the linking number remains fixed at 0 turns 44,63.

The hard-cylinder potential, , is evaluated by regarding each subunit rod as a cylinder with radius a=1.20nm about its bond vector. is infinite, if the cylinders overlap, and zero otherwise. Any overlaps are detected by a previously described algorithm 60, and any overlapped configuration is immediately rejected, because it lies outside the accessible region of configuration space, and a new move is attempted from the same previous configuration.

The electrostatic potential between subunit rods is evaluated by first placing charged spheres of radius, a=1.20nm, separated by 3.18nm along the entire chain of bond vectors, so exactly three spheres are found at identical positions (0, b/3, and 2b/3) along every bond vector. The electrostatic interaction between the ith and jth charged spheres is taken to be 64

(6)

where

is the distance between charges, e is the electronic charge, ɛ is the dielectric constant (ɛ=78.54 at 298K), κ is the Debye screening parameter due to small ions in the solution, and Z is the effective charge. Z is adjusted so that the electrostatic potential surrounding the middle of a linear array of 2001 such charges closely matches the solution of the nonlinear Poisson-Boltzmann equation for an infinitely long cylinder with the same charge per unit length as DNA at all distances >2a57,63,64. This has been done previously for uni-univalent salt concentrations in the range 0.01–1.0M, and for mixed salt solutions containing 0.10M NaCl+0.01M MgCl₂64. For a given value of the effective charge, Eq. (6) is the well-known Derjaguin-Landau-Verwey-Overbeek electrostatic potential between charged spheres in an ionic solution. The electrostatic energy of the supercoiled model DNA is reckoned by summing

over all pairs of spherical charges that reside on different subunit rods. Interactions between those pairs of charges that are separated by more than a specified cutoff distance are neglected. The cut-off distance is chosen so that

/kT<8×10⁻⁴ when

exceeds the cut-off distance. For simulations of the supercoiling free energies of model DNAs in ∼0.1M uni-univalent ionic strength, all electrostatic interactions between adjacent rods (along the chain) are ignored, so the local resistance to bending is determined almost exclusively by the bending torque constant,

. However, in simulations of the transfer of model DNAs from solution to the surface in both 0.161 and 0.01M uni-univalent ionic strength, all electrostatic interactions between adjacent subunit rods are included. Such interactions significantly affect the resistance to bending and the persistence length, and consequently also the choice of

that gives the desired persistence length, as described previously 57,63,64 and discussed further below.

When used in simulations with C=200 or 210fJ fm, the electrostatic potential in Eq. (6) yielded results in good agreement with the measured supercoiling free energy of pBR322 DNA in 0.02M ionic strength 63,68 and with AFM images of a native supercoiled pSA509 DNA in 0.161 and 0.010M ionic strength 56,57.

For simulations of the supercoiling free energies of model DNAs in 0.1M ionic strength, 1/κ=0.962nm, Z=−14.4 charges/sphere, and the cutoff distance is 12.0nm. For the simulations of surface-confined DNAs in 0.161M ionic strength, 1/κ=0.758nm, Z=−16.7 charges/sphere, and the cut-off distance is 12.0nm, and for surface-confined DNAs in 0.01M ionic strength, 1/κ=3.04nm Å, Z=−7.82 charges/sphere, and the cut-off distance is 24.0nm.

The bending potential energy is given by

(7)

where κ_B is the torque constant for bending of the DNA, and the bending angle,

, is the second rotation in the composite Euler rotation,

, that orients the frame of the (j+1)th subunit in that of the jth. The value of κ_B was chosen to yield a persistence length, P=50nm, in the following manner. Simulations of both 10- and 20-subunit linear chains were performed (10million moves/simulation) for several different values of κ_B, but using always the same appropriate values of the parameters in

and

that were discussed above. For each simulation, the persistence length was calculated according to

(8)

where the average,

, was taken over all subunits, j=1, …, M-1, in all configurations. The values,

=1.946×10⁻¹³, 1.806×10⁻¹³, and 1.060×10⁻¹³ dyne cm in, respectively, 100mM, 161mM, and 10mM ionic strength yielded P=50nm for both the 10- and 20-subunit chains within the simulation errors. The

-values in 161 and 10mM ionic strength are lower than that in 100mM, because the simulations in 161 and 10mM include electrostatic interactions between adjacent subunits, whereas the simulations in 100mM did not.

Each DNA in solution is simulated in the absence of via previously described algorithms and protocols 44,60,63.

For simulations of surface-confined DNAs, a surface potential of mean force, , is applied along the laboratory z axis, normal to the surface, to force the DNA to lie in the xy plane. This potential is applied separately to each subunit in the DNA, and takes the form

(9)

where z is the position of the subunit in Å and λ is a constant that is slowly increased from 0 to 1 (to turn on the attractive part). Because the force constants in Eq. (9) apply for z values in angstroms (Å), the unit of length in this and all subsequent discussions pertaining to surface-confined DNAs is taken to be 1Å=0.1nm. For the purpose of calculating

, the location of each subunit is taken to be at the joint between adjacent subunits. The value of the force constant for z≤0 was chosen so that the potential energy of a single subunit equals kT at −50Å and has units of erg/Å¹⁰. The value of the force constant for z>0 was chosen so that when λ=1, the potential energy equals kT at 150Å, and has units of erg/Å²57. As explained previously 57, this is an ad hoc potential that 1), hopefully reproduces certain features of the actual force field in the region around the potential minimum at z=0; 2), achieves a suitable degree of flattening of the DNA in the z direction normal to the surface; and 3) allows considerable equilibration of the surface-bound DNA within a reasonable simulation time. The force constant of the hard wall potential for z<0 was chosen to enable DNAs to equilibrate by passing DNA sections over each other, even when pressed onto the surface with mean positions near z∼0. The force constant of the harmonic potential (when λ=1.0) in the region z>0 was chosen so that the mean z thickness in 161mM ionic strength is ∼½ the diameter of an unperturbed interwound superhelix. At that point, the induced deformation is already quite substantial. A practical consideration is that any significant increase in this harmonic force constant would not only increase the degree of flattening, but would also greatly slow the equilibration process, and prohibitively increase the required simulation time. This ad hoc potential provides a force field with suitable properties for our purposes, which are to assess the effects of near-equilibrium flattening on the geometric, topological, and thermodynamic properties associated with the internal coordinates of the DNA. However, the value of this surface potential at any particular z coordinate relative to that of the corresponding solution DNA is certainly incorrect, so that the total work of transfer from solution to the surface is also incorrect, and cannot be used to estimate the equilibrium constant for surface binding.

Neglect of the hydration potential

Osmotic compression experiments revealed short-range repulsive forces between DNA duplexes that substantially exceed the expected electrostatic forces at short separations between duplex axes, d≤3.0nm 69. These short-range repulsions, which were attributed to hydration forces, are omitted from the interaction potential in Eq. (4). The resulting underestimate of is extremely slight, since the DNA subunits practically never approach to such close distances. The hard-core repulsion in precludes smaller separations, d<2.4nm, in any case.

Monte Carlo moves and simulation algorithms

DNAs in solution

In the basic Monte Carlo move (henceforth called a type I move), a different small random rotation about every local coordinate axis, x_j, y_j, and z_j, j=1,…,N, is applied simultaneously to every subunit in the chain 44,60, after which the Euler angles, , that orient each subunit in the laboratory frame are updated via the small-angle formulas 65. The maximum angle of rotation is very small, 0.008rad or less, so the rotations of each subunit around its three Cartesian axes commute to very high accuracy. This protocol samples configuration space in a uniform manner 60. Erratic behavior that occurs when is too close to either 0 or π, is removed by a π/2 rotation of the coordinate frame of the jth subunit around its own y_j axis, whenever or , where Θ is a suitable threshold angle. Details of this protocol were provided previously 44,60. In our early work 2,44,60,61,62, Θ was assigned a value corresponding to 17°, but in our recent simulations 57,63, and in the present study, this value was increased to 28°. The Euler angles, , are used to express the unit vectors, , of the jth coordinate frame in terms of the laboratory coordinates 44,60. At this point the ring-closure algorithm 60 is applied, which may result in additional small rotations of each subunit around its transverse axes, which in turn necessitates an update of the and the unit vectors. Next, the algorithm to detect hard-cylinder overlaps 60 is applied. If the new configuration exhibits no overlaps, then the knot-detection algorithm 60 is applied. If the new configuration is unknotted, then the frame-to-frame Euler rotations, , are reckoned from the and , as described previously 44. Thus, each type I move comprises 3N random subrotations, plus any small rotations that result from the ring-closure algorithm. Subsequently, the potential energy of the new configuration is evaluated and the Metropolis criterion applied to either accept the new configuration or reject it in favor of its predecessor. In the simulation used here, the potential function is evaluated directly from the frame-to-frame Euler angles on every step, so that it is not necessary to evaluate the writhe on every step, as was done in some of our earlier simulations, but only once every 2000 type I moves.

Surface-confined DNAs

To simulate the transfer of a supercoiled DNA from solution to the surface, a circular chain of 134 subunits is first simulated in the absence of any external potential (), while its linking difference is varied stepwise from 0 to −17 turns, which corresponds to native superhelix density (σ=−0.047). This model DNA with −17 turns is simulated for 35million moves, with single configurations selected at 15, 25, and 35million moves for transfer to the surface. Each of those three configurations is used as the initial configuration for a new simulation in a very weak surface potential (λ=1.407×10⁻⁵). When this simulation has equilibrated 57, its final configuration is used as the first configuration of a second simulation, which uses a larger value of λ. This process is iterated and, as λ is gradually increased in 10 steps from λ=1.407×10⁻⁵ to λ=1, the supercoiled DNA is slowly pressed into relatively flat configurations 57. Simulated transfers were performed using the same 11 values of λ used in our previous work (λ=1.407×10⁻⁵, 5.625×10⁻⁵, 2.250×10⁻⁴, 9.000×10⁻⁴, 3.600×10⁻³, 1.440×10⁻², 4.5918×10⁻², 1.4062×10⁻¹, 2.500×10⁻¹, 5.102×10⁻¹, and 1.000) 57. These 11 λ-values correspond to 11 different z-values at which , which are as follows (in Å): 40,000, 20,000, 10,000, 5000, 2500, 1250, 700, 400, 300, 210, and 150. Of course, for all of these potentials, at −50Å.

For simulations in the presence of , an additional type of Monte Carlo move must be employed. In a type II move, every subunit is translated by the same small random distance along a direction normal to the surface, as described previously 57. When is present, the program alternates between type I and type II moves. The maximum step size for type II (purely translational) Monte Carlo moves was 25Å for λ=1.407×10⁻⁵. As λ was increased, the maximum step size was decreased stepwise to 10Å to ensure that at least 50% of the type II moves were accepted. This entire surface transfer process is repeated for each of the three different initial configurations drawn from the original Monte Carlo simulation in the absence of . Three transfers using the three different values of C (200, 305, and 410fJ·fm) were simulated for each of the two ionic strengths, μ=0.01M and 0.161M, for which AFM images were reported 56. A total of 18 transfers, three for each pair of C- and μ-values, were simulated.

When λ=1.0, the surface potential in Eq. (10) flattens the distribution of duplex DNA centers to a root mean-squared (RMS) thickness in the z-direction of ∼45Å in 0.161M ionic strength, and ∼58Å in 0.01M ionic strength. In contrast, Lyubchenko and Shlyakhtenko reported a vertical image height of 17±3Å for single duplex strands, but performed no systematic AFM study of the height at points of duplex crossings 56. Such low heights imply that the z-distribution of duplex centers has negligible thickness compared to the duplex diameter. This observation raises the question of whether the DNA is so greatly flattened by the surface potential alone, or instead adopts such structures only under the additional force of the AFM tip, which is operated in tapping mode. When the amplitude of tip oscillation is reduced, the apparent height of the DNA is observed to increase by as much as twofold in some cases (personal communication from Y. Lyubchenko, Arizona State University, 2001). This suggests that, under the reduced force associated with the reduced amplitude of the tip oscillation, the DNA may sometimes be encountered and detected at a height significantly above the surface.

No useful information about additional flattening could be obtained by raising the value of λ above 1.0, because the equilibration slowed dramatically for significantly higher values of λ. Nevertheless, two factors suggest that sufficient flattening has been achieved, when λ=1.0. 1), The RMS thickness of the distribution of centers of the surface-confined DNA at λ=1.0 is ∼1/2; the diameter of a normal straight interwound superhelix in 161mM ionic strength, so deformation of the superhelix “structure” is already considerable. 2), It is likely that the AFM tip in tapping mode drives segments of the DNA that lie significantly “above” the surface down “onto” the surface, as noted above. The somewhat smaller than expected vertical height of the flat-lying DNA in AFM images suggests that the applied force is sufficient to compress even the diameter of the duplex and/or the underlying surface. It would require considerably less force to drive high-lying DNA segments, including those atop a plectonemic helix, down onto the surface. For such reasons, we suspect that in the absence of the tip force some of the duplex segments lie significantly above the surface, as is the case in these simulations.

The supercoiling free energy

According to experiments, the change in supercoiling free energy upon varying the linking difference from to is well represented by

(10)

where

is the twist energy parameter and N the number of base pairs. For DNAs in 0.1M ionic strength, both experiments and simulations indicate that

is independent of N and |Δℓ,| provided that N exceeds ∼2500bp 8,11,33,38,55,70,71. Equation (9) was observed experimentally to hold for p30δ DNA in 0.1M NaCl over a range of linking differences from 0 to native (∼−23 turns) 33. Previous simulations with C=200fJ fm and an effective hard-cylinder potential to represent electrostatic interactions were also in reasonable accord with both Eq. (10) and the experimental data, but those with C=300fJ fm exhibited a significant decline in the apparent

with increasing

, and the simulated

-values substantially exceeded the experimental data for all |Δℓ| 60. The more precise simulations described here, using a more realistic subunit interaction potential, yield very similar results, as detailed below.

The supercoiling free energy is defined here as the free energy to increase the linking difference of a closed circular DNA from 0 to Δℓ, and for the model DNAs here is given by

(11)

where

is the average torque exerted on the jth subunit by the j,j+1 spring, which must be overcome by the external agent to change the linking difference by

rad 44. For an isotropic bending potential,

(12)

For this case of an isotropic bending potential, Eq. (11) can be transformed to

(13)

where 〈w〉 is the average writhe of the circular DNA.

was calculated by using both Eqs. (11) and (13) as a check on the mean torque and writhe calculations, and identical results were obtained in all cases. The apparent twist energy parameter is reckoned from

according to

(14)

Determination of via either Eq. (11) or Eq. (13) is referred to as the reversible work method.

Simulation parameters and details

DNAs in solution

Simulations of 170 subunit (4770bp) model solution DNAs were performed for linking differences that spanned the range from 0 to −24 turns, and included, Δℓ=-1, −2, −3, −4, −5, and −6 in the weakly supercoiled regime. The ionic strength was 0.1M (uni-univalent salt) and the temperature was T=310K. The dielectric constant employed, ɛ=78.5, slightly exceeds the prevailing value at 310K, and leads to a slight (0.98-fold) underestimate of the Debye screening parameter, and a very slight overestimate of the repulsive interactions and . This error in is in the opposite direction to that due to neglect of hydration forces, and is expected to be much smaller than statistical errors in the simulations, especially for the small linking differences of interest, where the subunits rarely approach one another closely, even in the complete absence of long-range electrostatic interactions 57,62, unpublished results from our lab).

Surface-confined DNAs

The surface-confined model DNAs comprised 134 subunits (3760bp), and the linking difference was assumed to be −17 turns in 0.161M ionic strength and −15.4 turns in 0.01M ionic strength 57. The simulation temperature was taken to be 298K.

Equilibration

DNAs in solution

Of the various properties calculated for our model DNAs in solution, the radius of gyration is the most slowly varying. For the 170-subunit model DNA with |Δℓ|=2, the autocorrelation function of the radius of gyration () as a function of the number of moves falls below its starting value by a factor of e⁻² at ∼1.1million moves and fluctuates about 0 from ∼3.5million moves up to at least 100million moves; ½ of all moves are normally accepted. Typically, a new simulation for a particular choice of Δℓ and C is initiated by either incrementing Δℓ or altering C from its respective value in a previous simulation, the last configuration of which is the new starting configuration. An equilibration run of 20million moves is performed with the new value of Δℓ or C, before a final run of 340million moves is performed to generate the simulation data. In special cases, discussed below, the data collection runs were extended to 520million moves. We believe that our solution DNAs are generally well-equilibrated at all |Δℓ|-values. However, sampling of the thermally accessible volume in configuration space is significantly sparser for the smallest |Δℓ|-values, for which the DNA is less constrained, and that volume is substantially greater.

Surface-confined DNAs

As the 134-subunit model DNA is transferred from solution (λ=0) to the surface (λ=1.0), its equilibration proceeds in the following way 57. After incrementing λ to a new value, an equilibration simulation of 8million moves is performed. This is followed by two successive simulations, each comprising 4million moves. (A “move” in this context comprises a type I plus a type II move.) Each set of 4million moves is subdivided into five groups (the first 20%, the second 20%, etc.), and the radius of gyration is calculated for each group. The five average values are combined to reckon the total mean () and standard deviation (s(j)) for the jth set (j=1,2) of 4million moves. If lies within the range and also lies with the range , then the simulation is assumed to be satisfactorily equilibrated. If this criterion is not met, then another equilibration run of 8million moves is performed, followed by two more successive sets of 4million moves, which are compared as described above. This process is iterated until the specified condition is met. As λ approaches 1.0, the fraction of accepted Type I moves declines, so the equilibration simulation is extended to 10million moves, and the two data simulations are extended to 5million moves. The same condition must be met before the two most recent blocks of 5million moves are kept as simulated data.

As λ approaches 1.0, the simulation slows down in the sense that the structure changes shape more slowly with increasing move number. Especially, the interchange between unbranched and branched interwound configurations becomes very infrequent and is likely not adequately equilibrated and/or sampled. Such behavior mimics that of pSA509 DNA in the AFM studies, where only rather small changes in molecular shape are observed in the time required to raster the AFM tip over the entire field of view, which includes many DNA molecules. Nonetheless, structural properties, such as the radius of gyration and the writhe, of either kind of separate structure may be more completely equilibrated and adequately sampled.

Results and discussion

Simulated values of E_T and comparison with experiment for solution DNAs

The computed -values from simulations of our 170-subunit model DNA for each different trial value of the torsional rigidity, C, are plotted versus linking difference in Fig. 1. The error bars in the simulation represent one standard deviation of the mean for each point. The much larger standard deviations associated with the smaller linking differences, |Δℓ|≤3, reflect the much greater conformational space that must be explored in those cases to ensure that all of the relevant fluctuations are adequately sampled. Although substantially longer simulations of model DNAs with such small linking differences would be required to determine very accurate values of in this range, the present level of accuracy suffices for our purpose, which is to distinguish between DNAs simulated with C-values that differ by >∼15%.

Display large version of this figure

Figure 1

versus linking difference for 170 subunit chains, as predicted by Monte Carlo simulations. The values of C used in the simulations are (top to bottom, |Δℓ|=5) 410, 305, 269, 231, 200, and 170fJ fm. The bending persistence length of the chain is 50nm, the ionic strength is 100mM, and the temperature is 310K. The shaded area represents the range of measured values of

for p30δ from three separate sets of measurements 33,38,55. The lines simply connect the points from |Δℓ|=24 to |Δℓ|=5, and straight lines are drawn from |Δℓ|=5 to |Δℓ|=0 simply to guide the eye. In reality, these lines are most likely horizontal, but are sloped here to follow the data, which grow increasingly noisy with decreasing magnitude of the linking difference.

The proximity of the simulated -values for C=200 and 231fJ·fm after 340million moves led us to increase the size of the simulations for those two values of C for Δℓ=1, 2, 4, and 6. For each of those eight choices of (C,Δℓ), an additional 180million moves were performed. The average -values for the entire 520million moves in each of these cases appear in Fig. 1. Comparisons between the average -values for the initial simulations of 340million moves (not shown) and those for the entire simulations of 520million moves indicate that the larger number of moves did not increase the separation in between the C=200 and C=231fJ·fm simulations.

The shaded area in Fig. 1 represents the range of the experimental values, =900–1030, that were obtained for p30δ at 310K by three different groups of researchers over a period of more than a decade 33,38,55. All of these experimental -values were obtained by the topoisomer distribution (TD) method, wherein the relative intensities of several (7 to 9) topoisomer bands were simultaneously fitted by Eq. (10) to obtain and . The shaded region extends only to Δℓ=4, because topoisomers with linking differences >4 are either undetectable or present at such low concentration that the relative fluorescence intensities of their bands in the gel cannot be determined with acceptable precision. -values for other DNAs containing M≥2500bp were also obtained via the TD method at 310K 11,70,71. Practically all of the reported values fell in the range, =1000±100, and the majority also fell in the shaded region of Fig. 1.

The simulated -values for C=170fJ·fm clearly fall into the experimentally observed range. Although the -values for C=200 and 231fJ·fm lie just slightly above the shaded region, t

Article Outline document.write('<a id="close_link" onclick="javascript:document.getElementById(\'fulltext_body\').className=\'outline_closed\';document.getElementById(\'open_link\').style.display=\'block\';return false;" href="#" close="Close Article Outline">X</a>');

Article Information

PubMed

Related Articles

Torsional Rigidities of Weakly Strained DNAs

Abstract

Introduction

Variations in secondary structure and torsional rigidity

Methods to measure the torsional rigidity

Reported values of torsional rigidity

Writhe as a probe of the torsional rigidity

Objectives of this work

Theory

Topological and geometrical aspects of supercoiling

The mesoscopic model

Monte Carlo simulation models and methods

Potential functions

Neglect of the hydration potential

Monte Carlo moves and simulation algorithms

DNAs in solution

Surface-confined DNAs

The supercoiling free energy

Simulation parameters and details

DNAs in solution

Surface-confined DNAs

Equilibration

DNAs in solution

Surface-confined DNAs

Results and discussion

Simulated values of ET and comparison with experiment for solution DNAs

Article Outline

Simulated values of E_T and comparison with experiment for solution DNAs