Search

The RNA hairpin loops represent important RNA topologies with indispensable biological functions in RNA folding and tertiary interactions. 5′-UNCG-3′ and 5′-GNRA-3′ RNA tetraloops are the most important classes of RNA hairpin loops. Both tetraloops are highly structured with characteristic signature three-dimensional features and are recurrently seen in functional RNAs and ribonucleoprotein particles. Explicit solvent molecular dynamics (MD) simulation is a computational technique which can efficiently complement the experimental data and provide unique structural dynamics information on the atomic scale. Nevertheless, the outcome of simulations is often compromised by imperfections in the parametrization of simplified pairwise additive empirical potentials referred to also as force fields. We have pointed out in several recent studies that a force field description of single-stranded hairpin segments of nucleic acids may be particularly challenging for the force fields. In this paper, we report a critical assessment of a broad set of MD simulations of UUCG, GAGA, and GAAA tetraloops using various force fields. First, we utilized the three widely used variants of Cornell et al. (AMBER) force fields known as ff94, ff99, and ff99bsc0. Some simulations were also carried out with CHARMM27. The simulations reveal several problems which show that these force fields are not able to retain all characteristic structural features (structural signature) of the studied tetraloops. Then we tested four recent reparameterizations of glycosidic torsion of the Cornell et al. force field (two of them being currently parametrized in our laboratories). We show that at least some of the new versions show an improved description of the tetraloops, mainly in the syn glycosidic torsion region of the UNCG tetraloop. The best performance is achieved in combination with the bsc0 parametrization of the α/γ angles. Another critically important region to properly describe RNA molecules is the anti/high-anti region of the glycosidic torsion, where there are significant differences among the tested force fields. The tetraloop simulations are complemented by simulations of short A-RNA stems, which are especially sensitive to an appropriate description of the anti/high-anti region. While excessive accessibility of the high-anti region converts the A-RNA into a senseless “ladder-like” geometry, excessive penalization of the high-anti region shifts the simulated structures away from typical A-RNA geometry to structures with a visibly underestimated inclination of base pairs with respect to the helical axis.

When DNA interacts with a protein, its structure often undergoes significant conformational adaptation. Perhaps the most common is the transition from canonical B-DNA towards the A-DNA form, which is not a two-state, but rather a continuous transition. The A- and B-forms differ mainly in sugar pucker P (north/south) and glycosidic torsion χ (high-anti/anti). The combination of A-like P and B-like χ (and vice versa) represents the nature of the intermediate states lying between the pure A- and B- forms. In this work, we study how the A/B equilibrium and in particular the A/B intermediate states, which are known to be over-represented at protein-DNA interfaces, are modeled by current AMBER force fields. Eight protein-DNA complexes and their naked (unbound) DNAs were simulated with OL15 and bsc1 force fields as well as an experimental combination OL15χOL3. We found that while the geometries of the A-like intermediate states in the molecular dynamics (MD) simulations agree well with the native X-ray geometries found in the protein-DNA complexes, their populations (stabilities) are significantly underestimated. Different force fields predict different propensities for A-like states growing in the order OL15 < bsc1 < OL15χOL3, but the overall populations of the A-like form are too low in all of them. Interestingly, the force fields seem to predict the correct sequence-dependent A-form propensity, as they predict larger populations of the A-like form in naked (unbound) DNA in those steps that acquire A-like conformations in protein-DNA complexes. The instability of A-like geometries in current force fields may significantly alter the geometry of the simulated protein-DNA complex, destabilize the binding motif, and reduce the binding energy, suggesting that refinement is needed to improve description of protein-DNA interactions in AMBER force fields.

A method for averaging of NMR parameters by molecular dynamics (MD) has been derived from the method of statistical averaging in MD snapshots, benchmarked and applied to structurally dynamic interpretation of the

Although current AMBER force fields are relatively accurate for canonical B-DNA, many non-canonical structures are still described incorrectly. As non-canonical motifs are attracting increasing attention due to the role they play in living organisms, further improvement is desirable. Here, we have chosen Z-DNA molecule, can be considered a touchstone of the universality of empirical force fields, since the non-canonical α and γ backbone conformations native to Z-DNA are also found in protein-DNA complexes, i-motif DNA and other non-canonical DNAs. We show that spurious α/γ conformations occurring in simulations with current AMBER force fields, OL15 and bsc1, are largely due to inaccurate α/γ parameterization. Moreover, stabilization of native Z-DNA substates involving γ = trans conformations appears to be in conflict with the correct description of the canonical B-DNA structure. Because the balance of the native and spurious conformations is influenced by non-additive effects, this is a difficult case for an additive dihedral energy scheme such as AMBER. We propose new α/γ parameters, denoted OL21, and show that they improve the stability of native α/γ Z-DNA substates while keeping the canonical DNA description virtually unchanged, and thus represent a reasonable compromise within the additive force field framework. Although further extensive testing is needed, the new modification appears to be a promising step towards a more reliable description of non-canonical DNA motifs and provides the best performance for Z-DNA molecules among current AMBER force fields.

