This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.091280v1 91/12/4507    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 HighWire Citing Articles via Google Scholar Google Scholar Articles by Mazur, A. K. Search for Related Content PubMed PubMed Citation Articles by Mazur, A. K.

Evaluation of Elastic Properties of Atomistic DNA Models

Centre National de la Recherche Scientifique, UPR9080, Institut de Biologie Physico-Chimique, Paris, France

Correspondence: Address reprint requests to A. K. Mazur, Centre National de la Recherche Scientifique, UPR9080, Institut de Biologie Physico-Chimique, 13 rue Pierre et Marie Curie, Paris, 75005, France. E-mail: alexey{at}ibpc.fr.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
A number of intriguing aspects in dynamics of double-helical DNA is related to the coupling between its macroscopic and microscopic states. A link between the elastic properties of long DNA chains and their atom-level dynamics can be established by comparing the worm-like chain model of polymer DNA with the conformational ensembles produced by molecular dynamics simulations. This problem is complicated by the complexity of the DNA structure, the small size of DNA fragments, and relatively short trajectory durations accessible in computer simulations of microscopic DNA dynamics. A careful study of all these aspects has been performed by using longer DNA fragments and increased durations of MD trajectories as compared to earlier such investigations. Special attention is paid to the necessary conditions and criteria of time convergence, and the possibility to increase the sampling by using constrained DNA models and simplified simulation conditions. It is found that dynamics of 25-mer duplexes with regular sequences agrees well with the worm-like chain theory and that accurate evaluation of DNA elastic parameters requires at least two turns of the double helix and 20-ns duration of trajectories. Bond length and bond-angle constraints affect the estimates within numerical errors. In contrast, simplified treatment of solvation can strongly change the observed elastic parameters of DNA. The elastic parameters evaluated for AT- and GC-alternating duplexes reasonably agree with experimental data and suggest that, in different basepair sequences, the torsional and stretching elasticities vary stronger than the bending stiffness.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
The flexibility and internal dynamics of DNA are known to play an important role in its biological function. Biologically relevant dynamics of DNA spans a broad range of length scales starting from a subnanometer level of single basepairs and continues to macroscopic lengths where the double helix can be considered as a continuous flexible rod. All these motions are evidently coupled. For instance, the macroscopic elastic properties of DNA depend upon the basepair sequence, while external torsional stress applied to long chains affects the local rates of basepair opening. The corresponding molecular mechanisms are only partially understood and disclosing them represents a challenging task.

Long DNA double-helix behaves as a continuous elastic rod with harmonic bending, torsional, and stretching deformability. Its equilibrium shape is described by the worm-like chain (WLC) model (1,2), whereas the torsional and stretching fluctuations can be treated with the standard formalism of the classical statistical mechanics. This model accurately describes experimental data by using only a few adjustable parameters, which have been measured since the 1960s with progressively improved accuracy (3). More recently, single molecule nanomanipulations (4–10) became the major source of information on the DNA elasticity, and the WLC theory is successfully used for extracting information from such data (11–18). The validity of the WLC description of long DNA is corroborated by the success of Monte Carlo simulations of the WLC model with experimentally measured parameters. In many cases, such simulations demonstrated remarkable agreement with experiment and made predictions confirmed later (19–22).

Owing to important methodological advances made approximately 10 years ago (23–25), the local atom-level DNA properties are now reasonably well reproduced in classical molecular dynamics (MD) simulations of small duplexes in explicit aqueous environment (26,27). MD simulations are based upon empirical force fields that are parameterized by using experimental data for small molecules as well as quantum mechanics calculations (26,28). The experimental information about DNA elasticity is not used for parameterization; therefore, comparison of atom-level MD models with the WLC theory present significant fundamental interest. To establish this link between the microscopic MD and the WLC theory, we need to extract the effective WLC parameters from atom motions in short DNA fragments observed during relatively short times. This problem is not simple, because the DNA structure is complex, the duration of MD trajectories is limited by computer resources, and because we cannot exclude that the WLC theory works well for DNA only starting from prohibitively long fragments.

The first attempt to evaluate the elastic parameters of DNA from atom-level data was made by Olson et al. (29), who studied statistics of fluctuations in x-ray DNA structures assuming that the perturbations due to the crystal environment and proteins bound to DNA are equivalent to the heat bath effect. MD simulations of DNA were first analyzed from the same perspective by Bruant et al. (30), and later by Lankas et al. (31), who specifically targeted the sequence-dependent DNA elasticity. The obtained estimates of the WLC parameters reasonably agreed with other data. These studies, however, did not consider important issues concerning the overall agreement of DNA dynamics with the WLC theory, statistical convergence, and the accuracy of the final estimates. These aspects are crucial for probing relatively small modulations of DNA elasticity, for instance, due to small external stress. Here we study in detail the problems involved in the interpretation of MD trajectories of DNA in terms of the WLC theory by using longer DNA fragments and increased durations of MD trajectories. Special attention is paid to the necessary conditions and criteria of time convergence, and the possibility to increase the sampling by using constrained DNA models and simplified simulation conditions. It is found that dynamics of 25-mer duplexes with regular sequences agrees well with the WLC theory and that accurate evaluation of DNA elastic parameters requires at least two turns of the double helix and 20-ns duration of trajectories. Bond length and bond-angle constraints affect the estimates within numerical errors. In contrast, simplified treatment of solvation can strongly change the observed elastic parameters of DNA. The elastic parameters evaluated for AT- and GC-alternating duplexes suggest that, in different basepair sequences, the torsional and stretching elasticities vary stronger than the bending stiffness.

	THEORETICAL BACKGROUND AND METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
