Molecular dynamics and quantum mechanics of RNA: conformational and chemical change we can believe in.

Structure and dynamics are both critical to RNA's vital functions in biology. Numerous techniques can elucidate the structural dynamics of RNA, but computational approaches based on experimental data arguably hold the promise of providing the most detail. In this Account, we highlight areas wherein molecular dynamics (MD) and quantum mechanical (QM) techniques are applied to RNA, particularly in relation to complementary experimental studies. We have expanded on atomic-resolution crystal structures of RNAs in functionally relevant states by applying explicit solvent MD simulations to explore their dynamics and conformational changes on the submicrosecond time scale. MD relies on simplified atomistic, pairwise additive interaction potentials (force fields). Because of limited sampling, due to the finite accessible simulation time scale and the approximated force field, high-quality starting structures are required. Despite their imperfection, we find that currently available force fields empower MD to provide meaningful and predictive information on RNA dynamics around a crystallographically defined energy minimum. The performance of force fields can be estimated by precise QM calculations on small model systems. Such calculations agree reasonably well with the Cornell et al. AMBER force field, particularly for stacking and hydrogen-bonding interactions. A final verification of any force field is accomplished by simulations of complex nucleic acid structures. The performance of the Cornell et al. AMBER force field generally corresponds well with and augments experimental data, but one notable exception could be the capping loops of double-helical stems. In addition, the performance of pairwise additive force fields is obviously unsatisfactory for inclusion of divalent cations, because their interactions lead to major polarization and charge-transfer effects neglected by the force field. Neglect of polarization also limits, albeit to a lesser extent, the description accuracy of other contributions, such as interactions with monovalent ions, conformational flexibility of the anionic sugar-phosphate backbone, hydrogen bonding, and solute polarization by solvent. Still, despite limitations, MD simulations are a valid tool for analyzing the structural dynamics of existing experimental structures. Careful analysis of MD simulations can identify problematic aspects of an experimental RNA structure, unveil structural characteristics masked by experimental constraints, reveal functionally significant stochastic fluctuations, evaluate the structural role of base ionization, and predict structurally and potentially functionally important details of the solvent behavior, including the presence of tightly bound water molecules. Moreover, combining classical MD simulations with QM calculations in hybrid QM/MM approaches helps in the assessment of the plausibility of chemical mechanisms of catalytic RNAs (ribozymes). In contrast, the reliable prediction of structure from sequence information is beyond the applicability of MD tools. The ultimate utility of computational studies in understanding RNA function thus requires that the results are neither blindly accepted nor flatly rejected, but rather considered in the context of all available experimental data, with great care given to assessing limitations through the available starting structures, force field approximations, and sampling limitations. The examples given in this Account showcase how the judicious use of basic MD simulations has already served as a powerful tool to help evaluate the role of structural dynamics in biological function of RNA.

Why Use Computational Techniques on RNA?

The central role of RNA in numerous biological processes including translation,(1) protein localization,(2) gene regulation,(3) RNA processing,(4) and viral replication(5) calls for a detailed understanding of RNA function, structure, and conformational dynamics.(6) Accompanying and enhancing our increasing appreciation of RNA is the rapidly expanding availability of high-resolution structures of RNAs and RNA−protein (RNP) complexes. These atomic resolution snapshots provide detailed rationalization for existing biochemical data; however, biological function depends on the dynamic evolution of structures along functional pathways. A complete understanding of the relevant structural dynamics exhibited by RNA requires monitoring time scales from picoseconds to hours through the application of a correspondingly broad range of techniques (Figure 1),(6) with careful consideration given to the scope and limitation of each approach.

Figure 1

System size and time scales of techniques commonly used to evaluate RNA dynamics.

Provided they are judiciously applied, computational methods provide insights that are not fully accessible through experimental techniques. Reproduction of experimental data is desirable for assessing accuracy, but modern in silico approaches have matured beyond this modest initial goal. Molecular dynamics (MD) simulations can identify problematic aspects of experimental structures,[7−10] reveal functionally significant stochastic fluctuations,[11,12] predict the impact of base ionization on RNA structure and dynamics,[10,13] and characterize solvent behavior.[7,10,14,15] Combining simulations with quantum mechanical (QM) calculations in QM/MM approaches expands the repertoire of applications to mechanistic questions concerning the reaction chemistry of catalytic RNAs, or ribozymes.[16−20] The primary goal of this Account is to highlight areas of applicability of MD and QM techniques to RNA and their relation to complementary experimental studies.

