INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Sequence-dependent local conformational variability of DNA was first evidenced at the atomic resolution more than two decades ago (Wing et al., 1980). The structural basis of local variations, however, remains a matter of discussion. It is widely assumed that DNA local conformational variability and deformability play a considerable role in many biochemical and biological processes (Dickerson, 1999).

One of the most striking illustrations of sequence effects is the macroscopic curvature of DNA (Crothers et al., 1990; Marini et al., 1982; Diekmann, 1986; Olson and Zhurkin, 1996; Wu and Crothers, 1984). Phased runs of four to six consecutive adenines (a structure often termed as A-tract) induce a significant curvature. Despite intense research, the atomic-resolution picture of the A-tract-induced curvature remains elusive, and several models with distinct stereochemistry have been proposed (De Santis et al., 1990; Goodsell and Dickerson, 1994; Koo et al., 1986; Maroun and Olson, 1988; Olson et al., 1993). These include straight and intrinsically curved models of A-tract structures and suggest an important role played by the junctions between the A-tract and non-A-tract regions.

Fiber diffraction studies suggest that poly(dA)·poly(dT) adopts a distinct geometry that is characterized by a large propeller twist (around 26°), very narrow minor groove (~3.4 Å), and negative inclination of the basepairs with respect to the helical axis (Alexeev et al., 1987; Arnott et al., 1983; Coll et al., 1987). The high-propeller twist structure of poly(dA)·poly(dT) differs from the fiber B-DNA form adopted by a common sequence. Interestingly, the same architecture as for poly(dA)·poly(dT) was suggested for poly(dI)·poly(dC) and poly(dA-dI)·poly(dT-dC) (Leslie et al., 1980).

The fiber structural characteristics of poly(dA)·poly(dT) are essentially supported by oligonucleotide crystal structures of tracts with five or six consecutive adenines, showing a propeller twist around 20° and a minor groove width of 3.5-4.0 Å (Crothers and Shakked, 1999; DiGabriele et al., 1989; DiGabriele and Steitz, 1993; Nelson et al., 1987). Crystallographic studies bolster a view that A-tracts are straight with a clearly defined spine of hydration in their narrow minor groove (DiGabriele et al., 1989; DiGabriele and Steitz, 1993; Nelson et al., 1987). Hydroxyl radical footprinting experiments (Burkhoff and Tullius, 1987) show a progressive narrowing of the minor groove in the 5' to 3' direction, an observation which appears to be supported by NMR studies (Hud et al., 1999; Chuprina et al., 1991; Katahira et al., 1990). It has also been suggested, based on NMR experiments in the presence of NH₄⁺, that the minor groove of A-tracts is partly occupied by monovalent cations (Hud et al., 1999). A very recent NMR structure of d(GGCA₆CGG)·d(CCGT₆GCC) determined with residual dipolar coupling suggests a slightly reduced propeller twist of ~16° (MacDonald et al., 2001). Interestingly, A-tract-induced curvature and likely also the high-propeller A-tract structure persist only in the low-temperature region and are abolished under pre-melting conditions (Diekmann et al., 1987; Haran and Crothers, 1989; Herrera and Chaires, 1989; Chan et al., 1993; Marini et al., 1984; Jerkovic and Bolton, 2000). This indicates that the A-tract structure is energetically close to the common B-DNA geometry.

Recent studies brought attention to unusual structural and dynamical aspects of another purine tract motif, a sequence of consecutive guanines in one strand, namely the G-tract (Dornberger et al., 1999; Ng et al., 2000; Trantirek et al., 2000). G-tracts appear to have mechanical properties contrasting the A-tract (Lankas et al., 2000). G-tracts show an increasing ability to adopt an A-DNA-like structure (Cheatham et al., 1998; Lankas et al., 2000; Ng et al., 1999, 2000; Trantirek et al., 2000) depending on the tract length and environment. G-tracts are characterized by anomalously weak intrinsic base stacking due to unfavorable electrostatic and polarization contributions (Sponer et al., 2000b).

To obtain more information about the structure and dynamics of A-tracts and their influence on the curvature, experiments have been carried out with selective substitutions of certain A-T basepairs by a 6-oxopurine/cytosine (or 5-methyl cytosine) basepair, I-C (metC) (Koo and Crothers, 1987). The I-C (mC) basepair has two H-bonds similar to the A-T basepair. However, its H-bonding is substantially stronger because the electronic distribution and molecular dipole moments of I are very similar to G (Hobza and Sponer, 1999; Sponer et al., 1996a, 2000a). Thus, the I-mC basepair is sterically isostructural with the A-T basepair (see Fig. 1, A and B), while electrostatically this basepair strongly resembles the G-C one (Hobza and Sponer, 1999; Sponer et al., 1996a, 2000a). The latter fact is often ignored. The intrinsic stacking in an I-tract is weak due to the repulsive intrastrand electrostatics and polarization effects similar to a G-tract (Koo and Crothers, 1987; Shum and Crothers, 1983; Sponer et al., 2000a).

View larger version (18K): [in this window] [in a new window]   FIGURE 1   (A-D) Schematic representation of the four basepairs studied. The bottom part of each formula points to the minor groove. Note that the I-mC basepair is sterically similar to A-T while the D-T pair is isostructural with G-C. In contrast, the charge distribution and dipole moments of I-mC resemble G-C and those of D-T are similar to A-T. However, our results show that the presence or absence of the minor groove amino group rather than the electrostatic similarity decides the conformation of the tracts.

