Free Energy Reconstruction from Biased Experiments (FERBE) version 0.2

The calculation of free energy differences is a major focus in chemistry. Chemists are interested in free energy differences between thermodynamic states. For example, we would like to know how likely it is for a ligand to bind a particular receptor. We are also interested the free energy surface along a system order parameter, such as a protein dihedral angle or the length of a molecule; this is known as the potential of mean force (PMF). A difficulty with rigorous free energy calculations, however, is that they require thorough search of configuration space. This often includes the sampling of rare events. Fortunately, it is possible to improve the sampling of these events by modifying the search space - performing a biased experiment.

I've written several MATLAB programs to calculate free energy differences from different types of biased experiments.

1. Equilibrium Umbrella Sampling -> Weighted Histogram Analysis Method (WHAM). A common way to calculate the potential of mean force is to run several simulations with a harmonic bias that favors different values of a particular order parameter. After removing the effect of each bias, one is left with multiple short segments of the PMF. WHAM has been shown to be the optimal way to combine these different segments to estimate the overall PMF [1-3].
2. Unidirectional Single-Molecule Pulling -> the Hummer-Szabo method. A single-molecule pulling experiment is performed by adding a time-dependent bias to the system, i.e. using single-molecule force spectroscopy or a steered molecular dynamics simulation. The free energy difference between these states (defined by the Hamiltonian at a particular time) can be calculated from the path-ensemble average of the exponentially weighted work using Jarzynski's equality [4-5]. Hummer and Szabo's method [6-7] can be used to calculate the potential of mean force along the pulling coordinate. I've extended Hummer and Szabo's method to PMF reconstruction from pulling experiments performed under different protocols [8] and for multidimensional PMF reconstruction [9].
3. Bidirectional Single-Molecule Pulling -> the Bennett Acceptance Ratio (BAR) and the Conjugate Twin method. The Bennett Acceptance Ratio [10] is a classic method to calculate the free energy difference between two states based on the potential energy differences of individual configurations sampled from both states. More recently, Crooks extended it to nonequilibrium processes, replacing the potential energy difference with the work required to switch between two states [11]. The conjugate twin method [12] is an optimized way to calculate path-ensemble averages from bidirectional pulling. It can be used to calculate the free energy difference between every state in the switching process and estimate the PMF. I also have written an implementation of another bidirectional method developed by Chelii and coworkers [13].

The source codes for these programs are available for free under a general public license (GPL). I recognize that as a proprietary commerical software, MATLAB has the disadvantage of being less widely available than languages with open-source compilers, such as FORTRAN and C, but I have chosen to use it as a platform for algorithmic development because its extensive library of mathematical and matrix manipulation routines. These codes should at least serve a pedagogical purpose to computational scientists without MATLAB. All of the files below are available in FERBE.zip.

If you use these in a scientific publication I ask that you cite the appropriate work. Feel free to send me comments and suggestions.

About the versions...

version 0.2 (4/23/2008) - Analysis of bidirectional data added.
version 0.1 alpha (3/22/2007) - Original version. Includes WHAM and unidirectional pulling analysis.

1. Equilibrium Umbrella Sampling

This implementation of the WHAM method [1-3] refers to equations in Roux's paper [3]. I referred to Alan Grossfield's C++ implementation while writing it. I used it my multidimensional PMF paper [9] to compare umbrella sampling and nonequilibrium pulling experiments.

Simulation

potential_twoD.m - Calculate the potential and its partial derivatives.
umbrella_twoD.m - Run umbrella sampling along a two-dimensional surface.

Analysis

probability_twoD.m - Calculate the probability density as a function of X and Y.
hist_twoD.m - Generate a two dimensional matrix of counts.
WHAM.m - PMF reconstruction using WHAM.
WHAM_twoD.m - Two-dimensional PMF reconstruction using WHAM.
plot_WHAM_twoD.m - Creates a figure similar to Figure 1 in reference 9.

To generate Figure 1 in reference 9, enter the following commands in MATLAB:

a = 0;
[Xs,Ys,Bxs] = umbrella_twoD(a);
plot_WHAM_twoD(Xs,Ys,Bxs,a);

2. Unidirectional Single-Molecule Pulling

These codes implement Hummer and Szabo's free energy reconstruction method [6-7], as well as several extensions [8-9].

Simulation

potential_twoD.m - Calculate the potential and its partial derivatives.
probability_twoD.m - Calculate the probability density as a function of X and Y.
slowpull_DBwell.m - Slow pulling protocol.
mspull_DBwell.m - Dynamic force spectroscopy.
randpull_DBwell.m - Random biasing protocol.

Analysis

