Molecular dynamics studies of DNA

Coarse-grained model simulation-guided localized DNA signal amplification probe for miRNA detection

DNA-based enzyme-free signal amplification strategies are widely employed to detect biomarkers in low abundance. To enhance signal amplification, localized DNA reaction units which increases molecular collision probability is commonly utilized. However, the current understanding of the structure-function relationships in localized DNA signal amplification probes is limited, leading to unsatisfied performance. In this study, we introduced a coarse-grained molecular model to simulate the dynamic behavior of two DNA reaction units within a DNA enzyme-free signal amplification circuit called Localized Catalytic Hairpin Assembly (LCHA). We investigated the impact of localized distance and flexibility on reaction performance. The most efficient LCHA probe guided by simulation exhibits sensitivity 28 times greater that of free CHA, with a detection limit of miR-21 reaching 16 pM, while the least effective LCHA probe demonstrated a modest improvement of only 7 times. We successfully employed the optimized probe to differentiate cancer cells from normal cells based on their miR-21 expression levels, showcasing its quantification ability. By elucidating the mechanistic insights and structure-function relationship in our work, we aim to contribute valuable information that can save users' time and reduce costs when designing localized DNA probes. With a comprehensive understanding of how the localization affects probe performance, researchers can now make more informed and efficient decisions during the design process. This work would find broad applications of DNA nanotechnology in biosensing, biocomputing, and bionic robots.

Molecular mechanism of ATP and RNA binding to Zika virus NS3 helicase and identification of repurposed drugs using molecular dynamics simulations

Congenital Zika virus syndrome has caused a public health emergency of international concern. So far, there are no drugs available to prevent or treat the infection caused by Zika virus. The Zika virus NS3 helicase is a potential protein target for drug discovery due to its vital role in viral genome replication. NS3 helicase unwinds the viral RNA to enable the reproduction of the viral genome by the NS5 protein. NS3 helicase has two crucial binding sites; the ATP binding site and the RNA binding site. We used molecular docking and molecular dynamics (MD) simulations to study the structural behavior of Zika virus NS3 helicase in its apo form and in the presence of ATP, single-stranded RNA, and both ATP-RNA to understand their potential implications in NS3 helicase activity. Further, we have carried out virtual screening of FDA approved drugs, followed by molecular docking to identify the ATP-competitive hit molecules as probable Zika virus NS3 helicase inhibitors. The MD simulations trajectories were analyzed using normal mode analysis and principal component analysis that reveals fluctuations in the R-loop. These findings aid in understanding the molecular mechanisms of the simultaneous binding of ATP and RNA, and guide the design and discovery of new inhibitors of the Zika virus NS3 helicase as a promising drug target to treat the Zika virus infection. Communicated by Ramaswamy H. Sarma.

TacoxDNA: A user-friendly web server for simulations of complex DNA structures, from single strands to origami

Simulations of nucleic acids at different levels of structural details are increasingly used to complement and interpret experiments in different fields, from biophysics to medicine and materials science. However, the various structural models currently available for DNA and RNA and their accompanying suites of computational tools can be very rarely used in a synergistic fashion. The tacoxDNA webserver and standalone software package presented here are a step toward a long-sought interoperability of nucleic acids models. The webserver offers a simple interface for converting various common input formats of DNA structures and setting up molecular dynamics (MD) simulations. Users can, for instance, design DNA rings with different topologies, such as knots, with and without supercoiling, by simply providing an XYZ coordinate file of the DNA centre-line. More complex DNA geometries, as designable in the cadnano, CanDo and Tiamat tools, can also be converted to all-atom or oxDNA representations, which can then be used to run MD simulations. Though the latter are currently geared toward the native and LAMMPS oxDNA representations, the open-source package is designed to be further expandable. TacoxDNA is available at http://tacoxdna.sissa.it. © 2019 Wiley Periodicals, Inc.

Molecular Dynamics Simulations of Dna and Its Complexes

This article describes how classical molecular simulation methods are being used to gain a molecular-level understanding of the interaction mechanisms responsible for DNA–ligand recognition, and that govern the response of DNA to ligand binding. Case studies using a variety of different ligands—including small pharmaceutical drugs, proteins and lipids—are used to illustrate the power of modern molecular dynamics simulation methods for understanding how we may control the function and structure of DNA.

NAMD: a Parallel, Object-Oriented Molecular Dynamics Program

NAMD is a molecular dynamics program designed for high performance simulations of large biomolecular systems on parallel computers. An object-oriented design imple mented using C++ facilitates the incorporation of new algorithms into the program. NAMD uses spatial decom position coupled with a multithreaded, message-driven design, which is shown to scale efficiently to multiple processors. Also, NAMD incorporates the distributed par allel multipole tree algorithm for full electrostatic force evaluation in O( N) time. NAMD can be connected via a communication system to a molecular graphics program in order to provide an interactive modeling tool for viewing and modifying a running simulation. The application of NAMD to a protein-water system of 32,867 atoms illus trates the performance of NAMD.

Base damage, local sequence context and TP53 mutation hotspots: a molecular dynamics study of benzo[a]pyrene induced DNA distortion and mutability