What General Scope and Limitations Do MD Simulations Have?

Explicit solvent MD is an atomistic technique dealing with one or more primary solute molecules of defined geometry surrounded by an environment of water and ions that jointly undergo 1 to >1000 ns of dynamics simulated at room temperature. MD provides an unsurpassed level of detail of all aspects of the time evolution (with subpicosecond time resolution) of the three-dimensional structure, including the positions of all water molecules and ions. However, MD simulations are faced with two significant limitations. First, the sampling of conformational space reflects the short time scale of MD compared with actual biological processes. This limitation is slowly being overcome with the continuing emergence of more powerful computers and codes. Second, a fundamental limitation not waning with faster computers is the approximate nature of force fields, which are simple analytic atomistic functions relating structure with potential energy.[21,22] Despite sophisticated parametrization, the force field simply consists of sets of harmonic springs for both bond lengths and valence angles, supplemented by torsion profiles for dihedral angles.[21−24] Atoms are approximated as Lennard-Jones spheres with constant point charges localized at the atomic centers.[23,25]

There are two widely used nucleic acid (NA) force fields, Cornell et al.(21) and CHARMM27,(22) which share similar functional form but differ in parametrization.(26) Variants of the Cornell et al. (AMBER) force field (parm94-99, bsc0) have been extensively tested for many folded RNAs and noncanonical DNAs.[7−9,12−15,24,27−32] CHARMM27 describes B-DNA well(33) but has not yet been systematically tested for either folded RNAs or noncanonical DNAs. Both AMBER and CHARMM offer high-quality protein force fields for consistent description of NA−protein complexes.

The Cornell et al. parametrizations describe the key electrostatic terms using atomic charges derived to reproduce the electrostatic potential around the NA building blocks.(23) QM calculations show that base stacking is the best approximated term in NA simulations,(34) followed by base pairing, including non-Watson−Crick interactions utilizing the 2′-OH group.[28,35] The description of the flexible backbone is less straightforward.(24) In particular, the phosphate group is highly polarizable and the many dihedral angles of the backbone each possess multiple substates. The backbone description would therefore benefit from geometry-dependent electrostatic and polarization terms. Monovalent ions and solute−solvent interactions are thought to be reasonably well described, while the description of divalent ions is outside the applicability of force fields. Ions are simplified as Lennard-Jones spheres with constant point charges in their center. In reality, the first ligand shell of divalent metal ions is highly activated through polarization and charge-transfer so that constituent water molecules form very strong hydrogen bonds that propagate these effects beyond the first ligand shell.(36) All these contributions are neglected by the force field.

The quality of a force field’s performance inherently relies on the mutual compensation of errors, which in turn depends on a balance between forces in the system under study and the accuracy of the parametrization. There are two basic possibilities: (i) The compensation of errors is sufficient, and the force field finds the correct global minimum of the simulated system.(32) In this case not necessarily all details are correct, but the overall description is meaningful. The more qualitative the computational task, the more likely the force field description is sufficient. (ii) The force field does not give the correct global minimum and the simulated system eventually degrades.[24,37]

Biomolecular force fields are intentionally parametrized as multipurpose, with delicate trade-offs in parametrization. Tremendously challenging efforts are being expended to develop more physically accurate multipurpose polarization force fields.[38−43] A major problem in parametrization of these sophisticated biomolecular force fields may be to achieve a satisfactory overall balance among all their parameters.

Local conformational traps associated with standard MD can be overcome by enhanced sampling techniques (locally enhanced sampling, replica exchange, or targeted MD).[27,37] Robust sampling is also critical to obtain reliable results from free energy calculations that can provide useful information on RNA conformations,[17−20,27,44,45] but this goes beyond the scope of this Account.

What General Scope and Limitations Do Quantum Mechanical Calculations Have?

