Calibrating nucleic acids torsional energetics in force-field: insights from model compounds

Simulations Meet Experiment to Reveal New Insights into DNA Intrinsic Mechanics

The accurate prediction of the structure and dynamics of DNA remains a major challenge in computational biology due to the dearth of precise experimental information on DNA free in solution and limitations in the DNA force-fields underpinning the simulations. A new generation of force-fields has been developed to better represent the sequence-dependent B-DNA intrinsic mechanics, in particular with respect to the BI ↔ BII backbone equilibrium, which is essential to understand the B-DNA properties. Here, the performance of MD simulations with the newly updated force-fields Parmbsc0εζOLI and CHARMM36 was tested against a large ensemble of recent NMR data collected on four DNA dodecamers involved in nucleosome positioning. We find impressive progress towards a coherent, realistic representation of B-DNA in solution, despite residual shortcomings. This improved representation allows new and deeper interpretation of the experimental observables, including regarding the behavior of facing phosphate groups in complementary dinucleotides, and their modulation by the sequence. It also provides the opportunity to extensively revisit and refine the coupling between backbone states and inter base pair parameters, which emerges as a common theme across all the complementary dinucleotides. In sum, the global agreement between simulations and experiment reveals new aspects of intrinsic DNA mechanics, a key component of DNA-protein recognition.

Energies and 2'-Hydroxyl Group Orientations of RNA Backbone Conformations. Benchmark CCSD(T)/CBS Database, Electronic Analysis, and Assessment of DFT Methods and MD Simulations

Sugar-phosphate backbone is an electronically complex molecular segment imparting RNA molecules high flexibility and architectonic heterogeneity necessary for their biological functions. The structural variability of RNA molecules is amplified by the presence of the 2'-hydroxyl group, capable of forming multitude of intra- and intermolecular interactions. Bioinformatics studies based on X-ray structure database revealed that RNA backbone samples at least 46 substates known as rotameric families. The present study provides a comprehensive analysis of RNA backbone conformational preferences and 2'-hydroxyl group orientations. First, we create a benchmark database of estimated CCSD(T)/CBS relative energies of all rotameric families and test performance of dispersion-corrected DFT-D3 methods and molecular mechanics in vacuum and in continuum solvent. The performance of the DFT-D3 methods is in general quite satisfactory. The B-LYP-D3 method provides the best trade-off between accuracy and computational demands. B3-LYP-D3 slightly outperforms the new PW6B95-D3 and MPW1B95-D3 and is the second most accurate density functional of the study. The best agreement with CCSD(T)/CBS is provided by DSD-B-LYP-D3 double-hybrid functional, although its large-scale applications may be limited by high computational costs. Molecular mechanics does not reproduce the fine energy differences between the RNA backbone substates. We also demonstrate that the differences in the magnitude of the hyperconjugation effect do not correlate with the energy ranking of the backbone conformations. Further, we investigated the 2'-hydroxyl group orientation preferences. For all families, we conducted a QM and MM hydroxyl group rigid scan in gas phase and solvent. We then carried out set of explicit solvent MD simulations of folded RNAs and analyze 2'-hydroxyl group orientations of different backbone families in MD. The solvent energy profiles determined primarily by the sugar pucker match well with the distribution data derived from the simulations. The QM and MM energy profiles predict the same 2'-hydroxyl group orientation preferences. Finally, we demonstrate that the high energy of unfavorable and rarely sampled 2'-hydroxyl group orientations can be attributed to clashes between occupied orbitals.

Simulating DNA by molecular dynamics: aims, methods, and validation

The structure and dynamics of the B-DNA double helix involves subtle sequence-dependent effects which are decisive for its function, but difficult to characterize. These structural and dynamic effects can be addressed by simulations of DNA sequences in explicit solvent. Here, we present and discuss the state-of-art of B-DNA molecular dynamics simulations with the major force fields in use today. We explain why a critical analysis of the MD trajectories is required to assess their reliability, and estimate the value and limitations of these models. Overall, simulations of DNA bear great promise towards deciphering the structural and physical subtleties of this biopolymer, where much remains to be understood.

Understanding the Sequence Preference of Recurrent RNA Building Blocks using Quantum Chemistry: The Intrastrand RNA Dinucleotide Platform

