2017 AMBER tutorial with 4qmz

From Rizzo_Lab
Jump to: navigation, search

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

Introduction

Amber -Assisted Model Building with Energy Refinement - is a multi-program suite for macromolecular simulations.Amber is distributed in two parts: AmberTools15 and Amber14.Amber14 is the most recent version of the software and it includes new force fields such as ff14SB. In addition, in this release, more features from sander have been added to pmemd for both CPU and GPU platforms, including performance improvements, and support for extra points, multi-dimension replica exchange, a Monte Carlo barostat, ScaledMD, Jarzynski sampling, explicit solvent constant pH, GBSA, and hydrogen mass repartitioning. Support is also included for the latest Kepler, Titan and GTX7xx GPUs expanded options for Poisson-Boltzmann solvation calculations, accelerated molecular dynamics, additional features in sander pmemd code, and expanded methods for free energy calculations. Our lab is set up with Ambe r14 and the latest update of AmberTools15 which contains the programs such as antechamber and tleap to set up your simulation. The Amber 14 Manualis available to get started with using Amber14. You can search the document for keywords such as "tleap" if you use Adobe Acrobat to view the file. Additionally, AmberTools Reference Manualis another reference for the programs available under Amber tools.

Here below are some of the programs available in both Amber and AmberTools:

1.LEaP: LEaP is an X-windows-based program that provides for basic model building and Amber coordinate parameter/topology input file creation. It includes a molecular editor which allows for building residues and manipulating molecules.

2.ANTECHAMBER: This program suite automates the process of developing force field descriptors for most organic molecules. It starts with structures (usually in PDB format), and generates files that can be read into LEaP for use in molecular modeling. The force field description that is generated is designed to be compatible with the usual Amber force fields for proteins and nucleic acids.

3.SANDER: Sander is short for Simulated annealing with NMR-derived energy restraints. This allows for NMR refinement based on NOE-derived distance restraints, torsion angle restraints, and penalty functions based on chemical shifts and NOESY volumes. Sander is also the "main" program used for molecular dynamics simulations, and is also used for replica-exchange, thermodynamic integration, and potential of mean force (PMF) calculations. Sander also includes QM/MM capability.

4.PMEMD: This is an extensively-modified version (originally by Bob Duke) of the sander program, optimized for periodic, PME simulations, and for GB simulations. It is faster than sander and scales better on parallel machines.

5.PTRAJandCPPTRAJ: These are used to analyze MD trajectories, computing a variety of things, like RMS deviation from a reference structure, hydrogen bonding analysis, time-correlation functions, diffusional behavior, and so on.

6.MM_PBSA andMM_PBSA.py: These are scripts that automate post-processing of MD trajectories, to analyze energetics using continuum solvent ideas. It can be used to break energies energies into "pieces" arising from different residues, and to estimate free energy differences between conformational basins.

7.NAB: Originally named as "nucleic acid builder", NAB is a specialized language for writing programs that manipulate molecules and carry out molecular mechanics or distance-geometry based modeling. NAB provides and interface to Poisson-Boltzmann and RISM integral-equation solvent models. The "amberlite" package uses NAB to study protein-ligand interaction energetics. There is also a mailing list available as an additional resource. What you can do with it is: you document your questions and sent to this mail address, some specialists of Amber will be assigned to reply your email and help you.

4qmz:

Organizing Directories It makes things easier to organize your files in a clean and logical way. The following directory structure and naming scheme is a convenient way to organize your files. We could make these directories first before doing anything further

  ~username/AMS536-Spring2016/Amber_Tutorial/000.parameters/
                                             001.tleap_build/ 
                                             002.equilibration/       
                                             003.production/
                                             004.analysis/

II. Structural Preparation