In contrast to force fields, QM can achieve a physically more rigorous description of chemical systems. Ab initio QM methods are free of empirical parameters and offer a systematic (and controllable) tuning of their quality by improving the underlying basis sets of atomic orbitals together with a balanced inclusion of electronic correlation effects. Accurate QM calculations are, however, currently limited to ∼30−50 atoms and are carried out in the gas phase.(27) QM allows for reliable evaluation of intrinsic (gas-phase) interaction energies, defined as differences between the electronic energies of a dimer and its component noninteracting monomers. This direct structure−energy relationship can be accurately calculated for any single geometry of a stacking or base-pairing interaction to map the complete potential energy surface.[28,34,35] Such energies unambiguously reflect direct forces between the interacting partners with no influence of the environment, making QM a genuine reference tool to parametrize and verify other computational approaches.[24,28,35] When electron correlation calculations are expanded to complete basis sets of atomic orbitals and include corrections for higher-level electron correlation effects, QM calculations reach similar accuracies for both base pairing and stacking and are effectively converged.[34,35] However, it is not straightforward to extend QM calculations to biomolecules. NA conformations in particular result from a highly variable mixture of mutually compensating interactions, the balance of which may vary for distinct architectures. In addition, the strong electrostatic forces in NA are substantially modulated by solvent screening effects. Accurate inclusion of solvent effects is beyond the capability of modern QM approaches. Special care is needed when including the NA backbone in QM studies. Isolated small model systems (even as small as a single nucleotide) favor geometries that are biased by gas-phase specific intramolecular hydrogen bonds, where electrostatic effects of the phosphates dominate the energetics.(46) The problem is not so much the quality of the QM methods, but the incompleteness of the model system.

What Practical Considerations Arise during Implementation of MD ?

MD Relies on the Availability of Accurate High-Resolution Structures

If a reasonably accurate starting structure is available, MD in many cases can locally improve molecular interactions and backbone conformations in the experimental structure.(47) Due to force field and sampling limitations MD is unable to predict RNA structures without experimental input.(48) Additionally, unrealistic models often become swiftly distorted during a simulation, revealing their inadequacies. If the starting structure is in an incorrect conformation confined by large (>5 kcal/mol) energy barriers, an MD simulation cannot easily move it away from its starting geometry.(29) Reliable characterization of the molecular dynamics of an RNA therefore requires the use of high-resolution experimental starting structures.

Simulation Performance Needs To Be Judged by Comparison with Experimental Data

After collecting enough simulations (current state-of-the-art is multiple simulations with ∼20−50 ns duration each),(29) a careful comparison of the MD time trajectories with the experimental structure is needed. A simplified analysis of only a few heteroatom distances of interest accompanied by generally uninformative root-mean-square distance (rmsd) plots may mask considerable problems. The structural evolution results from a mixture of factors, including the actual stochastic flexibility of the RNA, experimental artifacts introduced through crystal contacts,(8) disorder or chemical modification, and force field artifacts. If this mixture is properly resolved the analysis of MD simulations can be very insightful.

Evaluation of RNA backbone conformations is difficult for both computational and experimental approaches. The flexibility and polarizability of the backbone is challenging for nonpolarizable force fields based on constant point charges. However, we usually observe surprisingly good agreement between experiment and RNA simulations due to compensating errors, base-pairing and stacking constraints, and the accuracy of the starting structures. For example, the dihedral backbone angles around an S-turn motif in a 2.05 Å resolution structure of the hairpin ribozyme with a single-nucleotide U39C mutation(49) differ from those of lower resolution (2.65 Å) structures carrying wild-type U39.(50) This difference could be due to either the distinct crystallography constructs used in the two studies or an artifact of the more limited resolution of the second structure. MD simulations resolved this ambiguity.(10) Starting from the lower resolution crystal structure, the backbone dihedrals switched to those observed in the higher resolution structure (Figure 2A). Structural bioinformatics using a recently developed backbone nomenclature assisted in the rapid evaluation of the backbone behavior.[10,51]

Figure 2

The good and the bad of MD simulations. (A) Overlay of the S-turn from a lower resolution structure (color) with a higher resolution structure (gray) before and after MD simulation. The backbone trace is highlighted (thick sticks). The dihedral angles of the backbone become more like those observed in the higher resolution structure. The angle of a δ dihedral is indicated (spheres) as an example of a large improvement in angle following MD. The change in the backbone suites of U39 following simulation are indicated in the figure, where !! signifies an atypical backbone conformation. (B) Apparent failure of the force field in describing guanine quadruplex DNA. The stem−loop shifts away from the experimental starting conformation during simulations, including loss of a critical potassium ion (spheres).

