key: cord-0888592-x28cxlm3
authors: Hori, Naoto; Denesyuk, Natalia A.; Thirumalai, D.
title: Salt Effects on the Thermodynamics of a Frameshifting RNA Pseudoknot under Tension
date: 2016-04-28
journal: Journal of molecular biology
DOI: 10.1016/j.jmb.2016.06.002
sha: ef4472aa4e66a2bfa77b1283679511564e523859
doc_id: 888592
cord_uid: x28cxlm3

Because of the potential link between -1 programmed ribosomal frameshifting and response of a pseudoknot (PK) RNA to force, a number of single molecule pulling experiments have been performed on PKs to decipher the mechanism of programmed ribosomal frameshifting. Motivated in part by these experiments, we performed simulations using a coarse-grained model of RNA to describe the response of a PK over a range of mechanical forces ($f$s) and monovalent salt concentrations ($C$s). The coarse-grained simulations quantitatively reproduce the multistep thermal melting observed in experiments, thus validating our model. The free energy changes obtained in simulations are in excellent agreement with experiments. By varying $f$ and $C$, we calculated the phase diagram that shows a sequence of structural transitions, populating distinct intermediate states. As $f$ and $C$ are changed, the stem-loop tertiary interactions rupture first, followed by unfolding of the $3^{prime}$-end hairpin ($textrm{I}rightleftharpoonstextrm{F}$). Finally, the $5^{prime}$-end hairpin unravels, producing an extended state ($textrm{E}rightleftharpoonstextrm{I}$). A theoretical analysis of the phase boundaries shows that the critical force for rupture scales as $left(log C_{textrm{m}}right)^{alpha}$ with $alpha=1,(0.5)$ for $textrm{E}rightleftharpoonstextrm{I}$ ($textrm{I}rightleftharpoonstextrm{F}$) transition. This relation is used to obtain the preferential ion-RNA interaction coefficient, which can be quantitatively measured in single-molecule experiments, as done previously for DNA hairpins. A by-product of our work is the suggestion that the frameshift efficiency is likely determined by the stability of the $5^{prime}$-end hairpin that the ribosome first encounters during translation.

The multifarious roles RNA molecules play in controlling a myriad of cellular functions [1] have made it important to understand in quantitative terms their folding [2] [3] [4] [5] [6] and how they respond to external stresses [7] . Among them, one of the simplest structural motifs is the RNA pseudoknot (PK), which is involved in many biological functions. The simplest type, referred to as H-type PK, consists of two stems connected by two loops in which one of the stems forms base pairs with the loop of the other. The PK motif, found in many RNA molecules such as telomerase [8] , mRNA [9] , ribosomal RNA [10] , transfer-messenger RNA, and viral RNA [11] , is functionally important [12] . For instance, a PK in the human telomerase RNA is essential for enzyme (a ribonucleoprotein complex) activity [13] . Many pathogenic mutations in the human telomerase RNA are found in the wellconserved PK region of the RNA, further underscoring the importance of the PK [8, 13] . The presence of PKs is also important in ribozyme catalysis and inhibition of ribosome translation.

