Structural equilibrium of DNA represented with different force fields

Formation pathways of a guanine-quadruplex DNA revealed by molecular dynamics and thermodynamic analysis of the substates

The formation of a cation-stabilized guanine quadruplex (G-DNA) stem is an exceptionally slow process involving complex kinetics that has not yet been characterized at atomic resolution. Here, we investigate the formation of a parallel stranded G-DNA stem consisting of four strands of d(GGGG) using molecular dynamics simulations with explicit inclusion of counterions and solvent. Due to the limitations imposed by the nanosecond timescale of the simulations, rather than watching for the spontaneous formation of G-DNA, our approach probes the stability of possible supramolecular intermediates (including two-, three-, and four-stranded assemblies with out-of-register base pairing between guanines) on the formation pathway. The simulations suggest that "cross-like" two-stranded assemblies may serve as nucleation centers in the initial formation of parallel stranded G-DNA quadruplexes, proceeding through a series of rearrangements involving trapping of cations, association of additional strands, and progressive slippage of strands toward the full stem. To supplement the analysis, approximate free energies of the models are obtained with explicit consideration of the integral cations. The approach applied here serves as a prototype for qualitatively investigating other G-DNA molecules using molecular dynamics simulation and free-energy analysis.

Titration in silico of reversible B <=> A transitions in DNA

Reversible transitions between the A- and B-forms of DNA are obtained in free molecular dynamics simulations of a single double helix immersed in a water drop with Na(+) counterions. The dynamics of the transitions agrees with their supposed cooperative character. In silico titration of the transitions was carried out by smooth variation of the drop size. The estimated range of hydration numbers corresponding to the transition roughly agrees with experimental data. The chain length dependence was studied for double helices from 6 to 16 base pairs. It appeared that the B --> A transition is hindered for DNA shorter than one helical turn. With increased NaCl concentration in the drop, stabilization of the B-form is observed accompanied by the salt crystallization. The results strongly suggest that the B --> A transition at low hydration is caused by Na(+) ions sandwiched between phosphate strands in the major groove and is driven by direct medium range electrostatic interactions. The role of the reduced water shell apparently consists of increasing the counterion concentration in the opening of the major groove. Analysis of the available experimental data suggests that this mechanism is perhaps generally responsible for the A/B polymorphism in DNA.

DNA polymorphism: a comparison of force fields for nucleic acids

The improvements of the force fields and the more accurate treatment of long-range interactions are providing more reliable molecular dynamics simulations of nucleic acids. The abilities of certain nucleic acid force fields to represent the structural and conformational properties of nucleic acids in solution are compared. The force fields are AMBER 4.1, BMS, CHARMM22, and CHARMM27; the comparison of the latter two is the primary focus of this paper. The performance of each force field is evaluated first on its ability to reproduce the B-DNA decamer d(CGATTAATCG)(2) in solution with simulations in which the long-range electrostatics were treated by the particle mesh Ewald method; the crystal structure determined by Quintana et al. (1992) is used as the starting point for all simulations. A detailed analysis of the structural and solvation properties shows how well the different force fields can reproduce sequence-specific features. The results are compared with data from experimental and previous theoretical studies.

Force field dependence of NMR-Based, restrained molecular dynamics DNA structure calculations including an analysis of the influence of residual dipolar coupling restraints

Restrained molecular dynamics is widely used to calculate DNA structures from NMR data. Here, results of an in silico experiment show that the force field can be significant compared to the NMR restraints in driving the final structures to converge. Specifically, we observed that i) the influence of the force field leads to artificially tight convergence within final families of structures and ii) the precision and character of resulting structures depend on the choice of force field used in the calculations. A canonical B-DNA model was used as a target structure. Distances, dihedral angles, and simulated residual dipolar couplings were measured in the target structure and used as restraints. X-PLOR and Discover, which use force fields developed for CHARMM and AMBER programs, respectively, were tested and found to produce different final structures despite the use of identical distance and dihedral restraints. Incorporation of residual dipolar coupling restraints in X-PLOR improves convergence with the target structure and between families of structures indicating that the force field dependence can potentially be overcome if residual dipolar coupling restraints are employed.

Theoretical methods for the simulation of nucleic acids

Different theoretical methods for the description of nucleic acid structures are reviewed. Firstly, we introduce the concept of classical force-field in the context of nucleic acid structures, discussing their accuracy. We then examine theoretical approaches to the description of nucleic acids based on: i) a rigid or quasi-rigid description of the molecule, ii) molecular mechanics optimization, and iii) molecular dynamics. Special emphasis is made ion current state of the art molecular dynamics simulations of nucleic acids structures.

Effect of a neutralized phosphate backbone on the minor groove of B-DNA: molecular dynamics simulation studies