Before we can submit our job for equilibration / production, we must prepare it for simulation. Ensure you have your environmental variables set to include your AMBERHOME, CUDA_HOME, and anaconda. The first step is to parametrize your ligand using antechamber and check for any missing force field parameters; run using this command:

  /path/to/amber/antechamber -i /path/to/file/sunitinib_sybyl.lig.mol2 -fi mol2 -o sunitinib_lig.am1bcc.mol2 -fo mol2 -at gaff -c bcc -rn LIG -nc 1
  /path/to/amber/parmchk -i sunitinib_lig.am1bcc.mol2 -f mol2 -o sunitinib_lig.am1bcc.frcmod

The next step is to create the topology and coordinate files for the simulation through tleap. First create a file called tleap.in with the following:

source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
set default PBradii mbondi3

gascomplex = combine {rec lig}
savepdb gascomplex 4qmz_sunitinib.gas.complex.pdb
saveamberparm gascomplex 4qmz_sunitinib.gas.complex.prmtop 4qmz_sunitinib.gas.complex.rst7
saveamberparm rec 4qmz_sunitinib.gas.receptor.prmtop 4qmz_sunitinib.gas.receptor.rst7
saveamberparm lig 4qmz_sunitinib.gas.ligand.prmtop 4qmz_sunitinib.gas.ligand.rst7
solvcomplex = combine {rec lig}
solvateoct solvcomplex TIP3PBOX 12.0

This will load in the pertinent libraries, solvate your system in a truncated octahedron 12.0 A out, and will give you topologies/rst/pdb's for the non-solvated (gas) ligand, gas receptor, solvated ligand, solvated receptor, gas-phase complex, and solvated complex. VISUALIZE THESE IN VMD TO ENSURE THEY LOOK REALISTIC! Many times errors down the line are resultant from skipping the inspection after running through tleap.

III. Structural Equilibration

We first start by performing a series of minimization and short MD simulation steps for the receptor and ligand. The purpose of these steps is to remove steric clashes (minimization) and allow the receptor:ligand to equilibrate into a low energy state (MD simulation). After minimization, we can perform MD simulations in which we could record the trajectories to predict the folded state of protein.

1) In the first minimization step, we use a large restraint_wt on all heavy atoms (5.0 kcal/molA^2). Thus, only hydrogens are allowed to relax and only steric clashes involving hydrogens are removed. Create the input file 01.min.in:

Minmize all the hydrogens
 &cntrl
  imin=1,           ! Minimize the initial structure
  maxcyc=5000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-288 & !@H=", ! atoms to be restrained
  restraint_wt=5.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format

imin = 1: flag to do minimization only

ntx = 1: read coordinates only from rst file

ntc = 1: flag for shake to perform bond length constraints. ntc =1 is the default, meaning SHAKE is not performed. In other words, minimization is performed on every bond.

ntb = 1: flag for constant volume

ntp = 0: NO constant pressure

ntwx = 1000: saving snapshot to a trajectory file every 1000 steps.

ntwe = 0: human readable energy file is appended every 0 steps.

ntpr = 1000: newest calculated energy is stored in .info file every 1000 steps.

cut: cutoff for energy evaluation (Angstroms)

ntr = 1: flag for restraining specified atoms

restraintmask = ':1-224 & !@H='; String that specifies restrained atoms (all atoms except hydrogens are restrained)

restraint_wt = 5.0; restraint strength (kcal/molA^2)


2) The second step is a short MD simulation to bring the receptor:ligand system into a more equilibrated state. Create an input file 02.equil.mdin by copying and editing 01.min.mdin. cp 01.min.mdin 02.equil.mdin

MD simualation
&cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntb=2,            ! Constant Pressure
 ntc=1,            ! No SHAKE on bonds between hydrogens
 dt=0.001,         ! Timestep (ps)
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen 
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-288 & !@H=", ! atoms to be restrained
 restraint_wt=5.0, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on


imin = 0: flag to do MD simulation

ntc = 2: flag for SHAKE to perform bond length constraints involving hydrogen

ntb = 2: flag for constant pressure

