Difference between revisions of "2016 AMBER tutorial with Beta Trypsin"

From Rizzo_Lab
Jump to: navigation, search
(→‎III. Simulation using pmemd)
(→‎MM-GBSA Energy Calculation)
Line 78: Line 78:
 
===MM-GBSA Energy Calculation===
 
===MM-GBSA Energy Calculation===
  
−
Monaf
+
begin by creating a new directory
 +
 
 +
  mkdir 005.MMGBSA/
 +
 
 +
as well as an input file
 +
 
 +
  vim.gbrescore.in
 +
 
 +
the following parameters will provide a good estimate of the energy of binding based on solvation enthalpy and binding enthalpy
 +
 
 +
  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,
 +
/
 +
 
 +
We can now submit a job to the cluster we are working on. This script, inside a file named run.sander.rescore.csh , is designed for LIRED, however it can be modified for whatever supercomputer you are working on
 +
 
 +
  #! /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/guest46/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
 +
 
 +
Then submit the job as a .csh
 +
 
 +
  qsub run.sander.rescore.csh
 +
 
 +
Your ligand, receptor, and receptor-ligand complex will all have .out files. The format should look
 +
 
 +
  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
 +
 
 +
We can now use a system of equations to elucidate various energies related to our system.
 +
VDWAALS = ΔGvdw
 +
E(ELS) = ΔGcoulombic
 +
E(GB) = ΔGpolar
 +
SASA = ESURF
 +
 
 +
 
 +
ΔGnonpolar can be estimated by:
 +
 
 +
ΔGnonpolar = (Solvent Accessible Surface Area)*0.00542 + 0.92
 +
 
 +
This equation is based on empirical observations relating a molecules surface area and solvation energy
 +
 
 +
ΔGnonpolar and ΔGmmgbsa are related by anther equation:
 +
 
 +
ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar
 +
 
 +
all the information needed is in our output files.
 +
 
 +
ΔΔGbind can now be solved via the equation:
 +
 
 +
ΔΔGbind = ΔG(mmgbsa complex) – (ΔG(mmgbsa ligand) + ΔG(mmgbsa receptor) (3)
 +
 
 +
ΔΔGbind is a function of ligand position, so it is useful to plot it. For publications, it is standard protocol to include the mean and standard deviation for your ΔΔGbind.
 +
 
 +
Lets make a script named get.mmgbsa.csh to extract this sophisticated energetic estimation from the three initial output files of the calculations obtained solve for ΔΔGbind:
 +
 
 +
    bash get.mmgbsa.csh
 +
 
 +
The output file will be .dat and we can easily use the program xmgrace to plot the information visually
 +
   
 +
  xmgrace MMGBSA_vs_time.dat
  
 
==V. Frequently Encountered Problems==
 
==V. Frequently Encountered Problems==

Revision as of 07:47, 10 April 2016

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

Yaping

II. Structural Preparation

Antechamber, Parmchk, tLeap

Omar, Katie

III. Simulation using pmemd

PMEMD

Agatha, Beilei

MD in put file1 01md.in: 

mkdir 01md.in

01mi.in: equilibration
&cntrl
imin = 1, maxcyc = 1000, ntmin = 2,
ntx = 1, ntc = 1, ntf = 1,
ntb = 1, ntp = 0,
ntwx = 1000, ntwe = 0, ntpr = 1000,
cut = 8.0,
ntr = 1,
restraintmask = ':1-224 & !@H=',
restraint_wt = 5.0,


Copy 01md.in and rename it to 02md.in:

IV. Simulation Analysis

Ptraj

Lauren, Haoyue

ptraj.strip.wat.in: 

trajin ../003.PMEMD/10md.trj 1 1000 1
trajin ../003.PMEMD/11md.trj 1 1000 1
trajout 1BJU.trj.strip nobox
strip :WAT

ptraj.rec.rmsd.in: 

trajin 1BJU.trj.strip
trajout 1BJU.com.trj.stripfit 
reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7
rms reference out 1BJU.rmsd.CA.dat :1-223@CA

ptraj.lig.rmsd.in: 

trajin 1BJU.com.trj.stripfit
reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7
rms reference out 1BJU.lig.rmsd.dat :224@C*,N*,O*,S* nofit

ptraj.strip.rec.in: 

trajin 1BJU.com.trj.stripfit 1 2000 1
trajout 1BJU.lig.trj.stripfit
strip :1-223

ptraj.strip.lig.in: 

trajin 1BJU.com.trj.stripfit 1 2000 1
trajout 1BJU.rec.trj.stripfit
strip :224

ptraj.dist.in: 

trajin 1BJU.com.trj.stripfit
distance 171OD2_224N1 :171@OD2 :224@N1 out 171OD2_224N1.dat
distance 171OD1_224N2 :171@OD1 :224@N2 out 171OD1_224N2.dat

MM-GBSA Energy Calculation

begin by creating a new directory 

  mkdir 005.MMGBSA/

as well as an input file 

  vim.gbrescore.in

the following parameters will provide a good estimate of the energy of binding based on solvation enthalpy and binding enthalpy 

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

We can now submit a job to the cluster we are working on. This script, inside a file named run.sander.rescore.csh , is designed for LIRED, however it can be modified for whatever supercomputer you are working on 

  #! /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/guest46/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

Then submit the job as a .csh 

  qsub run.sander.rescore.csh

Your ligand, receptor, and receptor-ligand complex will all have .out files. The format should look 

  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

We can now use a system of equations to elucidate various energies related to our system. VDWAALS = ΔGvdw E(ELS) = ΔGcoulombic E(GB) = ΔGpolar SASA = ESURF


ΔGnonpolar can be estimated by:

ΔGnonpolar = (Solvent Accessible Surface Area)*0.00542 + 0.92

This equation is based on empirical observations relating a molecules surface area and solvation energy

ΔGnonpolar and ΔGmmgbsa are related by anther equation:

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

all the information needed is in our output files.

ΔΔGbind can now be solved via the equation:

ΔΔGbind = ΔG(mmgbsa complex) – (ΔG(mmgbsa ligand) + ΔG(mmgbsa receptor) (3)

ΔΔGbind is a function of ligand position, so it is useful to plot it. For publications, it is standard protocol to include the mean and standard deviation for your ΔΔGbind.

Lets make a script named get.mmgbsa.csh to extract this sophisticated energetic estimation from the three initial output files of the calculations obtained solve for ΔΔGbind: 

   bash get.mmgbsa.csh

The output file will be .dat and we can easily use the program xmgrace to plot the information visually 

  xmgrace MMGBSA_vs_time.dat

V. Frequently Encountered Problems

Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2016_AMBER_tutorial_with_Beta_Trypsin&oldid=9768"