Kinetic and Thermodynamic DNA Elasticity at Micro- and Mesoscopic

J. Phys. Chem. B 2009, 113, 2077–2089

2077

Kinetic and Thermodynamic DNA Elasticity at Micro- and Mesoscopic Scales Alexey K. Mazur CNRS UPR9080, Institut de Biologie Physico-Chimique 13, rue Pierre et Marie Curie, Paris 75005, France ReceiVed: NoVember 10, 2008

Recent theoretical and experimental studies have suggested that the elastic behavior of the small-length doublehelical DNA does not correspond to the simple harmonic model. This article presents a thorough comparison of classical atom-level molecular dynamics (MD) and coarse-grained harmonic approximations. It is shown that the previously predicted duration of MD trajectories necessary for accurate assessment of DNA elasticity was significantly overestimated and that reliable estimates of elastic parameters can be obtained after a few tens of nanoseconds. The all-atom and coarse-grained ensembles were compared head-to-head, including the amplitudes and relaxation rates of internal fluctuations as well as translational diffusion. The computed diffusion rates were found to be similar, with good correspondence to experimental data. The torsional persistence length (PL) in MD agrees reasonably well with experiment, with the relaxation rate of twisting fluctuations corresponding well to the harmonic model. The bending PL also agrees reasonably well with experiment, but the corresponding relaxation rate is much higher than the harmonic approximation. For a tetradecamer DNA, the difference reaches 1 order of magnitude, with the bending dynamics faster than the twisting dynamics, in qualitative contrast to the coarse-grained model. The possible mechanisms of this anomalous behavior are discussed. Introduction The mechanical properties of the DNA double helix are important for many biological functions. The strand separation necessary for DNA transcription, replication, and repair is preceded, controlled, and accompanied by complex mechanical deformations of the intact double helix at length scales from a few to hundreds of base pair steps (bps). Long double helixes behave as continuous elastic rods described by the worm-like chain (WLC) theory;1,2 however, for a stiff polymer with a macroscopic persistence length (PL) of 50 nm, DNA makes surprisingly small loops,3 and the validity of the WLC model for short and medium-length segments is not certain.4-9 In living cells, DNA is always subjected to external macroscopic stress, which can affect its local elastic properties and cause conformational transitions.10,11 In addition, biological functions of DNA are somehow controlled by the sequence, which remains entirely puzzling even in extensively studied domains such as intrinsic sequence-directed curvature.12,13 All these and related problems represent major challenges for atom-level molecular modeling. Continuous improvements in force fields14-16 and simulation techniques17,18 now make possible free molecular dynamics (MD) simulations that reproduce conformational ensembles of DNA in reasonable agreement with experimental data.19,20 It is not known, however, whether the microscopic DNA dynamics predicted by MD agrees with experimental mechanical properties. The latter are very well described by the coarse-grained WLC model parametrized by direct fitting of the rheological and elastic properties to experimental data, with harmonic potentials assumed for bending, torsional, and stretching deformations.21-24 In contrast, the MD force fields are parametrized by using only spectroscopic and structural data for small molecules, as well as quantum mechanics calculations.19,25 The consistency relation between the two representations is difficult to check. The elastic parameters can be evaluated a posteriori, from ensembles of computed DNA conformations,26-28 but this task appears to be more complicated than it initially seemed.29,30

Reasonably good sampling can be obtained in MD only for short oligomers where the amplitudes of WLC fluctuations are small, and the DNA length, bend, and twist angles are needed with high precision. Taking such measurements in short fragments is problematic and prone to complex artifacts.29,31,32 Also, even for short DNA, it takes a relatively long time to visit the number of distinct conformational states necessary for the evaluation of elastic parameters.30 All of these issues and the possible breakdown of the WLC theory for medium-length DNA are interconnected and should be considered together. A recent report30 demonstrated that new insights in this direction can be obtained by directly comparing the atomistic MD model of DNA with its coarse-grained representations parametrized to reproduce the macroscopic elasticity. Brownian dynamics (BD) of discrete WLC models provides the best possible extrapolation of macroscopic elastic properties to small chain lengths where relevant experimental data are not available. These models were originally introduced by Allison et al.33,34 and were further developed by different authors.35,36 The BD method reproduces thermodynamic and kinetic parameters of long double helixes in good agreement with experiment.37,38 Analysis of MD trajectories usually gives complex patterns in which genuine WLC behavior is obscured by artifacts that might be due to either poor statistics or flaws in the methods of data processing. When the same methods are benchmarked against BD models with known elastic properties, methodological artifacts can be distinguished and corrected. One can also reveal the effects of poor statistics and estimate the required duration of MD trajectories. Earlier studies demonstrated that the statistics of bending MD fluctuations in short DNA qualitatively agree with the WLC theory.8,29,30 The apparent bending PL is somewhat larger than the experimental value, but the difference can be attributed to imperfect modeling of ionic conditions in MD. Unexpectedly, it was found that correct canonical ensembles of bent conformations are obtained in MD with trajectory durations significantly

10.1021/jp8098945 CCC: $40.75  2009 American Chemical Society Published on Web 01/23/2009

2078 J. Phys. Chem. B, Vol. 113, No. 7, 2009

Figure 1. Composite-bead model of a DNA double helix. A composite bead consists of four particles that define the local Cartesian frame. The particles are labeled as Pji, with the upper and lower indexes showing the number of the composite bead in the chain and the particle number within the bead, respectively.

below the estimates provided by BD models. The difference is due to surprisingly fast relaxation of bending fluctuations in atomistic MD. This condition is fortunate for convergence of computer simulations, but its physical origin is intriguing. Simple harmonic models can account for these data only by assuming nonphysically small values of DNA radius and/or solvent viscosity, which can be viewed as a divergence of kinetic and thermodynamic elasticities.30 On the other hand, indirect experimental observations by time-resolved Stokes shift (TRSS) spectroscopy revealed the presence of very slow relaxation in the dynamics of short DNA.39,40 These slow modes were originally attributed to solvent, but recent investigations suggest that they instead should be intrinsic in the DNA structure.41-43 The only reasonable candidates are bending and/or twisting, but this would correspond to WLC relaxation with an unphysically large DNA radius and/or solvent viscosity. The most intriguing are the physical mechanisms that could, in principle, produce very fast or very slow relaxation in the double helix. To gain further insight into these problems, we present a thorough re-examination of the correspondence between the apparent kinetic and thermodynamic elasticities in double-helical DNA as modeled by all-atom MD. To this end, we provide a head-to-head comparison of the atomistic and coarse-grained dynamics, including the amplitudes and the relaxation rates of all types of WLC fluctuations, as well as translational diffusion. We show that the previously predicted duration of MD trajectories necessary for accurate assessment of the macroscopic DNA elasticity was significantly overestimated and that accurate estimates of elastic parameters can be obtained after a few tens of nanoseconds. As expected, BD simulations reproduce the rate of translational diffusion in excellent agreement with experiment. MD simulations agree with these data reasonably well, suggesting that the effective solvent viscosity cannot be responsible for any apparent divergence of kinetic and thermodynamic elasticities. For twisting dynamics, we observe good agreement of the MD-measured and experimental torsional PLs and good correspondence of the twisting relaxation rate with the harmonic approximation. The measured bending PL also agrees reasonably well with experiment, but the corresponding relaxation rate is much higher than the harmonic approximation. For a tetradecamer DNA, the difference reaches 1 order of magnitude, with bending faster than twisting, in qualitative difference from the coarse-grained BD. We discuss some possible mechanisms of this anomalous behavior. Theory and Methods Composite-Bead Model of DNA. The DNA model used in BD simulations is displayed in Figure 1. A composite bead consists of four particles labeled P0i -P3i , where the upper index refers to the bead number in the chain. The beads are rigid, with 1-Å interparticle distances P0i P1i , P0i P2i , and P0i P3i and right

Mazur angles between the corresponding vectors. Particle P0i is used as the center of the local Cartesian frame; particles P1i , P2i , and P3i define the corresponding three orthogonal axes. Random forces and hydrodynamic interactions are applied only to P0i ; these particles are considered spherical with a constant radius a. Particles P1i -P3i are virtual; they provide the beads with a definite three-dimensional orientation and are used for applying angular potentials. The total energy of an N-bead chain is composed of bond stretching, torsion and bending terms as

h U) 2

N-1

∑ i)1

q (li - l0) + 2 2

N-1

∑ (φi - φ0)2 + Ub

(1)

i)1

