Kinetic Mechanism of RNA Helix-Terminal Basepairing—A Kinetic Minima Network Analysis

Article

Kinetic Mechanism of RNA Helix-Terminal Basepairing—A Kinetic Minima Network Analysis Fengfei Wang,1 Li-Zhen Sun,2 Pinggen Cai,2 Shi-Jie Chen,3,* and Xiaojun Xu1,* 1 Institute of Bioinformatics and Medical Engineering, School of Mathematics and Physics, Jiangsu University of Technology, Changzhou, Jiangsu, China; 2Department of Applied Physics, Zhejiang University of Technology, Hangzhou, Zhejiang, China; and 3Department of Physics, Department of Biochemistry, and Informatics Institute, University of Missouri, Columbia, Missouri

ABSTRACT RNA functions are often kinetically controlled. The folding kinetics of RNAs involves global structural changes and local nucleotide movement, such as base flipping. The most elementary step in RNA folding is the closing and opening of a basepair. By integrating molecular dynamics simulation, master equation, and kinetic Monte Carlo simulation, we investigate the kinetics mechanism of RNA helix-terminal basepairing. The study reveals a six-state folding scheme with three dominant folding pathways of tens, hundreds, and thousands of nanoseconds of folding timescales, respectively. The overall kinetics is rate limited by the detrapping of a misfolded state with the overall folding time of 10
5 s. Moreover, the analysis examines the different roles of the various driving forces, such as the basepairing and stacking interactions and the ion binding/dissociation effects on structural changes. The results may provide useful insights for developing a basepair opening/closing rate model and further kinetics models of large RNAs.

SIGNIFICANCE This article addresses the kinetic mechanism for the most fundamental events in the process of RNA folding, namely RNA helix-terminal basepairing. We perform direct folding study by the kinetic minima network approach integrated with molecular dynamics simulation, kinetic Monte Carlo simulation, and the master equation method. Our study reveals a new kinetic mechanism for RNA helix-terminal basepairing: a six-state folding scheme with two dominant on pathways and one rate-limited off pathway. Different roles of various driving forces, such as basepairing/stacking and ion binding/dissociation, on structural changes are also examined in detail. These new kinetic rules would be highly useful for understanding RNA function such as RNA enzymatic reactions and for rational design of RNA-based molecular systems.

INTRODUCTION RNA is a polymeric molecule essential in coding, decoding, regulation, and expression (1–4). The functional structures can be anything from folded architectures to dynamically heterogeneous ensembles, and the folding timescales can range from microseconds for simple basepair opening and closing to seconds or even longer for complex structural refolding events (5–21). The elucidation of RNA folding landscapes, especially the knowledge of alternative stable structures and the structural rearrangements between functional conformations, is critical for uncovering the physical mechanisms of RNA function. Compared with the protein backbone, RNA backbone contains more rotatable bonds per residue, contributing to substantial conformational flexSubmitted May 31, 2019, and accepted for publication September 17, 2019. *Correspondence: [email protected] or [email protected] Editor: Alan Grossfield. https://doi.org/10.1016/j.bpj.2019.09.017 Ó 2019 Biophysical Society.

1674 Biophysical Journal 117, 1674–1683, November 5, 2019

ibility. The complexity of RNA folding energy landscape (22–30) and the limitations of experimental tools (31–41) have hampered our ability to understand RNA folding kinetics, such as RNA basepair formation and fraying, and larger structure folding kinetics. Many efforts have been made to predict RNA folding kinetics at the secondary structural level (42–52). The development of the kinetic rate model for basepair opening and closing (49–51), the most elementary kinetic steps in RNA folding, however, remains a challenge. To understand the kinetic mechanism of nucleic acids folding at single-nucleotide level, several theoretical studies have been carried out (53–62). For example, Colizzi and Bussi (54) applied external force-induced unfolding simulations and constructed the free energy landscapes for RNA duplexes, providing a molecular motion picture of helix opening. Wang et al. (56) obtained the thermodynamic and kinetic parameters of an RNA basepair through a long-time molecular dynamics (MD) simulation of the opening-closing

RNA Basepairing Kinetics

switch process of the basepair near its melting temperature. The predicted diffusion constant exhibited a super-Arrhenius behavior with strong and weak temperature-dependent opening and closing rates, respectively. Recently, Pinamonti et al. (57) used MD simulations and the Markov state model to characterize the kinetics of RNA fraying and its sequence and direction dependence. The results showed an intriguing interplay between the stability of intermediate states and the overall fraying kinetics. These previous studies, which focused mostly on the unfolding of basepairs, offered significant insights into the RNA basepair folding kinetics. Recently, by integrating MD simulation, kinetic Monte Carlo (KMC), and master equation (ME) methods, we quantitatively investigated the singlenucleotide kinetics in a basepairing process (58). Through direct folding simulation, our results revealed a four-state kinetic mechanism involving a trapped state and two dominant kinetic barriers for the overall single-nucleotide folding process. As suggested by the terminal mismatch energies of the Turner parameters from experiments (63) as well as the previous pulling simulations (54,55), the 30 -nucleotide is generally more stable than its 50 -nucleotide basepairing partner, and the basepair formation is expected to fold the 30 -nucleotide first, followed by the folding of the 50 -nucleotide. Therefore, our previous study focused on the second step of the basepair formation, namely the folding of the 50 -nucleotide after the 30 -nucleotide is folded to the native state. In this article, we employ a similar hybrid computational method to further investigate the kinetic mechanism of RNA basepair formation, which is the fundamental event in the process of helix dissociation, elongation, and structural rearrangements thereafter. In contrast to the previous simulations for temperature- or force-induced unfolding reactions, here, we directly simulate the basepair folding process, starting from the structures with both the 50 -nucleotide and the 30 -nucleotide unstacked and extract the detailed kinetics information from a novel kinetic minima network (KMN) model. The results would help elucidate the interplay between the two nucleotides, folding pathways, and rate-limiting steps as well as the roles of water and ions. Furthermore, the mechanism of RNA basepair folding kinetics may provide useful insights into the development of a rate model for basepairing/stacking.

