The conformational flexibility of DNA is central to its many biological functions including recognition by proteins during gene regulation, DNA repair, packaging in the cell, and transient melting during transcription and replication (1–6). Specific binding by proteins is not only determined by specific interactions between DNA and proteins but also by the sequence-dependent structure and deformability of the DNA helix (1–4,7–16). In vivo DNA is mostly found in circularly closed form including supercoiling that causes a mechanical stress (torque) on the DNA twist and can strongly affect recognition by proteins and influence transcription and replication fidelity (6,15). Therefore, better understanding of the sequence-dependence of the DNA twist flexibility is of particular interest. The total twist and to some degree also the twist elasticity of DNA can be measured using DNA ring-closure ligation (cyclization) experiments with DNA molecules of different length and including a helix phasing analysis (17–19). This technique has been used to study the twist flexibility of DNA fragments of various length and sequence in solution. Alternatively, time-resolved fluorescence polarization anisotropy (on DNA with intercalated ethidium-bromide) can also be used to measure torsional DNA flexibility near the equilibrium state (20,21). These studies indicate that the overall base composition has only a modest influence on the average twist flexibility of DNA. Recently, it has also become possible to study the elastic response of single DNA molecules upon application of external torques to un- or over-twisted DNA (22). Such experiments provide valuable insights into the physical and elastic properties of DNA molecules. However, only relatively long DNA molecules (>1 kbp) can be studied that can react on the external stress by local as well as global conformational changes (e.g., supercoiling that relaxes part of the twist-stress). Similar to the DNA ring closure experiments, an overall twist elasticity can be obtained—but with only modest insight provided into the conformational changes taking place at the molecular level.

High-resolution experimental structures of isolated DNA molecules and complexes with organic ligands and proteins allow us to study the fine structure of DNA. Under the assumption that observed structural variations reflect the intrinsic deformability of B-DNA the analysis provides insights into the flexibility of DNA near the native state at atomic resolution including also DNA twist elasticity (3,4,7–15). Olson and co-workers (4) have extensively analyzed structural variations in available DNA oligonucleotide structures and complexes with proteins and were able to derive a sequence-dependent empirical energy function for the DNA helical elasticity. The approach also allows to investigate possible correlations between helical parameters that describe the flexibility of nucleic acids. However, it is not clear how well a set of crystal structures reflects the structural flexibility of DNA in solution.

Alternatively, the progress in molecular mechanics force fields and simulation methodology has made it possible to investigate DNA helical flexibility using molecular dynamics (MD) simulations including surrounding water molecules and ions (reviewed in (23–26)). Such MD simulations of DNA result in stable structures close to the experimental DNA conformations on the nanosecond timescale and can be used to characterize the equilibrium fluctuations of helical parameters (27–31). Computational large-scale studies on many different DNA molecules have been used to study the sequence-dependence of DNA flexibility and are also helpful to improve the molecular mechanics force fields (32,33). However, unrestrained MD simulations may allow only a limited sampling of possible substates due to energy barriers. To better understand the molecular mechanism of elastic twist deformations in DNA, in the current study, MD simulations combined with a twistlike restraining potential and the umbrella sampling method were used. Simulations were performed on several 12-bp DNA oligo-nucleotides with different sequences and for a range of total twist angles of the central 9-bp steps. The simulations allowed us to calculate the total free energy change, to characterize the twist change of individual basepair steps and to analyze the change of other DNA helical parameters in response to external twist deformations. In addition, the onset of an untwisting transition starting from one end of a DNA helix was observed.