Alternative models have been presented to provide explanations for the sequence-dependent variation of the DNA minor groove width. In a structural model groove narrowing in A-tracts results from direct, short-range interactions among DNA bases. In an electrostatic model, the narrow minor groove of A-tracts is proposed to respond to sequence-dependent localization of water and cations. Molecular dynamics simulations on partially methylphosphonate substituted helical chains of d(TATAGGCCTATA) and d(CGCGAATTCGCG) duplexes have been carried out to help evaluate the effects of neutralizing DNA phosphate groups on the minor groove width. The results show that the time-average minor groove width of the GGCC duplex becomes significantly more narrow on neutralizing the phosphate backbone with methylphosphonates. The minor groove of the AATT sequence is normally narrow and the methylphosphonate substitutions have a smaller but measurable affect on this sequence. These results and models provide a system that can be tested by experiment and they support the hypothesis that the electrostatic environment around the minor groove affects the groove width in a sequence-dependent dynamic and time-average manner.

Simulations of nucleic acids and their complexes

Recent years have seen considerable progress in simulations of nucleic acids. Improvements in force fields, simulation techniques and protocols, and increasing computer power have all contributed to making nanosecond-scale simulations of both DNA and RNA commonplace. The results are already helping to explain how nucleic acids respond to their environment and to their base sequence and to reveal the factors underlying recognition processes by probing biologically important nucleic acid-protein interactions and medically important nucleic acid-drug complexation. This Account summarizes methodological progress and applications of molecular dynamics to nucleic acids over the past few years and tries to identify remaining challenges.

Solvation and hydration of proteins and nucleic acids: a theoretical view of simulation and experiment

Many theoretical, computational, and experimental techniques recently have been successfully used for description of the solvent distribution around macromolecules. In this Account, we consider recent developments in the areas of protein and nucleic acid solvation and hydration as seen by experiment, theory, and simulations. We find that in most cases not only the general phenomena of solvation but even local hydration patterns are more accurately discussed in the context of water distributions rather than individual molecules of water. While a few localized or high-residency waters are often associated with macromolecules in solution (or crystals from aqueous liquors), these are readily and accurately included in this more general description. The goal of this Account is to review the theoretical models used for the description of the interfacial solvent structure on the border near DNA and protein molecules. In particular, we hope to highlight the progress in this field over the past five years with a focus on comparison of simulation and experimental results.

Critical effect of the N2 amino group on structure, dynamics, and elasticity of DNA polypurine tracts

Unrestrained 5-20-ns explicit-solvent molecular dynamics simulations using the Cornell et al. force field have been carried out for d[GCG(N)11GCG]2 (N, purine base) considering guanine*cytosine (G*C), adenine*thymine (A*T), inosine*5-methyl-cytosine (I*mC), and 2-amino-adenine*thymine (D*T) basepairs. The simulations unambiguously show that the structure and elasticity of N-tracts is primarily determined by the presence of the amino group in the minor groove. Simulated A-, I-, and AI-tracts show almost identical structures, with high propeller twist and minor groove narrowing. G- and D-tracts have small propeller twisting and are partly shifted toward the A-form. The elastic properties also differ between the two groups. The sequence-dependent electrostatic component of base stacking seems to play a minor role. Our conclusions are entirely consistent with available experimental data. Nevertheless, the propeller twist and helical twist in the simulated A-tract appear to be underestimated compared to crystallographic studies. To obtain further insight into the possible force field deficiencies, additional multiple simulations have been made for d(A)10, systematically comparing four major force fields currently used in DNA simulations and utilizing B and A-DNA forms as the starting structure. This comparison shows that the conclusions of the present work are not influenced by the force field choice.

Molecular dynamics simulation of nucleic acids: successes, limitations, and promise

In the last five years we have witnessed a significant increase in the number publications describing accurate and reliable all-atom molecular dynamics simulations of nucleic acids. This increase has been facilitated by the development of fast and efficient methods for treating the long-range electrostatic interactions, the availability of faster parallel computers, and the development of well-validated empirical molecular mechanical force fields. With these technologies, it has been demonstrated that simulation is not only capable of consistently reproducing experimental observations of sequence specific fine structure of DNA, but also can give detailed insight into prevalent problems in nucleic acid structure, ion association and specific hydration of nucleic acids, polyadenine tract bending, and the subtle environmental dependence of the A-DNA-B-DNA duplex equilibrium. Despite the advances, there are still issues with the methods that need to be resolved through rigorous controlled testing. In general, these relate to deficiencies of the underlying molecular mechanical potentials or applied methods (such as the imposition of true periodicity in Ewald simulations and the need for energy conservation), and significant limits in effective conformational sampling. In this perspective, we provide an overview of our experiences, provide some cautionary notes, and provide recommendations for further study in molecular dynamics simulation of nucleic acids.

Development and current status of the CHARMM force field for nucleic acids

