Equilibrium binding energies from fluctuation theorems and force spectroscopy simulations

Brownian dynamics simulations are used to study the detachment of a particle from a substrate. Although the model is simple and generic, we attempt to map its energy, length and time scales onto a specific experimental system, namely a bead that is weakly bound to a cell and then removed by an optical tweezer. The external driving force arises from the combined optical tweezer and substrate potentials, and thermal fluctuations are taken into account by a Brownian force. The Jarzynski equality and Crooks' fluctuation theorem are applied to obtain the equilibrium free energy difference between the final and initial states. To this end, we sample non--equilibrium work trajectories for various tweezer pulling rates. We argue that this methodology should also be feasible experimentally for the envisioned system. Furthermore, we outline how the measurement of a whole free energy profile would allow the experimentalist to retrieve the unknown substrate potential by means of a suitable deconvolution. The influence of the pulling rate on the accuracy of the results is investigated, and umbrella sampling is used to obtain the equilibrium probability of particle escape for a variety of trap potentials.

Brownian dynamics simulations are used to study the detachment of a particle from a substrate. Although the model is simple and generic, we attempt to map its energy, length and time scales onto a specific experimental system, namely a bead that is weakly bound to a cell and then removed by an optical tweezer. The external driving force arises from the combined optical tweezer and substrate potentials, and thermal fluctuations are taken into account by a Brownian force. The Jarzynski equality and Crooks fluctuation theorem are applied to obtain the equilibrium free energy difference between the final and initial states. To this end, we sample non-equilibrium work trajectories for various tweezer pulling rates. We argue that this methodology should also be feasible experimentally for the envisioned system. Furthermore, we outline how the measurement of a whole free energy profile would allow the experimentalist to retrieve the unknown substrate potential by means of a suitable deconvolution. The influence of the pulling rate on the accuracy of the results is investigated, and umbrella sampling is used to obtain the equilibrium probability of particle escape for a variety of trap potentials.

