(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2006

AMBER ADVANCED TUTORIALS
TUTORIAL 3 - SECTION 3

Perl Script mm_pbsa.pl

By Ross Walker & Thomas Steinbrecher

3) Calculate the binding free energy and analyse the results. (All calculations performed using mm_pbsa.pl here)

We now need to extract snapshots (without the water) from our production runs for use in the MM-PBSA calculation. The mm_pbsa.pl script (in the $AMBERHOME/bin directory) automates this extraction process for us. Note if you don't have gzcat installed you will need to unzip the trajectory files before running mm_pbsa. We have to provide the following input file (see the linked file for an explanation of each of the various terms):

extract_coords.mmpbsa

@GENERAL
PREFIX snapshot PATH ./ COMPLEX 1 RECEPTOR 1 LIGAND 1 COMPT ./ras-raf.prmtop RECPT ./ras.prmtop LIGPT ./raf.prmtop GC 1 AS 0 DC 0 MM 0 GB 0 PB 0 MS 0 NM 0 @MAKECRD BOX YES NTOTAL 42193 NSTART 1 NSTOP 200 NFREQ 1 NUMBER_LIG_GROUPS 1 LSTART 2622 LSTOP 3862 NUMBER_REC_GROUPS 1 RSTART 1 RSTOP 2621 @TRAJECTORY TRAJECTORY ./prod1.mdcrd TRAJECTORY ./prod2.mdcrd TRAJECTORY ./prod3.mdcrd TRAJECTORY ./prod4.mdcrd @PROGRAMS

This just specifies which atoms are part of the receptor, ligand and complex as well as specifying the prmtop files corresponding to the unsolvated structures, the total number of snapshots in the trajectores, the stride length and the names of the trajectory files. We have specified that each of the output files should be written with the prefix 'snapshot'. In this setup we have defined RAS to be the receptor and RAF to be the ligand. This is purely a naming convention.

$AMBERHOME/bin/mm_pbsa.pl extract_coords.mmpbsa > extract_coords.log

This will take a few minutes to run. Here are the relevant output files: extract_coords.tar.gz (14 MB)

You should check the log file for any errors. Also make sure the box sizes look reasonable. If there are any strange box sizes this is likely a problem with the atom numbers you selected or a corruption in one of your trajectory files.

Starting with the snapshots we just extracted, we will now calculate the interaction energy and solvation free energy for complex, receptor, and ligand and average the results to obtain an estimate of the binding free energy. Please note that we will not perform a calculation of the entropy contribution to binding in this tutorial and so strictly speaking our result will not be a true free energy but could be used to compare against similar systems. E.g. one could carry out an analysis of the effect of amino acid point mutations along the binding interface. A common approach to this is referred to as alanine scanning.

For demonstration purposes we will carry out the binding energy calculation using both the MM_PBSA method and the MM_GBSA method.  This is accomplished with the following input file for mm_pbsa.pl (see the linked file for comments):

binding_energy.mmpbsa

VERBOSE               0
PARALLEL              0
PREFIX                snapshot
PATH                  ./
START                 1
STOP                  200
OFFSET                1
COMPLEX               1
RECEPTOR              1
LIGAND                1
COMPT                 ./ras-raf.prmtop
RECPT                 ./ras.prmtop
LIGPT                 ./raf.prmtop
GC                    0
AS                    0
DC                    0
MM                    1
GB                    1
PB                    1
MS                    1
NM                    0
@PB
PROC                  2
REFE                  0
INDI                  1.0
EXDI                  80.0
SCALE                 2
LINIT                 1000
ISTRNG                0.0
RADIOPT               0
ARCRES                0.0625
INP                   1
SURFTEN               0.005
SURFOFF               0.00
IVCAP                 0
CUTCAP                -1.0
XCAP                  0.0
YCAP                  0.0
ZCAP                  0.0
@MM
DIELC                 1.0
@GB
IGB                   2
GBSA                  1
SALTCON               0.00
EXTDIEL               80.0
INTDIEL               1.0
SURFTEN               0.005
SURFOFF               0.00
@MS
PROBE                 0.0
@PROGRAMS

The various portions of this input file specify which calculations to run. On which files to run them and any special parameters necessary to calculate the different contributions to the binding free energy. If you open the linked file above you will find explanations of the different terms. Values for the different parameters are determined based on empirical data and are subject to on-going research. Please find here a list of the currently recommended settings for the calculation of the non-polar solvation free energy contribution. Example input files for binding free energy calculations with common parameter settings can be found in the $AMBERHOME/AmberTools/src/mm_pbsa/Examples/TEMPLATE_INPUT_SCRIPTS directory. Refer to publications in the field for more information. Note early versions of AMBER required an external poisson Boltzmann solver such as Delphi but as of AMBER 8 the internal PBSA program can be used. This program provides a significant speedup and easier integration into the AMBER framework. You can run the calculation with:

$AMBERHOME/bin/mm_pbsa.pl binding_energy.mmpbsa > binding_energy.log

You can watch the progress of your calculation by:

tail -f binding_energy.log

This calculation will take around 2 hours to run (P4 3.2 GHz). The PBSA part of the calculation on each of the 600 snapshots is what takes the majority of the time. The GBSA part of the calculation is done in seconds. To speedup the calculation you can also request parallel execution of mm_pbsa analyses for more than one snapshot at a time by specifying the number of available processors under PARALLEL.

When finished you should find the following output files: binding_energy.log, snapshot_statistics.out, snapshot_com.all.out, snapshot_rec.all.out, snapshot_lig.all.out

The log file just shows if the calculation completed successfully. The all.out files give the individual energy contributions for each of the snapshots for each of the species while the statistics.out file contains the final averaged binding free energy results:

snapshot_statistics.out

#                  COMPLEX                RECEPTOR                  LIGAND        
#          ----------------------- ----------------------- -----------------------
#                  MEAN        STD         MEAN        STD         MEAN        STD
#          ======================= ======================= =======================
ELE            -8656.78      70.18     -5602.09      63.10     -2102.25      52.57
VDW             -984.99      24.34      -661.18      20.33      -256.02      12.93
INT             5085.33      50.22      3449.57      38.65      1635.76      29.42
GAS            -4556.44      75.96     -2813.70      65.21      -722.52      53.50
PBSUR             65.09       1.05        45.25       0.64        27.24       0.46
PBCAL          -3223.64      58.68     -2490.86      48.73     -1671.27      47.46
PBSOL          -3158.55      58.26     -2445.62      48.45     -1644.03      47.31
PBELE         -11880.42      34.25     -8092.96      29.34     -3773.52      17.30
PBTOT          -7714.99      48.25     -5259.32      36.97     -2366.55      26.61
GBSUR             65.09       1.05        45.25       0.64        27.24       0.46
GB             -3407.82      58.49     -2631.83      50.08     -1731.06      47.68
GBSOL          -3342.74      58.15     -2586.58      49.83     -1703.82      47.55
GBELE         -12064.60      26.94     -8233.92      23.57     -3833.31      13.40
GBTOT          -7899.17      47.07     -5400.28      35.65     -2426.34      26.80

#                    DELTA        
#          -----------------------
#                  MEAN        STD
#          =======================
ELE             -952.43      44.10
VDW              -67.79       5.18
INT               -0.00       0.00
GAS            -1020.22      44.58
PBSUR             -7.40       0.41
PBCAL            938.50      42.51
PBSOL            931.09      42.31
PBELE            -13.94       9.43
PBTOT            -89.13       7.94
GBSUR             -7.40       0.41
GB               955.07      41.30
GBSOL            947.66      41.10
GBELE              2.63       7.41
GBTOT            -72.56       6.40

The meaning of the different terms in this file is as follows:

ELE = electrostatic energy as calculated by the MM force field.
VDW = van der Waals contribution from MM.
INT = internal energy arising from bond, angle, and dihedral terms in the MM force field. (this term always amounts to zero in the single trajectory approach).
GAS = total gas phase energy (sum of ELE, VDW, and INT).
PBSUR/GBSUR = non-polar contribution to the solvation free energy calculated by an empirical model.
PBCAL/GB = the electrostatic contribution to the solvation free energy calculated by PB or GB, respectively.
PBSOL/GBSOL = sum of non-polar and polar contributions to solvation.
PBELE/GBELE = sum of the electrostatic solvation free energy and MM electrostatic energy.
PBTOT/GBTOT = final estimated binding free energy calculated from the terms above. (kcal/mol)

There are several things to note here. Normally the ELE, PBCAL/GBCAL, and VDW parts make up the major contributions to the binding free energy and the first two of these terms tend to approximately cancel each other out which can be checked by the value of PBELE/GBELE which should be much smaller than the contributions to it.

One would typically expect to find an extremely favourable electrostatic energy and a dis-favourable solvation free energy. This symbolises the energy that ones has to use to de-solvate the binding particles and to align their binding interfaces.

From the negative total binding free energy -89.13 kcal/mol we clearly see that this is a favourable protein-protein complex in pure water but keep in mind that the result does not equal the real binding free energy since we did not estimate the (dis-favourable) entropy contribution to binding. Note that the GB approach gives a slightly lower binding energy but still suggests that this is a favourable bound state.

To extend this tutorial you could investigate what happens if you change the salt content of the solvent, modify specific residues or choose different starting geometries. You could also look into running the nmode calculations for finding the entropy but be aware that for a complex of this size this will be extremely memory and computationally expensive.


RETURN TO INTRODUCTION


(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2006