The mutational pattern for the TP53 tumour suppressor gene in lung tumours differs to other cancer types by having a higher frequency of G:C>T:A transversions. The aetiology of this differing mutation pattern is still unknown. Benzo[a]pyrene,diol epoxide (BPDE) is a potent cigarette smoke carcinogen that forms guanine adducts at TP53 CpG mutation hotspot sites including codons 157, 158, 245, 248 and 273. We performed molecular modelling of BPDE-adducted TP53 duplex sequences to determine the degree of local distortion caused by adducts which could influence the ability of nucleotide excision repair. We show that BPDE adducted codon 157 has greater structural distortion than other TP53 G:C>T:A hotspot sites and that sequence context more distal to adjacent bases must influence local distortion. Using TP53 trinucleotide mutation signatures for lung cancer in smokers and non-smokers we further show that codons 157 and 273 have the highest mutation probability in smokers. Combining this information with adduct structural data we predict that G:C>T:A mutations at codon 157 in lung tumours of smokers are predominantly caused by BPDE. Our results provide insight into how different DNA sequence contexts show variability in DNA distortion at mutagen adduct sites that could compromise DNA repair at well characterized cancer related mutation hotspots.

Close encounters with DNA

Over the past ten years, the all-atom molecular dynamics method has grown in the scale of both systems and processes amenable to it and in its ability to make quantitative predictions about the behavior of experimental systems. The field of computational DNA research is no exception, witnessing a dramatic increase in the size of systems simulated with atomic resolution, the duration of individual simulations and the realism of the simulation outcomes. In this topical review, we describe the hallmark physical properties of DNA from the perspective of all-atom simulations. We demonstrate the amazing ability of such simulations to reveal the microscopic physical origins of experimentally observed phenomena. We also discuss the frustrating limitations associated with imperfections of present atomic force fields and inadequate sampling. The review is focused on the following four physical properties of DNA: effective electric charge, response to an external mechanical force, interaction with other DNA molecules and behavior in an external electric field.

HYDRATE MODEL OF THE EQUILIBRIUM DNA–WATER SYSTEMS

A hydrate model for the DNA– H2O system is proposed. For the notions developed, available data on a tendency of H2O molecules to form structured hydrates containing atomic groups, molecules, or atoms housed within the structural cavities formed by H-bonded H2O molecules are used. It is shown that the large and small cavities of the hydrate structure II are in close geometric agreement with N-bases and deoxyribose and with other atomic groups of DNA molecules, respectively. On the basis of the model proposed and with the Watson–Crick base-pairing scheme, the number of pairs of N-bases per helix turn (11.25), the density (1.161 g/cm3), and the helix step (0.567 nm) for the quasi-equilibrium DNA (RNA)–water system are computed, the last value nearly coinciding with the well-known alpha-helix step in the protein secondary structure (5.44 nm). Assumptions on the phenomenology of some stages of mitosis are presented.

Analysis of a DNA simulation model through hairpin melting experiments

We compare the predictions of a two-bead Brownian dynamics simulation model to melting experiments of DNA hairpins with complementary AT or GC stems and noninteracting loops in buffer A. This system emphasizes the role of stacking and hydrogen bonding energies, which are characteristics of DNA, rather than backbone bending, stiffness, and excluded volume interactions, which are generic characteristics of semiflexible polymers. By comparing high throughput data on the open-close transition of various DNA hairpins to the corresponding simulation data, we (1) establish a suitable metric to compare the simulations to experiments, (2) find a conversion between the simulation and experimental temperatures, and (3) point out several limitations of the model, including the lack of G-quartets and cross stacking effects. Our approach and experimental data can be used to validate similar coarse-grained simulation models.

Molecular simulations of carbohydrates and protein-carbohydrate interactions: motivation, issues and prospects

The characterization of the 3D structure of oligosaccharides, their conjugates and analogs is particularly challenging for traditional experimental methods. Molecular simulation methods provide a basis for interpreting sparse experimental data and for independently predicting conformational and dynamic properties of glycans. Here, we summarize and analyze the issues associated with modeling carbohydrates, with a detailed discussion of four of the most recently developed carbohydrate force fields, reviewed in terms of applicability to natural glycans, carbohydrate-protein complexes and the emerging area of glycomimetic drugs. In addition, we discuss prospectives and new applications of carbohydrate modeling in drug discovery.

Structural properties of polymeric DNA from molecular dynamics simulations

Most of the reported DNA structural studies are based on oligonucleotide structures, which have artifacts due to unstable terminal base pairs (bps). We have carried out molecular dynamics simulation of DNA oligonucleotides in such a manner that gives rise to properties of polymeric DNA of infinite length. Molecular dynamics simulation studies of six homo- and heteropolymeric DNA sequences are reported here to understand structural features of all ten unique dinucleotide sequences. We observe that each of these dinucleotide sequences has unique features in agreement with Calladine's rule [C. R. Calladine, J. Mol. Biol. 161, 343 (1982)]. We noticed significant structural alternation between B(I) and B(II) forms for d(CA).d(TG) dinucleotide, where one of the strands showed frequent transitions between usual and unusual epsilon and zeta torsion angles associated with bp stacking geometry. In terms of the calculated bending rigidity and persistence length, pyrimidine-purine bp steps, namely, d(TA).d(TA), d(CA).d(TG), and d(CG).d(CG) are the most flexible dinucleotide bp steps. We estimated the major groove widths from our simulations. We did not observe much variation in major and minor groove widths depending on the base sequence. However, the distribution of water molecules in the minor groove shows sensitivity to the DNA sequence.

