Improving free-energy estimates from unidirectional work measurements: theory and experiment
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).
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.
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.
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.
supplemental/10.1103/PhysRevLett.107.060601 for the appendices.