Folded RNA molecules are shaped by an astonishing variety of highly conserved noncanonical molecular interactions and backbone topologies. The dinucleotide platform is a widespread recurrent RNA modular building submotif formed by the side-by-side pairing of bases from two consecutive nucleotides within a single strand, with highly specific sequence preferences. This unique arrangement of bases is cemented by an intricate network of noncanonical hydrogen bonds and facilitated by a distinctive backbone topology. The present study investigates the gas-phase intrinsic stabilities of the three most common RNA dinucleotide platforms - 5'-GpU-3', ApA, and UpC - via state-of-the-art quantum-chemical (QM) techniques. The mean stability of base-base interactions decreases with sequence in the order GpU > ApA > UpC. Bader's atoms-in-molecules analysis reveals that the N2(G)…O4(U) hydrogen bond of the GpU platform is stronger than the corresponding hydrogen bonds in the other two platforms. The mixed-pucker sugar-phosphate backbone conformation found in most GpU platforms, in which the 5'-ribose sugar (G) is in the C2'-endo form and the 3'-sugar (U) in the C3'-endo form, is intrinsically more stable than the standard A-RNA backbone arrangement, partially as a result of a favorable O2'…O2P intra-platform interaction. Our results thus validate the hypothesis of Lu et al. (Lu Xiang-Jun, et al. Nucleic Acids Res. 2010, 38, 4868-4876), that the superior stability of GpU platforms is partially mediated by the strong O2'…O2P hydrogen bond. In contrast, ApA and especially UpC platform-compatible backbone conformations are rather diverse and do not display any characteristic structural features. The average stabilities of ApA and UpC derived backbone conformers are also lower than those of GpU platforms. Thus, the observed structural and evolutionary patterns of the dinucleotide platforms can be accounted for, to a large extent, by their intrinsic properties as described by modern QM calculations. In contrast, we show that the dinucleotide platform is not properly described in the course of atomistic explicit-solvent simulations. Our work also gives methodological insights into QM calculations of experimental RNA backbone geometries. Such calculations are inherently complicated by rather large data and refinement uncertainties in the available RNA experimental structures, which often preclude reliable energy computations.

Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles

We report a reparameterization of the glycosidic torsion χ of the Cornell et al. AMBER force field for RNA, χ(OL). The parameters remove destabilization of the anti region found in the ff99 force field and thus prevent formation of spurious ladder-like structural distortions in RNA simulations. They also improve the description of the syn region and the syn-anti balance as well as enhance MD simulations of various RNA structures. Although χ(OL) can be combined with both ff99 and ff99bsc0, we recommend the latter. We do not recommend using χ(OL) for B-DNA because it does not improve upon ff99bsc0 for canonical structures. However, it might be useful in simulations of DNA molecules containing syn nucleotides. Our parametrization is based on high-level QM calculations and differs from conventional parametrization approaches in that it incorporates some previously neglected solvation-related effects (which appear to be essential for obtaining correct anti/high-anti balance). Our χ(OL) force field is compared with several previous glycosidic torsion parametrizations.

Development of CHARMM polarizable force field for nucleic acid bases based on the classical Drude oscillator model

A polarizable force field for nucleic acid bases based on the classical Drude oscillator model is presented. Parameter optimization was performed to reproduce crystallographic geometries, crystal unit cell parameters, heats of sublimation, vibrational frequencies and assignments, dipole moments, molecular polarizabilities and quantum mechanical base-base and base-water interaction energies. The training and validation data included crystals of unsubstituted and alkyl-substituted adenine, guanine, cytosine, uracil, and thymine bases, hydrated crystals, and hydrogen bonded base pairs. Across all compounds, the RMSD in the calculated heats of sublimation is 4.1%. This equates to an improvement of more than 2.5 kcal/mol in accuracy compared to the nonpolarizable CHARMM27 force field. However, the level of agreement with experimental molecular volume decreased from 1.7% to 2.1% upon moving from the nonpolarizable to the polarizable model. The representation of dipole moments is significantly improved with the Drude polarizable force field. Unlike in additive force fields, there is no requirement for the gas-phase dipole moments to be overestimated, illustrating the ability of the Drude polarizable force field to treat accurately differently dielectric environments and indicating the improvements in the electrostatic model. Validation of the model was performed on the basis of the calculation of the gas-phase binding enthalpies of base pairs obtained via potential of mean force calculations; the additive and polarizable models both performed satisfactorily with average differences of 0.2 and 0.9 kcal/mol, respectively, and rms differences of 1.3 and 1.7 kcal/mol, respectively. Overall, considering the number of significant improvements versus the additive CHARMM force field, the incorporation of explicit polarizability into the force field for nucleic acid bases represents an additional step toward accurate computational modeling of biological systems.

