DensityDependent Analysis of Nonequilibrium Paths Improves Free Energy Estimates
Abstract
When a system is driven out of equilibrium by a timedependent protocol that modifies the Hamiltonian, it follows a nonequilibrium path. Samples of these paths can be used in nonequilibrium work theorems to estimate equilibrium quantities, such as free energy differences. Here, we consider analyzing paths generated with one protocol using another one. It is posited that analysis protocols which minimize the lag, the difference between the nonequilibrium and the instantaneous equilibrium densities, will reduce the dissipation of reprocessed trajectories and lead to better free energy estimates. Indeed, when minimal lag analysis protocols based on exactly soluble propagators or relative entropies are applied to several test cases, substantial gains in the accuracy and precision of estimated free energy differences are observed.
I Introduction
The accurate and efficient estimation of free energy differences is an important goal in chemical physics and remains an active area of research. One promising approach to free energy estimation entails measuring the work done on a system over repetitions of an irreversible process. According to the second law of thermodynamics, the mean work is greater than the free energy difference between the end states of the process, . Nonequilibrium work theorems Jarzynski (1997a, b); Crooks (1998, 1999, 2000) supplement this upper bound by rigorously equating with other averaged functions of the work. These theorems have been empirically validated in singlemolecule pulling experiments Liphardt et al. (2002); Collin et al. (2005) and computer simulations (e.g. Ref. Hummer (2001)).
Jarzynski’s equality, Jarzynski (1997a, b) a unidirectional nonequilibrium work theorem, relates the free energy difference to an exponential average of the work. Unfortunately, because it uses a nonlinear (specifically, a logarithmic) function of the average, the free energy estimator based on this equality suffers from a systematic finitesampling bias. Zuckerman and Woolf (2002); Gore et al. (2003); Zuckerman and Woolf (2004) While accurate in the limit of infinite sampling, this estimator is usually dominated by rare events where the work is less than the free energy difference, and thereby converges slowly. Jarzynski (2006)
If the average amount of work dissipated as heat is reduced, these lowwork events will be more frequent and accurate free energy estimation will usually require fewer work samples. The most straightforward way to reduce heat dissipation is to slow the rate of the process; in the limit of infinitely slow switching, the process is reversible and the work is equal to the free energy difference. Unfortunately, reducing the switching rate requires additional time and lowers the signaltonoise ratio in singlemolecule pulling experiments. Maragakis et al. (2008) Under the constraint of constant experiment length, it is possible to reduce heat dissipation by optimizing the switching protocol that controls how the thermodynamic state changes with time. Protocol variation predates Jarzynski’s equality, having been applied to tightening free energy bounds from the second law of thermodynamics. Mark et al. (1990); Reinhardt and Hunter (1992); Hunter et al. (1993); Schon (1996); Jarque and Tidor (1997) More recently, variational calculus has been applied to find optimal protocols that minimize the mean work. Schmiedl and Seifert (2007); Then and Engel (2008); GomezMarin et al. (2008)
While protocol variation is, in principle, feasible in laboratory experiments, many more approaches to improving nonequilibriumbased free energy estimation are possible in computer simulations. For example, Wu and Kofke were inspired by the Rosenbluth chain sampling scheme to develop methods for generating lowwork nonequilibrium paths. Wu and Kofke (2005) Vaikuntanathan and Jarzynski took another approach, altering the system dynamics, to reduce heat dissipation and improve free energy estimates. Vaikuntanathan and Jarzynski (2008) The approach most mathematically similar to this work, however, is importance sampling in nonequilibrium path space. Sun (2003); Atilgan and Sun (2004); Ytreberg and Zuckerman (2004); Oberhofer et al. (2005); Oberhofer and Dellago (2008)
In importance sampling, samples from one distribution are used to estimate expectations in another. The technique is often applied to Markov chain Monte Carlo and molecular dynamics simulations (where it is usually called umbrella sampling): after applying a configurational bias to overcome energy barriers and promote ergodicity, expectations are calculated for the unbiased ensemble. Importance sampling has been extended to transition path sampling Pratt (1986); Dellago et al. (1998) with nonequilibrium trajectories. In this algorithm, a biasing function modifies the Monte Carlo acceptance criteria of proposed paths in a way that improves the convergence of free energy estimates. Sun (2003); Atilgan and Sun (2004); Ytreberg and Zuckerman (2004); Oberhofer et al. (2005); Oberhofer and Dellago (2008)
Here, we apply the importance sampling formalism in a completely different way. Instead of sampling nonequilibrium trajectories in a biased manner, we focus on the analysis of previously generated paths. Instead of asking which pathensemble we would like to sample from, we ask which pathensemble average we would like to evaluate. This is accomplished by processing paths generated using one protocol  the sampling protocol,  using another  the analysis protocol, .
While we have infinite freedom in selecting an analysis protocol, not all choices will improve the convergence of free energy estimates. One reasonable strategy for choosing is to minimize the lag, the difference between the nonequilibrium and instantaneous equilibrium densities; Vaikuntanathan and Jarzynski found that under certain dynamics, dissipation is eliminated if there is no lag, leading to a zerovariance estimator of . Vaikuntanathan and Jarzynski (2008) To reduce the lag, they modified their equation of motion with an additional flowfield term that “escorts” the system along a nearequilibrium path. Essentially, this strategy modifies the nonequilibrium density. In this paper, we take the opposite approach: using the analysis protocol to choose an instantaneous equilibrium density that closely matches the sampled nonequilibrium density.
As an illustrative case, consider a Brownian particle in a harmonic oscillator, or spring, which moves at a constant velocity (Fig. 1). If the system starts in thermal equilibrium, its density is a Gaussian about the initial spring position. When the spring starts moving, the density remains a Gaussian with the same variance, but its mean position, , lags behind the spring position. Mazonka and Jarzynski (1999); Minh and Adib (2009) For this particular system, an analysis protocol based on will have no lag. We shall further explore this system in Section III.
One complication with using an analysis protocol that minimizes the lag is that its end state is usually not the same as in the sampling protocol. Thus, the free energy difference being estimated differs. To estimate the same with a minimal lag analysis protocol, it may be necessary to extend or otherwise modify the sampling. To distinguish the two situations, we shall refer to the former as protocol postprocessing and the latter as nonequilibrium densitydependent sampling (NEDDS). Both fall under the aegis of densitydependent analysis.
The structure of this paper is as follows: in Section II, the importance sampling form of Jarzynski’s equality is detailed; in Section III, densitydependent analysis is demonstrated on two cases in which the propagator is analytically known; in Section IV, a general method for finding minimal lag analysis protocols is described, applied to an adaptive algorithm for NEDDS, and tested on the model system; and lastly, implications of this method and possible future directions are discussed.
Ii Free Energy Formalism
Consider a system whose Hamiltonian, , depends on , its position in phase space (or configuration space), and a control parameter, . Initially, the system is prepared in thermal equilibrium at . The parameter is perturbed according to a protocol until it reaches a final state at . Jarznyski’s equality, Jarzynski (1997a, b)
(1) 
relates the free energy difference between the initial and final states of the protocol, , to an average over all possible paths, , resulting from this nonequilibrium procedure. Specifically, this expectation (denoted by the angled brackets ), is a path integral over infinitesimal elements with the protocoldependent density . During each process, the work done on the system is . (In this paper, all energies will be expressed in units of .)
Suppose that instead of , we consider an alternate density of paths, . The free energy difference can be calculated by applying a reweighed form of Jarzynski’s equality, Ytreberg and Zuckerman (2004)
(2) 
where is the ratio of probabilities of observing the trajectory given the densities. To analyze a finite sample of paths drawn from , we replace the expectations with sample mean estimators, obtaining, Ytreberg and Zuckerman (2004)
(3) 
where is the sample size. In a standard Jarzynski estimate, .
Previous workers have improved the convergence properties of Eq. (3) by choosing to be various workweighted functionals of the original density . Ytreberg and Zuckerman (2004); Oberhofer et al. (2005); Oberhofer and Dellago (2008) When introducing the singleensemble biased path sampling approach, Ytreberg and Zuckerman picked , such that . Ytreberg and Zuckerman (2004) In a paper comparing the method with conventional equilibrium procedures, Oberhofer et. al. considered . Oberhofer et al. (2005) By variation of the asymptotic variance with respect to the sampling bias, Oberhofer and Dellago found that optimal workweighted sampling is given by . Oberhofer and Dellago (2008) Unfortunately, this optimal choice is impractical because it includes the sought quantity .
In the present method, which applies Eq. (3) in a novel manner, depends on the sampling protocol and differs from unity when the analysis protocol varies from . Notably, the relevant work is , not , meaning that different choices of will result in various work distributions and convergence properties. This new way of applying importance sampling leads to different, albeit analogous, asymptotic variance expressions. Min () The present approach is more general than previous applications of Eq. (3), which require transition path sampling, because it does not require biased sampling and paths can be generated by ordinary dynamical equations. Indeed, under certain assumptions, such as those suggested by Nummela and Andricioaei, Nummela and Andricioaei (2007) it should be possible to apply the present method to laboratory experiments.
We note, as a caveat, that the importance sampling form of Jarzynski’s equality will only be useful for stochastic dynamics where can be computed. Under deterministic dynamics, is a delta function and having different sampling and analysis protocols will not improve free energy estimates.
Iii Cases with an Analytical Propagator
As mentioned earlier, we would like to choose an analysis protocol that minimizes the lag, such that the instantaneous equilibrium density corresponds with the sampled nonequilibrium density. This is particularly tractable when the propagator is exactly known. Here, we demonstrate Eq. (3) on two such cases: a Brownian particle in a harmonic oscillator (i) moving at a constant velocity or (ii) with a timedependent natural frequency. With both, the potential energy has the general form and the nonequilibrium density is
(4) 
where and are the most typical paths and spring coefficients, respectively. As these propagators can be obtained by close analogy to the path integral derivation of workweighted propagators, Minh and Adib (2009) their derivations are not detailed here. In case (i), is constant and moves the spring position according to , such that , , and
(5) 
In case (ii), is zero and controls the spring coefficient, , such that . In the corresponding nonequilibrium density, , and
(6) 
Based on these expressions, it is evident that, for case (i), the minimal lag analysis protocol is from Eq. (5), and for case (ii), it is from Eq. (6). In these special cases, the nonequilibrium density is exactly the equilibrium density corresponding to and there is no lag.
To test whether densitydependent analysis leads to improved free energy estimates, onedimensional Brownian dynamics simulations were run with the equation of motion,
(7) 
where is the position at time , is the time step, and is a standard normal random variable. The primes denote spatial derivatives such that and . For a discrete trajectory sampled with the protocol , where is the total number of steps, the probability ratio is , where and is the stochastic action (discretized from Ref. Minh and Adib (2009)),
(8)  
The work was evaluted with the discrete formula, . This action is valid in the continuum limit, as and . To approach this limit, we chose and a time step of .
The simulations were performed over steps (truncated to be an integer), where m refers to 7 evenly spaced numbers between 1.5 and 3. In case (i), was set to 25 and was chosen to start from and linearly progress to the target state at . With case (ii), is a linear interpolation between 1 and 100. Afterwards, the trajectories both analyzed with the standard Jarzynski estimate and subjected to protocol postprocessing with .
For comparison, NEDDS was implemented by switching at a faster rate such that the final state went beyond and the final nonequilibrium density, according to the propagators, corresponded to the target state. This is illustrated in Fig. 1d, where moving the harmonic oscillator at a faster rate than in Fig. 1c allows for the final density to correspond to the equilibrium state with . These trajectories, which took the same amount of simulation time for the same number of steps, were then reanalyzed with the appropriate .
In case (i), we find that protocol postprocessing with leads to a desirable result: most work values are reduced such that a larger fraction of them are less than the free energy difference (Fig. 2). Of these negative dissipation trajectories, most have a probability ratio less than one. Conversely, several positive dissipation trajectories have a probability ratio greater than one. For this set of trajectories, the modified work distribution leads to a more accurate free energy estimate.
Over a large number of repetitions and range of switching speeds, we find that free energy estimates based on are vastly improved over the standard procedure, , having significantly less variance and systematic bias (Fig. 3). The standard estimator only approaches the accuracy and precision of protocol processing for slow switches. Clearly, these trajectories are much better at estimating the end state free energy differences for than for . The estimates of from NEDDS also require considerably less sampling than the standard procedure, although the effect is somewhat less dramatic.
Similarly, in case (ii), densitydependent methods also show improvement over the standard Jarzynski estimate. For the timedependent natural frequency, the systematic bias of the standard estimate is relatively small but nonetheless evident at all sampled switching rates (Fig. 4). Estimates from both densitydependent methods have reduced bias and variance, and are found to be of similar quality to each other.
Iv General Case
In most practical situations, unfortunately, the propagator is not known ahead of time. Thus, prior to simulations, it is unclear how long paths need to be generated before the nonequilibrium density matches a density characteristic of the target state. While paths are being generated, however, it is possible to estimate the difference between the sampled density and arbitrary equilibrium states. States which minimize this difference can be be collected in an analysis protocol with minimal lag.
One measure of the distance between two probability distributions is the KullbackLeibler divergence, or the relative entropy. The relative entropy between the nonequilibrium density and an arbitrary equilibrium state is,
(9) 
When the integral is separated into two at the logarithm, one part is a constant with respect to . The divergence is minimized by finding a state where the other, , is the least. Using sampled discrete paths, this integral can be estimated by , where is the position at step of path . For a state , the equilibrium density is , where is the test state Hamiltonian and is its free energy. Thus, the relative entropy is minimized by the smallest value of,
(10) 
among different states . Generally, the free energy, , is unknown, but for states which occur along the switching protocol, (where is the free energy at ) can be estimated using the standard form of Jarzynski’s equality. These states constitute our search space for minimizing the lag.
Suppose we are interested in the free energy difference between the states defined by and . We can use to estimate on the fly and determine when to stop sampling via the following adaptive algorithm:

Start with and the work . For each of paths, obtain by drawing samples from the equilibrium ensemble at .

Propagate each path, calculating using a dynamical equation such as Eq. 7. To obtain , calculate the work done on the system during the time step and add it to . The next step in the sampling protocol, , is found by adding a predetermined value, , to . The sign of must be the same as . Increment by one.

Using values in the standard form of Jarzynski’s equality, estimate , the free energy difference between the states with and .

For each state corresponding to , , … , use and the free energy difference estimated in the previous algorithm step to calculate . The which minimizes is . Add to the minimal lag protocol .

If hasn’t crossed , repeat from algorithm step 2. Otherwise, set the final value in to .

Estimate the free energy difference using Eq. (3) with .
This algorithm was tested on Sun’s system, Sun (2003) where the potential energy is . Using Eq. (3), the free energy difference was estimated between the initial state with , where the potential is a single well, and the target state , where it is a double well, such that . Oberhofer et al. (2005) Brownian dynamics simulations were performed with the same diffusion coefficient, time step, and equation of motion as in Section III. The increment of at each time step was , where and refers to 9 evenly spaced values between 0 and 2. For comparison, the standard Jarzynski estimate was applied to simulations where is switched between 0 and 1 at a slower velocity, taking the same total time as in the corresponding NEDDS simulations.
In a representative set of simulations, the density most noticeably lags behind the sampling state at the beginning (Fig. 5). Around the state defined by , the lag quickly diminishes. However, the minima of does not reach the target state until the sampling is beyond 1.
Based on many repetitions of this procedure at different pulling speeds, we find that our NEDDS algorithm converges much more quickly than the standard Jarzynski estimate (Fig. 6). The systematic bias is largely eliminated with simulations that are switched nearly an order of magnitude faster. At the fastest switching rates, NEDDS remains biased but still outperforms the standard Jarzynski estimate. With these fast switchings, it is possible that the nonequilibrium density does not correspond well to any traversed equilibrium state.
V Discussion and Conclusion
With the goal of minimizing the lag via the choice of analysis protocol, we have developed densitydependent methods to analyze nonequilbrium paths, to estimate which states may constitute a protocol that minimizes the lag, and to adaptively sample paths until the desired density is achieved. Our promising results validate the strategy and provide further evidence for the link between lag and heat dissipation. They also hint that the accurate estimation of free energy differences may require adequate sampling in the important regions of both end states.
Analysis protocols provide another degree of freedom for lag reduction, and can be used in conjunction with other methods, such as sampling protocol optimization or biased path sampling. Furthermore, their use should extend beyond Jarzynski’s equality; they can potentially be applied in bidirectional nonequilibrium work expressions Min () or any relationship between a nonequilibrium process and a state function, such as Hummer and Szabo’s expression for the potential of mean force. Hummer and Szabo (2001) Quite possibly, our results are just the tip of an iceberg and this paper will open up new research directions for sampling and analyzing nonequilibrium trajectories.
Vi Acknowledgments
The author thanks Artur Adib, Christopher Jarzynski, Attila Szabo, and Suriyanarayanan Vaikuntanathan for pertinent discussions, and Gerhard Hummer for suggesting that he considers the lag. He also thanks Artur Adib for supporting a postdoctoral fellowship. This research was supported by the Intramural Research Program of the NIH, NIDDK.
In these appendices, we derive asymptotic variance and bias expressions for free energies estimated using protocol postprocessing. Our derivations are similar to that of Oberhofer et. al. for biased sampling of nonequilibrium trajectories. Oberhofer et al. (2005); Oberhofer and Dellago (2008) Although these expressions are not directly relevant to the main text, they support the line of inquiry explored in the paper and may be useful to those who wish to expand upon the results. Both unidirectional and bidirectional expressions are considered.
Appendix A Unidirectional
Expressed in importance sampling form, Jarzynski’s equality is,
(11) 
This is Eq. (3) in the main text, reproduced here for convenience.
Towards deriving the asymptotic expressions, we first define,
(12)  
(13) 
The sample mean estimators for the expectations of and are,
(14)  
(15) 
These form an estimator for the free energy,
(16) 
However, the sample mean estimators deviate from their true expectations,
(17)  
(18) 
Here, subscripts on the angle brackets are omitted for notational simplicity.
Deviations in these expectations lead to variance and bias in . The magnitude of the error can be estimated by a Taylor series expansion about , which becomes an increasingly reasonable approximation in the large sampling, or asymptotic, limit. To the first order, this expansion is,
(19)  
(20)  
(21) 
The partial derivatives are evaluated at their mean values.
The variance is defined as . Using the firstorder Taylor series expansion, this is,
(22)  
(23)  
(24)  
(25) 
The variance of a sample mean is the variance of the variable over the number of samples. Since and are functions of the same data points, the covariance of their sample means is, similarly, . Thus, the asymptotic variance is,
(26)  
(27)  
(28)  
(29) 
The bias is defined as . If we approximate this error with a firstorder Taylor series expansion, it is always zero. In order to obtain a nonzero bias expression, we use a secondorder Taylor series expansion about ,
(30)  
(31)  
(32) 
Using this approximation, the bias is found to be,
(33)  
(34)  
(35)  
(36) 
Notably, when the dissipated work is zero, , expressions for both the variance and bias are likewise zero.
Appendix B Bidirectional
Here, we consider the possibility of reanalyzing bidirectional data, collected using both a protocol and its time reversal . The results derived in this section suggest that for bidirectional data, the optimal analysis protocol is actually the sampling protocol. Thus, protocol postprocessing is less promising when applied to bidirectional data than to unidirectional data.
For notational consistency, we start our discussion with the Crooks Fluctuation Theorem,Crooks (1998, 1999)
(37) 
As in the main text, is the probability of observing trajectory , given the protocol . Analogously, is the probability of observing the time reversal, or conjugate twin, of , using the reverse protocol .
This theorem can be used to derive a relationship between forward and reverse pathensemble averages, Crooks (2000)
(38)  
(39)  
(40) 
In the above, is an arbitrary functional. As the last pathensemble average is over trajectories , trajectories sampled from must be reversed prior to being evaluated with the functional.
Next, we rearrange the pathensemble average theorem into an expression for the free energy. Crooks (2000)
(41) 
As in the unidirectional case, these pathensemble averages can be written as reweighed samples from other sampling densities.
(42)  
(43) 
The probability ratio is defined similarly to ,
(44) 
Using Eq. (37), we can show that it is related to by
(45) 
The ratio of can be shown to be unity by converting into a forward pathensemble average using Eqs. (40) and (45), and applying the importance sampling form of Jarzynski’s equality, Eq. (11).
The asymptotic variance of Eq. (43) can be calculated by a similar procedure to the unidirectional case. We start by defining,
(46)  
(47) 
Replacing and with and , we follow the same logic as in the unidirectional case from Eqs. (16) to (24). Next, we note that C and D are independent samples drawn from different ensembles and their correlation is zero. Thus, the variance estimate is,
(48)  
(49) 
We would like to combine the two terms including C and D in a single pathensemble average. In order to do that we need to convert the ensemble averages containing D to the forward direction. For , this is,
(50)  
(51)  
(52)  
(53)  
(54) 
For , this is
(55)  
(56)  
(57)  
(58)  
(59) 
Now if we choose the functionals,
(61)  
(62) 
then we obtain a generalized form of the Bennett Acceptance Ratio, as derived by Crooks. Crooks (2000)
Variational optimization of Eq. (60), however, leads to the functionals,
(63)  
(64) 
This expression is analogous to Bennett’s original expression. It can be solved selfconsistently or by rearrangement into,
(68) 
Using the sample mean estimator for the expectations, this is,
(69) 
which is similar to the expression of Shirts et. al. for the Bennett Acceptance Ratio Shirts et al. (2003). This is an implicit function of which is solved by finding the zero of the equation.
The variance of this expression can be found by plugging the optimal functionals into Eq. (60),