2014 AMBER tutorial with HIV Protease

For additional Rizzo Lab tutorials see AMBER Tutorials.

In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).

I. Introduction

Mike and Ivan

IV. Simulation Analysis

MM-GBSA Energy Calculation

MM/GBSA is the acronym for Molecular Mechanics/Generalized Born Surface Area. This part of AMBER combines molecular mechanics (MM) with both the electrostatic (PB) and nonpolar (SA) contribution to solvation . Topology files are needed for the receptor, ligand, and receptor-ligand complex. The trajectory files generate the coordinates. Therefore, molecular dynamics is used to generate a set of snapshots taken at fixed intervals from the trajectories. These snapshots are processed to remove solvent and generate the free energy of binding.

In the AMBER_Tutorial directory, create a new directory:

 mkdir 005.MMGBSA

In this new directory, create the file gb.rescore.in containing:

 Single point GB energy calc
&cntrl
 ntf    = 1,        ntb    = 0,        ntc    = 2,
 idecomp= 0,
 igb    = 5,        saltcon= 0.00,
 gbsa   = 2,        surften= 1.0,
 offset = 0.09,     extdiel= 78.5,
 cut    = 99999.0,  nsnb   = 99999,
 imin   = 5,        maxcyc = 1,        ncyc   = 0,
 /

Then create a csh script, run.sander.rescore.csh, that contains the following lines of command:

 #! /bin/csh
 #PBS -l nodes=1:ppn=2
 #PBS -l walltime=24:00:00
 #PBS -o zzz.qsub.out
 #PBS -e zzz.qsub.err
 #PBS -V
 set workdir = "/nfs/user03/yechen1/AMBER-Tutorial/005.MMGBSA"
 set sander  = "/nfs/user03/wjallen/local/amber12/bin/sander"
 cd ${workdir}
 $sander -O -i gb.rescore.in -o gb.rescore.out.com -p ../002.ANTE.TLEAP  /1LOQ.com.gas.leap.prm7 \
 -c ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.com.trj.stripfit \
 -r restrt.com -ref ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7 -x mdcrd.com -inf mdinfo.com
 $sander -O -i gb.rescore.in -o gb.rescore.out.lig -p ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.prm7 \
 -c ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.lig.trj.stripfit \
 -r restrt.lig -ref ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.rst7 -x mdcrd.lig -inf mdinfo.lig
 $sander -O -i gb.rescore.in -o gb.rescore.out.rec -p ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.prm7 \
 -c ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.rec.trj.stripfit \
 -r restrt.rec -ref ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.rst7 -x mdcrd.rec -inf mdinfo.rec
 exit

~ ~ ~ ~

Then this script should be sent to the queue, i.e., qsub the script using the commands:

 qsub run.sander.rescore.csh

You can monitor your progress by typing

 qstat -u username

When the job is complete, you will obtain the following output files: gb.rescore.out.com, gb.rescore.out.lig, and gb.rescore.out.rec In these files, the single point energy calculation results will be written for each individual frame. It will be found in the results section and the output file will have an infrastrucutre that is similar to the following:

     FINAL RESULTS
   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
     1       3.6269E+03     1.8737E+01     1.0472E+02     CB        585
  BOND    =      580.2786  ANGLE   =     1563.7704  DIHED      =     2161.5659
  VDWAALS =    -1684.2762  EEL     =   -13809.8494  EGB        =    -2953.6681
  1-4 VDW =      756.7767  1-4 EEL =     7260.2823  RESTRAINT  =        0.0000
  ESURF   =     9752.0291
 minimization completed, ENE= 0.36269092E+04 RMS= 0.187371E+02
 TRAJENE: Trajectory file ended
 TRAJENE: Trajene complete.

In the command line, type:

 grep VDWAALS gb.rescore.out.com > vdw.com.txt.
 grep ESURF gb.rescore.out.com > surf.com.txt.   

You can take these text files, import them into Excel, and do the rest of your modifications there.

Equations for analysis

Remember that to obtain the Gvdw term, you need to take the SASA (which is ESURF) and input into equation 1:

ΔGnonpolar = SASA*0.00542 + 0.92

Also, the mmgbsa of a given system can be determined by equation 2:

ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar

From the output file:

VDWAALS = ΔGvdw

EELS = ΔGcoul

EGB = ΔGpolar

You can then easily calculate the ΔΔGbind by using equation 3:

ΔΔGbind = ΔGmmgbsa,complex – (ΔGmmgbsa,lig + ΔGmmgbsa,rec) You will want to careful when doing your analysis that the results from frame 1 for the receptor and ligand are subtracted from the results from frame 1 for your complex. By doing this in excel, you should have 2000 frames for each, and the values should cleanly line up. Finally, you will want to plot your ΔΔGbind and examine if you see corresponding changes in the ligand position and the ΔΔGbind. Also, you should determine the mean and standard deviation for your ΔΔGbind.

Plotting Energy

When your rescoring calculation finishes, you have the following three output files: "gb.rescore.out.com", "gb.rescore.out.lig", and "gb.rescore.out.rec".

Use the following script, entitled get.mmgbsa.bash, to extract data and calculate MMGBSA energy for each snap shot.

 #! /bin/bash
 # by Haoquan
 echo com lig rec > namelist
 LIST=`cat namelist`
 for i in $LIST ; do
 grep VDWAALS gb.rescore.out.$i | awk '{print $3}' > $i.vdw
 grep EGB     gb.rescore.out.$i | awk '{print $9}' > $i.polar
 grep EELS    gb.rescore.out.$i | awk '{print $6}' > $i.coul
 grep ESURF   gb.rescore.out.$i | awk '{print $3 * 0.00542 + 0.92}' > $i.surf
 paste -d " " $i.vdw $i.polar $i.surf $i.coul | awk '{print $1 + $2 + $3 + $4}' > data.$i
 rm $i.*
 done
 paste -d " " data.com data.lig data.rec | awk '{print $1 - $2 - $3}' > data.all 
 for ((j=1; j<=`wc -l data.all | awk '{print $1}'`; j+=1)) do
 echo $j , >> time
 done
 paste -d " " time data.all > MMGBSA_vs_time.dat  
 rm namelist time data.*

To run this script do:

 bash get.mmgbsa.sh

This will create a text file called MMGBSA_vs_time.dat with x and y values separated by a space and comma. These values can be imported to Excel or Origin or to XMGRACE if you are using Linux:

 xmgrace MMGBSA_vs_time.dat

2013 AMBER MMGBSA plot.jpg