build the initial folded structure of the 3-bp duplex with the sequence of 50 GGG30 -50 UCC30 , as shown in Fig. 1 A. Based on the 3-bp folded duplex structure, we rotate the torsional angles of the four bonds (P(G2)-O30 (G1) and C10 (G1)-N9(G1) for G1, and O30 (C5)-P(C6) and C10 (C6)-N1(C6) for C6, respectively; shown in Fig. S1) to generate the (complete) unfolded conformational ensemble for the G1-C6 basepair. Using fivefold uniform rotational angles for each bond, we generate an ensemble of 5 (4) ¼ 625 initial structures. Excluding structures disallowed by steric clashes, we obtain 600 viable initial unfolded structures (see stage 1 in Fig. 1 B).

MD simulations For each initial unfolded structure, after a 5000-step minimization for the whole system, we perform MD simulations at 310 K with the constraints of keeping the G2-C5 and G3-U4 basepairs fixed in the folded state (applying the position constraints on the C2, C4, and C6 atoms of the nucleotides G2, G3, U4, and C5 with the constraintScaling of 50.0 in the MD simulations), such that we focus only on the folding process of G1C6 basepair (Fig. 1 A). For each unfolded structure, we run a 100-ns MD simulation. The coordinates of all the atoms are written to the NAMD dcd file every 20 ps. In total, we simulate 600 MD trajectories. The stage

MATERIALS AND METHODS Initial structures The initial coordinates for the 3-bp helix are extracted from the experimentally determined hairpin structure (Protein Data Bank [PDB]: 1c0o (64)). The hairpin is embedded in a TIP3 water box with the water shell of ˚ . Sodium and chloride ions are used to neutralize the system and to 15 A keep 1 M sodium concentration of the system. Three different steps are used to minimize and relax the whole system: 1) fix all the atoms in RNA for 1000-step minimization; 2) fix the heavy atoms of RNA for 5 ns-long MD simulations at 310 K; and 3) fix the phosphorus atoms for 5 ns-long MD simulations at 310 K. We retain the first three basepairs to

FIGURE 1 (A) Folding of a G-C basepair from the unfolded to the native folded state. (B) Stage 1 includes exhaustive sampling for the discrete unfolded structures. The figure shows only 10 initial structures for clarity. Stage 2 includes enhanced conformational sampling with short-time MD simulations from the unfolded conformational ensemble; and stage 3 includes construction of the KMN based on the seven minima extracted from the free energy landscape (depicted with two reaction coordinates for clarity). To see this figure in color, go online.

Biophysical Journal 117, 1674–1683, November 5, 2019 1675

Wang et al. 2 in Fig. 1 B shows one of the sampled trajectories. All the simulations are performed using the NAMD 2.10 simulation package with the CHARMM27 force field (65). We keep the temperature at 310 K for the system through coupling to a Langevin heat bath. An integration time step of 2 fs is employed with the rigidBonds ‘‘on’’ option. Long-range electrostatic interactions are evaluated using the smooth particle mesh Ewald ˚ . The simulation is manually terminated method with a grid spacing of 1 A before 100 ns if the stable native folded structure is reached (root mean˚ ). square deviation (RMSD) < 1.5 A

KMN A major issue in MD simulations is that during the simulation, the sampled conformations often remain close to the initial structure, and the computational time is mostly spent on local sampling around the nearby free energy minima, as shown in stage 3 of Fig. 1 B. To circumvent this problem, we start from each and every (total 600) unfolded conformation in the complete conformational ensemble as the initial state for (short-time) MD simulations. The use of the complete ensemble of the unfolded states is expected to enhance the completeness of the conformation sampling. We use the following steps to construct a KMN network from the sampled short-time MD trajectories and predict the detailed long-time folding kinetics from the network: 1) To determine the KMN from the simulated trajectories. We identify the segments in the sampled trajectories (Fig. S2) such that the RMSDs between any two snapshots in the segment are less than a cutoff value ˚ ). We define the centroid of each such segment as a structure (1.5 A (node) of the KMN. After removing all the possible redundant structures, we find the KMN as defined by the respective centroid structures. 2) To extract kinetic connectivity and residence time. As illustrated in stage 2 of Figs. 1 B and S3, the RMSD values of the snapshots along a trajectory provide the information about the transitions between the minima and the residence time of each minimum. For the case shown in stage 2 of Fig. 1 B, the folding trajectory follows the sequential transitions of minimal M4 / M2 / M3 / M1. Therefore, we can obtain the connectivities M4 / M2, M2 / M3, and M3 / M1 and the corresponding residence times of the minima. 3) To calculate transition rate constants. From the kinetic connectivities between the minima extracted in step 2, we calculate the overall transition probability between a pair of minima: PMM’ ¼ NMM0 /NtotalM, where NMM0 and NtotalM are the numbers of occurrences for the M / M0 transition and for all the exit transitions from minimal M, respectively. To determine the rate constants for a transition M / M0 , we first compute the mean residence time for M from the MD trajectories (step 2). The inverse of gives the total rate ktotalM for all the transitions leaving M. The rate constant kMM0 for each individual transition M / M0 can be computed from the transition partition: kMM0 ¼ PMM0  ktotalM. 4) To construct KMN network and extract folding kinetics. Based on the kinetic minima and the interminimal rate constants (see stage 3 of Fig. 1 B), we apply KMC and ME methods to the KMN network and predict the folding kinetics of the system. In general, the ME method gives the overall relaxation kinetics (relaxation time and intermediate states) with a given initial condition, whereas the KMC simulation offers the detailed folding pathway information by tracking the first passage time of each state in the network (58).

RESULTS AND DISCUSSION Characterization of kinetic minima From the structural similarity as measured by the RMSD, we determine the kinetic minima from MD trajectories

1676 Biophysical Journal 117, 1674–1683, November 5, 2019

