Two dimensional Replica Exchange

Table of Contents

1 Introduction

This example shows how to use SWHAM to analyze the data generated by two dimensional (Hamiltonian and temperature) replica exchange simulations.

1.1 System

We study the binding affinity of a guest molecule (heptanoate) to a host molecule (\(\beta\)-cyclodextrin) in implicit solvent (OPLA-AA/AGBNP2) at \(15\) temperatures. Here the binding energy distribution analysis method (BEDAM) was applied to obtain the binding free energy and binding energy distributions of this complex at different temperatures.

BEDAM is a free energy method based on the Hamiltonian replica exchange (RE) algorithm in which the interaction between the ligand and acceptor is scaled by the factor \(\lambda\) changing gradually from \(0\) to \(1\). Suppose there are \(M\) parallel simulations in BEDAM, the Hamiltonian (potential) function of the \(i^{th}\) \(\lambda\)-state is

\begin{equation} H[i] = V_0 + \lambda[i] u\,, \end{equation}

where \(V_0\) is the effective potential energy of the complex without the direct and solvent-mediated ligand-receptor interactions, and \(u\) is the binding energy.

The raw data used in this example were generated by \(15\) separated BEDAM simulations at temperatures \(200K\), \(206K\), \(212K\), \(218K\), \(225K\), \(231K\), \(238K\), \(245K\), \(252K\), \(260K\), \(267K\), \(275K\), \(283K\), \(291K\), \(300K\). Each BEDAM simulation was run by using \(16\) \(\lambda\)-states. The chosen \(\lambda\) values are \(0.0\), \(0.001\), \(0.002\), \(0.004\), \(0.01\), \(0.04\), \(0.07\), \(0.1\), \(0.2\), \(0.4\), \(0.6\), \(0.7\), \(0.8\), \(0.9\), \(0.95\), \(1.0\). All the BEDAM simulations lasted \(72\ ns\). And the values of \(V_0\) and \(u\) were recorded very \(0.5\ ps\). Therefore, in this system there are totally \(16 \times 15 = 240\) states, and each state has \(144,000\) data points. Although there are no exchanges between replicas at different temperatures, the procedure described in this tutorial can be applied to two dimensional (Hamiltonian and temperature) replica exchange simulations without any alteration. The goal of this practice is to obtain the best estimates of the binding affinity at \(200\ K\), which is the most difficult one to converge because of the low temperature.

2 Files

Please download all these files into a same directory.

2.1 Data Files

  • data.tgz: Unzipping this file generates a directory called "Data" which contains \(240\) data files. Each data file corresponds to one \(\lambda\)-state and contains two columns of data. the first column is the effective potential energy \(V_0\) of each observation; the second column is the binding energy \(u\). Both are in units of "\(kcal/mol\)".
  • stateinfo.dat: This file contains the information of \(\lambda\)-states. There are \(240\) lines in this file. Each line corresponds to a \(\lambda\)-state. The first column is the temperature of each \(\lambda\)-state, and the second column is the \(\lambda\) value of each \(\lambda\)-state.

2.2 Python Scripts:

  • Histogram.py: Python script constructs a one dimensional histogram.

2.3 Others

3 RE-SWHAM Analysis