Simulations were performed on six 12-bp B-DNA molecules (Table 1) with the general sequence (5′-CGCGNNNNCGCG)₂ and six different central sequences (AATT, TATA, CGCG, GATC, CATG, CTAG). The central sequences cover each of the 10 possible basepair steps in DNA (AA, GG, GC, AG, TA, CG, CA, GC, AT, GT) at least once. Standard B-DNA start structures were generated using the Nucgen program of the AMBER8 package (34). Each system was neutralized by adding 22 K⁺ counter ions and solvated with ~4500 TIP3P water molecules (35) in a rectangular box using the xleap module of AMBER8. The simulation systems were subjected to energy minimization (1000 steps) using the Sander module. The PARM99 force field (36,37) was used for all simulations. During MD each DNA was initially harmonically restrained (25 kcal mol⁻¹ Å⁻²) to the energy-minimized start coordinates and the system was heated up to 300 K in steps of 100 K followed by gradual removal of the positional restraints and 1-ns unrestraint equilibration of each system at 300 K. During MD the long-range electrostatic interactions were treated with the particle-mesh Ewald method (38) using a real-space cutoff distance of r_cuttoff = 9 Å. The RATTLE algorithm (39) was used to constrain bond vibrations involving hydrogen atoms, which allowed a time step of 2 fs. To introduce the external twist stress, a modified version of the tornrg subroutine (which calculates the torsion-angle restraint energy in Sander) was used to harmonically restrain a twist angle between the second and eleventh basepair of each 12-bp duplex DNA. This twist angle (in the following termed τ) is formed by the distance vector between C₁′ atoms of the 11th basepair projected onto the plane defined by the distance vector of C₁′ atoms of the second basepair and the axis connecting the midpoints of the two pairs of C₁′ atoms (see Fig. 1). It is important to note that the above-defined twist-angle is only an approximation to the “real” twist because it is calculated with respect to a linear helical axis along the DNA. To compare the approximate twist restraining coordinate with a more accurate definition that also accounts for a possible curved helical axis, it was compared to the twist obtained using the program Curves (40,41). For the range of twist angles considered in this study, a good correspondence was observed (Fig. 2). During restraint MD simulations a quadratic penalty function for deviations of the 9-bp twist with respect to a reference twist value was used. To obtain start structures for the umbrella sampling simulations the reference τ-angles was changed in steps of 5° with an equilibration run of 0.2 ns for each τ. A total τ-range of ~225–360° per 9-bp steps was covered representing a range of ~25–40° per basepair step. The twist along the central 9-bp steps after the unrestrained equilibrium simulations was ~292° (corresponding to ~32.5° per bp step and depending slightly on the sequence). The total τ-range of ~225–360° per 9-bp steps represents a twist deformation of ~±70° per 9 bp steps (or ±7.5° per bp step) with respect to the equilibrium twist. During the equilibration the central basepair was weakly positionally restrained (force constant 0.01 kcal mol⁻¹ Å⁻²). This does not affect the internal dynamics of the basepair to any significant extent but avoids any overall rotation of the DNA and keeps the helical axis approximately aligned with the long axis of the rectangular box. Subsequently, however, during 1-ns data gathering for each reference τ-angle, the DNA was completely free to move except for the twist restraint. The actual restraining twist angle was recorded every 10 simulation steps. Complete coordinate frames were stored every 2 ps. The potential of mean force (PMF) along the reaction coordinate was calculated using the weighted histogram analysis (WHAM) method (42,43). Helical and dihedral angle parameters of recorded structures were calculated using Curves (40,41).

TABLE 1 Optimal twist and fourth-order polynomial fit of the calculated PMF for twist deformations

FIGURE 1 Definition of the twist-restraining coordinate. The twist-restraining angle (termed τ) is formed by the distance vector (upper double arrow) between C1′ atoms (blue spheres in the left panel) of the 11th basepair projected onto the plane (more ...)

FIGURE 2 Correlation of twist (per basepair step) calculated using the program Curves (40,41) and the twist as defined in Fig. 1 (given as average per-bp step twist instead of per 9-bp step as indicated in Fig. 1) that has been used as twist-restraining coordinate (more ...)