and rank them by the average residence time . Larger corresponds to a higher transition barrier between M and its kinetic neighbors, suggesting a higher possibility of a kinetic trap in folding kinetics. We use 0.05 ns as the threshold for to build the KMN network. Because the MD simulations are terminated once the stable native structure is reached, we manually add the folded state Mfolded as a kinetic minimal and treat it as a ‘‘black hole’’ with infinite residence time (i.e., no exit transitions from Mfolded). We also add the unfolded state, Munfolded, which is the conformational ensemble of all structures with both G1 and C6 bases flipping out. Moreover, we consider all the sampled structures, which are not in any energy basin of the extracted kinetic minima and have RMSD to the ˚ as the unfolded state native folded state Mfolded > 8.0 A Munfolded. In total, we have 36 minima, including Mfolded and Munfolded, in the KMN network. All structures (except Munfolded) are shown in Fig. S4. Figs. S5 and S6 show 10 independent 100 ns-long MD simulations at 310 K starting from the native folded (stacked) structure. Without applying enhanced conformational sampling technique, such as external force- and high temperature-induced enhanced sampling approaches, the regular MD simulations for quantitative study of the folding stability are usually computationally inefficient, especially at low temperatures. On the other hand, the stacking energy from the Turner parameters (63) for the base stacking between G1-C6 and G2-C5 is
3.3 kcal/mol, corresponding to a ratio of 200:1 between the stacked and unstacked populations at 310 K. Therefore, given the high stability of the stacked state, by assuming an absorbing boundary for the folded state, our kinetic model from direct folding simulations may provide a reliable description for the nonequilibrium process of helix-terminal basepairing. To better characterize the structural and energetic properties of the extracted kinetic minima, we define several order parameters. dxy is the distance between the geometric centers of all the heavy base atoms in the corresponding nucleotides. qxy is the angle between two base planes (defined by the atoms of C2, C4, and C6) in the corresponding nucleotides. Because the all-atom energies as determined by an MD force field are much more sensitive to the small changes in bond lengths, bond angles, and torsional angles than the RMSD values, we identify the conformations belonging to ˚ as a each minima by the RMSD calculations (1.5 A threshold) and use the averaged values of each energy parameter to illustrate the energetic properties of each centroid structure. Here, Exy is the nonbonded interaction energies (i.e., van der Waals and electrostatic energies analyzed by the VMD package with the CHARMM27 force field) between the corresponding nucleotides (see Supporting Materials and Methods for the detailed VMD scripts). The detailed information of all the order parameters is listed in Tables S1 and S2.

RNA Basepairing Kinetics

From the distributions of d12 and d56 in Fig. 2 A, we can see that all the extracted kinetic minima have at least one (of two, G1 or C6) nucleotide staying close to its adjacent nucleotide with a small value of either d12 or d56. This ‘‘L’’-shaped relationship between d12 and d56 suggests that 1) the energy landscape near the basepaired state Mfolded is more rugged than that for the unfolded states with both bases flipping out; 2) we may consider all the unfolded states as a macrostate, Munfolded, because of the relative smooth energy landscape with fast interconversions to reach a pre-equilibrium state; 3) and the exploration on the energy

landscape near Mfolded is mainly enthalpic to overcome the folding barriers on the rugged local energy landscape, whereas the entropic barrier dominates the folding process near Munfolded, which involves the conformational sampling of unfolded conformations on the smooth local folding landscape near Munfolded. From the energetic point of view as shown in Fig. 2 B, the single-stranded base stacking interactions E56 between the sequentially neighboring nucleotides C5 and C6 shows the dominant contributions compared with the base stacking interactions E12 between G1 and G2 and the basepairing interactions E16 between G1 and C6. Earlier stretching experiments suggested that RNA poly(C) is more stable than other single-stranded RNA homopolymers (66,67), which is consistent with the energy terms in Fig. 2 B. The most stable E16 of Mfolded indicates that the basepairing interactions also play a dominant role to stabilize the native folded state. On the other hand, the 30 dangling ends in RNA are generally more stable than the 50 dangling ends. From the Turner parameters in the NNDB (63), the 30 dangling C6 adjacent to C5-G2 has a stability of
0.8 kcal/mol, whereas the free energy for the 50 dangling G1 adjacent to G2-C5 is 0 kcal/mol. It should be noted that the overall energy landscape is determined by all the interactions, including the basepairing, and stacking contributions between nucleotides and the interactions with ions and water molecules. Different structures involve the competition between the different interactions. Furthermore, many extracted kinetic minima have parallel or antiparallel base orientations to Mfolded, which are revealed by the base-base angle qxy in Tables S1 and S2. Nucleotide C6 prefers to have the native-like base configuration with 17 out of 34 extracted kinetic minima in parallel orientation (q56 < 30 ). Nucleotide G1, however, favors the antiparallel orientations with 20 of 34 minima having q12 > 150 . As a result, to fold to the native folded state, these minima first flip out their bases so they can make necessary adjustments for their backbone and base configurations. After the rearrangement of the orientation, the bases flip back to the native configuration. Depending on the interactions between G1 and C6, the adjustment may be (partially) cooperative or independent. Validation of network

FIGURE 2 Selected order parameters to characterize the extracted 34 kinetic minima and the folded state. (A) d12 and d56 are the distances between the geometric centers of heavy base atoms between the corresponding nucleotides. (B) E12, E56, and E16 are the nonbonded interactions (van der Waals and electrostatic energies given by VMD) between the corresponding nucleotides. E12 and E56 characterize the single-stranded base stacking interactions. E16 represents the basepairing interactions. The stacking interaction E12, with the given force field, is much weaker than the other two interactions. Other order parameters are shown in Tables S1 and S2. To see this figure in color, go online.

With the connectivity between kinetic minima and the interminimal transition rate constants extracted from MD trajectories, we build a 36-minimal KMN network. To validate the KMN, we first bin all the conformations (snapshots) in the sampled trajectories into their respective minima, according to the structural similarity from RMSD calculations, to compute the time-dependent populations for each minimum. We then use the population of each minimum at time t ¼ 0 as the initial population and use the ME to predict the subsequent populational kinetics based on KMN. Note