The A/B transition is a basic element of DNA conformational change. Because of its involvement in the sensing of the ionic conditions by DNA and in specific protein-DNA interactions, this transition is important for biological functions of DNA. Therefore, accurate modeling of the A/B equilibrium by means of empirical force fields is of utmost interest. In this work, we examine the A/B equilibrium in three AMBER force fields, including the recent bsc1 and OL15 modifications, using much longer MD simulations than attempted before. Special attention is paid to the coupling of the A/B equilibrium with the south/north (S/N) transition of the sugar pucker. We found that none of the tested force fields provided a satisfactory description of the A/B equilibrium because the B-form was predicted to be much too stable and the A-form was predicted to be almost absent even in concentrated trifluoroethanol solutions. Based on comparison with NMR data for duplexes and single nucleosides, we hypothesize that this problem arose from the incorrect description of the S/N equilibrium of sugar pucker, where the south conformation is much too stable, thus stabilizing the B-form. Because neither the A/B equilibrium in duplexes nor the S/N equilibrium in nucleosides was described accurately, further refinements of the AMBER DNA force fields are needed.

The sugar–phosphate backbone of RNA can exist in diverse rotameric substates, giving RNA molecules enormous conformational variability. The most frequent noncanonical backbone conformation in RNA is α/γ = t/t, which is derived from the canonical backbone by a crankshaft motion and largely preserves the standard geometry of the RNA duplex. A similar conformation also exists in DNA, where it has been extensively studied and shown to be involved in DNA–protein interactions. However, the function of the α/γ = t/t conformation in RNA is poorly understood. Here, we present molecular dynamics simulations of several prototypical RNA structures obtained from X-ray and NMR experiments, including canonical and mismatched RNA duplexes, UUCG and GAGA tetraloops, Loop E, the sarcin–ricin loop, a parallel guanine quadruplex, and a viral pseudoknot. The stability of various noncanonical α/γ backbone conformations was analyzed with two AMBER force fields, ff99bsc0χOL3 and ff99bsc0χOL3 with the recent eζOL1 and βOL1 correc...

Molecular dynamics (MD) simulations became a leading tool for investigation of structural dynamics of nucleic acids. Despite recent efforts to improve the empirical potentials (force fields,ffs), RNAffshave persisting deficiencies, which hamper their utilization in quantitatively accurate simulations. Previous studies have shown that at least two salient problems contribute to difficulties in description of free-energy landscapes of small RNA motifs: (i) excessive stabilization of the unfolded single-stranded RNA ensemble by intramolecular base-phosphate and sugar-phosphate interactions, and (ii) destabilization of the native folded state by underestimation of stability of base pairing. Here, we introduce a generalffterm (gHBfix) that can selectively fine-tune non-bonding interaction terms in RNAffs, in particular the H-bonds. gHBfix potential affects the pair-wise interactions between all possible pairs of the specific atom types, while all other interactions remain intact, i.e., it is not a structure-based model. In order to probe the ability of the gHBfix potential to refine theffnon-bonded terms, we performed an extensive set of folding simulations of RNA tetranucleotides and tetraloops. Based on these data we propose particular gHBfix parameters to modify the AMBER RNAff. The suggested parametrization significantly improves the agreement between experimental data and the simulation conformational ensembles, although our currentffversion still remains far from being flawless. While attempts to tune the RNAffsby conventional reparametrizations of dihedral potentials or non-bonded terms can lead to major undesired side effects as we demonstrate for some recently publishedffs, gHBfix has a clear promising potential to improve theffperformance while avoiding introduction of major new imbalances.