All other interactions between the beads can be neglected because we consider only short stiff chains that cannot loop. Here, li is the distance between particles P0i and P0i+1, and l0 is its equilibrium value. In our calculations, l0 was either 5 or 0.33 nm. The first value is commonly used in BD simulations of long DNA.37 For the second value, one bead corresponds to one base pair. The torsion potential is applied to dihedral angles φi ) P1i P0i P0i+1P1i+1, with the equilibrium helical twist angle φ0. The bending energy, Ub, is defined in two alternative ways. In the first case, the bending is isotropic, and Ub is computed similarly to earlier simulations as

g 2

Ub )

N-2

∑ θi2

(2)

i)1

where the angular displacement is measured for angles θi ) P0i P0i+1P0i+2. The second case takes into account the fact that a single base pair step of the DNA double helix is stiff for bending in the plane of the backbone strands, but much softer in the perpendicular direction (bending toward the grooves).26,44 In addition, bending in the soft plane is often easier toward one of the grooves than in the opposite direction. All of these features are included in the expression N-1

Ub )

∑ i)1

with

gR π R 2 i 2

(

2

)

{

N-1

+

∑ i)1

gβ(βi) π βi 2 2

π g+, β > 2 gβ(β) ) π g-, β < 2

(

2

)

(3)

(4)

i+1 The “soft” and “stiff” bend angles are defined as βi ) Pi0Pi+1 0 P1 i i+1 i+1 and Ri ) P0P0 P2 , respectively. The chain is built so that vector P0i P3i of the ith bead always points toward P0i+1. To this end, after every translational BD displacement, each bead is properly rotated around its center. Such correcting rotations are conceptually similar to using projection operators for the treatment of holonomic constraints in MD.45 Projection operators replace reaction forces that act in linear subspaces orthogonal to particle movements without producing mechanical work. In the present algorithm, this condition is fulfilled only approximately because the correcting rotations might have nonzero projections on the gradients of

Kinetic and Thermodynamic DNA Elasticity

J. Phys. Chem. B, Vol. 113, No. 7, 2009 2079

torsional and bending potentials. These side effects are reduced by introducing explicit mechanical coupling between the translational and rotational movements of the beads. Particle P0i+1 is bound to vector P0i P3i , and it effectively belongs to the rigid body of the ith bead. Therefore, a pair of forces applied to particles P0i and P0i+1 should create an additional torque on the ith bead computed as

ti ) r × (Fi+1 - Fi)

1 Ti × r r2

(6)

and an opposite force on P0i . This coupling reduces the discrepancy between the theoretical and observed values of the torsional PLs, especially with the second form of the bending potentials. Interactions defined by eq 1 are treated as usual in molecular mechanics, that is, energy gradients are transformed to Cartesian forces applied to particles. This makes possible the use of conventional molecular mechanics codes and facilitates the addition of complex energy terms such as those in eqs 3 and 4. In the earlier implementations of this method, only the type of bending potential shown in eq 2 was used, and the torsional gradients were computed with respect to Euler angles.46-48 This difference is a matter of convenience and is not physically significant. The last composite bead in the linear chain differs from the rest in that vector P0NP3N is not constrained, and its orientation is arbitrarily chosen to be collinear to vector P0N-1P3N-1. The properties of the Nth bead are slightly perturbed as a result, but this end effect is noticeable only for very small fragments. BD Simulations. BD simulations were carried out using an algorithm developed previously and adapted for DNA by several authors.35,36,49,50 For translational movements, only particles P0i are considered. In the original (first-order) procedure by Ermak and McCammon,49 translational displacements on time step ∆t are computed as

∆ri )

∑ j

Dij0 Fj0 ∆t + Ri(∆t) kT

(7)

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 kT is the Boltzmann factor. The zero superscript refers to the state before the time step. The Gaussian random force, Ri(∆t), is sampled from a multivariate distribution with

〈Ri(∆t) Rj(∆t)〉 ) 2Dij0 ∆t

{

[(

) (

[(

)

]

)]

(5)

where r ) P0i P0i+1. Similarly, a torque Ti on the ith bead creates an additional force applied to particle P0i+1

fi+1 )

kT I 6πηa rijrij rijrij kT 2a2 1 I + 2 + 2 I - 2 , if rij > 2a 8πηrij rij rij 3 rij Dij ) r r r kT 9 ij 3 ij ij if rij e 2a 1I+ , 6πηa 32 a 32 arij (9)

Dii )

(8)

simulated by 3N independent standard Gaussian random generators as proposed by Ermak and McCammon.49 As in earlier DNA simulations,36 the tensor Dij is computed from interparticle vectors using the Rotne-Prager formulas51

In contrast to eq 7, here, subscripts i and j refer to particle numbers, with vector rij ) P0i P0j and rij ) |rij|. The bold capital letters denote three-dimensional tensors, notably, I is a unit three-dimensional tensor, and tensor multiplication of vectors is assumed throughout. The solvent viscosity, η, is that of water near 25 °C (0.89 mPa · s), which is close to the conditions used in MD simulations. In the previous report,30 a slightly larger value was used (1 mPa · s). The particle radius, a ) 1.78 nm, corresponds to the value used in earlier BD simulations of DNA.37 The torsional BD movements were computed according to the original version of the algorithm developed by Allison et al.46,52 In this implementation, rotation of composite beads occurs around vectors P0i P3i according to the formula

∆φi )

∆t R 0 D Ti + Si(∆t) kT

(10)

where ∆φi denotes the rotation increment at one time step and Ti is the component of the total torque applied around vector P0i P3i . The scalar character of the torque implies that all other components are compensated by reaction forces. The random rotations, Si, are computed as

Si ) √2DR∆tXi where DR is the rotational diffusion coefficient and Xi is a Gaussian random value with zero mean and variance equal to 1. The diffusion coefficient, DR, is a scalar, which means that hydrodynamic interactions are negligible.46 In the simplest case, DR is computed assuming that the beads are joined by cylindric sticks of radius b articulated at particles P0, which gives46

DR )

kT 4πηl0b2

The value of b is chosen independently of the particle radius a in eq 9, which allows independent fitting of translational and rotational diffusion rates to experimental data. In the present calculations, we used b ) 1.2 nm as in the early study by Allison et al.46 To increase the time step, we used the second-order version of this algorithm35,50 simplified as proposed by Jian et al.36 At the beginning, the chain is found in a correct configuration, with i i all particles Pi+1 0 on vectors P0P3. In this state, all potential forces and torques are computed, and then the coupling relations in eqs 5 and 6 are applied. The resulting total torques are projected onto vectors P0i P3i , and the predictor BD movement is carried out by using eq 10, followed by translations according to eq 7 and adjustments of the rotations of the beads to align vectors P0i P3i to the new positions of P0i+1. At this state, the potential forces and torques are computed again, and the corrector BD step is made by using the average of the two sets of potential forces and torques. The random increments for the correction are taken from the predictor step. The most time-consuming


Mazur

phase of the algorithm is the evaluation of the hydrodynamic tensor Dij. Following earlier experience,36 the components of this tensor are updated once every 10 steps of BD. The time step can be limited by the bending, twisting, or stretching motions, depending on the force-field parameters and the link length l0. Under our conditions, it was mainly limited by twisting. With l0 ) 5 nm, the time step was 25 ps as in earlier simulations.30 With l0 reduced and coefficient q in eq 1 increased to reach agreement with MD simulations, the time step should be scaled approximately as l02/q. Atomistic MD Simulations. The classical MD simulations were carried out using previously described protocols.29,30 The DNA duplexes were modeled in explicit aqueous environment with neutralizing amounts of Na+ ions. The AMBER98 forcefield parameters14,53 were used with the rigid TIP3P water model54 under NVT ensemble conditions at T ≈ 293 K, with a water density of around 0.997. To increase the time step, MD simulations were carried out by the internal coordinate method (ICMD),45,55,56 with the internal DNA mobility limited to essential degrees of freedom and rotation of water molecules and internal DNA groups including only hydrogen atoms slowed down by weighting of the corresponding inertia tensors.57,58 The double-helical DNA was modeled with all backbone torsions, free bond angles in the sugar rings, and rigid bases and phosphate groups. The effect of these constraints is insignificant, as was previously checked through comparisons with standard Cartesian dynamics.29,57 The time step was 0.01 ps, and the DNA conformation were saved every 2.5 or 5 ps. As in the previous report,30 we considered DNA duplexes with the AT alternating sequence. These molecules are homopolymers of AT units, so they cannot have distinguished asymmetric structures such as static bends. For consistency, we considered the results of 28-ns MD simulations of a 25-mer DNA fragment29 previously analyzed by different methods. These data, however, did not allow accurate evaluation of diffusion coefficients. To improve sampling and increase the ratio of the trajectory duration to the maximum relaxation time, an additional 100-ns trajectory was computed for a tetradecamer duplex. These calculations were carried out in a 6.2-nm cubic cell with periodic boundaries using the same protocol. Evaluation of Elastic Parameters. The DNA elasticity is conveniently characterized by three persistence lengths (PLs) corresponding to bending, twisting, and stretching that we denote here as lb, lt, and ls, respectively. Parameters g, q, and h in eqs 1-3 were chosen so that the theoretical values of the PLs were close to either experimental data or MD simulations used in comparisons. In the limit of small deviations, the PLs can be computed using the general theory of fluctuations1 as