Interestingly, while a common base substitution within the A-tract usually leads to a drastic reduction of the curvature, the AI substitution affects the gel migration only slightly, and AAIAA and AIAIA sequences only modestly reduce the curvature (Diekmann et al., 1987; Koo and Crothers, 1987; Mollegaard et al., 1997). This indicates that the A-tract structure is not stabilized by major groove-bifurcated H-bonds, as these cannot be formed in the case of AIAIA tracts (Diekmann et al., 1992; Koo and Crothers, 1987). Nevertheless, the AI step is capable of forming a close contact (3.1-3.2 Å) of the two amino groups across the major groove (Luisi et al., 1998; Shatzky-Schwartz et al., 1997; Sponer and Kypr, 1994). The close amino group contact, which was first identified in AT steps (Sponer and Hobza, 1994; Sponer and Kypr, 1994), can be stabilized by asymmetrical pyramidalization of the two contacted amino groups (Hobza and Sponer, 1999; Luisi et al., 1998; Sponer and Hobza, 1994).

The I-tract itself does not lead to a marked curvature (Koo and Crothers, 1987). Nevertheless, solution experiments suggest, in line with fiber diffraction data (see above) that the I-tract is characterized by a very narrow minor groove similar to A-tract (Bailly et al., 1995; Diekmann et al., 1992). Therefore, it appears that the I-tract has inherently the same structure as the A-tract, while its 3' end junction is different from that of the A-tract. Here the very different electrostatic potentials of A-T and I-C basepairs seem to be important (Sponer et al., 2000a). In contrast, an oligonucleotide x-ray study shows a somewhat reduced propeller twist in an I₃mC₃ tract sequence compared to A₃T₃ and AIATmCT (Shatzky-Schwartz et al., 1997). However, more high-resolution data on longer I_n tract sequences would be necessary to achieve an unambiguous picture.

Another substitution useful in studies of A-tracts is an insertion of a 2,6-diaminopurine-uracil (thymine) basepair, D-U(T). The D-T basepair is sterically similar to the G-C basepair (see Fig. 1, C and D), and has three hydrogen bonds, but these are substantially weaker than in the G-C basepair due to the very low polarity of D (Sponer et al., 1996a). Phased D-tracts produce no curvature, the D-tract appears to have a wide minor groove, and a substitution of a single D into an A-tract has a detrimental effect on the curvature (Bailly and Waring, 1998; Koo and Crothers, 1987; Mollegaard et al., 1997). All the experimental data suggest that the amino group in the minor groove has a very significant effect on the structure of oligo and polypurine tracts. At the same time, the experiments indicate that the pyrimidine methyl group in the major groove does not have a substantial effect on the curvature (Koo and Crothers, 1987).

In the present paper we study properties of polypurine tracts (N-tracts) using large-scale molecular dynamics simulations with explicit inclusion of solvent and counterions and with an accurate treatment of electrostatic interactions (Beveridge and McConnell, 2000; Cheatham and Brooks, 1998; Cheatham and Kollman, 2000). This computational technique represents a powerful tool to study DNA structure and dynamics at the atomic resolution level. We have carried out simulations of d[GCG(G)₁₁GCG]₂, d[GCG(I)₁₁GCG]₂, d[GCG(A)₁₁GCG]₂, d[GCG(D)₁₁GCG]₂, and d[GCG(AI)₅AGCG]₂ duplexes utilizing I-mC and D-T basepairs. The five structures are designated as G-tract, I-tract, A-tract, D-tract, and AI-tract throughout this paper. Each trajectory has been extended to a comparably long time of 5 ns, and the A-, I-, and G-tract trajectories were prolonged up to 20 ns. The simulated structures contain an 11-bp N-tract (a full repeat with 10-bp steps) in its center capped by two GCG duplex segments to protect the N-tracts from end effects. The primary aim of our study is to investigate intrinsic properties of N-tracts, and mainly the effect of the molecular interactions of nucleobases on the N-tract architecture. Our simulations are not directly designed to study DNA bending, as we did not investigate either the tract phasing or tract-nontract junctions. However, a proper understanding of the intrinsic structure and dynamics of N-tracts themselves is a key step toward understanding their role in DNA bending. As will be shown below, simulations of long tracts as utilized in this study point out their intrinsic properties, allow for a substantial increase of the effective time scale of the study, and improve convergence of the simulated trajectories.

One of the main tasks of our study was to investigate the capability of contemporary MD techniques to reproduce the intrinsic properties of the N-tracts as compared to data from x-ray, fiber diffraction, footprinting, and other experiments. Despite the great recent advances, the accuracy of current simulation techniques is still limited by the force field approximations and the nanosecond time scale of the simulations. Because the formation of the distinct architecture of the A-tract substate is likely due to a very subtle balance of different contributions, the outcome of the simulations can be affected by force field approximations, which again justifies our choice of long N-tracts. The long N-tracts allow us to examine their intrinsic properties without perturbation from proximal tract-nontract junctions, which may represent a separate difficult task for the simulation technique. Nevertheless, to address the force field issue directly, we have carried out additional simulations of d(A)₁₀ utilizing four major force fields currently used in DNA simulations (Cornell et al., 1995; Foloppe and MacKerell, Jr., 2000; Cheatham et al., 1999; Langley, 1998; Wang et al., 2000). The results confirm that conclusions of the present study do not strongly depend on the force field choice.