Molecular dynamics (MD) simulations became a leading tool for investigation of structural dynamics of nucleic acids. Despite recent efforts to improve the empirical potentials (force fields, ffs), RNA ffs have persisting deficiencies, which hamper their utilization in quantitatively accurate simulations. Previous studies have shown that at least two salient problems contribute to difficulties in description of free-energy landscapes of small RNA motifs: (i) excessive stabilization of the unfolded single-stranded RNA ensemble by intramolecular base-phosphate and sugar-phosphate interactions, and (ii) destabilization of the native folded state by underestimation of stability of base pairing. Here, we introduce a general ff term (gHBfix) that can selectively fine-tune non-bonding interaction terms in RNA ffs, in particular the H-bonds. gHBfix potential affects the pair-wise interactions between all possible pairs of the specific atom types, while all other interactions remain intact, i.e., it is not a structure-based model. In order to probe the ability of the gHBfix potential to refine the ff non-bonded terms, we performed an extensive set of folding simulations of RNA tetranucleotides and tetraloops. Based on these data we propose particular gHBfix parameters to modify the AMBER RNA ff. The suggested parametrization significantly improves the agreement between experimental data and the simulation conformational ensembles, although our current ff version still remains far from being flawless. While attempts to tune the RNA ffs by conventional reparametrizations of dihedral potentials or non-bonded terms can lead to major undesired side effects as we demonstrate for some recently published ffs, gHBfix has a clear promising potential to improve the ff performance while avoiding introduction of major new imbalances.

We report an extensive set of explicit solvent molecular dynamics (MD) simulations (∼25 μs of accumulated simulation time) of the RNA kissing-loop complex of the HIV-1 virus initiation dimerization site. Despite many structural investigations by X-ray, NMR, and MD techniques, the position of the bulged purines of the kissing complex has not been unambiguously resolved. The X-ray structures consistently show bulged-out positions of the unpaired bases, while several NMR studies show bulged-in conformations. The NMR studies are, however, mutually inconsistent regarding the exact orientations of the bases. The earlier simulation studies predicted the bulged-out conformation; however, this finding could have been biased by the short simulation time scales. Our microsecond-long simulations reveal that all unpaired bases of the kissing-loop complex stay preferably in the interior of the kissing-loop complex. The MD results are discussed in the context of the available experimental data and we suggest that both conformations are biochemically relevant. We also show that MD provides a quite satisfactory description of this RNA system, contrasting recent reports of unsatisfactory performance of the RNA force fields for smaller systems such as tetranucleotides and tetraloops. We explain this by the fact that the kissing complex is primarily stabilized by an extensive network of Watson-Crick interactions which are rather well described by the force fields. We tested several different sets of water/ion parameters but they all lead to consistent results. However, we demonstrate that a recently suggested modification of van der Waals interactions of the Cornell et al. force field deteriorates the description of the kissing complex by the loss of key stacking interactions stabilizing the interhelical junction and excessive hydrogen-bonding interactions.

Z-DNA duplexes are a particularly complicated test case for current force fields. We performed a set of explicit solvent molecular dynamics (MD) simulations with various AMBER force field parametrizations including our recent refinements of the ε/ζ and glycosidic torsions. None of these force fields described the ZI/ZII and other backbone substates correctly, and all of them underpredicted the population of the important ZI substate. We show that this underprediction can be attributed to an inaccurate potential for the sugar-phosphate backbone torsion angle β. We suggest a refinement of this potential, β(OL1), which was derived using our recently introduced methodology that includes conformation-dependent solvation effects. The new potential significantly increases the stability of the dominant ZI backbone substate and improves the overall description of the Z-DNA backbone. It also has a positive (albeit small) impact on another important DNA form, the antiparallel guanine quadruplex (G-DNA), and improves the description of the canonical B-DNA backbone by increasing the population of BII backbone substates, providing a better agreement with experiment. We recommend using β(OL1) in combination with our previously introduced corrections, εζ(OL1) and χ(OL4), (the combination being named OL15) as a possible alternative to the current β torsion potential for more accurate modeling of nucleic acids.

A recent study described an allosteric effect in which the binding of a protein to DNA is influenced by another protein bound nearby. The effect shows a periodicity of ∼10 basepairs and decays with increasing protein-protein distance. As a mechanistic explanation, the authors reported a similar periodic, decaying pattern of the correlation coefficient between major groove widths inferred from a shorter molecular dynamics simulation. Here we show that in a state-of-the-art, microsecond-long simulation of the same DNA sequence, the periodicity of the correlation coefficient is not observed. To study the problem further, we extend an earlier mechanical model of DNA allostery based on constrained minimization of effective quadratic deformation energy of the DNA. We demonstrate that, if the constraints mimicking the bound proteins are properly applied, the periodicity in the binding energy is indeed recovered.

