Molecular dynamics simulations of the 136 unique tetranucleotide sequences of DNA oligonucleotides. II: sequence context effects on the dynamical structures of the 10 unique dinucleotide steps

Molecular dynamics (MD) simulations including water and counterions on B-DNA oligomers containing all 136 unique tetranucleotide basepair steps are reported. The objective is to obtain the calculated dynamical structure for at least two copies of each case, use the results to examine issues with regard to convergence and dynamical stability of MD on DNA, and determine the significance of sequence context effects on all unique dinucleotide steps. This information is essential to understand sequence effects on DNA structure and has implications on diverse problems in the structural biology of DNA. Calculations were carried out on the 136 cases embedded in 39 DNA oligomers with repeating tetranucleotide sequences, capped on both ends by GC pairs and each having a total length of 15 nucleotide pairs. All simulations were carried out using a well-defined state-of-the-art MD protocol, the AMBER suite of programs, and the parm94 force field. In a previous article (Beveridge et al. 2004. Biophysical Journal. 87:3799-3813), the research design, details of the simulation protocol, and informatics issues were described. Preliminary results from 15 ns MD trajectories were presented for the d(CpG) step in all 10 unique sequence contexts. The results indicated the sequence context effects to be small for this step, but revealed that MD on DNA at this length of trajectory is subject to surprisingly persistent cooperative transitions of the sugar-phosphate backbone torsion angles alpha and gamma. In this article, we report detailed analysis of the entire trajectory database and occurrence of various conformational substates and its impact on studies of context effects. The analysis reveals a possible direct correspondence between the sequence-dependent dynamical tendencies of DNA structure and the tendency to undergo transitions that "trap" them in nonstandard conformational substates. The difference in mean of the observed basepair step helicoidal parameter distribution with different flanking sequence sometimes differs by as much as one standard deviation, indicating that the extent of sequence effects could be significant. The observations reveal that the impact of a flexible dinucleotide such as CpG could extend beyond the immediate basepair neighbors. The results in general provide new insight into MD on DNA and the sequence-dependent dynamical structural characteristics of DNA.     

Molecular dynamics simulations of the 136 unique tetranucleotide sequences of DNA oligonucleotides. I. Research design and results on d(CpG) steps

We describe herein a computationally intensive project aimed at carrying out molecular dynamics (MD) simulations including water and counterions on B-DNA oligomers containing all 136 unique tetranucleotide base sequences. This initiative was undertaken by an international collaborative effort involving nine research groups, the "Ascona B-DNA Consortium" (ABC). Calculations were carried out on the 136 cases imbedded in 39 DNA oligomers with repeating tetranucleotide sequences, capped on both ends by GC pairs and each having a total length of 15 nucleotide pairs. All MD simulations were carried out using a well-defined protocol, the AMBER suite of programs, and the parm94 force field. Phase I of the ABC project involves a total of approximately 0.6 mus of simulation for systems containing approximately 24,000 atoms. The resulting trajectories involve 600,000 coordinate sets and represent approximately 400 gigabytes of data. In this article, the research design, details of the simulation protocol, informatics issues, and the organization of the results into a web-accessible database are described. Preliminary results from 15-ns MD trajectories are presented for the d(CpG) step in its 10 unique sequence contexts, and issues of stability and convergence, the extent of quasiergodic problems, and the possibility of long-lived conformational substates are discussed.     

Studies of base pair sequence effects on DNA solvation based on all-atom molecular dynamics simulations

