Difference between revisions of "2021 AMBER tutorial 1 with PDBID 1HW9"

From Rizzo_Lab
Jump to: navigation, search
(→‎MMGBSA)
 
(4 intermediate revisions by the same user not shown)
Line 474: Line 474:
  
 
== '''MD Analysis''' ==
 
== '''MD Analysis''' ==
 +
Now that the MD trajectory and results have been generated, we can analyze them a bit further. Our analysis will cover three topics, RMSD,Hydrogen bonds, and MMGBSA.
 +
 +
 +
=== '''RMSD''' ===
 +
Here, we'll be determining how far the ligand and receptor move when interacting. This will require three input files. Our first input file needs to remove all water and ions from the production trajectory.
 +
 +
vi cpptraj_strip_wat.in
 +
 +
Place the following parameters in this input file:
 +
 +
#!/usr/bin/sh
 +
parm ../003_leap/1HW9.wet.complex.parm7
 +
#read in trajectory
 +
trajin ../004_production/10.prod.trj
 +
#read in reference
 +
reference ../003_leap/1HW9.wet.complex.rst7
 +
#compute RMSD + align CA to crystal structure
 +
rmsd rms1 reference :1-263@CA
 +
#strip solvent
 +
strip :WAT:Na+:Cl-
 +
#create gas-phase trajectory
 +
trajout 1HW9_stripfit_restrained_gas.trj nobox
 +
 +
Your reference mask should include both the ligand and receptor. In order to execute:
 +
 +
cpptraj -i cpptraj_strip_wat.in
 +
 +
Next, we'll measure the RMSD of the ligand, separate from the receptor. In order to do this, create the following input file:
 +
 +
vi cpptraj_rmsd_lig.in
 +
 +
The input file should have the following parameters:
 +
 +
#!/usr/bin/sh
 +
#trajin the trajectory
 +
trajin 1HW9_stripfit_restrained_gas.trj
 +
#read in reference
 +
reference ../003_leap/1HW9.gas.complex.rst7
 +
#compute RMSD (do not fit internal geometris first, included rigid body motions, convert frames to ns (framenum*.005)
 +
rmsd rms1 ":254&!(@H=)" nofit mass out 1HW9_lig_restrained_rmsd_nofit.dat time .005
 +
#histogram the nofit rmsd
 +
histogram rms1,*,*,.1,* norm out 1HW9_lig_restrained_rmsd_nofit_histogram.dat
 +
 +
Make sure that your reference mask only includes the ligand this time. In order to execute:
 +
 +
cpptraj -p ../003_leap/1HW9.complex.parm7 -i cpptraj_rmsd_lig.in
 +
 +
Finally, we'll measure the RMSD of the receptor, separate from the ligand. In order to do this. create the following input file:
 +
 +
vi cpptraj_rmsd_rec.in
 +
 +
The file should contain the following parameters:
 +
 +
#!/usr/bin/sh
 +
#trajin the trajectory
 +
trajin 1HW9_stripfit_restrained_gas.trj
 +
#read in reference
 +
reference ../003_leap/1HW9.gas.complex.rst7
 +
#compute RMSD (do not fit internal geometries first, included rigid body motions, convert frames to ns (framenum*.005)
 +
rmsd rms1 ":1-253&!(@H=)" nofit mass out 1HW9_rec_restrained_rmsd_nofit.dat time .005
 +
#histogram the nofit rmsd
 +
histogram rms1,*,*,.1,* norm out 1HW9_rec_restrained_rmsd_nofit_histogram.dat
 +
 +
Make sure that your reference mask only includes the receptor this time. In order to execute:
 +
 +
cpptraj -p ../003_leap/1HW9.gas.complex.parm7 -i cpptraj_rmsd_rec.in
 +
 +
 +
=== '''Hydrogen Bonding''' ===
 +
We can also determine if our ligand forms any hydrogen bonds with the receptor protein. In order to do this, we'll need to create the following input file:
 +
 +
vi cpptraj_hbond.in
 +
 +
The input file should have the following parameters:
 +
 +
!/usr/bin/sh
 +
#read in trajectory
 +
trajin ../004_production/10_prod.trj
 +
