(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 McGee, Miller, and Swails 2009
Python Script MMPBSA.py
Dwight McGee, Bill Miller III, and Jason Swails
The important files for calculating the binding free energy using MMPBSA.py are the topology files and the mdcrd file (ras-raf_top_mdcrd.tgz)
Calculate the binding free energy of Ras-Raf.
We will now calculate the interaction energy and solvation free energy for the 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 part of the tutorial and so strictly speaking our result will not be a true free energy but could be used to compare against similar systems. See Section 3.5 for an example of using Normal Mode Analysis (Nmode) to calculate the entropy contribution for a system or uncomment out the last line in the &general namelist of the input file below to perform a Quasi-Harmonic entropy calculation using the ptraj module in AMBER.
We will carry out the binding energy calculation using both the MM-GBSA method and the MM-PBSA method for comparison. This is accomplished with the following input file for MMPBSA.py:
| mmpbsa.in |
Input file for running PB and GB &general endframe=50, verbose=1, # entropy=1, / &gb igb=2, saltcon=0.100 / &pb istrng=0.100, / |
The input files for MMPBSA.py are designed to be similar to the setup of an mdin file used in the sander module of AMBER. The start of each namelist is designated by an ampersand (&) followed by the name of the namelist. Furthermore, a backslash (/) or '&end' can be used to end the namelist. For a complete list of all variables please see the Amber Manual or type $AMBERHOME/bin/MMPBSA.py --input-file-help. This input file is divided into three namelists: general, pb, and gb. The general namelist is designed to specify variables that are not specific to a particular part of the calculation, but to all parts. In this setup we have defined RAS to be the receptor and RAF to be the ligand. The 'endframe' variable sets what frame of the mdcrd to stop on. The '&gb' and '&pb' namelist markers let the script know to perform MM-GBSA and MM-PBSA calculations with the given values defined within those namelists. The 'verbose' variable allows the user to specify how much output is written to the output file.
$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp ras-raf_solvated.prmtop -cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y *.mdcrd
This will run the script interactively and print the progress of the calculation to STDOUT and any errors or warnings to STDERR. Finally, timings will be printed once the calculation has completed showing the time taken during each step of the calculation.
Command-line arguments can be given with shell-recognized wildcards (i.e. * and ? for bash). For example, the '-y *.mdcrd' on the command line tells the script to read in all files in the working directory that end in '.mdcrd' and use them as the trajectories to be analyzed.
Here are all the output files created by this script: pb_gb_output1.tgz.
The script creates three unsolvated mdcrd files (complex, receptor, and ligand) using ptraj that are the coordinates analyzed during the GB and PB calculations. The *.mdout files contain the energies for all frames specified. A PDB file of the average structure is created align (via RMS) all snapshots to prepare for a quasi-harmonic entropy calculation with ptraj if one is requested. All files created by MMPBSA.py should begin with the prefix '_MMPBSA_' except for the final output file, FINAL_RESULTS_MMPBSA.dat.
| FINAL_RESULTS_MMPBSA.dat |
| Run on Thu Feb 11 12:18:37 EST 2010 |Input file: |-------------------------------------------------------------- |Input file for running PB and GB |&general | endframe=50, verbose=1, |# entropy=1, |/ |&gb | igb=2, saltcon=0.100 |/ |&pb | istrng=0.100, |/ |-------------------------------------------------------------- |Solvated complex topology file: ras-raf_solvated.prmtop |Complex topology file: ras-raf.prmtop |Receptor topology file: ras.prmtop |Ligand topology file: raf.prmtop |Initial mdcrd(s): prod.mdcrd | |Best guess for receptor mask: ":1-166" |Best guess for ligand mask: ":167-242" |Calculations performed using 50 frames. |Poisson Boltzmann calculations performed using internal PBSA solver in sander. | |All units are reported in kcal/mole. ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- GENERALIZED BORN: Complex: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1863.7944 17.1704 2.4283 EEL -17200.7297 75.9366 10.7391 EGB -3249.6511 65.2075 9.2217 ESURF 91.3565 1.3938 0.1971 G gas -19064.5240 77.8536 11.0102 G solv -3158.2946 65.2224 9.2238 TOTAL -22222.8186 51.0216 7.2155 Receptor: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1268.1888 14.2342 2.0130 EEL -11557.0773 71.7127 10.1417 EGB -2532.0669 57.7003 8.1600 ESURF 64.2843 1.1143 0.1576 G gas -12825.2661 73.1118 10.3396 G solv -2467.7826 57.7110 8.1616 TOTAL -15293.0487 35.3527 4.9996 Ligand: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -529.3090 9.4198 1.3322 EEL -4684.4720 36.1449 5.1117 EGB -1688.9631 26.5353 3.7527 ESURF 37.0493 0.6185 0.0875 G gas -5213.7811 37.3522 5.2824 G solv -1651.9138 26.5425 3.7537 TOTAL -6865.6949 25.8878 3.6611 Differences (Complex - Receptor - Ligand): Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -66.2966 4.2751 0.6046 EEL -959.1803 34.9190 4.9383 EGB 971.3789 33.0497 4.6739 ESURF -9.9770 0.3759 0.0532 DELTA G gas -1025.4769 35.1797 4.9752 DELTA G solv 961.4018 33.0518 4.6742 DELTA G binding = -64.0750 +/- 6.3729 0.9013 ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- POISSON BOLTZMANN: Complex: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1863.7944 17.1704 2.4283 EEL -17200.7297 75.9366 10.7391 EPB -3207.7160 66.4023 9.3907 ECAVITY 67.8762 0.7818 0.1106 G gas -19064.5240 6061.1875 857.1813 G solv -3139.8399 66.4069 9.3914 TOTAL -7686.8660 52.5400 7.4303 Receptor: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1268.1888 14.2342 2.0130 EEL -11557.0773 71.7127 10.1417 EPB -2483.7242 56.4551 7.9840 ECAVITY 47.1495 0.4737 0.0670 G gas -12825.2661 5345.3320 755.9441 G solv -2436.5747 56.4571 7.9842 TOTAL -5250.2060 38.5188 5.4474 Ligand: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -529.3090 9.4198 1.3322 EEL -4684.4720 36.1449 5.1117 EPB -1670.4169 27.6694 3.9131 ECAVITY 28.0328 0.4133 0.0584 G gas -5213.7811 1395.1865 197.3092 G solv -1642.3841 27.6725 3.9135 TOTAL -2350.3020 25.1197 3.5525 Differences (Complex - Receptor - Ligand): Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -66.2966 4.2751 0.6046 EEL -959.1803 34.9190 4.9383 EPB 946.4251 34.5128 4.8808 ECAVITY -7.3062 0.3004 0.0425 DELTA G gas -1025.4769 1237.6138 175.0250 DELTA G solv 939.1189 34.5141 4.8810 DELTA G binding = -86.3579 +/- 8.3264 1.1775 ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- WARNINGS: igb=2 should be used with mbondi2 pbradii set. Yours are modified Bondi radii (mbondi) |
The beginning of the statistics file includes the date/time, any warnings based on the values and files given, the mmpbsa.in text, the files used by the script, the number of frames analyzed, and which PB solver (if any) was used. The rest of the statistics file includes all the average energies, standard deviations, and standard error of the mean for GB followed by PB. After each section, the ΔG of binding is given along with the error values. The meaning of the different terms in this file is as follows:
VDWAALS = van der Waals contribution from MM.
EEL = electrostatic energy as calculated by the MM force field.
EPB/EGB = the electrostatic contribution to the solvation free energy calculated by PB or GB respectively.
ECAVITY = nonpolar contribution to the solvation free energy calculated by an empirical model.
DELTA G binding = final estimated binding free energy calculated from the terms above. (kCal/mol)
Note that the total gas phase energy has not been reported because the values of the bonded potential terms for the receptor and ligand should exactly cancel those for the complex using the single trajectory approach. An error message will result if the energies do not cancel within an allowed tolerance.
One would typically expect to find an extremely favorable electrostatic energy and a unfavorable solvation free energy. This symbolises the energy that one has to use to de-solvate the binding particles and to align their binding interfaces.
From the negative total binding free energy -86.36 kcal/mol we clearly see that this is a favorable 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 (unfavorable) entropy contribution to binding. Note that the GB approach gives a slightly lower binding energy but still suggests that this is a favorable bound state.
CLICK HERE TO GO TO SECTION 3.2
CLICK HERE TO GO TO MMPBSA.py HOME
CLICK HERE TO GO BACK TO THE 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 McGee, Miller, and Swails 2009