Originally published as Biophys J. BioFAST on March 11, 2005.
doi:10.1529/biophysj.104.045344
Biophysical Journal 88:3224-3242 (2005)
© 2005 The Biophysical Society
Thermodynamic and Kinetic Aspects of RNA Pulling Experiments
M. Manosas and
F. Ritort
Departament de Fisica Fonamental, Universitat de Barcelona, Barcelona, Spain
Correspondence: Address reprint requests to Felix Ritort, Tel.: 0034-0-0034-934-035-869; E-mail: ritort{at}ffn.ub.es.
 |
ABSTRACT
|
|---|
Recent single-molecule pulling experiments have shown how it is possible to manipulate RNA molecules using laser tweezers. In this article we investigate a minimal model for the experimental setup which includes an RNA molecule connected to two polymers (handles) and a bead trapped in the optical potential and attached to one of the handles. We start by considering the case of small single-domain RNA molecules, which unfold in a cooperative way. The model qualitatively reproduces the experimental results and allows us to investigate the influence of the bead and handles on the unfolding reaction. A main ingredient of the model is to consider the appropriate statistical ensemble and the corresponding thermodynamic potential describing thermal fluctuations in the system. We then investigate several questions relevant to extract thermodynamic information from experimental data. The kinetics of unfolding is also studied by introducing a dynamical model. Finally, we apply the model to the more general problem of a multidomain RNA molecule with Mg2+ tertiary contacts that unfolds in a sequential way.
|
INTRODUCTION
|
|---|
The RNA molecule plays a central role in molecular biology, showing an enzymatic function during the translation and splicing processes (Doudna and Cech, 2002
; Moore and Steitz, 2002). Experiments based on the manipulation of single biomolecules, such as laser tweezers with force microscopy, allow scientists to investigate their mechanical properties. These give information about the structure, stability, and the interactions involved in the formation of such structures (Bustamante et al., 2000; Smith et al., 1992, 1996; Cluzel et al., 1996; Essevaz-Roulet et al., 1997; Russell et al., 2002a; Zhuang et al., 2002). In these experiments mechanical force is applied to the ends of an RNA molecule. The molecule is then pulled (Liphardt et al., 2001; Onoa et al., 2003) until a value of the force is reached such that the molecule unfolds. If the pulling process is reversed then the molecule refolds again. In these experiments the force exerted upon the system is recorded as a function of the end-to-end distance giving the so-called force-extension curve (FEC). The nature of the unfolding-refolding reaction is stochastic and therefore the values of the force at which the molecule unfolds-refolds change from experiment to experiment. Sometimes (e.g., in presence of Mg2+ tertiary contacts), it is not possible to pull the molecule in quasistatic conditions because the relaxation time is too large for the experimental possibilities, which are largely limited due to the presence of strong drift effects in the machine. Therefore, during the pulling process, the molecule is driven to a non-equilibrium state, which is characterized by strong irreversibility effects. The study of this pulling process might be useful to understand many biological processes where biomolecules are unfolded under locally applied force; for example, when the mRNA goes through the ribosome during the translation process.
To manipulate an RNA molecule some synthesized polymers, typically several hundred nanometers long (called handles), have to be chemically linked to the extremes of the RNA molecule. Two polystyrene beads are then chemically attached to the end of these handles and one bead is used to measure the force by reading its position inside the optical trap. These additional elements (beads and handles) are an inseparable part of any pulling experiment and they have an influence on the unfolding process. To characterize the thermal behavior of the pulled global system (bead, handles plus RNA molecule) it is important to identify the proper control parameter. This is an essential step toward the modelization of the experiment and has several consequences. For instance, the force acting on the extremes of the RNA molecule cannot be externally controlled but fluctuates, and its mean value depends in a nonlinear way on the value of the control parameter. The control parameter determines the relevant thermodynamic potential that defines the equilibrium state of the global system as well as the magnitude of the fluctuations around that state. A proper inclusion of these parts is necessary to accurately interpret the experimental data. Another important aspect to consider in a theoretical treatment is the model for the RNA molecule. In this article we treat the RNA molecule as composed by different domains, each one showing cooperative unfolding. Each domain is then modeled as a two-state system: the unfolded state (UF) and the folded one (F), which are separated by a kinetic barrier. A main effort throughout this article is to present, in the clearest way, the appropriate theoretical framework to understand pulling experiments—leaving aside further additional complications, nevertheless important, such as the detailed response of the laser tweezers machine or the microscopic structure of the RNA molecule.
The goal of this article is twofold:
- We show how to build a minimal model aiming to reproduce the experimental setup, including all the aforementioned elements (bead, handles, and the RNA molecule), and quantitatively reproducing various experimental results.
- We show how to analyze experimental data extracted from both quasistatic and out-of-equilibrium pulling experiments to obtain thermodynamic and kinetic information about the unfolding reaction.
The article is divided into three main parts. In the first part, we describe the model for the experimental setup and introduce the ensemble that is relevant to model the pulling experiment. Then we describe the two-states model convenient to reproduce the cooperative unfolding of the RNA molecule and the models used for the bead and handles. In the second part of the article, we analyze the unfolding-refolding behavior of a cooperative two-states RNA molecule in a pulling experiment for both equilibrium and non-equilibrium regimes. For the equilibrium regime, we compute the partition function in the ensemble that is experimentally relevant, and derive an expression for the quasistatic work exerted upon the system as the molecule unfolds. This expression relates the work measured in a quasistatic pulling process to the difference of free energy between the F and UF states at zero force, 
G0. We analyze in detail the different thermodynamic contributions to the total work, the influence of the parameters describing bead and handles on the FEC, and obtain an expression for the force at the midpoint of the transition.
For the non-equilibrium behavior we investigate in detail the fraction of molecules that unfold (refold) more than once during the unfolding (refolding) path, which is a quantity amenable to experimental checks. We find that this fraction is related to the mean dissipated work exerted upon the system, which gives us a way to extract the reversible work in non-equilibrium processes just by measuring the total work. We also identify an interesting symmetry property relating these fractions for the forward and reverse processes. To endorse most of our theoretical results we also consider a simulation of a pulling experiment that allow us to obtain the characteristic FEC, either in a situation where the transition occurs in equilibrium or in a situation where it does not. In the third part of the article (Unfolding of Domains Stabilized by Mg2+ Tertiary Contacts), we address the unfolding behavior of complex RNA molecules with more than one folded-domain and in the presence of Mg2+-dependent barriers. In this case, back-refolding during the unfolding path is not observed at the experimental conditions, and the distribution of the breakage force is a first-order Markov process (Evans and Richie, 1997, 1999). We focus our attention in the specific case of RNA molecules where domains unfold in a sequential fashion according to a reproducible path. This unfolding mechanism is generally a consequence of the topological connectivity of the different parts of the molecule and of the blockade of the force induced by the most external tertiary contacts on the interior domains. Finally, we present the Conclusions. Three Appendices are devoted to describing some analytical calculations.
|
MODEL FOR THE EXPERIMENTAL SETUP
|
|---|
We consider a minimal model to reproduce the experimental setup of a pulling experiment carried out using laser tweezers (Liphardt et al., 2001; Smith et al., 2003). The model (Fig. 1) is composed by a small RNA molecule connected to two polymers called handles, which are used to attach the small RNA molecule to two beads at each end. One bead (B1) of radius Rbead is confined in the optical trap potential Vb(x) generated by the laser beams. The other one (B2) is held fixed to the tip of a micropipette by air suction. A micromanipulator controls the position of the micropipette relative to the optical trap. The stiffness of such a glass micropipette is much higher than the stiffness of either the optical trap or the handles; hence we neglect the micropipette-bead fluctuations. The molecule is pulled by moving the micropipette along the x direction. The configurational variables of this simplified system are taken as the projections of the end-to-end distances of each element along the force axis (Fig. 1): 
for the end-to-end distances of the handles, 
for the RNA end-to-end distance, and xb for the position of the bead B1 in the trap. We use the position xb of the bead B1 to read the force f acting on the system, as
 | (1) |