An example of force field problems arises in simulations of guanine quadruplex DNA (G-DNA) that consists of four-stranded stems formed by cation-stabilized guanine quartets complemented by single-stranded hairpin loops. The parm99 force field provides a global minimum consistent with the experimental structures for the G-DNA stems but not for the loops (Figure 2B).(37) Thus, in a given simulation different parts of a molecule can be described with different success. Unexpectedly, problems also arose in long (15 to >50 ns) MD simulations of B-DNA, due to an accumulation of irreversible, experimentally unobserved backbone substates with concomitant progressive degradation of the whole structure. The γ backbone torsional profile was subsequently reparametrized (parmbsc0).(24) This improved force field allows stable microsecond time scale simulations of B-DNA and even repairs partially degraded B-DNA structures, indicating that B-DNA is now the global minimum.(32) The two examples above are the worst known problems of the Cornell et al. force fields. For most other systems, we do not encounter inaccuracies of this magnitude. Considering its simplicity, the performance of the AMBER Cornell et al. force fields for RNA so far has been remarkable. It is nevertheless important to note that the force field is still, even after recent adjustment,(24) not satisfactory for G-DNA loops. Similarly, replica exchange MD (REMD) reported basic folding of an RNA stem−loop; however, the functional tetraloop structure was not sampled, likely due to both REMD and force field limitations.(52)

What Are Specific Questions That MD simulations of RNA Can Currently Address?

Resolving Experimental Artifacts

Although simulations cannot predict the overall folding of an RNA, they can locally resolve regions of limited resolution in known experimental structures and reveal structural defects due to crystal packing. For example, a local region of lower resolution is observed near the conformationally dynamic active site in precursor crystal structures of the hepatitis delta virus (HDV) ribozyme. An unusual set of backbone dihedral angles at the active site become more canonical during MD simulations, ultimately adopting a common U-turn motif (Figure 3A).(9) Crystal structures of the HDV ribozyme also show an extruded guanine (G76) that participates in crystal packing.(53) Multiple simulations reveal the rapid loss of this particular conformation and suggest a possible role of G76 in promoting catalysis through novel hydrogen-bonding interactions with stem P1.(8) Similarly, inactivating 2′-deoxy or 2′-O-methyl backbone modifications or base mutations used to trap ribozymes in their precatalytic structures may distort the active site. Multiple MD simulations of the hairpin ribozyme consistently result in a change in the A-1 sugar pucker in the absence of the 2′-O-methyl modification present in the experimental structure, leading to significant repositioning of the catalytically important nucleotides G8 and A38 (Figure 3B).(10)

Figure 3

Improving on experimental RNA backbone structures. (A) Rotation of dihedrals toward values consistent with the U-turn motif from the hammerhead ribozyme (gray) occurs over the course of MD simulations of the HDV ribozyme (color), with a pronounced rotation of the γ dihedral (spheres). (B) An active site conformational change in the hairpin ribozyme during MD is a consequence of the removal of a catalytically inactivating 2′O-methyl modification necessary for crystallization. A green arrow indicates the cleavage site.

Flexibility of RNA Building Blocks

Stochastic flexibility is a key functional feature of RNA that is difficult to derive from experiment. MD fills the gap by achieving a qualitative, atomistic understanding of the stochastic dynamics and flexibility of RNA building blocks.[11,12] For example, simulations have revealed striking intrinsic dynamics of RNA kink-turns, which can act as anisotropic molecular elbows to facilitate functional dynamics of the ribosome (Figure 4).(11) The dynamics of the helix 42−44 segment of 23S rRNA, part of the GTPase associated center, is substantially nonharmonic and is not satisfactorily captured by a coarse-grained normal-mode analysis.(12)

Figure 4

MD can capture the intrinsic flexibility of rRNA building blocks. (A) Position of the helix (H) 42−44 segment (gray) within the ribosome (left and right shown with and without proteins, respectively). Kink-turn (Kt) 42 resides in the center of H42, adjacent to the GTPase associated center RNA (rGAC, H43−44). (B) Range of spontaneous thermal fluctuations sampled during simulations. (C) While Kt-42 is a genuine elbow (hinge) element, its flanking stems are relatively stiff. Hinge-like dynamics can also occur between the NC-stem and rGAC.