Please read the manual of the SWHAM program before practicing.

  1. First we construct input file for RE-SWHAM by running

    for i in {01..15}
    do
    for j in {01..16}
    do
    echo "# state $i $j" >> input4swham.dat
    cat Data/T_${i}_lambda_$j.dat >> input4swham.dat
    done
    done
    

    This shell script simply combines the raw data of each \(\lambda\)-state into one file "input4swham.dat". The data points of different \(\lambda\)-states are separated by a line beginning with "#".

  2. Run RE-SWHAM

    /path of SWHAM/swham.o -d input4swham.dat -f stateinfo.dat -u 3 -i 1 -o 1 -q 100 -n 1000 -x 2000 -s 1-16 -p 2 -k > reswham.out
    

    The explanations for options are:

    -d input4swham.dat read the file "input4swham.dat" for the information of observations
    -f stateinfo.dat read the file "stateinfo.dat" for the information of \(\lambda\)-states
    -u 3 use the type 3 potential function implemented in the SWHAM program
    -i 1 -o 1 the input and out data are in units of \(kcal/mol\)
    -q 100 -n 1000 -k equilibrate the system for 100,000 cycles first; then run swham for 1,000,000 cycles
    -x 2000 2000 exchange attempts per RE cycle to reach the infinite swapping limit
    -s 1-16 -p 2 print out the second property (binding energy) of each observation sampled at the first \(16\) \(\lambda\)-states (\(200\ K\))
  3. Next we plot the RE-SWHAM estimates of the equilibrium distribution of the interested \(\lambda\)-state (\(T=200\ K\); \(\lambda=1.0\)), and compare the RE-SWHAM estimate with the raw data. To construct the histogram based on the RE-SWHAM estimates, run the command

    ./Histogram.py -f reswham.out -i 16 -b 100 -n -45 -x -10 -a > bindingE_swham.hist
    
    -f reswham.out -i 16 read the sixteenth column of the file "reswham.out" as the data points
    -b 100 totally 100 bins in the histogram
    -n -45 -x -10 the minimum of the histogram is \(-45\ kcal/mol\) and the maximum is \(-10\ kcal/mol\)
    -a the probability of the data points \(<\ -45.0\ kcal/mol\) and \(> 5.0\ kcal/mol\) are included in the first and the last bin respectively.

    To construct the histogram based on the raw data sampled at the \(\lambda\)-state (\(T=200\ K\); \(\lambda=1.0\)) by BEDAM, run the command

    ./Histogram.py -f Data/T_01_lambda_16.dat -i 2 -b 100 -n -45 -x -10 -a > bindingE_rawdata.hist 
    

    Then we compare these two histograms. This can be done by running gnuplot (>v5.0) and loading the file "bindingE.gnu". The resulting plot should look like this

    bindingE200.jpg

    As can be seen, (1D) replica exchange had not converged at \(200\ K\) after \(72\ ns\) simulations. Applying RE-SWHAM to this unconverged data set led to significant improvements for the binding energy distribution at \(200\ K\) by using information from the simulations at higher temperatures.1

  4. The output of the RE-SWHAM analysis, "reswham.out", contains the equilibrium distributions of binding energies of the \(16\) \(\lambda\)-states at \(200\ K\). It can be used to estimate the binding free energy at \(200\ K\) by using UWHAM. The UWHAM estimate of the binding free energy at \(200\ K\) is about \(-6.3\ kcal/mol\).2 We can separate the binding energies of each \(\lambda\)-state by running the command lines

    mkdir -p UWHAM/Data
    for i in {0..15}
    do
    let j=i+1
    filename=$(printf 'lambda_%02d.dat' $j)
    perl -wnal -e " /#/ or printf(\"%12d %12g\\n\", \$., \$F[$i]) " reswham.out > UWHAM/Data/$filename
    done
    

    Then the input file for UWHAM can be generated by following the steps explained in the example "Hamiltonian Replica Exchange".

  5. Another option to obtain the estimate for the binding free energy at \(200\ K\) is to run ST-SWHAM using the local jump algorithm

    /path of SWHAM/swham.o -d input4swham.dat -f stateinfo.dat -u 3 -i 1 -o 1 -q 50 -n 200 -g 50 -m -s 1-16 -x 12
    

    The explanations for options are:

    -d input4swham.dat read the file "input4swham.dat" for the information of observations
    -f stateinfo.dat read the file "stateinfo.dat" for the information of \(\lambda\)-states
    -u 3 use the type 3 potential function implemented in the SWHAM program
    -i 1 -o 1 the input and out data are in units of \(kcal/mol\)
    -g 50 -m run ST-SWHAM and print out the estimates of free energies every 50,000,000 cycles
    -q 50 -n 200 -m equilibrate the system for 50,000,000 cycles first; then run swham for 200,000,000 cycles
    -x 12 12 local jump attempts per ST cycle
    -s 1-16 print out the free energy estimates for the first \(16\) \(\lambda\)-states (\(200\ K\))

    Does the ST-SWHAM estimate for the binding free energy at \(200\ K\) agree with the UWHAM estimate?

Footnotes:

1

A Stochastic Solution to the Unbinned WHAM Equations (link)

2

Locally Weighted Histogram Analysis and Stochastic Solution for Large-Scale Multi-State Free Energy Estimation (link)

Author: Bin Zhang

Created: 2018-05-03 Thu 16:08

Validate