The persistence lengths of the WLC model
The persistence length (PL) is a convenient measure of statistical flexibility of polymer chains that can be estimated experimentally. The advantage of the WLC model is that all its properties, including the local ones, are fully determined by the bending, torsional, and stretching PLs, which gives a theoretical possibility of evaluating the DNA elasticity from the MD of relatively short chains.

The standard definition of the bending PL is as follows. Consider a chain fragment of the contour length L placed in a thermal bath. We connect the two chain-ends by a vector and compute its projection upon the local chain direction measured at its origin. The Boltzmann average of this projection characterizes the internal molecular flexibility. The bending PL is defined as the limit of this projection with L .

We assume that the polymer can be approximated by a continuous flexible rod and try to compute its bending PL from local elastic properties. Statistical mechanics of flexible rods was first analyzed by Bresler and Frenkel in the 1930th (2). For an isotropic rod of length L, its bending energy in the first approximation is

(1)

where (l) is the bend angle of fragment (0, l) and A is a constant. The Boltzmann average of the bending angle cos is computed by using the general theory of fluctuations (2)

(2)

and it follows that

(3)

Equation 2 can be rewritten in a linear form suitable for data processing

(4)

The D_a values estimated for DNA fragments of increasing lengths should grow linearly and the rate of this growth gives an estimate of the bending PL. Several similar functions will be introduced below and, for convenience, they are all referred to as "deviations." Subscript a indicates that D_a and A_a are obtained from average bend angles.

The WLC bending PL can be also estimated from the average end-to-end distance. Vector R joining the chain-ends is

(5)

where t(l) is a unit tangential vector. The average square distance R² is obtained by straightforward integration using Eq. 5 together with Eq. 2, which gives (2)

(6)

For DNA fragments used in MD, we always have L << A, therefore Eq. 6 is approximated as

(7)

which gives a linear relationship similar to Eq. 4,

(8)

The D_r value can be computed from MD trajectories and used for evaluating A_r.

If the DNA molecule behaves as an elastic rod, the values of A_a and A_r must be equal. Even though Eqs. 4 and 8 are not independent, comparison of A_a and A_r appears useful in practice. A serious difficulty inherent in this type of computation is in the unknown rate of convergence. Both cos and R² are positive and there is no other parameter that would have a standard limiting value. Therefore, it is not easy to judge how complete the sampling is at a given stage of simulation. In the course of MD, the values of A_a and A_r can approach the same limiting value differently. For instance, local bending dynamics may be fast with respect to translational motion of the chain ends. In this case the rapid bending dynamics should be dominated by correlated fluctuations that leave the chain-ends fixed. As a result, A_a and A_r measured for finite trajectories may differ significantly. Comparison of A_a and A_r provides a valuable additional check of convergence and consistency of the data. Moreover, it is desirable to select a procedure that allows obtaining accurate estimates from possibly short MD trajectories. From this prospect, an appealing possibility for testing resides in using fluctuations of R² and cos rather than their averages. Here, we test this idea for cos. Its standard deviation cos is

(9)

The WLC theory does not offer an analytical relationship for cos(L), but the angle distributions obtained by MC simulations of a genuine WLC model reasonably well agree with

(10)

For WLC chains of 10-nm length and the bending PL of 20–60 nm, this equation gives an estimate of the persistence length with an error 2%.

An additional complication in the DNA bending dynamics lies in the possibility of static curvature. Certain DNA sequences are known to induce stable bends and, probably, many more such sequences exist but are not known. Moreover, static and metastable bending components are indistinguishable in the course of a finite MD simulation; therefore it is important to have a monitoring procedure to check whether the DNA bending observed in an MD trajectory contains such components. Here we apply the following approach. A reference Cartesian frame is placed at the chain origin, with the unit vector e_z looking along the helical axis, and vectors e_x and e_y being perpendicular to it. The unit vector of the helical axis at the opposite chain-end is denoted h. For true WLC dynamics we should have h_x = h_y = 0, and h_z = cos > 0. If h_x 0, or h_y 0, the DNA fragment under consideration is curved, on average. Let us define a unit vector with components

(11)

Vector gives an approximate direction of a helical axis of the average DNA conformation. Now we can measure the angle between vectors h and and compute its average cosine as

(12)

Just by analogy with Eq. 4, we can write

(13)

and

(14)

The last two equations do not have a rigorous justification. Intuitively, one can expect that the left-hand values in Eqs. 13 and 14 both grow with L; therefore, these formulae can be considered as first approximations of more complex functions. Below, parameters A_d and A_s are referred to as dynamic and static PL, respectively. In the course of MD of a true WLC model one should observe that A_s and A_d A_a; therefore, the closeness of A_d and A_s to the limiting values serve as additional criteria of convergence. The distinction between the static and dynamic PL is usually used in the framework of the wedge theory of DNA curvature (32–34). Here, we use the same terms in the context of dynamics of a single DNA molecule; therefore, the above definition of A_s and A_d is adapted for processing MD data. The physical meaning of A_s and A_d is not the same as in the wedge theory. Notably, no simple additivity rules exist that connect A_s and A_d with the total PL.

The classical WLC model can be extended to take into account fluctuations of the overall winding angle and the chain length. To this end, the total elastic energy is written as

(15)

(16)