Detailed analyses of the sequence-dependent solvation and ion atmosphere of DNA are presented based on molecular dynamics (MD) simulations on all the 136 unique tetranucleotide steps obtained by the ABC consortium using the AMBER suite of programs. Significant sequence effects on solvation and ion localization were observed in these simulations. The results were compared to essentially all known experimental data on the subject. Proximity analysis was employed to highlight the sequence dependent differences in solvation and ion localization properties in the grooves of DNA. Comparison of the MD-calculated DNA structure with canonical A- and B-forms supports the idea that the G/C-rich sequences are closer to canonical A- than B-form structures, while the reverse is true for the poly A sequences, with the exception of the alternating ATAT sequence. Analysis of hydration density maps reveals that the flexibility of solute molecule has a significant effect on the nature of observed hydration. Energetic analysis of solute-solvent interactions based on proximity analysis of solvent reveals that the GC or CG base pairs interact more strongly with water molecules in the minor groove of DNA that the AT or TA base pairs, while the interactions of the AT or TA pairs in the major groove are stronger than those of the GC or CG pairs. Computation of solvent-accessible surface area of the nucleotide units in the simulated trajectories reveals that the similarity with results derived from analysis of a database of crystallographic structures is excellent. The MD trajectories tend to follow Manning's counterion condensation theory, presenting a region of condensed counterions within a radius of about 17 A from the DNA surface independent of sequence. The GC and CG pairs tend to associate with cations in the major groove of the DNA structure to a greater extent than the AT and TA pairs. Cation association is more frequent in the minor groove of AT than the GC pairs. In general, the observed water and ion atmosphere around the DNA sequences is the MD simulation is in good agreement with experimental observations.     

The ABCs of molecular dynamics simulations on B-DNA, circa 2012

This article provides a retrospective on the ABC initiative in the area of all-atom molecular dynamics (MD) simulations including explicit solvent on all tetranucleotide steps of duplex B-form DNA duplex, ca. 2012. The ABC consortium has completed two phases of simulations, the most current being a set of 50-100 trajectories based on the AMBER ff99 force field together with the parmbsc0 modification. Some general perspectives on the field of MD on DNA and sequence effects on DNA structure are provided, followed by an overview our MD results, including a detailed comparison of the ff99/parmbsc0 results with crystal and NMR structures available for d(CGCGAATTCGCG). Some projects inspired by or related to the ABC initiative and database are also reviewed, including methods for the trajectory analyses, informatics of dealing with the large database of results, compressions of trajectories for efficacy of distribution, DNA solvation by water and ions, parameterization of coarse-grained models with applications and gene finding and genome annotation.     

Exploring polymorphisms in B-DNA helical conformations

The traditional mesoscopic paradigm represents DNA as a series of base-pair steps whose energy response to equilibrium perturbations is elastic, with harmonic oscillations (defining local stiffness) around a single equilibrium conformation. In addition, base sequence effects are often analysed as a succession of independent XpY base-pair steps (i.e. a nearest-neighbour (NN) model with only 10 unique cases). Unfortunately, recent massive simulations carried out by the ABC consortium suggest that the real picture of DNA flexibility may be much more complex. The paradigm of DNA flexibility therefore needs to be revisited. In this article, we explore in detail one of the most obvious violations of the elastic NN model of flexibility: the bimodal distributions of some helical parameters. We perform here an in-depth statistical analysis of a very large set of MD trajectories and also of experimental structures, which lead to very solid evidence of bimodality. We then suggest ways to improve mesoscopic models to account for this deviation from the elastic regime.     

DNA reference alignment benchmarks based on tertiary structure of encoded proteins

MOTIVATION: Multiple sequence alignments (MSAs) are at the heart of bioinformatics analysis. Recently, a number of multiple protein sequence alignment benchmarks (i.e. BAliBASE, OXBench, PREFAB and SMART) have been released to evaluate new and existing MSA applications. These databases have been well received by researchers and help to quantitatively evaluate MSA programs on protein sequences. Unfortunately, analogous DNA benchmarks are not available, making evaluation of MSA programs difficult for DNA sequences. RESULTS: This work presents the first known multiple DNA sequence alignment benchmarks that are (1) comprised of protein-coding portions of DNA (2) based on biological features such as the tertiary structure of encoded proteins. These reference DNA databases contain a total of 3545 alignments, comprising of 68 581 sequences. Two versions of the database are available: mdsa_100s and mdsa_all. The mdsa_100s version contains the alignments of the data sets that TBLASTN found 100% sequence identity for each sequence. The mdsa_all version includes all hits with an E-value score above the threshold of 0.001. A primary use of these databases is to benchmark the performance of MSA applications on DNA data sets. The first such case study is included in the Supplementary Material.     

