Originally published as Biophys J. BioFAST on February 10, 2006.
doi:10.1529/biophysj.105.078030
Biophysical Journal 90:3410-3427 (2006)
© 2006 The Biophysical Society
Forced-Unfolding and Force-Quench Refolding of RNA Hairpins
Changbong Hyeon * and
D. Thirumalai * 
* Biophysics Program Institute for Physical Science and Technology and Department of Chemistry and Biochemistry, University of Maryland, College Park, Maryland
Correspondence: Address reprint requests to D. Thirumalai, Tel.: 301-405-4803; E-mail: thirum{at}glue.umd.edu.
 |
ABSTRACT
|
|---|
Nanomanipulation of individual RNA molecules, using laser optical tweezers, has made it possible to infer the major features of their energy landscape. Time-dependent mechanical unfolding trajectories, measured at a constant stretching force (fS) of simple RNA structures (hairpins and three-helix junctions) sandwiched between RNA/DNA hybrid handles show that they unfold in a reversible all-or-none manner. To provide a molecular interpretation of the experiments we use a general coarse-grained off-lattice G
-like model, in which each nucleotide is represented using three interaction sites. Using the coarse-grained model we have explored forced-unfolding of RNA hairpin as a function of fS and the loading rate (rf). The simulations and theoretical analysis have been done both with and without the handles that are explicitly modeled by semiflexible polymer chains. The mechanisms and timescales for denaturation by temperature jump and mechanical unfolding are vastly different. The directed perturbation of the native state by fS results in a sequential unfolding of the hairpin starting from their ends, whereas thermal denaturation occurs stochastically. From the dependence of the unfolding rates on rf and fS we show that the position of the unfolding transition state is not a constant but moves dramatically as either rf or fS is changed. The transition-state movements are interpreted by adopting the Hammond postulate for forced-unfolding. Forced-unfolding simulations of RNA, with handles attached to the two ends, show that the value of the unfolding force increases (especially at high pulling speeds) as the length of the handles increases. The pathways for refolding of RNA from stretched initial conformation, upon quenching fS to the quench force fQ, are highly heterogeneous. The refolding times, upon force-quench, are at least an order-of-magnitude greater than those obtained by temperature-quench. The long fQ-dependent refolding times starting from fully stretched states are analyzed using a model that accounts for the microscopic steps in the rate-limiting step, which involves the trans to gauche transitions of the dihedral angles in the GAAA tetraloop. The simulations with explicit molecular model for the handles show that the dynamics of force-quench refolding is strongly dependent on the interplay of their contour length and persistence length and the RNA persistence length. Using the generality of our results, we also make a number of precise experimentally testable predictions.
|
INTRODUCTION
|
|---|
RNA molecules adopt precisely defined three-dimensional structures to perform specific functions (1
). To reveal the folding pathways navigated by RNA en route to their native conformations requires exhaustive exploration of the complex underlying energy landscape over a wide range of external conditions. In recent years, mechanical force has been used to probe the unfolding of a number of RNA molecules (2–4). Force is a novel way of probing regions of the energy landscape that cannot be accessed by conventional methods (temperature changes or variations in counterion concentrations). In addition, response of RNA to force is relevant in a number of cellular processes such as mRNA translocation through the ribosome and the activity of RNA-dependent RNA polymerases. Indeed, many dynamical processes are controlled by deformation of biomolecules by mechanical force.
By exploiting the ability of single molecule laser optical tweezer (LOT) to control the magnitude of the applied force, Bustamante and co-workers have generated mechanical unfolding trajectories for RNA hairpins and the Tetrahymena thermophila ribozyme (5,6). The unfolding of the ribozyme shows multiple routes with great heterogeneity in the unfolding pathways (6). In their first study (5), they showed that simple RNA constructs (P5ab RNA hairpins or a three-helix junction) unfold reversibly at equilibrium. From the time traces of the end-to-end distance (R) of P5ab, for a number of force values, Liphardt et al. (5) showed that the hairpins unfold in a two-state manner. The histograms of time-dependent R (and assuming ergodicity) were used to calculate the free energy difference between the folded and unfolded states. Unfolding kinetics, as a function of the stretching force fS, was used to identify the position of the transition states (5,7,8). These experiments and subsequent studies have established force as a viable way of quantitatively probing the RNA energy landscape with R serving as a suitable reaction coordinate.
The experiments by Bustamante and co-workers have led to a number of theoretical and computational studies using a variety of different methods (9–14). These studies have provided additional insights into the mechanical unfolding of RNA hairpins and ribozymes. In this article, we build on our previous work (14) and new theoretical analysis to address a number of questions that pertain to mechanical unfolding of RNA hairpins. In addition, we also provide the first report on force-quench refolding of RNA hairpins. This article addresses the following major questions:
- Are there differences in the mechanism of thermal and mechanical unfolding? We expect these two processes to proceed by different pathways because denaturation induced by temperature jump results in a stochastic perturbation of the native state while destabilization by force is due to a directed perturbation. The molecular model for P5GA gives a microscopic picture of these profound differences.
- For a given sequence, does the position of the transition state move in response to changes in the loading rate (rf) or the stretching force fS? Based on the analysis of experiments over a narrow range of conditions (fixed temperature and loading rate), it has been suggested that the location of the sequence-dependent unfolding transition state (TS) for secondary structure is midway between the folded state while the TS for the ribozyme is close to the native conformation (5). Explicit simulations show that the TS moves dramatically, especially at high values of fS and rf. As a consequence, the dependence of the unfolding rate on fS deviates from the predictions of the Bell model (15).
- What is the origin of the dramatic differences in refolding by force-quench from stretched conformations and temperature-quench refolding? Experiments by Fernandez and Li (8) on refolding initiated by force-quench on polyubiquitin construct suggest similar differences in the refolding timescales. For RNA hairpins, we show that the incompatibility of the local loop structures in the stretched state and the folded conformations lead to extremely long refolding times upon force-quench.
- What are the effects of linker dynamics on forced-unfolding and force-quench refolding of RNA hairpins? The manipulation of RNA is done by attaching a handle or linker to the 3' and 5' ends of RNA. Linkers in the LOT experiments, which are done under near-equilibrium conditions, are RNA/DNA hybrid handles. These are appropriately modeled as wormlike chains (WLCs). By adopting an explicit polymer model for semiflexible chains we show that, under certain circumstances, nonequilibrium response of the handle (which is not relevant for the LOT experiments) can alter the forced-unfolding dynamics of RNA. We probe the effects of varying the linker characteristics on the dynamics of folding/unfolding of RNA for the model of the handle-RNA-handle construct. In certain ranges of rf, nonequilibrium effects on the dynamics of linkers can affect the force-extension profiles.
|
METHODS
|
|---|
Hairpin sequence
To probe the dynamics of unfolding and force-quench refolding we have studied in detail the 22-nucleotide hairpin, P5GA, whose solution structure has been determined by NMR (Protein Data Bank (PDB) id:1eor). In many respects P5GA is similar to P5ab in the P5abc domain of group I intron (16). Both of these structures have GA mismatches, and are characterized by the presence of the GAAA tetraloop. The sequence of P5GA is GGCGAAGUCGAAAGAUGGCGCC (17).
RNA model
Because it is difficult to explore the unfolding and refolding of RNA hairpins over a wide range of external conditions such as temperature (T), stretching force (fS), and quench force (fQ) using all-atom models of biomolecules in explicit water, we have introduced a minimal off-lattice coarse-grained model (14). (Note that throughout the article, while fS and fQ each have a specific meaning, we use the notation f for referring to mechanical force in general terms.) In this model, each nucleotide is represented by three beads with interaction sites corresponding to the phosphate (P) group, the ribose (S) group, and the base (B) (14). Thus, the RNA backbone is reduced to the polymeric structure ((P – S)n) with the base that is covalently attached to the ribose center. In the minimal model, the RNA molecule with N nucleotides corresponds to 3N interaction centers. The secondary structure and the lowest energy structure using the minimal model are shown in Fig. 1.
View larger version (33K):
[in this window]
[in a new window]
|
FIGURE 1 Coarse-grained representation of a molecule of RNA using three interaction sites per nucleotide—phosphate (P), sugar (S), and base (B). On the left we present the secondary structure of the 22-nt P5GA hairpin, in which the bonds formed between basepairs are labeled from 1 to 9. The PDB structure (17) and the lowest energy structure obtained with the coarse-grained model are shown on the right.
|
|
Energy function
The total energy of a RNA conformation, specified by the coordinates of the 3N sites, is written as VT = VBL + VBA + VDIH + VSTACK + VNON + VELEC. Harmonic potentials are used to enforce structural connectivity and backbone rigidity. The connectivity between two beads (PiSi, SiPi+1, and BiSi) is maintained using
 | (1) |
