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

Revision as of 13:07, 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.mol.prep folder, add the following files from your docktutorial directory:

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

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

Open 1BJU.rec.noH.pdb

Before running tleap, make a new directory: mkdir 002.ANTE.TLEAP

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

tleap.lig.in

source leaprc.ff14SB                                               #Load a force field
source leaprc.gaff                  
LIG = loadmol2 1BJU.lig.mol2                                       #Conformation file
saveamberparm LIG 1BJU.lig.gas.leap.prm7 1BJU.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 1BJU.lig.wat.leap.prm7 1BJU.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 '1BJU.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 1BJU.lig.mol2 –fi mol2 –o 1BJU.lig.gaff.mol2 –fo mol2

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

Re-run tleap:

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

So now run parmchk to fix the missing parameters

 parmchk –i 1BJU.lig.gaff.mol2 –f mol2 –o 1BJU.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/1BJU.rec.noH.pdb                         #Conformation file for receptor 
saveamberparm REC 1BJU.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 1BJU.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

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

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

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

@@ Line 22: / Line 22: @@
 BJU.lig.mol2
 BJU.rec.noH.pdb
-BJU.rec.noH.pdb
@@ Line 48: / Line 48: @@
 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 1BJU.lig.mol2 –fi mol2 –o 4TKG.lig.gaff.mol2 –fo mol2
+   antechamber –i 1BJU.lig.mol2 –fi mol2 –o 1BJU.lig.gaff.mol2 –fo mol2
 Change the content of tleap.lig.in to
 (gaff file)
@@ Line 65: / Line 65: @@
 *tleap.rec.in
   source leaprc.ff14SB                                                  #Load a force field
-  REC = loadpdb ../001.chimera/1BJU.rec.noH.ZIN.pdb                     #Conformation file for receptor
+  REC = loadpdb ../001.chimera/1BJU.rec.noH.pdb                         #Conformation file for receptor
   saveamberparm REC 1BJU.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
@@ Line 73: / Line 73: @@
 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==

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

Revision as of 13:07, 12 April 2016

Contents

I. Introduction

II. Structural Preparation

Antechamber, Parmchk, tLeap

Antechamber, Parmchk, tLeap

III. Simulation using pmemd

PMEMD

IV. Simulation Analysis

Ptraj

MM-GBSA Energy Calculation

V. Frequently Encountered Problems

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Rizzo Lab

Courses

Toolbox