Performance of Molecular Mechanics Force Fields for RNA Simulations: Stability of UUCG and GNRA Hairpins

Abstract

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.

Figure 1.

(Left) Secondary structures of the studied systems with base pairing and base–phosphate interactions annotated according to the standard classifications.^, G_L4 of the UUCG tetraloop having syn orientation is highlighted in red. The modeled GC pairs in the GAGA system are shown in gray. The loop residues are labeled as L1–L4 to avoid context numbering. For instance, U6 of UUCG is labeled as U_L1. (Right) Three-dimensional structures of studied systems. The A-RNA stem part is shown in red, while the tetraloop nucleotides are in blue.

Figure 2.

Signature H-bonds (black dashed lines) of UUCG and GAAA tetraloops on the left and right sides, respectively.

Figure 3.

The MD snapshots (colored by atom types) compared with high-resolution NMR structure (in red) showing structural problems seen in simulations of the UUCG tetraloop (some atoms are not shown for clarity, and important parts are shown as sticks). (A) The disruption of the U_L1(O2′)•••G_L4(O6) H-bond and formation of a new U_L1(O2′)•••U_L2(O5′) H-bond observed in all MD simulations with standard χ profiles are highlighted by the blue arrow, while the simultaneous decrease of χ of G_L4 leading to a change in the U_L1/G_L4 propeller is shown by the red arrow. (B) The U_L2 phosphate α/γ flip is depicted by the black arrow. See the text for full details.

Figure 4.

Structures of UUCG TL at the beginning of the UUCG_bsc0 simulation, at 50 ns and at 90 ns, showing the distortion of the UUCG tetraloop. C_L3 is not shown, for clarity.

Figure 5.

Ladder-like conformer as observed in a simulation of the GNRA tetraloop with AMBER force fields, unless the χ torsion profile is appropriately modified. Initial geometry is on the left, and “ladder-like” conformer is on the right.

Figure 6.

Typical progression of GNRA tetraloop AMBER simulations. Left: simulations without the χ correction or with Ode et al.’s correction illustrated by a 25 ns GAGA_bsc0 simulation. Right: simulations with χ_OL-DFT, χ_OL, and χ_YIL variants illustrated by a 1 μs GAGA_bsc0χ_OL simulation. The upper graphs present a time evolution of the β(G_L3) torsion (black line), γ(G_L3) torsion (red line), and mean χ torsion averaged over either stem nucleobases (green line) or the GNRA tetraloop (blue line). The middle graph shows the G_L1/A_L4 buckle, and the lower graph presents the GNRA tetraloop signature H-bonds: G_L1(N2)•••A_L4(N7), G_L1(N2)•••A_L4(pro-R_p), and G_L1(O2′)•••G_L3(N7) in black, red, and green, respectively.