One-nanosecond-length simulations of a long A- and G-tract (10 bp) have been reported before as a part of an effort to investigate the B-A DNA equilibrium (Cheatham et al., 1998). Other groups attempted to characterize properties of short A-, G-, and I-tracts (Feig and Pettitt, 1998; Pastor et al., 1999; Sherer et al., 1999; Sprous et al., 1999; Trantirek et al., 2000; Young and Beveridge, 1998; Strahs and Schlick, 2000, Sponer et al., 2000b). Our study for the first time considers the D-tracts. Further, the extent of our simulation is considerably larger than in the preceding studies.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Isothermal-isobaric (NPT) ensemble molecular dynamics were performed on each of the five 17-bp duplex oligonucleotides together with water and neutralizing counterions. The AMBER 5.0 package with the Cornell et al. (1995) force field was used. The modified bases were parametrized by means of the RESP methodology (Bayly et al., 1993; Cornell et al., 1993). At first, canonical B-DNA was built using the Nucgen module in the AMBER package; 32 sodium ions were then added to maintain charge neutrality. For this purpose, the Cion procedure of the AMBER Edit module was used, which constructs a crude 1 Å grid and then iteratively places ions at the points of least electrostatic potential. After adding the ions the whole system was immersed in a rectangular box of TIP3P water molecules (Jorgensen et al., 1983). The box dimensions were chosen to achieve a minimum distance of 10 Å from all DNA atoms to the walls, resulting in a box of ~41 Å × 41 Å × 78 Å with ~4200 water molecules. This leads to a counterion concentration of ~0.4 M (based on box dimensions). The system was then equilibrated through a series of energy minimizations and short MD runs, and heated to 300 K. The particle mesh Ewald method was used to treat electrostatic interactions. The total length of each trajectory was 5 ns (20 ns for the A-, I-, and G-tract); snapshots were taken every 2 ps. The simulations were run on the IBM SP2 at the German Cancer Research Centre.

Due to incomplete energy conservation and small losses of energy (due to application of constant pressure, finite SHAKE tolerances, small truncations in the direct space, interaction and spline-fitting of the reciprocal interactions to a grid), when coupled with velocity scaling to maintain temperature, a small growth in center of mass translation in periodic systems can appear (Harvey et al., 1998, Chiu et al., 2000). To avoid this, it is necessary to periodically remove center of mass translation, and we did remove it every 10 ps of our simulations. However, it is not necessary to remove the rotation (and it is in fact ill-defined because it will add a net-torque to the overall system). Although rotation can be removed through the addition of weak restraints, we did not remove the center of mass rotation. This could, in principle, lead to a close approach of the periodic images. However, although our DNA rotated somewhat in the box during dynamics, as is expected (due to the relatively fast reorientation times of small duplexes in solution), not more than one or two terminal basepairs left the box for several tens of picoseconds. Since the three terminal basepairs were excluded from the analysis and their structure remained intact, we conclude that the fragment rotation did not unduly influence the structure or fluctuations.

The average structures were calculated from each trajectory (with the exclusion of the first nanosecond) using the Carnal module of the AMBER package. The SCHNAaP code (Lu et al., 1997), which implements the CEHS helicoidal analysis scheme (El Hassan and Calladine, 1995), was used to obtain conformational parameters of the averaged structures. The minor groove width was calculated as the minimum distance from the phosphorus atom in the first strand to a P atom in the second strand reduced by 5.8 Å for phosphate van der Waals radii (Lu et al., 1997). Only the central 11-bp part of each oligonucleotide was analyzed; the terminal GCG fragments were omitted from the analyses.

The elastic properties of the oligonucleotides were evaluated by analyzing the correlations of fluctuations of structural variables along the trajectories. The method is described in detail elsewhere (Lankas et al., 2000). Each snapshot was first analyzed by the CURVES conformational analysis program (Lavery and Sklenar, 1988, 1989), which has the unique property of constructing an optimal, curvilinear global helical axis of the DNA fragment. Four deformation variables were defined: the fragment length (taken as the global axis length), the total twist (calculated as the sum of twists of the individual basepair steps), and two bending angles describing the bending of the helical axis into the grooves of the central basepair and to the perpendicular direction, respectively. Ensemble correlations of the deformation variables were calculated, and the known relation between the correlations of quantities with Gaussian distribution and elements of the stiffness matrix (Landau and Lifshitz, 1980) was utilized to obtain the harmonic elastic constants, including anisotropic bending and all the coupling terms. The isotropic bending constant, A_iso, was calculated as the harmonic mean of the two anisotropic ones, A₁ and A₂: 2/A_iso = 1/A₁ + 1/A₂. Again, only the central 11-bp portion of the oligomers was analyzed.