where (L) is the total winding angle while ₀ = and L₀ = L denote minimum energy values. Parameter C in Eq. 16 is called torsional PL (by analogy with A in Eq. 1) and Y_f is Young's stretching modulus. Equation 15 is the first approximation of the elastic energy, with the possible coupling between bending, torsional, and stretching motions neglected. In contrast to the bending-angle , the winding and the length of consecutive DNA stretches are simply additive; therefore, expressions similar to Eq. 1 for torsional and stretching energies readily give Eq. 16. The Boltzmann averaging of Eq. 16 gives expressions similar to Eq. 4,

(17)

and

(18)

used for estimating the values of C and Y_f from MD results. The stretching PL, B, is usually defined as (15)

(19)

where l₀/2 is a scaling coefficient, with l₀ = 3.4 nm (the length of one helical turn in B-DNA).

Despite the apparent similarity of Eqs. 4, 17, and 18, the bending PL is clearly distinguished from the other two because it is determined from the average bending angle rather than the standard deviation. This difference is due to the nontrivial additivity of bending angles along the chain that requires a special treatment (the WLC model), and it does not depend upon the chain length. For instance, if both L and are very small, the bending energy of a flexible rod can be approximated as

(20)

and expansion of both sides of Eq. 2 gives

(21)

If we replace ² on the left in Eq. 21 by (by analogy with Eq. 17), the bending PL A on the right must be replaced by a larger value because 0 for the WLC model. Parameter also characterizes the bending rigidity, but it can significantly differ from A even with L 0. Therefore, evaluation of the bending PL from standard deviations rather than averages of bending angles is possible only by approximate empirical relationships like Eq. 10.

Construction of helical axes
To compare an ensemble of DNA conformations with the WLC model, every DNA structure should be replaced by an oriented space line. This issue would be unimportant if one could run MD simulations for very long double helices, which is hardly possible even in long-term prospect. Following earlier similar studies (30,31), we decided to employ the Curves algorithm by Lavery and Sklenar (35) to construct the DNA helical axes and measure bending, torsional, and stretching fluctuations. To this end, the algorithm was reimplemented and verified against the original Curves program for selected DNA conformations and also by processing intervals of MD trajectories by the two implementations in parallel. The Curves algorithm starts by constructing a local Cartesian frame at every base, with the frame orientation depending only upon the coordinates of N1/N9 atoms and the base plane. The optimal helical axis consists of a new sequence of Cartesian frames, one for each basepair, that we call axis frames. They are positioned in space by numerical optimization of a target function that can be written as

(22)

The exact definition of the four terms denoted A1, A2, B1, and B2 are given in the original article (35). Qualitatively, conditions A1 and A2 require that orientations of base coordinate frames with respect to the axis frames of the same level were possibly similar for neighboring levels. Conditions B1 and B2 do not consider the DNA structure and only require that the optimal helical axis be locally close to a straight line. Essentially, the axis is modeled as a flexible rod with elastic properties determined by conditions B1 and B2, and the numerical optimization of function (22) is used to fit it to the DNA structure. In the original version of the algorithm, factor Q_c was absent (Q_c = 1). Here it was added to vary the bending elasticity of the axis rod.

MD simulations
We study dynamics of 25-mer double-helical DNA fragments with AT- and GC-alternating sequences (AT25 and GC25, respectively). AT25 is modeled in three simulation regimes that differ by hydration conditions as well as the number of degrees of freedom in the DNA duplexes. In all cases, the AMBER98 force-field parameters (23,36) were used with the rigid TIP3P water model. Molecular dynamics simulations of rigid and partially fixed molecules were carried out with the ICMD method (37,38) adapted for DNA (39,40). In the first regime (trajectory AT25a, 16 ns, time-step 0.002 ps), the duplex was modeled with all degrees of freedom in a rectangular water box with a neutralizing number of sodium ions. In the second regime (trajectory AT25b, 28 ns, time-step 0.01 ps), the hydration conditions were the same as in AT25a, but the duplex was modeled with all backbone torsion degrees of freedom, free bond angles in sugar rings, and rigid bases and phosphate groups. In the third regime (trajectory AT25c, 120 ns, time-step 0.01 ps), the minimal B-DNA model was used, with semi-implicit treatment of solvent as described earlier (39,41,42). Duplex GC25 was modeled in conditions corresponding to AT25b (trajectory GC25 20 ns, time-step 0.01 ps). Trajectories AT25a, AT25b, and GC25 were run with periodical boundaries, in NVT ensemble conditions with water density at 0.997. The electrostatic interactions were treated by the SPME method (25), with the common values of Ewald parameters; that is, 9 Å truncation for the real-space sum and ß 0.35. In all cases, the fiber canonical B-DNA model (43) was used as the starting state. The starting states of trajectories AT25a and AT25b were prepared as follows. The DNA fragment was immersed in a rectangular water box 110 x 46 x 46 Å, with a higher water density of 1.04. The box was neutralized by placing Na⁺ ions at random water positions at least 5 Å from the solute. The system was energy-minimized and dynamics were initiated with the Maxwell distribution of generalized momenta at low temperature. The system was next slowly heated to 293 K and equilibrated during 0.6 ns. After that, the water density was adjusted to 0.997 by removing the necessary number of water molecules selected randomly at least 5 Å from DNA and ions, and the simulations were continued with NVT ensemble conditions. In all simulations, the temperature was maintained by the Berendsen algorithm (44) with a relaxation time of 10 ps. In all trajectories the DNA structures were saved every 2.5 ps. These conformations were later processed with the Curves algorithm to accumulate the data for the statistical analysis according to the WLC model. To increase the sampling, all possible internal fragments of the 25-mer duplexes were considered; for example, averaging for dimers involved 24 times more conformations than for 25-mer.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Artifacts of elasticity of helical axis
Fig. 1 displays the apparent chain-length dependence of different bending deviations for trajectory AT25b when the helical axes of DNA structures are computed with the original parameterization of the Curves algorithm; that is, with Q_c = 1 in Eq. 22. According to the WLC theory described above, all plots in this figure except D_d and D_s must be linear with the same slope. Instead, they all display strong deviations. The most striking are evident regular oscillations of D_a with a period close to that of DNA. Their origin becomes evident if we compare the traces of D_s and D_d. It is apparent that the helical axis exhibits static bending that regularly changes with the chain length. The corresponding DNA conformations, however, are quite straight and show no deformations that could have been responsible for such effects. The apparent static bending reaches maxima for DNA lengths of 0.5 and 1.5 helical turns, but is close to zero for full turns, which means that the axis is helical itself, with the winding period identical to that of the double helix.