I. INTRODUCTION
The adhesion of a cell to a substrate 1-3 occurs in a number of biophysical contexts, and is hence a very important phenomenon to study. Beyond its relevance for understanding biological phenomena in general, many clinical applications in both diagnostics and therapeutics fundamentally involve adhesion. Examples include: (i) the sequestration of red blood cells in small blood vessels due to infection with malaria, [4][5][6][7] (ii) the growth of metastases in cancer, [8][9][10] and (iii) the formation of platelets at the site of a vascular injury. 11 A variety of experimental techniques have been developed 12 to measure the adhesive properties of a single cell, such as atomicforce microscopy, [13][14][15][16] surface-force apparatus measurements, 17 micropipette manipulation, [18][19][20] as well as magnetic 21 and optical [22][23][24] tweezers. All these methods subject the cell to external time-dependent forces, with the aim of quantifying the energetics of the binding.
The theoretical framework to analyse such experiments are the recently developed non-equilibrium work theorems, 25,26 most notably the Jarzynski theorem, 27,28 and a) Electronic mail: ravi.jagadeeshan@monash.edu Crooks fluctuation theorem, [29][30][31] which have been used with great success to interpret data from both computer simulations and experiments. [32][33][34][35][36][37] These theorems combine in a coherent fashion the three salient aspects of the experiments, which are (i) the system's equilibrium statistical physics (in particular the binding enthalpy), (ii) the fact that time-dependent manipulation necessarily implies non-equilibrium statistical physics (where the degree of deviation from equilibrium is determined by the pulling speed or a similar parameter), and (iii) the influence of thermal fluctuations. The central quantity of the theorems is the non-equilibrium work that the external forces do on the system. As soon as the external driving happens on a time scale that is faster than the typical relaxation times of the system, the non-equilibrium work is no longer simply given by the free energy difference between final and initial state (as would be the case for infinitely slow or quasi-static driving), but rather acquires a dissipative contribution, which, as a result of thermal fluctuations, has a statistical distribution of values. The theorems make detailed statements on the relation between the probability distribution of the nonequilibrium work and the underlying equilibrium free energies, and are hence immensely useful to obtain the latter under experimental conditions that cannot be considered as quasi-static. Essentially the extraction of equi-librium properties from the non-equilibrium work distribution is tantamount to reweighting the latter. Therefore, the theorems, although in theory being applicable to a large class of physical situations, have limitations in practice, since the equilibrium free energy difference should not differ from the mean non-equilibrium work by more than a few standard deviations -and this becomes more and more unfavourable both with increasing dissipation and increasing system size. In practice, this means that a reliable acquisition of equilibrium properties requires more and more trajectories over which one needs to average. 27,28 In this context, it should be noted that the theorems always consider transitions from an equilibrium initial state to a final state, which is typically out of equilibrium. These states are not given by some reaction coordinate of the system, but rather by the external driving. Furthermore, we would like to mention that not only free energies, but also other equilibrium properties (like e. g. the probability of attachment) can be obtained in an analogous fashion by a suitable reweighting (or "umbrella sampling") procedure.
Binding between cells is complex and involves a slew of interactions, which are both specific and nonspecific. [38][39][40] The most important ingredient, however, are bonds that arise from receptor-ligand pairs. Typically, a single receptor-ligand interaction is fairly strong, i. e., of order of few k B T to 100 k B T , where k B is Boltzmann's constant and T the absolute temperature at ambient conditions, 17,41 i. e. k B T 4pN nm. Moreover, cell adhesion will in most cases involve many ligands, giving rise to a net total interaction of typically several hundred k B T . The mechanical detachment of a cell from a "substrate" to which it is bound via receptorligand pairs (the latter can for example be another cell, or a ligand-coated bead) is thus a very complex process. 38,42 In a highly simplified picture, we envision it to be roughly analogous to the pulling-off of a plaster from skin, or to the pinch-off of a water droplet from a dripping faucet. In optical tweezer experiments 43 we have observed that the same external force can be sufficient to break some cell-substrate pairs but insufficient to break others of the same type. In our opinion, this provides an indication that the underlying dynamics matter, and this will depend on details of variables such as the number of receptor-ligand pairs, their density, and their geometrical arrangement. At any rate, this means that a faithful modeling of cell-substrate detachment or attachment would need to take into account a large arrangement of receptor-ligand bonds, and their (elastic) interactions. The single pair, in turn, is weak enough that thermal fluctuations crucially contribute to its formation and breaking.
As a first step in the modeling of micromechanical manipulation of cell attachment and detachment, we focus in the present paper on the case of just a single ligandreceptor pair. This is clearly the easiest situation, since in principle this allows us to just consider a single coordinate x as a degree of freedom, which may be viewed as the cell-substrate distance. This degree of freedom can then be viewed as subject to (i) forces from the cell-substrate interaction, (ii) forces from the time-dependent external pulling, and (iii) thermal agitation. This situation is less artificial than one might think at first glance, since it is experimentally possible to modify the adhesive properties of cells through gene-knockout techniques and/or inhibitors, [44][45][46][47][48] such that receptor-ligand interactions are systematically turned off. The aim of the present theoretical study is to demonstrate that in this weak-binding situation the theorems can actually be applied practically to obtain reliable results on free energies, and, as a consequence, on the binding energetics. To do this, we study the attachment or detachment process within the framework of a very simple theoretical model, whose dynamics is simulated by means of Brownian Dynamics. An important aspect here is the fact that the simulation parameters (strength and range of interactions, pulling speed) roughly match those of real experiments. In the subsequent sections we will provide details on the choice of parameters, and discuss the relation between the free energies from the fluctuation theorems on the one hand, and the binding forces on the other.
It should be emphasised that our numerical model is fairly generic and therefore in principle applicable to any micromechanical manipulation that detaches one object from another (or attaches it to it), as long as this process can be described by a single reaction coordinate, and involves energies that are roughly comparable with k B T . However, what we have principally in mind are experiments with optical tweezers. We believe this technique has a great potential in the future, since it is fairly noninvasive, and provides good quantitative control over the external forces involved. For this reason, we choose our parameters in rough accordance with a typical tweezer experiment, and also use a nomenclature that refers to this situation. More precisely, we think of a cell tightly "glued" to a glass surface, 22 while a ligand-coated bead is moved due to the influence of a time-dependent (harmonic) tweezer potential. The forces that the cell exerts on the bead are then described by a fixed (not timedependent) "membrane potential".
It is worth noting that fluctuation theorems have already been used to computationally calculate binding free energies in drug-receptor systems. 32,33 These computations involve deterministic nonequilibrium molecular dynamics of ligand-receptor pairs whose molecular properties, such as Lennard-Jones parameters and force fields are known. In this paper, the analysis of single cell detachment events will be described and the usefulness of fluctuation theorems demonstrated, using data generated by stochastic simulation of a model cell and substrate. Since the situation in the numerical study is fairly similar to a typical experiment, we believe that this also demonstrates the usefulness of the approach to experimentally estimate the strength of binding -with the caveat that the experiments will be less accurate, since it is experimentally not possible to study O(10 6 ) trajec-tories, as was done in the present investigation.
The remainder of the paper is organised in the following manner: First, details of the Langevin simulation will be presented, including code validation. Second, the Jarzynski and Crooks fluctuation theorems are shown to be valid for this two state system. As a result, nonequilibrium work trajectories, calculated for the different trap velocities, can be used to obtain the equilibrium free energy difference between the final and the initial state. We will also briefly outline (although this has not been done in the present work) how this information can in principle be used to retrieve the membrane potential, which in an experiment is of course unknown. Third, limitations of numerical calculations using the fluctuation theorems will be discussed and illustrated with the use of cumulants. Finally, umbrella sampling will be used to derive equilibrium values such as the probability of detachment or adhesion for a variety of different trap potentials.