We present a refinement of the backbone torsion parameters ε and ζ of the Cornell et al. AMBER force field for DNA simulations. The new parameters, denoted as εζOL1, were derived from quantum-mechanical calculations with inclusion of conformation-dependent solvation effects according to the recently reported methodology (J. Chem. Theory Comput. 2012, 7(9), 2886-2902). The performance of the refined parameters was analyzed by means of extended molecular dynamics (MD) simulations for several representative systems. The results showed that the εζOL1 refinement improves the backbone description of B-DNA double helices and G-DNA stem. In B-DNA simulations, we observed an average increase of the helical twist and narrowing of the major groove, thus achieving better agreement with X-ray and solution NMR data. The balance between populations of BI and BII backbone substates was shifted towards the BII state, in better agreement with ensemble-refined solution experimental results. Furthermore, the refined parameters decreased the backbone RMS deviations in B-DNA MD simulations. In the antiparallel guanine quadruplex (G-DNA) the εζOL1 modification improved the description of non-canonical α/γ backbone substates, which were shown to be coupled to the ε/ζ torsion potential. Thus, the refinement is suggested as a possible alternative to the current ε/ζ torsion potential, which may enable more accurate modeling of nucleic acids. However, long-term testing is recommended before its routine application in DNA simulations.

We provide a critical assessment of explicit-solvent atomistic molecular dynamics (MD) simulations of RNA and protein/RNA complexes, written primarily for non-specialists with an emphasis to explain the limitations of MD. MD simulations can be likened to hypothetical single-molecule experiments starting from single atomistic conformations and investigating genuine thermal sampling of the biomolecules. The main advantage of MD is the unlimited temporal and spatial resolution of positions of all atoms in the simulated systems. Fundamental limitations are the short physical time-scale of simulations, which can be partially alleviated by enhanced-sampling techniques, and the highly approximate atomistic force fields describing the simulated molecules. The applicability and present limitations of MD are demonstrated on studies of tetranucleotides, tetraloops, ribozymes, riboswitches and protein/RNA complexes. Wisely applied simulations respecting the approximations of the model can successfully complement structural and biochemical experiments. WIREs RNA 2017, 8:e1405. doi: 10.1002/wrna.1405 For further resources related to this article, please visit the WIREs website.

The utility of molecular dynamics (MD) simulations to model biomolecular structure, dynamics, and interactions has witnessed enormous advances in recent years due to the availability of optimized MD software and access to significant computational power, including GPU multicore computing engines and other specialized hardware. This has led researchers to routinely extend conformational sampling times to the microsecond level and beyond. The extended sampling time has allowed the community not only to converge conformational ensembles through complete sampling but also to discover deficiencies and overcome problems with the force fields. Accuracy of the force fields is a key component, along with sampling, toward being able to generate accurate and stable structures of biopolymers. The Amber force field for nucleic acids has been used extensively since the 1990s, and multiple artifacts have been discovered, corrected, and reassessed by different research groups. We present a direct comparison of two of the most recent and state-of-the-art Amber force field modifications, bsc1 and OL15, that focus on accurate modeling of double-stranded DNA. After extensive MD simulations with five test cases and two different water models, we conclude that both modifications are a remarkable improvement over the previous bsc0 force field. Both force field modifications show better agreement when compared to experimental structures. To ensure convergence, the Drew–Dickerson dodecamer (DDD) system was simulated using 100 independent MD simulations, each extended to at least 10 μs, and the independent MD simulations were concatenated into a single 1 ms long trajectory for each combination of force field and water model. This is significantly beyond the time scale needed to converge the conformational ensemble of the internal portions of a DNA helix absent internal base pair opening. Considering all of the simulations discussed in the current work, the MD simulations performed to assess and validate the current force fields and water models aggregate over 14 ms of simulation time. The results suggest that both the bsc1 and OL15 force fields render average structures that deviate significantly less than 1 Å from the average experimental structures. This can be compared to similar but less exhaustive simulations with the CHARMM 36 force field that aggregate to the ∼90 μs time scale and also perform well but do not produce structures as close to the DDD NMR average structures (with root-mean-square deviations of 1.3 Å) as the newer Amber force fields. On the basis of these analyses, any future research involving double-stranded DNA simulations using the Amber force fields should employ the bsc1 or OL15 modification.