#wrap everything in one periodic cell - could cause problems, may comment out #autoimage if problems later
 +
autoimage
 +
#compute intra + water-mediated H-bonds
 +
hbond hb1 :1-263 out 1HW9_calcipotriol_hbond.out avgout 1HW9_calcipotriol_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
 +
nointramol brid\
 +
geout 1HW9_calcipotriol_bridge-water.dat dist 3.0 angle 140
 +
 +
The hbobnd hb1 line should contain the residue numbers corresponding to your protein and ligand. In order to execute the script:
 +
 +
cpptraj -p ../003_leap/1HW9_wetcomplex.parm7 -i cpptraj_hbond.in
 +
 +
 +
=== '''MMGBSA''' ===
 +
In this section, we'll be estimating the free energy of binding of our ligand and receptor. This is performed by analyzing both the unsolvated and solvated ligand and receptor separately. In order to get started, we'll need the following input file:
 +
 +
vi mmgbsa.in
 +
 +
Your input file should contain the following parameters:
 +
 +
mmgbsa 1HW9 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,
 +
  /
 +
 +
In order to run the input file, we'll need to create a separate script.
 +
 +
vi mmgbsa.sh
 +
 +
This file should have the following parameters:
 +
 +
#bin/bash
 +
#SBATCH --time=2-00:00:00
 +
#SBATCH --nodes=2
 +
#SBATCH --ntasks=28
 +
#SBATCH --job-name=1HW9_MMGBSA
 +
#SBATCH --output=1HW9_MMGBSA.log
 +
#SBATCH -p extended-28core
 +
#Define topology files
 +
solv_parm="../003_leap/1HW9.wet.complex.parm7"
 +
complex_parm="../003_leap/1HW9.complex.parm7"
 +
receptor_parm="../003_leap/1HW9.gas.receptor.parm7"
 +
lig_parm="../003_leap/1HW9.gas.ligand.parm7"
 +
trajectory="10_prod.trj"
 +
  MMPBSA.py -O -i mmgbsa.in \
 +
      -o 1S19_mmgbsa_results.dat \
 +
      -eo 1S19_mmgbsa_perframe.dat \
 +
      -sp ${solv_parm} \
 +
      -cp ${complex_parm} \
 +
      -rp ${receptor_parm} \
 +
      -lp ${lig_parm} \
 +
      -y ${trajectory}
 +
 +
You'll need to submit this to the queue:
 +
 +
sbatch mmgbsa.sh

Latest revision as of 00:59, 7 May 2021

Introduction

In this tutorial, we'll be exploring the wonders of molecular dynamics through the use of AMBER 16. AMBER 16 we'll help study the binding characteristics of the ligand to the target protein!


Directory Setup

It's important to stay organized during this process as we'll be generating a number of files during the process. Overall, we'll be working within five directories. They should be created as follows: 

mkdir 001_structure
mkdir 002_parameters
mkdir 003_leap
mkdir 004_equil
mkdir 005_production

Before proceeding, double-check that all directories have correctly been made.


1HW9 Structures

Let's move into our directory titled 001_structure. 

cd 001_structure

Receptor

We need to retrieve a fresh copy of the receptor before proceeding forward. Even if you've done the DOCK tutorial, you should still retrieve a fresh .pdb just in case Chimera did something wrong during the preparation process. In order to do so, head over to the PDB website and type in 1HW9 into the header. It should take you straight to the webpage. Download the system as a pdb and open it in Chimera. Our focus will be on Chain A of the complex. In order to visualize the chain: 

 Select -> Chain -> A

Notice, however, that this does not remove Chains B, C, and D from our view. We need to perform this action manually. 

 Select -> Invert (Selected Molecules) -> Actions -> Atoms/Bonds -> Delete

After deletion, the only thing remaining should be Chain A of the original file. It is within your best interest to save this session as you'll need to come back to this file when preparing the ligand. Doing so will prevent you from having to do the above deletion step again. 

 File -> Save Session As ->

Now that the session has been saved, we can go ahead and isolate the receptor. While we've isolated Chain A, you'll notice that Simvastatin is still located on the molecule. In order to remove everything besides HMGR, we need to do the following: 

 Select -> Residue -> ADP -> Actions -> Atoms/Bonds -> Delete
 Select -> Residue -> HOH -> Actions -> Atoms/Bonds -> Delete
 Select -> Residue -> SIM -> Actions -> Atoms/Bonds -> Delete

