Improving free-energy estimates from unidirectional work measurements: theory and experiment

Improving free-energy estimates from unidirectional work measurements: theory and experiment

Matteo Palassini    Felix Ritort Departament de Física Fonamental, Universitat de Barcelona, Diagonal 647, E–08028 Barcelona, Spain
CIBER-BBN de Bioingeniería, Biomateriales y Nanomedicina, Instituto de Salud Carlos III, Madrid, Spain
July 9, 2019

We derive analytical expressions for the bias of the Jarzynski free-energy estimator from nonequilibrium work measurements, for a generic work distribution. To achieve this, we map the estimator onto the Random Energy Model in a suitable scaling limit parametrized by , where measures the width of the lower tail of the work distribution, and then compute the finite- corrections to this limit with different approaches for different regimes of . We show that these expressions describe accurately the bias for a wide class of work distributions, and exploit them to build an improved free-energy estimator from unidirectional work measurements. We apply the method to optical tweezers unfolding/refolding experiments on DNA hairpins of varying loop size and dissipation, displaying both near-Gaussian and non-Gaussian work distributions.


The accurate measurement of free-energy changes has important applications in physics, chemistry, and biology. Traditional measurement methods rely on reversible, near-equilibrium transformations, which however are often unfeasible. In recent years, new results in nonequilibrium statistical mechanics have suggested ways to measure free-energy changes from experiments (and simulations) far from equilibrium (see reviewJ () for review). The Crooks fluctuation theorem (CFT) crooks () states that the probability distribution of the work done on a system driven out of equilibrium following an arbitrary finite-time protocol obeys the relation . Here, is the work distribution (WD) for the corresponding time-reversed protocol, is the free-energy difference between the final and initial equilibrium states JE (), and is the temperature. Hence, can be estimated in bidirectional experiments by repeating many times the forward and reverse protocol, as demonstrated using single-molecule manipulation techniques bustamante (); collin (). An asymptotically unbiased estimator based on the CFT is the acceptance ratio (AR) estimator shirts2004 ().

In many experimental settings, which we shall call unidirectional, the reverse work cannot be measured. Examples are found in AFM pulling of biopolymers hummer-szabo (); kiang (), steered simulations schulten (), free-energy landscape reconstruction woodside (), and single-molecule experiments on protein unbinding, intercalation, specific cation binding, antigen-antibody interactions, and non-native protein conformations. In these cases, an alternative method is provided by a corollary of the CFT, the Jarzynski equality (JE) , where is the expectation over JE (). Given work measurements under the same protocol, the Jarzynski estimator


converges to from above as (here and henceforth we set and express all work values in units of at room temperature). In practice, convergence of requires that rare trajectories with be sufficiently represented, which in turn requires , where is the typical value of the dissipated work, jarzynski2006 (). Therefore, is a reliable estimator of only when is not much larger than . It is thus important to have a quantitative estimate of the bias . The mathematical problem faced is that of calculating the distribution of a (log)sum of exponentials of i.i.d. random variables, Eq.(1), which depends on the system- and protocol-specific WD. No closed solution to this problem is available romeo (), even for a Gaussian WD (GWD). Expansions in ZW (); gore () are only applicable when the bias is of order , i.e. smaller than the statistical error and thus negligible. In the relevant regime , power-law interpolations in gore () and other approximations romeo () have been discussed, but no reliable analytical theory exists.

In this Letter, we derive analytical expressions for the bias expectation for a wide class of WD’s and validate them by comparison with exact numerical simulations, also in the regime of large bias. We use these results to build an improved unidirectional free-energy estimator by correcting for the bias of Eq.(1). We then discuss unfolding/refolding experiments on DNA hairpins, which allow us to test our method against the bidirectional AR estimator.

The experimental setup is shown in Fig. 1(a). We synthesized five hairpins (A,B,C,D,E) with identical stem and (GAAA…) loops of 4,6,12,16,20 bases, respectively [Fig.1b]. The hairpins are inserted between two short (29bp) dsDNA handles to improve signal-to-noise resolution ForLorManHayHugRit11 (). The construct is tethered to two beads, one held by a pipette, the other by an optical trap created by counterpropagating laser beams Huguet10 (). The light deflected by the trapped bead provides a direct measurement of the force acting on the molecule. By moving the trap away from the pipette at constant velocity, the hairpin is stretched until it unfolds. Subsequent reversal of the velocity causes the hairpin to refold. By repeating this cycle ( times per experiment) we collect the histogram of the WD’s for the work to unfold (U) and refold (R) the hairpin, measured by integrating the force-distance curves (Fig.1c) for the forward and reverse part of each cycle (see Sec. 1 in Supp () for details).