(Note that this is not the way the force is usually measured in dual beam optical tweezers where two photosensitive detectors located at opposite sides of the chamber are used to collect the total amount of deflected light, which is then converted into force after calibration of the machine; see Smith et al., 2003.) To a very good approximation, the optical trap is harmonic. Therefore,
 | (2) |
where kb is the stiffness of the optical trap. We define the subsystem S as that composed by the two handles and the small RNA molecule. The end-to-end distance for the subsystem S is then given by 
(Fig. 1). The total distance between the center of the trap and the tip of the micropipette is given by XT + Rbead = xb + x + Rbead. Pulling experiments give FECs, f(x), corresponding to the force (Eq. 1) as a function of the end-to-end distance of the subsystem S.
Ensembles
It is experimentally possible to consider two different ensembles depending on which variable is used as the externally imposed nonfluctuating parameter.
Mixed ensemble
The total distance between the center of the trap and the tip of the micropipette is held fixed, hence XT is the externally controlled parameter. In this ensemble there are fluctuations in x and f given by Gerland et al. (2003, 2001) as
 | (3) |
where 
...
stands for thermal average, kB is the Boltzmann constant, T is the temperature of the bath, kb is the stiffness of the optical trap (Eq. 2), and kx(XT) is the effective rigidity of subsystem S. The latter is determined by the serial compliance,
 | (4) |