The CHARMM27 all-atom force field for nucleic acids represents a highly optimized model for investigations of nucleic acids via empirical force field calculations. The force field satisfactorily treats the A, B, and Z forms of DNA as well as RNA, and it also useful for nucleosides and nucleotides. In addition, it is compatible with the CHARMM force fields for proteins and lipids, allowing for simulation studies of heterogeneous systems.

Sequence-dependent B<-->A transition in DNA evaluated with dimeric and trimeric scales

Experimental data on the sequence-dependent B<-->A conformational transition in 24 oligo- and polymeric duplexes yield optimal dimeric and trimeric scales for this transition. The 10 sequence dimers and the 32 trimers of the DNA duplex were characterized by the free energy differences between the B and A forms in water solution. In general, the trimeric scale describes the sequence-dependent DNA conformational propensities more accurately than the dimeric scale, which is likely related to the trimeric model accounting for the two interfaces between adjacent base pairs on both sides (rather than only one interface in the dimeric model). The exceptional preference of the B form for the AA:TT dimers and AAN:N'TT trimers is consistent with the cooperative interactions in both grooves. In the minor groove, this is the hydration spine that stabilizes adenine runs in B form. In the major groove, these are hydrophobic interactions between the thymine methyls and the sugar methylene groups from the preceding nucleotides, occurring in B form. This interpretation is in accord with the key role played by hydration in the B<-->A transition in DNA. Importantly, our trimeric scale is consistent with the relative occurrences of the DNA trimers in A form in protein-DNA cocrystals. Thus, we suggest that the B/A scales developed here can be used for analyzing genome sequences in search for A-philic motifs, putatively operative in the protein-DNA recognition.

Molecular dynamics simulations of the denaturation and refolding of an RNA tetraloop

Tetraloops are very abundant structural elements of RNA that are formed by four nucleotides in a hairpin loop which is closed by a double stranded helical stem with some Watson-Crick base pairs. A tetraloop r(GCGAAGGC) was identified from the crystal structure of the central domain of 16S rRNA (727-730) in the Thermus thermophilus 30S ribosomal complex. The crystal structure of the 30S complex includes a total of 104 nucleotides from the central domain of the 16S rRNA and three ribosomal proteins S6, S15 and S18. Independent biochemical experiments have demonstrated that protein S15 plays the role in initiating the formation of the central domain of this complex. In the crystal, the tetraloop interacts with the protein S15 at two sites: one of them is associated with hydrogen bond interactions between residue His50 and nucleotide G730, and the other is associated with the occurrence of residue Arg53 beside A728. This paper uses molecular dynamics (MD) simulations to investigate the protein-dependent structural stability of the tetraloop and demonstrates the folding pathway of this tetraloop via melting denaturation and its subsequent refolding. Three important results are derived from these simulations: (i) The stability of nucleotide A728 appears to be protein dependent. Without the interaction with S15, A728 flips away from stacking with A729. (ii) The melting temperature demonstrated by the simulations is analogous to the results of thermodynamic experiments. In addition, the simulated folding of the tetraloop is stepwise: the native shape of the backbone is formed first; this is then followed by the formation of the Watson- Crick base pairs in the stem; and finally the hydrogen bonds and base stacking in the loop are formed. (iii) The tetraloop structure is similar to the crystal structure at salt concentrations of 0.1 M and 1.0 M used for the simulations, but the refolded structure at 0.1 M salt is more comparable to the crystal structure than at 1.0 M. The results from the simulations using both the Generalized Born continuum model and explicit solvent model (Particle Mesh Ewald) generate a similar pathway for unfolding/refolding of the tetraloop.

Molecular dynamics studies of trinucleotide repeat DNA involved in neurodegenerative disorders

Expansion of trinucleotide repeat DNA of the classes CAG-CTG, CGG-CCG and GAA-TTC are found to be associated with several neurodegenerative disorders. Different mechanisms have been attributed to the expansion of triplets, mainly involving the formation of alternate secondary structures by such repeats. This paper reports the molecular dynamics simulation of triplet repeat DNA sequences to study the basic structural features of DNA that are responsible for the formation of structures such as hairpins and slip-strand DNA leading to expansion. All the triplet repeat sequences studied were found to be more flexible compared to the control sequence unassociated with disease. Moreover, flexibility was found to be in the order CAG-CTG > CGG-CCG approximately GAA-TTC, the highly flexible CAG-CTG repeat being the most common cause of neurodegenerative disorders. In another simulation, a single G-C to T-A mutation at the 9th position of the CAG-CTG repeat exhibited a reduction in bending compared to the pure 15-mer CAG-CTG repeat. EPM1 dodecamer repeat associated with the pathogenesis of progressive myoclonus epilepsy was also simulated and showed flexible nature suggesting a similar expansion mechanism.