Prokaryotic gene finding based on physicochemical characteristics of codons calculated from molecular dynamics simulations

An ab initio model for gene prediction in prokaryotic genomes is proposed based on physicochemical characteristics of codons calculated from molecular dynamics (MD) simulations. The model requires a specification of three calculated quantities for each codon: the double-helical trinucleotide base pairing energy, the base pair stacking energy, and an index of the propensity of a codon for protein-nucleic acid interactions. The base pairing and stacking energies for each codon are obtained from recently reported MD simulations on all unique tetranucleotide steps, and the third parameter is assigned based on the conjugate rule previously proposed to account for the wobble hypothesis with respect to degeneracies in the genetic code. The third interaction propensity parameter values correlate well with ab initio MD calculated solvation energies and flexibility of codon sequences as well as codon usage in genes and amino acid composition frequencies in ∼175,000 protein sequences in the Swissprot database. Assignment of these three parameters for each codon enables the calculation of the magnitude and orientation of a cumulative three-dimensional vector for a DNA sequence of any length in each of the six genomic reading frames. Analysis of 372 genomes comprising ∼350,000 genes shows that the orientations of the gene and nongene vectors are well differentiated and make a clear distinction feasible between genic and nongenic sequences at a level equivalent to or better than currently available knowledge-based models trained on the basis of empirical data, presenting a strong support for the possibility of a unique and useful physicochemical characterization of DNA sequences from codons to genomes.     

EuSplice: a unified resource for the analysis of splice signals and alternative splicing in eukaryotic genes

MOTIVATION: Despite increased availability of genome annotation data, a comprehensive resource for in-depth analysis of splice signal distributions and alternative splicing (AS) patterns in eukaryote genomes is still lacking. To meet this need, we have developed EuSplice--a unique splice-centric database which provides reliable splice signal and AS information for 23 eukaryotes. RESULTS: The EuSplice database contains 95,822 AS events and 2.1 million splice signals associated with over 270,000 protein-coding genes. The intuitive, user-friendly EuSplice web interface has powerful data mining and graphics capabilities for inter-genomic comparative analysis of splice signals, putative cryptic splice sites and AS events. Moreover, the seamless integration of splicing data to extensive gene-specific annotations, such as homolog annotations, functional information, mutations and sequence details makes EuSplice a powerful one-stop information resource for investigating the molecular mechanisms of complex splicing events, disease associations and the evolution of splicing in eukaryotes. AVAILABILITY: http://66.170.16.154/EuSplice. SUPPLEMENTARY INFORMATION: Supplementary tables and figures at Bioinfo online.     

Molecular dynamics simulations of DNA curvature and flexibility: helix phasing and premelting

Recent studies of DNA axis curvature and flexibility based on molecular dynamics (MD) simulations on DNA are reviewed. The MD simulations are on DNA sequences up to 25 base pairs in length, including explicit consideration of counterions and waters in the computational model. MD studies are described for ApA steps, A-tracts, for sequences of A-tracts with helix phasing. In MD modeling, ApA steps and A-tracts in aqueous solution are essentially straight, relatively rigid, and exhibit the characteristic features associated with the B'-form of DNA. The results of MD modeling of A-tract oligonucleotides are validated by close accord with corresponding crystal structure results and nuclear magnetic resonance (NMR) nuclear Overhauser effect (NOE) and residual dipolar coupling (RDC) structures of d(CGCGAATTCGCG) and d(GGCAAAAAACGG). MD simulation successfully accounts for enhanced axis curvature in a set of three sequences with phased A-tracts studied to date. The primary origin of the axis curvature in the MD model is found at those pyrimidine/purine YpR "flexible hinge points" in a high roll, open hinge conformational substate. In the MD model of axis curvature in a DNA sequence with both phased A-tracts and YpR steps, the A-tracts appear to act as positioning elements that make the helix phasing more precise, and key YpR steps in the open hinge state serve as curvature elements. Our simulations on a phased A-tract sequence as a function of temperature show that the MD simulations exhibit a premelting transition in close accord with experiment, and predict that the mechanism involves a B'-to-B transition within A-tracts coupled with the prediction of a transition in key YpR steps from the high roll, open hinge, to a low roll, closed hinge substate. Diverse experimental observations on DNA curvature phenomena are examined in light of the MD model with no serious discrepancies. The collected MD results provide independent support for the "non-A-tract model" of DNA curvature. The "junction model" is indicated to be a special case of the non-A-tract model when there is a Y base at the 5' end of an A-tract. In accord with crystallography, the "ApA wedge model" is not supported by MD.     

