Amber masthead
Filler image AmberTools20 Amber20 Manuals Tutorials Force Fields Contacts History
Filler image

Useful links:

Amber Home
Download Amber
Installation
News
Amber Citations
GPU Support
Intel Support
Updates
Mailing Lists
For Educators
File Formats
7 Free Energies
 

(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

1) Build the starting structure and run a simulation to obtain an equilibrated system.

The system we will model in this simulation is the protein-ligand complex between the Estrogen Receptor protein and Raloxifene ligand. Here is a pre-prepared pdb file of the complex.

Estrogen_Receptor-Raloxifene.pdb

This structure contains the Estrogen Receptor protein along with a ligand called Raloxifene that has been previously docked to the protein as illustrated in the figure below:

This system was constructed in a similar manner to the Ras-Raf system in Section 1. For directions on how to build the starting structure and run a simulation to obtain an equilibrated system, please refer to Section 1 and Section 2 Note that you must also use antechamber to get the correct parameters for raloxifene. See the Sustiva Tutorial (Basic 4) for detailed instructions.

The important files from the MD simuation for calculating the binding free energy using MMPBSA.py are the topology files and the mdcrd file (Est_Rec_top_mdcrd.tgz)

2) Calculate the binding free energy of the Estrogen Receptor and Raloxifene

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 section 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, keep_files=2,
/
&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. 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 rather to all parts. In this setup the Estrogen Receptor is the receptor and Raloxifene is the ligand. The 'endframe' variable sets which frame of the mdcrd to stop on. The presence of the '&gb' and '&pb' namelists let the script know to perform MM-GBSA and MM-PBSA calculations with the given input values.

$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp 1err.solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.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.

With 'keep_files=2', here are all the output files: pb_gb_output2.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:44:26 EST 2010

|Input file:
|--------------------------------------------------------------
|Input file for running PB and GB
|&general
|   endframe=50, keep_files=2,
|/
|&gb
|  igb=2, saltcon=0.100,
|/
|&pb
|  istrng=0.100, 
|/
|--------------------------------------------------------------
|Solvated complex topology file:  1err.solvated.prmtop
|Complex topology file:           complex.prmtop
|Receptor topology file:          receptor.prmtop
|Ligand topology file:            ligand.prmtop
|Initial mdcrd(s):                1err_prod.mdcrd
|
|Best guess for receptor mask:   ":1-240"
|Best guess for  ligand  mask:   ":241"
|Ligand residue name is "RAL"
|
|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                  -2013.3801               20.3021              2.8712
EEL                     -16938.6450               85.7631             12.1287
EGB                      -3507.0086               67.7839              9.5861
ESURF                       97.5448                1.3301              0.1881

G gas                   -18952.0251               88.1333             12.4639
G solv                   -3409.4639               67.7969              9.5879

TOTAL                   -22361.4889               47.1982              6.6748


Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1955.2272               19.2311              2.7197
EEL                     -16895.0354               85.5797             12.1028
EGB                      -3528.7276               68.3585              9.6673
ESURF                      101.2613                1.3071              0.1849


G gas                   -18850.2626               87.7138             12.4046
G solv                   -3427.4663               68.3710              9.6691

TOTAL                   -22277.7288               48.1057              6.8032


Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                     -1.8595                2.0516              0.2901
EEL                         -5.5796                2.0333              0.2876
EGB                        -28.4863                0.6040              0.0854
ESURF                        4.4326                0.0462              0.0065


G gas                       -7.4391                2.8885              0.4085
G solv                     -24.0538                0.6058              0.0857

TOTAL                      -31.4929                5.0748              0.7177


Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -56.2934                2.9265              0.4139
EEL                        -38.0300                3.2114              0.4542
EGB                         50.2053                2.5869              0.3658
ESURF                       -8.1491                0.2589              0.0366


DELTA G gas                -94.3234                4.3449              0.6145
DELTA G solv                42.0562                2.5999              0.3677


 DELTA G binding =        -52.2672     +/-      2.4568                 0.3475
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

POISSON BOLTZMANN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -2013.3801               20.3021              2.8712
EEL                     -16938.6450               85.7631             12.1287
EPB                      -3329.1708               67.0354              9.4802
ECAVITY                     68.2656                0.5195              0.0735

G gas                   -18952.0251             7767.4837           1098.4881
G solv                   -3260.9052               67.0374              9.4805

TOTAL                    -5265.0831               49.0426              6.9357


Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1955.2272               19.2311              2.7197
EEL                     -16895.0354               85.5797             12.1028
EPB                      -3355.4746               67.3299              9.5219
ECAVITY                     70.1184                0.5285              0.0747

G gas                   -18850.2626             7693.7163           1088.0558
G solv                   -3285.3562               67.3320              9.5222

TOTAL                    -5279.4509               50.4067              7.1286


Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                     -1.8595                2.0516              0.2901
EEL                         -5.5796                2.0333              0.2876
EPB                        -31.3364                0.6953              0.0983
ECAVITY                      3.1896                0.0288              0.0041

G gas                       -7.4391                8.3434              1.1799
G solv                     -28.1468                0.6959              0.0984

TOTAL                       56.0934                5.0476              0.7138


Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -56.2934                2.9265              0.4139
EEL                        -38.0300                3.2114              0.4542
EPB                         57.6402                3.0642              0.4333
ECAVITY                     -5.0423                0.1683              0.0238

DELTA G gas                -94.3234               18.8778              2.6697
DELTA G solv                52.5978                3.0688              0.4340



 DELTA G binding =        -41.7256     +/-      2.9618                 0.4189
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

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 ones has to use to de-solvate the binding particles and to align their binding interfaces.

From the negative total binding free energy -41.73 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 PB approach gives a slightly lower binding energy but still suggests that this is a favorable bound state.


CLICK HERE TO GO TO SECTION 3.3

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