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 N 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 (log N)/m, where m measures the width of the lower tail of the work distribution, and then compute the finite-N corrections to this limit with different approaches for different regimes of (log N)/m. 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 [1] for review). The Crooks fluctuation theorem (CFT) [2] states that the probability distribution p(W ) of the work W done on a system driven out of equilibrium following an arbitrary finite-time protocol obeys the relation p(W )/p R (−W ) = e (W −∆F )/kB T . Here, p R (W ) is the work distribution (WD) for the corresponding time-reversed protocol, ∆F is the free-energy difference between the final and initial equilibrium states [3], and T is the temperature. Hence, ∆F can be estimated in bidirectional experiments by repeating many times the forward and reverse protocol, as demonstrated using single-molecule manipulation techniques [4,5]. An asymptotically unbiased estimator based on the CFT is the acceptance ratio (AR) estimator [6].
In many experimental settings, which we shall call unidirectional, the reverse work cannot be measured. Examples are found in AFM pulling of biopolymers [7,8], steered simulations [9], free-energy landscape reconstruction [10], and single-molecule experiments on protein unbinding, intercalation, specific cation binding, antigenantibody interactions, and non-native protein conformations. In these cases, an alternative method is provided by a corollary of the CFT, the Jarzynski equality (JE) e −W/kBT = e −∆F/kBT , where · is the expectation over p(W ) [3]. Given N work measurements W 1 , . . . , W N under the same protocol, the Jarzynski estimator converges to ∆F from above as N → ∞ (here and henceforth we set k B T = 1 and express all work values in units of k B T at room temperature). In practice, convergence of ∆F N requires that rare trajectories with W i < ∆F be sufficiently represented, which in turn requires N ≫ exp(D typ ), where D typ is the typical value of the dissipated work, D = W − ∆F [11]. Therefore, ∆F N is a reliable estimator of ∆F only when D typ is not much larger than k B T . It is thus important to have a quantitative estimate of the bias B N = ∆F N − ∆F . 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 [12], even for a Gaussian WD (GWD). Expansions in N −1 [13,14] are only applicable when the bias is of order N −1 , i.e. smaller than the O(N − 1 2 ) statistical error and thus negligible. In the relevant regime B N ≫ O(N −1 ), power-law interpolations in N [14] and other approximations [12] have been discussed, but no reliable analytical theory exists.
In this Letter, we derive analytical expressions for the bias expectation B N 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 [16]. The construct is tethered to two beads, one held by a pipette, the other by an optical trap created by counterpropagating laser beams [17]. 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 (≈ 200 − 1000 times per experiment) we collect the histogram of the WD's p U,R (W ) 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 [18] for details). We divide the data in blocks of N cycles, compute ∆F N  for each block, and average over the blocks to estimate ∆F N for U and R separately. As shown in Fig.1(d,e) for hairpins A and C, ∆F N tends to the AR estimate ∆F AR for large N , 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 W ≪ W c , where W c is a characteristic work value, Ω > 0 measures the tail width and q is a normalization constant. For the JE to hold, generally one must have δ > 1 [19]. Two key parameters in the following are A saddle point calculation gives Hence the JE implies µ → D c in this limit. An example of a WD obeying Eq.(2) is a GWD with mean W = W c and variance σ 2 W = Ω 2 /2 (i.e. α = 0, δ = 2, q = π −1/2 ), for which e −kD = exp(µk 2 − k D ) and thus µ = D = D c = σ 2 W /2 exactly for all Ω. This relation allows one to define another unidirectional [20], since ∆F v = ∆F for the GWD.
Scaling limit -Our strategy consists in computing first B N in a suitable scaling limit, and then the finite-N corrections to this limit.
We obtain the scaling limit by mapping the problem onto the Random Energy Model (REM) [21] as [21,22]. From the known limit of N −1 log Z N (β) for N → ∞ [22,23] we obtain B N → B REM in the scaling limit (N, Ω) → ∞ with λ finite. For λ > 1, all terms in Z N give a finite contribution and we find B REM = D c − µ. For λ < 1, corresponding to the glass phase of the REM [22], Z N is dominated by a finite number of terms and we obtain (4) 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-N corrections is apparent.
Finite-N corrections -One must distinguish three regimes, which require different analytical approaches: λ > 1, λ ≪ 1, and λ 1. For λ > 1, by partially resumming the 1/N expansion we obtain a closed expression for B N that improves considerably over the truncated expansions previously considered [13], which are valid only for λ ≫ 1 (see Sec. 2.1 in [18]). However, the most relevant regime in applications of the JE is λ < 1, since in practice one usually has N ≪ exp(D typ ). In the limit λ ≪ 1, using an extreme-value approach (Sec. 2.2 in [18]) we obtain to leading order where γ E is the Euler-Mascheroni constant. Cook and Derrida [24] were able to compute the finite-N corrections in the critical region λ 1 for the GWD case in the context of the REM. We have extended their travelingwave approach to the more general case of Eq.(2). In this way (see Sec. 2.3 in [18] for details) we recover Eq.(5) for λ ≪ 1, while for λ 1 we obtain where erfc is the complementary error function and θ = (λ (5) and (6) provide, via Eqs.(3) and (4), explicit expressions for B N as a function of D c , log N , and the shape parameters δ, Ω, α, q of the WD tail. To illustrate the validity of these expressions, we computed numerically B N by sampling W i from the GWD and from the Weibull WD (WWD): where W c is fixed numerically by imposing the JE. The WWD satisfies Eq.(2) (α = 1 − δ, q = δ) and allows us to model tails falling faster (δ > 2) or more slowly (δ < 2) than a GWD. In this In their respective range of validity in λ, Eqs.(5) and (6) agree very well with the numerical data for the entire range tested (1.1 ≤ δ ≤ 3, 1.41 ≤ Ω ≤ 16) also for large bias, as shown for some cases in Fig. 3(a) (for more examples see Fig.S4 in [18]). Furthermore, substituting D c with µ in Eq.(4) worsens only slightly the agreement, an important observation for the following (see Sec. 2.5 in [18]). In the special case of the GWD, Eq. (6) gives B N = log 2 at λ = 1. The empirical one-parameter power law (7) interpolating between λ = 1 and λ = 0 (for which B N = D ) 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 [18]). A power-law fit of the bias in N was proposed in Ref. [14]. 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 D c = µ suggests an improved estimator ∆F * N applicable to any problem involving the logarithm of an exponential average: i) given N measurements W i , compute ∆F N and the histogram of the WD; ii) estimate the shape parameters Ω, δ, α, q (and thus µ) by fitting the histogram tail to Eq.(2), taking for instance the maximum of the WD for W c [15]; iii) define ∆F * N = ∆F N − B N taking for B N either Eq.(5) or (6) (depending on the value of λ) and setting D c = µ.
In the special case of the Jarzynski estimator, we must take into account the stronger constraint imposed by the CFT, which implies that e −W is dominated by values of W near the maximum W † of the reverse WD [25]. Sampling these very rare events is usually unfeasible, so there is no guarantee that the fit to the measured p(W ) will continue to hold near W † . 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 δ ≈ 2), it is reasonable to assume that the fit will hold near W † , hence, ∆F * N should give a good estimate. The distance from a GWD can be self-consistently quantified a posteriori by the ratio We return now to our DNA experiments, for which we can compare the unidirectional estimators ∆F N , ∆F v , ∆F * N separately for U and R with the bidirectional estimator ∆F AR . Figure 1(d) shows a case with mild dissipation, for which the bias of ∆F N is small and all estimators converge. Figure 1(e) shows an experiment with intermediate dissipation. In this case sampling the tail of p U near the maximum of p R (or viceversa) would require an unfeasible number of cycles. Nevertheless, both tails are well fitted with δ = 2, α = 0.
[They are not perfect GWD, being slightly asymmetric. Using D = W − ∆F AR we obtain r = 0.67 (U), 0.81 (R).] The curves in Fig.1(e) represent ∆F AR + B N with B N given by Eq.(5), (6), (7). We find that ∆F * N agrees with ∆F AR within its statistical error for N ≥ 5, for all three expressions. This represents a significant improvement over the variance estimator, which has a bias ∆F v − ∆F AR = 4.3 ± 0.8 (U), −2.9 ± 0.8 (R), and over the uncorrected Jarzynski estimator. For instance, we have ∆F N − ∆F AR = 5.8 ± 3.0 (U), −6.6 ± 2.0 (R) for N = 36, and ∆F N − ∆F AR = 5.1 ± 0.6 (R) for N = 289. Furthermore, the fitted WD satisfy fairly well the CFT (see Sec. 4 in [18]), 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 ∆F N converges very slowly [ (Fig.3(c)]. Equation (2) fits reasonably well the WD tails for a fairly broad range of δ values (we take α = 0 for simplicity [15]). The estimator ∆F * N 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 [18] for other examples). Note also that Eqs. (5,6) are ill-defined as the exponent δ approaches 1 [19]. As δ decreases the fitted tails satisfy less and less the CFT (see Sec. 4 in [18]), signaling that δ must increase further out in the tails. It is an open problem to generalize our analytical approach to an effective δ varying with W .
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.