gl0 kT ql0 lt ) kT hl0 3.4 nm ls ) kT 2π lb )

(

2

)

(11)

Dx(L) )

L lx

where L is the DNA length and x stands for b, t, or s. The WLC deviations Dx(L) can be computed from appropriate canonical averages as

Db(L) ) - ln(〈cos[θ(L)]〉) Dt(L) ) var[φ(L)] 2π 2 var[L] Ds(L) ) 3.4 nm

(

)

(13)

where θ(L) and φ(L) are the bend angle and the overall winding of DNA, respectively; the angular brackets denote canonical averaging; and var[ · ] represents the variance of the variable in brackets. For various reasons as discussed in the text, the bending and torsional PLs measured a posteriori for BD simulations can differ slightly from the theoretical predictions provided by eq 11. For the analysis of MD results, equivalent BD models were finetuned by empirically adjusting parameters q and g to make the measured PLs close to the corresponding MD values. For the composite-bead DNA models, the bend angles are measured directly, whereas the DNA length and the overall twisting are obtained by summing the instantaneous link lengths li and twist angles φi along the chain. For DNA conformations sampled by atomistic MD, the corresponding measures were taken using the program 3DNA59 and an in-house implementation29 of the algorithm Curves.31 Unlike in the previous report,30 the latter method was used with its standard parameters. Time Autocorrelation Functions and Diffusion Coefficients. The relaxation rates are characterized by normalized time autocorrelation functions Cb(t), Ct(t), and Cs(t) and the corresponding effective relaxation times τb, τt, and τs for bending, twisting, and stretching, respectively. The autocorrelation functions are computed from the time series of appropriate parameters as

C(t) )

C˜(t) , C˜(0)

˜ (t) ) 〈x(τ + t) x(τ)〉τ - 〈x〉2 C

where x stands for the overall length, bend angle, or winding of a DNA fragment. The effective relaxation times are measured at the 1/e level, that is, by using the condition C(τ) ) 1/e. The statistical analysis of translational diffusion was carried out by averaging over all relevant subtrajectories of a given duration starting at different moments. The diffusion coefficients, D, were evaluated from linear fits of 〈∆r2〉(t) to the Einstein-Smoluchowski equation for small t, where ∆r is the displacement vector measured for the DNA center of mass. The errors were evaluated by dividing every trajectory into three equal parts and taking the average and standard deviation of the three independent measures. The diffusion coefficients of DNA and water in MD were obtained from the 100ns trajectory of the tetradecamer double helix. The duration of the BD trajectories was at least 25 µs. The experimental values of diffusion coefficients were estimated using the empirical relationship

DdsDNA ) 7.73 × 10-6 × N-0.672 (cm-2 s-1) The same parameters can be extracted from canonical conformational ensembles produced by BD or MD simulations using the WLC theory, which provides linear relationships of the following form

(12)

(14)

where N is the number of bps in a DNA fragment and dsDNA is double-stranded DNA. Equation 14 was obtained by Stellwagen et al.60 from a fit to experimental data. The values of N



Figure 2. Selection of an optimal sampling interval for evaluating the torsional PL. Data from 100 independent BD trajectories of identical duration are analyzed, gradually reducing the time step between the saved configurations. Each trajectory gives one count of the measured PL, and 100 independent counts are distributed according to the probability distributions shown in the figure. The vertical gray bars show the scattering of the corresponding average values. A minimal chain of three beads (two torsions) was used. The duration of trajectories was 12.5 µs. The sampling intervals are indicated in the figure in nanoseconds. The theoretical value of the torsion PL is 75 nm.

for eq 14 were computed from the DNA length assuming a bps length of 0.33 nm. Results Coarse-Grained Torsional Dynamics. We begin with a separate analysis of torsional fluctuations, which were neglected for simplicity in our previous report.30 The standard compositebead model is used with the link length l0 ) 5 nm and the bending potential defined by eq 2.48 The PL values estimated by using eq 11 are lb ) 50 nm, lt ) 75 nm, and ls ) 100 nm. Small chains of three to eight composite beads are considered. Calculation of a PL value from data accumulated during a single trajectory can be viewed as a stochastic process and characterized by the probability of hitting the thermodynamic limit within a certain error. The corresponding probability distributions are obtained by running a large number of BD trajectories for the same system. The distributions depend physically on the trajectory duration and the corresponding relaxation time, but in practice, they are also affected by the interval between the saved chain configurations (sampling interval). To avoid spurious effects, the ratio of the sampling interval to the minimum relaxation time should be sufficiently small. Figure 2 checks this condition for torsional dynamics. The distributions are bellshaped, with the average PL values reasonably close to the analytical prediction, which indicates that the BD simulations correctly sample from the canonical ensemble of configurations. The width of the distributions decreases with reduced sampling interval until its value becomes smaller than 0.2 ns. For bending and stretching, the appropriate intervals are 2.5 and 0.1 ns, respectively.30 These intervals are used for other tests with this model. The accuracy of the measured PLs naturally improves with the duration of the trajectories. This dependence is critical for atomistic MD simulations, and we consider it in the following section by using a refined discrete WLC model. For the present model, such tests were also run, and the results are shown in Figure 1S of the Supporting Information. For bending and stretching, the results do not differ from those in the previous report.30 The rate of torsional dynamics is intermediate between bending and stretching; notably, a 90% probability of getting lt with 10% error is reached with a trajectory length of about 0.25 µs. The measured values of bending and stretching PLs always slightly exceed the theoretical estimates.30 This is attributed to

Figure 3. Comparison of the effects of increased fragment length and increased trajectory duration on the accuracy of the evaluated PLs. The probability distributions were obtained for chains with two and six elements, with one element corresponding to one angle, torsion, or link for bending, twisting, or stretching, respectively. The plots are marked by two-number codes of the form P1/P2, where P1 gives the number of elements and P2 indicates the length of the trajectory in nanoseconds. The primed distributions were obtained by considering only the first element in the chain. In all other cases, all internal elements were analyzed. The vertical gray bar shows the scattering of the average PL values.

a small shift in the effective temperature caused by the discretization error of the BD algorithm. For both lb and ls, the shift is about 2% and can be reduced by decreasing the time step. In contrast, the torsional PL appears to be underestimated with respect to the theoretical prediction (see Figures 2 and 3 and Figure 1S in the Supporting Information). The discrepancy can be reduced by using smaller time steps, but this effect is also due to the coupling between twisting and bending, which is intrinsic in the composite-bead chain model. In Figure 1, note that some translations of particle P03 followed by adjusting rotations of vector P20P23 change the torsion of P11P10P20P21; therefore, bending fluctuations contribute to twisting and reduce the effective torsional PL. The discrepancy is reduced when the bending stiffness is increased, but with the bending flexibility close to experiment, accurate prediction of the torsion PL with eq 11 is not always possible, and the value of q in eq 1 should be fine-tuned empirically. The previous report30 revealed a curious qualitative difference between bending and stretching dynamics of internal fragments in medium-length DNA. In time-limited ensembles, local stretching fluctuations appear to be statistically independent. In contrast, local bending fluctuations are constrained by the global bend of the whole molecule. This effect is attributed to external friction, and it dictates different optimal strategies for measuring