Conformations of an adenine bulge in a DNA octamer and its influence on DNA structure from molecular dynamics simulations

Molecular dynamics simulations have been applied to the DNA octamer d(GCGCA-GAAC). d(GTTCGCGC), which has an adenine bulge at the center to determine the pathway for interconversion between the stacked and extended forms. These forms are known to be important in the molecular recognition of bulges. From a total of ~35 ns of simulation time with the most recent CHARMM27 force field a variety of distinct conformations and subconformations are found. Stacked and fully looped-out forms are in excellent agreement with experimental data from NMR and x-ray crystallography. Furthermore, in a number of conformations the bulge base associates with the minor groove to varying degrees. Transitions between many of the conformations are observed in the simulations and used to propose a complete transition pathway between the stacked and fully extended conformations. The effect on the surrounding DNA sequence is investigated and biological implications of the accessible conformational space and the suggested transition pathway are discussed, in particular for the interaction of the MS2 replicase operator RNA with its coat protein.

Improving the accuracy of NMR structures of DNA by means of a database potential of mean force describing base-base positional interactions

NMR structure determination of nucleic acids presents an intrinsically difficult problem since the density of short interproton distance contacts is relatively low and limited to adjacent base pairs. Although residual dipolar couplings provide orientational information that is clearly helpful, they do not provide translational information of either a short-range (with the exception of proton-proton dipolar couplings) or long-range nature. As a consequence, the description of the nonbonded contacts has a major impact on the structures of nucleic acids generated from NMR data. In this paper, we describe the derivation of a potential of mean force derived from all high-resolution (2 A or better) DNA crystal structures available in the Nucleic Acid Database (NDB) as of May 2000 that provides a statistical description, in simple geometric terms, of the relative positions of pairs of neighboring bases (both intra- and interstrand) in Cartesian space. The purpose of this pseudopotential, which we term a DELPHIC base-base positioning potential, is to bias sampling during simulated annealing refinement to physically reasonable regions of conformational space within the range of possibilities that are consistent with the experimental NMR restraints. We illustrate the application of the DELPHIC base-base positioning potential to the structure refinement of a DNA dodecamer, d(CGCGAATTCGCG)(2), for which NOE and dipolar coupling data have been measured in solution and for which crystal structures have been determined. We demonstrate by cross-validation against independent NMR observables (that is, both residual dipolar couplings and NOE-derived intereproton distance restraints) that the DELPHIC base-base positioning potential results in a significant increase in accuracy and obviates artifactual distortions in the structures arising from the limitations of conventional descriptions of the nonbonded contacts in terms of either Lennard-Jones van der Waals and electrostatic potentials or a simple van der Waals repulsion potential. We also demonstrate, using experimental NMR data for a complex of the male sex determining factor SRY with a duplex DNA 14mer, which includes a region of highly unusual and distorted DNA, that the DELPHIC base-base positioning potential does not in any way hinder unusual interactions and conformations from being satisfactorily sampled and reproduced. We expect that the methodology described in this paper for DNA can be equally applied to RNA, as well as side chain-side chain interactions in proteins and protein-protein complexes, and side chain-nucleic acid interactions in protein-nucleic acid complexes. Further, this approach should be useful not only for NMR structure determination but also for refinement of low-resolution (3-3.5 A) X-ray data.

Molecular dynamics simulation of nucleic acids

We review molecular dynamics simulations of nucleic acids, including those completed from 1995 to 2000, with a focus on the applications and results rather than the methods. After the introduction, which discusses recent advances in the simulation of nucleic acids in solution, we describe force fields for nucleic acids and then provide a detailed summary of the published literature. We emphasize simulations of small nucleic acids ( approximately 6 to 24 mer) in explicit solvent with counterions, using reliable force fields and modern simulation protocols that properly represent the long-range electrostatic interactions. We also provide some limited discussion of simulation in the absence of explicit solvent. Absent from this discussion are results from simulations of protein-nucleic acid complexes and modified DNA analogs. Highlights from the molecular dynamics simulation are the spontaneous observation of A B transitions in duplex DNA in response to the environment, specific ion binding and hydration, and reliable representation of protein-nucleic acid interactions. We close by examining major issues and the future promise for these methods.

Stability and structure of RNA duplexes containing isoguanosine and isocytidine