Probing the interactions of the solvated electron with DNA by molecular dynamics simulations: bromodeoxyuridine substituted DNA

Solvated electrons ((e-)(aq)) are produced during water radiolysis and can interact with biological substrates, including DNA. To augment DNA damage, radiosensitizers such as bromo-deoxyuridine (BUdR), often referred to as an "electron affinic radiosensitizer", are incorporated in place of isosteric thymidine. However, little is known about the primary interactions of (e-)(aq) with DNA. In the present study we addressed this problem by applying molecular modeling and molecular dynamics (MD) simulations to a system of normal (BUdR.A)-DNA and a hydrated electron, where the excess electron was modeled as a localized (e-)(H2O)6 anionic cluster. Our goals were to evaluate the suitability of the MD simulations for this application; to characterize the motion of (e-)(aq) around DNA (e.g., diffusion coefficients); to identify and describe configurational states of close (e-)(aq) localization to DNA; and to evaluate the structural dynamics of DNA in the presence of (e-)(aq). The results indicate that (e-)(aq) has distinct space-preferences for forming close contacts with DNA and is more likely to interact directly with nucleotides other than BUdR. Several classes of DNA - (e-)(aq) contact sites, all within the major groove, were distinguished depending on the structure of the intermediate water layer H-bonding pattern (or its absence, i.e., a direct H-bonding of (e-)(aq) with DNA bases). Large-scale structural perturbations were identified during and after the (e-)(aq) approached the DNA from the major groove side, coupled with deeper penetration of sodium counterions in the minor groove.

Probing potential binding modes of the p53 tetramer to DNA based on the symmetries encoded in p53 response elements

Symmetries in the p53 response-element (p53RE) encode binding modes for p53 tetramer to recognize DNA. We investigated the molecular mechanisms and biological implications of the possible binding modes. The probabilities evaluated with molecular dynamics simulations and DNA sequence analyses were found to be correlated, indicating that p53 tetramer models studied here are able to read DNA sequence information. The traditionally believed mode with four p53 monomers binding at all four DNA quarter-sites does not cause linear DNA to bend. Alternatively, p53 tetramer can use only two monomers to recognize DNA sequence and induce DNA bending. With an arrangement of dimer of AB dimer observed in p53 trimer-DNA complex crystal, p53 can recognize supercoiled DNA sequence-specifically by binding to quarter-sites one and four (H14 mode) and recognize Holliday junction geometry-specifically. Examining R273H mutation and p53-DNA interactions, we found that at least three R273H monomers are needed to disable the p53 tetramer, consistent with experiments. But just one R273H monomer may greatly shift the binding mode probabilities. Our work suggests that p53 needs balanced binding modes to maintain genome stability. Inverse repeat p53REs favor the H14 mode and direct repeat p53REs may have high possibilities of other modes.

Theoretical study of large conformational transitions in DNA: the B<-->A conformational change in water and ethanol/water

We explore here the possibility of determining theoretically the free energy change associated with large conformational transitions in DNA, like the solvent-induced B<-->A conformational change. We find that a combination of targeted molecular dynamics (tMD) and the weighted histogram analysis method (WHAM) can be used to trace this transition in both water and ethanol/water mixture. The pathway of the transition in the A-->B direction mirrors the B-->A pathway, and is dominated by two processes that occur somewhat independently: local changes in sugar puckering and global rearrangements (particularly twist and roll) in the structure. The B-->A transition is found to be a quasi-harmonic process, which follows closely the first spontaneous deformation mode of B-DNA, showing that a physiologically-relevant deformation is in coded in the flexibility pattern of DNA.

RNA unrestrained molecular dynamics ensemble improves agreement with experimental NMR data compared to single static structure: a test case

Nuclear magnetic resonance (NMR) provides structural and dynamic information reflecting an average, often non-linear, of multiple solution-state conformations. Therefore, a single optimized structure derived from NMR refinement may be misleading if the NMR data actually result from averaging of distinct conformers. It is hypothesized that a conformational ensemble generated by a valid molecular dynamics (MD) simulation should be able to improve agreement with the NMR data set compared with the single optimized starting structure. Using a model system consisting of two sequence-related self-complementary ribonucleotide octamers for which NMR data was available, 0.3 ns particle mesh Ewald MD simulations were performed in the AMBER force field in the presence of explicit water and counterions. Agreement of the averaged properties of the molecular dynamics ensembles with NMR data such as homonuclear proton nuclear Overhauser effect (NOE)-based distance constraints, homonuclear proton and heteronuclear (1)H-(31)P coupling constant (J) data, and qualitative NMR information on hydrogen bond occupancy, was systematically assessed. Despite the short length of the simulation, the ensemble generated from it agreed with the NMR experimental constraints more completely than the single optimized NMR structure. This suggests that short unrestrained MD simulations may be of utility in interpreting NMR results. As expected, a 0.5 ns simulation utilizing a distance dependent dielectric did not improve agreement with the NMR data, consistent with its inferior exploration of conformational space as assessed by 2-D RMSD plots. Thus, ability to rapidly improve agreement with NMR constraints may be a sensitive diagnostic of the MD methods themselves.