2082 J. Phys. Chem. B, Vol. 113, No. 7, 2009 the corresponding PLs. For stretching, good sampling can be achieved by running short trajectories of long DNA and then counting fluctuations in all internal fragments. In contrast, for bending PL, the accuracy can be improved only by increasing the ratio of the trajectory length to the relaxation time, that is, by taking possibly short DNA fragments and simulating their free dynamics for maximal time. These properties are checked for torsional fluctuations in Figure 3, with bending and stretching data also included for comparison. Plots shown by open symbols demonstrate the narrowing of the probability distributions obtained by increasing the trajectory length by a factor of 3. For instance, the notation 2/250 in the middle panel means that the corresponding data were obtained from simulations of a chain with two torsions (three beads) run for 250 ns. This distribution is clearly wider than that labeled 2/750. The solid symbols mark the probability distributions obtained for 3-times-longer chains. Distribution 6/250 was obtained by considering all internal fragments, i.e., the number of chain configurations is identical to that in the narrow distribution 2/750. Qualitatively, torsional fluctuations appear to be intermediate between stretching and bending. With the number of configurations increased by elongating the chain, the accuracy is improved, but the effect is smaller than with the equivalent increase obtained by continuing the trajectory. The primed probability distribution 6/3′ obtained by counting only the first torsion is wider than the 6/3 distribution, which is qualitatively similar to stretching but differs from bending. Fine-Grained Discrete WLC Model. The discrete WLC models are designed for simulations of long DNA and employ large l0 values because the overall computational load grows as the square of the number of beads. These models behave as true WLCs for L . l0, but for the short fragments required for comparison with MD, it is preferable to reduce the link length to one base pair step. For convenience, in the subsequent discussion, the coarse-grained model with l0 ) 5 nm is called BD1, and the fine-grained representation with reduced l0 is referred to as BD2. The BD2 model used in this section has l0 ) 0.33 nm and φ0 ) 36°. Theoretically, a smaller link length can be used provided that the bead radius, a, in eq 9 is properly adjusted. The value of 1.78 nm for BD1 was obtained by fitting the computed and experimental sedimentation coefficients for a random-coil ensemble of chain configurations.36,61 We need the BD2 model for comparisons with BD1 and atomistic MD, and we will apply it to short DNA fragments. For these purposes, it is reasonable to adjust the bead radius by fitting the diffusion rates of identical DNA fragments modeled by BD1 and BD2. To begin, we compared the diffusion coefficients of the two models with the beads unchanged. The translational mobility was analyzed for the centers of mass of the DNA fragments and compared with experimental data as described in Theory and Methods. The results are reported in Table 1. It can be seen that the difference between the coarse- and fine-grained models is surprisingly small. Moreover, the latter is more mobile for small lengths, which was not expected. The higher mobility of BD2 can be explained by the properties of the Rotne-Prager diffusion tensor equations (eqs 9). With an increase in the number of beads on the same string, the friction coefficient grows until the beads begin to overlap. After that, it should decrease because the surface area of the limiting cylinder is smaller than that of the contacting beads. In BD2, the beads are very close, and the enveloping surface is almost cylindrical. The contact surface for BD1 evaluated directly is larger than the cylindrical limit, which explains its lower mobility. The data

Mazur TABLE 1: Diffusion Coefficients (D) of a Few Short DNA Fragments According to Different Simulation Modes and Experimental Dataa L ≈ 5 nm b

L ≈ 10 nm

L ≈ 15 nm

L ≈ 20 nm

N

N

method

N

D

N

BD1 BD2 MD exp

2 16 14 15.1

9.4 (0.0) 10.3 (0.0) 7.4 (0.7) 12.4

3 31 30.3

D

7.4 (0.0) 4 7.7 (0.0) 46 7.8 45.5

D

6.1 (0.0) 5 6.4 (0.0) 62 5.9 60.6

D 5.4 (0.0) 5.4 (0.0) 4.9

a Diffusion coefficients measured in units of 10-7 cm2/s. Values obtained from simulations are accompanied by standard deviations shown in parentheses. b Method denoted as follows: BD1, discrete WLC model with l0 ) 5 nm; BD2, fine-grained model with l0 ) 0.33 nm; MD, molecular dynamics simulations as described in Theory and Methods; exp, experimental estimates with the empirical formula of Stellwagen et al.60

in Table 1 indicate that, for our purposes, the adjustment of the bead radius is not essential, and to simplify comparisons, we prefer to use the same value in both models. Table 1 also reveals the excellent agreement of BD simulations with experimental diffusion coefficients obtained from electrophoretic mobility of small DNA. This means that the sedimentation coefficients originally used for fitting the bead radius61 are consistent with the recent data from Stellwagen et al.60 used in Table 1. The agreement is even better than shown here because the empirical power law given by eq 14 slightly overestimates the diffusion rates for N < 20 and underestimates them for 20 < N < 100.60 The time autocorrelation functions of the two models are compared in Figure 4. The stretching relaxation is very close to exponential, with the reduced link length producing virtually no effect on the relaxation times. For twisting and bending, the dynamics are significantly accelerated due to local modes, which is seen as deviations of the dashed plots from simple exponents. The amplitude of this effect is stronger for bending, and its relative weight is larger for small chain lengths. For 10-nm chains, the reduction in τb reaches a factor of 2. This is probably close to the theoretical maximum, because l0 ) 0.33 nm is already at the limit of physical validity of such modeling49 and reduction to 0.2 and 0.1 nm gives virtually no acceleration. All relaxation times fall with the chain length, but for stretching the decrease is much slower, which is additionally accentuated by the effect of local modes. As a result, for short fragments, the torsional dynamics becomes the fastest, with τt ≈ τs and τb only a few times larger. The accelerated relaxation should favorably affect the accuracy of elastic parameters evaluated from short trajectories. The estimates shown in Figure 5 were obtained from BD2 simulations of a tetradecamer DNA fragment. This length is appropriate for comparison with MD simulations because it is just sufficient for probing deformations in one helical turn. The results in Figure 5 were obtained similarly to the BD1 calculations described above and in ref 30. First, bell-shaped probability distributions were obtained by running a large number of trajectories of the same duration. The resulting distributions were integrated over appropriate error intervals to produce the plots shown in Figure 5. The durations of trajectories selected for this figure were chosen so that the corresponding 10% error intervals were hit with a probability of 90% or larger. The top panel of Figure 5 shows that, if the WLC fluctuations are measured for 1 bps, 10-ns trajectories should be largely sufficient for computing all elastic parameters with very high accuracy. This explains the good apparent convergence of



Figure 5. Expected accuracy of evaluation of the bending, torsional, and stretching PLs for simulations of a tetradecamer DNA. The plots show the probability that a PL value computed from a single trajectory appears within a given error interval from the average. The probability distributions were computed similarly to Figure 1S, that is, by running a large number of trajectories of a few selected durations with a fixed small sampling interval. The PL values were estimated with eq 13 by using the fluctuations in 1 bps (top panel), the initial slopes of the D(L) plots (middle panel), or the fluctuations in one helical turn (11 bps, bottom panel). Figure 4. Comparison of time autocorrelation functions for the BD1 (red solid lines) and BD2 (blue dashed lines) WLC models of DNA. The upper red numbers indicate the chain lengths for the red curves. The lower blue numbers indicate the corresponding lengths of DNA fragments in nanometers. In all cases, the apparent PLs were 52.5 ( 0.5 nm for bending, 71.5 ( 1 nm for twisting, and 102 ( 1 nm for stretching.

averages in early MD simulations.28,29 We will see below that it is sometimes necessary to measure the PLs from the slopes of the D(L) plots. The middle panel of Figure 5 demonstrates that, in this case, the accuracy is less good, but a few tens of nanoseconds should be sufficient. In practice, however, the bps parameters extracted from real conformations are strongly affected by small-scale atom motions, and their fluctuations are not always relevant to the WLC behavior of long chains. As shown in the next section, one helical turn is the most appropriate DNA length for measuring elastic parameters in MD. The bottom panel of Figure 5 shows that, in this case, about 50 ns is necessary for evaluating the torsional and stretching PLs, and 200 ns should be sufficient for bending. These estimates do not take into account the possible difference in the corresponding relaxation rates in BD and MD. For stretching, the data in Figure 5 agree reasonably well with the earlier report,30 but for bending, the necessary duration of trajectories appears to be about an order of magnitude shorter. Figure 4 indicates that the acceleration due to the fine-grained discretization can account for only part of this difference. The

rest of the discrepancy is caused by the earlier erroneous assignment of the true DNA length corresponding to the minimal chains in the BD1 model. The correct assignment follows from Figure 4 because, in the long-time limit, the slowest exponential components dominate and the corresponding lines for BD1 and BD2 become nearly parallel. It is seen that the minimal threebead chain in BD1 approximates a 15-nm DNA for twisting and a 10-nm DNA for stretching and bending. Only for stretching is this correspondence intuitively obvious. A threebead chain has two torsions, but its twisting friction corresponds to three linked cylinders. In our previous report,30 we took for granted that, in terms of bending dynamics, a three-bead BD1 chain approximated a 5-nm DNA fragment because only in this case are the chain length dependences smoothly extrapolated to theoretical limits at N ) 0 (this model has one angle per 5 nm on average). This logic appears to be wrong. In fact, there is no BD1 representation for bending in 5-nm DNA. A threebead BD1 chain approximates a 10-nm double helix, and it should be compared with atomistic MD of 30- rather than 15bp fragments, which additionally reduces the earlier difference in τb by a factor of 3. However, even the combined effect of the misassigned length and the local bending modes is not sufficient to obtain τb values close to those from atomistic MD. The bending relaxation in MD remains much faster, and we try to evaluate this difference more accurately in the following sections.


Figure 6. Effects of anisotropic bending on characteristic chain length dependences. BD simulations were carried out for two chains with anisotropic bending potentials defined by eqs 3 and 4 with gR/gβ ) 20. In the first case (open circles), the rolling is symmetric (g+ ) g-), and in the second case (closed circles), it is strongly asymmetric (g+/g- ) 25). Each chain consists of 36 beads, with l0 ) 0.33 nm and φ0 ) 36°. The chain length dependences were obtained by averaging over all internal fragments. The bending PLs measured for one helical turn are 48.4 and 50.8 nm for the symmetric and asymmetric cases, respectively. (Upper panel) Bending WLC deviations. The solid line corresponds to an ideal WLC with lb ) 50.8 nm. (Lower panel) Bending relaxation times.