View larger version (19K): [in this window] [in a new window]   FIGURE 1  Chain-length dependences of bending deviations measured by different methods. Results of trajectory AT25b were processed by the Curves algorithm with Qc = 1.

In the earlier studies based upon the Curves analysis of DNA structures these artifacts were not noticed (30,31). To check if they are specific to our implementation of Curves the same data were processed with the original program, but the results did not change. Moreover, this effect is not specific to the Curves algorithm. Similar features persist with the alternative procedure (41,45) that also uses numerical optimization, but finds a common axis of coaxial cylindrical surfaces passing through different atoms of DNA backbone. Therefore, the origin of the artifact helicity is in the very principle of fitting an elastic axis rod to a DNA structure. This may be a problem of local minima because the initial axis trace usually has some degree of helicity. It is also possible that the artificial helicity appears spontaneously during fitting because this is the lowest-energy deformation of a straight rod. In any case, evaluation of the WLC parameters of DNA from the data shown Fig. 1 is hardly justified. Fortunately, it appears that these artifacts can be reduced by adjusting the stiffness of the axis rod.

Fig. 2, a and b, shows how the artificial stiffness of the axis rod affects the measured bending deviations. The original Curves parameterization (Q_c = 1) gives a relatively soft axis that, for short DNA fragments, bends excessively in response to local conformational fluctuations. With the stiffness increased, the artifact oscillations in the D_a traces are reduced as well as the short length burst in the D_r traces, and for Q_c = 10 the chain length dependences demonstrate good agreement with the WLC theory. With the stiffness further increased, the initial artifacts change their sign. Notably, the bending looks smaller for short DNA fragments as well as near chain-ends, because in these cases, the fitted axis is less loaded by the DNA structure and it is easier for it to straighten.

View larger version (12K): [in this window] [in a new window]   FIGURE 2  Chain-length dependences of deviations Da, Dr, Dc, and Db with increased stiffness of the axis rod in Curves. The multiple plots in all plates were obtained by increasing Qc in Eq. 22, Qc = 1, 5, 10, 15, 20, 25, 30, 35, and 40. Traces for Qc = 1 and 10 are shown by and , respectively. Other traces are shown by •. The initial data were the same as in Fig. 1.

The results in Fig. 2 suggest that reasonably good estimates of DNA elastic parameters can be obtained from the slopes of linear fits of D(L) obtained with increased artificial stiffness of the helical axis. It is understood that as Q_c in Eq. 22 is gradually increased, the helical axis becomes more and more straight, and the question arises as to whether the elastic parameters we measure refer to DNA rather than to the artificial stiffness of the axis rod. Fig. 3 shows, however, the contribution of the artificial stiffness is not large and can be accounted for. The corresponding persistence lengths were obtained by linear regression analysis of D(L) dependences evaluated for trajectory AT25b. It is seen that, as expected, different estimates of the bending PL in Fig. 3 a all grow with Q_c. With Q_c > 5 this growth is small, so that an order-of-magnitude increase in Q_c gives only 10% increase in A_a. Moreover, the plots in Fig. 3 a are nearly linear when Q_c > 10; therefore, they can be extrapolated to Q_c = 0 to give asymptotic values. Interestingly, all estimates of the bending PL linearly extrapolated to Q_c 0 from Q_c > 10 converge to 80 nm.

View larger version (7K): [in this window] [in a new window]   FIGURE 3  Dependence of the estimated DNA elastic parameters upon the artificial stiffness of the helical axis. (a) Bending PL. Aa, ; Ad, •; Ar, ; Af, ; and A'a, +. (b) Torsional PL. (c) Stretching PL.

Fig. 2 indicates that the artificial stiffness Q_c = 10 represents a reasonable compromise because, on the one hand, it scales down the artifacts displayed in Fig. 1, and, on the other hand, it is relatively low and gives estimates of elastic parameters close to the Q_c = 0 asymptotic intercepts.