Effect of DNA structural flexibility on promoter strength--molecular dynamics studies of E. coli promoter sequences

To study possible correlations between promoter activity and the structural flexibility of the DNA helix, we have carried out unrestrained molecular dynamics simulations of the -10 consensus region sequence and five variants forming the -10 region of various Escherichia coli promoter sequences. Analyses of the trajectories obtained from the simulations show that the consensus sequence has a pattern of two structurally flexible nucleotide steps sandwiched between two stiff steps. In the other sequences, this pattern varies in consonance with the change in the sequence. The variations in the patterns show correlation with the promoter strength.

M.TaqI facilitates the base flipping via an unusual DNA backbone conformation

MD simulations have been carried out to understand the dynamical behavior of the DNA substrate of the Thermus aquaticus DNA methyltransferase (M.TaqI) in the methylation process at N6 of adenine. As starting structures, an x-ray structure of M.TaqI in complex with DNA and cofactor analogue (PDB code: 1G 38) and free decamer d(GTTCGATGTC)(2) were taken. The x-ray structure shows two consecutive BII substates that are not observed in the free decamer. These consecutive BII substates are also observed during our simulation. Additionally, their facing backbones adopt the same conformations. These double facing BII substates are stable during the last 9 ns of the trajectories and result in a stretched DNA structure. On the other hand, protein-DNA contacts on 5' and 3' phosphodiester groups of the partner thymine of flipped adenine have changed. The sugar and phosphate parts of thymine have moved further into the empty space left by the flipping base without the influence of protein. Furthermore, readily high populated BII substates at the GpA step of palindromic tetrad TCGA rather than CpG step are observed in the free decamer. On the contrary, the BI substate at the GpA step is observed on the flipped adenine strand. A restrained MD simulation, reproducing the BI/BII pattern in the complex, demonstrated the influence of the unusual backbone conformation on the dynamical behavior of the target base. This finding along with the increased nearby interstrand phosphate distance is supportive to the N6-methylation mechanism.

Microscopic Kinetics of DNA Translocation through synthetic nanopores

We have previously demonstrated that a nanometer-diameter pore in a nanometer-thick metal-oxide-semiconductor-compatible membrane can be used as a molecular sensor for detecting DNA. The prospects for using this type of device for sequencing DNA are avidly being pursued. The key attribute of the sensor is the electric field-induced (voltage-driven) translocation of the DNA molecule in an electrolytic solution across the membrane through the nanopore. To complement ongoing experimental studies developing such pores and measuring signals in response to the presence of DNA, we conducted molecular dynamics simulations of DNA translocation through the nanopore. A typical simulated system included a patch of a silicon nitride membrane dividing water solution of potassium chloride into two compartments connected by the nanopore. External electrical fields induced capturing of the DNA molecules by the pore from the solution and subsequent translocation. Molecular dynamics simulations suggest that 20-basepair segments of double-stranded DNA can transit a nanopore of 2.2 x 2.6 nm(2) cross section in a few microseconds at typical electrical fields. Hydrophobic interactions between DNA bases and the pore surface can slow down translocation of single-stranded DNA and might favor unzipping of double-stranded DNA inside the pore. DNA occluding the pore mouth blocks the electrolytic current through the pore; these current blockades were found to have the same magnitude as the blockade observed when DNA transits the pore. The feasibility of using molecular dynamics simulations to relate the level of the blocked ionic current to the sequence of DNA was investigated.

Water-mediated contacts in thetrp-repressor operator complex recognition process

Water-mediated contacts are known as an important recognition tool in trp-repressor operator systems. One of these contacts involves two conserved base pairs (G(6).C(-6) and A(5). T(-5)) and three amino acids (Lys 72, Ile 79, and Ala 80). To investigate the nature of these contacts, we analyzed the X-ray structure (PDB code: 1TRO) of the trp-repressor operator complex by means of molecular dynamics simulations. This X-ray structure contains two dimers that exhibit structural differences. From these two different starting structures, two 10 ns molecular dynamics simulations have been performed. Both of our simulations show an increase of water molecules in the major groove at one side of the dimer, while the other side remains unchanged compared to the X-ray structure. Though the maximum residence time of the concerned water molecules decreases with an increase of solvent at the interface, these water molecules continue to play an important role in mediating DNA-protein contacts. This is shown by new stable amino acids-DNA distances and a long water residence time compared to free DNA simulation. To maintain stability of the new contacts, the preferential water binding site on O6(G6) is extended. This extension agrees with mutation experiment data on A5 and G6, which shows different relative affinity due to mutation on these bases [A. Joachimiak, T. E. Haran, P. B. Sigler, EMBO Journal 1994, Vol. 13, No. (2) pp. 367-372]. Due to the rearrangements in the system, the phosphate of the base G6 is able to interconvert to the B(II) substate, which is not observed on the other half side of the complex. The decrease of the number of hydrogen bonds between protein and DNA backbone could be the initial step of the dissociation process of the complex, or in other words an intermediate complex conformation of the association process. Thus, we surmise that these features show the importance of water-mediated contacts in the trp-repressor operator recognition process.