where 
(i = 1, 2) and kr are the rigidities of the handles 1, 2, and the RNA, respectively. These rigidities are XT-dependent and so are the fluctuations (Eq. 3).
Force ensemble
In this case a piezo actuator controls the force (and therefore the position of the bead B1). In this ensemble XT and x are fluctuating variables, 
XT2 = x2 = kBT/kx(f), where kx(f) is the stiffness of the subsystem S when the force is held fixed, 
.
Most of the theoretical work for the force denaturation of RNA in pulling experiments considers the force ensemble. It might be possible to control the force using magnetic tweezers, which allows us to stretch and twist molecules by exerting forces in the range [1fN–10pN]. However, using optical tweezers it is experimentally very difficult to work in the force ensemble where either the force or the variable xb must be controlled and XT is a fluctuating variable. To compensate the fluctuations in the force, the distance XT should be corrected by a feedback mechanism that is difficult to implement. Therefore the most natural ensemble is that where XT is constant. Indeed, this is the most relevant ensemble for the experiments and therefore we will work in the mixed ensemble throughout this article.
|
MODELING THE DIFFERENT PARTS OF THE SETUP
|
|---|
Two-states model for a single RNA domain under mechanical load
The unfolding of some biomolecules under the effect of a mechanical force is a highly cooperative process that can be qualitatively described by a two-states model. The two-states model has a long tradition in physics, and has been applied previously by several authors to explain the unfolding behavior of single domains of proteins and RNA hairpins (Liphardt et al., 2001; Ritort et al., 2002; Fernandez et al., 2001; Muñoz et al., 1997; Bokinsky et al., 2003; Zhuang et al., 2000a). Recently, it has been shown how such a simple phenomenological description, with Kramer transition-rates, does not fully reproduce the kinetics observed in pulling experiments of the protein Titin, and more realistic descriptions have been proposed (Hummer and Szabo, 2003).
Let us consider an individual RNA molecule in thermal equilibrium with water solvent (at physiological conditions) at constant temperature, pressure, and zero force. In the simplest description, both states (hereafter denoted by UF, unfolded; and F, folded) are characterized by their Gibbs free energy 
and 
respectively, and the RNA molecule occupies each state with a probability given by the Boltzmann distribution. In a more realistic description the molecule can also occupy intermediate configurations, depending on the number n of the first-opened, or denaturated, basepairs (Cocco et al., 2003, Marinari et al., 2002).
When an externally controlled force f is applied to the ends of the RNA molecule, the adequate thermodynamic potential to consider is the Legendre transform of the Gibbs free energy G'(n) = G0(n) – fxr(n) (Tinoco and Bustamante, 2002), where G0(n) and xr(n) stand for the free energy and the projection of the end-to-end distance in the axis force of a hairpin with the first n basepairs opened, respectively. The free-energy landscape G' is then tilted along the reaction coordinate xr, which explicitly depends on the number of opened basepairs n. Since we work in the ensemble where neither f nor xr are control parameters, the nonfluctuating parameter XT determines the adequate thermodynamic potential 
The free-energy 
of the system shown in Fig. 1 is a potential of mean force that characterizes the equilibrium state of the whole system, including the handles, the bead, and the RNA molecule at a fixed value of XT. The potential free-energy landscape associated to is shown in Fig. 2, where we represent it as a function of the end-to-end distance of the subsystem S, 
. The shape of the potential shows two pronounced minima corresponding to the F and UF states. The discrete variable 
stands for the state of the domain: the value = 0 denotes the F state and = 1 the UF state. The relative thermodynamic stability of these states depends on the difference of free energy between them, G(XT). Moreover, we will consider the existence of a transition state along the reaction path from the F to the UF state and vice versa. This transition state is the RNA configuration with highest free energy connecting the F and the UF states along the reaction path. It may correspond to an RNA configuration where the first n = n* basepairs are opened. (We stress that the shape of the free-energy landscape depends on XT as well as the location of the barrier corresponding to the transition state. However, for the sake of simplicity, we will assume n* independent of XT.) In the simplest scenario the transition state can be assumed to have a very short lifetime. Therefore it can be represented by an activation barrier whose main effect is to hinder transitions between the F and UF states. This is the model we will adopt throughout the article. The F and UF states are separated by a barrier of height B(XT) measured relative to the F state. The barrier is located at a distance x1(XT) from the F state and x2(XT) from the UF state. The distance between the two states is xm(XT) = x1(XT) + x2(XT). Since the rigidity of the RNA molecule in the F state is very large, we can assume this state to be characterized by a single configuration corresponding to the value xr = 0 of the reaction coordinate. The RNA in the UF state has a finite rigidity, hence it is represented by a set of configurations within a continuous range of values of xr (Fig. 2).
Modeling the bead, handles, and the ssRNA
In this section we specify the models for the different elements of the system: the bead trapped in the optical-tweezers potential, the two handles, and the single-stranded RNA (ssRNA) molecule.
Model for the optical tweezers: a bead matched to a spring
We model the optical potential as an harmonic potential of stiffness kb (Eq. 2); hence the bead in the optical trap can be considered as a bead matched to a spring. We consider that the bead follows a Langevin dynamics of an overdamped particle (i.e., without inertial term),
 | (5) |