Isoguanosine (iG) and isocytidine (iC) differ from guanosine (G) and cytidine (C), respectively, in that the amino and carbonyl groups are transposed. The thermodynamic properties of a set of iG, iC containing RNA duplexes have been measured by UV optical melting. It is found that iG-iC replacements usually stabilize duplexes, and the stabilization per iG-iC pair is sequence-dependent. The sequence dependence can be fit to a nearest-neighbor model in which the stabilities of iG--iC pairs depend on the adjacent iG--iC or G--C pairs. For 5'-CG-3'/3'-GC-5' and 5'-GG-3'/3'-CC-5' nearest neighbors, the free energy differences upon iG-iC replacement are smaller than 0.2 kcal/mol at 37 degrees C, regardless of the number of replacements. For 5'-GC-3'/3'-CG-5', however, each iG--iC replacement adds 0.6 kcal/mol stabilizing free energy at 37 degrees C. Stacking propensities of iG and iC as unpaired nucleotides at the end of a duplex are similar to those of G and C. An NMR structure is reported for r(CiGCGiCG)(2) and found to belong to the A-form family. The structure has substantial deviations from standard A-form but is similar to published NMR and/or crystal structures for r(CGCGCG)(2) and 2'-O-methyl (CGCGCG)(2). These results provide benchmarks for theoretical calculations aimed at understanding the fundamental physical basis for the thermodynamic stabilities of nucleic acid duplexes.

Water and ion binding around r(UpA)12 and d(TpA)12 oligomers--comparison with RNA and DNA (CpG)12 duplexes

The structural and dynamic properties of the water and ion first coordination shell of the r(A-U) and d(A-T) base-pairs embedded within the r(UpA)12 and d(TpA)12 duplexes are described on the basis of two 2.4 ns molecular dynamics simulations performed in a neutralizing aqueous environment with 0.25 M added KCl. The results are compared to previous molecular dynamics simulations of the r(CpG)12 and d(CpG)12 structures performed under similar conditions. It can be concluded that: (i) RNA helices are more rigid than DNA helices of identical sequence, as reflected by the fact that RNA duplexes keep their initial A-form shape while DNA duplexes adopt more sequence-specific shapes. (ii) Around these base-pairs, the water molecules occupy 21 to 22 well-defined hydration sites, some of which are partially occupied by potassium ions. (iii) These hydration sites are occupied by an average of 21.9, 21.0, 20.1, and 19.8 solvent molecules (water and ions) around the r(G=C), r(A-U), d(G=C), and d(A-T) pairs, respectively. (iv) From a dynamic point of view, the stability of the hydration shell is the strongest for the r(G=C) pairs and the weakest for the d(A-T) pairs. (v) For RNA, the observed long-lived hydration patterns are essentially non-sequence dependent and involve water bridges located in the deep groove and linking OR atoms of adjacent phosphate groups. Maximum lifetimes are close to 400 ps. (vi) In contrast, for DNA, long-lived hydration patterns are sequence dependent and located in the minor groove. For d(CpG)12, water bridges linking the (G)N3 and (C)O2 with the O4' atoms of adjacent nucleotides with 400 ps maximum lifetimes are characterized while no such bridges are observed for d(TpA)12. (vii) Potassium ions are observed to bind preferentially to deep/major groove atoms at RpY steps, essentially d(GpC), r(GpC), and r(ApU), by forming ion-bridges between electronegative atoms of adjacent base-pairs. On average, about half an ion is observed per base-pair. Positive ion-binding determinants are related to the proximity of two or more electronegative atoms. Negative binding determinants are associated with the electrostatic and steric hindrance due to the proximity of electropositive amino groups and neutral methyl groups. Potassium ions form only transient contacts with phosphate groups.

Molecular dynamics of DNA quadruplex molecules containing inosine, 6-thioguanine and 6-thiopurine

The ability of the four-stranded guanine (G)-DNA motif to incorporate nonstandard guanine analogue bases 6-oxopurine (inosine, I), 6-thioguanine (tG), and 6-thiopurine (tI) has been investigated using large-scale molecular dynamics simulations. The simulations suggest that a G-DNA stem can incorporate inosines without any marked effect on its structure and dynamics. The all-inosine quadruplex stem d(IIII)(4) shows identical dynamical properties as d(GGGG)(4) on the nanosecond time scale, with both molecular assemblies being stabilized by monovalent cations residing in the channel of the stem. However, simulations carried out in the absence of these cations show dramatic differences in the behavior of d(GGGG)(4) and d(IIII)(4). Whereas vacant d(GGGG)(4) shows large fluctuations but does not disintegrate, vacant d(IIII)(4) is completely disrupted within the first nanosecond. This is a consequence of the lack of the H-bonds involving the N2 amino group that is not present in inosine. This indicates that formation of the inosine quadruplex could involve entirely different intermediate structures than formation of the guanosine quadruplex, and early association of cations in this process appears to be inevitable. In the simulations, the incorporation of 6-thioguanine and 6-thiopurine sharply destabilizes four-stranded G-DNA structures, in close agreement with experimental data. The main reason is the size of the thiogroup leading to considerable steric conflicts and expelling the cations out of the channel of the quadruplex stem. The G-DNA stem can accommodate a single thioguanine base with minor perturbations. Incorporation of a thioguanine quartet layer is associated with a large destabilization of the G-DNA stem whereas the all-thioguanine quadruplex immediately collapses.