Effects of Bending Anisotropy. The BD2 model can be made more realistic by taking into account the bending anisotropy of individual base pair steps. A stack of two base pairs has two qualitatively different orthogonal planes parallel to the helical axis. The first plane goes along the base pair approximately parallel to the hydrogen bonds. Bending in this plane is difficult because such deformations (tilting) are hindered by the DNA backbone. Bends in the perpendicular plane (rolling) deviate the helical axis toward the grooves, which costs much less energy. Depending on the sequence, bending toward the major groove (roll > 0) can be easier than in the opposite direction (roll < 0) or vice versa.26,44 The bending potential defined by eqs 3 and 4 describes this type of bending anisotropy, with gR and gβ terms corresponding to tilting and rolling, respectively. The effect of such anisotropy on the WLC behavior is analyzed in Figures 6 and 7. The results are shown only for bending because torsional and stretching dynamics are not affected. It is seen in Figure 6 that the strongly dissimilar energies of tilting and rolling have virtually no effect on the correspondence of this model with the WLC theory. In contrast, the asymmetry of rolling induces strong perturbations seen as regular oscillations of the Db(L) plots in phase with the helical winding defined by the twist angle φ0. This observation is important because very similar behavior was previously observed in all-atom MD.29 The oscillating Db(L) plot indicates that the chain of P0P3 vectors on average forms a helix, which is readily rationalized by noting that, with g+ * g- in eq 4, the canonical average 〈β〉 * π/2 at any nonzero temperature. In fragments including integral numbers of helical turns, the asymmetry of rolling appears to be approximately compensated by helical winding. As shown in the upper panel of Figure 6, fragments with integral numbers of helical turns produce plot points found approximately on a single straight line with zero intercept. This line corresponds

Mazur

Figure 7. Probability distributions of bend angles for different internal fragments of the two anisotropic chains considered in Figure 6. Data counts were accumulated in 80 windows evenly spaced on 0.6 < cos θ < 1. The left and middle panels show the distributions for N ) 3, 5, 7, 9, 11, and 13 (bottom to top) for the chains with symmetric and asymmetric rolling, respectively. The right panel shows the distributions for N ) 11, 21, and 31 (integral number of helical turns) for the model with asymmetric rolling. The coordinates were chosen to linearize the small-angle WLC distributions.8 Each distribution is accompanied by a theoretical straight line corresponding to the ideal WLC with an equivalent value of the bending PL, that is, lb ) 48.4 and lb ) 50.8 nm for symmetric and asymmetric rolling, respectively. All plots are consecutively shifted by 3.0 for clarity.

to an ideal WLC with the value of the bending PL equal to that of an infinitely long chain with such anisotropy. The lower panel in Figure 6 shows that the bending relaxation times oscillate with the same phase, which means that not only the average shapes but also the dynamic properties are affected. Notably, fragments including half-integer numbers of helical turns such as pentamers or pentadecamers appear to be considerably softer than they should according to the measured bending PL. The pattern in the upper panel of Figure 6 demonstrates that the PL values measured for 11-, 21-, and 31-mers, i.e., fragments with integral numbers of helical turns, provide increasingly accurate lower estimates of the asymptotic PL value. Such fragments can be used for the evaluation of elastic parameters even in the presence of strong asymmetric anisotropy. The accuracy is better for longer chains because the directional compensation improves with helical rotation, but for practical purposes, one helical turn should already be sufficient. The effect of directional compensation is confirmed in Figure 7, which compares the bend-angle distributions in the two models with that of the WLC theory. The asymmetric rolling causes various and sometimes very strong perturbations in the shapes of the bend-angle distributions. However, with integral numbers of helical turns, the distributions agree well with the WLC theory both qualitatively and quantitatively. Notably, the plots for 11-, 21-, and 31-mer fragments shown in the right panel overlap with the straight lines corresponding to ideal WLC distributions with the bending PL equal to that measured for one helical turn. The left panel demonstrates that, with symmetric rolling, the anisotropy of bending also perturbs the shapes of the bendangle distributions of short chains. These perturbations disappear for N > 5, that is, one-half of the helical turn. Atomistic MD. Comparison of MD trajectories with the properties of long DNA observed in experiments involves two steps. First, the experimental observations are extrapolated to short chain lengths accessible in atomistic MD. This is achieved by using the WLC theory and BD simulations considered in the previous sections. Second, parameters corresponding to the coarse-grained DNA representation (that is, bending, twisting, and stretching) should be evaluated for DNA conformations produced by MD. Such analysis of DNA structures involves a


Figure 8. Comparison bending dynamics of a 25-mer DNA fragment according to atomistic and coarse-grained simulations. Open squares show the results of BD simulations of a 25-mer composite-bead chain with anisotropic bending. The relative anisotropy of bending is as in Figure 6. Open and closed circles show the results of all-atom MD simulations analyzed by programs 3DNA and Curves, respectively. The measured bending PL for the BD model is lb ) 239 bps. (Upper panel) Bending WLC deviations. (Lower panel) Bending relaxation times. The left and right Y axes show the scales for the atomistic and coarsegrained models, respectively. The arrow points to the axis corresponding to the BD plot.

few controversial issues (reviewed in ref 32) that go beyond our study. For reliability, we apply two different methods most often used in recent years and concentrate on the aspects where the two methods are convergent. The pattern of bending fluctuations in 25-mer AT-alternating DNA fragments is analyzed in the upper panel of Figure 8. For convenience, two plots of bending WLC deviations produced from MD data by the programs 3DNA and Curves are supplemented by the results of BD simulations of a fine-grained BD model with a similar bending PL. Note the remarkable qualitative similarity between this picture and that shown in Figure 6 for anisotropic rolling. All plots exhibit regular oscillations corresponding to superhelical winding. The algorithm of Curves is designed to suppress any supercoiling of the same period as the double helix.29,31 It can be seen that the supercoiling is reduced but not canceled. Nevertheless, the bending PLs measured for one helical turn by the two methods appear to be remarkably similar. A small discrepancy seen for two helical turns is attributable to insufficient sampling. These conclusions are confirmed in Figure 9, where the bend-angle distributions in fragments of different length are compared with those obtained from the WLC theory. The pattern here is also similar to that in Figure 7. With both methods, the distributions obtained for one and two helical turns agree with the theory much better than with other lengths. The best-fit lb value corresponding to these plots appears to be somewhat higher than that of the BD model used as a reference in Figures 8 and 9, but the difference is not large. The lower panel in Figure 8 displays the chain length profiles of bending relaxation times. Qualitatively, this pattern resembles the lower panel of Figure 6. The relaxation times measured with Curves and 3DNA are similar, although the amplitudes of the oscillations differ in correspondence with the upper panel. At the same time, the absolute rate of bending fluctuations in MD is much larger than in the equivalent fine-grained BD model.