ntp = 1: constant pressure

ntr = 1: flag for restraining specified atoms

restraintmask = ':1-224 & !@H='; String that specifies restrained atoms (all atoms except hydrogens are restrained)

restraint_wt = 5.0; restraint strength (kcal/molA^2)

3) Create more minimization input files with decreasing restraint_wt (03.min.mdin, 04.min.mdin, 05.min.mdin). Create more MD simulation input files with different restraints (06.equil.mdin, 07.equil.mdin, 08.equil.mdin, 09.equil.mdin:

03.min.mdin:

Minmize all the hydrogens
&cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-288 & !@H=", ! atoms to be restrained
 restraint_wt=2.0, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format

04.min.mdin:

Minmize all the hydrogens
&cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-288 & !@H=", ! atoms to be restrained
 restraint_wt=0.1, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format

05.min.mdin:

 Minmize all the hydrogens
&cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-288 & !@H=", ! atoms to be restrained
 restraint_wt=0.05, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format


06.equil.mdin:

MD simualation
&cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntb=2,            ! Constant Pressure
 ntc=1,            ! No SHAKE on bonds between hydrogens
 dt=0.001,         ! Timestep (ps)
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen 
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-288 & !@H=", ! atoms to be restrained
 restraint_wt=1.0, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on


07.equil.mdin:

MD simulation
&cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen 
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-288 & !@H=", ! atoms to be restrained
 restraint_wt=0.5, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on


08.equil.mdin:

MD simulations
&cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen 
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-287@CA,C,N", ! atoms to be restrained
 restraint_wt=0.1, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on


09.equil.mdin:

MD simulations
&cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen 
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-287@CA,C,N", ! atoms to be restrained
 restraint_wt=0.1, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on

It is ideal to write a script that will send these 9 files to the queue as they are finished. An example of a script that can do this is md.equilibration.sh:

#!/bin/sh
#PBS -N 4qmz_equilibration
#PBS -l walltime=02:00:00
#PBS -l nodes=2:ppn=28
#PBS -j oe
#PBS -q long
cd $PBS_O_WORKDIR
echo "Started Equilibration on `date` "
do_parallel="mpirun pmemd.MPI"
prmtop="../001.tleap_build/4qmz_sunitinib.wet.complex.prmtop"
coords="../001.tleap_build/4qmz_sunitinib.wet.complex"
MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)
for input in ${MDINPUTS[@]}; do
  $do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf 
${input}.info -r ${input}.rst7
  coords=$input
done
echo "Finished Equilibration on `date` "


At this point, after the script has been submitted to the queue via:

qsub md.equilibration.sh

There will be 4 new files for each mdin submitted. For the first equilibration step, the files output are as follows:

01.min.info
01.min.mdout
01.min.trj
01.min.rst7

It is always a good idea to scroll through the outputs to ensure there are not any obvious error messages. Additionally, it is a good idea to visualize the last equilibration trajectory using VMD to ensure your structure is reasonable before beginning the production MD run.

IV. Simulation Production

Change your directory to 003.production. Once here, create two new directories, 001.restrained and 002.unrestrained. The simulation will be run twice, once with a restraint on the backbone and once with no restraint.

Move to the 001.restrained directory and create an input file called 10.prod.mdin using the command

  vi 10.prod.mdin