Molecular dynamics simulation of the RNA complex of a double-stranded RNA-binding domain reveals dynamic features of the intermolecular interface and its hydration

The interaction between double-stranded RNA (dsRNA) and the third double-stranded domain (dsRBD) from Drosophila Staufen protein represents a paradigm to understand how the dsRBD protein family, one of the most common RNA-binding protein units, binds dsRNA. The nuclear magnetic resonance (NMR) structure of this complex and the x-ray structure of another family member revealed the stereochemical basis for recognition, but also raised new questions. Although the crystallographic studies revealed a highly ordered interface containing numerous water-mediated contacts, NMR suggested extensive residual motion at the interface. To address how interfacial motion contributes to molecular recognition in the dsRBD-dsRNA system, we conducted a 2-ns molecular dynamics simulation of the complex derived from Staufen protein and of the separate protein and RNA components. The results support the observation that a high degree of conformational flexibility is retained upon complex formation and that this involves interfacial residues that are critical for dsRBD-dsRNA binding. The structural origin of this residual flexibility is revealed by the analysis of the trajectory of motion. Individual basic side chains switch continuously from one RNA polar group to another with a residence time seldom exceeding 100 ps, while retaining favorable interaction with RNA throughout much of the simulation. Short-lived water molecules mediate some of these interactions for a large fraction of the trajectory studied here. This result indicates that water molecules are not statically associated with the interface, but continuously exchange with the bulk solvent on a 1-10-ps time scale. This work provides new insight into dsRBD-dsRNA recognition and builds upon a growing body of evidence, suggesting that short-lived dynamic interactions play important roles in protein-nucleic acid interactions.

PvuII-endonuclease induces structural alterations at the scissile phosphate group of its cognate DNA

We investigated the PvuII endonuclease with its cognate DNA by means of molecular dynamics simulations. Comparing the complexed DNA with a reference simulation of free DNA, we saw structural changes at the scissile phosphodiester bond. At this GpC step, the enzyme induces the highest twist and axial rise, inclination is increased and the minor groove widened. The distance between the scissile phosphate group and the phosphate group of the following thymine base is shortened significantly, indicating a substrate-assisted catalysis. A feasible reason for this vicinity is the catalytically important amino acid residue lysine 70, which bridges the free oxygen atoms of the successive phosphate groups. Due to this geometry, a compact reaction pocket is formed where a water molecule can be held, thus bringing the reaction partners for hydrolysis into contact. The O1-P-O2 angle of the scissile nucleotide is decreased, probably due to a complexation of the negative oxygen atoms through protein and solvent contacts.

Influence of netropsin's charges on the minor groove width of d(CGCGAATTCGCG)2

The exact understanding of the interaction of minor groove binding drugs with DNA is of interest due to their importance as transcription controlling drugs. In this study we performed four molecular dynamics simulations, one of the uncomplexed d(CGCGAATTCGCG)(2) dodecamer and three simulations of the DNA complexed with the minor groove binder netropsin. The charged guanidinium and amidinium ends of the small ligand were in one simulation formally uncharged, in the second one normally charged, and in the third simulation we doubled the charges of the two ends. So we are able to filter out the influence the charges exert on the DNA structure. The positive charges reduce the width of the minor groove showing that charges are able to modify the groove width by charge neutralization of the negative phosphate groups. The quality of the used force field was successfully tested by comparing the results of the uncomplexed dodecamer with already reported NMR and x-ray studies. Thus our simulations should be able to describe the minor groove width of DNA in a correct manner underlying the validity of the results.

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.

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.

Structure and hydration of BamHI DNA recognition site: a molecular dynamics investigation

The results of a 3-ns molecular dynamics simulation of the dodecamer duplex d(TATGGATCCATA)(2) recognized by the BamHI endonuclease are presented here. The DNA has been simulated as a flexible molecule using an AMBER force field and the Ewald summation method, which eliminates the undesired effects of truncation and permits evaluation of the full effects of electrostatic forces. The starting B conformation evolves toward a configuration quite close to that observed through x-ray diffraction in its complex with BamHI. This configuration is fairly stable and the Watson-Crick hydrogen bonds are well maintained over the simulation trajectory. Hydration analysis indicates a preferential hydration for the phosphate rather than for the ester oxygens. Hydration shells in both the major and minor groove were observed. In both grooves the C-G pairs were found to be more hydrated than A-T pairs. The "spine of hydration" in the minor groove was clear. Water residence times are longer in the minor groove than in the major groove, although relatively short in both cases. No special long values are observed for sites where water molecules were observed by x-ray diffraction, indicating that water molecules having a high probability of being located in a specific site are also fast-exchanging.

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.

Probing the role of structural water in a duplex oligodeoxyribonucleotide containing a water-mimicking base analog