where kr = 20 kcal/(mol x Å2), and (
)i and (
)i are the distances between covalently bonded beads in the PDB structure. The notation (SP)i, denotes the ith backbone bead S or P. The angle 
formed between three successive beads (Pi – Si – Pi+1 or Si–1 – Pi – Si) along sugar-phosphate backbone is subject to the bond-angle potential,
 | (2) |
where k = 20 kcal/(mol x rad2), and 
is the value in the PDB structure.
Dihedral angle potential
The dihedral angle potential (VDIH) describes the ease of rotation around the angle formed between four successive beads along the sugar-phosphate backbone (Si–1PiSiPi+1 or PiSiPi+1Si+1). The ith dihedral angle 
i, which is the angle formed between the two planes defined by four successive beads i to i + 3, is defined by 
. In the coarse-grained model, the right-handed chirality of RNA is realized by appropriate choices of the parameters in the dihedral potential. Based on the angles in the PDB structure (
), one of the three types of dihedral potentials, trans (t, 0 < < 2
/3), gauche (+) (g+, 2/3 < < 4/3), and gauche (–) (g–, 4/3 < < 2), is assigned to each of the four successive beads along the backbone. The total dihedral potential of the hairpin is
 | (3) |
where the parameters (in kcal/mol) defined for t, g+, and g– are
To account for the flexibility in the loop region we reduce the dihedral angle barrier by halving the parameter values in 19 
i 24.
Stacking interactions
Simple RNA secondary structures, such as hairpins, are largely stabilized by stacking interactions whose context-dependent values are known (18–20). The folded P5GA RNA hairpin is stabilized by nine hydrogen bonds between the basepairs (see Fig. 1 B), including two GA mismatch pairs (17). The stacking interactions that stabilize a hairpin can be written as 
(nmax = 8 in P5GA). We assume that the orientation-dependent Vi is
 | (4) |