Figure 9. Probability distributions of bend angles in DNA fragments of different lengths obtained by atomistic MD. For each chain length, the distributions result from analysis of all internal fragments of a 25mer double helix. The DNA structures were processed by the programs Curves (left panel) and 3DNA (right panel). The distributions are shown for N ) 5, 9, 13, 17, and 25 (bottom to top). Data counts were accumulated in 80 windows evenly spaced on 0.6 < cos θ < 1. The coordinates were chosen to linearize the small-angle WLC distributions.8 The straight lines show the theoretical WLC distributions for lb ) 239 bps corresponding to the BD model in Figure 8. All plots are consecutively shifted by 3.0 for clarity.

Figure 10. Comparison torsional dynamics of a 25-mer DNA fragment according to atomistic and coarse-grained simulations. The three model systems and the notation are described in the caption of Figure 8. The measured value of the torsional PL for the BD model is lt ) 369 bps. (Upper panel) Torsional WLC deviations. (Lower panel) Torsional relaxation times.

This difference is a factor of 6-8, which is smaller than the earlier estimate30 but still significant. A similar analysis of torsional fluctuations is presented in Figure 10. In the upper panel, the plots of torsional WLC deviations obtained with two methods both exhibit nonzero intercepts in qualitative disagreement with the WLC theory. At the same time, they are nearly linear and parallel for small lengths where the sampling is best. The upward deviations observed for longer chains are partially attributable to reduced sampling. The initial slope of these plots is a convergent feature that is well reproduced with MD of a tetradecamer ATalternating duplex considered below. In all of these cases, linear regression over the first five points gives lt of about 370 base pair steps (bps) or 130 nm. We take this value as the best current estimate of the asymptotic lt value for long chains. A more


Figure 11. Comparison stretching dynamics of a 25-mer DNA fragment according to atomistic and coarse-grained simulations. The three model systems and the notation are described in the caption of Figure 8. The measured value of the stretching PL for the BD model is ls ) 303 bps.

detailed analysis of the Dt(L) plots in Figure 10 is beyond the scope of this study. The discrepancy from the WLC theory might be partially due to physical reasons; however, the appreciable difference between the Dt(L) plots of Curves and 3DNA indicates that technical details of the two methods contribute as well. A qualitatively similar difference is observed for the longer MD of a tetradecamer (see Figure 2S). The lower panel in Figure 10 compares the torsional relaxation times. It can be seen that the two methods produce very similar saturating profiles, with good quantitative agreement. In contrast to bending, the rate of torsional relaxation in MD appears to be similar to that of the equivalent fine-grained BD model. This observation indicates that the strikingly fast bending dynamics in MD does not result from overall acceleration due to, for instance, underestimated viscosity of the water model. Finally, Figure 11 exhibits a similar comparative analysis of fluctuations of the chain length. In this case, Curves and 3DNA give strongly divergent results. The length of 1 bps measured by 3DNA is about 10% larger than with Curves, which is not surprising in the context of earlier debates concerning these two methods.32 This large difference gives a systematic discrepancy in the measured values of PLs, and that is why the DNA length in Figures 8-11 had to be expressed in base pair steps. At the same time, the fluctuations measured with Curves appear to be much larger, and accordingly, the apparent stretching PL is smaller. The origin of these intriguing observations is rooted in properties of the corresponding methods and is outside the scope of this study. In Figures 8-11, we use raw MD data obtained some time ago,29 which is advantageous for comparisons with earlier conclusions. However, these data did not allow evaluation of the translational mobility of DNA in order to take into account the possible discrepancy in water viscosity as modeled in MD and BD. The foregoing results also suggest that one helical turn is already sufficient for probing the WLC properties of long DNA and that one can take a shorter DNA fragment in MD to reduce the computational load and provide a better ratio of the trajectory length to relevant relaxation times. Based on these considerations, a 100-ns trajectory was computed for a tetradecamer AT-alternating DNA duplex as described in Theory and Methods. It was expected that the rate of translational diffusion of DNA was overestimated in MD as a result of the known low viscosity of the TIP3P water model,62 which could partially account for the strikingly fast bending DNA dynamics. Surprisingly, the translational diffusion coefficient of the DNA center of mass shown in Table 1 appears below that in BD and experiment.

Mazur

Figure 12. Comparison of time autocorrelation functions for atomistic MD and fine-grained BD simulations of a tetradecamer DNA fragment. The plots corresponding to bending and twisting fluctuations in one helical turn are shown in red and blue, respectively. The curves shown by solid circles feature the results of MD simulations (processed with the program 3DNA). The results of BD simulations with standard viscosity are shown by solid lines. The dashed lines display the results of BD simulations with viscosity increased to make the translational diffusion coefficients similar in MD and BD.

This unexpected result is attributed to the specific conditions employed in our simulations. Periodic boundary conditions, with Ewald summation of electrostatic interactions, reduce the rate of translational diffusion of water and solutes through very longrange hydrodynamic interactions between periodic images.63 In addition, the inertia weighting used for increasing the time step58 slows water rotation and indirectly affects the rate of translational diffusion of water as well as ions. All of these details are not essential for our present purposes. The results reported in Table 1 clearly show that the rate of translational DNA diffusion in our MD simulations is somewhat underestimated with respect to BD and experiment; therefore, the effective viscosity of water cannot be responsible for the anomalously rapid bending dynamics in MD. The chain length dependences of the WLC deviations are qualitatively similar to those of a 25-mer (see Figure 2S). The lb value measured for one helical turn is about 205 bps (72 nm), whereas the asymptotic lt value estimated from the slope of the Dt(L) plot is about 385 bps (135 nm). The two values are roughly 35% larger than the corresponding experimental estimates,22,64 which is highly satisfactory for all-atom simulations using a universal force field. The time autocorrelation functions for bending and twisting are analyzed in Figure 12. The anomalous character of bending in MD is evident from two observations. First, MD relaxation is faster for bending than for twisting, which is qualitatively different from the BD model with the same elasticity. Second, the absolute rate of bending relaxation in MD is very high, whereas for twisting, it is even lower than in BD. The discrepancy in the rates of torsional relaxation can be rationalized. The lower diffusion coefficient (see Table 1) suggests that all relaxation rates in MD should be lower by a factor of 1.4 because of the higher effective viscosity. Moreover, the asymptotic lt value underestimates the magnitude of torsional fluctuations in the tetradecamer by about 10% because of the vertical shift of the Dt(L) plots (see Figure 10 and Figure 2S). When these two factors are taken into account in BD simulations, the difference in τt essentially disappears, and the corresponding time autocorrelation functions almost overlap for t e τt (see dashed plots in Figure 12). In contrast, for bending, the difference is further increased as should have been expected. Discussion The long-term objective of the present research line is to understand the mechanisms of coupling between macroscopic

Kinetic and Thermodynamic DNA Elasticity and microscopic states of the double-helical DNA. A number of interesting issues in this field can potentially be clarified by using computer molecular modeling. As a step toward this goal, we and others are looking for practical methods of probing the mesoscopic polymer properties of DNA by atomistic MD simulations.28-30 Several problems in this approach have been defined and partially solved, gradually opening a possibility for more extensive investigations. We try to gain further insight in this direction by directly comparing all-atom MD with BD of coarse-grained models. The latter are parametrized to reproduce the experimental properties of polymer DNA and provide the best extrapolation to small chain lengths used in MD. This approach allowed us to assess the accuracy and validate some of the earlier MD results. In addition, two interesting new effects have been documented, namely, anomalously fast bending dynamics in MD and an unusual pattern of relaxation of bending in internal fragments of long chains. Our BD simulations employ an improved discrete WLC model with bending, stretching, and twisting degrees of freedom. The modifications consist of (i) using point particles for constructing coordinate frames and (ii) incorporating explicit coupling of bending and twisting displacements. These features facilitate calculations with complex local potentials, notably, including directional anisotropy of local bending. Further development is possible toward detailed parametrizations used for studying sequence-dependent DNA elasticity.65 The torsional degree of freedom does not affect the earlier conclusions concerning evaluation of bending and stretching PLs from timelimited ensembles of chain configurations.30 Qualitatively, torsional dynamics look similar to stretching. Twisting and stretching fluctuations in internal fragments of long chains are almost uncoupled as if these fragments were not united. In contrast, local bending fluctuations appear to be constrained by the overall bend. The last effect was first observed and discussed previously30 and is confirmed here. Additional studies are necessary to check whether it persists with chain lengths exceeding the bending PL. The relative rate of torsional dynamics depends on the chain length. In long DNA, it is intermediate between bending and stretching, but in very short fragments, the twisting fluctuations are perhaps faster than the others. To make the composite-bead WLC model of DNA closer to the all-atom representation, we made the link length equal to the interbase pair distance. For distinction, we call such chains fine-grained. Similar models are employed in the literature for studying sequence-dependent DNA elasticity,65 and they are sometimes referred to as rigid-base-pair chains.66 To the best of my knowledge, however, such discretization has not been applied in BD simulations. Validation tests show that, for short DNA, the rate of translational diffusion is close to experiment and that the fine- and coarse-grained representations appear to be mutually consistent without adjustment of the bead radius. This condition is fortuitous, but it does not contradict the physical intuition and probably results from the properties of the Rotne-Prager diffusion tensor. The fine-grained discretization allows examination of the possible effects of the well-known mechanical anisotropy of single base-pair steps. We found that the tilt/roll anisotropy does not cause noticeable perturbations of the WLC behavior. In contrast, the minor/major groove asymmetry of rolling causes strong specific perturbations previously documented for MD simulations.29 Regular oscillations of Db(L) in MD can have different origins; notably, they would be produced by any naive algorithm that measures DNA bend angles between vectors