where 
(with = 6
Rbead, , and Rbead being the viscosity of the water and the radius of the bead, respectively) is the friction coefficient and FR is the resultant force applied to the bead. (In Eq. 5, we are neglecting the drag force felt by the bead, equal to –v, as the chamber is moved and the water dragged relative to the lab frame at a certain pulling speed 
For the range of pulling speeds used in the experiments this contribution is negligible, of the order of 0.1 pN.) The stochastic term 
(t) is a white noise with mean value (t) = 0 and variance (t)(t') = 2 kBT(t – t'). The force FR has two contributions, FR = fx – f: the force generated by the optical trap potential, f, given by Eq. 2, and the tension exerted by the subsystem S, fx. Using the equilibrium condition FR = 0 or f = fx, and doing an expansion around the equilibrium position of the bead, xeq, we get
 | (6) |
where kR is the effective spring constant applied to the bead, kR = kx + kb, with kx given by Eq. 4. The relaxation time 
b of the system (i.e., the typical time during which the position of the bead de-correlates) is given by b = /kR.
Polymer model for the handles and the ssRNA
To model the handles and the single-stranded RNA (ssRNA) we use the worm-like-chain (WLC) model. The thermodynamic properties of this model cannot be exactly computed, yet there are useful extrapolation formulas. A simple expression has been proposed (Bustamante et al., 1994) for the force as a function of mean end-to-end distance of the polymer x,
 | (7) |
where Lo and P are the contour and persistence lengths of the polymer, respectively. Equation 7 converges asymptotically to the exact solution as x approaches either zero or Lo and is accurate at least up to 90% in between. Bouchiat et al. (1999) have given an expression with an accuracy of 99% by adding a polynomial of seventh order to Eq. 7. The WLC model works well only at low forces, in the so-called entropic regime, where the molecule behaves as an entropic spring. At high forces there is an enthalpic correction due to the fact that the phosphodiester bonds along the backbone are stretched and the contour length Lo increases. To incorporate this effect it is common to replace x/Lo by x/Lo – f/Ey in Eq. 7, where Ey is the Young modulus of the polymer.
|
THERMODYNAMIC ANALYSIS
|
|---|
In this section we use the tools of statistical mechanics to analyze the thermodynamics of the system represented in Fig. 1. Most of the analytical treatment is described in the Appendix A. In what follows, we review the main results of these calculations.
Definitions
In equilibrium the observables x
and their conjugated forces f with = h1, h2, r, b (referring to the different elements: handle 1 and handle 2, RNA and bead, respectively) are fluctuating quantities. However, the thermodynamic free energy is only a function of the mean values of these observables that we denote by x, f. A representation of f versus x gives what we call the thermodynamic force extension curve (TFEC) for the element in the mixed ensemble. If refers to the whole subsystem S, then the TFEC corresponds to the usual force-extension curve (FEC) recorded in RNA pulling experiments, assuming that the pulling process is carried out reversibly. Throughout the thermodynamic analysis, and to simplify the notation, we will use indistinctly f(x) or f(x) to denote the TFEC. We can also define the restricted average O (XT) as the mean value of the observable 
when the RNA molecule is in the state (i.e., folded or unfolded) for a fixed total end-to-end distance XT. From now on, all the dependencies of the observables on the variable XT will not be explicitly written, hence (XT) 
. In Appendix A, we derive an expression for the partition function Z(XT) corresponding to the system schematically represented in Fig. 1. Applying the saddle point technique, and separating the contributions from the F ( = 0) and the UF ( = 1) state we get
 | (8) |
where
 | (9) |
 | (10) |
with 
Here Vb represents the optical trap potential and G0 is the free-energy difference between the F and the UF states at zero force. The function W(x) corresponding to the reversible work performed by adiabatically stretching the element from x = 0 to x = x, when the molecule is in the state , reads
 | (11) |
where f(y) is the TFEC for the element . The thermodynamic value of any observable can be expressed as
 | (12) |