II. PROBLEM FORMULATION
A. The model unbinding experiment Truncated harmonic potentials, as shown in Figs. 1 (a) and (b), are used to describe the interaction of the bead with both the membrane and the optical trap. These potentials are made dimensionless by scaling with the natural energy scale k B T , and defined by the expressions and where U M and U OT are the dimensionless membrane and optical trap potential energies, respectively. The distance x, measured from the fixed location of the minimum of the membrane potential, is made dimensionless by scaling with a length k B T /k s , where k s is a typical spring constant. We now choose the dimensionless parameters M and OT of order unity, which means that the involved energy scales are O(k B T ), as in the envisioned experiments. Furthermore, we assume that the spring constant k s is a value that corresponds to a typical optical trap strength of O(10 −3 pN/nm), 17 implying that k OT is a dimensionless parameter of order unity. At ambient conditions, k B T 4pN nm, meaning that a typical thermal displacement within the trap (which is our unit of length) is several tens of nanometers. The typical displacements that we observe for cell detachment 43 are of similar order, and therefore we set k M as a parameter of order unity as well.
The repulsive segment of the membrane potential The total potential, U = UM + UOT, experienced by the bead at some time t > 0. In order to detach from the membrane the bead needs an energy greater than M, while in order for the bead to go from being unattached to attached, it would require an energy of order OT or greater.
(−∞ < x ≤ 0) accounts for the impenetrability of the membrane to the bead, while the attractive segment (0 < x ≤ x ub M ) represents the adhesive force exerted by the membrane on the bead ( Fig. 1(a)). Beyond this distance, the bead detaches from the membrane and the influence on the bead by the membrane potential becomes negligible. Note that the minimum of the potential is held fixed at the origin (x = 0) for all time. Traditionally optical tweezer potentials are represented by harmonic wells. 37,49 However, for investigations of detachment or attachments one should take into account that the optical trap has a finite range of attraction as well, such that a truncated harmonic potential is more reasonable. In principle this consideration holds for both branches x < x OT and x > x OT , where x OT is the (timedependent) location of the minimum of U OT . However, it is crucially important only for x < x OT because this controls the energy barrier between the membrane and the trap potential. For x > x OT we do not truncate the tweezer potential, in order to obtain finite expressions in the equilibrium statistical mechanics of the system: If the total potential would exhibit an infinite range of vanishing potential, then this region would correspond to an infinite translational entropy, meaning that at any finite temperature there could be no equilibrium adsorption of the bead. Dynamically, this behavior would correspond to "evaporation" of the bead at sufficiently long times. It is therefore reasonable to study the particle in a potential that results in a converging partition function, and by this to strictly disregard such "evaporation" events (which, in a typical experiment, are anyway not observed). These considerations lead us to assume a model tweezer potential U OT as depicted in Fig. 1 (b). The total potential, U (x) = U M (x) + U OT (x), at some time t > 0, is shown schematically in Fig. 1 (c).
The optical trap potential minimum is located at the origin at time t = 0, i. e., x OT (t = 0) = 0. At later times, the optical trap is translated horizontally linearly with time, at varying speeds v OT (i. e., x OT (t) = v OT t), in order to simulate the process of bead detachment by the optical trap. The final position of the trap minimum is always at a fixed location, x final OT = 6, regardless of the value of v OT . The summed potential U is time dependent because of the time dependence of the optical potential. For the purpose of illustration, the shapes of the membrane and optical trap potentials, along with the summed potential, during the course of the simulation, at three different locations of the optical trap minimum are shown in Fig. 2.
The relative ease of attachment and detachment is controlled by the magnitudes of the barrier heights for the membrane ( M ) and the optical tweezer ( OT ) potentials, respectively, and also by their respective strengths k M and k OT . In order to model different adhesive interactions between the bead and the membrane, the barrier heights and spring constants can be changed appropriately. In the present work, we choose three different sets of values for these parameters (given in Table I), allowing different scenarios to be tested, as illustrated in Fig. 3. In Fig. 3 (a), the membrane potential is weaker than the optical trap in both strength and depth. In Fig. 3 (b), both the potentials have the same strength and depth, with the dimensional depth being of order 10 k B T , while in Fig. 3 (c), their dimensional depths are of order 1k B T . As will be seen subsequently, these three different scenarios lead to considerably different adhesive behaviour.

B. The Langevin equation
In the absence of inertia, the time evolution of the particle's position x(t), subject to an external force due to the presence of the membrane and optical potentials, and subject to thermal fluctuations, is described by a Langevin equation where the coordinate x is dimensionless as described above, and time is also made dimensionless by scaling with the typical time scale ζ/k s , ζ being the friction coefficient of the particle. F ext is the dimensionless external force due to the combined potential, given by F ext = −∂U/∂x, while F rand is the dimensionless random force (Gaussian white noise) with mean and variance We use an Euler algorithm with a time step ∆t, to numerically integrate the Langevin equation. Here r is a random number with r = 0 and r 2 = 1. We use Gaussian random numbers, applying the standard Box-Muller method. Details of time step sizes and the number of trajectories used in the simulations are given in the context of the various results discussed below.
Assuming a typical bead radius of 4 µm, and an aqueous environment with viscosity 10 −3 Pa s, we find a  Table I.
Stokes friction coefficient of 0.075 × 10 −3 pN s/nm, meaning that for a spring constant of 10 −3 pN/nm our unit of time is 0.075 seconds.
The non-equilibrium aspect of the computer experiment comes in through the finite pulling rate v (the velocity at which the location of the tweezer potential travels). For this we choose dimensionless values between 0.01 and 1. In experimental units, this means that even for the fastest process we pull the bead on a time scale of not much less than roughly 0.1 seconds, over a length scale of a few ten nanometers, which means that the simu-lated process corresponds well to experimentally feasible scales.

