2016 AMBER tutorial with Thrombin

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).

IV. Simulation Analysis

MM-GBSA Energy Calculation

Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) is a great method to calculate or estimate relative binding affinity of a ligand(s) to a receptor. The binding energy calculated from this method are also known as free energies of binding, where the more negative values indicate stronger binding. For this section, the topology files for the ligand, receptor and complex are needed.

Create a new directory:

 mkdir 005.MMGBSA

Create an input file name

 vim gb.rescore.in

Enter the following into the input file:

  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,
 /

Create a tcsh/bash/csh script (run.sander.rescore.csh) with the following information:

 #! /bin/tcsh
 #PBS -l nodes=1:ppn=1
 #PBS -l walltime=48:00:00
 #PBS -o zzz.qsub.out
 #PBS -e zzz.qsub.err
 #PBS -V
 #PBS -N mmgbsa
 
 set workdir = /nfs/user03/kbelfon/amber_tutorial/005.mmgbsa
 
 cd $workdir
 
 sander -O -i gb.rescore.in \
 -o gb.rescore.out.com \
 -p ../002.tleap/4TKG.com.gas.leap.prm7 \
 -c ../002.tleap/4TKG.com.gas.leap.rst7 \
 -y ../004.ptraj/4TKG.com.trj.stripfit \
 -r restrt.com \
 -ref ../002.tleap/4TKG.com.gas.leap.rst7 \
 -x mdcrd.com \
 -inf mdinfo.com

 sander -O -i gb.rescore.in \
 -o gb.rescore.out.lig \
 -p ../002.tleap/4TKG.lig.gas.leap.prm7 \
 -c ../002.tleap/4TKG.lig.gas.leap.rst7 \
 -y ../004.ptraj/4TKG.lig.trj.stripfit \
 -r restrt.lig \
 -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \
 -x mdcrd.lig \
 -inf mdinfo.lig 
 
 sander -O -i gb.rescore.in \
 -o gb.rescore.out.test.rec \
 -p ../002.tleap/4TKG.rec.gas.leap.prm7 \
 -c ../002.tleap/4TKG.rec.gas.leap.rst7 \
 -y ../004.ptraj/4TKG.rec.trj.stripfit \
 -r restrt.rec \
 -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \
 -x mdcrd.rec \
 -inf mdinfo.rec 
 
 exit

Execute this script on the seawulf cluster or machine(s) of your preference

 qsub run.sander.rescore.csh

Three output files will be generated once the job is completed: gb.rescore.out.com, gb.rescore.out.lig, and gb.rescore.out.rec These files represent the single point energy calculation results for the complex (.com), the ligand (.lig) and the receptor (.rec). The energy will be output by the program Sander for each frame specified in the input file. The final results for one frame in one of the three files should look as the following:

                                   FINAL RESULTS
  
  
   
   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1       5.9132E+03     2.0005E+01     1.2640E+02     C         159
  
 BOND    =      661.8980  ANGLE   =     1751.7992  DIHED      =     2581.7692
 VDWAALS =    -1696.6585  EEL     =   -13958.9335  EGB        =    -3125.9524
 1-4 VDW =      747.0185  1-4 EEL =     7750.8118  RESTRAINT  =        0.0000
 ESURF   =    11201.4791
minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02

Extracting Data from MM-GBSA calculation and calculating Free energy of Binding From the output files above:

VDWAALS = ΔGvdw

EELS = ΔGcoul

EGB = ΔGpolar

SASA = ESURF

With this information ΔGnonpolar can be solved using equation(1):

ΔGnonpolar = SASA*0.00542 + 0.92 (1)

Once ΔGnonpolar is solved then ΔGmmgbsa can be determined by equation(2):

ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar (2)

Solve equation 2 and 3 using the extracted information from all three output files. So therefore you should have ΔGmmgbsa for the complex, ligand and receptor

Finally ΔΔGbind can be calculated using equation (3):

ΔΔGbind = ΔGmmgbsa,complex – (ΔGmmgbsa,lig + ΔGmmgbsa,rec) (3)

Plot your ΔΔGbind and examine the plot for changes in the ligand position and the ΔΔGbind. Also, you should calculate the mean and standard deviation for your ΔΔGbind.

The following script (get.mmgbsa.sh) can be used to extract the energy from the three output files obtained above and calculate ΔΔGbind:

 #! /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.*

Execute this script:

 bash get.mmgbsa.sh

A text file called MMGBSA_vs_time.dat with x and y values separated by a space and comma should be created. Use XMGRACE to plot this dat file using the following command in Linux:

 xmgrace MMGBSA_vs_time.dat