where p0 and p1 are the probabilities for the RNA molecule to be in the F and UF states, respectively,
 | (13) |
At the transition midpoint both states are equally probable,
 | (14) |
where these functions have been defined in Eqs. 9, 10, and 13. Hence, the transition midpoint in the mixed-ensemble is defined by the value of the control parameter 
that verifies Eq. 14.
Computation of the transition force Fc, the TFEC, and the different contributions to the reversible work
The force at the transition, Fc, is computed as the mean value of the force at given by Eq. 14. To reproduce the experimental results obtained for the P5ab RNA molecule in 10 mM Mg2+ (Liphardt et al., 2001) we use the parameters given by Tables 1 and 2 getting Fc = 15.2 pN. This value is close to the one reported from the experiments 
(Liphardt et al., 2001). We also verify that the value of the computed force at the transition, Fc, is quite stable with respect to changes in the parameters of the problem used to model the handles and the bead, such as the persistence and contour lengths of the handles, the spring constant, and the bead radius. However, because the value of Fc is highly influenced by the characteristics of the RNA molecule, we conclude that the dependence of the value of Fc with the system is basically through the quantities G0, Lr, and Pr.
Another interesting magnitude to measure is the reversible work 
done upon the system when pulling from an initial value XT = XT0 to a final value of XT. This work is given by
 | (15) |
where we used Eq. 8. The total reversible work in Eq. 15 defines the change in the free energy of the system. In Fig. 3 A we show the total work and its different contributions, 
as a function of XT as derived from the numerical computation of Z(XT), where the reversible work exerted on each element (handles 1 and 2, bead, and RNA molecule) is defined as
 | (16) |
where
 | (17) |
 | (18) |
 | (19) |
View larger version (16K):
[in this window]
[in a new window]
|
FIGURE 3 (A) Different contributions to the reversible work obtained from the partition function analysis:  as a function of XT. Note that the smallest contribution to the total work comes from the RNA molecule. (B) The continuous line corresponds to the results obtained from the numerical computation of the TFEC. It is also shown that the TFEC is obtained by averaging over 1000 different trajectories, as explained in Simulation of a Pulling Experiment. The pulling is carried out at an approximate loading rate of 0.5 pN/s, slow enough to generate a quasistatic process. One can observe that both curves agree.
|
|
The functions Vb, Wh, and Wr correspond to the change in the potential energy of the bead in the optical trap and the work exerted upon the handles and the RNA molecule by moving the total end-to-end distance from the initial to the final value of XT, respectively. Finally in Fig. 3 B we represent the TFEC for the subsystem S, f versus x. This is obtained by numerical computation of the partition function using the relation
 | (20) |
To calculate f we use the expression (Eq. A1) for the partition function of the system Z(XT). Integrating this expression over the bead position xb gives
 | (21) |
where Z(x) is the partition function of the element , with = h1, h2, r, and b. The partition function for the bead (Eq. A2) satisfies
 | (22) |
where we used Eq. 2. Therefore the absolute value of the mean force, Eq. 20, can be computed as
 | (23) |
where 
The generalized force, Eq. 20, is the average force measured in the optical trap.
Reversible work across the transition
The quasistatic work 
exerted upon the subsystem S across the transition is the area under the TFEC (Fig. 4), f(x), from the folded branch x = xc0 ( = 0) to the unfolded branch x = xc1 ( = 1) (see Appendix A), where the super-index c indicates that the system is at the transition midpoint, 
(Eq. 14),
 | (24) |
View larger version (15K):
[in this window]
[in a new window]
|
FIGURE 4 The shadow area under the TFEC along the transition corresponds to the quasistatic work (schematic representation).
|
|
At the transition midpoint, both states are equally populated and Eq. 14 holds. Therefore identifying Eqs. 9 and 10, we can write Eq. 24 as
 | (25) |
where the functions with a super-index c are evaluated at the mean value of their variables at the critical extension 
The value Wr, given by Eq. 11, is the loss of entropy of the RNA molecule along the transition due to the stretching. The value Wh is the free-energy change of the handles between the folded and unfolded branches, and is given by
 | (26) |
