Multidimensional Replica Exchange Molecular Dynamics Yields a Converged Ensemble of an RNA Tetranucleotide

Abstract

A necessary step to properly assess and validate the performance of force fields for biomolecules is to exhaustively sample the accessible conformational space, which is challenging for large RNA structures. Given questions regarding the reliability of modeling RNA structure and dynamics with current methods, we have begun to use RNA tetranucleotides to evaluate force fields. These systems, though small, display considerable conformational variability and complete sampling with standard simulation methods remains challenging. Here we compare and discuss the performance of known variations of replica exchange molecular dynamics (REMD) methods, specifically temperature REMD (T-REMD), Hamiltonian REMD (H-REMD), and multidimensional REMD (M-REMD) methods, which have been implemented in Amber's accelerated GPU code. Using two independent simulations, we show that M-REMD not only makes very efficient use of emerging large-scale GPU clusters, like Blue Waters at the University of Illinois, but also is critically important in generating the converged ensemble more efficiently than either T-REMD or H-REMD. With 57.6 μs aggregate sampling of a conformational ensemble with M-REMD methods, the populations can be compared to NMR data to evaluate force field reliability and further understand how putative changes to the force field may alter populations to be in more consistent agreement with experiment.

Figure 1

Schematic diagram of the M-REMD simulation. Each plane represents a different Hamiltonian, and each arrow represents a different temperature. The unbiased Hamiltonian, shown in yellow, can be recovered at the Temperature of interest. Exchange in temperature dimension between 285.7 and 290.2 K is followed by an exchange in Hamiltonian space between H₀ and H₁. This figure represents four out of 24 temperatures and three out of eight Hamiltonian dimensions used in this work.

Figure 2

Population analysis showing the number of structures at specific RMSD values from an A-form reference structure. Mass weighted RMSD histograms of the unbiased replicas at 300 K are averaged (shown in black) between two runs (shown in red and blue). Error bars represent standard deviation between two runs. Independent sets of starting structures were used for each simulation.

Figure 3

Representative structures from DBscan cluster analysis at 300 K. The top six most populated clusters in most REMD simulations are shown above. For the r(GACC) sequence the coloring is G1 (red), A2 (green), C3 (cyan), and C4 (magenta). Cluster populations for M-REMD are shown in black, H-REMD in red, and T-REMD in blue. M-REMD and H-REMD cluster populations represent the average of two independent simulations. Cluster names were assigned on the basis of analysis of the 277 K T-REMD trajectory.

Figure 4

Correlation of cluster populations between independent H-REMD runs and independent M-REMD runs. Linear fit and correlation coefficients for H-REMD and M-REMD are shown, where red fits all clusters. The shaded area represents a 95% confidence interval where the slope and y-intercept bounds are denoted by α and β, respectively.

Figure 5

Principal component projection of top five modes onto (left) H-REMD run 1 (solid lines) and H-REMD run 2 (dashed lines) show little overlap in the low frequency modes, whereas (right) M-REMD run 1 and M-REMD run 2 show the overlap in the low frequency modes.

Figure 6

Kullback–Leibler divergence analysis. Time dependent Kullback–Leibler divergence. Left: first three principal components from the H-REMD and M-REMD simulations. H-REMD principal components 1, 2, and 3 are maroon, red, and orange. M-REMD principal components 1, 2, and 3 are blue, purple, and cyan. Right: mass-weighted heavy atom RMSD to an A-form reference from the M-REMD (blue) and H-REMD (red) simulations.

Figure 7

RMSD population histograms at 300 K from the various REMD simulations. Shown are the normalized populations (y-axis) of structures at particular mass-weighted RMSD values (x-axis) of all atoms in residues 1–4 to a reference structure (A-form RNA) for the various REMD simulations.