Structure and axis curvature in two dA6 x dT6 DNA oligonucleotides: comparison of molecular dynamics simulations with results from crystallography and NMR spectroscopy

Molecular dynamics (MD) simulations have been performed on the A6 containing DNA dodecamers d(GGCAAAAAACGG) solved by NMR and d(CGCAAAAAAGCG) solved by crystallography. The experimental structures differ in the direction of axis bending and in other small but important aspects relevant to the DNA curvature problem. Five nanosecond MD simulations of each sequence have been performed, beginning with both the NMR and crystal forms as well as canonical B-form DNA. The results show that all simulations converge to a common form in close proximity to the observed NMR structure, indicating that the structure obtained in the crystal is likely a strained form due to packing effects. A-tracts in the MD model are essentially straight. The origin of axis curvature is found at pyrimidine-purine steps in the flanking sequences.     

Conformational flexibility of two RNA trimers explored by computational tools and database search

Two RNA sequences, AAA and AUG, were studied by the conformational search program CICADA and by molecular dynamics (MD) in the framework of the AMBER force field, and also via thorough PDB database search. CICADA was used to provide detailed information about conformers and conformational interconversions on the energy surfaces of the above molecules. Several conformational families were found for both sequences. Analysis of the results shows differences, especially between the energy of the single families, and also in flexibility and concerted conformational movement. Therefore, several MD trajectories (altogether 16 ns) were run to obtain more details about both the stability of conformers belonging to different conformational families and about the dynamics of the two systems. Results show that the trajectories strongly depend on the starting structure. When the MD start from the global minimum found by CICADA, they provide a stable run, while MD starting from another conformational family generates a trajectory where several different conformational families are visited. The results obtained by theoretical methods are compared with the thorough database search data. It is concluded that all except for the highest energy conformational families found in theoretical result also appear in experimental data. Registry numbers: adenylyl-(3' --> 5')-adenylyl-(3' --> 5')-adenosine [917-44-2] adenylyl-(3' --> 5')-uridylyl-(3' --> 5')-guanosine [3494-35-7].     

Classical molecular interaction potentials: improved setup procedure in molecular dynamics simulations of proteins.

The latest version of the classical molecular interaction potential (CMIP) has the ability to predict the position of crystallographic waters in several proteins with great accuracy. This article analyzes the ability of the CMIP functional to improve the setup procedure of the molecular system in molecular dynamics (MD) simulations of proteins. To this end, the CMIP strategy is used to include both water molecules and counterions in different protein systems. The structural details of the configurations sampled from trajectories obtained using the CMIP setup procedure are compared with those obtained from trajectories derived from a standard equilibration process. The results show that standard MD simulations can lead to artifactual results, which are avoided using the CMIP setup procedure. Because the CMIP is easy to implement at a low computational cost, it can be very useful in obtaining reliable MD trajectories.     

Conformational Flexibility of Two RNA Trimers Explored by Computational Tools and Database Search

Abstract