Equation 25 tells us that the quasistatic work coincides with the change of free energy of the different elements that form the subsystem S across the transition. The value is experimentally measurable as the area under the rip observed in the TFEC corresponding to the F–UF transition (Fig. 4). Therefore Eq. 25 provides a way to estimate the unfolding free energy of the molecule G0 from the TFEC, which is a quantity biologically relevant as it determines the direction of biochemical reactions. This free energy G0 is equal to the Gibbs free energy measured by thermal denaturation in bulk experiments extrapolated to the working temperature.
In Fig. 5 we show two TFECs obtained from the partition function analysis corresponding to two systems with different kb but with the same handles and RNA molecule with parameters given in Tables 1 and 2, respectively. We use Eq. 25 to extract the value of G0 by computing as the area under the rip in the TFEC (Fig. 4). As expected for an harmonic trap (Eq. 2), the TFEC in Fig. 5 shows an slope at the transition (rip) proportional to –kb. To obtain the different contributions to Eq. 25 we first use the WLC model (Bouchiat et al., 1999) to estimate 
and 
given by Eqs. 11 and 26. Finally, we compute the area under the TFEC across the transition (rip) to obtain and use Eq. 25 to extract G0. The results are given in Table 3. Note that the contribution is negative because when the RNA molecule opens the force relaxes and the handles contract, hence the free energy of the handles across the transition decreases. Neglecting the contribution that comes from the handles across the transition is a typical approximation often applied to experimental results. However, this is not always accurate as this contribution can be large. In the previous example, even in the case of small kb, we would lose 8 kBT in the balance equation (Eq. 25). In Fig. 6 we show, for a small value of kb (kb = 0.1 pN/nm), how the different contributions to Eq. 25 change when considering systems with different values for the ratio Lh/Ph. The stretching contribution to the UF state of the RNA, 
does not change when modifying the magnitude Lh/Ph, because the forces at which the transition occurs are quite stable under changes of Lh/Ph. However, the magnitude of the contribution tends to notably increase as Lh/Ph becomes larger.
View larger version (15K):
[in this window]
[in a new window]
|
FIGURE 5 TFEC corresponding to two systems with handles and RNA characterized by the parameters given in Tables 1 and 2 and with an optical trap stiffness kb = 0.1 pN/nm and kb = 1 pN/nm, respectively. Note that the slope of the TFEC at the transition (rip) is proportional to –kb.
|
|
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 6 The different contributions to the free-energy change across the transition presented as a function of the ratio Lh/Ph. Note that the value of G0 is independent of that ratio (see also Table 3).
|
|
|
SIMULATION OF A PULLING EXPERIMENT
|
|---|
To simulate a pulling experiment it is important to distinguish the different timescales involved in the problem. Typically the bead has a much bigger size than the other components of the system (handles and ssRNA), therefore the bead is the element with largest dissipation and slowest relaxation compared to the elastic and bending modes of the handles and the ssRNA: b >> handles, ssRNA. The characteristic time b at which the bead relaxes to its equilibrium position can be computed as the ratio between the friction coefficient of the bead and the effective spring constant applied to the bead kR = kx + kb (see Eq. 6), with kx given by Eq. 4, i.e., b = /kR. For the typical experimental values for the trap stiffness and the radius of the beads, kb 
0.05–0.15 pN/nm and Rbead 1–3 µm, the time b lies in the range [10–3 s–10–6 s]; its value depends upon the value of the control parameter XT. The characteristic time F–UF at which the RNA hairpin folds and unfolds depends on the sequence of bases and also on the presence of tertiary contacts that slow down the kinetics of the unfolding reaction. Typical values are of the order of seconds-to-milliseconds. Hence the dynamics of the system shows the following separation of timescales: F–UF >> b >> handles, ssRNA. Therefore we can consider an instantaneous relaxation for the handles and the bead to solve the dynamical equations that describe the folding-unfolding kinetics of the RNA molecule. This hypothesis is valid, as long as the data is collected at frequencies smaller than the relaxational frequency of the bead, which is the element with the largest relaxation time. The dynamics for the RNA molecule is governed by a master equation for the probability p (Eq. 13),
 | (27) |
The functions k
and k
are the unfolding and folding rates corresponding to the activated process schematically represented in Fig. 2,
 | (28) |
where k0 is an attempt frequency. These rates satisfy the detailed balance condition,
 | (29) |
The expressions of G(XT) and B(XT) are derived in Appendix B using the partition function analysis.
To simulate a pulling experiment we use an adiabatic approximation by taking advantage of the great separation of timescales between the folding-unfolding kinetics and the relaxational dynamics of the different elements of the system. At each value of the extension XT and for a given state of the RNA molecule ( = 0, 1) we determine the values of the mean extension and force for the bead, handles, and ssRNA using the equilibrium equations. At the same time we numerically solve the dynamics for the RNA molecule (Eq. 27). In what follows, we describe the steps of the algorithm:
Step 1. We increase XT by vt, where v is the pulling speed, i.e., the velocity at which the micropipette is pulled, 
and t is the iteration time, hence 
is the frequency at which data is collected. Note that the relation between the pulling speed v and the loading rate r, i.e., the velocity at which the force increases, can be found using the relation between the force and displacement increments, f = keff(f)XT, as
 | (30) |