With the receptor isolated, go ahead and save the file as 

1HW9_fresh.pdb

Ligand

Open the session in Chimera we previously saved that consisted of Chain A with both the receptor and Simvastatin ligand. This time around, we'll be using the session file to isolate the ligand instead of the receptor. 

 Select -> Residue -> SIM -> Select -> Invert (Selected Molecules) -> Actions -> Atoms/Bonds -> Delete

We'll need to properly protonate the ligand in order for the AMBER calculations to function: 

Tools -> Structure Editing -> Add H
Tools -> Structure Editing -> Add Charge -> (have Amber ff14SB and AM1-BCC selected) -> Ok

Once done, make sure to save your file as: 

1HW9_lig_wH.mol2

NOTE: Before proceeding, make sure both the ligand and receptor file you just created have been copied from your local computer to the 001_structure directory in your Seawulf account.


Simulation Parameters

Let's move into the second directory: 

cd 002_parameters

Here, we'll be generating force field parameters for our ligand so that AMBER can use them in later calculations. We'll be doing this using antechamber and parmchk2: 

antechamber -i ../001_structure/1HW9_lig_wH.mol2 -fi mol2 -o 1HW9_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -1

Notice the -nc flag at the end of the command line. Make sure that this value corresponds to the protonation state of your ligand! Once the output file is generated, we need to run parmch2 to generate a frcmod file: 

parmchk2 -i 1HW9_ligand_antechamber.mol2 -f mol2 -o 1HW9_ligand.am1bcc.frcmod

TLeap System

Let's now move into the third directory: 

cd 003_leap

In this stage, we'll be creating files to simulate the ligand and receptor as one system. This way, calculations like binding affinities and energies can be performed on the system as a whole. Fortunately, tleap will do this for us by generating parameter (parm7) and restart (rst7) files. In order to get started, we'll need to create the input file. 

vi leap.in

Copy the following code into your new input file: 

#!/usr/bin/sh 
###load protein force field
source leaprc.protein.ff14SB
###load GAFF force field (for our ligand)
source leaprc.gaff
###load TIP3P (water) force field
source leaprc.water.tip3p
###load ions frcmod for the tip3p model
loadamberparams frcmod.ionsjc_tip3p
###needed so we can use igb=8 model
set default PBradii mbondi3
###load protein pdb file
rec=loadpdb ../001_structure/1HW9_fresh.pdb
###load ligand frcmod/mol2
loadamberparams ../002_parameters/1HW9_ligand.am1bcc.frcmod
lig=loadmol2 ../002_parameters/1HW9_ligand_antechamber.mol2
###create gase-phase complex
gascomplex= combine {rec lig}
###write gas-phase pdb
savepdb gascomplex 1HW9.gas.complex.pdb
###write gase-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 1HW9.complex.parm7 1HW9.gas.complex.rst7
saveamberparm rec 1HW9.gas.receptor.parm7 1HW9.gas.receptor.rst7
saveamberparm lig 1HW9.gas.ligand.parm7 1HW9.gas.ligand.rst7
###create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###solvate the system
solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0
#write solvated pdb file
savepdb solvcomplex 1HW9.wet.complex.pdb
###check the system
charge solvcomplex
check solvcomplex
###write solvated toplogy and coordinate file
saveamberparm solvcomplex 1HW9.wet.complex.parm7 1HW9.wet.complex.rst7
quit

In order to now execute: 

tleap -f leap.in

Once this is done, you should copy your parameter and coordinate files to your local computer and view them in Chimera. You want to make sure that everything looks correct before proceeding forward. For example, your protein ligand complex should not be outside the water molecule box that you just created. It's little things like this that can create headaches later on. In order to view in Chimera: 

Tools -> MD/Ensemble Analysis -> MD Movie --> Choose appropriate parameter and coordinate files to load

Equilibration

Let's now move into our fourth directory: 

cd 004_equil