Chain-end effects
In analysis of DNA structures, a few terminal basepairs are commonly excluded. This is reasonable because the DNA ends may be physically more flexible and also because fitting of helical axis near chain ends is less reliable. Fig. 4 shows how the estimated persistence length changes when terminal basepairs are excluded from analysis. With 10 basepairs excluded, only the central 5-mer is considered, and so forth. The four different estimates of the bending PL shown in Fig. 4 a were obtained by linear regression analysis of the corresponding D(L) plots, with the helical stiffness Q_c = 10. It is seen that the end effects are small. A descending trend in A_a and A_r is larger for short trajectories, but it reduces with time. Certainly, long chains need longer time for sampling configurations corresponding to a given persistence length. However, a much more noticeable trend observed in Fig. 4 a consists of the divergence of the four estimates of the bending PL, which becomes dramatic when the length of the remaining part of DNA approaches one helical turn. With the reduced DNA length, the amplitudes of bending fluctuations are reduced and the above-discussed artifacts of Curves become more pronounced. In addition, the sampling is reduced due to the smaller number of DNA conformations analyzed and also due to a smaller number of points for the linear regression analysis. As seen in Fig. 4, b and c, these factors affect the estimates of the torsional and stretching PLs only when <1 helical turn remains for analysis. In contrast to bending, however, here the expected strong chain-end effects are evident. The DNA ends look softer than its inner part as regards both winding and stretching. As noted already, factor Q_c directly affects only the bending stiffness. With Q_c = 10 the chain-end effects seem to be significantly reduced for bending, whereas for stretching and torsional deformations they are still distinguishable.

View larger version (12K): [in this window] [in a new window]   FIGURE 4  Bending, torsional, and stretching PL measured for the central part of the 25-mer fragment. The number of basepairs excluded from both ends is shown on the x axis. (a) Different estimates of the bending PL. Aa, ; Ad, •; Ar, ; and Af, . (b) Torsional PL. (c) Stretching PL.

1. The helical stiffness Q_c = 10 is used as optimal, but the analysis is repeated for a series of Q_c values to check its effect upon the results.
   
2. Three basepairs from both ends are excluded from analysis in correspondence with the earlier studies (31).
   
3. Based upon the results shown in Fig. 2, the persistence lengths are estimated by linear regression analysis of the D(L) plots over L > 6 for bending-angle deviations (except D_s(L)), L > 12 for D_r, and L > 3 for torsional and stretching deviations. The D_s(L) plots remain oscillating with high Q_c values; therefore, the static PL is estimated by simple averaging of A_s(L). Every point on D(L) plots is weighted by the corresponding number of sampled DNA conformations.
   

Time-convergence
Perhaps the most difficult issue concerns the duration of MD trajectories necessary for reliable estimates of the DNA elastic parameters. Recent experimental reports suggest that the conformational relaxation times can be hopelessly long even in short DNA fragments (46–48). In contrast, certain structural parameters of DNA seem to converge very rapidly during MD; for example, the average helical parameters and the average fraction of B_II conformations (49). Fig. 5 shows the time-dependence of the measured DNA elastic parameters for trajectory AT25b. When MD trajectories start from canonical structures, DNA commonly exhibits a conformational drift at the beginning. The duration of this initial phase cannot be known a priori because it depends upon the specific conformational parameter considered and because it is difficult to distinguish between slow drift and a slow equilibrium fluctuation. If the data with such drift are interpreted as fluctuations, statistical parameters measured during the first nanoseconds can be strongly biased, and it takes a long time before they reach their true values. For torsional and stretching PL, this effect is always seen as very slow growth, because these are computed from standard deviations of the winding angle and the length of DNA, respectively. The corresponding manifestations in the measured bending PL may vary.

View larger version (17K): [in this window] [in a new window]   FIGURE 5  Time-convergence of the measured DNA elastic parameters in trajectory AT25b. (a,d) Bending PL with normal and inverted time directions, respectively. Aa, ; Ad, •; Ar, ; Af, ; and A'a, +. (b,e) Torsional PL with normal and inverted time directions, respectively. (c,f) Stretching PL with normal and inverted time directions, respectively.

These figures show also that different estimates of the bending PL become closer with time, which is very important because this feature can be used as a unique necessary criterion of convergence. The rate of convergence is similar for all five estimates of the bending PL so that none of them generally give better estimates for short trajectories. Nevertheless, all of them are useful because, as seen in Fig. 5, a and d, only comparison of all these values gives a correct appreciation of convergence. The discrepancy between A_d and the total PL falls down initially, but eventually stabilizes at a non-zero level. Simultaneously, the static PL, A_s, grows and stabilizes around 10^3.5 (data not shown). Even though A_s 10^3.5 seems large, it can result in a nonnegligible deviation of A_d from the total PL due to specific combinations of bending directions. This residual divergence between A_d and the total PL is mainly due to insufficient sampling of curved conformations with alternative bending directions, although it also involves a contribution from artifacts in the construction of the helical axis that are reduced, but not completely suppressed, with Q_c = 10.

Fig. 5 shows that all the measured elastic parameters are characterized by similar relaxation times. The convergence of the torsional and stretching PLs was expected to be somewhat faster because the corresponding fluctuations require relatively small atom displacements and they should be also much less hindered by water. Moreover, in the course of MD simulations, the average twist and rise normally reach stationary values after 1–2 ns. Surprisingly, Fig. 5 clearly shows that the convergence is similar for all elastic parameters. This may indicate that all these motions are significantly coupled, which should have been taken into account in Eq. 15, as earlier suggested by some groups (31). Accurate analysis of this issue requires further studies.