Revealing Solvent Dynamics

MD simulations can detect long-residency water molecules that occupy prominent hydration sites and stay bound for many nanoseconds, contrasting with common water binding events of only ∼50−500 ps. Long-residency water molecules can serve structural, functional, and possibly catalytic roles.[7,11,14,54] A structurally relevant long-residency hydration site was detected in the A-minor I tertiary interactions of kink-turns 38 and 42 in 23S rRNA. Their A:C base pairs dynamically oscillate between a direct and water-mediated hydrogen bond whose interconversion significantly contributes to the elbow-like flexibility of the kink-turns. The static X-ray structures show both geometries, where the A:C interaction of kink-turn 38 is water-mediated and that of kink-turn 42 is direct.(11) Simulations based on lower resolution crystal structures of the hairpin ribozyme predicted the presence of interdomain long residency water molecules that were ultimately verified by the emergence of higher resolution structures.[7,55]

MD simulations can qualitatively characterize major binding sites of monovalent ions that are primarily determined by electrostatic interactions. Simulations in monovalents alone have revealed ion density in known multivalent ion binding sites.[7,13−15,29,30,56] For example, simulations of the HDV ribozyme predict monovalent cation binding at the cleavage site in a crystallographically resolved divalent metal ion binding site proximal to the 5′-O leaving group. Two Na+ ions and their accompanying first hydration spheres fill the catalytic pocket and may contribute to catalysis in the absence of divalents.(15) This MD prediction was later verified by crystal structures solved in the presence of Tl+ that reveal two Tl+ ions at the active site.(53) Simulations further predict a competition between ion binding and protonation of C75,(15) a feature not evident from the crystal structures but supported by mechanistic studies.(57) An additional binding site was predicted near the 2′-OH nucleophile. This site is again verified by crystal structures; however the exact coordination geometry differs between experiment and simulation, likely reflecting a combination of differences between the ions used, force field approximations, and crystallographic ambiguities.(15) Ion binding sites may also elude experimental detection, due to either low resolution or ion delocalization in the pocket, as observed for 5S rRNA loop E and the human immunodeficiency virus (HIV) dimerization initiation signal (DIS) kissing complex.[14,29,56]

Probing the Structural Effects of Base Substitutions and Ionizations

MD can assess the effects of base substitutions at atomic resolution to complement experimental mutagenesis studies. Thus, several base substitutions were modeled into the experimentally determined crystal structures of the hairpin ribozyme,(7) with good agreement between the MD predicted and experimentally determined stability of the tertiary structure.(7) The same simulations also revealed the importance of coupled networks of hydrogen bonds involving long residency water molecules for tertiary structure stability, whereby mutations exert significant long-range effects.(7) In simulations of the HDV ribozyme, each of the four standard nucleobases was separately modeled into the position immediately 5′ of the cleavage site. The wild-type U-1 was found to have the most tightly folded catalytic core, consistent with experimental footprinting data.(9) The same simulations revealed that a hydrogen bond characteristic of U-turn motifs from U1 to the phosphate of C3 is only transiently sampled (Figure 3A), reflecting a local flexibility at the cleavage site that correlates with increased catalytic activity. Simulations of ribozymes thus reveal additional structural and functional features that expand on experimental structures.

MD simulations of the HDV and hairpin ribozymes were used to assess the impact of protonation state on catalytically relevant structures. In simulations of the HDV ribozyme in which C75 is in its neutral (unprotonated) form, a geometry is adopted that is suitable for general base catalysis,[13,16] while simulations with a protonated C75H+ do not predict a reasonable geometry for C75 to act as a general acid.(13) (Note that this observation may indicate that the crystal structure fails to capture the catalytically relevant geometry.) For hairpin ribozyme catalysis, compelling mechanistic evidence likewise suggests a direct catalytic role for A38. Similar to the HDV ribozyme, MD simulations of the hairpin ribozyme with an unprotonated A38 lead to a geometry compatible with A38 acting as a general base.(10) Such a role is consistent with the available biochemical evidence but was previously discounted based on heteroatom distances likely influenced by the backbone 2′-O-methyl modification of the cleavage site in crystal structures (Figure 3B).[50,55] Unlike the HDV ribozyme, hairpin ribozyme simulations with a protonated A38H+ provide a geometry suitable for general acid catalysis as well.(10) Obviously, assessing a catalytic mechanism is limited in classical MD since bond breaking and making are by definition unobservable, thus warranting the use of QM to further evaluate the feasibility of a specific mechanism suggested by classical MD and biochemical data.

