2016 AMBER tutorial with Beta Trypsin

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 input 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:

cp 01md.in 02md.in

02mi.in: equilibration
&cntrl
  imin = 0, ntx = 1, irest = 0, nstlim = 50000,
  temp0 = 298.15, tempi = 298.15, ig = 71287,
  ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
  ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
  ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
  cut = 8.0, iwrap = 1,
  ntr = 1, nscm = 100,
  restraintmask = ':1-224 & !@H=', restraint_wt = 5.0,

Using the same way to make 03md.in, 04md.in, 05md.in, 06md.in, 07md.in, 08md.in, 09md.in, 10md.in, 11md.in:

03mi.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 = 2.0,

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