Alternative DNA models
Figs. 6 and 7 show the time dependence of the measured elastic parameters for the two alternative representations of the AT-alternating 25-mer duplex. The three models we consider strongly differ in performance. Model AT25a (Fig. 6) is most detailed, but also most computationally demanding. Model AT25b (Fig. 5) is five times less expensive computationally and may represent a reasonable compromise for practical computations. Model AT25c (Fig. 7) is computationally much faster than the former two, but its approximations may affect DNA dynamics. Trajectory AT25a was continued to 16 ns, which, according to Fig. 5, is the necessary minimum for convergence. Comparison of Figs. 5 and 6 suggests that the standard geometry fixations used in model AT25b do not significantly affect the elastic properties of DNA. The only significant difference concerns the lower apparent torsional PL for AT25b, but it appears to be due to a more pronounced chain-end effect. The behavior of the torsional PL for AT25a is qualitatively similar to that shown in Fig. 4, but the plateau at 150 nm is reached only when six basepairs from both ends are excluded. We conclude, therefore, the elastic parameters can be measured with reasonable accuracy by using the partially fixed DNA model. MD trajectories for the third model, AT25c, can be run for much longer durations, but as seen in Fig. 7, the absolute rate of convergence in this case is also the lowest. This was expected since the distance scaling of electrostatic interactions is known to increase artificially the strength of all polar contacts including hydrogen bonds (50). The same effect can also be responsible for the generally higher stiffness of the AT25c model compared to the other two.

View larger version (14K): [in this window] [in a new window]   FIGURE 6  Time-convergence of the measured DNA elastic parameters in trajectory AT25a. The same notations as in Fig. 5.

View larger version (17K): [in this window] [in a new window]   FIGURE 7  Time-convergence of the measured DNA elastic parameters in trajectory AT25c. The same notations as in Fig. 5.

View larger version (16K): [in this window] [in a new window]   FIGURE 8  Time-convergence of the measured DNA elastic parameters in trajectory GC25. The same notations as in Fig. 5.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

Much of our efforts were spent on sorting out different factors that hinder interpretation of the atomic-level DNA dynamics in terms of the WLC theory. The major such factor is related with the necessary reduction of a detailed DNA conformation to an oriented axis rod (optimal helical axis). The Curves algorithm was proposed approximately 20 years ago (35), and since then became a standard instrument for treating this problem (30,31). Here it was found that the helical axes computed by Curves generate systematic artifact deviations from the WLC theory that can be reduced by adjusting the internal parameters of the algorithm. This does not mean that the earlier usage of Curves was biased by artifacts. Real double-helical DNA conformations do not have strict helical symmetry; therefore, the concept of optimal helical axis is relative and valid for a certain choice of criteria and parameters. Any such axis is good for comparisons if it is chosen as a conventional standard. The axes computed with the standard Curves parameterization deviate from the WLC theory, but those computed with increased axis stiffness may be inadequate for other purposes.

One might suggest that the Curves algorithm should be corrected to suppress the observed deviations from the WLC theory without increasing the artificial axis stiffness. In our view, this can only hide the problem and complicate interpretation of the results. The oscillating character of angular deviations in Fig. 1 varies with the basepair sequence of the DNA fragment, and sometimes it is not that obvious. In contrast, our earlier experience shows that, with the standard Curves algorithm, the A_r trace always has a concave shape. Such pattern can be misinterpreted as if short DNA fragments experience strong local bends that are anticorrelated along the sequence and cancel out in longer stretches. The oscillations shown in Fig. 1 in fact represent a fortunate feature that immediately reveals their nonphysical nature. Any improved axis-fitting algorithm will still be prone to local artifacts because the problem is inherent in the very fitting principle. For instance, the Curves algorithm assumes that base orientations with respect to the optimal helical axis change smoothly along DNA. This certainly contradicts the zigzag character of DNA with alternating sequences and perhaps results in the regular oscillations in Fig. 1. Unfortunately, it is hardly possible to propose a simple alternative algorithm free from such defects. We believe that accurate evaluation of DNA elasticity in MD simulations must, in any case, include systematic variation of the parameters employed by the axis-fitting algorithm and extrapolation to zero artificial stiffness. From this perspective, the Curves method is advantageous because it is rather rapid.

Earlier estimates of elastic parameters of MD models of DNA were made with the same force field and similar simulation conditions (30, 31), but with double helices of 15–16 basepairs and the duration of trajectories of 1–5 ns. We believe that the relatively small DNA length in these earlier studies explains why the Curves effect discussed above was not noticed. Instead, the apparent chain-length dependences of elastic parameters were interpreted as physical effects (31). In addition, unlike these earlier reports, we did not take into account the coupling between different elastic deformations because we were not sure if the accuracy of our data justifies a more sophisticated analysis. It is not surprising, therefore, that some of the values reported in Table 1 significantly differ from earlier estimates. Nevertheless, the qualitative trends are well reproduced—notably, the AT-alternating duplex exhibits higher torsional, but lower stretching stiffness than the GC-alternating DNA, which was also observed by Lankas et al. (31).

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Alexander Vologodskii has generously provided MC data for comparisons and fitting of Eq. 10. He is also gratefully acknowledged for many fruitful discussions during the course of these studies.

Submitted on June 14, 2006; accepted for publication September 13, 2006.

	REFERENCES

TOP ABSTRACT INTRODUCTION THEORETICAL BACKGROUND AND... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
1. Cantor, C. R., and P. R. Schimmel. 1980. Biophysical Chemistry, Part III: The Behavior of Biological Macromolecules. W. H. Freeman, San Francisco, CA.

2. Landau, L. D., and E. M. Lifshitz. 1976. Statistical Physics, Part 1. Nauka, Moscow, Russia.

3. Hagerman, P. J. 1988. Flexibility of DNA. Annu. Rev. Biophys. Biophys. Chem. 17:265–286.[CrossRef][Medline]