J. Phys. Chem. B, Vol. 113, No. 7, 2009 2087 rigidly bound to nucleobases. However, a remarkable qualitative similarity between MD and BD, notably, the fact that coherent oscillations are also observed for the relaxation times, strongly suggests that the asymmetric rolling contributes to this behavior in MD. This observation is essential for the issue of the existence of an optimal helical axis in any instantaneous DNA conformation.31 Bending via asymmetric rolling might be inherent in double-helical DNA. Therefore, the optimal axis can legitimately look like the corresponding BD2 model, and construction of an axis without helical winding might be fundamentally incorrect. It is not evident that such an axis can be constructed for an instantaneous configuration of the BD2 model, and nor is it clear what the base-pair parameters will be in this case. The local flexibility introduced with the fine graining only weakly affects stretching dynamics, but it accelerates the relaxation of bending and twisting fluctuations. In addition, comparison of time autocorrelation functions revealed a flaw in the previously assumed correspondence of DNA lengths in the coarse-grained and atomistic models. These two factors result in appreciable downward re-evaluation of the duration of MD trajectories necessary for accurate assessment of DNA elasticity. We believe that these new estimates are more accurate and that they correctly account for the kinetic constraints imposed by the rheological conditions of hydrated DNA. It appears that 10-ns MD of a tetradecamer DNA already assures very high accuracy if the WLC parameters are evaluated from thermal fluctuations of 1 bps. This approach is well-founded;26,28,66 however, a closer view of MD data shows that the accuracy of such calculations is doubtful. First, these fluctuations are small, so the relative errors in the measured values are high. Second, both experiments and simulations suggest that fluctuations in neighboring steps commonly are correlated.67,68 The nonzero intercepts and nonlinearity of the D(L) plots in Figures 10 and 11 are formal signs of strong anticorrelations because the corresponding linear relations appear as a summation of variances of similarly distributed and supposedly independent random values (see eq 13). The correlations might be due, for instance, to long-range interactions between non-neighbor base pairs. They might also be artificial because the base-pair geometry is affected by many local atom movements that do not change the overall DNA length, bending, or twisting. For instance, a small displacement of one atom in a nucleobase slightly moves the average base-pair plane, which concertedly changes the measured distances from the two neighbor planes. The corresponding energy change might be independent of the instantaneous DNA length, so that this local motion is not coupled to the overall stretching. In this case, we would obtain spurious anticorrelations of consecutive interbase pair distances, and the fluctuations of the overall length would be smaller than the sum of the base-pair-step fluctuations. Finally, the evaluation of elastic parameters from single base-pair steps is incorrect if the D(L) dependences oscillate as in Figure 8. The foregoing problems are reduced if the WLC fluctuations are measured for DNA stretches of several base-pair steps. Our results demonstrate that one helical turn is a compromise length that allows radical improvement in the accuracy with tolerable increase in computational costs. BD simulations show that the perturbations caused by anisotropic local bending do not affect the WLC behavior of fragments corresponding to integral numbers of helical turns. This finding also appears true for MD simulations, and one turn is the minimum length that allows relatively accurate assessments. Longer fragments might certainly be necessary for complex base pair sequences with longrange correlations. Our data show that, for double helixes of

2088 J. Phys. Chem. B, Vol. 113, No. 7, 2009 11-15 base pairs, MD trajectories of 50-100 ns are sufficient, which is quite ordinary for the modern all-atom simulations. Remarkably, similar values for the bending PL are obtained using the programs Curves and 3DNA, suggesting that construction of an optimal helical axis might be unnecessary. Such axes are constructed by optimization algorithms that strongly depend on built-in weighting factors that might require ad hoc adjustment.29 This causes side effects that are difficult to trace, and getting rid of this difficulty would significantly simplify the calculations. For twisting and stretching, the convergence between the two methods is less evident, however, and this issue requires further examination. Our results refine and confirm the earlier finding concerning the very high absolute rate of bending fluctuations in MD. Here, we presented a thorough verification including parallel examination of the rates of twisting and translational diffusion. It appears that, for a tetradecamer double helix with identical rates of translational diffusion, the bending fluctuations in MD relax an order of magnitude faster than predicted by BD (see dashed traces in Figure 12). The analytical estimate of the absolute rate28 agrees with BD simulations. In contrast, good quantitative agreement is demonstrated for the rates of twisting fluctuations. As a result, whereas bending is significantly slower than twisting in the analytical and discrete WLC models, in the atomistic representation, it is much faster. The torsional dynamics in MD is relatively well described by the harmonic law assumed in the analytical and discrete WLC models; that is, both the amplitudes and the decay rates of fluctuations correspond to approximately the same harmonic stiffness. The bending dynamics is evidently anomalous because the fluctuations decay 10 times faster than expected for overdamped harmonic regime. The foregoing anomaly is formally similar to the discrepancy between the kinetic and thermodynamic DNA stiffnesses invoked in the earlier literature in relation to experiments that probed relaxation rates of short-time bending fluctuations.52,69,70 In the submicrosecond range accessible to these methods, the measured relaxation rates appear to be higher than predicted by harmonic models with a canonical bending PL of 50 nm. This is explained by using the concepts of static and dynamic contributions to the bending PL.71 The static contribution involves all motions that occur beyond the experimental time range; it is assumed that these are conformational transitions between strong “static” bends. This discrepancy is an experimental limitation, and it should disappear if the measured time range were extended to infinity. Our results are qualitatively different because we probe the apparent kinetic and thermodynamic bending stiffnesses in the same time range. The argument of the limited duration of MD trajectories cannot save this situation even in principle. As we showed, the duration of 100 ns is sufficient for sampling the corresponding bending fluctuations in BD. The existence in MD of undetectable slower motions cannot be rationalized in terms the current coarsegrained models. Moreover, to increase τb by a factor of 10, these hypothetical motions must have large amplitudes and reduce the bending PL accordingly, that is, to unphysically low values. The observed discrepancy is too large to be attributable to the unharmonic contributions of the MD force field. It can be explained only by a special mechanism of bending relaxation that decouples the kinetic and thermodynamic stiffnesses. As suggested previously,30 this mechanism might be due to the postulated intrinsic compression of the DNA backbone in the double helix.12,13,72 The compressed backbone hypothesis (CBH) predicts metastable local bends in a homopolymer double helix due to regular modulations of the minor and major DNA