The range of twist values covered during the simulations was limited to an average 25–40° per bp step (Fig. 3). Further over-twisting of the DNA beyond these limits led in all cases to a steep rise of the calculated PMF due to sterical strain (see Fig. 3). Untwisting below ~25° per bp resulted in the onset of significant structural changes (stable transitions) of the DNA discussed in the last paragraph of Results and Discussion. These structural changes occurred spontaneously and asymmetrically starting at one end of the DNA and require much longer simulations to achieve convergence of the associated free-energy changes. The limitation to an average range of ~25–40° per bp step allowed a very good convergence of the free energy curves (compare dotted, dashed, and bold lines in Fig. 3 that correspond to different simulation times for each twist reference value). For all sequences and most parts of the plots, the final drift of the free energy curves was smaller than 0.5 kcal mol⁻¹. The basic shape of the free energy curves was similar for all sequences (Fig. 3) with a minimum in the free energy curves in the range of ~31–33°, which is ~3–5° smaller than the experimentally observed 36° per bp step. Such smaller average twist angles compared to standard B-DNA were also observed in unrestrained MD simulations using the same force field (31.5–33° per bp step, data not shown). It agrees also with previous MD-simulation studies on B-DNA using the AMBER force field, indicating that both the PARM94 (44) and, to a lesser degree, the current PARM99 force-field parameters (36,37), result in slightly under-twisted average structures during MD simulations (32,33,36)—which also implies that these force fields do not reproduce the sequence-dependent equilibrium twist of B-DNA.

FIGURE 3 Potential of mean force (PMF) of twist deformations obtained during umbrella-sampling molecular dynamics simulations (central sequences are indicated in each panel). The calculated PMF (using the WHAM approach (42)) is plotted versus the applied twist-restraining (more ...)

The optimal twist differs slightly between the six sequences (Table 1). As expected, the calculated free energy curves are not symmetric with respect to the twist minimum. Over-twisting of the DNA causes generally a stronger increase of the free energy compared to untwisting of the DNAs. The calculated PMFs can be approximated very well by a fourth-order polynomial (Table 1). The quadratic term which dominates for small twist deformations indicates a range of twist elastic constants between ~0.03 and ~0.07 kcal mol⁻¹ deg⁻². The calculated range of twist elastic constants is in very good agreement with available experimental data (18–21) and results of unrestrained DNA simulations that obtain effective force constants for twist deformations from twist distribution functions (29–31). It translates to a range of twist fluctuations per bp step (= sqrt (R T/ c), where R is the gas constant, T is the temperature, and c is the elastic constant) at room temperature of 3–5°. In terms of overall twist flexibility near the optimal twist, the central TATA sequence was the most flexible, followed by CATG, CTAG, and GATC (Table 1). The stiffest sequences, with respect to twist, are the DNAs with central AATT and central GCGC sequence. The overall modest variation of the effective twist flexibility among the various sequences agrees with the experimental observation based on DNA cyclization experiments (17–19) and time-resolved fluorescence polarization anisotropy (20,21).