reconstruct_time.m - PMF reconstruction using time slices [6].
reconstruct_bias.m - PMF reconstruction using bias slicing [8].
reconstruct_time_twoD.m - Two-dimensional PMF reconstruction using time slicing [9].
plot_pullreconstruct_twoD.m - Creates a figure similar to Figures 2 and 4 in reference 9.

To generate Figure 1 in reference [8] using an reconstruction algorithm that assumes a fixed biasing protocol, enter the following commands in MATLAB:

[Xs,Bxs,Wt,Fs] = slowpull_DBwell;
[HSJE,meanW,FD,xspace] = reconstruct_time(Xs,Bxs(1,:)',Wt,2.0);
plot(xspace,[HSJE;FD;meanW]);

Notice that the input for reconstruct_time includes Bxs(1,:)'. This is because the input for the pulling protocol of reconstruct_time is a Nx1 vector, whereas the output from slowpull_DBwell is a matrix. The reason it is designed to be a matrix is to make it compatible with reconstruct_bias, in which the pulling protocol for each run can be different.

To generate the same figure using a reconstruction algorithm that allows you to input the biasing protocol for each trajectory, enter the following commands in MATLAB:

[Xs,Bxs,Wt,Fs] = slowpull_DBwell;
[HSJE,meanW,FD,xspace] = reconstruct_bias(Xs,Bxs,Wt,2.0);
plot(xspace,[HSJE;FD;meanW]);

To generate Figure 2 in reference [9], enter the following commands in MATLAB:

a = 0.3;
[Xs,Ys,Bxs,Wt] = pull_twoD(a);
plot_pullreconstruct_twoD(Xs,Ys,Bxs,Wt,a);

3. Bidirectional Single-Molecule Pulling

These codes implement the Bennett Acceptance Ratio, the Conjugate Twin Estimator, and Chelli et al.'s Maximum Likelihood Method.

Simulation

potential.m - Calculate the potential and its partial derivatives.
pull.m - Pulling with Brownian dynamics.
pull_vv.m - Pulling with velocity verlet.
test_pullFR.m - Performs the test simulations in reference [12].

Analysis

histw.m - Calculates a weighted histogram.
BAR.m and BAR_sc.m - The Bennett Acceptance Ratio.
Chelli.m and Chelli_sc.m - Chelli et al.'s method.
pmf_pullF.m - PMF reconstruction using time slices of unidirectional pulling.
pmf_pullFR.m - PMF reconstruction using time slices of bidirectional pulling.

To generate Figures 1 and 2 in reference [12], run test_pullFR.m

Citations

[1] A. Ferrenberg and R. Swendsen. Optimized Monte Carlo data analysis. Physical Review Letters 63, 1195 (1989).
[2] S. Kumar, J.M. Rosenberg, D. Bouzida, R.H. Swendsen, and P.A. Kollman. Multidimensional free-energy calculations using the weighted histogram analysis method. Journal of Computational Chemistry 16(11):1339-1350 (1995).
[3] B. Roux. The calculation of the potential of mean force using computer simulations. Computer Physics Communications 91, 275 (1995).
[4] C. Jarzynski. Nonequilibrium Equality for Free Energy Differences. Physical Review Letters 78, 2690 (1997).
[5] C. Jarzynski. Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach. Physical Review E 56, 5018-5035 (1997).
[6] G. Hummer and A. Szabo. Free energy reconstruction from nonequilibrium single-molecule pulling experiments. Proceedings of the National Academy of Sciences U.S.A. 98, 3658-3661 (2001).
[7] G. Hummer and A. Szabo. Free Energy Surfaces from Single-Molecule Force Spectroscopy. Accounts of Chemical Research 38(7), 504-513 (2005).
[8] D. Minh. Free-energy reconstruction from experiments performed under different biasing programs. Physical Review E 74, 061120 (2006).
[9] D. Minh. Multidimensional Potentials of Mean Force from Biased Experiments Along a Single Coordinate. Journal of Physical Chemistry B 111(16), 4137-4140 (2007).
[10] C. Bennett. Efficient Estimation of Free Energy Differences from Monte Carlo Data. Journal of Computational Physics 22, 245-268 (1976).
[11] G. Crooks. Path-ensemble averages in systems driven far from equilibrium. Physical Review E 61, 2361-2366 (2000).
[12] D. Minh and A. Adib. Optimized free energies from bidirectional single-molecule force spectroscopy. Physical Review Letters 100, 180602 (2008).
[13] R. Chelli, S. Marsili, and P. Procacci. Calculation of the potential of mean force from nonequilibrium measurements via maximum likelihood estimators. Physical Review E 77, 031104 (2008).