A procedure for deriving force field torsion parameters including certain previously neglected solvation effects is suggested. In contrast to the conventional in vacuo approaches, the dihedral parameters are obtained from the difference between the quantum-mechanical self-consistent reaction field and Poisson-Boltzmann continuum solvation models. An analysis of the solvation contributions shows that two major effects neglected when torsion parameters are derived in vacuo are (i) conformation-dependent solute polarization and (ii) solvation of conformation-dependent charge distribution. Using the glycosidic torsion as an example, we demonstrate that the corresponding correction for the torsion potential is substantial and important. Our approach avoids double counting of solvation effects and provides parameters that may be used in combination with any of the widely used nonpolarizable discrete solvent models, such as TIPnP or SPC/E, or with continuum solvent models. Differences between our model and the previously suggested solvation models are discussed. Improvements were demonstrated for the latest AMBER RNA χOL3 parameters derived with inclusion of solvent effects in a previous publication (Zgarbova et al. J. Chem. Theory Comput.2011, 7, 2886). The described procedure may help to provide consistently better force field parameters than the currently used parametrization approaches.

Refinement of empirical force fields for nucleic acids requires their extensive testing using as wide range of systems as possible. However, finding unambiguous reference data is not easy. In this paper, we analyze four systems which we suggest should be included in standard portfolio of molecules to test nucleic acids force fields, namely, parallel and antiparallel stranded DNA guanine quadruplex stems, RNA quadruplex stem, and Z-DNA. We highlight parameters that should be monitored to assess the force field performance. The work is primarily based on 8.4 μs of 100-250 ns trajectories analyzed in detail followed by 9.6 μs of additional selected back up trajectories that were monitored to verify that the results of the initial analyses are correct. Four versions of the Cornell et al. AMBER force field are tested, including an entirely new parmχ(OL4) variant with χ dihedral specifically reparametrized for DNA molecules containing syn nucleotides. We test also different water models and ion conditions. While improvement for DNA quadruplexes is visible, the force fields still do not fully represent the intricate Z-DNA backbone conformation.

Sequence-dependent local variations of helical parameters, structure, and flexibility are crucial for molecular recognition processes involving B-DNA. A-tracts, i.e., stretches of several consecutive adenines in one strand that are in phase with the DNA helical repeat, mediate significant DNA bending. During the past few decades, there have been intense efforts to understand the sequence dependence of helical parameters in DNA. Molecular dynamics (MD) simulations can provide valuable insights into the molecular mechanism behind the relationship between sequence and structure. However, although recent improvements in empirical force fields have helped to capture many sequence-dependent B-DNA properties, several problems remain, such as underestimation of the helical twist and suspected underestimation of the propeller twist in A-tracts. Here, we employ reference quantum mechanical (QM) calculations, explicit solvent MD, and bioinformatics to analyze the underestimation of propeller twisting of A-tracts in simulations. Although we did not identify a straightforward explanation, we discovered two imbalances in the empirical force fields. The first was overestimation of stacking interactions accompanied by underestimation of base-pairing energy, which we attribute to anisotropic polarizabilities that are not reflected by the isotropic force fields. This may lead to overstacking with potentially important consequences for MD simulations of nucleic acids. The second observed imbalance was steric clash between A(N1) and T(N3) nitrogens of AT base pairs in force-field descriptions, resulting in overestimation of the AT pair stretch in MD simulations. We also substantially extend the available set of benchmark estimated CCSD(T)/CBS data for B-DNA base stacking and provide a code that allows the generation of diverse base-stacking geometries suitable for QM computations with predefined intra- and interbase pair parameters.

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.

The folded structures of RNA molecules and large ribonucleoprotein particles are stabilized by a wide range of base pairs that actively utilize the 2'-OH groups of ribose for base pairing. Such base pairing does not occur in DNA and is essential for functional RNAs. We report reference quantum chemical calculations of base pairing energies for a representative selection of 25 RNA base pairs utilizing the ribose moiety for base pairing, including structures with amino acceptor interactions. All base pairs are evaluated at the MP2 level with extrapolation to the complete basis set (CBS) of atomic orbitals. CCSD(T) correction terms were obtained for four base pairs. In addition, the base pairing is evaluated using the DFT-SAPT perturbational procedure along with the aug-cc-pVDZ basis set, which allows for the decomposition of the interaction energies into separate, physically meaningful, components. These calculations confirm that, compared to canonical base pairs, many RNA base pairs exhibit a modestly increased role of dispersion attraction compared to canonical base pairs. However, the effect is smaller than one would assume based on assessment of the ratio of HF and correlation components of the interaction energies. Interaction energies are further calculated using the SCS(MI)-MP2 and DFT-D methods. Finally, we estimate the effect of aqueous solvent screening on the base pairing stability using the continuum solvent approach.