C. Fluctuation theorems
The initial and final states of our system are respectively defined as (i) x OT = 0, a situation where the tweezer potential keeps the bead at a location close to the membrane, and (ii) OT where it has moved the bead quite far away from it, such that it feels only the force from the optical trap. The fluctuation theorems are concerned with the free energy difference ∆F between these two states.
If the unbinding is carried out isothermally and infinitesimally slowly, then ∆F is equal to the work W performed during the process. On the other hand, if the unbinding experiment is carried out at a finite rate over a period of time t D , the work performed will not be unique. Rather, an ensemble of such unbinding experiments will lead to a distribution of work values, P F (W ) (where the subscript 'F' indicates the experiment is carried out in the forward direction, from the cell and bead being bound together to being unbound). Note that in this scenario, it is possible that at the end of the experiment, the bead remains close to the cell, even though work has been performed. In the quasi-static limit t D → ∞, . For finite rates of detachment, however, The great advance that has been made with the recently developed fluctuation theorems is that, contrary to the suggestion of Eqn. (6), a knowledge of the nonequilibrium work distribution is sufficient to determine the equilibrium free energy ∆F exactly.
The two fluctuation theorems that are primarily used in this work are the Crooks fluctuation theorem, [29][30][31] and the Jarzynski equality. 27,28 Both these theorems are based on the following set of assumptions. The system, whose dynamics is in our case stochastic and Markovian, is driven by an external perturbation from an initial equilibrium state, to a final state that is not necessarily at equilibrium. The external parameter driving the perturbation at a finite rate from the initial to the final state is denoted by λ, with values λ 0 in the initial equilibrium state, and λ f in the final state.
The Crooks fluctuation theorem states that [29][30][31] where both the work and the free energy have been made dimensionless by scaling with our energy unit k B T . The distribution P F (W ) is the probability that the work of magnitude W is performed in perturbing the system from an initial equilibrium state with λ = λ 0 to a final state with λ = λ f in a finite time t D , while P R (−W ) is the probability that work of the same magnitude but opposite sign will be performed on perturbing the system in the reverse path, from an equilibrium state with λ = λ f to a state with λ = λ 0 , over the same length of time.
Equation 7 clearly suggests that the value of work W * at which P F (W * ) = P R (−W * ), is nothing but the equilibrium free energy difference between the initial and final states. We use this result subsequently in order to estimate the free energy of binding.
The Jarzynski equality in its original form 27,28 only considers perturbations from λ 0 to λ f , and states that where the subscript 'F' on the ensemble average on the left hand side indicates an average over forward trajectories. While the ensemble average of the non-equilibrium work is always greater than the equilibrium free energy for finite rates of system perturbation, Jarzynski's equality states that an ensemble average of the exponential of (−W ) can be used to directly evaluate the equilibrium free energy. As will be seen subsequently, however, driving the system from λ 0 to λ f at increasingly rapid rates leads to a widening of the distribution P F , and consequently requires larger and larger ensembles to obtain an accurate estimate of ∆F .
The experimental and practical relevance of these relations becomes clear when considering the defining relation for the free energy, where we emphasise that the tweezer potential depends on the difference x − x OT . Now, the fluctuation theorems permit us to determine the free energy not only for the final state of the tweezer potential, but also for any intermediate state x interm OT . We thus find Defining which we can assume to be known since the properties of the optical trap are known, and which is not known, we can write In other words, the exponential of the free energy profile, which is experimentally accessible via the fluctuation theorems, is nothing but the convolution of the known Boltzmann factor of the tweezer potential with the unknown Boltzmann factor of the membrane potential. Therefore, it should be possible to retrieve the latter by just a numerical deconvolution, assuming that the free energy profile is known with sufficient accuracy. More precisely, the procedure yields U M up to an unknown constant, which is however obviously irrelevant. Mapping out the membrane potential is, in our opinion, the ideal goal of such experiments. In the present work, we do not perform this program, but rather confine ourselves to the simpler task of just determining ∆F for a single final state.

D. Non-equilibrium work
The application of the fluctuation theorems requires the determination of the distribution of work P F (W ) when the system is driven from λ 0 to λ f in the forward path, and the distribution P R (W ) when the path is reversed. Following the arguments of Jarzynski, 28 we introduce the function H λ (x), as the energy of the system for any fixed value of λ, where x(t) is the stochastic phasespace trajectory that describes the time evolution of the system, which depends on the time dependence of the external parameter λ. The total work performed on the system, when it evolves from λ = λ 0 to λ = λ f , in a time period t D , is 28 whereλ = dλ/dt. The stochastic phase-space trajectory x(t) of the bead is determined here by solving the Langevin equation (3). In the model system considered here, the only component of the system's energy that depends on the external driving parameter λ (= x OT ), is the potential energy of the trap, U OT . As a result, Equations (15) and (16) are used here to calculate the work done on the bead when the optical trap is translated from x OT = 0 to x OT = x final OT , at all times t at which the bead's location satisfies x(t) ≥ x lb OT . At other times, when the force of the optical trap on the bead is zero, the contribution to the work is zero. At any time t during the course of the Langevin simulation, the accumulated work until time t is calculated by numerically evaluating the integral in Eqn. (16) from t = 0 to t = t. Since the typical time steps used in the simulation are very small (∆t = 10 −4 to ∆t = 10 −3 ), a simple rectangular method was used to carry out the quadrature, where at each time step, the accumulated work at the end of the previous time step is augmented by the product of the value of the integrand at the beginning of the time step with ∆t.

E. Analytical evaluation of the free energy
For the simple model considered here, the free energy difference between the initial and final states can be evaluated analytically exactly, and is given by where the respective partition functions are given by the expressions The bounds on the integrals in the expressions above can be understood from the schematic representations of the potentials in Figs. 1 and 2. These integrals can be evaluated analytically, and give rise to the following expressions for the partition functions of the initial and final states, respectively, Equations (20) and (21) can be used along with Eqn. (17) to obtain the exact value of the free energy difference between the initial and final state for any choice of parameter values in the potentials U M (x) and U OT (x). Free energy differences for the particular choice of values listed in Table I as parameter sets 1, 2 and 3, are given in Table II. They are used to evaluate the accuracy of the free energy differences predicted by the Crooks and Jarzynski fluctuation theorems.