A separate study was performed to investigate the force field dependence of the results. In this case, 7-ns length simulations of a model of polyA-polyT consisting of 10 consecutive basepairs were run with the Cornell et al. (Cornell et al., 1995), MacKerell all_27 (Foloppe and MacKerell, Jr., 2000), and Langley BMS (Langley, 1998) force fields using CHARMM (Brooks et al., 1983), version c26n1. Canonical models were built using the Nucgen module of AMBER 5.0. Both A- and B-forms were considered as the starting structures. Net-neutralizing ions were initially placed 4 Å perpendicularly off the phosphate oxygen bisector and then relaxed via a gas phase 1000-step adopted-basis set Newton-Raphson minimization (with the DNA held fixed and the radii of the ions artificially increased to 5.0 Å to mimic the hydration shell). This complex was then solvated with TIP3P water (Jorgensen et al., 1983) in a rhombic dodecahedron (angles 60°, 90°, 60°) with side lengths of ~50 Å. In all of the CHARMM calculations described, the particle mesh Ewald method (Darden et al., 1993) was applied throughout with fairly stringent accuracy, including a 54 × 54 × 54 FFT charge grid, an order of six for the B-spline interpolation,  = 0.34, and a 10 Å cutoff (with the pairlist built out to 12 Å and updated heuristically). The van der Waals interactions were smoothly shifted to zero at the cutoff. During all the molecular dynamics simulations with CHARMM, a 2-fs time step, SHAKE (Ryckaert et al., 1977), and constant temperature and pressure were applied with the Langevin Piston (Feller et al., 1995) (pressure = 1.0 atm, piston mass = 500, piston gamma = 20.0) and Nosé-Hoover (mass = 1000) methods (Nosé, 1984; Hoover, 1995). The calculations were run on multiple processors of a cluster of Pentium processors (on the National Institutes of Health LoBoS I and LoBoS II computers, http://www.lobos.nih.gov) using a MPI parallelized version of CHARMM.

The first round of equilibration was performed with the Cornell et al. force field and involved the imposition of harmonic best-fit restraints (with a force constant of 25.0 kcal mol¹ Ų) on all the DNA atoms. All other atoms were free to move. Minimization was performed to relax the water and ions and applied 250 steepest descent (sd) steps followed by 750 steps applying the adopted-base Newton Raphson (abnr) algorithms implemented in CHARMM. After this, 5 ps of dynamics were performed ramping the temperature from 50 to 300 K in 1-ps intervals, followed by an additional 20 ps of dynamics at 300 K at constant pressure and temperature. This was followed by 500 steps of abnr minimization and an additional 50 ps of dynamics at 300 K. This final set of coordinates was used as the starting point for each of the force fields and followed by further equilibration that involved rounds of minimization (250 abnr steps) and 5 ps of dynamics at 300 K, with the harmonic best-fit restraint force constant set at 20.0, 15.0, 10.0, and 5.0 kcal/mol/Ų and then one final round without any restraints applied. After this, unrestrained production dynamics were started under the same conditions as described above. Center of mass translation was removed every 5 ps. Helicoidal parameters were calculated for the average structures corresponding to five 2-ns portions of each trajectory (1-3, 2-4, 3-5, 4-6, and 5-7 ns). An average over all the basepairs in the five structures was then calculated. Besides that, the standard deviations of helicoidal parameters were computed using the 1-7-ns portion of each trajectory from each configuration at 1-ps intervals. Dials and Windows were used for the calculation of helicoidal parameters; these are very similar to the values calculated via SCHNAaP. The reported minor groove widths are averages over all frames, at 1-ps intervals, for the five internal closest phosphate-phosphate distances across the minor groove, minus 5.8 Å for the phosphate radii.

In addition, 7-ns-length simulations of the 10-mer polyA-polyT were run in AMBER version 6.0 using the parm99.dat force field starting from Nucgen-generated models of A-DNA and B-DNA. Ions were placed via LEaP using a crude electrostatic potential estimation and the DNA was solvated into truncated octahedral unit cells with sides of ~55 Å. A less stringent equilibration protocol was used that involved 1000 steps of steepest-descent minimization and 100 ps of dynamics (ramping the temperature from 100 to 300 K in the first 10 ps) with the DNA held fixed. This was followed by unrestrained minimization of 1000 steps and then production dynamics at 300 K. In all the simulations, the particle mesh Ewald method was applied with a 54 × 54 × 54 FFT charge grid, order = 4,  = ~0.30768, and a 9 Å cutoff. The pairlist was updated heuristically and build out to 10 Å. A 2-fs time step was applied and temperature and pressure were maintained with Berendsen coupling (Berendsen et al., 1984), with both coupling times of 1 ps. Center of mass translation was removed every 5 ps (Harvey et al., 1998, Chiu et al., 2000).

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

RMSD and geometries of the N-tracts

Analysis of root-mean-square deviations

View larger version (51K): [in this window] [in a new window]   FIGURE 2   (A-E) Root-mean-square deviations (RMSD) of the atomic positions along the molecular dynamics trajectories. Only the central 11-bp portions of the oligomers, i.e., the tracts under study, were analyzed. The GCG ends were excluded from the analysis. RMSD with respect to canonical B-DNA (thick line), canonical A-DNA (thin line), and with respect to the average structure from the last 4 ns of the simulation (dashed line) are shown.

Conformations of the N-tracts

View larger version (26K): [in this window] [in a new window]   FIGURE 3   Stereoviews of the average structures for the five tracts. Average structures were calculated from the last 4 ns of the simulations. Only the central 9 bp from each tract are shown.

Variation of the structural parameters along the tracts

Minor-groove width and sugar puckers

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Minor groove widths of the five tracts at individual basepair steps, shown in the 5' to 3' direction of the first strand. The widths were calculated as the minimum distance from the P atom in the first strand to a P atom in the second strand, reduced by 5.8 Å to account for the phosphorus van der Waals radii. A-tract, circles; I-tract, squares; AI-tract, diamonds; G-tract, plus sign; D-tract, ×.

Base stacking energies

Table 3 summarizes intrinsic stacking energies of the NpN steps in the simulated polypurine tracts. The stacking energies are calculated using the Cornell et al. (1995) force field utilized in the simulations. The sugar-phosphate units are replaced by hydrogen atoms while charges of these hydrogens are adjusted to keep the bases neutral. The table shows mean values of base stacking energies of the inner three steps in the averaged simulated structures. The stacking energies are calculated assuming a dielectric constant of 1 and thus show the intrinsic stacking of bases in these geometries with no solvent screening (Sponer et al., 2000b). The base stacking energy consists of two parts. The van der Waals component of stacking (not shown) is overlap-dependent and includes the dominating dispersion attraction as well as steric effects. The electrostatic component of stacking (values in parentheses) reflects the mutual interaction of molecular electrostatic potentials of bases (Hobza and Sponer, 1999; Sponer et al., 1997, 2000a). Note that the Cornell et al. force field provides an excellent approximation of high-level ab initio calculations of base stacking (Hobza et al., 1997; Sponer et al., 1997), although the force field obviously neglects polarization effects (Sponer et al., 1997, 2000b).

Table 3 clearly shows that the GpG stacking is rather weak and is characterized by a substantial intrastrand electrostatic repulsion (Sponer et al., 1997, 2000b). It has been suggested that this salient stacking might play a role in establishing the distinct features of G-tracts (Sponer et al., 1997, 200b; Ng et al., 2000; Lankas et al., 2000). However, the IpI step has an almost identical balance of individual intrinsic stacking energy contributions as the GpG step, while the ApA step has entirely different distribution and magnitude of the electrostatic stacking energy terms. Considering the fact that the simulated IpI steps adopt almost identical geometries as the ApA steps, we can conclude that the geometry of NpN tracts is not primarily determined by the electrostatic component of stacking. The DpD stacking decomposition is very similar to the ApA one. The structure of DpD steps is, however, much closer to the GpG steps, and this confirms the role of the amino group protruding into the minor groove as well as the likely limited effect exerted by the pyrimidine methyl group. Finally, the ApI and IpA steps of the AI-tract show rather significant difference of their individual electrostatic stacking energy terms, resulting also in different total stacking energies. This is accompanied by rather subtle local conformational variations (see above). Thus our results indicate that the presence of efficient solvent screening reduces the role of the electrostatic components of stacking (Florian et al., 1999; Sponer et al., 2000b; Friedman and Honig, 1992; Luo et al., 2001). Note that solvent screening is capable of eliminating electrostatic repulsion even between two consecutive protonated basepairs, at least in some DNA architectures (Gallego et al., 1999; Spackova et al., 1998).

In contrast to the variable electrostatic stacking terms, the van der Waals contributions (not shown in the table) are essentially independent of the sequence. Their total values are highly conserved, ranging from 16.6 to 16.9 kcal/mol for the six different NpN steps. There are moderate variations in the Pyr-Pyr and Pur-Pyr(3') terms, correlated with the presence or absence of the N2-amino group.

In addition to the base-base interactions it would be possible to analyze other contributions to the potential energy as well, e.g., sugar-base interactions. However, the interbase interactions most probably contain the majority of sequence-dependent effects.

Interaction of the tracts with monovalent ions

Several preceding MD studies have revealed that monovalent cations are capable of penetrating into the DNA grooves to establish direct interaction with nucleobases and even to affect the fine aspects of the DNA local architecture (Hamelberg et al., 2000, 2001; Spackova et al., 2000; Stefl and Koca, 2000; Young et al., 1997a, b; McConnell and Beveridge, 2001; Cheatham and Young, 2001). Specifically, the penetration of monovalent cations into the minor groove at ApT B-DNA steps has been studied in detail (Hamelberg et al., 2000, 2001; Stefl and Koca, 2000; Young et al., 1997a, b; McConnell and Beveridge, 2000). Of special importance also is the stabilization of DNA quadruplex molecules by direct interaction of monovalent cations with the guanine quartets (Spackova et al., 1999, 2001; Stefl et al., 2001). Thus, we have systematically monitored the direct interactions of sodium cations with nucleobases in the course of our simulations in the period 3-5 ns (we omitted the first three nanoseconds of simulations to provide a little more time for the ions to equilibrate before analyzing their positions). To our surprise, we found almost no interaction of sodium cations with the minor groove edges of bases in our simulated N-tracts. We have found only a single such event of a cation coordination lasting ~0.1 ns (with a cation interacting directly with O₂ atoms of cytosines during the G-tract simulation). Instead, we have found more frequent penetration of cations into the major groove. For the A-tract we found (in the period of 3-5 ns) six major contacts lasting 0.1-0.6 ns (~0.25 ns on average) distributed over the whole tract and involving six distinct cations. For the D-tract, we have evidenced four major contacts (0.1-0.35 ns) involving three distinct cations, and for the I-tract, four major contacts (0.1-0.2 ns) involving three different cations. More major-groove direct cation-base contacts have been seen for the G-tract (13 contacts on a time scale of 0.1-0.55 ns involving 7 cations) and AI-tract (12 contacts lasting again 0.1-0.55 ns and involving 7 different cations). The direct major-groove cation-base contacts preferably involve the N7 position of purine bases and less frequently the carbonyl oxygen atoms, when available. Despite the fact that direct cation-base contacts in the major groove are more frequent, we do not believe they critically influence the structure of the simulated N-tracts, because for the major portion of time the bases remain uncontacted with the cations. Furthermore, while the cation-base contacts are most abundant in case of the AI- and G-tracts, these molecules adopt substantially different architectures.

To further investigate the interaction of the ions, the A-tract, G-tract, and I-tract simulations were extended to 20 ns. Analysis of the ion interaction over the time range of 5-20 ns confirmed our previous observations. Specifically, very little direct interaction of sodium ions was observed directly to bases in the minor groove (defined as interaction of 3.5 Å or less to N3, O₂, and N2-H22 atoms) for each of the sequences. The only consistent set of interactions in the minor groove observed are with the phosphate O1P atoms (which point into the minor groove) where direct ion interactions are observed throughout the sequence with occupancies in the 5-15% range (with lifetimes on the order of 50 ps or less and maximum lifetimes of 300-400 ps). Consistent with the analysis over shorter time frames, most of the ion interaction (beyond direct interaction with the phosphates) is transient interaction with bases (primarily N7 atoms) in the major groove. Occupancies in the range of 5-10% are common throughout the sequence with measured lifetimes of typically <200 ps with frequent exchange of ions.

The absence of cations in the minor groove of the A-tract might appear surprising, especially given the observations of direct ion interaction seen by others with the same force field and equivalent simulation protocol (Hamelberg et al., 2000; Young et al., 1997a; McConnell and Beveridge, 2001; Cheatham and Young, 2001). However, ion localization is likely influenced by issues including the starting geometry (which in our case involved placing the ions initially at favorable electrostatic potential regions), the cation concentration (around 350 mM), the ion parameters, and the time-length of the simulation. To investigate issues related to the time scale and initial ion placement we extended the simulations to 20 ns. In the longer simulations we continued to observe very low direct occupancy of Na⁺ ions in the minor groove, along with significant exchange of the ions, suggesting that the initial ion placement and short (3-5 ns) time scale was not an issue. This could be a factor in the earlier ApT studies, where ion localization at the ApT steps was seen very early in the simulations (Hamelberg et al., 2000; Young et al., 1997a); however, a higher occupancy at these steps is expected due to its low negative electrostatic potential (Lavery and Pullman, 1985). A more likely cause for low ion occupancy in the current set of simulations is the low cation concentration. Although the cation concentration in the simulations is on the order of 350 mM, this only represents enough salt to directly neutralize the DNA and more realistically represents the case of "no added salt." Without additional salt, direct interaction with the minor groove bases may be disfavored. This is consistent with recent work investigating similar A-tract sequences (McConnell and Beveridge, 2001), where in significantly shorter (3 ns) simulations low (<2.5%) occupancies were observed at ApA steps in the minor groove under similar, net-neutralizing, salt conditions. It is also consistent with observations made by Strahs and Schlick for their A6 tract simulations (Strahs and Schlick, 2000). At higher added salt concentrations, on the order of ~200 mM added Na⁺/Cl, significantly more ion occupancy in A-tract regions is observed (~25% or more) (Cheatham and Young, 2001). Based on our observations, and under a low salt limit, we can nevertheless conclude that the structural and elastic properties of the N-tracts, as revealed by the present simulations, are not affected by direct cation-base contacts in the minor groove. Despite less direct association in the minor groove than expected, it should be nevertheless noted that the lack of Na⁺ association with the bases in the minor groove is not in disagreement with NMR studies that show a considerable presence of NH₄⁺ there. As noted by Denisov and Halle (2000), NH₄⁺ appears to bind to the A-tracts by an order of magnitude more strongly than Na⁺. Note also that the current experimental evidence concerning direct association with nucleobases in the minor groove of B-DNA is somewhat conflicting, and it may well be that the degree of cation association is highly dependent on the type of cation, concentration, temperature, and other factors (Denisov and Halle, 2000; McConnell and Beveridge, 2000; Stellwagen et al., 2001).

Elastic properties of the tracts

There are important biological processes in which not only the structure, but also the mechanical deformability of DNA plays a crucial role. These include specific DNA-binding to proteins or processes related to long-range interactions and higher-order organization of DNA in the nucleus. Many such biological functions occur within the linear elasticity regime. In that case the elastic energy can be expressed as a quadratic function of appropriate deformation variables.

Various methods have been devised to theoretically calculate the elastic constants. One of the most promising is the "method of fluctuations": if an ensemble of deformation variables with Gaussian distribution is available, then the correlations of these variables are directly related to the harmonic force constants (Landau and Lifshitz, 1980). Olson et al. (1998) used this approach to calculate empirical deformabilities of dinucleotide steps based on an ensemble of x-ray protein-DNA complexes. They obtained force constants within an unknown multiplicative factor, the effective temperature of the ensemble.

It is, however, also desirable to know elastic properties of much larger portions of DNA than basepair steps. In a recent attempt, Bruant et al. (1999) used ensembles of snapshots from 1-ns molecular dynamics trajectories to calculate the isotropic bending, twisting, and stretching rigidities of two 15-bp oligomers. While their bending and twisting constants were in reasonable agreement with experiment, the stretch modulus was much higher than expected from the macroscopic experiments.

In our previous work (Lankas et al., 2000) we performed 5-ns molecular dynamics runs on four 17-bp DNA oligomers with different sequences and calculated the elastic constants including anisotropic bending and all the coupling terms. In particular, the value of twist-stretch coupling was established, and the twist-bend coupling constant, whose existence was theoretically predicted some time ago (Marko and Siggia, 1994; O'Hern et al., 1998), was calculated for the first time. Although the values of the force constants were in very good overall agreement with experimental values for random sequences, a pronounced sequence-dependence of elastic properties was observed, especially regarding the stretch modulus.

Large-scale molecular dynamics simulations thus represent a unique source of data to study DNA elasticity, including its sequence dependence. Here, we use the method from our previous work to compare the deformabilities of the five tracts under study.

The elastic constants are shown in Table 4. The values of the stretch modulus, Y_f, strongly confirm the similarity of A-, I-, and AI-tracts versus G- and D-tracts. The former are much stiffer than the latter, with AI and D as extreme examples. Just the opposite, yet not so markedly pronounced, holds for the twist persistence length C: the former group is in general softer with respect to twisting than the latter one. No clear tendency has been observed in the case of the isotropic bending persistence length A_iso; here the AI-tract alone stands out as the most flexible. It thus seems that the bending stiffness is not a simple consequence of the presence or absence of the C2 amino group. The same applies to the twist-stretch coupling D.

The coupling between twist and bending into the major groove of the tract's central basepair, A_1t, is the most important coupling term for fragments shorter than one helical turn. On the dinucleotide level it corresponds to the anticorrelation between twist and roll, a phenomenon observed in the crystallographic data (Gorin et al., 1995; Olson et al., 1998). This coupling diminishes as the fragment length approaches one helical turn; therefore we report in the table the values for 6-bp fragments (averaged over all possible 6-mers of the tract). All the other coupling terms are small and are not reported in the Table.

As far as the anisotropic bending is concerned, it is markedly pronounced for small numbers of basepairs in a fragment; our calculations show that about two times less energy is needed to bend a trimer into the grooves of its central basepair than in the perpendicular direction. This is in line with the view that changes in roll should be much easier to perform than changes in tilt (Olson et al., 1998; Zhurkin et al., 1979). This anisotropy vanishes at the half of the helical turn, and is quite small at the full length of our tracts, in accordance with our previous results (Lankas et al., 2000).

Besides the global elastic properties mentioned above, one might be interested in the DNA deformability on the basepair or basepair step level. The corresponding force constants may be inferred from the molecular dynamics trajectories as well (Lankas et al., manuscript in preparation).

Structure of d(A)₁₀: Comparison of four force fields

The lower than experimentally observed propeller twist for the A-tract sequence studied (see above) led us to extend our study and to investigate the force field dependence and possible bias in the observed values. To characterize this, simulations of a 10-mer polyA-polyT were performed with four of the most reliable and currently available empirical, all-atom nucleic acid force fields. These include the Cornell et al. (1995), Cheatham et al. (parm98/99) (Cheatham et al., 1999; Wang et al., 2000), the MacKerell all_27 nucleic acid force field (Foloppe and MacKerell, Jr., 2000), and Langley's BMS force field (Langley, 1998). All simulations were carried out for 7 ns. Furthermore, to verify the effect of the starting structure, the simulations were carried out taking both A- and B-DNA structures as starting points. Thus, we have carried out altogether ~60 ns of simulations of the 10-mer polyA-polyT. With each of the force fields, a conversion from A-DNA to B-DNA is seen (data not shown) with convergence to a common set of structures for each applied force field. Thus, in the following, we present only data obtained for the structures initiated in the B-DNA form geometry, and these are summarized in Table 5. Comparing the B-DNA simulation results across the force fields suggests that all of these force fields tend to underestimate the average propeller twist compared to crystallographic data. The best behavior, arguably, comes from the BMS force field, which shows propeller twisting even slightly larger than suggested recently by NMR spectroscopy (MacDonald et al., 2001). It should be noted that in the simulated 10-mer polyA-polyT structures a certain variation is seen in the values at each basepair, suggesting that there is still somehow insufficient sampling in the current set of 7-ns simulations. As has been noted by us and others (Feig and Pettitt, 1998; Spackova et al., 2001; Sponer et al., 2000b), nanosecond-length simulations may not be sufficient to fully converge the structures. There are relatively large variations (standard deviations, see Table 5) observed during the simulation for the propeller twist, compared to helical twist, suggesting a wide range of motion of the basepairs. Note, however, that the structures of the N-tracts in our 17-bp duplex simulations studied herein and discussed previously are considerably more uniform. This may be related to the presence of terminal GCG segments shielding the core N-tract from end-effects. Nevertheless, convergence issues aside, all the force fields show similar trends.

It is difficult to rationalize at this point whether the lower than expected average propeller twist in the simulated A-tracts is due to lack of convergence or insufficient sampling, a potential bias in the experimental structures, or a systematic bias of the force field. The latter is distinctly possible due to the neglect of hydrogen bond directionality, polarization effects, or other force field approximations such as neglect of anisotropy of the short-range repulsion (Sponer et al., 1996b). However, despite the observation of lower propeller twist values, it is reassuring that all three different force field philosophies tend to give comparable results and show similar trends. A comparison of the Cornell et al. results to the other force fields suggests that the root of the "lower than expected" propeller twist in the Cornell et al. force field is not simply due to the known deficiencies of this force field. Although the Cornell et al. force underestimates the average sugar pucker phases and helical repeat, the MacKerell all_27 force field, which does a much better job in this regard, shows larger than expected minor groove widths and lower than expected propeller twist. As the parm98/99 results suggest a larger magnitude for the propeller twist, clearly the sugar pucker phases and helical repeat (which are improved in this force field) are involved. The increase in propeller twisting compared to the Cornell et al. force field, however, is not substantial. Despite the deficiencies and absence of specific hydrogen bond directionality and explicit polarization, it is reassuring that the BMS force field, which uses the Cornell et al. charge model and all_22 internal (bond/angle) parameters and an iteratively improved dihedral potential, does better. This suggests that with further effort to balance and tweak these simple empirical potentials it will be possible to generate an even more reliable (and general) nucleic acid force field.

Nevertheless, the most important conclusion regarding the purpose of the present study is that all the applied force fields show qualitatively similar results. Thus, the main results with the Cornell et al. force field described herein, regarding the intrinsic flexibility and variations in helicoidal parameters and minor groove widths across the polypurine sequences studied, seem to be general.

	DISCUSSION AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

We have performed 5-20-ns unrestrained simulations of d[GCG(N)₁₁GCG]₂ duplex structures. The simulations have been carried out with explicit inclusion of solvent and counterions and with the particle mesh Ewald method for an accurate treatment of long-range electrostatic forces. This technique is currently the most robust computational tool available for studies of structure and dynamics of oligonucleotides. We have used the Cornell et al. force field, which is known to provide a realistic description of base stacking interactions (Hobza et al., 1997).

The simulation reproduces many essential features of the A-tracts. Mainly, it shows that the A-tract structure is substantially propeller-twisted, has a narrow minor groove, and the width of the minor groove systematically decreases along the tract in the 3' to 5' direction. The narrowing of the groove appears to be saturated after approximately six consecutive adenines. The A-tract is found to be rigid with respect to stretching and relatively flexible with respect to twisting, while its bending rigidity does not seem to differ much from that of the other tracts.

However, the simulation also shows rather substantial differences compared to the experimental findings. The propeller twist is (in absolute value) underestimated, being 12° in the simulations, while oligonucleotide crystals and fibers show a value around 20° and 26°, respectively. The calculated helical twist of the A-tract is underestimated by ~3-4°. Also, the narrowing of the minor groove in the simulated structure is not as significant as in the experiments. These observations are in agreement with literature data (Cheatham and Kollman, 2000). These data indicate that while qualitatively correct, the force field is not capable of quantitatively reproducing some fine structural aspects of the A-tract. The most likely reasons might be the imbalance of the parameters for the sugar-phosphate backbone, imbalance of solute-solvent interactions, possible nonplanarity of the N2 amino group, or the absence of polarizabilities in the force field. Nevertheless, the simulation shows key features of the A-tracts and, therefore, it is justified to carry out simulations with different sequences. The relative differences among various sequences should be affected by the force field inaccuracies to a substantially smaller degree, as all structures should be biased in the same direction, leaving sequence differences fully expressed. Furthermore, it has been explicitly demonstrated that the Cornell et al. force field provides a very good and well-balanced description of the base-base interactions (Hobza et al., 1997).

It should be noted that our simulations have been carried out with the original version of the Cornell et al. force field (Cornell et al., 1995). Recently, new variants of the force field, parm98 (Cheatham et al., 1999) and parm99 (Wang et al., 2000) have been released, with modified sugar pucker parameters. This influences the B-A equilibria provided by the force fields, helical twist parameters, and some other features of the simulated structures. Furthermore, there exist two other force fields designed for DNA simulations, BMS (Langley, 1998) and CHARMM (Foloppe and MacKerell, Jr., 2000). This has inspired us to carry out simulations of a long adenine d(A)₁₀ tract using all four parametrizations. These additional simulations show that, at least considering the presently available generation of force fields, our results obtained with the original Cornell et al. force field are general. These results are also in line with our preceding paper (Spackova et al., 2000), where we have found almost no difference in simulations of an adenine-zipper duplex with parm94 and parm98. This view is also in agreement with results obtained by others (Cheatham and Kollman, 2000). It is to be noted that all three AMBER force field variants have identical parametrization of DNA bases, which is the most important aspect of the force field for the purpose of the present study. Thus the choice of the force field (considering presently available parametrizations) should have no substantial influence on the outcome and conclusions of the present paper.

Simulations of d[GCG(I)₁₁GCG]₂ and d[GCG(AI)₅AGCG]₂ duplexes show unambiguously that both I- and AI-tracts adopt a high-propeller structure with a narrow minor groove width very similar to that of A-tract. This finding is in full agreement with all available experimental data (Bailly et al., 1995; Diekmann et al., 1992; Leslie et al., 1980). The only experimental result that might be in disagreement with our data is a crystal structure of a short medium-resolution I₃mC₃ tract (Shatzky-Schwartz et al., 1997). An analogous crystal structure with AIATmCT sequence shows a high propeller twist (Shatzky-Schwartz et al., 1997), entirely in agreement with our present data. The mechanical deformability of the I- and AI-tract is characterized by a high stretch modulus and rather low twist persistence length, again entirely in accordance with the properties of the A-tract.

The A-T and I-mC basepairs have the same van der Waals contour; however, they differ substantially in their molecular electrostatic potentials. As a result, AA and II steps differ considerably in the balance of intrastrand and interstrand contributions to stacking and in the interplay between van der Waals and electrostatic contributions (Sponer et al., 1997, 2000a). The base stacking in the II step is weak and is characterized by an intrastrand electrostatic repulsion. The electrostatic component of stacking in an I-tract is considerably closer to stacking in a G-tract than to stacking in an A-tract. The alternating AI structure, where rather nonpolar AT basepairs are interdispersed between very polar IC basepairs, has electrostatic interactions obviously closer to the A-tract sequence. In our simulations, the structure of the I-tract remains similar to the A-tract, and the same is true for the AI-tract. Thus, the simulations clearly show that the steric effects associated with the presence or absence of the minor-groove amino group within the tract is the critical factor governing the formation of A-tract-like structure, while the electrostatic differences of the stacking remain unexpressed. It does not rule out, however, that the substantial difference in molecular electrostatic potentials of A-T and