Terminal base pairs of DNA and RNA molecules in solution are known to undergo frequent transient opening events (fraying). Accurate modeling of this process is important because of its involvement in nucleic acid end recognition and enzymatic catalysis. In this article, we describe fraying in molecular dynamics simulations with the ff99bsc0, ff99bsc0χOL3, and ff99bsc0χOL4 force fields, both for DNA and RNA molecules. Comparison with the experiment showed that while some features of fraying are consistent with the available data, others indicate potential problems with the force field description. In particular, multiple noncanonical structures are formed at the ends of the DNA and RNA duplexes. Among them are tWC/sugar edge pair, C-H edge/Watson-Crick pair, and stacked geometries, in which the terminal bases are stacked above each other. These structures usually appear within the first tens to hundreds of nanoseconds and substantially limit the usefulness of the remaining part of the simulation due to geometry distortions that are transferred to several neighboring base pairs ("end effects"). We show that stability of the noncanonical structures in ff99bsc0 may be partly linked to inaccurate glycosidic (χ) torsion potentials that overstabilize the syn region and allow for rapid anti to syn transitions. The RNA refined glycosidic torsion potential χOL3 provides an improved description and substantially more stable MD simulations of RNA molecules. In the case of DNA, the χOL4 correction gives only partial improvement. None of the tested force fields provide a satisfactory description of the terminal regions, indicating that further improvement is needed to achieve realistic modeling of fraying in DNA and RNA molecules.

The importance of allosteric effects in DNA is becoming increasingly appreciated, but the underlying mechanisms remain poorly understood. In this work, we propose a general modeling framework to study DNA allostery. We describe DNA in a coarse-grained manner by intra-base pair and base pair step coordinates, complemented by groove widths. Quadratic deformation energy is assumed, yielding linear relations between the constraints and their effect. Model parameters are inferred from standard unrestrained, explicit-solvent molecular dynamics simulations of naked DNA. We applied the approach to study minor groove binding of diamidines and pyrrole-imidazole polyamides. The predicted DNA bending is in quantitative agreement with experiment and suggests that diamidine binding to the alternating TA sequence brings the DNA closer to the A-tract conformation, with potentially important functional consequences. The approach can be readily applied to other allosteric effects in DNA and generalized to model allostery in various molecular systems.

We present a brief overview of explicit solvent molecular dynamics (MD) simulations of nucleic acids. We explain physical chemistry limitations of the simulations, namely, the molecular mechanics (MM) force field (FF) approximation and limited time scale. Further, we discuss relations and differences between simulations and experiments, compare standard and enhanced sampling simulations, discuss the role of starting structures, comment on different versions of nucleic acid FFs, and relate MM computations with contemporary quantum chemistry. Despite its limitations, we show that MD is a powerful technique for studying the structural dynamics of nucleic acids with a fast growing potential that substantially complements experimental results and aids their interpretation.

Recent advances in polarizable force fields have revealed that major reparameterization is necessary when the polarization energy is treated explicitly. This study is focused on the torsional parameters, which are crucial for the accurate description of conformational equilibria in biomolecules. In particular, attention is paid to the influence of polarization on the (i) transferability of dihedral terms between molecules, (ii) transferability between different environments, and (iii) additivity of dihedral energies. To this end, three polarizable force fields based on the induced point dipole model designed for use in AMBER are tested, including two recent ff02 reparameterizations. Attention is paid to the contributions due to short range interactions (1-2, 1-3, and 1-4) within the four atoms defining the dihedral angle. The results show that when short range 1-2 and 1-3 polarization interactions are omitted, as for instance in ff02, the 1-4 polarization contribution is rather small and unlikely to improve the description of the torsional energy. Conversely, when screened 1-2 and 1-3 interactions are included, the polarization contribution is sizeable and shows potential to improve the transferability of parameters between different molecules and environments as well as the additivity of dihedral terms. However, to reproduce intramolecular polarization effects accurately, further fine-tuning of the short range damping of polarization is necessary. © 2015 Wiley Periodicals, Inc.