Mazur grooves. The bending fluctuations occur as a result of 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 to be slightly strained, and this internal strain drives its trace to a new configuration. In such a system, the rate of relaxation might not correspond to the value of the bending PL. It might be faster, as in our simulations, which use a homogeneous base pair sequence where the groove profiles might have an infinite number of configurations with the same energy (strong frustration). The relaxation can also be much slower if the base pair sequence provides long-lived metastable states that can temporarily fix the phase of the groove modulations. This situation might be relevant to the time-resolved Stokes shift experiments that revealed very slow relaxation processes in short DNA that also cannot be explained by the coarse-grained WLC models.39,40 The foregoing issues will be clarified in future studies. An immediate practical conclusion following from our results concerns the possibility of using free MD simulation for probing the DNA elasticity. Reliable and accurate estimates of the elastic parameters can be obtained with reasonable efforts by properly choosing the simulation condition and the methods of analysis. This opens way to more extensive applications including the effects of the environmental factors as well as the weak external stress comparable with experimental and in vivo conditions. This work is now in progress. Supporting Information Available: Evaluation of three PLs from trajectories of increased durations with constant sampling intervals and comparison of bending and twisting dynamics of a tetradecamer DNA fragment according to atomistic and coarsegrained simulations. This information is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Landau, L. D.; Lifshitz, E. M. Statistical Physics, Part 1; Nauka: Moscow, U.S.S.R., 1976. (2) Cantor, C. R.; Schimmel, P. R. Biophysical Chemistry, Part III: The BehaVior of Biological Macromolecules; W. H. Freeman: San Francisco, 1980. (3) Garcia, H. G.; Grayson, P.; Han, L.; Inamdar, M.; Kondev, J.; Nelson, P. C.; Phillips, R.; Widom, J.; Wiggins, P. A. Biopolymers 2007, 85, 115–130. (4) Cloutier, T. E.; Widom, J. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 3645–3650. (5) Du, Q.; Smith, C.; Shiffeldrim, N.; Vologodskaia, M.; Vologodskii, A. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 5397–5402. (6) Wiggins, P. A.; Heijden, T. V. D.; Moreno-Herrero, F.; Spakowitz, A.; Phillips, R.; Widom, J.; Dekker, C.; Nelson, P. C. Nature Nanotechnol. 2006, 1, 137–141. (7) Wiggins, P. A.; Nelson, P. C. Phys. ReV. E 2006, 73, 031906. (8) Mazur, A. K. Phys. ReV. Lett. 2007, 98, 218102. (9) Yuan, C.; Chen, H.; Lou, X. W.; Archer, L. A. Phys. ReV. Lett. 2008, 100, 018102. (10) Travers, A.; Muskhelishvili, G. Nat. ReV. Microbiol. 2005, 3, 157– 169. (11) Fujimoto, B. S.; Brewood, G. P.; Schurr, J. M. Biophys. J. 2006, 91, 4166–4179. (12) Mazur, A. K.; Kamashev, D. E. Phys. ReV. E 2002, 66, 011917. (13) Kamashev, D. E.; Mazur, A. K. Biochemistry 2004, 43, 8160–8168. (14) 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. (15) MacKerell, A. D., Jr.; Wi′orkiewicz-Kuczera, J.; Karplus, M. J. Am. Chem. Soc. 1995, 117, 11946–11975. (16) Perez, A.; Marchan, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Biophys. J. 2007, 92, 3817–3829. (17) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089– 10092. (18) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (19) Cheatham, T. E., III; Kollman, P. A. Annu. ReV. Phys. Chem. 2000, 51, 435–471.

Kinetic and Thermodynamic DNA Elasticity (20) Perez, A.; Luque, F. J.; Orozco, M. J. Am. Chem. Soc. 2007, 129, 14739–14745. (21) Hagerman, P. J. Annu. ReV. Biophys. Biophys. Chem. 1988, 17, 265–286. (22) Bustamante, C.; Marko, J. F.; Siggia, E. D.; Smith, S. Science 1994, 265, 1599–1600. (23) Vologodskii, A. V. Macromolecules 1994, 27, 5623–5625. (24) Bouchiat, C.; Wang, M. D.; Allemand, J.; Strick, T.; Block, S. M.; Croquette, V. Biophys. J. 1999, 76, 409–413. (25) MacKerell, A. E., Jr.; Banavali, N.; Foloppe, N. Biopolymers 2001, 56, 257–265. (26) 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–11168. (27) Bruant, N.; Flatters, D.; Lavery, R.; Genest, D. Biophys. J. 1999, 77, 2366–2376. (28) Lankas, F.; Sponer, J.; Hobza, P.; Langowski, J. J. Mol. Biol. 2000, 299, 695–709. (29) Mazur, A. K. Biophys. J. 2006, 91, 4507–4518. (30) Mazur, A. K. J. Phys. Chem. B 2008, 112, 4975–4982. (31) Lavery, R.; Sklenar, H. J. Biomol. Struct. Dyn. 1988, 6, 63–91. (32) Lu, X.-J.; Olson, W. K. J. Mol. Biol. 1999, 285, 1563–1575. (33) Allison, S. A.; McCammon, J. A. Biopolymers 1984, 23, 363–375. (34) Allison, S. A. Macromolecules 1986, 19, 118–124. (35) Chirico, G.; Langowski, J. Macromolecules 1992, 25, 769–775. (36) Jian, H.; Vologodskii, A. V.; Schlick, T. J. Comput. Phys. 1997, 136, 168–179. (37) Vologodskii, A. Biophys. J. 2006, 90, 1594–1597. (38) Vologodskii, A. In Computational Studies of DNA and RNA; Lankas, F., Sponer, J., Ed.; Springer: Amsterdam, 2006; pp 579-604. (39) Brauns, E. B.; Madaras, M. L.; Coleman, R. S.; Murphy, C. J.; Berg, M. A. Phys. ReV. Lett. 2002, 88, 158101. (40) Berg, M. A.; Colemanb, R. S.; Murphy, C. J. Phys. Chem. Chem. Phys. 2000, 10, 1229–1242. (41) Halle, B. Philos. Trans. R. Soc. London B: Biol. Sci. 2004, 359, 1207–1224. (42) Nilsson, L.; Halle, B. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 13867–13872. (43) Furse, K. E.; Lindquist, B. A.; Corcelli, S. A. J. Phys. Chem. B. 2008, 112, 3231–3239. (44) Zhurkin, V. B.; Lysov, Y. P.; Ivanov, V. I. Nucleic Acids Res. 1979, 6, 1081–1096. (45) Mazur, A. K. J. Chem. Phys. 1999, 111, 1407–1414. (46) Allison, S.; Austin, R.; Hogan, M. J. Chem. Phys. 1989, 90, 3843– 3854.

J. Phys. Chem. B, Vol. 113, No. 7, 2009 2089 (47) Chirico, G.; Langowski, J. Biopolymers 1994, 34, 415–433. (48) Jian, H.; Schlick, T.; Vologodskii, A. J. Mol. Biol. 1998, 284, 287– 296. (49) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352– 1360. (50) Iniesta, A.; de la Torre, J. G. J. Chem. Phys. 1990, 92, 2015–2018. (51) Rotne, J.; Prager, S. J. Chem. Phys. 1969, 50, 4831–4837. (52) Allison, S. A.; Schurr, J. M. Macromolecules 1997, 30, 7131–7142. (53) Cheatham, T. E., III; Cieplak, P.; Kollman, P. A. J. Biomol. Struct. Dyn. 1999, 16, 845–862. (54) Jorgensen, W. L.; Chandreskhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (55) Mazur, A. K. J. Comput. Chem. 1997, 18, 1354–1364. (56) Mazur, A. K. In Computational Biochemistry and Biophysics; Becker, O. M., MacKerell, A. D., Jr., Roux, B., Watanabe, M., Ed.; Marcel Dekker: New York, 2001; pp 115-131. (57) Mazur, A. K. J. Am. Chem. Soc. 1998, 120, 10928–10937. (58) Mazur, A. K. J. Phys. Chem. B 1998, 102, 473–479. (59) Lu, X.-J.; Olson, W. K. Nucleiic Acids Res. 2003, 31, 5108–5121. (60) Stellwagen, E.; Lu, Y.; Stellwagen, N. C. Biochemistry 2003, 42, 11745–11750. (61) Hagerman, P. J.; Zimm, B. H. Biopolymers 1981, 20, 1481–1502. (62) Yeh, I.-C.; Hummer, G. J. Am. Chem. Soc. 2002, 124, 6563–6568. (63) Yeh, I.-C.; Hummer, G. Biophys. J. 2004, 86, 681–689. (64) Bryant, Z.; Stone, M. D.; Gore, J.; Smith, S. B.; Cozzarelli, N. R.; Bustamante, C. Nature 2003, 424, 338–341. (65) Swigon, D.; Coleman, B. D.; Olson, W. K. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 9879–9884. (66) Becker, N. B.; Everaers, R. Phys. ReV. E 2007, 76, 021923. (67) Yanagi, K.; Prive′, G. G.; Dickerson, R. E. J. Mol. Biol. 1991, 217, 201–214. (68) Lankas, F.; Sponer, J.; Langowski, J.; Cheatham, T. E., III Biophys. J. 2003, 85, 2872–2883. (69) Okonogi, T. M.; Reese, A. W.; Alley, S. C.; Hopkins, P. B.; Robinson, R. H. Biophys. J. 1999, 77, 3256–3276. (70) Naimushin, A. N.; Fujimoto, B. S.; Schurr, J. M. Biophys. J. 2000, 78, 1498–1518. (71) 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. (72) Mazur, A. K. J. Am. Chem. Soc. 2000, 122, 12778–12785.

JP8098945

Kinetic and Thermodynamic DNA Elasticity at Micro- and Mesoscopic

Recommend Documents