Molecular dynamics simulations were performed on models of the dodecamer DNA double-stranded segment, [d(CGCGAATTCGCG)](2), in which each of the adenine residues, individually or jointly, was replaced by the water-mimicking analog 2'-deoxy-7-(hydroxy-methyl)-7-deazaadenosine (hm(7)c(7)dA) [Rockhill, J.K., Wilson,S.R. and Gumport,R.I. (1996) J. Am. Chem. Soc.,118, 10065-10068]. The simulations, when compared with those of the dodecamer itself, show that incorporation of the analog affects neither the overall DNA structure nor its hydrogen-bonding and stacking interactions when it replaces a single individual base. Furthermore, the water molecules near the bases in the singly-substituted oligonucleotides are similarly unaffected. Double substitutions lead to differences in all the aforementioned parameters with respect to the reference sequence. The results suggest that the analog provides a good mimic of specific 'ordered' water molecules observed in contact with DNA itself and at the interface between protein and DNA in specific complexes.

Unrestrained stochastic dynamics simulations of the UUCG tetraloop using an implicit solvation model

Three unrestrained stochastic dynamics simulations have been carried out on the RNA hairpin GGAC[UUCG] GUCC, using the AMBER94 force field (Cornell et al., 1995. J. Am. Chem. Soc. 117:5179-5197) in MacroModel 5.5 (Mohamadi et al., 1990. J. Comp. Chem. 11:440-467) and either the GB/SA continuum solvation model (Still et al., 1990. J. Am. Chem. Soc. 112:6127-6129) or a linear distance-dependent dielectric (1/R) treatment. The linear distance-dependent treatment results in severe distortion of the nucleic acid structure, restriction of all hydroxyl dihedrals, and collapse of the counterion atmosphere over the course of a 5-ns simulation. An additional vacuum simulation without counterions shows somewhat improved behavior. In contrast, the two GB/SA simulations (1.149 and 3.060 ns in length) give average structures within 1.2 A of the initial NMR structure and in excellent agreement with results of an earlier explicit solvent simulation (Miller and Kollman, 1997. J. Mol. Biol. 270:436-450). In a 3-ns GB/SA simulation starting with the incorrect UUCG tetraloop structure (Cheong et al., 1990. Nature. 346:680-682), this loop conformation converts to the correct loop geometry (Allain and Varani, 1995. J. Mol. Biol. 250:333-353), suggesting enhanced sampling relative to the previous explicit solvent simulation. Thermodynamic effects of 2'-deoxyribose substitutions of loop nucleotides were experimentally determined and are found to correlate with the fraction of time the ribose 2'-OH is hydrogen bonded and the distribution of the hydroxyl dihedral is observed in the GB/SA simulations. The GB/SA simulations thus appear to faithfully represent structural features of the RNA without the computational expense of explicit solvent.

Molecular dynamic simulations of environment and sequence dependent DNA conformations: the development of the BMS nucleic acid force field and comparison with experimental results

Molecular dynamic (MD) simulations using the BMS nucleic acid force field produce environment and sequence dependent DNA conformations that closely mimic experimentally derived structures. The parameters were initially developed to reproduce the potential energy surface, as defined by quantum mechanics, for a set of small molecules that can be used as the building blocks for nucleic acid macromolecules (dimethyl phosphate, cyclopentane, tetrahydrofuran, etc.). Then the dihedral parameters were fine tuned using a series of condensed phase MD simulations of DNA and RNA (in zero added salt, 4M NaCl, and 75% ethanol solutions). In the tuning process the free energy surface for each dihedral was derived from the MD ensemble and fitted to the conformational distributions and populations observed in 87 A- and B-DNA x-ray and 17 B-DNA NMR structures. Over 41 nanoseconds of MD simulations are presented which demonstrate that the force field is capable of producing stable trajectories, in the correct environments, of A-DNA, double stranded A-form RNA, B-DNA, Z-DNA, and a netropsin-DNA complex that closely reproduce the experimentally determined and/or canonical DNA conformations. Frequently the MD averaged structure is closer to the experimentally determined structure than to the canonical DNA conformation. MD simulations of A- to B- and B- to A-DNA transitions are also shown. A-DNA simulations in a low salt environment cleanly convert into the B-DNA conformation and converge into the RMS space sampled by a low salt simulation of the same sequence starting from B-DNA. In MD simulations using the BMS force field the B-form of d(GGGCCC)2 in a 75% ethanol solution converts into the A-form. Using the same methodology, parameters, and conditions the A-form of d(AAATTT)2 correctly converts into the B-DNA conformation. These studies demonstrate that the force field is capable of reproducing both environment and sequence dependent DNA structures. The 41 nanoseconds (nsec) of MD simulations presented in this paper paint a global picture which suggests that the DNA structures observed in low salt solutions are largely due to the favorable internal energy brought about by the nearly uniform screening of the DNA electrostatics. While the conformations sampled in high salt or mixed solvent environments occur from selective and asymmetric screening of the phosphate groups and DNA grooves, respectively, brought about by sequence induced ion and solvent packing.

The biological implications of damage to DNA incorporating an 8-oxodeoxyguanine:cytosine basepair