4. Cluzel, P., A. Lebrun, C. Heller, R. Lavery, J. L. Viovy, D. Chatenay, and F. Caron. 1996. DNA: an extensible molecule. Science. 271:792–794.[Abstract]

5. Lavery, R., A. Lebrun, J.-F. Allemand, D. Bensimon, and V. Croquette. 2002. Structure and mechanisms of single biomolecules: experiment and simulation. J. Phys. Condens. Matter. 14:R383–R414.[CrossRef]

6. Perkins, T. T., S. R. Quake, D. E. Smith, and S. Chu. 1994. Relaxation of a single DNA molecule observed by optical microscopy. Science. 264:822–826.[Abstract/Free Full Text]

7. Smith, S. B., Y. Cui, and C. Bustamante. 1996. Overstretching B-DNA: the elastic response of individual double-stranded and single-stranded DNA molecules. Science. 271:795–799.[Abstract]

8. Smith, S. B., L. Finzi, and C. Bustamante. 1992. Direct mechanical measurements of the elasticity of single DNA molecules by using magnetic beads. Science. 258:1122–1126.[Abstract/Free Full Text]

9. Strick, T. R., J. F. Allemand, D. Bensimon, A. Bensimon, and V. Croquette. 1996. The elasticity of a single supercoiled DNA molecule. Science. 271:1835–1837.[Abstract]

10. Strick, T. R., M.-N. Dessignes, G. Charvin, N. H. Dekker, J.-F. Allemand, D. Bensimon, and V. Croquette. 2003. Stretching of macromolecules and proteins. Rep. Prog. Phys. 66:1–45.[CrossRef]

11. Baumann, C. G., S. B. Smith, V. A. Bloomfield, and C. Bustamante. 1997. Ionic effects on the elasticity of single DNA molecules. Proc. Natl. Acad. Sci. USA. 94:6185–6190.[Abstract/Free Full Text]

12. Bouchiat, C., M. D. Wang, J. Allemand, T. Strick, S. M. Block, and V. Croquette. 1999. Estimating the persistence length of a worm-like chain molecule from force-extension measurements. Biophys. J. 76:409–413.[Abstract/Free Full Text]

13. Bustamante, C., J. F. Marko, E. D. Siggia, and S. Smith. 1994. Entropic elasticity of lambda-phage DNA. Science. 265:1599–1600.[Free Full Text]

14. Moroz, J. D., and P. Nelson. 1997. Torsional directed walks, entropic elasticity, and DNA twist stiffness. Proc. Natl. Acad. Sci. USA. 94:14418–14422.[Abstract/Free Full Text]

15. Nelson, P. 1998. New measurements of DNA twist elasticity. Biophys. J. 74:2501–2503.[Abstract/Free Full Text]

16. Vologodskii, A. V. 1994. DNA extension under the action of an external force. Macromolecules. 27:5623–5625.[CrossRef]

17. Vologodskii, A. V., and J. F. Marko. 1997. Extension of torsionally stressed DNA by external force. Biophys. J. 73:123–132.[Abstract/Free Full Text]

18. Wang, M. D., H. Yin, R. Landick, J. Gelles, and S. M. Block. 1997. Stretching DNA with optical tweezers. Biophys. J. 72:1335–1346.[Abstract/Free Full Text]

19. Huang, J., T. Schlick, and A. Vologodskii. 2001. Dynamics of site juxtaposition in supercoiled DNA. Proc. Natl. Acad. Sci. USA. 98:968–973.[Abstract/Free Full Text]

20. Klenin, K., J. Langowski, and A. Vologodskii. 2002. Computational analysis of the chiral action of type II DNA topoisomerase. J. Mol. Biol. 320:359–367.[CrossRef][Medline]

21. Podtelezhnikov, A. A., N. R. Cozzarelli, and A. V. Vologodskii. 1999. Equilibrium distributions of topological states in circular DNA: interplay of supercoiling and knotting. Proc. Natl. Acad. Sci. USA. 96:12974–12979.[Abstract/Free Full Text]

22. Vologodskii, A. 2006. Brownian dynamics simulation of knot diffusion along a stretched DNA molecule. Biophys. J. 90:1594–1597.[Abstract/Free Full Text]

23. Cornell, W. D., P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman. 1995. A second generation force field for the simulation of proteins, nucleic acids and organic molecules. J. Am. Chem. Soc. 117:5179–5197.[CrossRef]

24. Darden, T., D. York, and L. Pedersen. 1993. Particle mesh Ewald. An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 98:10089–10092.[CrossRef]

25. Essmann, U., L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen. 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103:8577–8593.[CrossRef]

26. Cheatham III, T. E., and P. A. Kollman. 2000. Molecular dynamics simulations of nucleic acids. Annu. Rev. Phys. Chem. 51:435–471.[CrossRef][Medline]

27. Miller, J. L., T. E. Cheatham III, and P. A. Kollman. 1999. Simulation of nucleic acid structure. In Oxford Handbook of Nucleic Acid Structure. S. Neidle, editor. Oxford University Press, New York. 95–115.

28. MacKerell, A. E., Jr., N. Banavali, and N. Foloppe. 2000–2001. Development and current status of the CHARMM force field for nucleic acids. Biopolymers. 56:257–265.[CrossRef]

29. Olson, W. K., A. A. Gorin, X. J. Lu, L. M. Hock, and V. B. Zhurkin. 1998. DNA sequence-dependent deformability deduced from protein-DNA crystal complexes. Proc. Natl. Acad. Sci. USA. 95:11163–11168.[Abstract/Free Full Text]

