(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

MM-PBSA

By Ross Walker & Thomas Steinbrecher

3) Calculate the binding free energy and analyse the results.

Starting with the snapshots we extracted at the end of the last section 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 calculate 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.

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

@GENERAL
PREFIX                snapshot
PATH                  ./
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
PRBRAD                1.4
ISTRNG                0.0
RADIOPT               0
NPOPT                 1
CAVITY_SURFTEN        0.0072
CAVITY_OFFSET         0.00
SURFTEN               0.0072
SURFOFF               0.00
@MM
DIELC                 1.0
@GB
IGB                   2
GBSA                  1
SALTCON               0.00
EXTDIEL               80.0
INTDIEL               1.0
SURFTEN               0.0072
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 explanations of the different terms. Values for the different parameters are determined based on empirical data and are subject to on-going research. 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/exe/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. At present the mm_pbsa.pl script is not parallel but one could easily modify things to split the trajectory and run several mm_pbsa analysis scripts in parallel on a cluster for example.

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:

#                  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             93.72       1.51        65.15       0.93        39.23       0.67
PBCAL          -3244.39      59.67     -2516.48      50.04     -1672.39      47.87
PBSOL          -3150.66      59.09     -2451.33      49.64     -1633.16      47.67
PBELE         -11901.17      33.10     -8118.58      28.08     -3774.64      16.64
PBTOT          -7707.10      47.89     -5265.03      36.65     -2355.68      26.66
GBSUR             93.72       1.51        65.15       0.93        39.23       0.67
GB             -3250.03      59.24     -2540.49      50.86     -1672.13      48.22
GBSOL          -3156.30      58.77     -2475.34      50.52     -1632.90      48.04
GBELE         -11906.80      27.15     -8142.59      23.52     -3774.38      13.43
GBTOT          -7712.74      47.25     -5289.04      35.96     -2355.41      27.07

#                    DELTA        
#          -----------------------
#                  MEAN        STD
#          =======================
ELE             -952.43      44.10
VDW              -67.79       5.18
INT               -0.00       0.00
GAS            -1020.22      44.58
PBSUR            -10.66       0.58
PBCAL            944.49      43.14
PBSOL            933.83      42.84
PBELE             -7.95       9.81
PBTOT            -86.39       8.24
GBSUR            -10.66       0.58
GB               962.59      41.40
GBSOL            951.93      41.11
GBELE             10.16       7.77
GBTOT            -68.29       6.58

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 = nonpolar 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 nonpolar 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-avourable 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 -86.39 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