Recently, there is heightened interest in the biophysics of PK folding because it plays a crucial role in affecting the efficiency of −1 programmed ribosomal frameshifting (−1 PRF) [14] [15] [16] . Usually, proteins are synthesized when the ribosome reads the mRNA code in three nucleotide steps until the stop codon is reached. In −1 PRF, the open reading frame of mRNA being translated within the ribosome is programmed to be shifted by one nucleotide, and consequently, the mRNA becomes nonfunctional or produces an entirely different protein [5, [17] [18] [19] [20] . The downstream PK of the mRNA is a roadblock at the entrance to the mRNA channel on the ribosome, which impedes the translocation of the ribosome. The ribosome must unfold the PK, presumed to occur by exertion of mechanical force, to complete translocation. Because frameshifting efficiency could depend on how the PK responds to force, a number of single-molecule pulling experiments have focused on PKs [15, 21, 22] . Several factors could determine −1 PRF, as evidenced by the multitude of proposals based on many experimental studies [14, 15, [23] [24] [25] [26] [27] [28] . Nevertheless, given that ribosome translocation exerts mechanical force on the downstream mRNA, it is physically reasonable to assume that resistance of PK to force should at least be one factor determining the frameshifting efficiency [14, 29] .

The considerations given above and the availability of both ensemble and single-molecule pulling experimental data prompted us to focus on the effects of temperature, mechanical force, and salt concentration on a 28nucleotide H-type RNA PK from the beet western yellow virus (BWYV) [19] . This PK has been the subject of several ensemble experiments, which have characterized the folding thermodynamics as a function of monovalent and divalent cations [30, 31] . The BWYV PK has two stems (S1 and S2) connected by two loops (L1 and L2) as shown in Fig. 1 . The crystal structure reveals that L2 forms triplex-type tertiary interactions with S1 [32] . In another context, it has been shown that such interactions and mutations affecting them dramatically change the efficiency of frameshifting [21] . [32] color-coded according to secondary structures. The black arrows indicate positions at which the external mechanical force is applied. (B) Schematic representation of the secondary structure. Unlike the typical H-type pseudoknot, the two helix stems (S1 and S2) are not coaxially stacked in the BWYV pseudoknot [25] . Tertiary interactions are indicated by red dashed lines.

variations in the frameshift efficiency, which, for some mutants, is attributed to changes in the stem-loop interactions affecting the stability of PK. For others, direct involvement of the ribosome is invoked [25] . In addition to mechanical force (f ), it is known from several ensemble experiments that the stability of BWYV PK, like other RNA molecules, depends on monovalent and divalent cations [33] . Thus, it is important to understand in quantitative detail how this crucial structural motif respond to f and changes in ion concentration.

Although single-molecule pulling experiments suggest how force affects the stability of RNA [19] , structural information obtained from such experiments is limited only to changes in molecular extensions. In order to provide a theoretical basis for the changes in the stability and structural transitions as f and monovalent salt concentration (C) are altered, we performed simulations using a coarse-grained three-site-interaction model (TIS model) [34] of the PK. After demonstrating that our simulations reproduce the temperature-induced structural transitions observed in ensemble experiments, we calculated the global phase diagrams in the [C, f ] plane. The simulations, performed in the same force range as used in Laser Optical Tweezer (LOT) experiments [19] , quantitatively reproduce the experimental observations. By combining our simulations and a theory based on ion-RNA interactions expressed in terms of preferential interaction coefficients [31, 35] , we produce several key predictions: (1) The phase diagram in the [C, f ] plane shows that the folded state (F) ruptures in stages populating dis-tinct intermediates as f increases. Structurally, the low force intermediate is similar to the one populated during thermal melting. The sequence of f -induced structural transitions can be described as F I E, where E is an extended state; (2) The critical force, f c , to rupture the PK scales as (log C m ) α with α = 1 for the I E transition and α = 0.5 for the I F transition. We expect this result to be valid for H-type PKs composed of stems with vastly different stability, such as the BWYV PK. This result is of general validity and is applicable to other PKs as well.

The slope of the f c versus log C m is proportional to the changes in the preferential interaction coefficient between an intermediate and unfolded states. This prediction can be used in conjunction with single-molecule experiments allowing for direct measurement of a key quantity in ion-RNA interactions; and (3) To the extent the downstream PK is a roadblock for translation, our work suggests that the ribosome should exert a force ∼ (10 − 15) pN [36] for unimpeded translocation along mRNA, although the actual force value could vary depending on factors such as interaction between PK and ribosome. Based in part on our work we propose that there should be a link between the stimulatory effect of frameshifting and the stability of the 5 -end of the PK. Our predictions are amenable to validations using single-molecule experiments.

In order to validate our model, we first performed temperature replica-exchange molecular dynamics (REMD) simulations at C = 500 mM with f = 0. Formations of stems (S1 and S2) and tertiary (stem-loop) interactions (L1 and L2) are characterized by hydrogen bond (HB) energies U S HB and U L HB , respectively. Changes in the distributions of HB energies as the temperature is increased are shown in Fig. S1 in the Supporting Information. Fig. S1 shows that there are two minima corresponding to two distinct states. The minimum that dominates at T = 20 • C corresponds to the F state in which both S1 and S2 are formed (U S HB ≈ −45 kcal/mol). The minimum with U S HB ≈ −30 kcal/mol at 80 • C corresponds to a state in which S1 is intact (Fig. S1 ). The free energy profile of U L HB shows only one possible minimum around U L HB ≈ −20 kcal/mol at T = 20 • C, indicating that L1 and L2 fold/unfold cooperatively as temperature changes.

We calculated the free energy profiles G(R α ) = −k B T log P (R α ) (R α = R ee , end-to-end distance, or R g , radius of gyration) from the distributions of P (R α ), which are given in Fig. S2 in the Supporting Information. The locations and the number of minima report on the thermodynamic transitions between the states (Fig. 2) .

(1) At a low temperature, T = 20 • C, the PK is completely folded (F state) including all the structural motifs, S1, S2, L1 and L2. The existence of F state is also confirmed by the presence of a minimum at R ee ≈ 3 nm ( Fig. 2A) and at R g ≈ 1.3 nm (Fig. 2B) . These values are similar to those obtained directly from the crystal structure (R native ee = 2.8 nm, R native g = 1.3 nm). (2) At T = 60 • C, the free energy profile for U L HB shows that some of the tertiary interactions, involving L1 and L2, are no longer stably formed (Fig. S1B ). On the other hand, G(R g ) shows two stable minima, indicating that there are two possible states (Fig. 2B) . In one state, both stems S1 and S2 are formed but tertiary interactions are absent, corresponding to the I1 state identified using UV absorbance and differential scanning calorimetry experiments (Fig. 3) [30, 31] . The other conformation, in which only S1 is formed, has been also observed in experiments and is referred to as the I2 state. It should be noted that the distribution of U S HB also has two minima at 60 • C (Fig. S1A) whereas G(R ee ) has only one minimum ( Fig. 2A). (3) At a still higher temperature T = 80 • C, G(R ee ) and G(R g ) each have a minimum with S1 stably formed. This is also reflected in Fig. S1 , which shows a minimum at U S HB ≈ −30 kcal/mol and a minimum at U S HB ≈ 0. Thus, completely unfolded conformations, U state, and the I2 state are populated. (4) At T = 110 • C, both U S HB and U L HB are 0 , indicating that only the U state exists. This is also reflected in G(R ee ) and G(R g ), which show that the PK is completely unfolded. In both profiles (Figs. 2A and 2B), the center of minimum is located at larger values than in the F state (R ee ≈ 5 nm and R g ≈ 2 nm).

The simulation results in Fig. 2 and Fig. S1 show that melting of BWYV PK undergoes through two intermediate states, F I1 I2 U where I1 has the topology of the F state but without significant tertiary interactions (Fig. 3) . Although simulations can be used to distinguish I1 and F states, the intermediate does not have significant population. The identification of I1 in our simulations provides a structural support to experiments in which I1 was inferred from three broad overlapping peaks in the heat capacity curves [30, 31] . It should be stressed that the experimental heat capacity curves exhibit only two discernible peaks, which would normally imply that there is only one intermediate present. The range of melting temperatures for the three transitions can be obtained from the temperature dependence of the various states f α (T ) (α=F, I1, I2 or U) given in Fig. 2C . The I1 state has a maximum population around T ≈ 50 • C. We equate this with the melting temperature T m1 associated with the F I1 transition. Similarly, T m2 for the I1 I2 transition is ≈ 70 • C. The T m3 for the I2 U transition occurs around 90

. The values of the three melting temperatures obtained in our simulations are in remarkable agreement with experiments (T exp. m1 =59.4 • C, T exp. m2 =69.4 • C, and T exp. m3 =91.2 • C taken from results of a differential scanning calorimetry experiment at pH=7 in Table 1 of Nixon and Giedroc [30] ).

The dependence of the stability of the F state with respect to the U state, ∆G UF (T ) is calculated using a procedure that does not rely on any order parameter [30, 31] . We also predict that ∆G UF (T = 25 • C) is −19.0 kcal/mol. The excellent agreement between simulations and experiments for the temperatureinduced structural transitions and the ∆G UF (T ) validates the model allowing us to predict the effect of f on BWYV PK.

In order to obtain the phase diagram in the [C, f ] plane of the BWYV PK, we performed a series of low friction Langevin dynamics simulations by varying the salt concentration from 5 to 1200 mM and the mechanical force from 0 to 20 pN, at a constant temperature of 50 • C. We varied C by changing the Debye length in the simulations, and the mechanical force was externally applied to the terminal nucleotides (Fig. 1A) .

Determination of the phase diagram in the [C, f ] plane requires an appropriate order parameter that can distinguish between the various states. In single-molecule pulling experiments, the variable that can be measured is the molecular extension, R ee , which is conjugated to the applied force. The overall size of the RNA molecule is assessed using the radius of gyration, R g . Therefore, we characterized the states of the RNA using both R g and R ee as order parameters. The values of the predicted R g at f = 0 can be measured using scattering experiments. Using these two parameters, we obtained the [C, f ] phase diagram (Fig. 4) . Comparison of the diagram of states in Figs. 4D and 4E reveals common features and some differences. From Figs. 4D and 4E, we infer that at f > 12.5 pN, extended (E) conformations are dominant at all values of C. As the force decreases, the PK forms compact structures. The boundary separating the extended and relatively compact phases depends on the salt concentration and the value of the mechanical force. The critical force to rupture the compact structures increases linearly as a function of logarithm of salt concentration (Fig. 4D ; boundary between red and green regions). At low forces (f < 2.5 pN), the diagram of states based on R ee shows that the extension is relatively small as C changes from a low to a high value. From this finding, one might be tempted to infer that the PK is always folded, which is not physical especially at low (C ≈ 10 mM) ion concentrations. In contrast, Fig. 4E shows that below 5 pN, there is a transition from compact structures (R g ≈ 1.3 nm in the blue region) at high C to an intermediate state (R g > 2.2 nm in the green region) at C ≈ 100 mM. The differences between the diagrams of states in the [C, f ] plane using R ee and R g as order parameters are more vividly illustrated in terms of the free energy profiles G(R α ) = −k B T log(R α ) where R α is R ee or R g (Fig. 5 ). The profiles G(R ee ), at three values of C and different f values, show two minima at most. At f = 0, there is only one discernible minimum at R ee ≈ 4.2 nm at C = 10 mM. The minimum moves to R ee ≈ 3 nm at C = 1200 mM corresponding to a folded PK. At f = 5 pN there are two minima at C = 10 mM corresponding to a compact structure and a stretched conformation (see the cyan profile in Fig. 5A ). As f exceeds 5 pN, there is essentially one minimum whose position increases as f increases. In the force regime (f > 5 pN), only the E state is visible at all C.

In order to compare the thermal free energy profiles generated at C = 500 mM (Fig. 2) , we calculated the G(R g ) and G(R ee ) at different forces with C fixed at 500 mM. Comparison of Fig. 2 and Fig. S3 shows that at f = 0 the free energy profiles are similar. At f = 0 the folded state is destabilized, and at f = 15 pN, the unfolded state is more stable than the folded PK.

A much richer diagram of states emerges when the G(R g ) = −k B T log P (R g ) is examined. The G(R g ) profile at C = 10 mM is similar to G(R ee ). However, at higher concentrations, there is clear signature of three states (especially at f = 12.5 pN) corresponding to the folded PK, an intermediate state, and the E state. The free energy profiles show that under force, there are minimally three states, which are F, I, and E. Because R ee cannot capture the complexity of the states populated as C and f are varied, it is important to consider R g as well to fully capture the phase diagram of the PK. We also computed the free energy profiles as a function of root-mean-square-deviations (RMSD). In terms of the number and relative positions of basins, the profiles based on RMSD are qualitatively similar to one of R g (Fig. S4 ).

In order to highlight the structural changes that occur as C and f are varied, we used, Q, the fraction of native HB interactions as an order parameter. The [C, f ] phase diagram (Fig. 6A ) calculated using the average Q is similar to the R g -phase diagram (Fig. 4E) , indicating the presence of three states (Fig. 6A ). Using Q, we can also quantitatively assess the contributions from different structural elements to the stability of the PK and correctly infer the order in which the structural elements of the PK rupture as f is increased. In order to determine the structural details of each state, we calculated the individual Qs for the two stems (Q S1 , Q S2 ) and the two loops (Q L1 Q L2 ) (Fig. 6B) . The dependence of Q S1 as a function of C and f shows that Stem 1 (S1) is extremely stable and remains intact at all salt concentrations, rupturing only at f ≈ 10 pN (Fig. 6B, upper left panel) . In contrast, the upper right panel in Fig. 6B shows that Stem 2 (S2) is only stably folded above a moderate salt concentration (C ≈ 80 mM) and low f . The stability difference between S1 and S2 can be explained by the number of canonical G-C base pairs; S1 has five G-C base pairs, whereas only three such pairs are in S2. Consequently, S1 cannot be destabilized even at low C, but its rupture to produce extended states requires the use of mechanical force. Above C > 100 mM, the fraction of contacts associated with S2 and the two loops (Q S2 , Q L1 and Q L2 ) simultaneously increase. All the interactions are, however, still not completely formed (Q ≈ 0.8) even at the highest salt concentration, C = 1200 mM. In particular, the contacts involving L2, Q L2 ≈ 0.6, implying that tertiary interactions are only marginally stable at T = 50 • C.

The results in Fig. 6 provide a structural explanation of the transitions that take place as C and f are changed.

(1) In the high force regime, the two stems are fully ruptured and native interactions are absent creating the E state. (2) Below C ≈ 100 mM and f 10 pN (left bottom on the [C, f ] plane), only S1 is intact, and thus, a hairpin conformation is dominant. This is the interme-diate state, I2, identified by Soto et al. who found that S1 remains stable even at salt concentration as less as 20 mM NaCl [31] . (3) Above C 100 mM and f 10 pN (right bottom), both S1 and S2 stems are formed, and the tertiary interactions involving the two loops are at least partly formed, ensuring that the topology of the PK is native-like. At f = 0 and temperature as a perturbation, it has been shown that another intermediate state, I1, exists between I2 and folded states [30, 31] . The I1 state may correspond to conformations in which two stems are fully folded with the entire topology of the native conformation intact, but tertiary interactions between stems and loops are incomplete. Although we have accounted for the presence of such a state in our simulations (see Fig. 2 ), from here on, we do not distinguish between I1 and the completely folded state in the phase diagram since the increases in Q S2 , Q L1 and Q L2 on the [C, f ] plane almost overlap. Therefore, we classify the right bottom region in the [C, f ] phase diagram as (I1 + F) state.

The [C, f ] phase diagrams provide a thermodynamic basis for predicting the structural changes that occur as f and C are varied. However, they capture only the average property of the ensemble within each state without the benefit of providing a molecular picture of the distinct states, which can only be described using simulations done at force values close to those used in LOT experiments. Therefore, we investigated the distribution of interactions to ascertain the factors that stabilize the three states as C and f are varied. Our analysis is based on the HB energy, U HB . The energy function for HB captures the formation of base pairs and loop-stem tertiary interactions, taking into account both distance and angular positions of the donor and the acceptor [37] . Fig. 7 shows that the probability densities of the total HB energy have three peaks at U HB ≈ 0, −30, and −60 kcal/mol. In the high-force regime (f > 12.5 pN), there is only one population at all C around U HB ≈ 0, which obviously corresponds to the fully extended state. The total HB energy can be decomposed into contributions arising from four distinct structural motifs in the PK (Fig. S5) . The peak at U HB ≈ −30 kcal/mol (Fig. 7) arises due to the formation of S1 (compare Figs. 7 and S5A). The broad peak at U HB ≈ −60 kcal/mol (Fig. 7A) is due to the sum of contributions from S2 and the tertiary interactions.

Parsing the U HB due to contributions from interactions among S1, S2, L1, and L2 produces a picture of thermodynamic ordering as C (f ) is increased (decreased) from low (high) value to high (low) value. Low C and high f correspond to the top left region in the phase diagrams (Figs. 4 and 6) . In all conditions, S1 is the first to form as indicated by the presence of distinct peaks (Figs. S5 B and C) . Only at the highest value of C and f 5 pN, tertiary interactions involving L2 emerge. Interestingly, the order of formation of the various structural elements en route to the formation of their folded state follows the stabilities of the individual elements with the most stable motif forming first. The PK tertiary interactions are consolidated only after the assembly of the stable structural elements. This finding follows from our earlier work establishing that the hierarchical assembly of RNA PKs is determined by the stabilities of individual structural elements [38] .

We did not find any condition in which S1 and S2 are completely folded, but there are no tertiary interactions, U L HB ≈ 0 (i.e., isolated I1 state). Our simulations show that the tertiary interactions contribute to the loop energy only after the stems are formed. Thus, the actual population is distributed between I1 and F, and the tertiary interactions are only marginally stable in the conditions examined here. Soto et al. suggested that diffuse Mg 2+ ions play an important role in stabilizing such tertiary interactions [31] . Because Mg 2+ is not considered here, we cannot determine the precise interactions stabilizing the I1 state, especially at f = 0. In the discussion hereafter, for simplicity, we refer to (I1 + F) state as F.

In order to delineate the phase boundaries quantitatively, we calculated the free energy differences between the three major states in the [C, f ] plane using ∆G EI (C, f ) = −k B T log P (I;C,f ) P (E;C,f ) , where the classification of the states, E, I (=I2), or F, is based on the HB interaction energies of stems, U S1 HB and U S2 HB (Fig. 8) . A similar equation is used to calculate ∆G IF . Based on the distribution of U HB , we determined the threshold values for S1 formation as U S1 HB < −15 kcal/mol (Fig. S5A) . Similarly, S2 formation is deemed to occur if U S2 HB < −10 kcal/mol (Fig. S5B) . We classified each structure depending on whether both stems are formed (F) with both U S1 HB < −15 kcal/mol and U S2 HB < −10 kcal/mol, only S1 is formed (I (=I2)), or fully extended (E).

The classification of the diagram of states based on the calculation of ∆G EI and ∆G IF is consistent with the R g and Q phase diagrams (Figs. 8C and 8D ). At all values of f , ∆G depends linearly on log C over a broad range of C as shown in Figs. 8A In what follows, we provide a theoretical interpretation of these results, which automatically shows that the difference in preferential ion interaction coefficients between states can be directly measured in LOT experiments [39] [40] [41] [42] [43] .

The salt concentration dependence of RNA stability can be determined using the preferential interaction coefficient, defined as Γ = ∂C ∂CRNA µ , where C RNA is the concentration of RNA, and the chemical potential of the ions, µ, is a constant [35, 44] . The free energy difference, ∆G αβ (=G β − G α ), between two states α and β (such as E and I or I and F) can be expressed as:

where C 0 is an arbitrary reference value of the salt concentration [45] . Note that we consider only 1:1 monovalent salt such as KCl or NaCl for which ion activity can be related to salt concentration C. The factor of 2 in Eq. (1) arises from charge neutrality. The difference, ∆Γ αβ (=Γ β −Γ α ), is interpreted as an effective number of bound or excluded ions when the RNA molecule changes its conformation between the two states [35, 44] . The free energy change upon application of an external force f is ∆G αβ (f ) = ∆G αβ (0)+∆R αβ ee ·f where ∆R αβ ee = R α ee −R β ee [29] . Thus,

In the [C, f ] phase diagram, ∆G αβ is 0 along the phase boundaries. Since the reference salt concentration C 0 is arbitrary, we determined its value using ∆G αβ (C 0 , f = 0) = 0. Thus, C 0 = C αβ m,f =0 is the salt concentration at the midpoint of the transition at zero force. The determination of ∆Γ αβ using this procedure from single-molecule pulling data has been described in several key papers [39] [40] [41] [42] . By adopting the procedure outlined in these studies and with our choice of C 0 , we rearrange Eq.(2) to obtain the phase boundary at f = 0,

where C αβ m is the midpoint of salt concentration, and f αβ c is the critical force associated with the transition between states α and β. Measurement of f αβ c using data from single-molecule pulling experiments has been used to obtain ∆Γ αβ for few systems that exhibit two-state transitions [39] [40] [41] [42] . Here, we have adopted it for a PK exhibiting more complex unfolding transitions. The connection between the relation in Eq. (3), derived elsewhere, to Clausius-Clapeyron equation was established recently by Saleh and coworkers [41] .

It follows from Eq.(3) that the critical force, The observed nonlinearity of the phase boundary separating I and F over the entire concentration range (Fig. 8D ) arises because ∆R IF ee is not a constant but varies linearly with f IF c (Figs. 9B and 9C). In giving from F↔I, there is only a rupture of few tertiary interactions and unfolding of the relatively unstable S2. Because this transition is not cooperative (interactions break in a sequential manner), we believe that a linear behavior over a limited force range is not unreasonable. In contrast, the unfolding of S1 is cooperative, like some hairpins studied using pulling experiments, in the E ↔ I transition, resulting in the constant value for ∆R ee over the force range probed. Thus, using ∆R αβ ee = af αβ c + b and Eq.(4), we find that f αβ c satisfies,

Note that if a is zero, then Eq.(5) reduces to Eq.(4) with b = ∆R αβ ee = const. Solving Eq.(5), the nonlinear dependence of f IF c is expressed as, 

The simulation data can be quantitatively fit using Eq.(6) (Fig. 8D) . In general, we expect ∆R αβ ee depends on f αβ c , and hence, f αβ c ∼ log C αβ m . It is worth noting that the estimation of ∆Γ αβ requires only the knowledge of how critical force varies with salt concentrations (∆G = 0). For any given salt concentration, the best statistics is obtained at the critical force because at f c , multiple transitions between the two states can be observed.

From the coefficients in the dependence of f c on log C m , we can estimate ∆Γ which provides a quantitative estimation of the extent of ion-RNA interactions [45] . For the E I transition, ∆R ee ≈ 4.9 nm and the slope of the linear function is 1.7 pN, which leads to ∆Γ EI = 0.96 (Fig. 8C ). This indicates that upon the conformational change from unfolded to the intermediate state, one ion pair enters the atmosphere of the PK. For the transition between I and F, fitting the dependence of ∆R IF ee on f IF c with a = 0.14 nm/pN and b = 0.27 nm (Fig. 9C) , we estimate ∆Γ IF = 1.8 using Eq.(6) (Fig. 8D) . Consequently, the total number of excess cations effectively bound to folded RNA is expected to be around 2.8, a value that is in rough accord with the experimental estimate of 2.2 [30] . More importantly, ∆Γ αβ can be measured using LOT and magnetic tweezer experiments [39] [40] [41] [42] [43] , thus expanding the scope of such experiments to measure ion-RNA interactions. A recent illustration of the efficacy of single-molecule experiments is the study by Jacobson and Saleh who quantified ∆Γ from force measurements for an RNA hairpin [43] .

Our analysis is based on the assumption that ∆Γ is a constant independent of the salt concentration. However, it is known experimentally that ∆Γ deviates from a constant value at high salt concentration [41, 43, 46] . In Fig. 8C and D, theoretical lines (blue dotted) show such a deviation from ∆G = 0 of the simulation data (white region) at high salt concentration (C 500 mM). Therefore, we believe that our theoretical analysis is most accurate for C 500 mM. Our coarse-grained simulation reflects the nonlinearity of ∆Γ observed in previous experiments. Understanding the origin of the nonlinearity requires a microscopic theory that should take into account the effects of ion-ion correlations. For complex architectures such as PK or ribozymes, this does not appear straightforward.

The three-interaction-site (TIS) model

We employed a variant of the three-interaction-site (TIS) model, a coarse-grained model first introduced by Hyeon and Thirumalai for simulating nucleic acids [34] . The TIS model has been previously used to make several quantitative predictions for RNA molecules, ranging from hairpins to ribozymes with particular focus on folding and response to f [34, 38, [47] [48] [49] . More recently, we have incorporated the consequences of counter ion condensation into the TIS model, allowing us to predict thermodynamic properties of RNA hairpins and PK that are in remarkable agreement with experiments [37] . In the TIS model, each nucleotide is represented by three coarse-grained spherical beads corresponding to phosphate (P), ribose sugar (S), and a base (B). Briefly, the effective potential energy (for details, see Ref. [37] ) of a given RNA conformation is U TIS = U L + U EV + U ST + U HB +U EL , where U L accounts for chain connectivity and angular rotation of the polynucleic acids, U EV accounts for excluded volume interactions of each chemical group, and U ST and U HB are the base-stacking and HB interactions, respectively. Finally, the term U EL corresponds to electrostatic interactions between phosphate groups.

The determination of the parameters for use in simulations of the coarse-grained model of RNA was described in detail in a previous publication [37] . Briefly, we determined the stacking interactions by a learning procedure, which amounts to reproducing the measured melting temperature of dimers. We showed that a single choice of model parameters (stacking interactions, the Debye-Hückel potential for electrostatic interactions, and structure-specific choice for hydrogen bonds) is sufficient to obtain detailed agreements with available experimental data for three different RNA molecules. Thus, the parameters are transferable and we use them verbatim in the present simulations. The determination of hydrogen bond interaction is predicated on the structure. Most of these are taken from the A-form RNA structure and hence can be readily used for any RNA molecule. Only the interaction lengths and angles in hydrogen-bonding interaction of noncanonical base pairs are based on the specific PDB structure, which in our case is the BWYV pseudoknot in this study. We list all the BYWV-PKspecific values used in our simulations in Table S1 .

The repulsive electrostatic interactions between the phosphate groups is taken into account through the Debye-Hückel theory, U EL = i,j q * 2 e 2 4πε0ε(T )rij exp − rij λD , where the Debye length is λ D = ε(T )kBT 4πe 2 I . In the simulations, salt concentration can be varied by changing the ionic strength I = Σq 2 n ρ n , where q n is the charge of ion of type n, and ρ n is its number density. Following our earlier study [37] , we use an experimentally fitted function for the temperature-dependent dielectric constant ε(T ) [50] . The renormalized charge on the phosphate group is −q * e (q * < 1). Because of the highly charged nature of the polyanion, counterions condense onto the RNA, thus effectively reducing the effective charge per phosphate. The extent of charge reduction can be calculated using the Oosawa-Manning theory. Charge renormalization, first used in the context of RNA folding by Heilman-Miller et al. [51] and more recently incorporated into CG models for RNA by us [37] and others [52, 53] , is needed to predict accurately the thermodynamics and kinetics of RNA folding [54] . The renormalized value of the charge on the P group is approximately given by −q * (T )e = −be lB(T ) , where the Bjerrum length is l B (T ) = e 2 ε(T )kBT and b is the mean distance between the charges on the phosphate groups [55] . The mean distance (b) between charges on RNA is difficult to estimate (except for rod-like polyelectrolytes) because of its shape fluctuations. Therefore, it should be treated as an adjustable parameter. In our previous study, we showed that b = 4.4Å provides an excellent fit of the thermodynamics of RNA hairpins and the MMTV pseudoknot [37] . We use the same value in the present study of BWYV as well. Thus, the force field used in this study is the same as in our previous study attesting to its transferability for simulating small RNA molecules.

The native structure of the BWYV PK is taken from the PDB (437D) [32] . The structure contains 28 nucleotides. We added an additional guanosine monophosphate to both the 5 and 3 terminus. Thus, the simulated PK has 30 nucleotides. The mechanical force is applied to the ribose sugars of both the 5 -and 3 -terminus (Fig. 1A) .

We performed Langevin dynamics simulations by solving the equation of motion, mẍ = − ∂UTIS ∂x − γẋ + R, where m is mass; γ, the friction coefficient, depends on the size of the particle type (P, S, or B); and R is a Gaussian random force, which satisfies the fluctuationdissipation relation. The [C, f ] phase diagram is calculated at a somewhat elevated temperature, T =50 • C, in order to observe multiple folding-unfolding transitions when the monovalent salt concentration is varied over a wide range, from 5 to 1200 mM. The numerical integration is performed using the leap-frog algorithm with time step length δt L = 0.05τ , where τ is the unit of time. We performed low friction Langevin dynamics simulations to enhance the rate of conformational sampling [56] . For each condition (C and f ), we performed 50 independent simulations for 6 × 10 8 time steps and the last two-third of trajectories are used for data analyses. In addition, we also performed a set of simulations using the replicaexchange molecular dynamics (REMD) to confirm the validity of the results obtained by the low friction Langevin dynamics simulations [57] . In the REMD simulations, both C and f are treated as replica variables and are exchanged between 256 replicas, by covering C from 1 to 1200 mM and f from 0 to 20 pN. The replicas are distributed over 16 × 16 grid. The combined set of simulations using different protocols assures us that the simulations are converged at all conditions. Our simulations are done at constant forces under equilibrium conditions in which we observe multiple transitions between the various states. We have made comparison to LOT experiments, which are performed by pulling the RNA at a constant velocity. In general, if RNA is stretched at a constant velocity, then one has to be concerned about non-equilibrium effects. However, in optical tweezer experiments, the pulling speed is low enough so that the tension propagates uniformly throughout the RNA prior to initiation of unfolding [58] . Thus, the unfolding of PK studied using LOT experiments [19] occurs at equilibrium, justifying comparisons to our explicit equilibrium simulations.

To analyze the structural transitions from which the phase diagrams are determined, we stored snapshots every 10 4 time steps from all the trajectories. We computed the radius of gyration (R g ) and the end-to-end distance (R ee ) using the stored snapshots. The fraction of native contacts, Q, is calculated by counting the number of hydrogen bond (HB) interaction pairs. The assessment of HB formation is based on the instantaneous value of the HB interaction energy, U HB . In our model, each HB interaction pair contributes up to nU 0 HB toward stability, where U 0 HB is -2.43 kcal/mol, corresponding to the stability of one HB, and n represents number of HB associated in the interaction [37] . A cutoff value nU C HB is defined to pick out contacts that are established. We use a value, U C HB = −1.15 kcal/mol to obtain the diagram using Q as an order parameter. Modest changes in U C HB do not affect the results qualitatively.

Motivated by the relevance of force-induced transitions in PK to PRF, we have conducted extensive simulations of BWYV PK using a model of RNA, which predicts with near quantitative accuracy the thermodynamics at f = 0. The phase diagram in the [C, f ] plane, which can be measured using single-molecule LOT pulling experiments, shows that the thermodynamics of rupture as f is increased, or the assembly of the PK as f is decreased, involves transitions between the extended, intermediate, and the folded states. The predicted linear relationship between f c and log C m [Eq.4] or f c ∼ √ log C m [Eq.6] shows that ∆Γ can be directly measured in LOT experiments. The theoretical arguments leading to Eqs.(4) and (6) are general. Thus, f can be used as a perturbant to probe ion-RNA interactions through direct measurement of ∆Γ.

This work was completed when the authors were at the University of Maryland. We are grateful to Jon Dinman, Xin Li, Pavel Zhuravlev, Huong Vu, and Mauro Mugnai for comments on the manuscript. This work was supported in part by a grant from the National Science Foundation (CHE-1636424). S1. List of hydrogen bonding interactions in BWYV PK. First 8 pairs correspond canonical base pairs forming stems S1 and S2. Other 15 pairs are all tertially interactions forming either L1 or L2. Since we use the three-interaction-site model, a letter "s", "b" or "p" is put following each nucleotide number to indicate one of sugar, base, or phosphate, respectively.

Nucleotide pair Interaction type (atoms)

The noncoding RNA revolution-trashing old rules to forge new ones

How RNA folds

Structure and assembly of group I introns

RNA and protein folding: common themes and variations

RNA folding: conformational statistics, folding kinetics, and ion electrostatics

RNA folding pathways and the selfassembly of ribosomes

Determination of thermodynamics and kinetics of RNA reactions by force

Structure and function of telomerase RNA

Characterization of an efficient coronavirus ribosomal frameshifting signal: requirement for an RNA pseudoknot

A functional pseudoknot in 16S ribosomal RNA

Structural and functional aspects of RNA pseudoknots

Pseudoknots: RNA structures with diverse functions

The telomerase RNA pseudoknot is critical for the stable assembly of a catalytically active ribonucleoprotein

Characterization of the mechanical unfolding of RNA pseudoknots

Programmed −1 frameshifting efficiency correlates with RNA pseudoknot conformational plasticity, not resistance to mechanical unfolding

A frameshifting stimulatory stem loop destabilizes the hybrid state and impedes ribosomal translocation

Single-molecule mechanical unfolding and folding of a pseudoknot in human telomerase RNA

Frameshifting RNA pseudoknots: structure and mechanism

Mechanical unfolding of the beet western yellow virus -1 frameshift signal

Mechanisms and implications of programmed translational frameshifting

Triplex structures in an RNA pseudoknot enhance mechanical stability and increase efficiency of -1 ribosomal frameshifting

Singlemolecule measurements of the CCR5 mRNA unfolding pathways

The role of RNA pseudoknot stem 1 length in the promotion of efficient -1 ribosomal frameshifting

Ribosomal pausing at a frameshifter RNA pseudoknot is sensitive to reading phase but shows little correlation with frameshift efficiency

Specific mutations in a viral RNA pseudoknot drastically change ribosomal frameshifting efficiency

HIV-1 frameshift efficiency is primarily determined by the stability of base pairs positioned at the mRNA entrance channel of the ribosome

Ribosomal frameshifting in the CCR5 mRNA is regulated by miRNAs and the NMD pathway

Transactivation of programmed ribosomal frameshifting by a viral protein

The effect of force on thermodynamics and kinetics of single molecule reactions

Energetics of a strongly pH dependent RNA tertiary structure in a frameshifting pseudoknot

Tertiary structure of an RNA pseudoknot is stabilized by "diffuse" Mg 2+ ions

Minor groove RNA triplex in the crystal structure of a ribosomal frameshifting viral pseudoknot

Salt contribution to RNA tertiary structure folding stability

Mechanical unfolding of RNA hairpins

Salt dependence of oligoion-polyion binding: a thermodynamic description based on preferential interaction coefficients

Direct measurement of the mechanical work during translocation by the ribosome

Coarse-grained model for predicting RNA folding thermodynamics

Assembly mechanisms of RNA pseudoknots are determined by the stabilities of constituent secondary structures

Interplay of ion binding and attraction in DNA condensed by multivalent cations

Maxwell relations for single-DNA experiments: Monitoring protein binding and double-helix torque with force-extension measurements

Single-molecule methods for ligand counting: Linking ion uptake to DNA hairpin folding

Measuring the differential stoichiometry and energetics of ligand binding to macromolecules by single-molecule force spectroscopy: an extended theory

Quantifying the ion atmosphere of unfolded, single-stranded nucleic acids using equilibrium dialysis and single-molecule methods

Conformational transitions of duplex and triplex nucleic acid helices: thermodynamic analysis of effects of salt concentration on stability using preferential interaction coefficients

Analysis of effects of salts and uncharged solutes on protein and nucleic acid equilibria and processes: a practical guide to recognizing and interpreting polyelectrolyte effects, Hofmeister effects, and osmotic effects of salts

Thermodynamic analysis of ion effects on the binding and conformational equilibria of proteins and nucleic acids: the roles of ion association or release, screening, and ion effects on water activity

Folding of human telomerase RNA pseudoknot using ion-jump and temperature-quench simulations

Kinetics of allosteric transitions in S-adenosylmethionine riboswitch are accurately predicted from the folding landscape

How do metal ions direct ribozyme folding?

Dielectric constant of water from 0 • to 100 • C

Role of counterion condensation in folding of the Tetrahymena ribozyme. I. Equilibrium stabilization by cations

Reduced model captures Mg 2+ -RNA interaction free energy of riboswitches

Generalized Manning condensation model captures the RNA ion atmosphere

Role of counterion condensation in folding of the Tetrahymena ribozyme II. Counterion-dependence of folding kinetics

Limiting Laws and Counterion Condensation in Polyelectrolyte Solutions I

The nature of folded states of globular proteins

Replica-exchange molecular dynamics method for protein folding

Pathways and kinetic barriers in mechanical unfolding and refolding of RNA and proteins

Free energy profiles of the hydrogen bond energy of two stems (A; U S HB ) and two loops (stem-loop interactions

C = 500 mM and f = 0

Probability distributions of the molecular extension (A

Ree) and radius of gyration (B; Rg) at C = 500 mM and f = 0. FIG. S3. Free energy profiles as a function of Ree (A) and Rg (B) at 500 mM of monovalent salt

Free energy profiles as a function of RMSD at four different salt concentrations, 10 mM (A), 100 mM (B), 500 mM (C), and 1200 mM (D)

Probability distributions of hydrogen-bond energy. (A-D) Distributions of hydrogen bond energies for individual structural elements (A) Stem1; (B) Stem 2; (C) Loop 1; and (D) Loop2. The salt concentrations and force values are explicitly

G4b -A20b N2 -N3

G4s -A20b O2 -N1

C5b -A20s O2 -O2

G7b -A24b N2 -N1

C8b -G12b N4 -N7

C8b -C26b O2 -N4

C14s -A25b O2 -N1

C15s -A23b O2 -N1