where: 
G(T) = H – TS; the bond angles {} are 1i 
SiBiBj, 2i BiBjSj, 3i Si+1Bi+1Bj–1, and 4i Bi+1Bj–1Sj–1; the distance between two paired bases is rij = |Bi – Bj| and ri+1j–1 = |Bi+1 – Bj–1|; and 
1i and 2i are the dihedral angles formed by the four beads BiSiSi+1Bi+1 and Bj–1Sj–1SjBj, respectively. The superscript o refers to angles and distances in the PDB structure. The values of 
st, ßst, and 
st are 1.0, 0.3 Å–2, and 1.0, respectively. We take H and S from Turner's thermodynamic data set (18,19). There are no estimates for GA-related stacking interactions. Nucleotides G and A do not typically form a stable bond, and hence GA pairing is considered a mismatch. We use the energy associated with GU for the GA basepair.
Nonbonded interactions
We use the Lennard-Jones interactions between nonbonded interaction centers to account for the hydrophobicity of the purine/pyrimidine groups. The total nonbonded potential is
 | (5) |
where 
. The prime in the second term on the Eq. 5 denotes the condition m 
2i – 1. In our model, a native contact exists between two noncovalently bound beads, provided they are within a cutoff distance rc (= 7.0 Å). Two beads beyond rc are considered to be nonnative. For a native contact,
 | (6) |
where 
is the distance between beads in PDB structure and 
for all native contact pairs, except for the B10B13 basepair associated with the formation of the hairpin loop, for which 
. The additional stability for the basepair associated with loop formation is similar to the Turner's thermodynamic rule for the free energy gain in the tetraloop region. For beads beyond rc the interaction is
 | (7) |
with a = 3.4 Å and CR = 1 kcal/mol. The value of 
(= 1.5 kcal/mol) has been chosen so that the hairpin undergoes a first-order transition from unfolded states. Our results are not sensitive to minor variations in Ch.
Electrostatic interactions
The charges on the phosphate groups are efficiently screened by counterions so that in the folded state the destabilizing contribution of the electrostatic potential is relatively small. However, the nature of the RNA conformation (especially tertiary interactions) can be modulated by changing counterions. For simplicity, we assume that the electrostatic potential between the phosphate groups is pairwise-additive 
For 
we use Debye-Hückel potential, which accounts for screening by condensed counterions and hydration effects, and is given by
 | (8) |
where 
is the charge on the phosphate ion, 
r = /0, and the Debye length is 
with 
= 8.99 x 109 JmC. To calculate the ionic strength 
, we use the value ci = 200 mM-NaCl from the header of a PDB file (17). We use r = 10 in the simulation (21). Because of the Debye screening length, 
, the strength of electrostatic interactions between the phosphate groups are temperature-dependent even when we ignore the variations of with T. At room temperature (T 
300 K), the electrostatic repulsion between the phosphate groups at r 5.8 Å, which is the closest distance between phosphate groups, is 
Thus, VELEC between phosphate groups across the basepairing (r = 16 18 Å) is almost negligible. The Debye-Hückel interactions are most appropriate for monovalent counterions like Na+.
Models for linker or handle
In laser optical tweezer (LOT) experiments, the RNA molecules are attached to the polystyrene beads by an RNA/DNA hybrid handle or linker. A schematic illustration of the pulling simulations that mimic the experimental setups for the LOT (Fig. 2 A) is shown in Fig. 2 B. Globally, the principles involved in atomic force microscopy (AFM) and LOT for mechanical unfolding of biomolecules are essentially the same, except for the difference in the effective spring constant. The spring constant of the nearly harmonic potential of the optical trap is in the range of k = 0.01–0.1 pN/nm, whereas the cantilever spring in AFM experiments (Fig. 2 B) is much stiffer and varies from 1 to 10 pN/nm. To fully characterize the RNA energy landscape, it is necessary to explore a wide range of loading rates (22). To vary rf we have used different values of k in the simulations.
View larger version (40K):
[in this window]
[in a new window]
|
FIGURE 2 (A) Schematic illustration of laser optical tweezer (LOT) setup for RNA stretching. A single RNA molecule is held between two polystyrene beads via molecular handles with one of the polystyrene beads being optically trapped in the laser light. The location of the other bead is changed by micropipette manipulation. The extension of the molecule through the molecular handles induces the deviation in the position of the polystyrene bead held in the force-measuring optical trap. (B) Both LOT and AFM can be conceptualized as schematically shown. The RNA molecule is sandwiched between the linkers and one end of the linker is pulled. The spring constant of the harmonic trap in LOT (or the cantilever in AFM experiments) is given by k and v is the pulling speed.
|
|
We simulate mechanical unfolding of RNA hairpins either by applying a constant force on the 3'-end with the 5'-end being fixed (unfolding at constant force), or by pulling the 3'-end through a combination of linkers and harmonic spring (Fig. 2 B) at a constant speed in one direction (unfolding at constant loading rate). Comparison of the results allows us to test the role of the linker dynamics on the experimental outcome. If the linker is sufficiently stiff then it should not affect the dynamics of RNA. On the other hand, unfolding at constant loading rate (rf k x v, where v is the pulling speed and the stretching force (fS) is computed using f = k x 
z with z being the displacement of the spring) can be modulated either by changing k or v. Instead of k, an effective spring constant keff, with 
, should be used to compute the loading rate. Typically, 
and 
(where kmol is the stiffness associated with the model) are small, and hence, keff 
k. In this setup, the relevant variables are k, v, and the contour length (L) of the linker and the effective stiffness of the linker. We explore the effects of these factors by probing changes in the force extension curve (FEC). The variations in L are explored only using constant loading rate simulations. Manosas and Ritort (13) used an approximate method to model linker dynamics.
The energy function for the linker molecule is
 | (9) |
