Asymptotic work distributions in driven bistable systems

Traditionally, statistical mechanics is concerned with averages; the probability distributions for thermodynamic quantities of macroscopic systems are so exceedingly sharp that only their most probable values matter. These in turn are practically identical with the averages. Only near instabilities, deviations from averages become important. In static cases these fluctuations are well described by the second moments of the respective distributions; if the dynamics is of importance as well, they are characterized by the correlation functions.

When investigating nanoscopic systems from the point of view of thermodynamics, the situation changes. Fluctuations are now strong and ubiquitous, and correspondingly, the probability distributions of relevant quantities are broad and poorly characterized by their leading moments alone. Whereas this seems rather obvious, it came as a real surprise that work [1] and fluctuation [2, 3] theorems which form the cornerstones of the emerging field of stochastic thermodynamics [4–8] are relations that probe the very far tails of the respective probability distributions. Very unlikely events now carry significant information about the physics of the system under consideration. On the one hand, interest in mathematical investigations like large deviation theory (see [9] and references therein) is renewed, on the other hand, techniques are in demand that allow to determine the asymptotics of probability distributions. Since rare events are hard to get in experiments and numerical simulations, approximate analytical procedures have to be developed.

In the present paper we apply a recently proposed method for the analytical determination of the asymptotics of work distributions in driven Langevin systems [10] to a simple model for single-molecule force spectroscopy. In these experiments (for a recent review see [11, 12]) the free-energy difference ΔF between the folded and the unfolded state of a biomolecule is determined from the distribution of work W obtained in isothermal unfolding and refolding processes. If only one transition is monitored [13], the Jarzynski equation [1]

$$\begin{eqnarray} &&\left\langle \mathrm{e}^{-\beta W} \right\rangle = \mathrm{e}^{-\beta \Delta F} \end{eqnarray} \tag{ 1 }$$

is employed, where β denotes the inverse of the temperature. If histograms of work for both the forward and the reverse process are compiled [14], it may be advantageous to use the Crooks fluctuation theorem [15]

$$\begin{eqnarray} &&\frac{P(W)}{P_{\mathrm{r}}(-W)}=\mathrm{e}^{\beta(W-\Delta F)}. \end{eqnarray} \tag{ 2 }$$

Here, the free-energy difference ΔF is identified from the intersection of the probability density functions P(W) of the forward and P_r(−W) of the reverse process. Note that in both cases an accurate estimate for ΔF requires reliable information about the tails of the work probability distributions.

The paper is organized as follows. In section 2 we introduce the notation and review the general method [10]. Section 3 contains the application to an one-dimensional stochastic process with a time-dependent double-well potential and the comparison with numerical simulations of the system. In section 4 we discuss our results and in particular compare them with a recent proposal for the general shape of the tail of the work distribution in single-molecule experiments [19].

To pave the way for the analysis and to fix the notation we summarize in the present section in a very condensed way the main steps of our approach to determine the asymptotic tail of work distributions in driven Langevin systems. For more details the reader is referred to [10].

We investigate a system described by an one-dimensional, overdamped Langevin equation which in dimensionless units has the form

$$\begin{eqnarray} &&\dot{x}=-V'(x,t) +\sqrt{2/\beta}\; \xi(t). \end{eqnarray} \tag{ 3 }$$

Here x denotes the degree of freedom, V is a time-dependent potential modelling the external driving, and ξ(t) is Gaussian white noise with 〈ξ(t)〉 ≡ 0 and 〈ξ(t)ξ(t')〉 = δ(t − t'). We denote derivatives with respect to x by a prime, and those with respect to t by a dot. The initial state x(t = 0) = :x₀ of the process is sampled from the equilibrium distribution at t = 0 with inverse temperature β,

$$\begin{eqnarray} &&\rho_0(x_0)=\frac{1}{Z_0}\exp\left[-\beta V_0(x_0)\right]. \end{eqnarray} \tag{ 4 }$$

Accordingly, V₀(x): = V (x,t = 0) is the initial potential and

$$\begin{eqnarray} &&Z_0=\int \mathrm{d}x \exp\left[-\beta V_0(x)\right] \end{eqnarray} \tag{ 5 }$$

the corresponding partition function.

The work performed on the system for a particular trajectory x(·) is given by [16]

$$\begin{eqnarray} &&W[x(\cdot)]=\int_0^{T} \mathrm{d} t\dot{V}\big(x(t),t\big). \end{eqnarray} \tag{ 6 }$$

Due to the randomness inherent in x(·), the work W is itself a random quantity and its pdf can be written as

$$\begin{eqnarray} P(W)&=&\int \frac{\mathrm{d}x_0}{Z_0}\exp\left[-\beta V_0(x_0)\right] \int \mathrm{d} x_T \nonumber\\ & & \times\int\nolimits_{x(0)=x_0}^{x(T)=x_T} { \cal D}x(\cdot)\; p[x(\cdot)] \; \delta(W- W[x(\cdot)]). \end{eqnarray} \tag{ 7 }$$

For mid-point discretization we have [17]

$$\begin{eqnarray} &&p[x(\cdot)]={\cal N}[x(\cdot)]\exp\left[-\frac{\beta}{4}\int_0^{T} \mathrm{d} t\; \big(\dot{x}+V'(x,t)\big)^2\right]\; \end{eqnarray} \tag{ 8 }$$

with the normalization factor

$$\begin{eqnarray} &&{\cal N}[x(\cdot)]=\exp\left[\frac{1}{2}\int_0^T \mathrm{d} t\; V''(x(t),t)\right]. \end{eqnarray} \tag{ 9 }$$

Hence

$$\begin{eqnarray} P(W)&=&\int \frac{\mathrm{d} x_0}{Z_0}\int \mathrm{d} x_T \int \frac{\mathrm{d} q}{4\pi/\beta}\nonumber\\ & & \times\int\nolimits_{x(0)=x_0}^{x(T)=x_T}{ \cal D}x(\cdot) \,{\cal N}[x(\cdot)] \exp\left\lbrace-\beta S[x(\cdot),q]\right\rbrace \end{eqnarray} \tag{ 10 }$$

with the action

$$\begin{eqnarray} &&S[x(\cdot),q] = V_0(x_0) + \int_0^{T} \mathrm{d} t \left[\frac{1}{4}(\dot{x} + V')^2 + \frac{{\rm i}q}{2}\dot{V}\right] - \frac{{\rm i}q}{2}W.\nonumber\\ \end{eqnarray} \tag{ 11 }$$

We are interested in the asymptotic behaviour of P(W). Rare values of W are brought about by unlikely trajectories x(·). In the spirit of the contraction principle of large deviation theory [9], we expect that in the asymptotic tails of P(W) there is one trajectory for each value of W that dominates P(W). To find it, we have to maximize P[x(·)]: = ρ₀(x₀)p[x(·)] under the constraint W = W[x(·)]. This can be done by using a saddle-point approximation of the integrals in (10). The result is

$$\begin{eqnarray} &&P(W)=\frac{\bar{{\cal N}}\sqrt{2}}{Z_0}\,\frac{\exp\left[-\beta\bar{S}\right]}{\sqrt{\bar{R}\det A}} \big(1+{O}(1/\beta)\big) . \end{eqnarray} \tag{ 12 }$$

To use this expression in explicit examples, we first have to determine the most probable trajectory $\bar {x}(\cdot )$ satisfying the Euler–Lagrange equation (ELE)

$$\begin{eqnarray} &&\ddot{\bar{x}}+(1-{\rm i}\bar{q}){\dot{\bar{V}}}{}'-\bar{V}'\bar{V}''=0, \end{eqnarray} \tag{ 13 }$$

where $\bar {V}(t):=V(\bar {x}(t),t)$ and similarly for derivatives of V . The ELE is completed by the boundary conditions

$$\begin{eqnarray} &&\dot{\bar{x}}_0-\bar{V}'_0=0, \quad \dot{\bar{x}}_T+\bar{V}'_T=0 \end{eqnarray} \tag{ 14 }$$

and by the corresponding value $\bar {q}$ of the Lagrange parameter q ensuring $W[\bar {x}(\cdot )]=W$. Using $\bar {x}(\cdot )$ and $\bar {q}$, we calculate $\bar {S}:=S[\bar {x}(\cdot ),\bar {q}]$ and $\bar {{\cal N}}:={\cal N}[\bar {x}(\cdot )]$. Then all terms in (12) depending solely on the optimal trajectory are determined.

The denominator $\sqrt {\bar {R}\det A}$ in (12) comprises the contribution from the neighbourhood of the optimal path and stems from the integral over the Gaussian fluctuations around $\bar {x}(\cdot )$ and $\bar {q}$. Here

$$\begin{eqnarray} &&A:= -\frac{\mathrm{d}^2}{\mathrm{d} t^2}+(\bar{V}'')^2+\bar{V}'\bar{V}'''-(1-{\rm i}\bar{q})\dot{\bar{V}}{}'' \end{eqnarray} \tag{ 15 }$$

denotes the operator of quadratic fluctuations which acts on functions φ(t) on the interval 0 < t < T obeying the boundary conditions

$$\begin{eqnarray} &&\begin{array}{@{}c@{}} \bar{V}_0''\varphi(0)-\dot{\varphi}(0)=0,\\ \bar{V}_T''\varphi(T)+\dot{\varphi}(T) =0. \end{array} \end{eqnarray} \tag{ 16 }$$

A simple prescription to calculate $\det A$ is as follows [18]. Solve the initial value problem

$$\begin{eqnarray} A \chi(t)&=&0,\nonumber\\ \chi(0)&=&1, \qquad \dot{\chi}(0)=\bar{V}_0'', \end{eqnarray} \tag{ 17 }$$

then

$$\begin{eqnarray} &&\det A=2\big(\bar{V}_T''\chi(T)+\dot{\chi}(T)\big). \end{eqnarray} \tag{ 18 }$$

The factor $\bar {R}$ in (12) accounts for the influence of the constraint (6) on the fluctuations around $\bar {x}(\cdot )$. Since also the trajectories from the neighbourhood of $\bar {x}(\cdot )$ have to yield the very same value of W, fluctuations violating this constraint are suppressed. This gives rise to a correction to the fluctuation determinant of the form

$$\begin{eqnarray} &&\bar{R}=\int_0^T \mathrm{d}t\,\dot{\bar{V}}{}'(t) \; A^{-1} \dot{\bar{V}}{}'(t) , \end{eqnarray} \tag{ 19 }$$

where A⁻¹ denotes the inverse operator of A. To explicitly determine $\bar {R}$, it is convenient to solve the ordinary differential equation

$$\begin{eqnarray} &&A\, \psi(t)=\dot{V}'(\bar{x}(t),t) \end{eqnarray} \tag{ 20 }$$

with boundary conditions (16) and to use

$$\begin{eqnarray} &&\bar{R}=\int_0^T \mathrm{d} t\, \psi(t) \dot{\bar{V}}{}'(t). \end{eqnarray} \tag{ 21 }$$

The unfolding and refolding of single molecules can be modelled by a time-dependent double-well potential of the form

$$\begin{eqnarray} &&V(x,t) = ax^4 - bx^2 + r(c-t)x , \end{eqnarray} \tag{ 22 }$$

where x denotes the extension of the molecule in the direction of the force [13, 14, 19]. The parameters a and b characterize depth and separation of the two minima of V , c fixes the moment at which V is symmetric, V (x) = V (−x), and r denotes the transition rate. Choosing T > c for the final time, the two minima will interchange global stability during the process. We have used two exemplary sets of parameters for which the time evolution of V (x,t) is sketched in figure 1. The main differences between the two sets are that in (b) the minima are further apart and the transition rate r is smaller.

Zoom In Zoom Out Reset image size

We now apply the analytic method to obtain the asymptotics (12) of the work distribution of the dynamics defined by the potential V (x,t) (22). To clarify the procedure, we compile the necessary equations. We have from (5) and (9)

$$\begin{eqnarray} &&Z_0 = \int\, \mathrm{d} x\; \exp\left[-\beta(ax^4-bx^2+rcx)\right], \end{eqnarray} \noindent \tag{ 23 }$$

$$\begin{eqnarray} &&{\cal N}[x(\cdot)] = \exp\left[6a\int_0^T \mathrm{d} t\; x^2(t)-bT \right]. \end{eqnarray} \tag{ 24 }$$

The ELE (13) reads

$$\begin{eqnarray} \ddot{\bar{x}} &=& 48a^2\bar{x}^5 - 32a\bar{x}^3b + r(c-t)\big[12a\bar{x}^2-2b\big]\nonumber \\ & & + 4b^2\bar{x} + r(1-{\rm i}\bar{q}), \end{eqnarray} \tag{ 25 }$$

and its boundary conditions (14) are of the form

$$\begin{eqnarray} &&\begin{array}{@{}c@{}} 0 = \dot{\bar{x}}_0 - 4a\bar{x}_0^3+2b\bar{x}_0-rc,\\ 0 = \dot{\bar{x}}_T + 4a\bar{x}_T^3-2b\bar{x}_T+r(c-T). \end{array} \end{eqnarray} \tag{ 26 }$$

The constraint (6) is

$$\begin{eqnarray} &&W = - r \int_0^T \mathrm{d} t\;\bar{x}(t;\bar{q}) . \end{eqnarray} \tag{ 27 }$$

Using a standard relaxation algorithm, the numerically solution of equations (25)–(27) for desired work values W yields optimal trajectories $\bar {x}(t)$ depicted exemplarily in figure 2.

Zoom In Zoom Out Reset image size

The operator A from (15) acquires the form

$$\begin{eqnarray} &&A = -\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 240a^2\bar{x}^4 - 96ab\bar{x}^2 + 24ar(c-t)\bar{x} + 4b^2 \end{eqnarray} \tag{ 28 }$$

with the boundary conditions (16)

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} 0 = \big[12a\bar{x}_0^2-2b\big]\varphi(0)-\dot\varphi(0),\\ 0 = \big[12a\bar{x}_T^2-2b\big]\varphi(T)+\dot\varphi(T). \end{array} \end{eqnarray} \tag{ 29 }$$

To obtain $\det A$ according to (18), we determine

$$\begin{eqnarray} &&\det A = \big[24a\bar{x}_T^2-4b\big]\chi(T) + 2 \dot\chi(T) \end{eqnarray} \tag{ 30 }$$

by solving numerically the initial value problem (17)

$$\begin{eqnarray} \ddot\chi(t) &=& \big[240a^2\bar{x}^4 - 96ab\bar{x}^2+ 24ar(c-t)\bar{x} + 4b^2\big]\chi(t) = 0,\nonumber\\ \dot\chi(0) &=& 1,\quad \dot\chi(0) = 12a\bar{x}_0^2-2b. \end{eqnarray} \tag{ 31 }$$

This has to be done for each value of W separately by using the appropriate results for $\bar {x}(t;W)$ and $\bar {q}(W)$.

The last ingredient for the pre-exponential factor is $\bar {R}$ from (21). To determine it, we need to solve the boundary value problem (20) and (16)

$$\begin{eqnarray}\ \ddot\psi(t)&=&\bigl[240a^2\bar{x}^4 - 96ab\bar{x}^2+ 24ar(c-t)\bar{x} + 4b^2 \bigr]\psi(t) + r , \nonumber\\ 0 &=& \big[12a\bar{x}_0^2-2b\big]\psi(0)-\dot\psi(0),\\ 0 &= &\big[12a\bar{x}_T^2-2b\big]\psi(T)+\dot\psi(T)\nonumber \end{eqnarray} \tag{ 32 }$$

for each $\bar {x}(t;W)$ and $\bar {q}(W)$ and use the result in (21)

$$\begin{eqnarray} &&\bar{R} = -r\int_{0}^{T} \mathrm{d} t\;\psi(t) . \end{eqnarray} \tag{ 33 }$$

Plugging the numerical results for Z₀, $\bar {{\cal N}}$, $\bar {S}$, $\det A$ and $\bar {R}$ into (12), we obtain the final result for the asymptotic form of the work distribution. We carried out this program for the two processes characterized by the parameter sets given in the caption of figure 1, including for both cases the reversed processes defined by the substitution t → (T − t). From simulations of the Langevin equation (3), we also obtained from (6) the corresponding histograms of work distributions. The values of β and T are chosen such that most of the trajectories reach the right minimum of V at the end of the forward process, i.e. the molecules are stretched until virtually all of them unfold (cf figure 3). The results for the asymptotics are shown in figure 4, together with the outcome of our numerical simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In a recent paper [19], Palassini and Ritort propose an universal form for the tails of work distributions for single molecule stretching experiments given by

$$\begin{eqnarray} &&P(W) \sim n \frac{\Omega^{\alpha-1}}{|W-W_{\mathrm{c}}|^\alpha}\,\exp\left[-\frac{|W-W_{\mathrm{c}}|^\delta}{\Omega^\delta}\right] . \end{eqnarray} \tag{ 34 }$$

Here, W_c is a characteristic work value, Ω > 0 measures the tail width, and n, α and δ are further constants. The Jarzynski equation (1) stipulates δ > 1. Based on (34), Palassini and Ritort present in [19] three slightly different analytical methods to improve the estimation of free-energy differences ΔF from the Jarzynski equation (1). To decide which approach to use, they distinguish three regimes defined by the parameter combination

$$\begin{eqnarray} &&\lambda := \left(\delta/\Omega\right)^{\delta/(\delta-1)} \,\ln N . \end{eqnarray} \tag{ 35 }$$

The three regimes are then identified by λ > 1, λ ≪ 1 and λ ≲ 1 respectively. They test their method with experimental data from DNA stretching experiments. For more details regarding the improved estimation of ΔF and the experiments see [19].

We fitted the empirical asymptotic form (34) to the tails of the work distributions obtained in our simulations by standard least-square fits, starting with a Gaussian distribution specified by $W_{\mathrm {c}}=\left \langle W \right \rangle$, $\Omega ^2=2\left \langle (W-W_{\mathrm {c}})^2 \right \rangle$, α = 0, δ = 2 and $n=1/\sqrt {\pi }$. A subtle point in the procedure is to find the optimal interval of work values for the fit. The resulting parameter values are listed in table 1, the corresponding fits are included into figure 4.

As shown in figure 4 for both parameter sets (a) and (b), we achieve an excellent agreement between simulation and asymptotics, not only asymptotically for the tails, but also for the centre of the work distributions. The good reproduction of the centre of the distribution is presumably due to the fact that also the probabilities of typical work values are dominated by single trajectories and their neighbourhood. Note also the perfect match between the free-energy difference ΔF and the intersection of our asymptotics P(W) and P_r(−W), which demonstrates that the Crooks relation (2) holds in its exact form for our asymptotics, as has been shown in [10].

In addition, figure 2 illustrates the optimal trajectories $\bar {x}(\cdot )$ which dominate the asymptotics of the work distributions (12). In comparison with the trajectories obtained from our simulations shown in figure 3, the trajectories $\bar {x}(\cdot )$ contributing to the tails of the work distribution are of much broader variety. Interestingly, the probability of small work values is dominated by $\bar {x}(\cdot )$ that run into the evolving right minimum even before the minimum is shaped.

For the parameter set (a), the fit of (34) compares well with the analytic asymptotics. In this case a combination of a histogram from experimental values and (34) would therefore result in reliable estimates for the free-energy difference ΔF. For the parameter set (b) only the forward process is well described by the fit, for the reverse one the tail of the distribution is markedly overestimated. This shifts the intersection point between P(W) and P_r(−W) away from the correct value of ΔF as can be seen in figure 4(b). Also, some fit parameters for this case listed in table 1 clearly deviate from the other cases. To investigate that mismatch more closely, we also fitted (34) to the analytic asymptotics (not shown). This results into similar conspicuous fit parameters, but the fitting curve now is almost congruent with the analytic asymptotics. From that we conclude that the empirical asymptotics (34) is also valid in this case, but the number of work values obtained from the simulation is not enough to reliably fit (34). Note that rather than taking the intersection point between P(W) and P_r(−W), Palassini and Ritort use in [19] much more sophisticated methods to obtain estimates of ΔF, based on their empirical asymptotics (34). Our investigation aims only at validating (34).

In addition to the two parameter sets displayed in figure 1, we also investigated several other realizations of the potential (22) (not shown). Mostly, we found that fits of (34) to histograms of 10⁴–10⁵ work values extrapolate reliably to the asymptotic tail of the distribution. But as exemplified by the case shown in figure 4(b), there is no guarantee that a number of 10⁵ work values is sufficient.

The asymptotics of work distributions for driven Langevin systems can be determined by the method proposed in [10]. We employed this method to determine the asymptotics of a system defined by the potential (22) modelling the stretching of single molecules. The unfolding of the molecules corresponds to the forward process, the refolding to the reverse process. We obtained histograms of work by simulating the Langevin equation (3). The form of the work distributions are found to be near-Gaussian, similar to distributions measured in a DNA stretching experiment [19]. We observe excellent agreement between asymptotics and work distribution, not only for the asymptotic regime but also for the whole range of work values.

One aim of single molecule experiments is to obtain the free-energy difference between the folded and unfolded state of the molecule. If both the work distribution for the forward and reverse process is available, the Crooks relation (2) can be used to determine the free-energy difference. It is shown in [10] that the asymptotics (12) generally satisfies the Crooks relation exactly, which we demonstrated for two representative examples of the potential (22).

Finally, we tested the universal form (34) of the tails of work distributions proposed by M Palassini and F Ritort against our results for the asymptotics (12). For a broad range of parameters used for the model potential (22), we found a good agreement between (34) and our asymptotics. Only if the work distribution differs markedly from a Gaussian form, a reliable fit of (34) is likely to require more data points than usually acquired in single-molecule experiments. This might lead to a significant difference between the exact and estimated free-energy difference as illustrated by our examples.

We would like to thank Felix Ritort for interesting discussions and the organizers of the 2011 Nordita workshop Foundations and Applications of Non-equilibrium Statistical Mechanics where part of this work was done, for providing a stimulating atmosphere, and the ESF for financially supporting the workshop. DN acknowledges financial support from the Deutsche Forschungsgemeinschaft under grant EN 278/7.