DNA damage produced by free radicals is probably the most frequent lesion encountered by cells (Wallace, S.S., Environmental and Molecular Mutagenesis 12, 431-477, 1988 (1)). One of the most common effects is the formation of 7-hydro-8-oxodeoxyguanine due to oxygen radicals interacting with the normal guanine base. Such chemical changes appear to be important in mutagenesis, cancer and aging. We have used computer simulation techniques to model the effect of inclusion of such a modified base within a duplex strand of DNA. We find that such modifications can be stabilized within a normal sequence. The conformation of the modified base relative to the sugar residue depends on many local interactions not accessible to the isolated nucleoside. We have also studied the essential dynamics of both normal and modified sequences and show that there are only subtle changes to the dynamics on inclusion of such a modification.

Hydration of RNA base pairs

The hydration patterns around the RNA Watson-Crick and non-Watson-Crick base pairs in crystals are analyzed and described. The results indicate that (i) the base pair hydration is mostly "in-plane"; (ii) eight hydration sites surround the Watson-Crick G-C and A-U base pairs, with five in the deep and three in the shallow groove, an observation which extends the characteristic isostericity of Watson-Crick pairs; (iii) while the hydration around G-C base pairs is well defined, the hydration around A-U base pairs is more diffuse; (iv) the hydration sites close to the phosphate groups are the best defined and the most recurrent ones; (v) a string of water molecules links the two shallow groove 2'-hydroxyl groups, and (vi) the water molecules fit into notches, the size and accessibility of which are almost as important as the number and strength of the hydrophilic groups lining the cavity. Residence times of water molecules at specific hydration sites, inferred from molecular dynamics simulations, are discussed in the light of present data.

Electrostatic interaction between helical macromolecules in dense aggregates: an impetus for DNA poly- and meso-morphism

DNA exhibits a surprising multiplicity of structures when it is packed into dense aggregates. It undergoes various polymorphous transitions (e.g., from the B to A form) and mesomorphous transformations (from hexagonal to orthorhombic or monoclinic packing, changes in the mutual alignment of nearest neighbors, etc). In this report we show that such phenomena may have their origin in the specific helical symmetry of the charge distribution on DNA surface. Electrostatic interaction between neighboring DNA molecules exhibits strong dependence on the patterns of molecular surface groups and adsorbed counter-ions. As a result, it is affected by such structural parameters as the helical pitch, groove width, the number of base pairs per helical turn, etc. We derive expressions which relate the energy of electrostatic interaction with these parameters and with the packing variables characterizing the axial and azimuthal alignment between neighboring macromolecules. We show, in particular, that the structural changes upon the B-to-A transition reduce the electrostatic energy by approximately kcal/mol per base pair, at a random adsorption of counter ions. Ion binding into the narrow groove weakens or inverts this effect, stabilizing B-DNA, as it is presumably the case in Li+-DNA assemblies. The packing symmetry and molecular alignment in DNA aggregates are shown to be affected by the patterns of ion binding.

Hydration of the phosphate group in double-helical DNA

Water distributions around phosphate groups in 59 B-, A-, and Z-DNA crystal structures were analyzed. It is shown that the waters are concentrated in six hydration sites per phosphate and that the positions and occupancies of these sites are dependent on the conformation and type of nucleotide. The patterns of hydration that are characteristic of the backbone of the three DNA helical types can be attributed in part to the interactions of these hydration sites.

Water molecules in DNA recognition II: a molecular dynamics view of the structure and hydration of the trp operator

The structure and hydration of the DNA duplex d-(AGCGTACTAGTACGCT)2 corresponding to the trp operator fragment used in the crystal structure of the half site complex (PDB entry 1TRR) was studied by a 1.4 ns molecular dynamics simulation in water. The simulation, starting from a B-DNA conformation, used a non-bonded cutoff of 1.4 nm with a reaction field correction and resulted in a stable trajectory. The average DNA conformation obtained was closer to the ones found in the crystal structures of the complexes (PDB entries 1TRO and 1TRR) than to the crystal structure of unbound trp operator (Nucleic Acid Database entry BDJ061). The DNA hydration was characterized in terms of hydrogen bond percentages and corresponding residence times. The residence times of water molecules within 0.35 nm of the DNA non-exchangeable protons were calculated for comparison with NMR measurements of intermolecular water-DNA NOEs and nuclear magnetic relaxation dispersion measurements. No significant difference was found between major and minor groove hydration. The DNA donors and acceptors were hydrogen bonded to water molecules for 77(+/-19)% of the time on average. The average residence time of the hydrogen bonded water molecules was 11(+/-11) ps with a maximum of 223 ps. When all water molecules within NOE distance (0.35 nm) of non-exchangeable protons were considered, the average residence times increased to an average of 100(+/-4) ps and a maximum of 608 ps. These results agree with the experimental NMR results of Sunnerhagen et al. which did not show any evidence for water molecules bound with more than 1 ns residence time on the DNA surface. The exchange of hydration water from the DNA occurred in the major groove primarily through direct exchange with the bulk solvent, while access to and from the minor groove frequently proceeded via pathways involving ribose O3' and O4' and phosphate O2P oxygen atoms. The most common water diffusion pathways in the minor groove were perpendicular to the groove direction. In general, water molecules visited only a limited number of sites in the DNA grooves before exiting. The hydrogen bonding sites, where hydrogen bonds could be formed with donor and acceptor groups of the DNA, were filled with water molecules with an average B-factor value of 0.58 mn2. No special values were observed at any of the sites, where water molecules were observed both in the trp repressor/operator co-crystals and in the crystal structure of unbound DNA.