where keff is the effective stiffness of the system, computed as
 | (31) |
and where kx has been defined in Eq. 4 and kb is the stiffness of the optical trap. The F–UF transition for a small single RNA domain typically occurs at forces in the range 8–20 pN. At these forces the system verifies that kb is much smaller than the stiffness of the handles and the RNA molecule, 
and kr; and therefore we can safely take v = r/kb.
Step 2. We compute the new f and x iteratively, using the saddle point equations for the partition function. To these mean values we add Gaussian fluctuations of zero mean and variance given by Eq. 3. We then obtain the FEC, f(x), which should qualitatively reproduce the experimental one.
Step 3. The RNA molecule is then unfolded or folded, with a probability k(XT)t or k(XT)t, respectively, where t is the iteration time. For the rates, we use Eq. 28 with the functions G(XT) and B(XT) given by Eqs. B2 and B3, respectively.
Force-extension curve results (FEC)
In Fig. 7, A and B, we show the resulting FEC of our simulations for the values used in the experiment of Liphardt et al. (2001) shown in Tables 1, 2, and 4, corresponding to a P5ab RNA molecule and for a loading rate of r = 1 pN/s and of r = 50 pN/s, respectively. In these simulations we implement the dynamical algorithm previously described for the forward and reverse processes where XT increases and decreases in time, respectively. As shown in Fig. 7 A, at a loading rate of 1 pN/s, different transition jumps are observed along both the forward and reverse processes, because the pulling speed (v) is low enough. Comparing these simulation results with the experimental FEC (Liphardt et al., 2001) shown in Fig. 8, we find a qualitative agreement, and the shape of the curve around the transition region is qualitatively reproduced. However, we find some discrepancies:
- The simulated curve is shifted in the x direction in comparison with the experimental one. This is because, experimentally, the quantity measured is the relative change in x rather than its absolute value. Indeed, there may be some uncertainty (typically of the order of 100 nm) in the diameter of the bead used in the experiments. The value of the diameter is required to determine the distance x from the experimentally measured value of the distance between the centers of the two beads (equal to x + 2Rbead). Therefore, in Fig. 8, the extension represented in the x axis corresponds to changes in the value of x with respect to an initial extension of
100 nm.
- As the force increases, the experimental curve separates from the theoretical WLC prediction and therefore from the simulated results. The agreement can be improved by considering larger values for the Young modulus of the handles and of the ssRNA. Furthermore, by extending the RNA molecule model to include intermediate configurations, which depend on the number of opened basepairs n, we realize that the cooperative transition might not be between the F (n = 0) and UF (n = N) states, but between partially folded and partially unfolded states. For instance, for the P5ab RNA molecule, the cooperative folding-unfolding transition is between the state n = 3 and the state n = N (Cocco et al., 2003). This means that typically the first three basepairs open before the transition occurs, increasing the extension of the handles.
View larger version (17K):
[in this window]
[in a new window]
|
FIGURE 7 Results for the FEC obtained from the simulation of a pulling experiment with iteration time t = 10–2 s. (A) Pulling rate r = 1 pN/s. (B) Pulling rate r = 50 pN/s. At this pulling rate the molecule is driven out of equilibrium, and hysteresis is observed around the transition.
|
|
Fig. 7 B shows the FEC corresponding to a pulling process carried out at a loading rate of r = 50 pN/s. At this pulling speed, the process is not in equilibrium, and hysteresis effects are observed around the transition region.
Fraction of trajectories that have at least one refolding
We consider a system with a control parameter (generally denoted by y) that is pulled by changing y at certain speed 
. The forward (reverse) pulling process starts at a initial value of the control parameter yi (yf) where the RNA is in the F (UF) state and finishes at a final value of the control parameter yf (yi) where the RNA is in the UF (F) state. We then define NF and NR as the fractions of forward and reverse trajectories that have at least one refolding, respectively (Fig. 9). These fractions are given by
 | (32) |
 | (33) |