Figure 1: Bias measurements in DNA hairpins. a) Experimental setup. b) Hairpin B sequence. c) Examples of force-distance cycles for hairpins A-E. d) Left: histograms of (upper), (lower) for hairpin A pulled at 400 nm/s. The horizontal line is the AR estimate , giving =4.3 (U), 3.9 (R). The lines are GWDs fitting the lower (upper) tail of (). Right: Jarzynski estimator as a function of for U and R. Errors are estimated by jackknife. The lines represent with given by Eq.(6) (dashed line) and Eq.(7) (continuous line) for the GWD case, assuming , where is estimated from the GWD fit to the tails. e) Same as d) but for hairpin C pulled at 65 nm/s [=13.6 (U); 14.8 (R)]. Also shown is with given by Eq.(5) (dotted line).

We divide the data in blocks of cycles, compute for each block, and average over the blocks to estimate for U and R separately. As shown in Fig.1(d,e) for hairpins A and C, tends to the AR estimate for large , from opposite sides for U and R. Note that the dissipation increases with loop size and pulling speed.

We analyze theoretically the bias for a generic WD with finite mean and an unbounded lower tail which decays as


for , where is a characteristic work value, measures the tail width and is a normalization constant. For the JE to hold, generally one must have tp (). Two key parameters in the following are


A saddle point calculation gives for large , where . Hence the JE implies in this limit. An example of a WD obeying Eq.(2) is a GWD with mean and variance (i.e. , , ), for which and thus exactly for all . This relation allows one to define another unidirectional estimator wood (), since for the GWD.

Scaling limit – Our strategy consists in computing first in a suitable scaling limit, and then the finite- corrections to this limit. We obtain the scaling limit by mapping the problem onto the Random Energy Model (REM) derrida1 () as , where is the REM partition function and the i.i.d. variables have a distribution decaying as for derrida1 (); MB (). From the known limit of for MB (); benarous () we obtain in the scaling limit with finite. For , all terms in give a finite contribution and we find . For , corresponding to the glass phase of the REM MB (), is dominated by a finite number of terms and we obtain


Figure 2a shows the approach to the scaling limit as increases, for the GWD case. Significant deviations occur for moderate , from which the need to compute finite- corrections is apparent.

Figure 2: Finite- corrections to the bias for a GWD. a) Convergence of to its scaling limit. The points joined by lines are averages of over many sets () sampled from a GWD with variance and mean . The dashed line represents [(Eq.(4) for ]. Inset of a): The unscaled data. The continuous lines represent Eq.(7). The horizontal dashed line indicates the accuracy limit common in biophysical studies. b) Bias (estimated as ) for all experiments on hairpins A,B,C and , including both U and R. The pulling speed is in the range 25 - 300 nm/s and is in the range 1 - 20 . The lines show Eq.(6) for the GWD.

Finite- corrections – One must distinguish three regimes, which require different analytical approaches: , , and . For , by partially resumming the expansion we obtain a closed expression for that improves considerably over the truncated expansions previously considered ZW (), which are valid only for (see Sec. 2.1 in Supp ()). However, the most relevant regime in applications of the JE is , since in practice one usually has . In the limit , using an extreme-value approach (Sec. 2.2 in Supp ()) we obtain to leading order


where is the Euler-Mascheroni constant. Cook and Derrida cook () were able to compute the finite- corrections in the critical region for the GWD case in the context of the REM. We have extended their traveling-wave approach to the more general case of Eq.(2). In this way (see Sec. 2.3 in Supp () for details) we recover Eq.(5) for , while for we obtain


where erfc is the complementary error function and .

Numerical test – Equations (5) and (6) provide, via Eqs.(3) and (4), explicit expressions for as a function of , and the shape parameters of the WD tail. To illustrate the validity of these expressions, we computed numerically by sampling from the GWD and from the Weibull WD (WWD): for ; for , where is fixed numerically by imposing the JE. The WWD satisfies Eq.(2) (, ) and allows us to model tails falling faster () or more slowly () than a GWD. In this case , .