Local conformational variations observed in B-DNA crystals do not improve base stacking: computational analysis of base stacking in a d(CATGGGCCCATG)(2) B<-->A intermediate crystal structure

The crystal structure of d(CATGGGCCCATG)(2) shows unique stacking patterns of a stable B<-->A-DNA intermediate. We evaluated intrinsic base stacking energies in this crystal structure using an ab initio quantum mechanical method. We found that all crystal base pair steps have stacking energies close to their values in the standard and crystal B-DNA geometries. Thus, naturally occurring stacking geometries were essentially isoenergetic while individual base pair steps differed substantially in the balance of intra-strand and inter-strand stacking terms. Also, relative dispersion, electrostatic and polarization contributions to the stability of different base pair steps were very sensitive to base composition and sequence context. A large stacking flexibility is most apparent for the CpA step, while the GpG step is characterized by weak intra-strand stacking. Hydration effects were estimated using the Langevin dipoles solvation model. These calculations showed that an aqueous environment efficiently compensates for electrostatic stacking contributions. Finally, we have carried out explicit solvent molecular dynamics simulation of the d(CATGGGCCCATG)(2) duplex in water. Here the DNA conformation did not retain the initial crystal geometry, but moved from the B<-->A intermediate towards the B-DNA structure. The base stacking energy improved in the course of this simulation. Our findings indicate that intrinsic base stacking interactions are not sufficient to stabilize the local conformational variations in crystals.

On the truncation of long-range electrostatic interactions in DNA

Long-range interactions are known to play an important role in highly polar biomolecules like DNA. In molecular dynamics simulations of nucleic acids and proteins, an accurate treatment of the long-range interactions are crucial for achieving stable nanosecond trajectories. In this report, we evaluate the structural and dynamic effects on a highly charged oligonucleotide in aqueous solution from different long-range truncation methods. Two group-based truncation methods, one with a switching function and one with a force-switching function were found to fail to give accurate stable trajectories close to the crystal structure. For these group-based truncation methods, large root mean square (rms) deviations from the initial structure were obtained and severe distortions of the oligonucleotide were observed. Another group-based truncation scheme, which used an abrupt truncation at 8. 0 A or at 12.0 A was also investigated. For the short cutoff distance, the conformations deviated far away from the initial structure and were significantly distorted. However, for the longer cutoff, where all necessary electrostatic interactions were included, the trajectory was quite stable. For the particle mesh Ewald (PME) truncation method, a stable DNA simulation with a heavy atom rms deviation of 1.5 A was obtained. The atom-based truncation methods also resulted in stable trajectories, according to the rms deviation from the initial B-DNA structure, of between 1.5 and 1.7 A for the heavy atoms. In these stable simulations, the heavy atom rms deviations were approximately 0.6-1.0 A lower for the bases than for the backbone. An increase of the cutoff radius from 8 to 12 A decreased the rms deviation by approximately 0.2 A for the atom-based truncation method with a force-shifting function, but increased the computational time by a factor of 2. Increasing the cutoff from 12 to 18 A for the atom-based truncation method with a force-shifting function requires 2-3 times more computational time, but did not significantly change the rms deviation. Similar rms deviations from the initial structure were found for the atom-based method with a force-shifting function and for the PME method. The computational cost was longer for the PME method with a cutoff of 12. 0 A for the direct space nonbonded calculations than for the atom-based truncation method with a force-shifting function and a cutoff of 12.0 A. If a nonperiodic boundary, e.g., a spherical boundary, was used, a considerable speedup could be achieved. From the rms fluctuations, the terminal nucleotides and especially the cytidines were found to be more flexible than the nonterminal nucleotides. The B-DNA form of the oligonucleotide was maintained throughout the simulations and is judged to depend on the parameters of the energy function and not on the truncation method used to handle the long-range electrostatic interactions. To perform accurate and stable simulations of highly charged biological macromolecules, we recommend that the atom-based force-shift method or the PME method should be used for the long-range electrostatics interactions.

Effect of coordinated ions on structure and flexiblity of parallel G-quandruplexes: a molecular dynamics study

