(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
MM-PBSA - 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 |
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.
(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