where ri, i+1 and 
are the distance and unit vector connecting i and i + 1 residue, respectively. For the bond potential we set kB = 20 kcal/(mol · Å2) and b = 5 Å. This form of the energy function describes the wormlike chain (WLC) (23) (appropriate for RNA/DNA hybrid handles used by Liphardt et al. (5)) when kA is large. The linker becomes more flexible as the parameter describing the bending potential, kA, is reduced. Thus, by varying kA the changes in the entropic elasticity of the linker on RNA hairpin dynamics can be examined. We use kA = 80 kcal/mol or 20 kcal/mol to change the flexibility of the linker. For the purpose of computational efficiency, we did not include excluded volume interactions between linkers or between the linker and RNA. When linkers are under tension, the chains do not cross unless thermal fluctuations are larger than the energies associated with force. To study the linker effect on force extension curves (FEC), we attach two linker polymers, each with the contour length L/2, to the ends of the hairpin and stretch the molecule using the single pulling speed 0.86 x 102 µm/s with spring constant k = 0.7 pN/nm. The total linker length is varied from L = (10–50) nm.
Simulations
We assume that the dynamics of the molecules (RNA hairpins and the linkers) can be described by the Langevin equation. The system of Langevin equations is integrated as described before (14,24). Using typical values for the mass of a bead in a nucleotide (Bi, Si, or Pi), m = 100 g/mol 160 g/mol, the average distance between the adjacent beads a = 4.6 Å, the energy scale h = 1 2 kcal/mol, the natural time is 
= 2.6 2.8 ps. We use 
L = 2.0 ps to convert the simulation times into real times. To estimate the timescale for thermal and mechanical unfolding dynamics, we use a Brownian dynamics algorithm (25,26) for which the natural time for the overdamped motion is 
. We use 
= 50 
which approximately corresponds to friction constant in water, in the overdamped limit. To probe the thermodynamics and kinetics of folding we used a number of physical quantities (end-to-end distance (R), fraction of native contacts (Q), the structural overlap function (
), number of hydrogen bonds nbond, etc.) to monitor the structural change in the hairpin.
Computation of free energy profiles
We adopted the multiple histogram technique (27,28) to compute the thermodynamic averages of all the observables at various values of T and f. For example, the thermodynamic average of the fraction of native contact, Q, can be obtained at arbitrary values of T and f if the conformational states are well sampled over a range of T- and f-values. The thermodynamic average of Q is given by
 | (10) |
where K is the number of histograms, hk(E, R, Q) is the number of states between E and E + E, R and R + R, Q and Q + Q in the kth histogram, 
, and Tk and fk are the temperature and the force in the simulations used to generate the kth histogram, respectively. The free energy, Fk, which is calculated self-consistently, satisfies
 | (11) |
Using the low friction Langevin dynamics, we sampled the conformational states in the (T,f) in the range {0 K < T < 500 K, f = 0.0 pN} and {0.0 pN < f < 20.0 pN, T = 305 K}. Exhaustive samplings around the transition regions at {305 K T 356 K, f = 0.0 pN} and {5.0 pN f 7.0 pN, T = 305 K} is required to obtain reliable estimates of the thermodynamic quantities. The free energy profile, with Q as an order parameter, is given by
 | (12) |
where Fo(T, f) = –kBT log Z(T, f); 
; and P(Q(T, f)), as defined in Eq. 10. The free energy profile F(R) with R as a reaction coordinate can be obtained using a similar expression.
|
RESULTS
|
|---|
Mechanisms of thermal denaturation and forced-unfolding are different
We had previously reported (14) the thermodynamic characteristics of the P5GA hairpin as a function of T and f. The native structure of the hairpin, determined using a combination of multiple slow cooling, simulated annealing, and steepest-descent quench, yielded conformations whose average root mean-square deviation with respect to the PDB structure is 0.1 Å (14). The use of a G-like model leads to a small value of root mean-square deviation. From the equilibrium (T,f) phase diagram in Hyeon and Thirumalai (14), the melting temperature, Tm 341 K, at fS = 0. Just as in thermal denaturation at Tm at zero f, the hairpin unfolds by a weak first-order transition at an equilibrium critical force fc. Above fc, which is temperature-dependent, the folded state is unstable.
To monitor the pathways explored in the thermal denaturation, we initially equilibrated the conformations at T = 100 K at which the hairpin is stable. Thermal unfolding (f = 0) was initiated by a temperature jump to T = 346 K > Tm. Similarly, forced-unfolding is induced by applying a constant force fS = 42 pN to thermally equilibrated initial conformation at T = 254 K (14). The value of fS = 42 pN far exceeds fc = 15 pN at T = 254 K. Upon thermal denaturation, the nine bonds fluctuate (in time) stochastically, in a manner independent of one another until the hairpin melts (Fig. 3 A). Forced-unfolding, on the other hand, occurs in a directed manner. Mechanical unfolding occurs by sequential unzipping with force unfolding the bond, from the ends of the hairpin (beginning at bond 1) to the loop (Fig. 3 B).
The differences in the folding pathways are also mirrored in the free energy profiles. Assuming that Q is an adequate reaction coordinate for thermal unfolding (29) we find, by comparing F(Q) at T = 100 K and T = 346 K, that thermal melting occurs by crossing a barrier. The native basin of attraction (NBA) at Q = 0.9 at T = 100 K is unstable at the higher temperature and the new equilibrium at Q 0.2 is reached in an apparent two-state manner (Fig. 3 C). Upon directed mechanical unfolding, the free energy profile F(R) is tilted from the NBA at R 1.5 nm to R 12 nm, at which the stretched states are favored at f = 42 pN. The forced-unfolding transition also occurs abruptly once the activation barrier at R 1.5 nm (close to the folded state) is crossed (Fig. 3 D).
Free energy profiles and transition state (TS) movements
Based on the proximity of the average transition state location, 
, it has been suggested that folded states of RNA (6) and proteins (7) are brittle. If the experiments are performed by stretching at a constant loading rate, then is calculated using 
log rf (30), where f* is the most probable unfolding force and rf, the loading rate, is rf = df/dt = kv. Substantial curvatures in the dependence of f* on log rf ((f*, log rf) plot) have been observed, especially if rf is varied over a wide range (22). Similarly, in constant force unfolding experiments, is obtained from the Bell equation (15), which relates the unfolding rate to the applied force, 
where 
is the unfolding rate in the absence of force. In the presence of curvature in the (f*, log rf) plots or when the Bell relationship is violated (31), it is difficult to extract meaningful values of or by a simple linear extrapolation. By carefully examining the origin of curvature in the (f*, log rf) plot or in kU, we show that in the unfolding of hairpin the observed nonlinearities are due to movements in the transition state ensemble, i.e., depends on fS and rf.
Unfolding at constant loading rate
We performed forced unfolding simulations by varying both the pulling speed v and the spring constant k so that a broad range of loading rates can be covered. The unfolding forces at which all the hydrogen bonds are ruptured are broadly distributed, with the average and the dispersion that increase with growing loading rates (Fig. 4 A). The plot of f* as a function of log rf (Fig. 4 B) shows marked departure from linearity. The slope of the plot ((f*, log rf)) increases sharply as rf increases. There are two possible reasons for the increasing tangent. One is the reduction of , which would lead to an increase in the slope (
) of f* vs. log rf. The other is the increase of curvature at the transition state region, i.e., the barrier top of the free energy landscape. Regardless of the precise reason, it is clear that the standard way of estimating using f* at large loading rates results in a very small value of . We estimate from (f*, log rf) plot to be 4 Å at rf 105 pN/s. The estimated value of is unphysical, because 4 Å is less than the average distance between neighboring P atoms. The minimum pulling speed used in our simulations is nearly five orders-of-magnitude greater than in experiments. The use of large loading rate results in small values of . If the simulations can be performed at small values of rf we expect the slope of f* vs. log rf to decrease, which would then give rise to physically reasonable values of at low rf. Our simulations suggest that the curvature in the plot of f* as a function of log rf is due to the dependence of on rf and not due to the presence of multiple transition states (22). As a result, extrapolation to low rf values can give erroneous results (Fig. 4 B).
View larger version (24K):
[in this window]
[in a new window]
|
FIGURE 4 Constant loading rate force unfolding. (A) The unfolding force distributions at different pulling speeds with hard (k = 70 pN/nm, top) and soft springs (k = 0.7 pN/nm, bottom). For the hard spring, the pulling speeds from right to left are v = 8.6 x 104, 8.6 x 103, and 8.6 x 102 µm/s. For the soft spring, the pulling speeds are v = 8.6 x 104, 1.7 x 104, 8.6 x 103, 8.6 x 102, and 8.6 x 101 µm/s from right to left force peaks. The peak in the distributions which are fit to a Gaussian is the most probable force f*. (B) The dependence of f* as a function of the loading rates, rf. The results from the hard spring and soft spring are combined using the loading rate as the relevant variable. The inset illustrates the potential difficulties in extrapolating from simulations at large rf to small values of rf.
|
|
Unfolding at constant force
To monitor the transition state movements, we performed a number of unfolding simulations by applying fS > 20 pN at T = 290 K. The unfolding rates are too slow at fS < 20 pN to be simulated. Nevertheless, the simulations give strong evidence for force-dependent movement of . For a number of values of fS in the range 20 pN < fS < 150 pN, we computed the distribution of first-passage times for unfolding. The first-passage time for the ith molecule is reached if R becomes R = 5 nm for the first time. From the distribution of first-passage times (for 50–100 molecules at each fS) we calculated the mean unfolding time. Just as for unfolding at constant rf, the dependence of log U on fS shows curvature (Fig. 5 A)—hence deviating from the often-used Bell model (31). By fitting U to the Bell formula (
), over a narrow range of fS we obtain 4 Å, which is too small to be physically meaningful at fS = 0.
Insights into the shift of as fS increases can be gleaned from the equilibrium force-dependent free energy F(R) as a function of R. The one-dimensional free energy profiles F(R) show significant movements in as fS changes (Fig. 5 C). As fS increases, decreases sharply, which implies that the unfolding TS is close to the folded state. At smaller values of fS, the TS moves away from the native state. At the midpoint of the transition, 
, which is approximately half-way to the native state. The result for 
(RU is the average value of R in the unfolded state) is in accord with experiments (5) done at forces that are not too far from the equilibrium unfolding force. Although is dependent on the RNA sequence, it is likely that, for simple hairpins, . The prediction that is dependent on fS is amenable to experimental test.
Force-quench refolding times depend (approximately) exponentially on fQ
One of the great advantages of force-quench refolding experiments (8) is that the ensemble of conformations with a predetermined value of R can be prepared by suitably adjusting the value of fS (8). Because force-quench refolding can be initiated from conformations with arbitrary values of R, regions of the energy landscape that are completely inaccessible in conventional experiments can be probed. To initiate refolding by force-quench, we generated extended conformations with R = 13.5 nm, using fS = 90 pN at T = 290 K. Subsequently, we reduced the force to fQ in the range 0.5 pN < fQ < 4 pN. For these values of fQ, the hairpin conformation is preferentially populated at equilibrium. The probability, PU(t), that the RNA hairpin remains unfolded at six fQ values, decays nonexponentially (Fig. 6). The mean refolding times, F(fQ), upon force-quench, which are computed using PU(t) (
), show that in the range 0.5 pN < fQ < 4 pN (Fig. 5 B),
 | (13) |
where 
is the distance from the unfolded basin of attraction to the refolding transition state, and 
is the refolding time in the absence of force. Linear regression gives 
and 
. The value of , which is obtained from kinetic simulations, is in accord with the location of obtained directly from the equilibrium free energy profiles F(R) (Fig. 5 C). In the 0.5 pN < fQ < 4 pN range, the distance from UBA to the TS is 1 nm, and is a constant. Thus, kinetic and thermodynamic data lead to a consistent picture of force-quench refolding.
View larger version (32K):
[in this window]
[in a new window]
|
FIGURE 6 Time-dependence of the probability that RNA is unfolded upon force-quench. In these simulations, T = 290 K, and the initial stretching force, fS = 90 pN, and fQ, the quench force, is varied (values are given in each panel). The simulation results are fit using Eq. 14, which is obtained using the kinetic scheme  . Here, I represents conformations within certain fractions of incorrect dihedral angles. The time constants (1,2), in µs, at each force are: (81.2, 101.3) at fQ = 0 pN; (159.5, 160.8) at fQ = 0.5 pN; (180.0, 174.8) at fQ = 1 pN; (237.6, 240.5) at fQ = 2 pN; (326.8, 335.6) at fQ = 3 pN; and (347.7, 329.7) at fQ = 4 pN.
|
|
If the simulations are done with fQ 0 pN, then we find that F(fQ = 0) 191 µs (Fig. 5 B), which differs from obtained using Eq. 13. When fQ is set to zero, the 3' end fluctuates, whereas when fQ 0, the 3' end remains fixed. The rate-limiting step in the hairpin formation is the trans
gauche transitions in the dihedral angles in the GAAA tetraloop region (see below). When one end is free to fluctuate, as is the case when fQ = 0, the transgauche occurs more rapidly than fQ 0. The difference between F(fQ = 0) and is, perhaps, related to the restriction in the conformational search among the compact structures, which occurs when one end of the chain is fixed.
Metastable intermediates lead to a lag phase in the refolding kinetics
The presence of long-lived conformations (see below), with several incorrect dihedral angles in the GAAA tetraloop, is also reflected in the refolding kinetics as monitored by PU(t). If there are parallel routes to the folded state, then PU(t) can be described using a sum of exponentials. The lag kinetics, which is more pronounced as fQ increases (see especially fQ = 4 pN in Fig. 6), can be rationalized using a kinetic scheme where S is the initial stretched state, I is the intermediate state, and F is the folded hairpin. Setting PU(t) PS(t) + PI(t) = 1 – PF(t), we obtain
 | (14) |