where the first integral in the right-hand side of both equations accounts for the probability of unfolding (folding) before a certain value of the control parameter y is reached and the second integral accounts for the probability of refolding once the RNA molecule has been unfolded (folded). The function 
F(R)(z, z') is the probability that the RNA molecule remains at the state until y = z' starting at y = z in the forward (reverse) process. The term is the solution of the master equation:
 | (34) |
 | (35) |
with initial condition 
In Appendix C we prove that the fraction NF is equal to NR if the perturbation protocol for the control parameter is symmetric, i.e., if the velocities along the forward and reverse process verify vF(y) = – vR(y). In our analysis the control parameter y corresponds to the total distance XT and the folding-unfolding rates are given by Eq. 28. The detailed analytical expressions for the rates have been given in Eqs. B5 and B6. Analytical computations with such rates appear quite cumbersome and it is preferable to simplify them. For analytical purposes, we will consider effective rates where the functions B1 and G1 given by Eq. B7 and x1 and x2 (the distances from the F and UF states to the transition state along the x axis, see Fig. 2) are effective parameters independent of XT. We call these 
obtaining
 | (36) |
where the force f ( = 0, 1) is the force acting upon the system at a given value of XT when the RNA is in the state . (The approximation expression in Eq. 36, where force does not fluctuate near the transition, is well justified. In fact, when the RNA is in a given state, i.e., folded or unfolded, the magnitude-of-force fluctuations is negligible, with the RMS in the range 0.03–0.1 pN, so one can consider the instantaneous force equal to the mean force. Hence the fluctuations in force near the transition arise solely from the force jump between the F and UF states.) The forces f0 and f1 in Eq. 36 correspond to the two branches (Eq. A16): f1 = kb(XT – x1) and f0 = kb(XT – x0), where we used Eq. 23. Therefore, the relation between f0 and f1 reads as
 | (37) |
where 
is the distance between the F and UF states along the x axis, 
. Using Eq. 37, it is straightforward to see that the effective rates (Eq. 36) satisfy the detailed balance condition (Eq. 29). We can now compute the fractions (Eqs. 32 and 33) as a function of the loading rate r. In Fig. 10, we show the results obtained for the fractions NF and NR from the numerical computation of Eqs. 32 and 33, using the effective rates (Eq. 36) with the definitions in Eqs. B5–B7. We also show the results obtained from the simulations for the fractions NF and NR as a function of the loading rate r, and they agree fairly well.
View larger version (7K):
[in this window]
[in a new window]
|
FIGURE 9 Different trajectories that have at least one refolding. The ratio between the sum of these trajectories and the total number of trajectories gives the fraction NF for the forward process and NR for the reverse process.
|
|
View larger version (26K):
[in this window]
[in a new window]
|
FIGURE 10 The fractions NF and NR as a function of r. Simulation results correspond to 5000 realizations of a pulling experiment. We also show the numerical integration of Eq. 32 (equal to Eq. 33; see Appendix C) using the rates (Eq. 36) with the following parameters: 
|
|
From these simulations we can also compute the mean work exerted upon the system as a function of r,
 | (38) |
where fi is the force acting on the system (Eq. 2), XT is the uniform increase in the total end-to-end distance at each iteration, and n is the total number of iterations. The average is over different realizations of the simulation of the pulling process. The total mean work is the sum of the reversible work (i.e., the work measured in a quasistatic process, r 0), and the mean dissipated work, 
We then consider the fraction NF(r) for three different RNA molecules characterized by different parameters, i.e., G0, Lr, N (total number of basepairs), n*, and B0 ln k0; the results are shown in Fig. 11 A. Plotting these fractions NF as a function of the mean dissipated work Wdis exerted upon the system, we find that the three curves corresponding to the three RNA molecules show the same kind of dependence (Fig. 11 B). This dependence is not surprising as the average dissipated work has been already shown (Ritort et al., 2002) to be a useful quantity to characterize the non-equilibrium regime. In particular, in the linear response regime, the average dissipated work depends linearly on the loading rate r, the proportionality constant being a function of the relaxation time of the molecule, the unfolding free energy, and the transition force (Ritort et al., 2002). The collapse of all curves in Fig. 11 is, however, not restricted to the linear response regime. Indeed, we have verified that in the regime 2 kBT < Wdis < 5 kBT, where deviations from the linear response regime are observable (Fig. 12), there is still a good collapse in Fig. 11 B of the curves corresponding to the three molecules. Note that by measuring the fraction NF we can obtain information about the value Wdis, and from the knowledge of the total work we can extract the reversible work exerted upon the system. This provides an alternative way to derive equilibrium information from non-equilibrium experiments (Liphardt et al., 2002; Ritort, 2003).
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 11 (A) The fraction NF as a function of r for three different RNA molecules characterized by Molecule 1, G0 = 59 kBT; Lr = 28.9 nm; N = 24; n* = 12; and B0 ln k0 = 29 kBT. Molecule 2, | |
TOP
ABSTRACT
INTRODUCTION
and 
which are the projections of the end-to-end distance of each element along the reaction coordinate axis (i.e., the axis along which the force is applied). The potential Vb(xb) is well described by an harmonic potential for a one-dimensional spring with rest position at xb = 0.
for the whole system at temperature below the melting temperature, for 
(where 