DNA structure and fluctuations sensed from a 1.1ns molecular dynamics trajectory of a fully charged Zif268-DNA complex in water

Molecular dynamics simulations of the zinc finger domain of protein Zif268, in a complex with a high affinity DNA sequence, yields a globally stable system with small yet significant readjustments with persistence time of the order of 1.1ns. The results confirm the quality of the standard GROMOS87 force field with a corrected solvent-to-solute interaction that does not affect the water-water SPC interactions nor the intra-molecular cohesive forces. Specificity determinants are discussed. The simulations of DNA alone, with the same force field, showed the important role played by the solvent and the symmetry of the counterion distribution. (Tapia & Velázquez, J. Am. Chem. Soc., 119, 5934, 1997) In the present work, this feature was retained when appropriate. The results for root mean square deviations and temperature B-factors illustrate the reliability of this approach. The structure of DNA is held by its interactions with the zinc finger protein. This behavior is not much affected by the slow whithering away of finger-1 from DNA. The factors contributing to the molecular stability found in GROMOS' potential energy function appear to be sufficient to yield stable fluctuation patterns when surrounding medium effects are properly included.

Molecular dynamics simulations of an oligonucleotide duplex with adenine tracts phased by a full helix turn

A theoretical model of a DNA oligonucleotide duplex featuring A-tracts phased by a full helix turn is developed based on molecular dynamics computer simulation. The extent to which this model agrees with relevant experimental data on axis bending and the relationship of A-tracts to bending and other aspects of helix morphology is investigated. Specifically, a series of nanosecond-level molecular dynamics (MD) simulations have been carried out for the 25 bp duplex d(ATAGGCAAAAAATAGGCAAAAATGG) at various concentrations of saline solution. A 30 base-pair sequence composed of three 10 bp repeats of the BamHI recognition sequence ligated together, d(CGGGATCCCG. CGGGATCCCG.CGGGATCCCG), was simulated as a control. The MD was carried out using the AMBER 4.1 suite of programs, and utilized the Cornell et al. force-field with the electrostatic boundary conditions treated by the particle-mesh Ewald summation protocol. The MD results show that at a concentration of 60 mM KCl, 10 mM MgCl2 added salt plus minimal neutralizing cations, the MD model exhibits concerted axis bending to the extent of 15.5 degrees per A-tract. This compares favorably with the bending per turn of 17 to 21 degrees inferred from cyclization experiments. The MD model also exhibits a progressive 5' to 3' narrowing of the minor-groove region of A-tracts, a feature inferred from DNA footprinting experiments. Analysis of the dynamic structure of the MD models shows that the origin of the bending follows a junction-type bending model with an admixture of mixed sequence effects, with A-tracts relatively straight, as in oligonucleotide crystal structures of sequences containing A-tracts. The results are shown to be sensitive to environmental conditions: MD on d(ATAGGCAAAAAATAGGCAAAAATGG) in neutralizing Na+ buffer results in markedly reduced curvature, and the removal of Mg2+ measurably affects bending. Carrying out the simulations at experimental salt conditions appears to be essential to obtain an accurate account of the experimentally observed bending.

The intrinsic curvature of a 51 bp K-DNA fragment of Leishmania tarentolae: a molecular model

DNA intrinsic structure and curvature is a subject of debate because of the importance of these attributes in processes such as DNA packaging, transcription, and gene regulation. X-ray crystallography of DNA single crystals has provided a wealth of information about the local, short range conformational features of DNA. On the other hand, gel electrophoresis analysis of DNA has not only uncovered the macroscopic curvature of DNA but it also provides most of the available data on DNA intrinsic curvature. However, gel electrophoresis can not identify features of DNA structure at the nucleotide or atomic level. In order to address the problem of DNA intrinsic curvature in an attempt to bridge the gap between X-ray crystallography and gel electrophoresis, we use the computational method of molecular dynamics (MD). In this study, we report the results of 2.0 ns MD simulations on a 51 bp fragment of the K-DNA of Leishmania tarentolae containing several A-tracts. The K-DNA double helix is very stable and remains in an intermediate state between the canonical A and B forms of the duplex. The magnitude of global curvature (75 degrees) agrees well with the experimental estimate (72 degrees) available. Analysis of local (every base triplet) and sublocal (every helix turn) curvature shows that the 51 bp K-DNA fragment has curvature features also present in the Wedge, Junction and Calladine's models of DNA intrinsic curvature. We further characterize the flexibility of individual nucleotides in the molecule and find the sugar flexibility within the A-tracts to be strongly correlated with the pattern of A-tract cleavage by the hydroxyl radical. Differential curvature and flexibility at the 5' and 3'junctions between A-tracts and general-sequence DNA are found to modulate the global curvature of the K-DNA fragment.