X-ray crystallography can provide important insights into the structure of RNA enzymes (ribozymes). However, the details of a ribozyme's active site architecture are often altered by the inactivating chemical modifications necessary to inhibit self-cleavage. Molecular dynamics (MD) simulations are able to complement crystallographic data and model the conformation of the ribozyme's active site in its native form. However, the performance of MD simulations is driven by the quality of the force field used. Force fields are primarily parametrized and tested for a description of canonical structures and thus may be less accurate for noncanonical RNA elements, including ribozyme catalytic cores. Here, we show that our recent reparametrization of ε/ζ torsions significantly improves the description of the hairpin ribozyme's scissile phosphate conformational behavior. In addition, we find that an imbalance in the force field description of the nonbonded interactions of the ribose 2'-OH contributes to the conformational behavior observed for the scissile phosphate in the presence of a deprotonated G8(-). On the basis of the new force field, we obtain a reactive conformation for the hairpin ribozyme active site that is consistent with the most recent mechanistic and structural data.

A-tracts are functionally important DNA sequences which induce helix bending and have peculiar structural properties. While A-tract structure has been qualitatively well characterized, their mechanical properties remain controversial. A-tracts appear structurally rigid and resist nucleosome formation, but seem flexible in DNA looping. In this work, we investigate mechanical properties of symmetric AnTn and asymmetric A2n tracts for n = 3, 4, 5 using two types of coarse-grained models. The first model represents DNA as an ensemble of interacting rigid bases with non-local quadratic deformation energy, the second one treats DNA as an anisotropically bendable and twistable elastic rod. Parameters for both models are inferred from microsecond long, atomic-resolution molecular dynamics simulations. We find that asymmetric A-tracts are more rigid than the control G/C-rich sequence in localized distortions relevant for nucleosome formation, but are more flexible in global bending and twisting relevant for looping. The symmetric tracts, in contrast, are more rigid than asymmetric tracts and the control, both locally and globally. Our results can reconcile the contradictory stiffness data on A-tracts and suggest symmetric A-tracts to be more efficient in nucleosome exclusion than the asymmetric ones. This would open a new possibility of gene expression manipulation using A-tracts.

The DNA sugar–phosphate backbone has a substantial influence on the DNA structural dynamics. Structural biology and bioinformatics studies revealed that the DNA backbone in experimental structures samples a wide range of distinct conformational substates, known as rotameric DNA backbone conformational families. Their correct description is essential for methods used to model nucleic acids and is known to be the Achilles heel of force field computations. In this study we report the benchmark database of MP2 calculations extrapolated to the complete basis set of atomic orbitals with aug-cc-pVTZ and aug-cc-pVQZ basis sets, MP2(T,Q), augmented by ΔCCSD(T)/aug-cc-pVDZ corrections. The calculations are performed in the gas phase as well as using a COSMO solvent model. This study includes a complete set of 18 established and biochemically most important families of DNA backbone conformations and several other salient conformations that we identified in experimental structures. We utilize an electronically sufficiently complete DNA sugar–phosphate–sugar (SPS) backbone model system truncated to prevent undesired intramolecular interactions. The calculations are then compared with other QM methods. The BLYP and TPSS functionals supplemented with Grimme's D3(BJ) dispersion term provide the best tradeoff between computational demands and accuracy and can be recommended for preliminary conformational searches as well as calculations on large model systems. Among the tested methods, the best agreement with the benchmark database has been obtained for the double-hybrid DSD-BLYP functional in combination with a quadruple-ζ basis set, which is, however, computationally very demanding. The new hybrid density functionals PW6B95-D3 and MPW1B95-D3 yield outstanding results and even slightly outperform the computationally more demanding PWPB95 double-hybrid functional. B3LYP-D3 is somewhat less accurate compared to the other hybrids. Extrapolated MP2(D,T) calculations are not as accurate as the less demanding DFT-D3 methods. Preliminary force field tests using several charge sets reveal an almost order of magnitude larger deviations from the reference QM data compared to modern DFT-D3, underlining the challenges facing force field simulations of nucleic acids. As expected, inclusion of the solvent environment approximated by a continuum approach has a large impact on the relative stabilities of different backbone substates and is important when comparing the QM data with structural bioinformatics and other experimental data.

Knowledge of geometrical and physico-chemical properties of the sugar–phosphate backbone substantially contributes to the comprehension of the structural dynamics, function and evolution of nucleic acids. We provide a side by side overview of structural biology/bioinformatics, quantum chemical and molecular mechanical/simulation studies of the nucleic acids backbone. We highlight main features, advantages and limitations of these techniques, with a special emphasis given to their synergy. The present status of the research is then illustrated by selected examples which include classification of DNA and RNA backbone families, benchmark structure–energy quantum chemical calculations, parameterization of the dihedral space of simulation force fields, incorporation of arsenate into DNA, sugar–phosphate backbone self-cleavage in small RNA enzymes, and intricate geometries of the backbone in recurrent RNA building blocks. Although not apparent from the current literature showing limited overlaps between the QM, simulation and bioinformatics studies of the nucleic acids backbone, there in fact should be a major cooperative interaction between these three approaches in studies of the sugar–phosphate backbone.