From Fig. 6, which shows the fits of the simulated PU(t) to Eq. 14, we obtain the parameters (1, 2) at each fQ (see caption to Fig. 6 for the values). If folding is initiated by temperature-quench, 1 << 2, so that 
. Explicit simulations show that thermal refolding occurs in a two-state manner (data not shown). In force-quench refolding, both 1 and 2 increase as fQ increases, and 1/2 O(1) at the higher values of fQ. There are considerable variations in the structures of the metastable intermediate depending on fQ. The variations in the conformations are due to the differences in the number of incorrect or nonnative dihedral angles. As a consequence, there are multiple steps in force-quench refolding in contrast to forced-unfolding, which occurs in an all-or-none manner.
Transgauche transitions in the GAAA tetraloop dihedral angles lead to long refolding times
It is of interest to compare the refolding times obtained from stretched ensemble (F(fQ)) and the refolding time (F(T)) from thermally denatured ensemble. In a previous article (14), we showed that F(fQ = 0) = 15F(T) (Fig. 5 B). The large difference in refolding times may be because the initial conditions from which folding commences are vastly different upon force and temperature quench (14,32). The fully stretched conformations, with R = 13.5 nm, are very unlikely to occur in an equilibrated ensemble even at elevated temperatures. The canonical distribution of thermally denatured conformations even at T = 1500 K (>> TF) shows that the probability of populating conformations with R = 13.5 nm (Fig. 7 A) is practically zero. Thus, folding from thermally denatured ensemble starts from relatively compact conformations. In contrast, the initial condition for force-quench refolding can begin (as in our simulations) from fully stretched conformations upon force-quench. Both R and the radius of gyration (Rg) undergo substantial changes en route to the NBA. Indeed, the refolding trajectories from extended conformations exhibit broad fluctuations in R(t) in the order of 25–75 Å) for long time periods, followed by a cooperative reduction in R at the final stage (Fig. 7 B).
View larger version (28K):
[in this window]
[in a new window]
|
FIGURE 7 (A) The equilibrium distribution of the end-to-end distance at extremely high temperature (T = 1500 K). Even at this elevated temperature, the fully stretched conformations of R = 13.5 nm (arrow) are not found in the ensemble of thermally denatured conformations. (B) Refolding is initiated by a force-quench from the initial value fS = 90 pN to fQ = 4 pN. The five time-traces show great variations in the relaxation to the hairpin conformation. However, in all trajectories, R decreases in at least three distinct stages that are explicitly labeled for the trajectory in green.
|
|
The long refolding times upon force-quench starting from fully stretched conformations may be generic for folding of globular proteins as well. Recent experiments on force-quench refolding of poly-ubiquitin (8) show that R(t) for proteins exhibit behavior similar to that shown in Fig. 7 B. The resulting fQ-dependent refolding times for poly-ubiquitin (0.1–10 s) are considerably larger compared to F ( 5 ms (33,34)) in the absence of force. Because the collapse of a single ubiquitin molecule in solution occurs in less than a millisecond, the origin of the long refolding times has drawn considerable attention (32,35). The microscopic model considered here for RNA can be used to shed light on the origin of the generic long refolding times upon force-quench.
From the phase diagram (see Fig. 2 in (14)), it is clear that the routes navigated by RNA hairpin upon thermal and force-quench have to be distinct. Although the distinct initial conditions do not affect the native state stability (as long as the final values of T and fQ are the same), they can profoundly alter the rates and pathways of folding. The major reason for the long force-quench refolding times in RNA hairpins is that, in the initial stretched state, there is a severe distortion (compared to its value in the native state) in one of the dihedral angles. The 19th dihedral angle (found in the GAAA loop region) along the sugar-phosphate backbone is in g+ conformation in the native state, whereas, in the initial stretched conformation, it is in the t state (Fig. 8 A). Thus, if all the molecules are in the fully stretched conformation, then the 19th dihedral angle in each of them has to, during the force-quench refolding process, undergo the t g+ transition in the GAAA tetraloop region in order to fold. The enthalpic barrier associated with the t g+ transition is 1 kBT (Fig. 8 B). However, this transition is coupled to the dynamics in the other degrees of freedom and such a cooperative event (see Fig. 7 B for examples of trajectories) makes the effective free energy barrier even higher. Although significant fluctuations are found in thermally denatured ensemble at T = 500 K, they are not large enough to produce nonnative dihedral angles in the GAAA tetraloop. The dihedral angles in thermally denatured conformations do not deviate significantly from their values in the native conformation (Fig. 8 C). In contrast, upon fully stretching, P5GA dihedral angles in the GAAA tetraloop adopt nonnative values (Fig. 8 D).
The timescale for the t g+ transitions can be inferred from the individual trajectories. Typically, there are large fluctuations due to g+ 
t g– transitions in the dihedral angles in the flexible loop region (i = 19–24). For the trajectory in Fig. 9 B, we observe the persistence of incorrect dihedral angle in the loop region for t 300 µs. At t > 300 µs, the nativelike dihedral angles form. Subsequently, the formation and propagation of interaction stabilizing the native RNA hairpin takes place. These dynamical transitions are clearly observed in Fig. 9, B and C. The observed mechanism is reminiscent of a nucleation process. We conclude that the formation of the flexible loop with all the dihedral angles achieving near-native values is the rate-limiting step in the refolding kinetics of RNA hairpins upon force-quench starting from the fully stretched state. It should be stressed that the rate-limiting step for thermal refolding of P5GA, or force-quench refolding starting from partially stretched conformations (see below), is different. The observation that the zipping of the hairpin takes place upon synchronous formation of all the nativelike dihedral angles suggests the presence of a high entropic barrier. The crossing of the entropic barrier results in slow refolding if P5GA is fully stretched.
Linker effects on RNA force-extension curves
In LOT experiments, the force-extension curves (FECs) are measured for the handle(H)-RNA-handle(H) construct (5,6). To unambiguously extract the FEC for RNA alone (from the measured FEC for H-RNA-H construct), the properties of the handle, namely, the contour length LH and the persistence length 
, cannot be chosen arbitrarily. In the experiments by Liphardt et al. (5), LH = 320 nm and 
for the RNA/DNA hybrid handle. We expect that both 
and will affect the FEC curves of the object of interest, namely, the RNA molecule. To discern the signature for the force-induced transition in the RNA hairpin alone from the FEC for H-RNA-H construct, 
has to be large. If we assume 
, then the experimental value of 50, which is large enough to extract the transitions in the RNA hairpin. If 1, then the entropic fluctuations in the handle can mask the signals in RNA (36). Similarly, the end-to-end distance fluctuation in the handle, R, should be smaller then the extension in RNA. Because R grows with LH (see below), it follows that, if very long LH is used, even with >> 1, the signal from RNA can be masked. The square of the fluctuations in the end-to-end distance R of the linker is given by (R)2 = 
R
/(ßfS), where R is the mean end-to-end distance. For WLC, which describes the linkers, we expect R (2LHlp)1/2 for small fS and R LH for large fS, provided LH is long. Thus,
 | (15) |
The mean fluctuation in the extension of the spring is, using equipartition theorem, given by
 | (16) |