We now need to perform energy minimization and equilibration similarly to how it was performed in the DOCK tutorial. Unfortunately, this also means we need to create 9 input files as this is a 9 step procedure: You're input files should be named and contain the following parameters: 

vi 01.min.mdin

Minimize all the hydrogens
&cntrl
imin=1,           ! Minimize the initial structure
ntmin=2,         ! Use steepest descent Ryota Added
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="!@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
/

vi 02.equil.mdin

MD simulation
&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=":!@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
/

vi 03.min.mdin

Minimize 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="!@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
/

vi 04.min.mdin

Minimize 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="!@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
/

vi 05.min.mdin

Minimize 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="!@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 
/

vi 06.equil.mdin

MD simulation
&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="!@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
/

vi 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="!@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
/

vi 08.equil.mdin *See below for specific comments on this input

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-433@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
/

vi 09.equil.mdin *See below for comments on this input

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-433@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
/

NOTE: The restraint mask lines in files 8 and 9 have to be specific to your system. These numbers correspond to residues in your system. These residue numbers can be found in the wet.complex.pdb file we generated back in step 3. File 8's restraint mask should include both the receptor and the ligand. You'll know the residue number of the ligand as the pdb will have is labeled as LIG. File 9's restraint mask is only the receptor so the value should actually be 1 less.

In order to run all 9 scripts, an input file has to be created: 

vi mdequilibration.sh

In this input file, you should have the following code: 

#!/bin/sh
#SBATCH --job-name=1HW9_equilibration
#SBATCH --ntasks-per-node=40
#SBATCH --nodes=2
#SBATCH --time=8:00:00
#SBATCH -p long-40core  
echo "started Equilibration on 'date' "
parm7="../003_leap/1HW9.wet.complex.parm7"
coords="../003_leap/1HW9.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

   pmemd -O -i ${input}.mdin -o ${input}.mdout -p $parm7 -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7 
coords=$input

done

echo " Finished equilibration on 'date' "

In order to submit this script: 

sbatch mdequilibration.sh


Production

We can now move into our fifth and final director: 

cd 005_production

Here, we'll be performing the production run and any further analysis. First, we need to generate the production input file to create our MD trajectory. 

vi 10.prod.mdin

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-432@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
/

NOTE: Rember to modify the restrain mask! This one should only include the receptor, not the ligand.

We now need to be able to run the script: 

vi mdproduction.sh

!/bin/sh
#SBATCH --job-name=1HW9_prod
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=2
#SBATCH --time=7-00:00:00
#SBATCH -p extended-24core
echo "started Equilibration on 'date' "
parm7="../003_leap/1HW9.wet.complex.parm7"
coords="../004_equil/09.equil"
MDINPUTS=( 10.prod)
for input in ${MDINPUTS[@]}; do
  pmemd -O -i ${input}.mdin -o ${input}.mdout -p ${parm7} -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
  coords=$input
done
echo "Finished Equilibration on `date` "

We'll need to submit this job as well. 

sbatch mdproduction.sh

Be prepared for this to take several days. It most certainly is NOT a fast process so be sure to plan accordingly.


MD Analysis

Now that the MD trajectory and results have been generated, we can analyze them a bit further. Our analysis will cover three topics, RMSD,Hydrogen bonds, and MMGBSA.


RMSD

Here, we'll be determining how far the ligand and receptor move when interacting. This will require three input files. Our first input file needs to remove all water and ions from the production trajectory. 

vi cpptraj_strip_wat.in

Place the following parameters in this input file: 

#!/usr/bin/sh
parm ../003_leap/1HW9.wet.complex.parm7
#read in trajectory
trajin ../004_production/10.prod.trj
#read in reference
reference ../003_leap/1HW9.wet.complex.rst7
#compute RMSD + align CA to crystal structure
rmsd rms1 reference :1-263@CA
#strip solvent
strip :WAT:Na+:Cl-
#create gas-phase trajectory
trajout 1HW9_stripfit_restrained_gas.trj nobox

Your reference mask should include both the ligand and receptor. In order to execute: 

cpptraj -i cpptraj_strip_wat.in

Next, we'll measure the RMSD of the ligand, separate from the receptor. In order to do this, create the following input file: 

vi cpptraj_rmsd_lig.in

The input file should have the following parameters: 