RNA molecules are stabilized by a wide range of noncanonical interactions that are not present in DNA. Among them, the recently classified base-phosphate (BPh) interactions belong to the most important ones. Twelve percent of nucleotides in the ribosomal crystal structures are involved in BPh interactions. BPh interactions are highly conserved and provide major constraints on RNA sequence evolution. Here we provide assessment of the energetics of BPh interactions using MP2 computations extrapolated to the complete basis set of atomic orbitals and corrected for higher-order electron correlation effects. The reference computations are compared with DFT-D and DFT-D3 approaches, the SAPT method, and the molecular mechanics force field. The computations, besides providing the basic benchmark for the BPh interactions, allow some refinements of the original classification, including identification of some potential doubly bonded BPh patterns. The reference computations are followed by analysis of some larger RNA fragments that consider the context of the BPh interactions. The computations demonstrate the complexity of interaction patterns utilizing the BPh interactions in real RNA structures. The BPh interactions are often involved in intricate interaction networks. We studied BPh interactions of protonated adenine that can contribute to catalysis of hairpin ribozyme, the key BPh interaction in the S-turn motif of the sarcin-ricin loop, which may predetermine the S-turn topology and complex BPh patterns from the glmS riboswitch. Finally, the structural stability of BPh interactions in explicit solvent molecular dynamics simulations is assessed. The simulations well preserve key BPh interactions and allow dissection of structurally/functionally important water-meditated BPh bridges, which could not be considered in earlier bioinformatics classification of BPh interactions.

The intermolecular interaction energy components for several molecular complexes were calculated using force fields available in the AMBER suite of programs and compared with Density Functional Theory-Symmetry Adapted Perturbation Theory (DFT-SAPT) values. The extent to which such comparison is meaningful is discussed. The comparability is shown to depend strongly on the intermolecular distance, which means that comparisons made at one distance only are of limited value. At large distances the coulombic and van der Waals 1/r(6) empirical terms correspond fairly well with the DFT-SAPT electrostatics and dispersion terms, respectively. At the onset of electronic overlap the empirical values deviate from the reference values considerably. However, the errors in the force fields tend to cancel out in a systematic manner at equilibrium distances. Thus, the overall performance of the force fields displays errors an order of magnitude smaller than those of the individual interaction energy components. The repulsive 1/r(12) component of the van der Waals expression seems to be responsible for a significant part of the deviation of the force field results from the reference values. We suggest that further improvement of the force fields for intermolecular interactions would require replacement of the nonphysical 1/r(12) term by an exponential function. Dispersion anisotropy and its effects are discussed. Our analysis is intended to show that although comparing the empirical and non-empirical interaction energy components is in general problematic, it might bring insights useful for the construction of new force fields. Our results are relevant to often performed force-field-based interaction energy decompositions.

Accurate representation of nucleic acids in molecular dynamics simulations depends critically on the quality of the applied empirical force field. Among force field terms, the torsion parameters are known to strongly influence the conformational equilibria and molecular structures. Unfortunately, past several years witnessed severe problems in describing the torsion space in nucleic acids by current force fields and more problems continue emerging. In an attempt for improvement, we suggested a novel parameterization procedure that incorporates some previously neglected solvation-related effects, which proved to be essential for obtaining accurate torsion profiles. The suggested approach avoids double counting of solvation effects and provides parameters that may be used in combination with any of the widely used nonpolarizable discrete solvent models or with the continuum solvent models. Improvements are demonstrated for the latest AMBER force field for RNA simulations, ff10, which incorporates parameters...

Dihedral angle parameters in force fields are often obtained with the help of quantum chemical calculations. These are usually performed in vacuum and it is assumed that the solvent model used in the molecular dynamics simulation (explicit or implicit) satisfactorily accounts for the major part of the solvation energy. However, there are certain nonspecific solvation effects that are not properly accounted for, in this approach. They relate to conformation-dependent solute polarization and solvation of conformation-dependent charge distribution. On an example of the glycosidic torsion, we show that the contribution resulting from these effects is substantial and provides important correction to the torsion potential. A parameterization procedure is suggested that incorporates the missing conformation-dependent solvation effects into the torsion parameters, based on the difference between the quantum mechanical self-consistent reaction field and Poisson–Boltzmann continuum solvation models. The suggested a...