Single tract guanine residues can associate to form stable parallel quadruplex structures in the presence of certain cations. Nanosecond scale molecular dynamics simulations have been performed on fully solvated fibre model of parallel d(G7) quadruplex structures with Na+ or K+ ions coordinated in the cavity formed by the 06 atoms of the guanine bases. The AMBER 4.1 force field and Particle Mesh Ewald technique for electrostatic interactions have been used in all simulations. These quadruplex structures are stable during the simulation, with the middle four base tetrads showing root mean square deviation values between 0.5 to 0.8 A from the initial structure as well the high resolution crystal structure. Even in the absence of any coordinated ion in the initial structure, the G-quadruplex structure remains intact throughout the simulation. During the 1.1 ns MD simulation, one Na+ counter ion from the solvent as well as several water molecules enter the central cavity to occupy the empty coordination sites within the parallel quadruplex and help stabilize the structure. Hydrogen bonding pattern depends on the nature of the coordinated ion, with the G-tetrad undergoing local structural variation to accommodate cations of different sizes. In the absence of any coordinated ion, due to strong mutual repulsion, 06 atoms within G-tetrad are forced farther apart from each other, which leads to a considerably different hydrogen bonding scheme within the G-tetrads and very favourable interaction energy between the guanine bases constituting a G-tetrad. However, a coordinated ion between G-tetrads provides extra stacking energy for the G-tetrads and makes the quadruplex structure more rigid. Na+ ions, within the quadruplex cavity, are more mobile than coordinated K+ ions. A number of hydrogen bonded water molecules are observed within the grooves of all quadruplex structures.

Effect of neighboring bases on base-pair stacking orientation: a molecular dynamics study

It is generally believed that base-pair stacking interaction in DNA double helix is one of the strongest interactions that governs sequence directed structural variability. However, X-ray crystal structures of some base-paired doublet sequences have been seen to adopt different structures when flanked by different base-pairs. DNA crystal database, however, is still too small to make good statistical inference about effect of such flanking residues. Influence of neighboring residue on the local helical geometry of a base-paired doublet in B-DNA has been investigated here using molecular dynamics simulation. We have generated ensembles of structures for d(CA).d(TG) and d(AA).d(TT) base-paired doublets located at the centers of d(CGCGCAAAGCG).d(CGCTTTGCGCG) and d(CGCGAAAACGCG).d(CGCGTTTTCGCG) sequences along with their analogs by varying the bases either at 5'- or 3'- position to the central doublet. Comparison of base paired doublet parameters for the ensembles of structures show that stacking geometry of d(CA).d(TG) doublet depends on some of the flanking base-pairs. On the other hand d(AA).d(TT) doublet remains nearly unperturbed when the flanking residues are altered.

Water and ion binding around RNA and DNA (C,G) oligomers

The dynamics, hydration, and ion-binding features of two duplexes, the A(r(CG)(12)) and the B(d(CG)(12)), in a neutralizing aqueous environment with 0.25 M added KCl have been investigated by molecular dynamics (MD) simulations. The regular repeats of the same C=G base-pair motif have been exploited as a statistical alternative to long MD simulations in order to extend the sampling of the conformational space. The trajectories demonstrate the larger flexibility of DNA compared to RNA helices. This flexibility results in less well defined hydration patterns around the DNA than around the RNA backbone atoms. Yet, 22 hydration sites are clearly characterized around both nucleic acid structures. With additional results from MD simulations, the following hydration scale for C=G pairs can be deduced: A-DNA<RNA (+3 H(2)O) and B-DNA<RNA (+2 H(2)O). The calculated residence times of water molecules in the first hydration shell of the helices range from 0.5 to 1 ns, in good agreement with available experimental data. Such water molecules are essentially found in the vicinity of the phosphate groups and in the DNA minor groove. The calculated number of ions that break into the first hydration shell of the nucleic acids is close to 0.5 per base-pair for both RNA and DNA. These ions form contacts essentially with the oxygen atoms of the phosphate groups and with the guanine N7 and O6 atoms; they display residence times in the deep/major groove approaching 500 ps. Further, a significant sequence-dependent effect on ion binding has been noted. Despite slight structural differences, K(+) binds essentially to GpC and not to CpG steps. These results may be of importance for understanding various sequence-dependent binding affinities. Additionally, the data help to rationalize the experimentally observed differences in gel electrophoretic mobility between RNA and DNA as due to the difference in hydration (two water molecules in favor of RNA) rather than to strong ion-binding features, which are largely similar for both nucleic acid structures.

Sequence-dependent elastic properties of DNA

Harmonic elastic constants of 3-11 bp duplex DNA fragments were evaluated using four 5 ns unrestrained molecular dynamics simulation trajectories of 17 bp duplexes with explicit inclusion of solvent and counterions. All simulations were carried out with the Cornell et al. force-field and particle mesh Ewald method for long-range electrostatic interactions. The elastic constants including anisotropic bending and all coupling terms were derived by analyzing the correlations of fluctuations of structural properties along the trajectories. The following sequences have been considered: homopolymer d(ApA)(n) and d(GpG)(n), and alternating d(GPC)(n) and d(APT)(n). The calculated values of elastic constants are in very good overall agreement with experimental values for random sequences. The atomic-resolution molecular dynamics approach, however, reveals a pronounced sequence-dependence of the stretching and torsional rigidity of DNA, while sequence-dependence of the bending rigidity is smaller for the sequences considered. The earlier predicted twist-bend coupling emerged as the most important cross-term for fragments shorter than one helical turn. The calculated hydrodynamic relaxation times suggest that damping of bending motions may play a role in molecular dynamics simulations of long DNA fragments. A comparison of elasticity calculations using global and local helicoidal analyses is reported. The calculations reveal the importance of the fragment length definition. The present work shows that large-scale molecular dynamics simulations represent a unique source of data to study various aspects of DNA elasticity including its sequence-dependence.