Biophysical Journal 117, 1674–1683, November 5, 2019 1677

Wang et al.

that the conformational fluctuations (RMSD distribution) induced by interactions with water and ion molecules for each minimal shape the local energy landscape. For the small system studied here, with only one basepair folding freely while keeping others in their folded states, the ˚ . Other average RMSD fluctuation for Mfolded is 1.0 A minima states are usually more flexible with larger conformational fluctuations than Mfolded. Fig. S3 gives the details of KMN network extraction from the MD trajectories with a given RMSD threshold. Different RMSD thresholds can lead to different distributions of the residence time and interminimal connectivity, resulting in different KMN-predicted kinetics. We use different RMSD thresholds to extract MD-based and KMN-predicted populational kinetics. As we can see from Fig. S7, a smaller RMSD threshold gives better agreement. Fig. 3 A shows the results ˚ . We find that such KMNwith the RMSD threshold of 1.5 A predicted populational kinetics agrees with the original populations predicted from MD trajectories. The result suggests that the KMN may be reliable. Folding kinetics from network To obtain the overall folding kinetics starting from the unfolded structure, we set the initial fractional population of Munfolded to be 100%. The ME solution, as shown in Fig. S8, predicts that through conversions mainly involving M1, M2, M6, M12, M27, and M28, Munfolded quickly relaxes, and its population decays away within 30 ns, resulting in a gradual increase in the population of Mfolded (PMfolded / 11%). After such a rapid initial relaxation, M2, M6, and M12 show a significant populational jump, with peak populations PM2 z 12%, PM6 z 10%, and PM12 z 35% at t z 70, 60, and 700 ns, respectively. In the meantime, Mfolded shows a second populational jump and reaches 100% at t z 104 ns through transitions from other minima. The overall folding time for the kinetic process of Munfolded / Mfolded is 10
5 s, which is consistent with the experimental results about time of helix elongation of DNA and RNA (18–21). To better capture the main features of the basepairing kinetics, especially the folding pathways, we manually classify the 36 minima into six states, based on the behavior of the time-dependent population (kinetic profile shown in Fig. S8) for each minimum. Specifically, we identify the minima that emerge at the same/similar time with the same/similar time dependence of the population and group them as an effective macrostate. Because these minima are roughly formed commensurately, they may belong to the same macrostate such that their populational kinetics follows the kinetics of the macrostate. S1 for the unfolded structure ensembles Munfolded; S6 for the folded structure ensembles Mfolded; S2, S4, and S5 contain the minima with their populations showing one dominant peak at t z 10, 100, and 1000 ns, respectively; and S3 contains the minima

1678 Biophysical Journal 117, 1674–1683, November 5, 2019

FIGURE 3 (A) Comparison between the kinetic results from MD simulations (black) and from the 36-minimal KMN (red). (B) The predicted single basepair folding kinetics starting from the unfolded state S1, based on the six-state network (structures are shown in (C)). (D) The six-state network with two dominant interstate transitions (S1 4 S2, and S3 4 S4) and one dominant one-way transition (S5 / S3). The detailed connectivities and transition partitions are shown in Fig. S9. (E) Three folding pathways categorized by the KMC studies with their folding partitions of (0.14, 0.11, and 0.75). To see this figure in color, go online.

with their populations showing two significant peaks in the time window of (10, 1000) ns. From the original 36-minimal kinetic network, we use the KMC method to extract the connectivity and the residence time for each state (detailed data shown in Fig. S9) and reconstruct a reduced six-state kinetic network and predict the overall folding process starting from the unfolded state (S1) (see Fig. 3 B). We find that the states of S2, S3, S4, and S5 connect the unfolded and folded states of S1 and S6, and the overall folding kinetics

RNA Basepairing Kinetics

is rate limited by the detrapping of misfolded state S5. From the structural point of view, as shown in Fig. 3 C and Tables S1 and S2, states S2, S3, and S4 have different origins (parallel, antiparallel, or in-between) for both G1 and C6 nucleotides; state S5 is a misfolded state with both G1 and C6 nucleotides in the native folded position with antiparallel orientations (q12 and q56 > 150 ). The detrapping of misfolded S5 involves base flip for both nucleotides so that the rearrangement of the base and backbone conformations can occur, resulting in the largest average residence time and the rate-limiting step for the overall folding kinetics. Folding pathways From the detailed partitioning of the kinetic transitions, as shown in Figs. 3 D and S9, we find that the transitions from state S1 are well distributed, ranging from 0.07 to 0.35, among transitions to all the states. The transitions from state S2 are also well distributed, ranging from 0.15 to 0.42 to other states, except for the transition to S4. The average residence times for states S1 and S2 are 10.6 and 9.9 ns, respectively. As a result, the relatively significant interstate transition S1 4 S2 with the partition of 24% for S1 / S2 and 29% for S2 / S1 contributes to the initial rise of the S6 population in the time window of (0, 30) ns. As shown in Fig. 3 E, the KMC study confirms that 14% of the population flows from S1 to S6 with the help of the macrostate (S1, S2) (without reaching other states), which is preequilibrated through the S1 4 S2 interstate transitions. The initial relaxation process also leads to the population flows to the states S3, S4, and S5. As shown in Fig. S9, state S3 has two dominant transition partitions: 0.76 for S3 / S4 and 0.13 for S3 / S5. Both the S4 and S5 states show only one dominant transition: 0.95 for S4 / S3 and 0.9 for S5 / S3. The average residence times for states S3, S4, and S5 are 24.1, 30.8, and 350.1 ns, respectively. Because of the dominant interstate transition S3 4 S4 with the partition of 76% for S3 / S4 and 95% for S4 / S3, the KMC simulation results in Fig. 3 E show that 11% of the population flows through the pre-equilibrium macrostate (S1, S2, S3, S4) into S6 in the time window of (30, 500) ns. Because of the large residence time of S5 (350.1 ns on average), S5 reaches its populational peak of 0.35 at t z 600 ns after the two preceding time windows. Through transitions to other states, with the dominant transition partition of 0.9 for S5 / S3, the rest population merges eventually with S6 during the time window of (500, 10,000) ns. From the KMC simulations (Fig. 3 E), we find that 75% of the population from S1 must go through state S5 (within the macrostate (S1, S2, S3, S4, S5)) before reaching the folded state S6. Because of the two dominant interstate transitions (S1 4 S2 and S3 4 S4) and one dominant one-way transition (S5 / S3), as shown in Fig. 3 D, the KMC trajectories of the six-state network can be categorized into two on pathways and one off pathway (shown in Fig. 3 E): (I)

