J. Phys. Chem. B 2008, 112, 4975-4982
4975
Statistics of Time-Limited Ensembles of Bent DNA Conformations Alexey K. Mazur* CNRS UPR9080, Institut de Biologie Physico-Chimique, 13, rue Pierre et Marie Curie, Paris, 75005, France ReceiVed: December 17, 2007; In Final Form: February 5, 2008
The paper considers statistical properties of ensembles of chain conformations obtained by short-time Brownian dynamics (BD) of a coarse-grained DNA model in order to find out if the conditions necessary for accurate evaluation of the polymer elasticity are attainable in atom-level molecular dynamics (MD) simulations. To measure the bending persistence length (PL) with a 10% error using data accumulated in a single trajectory of a double helix of 15 base pairs, dynamics should be continued for a few microseconds. However, these estimates should be scaled down by about 2 orders of magnitude because the bending dynamics of short double helices in MD features much smaller relaxation times. As a result, good qualitative agreement with the worm-like chain (WLC) theory is reached in MD after tens of nanoseconds. The presently accessible durations of MD trajectories provide reasonably accurate evaluation of DNA elasticity and allow modeling of its mesoscopic properties. The surprisingly fast bending dynamics of short double helices in MD suggests that the microscopic mechanisms of DNA flexibility differ from a simple harmonic model.
Introduction Internal dynamics and flexibility of DNA are very important for its biological functions. Biologically relevant dynamics of DNA spans a broad range of length scales starting from a subnanometer level of single base pairs and continuing to macroscopic lengths where the double helix can be considered as a continuous flexible rod. The coupling between the macroscopic and the microscopic DNA mechanics is very interesting and it can also be important for biology because macroscopic factors such as external bending or torsion strain applied to long stretches of the double helix are involved in regulation of its activity. DNA in vivo is always found under variable degrees of torsion stress, which is arguably vital for the living cell.1 Recent in vitro studies of DNA elasticity revealed intricate relationships between torsion and stretching deformations.2,3 Notably, there is a controversy in the literature concerning the possible effect of external bending and/or stretching strain upon the torsion stiffness of the double helix.4,5 The corresponding molecular mechanisms are puzzling and accurate atom-level molecular modeling can help in disclosing them. In order to consider the above problems in the context of atom-level modeling, one has to establish a correspondence between the atomistic and the coarse-grained DNA representations. 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) model6,7 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 that are measured since sixties with progressively improved accuracy.8-11 On the other hand, local atom level DNA properties are reasonably well-reproduced in classical molecular dynamics (MD) simulations of small duplexes in explicit aqueous environment.12,13 MD simulations are based upon empirical force fields * Corresponding author.
that are parametrized by using experimental data for small molecules as well as quantum mechanics calculations.13,14 The experimental DNA elasticity is not taken into account because there are no practical ways to do that. Evaluation of elastic properties of atomistic DNA models encounters many problems,15 and one of them is that we do not know how long it would take for real DNA to visit the number of distinct conformational states necessary for evaluation of its elastic properties. Accordingly, it is not known if the necessary sampling is theoretically attainable in MD. A few authors earlier tried to evaluate the elastic parameters of DNA from atom-level data. The first attempt was made by Olson et al. on the basis of fluctuations of X-ray DNA structures rather than simulations.16 Other authors later tried to extract this information from classical MD trajectories.15,17,18 All of these studies arrived at a reasonable agreement with experiment although the accuracy of the obtained estimates was not verifiable. The most recent study15 used qualitative agreement with the WLC theory and time convergence of elastic parameters measured by different techniques as an internal criterion of accuracy. It was found that 20 ns of atomistic MD of 25-mer duplexes is sufficient to reach a demonstrative agreement with the WLC theory and to get convergent estimates of elastic parameters. However, it is possible that the above criterion is fulfilled before the true statistical convergence is reached and that, with longer MD, the measured WLC parameters could drift very far, albeit slowly. This issue can be clarified by analyzing trajectories of Brownian dynamics (BD) of coarse-grained DNA models that reproduce kinetics of DNA fluctuations in solution. A suitable beads-on-string model was originally introduced by Allison et al.19,20 and later improved by different authors.21,22 The recent version of this model predicts the rate of diffusion of knots along a stretched DNA in good agreement with experiment.23,24 In the present paper, we consider an inverse problem of WLC dynamics, that is, extracting elastic parameters from time-limited ensembles of configurations. DNA fragments of a few helical turns are modeled as discrete WLCs with hydrodynamic
10.1021/jp711815x CCC: $40.75 © 2008 American Chemical Society Published on Web 03/29/2008
4976 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Mazur
parameters earlier shown to reproduce experimental data. Timelimited ensembles of chain configurations are obtained by BD simulations and compared with the results of atomistic MD. It appears that at least a few microseconds of BD are necessary in order to accumulate an ensemble of distinct bent conformations sufficient for a reasonably accurate evaluation of the bending persistence length (PL). A similarly long duration is necessary to arrive at a qualitative agreement with the WLC theory. In contrast, in atomistic MD, a comparable degree of qualitative agreement with the theory is already observed after tens of nanoseconds. This difference is due to much faster bending dynamics predicted by the atomistic MD, which can be confirmed by a direct comparison of the corresponding time autocorrelation functions. Interestingly, the small bending relaxation times do not correspond with the measured PLs, which indicates that the microscopic mechanisms of bending of short double helices significantly differ from simple harmonic flexibility. The absolute rates of bending dynamics in MD may not agree with experiment; nevertheless, the duration of MD trajectories accessible at present appears sufficient for a reasonably accurate evaluation of elastic parameters and the studying of mesoscopic properties of double helical DNA.
where Yf and B stand for the Young modulus and the stretching PL, respectively. Parameter B was introduced in order to use the same units for different types of DNA elasticity.26 We carried out all BD simulations under T ) 300 K, with parameters corresponding to A ) 50 nm and B ) 100 nm. BD Simulations. BD simulations were carried out by using the algorithm earlier developed and adapted for DNA by several authors.21,22,27,28 In the original (first-order) procedure by Ermak and McCammon,27 the one-step particle propagation is computed as
Theory and Methods
Subscript zero in eq 2 refers to the initial state. Similar to earlier DNA simulations,22 tensor Dij was computed from particle coordinates by using Rotne-Prager formulas for non-overlapping beads.29 The corresponding random force components were obtained from 3N independent standard Gaussian random generators as proposed by Ermak and McCammon.27 To increase the time step, we used the second-order version of this algorithm.21,22,28 The time step was 25 ps, which is smaller than that in earlier studies because of the increased stretching stiffness. Atomistic MD Simulations. The BD of the coarse-grained DNA model is compared with classical MD of atom-level models of 15-mer and 25-mer DNA duplexes with the ATalternating sequence. This choice is appropriate because the ATalternating duplexes are homopolymers of AT-units; therefore, they cannot have distinguished asymmetric structures like static bends. We mainly consider the data earlier used for evaluation of DNA elastic parameters15 supplemented by a 30 ns trajectory for the 15-mer. The simulations were carried out in explicit aqueous environment with a neutralizing amount of Na+ ions. The AMBER98 forcefield parameters30,31 were used with the rigid TIP3P water model32 in NVT ensemble conditions under T ≈ 293 K, with the water density around 0.997. To increase the time step, MD simulations were carried out by the internal coordinate method (ICMD),33-36 with the internal DNA mobility limited to essential degrees of freedom. The double helical DNA was modeled with all backbone torsions, free bond angles in sugar rings, and rigid bases and phosphate groups. The effect of these fixation is insignificant, which was checked by comparison with standard Cartesian dynamics of the same 25mer DNA fragment. The time step was 0.01 ps and the DNA conformation were saved every 2.5 ps. Bending in the saved DNA conformations was analyzed by the Curves algorithm,37 with one of the fitting parameters increased by a factor of 10. This modification was described in the previous paper15 where one can also find additional details of the simulation protocols. The normalized time autocorrelation function C(t) for bending in different fragments is computed from the time series of the corresponding bend angle θ(t) as
Beads-on-String DNA Model. Following earlier studies,20-22 a DNA molecule is modeled by N beads of radius a linked by N - 1 bonds. The total energy is composed of bond stretching and bending terms as
U)
h N-1 2
∑ i)1
(li - l0)2 +
g N-2 2
θ2i ∑ i)1
Here, li and l0 are the length of the ith link and its equilibrium value, respectively; θi is the angular displacement of bond i + 1 relative to bond i. All other interactions between the beads can be neglected because we consider only short stiff chains that cannot loop. We are mainly interested in the bending elasticity; the stretching dynamics is analyzed for comparison. For simplicity, torsion degrees of freedom are not considered. The choice of parameters a, h, l0, and g was considered in detail in the earlier literature.22,25 We use the equilibrium link length l0 ) 5 nm and the bead radius a ) 1.78 nm, which is close to the values employed in the recent BD simulations of DNA.23 In order for different chain length dependences to be smoothly extrapolated to zero, the minimal working length should also equal l0. Therefore, the equivalent DNA length is computed as l0(N - 2) and l0(N - 1) for bending and stretching, respectively. Parameters g and h are chosen so that the theoretical values of experimental observables correspond to the data available for DNA. In the limit of small angle deviations, the bending PL can be computed by using the general theory of fluctuations6 as
A)
gl0 kT
(1)
where k is the Boltzmann constant and T is the temperature. Equation 1 is strictly valid only for l0 , A, but the difference from the exact procedure22 is negligible in our case. Stretching observables are computed as
Yf ) bl0,
B)
bl0 3.4 nm 2 kT 2π
(
)
ri )
r0i
+
∑j
D0ij F0j kT
∆t + Ri(∆t)
(2)
where ri is one of the 3N Cartesian coordinates, Fj is the corresponding component of the regular force, Dij is the hydrodynamic diffusion tensor, and Ri(∆t) is the Gaussian random force sampled from a multivariate distribution with
〈Ri(∆t)Rj(∆t)〉 ) 2D0ij∆t
C(t) )
C ˜ (t) ; C ˜ (t) ) 〈θ(τ + t)θ(τ)〉τ - 〈θ〉2 C ˜ (0)
Ensembles of Bent DNA Conformations
J. Phys. Chem. B, Vol. 112, No. 16, 2008 4977
Evaluation of Elastic Parameters. The WLC theory offers a few alternative ways for evaluating the bending PL from a canonical ensemble of chain configurations.15 The average bending angle grows with the chain length L according to the following relationship
Da(L) ) -ln(〈cos θ〉) )
L Aa
(3)
This linear law, when checked in computations, provides a useful test for convergence and the correspondence with the WLC theory. The PL can be obtained either from a single measurement or from the slope of a linear regression corresponding to eq 3. Subscript a in eq 3 indicates that the corresponding values are computed from angle deviations and it is dropped below. The bending PL can also be obtained from fluctuations of end-to-end distances, but this approach is not considered here because the corresponding linear approximation15 additionally requires that l0 is small compared with L. Finally, the probability distribution of the bending angle can be written as
[
dP ∼ exp -
A (1 - cos θ) d(1 - cos θ) L
]
(4)
which gives descending straight lines in coordinates (ln P) versus (1 - cos θ).38 The predicted regular variation of the slopes with the chain length is used to check the qualitative agreement with the WLC theory. The value of the bending PL can also be obtained from fits of eq 4 to MD data. The stretching PL is evaluated from a relationship similar to eq 3 that can be written as15
Db(L) ) ∆2L )
L 3.4 nm 2 B 2π
(
)
(5)
where ∆L is the standard deviation of the chain length from its equilibrium value L. Results Effects of Trajectory Lengths and Sampling Intervals. Because of common computational limitations, the elastic parameters of DNA models have to be estimated from limited ensembles of configurations usually taken from a single MD trajectory. Such calculation can be viewed as a stochastic process and characterized by a probability of hitting a certain error interval around the thermodynamic limit. Accordingly, we will be interested in the probability distributions for the bending and stretching PLs computed by different methods from a single DNA trajectory of a given duration. Suppose we have a large set of trajectories, each providing one count of the PL. If the DNA configurations are saved with a large time interval, the PL probability distribution will be very broad because every count comes out from a few structures. With the saving interval reduced, the sampling volume increases and the accuracy will improve until the interval becomes so small that the consecutive DNA conformations are essentially identical. For our purposes, the sampling intervals should be small enough so that the accuracy of the measured PLs depends only on the duration of trajectories. The results of the foregoing analysis are shown in Figure 1 for a chain of 4 beads (total length 15 nm). The bending and stretching PLs were computed by averaging fluctuations of one angle and one chain link by using eqs 3 and 5, respectively. The average PL values are very close to the analytical
Figure 1. Selection of optimal sampling time steps for evaluating bending and stretching PLs. Data from a series of BD trajectories of identical duration were analyzed, gradually reducing the time interval between the saved configurations. Each trajectory gives one count of the measured PL. The probability distributions shown in the figure represent results of 100 independent counts. The vertical gray bars show the scattering of the corresponding average values. Upper panel: The bending PL; duration of trajectories, 12.5 µs. Lower panel: The stretching PL; duration of trajectories, 1.25 µs. The sampling intervals are shown in both panels in nanoseconds.
predictions suggesting that the method of Brownian dynamics we use correctly samples from the canonical ensemble of configurations. Both the bending and stretching PLs appear slightly overestimated, which can be attributed to the discretization errors of the algorithm. It is seen that, for the bending PL, the probability distribution stops evolving when the sampling interval becomes less than 5 ns. The corresponding value for stretching is 0.125 ns. The difference is certainly due to the slower bending dynamics. The time autocorrelation function of the bend angle in a chain of 3 beads exhibits simple exponential relaxation with the relaxation time τb ≈ 3 ns. The chain length relaxation in a chain of two beads is also exponential with the relaxation time τs ≈ 0.15 ns. In longer chains, the relaxation law deviates from exponential,22 but, as we see, the short-chain relaxation time gives a good estimate of the optimal sampling interval. For the following analysis, we fix the sampling interval at 2.5 and 0.125 ns for bending and stretching, respectively. The effect of the trajectory length upon the accuracy of the measured PLs is shown in Figure 2. The BD data were processed similarly to Figure 1; that is, dynamics of a 15 nm chain was considered, with all internal angles and links used for averaging in eqs 3 and 5. The probability distributions are bell-shaped and they narrow progressively as the trajectory length is increased. For durations below a microsecond, the probability distributions are very broad especially for the bending PL. The average PL values are reasonably close to the analytical predictions, with the difference decreasing as the trajectory length is increased. By integrating the distributions shown in Figure 2, one can evaluate the probability of getting a good PL estimate from a single trajectory (Figure 3). It is seen that in order to obtain with 90% probability the bending PL within 10% error one should have a trajectory of at least 1.25 µs. The corresponding duration for stretching is 62.5 ns. In both cases,
4978 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Mazur
Figure 4. Probability distributions of the bending PLs computed from BD trajectories of 3.75 µs by four different methods: by using eq 3 for one angle (filled squares); from the slope of the Da(L) plot at L ) l0 (open circles); from the average slope of the Da(L) plot (filled circles); from the slope of the linear regression straight line approximating all Da(L) data (open squares). The vertical gray bar shows the scattering of the corresponding average values.
Figure 2. Probability distributions of the bending and stretching PLs computed from BD trajectories of different duration. The lengths of trajectories are shown in the figure in microseconds (upper panel) and nanoseconds (lower panel). The average PLs are shown by dashed lines for the shortest trajectories and by gray bars for the rest.
Figure 3. Probability that a PL value computed from a single trajectory is found within a given error interval from the average. The lengths of trajectories are shown in microseconds for bending (upper panel) and in nanoseconds for stretching (lower panel).
these values are about 500 times larger than the short chain relaxation times. On the other hand, for a trajectory of 125 ns, which is close to the limit of modern MD simulations, only 40% accuracy of the bending PL is warranted with 90% probability while the chances to hit a 10% error interval are only 30%. As we see, the time convergence is very slow especially for the bending PL. In principle, the rate of convergence may vary between different methods of computing PLs,15 and it is worth checking different possibilities. One of them consists of using not only the first point of the Da(L) plots corresponding to L ) l0 (as for Figures 1-3) but also the slopes of such plots. This is necessary in practice for extrapolating to the long chain limit because bends in short DNA cannot be measured accurately. A
few alternative possibilities of measuring the bending PL are compared in Figure 4. Indeed, the widths of the probability distributions vary considerably indicating different rates of convergence. The method employed for Figures 1-3 is more accurate because it uses only single angle fluctuations that have the largest statistical weight as well as the highest relaxation rate. In contrast, a linear regression analysis of data for all chain lengths results in much slower convergence due to a contribution of slower bending fluctuations of long chains. Sampling of Internal Fragments of Long Chains. The sampling of configurations of short DNA is significantly enhanced by considering all internal fragments of longer molecules. For instance, addition of one base pair to a nanomer DNA fragment only slightly augments the computational load whereas the number of available nanomer configurations is doubled. This trick was used above as well as in the recent MD studies of DNA elasticity.15,18 By extrapolating to infinity, one concludes that the elastic parameters can be computed from a single (random) configuration of a very long chain. The last statement is perhaps correct for stretching, but not for bending dynamics. The effects of chain elongation are checked in Figure 5. Plots shown by open symbols demonstrate the narrowing of the PL probability distributions obtained by three times increasing the trajectory lengths. To this end, we use the same data as in Figure 3. In all cases, the PL values were computed from fluctuations of one angle/link by using eqs 3/5. Notation “2/ 1.25” means that the corresponding bending data were obtained from simulations of a chain with two angles (four beads) during 1.25 µs. Accordingly, code “3/62.5” refers to stretching data from 62.5 ns trajectories of the same chain (three links). The filled symbols mark the probability distributions obtained for three times longer chains. Distributions 6/1.25 and 9/62.5 were obtained by considering all internal fragments; that is, the number of chain configurations is identical to that in the narrow distributions 2/3.25 and 3/187.5, respectively. For stretching, the effect of chain elongation is approximately equivalent to that of increased trajectory length. In contrast, for bending PL, the probability distributions seem to depend only on the trajectory duration; that is, the enhancement of sampling due to chain elongation does not result in a noticeable narrowing of the distribution. The expected effect of chain elongation could be masked by that featured in Figure 1, that is, the correspondence between the rate of relaxation and the sampling time step. If in longer chains the bending dynamics of internal fragments is significantly slower, the earlier selected sampling time interval may be too small. In this case, the consecutive configurations are
Ensembles of Bent DNA Conformations
J. Phys. Chem. B, Vol. 112, No. 16, 2008 4979
Figure 6. Chain length dependences of the angular deviation defined by eq 3 obtained for a 10-bead WLC chain in three BD trajectories with increased duration. The durations of trajectories are shown in the upper left corner (in microseconds).
Figure 5. Probability distributions of the bending (upper panel) and stretching (lower panel) PLs computed from BD trajectories of WLC chains by varying both the chain lengths and the duration of trajectories. The plots are marked by two-number codes of the form P1/P2. Parameter P1 shows the number of angles (upper panel) and links (lower panel) in the modeled chain. Parameter P2 shows the length of trajectories in microseconds (upper panel) and nanoseconds (lower panel). The prime indicates that in this case only the first angle/link of the chain were considered. Otherwise all internal fragments of the chain were taken into account. The vertical gray bar shows the scattering of the corresponding average PL values.
similar and the effective number of configurations provided by three times longer chains would be smaller than expected. This possibility was checked by computing the PL probability distributions from fluctuations of a single angle/link in the longer chains (codes 6/1.25′ and 9/62.5′ in Figure 5). These distributions result from the smallest number of configurations; therefore, they should be wider than others in Figure 5, and the suggested slow bending relaxation should make the corresponding distribution in the upper panel still wider. Unlike that, the widening is seen only for stretching, whereas for bending PL, there is essentially no effect at all. It means that fluctuations of internal chain angles occur at similar rates regardless of the chain length. To explain the results shown in the upper panel of Figure 5, we have to admit that fluctuations of internal angles in long chains sampled during finite time intervals are not statistically independent. Analysis of cross-correlation functions between different angles and links of the same chain does not show any correlations. However, not all types of coupling can be revealed by such analysis. Suppose we fix the distance between the opposite ends of a relatively short chain at an unusually small value. Its internal fragments will exhibit dynamics corresponding to a very low PL similar for all of the fragments. The internal angles of such a chain are not independent, but they are not correlated either. This example possibly reproduces the essential qualitative feature of DNA dynamics modeled in BD. To sample all configurations corresponding to a given PL, the ends of a long chain must cover very long distances compared with those of its internal fragments. With strong friction, it essentially means that the internal fragments can equilibrate like in the above example before the chain ends move significantly. The foregoing considerations probably do not apply to stretching because according to eq 5 the distances corresponding to chain length fluctuations grow as xL, which is much slower than for
bending. The evident practical conclusion is that, in MD studies of the bending elasticity, the trajectory length cannot be traded for the chain length; consequently, the DNA length should be reduced as much as possible. WLC Consistency Checks. The foregoing observations and conclusions question the credibility of the recently measured elastic parameters of atomistic DNA models.15 Extrapolated to common conditions of MD, these observations suggest that trajectories of 30 ns and less can give essentially arbitrary values of the bending PL depending upon the initial conditions. On the other hand, it was shown earlier that in many cases MD data exhibited nontrivial qualitative agreement with the WLC theory.15,38 A question arises as to if the qualitative agreement with the WLC theory is a valid criterion of good sampling and statistical convergence. To check how this approach works in BD, we use eqs 3 and 4. A BD trajectory of a 10-bead chain was chosen from the set analyzed in Figure 5 so that the bending PL evaluated from fluctuations of a single angle significantly exceeded the true value. This trajectory was continued to 12.5 µs and the chain configurations from the first 0.125, 1.25, and 12.5 µs were analyzed separately. The Da(L) dependences corresponding to eq 3 are shown in Figure 6. The agreement with the linear form predicted by the WLC theory is rather poor for trajectories of 0.125 and 1.25 µs, and only with 12.5 µs does it become good. This is true for long as well as for short chains. Comparing it with Figure 3, we see that the qualitative agreement with the WLC theory becomes evident when a single trajectory has an almost 100% probability to hit the 10% error interval around the true bending PL. The last conclusion is further corroborated by Figure 7 where the shapes of the bend angle distributions are analyzed. A regular pattern with all computed distributions close to the theoretical straight lines emerges by 1.25 µs and further improves at the last panel. For comparison, a similarly nice regular pattern was obtained in MD by 20 ns.38 Note that in the first two panels the computed distributions sometimes display partial agreement with the theoretical straight lines. This is the case even when the PLs computed from the average bending angles are incorrect (the latter can be visually estimated from the plots in Figure 6). For instance, after 125 ns, the average bending angle of 5 nm fragments gives an apparent bending PL beyond 100 nm. The corresponding plot in Figure 7 (the bottom trace in the left panel) reveals that the angle distribution has a reasonably good WLC shape but with a much lower PL. The discrepancy is due to a steep initial slope of the plot in Figure 7 as well as the low populations of large angles. At the same time, no temporary WLC-like patterns corresponding to incorrect PLs can be seen at early stages of simulations.
4980 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Figure 7. Probability distributions of bend angles in fragments of L ) 5, 10, 15, 20, 25, 30, 35, and 50 nm observed in three BD trajectories of a 10-bead WLC chain. Different curves in the three panels refer to increasing chain lengths from bottom to top. Data counts were accumulated in 80 windows evenly spaced on -1 < cos θ < 1. The coordinates are chosen to linearize the small angle WLC distributions. The straight lines correspond to the bending PL of 51 nm. All plots are consecutively shifted by 3.0 for clarity. The durations of trajectories are shown in the upper right corners of the three panels (in microseconds).
Mazur
Figure 9. Initial rates of relaxation of bending and stretching fluctuations in DNA according to the discrete WLC model. The relaxation times were estimated from the initial slopes of time autocorrelation functions partially shown in Figure 8.
Figure 10. Relaxation of bending fluctuations in DNA according to classical MD. Dynamics of double helices with AT-alternating base pair sequences were modeled in explicit water environment with neutralizing amount of sodium ions. The solid traces show the time autocorrelation functions for internal fragments of a 25-mer DNA. The numbers near the curves show the corresponding chain lengths (in base pairs). The length of 15 bp fragments is close to 5 nm. The dashed line shows the autocorrelation function for bending of an isolated 15mer DNA with the same sequence modeled in identical conditions. Figure 8. Relaxation of bending fluctuations in DNA according to the discrete WLC model. Time autocorrelation functions are shown for bend angles in internal fragments of a 10-bead WLC chain. The curves are labeled by the corresponding DNA lengths (in nanometers). The dashed line shows the autocorrelation function for bending of an isolated chain of 5 nm.
Rates of Coarse Grained and Atomistic DNA Dynamics. Thus, in terms of the degree of the qualitative agreement with the WLC theory a few microseconds of BD of a coarse-grained DNA model correspond to a few tens of nanoseconds of all atom MD. To explain this contradiction, we should assume that the bending dynamics of equivalent DNA fragments in MD is 2 orders of magnitude faster. Figure 8 exhibits the time autocorrelation functions for bending of internal fragments in a discrete WLC model of DNA. For fragments longer than half of the full length, the plots nearly overlap and diverge only at longer times. The relaxation is multiexponential for all internal fragments because their bending dynamics is affected by the relaxation modes of both shorter and longer chains. For instance, the initial slope of the autocorrelation function of the full-length 40 nm DNA is larger than the long-time limit due to the faster motions of the chain ends. At the opposite end of the spectrum, the initial slope of the autocorrelation function of 5 nm chains gives exactly the relaxation time for bending of isolated 5 nm DNA. In the latter case, the correlation decay is exponential (the dashed line in Figure 8). For internal stretches of long DNA, however, this regime accounts for only a part of the whole process dominated by slower modes imposed by longer chains by the mechanism discussed above in relation to Figure 5. The
foregoing considerations apply to all other chain lengths; therefore, the overall pattern is strongly hierarchical. In the long time limit, all plots in Figure 8 become parallel to the top trace. The relaxation times corresponding to the initial rates of the correlation decay are shown in Figure 9 for bending and also for stretching fluctuations. The relaxation times naturally grow with the chain length and in both cases the growth is approximately linear, with some signs of saturation. (Note that these plots do not have to extrapolate to zero.) The bending relaxation times grow much faster, which partially explains the qualitative difference between bending and stretching revealed in Figure 5. In both cases, however, the growth is slower than the analytical laws corresponding to diffusion of rigid cylinders in viscous media, that is, L2/ln(L) and L4 for stretching and bending, respectively.18 A few representative time autocorrelation functions for bend angles in DNA modeled by atom level MD are displayed in Figure 10. It is readily seen that the relaxation is indeed markedly faster for the bending dynamics shown in Figure 10 compared with those in Figure 8. Nevertheless, the overall patterns are qualitatively similar. All plots in Figure 10 exhibit a rapid fall within the 2.5 ps interval between the saved configurations. This is due to small scale atom motions that cannot be unequivocally interpreted in terms of the overall bending of DNA. The apparent amplitude of this fast decay varies when the bend angles are evaluated with different internal parameters of the Curves algorithm as discussed elsewhere.15 This difficulty precludes comparison of initial slopes as in Figure 8; but for the chain lengths of 15 pb and longer, the amplitude
Ensembles of Bent DNA Conformations of the initial fall is negligible, and the relaxation is close to exponential. The limiting relaxation time is that of a 25-mer DNA and it is about 0.12 ns. The relaxation time corresponding to the coarse-grained model considered in Figures 1-3 is that of a 10 nm chain in Figure 8, that is, about 7 ns. Therefore, the bending in the atom level DNA proceeds 60 times faster, which means that, in terms of convergence, trajectories of 30 ns in MD and 1.7 µs in BD should be comparable. Referring to Figure 3, we conclude that the bending PL evaluated for the atomistic DNA models from 20 to 30 ns MD trajectories of 25-mers is probably accurate to 15% or better, which also explains a good qualitative agreement of MD data with the WLC theory.15,38 Figure 10 also shows the time autocorrelation function for bending of an isolated 15-mer DNA with the same sequence. In this case, the relaxation is much faster, with the difference from BD approaching 2 orders of magnitude. The correlation decay significantly deviates from the exponential form, but the effective relaxation rate is approximately 3 times higher than that for the free 25-mer duplex. The difference in lengths is by a factor of 1.7 only; therefore, the relaxation time grows at least as L2, which is faster than that in Figure 8, but still much slower than the L4 law predicted by the model of rigid cylinders.18 The most striking, however, are the very small absolute values of the bending relaxation times revealed by Figure 10. Both the discrete WLC model in Figure 8 and the analytical estimates18 predict relaxation times of a few nanoseconds for DNA fragments of 5-10 nm. The fast relaxation like in Figure 10 can hardly be reproduced by any coarse grained model of DNA with realistic radii of beads or cylinders, and the elasticity corresponding to the known value of the bending PL. The evident discrepancy between Figure 8 and Figure 10 requires further examination and additional studies. Discussion A number of intriguing aspects in dynamics of double helical DNA is related with the coupling between its macroscopic and microscopic states. It is well-known that the external stress applied to long chains causes conformational transitions, that it changes the local rates of base pair opening, that the macroscopic DNA elasticity depends upon the base pair sequence, and so forth. The corresponding molecular mechanisms are only partially understood. A link between the properties of long DNA chains and their atom level dynamics can be established by comparing the WLC model of polymer DNA with the conformational ensembles produced by MD simulations. As a first step toward this objective, a few recent studies probed the bending, stretching, and torsion elasticity of DNA in classical MD.15,18 A surprisingly good agreement with the WLC theory was observed, with the measured PLs reasonably close to experiment. At the same time, the durations of MD trajectories looked small taking into account the experimental estimates of the corresponding relaxation times. They were also much smaller than the common lengths of trajectories in BD simulations of coarse grained models with parameters directly fitted to hydrodynamic properties of long DNA chains. This was not a contradiction since the DNA length in MD is significantly smaller than in both experiments and BD simulations, and since it is not known how long it would take for real DNA to visit the number of distinct conformational states necessary for evaluation of its elastic properties by the specific methods employed in MD. Nevertheless, the accuracy of these MD results looked questionable. Here, we tried to shed light upon the foregoing controversy by processing BD trajectories of the discrete WLC model of
J. Phys. Chem. B, Vol. 112, No. 16, 2008 4981 DNA in the same way as MD data, using limited durations of trajectories. We take advantage of the fact that all of the PLs in this case can be directly computed from input parameters by well-known analytical formulas, which makes possible proper checks of the accuracy and time convergence. It is shown that reasonably accurate values of PLs can be obtained with high probability from a single trajectory only if its duration exceeds the corresponding relaxation time by at least 2 orders of magnitude. Moreover, the accuracy of the bending PL is always limited by the relaxation time of the whole chain, and it cannot be significantly improved by considering only the dynamics of shorter subfragments. MD simulations, however, suffer from these limitations less than expected because DNA bending in MD turns out to be much faster than in the discrete WLC model. The difference is about 2 orders of magnitude, which approximately corresponds to the difference in durations of trajectories in MD and BD. Because of the small bending relaxation time, the atom-level MD reaches good agreement with the WLC theory after a few tens of nanoseconds. Thus, the above apparent controversy is completely clarified. The surprisingly fast bending of DNA in the atom-level MD simulations looks intriguing from two different points of view. First, it is interesting how this observation compares with experiment. Second, it is not clear how this result can be explained in terms of the intrinsic mechanisms of the employed models. To my knowledge, there is only one published experimental estimate of the bending relaxation time in short DNA of comparable length.39 The reported value, about 0.4 ns, reasonably agrees with MD taking into account that the viscosity coefficient of TIP3P water is three times lower than the experimental value.40 Unfortunately, this estimate was obtained indirectly by using a model-dependent data treatment. Experimental data for long DNA usually are well-described by the discrete WLC model;41,42 however, the corresponding low relaxation rates39,43 cannot be extrapolated to small chain lengths because of the possible rapid nonlinear growth of the bending relaxation time with the chain length.18,39 Comparison with experiment, therefore, should be left aside for the moment, and we consider the second (simpler) issue, namely, if the difference in rates of the bending relaxation in MD and BD models of DNA can be rationalized in terms of internal parameters and mechanisms. In the discrete WLC model, the rate of bending relaxation is determined by the bending energy coefficient g given by eq 1, the solvent viscosity, and the bead diameter. The first two parameters are just experimental values. The third parameter was originally fitted to sedimentation data,44 and it is somewhat larger than the true DNA diameter, which is attributed to water and condensed ions that accompany the double helix in longrange translations. Bending in small DNA may occur without such translations; therefore, the effective diameter probably should be reduced to the true value, but this can give acceleration by a factor of 2, not more. Thus, a flexible body with dimensions of a 25-mer double helix and harmonic elasticity corresponding to the DNA bending PL will show in water a relaxation time in the nanosecond time range. Despite its complex structure, the atomistic model of the same fragment cannot have a hundred times smaller effective diameter, and one may wonder how this object arrives at that high bending rate. Evidently, this can be due to either very low solvent viscosity or very strong bending stiffness. The viscosity of the TIP3P water model is smaller than the experimental value by a factor of 3.40 Accordingly, the diffusion rates of small chain molecules are overestimated, but very good agreement with experiment is achieved by simple
4982 J. Phys. Chem. B, Vol. 112, No. 16, 2008 rescaling.45 Together with about two times overestimated bending PL,15 the low viscosity of the TIP3P water can account for acceleration of the bending dynamics by a factor of 3-5 maximum. Consequently, the fast relaxation in MD can only be due to a much higher effective bending stiffness of DNA. It means that the “kinetic” stiffness that governs bending relaxation differs from the “thermodynamic” value computed from the measured PL by eq 1, that is, in the harmonic approximation. To account for our data, the difference can reach a factor of 10. The possibility of divergence between the kinetic and the thermodynamic DNA stiffness was invoked in the earlier literature. For instance, the concepts of static and dynamic bending PLs were used for interpretation of electron paramagnetic resonance (EPR),39 as well as time-resolved fluorescence polarization anisotropy (FPA,)46 and transient polarization grating (TPG)43 studies of 14-200 bp DNA. According to this model,47 the total PL has a contribution from permanent static bends in DNA that increase the average bend angles but do not affect the sub-microsecond fluctuations probed by the spectroscopic methods. The estimated dynamic part of the stiffness was 3-4 times larger than the harmonic estimate corresponding to the total PL.43,46 This particular model is not applicable to our case because the value of the bending PL and the relaxation time are estimated from the same MD trajectory; that is, they correspond to the same time scale. Also, the choice of the ATalternating sequence ensures that the DNA fragments considered cannot have static bends. However, it is not the only model to be considered. The compressed backbone hypothesis (CBH) predicts metastable local bends in a homopolymer double helix due to regular modulations of the minor and major DNA grooves.48-50 In this case, bending fluctuations occur by migration of the phase of the groove modulations, which means that it is driven by concerted conformational transitions in the DNA backbone. After every transition, the DNA molecule appears slightly strained, and this internal strain drives its trace to a new configuration. In such a system, the rate of relaxation may not correspond to the value of the bending PL. The rate of DNA bending in MD and the mechanisms of relaxation should be checked in further studies. It will be interesting to see how the solvent viscosity affects the bending relaxation rates. To this end, one can employ water models with the viscosity coefficient closer to the experiment. It would also be interesting to compare the torsion elasticity of DNA observed in MD with those from coarse grained BD simulations. This issue was omitted here because including the torsion degree of freedom into the treatment of atomistic DNA models calls for analysis of bend directions and anisotropy of bending, which significantly complicates the problem and requires for comparison extended discrete WLC models that were not tested earlier. Qualitative and quantitative comparison of torsion and bending relaxation in atom-level DNA simulations can also be instructive in the context of the present discussion. This work is now in progress. Acknowledgment. The author is very grateful to Alexander Vologodskii for instructive advice, comments, and useful discussions. References and Notes (1) Travers, A.; Muskhelishvili, G. Nat. ReV. Microbiol. 2005, 3, 15769. (2) Gore, J.; Bryant, Z.; Nollmann, M.; Le, M. U.; Cozzarelli, N. R.; Bustamante, C. Nature 2006, 442, 836-9. (3) Lionnet, T.; Joubaud, S.; Lavery, R.; Bensimon, D.; Croquette, V. Phys. ReV. Lett. 2006, 96, 178102.
Mazur (4) Bryant, Z.; Stone, M. D.; Gore, J.; Smith, S. B.; Cozzarelli, N. R.; Bustamante, C. Nature 2003, 424, 338-41. (5) Fujimoto, B. S.; Brewood, G. P.; Schurr, J. M. Biophys. J. 2006, 91, 4166-79. (6) Landau, L. D.; Lifshitz, E. M. Statistical Physics, Part 1; Nauka: Moscow, 1976. (7) Cantor, C. R.; Schimmel, P. R. Biophys. Chem. Part III: The BehaVior of Biological Macromolecules; W. H. Freeman: San Francisco, 1980. (8) Hagerman, P. J. Ann. ReV. Biophys. Biophys. Chem. 1988, 17, 265286. (9) Bustamante, C.; Marko, J. F.; Siggia, E. D.; Smith, S. Science 1994, 265, 1599-600. (10) Vologodskii, A. V. Macromolecules 1994, 27, 5623-5625. (11) Bouchiat, C.; Wang, M. D.; Allemand, J.; Strick, T.; Block, S. M.; Croquette, V. Biophys. J. 1999, 76, 409-13. (12) Miller, J. L.; Cheatham, T. E., III; Kollman, P. A. In Oxford Handbook of Nucleic Acid Structure; Neidle, S., Ed.; Oxford University Press: New York, 1999; pp 95-115. (13) Cheatham, T. E., III; Kollman, P. A. Annu. ReV. Phys. Chem. 2000, 51, 435-471. (14) MacKerell, A. E., Jr.; Banavali, N.; Foloppe, N. Biopolymers 20002001, 56, 257-265. (15) Mazur, A. K. Biophys. J. 2006, 91, 4507-4518. (16) Olson, W. K.; Gorin, A. A.; Lu, X. J.; Hock, L. M.; Zhurkin, V. B. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 11163-8. (17) Bruant, N.; Flatters, D.; Lavery, R.; Genest, D. Biophys. J. 1999, 77, 2366-76. (18) Lankas, F.; Sponer, J.; Hobza, P.; Langowski, J. J. Mol. Biol. 2000, 299, 695-709. (19) Allison, S. A.; McCammon, J. A. Biopolymers 1984, 23, 363375. (20) Allison, S. A. Macromolecules 1986, 19, 118-124. (21) Chirico, G.; Langowski, J. Macromolecules 1992, 25, 769-775. (22) Jian, H.; Vologodskii, A. V.; Schlick, T. J. Comput. Phys. 1997, 136, 168-179. (23) Vologodskii, A. Biophys. J. 2006, 90, 1594-1597. (24) Vologodskii, A. In Computational Studies of DNA and RNA; Lankas, F., Sponer, J., Ed.; Springer: Netherlands, 2006; pp 579-604. (25) Jian, H.; Schlick, T.; Vologodskii, A. J. Mol. Biol. 1998, 284, 28796. (26) Nelson, P. Biophys. J. 1998, 74, 2501-2503. (27) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 13521360. (28) Iniesta, A.; de la Torre, J. G. J. Chem. Phys. 1990, 92, 20152018. (29) Rotne, J.; Prager, S. J. Chem. Phys. 1969, 50, 4831-4837. (30) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (31) Cheatham, T. E., III; Cieplak, P.; Kollman, P. A. J. Biomol. Struct. Dyn. 1999, 16, 845-862. (32) Jorgensen, W. L.; Chandreskhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (33) Mazur, A. K. J. Comput. Chem. 1997, 18, 1354-1364. (34) Mazur, A. K. J. Am. Chem. Soc. 1998, 120, 10928-10937. (35) Mazur, A. K. J. Chem. Phys. 1999, 111, 1407-1414. (36) Mazur, A. K. In Computational Biochemistry and Biophysics; Becker, O. M., MacKerell, A. D., Jr., Roux, B., Watanabe, M., Eds.; Marcel Dekker: New York, 2001; pp 115-131. (37) Lavery, R.; Sklenar, H. J. Biomol. Struct. Dyn. 1988, 6, 63-91. (38) Mazur, A. K. Phys. ReV. Lett. 2007, 98, 218102. (39) Okonogi, T. M.; Reese, A. W.; Alley, S. C.; Hopkins, P. B.; Robinson, R. H. Biophys. J. 1999, 77, 3256-3276. (40) Yeh, I.-C.; Hummer, G. J. Am. Chem. Soc. 2002, 124, 6563-6568. (41) Allison, S.; Austin, R.; Hogan, M. J. Chem. Phys. 1989, 90, 38433854. (42) Allison, S.; Sorlie, S. S.; Pecora, R. Macromolecules 1990, 23, 1110-1118. (43) Naimushin, A. N.; Fujimoto, B. S.; Schurr, J. M. Biophys. J. 2000, 78, 1498-1518. (44) Hagerman, P. J.; Zimm, B. H. Biopolymers 1981, 20, 1481-1502. (45) Yeh, I.-C.; Hummer, G. Biophys. J. 2004, 86, 681-9. (46) Allison, S. A.; Schurr, J. M. Macromolecules 1997, 30, 71317142. (47) Trifonov, E. N.; Tan, R. K. Z.; Harvey, S. C. In DNA Bending and CurVature; Olson, W. K., Sarma, M. H., Sarma, R. H., Sundaralingam, M., Eds.; Adenine Press: New York, 1988; pp 243-253. (48) Mazur, A. K. J. Am. Chem. Soc. 2000, 122, 12778-12785. (49) Mazur, A. K.; Kamashev, D. E. Phys. ReV. E 2002, 66, 011917. (50) Kamashev, D. E.; Mazur, A. K. Biochemistry 2004, 43, 81608168.