Modeling DNA deformations

Recent developments have been made in modeling double-helical DNA at four levels of three-dimensional structure: the all-atom level, whereby an oligonucleotide duplex is surrounded by a shroud of solvent molecules; the base-pair level, with explicit backbone atoms; the mesoscopic level, that is, a few hundred base pairs, with the local duplex conformation described by knowledge-based harmonic energy functions; and the scale of several thousand nucleotides, with the duplex described as an ideal elastic rod. Predictions of the sequence-dependent bending and twisting of the double helix, as well as solvent- and force-induced B-->A and over-stretching conformational transitions, are compared with experimental data. These subtle conformational changes are critical to the functioning of the double helix, including its packaging in the close confines of the cell, the mutual fit of DNA and protein in nucleoprotein complexes, and the effective recognition of base pairs in recombination and transcription.

Conformational deformability of RNA: a harmonic mode analysis

The harmonic mode analysis method was used to characterize the conformational deformability of regular Watson-Crick paired, mismatch- and bulge-containing RNA. Good agreement between atomic Debye-Waller factors derived from x-ray crystallography of a regular RNA oligonucleotide and calculated atomic fluctuations was obtained. Calculated helical coordinate fluctuations showed a small sequence dependence of up to approximately 30-50%. A negative correlation between motions at a given base pair step and neighboring steps was found for most helical coordinates. Only very few calculated modes contribute significantly to global motions such as bending, twisting, and stretching of the RNA molecules. With respect to a local helical description of the RNA helix our calculations suggest that RNA bending is mostly due to a periodic change in the base pair step descriptors slide and roll. The presence of single guanine:uridine or guanine:adenine mismatches had little influence on the calculated RNA flexibility. In contrast, for tandem guanine:adenine base pairs the harmonic mode approach predicts a significantly reduced conformational flexibility in the case of a sheared arrangement and slightly enhanced flexibility for a face-to-face (imino proton) pairing relative to regular RNA. The presence of a single extra adenine bulge nucleotide stacked between flanking sequences resulted in an increased local atomic mobility around the bulge site (approximately 40%) and a slightly enhanced global bending flexibility. For an adenine bulge nucleotide in a looped-out conformation a strongly enhanced bulge nucleotide mobility but no increased bending flexibility compared to regular RNA was found.

Nucleic acids: theory and computer simulation, Y2K

Molecular dynamics simulations on DNA and RNA that include solvent are now being performed under realistic environmental conditions of water activity and salt. Improvements to force-fields and treatments of long-range interactions have significantly increased the reliability of simulations. New studies of sequence effects, axis bending, solvation and conformational transitions have appeared.

From atomic to mesoscopic descriptions of the internal dynamics of DNA

An analysis of four 1-ns molecular dynamics trajectories for two different 15-bp oligonucleotides is presented. Our aim is to show which groups of atoms can be treated as rigid bodies within a bead representation of DNA, independently of the base sequence and for any conformations belonging to the A/B family. Five models with moderate intragroup deformations are proposed in which the groups are formed of atoms belonging to a single nucleotide or to a complementary nucleotide pair. The influence of group deformation in two of these models is studied using canonical correlation analysis, and it is shown that the internal DNA dynamics is indeed dominated by the rigid motion of the defined atom groups. Finally, using one of the models within a bead representation of duplex DNA makes it possible to obtain stretching, torsional, and bending rigidities in reasonable agreement with experiment but points to strongly correlated stretching motions.

Sodium and chlorine ions as part of the DNA solvation shell

The distribution of sodium and chlorine ions around DNA is presented from two molecular dynamics simulations of the DNA fragment d(C(5)T(5)). (A(5)G(5)) in explicit solvent with 0.8 M additional NaCl salt. One simulation was carried out for 10 ns with the CHARMM force field that keeps the DNA structure close to A-DNA, the other for 12 ns with the AMBER force field that preferentially stabilizes B-DNA conformations (, Biophys. J. 75:134-149). From radial distributions of sodium and chlorine ions a primary ion shell is defined. The ion counts and residence times of ions within this shell are compared between conformations and with experiment. Ordered sodium ion sites were found in minor and major grooves around both A and B-DNA conformations. Changes in the surrounding hydration structure are analyzed and implications for the stabilization of A-DNA and B-DNA conformations are discussed.