In their respective range of validity in , Eqs.(5) and (6) agree very well with the numerical data for the entire range tested (, ) also for large bias, as shown for some cases in Fig. 3(a) (for more examples see Fig.S4 in Supp ()). Furthermore, substituting with in Eq.(4) worsens only slightly the agreement, an important observation for the following (see Sec. 2.5 in Supp ()). In the special case of the GWD, Eq.(6) gives at . The empirical one-parameter power law


interpolating between and (for which ) also fits fairly well the GWD data, as shown in Fig.2, although less so than Eq.(6) for large (see also Fig.S4 in Supp ()). A power-law fit of the bias in was proposed in Ref.gore (). Figure 2(b) shows that the bias of all our experiments with mild dissipation is well described by Eq.(6) for a GWD.

Figure 3: free-energy recovery for non-Gaussian WDs. a) Analytical estimates of as given in Eq.(5) (dashed lines) and Eq.(6) (dotted lines) with replaced by , tested against simulated data for the WWD (symbols). From top to bottom, = (1.1, 1.83, 31.0), (1.1, 1.68, 13.5), (1.5, 6, 30.1), (1.5, 4, 8.8). b) WD for hairpin E pulled at 130 nm/s and fitted WD tails assuming (dashed lines). The vertical line represents . c) for U (circles) and R (triangles). Note the very slow convergence. The improved estimator obtained by correcting for the bias with Eq.(5) (squares) and Eq.(6) (stars) using converges quickly to a value consistent with (dashed horizontal line) d) recovered from U and R as a function of for . Using Eq.(5) (dashed lines) or Eq.(6) (continuous lines) gives nearly the same result.

Free-energy recovery method – The fact that Eqs.(5) and (6) describe well the bias when suggests an improved estimator applicable to any problem involving the logarithm of an exponential average: i) given measurements , compute and the histogram of the WD; ii) estimate the shape parameters (and thus ) by fitting the histogram tail to Eq.(2), taking for instance the maximum of the WD for on_alpha (); iii) define taking for either Eq.(5) or (6) (depending on the value of ) and setting .

In the special case of the Jarzynski estimator, we must take into account the stronger constraint imposed by the CFT, which implies that is dominated by values of near the maximum of the reverse WD Ritort04 (). Sampling these very rare events is usually unfeasible, so there is no guarantee that the fit to the measured will continue to hold near . Nevertheless, we argue that a distinction should be made between near-Gaussian and non-Gaussian tails. In the former case (i.e. when the tail can be fitted with ), it is reasonable to assume that the fit will hold near , hence, should give a good estimate. The distance from a GWD can be self-consistently quantified a posteriori by the ratio ( for a GWD), taking .

We return now to our DNA experiments, for which we can compare the unidirectional estimators separately for U and R with the bidirectional estimator . Figure 1(d) shows a case with mild dissipation, for which the bias of is small and all estimators converge. Figure 1(e) shows an experiment with intermediate dissipation. In this case sampling the tail of near the maximum of (or viceversa) would require an unfeasible number of cycles. Nevertheless, both tails are well fitted with . [They are not perfect GWD, being slightly asymmetric. Using we obtain 0.67 (U), 0.81 (R).] The curves in Fig.1(e) represent with given by Eq.(5), (6), (7). We find that agrees with within its statistical error for , for all three expressions. This represents a significant improvement over the variance estimator, which has a bias (U), (R), and over the uncorrected Jarzynski estimator. For instance, we have (U), (R) for , and (R) for . Furthermore, the fitted WD satisfy fairly well the CFT (see Sec. 4 in Supp ()), giving further confidence in the consistency of the method.

