Ab initio RNA folding by discrete molecular dynamics: from structure prediction to folding mechanisms

Abstract

RNA molecules with novel functions have revived interest in the accurate prediction of RNA three-dimensional (3D) structure and folding dynamics. However, existing methods are inefficient in automated 3D structure prediction. Here, we report a robust computational approach for rapid folding of RNA molecules. We develop a simplified RNA model for discrete molecular dynamics (DMD) simulations, incorporating base-pairing and base-stacking interactions. We demonstrate correct folding of 150 structurally diverse RNA sequences. The majority of DMD-predicted 3D structures have <4 A deviations from experimental structures. The secondary structures corresponding to the predicted 3D structures consist of 94% native base-pair interactions. Folding thermodynamics and kinetics of tRNA(Phe), pseudoknots, and mRNA fragments in DMD simulations are in agreement with previous experimental findings. Folding of RNA molecules features transient, non-native conformations, suggesting non-hierarchical RNA folding. Our method allows rapid conformational sampling of RNA folding, with computational time increasing linearly with RNA length. We envision this approach as a promising tool for RNA structural and functional analyses.

FIGURE 1.

Coarse-grained structural model of RNA employed in DMD simulations. (A) Three consecutive nucleotides, indexed i − 1, i, i + 1, are shown. Beads in the RNA: sugar (S), phosphate (P), and base (B). (Thick lines) Covalent interactions, (dashed lines) angular constraints, (dashed–dotted lines) dihedral constraints. Additional steric constraints are used to model base stacking. (B) Hydrogen bonding in RNA base pairing. (Dashed lines) The base-pairing contacts between bases B _{i − 1}:B_{j + 1} and B _i:B _j. A reaction algorithm is used (see Materials and Methods) for modeling the hydrogen bonding interaction between specific nucleotide base pairs.

FIGURE 2.

Ab initio RNA folding using DMD. (A) Fraction of native base pairs (Q-values) present in the predicted RNA 3D structure. The maximum Q-values during the course of simulations are also shown, which depict the conformational sampling efficiency of the DMD algorithm to reach the native states. We also show the Mfold predicted Q-values. (B) Scatter plots of RMSD for the final folded conformation with respect to the experimentally derived native structure as a function of RNA size. Large RNA molecules have increased fluctuations due to larger conformational freedom and consequently have greater RMSD from the native conformation. (C) Normalized histogram of predicted and least RMSD to the native RNA structure.

FIGURE 3.

Ab initio folding kinetics and energetics of a model pseudoknot RNA. (A) Superposition of experimental pseudoknot structure (NDB code: 1A60, ribbon) against DMD prediction (ribbon backbone trace with backbone spheres). Backbone ribbons are colored blue (N terminus) to red (C terminus). (B) Graph of specific heat of the pseudoknot molecule as a function of simulation temperature. (C) Two-dimensional potential of mean force 2D-PMF for pseudoknot folding at T* = 0.245 (corresponds to the major peak in the specific heat). (I₁, I₂) The two intermediate states, (N) native state. (D) The 2D-PMF plot at T* = 0.21. (E) Internucleotide base-pairing contact frequencies at the first folding intermediate (I₁) corresponding to the state where the 5′ hairpin is folded. (F) Internucleotide base-pairing contact frequencies at the second intermediate state (I₂) corresponding to the formation of the major groove helix stem of the 3′ pseudoknots. (G) Contact map of the native state (N) as observed in the experimental structure (NDB code: 1A60).

FIGURE 4.

Ab initio folding kinetics and energetics of a model tRNA. (A) Mg²⁺ binding site (sphere) in the tRNA. (B) Superposition of experimental tRNA structure (NDB code: 1EVV, ribbon) against DMD prediction (ribbon backbone trace with backbone spheres). Backbone ribbons are colored blue (N terminus) to red (C terminus). D loop, TΨC loop, anticodon loop, and acceptor loop are indicated with color representing their position in the tRNA secondary structure. (C) The specific heat of the tRNA molecule as the function of simulation temperature. (D) The 2D-PMF as the function of the total number of contacts and the number of native contacts at T* = 0.22. (I₁, I₂, I₃) Folding intermediates, (N) native conformation. (E,F) Folding events in the trajectories of tRNA replica exchange simulation. Two folding events in corresponding different replicas are observed out of eight replicas. (G–I) Internucleotide base-pairing contact frequencies at the threefolding intermediate states, I₁, I₂, and I₃, respectively. (J) Contact map of the native conformation (N) as observed in the experimental structure (NDB code: 1EVV).

FIGURE 5.

Thermodynamics of B-RNA and 72 RNA variants. (A) Specific heat: (circles) B-RNA, (diamonds) 72-RNA, (squares) 72-C RNA, (triangles) 72-14 RNA (shown in DMD units). (B) Superposition of experimental B-RNA structure (ribbon) against DMD prediction (ribbon with backbone spheres). Backbone ribbons are colored blue (N terminus) to red (C terminus). (C) 2D-PMF of B-RNA as the function of the number of total base pairs and native base pairs. We find that there are three major basins in the 2D-PMF corresponding to intermediate states I₁, I₂, and I₃. (I₄) Near-native intermediate conformation, (N) native conformation. (D) Internucleotide contact frequencies at the intermediate state with about zero native contacts (i.e., non-native state I₁). (E) Internucleotide contact frequencies at the B-RNA folding intermediate state with about five native contacts, (non-native state I₂). (F) Internucleotide contact frequencies at the B-RNA folding intermediate state I₃ with about nine native contacts. (G) Internucleotide contact frequencies at the B-RNA folding intermediate state I₄ at near-native conformation. (H) Contact map in the native state (N) observed in the experimental structure (NDB code: 1C2W).