Spin-oscillator model for the unzipping of biomolecules by mechanical force

PHYSICAL REVIEW E 86, 021919 (2012) Spin-oscillator model for the unzipping of biomolecules by mechanical force A. Prados,1A. Carpio,2and L. L. Bonilla3 1F´ ısica Te´ orica, Universidad de Sevilla Apartado de Correos 1065, E-41080 Sevilla, Spain 2Departamento de Matem´ atica Aplicada, Universidad Complutense de Madrid, 28040 Madrid, Spain 3G. Mill´ an Institute for Fluid Dynamics, Nanoscience and Industrial Mathematics, Universidad Carlos III de Madrid, 28911 Legan´ es, Spain (Received 11 May 2012; revised manuscript received 8 July 2012; published 21 August 2012) A spin-oscillator system models unzipping of biomolecules (such as DNA, RNA, or proteins) subject to an external force. The system comprises a macroscopic degree of freedom, represented by a one-dimensional oscillator, and internal degrees of freedom, represented by Glauber spins with nearest-neighbor interaction and a coupling constant proportional to the oscillator position. At a critical value Fcof an applied external force F, the oscillator rest position (order parameter) changes abruptly and the system undergoes a first-order phase transition. When the external force is cycled at different rates, the extension given by the oscillator position exhibits a hysteresis cycle at high loading rates, whereas it moves reversibly over the equilibrium force-extension curve at very low loading rates. Under constant force, the logarithm of the residence time at the stable and metastable oscillator rest position is proportional to F−Fcas in an Arrhenius law. DOI: 10.1103/PhysRevE.86.021919 PACS number(s): 87.15.Cc, 05.40.−a, 05.50.+q, 64.60.De I. INTRODUCTION Many physical situations can be modeled by a mechanical system coupled to a thermal bath or to spin systems. Examples abound: The collective Jahn-Teller effect has been analyzed by spin-phonon systems [1–3], mass spectrometry through a nanoelectromechanical oscillator whose resonant frequency decreases as single molecules are added thereto [4], decoherence of a spin representing a two-level system due to coupling to a boson bath (the spin-boson system) [5], a classical oscillator coupled to a spin causing wave function collapse thereof [6], a1/2-spin representing a nonlinear Josephson phase quantum bit coupled to an oscillator (superconducting resonator) and to a classical signal [7,8], rippling in clamped graphene sheets investigated by means of a spin-string system [9], etc. Very recently we have introduced a simple model in which a single oscillator is coupled to a chain of Ising spins undergoing Glauber dynamics in contact with a thermal bath [10–12]. In our model the spins in the chain are coupled only to their nearest neighbors, but their coupling constant is proportional to the oscillator position, which makes their interaction effectively long range. In equilibrium, elimination of the oscillator coordinates gives rise to an effective spin interaction equivalent to a one-dimensional Ising model with mean-field coupling [13]. There is a second-order phase transition at a finite temperature, with the oscillator rest position as its order parameter. Above the critical temperature, the oscillator rest position is zero, thereby coinciding with that of the uncoupled oscillator. Below the critical temperature, two symmetric nonzero rest positions issue forth symmetrically from zero as in the case of a pitchfork bifurcation. In the limit of fast relaxation of the spins compared to the natural period of the oscillator, the oscillator position satisfies an effective equation having both nonlinear force and nonlinear friction terms [10,12]. Interestingly, this nonlinear friction arises from the coupling of the macroscopic elastic mode with the internal degrees of freedom (modeled in our system by the spins). A related mechanism has been proposed to explain the “internal friction” observed in experiments with proteins or polymers in solution [14]. In recent years, technological development has allowed to manipulate or visualize individual molecules and to measure microscopic forces with high-precision instruments. These single-molecule experiments (SMEs) provide key information about the thermodynamic and kinetic properties of biomolecules, offering a complementary but different perspective to understand molecular processes. An extensive review of these techniques can be found in Ref. [15]. Using SMEs, distributions describing certain molecular properties can be measured, thereby allowing one to characterize the kinetics of biomolecular reactions and the observation of possible intermediate states. A typical outcome of SMEs are the force-extension curves for DNA, RNA, and other biomolecules like proteins. In a seminal paper, Liphardt et al. [16] pulled a RNA hairpin by an increasing applied force until it unfolds at a critical value of the force, Fc≃14.5 pN. Afterwards, the molecule is pushed back, by decreasing the force, until it refolds. At low pulling (pushing) rates, the stretching and relaxing force-extension curves are superimposed, and the molecule unfolds at the critical value of the force Fc. Taking a closer look at this transition, hopping between the two possible extension values is observed. This suggests that the system is bistable: There are two possible states of the molecule with stochastic transitions between them. This physical picture is confirmed by experiments carried out at constant load, in a narrow region around the critical force [see Fig. 2(c) of Ref. [16]]. The residence times in the zipped and unzipped states have an Arrhenius-like dependence. When cycles of pulling and pushing the molecule are carried out at high loading rates, the extension of the molecule occurs at a higher force F+>F c, whereas the hairpin folds at a lower force F−<F c. Thus a hysteresis cycle arises, and some authors have been claimed this to be a signature of irreversible nonequilibrium behavior [16–19]. More recently many works have tried to understand these unzipping experiments from a physical point of view [17–26]. The elastic response of DNA [27–36] and other biomolecules like proteins [37–41] exhibits similar behavior. Despite some differences in minor details, in all cases the 021919-1 1539-3755/2012/86(2)/021919(9) ©2012 American Physical Society A. PRADOS, A. CARPIO, AND L. L. BONILLA PHYSICAL REVIEW E 86, 021919 (2012) “basic” behavior is like that found in Ref. [16]: There is a critical value of the applied force Fc, at which the molecule unzips (or unfolds) and its length increases discontinuously. Some peculiarities arise when going into fine detail though. On the one hand, in experiments with DNA or RNA hairpins the molecule is attached to two beads and the distance between them is controlled. Thus the total length of the molecule is kept fixed in these experiments, and it is the control parameter. The stretching transition is accompanied by a drop in the measured force [16–19]. In very recent experiments, the applied force is the control parameter, and the unzipping transition occurs at a constant value of the force [25]. For proteins, on the other hand, a typical sawtooth pattern with regularly spaced force peaks is often observed at the transition. This behavior is usually interpreted as the successive unfolding of each of the protein domains in a single protein molecule [37–39]. Performing pulling and pushing back cycles, both a “simple” reversible behavior and hysteresis cycles similar to those of Ref. [16]have been observed [40,41]. Finally, in the case of DNA, there are different transitions that depend on the experiment details [30]. There is a so-called overstretching (or B-DNA to S-DNA) transition at Fc≃65 pN. The base pair distance increases by approximately 70% for this value of the force [26–28,36], resulting in a sudden elongation of the molecule. In addition, the unzipping transition occurs at a critical force Fcin the range 10–20 pN [24,28,32] (of the same order of magnitude of the critical force for RNA hairpins in Ref. [16]).Atthis value of the force, many base pairs of double-stranded DNA (dsDNA) are “opened,” and the length suddenly increases. As has been already pointed out [29,33,34], this behavior is a signature of a first-order transition, with its associated region of metastability [32]. In this paper we add an external force to our previous oscillator-spin model [10] and analyze the resulting forceextension curves. Qualitatively these curves have the same features as those of the force-extension curves measured in experiments with DNA, RNA, and other biopolymers described above. At subcritical temperatures, our spin-oscillator system has a first-order phase transition at a critical force Fc with the oscillator rest state as its order parameter. We find that the DNA force-extension curves correspond to cycling at different rates the curves of the first-order phase transition. As in the experiments, we find a region of metastability in a certain range of forces, close enough to Fc. Moreover, the residence time spent at the basin of attraction of both the stable and metastable states obey an Arrhenius law: Its logarithm is proportional to F−Fc. A note about terminology: We will use the terms folded or zipped to denote the shortest value of the length, and unfolded or unzipped to denote the longest one. These terms are equivalent in the model, but they should be properly “translated” to the concrete physical situation at hand. For instance, folded and unfolded are most appropriate in the context of protein folding, although they are also used in the unzipping of RNA hairpins, while zipped or unzipped is the most adequate for the transition from double-stranded to single-stranded DNA. The rest of the paper is as follows. In Sec. II the oscillatorspin model is motivated in a biological context, and its equilibrium properties analyzed. The dynamical behavior of the model is analyzed in Sec. III. The oscillator obeys Newton’s second law with a mean-field force due to the coupling with the spins. The latter flip stochastically following Glauber dynamics at temperature T[42]. This causes the oscillator position to become a stochastic process. In the limit of fast spins, the spin-oscillator coupling gives rise to a nonlinear friction term, which drives the oscillator to equilibrium. In Sec. IV, we present Monte Carlo simulations of the system and analyze them in the light of the effective potential acting on the oscillator. Section Vcontains the main conclusions of the present work. II. THE MODEL We consider a one-dimensional chain of length L0. There is a large number N+1 of “internal” degrees of freedom sitting at regularly spaced lattice sites that are modeled by Ising spins. Thus, the distance between spins is d0=L0/N (Fig. 1). Assume that we stretch the chain so that its length becomes L=L0+. For the sake of simplicity, we will assume that the spins are regularly spaced after the stretching, so that the distance between two neighboring spins changes to d=d0+/N. This assumption amounts to a “mean-field” approximation. The potential energy of the system is V(,σ)=1 2mω22+J() N+1  i=1 σiσi+1,(1a) where J()=J0−μ (1b) is a function of . The potential Vcontains a harmonic macroscopic elastic term mω2/2 and a spin energy arising from a nearest-neighbor interaction. The spin coupling constant J depends linearly on the separation between sites, and equals J0 for the initial chain length L0. This simple choice is reasonable for L. The interaction between nearest neighbor spins mimics (in a very simple way) the short-ranged interaction between the internal degrees of freedom of complex biological molecules like nucleic acids. We assume that both J0and μ are positive. Let us define x=−J0 μ,(2) such that J(x) vanishes for x=0. For x<0 (folded state) the interaction between the spins is ferromagnetic, whereas for x>0 (unfolded state) it is antiferromagnetic. In terms of x,V becomes V=1 2mω2x2+Fcx−μx N+1  i=1 σiσi+1,F c=mω2J0 μ(3) except for an irrelevant additive constant. The parameter Fc has the dimensions of a force. If an external load Fis applied d0 L0 FIG. 1. Sketch of the model described in the main text. The spin representing the internal degree of freedom at each lattice site is shown. 021919-2 SPIN-OSCILLATOR MODEL FOR THE UNZIPPING OF ... PHYSICAL REVIEW E 86, 021919 (2012) to the system, a new term −Fx is added to (3) so that the potential energy is now V=1 2mω2x2−Hx −μx N+1  i=1 σiσi+1,H=F−Fc.(4) This potential energy is the same as introduced in Ref. [10–12], except for the extra term −Hx. Interestingly, in the “zero field” case, our system has a second-order phase transition at a critical temperature Tc, given by [10–12] Tc=μ2(N+1) mω2kB ,(5) where kBis Boltzmann’s constant. The rest position xis the order parameter of the transition: For T>T c,x=0, whereas for T<T cthere are two equally probable equilibrium states with rest positions x=±x0(x0>0). The potential energy (4) has the key ingredients to model DNA/RNA behavior in unfolding or refolding experiments. For T<T cand applied force F<F c,H<0 stabilizes the solution with rest state −x0(folded state), whereas forF>F c, H>0 and the stable solution has rest state +x0(unfolded state). The main effect of the “external field” His that the system undergoes a first-order phase transition at H=0(F= Fc)[43]. Thus, Fcis the a critical value of the force F:Atany temperature T<T cthe oscillator rest position as a function of the force changes abruptly at F=Fc. Furthermore there is a region of metastability around Fc, as discussed in Sec. II A. To analyze our model, it is convenient to render its equations of motion dimensionless first. Let ω−1be the time unit. The elastic and the spin term in the potential are of the same order if the scale of the oscillator position is [x]= μ(N+1)/(mω2). The force and the spin term are of the same order provided the scale of Fis [F]=μ(N+1), and therefore the order of magnitude of the potential is [F][x]=μ2(N+ 1)2/(mω2)=(N+1)kBTc. This scaling is reasonable because makes Vextensive. For the critical temperature Tcto be size-independent, we must assume that μscales as (N+1)−1/2 in the limit of large system size [10–12]. Thus we can define nondimensional variables according to x∗=x/[x], t∗=t/[t], V∗=V/[V],..., where the units [x], [t], [V],... are as defined in Table I. The dimensionless potential is V∗=V (N+1)kBTc =x∗2 2−H∗x∗−x∗ N+1 N+1  i=1 σiσi+1, (6a) H∗=F∗−F∗ c,(6b) with F∗ c=Fc μ(N+1) =mω2J0 μ2(N+1) =J0 kBTc .(7) TABLE I. Nondimensional units and parameters. xtF Vθ μ(N+1) mω2 1 ωμ(N+1) (N+1)kBTcT Tc We will drop the asterisks in the following (so as not to clutter our formulas), and from now on every expression will be written in terms of the dimensionless variables and parameters. A. Equilibrium state: Effective potential In equilibrium, the joint probability distribution for the oscillator position xand the spin configuration σ= {σ1,...,σ N+1}is the canonical distribution, which, in nondimensional variables, is Peq(x,σ)=1 Zexp[−(N+1)V(x,σ)/θ].(8) Here Zis a normalization constant. Let us study the equilibrium values of the oscillator position x. Then we sum over the spin variables to obtain the marginal distribution probability Peq(x)= σ Peq(x,σ)=1 ˜ Zexp[−(N+1)Veff(x)/θ],(9) where ˜ Z=2NZ.InEq.(9),Veff is an effective potential for the xvariable, Veff(x)=x2 2−Hx −θln cosh x θ,(10) whose minima will be the stable equilibrium values of x. Therefore, xeq =H+tanh xeq θ(11) gives the oscillator rest position xeq in equilibrium as a function of the dimensionless external field H=F−Fcand temperature θ. For H=0, we recover the model analyzed in Refs. [10–12], in which xeq =0 is always a solution for any θ.Forθ>1, it is the only solution; it corresponds to a maximum of Peq and is therefore stable. At θ=1twonew stable equilibria corresponding to two different maxima of Peq bifurcate from that having xeq =0. For θ<1, the positions of these maxima are ±x(0) eq .Asθ→1−,wehave x(0) eq ∼3(1 −θ).(12) On the other hand, for θ→0+,wehavex(0) eq →1, which is the maximum value of x(0) eq >0. For H=0, the two (positive or negative) equilibrium rest positions ±x(0) eq are equiprobable because they correspond to equally deep minima of the effective potential (10), which is an even function of xfor H=0. For H= 0 and temperatures below critical, θ<1, the field term −Hxin Eq. (10) breaks the symmetry between x>0 and x<0, and the effective potential is no longer an even function of x.IfH<0, the field term favors the negative branch x<0, since it gives a negative contribution to the effective potential. Therefore, we expect to find the system in the “folded” state x<0 for low values of the loading force, FFc.Onthe contrary, for H>0, the energy of the positive branch x>0 will be lowered by the field term, stabilizing it. Thus, the system will be in the “unfolded” state x>0 for high values of the loading force, FFc. 021919-3 A. PRADOS, A. CARPIO, AND L. L. BONILLA PHYSICAL REVIEW E 86, 021919 (2012) Differentiating Eq. (11) with respect to H, we find ∂Hxeq =1−1 θsech2xeq θ−1 .(13) The right-hand side of this expression is the reciprocal of the second derivative of the effective potential, and therefore we can rewrite (13) as ∂2 xV|x=xeq ∂Hxeq =1.(14) Then the effective potential has a local minimum, and the corresponding equilibrium rest position is stable for ∂Hxeq > 0, while the equilibrium rest state is unstable for ∂Hxeq <0. Since ∂Hxeq =θ/(θ−1) for xeq =0, the zero rest position of the oscillator is stable at θ>1 and unstable at subcritical temperatures θ<1. The equilibrium position xeq is a monotonically increasing function of the applied load F(or the applied field H=F− Fc) for supercritical temperatures, θ>1. There is a unique value of xeq for each value of the applied force Fif θ>1, and it is stable. As shown in Fig. 2,His nonmonotonic as a function of xeq for subcritical temperatures θ<1: It has a local maximum at −xband a local minimum at +xb, with cosh2xb θ=1 θ.(15) For each value of Hbetween −Hband +Hb, with Hb=xb−tanh xb θ,(16) eq eq FIG. 2. (Color online) Plot of Hvs xeq for θ=0.5. The values ±Hbbetween which there are three equilibrium values of the oscillator position are shown. For a given value of H,−Hb<H< +Hb, the two locally stable equilibrium points xLand xR(big green circles) and the unstable one xU(small red circle) are indicated. The two symmetric stable equilibrium points ±x(0) eq corresponding to the zero field case are also shown. The qualitative shape of the curve is the same for all the subcritical temperatures θ<1. there are three possible values of xeq:xUbetween −xband xb is therefore unstable, while the other two, xL<−xband xR> xb, are locally stable. The absolute minimum of the potential corresponds to the value of xhaving the same sign as the applied field H. The other local minimum, with sgn(x)= sgn(H), is a metastable state in the thermodynamic sense. Then we expect to find bistability in the system for |H|<H b, i.e., for a given range of loadings |F−Fc|<H baround the “critical” force Fc. III. DYNAMICS The Hamilton equations of motion corresponding to the nondimensional Hamiltonian function H(x,p,σ)=p2 2+V(x,σ),(17) with potential energy given by Eq. (6a),are ˙ x=p, ˙ p=−∂xV(x,σ),(18) so that ¨ x=−x+H+1 N+1 N+1  i=1 σiσi+1.(19) According to the stochastic Glauber dynamics, the spins flip at a rate [42] Wi(x,σ)=α 21−γ(x) 2σi(σi−1+σi+1),(20a) γ(x)=tanh 2x θ.(20b) Here Wi(x,σ) is the transition rate from configuration σto Riσ, the same as σexcept for the sign of the ith spin. Since the oscillator evolution equation (19) includes a term that depends on the stochastically changing spin configuration, the oscillator position becomes an stochastic process. The average values of the spin correlations, Ci,n =σiσi+n,(21) satisfy the system of equations dCi,n dt =α−2Ci,n+1 2γ(x)(Ci,n−1+Ci,n+1 +Ci−1,n+1+Ci+1,n−1),(22) for n⩾1, with the boundary condition Ci,0=1. If the oscillator position xwere time independent, the spins would reach the equilibrium distribution corresponding to the constant value x. Then the average spin correlations Ci,n=[tanh(x/θ)]n would be independent of idue to the spatial translation invariance. Something similar occurs in the limit of large system size, N1. Both xand the spin correlations Ci,n become macroscopic self-averaging variables; i.e., they tend to their respective “macroscopic values,” xand  Cn(independent of i), which coincide with their averages and are the most probable values of the corresponding stochastic variables [44]. Splitting both xand the correlations Ci,n in their corresponding macroscopic and fluctuating parts, x=x+δx, Ci,n = Cn+δCi,n,(23) 021919-4 SPIN-OSCILLATOR MODEL FOR THE UNZIPPING OF ... PHYSICAL REVIEW E 86, 021919 (2012) where δx and δCi,n are O(N−1/2), Eqs. (19),(22), and (23) yield the evolution equations, ¨ x=−x+H+ C1,(24a) d Cn dt =α[−2 Cn+γ(x)( Cn−1+ Cn+1)],n⩾1,(24b) to be solved with  C0=1 and appropriate initial conditions. In Eq. (24b) we have neglected terms or order 1/N such as (δx)2,δxδCi,n,etc. As the spins represent the internal degrees of freedom of the molecule, we assume that they evolve rapidly compared with the time scale of the macroscopic degree of freedom modeled by the oscillator position. Thus, the dimensionless characteristic attempt rate satisfies α1 (recall that the unit of time has been chosen as ω−1). Thus, we can solve approximately the system of Eqs. (24b) using a power series in the small parameter α−1[45–47], a procedure akin to the Hilbert method in kinetic theory. If xwere time independent, the spin correlations  Cnwould reach the equilibrium values corresponding to x,  Cn,eq =[ C1,eq]n, C1,eq =tanh x θ,(25) in the long time limit. The leading order correction of this result is [10]  C1= C1,eq −τd C1,eq dt ,(26) where τis the spins average relaxation time [45,46] τ=1 2α 1+ C1 2 eq 1− C1 2 eq2.(27) Equation (26) does not depend on the initial condition for  C1 which is forgotten after a time much shorter than the oscillator natural period. Inserting Eqs. (25)–(27) into (24a), we obtain d2x dt2+1 2αθ 1+tanh2x θ 1−tanh2(x θ) dx dt +x−H−tanh x θ=0, (28) which can be rewritten in terms of the nondimensional effective potential (10) and the friction coefficient R(x)=1+tanh2(x θ) 1−tanh2(x θ),(29) as d2x dt2=−V eff(x)−1 2αθ R(x)dx dt .(30) This approximate evolution equation gives the dynamics of the macroscopic value xof the oscillator position for fast spins: The nonlinear friction term drives the system towards equilibrium, which corresponds to the minima of Veff.Both the “renormalization” of the potential to Veff and the nonlinear friction term are a consequence of the coupling between the oscillator and the internal degrees of freedom. Equation (30) ceases to hold as θ→0+because the spin relaxation time given by (27) diverges. A detailed discussion on this point can be found in Refs. [10,11]. In the remainder of the paper, we will restrict ourselves to a temperature range for which the spins change rapidly compared to the oscillator motion and Eq. (30) holds. A. Metastability region For H=0 and subcritical temperatures, θ<1, the effective potential has two equally deep minima at the symmetric positions ±x(0) eq of Sec. II, and a maximum at x=0. For |H|< Hb(see Fig. 2for a qualitative picture), the effective potential has two minima at xR>0 and xL<0 and a metastability region appears. The globally stable position satisfies xiH>0 (i=R,L), while the other one is a metastable state in the thermodynamic sense. It must be stressed that this bistability is present for all the subcritical temperatures θ<1, it is not limited to a region near the critical temperature. Let us analyze in more detail the situation for weak fields. Expansion of Eq. (11) in powers of H=F−Fc1gives xR,L =±x(0) eq +χH +O(H2),(31) where x(0) eq is given by the solution of Eq. (11) for zero field, and χis the zero-field “susceptibility” χ=∂HxeqH=0=1−1 θ+x(0) eq 2 θ−1 .(32) Therefore, xR−xL=2x(0) eq to lowest order in H. Close to the critical temperature, θ→1−, substitution of Eq. (12) into (32) yields χ∼θ/[2(1 −θ)]. On the other hand, the unstable position xUchanges from zero to xU=−1 θ−1H+O(H2).(33) According to (9), the transitions between the two minima are hindered by the presence of large energy barriers, BR,L =(N+1)[Veff(xU)−Veff(xR,L)],(34) which are proportional to the system size N+1. For H=0, both equilibrium states have the same energy barrier, B(0) =−(N+1)Veffx(0) eq >0.(35) For H= 0, the barrier corresponding to the stable equilibrium point verifying xiH>0(i=R,L) is larger than the one for the metastable equilibrium point. The barrier from the metastable state tends to zero as |H|→Hbbecause xUand the metastable equilibrium position (xLfor H>0, xRfor H<0) coalesce in that limit. The maximum effective potential at xU separates the basins of attraction of the minima at xLand xR. The system spends long periods of time oscillating in the vicinity of either xLor xRuntil it is able to hop to the other minimum via a thermally activated process. The residence times in each basin of attraction increase exponentially with N+1, and we have to consider a system of moderate size in order to see hopping between the two minima on a reasonable time scale. Let us estimate the barrier height for weak fields. From Eqs. (10) and (31),wehave Veff(xR,L)=Veffx(0) eq ∓Hx(0) eq +O(H2),(36) 021919-5 A. PRADOS, A. CARPIO, AND L. L. BONILLA PHYSICAL REVIEW E 86, 021919 (2012) while Veff(xU)=O(H2). In Eq. (36) and for the rest of this section, the upper sign corresponds to xRand the lower sign to xL. The respective barriers from the equilibrium states xR,L defined in Eq. (34) are BR,L ≃B(0) ±(N+1)Hx(0) eq .(37) The residence times in the respective basins of the minima should have the Arrhenius form [32], τR,L =τ0exp(BR,L/θ),(38) where τ0is some characteristic time. By inserting Eq. (37) into (38), we get τR,L =τcexp ±(N+1)x(0) eq θH,(39) where τc=τ0exp(B(0)/θ) is the residence time for zero field in each basin (the same for them both). Equation (39) is the main result of this section; we should stress that it is valid for weak fields H1, but (N+1)Hcan be of the order of unity or even a large number. For H>0, we have τR>τ L, because xR is the globally stable state, while for H<0itisτR<τ L, since xRis the metastable state in that case. Interestingly, the ratio of the average lifetimes τR/τLgives the so-called equilibrium constant Kfor folding and unfolding at that force value [16]. Therefore, for our model we arrive at K≡τR τL =exp 2(N+1)x(0) eq θH,H=F−Fc,(40) ln Kis a linear function of the applied force F. A completely analogous behavior has been observed experimentally [16]. The exponent in Eq. (40) is simply N+1 times the difference of values of the effective potential (which plays the role of the free energy per particle) between both states, as readily seen by making use of Eq. (36). IV. NUMERICAL RESULTS We have performed Monte Carlo simulations of the system dynamics introduced above. In all the cases presented here, the dimensionless temperature has been chosen to be θ=0.9. The qualitative shape of the equilibrium Hversus xcurve is similar to the one shown in Fig. 2.Forθ=0.9, the oscillator rest position at equilibrium and zero field is ±x(0) eq ≃±0.525, as given by Eq. (11). Interestingly, the approximation in Eq. (12), which is expected to be valid very close to the critical temperature, gives quite a good estimate x(0) eq ≃0.547. The points at which the curve Hversus xhas either a maximum or a minimum are ±xb=±0.295; the corresponding values of the field are ∓Hb=∓2.15 ×10−2. There is metastability for applied fields in the interval |H|<H b. A. Force-extension curves Let us analyze the force extension curves of the model. The pulling-pushing cycle is as follows. We start at the equilibrium configuration corresponding to xmin <0 (folded state), so that the initial value of the field is Hmin =xmin −tanh(xmin/θ) given by Eq. (11). We pull the system by a stepwise increment of the field: At each step, the field is increased by H, and then the system is allowed to evolve during a given time t.At the end of this period, we record the oscillator position x. Then we increase again the field by H and continue the process in the same vein. The pulling process ends when we reach a positive value of the field Hmax =−Hmin. Then we start to push back, decreasing the field by H at each step, until we reach again the minimum field value Hmin. As during the pulling process, we record the value of xat fixed Hafter the evolution time t. This process is completely analogous to that carried out in unzipping experiments with biomolecules. Figure 3shows some typical pulling-pushing cycles with different loading rates H/t.Wehaveusedasystemwith N+1=1000 spins, and a spin attempt rate α=4, large enough for the spins to be fast as compared to the oscillator [10]. For all the curves, the minimum value of the oscillator position is xmin =−1.5, and the field increment at each step is H =10−3, which is smaller than Hb, and it allows the system to visit the metastability region. The system behavior is qualitatively similar for other parameter values, as long as the temperature θ<1. The loading rate is changed by varying the amount of time t at each step. The red (solid) lines correspond to the unfolding process (t =5, 100, and 104from top to bottom), and the green (dashed) lines to the refolding process (t =5, 100, and 104from bottom to top). These numerical curves are qualitatively similar to those observed in unzipping experiments with nucleic acids [16]. There is always some hysteresis as the unfolding and folding curves are not superimposed. The area of the hysteresis cycle increases with the loading rate, being large for the largest loading rate considered and almost zero for the smallest one. The main difference with the experimental results is that, in our model, the extension of the molecule is not linked to a dropping of the loading force (in order to see this effect in a real experiment, see, for instance, Figs. 2(A) and 2(E) of Ref. [16]). In the experiments of Ref. [16], the total length between the beads localizing the molecule is controlled. This x H 10.50-0.5-1 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 FIG. 3. (Color online) Hysteresis cycles for different values of the loading rate H/t. The red (solid) lines correspond to the unfolding process (t =5, 100, and 104from top to bottom), and the green (dashed) lines to the refolding process (t =5, 100, and 104from bottom to top). The Hvs xcurve at equilibrium is plotted with blue stars. 021919-6 SPIN-OSCILLATOR MODEL FOR THE UNZIPPING OF ... PHYSICAL REVIEW E 86, 021919 (2012) corresponds to fixing the length Lin our model, not the load given by Has we have done. When the force is externally controlled in experiments, there is no drop of the loading force at the extension transition, and the hysteresis cycle described by the molecule extension is completely analogous to ours [25]. Our numerical results can be explained as follows. As depicted in Fig. 2and observed in the last paragraph of Sec. II A,for|H|>H bwe have a unique stable equilibrium point for |H|>H b, the zipped state xeq <−xb<0forH< −Hb<0 and the unzipped state xeq >+xb>0forH>H b. For |H|<H b, there are two (locally) stable rest positions xL<0 and xR>0 and an unstable position xUbetween them. The true thermodynamic equilibrium state of the system corresponding to the global minimum of the potential satisfies xiH>0(i=L,R), whereas the metastable state has rest position such that xiH<0. Let us consider again the pulling processes in Fig. 3, starting from the equilibrium configuration corresponding to a low value of the applied force, H<0, |H||Hb|. When we increase the force at a moderate rate, the oscillator follows the equilibrium curve for negative values of H, with xeq <0, because its relaxation time is small compared to t, and the friction term in Eq. (30) can drive it to equilibrium. When the field reaches H=0 (or a very small value), the stable equilibrium position of the oscillator changes discontinuously from −x(0) eq to x(0) eq . The important question is now whether t insufficiently long for the oscillator position to overcome the energy barrier B(0) at zero field, given by Eq. (35).Ifthe answer is positive, the system jumps during the time interval t to the other stable branch where x>0, and it stays on that branch when the force is further increased. In the pushing back experiment, the system reverses its path and the behavior is almost reversible. This behavior is observed for the largest value t =104, corresponding to a loading rate H/t = 10−7. For larger pulling rates, such as the other ones considered in the same figure, the system does not have enough time to jump over the energy barrier at H=0. Therefore, xmoves over the metastable branch with xiH<0, until the barrier decreases sufficiently for the oscillator position to jump to the most stable branch, with x>0. Of course, the actual part of the metastable region visited by the system depends on the loading rate. For the highest loading rate considered, corresponding to t =5, the system visits the whole metastable branch up to the maximum. A similar line of reasoning explains the behavior observed in the refolding curve. It is interesting to note that the hysteresis cycle found for high loading rates is not a nonequilibrium behavior, as previously suggested [16–19], but it arises from the sampling of the regions of metastability for subcritical temperatures. A pulling experiment corresponding to a rate even slower than the smallest loading rate in Fig. 3is plotted in Fig. 4. Again, the minimum value of the oscillator position has been chosen to be xmin =−1.5, and the field increment at each step H =10−3, but the time spent by the system at each value of the force is very large, namely, t =105. For the sake of clarity, the region around the critical force F=Fc(H=0) has been zoomed in. We observe several jumps between the folded (x<0) and unfolded (x>0) states for |H|=|F−Fc|< 5×10−3. This behavior is a clear signature of the bistability x H 10.50-0.5-1 0.02 0.015 0.01 0.005 0 -0.005 -0.01 -0.015 -0.02 FIG. 4. (Color online) Detail of the metastability region |H|⩽ Hbfor a pulling experiment with a very slow loading rate, H =10−3 and t =105. Hopping between the zipped and unzipped state is clearly seen for |H|=|F−Fc|<5×10−3. shown by our system in the region |H|<H b. At zero field, the energy barriers between the two stable states are identical, τR=τL=τcin (39), and the reverse transition is equally likely. For t τc, the oscillator position has enough time to surpass the energy barriers several times, and it may go back and forth from one state to the other, as shown in Fig. 4. B. Constant force experiments We have also carried out Monte Carlo simulations at constant force. We have chosen θ=0.9, a small value of the field, H=3×10−3<H b, and a smaller size, N+1= 500, in order to keep the simulation time under control. Equation (11) gives the two locally stable oscillator rest positions xL=−0.509 (metastable) and xR=0.540 (globally stable), separated by the unstable oscillator position xU= −0.027. These values agree with the weak field expressions (31)–(33). Figure 5shows a time trace of the system. The oscillator jumps stochastically between the two values xLand xRcorresponding to the two locally stable equilibrium points. The system spends more time in the stable state xR(recall H=3×10−3>0), because it has to surpass a larger energy barrier so as to escape therefrom. In order to check Eqs. (39) or (40) for the residence times in each basin of attraction, we have measured the average time spent in each basin for different values of the applied field in the metastability region |H|<H b. We find that ln τR,L increases linearly with H, with a slope that agrees with the theoretical prediction of Eq. (39). We have plotted the ratio of average lifetimes K=τR/τLdefined in Eq. (40) as a function of the applied field Hin Fig. 6. Therein, ln Kshows a linear behavior, similarly to that seen in actual experiments [16]. The slope m= dln K/dH obtained numerically, m=542.5 agrees well with the theoretical prediction calculated from Eq. (40),m=583.8. This result strongly supports the physical picture developed in Sec. III A for the behavior of the system in the metastability region, both from a qualitative and a quantitative point of view. 021919-7 A. PRADOS, A. CARPIO, AND L. L. BONILLA PHYSICAL REVIEW E 86, 021919 (2012) t x 2500020000150001000050000 1 0.5 0 -0.5 -1 FIG. 5. (Color online) Time trace of the oscillator position for a constant force experiment in the metastability region, namely, H= 3×10−3. Hopping between the two locally stable equilibrium points of the oscillator is observed. V. CONCLUSIONS We have modeled biomolecules’ folding and unfolding under an external load force by a macroscopic linear oscillator coupled with internal degrees of freedom represented by Ising spins that undergo Glauber dynamics. The simple mean-field character of the model prevents us from doing quantitative comparisons with the real experiments, and we have to settle for qualitative comparisons. We cannot simulate position-controlled experiments, only force-controlled ones. To overcome the limitations of this simple system, we should consider more realistic models. One possibility is to use spin-string models, in which a discretized string is coupled to Ising-like variables. These models do not have a mean-field character because the couplings between the string variables and the spins are local, thereby exhibiting a much richer and more complex phenomenology [9]. In the context of DNA unzipping, it is worth considering extensions of the H ln K 0.010.0050-0.005-0.01 6 4 2 0 -2 -4 -6 FIG. 6. (Color online) Logarithm of the ratio K≡τR/τL (squares) of the residence times for the two equilibrium states xR,L in the metastability region. The number of spins is N+1=500. Also plotted is the best fit to the Arrhenius law (solid line) [Eq. (40)]. well-known Poland-Scheraga and Peyrard-Bishop models [29,33], which take into account the base pairing between the different DNA strands. This ingredient, which is missing in the spin oscillator model presented here, is crucial for understanding the importance of the entropic effect due to bubble formation in DNA denaturation. Despite its shortcomings, the present spin-oscillator model is an useful minimal model to understand experiments in which biomolecules are unzipped by mechanical force. This simple model allows us to obtain exact analytical results from which a physically appealing and general picture arises. The hysteresis cycles show that the system exhibits a metastable equilibrium behavior in the unzipping experiments, not a true nonequilibrium behavior, as was suggested previously [16–19]. In this regard, the unfolding and refolding cycles are quite different from the truly nonequilibrium hysteresis cycles exhibited by glass formers in cooling and heating processes (see Refs. [46,47] and references therein). In the cooling process, the glass formers depart from the equilibrium curve and end in a far from equilibrium state at low temperatures. In the reheating process, they return to equilibrium following a curve different from the cooling one, which typically overshoots the equilibrium curve. In contrast with this behavior, the hysteresis cycles in our system arise because the metastable equilibrium branches of the xversus Hcurve are swept at high loading rate (such that the system does not have enough time to find the true minimum of the potential in the bistable field interval |H|<H b). Due to the shape of the equilibrium xversus H curve in Fig. 2, the system unzips at a higher value of the field than the one at which it rezips. We expect that this general picture remains valid for more realistic models and/or actual biomolecules. As in the experiments, the system hops between the zipped and the unzipped state at very small loading rates. It has then enough time to surpass the energy barrier separating the unfolded and the folded states. When the force is held constant with |H|<H b(bistable region), the system hops back and forth between the two possible equilibrium values of the oscillator position. The average lifetimes show an Arrhenius-like dependence on the applied field H=F−Fc, again in agreement with the behavior of real systems. For all subcritical temperatures T<T c, these behaviors occur due to the first-order transition and its associated region of metastability. The phase transition is a consequence of the coupling between the oscillator and the one-dimensional spin system, which introduces an effectively long-range interaction between the spins. Similar hidden one-dimensional long-range effective correlations enable phase transitions in biological systems [48]. ACKNOWLEDGMENTS This research has been supported by the Spanish Ministerio de Econom´ ıa y Competitividad through Grants FIS201124460 (A.P., partially financed by FEDER funds), FIS201128838-C02-02 (A.C.), FIS2011-28838-C02-01 (L.L.B.), and FIS2010-22438-E (Spanish National Network Physics of Outof-Equilibrium Systems), and by UCM/BSCH CM 910143 (A.C.). 021919-8 SPIN-OSCILLATOR MODEL FOR THE UNZIPPING OF ... PHYSICAL REVIEW E 86, 021919 (2012) [1] J. Feder and E. Pytte, Phys.Rev.B8, 3978 (1973). [2] P. A. Rikvold, Z. Phys. B 26, 195 (1977). [3] P. A. Rikvold, Z. Phys. B 30, 339 (1978). [4] A. Boisen, Nat. Nanotechnol. 4, 404 (2009). [5] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987). [6] L. L. Bonilla and F. Guinea, Phys. Rev. A 45, 7718 (1992). [7] M. Hofheinz, H. Wang, M. Ansmann, R. C. Bialczak, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, J. Wenner, J. M. Martinis, and A. N. Cleland, Nature (London) 459, 546 (2009). [8] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang, M. Weides, J. Wenner, J. M. Martinis, and A. N. Cleland, Nature (London) 464, 697 (2010). [9] L. L. Bonilla, A. Carpio, A. Prados, and R. R. Rosales, Phys. Rev. E 85, 031125 (2012). [10] A. Prados, L. L. Bonilla, and A. Carpio, J. Stat. Mech. (Theor. Exp.) (2010) P06016. [11] L. L. Bonilla, A. Prados, and A. Carpio, J. Stat. Mech. (Theor. Exp.) (2010) P09019. [12] A. Prados, L. L. Bonilla, and A. Carpio, AIP Conf. Proc. 1332, 232 (2011). [13] W. L. Bragg and E. J. Williams, Proc. R. Soc. London A 145, 699 (1934). [14] H.-X. Zhou and R. Zwanzig, J. Phys. Chem. A 106, 7562 (2002). [15] F. Ritort, J. Phys.: Condens. Matter 18, R531 (2006). [16] J. Liphardt, B. Onoa, S. B. Smith, I. Tinoco, Jr., and C. Bustamante, Science 292, 733 (2001). [17] J. Liphardt, S. Dumont, S. B. Smith, I. Tinoco, Jr., and C. Bustamante, Science 296, 1832 (2002). [18] F. Ritort, in Progress of Mathematical Physics vol. 38 (Poincar´ e Seminar 2003), J. Dalibard and V. Rivasseau eds. (Birkh¨ auser Verlag, Basel, 2004), p. 193–227. [19] D. Collin, F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco, Jr., and C. Bustamante, Nature (London) 437, 231 (2005). [20] M. C. Williams and I. Rouzina, Curr. Opin. Struct. Biol. 12, 330 (2002). [21] M. Manosas and F. Ritort, Biophys. J. 88, 3224 (2005). [22] J.-D. Wen, M. Manosas, P. T. X. Li, S. B. Smith, C. Bustamante, F. Ritort, and I. Tinoco, Jr., Biophys. J. 92, 2996 (2007). [23] M. Manosas, J.-D. Wen, P. T. X. Li, S. B. Smith, C. Bustamante, I. Tinoco, Jr., and F. Ritort, Biophys. J. 92, 3010 (2007). [24] J. M. Huguet, N. Forns, and F. Ritort, Phys.Rev.Lett.103, 248106 (2009). [25] J. M. Huguet, Ph.D. thesis, University of Barcelona, 2010. [26] P. Gross, N. Laurens, L. B. Oddershede, U. Bockelmann, E. J. G. Peterman, and G. J. L. Wuite, Nature Phys. 7, 731 (2011). [27] S. B. Smith, Y. Cui, and C. Bustamante, Science 271, 795 (1996). [28] M. Rief, H. Clausen-Schaumann, and H. E. Gaub, Nat. Struct. Biol. 6, 346 (1999). [29] D. K. Lubensky and D. R. Nelson, Phys.Rev.Lett.85, 1572 (2000). [30] C. Bustamante, S. B. Smith, J. Liphardt, and D. Smith, Curr. Opin. Struct. Biol. 10, 279 (2000). [31] D. Anselmetti, J. Fritz, B. Smith, and X. Fernandez-Busquets, Single Mol. 1, 53 (2000). [32] S. Cocco, R. Monasson, and J. F. Marko, Proc. Natl. Acad. Sci. USA 98, 8608 (2001). [33] D. Marenduzzo, S. M. Bhattacharjee, A. Maritan, E. Orlandini, and F. Seno, Phys.Rev.Lett.88, 028102 (2001). [34] D. Marenduzzo, A. Maritan, E. Orlandini, F. Seno, and A. Trovato, J. Stat. Mech. (Theor. Exp.) (2009) L04001. [35] J. van Mameren, P. Gross, G. Farge, P. Hooijman, M. Modesti, M. Falkenberg, G. J. L. Wuite, and E. J. G. Peterman, Proc. Natl. Acad. Sci. USA 106, 18231 (2009). [36] A. Fiasconaro and F. Falo, arXiv:1205.6126. [37] T. E. Fisher, A. F. Oberhauser, M. Carrion-Vazquez, P. E. Marszalez, and J. M. Fernandez, Trends Biochem. Sci. 24, 379 (1999). [38] T. E. Fisher, P. E. Marszalek, and J. M. Fernandez, Nat. Struct. Biol. 7, 719 (2000). [39] A. F. Oberhauser and M. Carrion-Vazquez, J. Biol. Chem. 283, 6617 (2008). [40] B. Bullard, T. Garcia, V. Benes, M. C. leake, W. A. Linke, and A. F. Oberhauser, Proc. Natl. Acad. Sci. USA 103, 4451 (2006). [41] C. Cecconi, E. A. Shank, C. Bustamante, and S. Marqusee, Science 309, 2057 (2005). [42] R. J. Glauber, J. Math. Phys. 4, 294 (1963). [43] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics Vol. 5, Statistical Physics Part 1 (Pergamon, Oxford, 1980). [44] N. G. Van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1997). [45] J. J. Brey and A. Prados, Phys. Rev. E 47, 1541 (1993). [46] J. J. Brey, A. Prados, and M. J. Ruiz-Montero, J. Non-Cryst. Solids 172–174, 371 (1994). [47] A. Prados, J. J. Brey, and B. S´ anchez-Rey, Phys.Rev.B55, 6343 (1997); A. Prados and J. J. Brey, Phys.Rev.E64, 041505 (2001). [48] M. Peyrard, Nature Phys. 2, 13 (2006). 021919-9