write the following in the file

  MD simulations
 &cntrl
  imin=0,           ! Perform MD
  nstlim=5000000,   ! Number of MD steps
  ntx=5,            ! Positions and velocities read formatted
  irest=1,          ! Restart calculation
  ntc=2,            ! SHAKE on for bonds with hydrogen
  dt=0.002,         ! Timestep (ps)
  ntb=2,            ! Constant Pressure
  ntp=1,            ! Isotropic pressure scaling
  barostat=1        ! Berendsen
  taup=0.5          ! Pressure relaxtion time (ps)
  ntf=2,            ! No force evaluation for bonds with hydrogen
  ntt=3,            ! Langevin thermostat
  gamma_ln=2.0      ! Collision Frequency for thermostat
  ig=-1,            ! Random seed for thermostat
  temp0=298.15      ! Simulation temperature (K)
  ntwx= 2500,       ! Write to trajectory file every ntwx steps
  ntpr= 2500,       ! Print to mdout every ntpr steps
  ntwr= 5000000,    ! Write a restart file every ntwr steps
  cut=8.0,          ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-287@CA,C,N", ! atoms to be restrained
  restraint_wt=0.1, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  iwrap=1,          ! iwrap is turned on
/

In the 002.unrestrained directory, a similar input file should be made. Use the command

 vi 10.prod.mdin

to make a new input file. In it, type

  MD simulations
 &cntrl
  imin=0,           ! Perform MD
  nstlim=1000000    ! Number of MD steps
  ntx=5,            ! Positions and velocities read formatted
  irest=1,          ! Restart calculation
  ntc=2,            ! SHAKE on for bonds with hydrogen
  dt=0.002,         ! Timestep (ps)
  ntb=2,            ! Constant Pressure
  ntp=1,            ! Isotropic pressure scaling
  barostat=1        ! Berendsen
  taup=0.5          ! Pressure relaxtion time (ps) 
  ntf=2,            ! No force evaluation for bonds with hydrogen
  ntt=3,            ! Langevin thermostat
  gamma_ln=2.0      ! Collision Frequency for thermostat
  ig=-1,            ! Random seed for thermostat
  temp0=298.15      ! Simulation temperature (K)
  ntwx= 2500,       ! Write to trajectory file every ntwx steps
  ntpr= 2500,       ! Print to mdout every ntpr steps
  ntwr= 5000000,    ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  iwrap=1,          ! iwrap is turned on
/

IV. Simulation Analysis

Once the production simulation finishes successfully, it's good to visually check the trajectory (in vmd) and apply some post processing. Towards that end create an analysis directory containing three subdirectories:

   mkdir 004.analysis
   cd 004.analysis/
   mkdir 001.rmsd  002.hbond  003.mmgbsa


We will now remove the water molecules and create a trajectory aligned to the CA atoms of the crystal structure. To do this create a cpptraj input file named "cpptraj.strip.wat.in" containing the following:

   #!/usr/bin/sh
   #read in trajectory
   trajin ../../003.production/001.restrained/10.prod.trj
   #read in reference
   reference ../../001.tleap_build/4qmz_sunitinib.wet.complex.rst7
   #compute rmsd and align CA to the crystal structure
   rmsd rms1 reference :1-287@CA
   #strip Solvent
   strip :WAT:Na+:Cl-
   #create gas-phase trajectory
   trajout 4qmz_sunitinib.stripfit.restrained.gas.trj nobox

and run it:

   cpptraj -p ../../001.tleap_build/4qmz_sunitinib.wet.complex.prmtop -i cpptraj.strip.wat.in

Take a look at the trajectory in vmd.

Next is to compute ligand rmsd versus the crystal structure pose and histogram the RMSD. Create a new cpptraj input file named "cpptraj.rmsd.lig.in":

   #!/usr/bin/sh
   #trajin the trajectory
   trajin 4qmz_sunitinib.stripfit.restrained.gas.trj
   #read in the reference
   reference ../../001.tleap_build/4qmz_sunitinib.gas.complex.rst7
   #compute the RMSD (do not fit the internal geometries first, included rigid body motions
   #and convert the frames to ns (framenum*.005)
   rmsd rms1 ":288&!(@H=)" nofit mass out 4qmz_sunitinib.lig.restrained.rmsd.nofit.dat time .005
   #histogram the nofit rmsd
   histogram rms1,*,*,.1,* norm out 4qmz_sunitinib.lig.restrained.rmsd.nofit.histogram.dat


