The Crooks theorem

Recent development on non equilibrium thermodynamics have clarified that the PMF along the given reaction coordinate $z$ can actually be reconstructed exactly using an ensemble of steered molecular dynamics simulations without resorting to any assumption on, or having any knowledge of the frictional behavior of the system along the reaction coordinate. These developments date back to a paper by Evans, Searls[146] where the first example of transient fluctuation theorem for a system driven out of equilibrium was formulated, demonstrating the connection between the time integral of the phase compression factor in Liouville space along an arbitrary time interval and the probability ratio of producing the entropy $A$ and $-A$ along a deterministic trajectory of a many particles non equilibrium steady state system. Gavin Crooks in his phd thesis proposed[64], in the context of Monte Carlo simulations in the canonical ensemble (NVT), a transient[146] fluctuation formula (from now on indicated with CT) involving the dissipative work for systems driven out of equilibrium by varying some arbitrary mechanical parameter. The CT is actually even more general than the Evans and Searls fluctuation theorem[146] since in the latter the driven $z$ coordinate has an underlying zero PMF (i.e. only entropy is produced in the non equilibrium process) while in the former the system can also cross different thermodynamics states (i.e. the underlying PMF can also be non zero such that thermodynamic work can also be done). The Crooks theorem (CT) reads

$$\frac{p(\Gamma(z_0) \rightarrow \Gamma(z_\tau))}{p(\Gamma^\ast(z_... ...)} = \exp [ \beta ( W_{\Gamma(z_0) \rightarrow \Gamma(z_\tau )} - \Delta F ) ],$$ (8.2)

where $\tau$ is the duration time of the driven non equilibrium process, $W_{\Gamma(z_0) \rightarrow \Gamma(z_\tau)}$ is the work done on the system during the driven trajectory $\Gamma(z_0) \rightarrow \Gamma(z_\tau)$; $p(\Gamma(z_0) \rightarrow \Gamma(z_\tau))$ is the joint probability of taking the microstate $\Gamma(z_0)$ from a canonical distribution with a given initial Hamiltonian $H(z=z_0)$ and of performing the forward transformation to the microstate $\Gamma(z_\tau)$ corresponding to a different Hamiltonian $H(z=z_\tau)$; $p(\Gamma^\ast(z_0) \leftarrow \Gamma^\ast(z_\tau))$ is the analogous joint probability for the time reversal path, producing the work $W_{\Gamma(z_\tau) \rightarrow \Gamma(z_0)}= -W_{\Gamma(z_0) \rightarrow \Gamma(z_\tau)}$. $\Delta F = F(z=z_\tau) - F(z=z_0)$ is the free energy difference between the thermodynamic states associated to the Hamiltonians $H(z=z_\tau)$ and $H(z=z_0)$. Although the CT can be stated in a more general formulation (see Gavin Crooks, $phd$ thesis), here the essential assumptions are that i) the system is deterministic and satisfies the time reversal symmetry and ii) the reverse trajectory is done following a reversed time schedule such that $W_{\Gamma(z_\tau) \rightarrow \Gamma(z_0)}= -W_{\Gamma(z_0) \rightarrow \Gamma(z_\tau)}$. The first assumption is satisfied by any kind of standard MD equation of motion (Newtonian, Nosé-Hoover. Parrinello-Rahman) while the second condition can be easily imposed in a SMD experiment. A very simple proof of Eq. 8.2 goes as follows: suppose the $z_0$ is drawn from a canonical distribution, and that the driven trajectory that brings the system to $z_\tau$ is done adiabatically, i.e. removing the thermal bath. For the reverse trajectory, drawing $z_\tau$ from a canonical distribution, due to the time reversal symmetry of the Hamilton equations, one ends up adiabatically in $z_0$. Under these assumptions, the ratio of the two probabilities on the left hand side of Eq. 8.2 can be written as

$$\frac{p(\Gamma(z_0) \rightarrow \Gamma(z_\tau))}{p(\Gamma^\ast(z_0) \leftarrow \Gamma^\ast(z_\tau))} = \frac{e^{-\beta H(z=z_0)} }{Z_0} \frac{Z_\tau}{e^{-\beta H(z=z_\tau)} }$$

$$= e^{\beta (H_{(z=z_\tau)} - H_{(z=z_0}) - \Delta F)}$$

$$= \exp \beta(W_{\Gamma(z_0) \rightarrow \Gamma(z_\tau )}-\Delta F)$$ (8.3)

equation where we have used the facts that $\beta F(z=z_0)= \ln Z_0$ , $\beta F(z=z_\tau)= \ln Z_\tau$ and that the energy difference $H_B-H_Z$ in the forward adiabatic trajectory equals to the external work done on the systems.

Figure: Physical significance of the Crooks theorem for a general driven process: for nearly reversible processes (left) the forward $P_f(W)$ and backward $P_b(-W)$ work distributions overlap significantly. The dotted line is the the backward work distribution for the inverse process, without changing the sing of the work. The crossing of the two solid distribution occurs at the free energy value for the forward process $\Delta F =18$. When the process is done faster (right panel), the dissipation $W_d$ both in the forward and in the backward process is larger, the overlap is negligible and the crossing point of the two solid distribution can no longer easily identified.

        

Equation 8.2 refers to the probability of a single forward or backward trajectory. Suppose now to perform a large number of forward trajectories all with a give time schedule, but each started from a different initial phase point sampled according to the canonical equilibrium distribution characterized by the Hamiltonian $H(z=z_0)$ and a large and not necessarily equal number of backward trajectories with reverse time schedule and starting from initial phase points this time sampled according to the canonical equilibrium distribution characterized by the Hamiltonian $H(z=z_t)$.^8.2. By collecting all trajectories yielding the work $W$ in (8.2), the CT may compactly be written as:

$$\frac{P_F(W)}{P_R(-W)} = \exp[ \beta (W - \Delta F) ],$$ (8.4)

where $P_F(W)$ and $P_r(W)$ are the normalized forward and backward distribution functions (note that, due to the time reversal symmetry, for the backward distribution the work is taken with the minus sign, i.e. $P_R(-W)$ is the mirror symmetric with respect to $P_R(W)$). According to Eq. 8.4, the $\Delta F$ may be thus evaluated constructing the two work distribution function: $\Delta F$ is the work value where the two distribution cross, i.e $P_F(W=\Delta) = P_B(-W=\Delta F)$. We point out in passing that, the famous Jarzynski identity[63] (JI),

$$< e^{-\beta W} > = e^{-\beta \Delta F},$$ (8.5)

is actually a trivial consequence of the CT, being derived from the latter by integrating out the work variable and using the fact that the work distribution function $P_F(W)$ and $P_R(-W)$ are normalized.

The physical meaning of the Crooks equation sounds indeed very reasonable and can be even be considered as a probabilistic restatement of the second law or of a generalization of the H-Boltzmann theorem: Given a forward deterministic non equilibrium trajectory starting form equilibrium and producing a work $W$, the probability to observe a trajectory for the reverse process again starting from equilibrium and producing the work $-W$ is $e^{\beta{W_d}}$ small than the former, where $W_d = W-\Delta F$ is the dissipated work in the forward process. When the dissipated work is zero, i.e. when the driven process is quasi-static and is done always at equilibrium, then the two probabilities are identical. With this regard, one important point to stress is that the CT and the JI hold for all systems and for any kind of arbitrary non equilibrium process, no matter how fast is performed. In particular, if the non equilibrium process is instantaneous, i.e. if it is done at infinite speed, then the work done on the system is simply equal to $W = (H_1-H_0)$, with $H_0$ and $H_1$ being the Hamiltonian of the initial and final state, respectively. The JI reduces in this case to the to famous free energy perturbation Zwanzig[119] formula $<e^{-\beta(H_1-H_0)}>_0 = e^{-\beta \Delta F}$ with the subscript 0 indicating that the canonical average must be taken according to the equilibrium distribution of the system with Hamiltonian $H_0$.

For fast non equilibrium experiments, a large amount of the work, rather than in advancing the reaction coordinate, is dissipated in heat that is in turn (only partly) assimilated by the thermal bath^8.3 A consequence of this is that the maxima of two work distributions $P_F(W)$ and $P_R(W)$ tend to get further apart from each other so that the determination of $\Delta F$ becomes less accurate. The faster are performed the non equilibrium experiments, the large is the average dissipation and the smaller is the overlap between the two work distributions (see Fig. 8.1) The reason why CT and JI can be so useful in evaluating the free energies along given reaction paths in the molecular dynamics simulation of complex biological system lies on the fact that this methodologies are inherently more accurate the smaller is the sample. Let's see why. As one can see form Fig. 8.1, $\Delta F$ can be determined with accuracy if the two work distributions overlap appreciably, or stated in other terms, if there are sufficient trajectories that in both directions transiently violate the second law, i.e trajectories for which $W < \Delta F$. This is clearly not in contrast with the second law which states that $\bar W \le \Delta F$ where $\bar W=\int P(W)WdW$ is the mean irreversible work. In general, the probability of an overlap of the two work distributions (i.e. the probability of transiently violating the second law) is clearly larger the smaller is the system.

Figure: Effect of the size of the system on the overlap of the forward and backward work distributions. In the left panel the non equilibrium processes are done in a given time $\tau$ on a single molecule. In the right panel the processes, as in left panel of duration $\tau$, are done independently on three identical molecules. This implies a factor 3 on energies and a factor $3^{1/2}$ on widths. As a result of the increased size, the overlap between $P_f(W)$ and $P_b(-W)$ decreases significantly.

        

Suppose to simultaneously and irreversibly unfold $N$ identical proteins in a dilute solution starting from their native states. In the assumption that the intraprotein interaction are negligible, the mean work for this system will be simply $N$ times the mean work done on a single molecule, while the width of the work distribution for the $N$ molecule systems will be only $N^{1/2}$ larger than that of the single molecule system. This effect is illustrated in Fig. 8.2. Now, biomolecular simulation of biosystems are usually done, for computational reasons, on a single solvated biomolecule, i.e. in the conditions where the non equilibrium techniques, for the reason explained above, are deemed to be more successful.