What Can QM/MM Reveal about the Chemical Change Catalyzed by Ribozymes?

A wide spectrum of both fast and accurate QM approaches has recently emerged, allowing for the inclusion of hundreds of atoms in a QM calculation.(58) Unfortunately, making a QM system larger but still incomplete will only exacerbate the errors resulting from the incompleteness of the system.(46) However, fast QM methods facilitate applications of QM/MM hybrid methods where a smaller segment of the system is treated quantum chemically while the remainder, including the solvent, is treated classically using force fields.(59) QM is particularly attractive for ribozymes, since QM, but not MD can describe the catalyzed reactions.(60) The main limitations of current QM/MM methods derive from insufficient sampling, inaccuracies of the QM or MM method, and treatment of the QM/MM boundary. To enhance QM/MM sampling, semiempirical (such as AM1, SCC-DFTB) and empirical (EVB) methods are used.[59,60]

Calculations using QM/MM methods have predicted specific roles for nucleobases, divalent ions, or electrostatic stabilization in catalyzing self-cleavage by the hammerhead, hairpin, and HDV ribozymes.[16−18] QM/MM methods have also been successfully applied to the elucidation of the mechanism of peptide bond formation and translation termination on the ribosome.[20,61] MD simulations of the HDV ribozyme provided a suitable starting geometry for a mechanism in which an unprotonated neutral C75 acts as a general base. QM/MM calculations are consistent with a role of C75 as the general base and Mg2+ as the general acid, predicting an energy barrier of ∼20 kcal/mol for the catalyzed reaction, in good agreement with the available experimental data.(16) For the QM scans, a region made up of 80 atoms in the active site was treated quantum-mechanically. Multiple starting positions of a specifically bound Mg2+ were sampled, establishing a hexacordinated Mg2+ ion with a single inner-sphere contact to a cleavage site nonbridging oxygen as the most likely conformation, with the Mg2+ acting as a Lewis acid in the reaction (Figure 5). Mechanisms in which C75 acts as the general acid instead, suggested by some recent biochemical studies,[62,63] could not be explored due to the paucity of suitable starting geometries. In contrast, MD simulations of the hairpin ribozyme with protonated and unprotonated A38 result in plausible catalytic geometries for A38 acting as general acid or general base, respectively.(10) These simulations reveal in part the complex impact of base ionization on the starting ground-state geometry and may explain the apparent insensitivity to base ionization of an initial QM/MM analysis of the ribozyme based on a crystal structure of a transition-state analog.[17,18] Large-scale classical simulations may be essential in establishing starting geometries for subsequent QM/MM calculations.

Figure 5

Mechanism of the HDV ribozyme. Starting from an experimental structure, the unprotonated and protonated forms of C75 were used in simulations to derive reasonable starting structures for mechanisms in which C75 serves as either a general base or a general acid, respectively. MD simulations result in a reasonable geometry for the general base but not the general acid mechanism. A range of Mg2+ ion positions were sampled, and the most likely position was determined. QM/MM predicts an energy barrier for the reaction that is in reasonable agreement with the available experimental data.

Conclusions

MD simulations of RNA are a powerful tool to expand on experimental structures and biochemical data, providing unique atomistic descriptions of the dynamic roles of nucleobases, the backbone, counterions, and individual water molecules in imparting biological function to RNA. Experiments benefit from a side-by-side comparison with simulations, where MD can serve to refine, interpret, and better understand existing experimental structures. In evaluating MD simulations, it is essential to consider that ensemble averaging and error margins of the underlying experimental structures have an impact and that force field artifacts are pervasive. In some instances, the available force field may not be sufficient to obtain meaningful results, in which case the limitations should be fully acknowledged(37) or even be addressed by improving the MD method.(24) QM calculations, often in the form of hybrid QM/MM approaches, can further build on MD simulations to access reaction chemistry.