and run it:

   cpptraj -p ../../001.tleap_build/4qmz_sunitinib.gas.complex.prmtop -i cpptraj.rmsd.lig.in

This will output a file with the frame number and ligand rmsd in that frame compared to the crystal structure ligand pose. If the rmsd does not diverge too much throughout the simulation, it means that the ligand is stable in the binding pocket.


Next is to analyze hydrogen bonding between the ligand and the enzyme residues.

  cd 004.analysis/002.hbond
  #Compute hydrogen bonds
  cpptraj -p ../../001.tleap_build/4qmz_sunitinib.wet.complex.prmtop -i cpptraj.hbond.in

The input file "cpptraj.hbond.in" contains the following lines:

   #!/usr/bin/sh
   #read in trajectory
   trajin ../../003.production/001.restrained/10.prod.trj
   #wrap everything into one periodic cell
   autoimage
   #compute intra and water mediated hydrogen bonds
   hbond hb1 :1-288 out 4qmz_sunitinib.hbond.out avgout 4qmz_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT nointramol bridgeout 4qmz_sunitinib.bridge-water.dat dist 3.0 angle 140


Next we will calculate the binding free energy (delta G binding) of the ligand to the enzyme using the MMGBSA method (in AMBER it is done via a python script MMGBSA.py).

Go into the directory "003.mmgbsa/"

Create the input file "mmgbsa.in" with the following commands:

   mmgbsa 4QMZ analysis
   &general
     interval=1, netcdf=1,
     keep_files=0,
   
   /
   &gb
     igb=8,
     saltcon=0.0, surften=0.0072,
     surfoff=0.0, molsurf=0,
   /
   &nmode
     drms=0.001, maxcyc=10000,
     nminterval=250, nmendframe=2000,
     nmode_igb=1,
   /

The MMGBSA simulation should be submitted to the cluster using the submission script:

   qsub mmgbsa.sh

The above submission script "mmgbsa.sh" would contain:

   #!/bin/bash
   #PBS -l walltime=35:00:00
   #PBS -l nodes=1:ppn=28
   #PBS -N 4qmz_mmgbsa
   #PBS -V
   #PBS -j oe
   #PBS -q long
   
   cd $PBS_O_WORKDIR
   
   #Define topology files 
   solv_prmtop="../../001.tleap_build/4qmz_sunitinib.wet.complex.prmtop"
   complex_prmtop="../../001.tleap_build/4qmz_sunitinib.gas.complex.prmtop"
   receptor_prmtop="../../001.tleap_build/4qmz_sunitinib.gas.receptor.prmtop"
   ligand_prmtop="../../001.tleap_build/4qmz_sunitinib.gas.ligand.prmtop"
   trajectory="../../003.production/001.restrained/10.prod.trj"
   
   #create mmgbsa input file
   cat >mmgbsa.in<<EOF
   mmgbsa HIVgp41 analysis
   &general
     interval=1, netcdf=1,
     keep_files=0,
   
   /
   &gb
     igb=8,
     saltcon=0.0, surften=0.0072,
     surfoff=0.0, molsurf=0,
   /
   &nmode
     drms=0.001, maxcyc=10000,
     nminterval=250, nmendframe=2000,
     nmode_igb=1, 
   /
   EOF
   
   MMPBSA.py -O -i mmgbsa.in \
              -o  4qmz_sunitinib.mmgbsa.results.dat \
              -eo 4qmz_sunitinib.mmgbsa.per-frame.dat \
              -sp ${solv_prmtop} \
              -cp ${complex_prmtop} \
              -rp ${receptor_prmtop} \
              -lp ${ligand_prmtop} \
               -y ${trajectory}


The output file "4qmz_sunitinib.mmgbsa.results.dat" shows the calculations and the resulting computed energy value; the binding energy (delta G binding) is given at the bottom of the file (values can range roughly between -15 and -25 kcal/mol +- 5 kcal/mol) for this system.

THE END