S1 / (S1, S2) / S6; (II) S1 / (S1, S2, S3, S4) / S6; and (III) S1 / (S1, S2, S3, S4, S5) / S6. Even though S4 prefers the antiparallel orientation of G1 and the parallel position for C6 with no directed transitions to the folded state S6, most of the folding processes would visit S4 before reaching S6. The fast exit kinetics (a shorter residence time of 30.8 ns compared to S5) through the dominant interstate transitions between S3 and S4 indicates that the folding pathway of II is mainly an on pathway. However, the largest residence time with one-way dominant transition partition of state S5 serves as a trapped state in an off pathway of III for the overall kinetics. On the other hand, we calculate the RMSD values for the base atoms of G1 and C6 with respect to the corresponding native configurations to better illustrate the helix-terminal basepair folding pathways. From the sampled 600 shorttime MD trajectories, we observed 123 events of folding into the native folded state. Among these trajectories, only 10 of them show the (almost) simultaneous folding of G1 and C6, as shown as an example in Fig. 4 A. The majority of the observed folding events from MD simulations have the alternative folding with one of the two nucleotides folding first, followed by the multiple adjustment of both nucleotides before reaching their native configurations, as shown as an example in Fig. 4 B (see Fig. S10 for more cases). Specifically, G1 (C6) folds first for 39 (74) of 123 events. Compared with the weak interaction E12 between G1 and G2, the much stronger interaction of E56 between C5 and C6, shown in Fig. 2 B, can drive the folding of C6 toward its native configuration efficiently, leading to the faster folding of nucleotide C6 than G1. Furthermore, the folding of the second nucleotide may require (partial) unfolding of the prefolded nucleotide to lower the kinetic

FIGURE 4 RMSD distributions for the base atoms of G1 (in blue) and C6 (in red) with respect to their native base conformations along two MD trajectories. (A) G1 and C6 fold simultaneously into the native configurations. (B) G1 and C6 fold alternately for the cooperative adjustment of base and backbone orientations. To see this figure in color, go online.

Biophysical Journal 117, 1674–1683, November 5, 2019 1679

Wang et al.

barrier for the folding of both nucleotides to proceed. The basepairing interactions between G1 and C6, defined as E16, along with the stacking interactions E12 and E56 and other interactions, serve as the major driving forces for the adjustment of base and backbone orientations. Ion effects To investigate the ion effect in the basepairing kinetics, we track the ion positions and the ion-RNA interactions in the sampled MD trajectories. To account for the dynamic nature ˚ (to the of MD simulations, we use the RMSD cutoff of 1.5 A centroid structures) to select sampled structures in each kinetic minimal basin and consider the averaged ion distributions and ion-RNA interactions as the features of each kinetic minimal. Specifically, we calculate the minimal distance of each ion (Naþ and Cl
) to non-H atoms in RNA (i.e., the distance to RNA surface) and bin the extracted dis˚ to tances into tiers centered at 1.5, 4.5, 7.5, 10.5, and 13.5 A RNA surface (see the inset of Fig. 5 A) to show the ion distribution for each kinetic minima. Fig. 5, A and B give the ions (Naþ, Cl
) and the net charge (NNaþ - NCl
) distributions of the Mfolded state (S6) and the 34 kinetic minima. Because RNA molecules are negatively charged with
1e for each nucleotide, the negatively charged RNA is (gradually) neutralized by the ˚ to net charge distributions of Naþ and Cl
within 10 A RNA surface. Specifically, the
6e (with the CHARMM force field) of this 3-bp RNA system is neutralized by the net charges of þ3.0e, þ2.5e, and þ0.5e in the spatial layer

˚ , respectively. The (small) strucof (0, 3), (3, 6), and (6, 9) A tural differences between different kinetic minima (with one basepair folding/unfolding, whereas the other two basepairs constraint) leads to fluctuations in ion distributions for both Naþ and Cl
, especially in the vicinity of the RNA surface. On the other hand, as we can see from Figs. 5 B and S11, the S4 and S5 states have different net charge distributions from the folded native state (S6, in red). However, some of the minima in S2 and S3 exhibit similar distributions with S6. From the interstate kinetic connectivities shown in Fig. 3 D, there is no apparent event for the S4 / S6 and S5 / S6 transitions. Therefore, the net charge flow, mainly through Naþ, from the nearest to the next-nearest region, may help to trigger the detrapping from the non-native structures with the misfolded base flipping out to reach the native folded state. The interplay between the atoms of G1-C6 basepair and the surrounding ions may play an important role in the overall folding kinetics. Fig. 5, C and D give the distributions of the distances between the Naþ ions and the nucleotides G1 and C6 and the bond lengths of three hydrogen bonds (O6N4, N1-N3, and N2-O2) involved in the G1-C6 basepair along an MD trajectory, respectively. The Naþ binding in the vicinity of the RNA surface during the time windows of (14, 21) ns and (24, 30) ns precedes the formation of the stable G1-C6 basepair. Therefore, the change of ion distribution along the trajectories (i.e., ion binding and dissociation), is directly coupled to the RNA folding events. Physically, ion binding can lower the electrostatic repulsion between nucleotides, thus inducing the close approach and