Toward a full characterization of nucleic acid components in aqueous solution: simulations of nucleosides

The eight nucleoside constituents of nucleic acids were simulated for 50 ns in explicit water with molecular dynamics. This provides equilibrium populations of the torsional degrees of freedom, their kinetics of interconversion, their couplings, and how they are influenced by water. This is important, given that a full and quantitative characterization of the nucleosides in aqueous solution by experimental means has been elusive, despite immense efforts in that direction. It is with the anti/syn equilibrium that the simulations are most complementary to experiment, by accessing directly the influence of the sugar type, sugar pucker, and base on the anti/syn populations. The glycosidic torsion distributions in the anti conformation are strongly affected by water and depart from the corresponding X-ray modal values and the associated energy minima in vacuo. Water also preferentially stabilizes some sugar conformations, showing that potential energies in vacuo are not sufficient to understand the nucleosides. Deoxythymidine (but not other pyrimidines) significantly populates the syn orientation. Guanine favors the syn orientation more than adenine. The ribose favors the syn orientation significantly more than the deoxyribose. The NORTH pucker coexists with the syn conformers. A hydrogen bond is frequently formed between the 5'-OH group and the syn bases, despite competition by water. The rate of the anti/syn transitions with purines is on the nanosecond time scale, confirming a long held assumption underpinning the interpretation of ultrasonic relaxation studies. Therefore, our knowledge of the structure and dynamics of nucleosides in solvent is only limited by the accuracy of the potential used to simulate them, and it is shown that such simulations provide a distinct and unique test of nucleic acid force fields. This confirmed that the widely distributed CHARMM27 force field is, overall, well-balanced with a particularly good representation of the ribose. Specific improvements, however, are suggested for the deoxyribose and torsion gamma.

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.

Conformational characteristics and correlations in crystal structures of nucleic acid oligonucleotides: evidence for sub-states

Sugar phosphate backbone conformations are a structural element inextricably involved in a complete understanding of specific recognition nucleic acid ligand interactions, from early stage discrimination of the correct target to complexation per se, including any structural adaptation on binding. The collective results of high resolution DNA, RNA and protein/DNA crystal structures provide an opportunity for an improved and enhanced statistical analysis of standard and unusual sugar-phosphate backbone conformations together with corresponding dinucleotide sequence effects as a basis for further exploration of conformational effects on binding. In this study, we have analyzed the conformations of all relevant crystal structures in the nucleic acids data base, determined the frequency distribution of all possible epsilon, zeta, alpha, beta and gamma backbone angle arrangements within four nucleic acid categories (A-RNA and A-DNA, free and bound B-DNA) and explored the relationships between backbone angles, sugar puckers and selected helical parameters. The trends in the correlations are found to be similar regardless of the nucleic acid category. It is interesting that specific structural effects exhibited by the different unusual backbone sub-states are in some cases contravariant. Certain alpha/gamma changes are accompanied by C3' endo (north) sugars, small twist angles and positive values of base pair roll, and favor a displacement of nucleotide bases towards the minor groove compared to that of canonical B form structures. Unusual epsilon/zeta combinations occur with C2' (south) sugars, high twist angles, negative values of base pair roll, and base displacements towards the major groove. Furthermore, any unusual backbone correlates with a reduced dispersion of equilibrium structural parameters of the whole double helix, as evidenced by the reduced standard deviations of almost all conformational parameters. Finally, a strong sequence effect is displayed in the free oligomers, but reduced somewhat in the ligand bound forms. The most variable steps are GpA and CpA, and, to a lesser extent, their partners TpC and TpG. The results provide a basis for considering if the variable and non-variable steps within a biological active sequence precisely determine morphological structural features as the curvature direction, the groove depth, and the accessibility of base pair for non covalent associations.