A. Code validation
In order to validate the predictions of the current algorithm, comparisons were carried out with the results of two earlier studies which demonstrated the Evans-Searles fluctuation theorems using experiments and simulations involving an optical trap. 37,49 The transient fluctuation theorem (TFT) of Evans and Searles 25,50,51 states that while the integrated form of the transient fluctuation theorem (ITFT) states that Here, Σ t is the dissipation function, which is a dimensionless measure of the total entropy production that occurs along the system's trajectory, over time t. It assumes different forms depending on the system under consideration. The TFT relates the probability of observing a trajectory with entropy production, Σ t = A, to the probability of observing a trajectory with the consumption of the same magnitude of entropy, Σ t = −A. On the other hand, the integrated version of the theorem specifies a relationship between the frequency of entropyconsuming trajectories to that of entropy-producing trajectories, with the average on the right hand side of Eqn. (23) carried out over only entropy-producing trajectories.
In the first study considered here, Wang et al. 37 examined the trajectory of a colloidal particle captured in an optical trap translated at a uniform velocity relative to the surrounding medium. They experimentally demonstrated the validity of the ITFT, and also carried out molecular dynamics simulations to show that the predictions of both the TFT and the ITFT were correct. In the second study, Carberry et al. 49 observed the timedependent relaxation of a colloidal particle subjected to a step change in the strength of a stationary optical trap. In this case, they were able to experimentally demonstrate the validity of both the TFT and the ITFT.
We have carried out Langevin simulations of these two previously studied applications of the Evans-Searles fluctuation theorems in order to ensure that our algorithm was implemented correctly. In both these examples, only a single optical trap is involved. As a consequence, the external force (in Eqn. (3)) on the colloidal particle due to the optical trap is given by, where k OT and x OT (t) assume different expressions in the two studies. As mentioned earlier, the dissipation function Σ t is also different in the two cases. The relevant expressions are listed below. Study 1 (Wang et al. 37 ): where F OT (x) is given by Eqn. (15). Study 2 (Carberry et al. 49 ): where H(t) is the Heaviside step function, and k 0 and k 1 are constants equal to the optical trap strength before  and after the step change, respectively.
The Langevin simulation of both these cases was carried out with 2 × 10 6 trajectories, using a time step of 10 −4 . In both cases, after an initial equilibration time of 10 4 time steps, the distribution of particle positions was checked to see if the respective equilibrium distribution functions were obeyed. In Study 1, after equilibration, the optical trap was translated with a constant velocity v OT = 0.5, from time t = 0 to t = 10, with a constant trap strength k OT = 1. In Study 2, after equilibration, the optical trap strength was changed discontinuously from k 0 = 1 to k 1 = 2 at time t = 0, and the simulation continued until t = 10. The position of the colloidal particle at time t = 0 is taken to be x(0). Figures 4 and 5 summarise the results of the validation studies.
In order to demonstrate the TFT a histogram of the values of the dissipation function Σ t at the end of the simulation was constructed over the 2 × 10 6 trajectories. If N i is the number of trajectories with dissipation function between Σ t,i ± ∆/2 (where ∆ = 0.1 is the size of the  The ITFT is demonstrated for the two studies in Figs. 5(a) and 5(b), respectively, by plotting the ratio of the number of entropy consuming trajectories (Σ t < 0) to the number of entropy producing (Σ t > 0) trajectories as a function of time, along with the time dependence of the entropy production averaged over the subset of 2 × 10 6 trajectories in which entropy is produced.

B. Crooks fluctuation theorem
Simulations were carried out with the three sets of parameter values listed in Table I for the membrane and optical trap potentials, with a time step size ∆t = 10 −3 . Rather than running the simulations for an initial equilibration period, the positions of the bead at time t = 0 were chosen such that they satisfied the known initial equilibrium distribution functions. Two kinds of simulations were carried out. The first kind, that generated forward trajectories, started at time t = 0 with the optical trap minimum at x OT = 0, followed by the trap minimum being translated with a uniform velocity v OT until it was located at x final OT at time t = t D . The set of optical trap velocities v OT = {0.01, 0.05, 0.1, 0.5, 1} was used. Note that t D depends on the value of v OT since the location x final OT is fixed and the same for all simulations. The second set of simulations, which generated reverse trajectories, started at time t = 0 with the optical trap minimum at x OT = x final OT , followed by the trap minimum being translated with the same set of velocities (but with opposite sign), until the minimum was located at x OT = 0 at time t = t D . Each simulation in the forward and reverse direction consisted of 10 5 trajecto-ries. Ten such simulations were carried out in each case. The work values obtained after each trajectory in both sets of forward and reverse simulations (calculated using Eqn. (16)), were sorted into bins of width equal to 0.01. The distributions of work values obtained in this manner are plotted in Fig. 6 for the various cases.
Panel A in Fig. 6 plots, for parameter set 1, the probability of work W being performed in the forward path (P F (W )) alongside the distribution of work values in the reverse path (P R (W )) for the various trap velocities v OT indicated in the figure legend. While the work is predominantly positive in the forward trajectories (with a positive mean value), the work is predominantly negative in the reverse trajectories (with a negative mean value). The widening of the distributions with increasing trap velocities is also apparent. As noted previously, in the limit of a quasistatic process (v OT → 0), P F (W ) → δ(W − ∆F ), and W F = ∆F . However, for increasing values of v OT , the mean value shifts towards the right with a wider range of work values, and with W F ≥ ∆F .
The usefulness of Crooks fluctuation theorem is best appreciated when P F (W ) is plotted alongside P R (−W ) as shown in panels B, C and D of Fig. 6. These three figure panels correspond to the three potential parameter sets listed in Table I, respectively. As noted before, according to Eqn. (7), the value of work W * at which P F (W * ) = P R (−W * ) is nothing but the equilibrium free energy difference. Consequently, ∆F is estimated from Fig. 6 by finding the point of intersection of the forward and reverse probability curves for each of the trap velocities, for the three sets of parameter values. The values of ∆F obtained in this way are listed in Table II, along with an estimate of the error in finding the point of intersection due to the relatively coarse interval used for binning the work values. The percentage relative error in the free energy predicted by the Crooks fluctuation theorem, defined by the expression is also listed in Table II. It is worth noting that the error in finding the point of intersection consistently increases with the trap velocities, but is roughly the same order of magnitude in all cases. On the other hand, the percentage relative error varies without a set pattern for the different values of v OT , depending on how close the predicted value is to the analytical value. Remarkably, for each parameter set, the intersection of the forward and reverse probability curves occurs at nearly identical values, with the error in the estimated free energy being at most 2% even for large trap velocities.
The increase in error with increasing trap velocity can be understood by considering panel B in Fig. 6. As the velocity increases, it causes the mean value of work to shift away from the free energy value, with a simultaneous increase in the standard derivation of the distribu-tion. As a result, the crossover occurs at the tails of the distributions, where errors are high and therefore require much larger populations to ensure adequate statistics. Figure 6 indicates that the velocities at which this could become an issue is sensitive to the choice of potential parameters. Parameter set 1 (panel B), where the optical trap strength was double that of the membrane, and the barrier height for detachment was much lower than that of re-attachment (see Fig. 3a), seems to have the most movement of the mean away from the exact free energy value. On the other hand parameter set 3 (panel D), where barrier heights are of O(k B T ) (see Fig. 3c), seems to be the least affected by increased velocity.

C. Jarzynski equality
The form of the Jarzynski equality given by Eqn. (8) corresponds to switching the system from an initial equilibrium state with λ = λ 0 to a final state with λ = λ f . When the system is switched from an initial equilibrium state with λ = λ f to a final state with λ = λ 0 , the Jarzynski equality takes the form, 52 where the subscript 'R' on the ensemble average on the left hand side indicates an average over reverse trajectories, and the change in free energy is still defined by The sets of forward and reverse simulations carried out to demonstrate the Crooks fluctuation theorem can also be used to examine the usefulness of the Jarzynski equality. The ensemble averages on the left hand sides of Eqns. (8) and (26) were calculated using the values of work accumulated at the end of each of the 10 5 trajectories corresponding to a particular simulation. The sets of forward and reverse simulations were repeated ten times each, so that we obtain ten estimates for the equilibrium free energy in each case, and the errors can be estimated. The mean of these 10 values, and the standard error in these mean values are displayed in Fig. 7 for all the cases considered here. Parameter sets 1, 2, and 3 are shown in rows 1, 2, and 3 respectively, with the left hand column showing results for the forward trajectories whilst the right hand column shows results for reverse trajectories. The mean value of ∆F and the standard error in the mean are also compared with exact analytical values in Table II for simulations carried out in the forward direction. Note that the percentage relative error reported in the Table is calculated using Eqn. (25) with the mean value of ∆F .
A feature of all approaches for determining free energy differences using ensemble averages, of which the Jarzynski equality is no exception, is their limitation due to sample size. As argued by Jarzynski, 28 for systems where the spread in the distributions P F (W ) and P R (W ) is large, the function exp(−W ) varies signifi-TABLE II. Comparison of equilibrium free energies calculated with the Crooks fluctuation theorem, the Jarzynski equality, and from a sum over the first six terms of the cumulant expansion, with exact analytical values, for the various trap velocities. The three sets of values for the membrane and optical trap potential parameters are given in Table I cantly over many standard deviations about the mean value of work. As a result, the numerically determined average exp(−W ) can be dominated by work values that are by their very nature statistically rare. Therefore an unreasonable number of measurements of the work would be required to get an accurate result. This results in a practical restriction on the rates at which the system can be switched between λ 0 and λ f . As can be seen from Fig. 7 and Table II, the accuracy in the estimation of the free energy decreases with the trap velocity in all cases. A comparison of the relative errors in the free energies predicted by the Crooks fluctuation theorem and the Jarzynski equality (in the case of forward trajectories) in Table II shows that they are roughly similar in magnitude for the various cases. As noted earlier, there is a reduction in accuracy with increasing trap velocity, which appears to be magnified when either one or both the potential well depths are high compared to k B T , which is the case for parameter sets 1 and 2 (displayed in Fig. 3). The dependence of the error on well depth is studied shortly below.
For slow rates of switching between λ 0 and λ f , the distributions P F (W ) and P R (W ) are expected to be approximately Gaussian. 52 In this case, retaining only the first two terms in the cumulant expansion for exp(−W ) (which is discussed in greater detail in the section below), one can write, 52 where σ 2 F and σ 2 R are the variances of the work distributions P F (W ) and P R (W ), respectively. Defining the mean dissipated work W d as the difference between the mean actual work of the process and the reversible work (which is equal to the equilibrium free energy), we can estimate the departure from the Gaussian approximation by evaluating the error estimates E F and E R defined by, The values of mean actual work, variances, mean dissi-  pated work and error estimates, for membrane and optical trap potential parameters corresponding to Set 1, are displayed in Table III for both the forward and reverse paths. Clearly, the Gaussian approximation leads to an error of less than 9% up to trap velocities v OT = 0.5. Interestingly, the variances of P F (W ) and P R (W ) and the mean dissipated work in the forward and reverse paths are roughly equal in magnitude for identical velocities in the forward and reverse paths. For distributions that are not Gaussian, the exponential average in Jarzynski's equality can be expanded in terms of cumulants, 52 and the convergence of ∆F can be studied as a function of the various potential parameters, as discussed in the section below. It is worth noting that it is also possible to obtain estimates for the free energy that are accurate to a higher order in the cumulant expansion than the Gaussian approximation by suitably combining the mean work and variance in the forward and reverse paths. 52

D. Cumulant expansion for the free energy of binding
The average of the exponential of work on the left hand sides of Eqns. (8) and (26) in Jarzynski's equality can be expanded in terms of cumulants. 52 In the case of forward paths, this leads to the following expression for the free energy change: where Here, the cumulants C n are defined by the expressions . . .
with µ n being the central moments of P F (W ), The recursive relationship between the cumulants and central moments in Eqn. (31) has been given by Smith. 53 In the case of reverse paths, the cumulant expansion on the right hand side of Eqn. (30) leads to the free energy change −∆F = F λ0 − F λ f , with µ n in the expressions for C n being the central moments of P R (W ). An analysis of the simulation results for the forward and reverse paths in terms of the cumulant expansion is displayed in Fig. 8, where the difference between the values of ∆F k (which represent the approximate estimate of the free energy change given by k terms of the cumulant expansion) and the analytical value ∆F anal , is plotted against the trap velocities v OT (for values of k in the range 2 ≤ k ≤ 6). Additionally, the particular values obtained for ∆F 6 in the case of forward trajectories, and the relative error compared to the exact values are listed in Table II. As expected, at low trap velocities where the system approaches a quasistatic process, the work distribution approaches a Gaussian, and quite accurate results are obtained with two cumulants. However as the trap velocity increases, higher cumulant numbers are required until, for v OT = 1, even at cumulant numbers of 6 the system has still not converged.
An alternative representation of the cumulant expansion data is given in Fig. 9, where ∆F k −∆F anal is plotted as a function of k (2 ≤ k ≤ 6), at the lowest and highest trap velocities (v OT = 0.01 and v OT = 1.0), for parameter values corresponding to set 3. Since the cumulant expansion is an approximation for the left hand sides of Eqns. (8) and (26), we expect that the free energy difference ∆F k should converge to the free energy difference predicted by Jarzynski's equality ∆F Jarzynski , for sufficiently large values of k. This can be seen to be clearly the case for v OT = 0.01, for both the forward and reverse trajectories, from the top row in Fig. 9, where the solid line corresponds to the difference ∆F Jarzynski − ∆F anal . The scale of the y-axis in both the subfigures in the bottom row of Fig. 9 (corresponding to v OT = 1.0) makes it difficult to distinguish ∆F Jarzynski −∆F anal from 0. While the values of ∆F k − ∆F anal appear to be getting smaller with increasing k, there are still large changes in ∆F k with increasing k, and convergence has not occurred by k = 6, as was observed previously at this value of trap velocity in Fig. 8.
The cumulant expansion can also be used to examine the influence of well depth. In order to do so, simulations in the forward direction were carried out for 10 6 trajectories with time step ∆t = 10 −4 , for trap velocities v OT = {0.01, 0.05, 0.1, 0.5, 1}. In all cases, the final location of the trap potential minimum was x final OT = 6. The membrane potential depth was held fixed at M = 4, whilst a parameter sweep from 1 to 8 was carried out for the optical trap potential depth, OT . The trap strengths k M and k OT for both the membrane and the optical trap potentials were held constant at a value of two. Results of the cumulant analysis are plotted in Fig. 10 for the difference ∆F k − ∆F anal , as a function of trap velocity, at the various values of k, with each subfigure representing a different value of OT . Since the exact analytical value ∆F anal is different for each value of trap well depth, the values are given in the caption to Fig. 10. The cumulant analysis suggests that convergence occurs quickly at the low velocities and becomes poorer and poorer at higher velocities. It is also evident that increasing optical trap well depth significantly increases the error in the estimate of the free energy for a given value of the number of terms k in the cumulant expansion (note the different scales of the y-axes in the different subfigures of Fig. 10).

E. Probabilities of attachment and detachment via umbrella sampling
An important quantity that is frequently the focus of experiments on cell adhesion is the probability of adhesion. Measurements of the adhesion probability are often used to determine the kinetics of the adhesion process through the calculation of on and off-rates of binding etc. The experiments, which typically monitor whether a binding event occurs or not when ligand and receptor bearing surfaces are brought into contact, are by their very nature carried out at finite rates. As a result, a true measure of the equilibrium probability of binding is difficult to obtain. In this context, the method of nonequilibrium umbrella sampling 31,36,54,55 provides a means of determining the equilibrium binding probability from non-equilibrium measurements. Here, we demonstrate how non-equilibrium umbrella sampling can be used to find, at the end of the unbinding experiment, the probability of either the bead being attached to the cell, or being detached from it and held in the optical trap.
At the end of the computer experiment, when t = t D and the optical trap minimum is located at x final OT , it makes sense to sub-divide the x axis into three intervals (cf. Fig. 2 (c)): First, there is the interval −∞ < x(t D ) ≤ x ub M , which we define as the set of states that correspond to the bead still being attached to the cell (membrane). The second interval is x ub M < x(t D ) < x lb OT , where the potential is flat, and which we define as corresponding to an intermediate state of the bead. Finally, the interval x lb OT ≤ x(t D ) < ∞ corresponds, according to our definition, to the detached (or optically trapped) state of the bead. In what follows, we will focus on the equilibrium probabilities for the attached state and the detached state; the probability for the intermediate state then follows trivially by subtracting the sum of these values from one.
To formalise these definitions, it is useful to introduce the indicator functions χ A and χ D , The equilibrium probabilities for the attached and the detached states are then simply the Boltzmann averages of χ A and χ D , respectively. Here of course the Boltzmann distribution corresponding to the final potential profile (λ = λ f ) must be used: For the choice of potentials in the present work, it is straightforward to determine these values analytically. Using arguments along the lines of those in section II E for the analytical determination of free energy differences, we can show that where the quantities Z A and Z D in the equations above are given by and Z(x OT = x final OT ) is given by Eqn. (21). These expressions are useful to evaluate the degree of success of the non-equilibrium umbrella sampling technique in determining the equilibrium probabilities p A and p D from the nonequilibrium computer experiment. This latter analysis is done as follows: We denote the total number of detachment simulations with N T . Similarly, N A denotes the number of runs where the bead ends up in the attached state (χ A (x(t = t D )) = 1). Analogously, N D is the number of runs where the bead is finally detached. If p λ f neq (x(t D )) is the non-equilibrium distribution of bead positions at the final time t D , then the non-equilibrium probabilities of attachment and detachment, defined by the following expressions, are easily estimated by simulations from the ratios N A /N T and N D /N T , respectively: We now use the technique of non-equilibrium umbrella sampling to obtain the equilibrium probabilities from the non-equilibrium simulations. Let us outline this method in general terms: For an observable B that has been sampled by a nonequilibrium (computer) experiment, i.e., using the probability distribution p λ f neq (x), we simply have to multiply each data point with the ratio p λ f eq (x)/p λ f neq (x) such that the data point is given the weight p λ f eq (x) rather than p λ f neq (x). It can be shown 31,54,55 that the ratio p λ f eq (x)/p λ f neq (x) is nothing but the factor e −W , except for normalisation. Therefore, the data need to be reweighted according to the formula Application of this general formula to our observables (χ A , χ D ) yields Simulation data generated previously for examining the influence of well depth in Sec. III D has been used here for evaluating the usefulness of nonequilibrium umbrella sampling, for trap velocities v OT = {0.01, 0.05, 0.1, 0.5, 1}. The potential parameters used in the simulations are as given in the caption to Fig. 10, along with x final OT = 6.
The symbols in Figs. 11 (a) and (b) are the nonequilibrium probabilities of attachment and detachment p neq A and p neq D , determined from Eqns. (41) and (42) for various trap velocities, while the symbols in Figs. 11 (c) and (d) are the equilibrium probabilities p A and p D , determined by applying the umbrella sampling procedure as expressed in Eqns. (44) and (45). Error bars estimated from the ten repeated simulations are smaller than the symbol size in Figs. 11 (c) and (d). The curves in the subfigures of Fig. 11 represent the analytical equilibrium probabilities p anal A and p anal D (as appropriate), calculated from Eqns. (37) and (38), respectively.
As expected, Figs. 11 (a) and (b) indicate that the non-equilibrium probabilities are nearly identical to the equilibrium probabilities at low trap velocities, but deviate from the latter more and more as the trap velocity increases. Interestingly, the greatest departure occurs for membrane and optical trap potential well depths that are roughly equal in magnitude. Not surprisingly, the probability of detachment is greatest for the largest optical trap well depth, while the likelihood of remaining in the membrane potential is high at low trap well depths. For nearly all the trap velocities, except perhaps at v OT = 1 (for roughly equal trap strengths), application of umbrella sampling recovers the equilibrium probabilities nearly perfectly.

IV. CONCLUSIONS
A simple model for the detachment of a ligand coated bead with the help of an optical tweezer, from receptors on the surface of a cell to which it is bound, has been used to examine if fluctuation theorems are useful in de-termining equilibrium free energies, which in turn provide information about the binding energetics. By using truncated harmonic potentials to represent the stationary cell membrane and the moving optical trap, and a Langevin equation to model the stochastic motion of the bead in these potentials, the distribution of work performed in driving the system from an initial equilibrium state to a final non-equilibrium state (at various finite rates) has been calculated by carrying out repeated simulations of the Langevin equation in the forward and reverse directions. The former corresponds to the membrane and trap potentials being superposed at time t = 0, followed by the optical trap being translated uniformly until the two potentials are sufficiently apart at the final time t = t D . The latter refers to the opposite situation.
The calculation of work distributions enables the determination of the equilibrium free energy change between the initial and final states of the system, using both the Crooks fluctuation theorem and the Jarzynsky equality. The simplicity of the model also permits a straight forward determination of the exact free energy change by analytical means. It is found that both fluctuation theorems lead to excellent predictions provided the rate of switching from the initial to the final state is sufficiently slow. For relatively rapid rates of trap translation, sampling problems (for the given sample size) lead to a decrease in accuracy. The reduction in accuracy is discussed both in terms of a Gaussian approximation for the work distributions, and a cumulant expansion for the average of the exponential of work.
The method of non-equilibrium umbrella sampling has been used to determine the equilibrium probability that, after translating the trap from its initial to its final location, the bead and cell are still attached (i.e., the bead lies only within the range of influence of the membrane potential), and the equilibrium probability that the bead and cell have been detached (the bead lies only within the range of influence of the optical trap potential), for a range of different values of the optical trap well depth. It is seen that by appropriately analysing the non-equilibrium simulation data, accurate estimates of the equilibrium probabilities of attachment and detachment can be found for all but the highest rates of trap translation.
In conclusion, although a very simple model has been used, the present work demonstrates that nonequilibrium fluctuation theorems can be applied without significant statistical problems to binding/unbinding experiments carried out with optical trap velocities that are realizable under experimental conditions. Combined with the theoretical procedure outlined in section II C, they could provide a reliable means of extracting unknown membrane potentials (see Eqn. (13)).