30. Bruant, N., D. Flatters, R. Lavery, and D. Genest. 1999. From atomic to mesoscopic descriptions of the internal dynamics of DNA. Biophys. J. 77:2366–2376.[Abstract/Free Full Text]

31. Lankas, F., J. Sponer, P. Hobza, and J. Langowski. 2000. Sequence-dependent elastic properties of DNA. J. Mol. Biol. 299:695–709.[CrossRef][Medline]

32. Schellman, J. A., and S. C. Harvey. 1995. Static contributions to the persistence length of DNA and dynamic contributions to DNA curvature. Biophys. Chem. 74:2191–2198.

33. Trifonov, E. N., R. K. Z. Tan, and S. C. Harvey. 1988. Static persistence length of DNA. In DNA Bending and Curvature. W. K. Olson, M. H. Sarma, R. H. Sarma, and M. Sundaralingam, editors. Adenine Press, New York. 243–253.

34. Vologodskaia, M., and A. Vologodskii. 2002. Contribution of the intrinsic curvature to measured DNA persistence length. J. Mol. Biol. 317:205–213.[CrossRef][Medline]

35. Lavery, R., and H. Sklenar. 1988. The definition of generalized helicoidal parameters and of axis curvature for irregular nucleic acids. J. Biomol. Struct. Dyn. 6:63–91.[Medline]

36. Cheatham III, T. E., P. Cieplak, and P. A. Kollman. 1999. A modified version of the Cornell et al. force field with improved sugar pucker phases and helical repeat. J. Biomol. Struct. Dyn. 16:845–862.[Medline]

37. Mazur, A. K. 1997. Quasi-Hamiltonian equations of motion for internal coordinate molecular dynamics of polymers. J. Comput. Chem. 18:1354–1364.[CrossRef]

38. Mazur, A. K. 2001. Internal coordinate simulation method. In Computational Biochemistry and Biophysics. O. M. Becker, A. D. MacKerell, Jr., B. Roux, and M. Watanabe, editors. Marcel Dekker, New York. 115–131.

39. Mazur, A. K. 1998. Accurate DNA dynamics without accurate long range electrostatics. J. Am. Chem. Soc. 120:10928–10937.[CrossRef]

40. Mazur, A. K. 1999. Symplectic integration of closed chain rigid body dynamics with internal coordinate equations of motion. J. Chem. Phys. 111:1407–1414.[CrossRef]

41. Mazur, A. K. 2000. Theoretical studies of the possible origin of intrinsic static bends in double helical DNA. J. Am. Chem. Soc. 122:12778–12785.[CrossRef]

42. Mazur, A. K. 2001. Molecular dynamics of minimal B-DNA. J. Comput. Chem. 22:457–467.[CrossRef]

43. Arnott, S., and D. W. L. Hukins. 1972. Optimised parameters for A-DNA and B-DNA. Biochem. Biophys. Res. Commun. 47:1504–1509.[CrossRef][Medline]

44. Berendsen, H. J. C., J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684–3690.[CrossRef]

45. Mazur, A. K., and D. E. Kamashev. 2002. Simulated and experimental bending dynamics in DNA with and without regularly repeated adenine tracts. Phys. Rev. E. 66:011917-1–011917-13.

46. Andreatta, D., J. L. P. Lustres, S. A. Kovalenko, N. P. Ernsting, C. J. Murphy, R. S. Coleman, and M. A. Berg. 2005. Power-law solvation dynamics in DNA over six decades in time. J. Am. Chem. Soc. 127:7270–7271.[CrossRef][Medline]

47. Brauns, E. B., M. L. Madaras, R. S. Coleman, C. J. Murphy, and M. A. Berg. 2002. Complex local dynamics in DNA on the picosecond and nanosecond time scales. Phys. Rev. Lett. 88:158101.[CrossRef][Medline]

48. Okonogi, T. M., A. W. Reese, S. C. Alley, P. B. Hopkins, and R. H. Robinson. 1999. Flexibility of duplex DNA on the submicrosecond timescale. Biophys. J. 77:3256–3276.[Abstract/Free Full Text]

49. Trieh, M., C. Rauch, B. Wellenzohn, F. Wibowo, T. Loerting, and K. R. Liedl. 2004. Dynamics of DNA: B_I and B_II phosphate backbone transitions. J. Phys. Chem. 108:2470–2476.

50. Guenot, J., and P. A. Kollman. 1992. Molecular dynamics studies of a DNA-binding protein. 2. An evaluation of implicit and explicit solvent models for the molecular dynamics simulation of the Escherichia coli Trp repressor. Protein Sci. 1:1185–1205.[Abstract]

51. Bryant, Z., M. D. Stone, J. Gore, S. B. Smith, N. R. Cozzarelli, and C. Bustamante. 2003. Structural transitions and elasticity from torque measurements on DNA. Nature. 424:338–341.[CrossRef][Medline]

52. Bednar, J., P. Furrer, V. Katritch, A. Z. Stasiak, J. Dubochet, and A. Stasiak. 1995. Determination of DNA persistence length by cryo-electron microscopy. Separation of the static and dynamic contributions to the apparent persistence length of DNA. J. Mol. Biol. 254:579–594.[CrossRef][Medline]

53. Olson, W. K., and V. B. Zhurkin. 2000. Modeling DNA deformations. Curr. Opin. Struct. Biol. 10:286–297.[CrossRef][Medline]

This article has been cited by other articles: