This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.087593v1 91/11/4166    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Fujimoto, B. S. Articles by Schurr, J. M. Search for Related Content PubMed PubMed Citation Articles by Fujimoto, B. S. Articles by Schurr, J. M.

Torsional Rigidities of Weakly Strained DNAs

Department of Chemistry, University of Washington, Seattle, Washington

Correspondence: 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; E-mail: schurr{at}chem.washington.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX: SYSTEMATIC ERRORS IN... ACKNOWLEDGEMENTS REFERENCES

 
Measurements on unstrained linear and weakly strained large (340 bp) circular DNAs yield torsional rigidities in the range C = 170–230 fJ fm. However, larger values, in the range C = 270–420 fJ fm, are typically obtained from measurements on sufficiently small (247 bp) circular DNAs, and values in the range C = 300–450 fJ 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–230 fJ fm are compatible with experimental observations, whereas values in the range C 269 fJ 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

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX: SYSTEMATIC ERRORS IN... ACKNOWLEDGEMENTS REFERENCES

 
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 310 K; 2), overall base composition from 34% to 100% GC; 3), NaCl concentration from 0.01 to 1.0 M; or 4), Mg²⁺concentration from 0 to 40 mM in 0.1 M NaCl are all <15%, so the large spread in reported values must stem from other causes (1–6). All of the larger reported C-values pertain to DNAs that are subject to substantial bending strain in small circles with N 254 bp (7–19), or substantial tensile stresses in single-molecule pulling experiments (20–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–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–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_d (4,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 4 mM ionic strength at 294 K, the values C = 188 fJ fm, = 200 nm, and = 50 nm, were obtained. Robinson and co-workers employed electron paramagnetic resonance spin-label methods to assess , and obtained comparably large values in the range = 150–170 nm for numerous synthetic DNAs in 0.1 M 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–200 nm, and all of the C-values quoted below are obtained by using P_d = 180 nm. 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 (205N250 bp) 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 (150N350 bp) DNAs (7,14–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.4 bp 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 C (17). 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–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–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–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–170 fJ fm, which were obtained for subsequences of 150–170 bp 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.

View this table: [in this window] [in a new window]   TABLE 2  Simulation results for model 134-subunit DNAs in solution ( = 0) and for the final simulations ( = 1) of their transfers to a surface for an ionic strength of 161 mM

All of the C-values determined for circular DNAs with N 247 bp via the FPA and TR methods lie in the range 300–430 fJ fm. Some of the C-values obtained by the CK method likewise exceed 300 fJ fm, although values as low as 220–240 fJ fm were also reported. All of the C-values for equilibrium unstrained linear DNAs of any length and for circular DNAs with N 340 bp lie in the range 150–230 fJ 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–450 fJ 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 450 fJ 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 (4932 bp) at low superhelix densities, || 0.008 (33,38,55), where the strain free energy associated with the coherent (superhelical) deformation is 0.7kT/1000 bp. 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 ( = 50 nm) 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 (3760 bp) 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

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX: SYSTEMATIC ERRORS IN... ACKNOWLEDGEMENTS REFERENCES

 
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–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_j (66). For the subunit rod length, = 9.54 nm, 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.54 nm, was chosen to be slightly <1/5 of a persistence length, P_eq = 50 nm. 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.06 bp. Two different model DNAs are simulated. The first model comprises M = 170 subunits, corresponding to 4770 bp. This is comparable in size to p30 (4932 bp), whose supercoiling free energy was measured in numerous experiments (33,38,55). The second model comprises M = 134 subunits, corresponding to 3760 bp. 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–410 fJ·fm. The bending torque constant is always chosen to yield a total bending persistence length P_eq = 50 nm.

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.20 nm 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.20 nm, separated by 3.18 nm 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 298 K), 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 >2a (57,63,64). This has been done previously for uni-univalent salt concentrations in the range 0.01–1.0 M, and for mixed salt solutions containing 0.10 M NaCl + 0.01 M 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 x 10^–4 when exceeds the cut-off distance. For simulations of the supercoiling free energies of model DNAs in 0.1 M 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.01 M 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 210 fJ fm, the electrostatic potential in Eq. 6 yielded results in good agreement with the measured supercoiling free energy of pBR322 DNA in 0.02 M ionic strength (63,68) and with AFM images of a native supercoiled pSA509 DNA in 0.161 and 0.010 M ionic strength (56,57).

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

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 = 50 nm, in the following manner. Simulations of both 10- and 20-subunit linear chains were performed (10 million 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 x 10^–13, 1.806 x 10^–13, and 1.060 x 10^–13 dyne cm in, respectively, 100 mM, 161 mM, and 10 mM ionic strength yielded P = 50 nm for both the 10- and 20-subunit chains within the simulation errors. The -values in 161 and 10 mM ionic strength are lower than that in 100 mM, because the simulations in 161 and 10 mM include electrostatic interactions between adjacent subunits, whereas the simulations in 100 mM 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.1 nm. 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 z0. The force constant of the harmonic potential (when = 1.0) in the region z > 0 was chosen so that the mean z thickness in 161 mM 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.0 nm (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.4 nm, 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.008 rad 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–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 35 million moves, with single configurations selected at 15, 25, and 35 million 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 x 10^–5). 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 x 10^–5 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 x 10^–5, 5.625 x 10^–5, 2.250 x 10^–4, 9.000 x 10^–4, 3.600 x 10^–3, 1.440 x 10^–2, 4.5918 x 10^–2, 1.4062 x 10^–1, 2.500 x 10^–1, 5.102 x 10^–1, 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 x 10^–5. 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 410 fJ·fm) were simulated for each of the two ionic strengths, µ = 0.01 M and 0.161 M, 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.161 M ionic strength, and 58 Å in 0.01 M 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 the diameter of a normal straight interwound superhelix in 161 mM 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.1 M ionic strength, both experiments and simulations indicate that is independent of N and |,| provided that N exceeds 2500 bp (8,11,33,38,55,70,71). Equation 9 was observed experimentally to hold for p30 DNA in 0.1 M NaCl over a range of linking differences from 0 to native (–23 turns) (33). Previous simulations with C = 200 fJ 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 = 300 fJ 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 (4770 bp) 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.1 M (uni-univalent salt) and the temperature was T = 310 K. The dielectric constant employed, = 78.5, slightly exceeds the prevailing value at 310 K, 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 (3760 bp), and the linking difference was assumed to be –17 turns in 0.161 M ionic strength and –15.4 turns in 0.01 M ionic strength (57). The simulation temperature was taken to be 298 K.

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^–2 at 1.1 million moves and fluctuates about 0 from 3.5 million moves up to at least 100 million 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 20 million moves is performed with the new value of or C, before a final run of 340 million moves is performed to generate the simulation data. In special cases, discussed below, the data collection runs were extended to 520 million 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 8 million moves is performed. This is followed by two successive simulations, each comprising 4 million moves. (A "move" in this context comprises a type I plus a type II move.) Each set of 4 million 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 4 million 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 8 million moves is performed, followed by two more successive sets of 4 million 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 10 million moves, and the two data simulations are extended to 5 million moves. The same condition must be met before the two most recent blocks of 5 million 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

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX: SYSTEMATIC ERRORS IN... ACKNOWLEDGEMENTS REFERENCES

 
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%.

View larger version (24K): [in this window] [in a new window]   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 170 fJ fm. The bending persistence length of the chain is 50 nm, the ionic strength is 100 mM, and the temperature is 310 K. 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 shaded area in Fig. 1 represents the range of the experimental values, = 900–1030, that were obtained for p30 at 310 K 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 2500 bp were also obtained via the TD method at 310 K (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 = 170 fJ·fm clearly fall into the experimentally observed range. Although the -values for C = 200 and 231 fJ·fm lie just slightly above the shaded region, the error bars for all points with = 1, 2, 3, and 4 extend well into that shaded region, so those C-values should also be regarded as consistent with the measured range of . From this comparison between simulated and experimental -values, we infer that the range of "valid" C-values for weakly strained large circular DNAs includes C = 170–230 fJ fm, but does not include higher values in the range C 269 fJ fm! This inferred range of "valid" C-values falls within the range of C-values (150–230 fJ fm) measured for linear and large circular DNAs by FPA and for 340- to 350-base-pair circles by CK (cf. Table 1).

Higher C-values, in the range 310–400 fJ fm were measured for small circular DNAs with N = 181 bp by the FPA method, and C-values in the range 320–410 fJ fm were measured for small circles with 205 N 247 bp by the TR method (cf. Table 1). Most C-values obtained for small circles with N 254 bp by the CK method lie in the range 270–340 fJ fm, although smaller C-values in the range 220–270 fJ fm were also occasionally observed. As noted above, values obtained by the CK method may be only lower bounds, rather than the actual values. Hence, it is conceivable that all circles with N 254 bp exhibit C-values that exceed 270 fJ fm. Values of C >300 fJ·fm were also obtained via different SMP experiments. Nevertheless, the use of any C 269 fJ fm in the simulations yields -values that significantly exceed the measured values.

In conclusion, all values of C 269 fJ fm appear to be incompatible with the -values measured for weakly supercoiled DNAs with || 0.008 at 310 K!

The reported values of C seem to fall largely into two groups, the lower of which (150–230 fJ fm) is typical of unstrained or weakly strained DNAs, and the higher of which (270–450 fJ fm) typifies DNAs that are subject to either coherent bending strain 1.45°/bp or tension 0.1 pN.

Simulated structures of surface-confined DNAs
As reported previously, the transfer of a supercoiled DNA from solution to a surface significantly affects its morphology (57). The surface-confined DNA could be classified as having either an unbranched interwound (UI) or a branched interwound (BI) configuration, usually with three arms of different length. In contrast, the solution DNAs are all interwound, but are not so readily classified into the two simple categories (branched and unbranched) as are the surface-confined DNAs. The solution DNAs often exhibit loops or partial branches that do not properly belong in either category. Equilibration among the UI and BI configurations of the surface-confined DNAs is very slow, so adequate sampling of that equilibrium was probably not achieved in the earlier simulations (57) and is likely not fully achieved here. However, the writhe and number of crossovers in the projection of the DNA onto the xy plane are fairly similar for both UI and BI configurations under a given set of conditions, so incomplete equilibration of the UI BI interchange probably has little effect on these properties.

Quantitative results for simulations of the 134-subunit model in 0.161 M ionic strength are presented in Table 2. Results obtained using C = 200 fJ·fm for both solution ( = 0) and surface-confined ( = 1.0) DNAs show trends that are very similar to those reported previously (57). The twist energies of the present and previous simulations cannot be directly compared, since local fluctuations about the uniform net twist were omitted from the previous simulations, but they are included here. The magnitude of the simulated writhe is 1 turn less in this work than in the earlier simulations. This is apparently due to eliminating errors in the computed writhe, which arise from use of the full subunit length, b = 9.54 nm, in Eq. 3. However, the writhe still becomes more negative by 0.6 turns as the DNA is transferred from solution ( = 0) to the surface ( = 1.0), as found previously. The bending, electrostatic, and external potential energies and also in Table 2 are all very similar to their respective values obtained in the previous simulations. A decrease in the internal energy, as the DNA is transferred from the solution to the surface, is also seen here, and is of similar size to the decrease seen in the previous simulations. The current and previous simulations of 134-subunit solution DNAs yield similar values for and also for , but modest (10%) differences appear in the case of surface-confined DNAs. Since and for a surface configuration are determined primarily by whether that configuration is UI or BI, and the UI BI equilibrium is likely not adequately sampled in either the previous or present simulations, such modest differences between the two simulations in regard to these quantities are not a source of concern. Despite some differences in details, results of the current simulations using C = 200 fJ fm for solution and surface-confined DNAs in 0.161 M ionic strength generally confirm the conclusions implied by the earlier simulations (57).

The simulated surface-confined DNAs in 0.161 M ionic strength with C = 200 fJ fm exhibit configurations that qualitatively resemble those observed in AFM studies (56). Although it is not possible to unambiguously count the number of crossovers in those AFM images, plausible estimates are in the range 10–12. The number of crossovers exhibited by the simulated surface-confined ( = 1.0) DNAs with C = 200 fJ fm lies in the range 12–13 and the writhe is –11.7, all in good agreement with the AFM images. However, the number of crossovers observed increases to 14.5 for C = 305 fJ fm and to 15.5 for C = 410 fJ fm. The magnitude of the writhe also increases to –13.5 for C = 305 fJ fm and to –14.5 for C = 410 fJ fm. The agreement between the simulated structures and AFM images in 0.161 M ionic strength is significantly poorer for these larger values of C. This provides another indication that large values of C in the range 300–450 fJ fm do not apply to weakly strained DNAs.

The simulation results for solution 134-subunit model DNAs in 0.01 M ionic strength are presented in Table 3. The data obtained using C = 200 fJ fm for solution ( = 0) and surface-confined ( = 1.0) DNAs show trends that are similar to those reported previously (57). The bending, electrostatic, and external energies and also in Table 3 are all very similar to those obtained by the previous simulations. In both this work and the previous study, very little change is observed in the internal energy, as the DNA is transferred fro