Theory of force-extension curves for modular proteins and DNA hairpins

PHYSICAL REVIEW E 91, 052712 (2015) Theory of force-extension curves for modular proteins and DNA hairpins L. L. Bonilla,1A. Carpio,2and A. Prados3 1G. Mill´ an Institute, Fluid Dynamics, Nanoscience and Industrial Mathematics, Universidad Carlos III de Madrid, 28911 Legan´ es, Spain 2Departamento de Matem´ atica Aplicada, Universidad Complutense de Madrid, 28040 Madrid, Spain 3F´ ısica Te´ orica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain (Received 22 July 2014; revised manuscript received 29 March 2015; published 26 May 2015) We study a model describing the force-extension curves of modular proteins, nucleic acids, and other biomolecules made out of several single units or modules. At a mesoscopic level of description, the configuration of the system is given by the elongations of each of the units. The system free energy includes a double-well potential for each unit and an elastic nearest-neighbor interaction between them. Minimizing the free energy yields the system equilibrium properties whereas its dynamics is given by (overdamped) Langevin equations for the elongations, in which friction and noise amplitude are related by the fluctuation-dissipation theorem. Our results, both for the equilibrium and the dynamical situations, include analytical and numerical descriptions of the system force-extension curves under force or length control and agree very well with actual experiments in biomolecules. Our conclusions also apply to other physical systems comprising a number of metastable units, such as storage systems or semiconductor superlattices. DOI: 10.1103/PhysRevE.91.052712 PACS number(s): 87.15.La,72.25.Dc,73.63.Hs,75.50.Pp I. INTRODUCTION Nowadays technological advances allow manipulation of single molecules with sufficient precision to study many mechanical, kinetic, and thermodynamic properties thereof. Recent reviews of techniques used and results obtained in single-molecule experiments (SMEs) can be found in Refs. [1– 3]. In these experiments, typically the force applied to pull the biomolecule is recorded as a function of its end-to-end distance, thereby producing a force-extension curve (FEC). This FEC characterizes the molecule elasticity and provides information about its processes of folding and unfolding [4–11]. In the following, the end-to-end distance of the biomolecule is referred to as the total length [12]. The force-extension curves are different depending on whether the total length or the force are controlled. When the total length of the protein is used as a control parameter (length control), the unfolding transition is accompanied by a drop in the measured force and a sawtooth pattern is the typical force-extension curve [5,7–9,13]. When the force is the control parameter (force-control), unfolding of several or all single protein domains may occur at a constant value of the force [14]. Other questions are related to the rate at which the control parameter (length or force) sweeps the force-extension curve: depending on the loading rate, stochastic jumps between folded and unfolded protein states may be observed [1,8,10,13,14]. The analysis of the force vs extension curves provides valuable information about the polyprotein, the DNA, or the RNA hairpin. Let us consider atomic force microscope (AFM) experiments in which a modular protein comprising a number of identical folds (modules or units) is pulled at a certain rate (length control) [5,7]. The typical value of the force Fcat which the unfolding takes place is related to the mechanical stability of the units: a larger value of the force is the signature of higher stability. Nevertheless, it should be stressed that the unraveling of a domain is a stochastic event and occurs for forces within a certain range. A second feature of the sawtooth FEC is the spacing between consecutive force peaks. This spacing is directly related to the difference of length between the folded and unfolded configurations of one unit. This is the reason that the peaks of the FEC of artificially engineered modular proteins are regularly spaced. A typical example is I278, composed of eight copies of immunoglobulin domain 27 from human cardiac titin. The spacing between peaks for this protein is 28.4±0.3 nm at an unfolding force of 204 ±26 pN [5,7]. This length increment is found by fitting several peaks of the FEC with the worm-like chain (WLC) model of polymer elasticity [15,16]. More recently, force-controlled AFM experiments with a I27 single-module protein have been reported [17,18]. These experiments provide data free from the module to module variations that even an artificially engineered polyprotein has. Berkovich et al. have interpreted their results using a simple Langevin equation model that includes an effective potential with two minima for a range of the applied force [17]. The thermodynamics of pulling experiments is well established under both force and length control. For controlled force, the relevant thermodynamic potential is a Gibbs-like free energy, whereas for controlled length it is a Helmholtz-like free energy [10,19,20]. Interestingly, the sawtooth structure of the FEC of biomolecules is already present at equilibrium, as shown very recently in a simple model with a Landau-like free energy [21]. However, the control parameter in real experiments with biomolecules (force or length) changes usually with time at a finite rate [1,3,7,8,10,13,14,22]. Knowledge of these dynamical situations is not as complete as in the equilibrium case. Under force control, we can write a Langevin equation (or the associated Fokker–Planck equation) in which noise amplitude and effective friction are linked by a fluctuation-dissipation relation, as done in Refs. [17,23]. On the other hand, under length control, the situation is more complex: the force is no longer a given function of time but an unknown that must be calculated by imposing the length constraint. This has lead to the proposal of simple dynamical algorithms such as the quasi-equilibrium algorithm of Ref. [10]. While being successful in reproducing experimentally observed behavior, these algorithms do not correspond to the integration of well-defined evolution equations. 1539-3755/2015/91(5)/052712(19) 052712-1 ©2015 American Physical Society L. L. BONILLA, A. CARPIO, AND A. PRADOS PHYSICAL REVIEW E 91, 052712 (2015) In some cycling experiments, the biomolecule is switched between the folded and unfolded configurations at a certain switching rate [1,8,10,13,14,22]. After Liphardt et al., we call such a process a stretching-relaxing or an unfolding-refolding cycle [8,13]. The unfolding typically occurs at a force Fu that is larger than the refolding force Fr. Therefore, some hysteresis is present and, moreover, the unfolding (refolding) force typically increases (decreases) with the pulling rate. A reversible curve in which Fu=Fr=Fcis only observed for a sufficiently small rate. Some authors have claimed that this is a signature of irreversible nonequilibrium behavior and thus used these experiments to test nonequilibrium fluctuation theorems [13,24–26]. On the other hand, for a simple model for which only the force-controlled situation could be analyzed [27], it has been found that the observed behavior in biomolecules can be understood as the system sweeping a certain part of the metastable equilibrium region of the FEC that surrounds Fc. In this way, the system is exploring metastable minima of the system free energy landscape. One of the main goals of this work is to determine if this physical picture also holds for length-controlled experiments. In this paper, we add two important ingredients of real biomolecule pulling experiments to a simple model with independent domains and Landau-like free energy whose equilibrium analysis is given in Ref. [21]. We add: (i) dynamical effects and (ii) interacting units. Dynamical effects are introduced by means of Langevin or Fokker–Planck equations, both under force and, most interestingly, length control. Therefrom, we can carry out a systematic investigation of the dynamical FEC when the control parameter (force or length) is varied at a finite rate. The simplest way to introduce interaction between modules is via a harmonic potential trying to drive them to global equilibrium. In this way, the creation of bubbles, that is, regions of unfolded modules inside regions of folded ones, has a free energy cost. This is expected to be most relevant for systems in which the unfolding and refolding of units is mainly sequential, as in the unzipping of DNA hairpins [1]. Interestingly, the complex and force-sensitive behavior of polyproteins observed in force-clamp experiments has been recently explained by sequential unfolding [28]. The main ingredients of our model are bistability of protein modules and, in the length-controlled case, a global constraint that introduces a long-range interaction among modules. These features are quite general in physics, as they appear in many different fields. For instance, many particle storage systems such as the storage of lithium in multiparticle electrodes of rechargeable lithium-ion batteries [29,30], air storage in interconnected systems of rubber balloons [31], or voltage biased weakly coupled semiconductor superlattices [32–36]. Throughout the paper, the analogies and differences that arise in these different physical situations will be discussed. The rest of the paper is as follows: The model we use is described in Sec. II, in which we write down both the Langevin and the Fokker–Planck equations in Secs. II A and II B, respectively. In Sec. III, we investigate an ideal modular protein comprising many identical, noninteracting, units. In Sec. III A, we show that the equilibrium FEC corresponding to our Landau-like double-well free energy has multiple branches. Statistical mechanics considerations determine the stability of the equilibrium branches for: (a) force control in Sec. III B, and (b) length control in Sec. III C. We also consider dynamical situations when the control parameter (either force or length) varies at a finite rate. Section IV deals with a real chain, in which the nearest-neighbor modules interact via an extra harmonic term. First, we study the equilibrium situation in Sec. IV A, in which we show that the size of the branches is reduced, as compared to the ideal case. Sections IV B–IV D analyze the changes that the dynamics brings to the equilibrium picture by considering deterministic dynamics, quenched disorder and finite temperature dynamics (thermal noise), respectively. Final remarks are made in Sec. V. Appendix Aexplains unfolding and refolding under length control using a more realistic potential, whereas Appendixes B and C deal with some technical aspects not covered in the main text. II. MODEL To be specific, let us consider AFM unfolding of modular proteins: They are stretched between the tip of the microscope cantilever and a flat, gold-covered substance (platform), whose position is externally controlled. The forces acting on the molecule bend the cantilever which, in turn, determines the applied force with pN precision. See Fig. 1 of Refs. [3]or [7] for an idealized situation. In force-controlled experiments with a single-module protein, the free energy of an extending protein comprises at least two distinct components: an entropic term that accounts for chain elasticity and an enthalpic component that includes the short-range interactions arising between the neighboring amino acids as the protein contracts [17,18]. In a certain force range, these two components cause the single-module free energy to have two minima, corresponding to the folded and unfolded states of the domain [17]. Let us consider a system comprising Nmodules. The jth module extends from xjto xj+1, so that its extension is ηj= xj+1−xj,j=1,...,N. The configuration η={ηj}defines the polyprotein state at a mesoscopic level of description. When isolated, the free energy of the jth unit is a(ηj;Y,δj), a double-well potential whose minima correspond to the folded and unfolded states discussed above. Yis the set of relevant intensive parameters, like the temperature Tand the pressure pof the fluid (thermal bath) surrounding our system. The parameter δjaccounts for the slight differences from unit to unit: δj=0∀j, if all units are identical, and thus no quenched disorder is present in the system. As part of the tertiary structure of the polyprotein, modules are weakly interconnected by linkers in a structure-dependent way [37]. It seems reasonable that this weak interaction acts on the unfolding and refolding timescale and tries to bring the extensions of the modules to a common value, corresponding to global mechanical equilibrium. For the sake of simplicity, we model the linkers as harmonic springs. Thus, the system free energy Afor a given configuration of module extensions ηis A(η;Y)= N  j=1 a(ηj;Y,δj)+ N  j=2 kj(Y) 2(ηj−ηj−1)2.(1) 052712-2 THEORY OF FORCE-EXTENSION CURVES FOR MODULAR . . . PHYSICAL REVIEW E 91, 052712 (2015) If all the linkers are identical, kj=kfor all j=2,...,N, and the elastic constants may depend on the intensive parameters. The length Lof a polyprotein in a configuration ηis L(η)= N  j=1 ηj.(2) The experiments are carried out at either force-controlled or length-controlled conditions. First, we analyze case (i), in which a certain external force F=F(t) is applied to the ends of the protein or DNA hairpin. For a detailed discussion of how this is achieved in real experiments, see, for example, Refs. [14] (Chap. 6) and [38] for the optical-tweezers case, and Refs. [39,40] for the AFM case. In our simplified theoretical approach, we only have to add the term Ufc (η;F)=−FL(η)=Fx 1−Fx N+1(3) to the free energy A(η). In this way, we obtain a Gibbs free energy G(η;Y,F)=A(η,Y )+Ufc =A(η,Y )−FL(η), G(η;Y,F)= N  j=1 g(ηj;Y,F,δj)+ N  j=2 kj(Y) 2(ηj−ηj−1)2, (4a) g(ηj;Y,F,δj)=a(ηj;Y,δj)−Fη j.(4b) Note that we are not taking into account the limited bandwidth of the feedback device that controls force in real experiments; we are assuming that the desired force program F(t)is perfectly implemented. Second, we investigate the lengthcontrolled situation, case (ii). For a schematic representation of the experimental situation, see, for instance, Fig. 1 of Ref. [3]: The length L(t) between the base of the cantilever and the platform is the externally controlled quantity. On the other hand, the cantilever tip deflects a certain distance x from its base, such that x +L(η)=L(t). If the stiffness (spring constant) of the cantilever is χlc, we have an extra harmonic term in the potential Ulc =χlc(x)2/2, that is, Ulc (η;L)=χlc 2[L(η)−L(t)]2.(5) Therefore, an extra force Flc =−χlc[L(η)−L(t)] acts over each unit, trying to keep the polyprotein length equal to L(t): The larger χlc, the better the length control is, as expected on intuitive grounds and explicitly shown in Ref. [24]. In this paper, with the exception of Appendix A, we assume perfect length control, that is, we consider the limit χlc →∞that implies L(η)−L(t)→0 over the time evolution of the units and a finite value of the corresponding extra force Flc.In other words, Flc tends to a limiting value Fthat depends on the prescribed length L(t). This unknown value is the force required to attain the total length L(t), and it has to be calculated by imposing the constraint iηi=L. The effect of this limit on the relevant thermodynamic potential for the length-controlled case A+Ulc shall be discussed in the section on Fokker–Planck description of the dynamics. Finally, we would like to stress that the present model has some similarities with the more complicated one proposed by Hummer and Szabo for the unfolding of polyproteins several years ago; see Appendix C of Ref. [41]. In addition to the module extensions ηj, these authors consider the module centers of mass rjas independent unknowns. These variables interact through a quadratic potential that yields a linear restoring force whenever rj+1−rjdeparts from ηj+1+ηj 2. The site potential for the module extensions is the sum of a WLC potential and a harmonic potential [41], instead of the double-well potential we consider in the main text or the asymmetric potential we consider in Appendix A. Moreover, Hummer and Szabo introduce a WLC linker that connects the polyprotein with the length-controlling device, which is absent in our model. A. Langevin dynamics The extensions ηjobey coupled Langevin equations with the appropriate thermodynamic potential. The friction coefficient and the amplitude of the white noise are related by a fluctuation-dissipation theorem. The source for both the friction and the stochastic force is the fluid that the modules are immersed in, which is assumed to remain in equilibrium at temperature T. We assume that the modules’ inertia can be neglected and thus their evolution equations are overdamped, γj˙ηj=F−∂ ∂ηj A(η;Y)+2Tγ jξj(t),(6a) ξj(t)=0,ξj(t)ξl(t)=δjlδ(t−t),j=1,...,N. (6b) Here γjis the friction coefficient for the jth module, and we measure the temperature in units of energy (kB=1). In general, the friction coefficients γjmay depend on the system configuration η, if hydrodynamic interactions play a significant role in the considered unfolding scenario. For the sake of simplicity, we do not consider this possibility in the present paper. In this respect, it is interesting to remark that, in the more complicated model of Ref. [41], module centers of mass and extensions satisfy Langevin equations with different extension-dependent diffusion coefficients. Our presentation of the model above implies that the Langevin equations (6) are valid both in force-controlled and length-controlled experiments, but (i) in force-controlled experiments, F=F(t) is the known force program, whereas (ii) in length-controlled experiments we have L(η)=L(t). We differ the discussion on the experimental situation with an “imperfect” length control (because of the finite value of the stiffness χlc of the device controlling the length) to the next section on the equivalent Fokker–Planck description of the dynamics. For perfect length control, F(t) is determined by imposing the constraint L(η)=L, which yields F=γ N⎛ ⎝ dL dt + N  j=1 1 γj ∂A(η;Y) ∂ηj− N  j=12T γj ξj⎞ ⎠,(7a) γ−1=1 N N  j=1 γ−1 j.(7b) The parameter γis an average friction coefficient. In the case of identical units, γj=γ∀j. We split Fin two terms, a 052712-3 L. L. BONILLA, A. CARPIO, AND A. PRADOS PHYSICAL REVIEW E 91, 052712 (2015) “macroscopic term” FFP and a “fluctuating term” F ,as follows: F=FFP +F, (8a) FFP =γ N⎡ ⎣ dL dt + N  j=1 1 γj ∂A(η;Y) ∂ηj⎤ ⎦,(8b) F =−γ N N  j=12T γj ξj.(8c) We prove in Sec. II B that FFP is the force appearing in the flux term of the Fokker–Plack equation. Note that for any N,F =0 and then F=FFP. Furthermore, F is a sum of Gaussian variables, and thus its statistical properties are completely given by its first two moments. It can be easily shown that F (t)F (t)=N−1γδ(t−t), its variance tends to zero as N−1, which is the typical behavior of fluctuating quantities in statistical mechanics. Even so, it should be noted that in biomolecules Nis not necessarily very large and certainly not of the order of Avogadro’s number, and thus fluctuations play a major role. In force-extension experiments, the length is usually uniformly increased or decreased with time t,dL/dt =μwith a constant μ. It is convenient to render our equations dimensionless. We set the length unit [η] equal to the difference between the extensions of the two free energy minima of a single unit for a certain applied force. It is natural to adopt the critical force, at which the two minima are equally deep, as the unit of force: [F]=Fc. The parameters [η] and [F] depend on the specific choice of the double-well potential a(η;Y,0). The free energy unit is then [F][η]. We select the timescale as [t]=γ[η]/[F], where γis the typical friction coefficient experienced by the units. The typical value of γcan be obtained from the value of the diffusion coefficient D=T/γ of a single module protein being stretched [18]. In principle, we introduce a new notation for the dimensionless variables, F∗=F/[F], etc., but, in order not to clutter our formulas, we drop the asterisks in the remainder of the paper. B. Fokker–Planck equation and equilibrium distributions In force-controlled experiments, F(t) is a given function of time, and the set of Langevin equations (6a) is equivalent to the following Fokker–Planck equation for the probability density P(η,t) of finding the system with extension values η={η1,...,η N}at time t: ∂ ∂tP= N  j=1 1 γj ∂ ∂ηj∂G ∂ηj P+T N  j=1 1 γj ∂2P ∂η2 j .(9) where G=A−FL, as given by Eq. (4a). If the force Fis kept constant, Eq. (9) has a stationary solution, which is the statistical mechanics prescription, Peq (η)∝e−G(η;Y,F)/T .(10) Therefore, the equilibrium values of the module extensions ηeq are the functions of Fthat maximize Por, equivalently, minimize G; that is, they verify ∂G ∂ηjY,F eq =0⇒ηj=ηeq j(Y,F),j=1,...,N.(11) If there is only one minimum, this is the equilibrium configuration. If there is more than one, the absolute minimum is the thermodynamically stable configuration, while the other minima correspond to metastable states in the thermodynamic sense. For each equilibrium configuration, either stable or metastable, the equilibrium value of the free energy Gis Geq(Y,F)=G(ηeq(Y,F); Y,F).(12) Taking into account Eq. (11), we have ∂Geq ∂F Y=∂G ∂F η,Y eq =− N  j=1 ηeq j(Y,F)=−Leq(Y,F), (13) which gives the equilibrium FEC under force control. Let us consider now the length-control situation. In the experiments, the device controlling the length of the system does not have an infinite stiffness and thus the length control is not perfect, as discussed above (see also Ref. [24] and Appendix A). Had we taken into account this finite value of the stiffness χlc, the Fokker–Planck equation would have been obtained by substituting the Gibbs free energy G≡ A+Ufc in Eq. (9) by the corresponding thermodynamic potential A+Ulc. Thus, the stationary solution of this Fokker–Planck equation would be the equilibrium distribution Peq(η;Y,L)∝exp{−[A(η;Y)+Ulc(η;L)]/T}. Of course, in the limit as χlc →∞, the variance of the Gaussian factor exp[−Ulc(η;L)/T] vanishes and this factor tends to a delta function δ(L(η)−L), giving perfect length control. In the case of perfect length control, the correct Fokker– Planck equation can be obtained by taking the limit as χlc → ∞, but here we follow an alternative route. We calculate the first two moments of the extensions η, taking into account that not all the extensions ηjare independent and that the force F is given by Eq. (7a), ∂ ∂tP= N  j=1 1 γj ∂ ∂ηj∂A ∂ηj−FFPP +T N  j=1 1 γj N  k=1δjk −γ Nγk∂2 ∂ηj∂ηk P.(14) Here FFP is given by Eq. (8b). If the length is kept constant, dL/dt =0, Eq. (14) has a stationary solution, Peq (η;Y,L)∝δ(L(η)−L)e−A(η;Y)/T ,(15) as can be easily verified by inserting Eq. (15) into Eq. (14). This means that Ais the relevant potential for the statistical mechanics description at equilibrium, as was expected. Equation (15) is consistent with the limit as χlc →∞of the equilibrium distribution for realistic length control, as already discussed above. To obtain the equilibrium values for the extensions, we look for the minima of Awith the constraint given by the delta function in Eq. (15), L(η)=L. We have to introduce a 052712-4 THEORY OF FORCE-EXTENSION CURVES FOR MODULAR . . . PHYSICAL REVIEW E 91, 052712 (2015) Lagrange multiplier Fand look for the minima of A−FL; that is, the same minimization as in the force-controlled case. However, the Lagrange multiplier is an unknown that must be calculated at the end of the process by imposing the constraint, F=F(L). This Lagrange multiplier is, from a physical point of view, the force that must be applied to the system in order to have the desired length. The equilibrium extensions ηeq j(L) are thus given by the solutions of ∂A ∂ηjYeq =F, j =1,...,N; N  j=1 ηeq j(Y,F)=L. (16) The last equation gives the FEC, L=L(Y,F)orF=F(Y,L), from which we obtain ηeq =ηeq(Y,L). The thermodynamic potential Aeq is the Legendre transform of Geq with respect to F. In fact, the equilibrium value of A,Aeq(Y,L)= A(ηeq(Y,L); Y), verifies that ∂Aeq (Y,L) ∂L Y=F. (17) The proper variables for Aeq are the set of intensive parameters Y(temperature T, pressure p,... of the fluid in which the polyprotein is immersed) and the extensive length L[42], while the proper variables for Geq are all intensive, Yand F. In this sense, Aeq plays the role of Helmholtz free energy, while Geq is the analogous to the Gibbs free energy. It should be stressed that, (i) however, different notations are found in the literature for these two thermodynamic potentials; and (ii) as in the case of magnetic systems [43], there is a difference of sign with respect to the usual free energy terms with the pressure pand the volume V. The fluctuation theorems for Markov processes described by the Langevin (or the equivalent Fokker–Planck) equations have been thoroughly analyzed in Ref. [44]. The results therein are directly applicable to the Fokker–Planck equations derived here for the force-controlled and the realistic (finite χlc) length-controlled cases. When the controlled parameter (either force or length) is kept constant (time-independent), detailed balance applies and the corresponding stationary distributions are equilibrium (canonical) ones (in the terminology of Sec. II of Ref. [44]). In the limit χlc →∞, we expect this result to be still valid on physical grounds, but further mathematical work would be necessary to establish it rigorously: Some of the matrices defined in Ref. [44] become singular and thus have no inverse. This is a point that certainly deserves further investigation, but it is out of the scope of the present paper. III. THE IDEAL CHAIN In this section, we analyze the case of an ideal chain, in which the identical units do not interact either among themselves or with the cantilever or platform, kj=0 and δj=0 for all j[21]. We analyze the equilibrium situation and thus solve the minimization problems for the force-controlled and length-controlled cases of the previous section. We also investigate the dynamical situation arising in processes in which the force or length varies in time at a finite rate and compare these dynamical FECs to the equilibrium ones. A. Double-well potential; equilibrium branches In order to keep the notation simple, we omit the dependence on the intensive parameters Yof the free energy parameters. As a minimal model, we consider the polynomial form, ` a la Landau, for the free energy [21]: A(η)= N  j=1 a(ηj),a(η)=Fcη−αη2+βη4.(18) The parameters Fc,α, and βare all positive functions of the intensive parameters Y. Specifically, Fcplays the role of the critical force, above (below) which the unfolded (folded) configuration is the most stable one, as shown in what follows. The possible equilibrium extensions ηeq are the minima of a(η)−Fη, a(η(i))=F, (19) or, equivalently, −2αη(i)+4β(η(i))3=ϕ, ϕ ≡F−Fc.(20) We have introduced the notation η(i)because Eq. (20) has three solutions in the metastability region, given by |ϕ|=|F− Fc|<ϕ 0=(2α/3)3/2β−1/2. We set the indexes by choosing η(1) <η (2) <η (3). They depend on the force Fthrough ϕ (and on the intensive variables Ythrough {α,β,Fc}). The extensions η(1)(ϕ) and η(3)(ϕ) are locally stable because they correspond to minima of aj−Fη j, while η(2)(ϕ) corresponds to a maximum and is therefore unstable. The curvatures at the folded and unfolded states are χ(i)(ϕ)=a(η(i)(ϕ))= 12β[η(i)(ϕ)]2−2α,i=1,3. Both curvatures (i) are positive in the metastability region and (ii) vanish at their limits of stability, χ(1) (χ(3))atϕ=ϕ0(ϕ=−ϕ0). The situation is similar to that analyzed by Landau [45]for a second-order phase transition under an external field, with η and ϕ=F−Fcplaying the role of the order parameter and the external field, respectively. At the critical force ϕ=0, the stable equilibrium values of the extensions are η(3) c=−η(1) c=α 2β1/2 .(21) They are equiprobable, since g(η)=a−Fη is an even function of ηfor F=Fc, and g(1) =a(1) −Fcη(1) =g(3) = a(3) −Fcη(3), where we have introduced the notation a(1) ≡ a(η(1)), a(3) ≡a(η(3)), g(1) ≡g(η(1)), and g(3) ≡g(η(3)). For F= Fc, the “field” ϕfavors the state with sgn(ϕ)=sgn(η). In fact, at the limit of stability we have that g(1) =−13α2(6β)−1 for ϕ=−ϕ0[or g(3) =−13α2(6β)−1for ϕ=ϕ0]. Therefore, in the metastability region |ϕ|<ϕ 0, we have the following picture: For F<F c, the thermodynamically stable state is the folded one η(1) <0 and the unfolded one η(3) >0is metastable. For F>F c, the situation is simply reversed. On the other hand, the folded η(1) (unfolded η(3)) state also exists for forces below (above) the metastability region ϕ<−ϕ0 (ϕ>ϕ 0). In their respective regions of existence, both locally stable extensions η(1) and η(3) are increasing functions of ϕ(or F), since Eq. (19)impliesthatχ(k)(ϕ)dη(k)/dϕ =1. At zero force, one module can be folded or unfolded if ϕ0>F c, while we have only the folded state if ϕ0<F c. 052712-5 L. L. BONILLA, A. CARPIO, AND A. PRADOS PHYSICAL REVIEW E 91, 052712 (2015) 0.6 0.4 0.2 0.0 0.2 0.4 0.6 1.5 1.0 0.5 0.0 0.5 1.0 1.5 L Lc 0 FIG. 1. (Color online) Normalized FECs for N=8 (solid black) and N=20 (dashed red). Zero length corresponds to having half of the units unfolded at F=Fc.ThereareN+1 branches in the metastability region |ϕ/ϕ0|<1, with the number of unfolded units, J, increasing from left to right. The first (J=0) and last (J=N) branches are independent of N. Note that the branches become denser as Nincreases, and also the up-down and left-right symmetry thereof. These symmetries stem from the simple form of the Landau-like free energy (18), and thus they are not present if a more realistic potential is considered; see Appendix A. Any module can be either folded or unfolded in the metastability region, and thus a FEC with N+1 branches shows up, as seen in Fig. 1.TheJth branch of the F−L curve corresponds to Junfolded modules and N−Jfolded ones, J=0,...,N. Since there is no coupling among the units, the equilibrium value of Aover the Jth branch is Aeq J=(N−J)a(1) +Ja(3).(22a) The corresponding length is LJ=(N−J)η(1) +Jη(3).(22b) Both Aeq Jand LJare functions of Fand the intensive parameters Ythrough the equilibrium extensions. Equation (22b)isthe FEC, both for the forceand length-controlled cases. In Fig. 1, we have normalized the lengths with Lc=LN(Fc)−L0(Fc)=Nη(3) c−η(1) c,(23) which is the difference of lengths between the completely unfolded branch (J=N) and the completely folded one (J=0) at the critical force. It is interesting to note that similar multistable equilibrium curves appear in quite different physical systems: from storage systems [29–31] to semiconductor superlattices [32–36]. For instance, see Fig. 3 of Ref. [29] and Fig. 6 of Ref. [30] for the chemical potential vs charge curve in storage systems, and Fig. 8.13 of Ref. [35]forthe current-voltage curve of a superlattice. As discussed in the previous section, we have chosen [F]= Fc=1 and [η]=η(3) c−η(1) c=1 as units of force and length. Using Eq. (21), β=2αand η(3) c=−η(1) c=1/2. Moreover, the folded state η(1) is the most stable one at zero force. This means that the unstable state η(2) is closer to the metastable state η(3) for the simple Landau potential we are using [46]. For the sake of concreteness, we take η(2) −η(1) =0.9(η(3) −η(1)) at zero force, which leads to α=2733/2/1672 ≈2.697 787 and ϕ0=913/2/836 =1.038 378. It should be stressed that all the normalized plots in this section are independent of this particular choice of parameters. A more conventional definition of protein length could be to select at zero force (a) zero extension for the folded modules (b) the difference between the unfolded and folded configurations as the length unit. This “physical” definition would give a nondimensional extension u=η−η(1) (F=0) η(3) (F=0)−η(1) (F=0),(24a) and a polyprotein length Lu=−Nη(1) (F=0)+N j=1ηj η(3) (F=0)−η(1) (F=0)=19N 30 +√273 15 L, (24b) respectively. At zero force, the module extensions are u(1) =0, u(2) =0.9, and u(3) =1. The length Luis typically positive for F>0, that is, for ϕ>−Fc. We expect that the simple Landau-like free energy given by Eq. (18) should be relevant to investigate qualitatively the FECs for forces or lengths close to the metastability region. In particular, this minimal choice does not account for the existence of a maximum length of the polymer, its so-called contour length [12], a fact that becomes significant for high forces. In order to study the whole range of forces and/or try to describe quantitatively the experiments, we should use a more realistic potential, such as that proposed by Berkovich et al. [17] or modifications thereof. We did this in Ref. [28] to understand the stepwise unfolding observed in force-clamp experiments. The simpler potential used in this paper suffices for (i) showing that the key aspects of the experimental behavior observed in the unfolding or refolding region can be understood within a minimal model, and (ii) establishing connections with other physical systems such as storage devices [29–31] or semiconductor superlattices [32–36] that have similar behavior in the metastability region. We briefly investigate an asymmetric potential in Appendix A to understand why the experimentally observed FEC corresponding to unfolding under length control is reproduced by the Landau-like potential whereas the FEC corresponding to refolding is not; see Sec. III C for details. B. Force control In force-controlled experiments, the Gibbs free energy is the relevant thermodynamic potential because it appears in the equilibrium distribution (10). As discussed in Sec. II B, the stable state corresponds to the absolute minimum of G. All the units in our ideal chain are independent under force control. Therefore, by increasing quasistatically the force, the equilibrium FEC (22b) is swept. Over the Jth branch with J unfolded modules, Geq J=(N−J)g(1) +Jg(3),g (i)=a(i)−Fη(i).(25) For F<F c=1(F>F c), the absolute minimum of G corresponds to the folded (unfolded) state η(1) (η(3)) and the system moves over the force-extension branch in which none (all) of the units are unfolded, J=0(J=N). 052712-6 THEORY OF FORCE-EXTENSION CURVES FOR MODULAR . . . PHYSICAL REVIEW E 91, 052712 (2015) Unfolding is a first-order phase transition between these states that occurs at the critical force Fc=1 defined by continuity of forces and of the Gibbs free energies, Geq 0|Fc=Geq N|Fc.AtFc=1, all the units unfold simultaneously. The length, which is a function of Fgiven by Eq. (13), has a discrete jump equal to Lc, given by Eq. (23). It is worth recalling η(3) c−η(1) c=1 in nondimensional units. The free energy (25) produces d dF Geq N−Geq 0=−N(η(3)(F)−η(1)(F)) <0∀F, (26) consistently with Eq. (13). Then the basin of attraction of the completely folded branch is the largest one for F<F c, whereas the completely unfolded branch has the largest basin of attraction for F>F c. All the intermediate metastable branches with J= 0,N are not “seen” by the system in a quasistatic process that takes infinite time to occur; see the top panel of Fig. 2. For a real, nonquasistatic process, the simple equilibrium picture above is not realized. Depending of the rate of variation of the force and the strength of the thermal fluctuations, the system will explore the metastable branches of the FEC. Then intermediate states between the completely folded and unfolded configurations will be seen [1,14,27,47]. This is shown in the bottom panel of Fig. 2by solving Eqs. (6a) and (6b) with kj=δj=0 (in nondimensional form) for a 20-module protein, with T=2×10−5and T=0.02. All the T=2×10−5curves are superimposed on each other because the considered rates |dF/dt|are small enough to lead to the adiabatic limit. For the upsweeping (downsweeping) process, the system moves over the completely folded (unfolded) branch until it reaches the end thereof, ϕ=ϕ0(ϕ=−ϕ0). Then it jumps to the completely unfolded (folded) branch. The temperature is so small that the activated processes over the free energy barriers take place over a much longer timescale. For the higher temperature, T=0.02, the system can jump between the different minima of the potential and the force at which the system jumps between branches depends on the rate of variation of the force. Also, the system partially explores some of the intermediate branches. This picture is consistent: quite close to the adiabatic limit, the hysteresis cycle is large for the highest rate of variation, whereas the cycle shrinks towards the straight line ϕ=0(Fc=1) as the rate tends to zero. C. Length control In length-controlled experiments, the length constraint introduces a long-range interaction between the protein modules. The equilibrium probability of any configuration ηis now givenbyEq.(15). Then the equilibrium configuration ηeq is found by minimizing Awith the constraint (2), and the difference between values of Aeq at adjacent branches in the F−Ldiagram governs the stability thereof. The length J at which there is a change in the relative stability of two consecutive branches, with J−1 and Junfolded units, is determined by the equality of their respective free energies Aeq. The corresponding forces f− J≡FJ−1(J) and f+ J=FJ(J) over the branches with J−1 and Junfolded units obey the 0.6 0.4 0.2 0.0 0.2 0.4 0.6 1.5 1.0 0.5 0.0 0.5 1.0 1.5 L Lc 0 0.6 0.4 0.2 0.0 0.2 0.4 0.6 1.5 1.0 0.5 0.0 0.5 1.0 1.5 L Lc 0 FIG. 2. (Color online) (top) First-order transition in the length for a quasistatic increase of the force. We use different colors for the stable parts of the branches (solid black) and the metastable parts (dashed red). The first branch J=0 is swept until the critical force ϕ=0is reached. Then all the modules unfold simultaneously and the system goes directly (arrow) to the completely unfolded branch J=N=20 (bottom). Hysteresis cycles under force-controlled conditions for a N=20 system. The lines correspond to simulations of the Langevin equations (6) for the temperature T=0.02 and different rates of variation of the force, namely |dF/dt|=3×10−k, with k=2 (dot-dashed red), k=3(dottedgreen),k=4(solidblue),and k=5 (dashed orange). The same rates of variation of the force are considered for the very low temperature T=2×10−5. All the curves are superimposed and thus they are plotted with the same symbols (black dots). For the higher temperature, the area of the hysteresis cycle decreases with the rate, approaching the behavior for a quasistatic process. system of two equations Aeq J−1f− J=Aeq Jf+ J ,Leq J−1f− J=Leq Jf+ J .(27) The force rips at L=Jare Nfirst-order equilibrium phase transitions because (i) the thermodynamic potential Aeq is continuous at the transition, (ii) F=(∂Aeq/∂L)Yhas a finite jump, from f− Jto f+ J<f− Jat the Jth transition. In the top-left panel of Fig. 3, we explicitly show f− 1and f+ 1.Wehavethe following picture: As observed in Fig. 1, the branches J−1 and Jcoexist on a certain range of lengths. Inside this range, 052712-7 L. L. BONILLA, A. CARPIO, AND A. PRADOS PHYSICAL REVIEW E 91, 052712 (2015) f 1 f 1 0.6 0.4 0.2 0.0 0.2 0.4 0.6 1.5 1.0 0.5 0.0 0.5 1.0 1.5 L Lc 0 − 0.6 − 0.4 − 0.2 0.0 0.2 0.4 0.6 − 1.5 − 1.0 − 0.5 0.0 0.5 1.0 1.5 L/Lc /0 − 0.6 − 0.4 − 0.2 0.0 0.2 0.4 0.6 − 1.5 − 1.0 − 0.5 0.0 0.5 1.0 1.5 L/Lc /0 − 0.6 − 0.4 − 0.2 0.0 0.2 0.4 0.6 − 1.5 − 1.0 − 0.5 0.0 0.5 1.0 1.5 L/Lc /0 FIG. 3. (Color online) (top left) Equilibrium force rips in the F−Lcurve for a system with N=8 domains. We use different colors for the stable parts of the branches (solid black), metastable parts (dotted red), and the force rips (black arrows). The system follows the solid black curve in a quasistatic pulling process, with a series of first-order transitions in the force (marked by the arrows). At the Jth transition, the force changes from f− J[over the (J−1)st branch] to f+ J(over the Jth branch). These forces f± Jincrease with the number of unfolded units, as observed in AFM experiments with modular proteins, even though all the units are perfectly identical in the model. (top right) Hysteresis cycle for a system composed of N=8 modules. The dimensionless temperature T=0.02, and the rate of variation of the length is |dL/dt|=1.2×10−3(˙ L>0, solid blue; ˙ L<0, dashed red). (bottom left) The same as in the top-right panel, but for a smaller rate |dL/dt|=1.2×10−6. Aside from thermal fluctuations, the system almost sweeps the equilibrium curve. (bottom right) The same plot as in the bottom-left panel, but for T=2×10−5. Thermal fluctuations are so small that the system approaches the T=0 behavior, in which the branches are swept up to the end of the metastability region. Eq. (17) implies ∂ ∂L Aeq J−Aeq J−1Y=FJ(L)−FJ−1(L)<0,(28) wherewehaveusedEq.(16). At equal length values L,the force is larger on the branch with a smaller number of folded units, FJ(L)<F J−1(L)∀J. Therefore, Aeq J−1<Aeq J, and then the branch J−1 is the stable one and Jis metastable for L< J. The situation reverses for L> J, and there are not more stability changes between these branches because Aeq J− Aeq J−1decreases monotonically as a function of L,asgiven by Eq. (28). Each intermediate branch (J=1,...,N −1) is thus stable between Jand J+1, that is, between f+ Jand f− J+1 (see top-left panel of Fig. 3). A sawtooth pattern arises in the F−Lcurve, with Ntransitions between the N+1 branches at lengths 1,..., N. Similarly to the analysis for the force-controlled case, let us investigate the behavior of the system when the length is first increased and afterwards decreased with the same rate. Depending on the rate and the value of the temperature, a region of the metastable part of the branches is explored, and the force rips do not take place at the equilibrium values J. In Fig. 3, apart from the equilibrium force-extension curve (top-left panel), we plot three unfolding or refolding cycles for an ideal eight-module protein. In the top-right panel and the bottom-left panels, the temperature is T=0.02 and the rates are |dL/dt|=1.2×10−3and 1.2×10−6, respectively. For the smallest rate, the system basically sweeps the equilibrium curve, aside from thermal fluctuations. Note that some of the transitions are “blurred” because of the hopping between the two possible forces at the transition length for this very small rate. On the other hand, for the highest rate, some hysteresis is present. Finally, in the bottom-right panel, the temperature is much lower, T=2×10−5, which results in the largest hysteresis cycle. This low-temperature dynamical FEC is basically the same for the two rates considered before, |dL/dt|=1.2×10−3and |dL/dt|=1.2×10−6,so only the latter is shown. The system sweeps each branch up to its limit of (meta)stability, |ϕ/ϕ0|=±1. Interestingly, this low-temperature behavior resembles that of the chemical potential in a recent investigation of the thermodynamic origin of hysteresis in insertion batteries [29,30]. This indicates that thermal fluctuations are less relevant for insertion batteries than for modular proteins, despite the similarities in their mathematical description. In our simulations, we have averaged the force over a unit time interval to mimic the experimental situation, in which the measuring devices have finite resolution [48]. 052712-8 THEORY OF FORCE-EXTENSION CURVES FOR MODULAR . . . PHYSICAL REVIEW E 91, 052712 (2015) 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.6 0.4 0.2 0.0 0.2 0.4 0.6 L Lc 0 20 40 60 80 100 120 9.0 9.5 10.0 10.5 11.0 N N F Fc FIG. 4. (Color online) (top panel) Zoom of the metastability region for two chains with N=8 (solid line) and N=20 (dashed line). The size of the rips decreases as the number of units N increases, and it vanishes as N→∞. (bottom panel) Decrease of the force rips with the number of units N. We plot the size of the rips F =f− J−f+ Jfor the “central” transition with J=(M+1)/2, scaled with the factor N/Fc(circles). The limiting value 6√3, which represents the power law behavior given Eq. (29), is shown with a solid line. It is observed that the system approaches rapidly this asymptotic behavior, being very close to it for N20. In the unfolding process, the transitions occur at forces or lengths that are displaced upward with respect to that of the refolding process, as usually observed in the experiments [1,8,10,14,49–52]. However, the unfolding or refolding curves for the simple Landau-like quartic potential we are using are much more symmetric than the experimental ones, reflecting the symmetry of the potential. In the experiments with modular proteins, the unfolding FEC exhibits large force rips similar to ours, but the refolding FEC does not present a sawtooth pattern [49–51]. Recently, however, several force rips have been observed in the refolding of the NI6C protein, both in AFM experiments and steered molecular dynamics simulations [52]. In Appendix A, we briefly analyze the predictions of our theory for the more realistic potential introduced in Ref. [17]. In this case, the unfolding and refolding curves are strongly asymmetric and closely resemble the experimental ones. For both, equilibrium and dynamical FECs (the latter being closer to the real experimental situation), (i) the size of the force rips decrease with the number of units N, and (ii) f± J increase with the number of unfolded units Jfor moderate values of N. The equilibrium case is illustrated by the top panel of Fig. 4. Interestingly, the increase with Jof the rips forces has been observed in modular proteins [3,7](N∼10), whereas the rips forces are basically independent of Jfor nucleic acids experiments (larger N)[1,14,22]. Also, it is worth noting that it has recently been shown that the length-controlled and the flow-controlled scenarios in polymer stretching are thermodynamically equivalent [53]. Let us investigate in more detail the dependence of the force rips size with the number of units Nin the equilibrium case; see left panels of Fig. 3.ForlargeN, the free energy Aover each branch is extensive (∝N), whereas the difference of free energies over consecutive branches for a given value of the length is independent of N. Therefore, the relative free energy change between consecutive branches scales as N−1. Making use of Eq. (27) and neglecting terms of order N−3, we obtain [54] f± J−Fc ϕ0=∓3√3 N1∓rJ N,r J=2LJ LN−L0Fc .(29) Both f− Jand f+ Jincrease linearly with J, and so do rJand LJ,[LJis given by Eq. (22b)]. This is necessary to fulfill the continuity condition for the free energy at the rips. On the other hand, for very large N, the term proportional to rJis proportional to N−2and, therefore, it is small compared with the first term on the right-hand side of Eq. (29), which is proportional to N−1. As a consequence, in this limit the force rips become independent of Jand symmetrical with respect to Fc, f± J−Fc ϕ0∼∓3√3 N.(30) This is consistent with the behavior observed in nucleic acids [1,14,22], in which the number of units is much larger than that typical of modular proteins. Moreover, it shows that the rip size in equilibrium follow a simple power law; it decays as N−1for large N. We show the tendency to this power law in the bottom panel of Fig. 4, in which we plot the size of the rip F =f− J−f+ Jfor a specific value of J. We have chosen J=(M+1)/2, that is, the transition in which the number of unfolded units become larger than the number of folded ones, Jincreases from (M−1)/2to(M+1)/2(Modd). Note that the fact that limN→∞ f± J=Fcimplies that all the units of the system unfold simultaneously at the critical force Fc in the infinite-size limit. This is the expected behavior, since in the thermodynamic limit as N→∞force fluctuations disappear and the collectives with controlled force and controlled length should be utterly equivalent. IV. CHAINS WITH ELASTIC INTERACTIONS BETWEEN IDENTICAL MODULES In this section, we investigate the effect of the harmonic potential in Eq. (1) [proportional to (ηj−ηj−1)2]onthe FECs. This term tends to minimize the number of “domain walls” separating regions with folded units from regions with unfolded units, as the domain walls give a positive contribution to the free energy that is proportional to their number. This elastic interaction is expected to be more relevant in experiments in which the unfolding or refolding of units is basically sequential, as in the case of unzipping and rezipping of DNA and RNA hairpins. The harmonic potential does not completely prevent the formation of “bubbles,” regions of unfolded units inside a domain of folded ones, but adds a free energy cost 052712-9 L. L. BONILLA, A. CARPIO, AND A. PRADOS PHYSICAL REVIEW E 91, 052712 (2015) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.2 0.4 0.6 0.8 L F FIG. 11. (Color online) FECs for the BGMKUF potential, with N=8 (solid red) and N=15 (dashed blue). There are N+1 branches in the metastability region Fm<F <F M, with the number of unfolded units Jincreasing from left to right. The first (J=0) and last (J=N) branches are independent of N,asinFig.1.Note the asymmetry of the branches with respect to the critical force Fc (dot-dashed line). to 1 ms) to mimic the finite resolution of the measuring devices in real experiments. Due to the asymmetry of the equilibrium branches, the unfolding and refolding curves are quite different, as seen in experiments. A clear sawtooth pattern is present in the unfolding curve: The molecule clearly sweeps a certain part of each equilibrium branch until it reaches a length at which it jumps to the neighboring branch. Similarly to experimental observations, this jump is associated with a decrease in the force (force rip) [3,5,7,37]. On the other hand, in the refolding process, the curve is much smoother and it is much more difficult to identify the intermediate branch that the system is sweeping, at least for the first stage of the relaxation curve (here, for L0.2). This is analogous to the usual experimental behavior in the refolding process 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 L F FIG. 12. (Color online) Unfolding-refolding cycle for a modular protein with eight units with free energies given by the BGMKUF potential. The dot-dashed lines correspond to the equilibrium branches of the FEC; see Fig. 11. There are sharp force rips in the unfolding process, each corresponding to the unfolding of one of the units (jumps between neighboring branches). In the refolding process, there are no sharp peaks until the length has almost completely relaxed, L0.2. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 L F FIG. 13. (Color online) Unfolding-refolding cycle for a modular protein with eight units with free energies given by the BGMKUF potential and a finite value of the cantilever stiffness. The force F is plotted against the end-to-end distance of the molecule L(η). The equilibrium branches of the FEC (see Fig. 11) are the dot-dashed lines. Both the unfolding and refolding curves are very similar to those in Fig. 12, except for the force rips in the unfolding process not being perfectly vertical as a consequence of the imperfect length control. [49–51]. However, there appear clearer traces of force peaks in the refolding FEC when the molecule has partially relaxed (L0.2). This behavior resembles the FECs obtained for the NI6C protein in Ref. [52]; see Figs. 1C, 1D, and S5 therein. Although the previous unfolding-refolding cycle is very similar to those observed in experiments, it may be argued that our ideal length-control device may have some impact on the observed behavior. Therefore, we consider now a more realistic length-control device, such as that depicted in Fig. 1 of Ref. [3], which leads to the length-control potential term in Eq. (5), where χlc is the (finite) spring constant of the cantilever. A typical value of the spring constant for an AFM experiment is 6pN/nm, which gives a dimensionless value χlc =1.8. First, it is important to stress that the equilibrium branches of the FEC are not changed by the finite stiffness of the length-controlling device. The equilibrium extensions ηiare given by Eq. (19), a(ηi)=F, but now F=−χlc[L(η)−L] is the force exerted by the finite-stiffness control device. Metastability appears in the same range of applied forces as in the case of ideal length control; the only difference is that the end-to-end distance L(η) does not equal L, instead, L(η)=L−F/χlc <L. In other words, the tip of the cantilever has an equilibrium deflection x =F/ξlc for each considered force F. Repeating the unfolding-refolding process in Fig. 12, with the only difference of the finite value of the stiffness, we have obtained the results shown in Fig. 13. The unfolding-refolding cycles in both figures are very similar, although they would not match perfectly when superimposed. To obtain complete agreement with the perfect length control situation shown in Fig. 12, we should have employed a larger value of the spring constant, around 150 pN/nm or χlc =45. In particular, the refolding curve is again much smoother than the unfolding sawtooth pattern found with either the quartic or the BGMKUF potential, but with some minor upward traces for L0.2. 052712-16 THEORY OF FORCE-EXTENSION CURVES FOR MODULAR . . . PHYSICAL REVIEW E 91, 052712 (2015) APPENDIX B: EQUIVALENT ISING MODEL FOR THE FREE ENERGY MINIMA We can write down the length and Gibbs free energy (37) in an Ising-like manner. Let us assign a spin-down variable to the folded units, so that σj=−1ifηeq j,0=η(1), and a spin-up σj=+1 to the unfolded ones, with ηeq j,0=η(3). The number of unfolded units and domain walls are J= N  j=1 1+σj 2,M= N−1  j=1 1−σjσj+1 2.(B1) Except for an additive constant, the free energy (37b) becomes Geq (σ)=−H N  j=1 σj− N−1  j=1 σjσj+1+O(k2),(B2) an Ising system with an external field Hand ferromagnetic nearest-neighbor coupling given by H=g(3) −g(1) 2,=k[η(3) −η(1)]2 4>0.(B3) Interestingly, a similar expression for the free energy was proposed in Ref. [71]. The sign of Hdetermines which minimum of the Gibbs free energy g(η) is deepest, η(1) or η(3); at the critical force Fc=1, that is, H=0, they are equally deep. The ferromagnetic coupling ∝kfavors the configurations with domains of parallel spins and thus a minimal number of domain walls for a given number of unfolded units J[72]. Then M=0, when all the units are either folded or unfolded, or M=1, when there are both folded and unfolded units, produce the minimum free energy (B2). Given Eq. (B1), the length of the system at equilibrium is Leq(σ)=N 2(η(1) +η(3))+(N−1) +η(3) −η(1) 2 N  j=1 σj− N−1  j=1 σjσj+1,(B4) where =k[χ(3) −χ(1)] 2χ(1)χ(3) .(B5) The parameter can be positive or negative. For the simple quartic potential we are considering, =0 at the critical force Fc=1, >0forF<F c, and <0forF>F c. We have not considered here the quadratic corrections, proportional to k2, which only affect sites at the domain walls and their nearest neighbors. In this equivalent Ising description, they (i) change the first-order coupling constants and , and (ii) introduce a second-nearest-neighbor interaction. Similarly, by taking into account higher-order corrections, up to order kn, we get an Ising model with longer-ranged interactions up to the nth-nearest neighbors. APPENDIX C: LYAPUNOV FUNCTION FOR THE DETERMINISTIC DYNAMICS IN THE LENGTH-CONTROLLED CASE Unlike the Gibbs free energy Gin the force-controlled case, the Helmholtz free energy A, as given by Eq. (1), is no longer a Lyapunov function of the zero-noise dynamics under lengthcontrolled conditions with a known length dependence L(t). However, ˜ A(η)=A(η)+ N  j=1k 2(ηj+1−ηj)2 −ηj NN  k=1 a(ηk)+dL dt −a(ηj)−ηja(ηj) N(C1) is a Lyapunov function in this case. In fact, the governing nondimensional equations can be written as dηj dt =− ∂ ∂ηj ˜ A(η),(C2) after eliminating Fby means of Eq. (7a). Then, d dt ˜ A(η)=− N  j=1∂ ∂ηj ˜ A(η)2 ⩽0. Also, ˜ A(η)>Nmin ua(u)−Fu−a(u)−ua(u) N for Fm<F =1 N N  j=1 a(ηj)+1 N dL dt <F M. [1] F. Ritort, J. Phys.: Condens. Matter 18,R531 (2006). [2] S. Kumar and M. S. Li, Phys. Rep. 486,1(2010). [3] P. E. Marszalek and Y. F. Dufrˆ ene, Chem.Soc.Rev.41,3523 (2012). [4] S.B.Smith,Y.Cui,andC.Bustamante,Science 271,795 (1996). [5] M. Carrion-V´ azquez, A. F. Oberhauser, S. B. Fowler, P. E. Marszalek, S. E. Broedel, J. Clarke, and J. M. Fernandez, Proc. Natl. Acad. Sci. USA 96,3694 (1999). [6] H. Lu and K. Schulten, Proteins: Struct., Funct., Genet. 35,453 (1999). [7] T. E. Fisher, P. E. Marszalek, and J. M. Fernandez, Nat. Struct. Biol. 7,719 (2000). [8] J. Liphardt, B. Onoa, S. B. Smith, I. Tinoco, and C. Bustamante, Science 292,733 (2001). [9] C. Bustamante, Z. Bryant, and S. B. Smith, Nature (London) 421,423 (2003). [10] M. Manosas and F. Ritort, Biophys. J. 88,3224 (2005). [11] Y. Cao, R. Kuske, and H. Li, Biophys. J. 95,782 (2008). [12] Thus, the total length should not be confused with the contour length of the biomolecule, which is a fixed quantity. 052712-17 L. L. BONILLA, A. CARPIO, AND A. PRADOS PHYSICAL REVIEW E 91, 052712 (2015) [13] J. Liphardt, S. Dumont, S. B. Smith, I. Tinoco, and C. Bustamante, Science 296,1832 (2002). [14] J. M. Huguet, Ph.D. thesis, Universitat de Barcelona, 2010 (unpublished). [15] C. Bustamante, J. F. Marko, E. D. Siggia, and S. Smith, Science 265,1599 (1994). [16] J. F. Marko and E. D. Siggia, Macromolecules (Washington, DC, US) 28,8759 (1995). [17] R. Berkovich, S. Garcia-Manyes, J. Klafter, M. Urbakh, and J. M. Fernandez, Biophys. J. 98,2692 (2010); Biochem. Biophys. Res. Commun. 403,133 (2010). [18] R. Berkovich, R. I. Hermans, I. Popa, G. Stirnemann, S. GarciaManyes, B. J. Bernes, and J. M. Fernandez, Proc. Natl. Acad. Sci. USA 109,14416 (2012). [19] D. Keller, D. Swigon, and C. Bustamante, Biophys. J. 84,733 (2003). [20] J. M. Rub´ ı, D. Bedeaux, and S. Kjelstrup, J. Phys. Chem. B 110, 12733 (2006). [21] A. Prados, A. Carpio, and L. L. Bonilla, Phys. Rev. E 88,012704 (2013). [22] J. M. Huguet, C. V. Bizarro, N. Forns, S. B. Smith, C. Bustamante, and F. Ritort, Proc. Natl. Acad. Sci. USA 107, 15431 (2010). [23] D. Reguera, J. M. Rub´ ı,andJ.M.G.Vilar,J. Phys. Chem. B 109,21502 (2005). [24] S. Park, F. Khalili-Araghi, E. Tajkhorshid, and K. Schulten, J. Chem. Phys. 119,3559 (2003). [25] D. Collin, F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco Jr., and C. Bustamante, Nature (London) 437,231 (2005). [26] G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. USA 107, 21441 (2010). [27] A. Prados, A. Carpio, and L. L. Bonilla, Phys. Rev. E 86,021919 (2012). [28] L. L. Bonilla, A. Carpio, and A. Prados, Europhys. Lett. 108, 28002 (2014). [29] W. Dreyer, J. Jamnik, C. Guhlke, R. Huth, and M. Gaberscek, Nat. Mater. 9,448 (2010). [30] W. Dreyer, C. Guhlke, and R. Huth, Phys. D 240,1008 (2011). [31] W. Dreyer, C. Guhlke, and M. Herrmann, Continuum Mech. Thermodyn. 23,211 (2011). [32] H. T. Grahn, R. J. Haug, W. Muller, and K. Ploog, Phys. Rev. Lett. 67,1618 (1991). [33] M. Rogozia, S. W. Teitsworth, H. T. Grahn, and K. H. Ploog, Phys. Rev. B 65,205303 (2002). [34] L. L. Bonilla and H. T. Grahn, Rep. Prog. Phys. 68,577 (2005). [35] L. L. Bonilla and S. W. Teitsworth, Nonlinear Wave Methods for Charge Transport (Wiley-VCH, Weinheim, 2010). [36] Yu. Bomze, R. Hey, H. T. Grahn, and S. W. Teitsworth, Phys. Rev. Lett. 109,026801 (2012). [37] T. Hoffmann and L. Dougan, Chem.Soc.Rev.41,4781 (2012). [38] N. Forns, S. de Lorenzo, M. Manosas, K. Hayashi, J. M. Huguet, and F. Ritort, Biophys. J. 100,1765 (2011). [39] A. F. Oberhauser, P. K. Hansma, M. Carrion-Vazquez, and J. M. Fernandez, Proc. Natl. Acad. Sci. USA 98,468 (2001). [40] J. M. Fernandez and H. Li, Science 303,1674 (2004). [41] G. Hummer and A. Szabo, Biophys. J. 85,5(2003). [42] The length Lis extensive in the sense that it is proportional to the number of units Nin the limit N1. [43] C. J. Thompson, Classical Equilibrium Statistical Mechanics (Oxford University Press, Oxford, 1988). [44] V. Y. Chernyak, M. Chertkov, and C. Jarzynski, J. Stat. Mech: Theor. Exp. (2006)P08001. [45] L. D. Landau and E. M. Lifshitz, Statistical Physics Part 1 (Course of Theoretical Physics) (Pergamon Press, Oxford, 1980), Vol 5. [46] In biomolecules at zero force, the folded state is much more localized than the unfolded state but this feature cannot be reproduced with the simple Landau-like potential we are using. Due to its simplicity, the barrier separating the two minima is always closer to the metastable state. In order to make the “width” of the folded state smaller than that of the unfolded state, more realistic potentials like the one in Refs. [17,18,28] must be used. See also Appendix A. [47] R. Kapri, Phys. Rev. E 86,041906 (2012). [48] This time interval is much longer than the time step used to integrate the Langevin equations, but much shorter than the total time for the stretching (or relaxing) process. As a result, the plotted force is close to Fexp,Eq.(8b), because the time average of F is very small. [49] A. F. Oberhauser, P. E. Marszalek, M. Carrion-Vazquez, and Julio M. Fernandez, Nat. Struct. Biol. 6,1025 (1999). [50] M. Rief, J. Pascual, M. Saraste, and H. E. Gaub, J. Mol. Biol. 286,553 (1999). [51] I. Schwaiger, C. Sattler, D. R. Hostetter, and M. Rief, Nat. Mater. 1,232 (2002). [52] W. Lee, X. Zeng, H.-X. Zhou, V. Bennet, W. Yang, and P. E. Marszalek, J. Biol. Chem. 285,38167 (2010). [53] F. Latinwo, K.-W. Hsiao, and C. M. Schroeder, J. Chem. Phys. 141,174903 (2014). [54] In Ref. [21], it was LJ−1+LJinstead of 2LJin the numerator of rJ. Equation (29) is more consistent, since LJ−1−LJis of the order of N−1. [55] A. Carpio and L. L. Bonilla, Phys.Rev.Lett.86,6034 (2001). [56] A. Carpio, L. L. Bonilla, and G. Dell’Acqua, Phys. Rev. E 64, 036204 (2001). [57] A. Carpio and L. L. Bonilla, SIAM J. Appl. Math. 63,1056 (2003). [58] N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1997). [59] D. Bensimon, D. Dohmi, and M. M´ ezard, Europhys. Lett. 42, 97 (1998). [60] G. Giacomin and F. L. Toninelli, Phys. Rev. Lett. 96,070602 (2006). [61] S. Ares and A. S´ anchez, Eur. Phys. J. B 56,253 (2007). [62] R. Merkel, P. Nassoy, A. Leung, K. Ritchie, and E. Evans, Nature (London) 397,6714 (1999). [63] J. Bruji´ c, R. I. Hermans, K. A. Walther, and J. M. Fernandez, Nat. Phys. 2,282 (2006). [64] J. J. Brey and A. Prados, Phys. Rev. E 47,1541 (1993). [65] J. J. Brey and A. Prados, Phys.Rev.B49,984 (1994); J. J. Brey, A. Prados, and M. J. Ruiz-Montero, J. Non-Cryst. Solids 172-174,371 (1994). [66] A. Prados, J. J. Brey, and B. S´ anchez-Rey, Phys. Rev. B 55, 6343 (1997); A. Prados and J. J. Brey, Phys. Rev. E 64,041505 (2001). [67] F. Ritort and P. Sollich, Adv. Phys. 52,219 (2003). [68] A. S. Keysa, J. P. Garrahan, and D. Chandler, Proc. Natl. Acad. Sci. USA 110,4482 (2013). 052712-18 THEORY OF FORCE-EXTENSION CURVES FOR MODULAR . . . PHYSICAL REVIEW E 91, 052712 (2015) [69] P. J. Elms, J. D. Chodera, C. Bustamante, and S. Marqusee, Proc. Natl. Acad. Sci. USA 109,3796 (2012). [70] H. Bai, J. E. Kath, F. M. Z¨ orgiebel, M. Sun, P. Ghosh, G. F. Hatfull, N. D. F. Grindley, and J. F. Marko, Proc. Natl. Acad. Sci. USA 109,16546 (2012). [71] D. E. Makharov, Biophys. J. 96,2160 (2009). [72] The boundary conditions (32) imply that the terms corresponding to j=0andj=Ndo not contribute to the interaction between nearest-neighbor spins. 052712-19