FIGURE 5 (A) Ion distributions of Naþ and Cl
, obtained by tracking their occurrences in the spatial ˚ in width) centered at the distances of layers (of 3 A ˚ to RNA surface, as 1.5, 4.5, 7.5, 10.5, and 13.5 A illustrated in inset, for all the kinetic minima (except Munfolded) from MD simulations. (B) The average net charge (NNaþ - NCl
) distributions in states S2, S3, S4, S5, and S6 is shown. (C) Distribution of the distances between Naþ and the atoms of the nucleotides G1 and C6 along an MD trajectory is shown. (D) Distribution of the bond lengths of three hydrogen bonds (O6-N4, N1-N3, and N2O2) involved in the G1-C6 basepair along an MD trajectory is shown. To see this figure in color, go online.

1680 Biophysical Journal 117, 1674–1683, November 5, 2019

RNA Basepairing Kinetics

interaction between RNA bases/backbone. Depending on the nucleotide type and the location of the bound ions, different nucleotides may experience different ion effects. Table S3 lists three ion-involved energies for each kinetic minimal, calculated from VMD. EG1
ion, EC6
ion, and ERNA
ion are the interaction energies between G1, C6, and the whole RNA and the ions (Naþ and Cl
) in solution, respectively. EG1
ion is in the region of (
284.5,
101.9) kcal/mol, whereas EC6
ion is within (
88.6,
22.8) kcal/mol. The relative weak interactions from ions to nucleotide C6 may also contribute to the faster folding of nucleotide C6 than G1 as observed from MD simulations. Furthermore, the averaged ERNA
ion for states of S2, S3, S4, and S5 are
468.1,
510.8,
569.9, and
542.9 kcal/mol (see details in Table S3). The release of Naþ from the nearest to the next-nearest region destabilizes the misfolded states (especially S4 and S5), thus lowering the free energy barrier for the detrapping transition from the misfolded state to reach the folded state S6 with ERNA
ion of
440.3 kcal/mol. CONCLUSION Using an integrated approach with all-atom MD simulations, the ME method, and KMC simulations, we perform direct folding studies for RNA helix-terminal basepairing. Based on the kinetic minima extracted from short-time MD trajectories, we reveal a six-state kinetic mechanism with three dominant folding pathways. Because of the two dominant interstate transitions of S1 4 S2 and S3 4 S4 and one dominant one-way transition of S5 / S3 from the unfolded state S1, 14% of the population flows into the folded state S6 during the initial relaxation process of states S1 and S2 in the time window of (0, 30) ns, followed by 11% of the population flowing into S6 through steps involving states S3 and S4 in the time window of (30, 500) ns. The overall folding kinetics is rate limited by the detrapping of the misfolded state S5 in the time window of (500, 10,000) ns, which promotes the folding of the rest 75% of the population with the overall folding time of 10
5 s. RNA folding is a result of the complex interplay of the different forces, such as basepairing and stacking, RNAion, and other interactions. For the folding of helix-terminal basepair, we find that compared with the base stacking interactions E12 between G1 and G2 and the basepairing interactions E16 between G1 and C6, the single-stranded base stacking interactions E56 between the sequentially neighboring nucleotides C5 and C6 shows the dominant contributions, resulting in a faster folding of nucleotide C6 than G1. Furthermore, the interactions EG1
ion between nucleotide G1 and ions in the solution is much stronger than EC6
ion, which may also lead to a faster folding of nucleotides C6 than G1. The release of Naþ from the nearest to the nextnearest region of RNA surface destabilizes the misfolded states, thus lowering the free energy barrier for the detrap-

ping transition from the misfolded states to reach the folded state. In our approach, the reduction from hundred-thousands of snapshots (MD trajectories) to 36 kinetic minima and to six states allows us to focus on the overall folding kinetics while leaving the microscopic details of the energy landscape to the transition partitions and residence time of each minimum. The method here can be readily extended to study the sequence-dependent basepair folding process and other more complicated biological phenomena, such as singlenucleotide flipping in the middle of RNA/DNA helices and ion-induced bulge/internal loop bending. SUPPORTING MATERIAL Supporting Material can be found online at https://doi.org/10.1016/j.bpj. 2019.09.017.

AUTHOR CONTRIBUTIONS X.X. and S.-J.C. designed the research. F.W. and P.C. performed the simulations. L.-Z.S. and X.X. performed the analysis. F.W., X.X., and S.-J.C. contributed to writing the article.

ACKNOWLEDGMENTS This work was supported by Zhejiang Natural Science Foundation for Youth LQ17C100002 to P.C., the National Science Foundation of China 11704333 to L.-Z.S., and the National Institutes of Health R01GM063732 and R01- GM117059 to S.-J.C.

REFERENCES 1. Esteller, M. 2011. Non-coding RNAs in human disease. Nat. Rev. Genet. 12:861–874. 2. Quinn, J. J., and H. Y. Chang. 2016. Unique features of long noncoding RNA biogenesis and function. Nat. Rev. Genet. 17:47–62. 3. Anastasiadou, E., L. S. Jacob, and F. J. Slack. 2018. Non-coding RNA networks in cancer. Nat. Rev. Cancer. 18:5–18. 4. Ransohoff, J. D., Y. Wei, and P. A. Khavari. 2018. The functions and unique features of long intergenic non-coding RNA. Nat. Rev. Mol. Cell Biol. 19:143–157. 5. Russell, R., I. S. Millett, ., L. Pollack. 2002. Rapid compaction during RNA folding. Proc. Natl. Acad. Sci. USA. 99:4266–4271. 6. Wenter, P., B. F€urtig, ., S. Pitsch. 2005. Kinetics of photoinduced RNA refolding by real-time NMR spectroscopy. Angew. Chem. Int.Engl. 44:2600–2603. 7. Wenter, P., B. F€urtig, ., S. Pitsch. 2006. A caged uridine for the selective preparation of an RNA fold and determination of its refolding kinetics by real-time NMR. Chembiochem. 7:417–420. 8. Al-Hashimi, H. M., and N. G. Walter. 2008. RNA dynamics: it is about time. Curr. Opin. Struct. Biol. 18:321–329. 9. F€urtig, B., P. Wenter, ., H. Schwalbe. 2010. Probing mechanism and transition state of RNA refolding. ACS Chem. Biol. 5:753–765. 10. Rinnenthal, J., J. Buck, ., H. Schwalbe. 2011. Mapping the landscape of RNA dynamics with NMR spectroscopy. Acc. Chem. Res. 44:1292– 1301.

Biophysical Journal 117, 1674–1683, November 5, 2019 1681

Wang et al. 11. Dethoff, E. A., J. Chugh, ., H. M. Al-Hashimi. 2012. Functional complexity and regulation through RNA dynamics. Nature. 482:322– 330.

34. Ha, T., X. Zhuang, ., S. Chu. 1999. Ligand-induced conformational changes observed in single RNA molecules. Proc. Natl. Acad. Sci. USA. 96:9077–9082.

12. Dupuis, N. F., E. D. Holmstrom, and D. J. Nesbitt. 2014. Molecularcrowding effects on single-molecule RNA folding/unfolding thermodynamics and kinetics. Proc. Natl. Acad. Sci. USA. 111:8464–8469.

35. Nagel, J. H., A. P. Gultyaev, ., C. W. Pleij. 2002. A pH-jump approach for investigating secondary structure refolding kinetics in RNA. Nucleic Acids Res. 30:e63.

13. Mustoe, A. M., C. L. Brooks, and H. M. Al-Hashimi. 2014. Hierarchy of RNA functional dynamics. Annu. Rev. Biochem. 83:441–466.

36. Jung, J., and A. Orden Van. 2006. A three-state mechanism for DNA hairpin folding characterized by multiparameter fluorescence fluctuation spectroscopy. J. Am. Chem. Soc. 128:1240–1249.

14. Zhang, X., X. Xu, ., L. Q. Gu. 2015. Mimicking ribosomal unfolding of RNA pseudoknot in a protein channel. J. Am. Chem. Soc. 137:15742–15752. 15. Herschlag, D., B. E. Allred, and S. Gowrishankar. 2015. From static to dynamic: the need for structural ensembles and a predictive model of RNA folding and function. Curr. Opin. Struct. Biol. 30:125–133. 16. Watters, K. E., E. J. Strobel, ., J. B. Lucks. 2016. Cotranscriptional folding of a riboswitch at nucleotide resolution. Nat. Struct. Mol. Biol. 23:1124–1131. 17. Espah Borujeni, A., and H. M. Salis. 2016. Translation initiation is controlled by RNA folding kinetics via a ribosome drafting mechanism. J. Am. Chem. Soc. 138:7016–7023.

37. F€urtig, B., J. Buck, ., H. Schwalbe. 2007. Time-resolved NMR studies of RNA folding. Biopolymers. 86:360–383. 38. Fogarty, K., J. T. McPhee, ., A. Van Orden. 2009. Probing the ionic atmosphere of single-stranded DNA using continuous flow capillary electrophoresis and fluorescence correlation spectroscopy. Anal. Chem. 81:465–472. 39. Gupta, A. N., A. Vincent, ., M. T. Woodside. 2011. Experimental validation of free-energylandscape reconstruction from non-equilibrium single-molecule force spectroscopy measurements. Nat. Phys. 7:631–634.

18. Po¨rschke, D. 1974. Thermodynamic and kinetic parameters of an oligonucleotide hairpin helix. Biophys. Chem. 1:381–386.

40. Ding, B., M. R. Hilaire, and F. Gai. 2016. Infrared and fluorescence assessment of protein dynamics: from folding to function. J. Phys. Chem. B. 120:5103–5113.

19. Po¨rchke, D. 1976. The nature of stacking interations in polynucleotides. Molecular states in Oligo- and polyribocytidylic acids by relaxation analysis. Biochemistry. 15:1495–1499.

41. Orden, A. V., and J. Jung. 2008. Review fluorescence correlation spectroscopy for probing the kinetics and mechanisms of DNA hairpin formation. Biopolymers. 89:1–16.

20. Dewey, T. G., and D. H. Turner. 1980. Laser temperature jump study of solvent effects of poly(adenylic acid) stacking. Biochemistry. 19:1681– 1685.

42. Chen, S. J., and K. A. Dill. 2000. RNA folding energy landscapes. Proc. Natl. Acad. Sci. USA. 97:646–651.

21. Gueron, M., and J. L. Leroy. 1995. Studies of base pair kinetics by NMR measurement of proton exchange. Methods Enzymol. 261:383–413. 22. Collin, D., F. Ritort, ., C. Bustamante. 2005. Verification of the Crooks fluctuation theorem and recovery of RNA folding free energies. Nature. 437:231–234.

43. Cao, S., and S. J. Chen. 2009. A new computational approach for mechanical folding kinetics of RNA hairpins. Biophys. J. 96:4024–4034. 44. Liu, F., H. Tong, and Z. C. Ou-Yang. 2006. Force unfolding single RNAs. Biophys. J. 90:1895–1902. 45. Xayaphoummine, A., V. Viasnoff, ., H. Isambert. 2007. Encoding folding paths of RNA switches. Nucleic Acids Res. 35:614–622.

23. Chen, S. J. 2008. RNA folding: conformational statistics, folding kinetics, and ion electrostatics. Annu. Rev. Biophys. 37:197–214.

46. Geis, M., C. Flamm, ., C. Thurner. 2008. Folding kinetics of large RNAs. J. Mol. Biol. 379:160–173.

24. Woodson, S. A. 2010. Compact intermediates in RNA folding. Annu. Rev. Biophys. 39:61–77.

47. Hofacker, I. L., C. Flamm, ., P. F. Stadler. 2010. BarMap: RNA folding on dynamic energy landscapes. RNA. 16:1308–1316.

25. Solomatin, S. V., M. Greenfeld, ., D. Herschlag. 2010. Multiple native states reveal persistent ruggedness of an RNA folding landscape. Nature. 463:681–684.

48. Dotu, I., W. A. Lorenz, ., P. Clote. 2010. Computing folding pathways between RNA secondary structures. Nucleic Acids Res. 38:1711–1722.

26. Yang, C., E. Kim, and Y. Pak. 2015. Free energy landscape and transition pathways from Watson-Crick to Hoogsteen base pairing in free duplex DNA. Nucleic Acids Res. 43:7769–7778. 27. Miner, J. C., A. A. Chen, and A. E. Garcı´a. 2016. Free-energy landscape of a hyperstable RNA tetraloop. Proc. Natl. Acad. Sci. USA. 113:6665–6670. 28. Manz, C., A. Y. Kobitski, ., G. U. Nienhaus. 2017. Single-molecule FRET reveals the energy landscape of the full-length SAM-I riboswitch. Nat. Chem. Biol. 13:1172–1178. 29. Chakraborty, D., and D. J. Wales. 2018. Energy landscape and pathways for transitions between Watson-Crick and Hoogsteen base pairing in DNA. J. Phys. Chem. Lett. 9:229–241. 30. Bello-Rivas, J. M., and R. Elber. 2016. Simulations of thermodynamics and kinetics on rough energy landscapes with milestoning. J. Comput. Chem. 37:602–613. 31. Lynch, D. C., and P. R. Schimmel. 1974. Effects of abnormal base ionizations on Mg2 plus binding to transfer ribonucleic acid as studied by a fluorescent probe. Biochemistry. 13:1852–1861.

49. Zhang, W., and S. J. Chen. 2002. RNA hairpin-folding kinetics. Proc. Natl. Acad. Sci. USA. 99:1931–1936. 50. Xu, X., and S. J. Chen. 2012. Kinetic mechanism of conformational switch between bistable RNA hairpins. J. Am. Chem. Soc. 134:12499–12507. 51. Srinivas, N., T. E. Ouldridge, ., E. Winfree. 2013. On the biophysics and kinetics of toehold-mediated DNA strand displacement. Nucleic Acids Res. 41:10641–10658. 52. Xu, X., D. Duan, and S. J. Chen. 2017. CRISPR-Cas9 cleavage efficiency correlates strongly with target-sgRNA folding stability: from physical mechanism to off-target assessment. Sci. Rep. 7:143. 53. Hagan, M. F., A. R. Dinner, ., A. K. Chakraborty. 2003. Atomistic understanding of kinetic pathways for single base-pair binding and unbinding in DNA. Proc. Natl. Acad. Sci. USA. 100:13922–13927. 54. Colizzi, F., and G. Bussi. 2012. RNA unwinding from reweighted pulling simulations. J. Am. Chem. Soc. 134:5173–5179. 55. Zgarbova´, M., M. Otyepka, ., P. Jurecka. 2014. Base pair fraying in molecular dynamics simulations of DNA and RNA. J. Chem. Theory Comput. 10:3177–3189.

32. Bina-Stein, M., and D. M. Crothers. 1975. Localization of the structural change induced in tRNA fMET (Escherichia coli) by acidic pH. Biochemistry. 14:4185–4191.

56. Wang, Y., S. Gong, ., W. Zhang. 2016. The thermodynamics and kinetics of a nucleotide base pair. J. Chem. Phys. 144:115101.

33. Sclavi, B., M. Sullivan, ., S. A. Woodson. 1998. RNA folding at millisecond intervals by synchrotron hydroxyl radical footprinting. Science. 279:1940–1943.

57. Pinamonti, G., F. Paul, ., G. Bussi. 2019. The mechanism of RNA base fraying: molecular dynamics simulations analyzed with core-set Markov state models. J. Chem. Phys. 150:154123.

1682 Biophysical Journal 117, 1674–1683, November 5, 2019

RNA Basepairing Kinetics 58. Xu, X., T. Yu, and S. J. Chen. 2016. Understanding the kinetic mechanism of RNA single base pair formation. Proc. Natl. Acad. Sci. USA. 113:116–121. 59. Pinamonti, G., J. Zhao, ., G. Bussi. 2017. Predicting the kinetics of RNA oligonucleotides using Markov state models. J. Chem. Theory Comput. 13:926–934. 60. Wang, Y., Z. Wang, ., W. Zhang. 2018. The nearest neighbor and next nearest neighbor effects on the thermodynamic and kinetic properties of RNA base pair. J. Chem. Phys. 148:045101. 61. Bottaro, S., G. Bussi, ., K. Lindorff-Larsen. 2018. Conformational ensembles of RNA oligonucleotides from integrating NMR and molecular simulations. Sci. Adv. 4:eaar8521.  62. Sponer, J., G. Bussi, ., M. Otyepka. 2018. RNA structural dynamics as captured by molecular simulations: a comprehensive overview. Chem. Rev. 118:4177–4338.

63. Turner, D. H., and D. H. Mathews. 2010. NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure. Nucleic Acids Res. 38:D280–D282. 64. Colmenarejo, G., and I. Tinoco, Jr. 1999. Structure and thermodynamics of metal binding in the P5 helix of a group I intron ribozyme. J. Mol. Biol. 290:119–135. 65. Kale, L., R. Skeel, ., K. Schulten. 1999. NAMD2: greater scalability for parallel molecular dynamics. J. Comput. Phys. 151:283–312. 66. Seol, Y., G. M. Skinner, and K. Visscher. 2004. Elastic properties of a single-stranded charged homopolymeric ribonucleotide. Phys. Rev. Lett. 93:118102. 67. Seol, Y., G. M. Skinner, ., A. Halperin. 2007. Stretching of homopolymeric RNA reveals single-stranded helices and base-stacking. Phys. Rev. Lett. 98:158103.

Biophysical Journal 117, 1674–1683, November 5, 2019 1683