Finally, we consider an experiment where the tails are far from a GWD and the dissipation is large (Fig.3). In this case the WDs are wide apart [Fig.3(b)] and converges very slowly [(Fig.3(c)]. Equation (2) fits reasonably well the WD tails for a fairly broad range of values (we take for simplicity on_alpha ()). The estimator is shown in Fig.3d: the pronounced dependence on and the discrepancy between U and R confirm that the predictive power of Eqs.(5,6) relies on accurately knowing (see Sec. 3 in Supp () for other examples). Note also that Eqs.(5,6) are ill-defined as the exponent approaches 1 tp (). As decreases the fitted tails satisfy less and less the CFT (see Sec. 4 in Supp ()), signaling that must increase further out in the tails. It is an open problem to generalize our analytical approach to an effective varying with .

In summary, we obtained analytical expressions for the bias of the Jarzynski estimator, and showed that they can be used to obtain improved unidirectional estimates of the free energy of mechanical unfolding of DNA hairpins, provided the WD tail is described by a compressed exponential over a wide enough range of work values. These results are applicable to many unidirectional experiments and simulations, and are relevant to other contexts involving sums of random exponentials.

We thank N. Skantzos for discussions in the initial stages of this work. MP thanks NORDITA for hospitality. Work supported by MICINN (FIS2006-13321-C02-01, FIS2007-3454), HFSP Grant No. RGP55-2008, and ICREA Academia grants.



  • (1) C. Jarzynski, Eur. Phys. J. B 64, 331 (2008).
  • (2) G.E. Crooks, Phys. Rev. E 60, 2721 (1999).
  • (3) C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
  • (4) J. Liphardt et al., Science 296, 1832 (2002).
  • (5) D. Collin et al., Nature 437, 198 (2005).
  • (6) M. R. Shirts, E. Bair, G. Hooker and V. S. Pande, Phys. Rev. Lett. 91, 140601 (2003).
  • (7) G. Hummer and A. Szabo, PNAS 98, 3658 (2001).
  • (8) N. C. Harris and C.-H. Kiang, Phys. Rev. E 79, 041912 (2009).
  • (9) M. O. Jensen, S. Park, E. Tajkhorshid and K. Schulten, PNAS 99 6731 (2002).
  • (10) A.N. Gupta et al.,, Nat. Phys. (2011) doi:10.1038/nphys2022.
  • (11) C. Jarzynski, Phys. Rev. E 73, 046105 (2006).
  • (12) M. Romeo, V. Da Costa, and F. Bardou, Eur. Phys. J. B 32, 513 (2003).
  • (13) D. M. Zuckerman and T. B. Woolf, Phys. Rev. Lett. 89, 180602 (2002); J. Stat. Phys. 114, 1303 (2004).
  • (14) J. Gore, F. Ritort, and C. Bustamante, PNAS 100 12564 (2003).
  • (15) Since is slowly varying compared to the exponential, in fitting Eq.(2) to experiments one can safely set .
  • (16) N. Forns et al., Biophys. J. 100, 1765 (2011).
  • (17) J. M. Huguet et al., PNAS 107, 15431 (2010).
  • (18) See Supplemental Material at
    supplemental/10.1103/PhysRevLett.107.060601 for the appendices.
  • (19) The JE can also be satisfied with . This marginal case requires a separate analysis.
  • (20) R.H. Wood, W.F.C. Muhlbauer, and P.T. Thompson, J. Phys. Chem. 95, 6670 (1991).
  • (21) B. Derrida, Phys. Rev. B 24, 2613 (1981).
  • (22) J.-P. Bouchaud and M. Mézard, J. Phys. A 30 7997 (1997).
  • (23) The scaling limit of the REM has been studied rigorously in G. Ben Arous, L.V. Bogachev, S.A. Molchanov, Prob. Th. Rel. Fields, 132, 579 (2005); A. Bovier, I. Kurkova, and M. Löwe, Ann. Probab., 30, 605 (2002).
  • (24) J. Cook and B. Derrida, J. Stat. Phys. 63, 505 (1991).
  • (25) F. Ritort, J. Stat. Mech. P10016 (2004).
  • Comments 0
    Request Comment
    You are adding the first comment!
    How to quickly get a good reply:
    • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
    • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
    • Your comment should inspire ideas to flow and help the author improves the paper.

    The better we are at sharing our knowledge with each other, the faster we move forward.
    The feedback must be of minimum 40 characters and the title a minimum of 5 characters
    Add comment
    Loading ...
    This is a comment super asjknd jkasnjk adsnkj
    The feedback must be of minumum 40 characters
    The feedback must be of minumum 40 characters

    You are asking your first question!
    How to quickly get a good answer:
    • Keep your question short and to the point
    • Check for grammar or spelling errors.
    • Phrase it like a question
    Test description