Highly sampled tetranucleotide and tetraloop motifs enable evaluation of common RNA force fields

Abstract

Recent modifications and improvements to standard nucleic acid force fields have attempted to fix problems and issues that have been observed as longer timescale simulations have become routine. Although previous work has shown the ability to fold the UUCG stem-loop structure, until now no group has attempted to quantify the performance of current force fields using highly converged structural populations of the tetraloop conformational ensemble. In this study, we report the use of multiple independent sets of multidimensional replica exchange molecular dynamics (M-REMD) simulations with different initial conditions to generate well-converged conformational ensembles for the tetranucleotides r(GACC) and r(CCCC), as well as the larger UUCG tetraloop motif. By generating what is to our knowledge the most complete RNA structure ensembles reported to date for these systems, we remove the coupling between force field errors and errors due to incomplete sampling, providing a comprehensive comparison between current top-performing MD force fields for RNA. Of the RNA force fields tested in this study, none demonstrate the ability to correctly identify the most thermodynamically stable structure for all three systems. We discuss the deficiencies present in each potential function and suggest areas where improvements can be made. The results imply that although "short" (nsec-μsec timescale) simulations may stay close to their respective experimental structures and may well reproduce experimental observables, inevitably the current force fields will populate alternative incorrect structures that are more stable than those observed via experiment.

Keywords: AMBER; CHARMM; RNA; enhanced sampling; force fields; molecular dynamics; replica exchange.

FIGURE 1.

Images of RNA systems simulated in this study include (A) the r(CCCC) tetranucleotide, shown in an A-form conformation with bases colored by residue number, (B) the r(GACC) tetranucleotide, shown in the NMR major conformation with bases colored by residue number, (C) the r(GGCAC-UUCG-GUGCC) 2KOC NMR ensemble, where the A:U base pair is shown in green, G:C base pairs are shown in light blue, and the UUCG loop sequence is shown in dark blue, (D) r(CAC-UUCG-GUG), a truncated version of structure 1 from the 2KOC NMR ensemble, used as a starting structure for M-REMD simulations, and (E) The UUCG tetraloop with alternate stem sequence, used as a starting structure for Anton simulations. The G:C base pairs are shown in light blue and the UUCG loop sequence in dark blue. We note that an alternative stem structure was used in the Anton simulations (in the 2010–2011 timeframe) to be consistent with our earlier tetraloop investigations. In this earlier study, representatives of all tetraloop structures found in the PDB were grafted onto a common stem structure; these relatively short MD simulations with even more poorly performing force fields were never published.

FIGURE 2.

Mass-weighted loop RMSD (UUCG loop and closing base pair are included) versus time. (A) Twenty 200-nsec simulations of UUCG, starting from each of the 20 members of the NMR ensemble. RMSDs are shown with all data in gray, and 1-nsec running averages of each simulation in black (if they maintain native contacts) or in color (if native loop contacts are lost). (B) Three >1-μsec simulations of UUCG, performed on a combination of CPUs and GPUs, starting from the first three members of the NMR ensemble. (C) Two >5-μsec simulations performed on Anton. Increases in the RMSD can be attributed to a loss of loop structure.

FIGURE 3.

(Left) M-REMD RMSD to A-form Reference for r(GACC). (A) ff12, (B) ff12 + vdW_all, (C) ff12 + vdW_bb, (D) ff99 + Chen–Garcia, (E) ff99 + χ_Yil, and (F) C36 force fields. The averages between two runs per force field are shown, with error bars shown as the standard deviation. RMSDs corresponding to NMR major and minor structures are ∼2.0 Å and 2.6 Å, respectively. (G) Populations of major conformation types seen in cluster analysis of each M-REMD simulation. Bars indicate the average and standard deviation between two independent runs and are colored by force field to match the RMSD histograms. Experimental values are shown on the table in a blue striped pattern. Representative structures and classifications (including discussion of the conformations classified as “Other”) are shown in Supplemental Table 1 and Supplemental Figure 2.

FIGURE 4.

M-REMD RMSD to A-form Reference for r(CCCC). (A) M-REMD RMSD histogram profiles for r(CCCC) in five force fields. The averages of two runs per force field are shown, with error bars as standard deviation between runs. (B) Close-up of AMBER force fields, 4.8–5.7 Å. (C) A-form reference structure. (D) Top two clusters from combined cluster analysis and their populations in each force field, colored to match the force field designations from the top histogram plot. The mass-weighted RMSD to the A-form reference is provided for each structure.

FIGURE 5.

Combined clustering results for UUCG tetraloop simulations. Cluster number is shown on the x-axis and percentage of each simulation's 277 K ensemble is shown on the y-axis. Representative structures are shown below the bar graph. The native structure found in cluster 5 is shown overlapped with the NMR structure (in green).

FIGURE 6.

Free energy versus χ dihedral angle distribution at 277 K for the A:U stem base pair (top) and G_L4 residue (bottom) for four force fields. Syn minimum is 60°, anti minimum is 180–240°, and high anti is 270°. Data are the average of two independent runs with error bars representing standard deviation.

FIGURE 7.

Overlap of Anton Run 2 final structure (metallic red) with representative structure from Cluster 3 of combined cluster analysis at 277 K (blue stem and colored UUCG loop bases). All atom RMSD is 0.7 Å.