#!/usr/bin/sh
#trajin the trajectory
trajin 1HW9_stripfit_restrained_gas.trj
#read in reference
reference ../003_leap/1HW9.gas.complex.rst7
#compute RMSD (do not fit internal geometris first, included rigid body motions, convert frames to ns (framenum*.005)
rmsd rms1 ":254&!(@H=)" nofit mass out 1HW9_lig_restrained_rmsd_nofit.dat time .005
#histogram the nofit rmsd
histogram rms1,*,*,.1,* norm out 1HW9_lig_restrained_rmsd_nofit_histogram.dat

Make sure that your reference mask only includes the ligand this time. In order to execute: 

cpptraj -p ../003_leap/1HW9.complex.parm7 -i cpptraj_rmsd_lig.in

Finally, we'll measure the RMSD of the receptor, separate from the ligand. In order to do this. create the following input file: 

vi cpptraj_rmsd_rec.in

The file should contain the following parameters: 

#!/usr/bin/sh
#trajin the trajectory
trajin 1HW9_stripfit_restrained_gas.trj
#read in reference
reference ../003_leap/1HW9.gas.complex.rst7
#compute RMSD (do not fit internal geometries first, included rigid body motions, convert frames to ns (framenum*.005)
rmsd rms1 ":1-253&!(@H=)" nofit mass out 1HW9_rec_restrained_rmsd_nofit.dat time .005
#histogram the nofit rmsd
histogram rms1,*,*,.1,* norm out 1HW9_rec_restrained_rmsd_nofit_histogram.dat

Make sure that your reference mask only includes the receptor this time. In order to execute: 

cpptraj -p ../003_leap/1HW9.gas.complex.parm7 -i cpptraj_rmsd_rec.in


Hydrogen Bonding

We can also determine if our ligand forms any hydrogen bonds with the receptor protein. In order to do this, we'll need to create the following input file: 

vi cpptraj_hbond.in

The input file should have the following parameters: 

!/usr/bin/sh
#read in trajectory
trajin ../004_production/10_prod.trj
#wrap everything in one periodic cell - could cause problems, may comment out #autoimage if problems later
autoimage
#compute intra + water-mediated H-bonds
hbond hb1 :1-263 out 1HW9_calcipotriol_hbond.out avgout 1HW9_calcipotriol_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
nointramol brid\
geout 1HW9_calcipotriol_bridge-water.dat dist 3.0 angle 140

The hbobnd hb1 line should contain the residue numbers corresponding to your protein and ligand. In order to execute the script: 

cpptraj -p ../003_leap/1HW9_wetcomplex.parm7 -i cpptraj_hbond.in


MMGBSA

In this section, we'll be estimating the free energy of binding of our ligand and receptor. This is performed by analyzing both the unsolvated and solvated ligand and receptor separately. In order to get started, we'll need the following input file: 

vi mmgbsa.in

Your input file should contain the following parameters: 

mmgbsa 1HW9 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,
 /

In order to run the input file, we'll need to create a separate script. 

vi mmgbsa.sh

This file should have the following parameters: 

#bin/bash
#SBATCH --time=2-00:00:00
#SBATCH --nodes=2
#SBATCH --ntasks=28
#SBATCH --job-name=1HW9_MMGBSA
#SBATCH --output=1HW9_MMGBSA.log
#SBATCH -p extended-28core
#Define topology files
solv_parm="../003_leap/1HW9.wet.complex.parm7"
complex_parm="../003_leap/1HW9.complex.parm7"
receptor_parm="../003_leap/1HW9.gas.receptor.parm7"
lig_parm="../003_leap/1HW9.gas.ligand.parm7"
trajectory="10_prod.trj"
 MMPBSA.py -O -i mmgbsa.in \
     -o 1S19_mmgbsa_results.dat \ 
     -eo 1S19_mmgbsa_perframe.dat \
     -sp ${solv_parm} \
     -cp ${complex_parm} \
     -rp ${receptor_parm} \
     -lp ${lig_parm} \
      -y ${trajectory}

You'll need to submit this to the queue: 

sbatch mmgbsa.sh
Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2021_AMBER_tutorial_1_with_PDBID_1HW9&oldid=14392"