Two RNA sequences, AAA and AUG, were studied by the conformational search program CICADA and by molecular dynamics (MD) in the framework of the AMBER force field, and also via thorough PDB database search. CICADA was used to provide detailed information about conformers and conformational interconversions on the energy surfaces of the above molecules. Several conformational families were found for both sequences. Analysis of the results shows differences, especially between the energy of the single families, and also in flexibility and concerted conformational movement. Therefore, several MD trajectories (altogether 16 ns) were run to obtain more details about both the stability of conformers belonging to different conformational families and about the dynamics of the two systems. Results show that the trajectories strongly depend on the starting structure. When the MD start from the global minimum found by CICADA, they provide a stable run, while MD starting from another conformational family generates a trajectory where several different conformational families are visited. The results obtained by theoretical methods are compared with the thorough database search data. It is concluded that all except for the highest energy conformational families found in theoretical result also appear in experimental data.

Registry numbers:

adenylyl-(3′ →5′)-adenylyl-(3′ →5′)-adenosine [917-44-2]

adenylyl-(3′ →5′)-uridylyl-(3′ →5′)-guanosine [3494-35-7]     

Construction, MD Simulation, and Hydrodynamic Validation of an All-Atom Model of a Monoclonal IgG Antibody

At 150 kDa, antibodies of the IgG class are too large for their structure to be determined with current NMR methodologies. Because of hinge-region flexibility, it is difficult to obtain atomic-level structural information from the crystal, and questions regarding antibody structure and dynamics in solution remain unaddressed. Here we describe the construction of a model of a human IgG1 monoclonal antibody (trastuzumab) from the crystal structures of fragments. We use a combination of molecular-dynamics (MD) simulation, continuum hydrodynamics modeling, and experimental diffusion measurements to explore antibody behavior in aqueous solution. Hydrodynamic modeling provides a link between the atomic-level details of MD simulation and the size- and shape-dependent data provided by hydrodynamic measurements. Eight independent 40 ns MD trajectories were obtained with the AMBER program suite. The ensemble average of the computed transport properties over all of the MD trajectories agrees remarkably well with the value of the translational diffusion coefficient obtained with dynamic light scattering at 20°C and 27°C, and the intrinsic viscosity measured at 20°C. Therefore, our MD results likely represent a realistic sampling of the conformational space that an antibody explores in aqueous solution.     

A systematic molecular dynamics study of nearest-neighbor effects on base pair and base pair step conformations and fluctuations in B-DNA

It is well recognized that base sequence exerts a significant influence on the properties of DNA and plays a significant role in protein–DNA interactions vital for cellular processes. Understanding and predicting base sequence effects requires an extensive structural and dynamic dataset which is currently unavailable from experiment. A consortium of laboratories was consequently formed to obtain this information using molecular simulations. This article describes results providing information not only on all 10 unique base pair steps, but also on all possible nearest-neighbor effects on these steps. These results are derived from simulations of 50–100 ns on 39 different DNA oligomers in explicit solvent and using a physiological salt concentration. We demonstrate that the simulations are converged in terms of helical and backbone parameters. The results show that nearest-neighbor effects on base pair steps are very significant, implying that dinucleotide models are insufficient for predicting sequence-dependent behavior. Flanking base sequences can notably lead to base pair step parameters in dynamic equilibrium between two conformational sub-states. Although this study only provides limited data on next-nearest-neighbor effects, we suggest that such effects should be analyzed before attempting to predict the sequence-dependent behavior of DNA.     

Determining a unique defining DNA sequence for yeast species using hashing techniques

MOTIVATION: Yeasts are often still identified with physiological growth tests, which are both time consuming and unsuitable for detection of a mixture of organisms. Hence, there is a need for molecular methods to identify yeast species. RESULTS: A hashing technique has been developed to search for unique DNA sequences in 702 26S rRNA genes. A unique DNA sequence has been found for almost every yeast species described to date. The locations of the unique defining sequences are in accordance with the variability map of large subunit ribosomal RNA and provide detail of the evolution of the D1/D2 region. This approach will be applicable to the rapid identification of unique sequences in other DNA sequence sets. AVAILABILITY: Freely available upon request from the authors. Supplementary information: Results are available at http://www.sys.uea.ac.uk/~jjw/project/paper     

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.