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

Revision as of 12:49, 12 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

Antechamber, Parmchk, tLeap

Before beginning the Molecular Dynamics protocol using AMBER, you must first set up your files. In your 001.chimera folder, you will add the following files:

 1BJU.lig.mol2
 1BJU.rec.noH.pdb
 1BJU.rec.noH.ZIN.pdb

To prepare the first two files, please see the 2025 DOCK tutorial at the following link: [1]

Note: delete any headers before the atoms/helix information.

In order to prepare 4TKG.rec.noH.ZIN.pdb, first open 4TKG.rec.noH.pdb and use the command "SHIFT + g" to reach the bottom of the pdb. The last lines of your files should look like this:

 ATOM   1648  OE1 GLU B 205      23.839 -23.190  57.747  1.00  0.00           O
 TER    1649      GLU B 205
 HETATM 1650 ZN   ZN  B 206      28.130  -3.467  55.482  1.00  0.00          Zn
 END

To make 4TKG.rec.noH.ZIN.pdb, you will need to change the "ZN" atom ID to "ZIN" so that AMBER can read the atom type.

 ATOM   1648  OE1 GLU B 205      23.839 -23.190  57.747  1.00  0.00           O
 TER    1649      GLU B 205
 ATOM   1650  ZIN ZIN B 206      28.130  -3.467  55.482  1.00  0.00          Zn
 END

To run a tleap, first generate a 002.tleap file, with the command:

 mkdir 002.tleap

 cp ~/../../tleap.lig.in

• tleap.lig.in

source leaprc.ff14SB                                               #Load a force field
source leaprc.gaff                  
LIG = loadmol2 4TKG.lig.mol2                                       #Conformation file
saveamberparm LIG 4TKG.lig.gas.leap.prm7 4TKG.lig.gas.leap.rst7    #Save the ligand gas phase AMBER topology and coordinate file
solvatebox LIG TIP3PBOX 10.0                                       #Solvate the lig using TIP3P, solvent box radii 10 angstroms
saveamberparm LIG 4TKG.lig.wat.leap.prm7 4TKG.lig.wat.leap.rst7    #Save the ligand water phase AMBER topology and coordinate files
quit

 tleap –s –f tleap.lig.in > tleap.lig.out

The atom types which used in '4TKG.lig.mol2' file are not recognized.

So the name of the atoms has to be redefined based on the gaff force field, antechamber can be used to solve the error.

 antechamber –i 4TKG.lig.mol2 –fi mol2 –o 4TKG.lig.gaff.mol2 –fo mol2

Change the content of tleap.lig.in to (gaff file)

Re-run tleap by:

 tleap –s –f tleap.lig.in > tleap.lig.out

So now run parmchk to fix the missing parameters

 parmchk –i 4TKG.lig.gaff.mol2 –f mol2 –o 4TKG.lig.ante.frcmod

Change the content of tleap.lig.in to (load amber params) run tleap.lig.in, and all the errors are fixed.

 cp ~/../../tleap.rec.in

• tleap.rec.in

source leaprc.ff14SB                                                  #Load a force field
REC = loadpdb ../001.chimera/4TKG.rec.noH.ZIN.pdb                     #Conformation file for receptor 
saveamberparm REC 4TKG.rec.gas.leap.prm7 4TKG.rec.gas.leap.rst7       #Save the receptor gas phase AMBER topology and coordinate file
solvateBox REC TIP3PBOX 10.0                                          #Solvate the lig using TIP3P, solvent box radii 10 angstroms
saveamberparm REC 4TKG.rec.wat.leap.prm7 4TKG.rec.wat.leap.rst7       #Save the receptor water phase AMBER topology and coordinate file
quit

tleap –s –f tleap.rec.in > tleap.rec.out

For the missing parameters, cp the two ionparam from ../000.amberfiles/

 cp ../000.amberfiles/ions.lib
 cp ../000.amberfiles/ions.frcmod

Change the content of the tleap.rec.in to (load param)

run tleap.rec.in, and all the errors are fixed. (explain)

 cp ~/../../tleap.com.in
 tleap –s –f tleap.com.in > tleap.com.out

This whole process can be run on a cluster by using the following script:

III. Simulation using pmemd

PMEMD

Agatha, Beilei

MD input file1 01mi.in:

mkdir 01mi.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 01mi.in 02md.in

02md.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 03mi.in, 04mi.in, 05mi.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,

04mi.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 = 0.1,

05mi.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 = 0.05,

06md.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 = 1.0,

07md.in: equilibration
&cntrl
  imin = 0, ntx = 5, irest = 1, 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 = 0.5,

08md.in: equilibration
&cntrl
  imin = 0, ntx = 5, irest = 1, 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-223 & @CA,C,N', restraint_wt = 0.1,

09md.in: equilibration
&cntrl
  imin = 0, ntx = 5, irest = 1, 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-223 & @CA,C,N', restraint_wt = 0.1,

10md.in: production
&cntrl
  imin = 0, ntx = 5, irest = 1, nstlim = 500000,
  temp0 = 298.15, tempi = 298.15, ig = 71287,
  ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
  ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
  ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
  cut = 8.0, iwrap = 1,
  ntr = 1, nscm = 100,
  restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,

11md.in: production
&cntrl
  imin = 0, ntx = 5, irest = 1, nstlim = 500000,
  temp0 = 298.15, tempi = 298.15, ig = 71287,
  ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
  ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
  ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
  cut = 8.0, iwrap = 1,
  ntr = 1, nscm = 100,
  restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,

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