FIGURE 4 Average twist per basepair step versus applied twist stress. Each bar indicates the average twist per basepair (for the basepairs indicated on the x axis) along the central five basepair steps of three example DNA oligo-nucleotides (central sequence given (more ...)

FIGURE 5 Local twist of each of the 10 possible basepair steps versus applied twist (per basepair step along each DNA). Each panel represents the average twist (per indicated basepair step) from all simulations (at the indicated applied twist) on the six different (more ...)

TABLE 2 Twist flexibility of basepair steps and coupling to slide and roll

FIGURE 6 Correlation of helical parameters roll (A) and slide (B) with the twist angle obtained using Curves (40,41) at each of the 10 possible basepair steps (data obtained from twist restraining simulations on all six oligonucleotides, see Table 1).

As expected from the observed coupling of twist changes to other helical parameters—in particular, roll—the induced twist deformation also lead to an overall slight average bending of the DNAs (Fig. 7 a) ranging from ~10° (with an average restraining twist of 40°) to ~40° in case of a restraining twist of ~25° (Fig. 7 a). Note that even in unrestrained simulations, a B-DNA molecule for entropic reasons is, on average, bent (even if the straight form is the energetically most favorable form). DNA untwisting also resulted in an increase of the minor groove size (Fig. 7 b). As an example, the observed structural changes are illustrated for the TATA case in Fig. 8. The over-twisted structure has a narrow minor groove with a relatively straight helical axis. The untwisted structure has a wider minor groove and a narrower major groove compared to standard B-DNA more typical for A-form structures (Fig. 8). As expected, overall the untwisting increases the A-form character of the DNAs. Another indicator of the increased A-form character upon untwisting of the DNA is the observed shift of the distribution of the deoxy-ribose δ-dihedral angles (of the central parts of the DNA) from values characteristic for B-DNA (at high over-twisting stress) to smaller values with contributions characteristic for A-form (δ-dihedral angle at ~85°) upon untwisting (Fig. 9). Interestingly, the γ-dihedral angle distribution stayed close to the distribution characteristic for B-DNA (+gauche, Fig. 9). It has been reported that during very long MD simulations of B-DNA using the AMBER force-field transitions to backbone structure with α/γ, flips (crank shift transitions) may occur (47). This was not observed in the present simulations, possibly because the (independent) simulations per each twist-restraining window were, overall, relatively short (<2 ns). However, as indicated in Fig. 9 (lower two panels), a number of B_II states (a correlated shift of nucleic acid backbone - and ζ-dihedral angles from t to −g and −g to t, respectively) were observed (~10% B_II). This is slightly larger than what has been found in unrestrained simulations (~6% B_II, (33)) and likely due to the fact that B_I-B_II transitions can affect the twist of DNA. However, no clear correlation between B_II-state frequency and twist stress in the present twist range was observed (Fig. 9).

FIGURE 7 Correlation of bending (A) and minor groove width (B) versus applied restraining twist (per basepair step along each DNA). Global bend angle and minor groove width were calculated using Curves (40,41).

FIGURE 8 Comparison of average DNA structures (central TATA case) during restrained MD simulations at a twist restraint corresponding to 40° (A), 33° (B), and 25° (C) per bp step, respectively. The view in the case of the upper set of panels (more ...)

FIGURE 9 Nucleic acid backbone dihedral angle distributions (TATA case) for the central-six basepair steps during the data-gathering simulation period. Continuous line indicates distribution at high twist stress, 37–40° (per basepair step); the (more ...)

FIGURE 10 (A) Average structure of a DNA oligonucleotide (central TATA sequence) during a twist-restraining simulation (1 ns data gathering time) with a twist-restraint reference of 20° per bp step. The two views in panel A differ by a 90° rotation (more ...)

Understanding DNA deformability in particular of the DNA-helical twist is of fundamental importance to understand processes like DNA recognition, packing, transcription, and replication. In this study, a molecular dynamics umbrella sampling approach was used to study systematically, on several different DNA molecules, the free energy change associated with an over- or untwisting of a segment of DNA. It was possible to extract the shape of the free energy curves versus applied twist restraint, and qualitatively similar curves for the six tested DNA molecules were observed. The approach goes beyond previous unrestrained simulations that do not allow us to directly extract the free energy change associated with DNA deformations. It also complements experimental single molecule studies (22) mostly performed on long DNAs that largely relax applied torsional stress through supercoiling and approaches that evaluate DNA flexibility indirectly from analyzing experimental DNA crystal structure (4–12). The current results are in very good agreement with available experimental data on the overall twist flexibility of DNA as well as on the local per-bp step flexibility and the coupling to other helical parameters. This result indicates that the structural variation seen in DNA crystal structures does indeed reflect to a significant degree the DNA flexibility in solution. The simulation study indicates that the strongly nonuniform local response to external twist stress can significantly alter the local arrangement of chemical groups recognized by ligands. In addition, twist stress can alter the overall shape of a DNA segment (Fig. 8), including DNA bending and groove deformation again strongly modulating the accessibility by proteins and other ligands. Interestingly, the onset of a structural transition of the double helix to a structure with a strongly reduced twist compared to B-DNA and opened minor and major grooves was observed during extensive DNA untwisting simulations. The significance of this transition will be subject of a future study.

This study was performed using the computational resources of the Computational Laboratories for Analysis, Modeling and Visualization at International University Bremen and supercomputer resources of the Environmental Molecular Science Laboratories at the Pacific Northwest National Laboratories, grant No. GC11-2002. S.K. is supported by the BIOlogical RECognition graduate program at International University Bremen and by a grant from the VolkswagenStiftung to M.Z.