For the signal from RNA to be easily discerned from the experimentally measured FEC, the expansion of end-to-end distance of the molecule at transition should be larger than R and x. Since R grows sublinearly with the linker length (Eq. 15), the attachment of large linker polymer can mask the transition signal. These arguments show that the characteristics of the linker can obscure the signals from RNA.
The FECs in the simulations can also be affected by nonequilibrium effects due to the linker dynamics. In the handle(H)-RNA-handle(H) construct considered here, the force exerted at one of the linkers depends on the angle between fS and the end of the linker. The initial event in force transmission along the contour of the H-RNA-H construct is the alignment of the molecule along the force direction. The characteristic time for force to reach RNA so that unfolding can occur is fc/rf, where fc is the critical force. Nonequilibrium effects due to linker dynamics become relevant if R, the timescale for alignment of H-RNA-H along the force direction, is R > fc/rf. This condition is not relevant in experiments conducted at small values of rf. However, it is important to consider nonequilibrium effects in simulations performed at high rf values. The characteristic time depends on LH and .
We validate these arguments by obtaining FEC for H-P5GA-H by varying LH and . Using the WLC for the linker we calculated FEC where the LH is varied from 10 to 50 nm. To observe rapid unfolding we have carried out our simulation at the pulling speed v = 0.86 x 102 µm/s and k = 0.7 pN/nm. Under these conditions nonequilibrium effects are relevant for the linker (23), which is not the case in experiments. The FECs show clearly a plateau in the range 20 pN < fS < 40 pN, which corresponds to the two-state hairpin opening (Fig. 10). For the experimentally relevant plot (Fig. 10 A) showing FEC for H-P5GA-H, the transition plateau is present at all values of LH. However, when LH reaches 50 nm, the signal from P5GA is masked. The FEC for P5GA alone (Fig. 10 B) shows modest increase in the unfolding force as LH increases. Similarly, we find the value of unfolding force also increases as the linker flexibility increases. These observations are due to nonequilibrium effects on the linker dynamics because of the relatively large values of rf used in the simulations.
View larger version (29K):
[in this window]
[in a new window]
|
FIGURE 10 The force extension curves (FECs) of RNA hairpin at constant pulling speed and varying linker lengths and flexibilities. Pulling speed is v = 0.86 x 102 µm/s; spring constant is k = 0.7 pN/nm. FECs from different linker lengths are plotted in black (10 nm), red (20 nm), green (30 nm), blue (40 nm), and orange (50 nm). (A) Shows the experimentally relevant plots, namely, FECs for H-RNA-H construct. The signature for the hairpin opening transition region is ambiguous at large LH values. (B) FECs only, for P5GA corresponding to different LH values. The FEC for the linker is subtracted from A. The gradual increase of rupture force is observed as L increases. The value of kA (Eq. 9) of linker polymer used in A and B is 80 kcal/(mol · Å). (C) Comparison between FECs with LH (= 40 nm) fixed, but at different kA-values. The red curve is the average of seven individual FECs, which are shown in orange (80 kcal/(mol · Å)). The black curve is the average of nine individual FECs in gray (20 kcal/(mol x Å)). A less-stiff linker leads to a slightly larger unfolding force. It should be stressed that, for both values of kA, the -ratio is large.
|
|
Our simulations show that, at high loading rates, the length of the handle is also important. This issue is not relevant in experiments in which loading rates are much smaller. However, they become important in interpreting simulation results. In our case, rf (= kv) is 6 x 104 times that used in experiments. At such high rf, nonequilibrium effects control the linker dynamics (23). Thus, to extract unfolding signatures from RNA alone, it is necessary to use a high value of and relatively short values of L. In other words: FEC, when rf is varied, may be a complicated function of 
and 
.
Force-quench refolding of P5GA with attached linkers
It is convenient to monitor force-quench refolding of RNA alone using simulations. A similar experiment can only be performed by attaching handles to RNA. In such an experiment, which has not yet been done (however, see Note in Proof added at the end of this article), the H-RNA-H would be stretched by a stretching force fS > fc so that RNA unfolds. By a feedback mechanism, the force is quenched to fQ 0. We have simulated this situation for the H-P5GA-H construct with LH = 15 nm, 
for each handle. We chose very stiff handles (
) so that the dynamics of RNA can be easily monitored. The end-to-end distance of the H-P5GA-H system as a function of t with fS = 90 pN and fQ = 2 pN shows a rapid decrease from Rsys = 44 nm to Rsys = 37 nm in <100 µs (Fig. 11). The use of large values of fS will not affect the results qualitatively. The value of Rsys fluctuates at 36 nm for a prolonged period and eventually Rsys attains its equilibrium value at 30 nm.
View larger version (34K):
[in this window]
[in a new window]
|
FIGURE 11 Refolding trajectory of RNA that is attached to the handles. Linkers with  and are attached to both sides of the RNA hairpin. The initial force (90 pN) stretched P5GA to 14 nm. The value of fQ = 2 pN. The top panel shows the dynamics of the end-to-end distance, Rsys, of H-P5GA-H. The middle panel corresponds to the end-to-end distance of P5GA. The end-to-end distances, RH, of the 3' and 5'-side handles, fluctuate around its equilibrium value of 14 nm after the initial rapid relaxation. The right panels show the decomposition of RH into the longitudinal and the transverse components. The value of  agrees well with the equilibrium value obtained by solving | |
TOP
ABSTRACT
INTRODUCTION
exp(
), in the range of 0.5 pN < fQ < 4 pN, we obtain 
and 
) for ith dihedral angle 