|Home | About | Journals | Submit | Contact Us | Français|
The microscopic mechanics of DNA stretching was characterized using extensive molecular dynamics simulations. By employing an anisotropic pressure control method, realistic force-extension dependences of effectively infinite DNA molecules were obtained. A coexistence of B- and S-DNA domains was observed during the overstretching transition. The simulations revealed that strain softening may occur in the process of stretching torsionally constrained DNA. The latter observation was qualitatively reconciled with available experimental data using a random-field Ising model.
The mechanical properties of double-stranded DNA (dsDNA) have been the subject of extensive studies [1–7] because of their fundamental importance to gene regulation processes in biological cells. Single molecule manipulation experiments revealed a characteristic plateau in the force-extension curve of dsDNA that signifies a highly cooperative transition from a canonical B-DNA structure to an overstretched S-DNA conformation, the so-called B to S transition . Despite extensive experimental and theoretical studies, the nature of the B to S transition remains controversial . Molecular dynamics (MD) simulations revealed that stretching DNA with a force transforms the canonical B-DNA structure into a ladder-like conformation [2, 6]. However, subsequent thermodynamic analysis of DNA stretching suggested that the “B–S” plateau indicates a melting transition, i.e., separation of the two DNA strands that occurs at the begining of (and throughout) the plateau . Furthermore, both theory  and experiment  suggest that depending on the twist of the DNA helix, i.e., the number of base pairs per turn of the helix, several DNA conformations may coexist during the transition.
In contrast to previous MD studies of short ds-DNA fragments [2, 6, 9–11], in this letter we report the force-extension dependence of an effectively infinite DNA molecule that was stretched using an anisotropic pressure-control method. Figure 1(a) illustrates the setup of our simulations. A fragment of dsDNA two helical turns in length was submerged in an aqueous solution of 1 M KCl. Each strand of the DNA helix was covalently linked to itself over the periodic boundary of the system. Thus, any deformation modes with wave lengths bigger than the system size were suppressed. We constructed three systems having the following nucleotide sequences: poly(dA)·poly(dT), poly(d(CGATATATCG))·poly(d(GCTATATAGC)), and poly(dC)·poly(dG), which are referred to as sequence A, B and C, respectively. Note that sequence A and C is homogeneous, and sequence B is heterogeneous. According to experiment, the systems were built to contain 20 base pairs (bps) in two helical turns for sequences A and B [12, 13], and 21 bps for sequence C .
Following assembly, each system was equilibrated in the NPT ensemble  (P = 1 bar, T=310 K) for 10 ns using the program NAMD , the parm94 force field for DNA , the TIP3P model of water , standard parameters for ions , particle-mesh Ewald (PME) electrostatics and multiple time stepping . Van der Waals interactions were calculated using a smooth (10-12 Å) cutoff. The temperature was kept constant by applying a Langevin thermostat  with damping rate 1 ps−1 to all non-hydrogen atoms. After equilibration, the systems measured about 78 × 78 × L0 Å3, where the average DNA contour length L0 was 65.8, 65.3 and 67.0 Å for sequence A, B, and C, respectively. During the equilibration, DNA remained in the B-form conformation.
To obtain the force-extension curve, each system was simulated under anisotropic pressure conditions that were maintained using the Nosé-Hoover Langevin piston pressure control . While the Pxx and Pyy components of the pressure tensor were kept at 1 bar, the Pzz component was stepwise decreased [Fig. 1(a)]. When Pzz was assigned a negative value, the system stretched in the z direction until the total internal stress in the simulated system balanced the applied pressure Pzz. The water density remained constant. Repeating such simulations at different (negative) values of Pzz yielded the stress-strain curves shown in Fig. 1(b). Each point in the curves was obtained by equilibrating the system at a preset value of Pzz, requiring from 5 to 40 ns. Such equilibration times are considerably larger than a typical relaxation time of the internal stress, 2Lz/υ ~ 10 ps, where υ (~ 15 Å/ps) is the speed of sound in water.
Each stress-extention curve shown in Fig. 1(b) has three characteristic regions: elastic, transition, and overstretched. The DNA was observed to deform elastically up to a normalized extension (Lz/L0) of 1.06 for sequence A and B and up to 1.12 for sequence C. A plateau indicating a transition from canonical to overstretched structure was observed in a very narrow range of negative pressures [Fig. 1(a)]. The overstretched region begins at the end of the transition plateau where DNA extends to 1.6L0, in excellent agreement with experiment [2, 21]. Small plateaus in the overstretched region are reminiscent of the experimentally observed melting transition that follows the B to S transition [22, 23].
The kinetics of the B to S transition is illustrated in Fig. 2. At Pzz = −24 bar, the DNA molecule is elastically stretched; its average length is 67.7 Å [Fig. 2(a)]. When Pzz was set to −32 bar, the DNA's length was about 69.1 Å for the most of the 20 ns trajectory but could transiently increase to 78 Å as a group of 7-10 bps transiently lost its canonical B-DNA structure [Fig. 2(b)]. These elastic instabilities indicated proximity of the transition. When temperature was raised to 350 K, the transition occurred within 5 ns [Fig. 2(b)]. Further decreasing Pzz to −34 bar (at 310 K) initiated transition to the overstretched conformation [Fig. 2(c)]. After about 40 ns, the length of the DNA molecule reached 100 Å, i.e.,1.54L0. During the transition, the system visited several metastable states, which were identified as steps (Lz = 74, 80, 88 and 94 Å) in the DNA extension curve [Fig. 2(c)]. When Pzz was decreased directly from −32 to −36 bar, the transition occurred much faster than when quenched to −34 bar [Fig. 2(d)]; the system visited only two metastabe states (Lz = 85 and 93 Å). Decreasing Pzz from −36 to −44 bar increased the DNA length by 3 Å only [Fig. 2(e)], indicating the end of the transition plateau.
In our simulations, the applied external pressure Pzz is balanced by the internal stress in DNA and water. Assuming the pressure in bulk water is isotropic, the zz component of the water stress tensor bar. Hence, the tensile force inside DNA can be computed as , where S and SDNA are the area of the simulation system and of DNA, respectively, in the x–y plane. Because pN, we use the following approximate formula to compute the tensile force in DNA: F = −(Pzz − 1) × S.
Using the above formula, the force-extension dependence of a DNA molecule, Fig. 3(a), was computed from the stress-extension dependence of the entire system [Fig. 1(b)]. Estimated from the slope of the force-extension curve, the elastic modulus of sequence C is much smaller than that of sequence A and B, which is consistent with previous MD studies . A striking feature of the force-extension curves is a monotonic decrease of the tensile force in the transition region after initial increase to the yield force fy, i.e., a strain softening effect. This effect was also observed using an 11 pb per turn model of sequence A, and in the simulations performed using the latest parm-bsc0 force field  (data not shown). Microscopically, the strain softening effect is caused by the break up of the base pairing and base stacking interactions. Thus, depending on the history of the stretching process, the same stretching force can drive dsDNA into different conformations.
Figure 4(a-c) illustrates three possible conformations of DNA under the tensile force of 170 pN. The local deformation of the DNA structure was quantified by plotting the averaged over six bps distance between the phosphorous atoms of the DNA backbone versus the bp number [Fig 4(d)].
The conformation of the elastically stretched ds-DNA [Fig 4(a)] is close to the canonical B-DNA form, in which the tensile force is partially borne by orderly stacked bases. Beyond the elastic region, as the hydrogen bonds and base stacking yield, the DNA backbone stretches to bear a higher fraction of the tensile force. In the metastable state [Fig 4(b)], domains of B-DNA and overstretched S-DNA coexist, as also indicated by a non-uniform deformation of the DNA backbone [Fig 4(d)]. An independent 10-ns MD simulation restarted from a metastable state revealed a stable coexistence of B- and S-DNA domains when the DNA length Lz was kept constant. The tensile force in that simulation was constant and lower than the yield force. The coexistence of the Band S-DNA states, observed for the first time in our simulations, was postulated in previous theoretical models of DNA deformation [2, 26]. In the overstretched conformation [Fig. 4(c)], most of the bps are broken. Large fluctuations of in this regime reflect a variation of the local DNA structure. The small value of at n = 5 is due to the locally overwound DNA (the angle between consecutive bases in the same strand is ~150°), while the large value of at n = 14 corresponds to a “bubble” (the angle is ≤20°) of broken bps [Fig. 4(c)]. The convergence of the three force-extension curves in the overstretched region [Fig. 3(a)] indicates that the tensile force in this region is likely borne by DNA backbone only.
It has been established that a force-extension curve of a long DNA fragment such as λ-DNA exhibits an overstretching plateau at a tension of about 70 pN and that the tension slightly increases in the plateau region [2, 21]. In our simulations, the DNA molecules were torsionally restrained via the periodic boundary condition. Under similar conditions, experiment has shown a significantly higher yield force of about 110 pN , which is close to the yield force of sequence C [Fig. 3(a)]. Allowing DNA to unwind in our simulations was observed to greatly reduce the yield force. Figure 3(b) shows the force-extension curves of DNA molecules that have a nick in one of their strands. As such nicked DNA unwinds, its tertiary structure continuously transforms from the B-DNA to ladder-like structure. Such a continuous transformation of the tertiary structure does not require the hydrogen bonds between the two strands to break; the strain softening effect is abolished, which is consistent with experiments on stretching DNA of homogeneous sequences without applying torsional restraints .
Even under torsional constraints, experiments using λ-DNA revealed a force-extension plateau , i.e., no strain softening effect. We believe that this apparent discrepancy is due to the highly heterogeneous sequence of λ-DNA. Our simulations have demonstrated a strong dependence of the yield force on the DNA sequence. Hence, the B to S transition in dsDNA of a heterogeneous sequence would start with the DNA fragments characterized by the lowest yield force. The initial metastable growth of the S-DNA domains that could be accompanied by strain-softening would be limited to the domains of DNA sequence having similar yield forces. Driving the domain boundary towards B-DNA domains would require a higher stretching force that could also nucleate S-DNA domains elsewhere. Overall, the force-extension dependence appears as the “B–S” plateau.
To provide a theoretical model for the above description, we treat the microscopic process of the B to S transition as the motion of a magnetic domain driven by an external magnetic field in the Ising model . As the yield force was found to be sequence dependent, we considered a mean field model similar to the random field Ising model . The hamiltonian H = −J Σi(lili−1 + lili+1) + Σi(f − fi)li with li = 0.5(LS + LB) − Li, where the length of a DNA fragment (e.g. one helical turn) Li equals LB or LS for B- or S-DNA conformation, respectively. The coupling term −J Σi(lili−1 + lili+1) describes the interaction between neighboring DNA fragments. The condition of J > 0 is consistent with our MD results: at the same tensile force, the DNA conformation is stable when the two neighboring DNA fragments are in the B-DNA state [Fig. 4a] but unstable when one of them is stretched [Fig. 4b]. The yield force of the i-th fragment fi is analogous to the random local field in the Ising model. We define fi = f0 + f′ρi, where ρi is a random number between -1 and 1. The hamiltonian can be rewritten as H = − Σi liFi, where the effective local field Fi = J(li−1 + li+1) − f + fi. In the case of a homogeneous DNA sequence, ρi (= ρ) is a constant and the yield force of the whole DNA fy = J(LS − LB) + f0 + f′ρ. Here, fy is the force to nucleate an S-DNA domain. The required force to drive the domain walls, however, can be as low as f0 + f′ρ, consistent with the strain softening effect. For a heterogeneous DNA sequence, this model predicts that the “hysteresis” (or strain softening) effect is reduced or even eliminated when the sequence is random. This phenomenon was investigated in detail for the hysteresis loop of magnetization .
To test our theory, we performed simulations of torsionally constrained DNA having the following heterogeneous sequence: poly(dA10dC11)·poly(dT10dG11). In the middle of the overstretching plateau [Fig. 3(b)], we observed a stable (over 35 ns) state in which only the poly(dC11)·poly(dG11) fragment was overstretched. A force higher than the yield force of the dC11 fragment was required to induce stretching of the dA10 fragment. Hence, strain softening was not observed. Note that sequence B is also heterogeneous. However, the CGCG segment of sequence B could not maintain the B-DNA structure after neighboring ATATAT segments became overstretched. Hence, strain softening was observed.
In conclusion, we investigated the B to S transition in stretched dsDNA using an anisotropic pressure control method. Our simulations suggest that the force-extension dependence of a torsionally constrained DNA molecule can exhibit a strain softening effect if no stable coexistence of B-DNA and S-DNA domains is possible at a tensile force higher than the yield force. We expect such behavior to emerge when stretching torsionally constrained simple repeat DNA sequences that are abundant in mammalian genomes . The strain softening effect can be exploited in DNA nanotechnology to construct a mechanical “all or nothing” switch that unravels a torsionally constrained DNA helix when the tensile force exceeds a certain value. Physically, the strain softening effect is similar to the hysteresis effect in the first-order phase transition. Removing torsional restrains eliminates the strain softening effect, as in that case the DNA structure can continuously transform from B-DNA to S-DNA state through unwinding. Although the strain softening effect has not been experimentally observed in DNA, we believe such observations are possible using current experimental methods [8, 22].
The authors acknowledge useful discussions with Nigel Goldenfeld. Supercomputer time was provided via Large Resources Allocation Committee grant MCA05S028. This work was supported by grants from the National Institutes of Health (R01-HG003713 & PHS 5 P41-RR05969).