Difference between revisions of "2022 AMBER tutorial 3 with PDBID 1X70"

Latest revision as of 21:04, 4 May 2022

Introduction

In this tutorial, we will be doing molecular dynamics (MD) simulations using Amber 16. This will allow us to obtain much more detailed information about the binding affinities of a ligand with a receptor. Like the DOCK tutorial, we should start with setting up the directory. This can be done with the following command.

 mkdir 001_structure 002_parameters 003_leap 004_equil 005_production

Making the structures

It is recommended to start fresh with the receptor and ligand files. Do the same process as in the DOCK tutorial. Start with importing pdb code 1X70, then isolate one of the proteins and its associated ligand. With an isolated protein, you can do the following command

 Select -> residue -> all nonstandard

then do

 actions -> atoms/bonds -> delete

This should isolate the protein and remove all of the nonstandard amino acids. With this, we can save this as a pdb and title it 1X70_fresh.pdb

Moving on to the ligand, we can repeat the same process as above but just isolate the ligand instead. So start with the fresh pdb code, then select the ligand and do

 Select -> Invert (selected molecules)

and

 actions -> atoms/bonds -> delete

This should isolate the ligand. Then we have to make sure it is charged properly and has hydrogens. We can do that with the following two commands.

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

Once this is done, save the file as 1x70_lig_wH.mol2

After doing this with both files, upload them to seawulf and place them in the 001_structure directory.

Simulation parameters

In this section, we are generating the force field parameters that will be used in future steps. The first step is to swap to the parameters directory, 002_parameters, then use antechamber.

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

The -nc tag at the end is the total charge of the ligand, so this may be different depending on the system. In this case, the net charge of the ligand is 1. After this step is completed, we need to parmch2 to generate another file.

parmchk2 -i 1x70_ligand_antechamber.mol2 -f mol2 -o 1x70_ligand.am1bcc.frcmod

TLeap

The next step is to create the files needed to simulate the system. Copy the following code into a file named leap.in

 #!/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/1x70_fresh.pdb
 #bond rec.649.SG rec.762.SG
 #bond rec.454.SG rec.472.SG
 #bond rec.444.SG rec.447.SG
 #bond rec.328.SG rec.339.SG
 #bond rec.385.SG rec.394.SG
 ###load ligand frcmod/mol2
 loadamberparams ../002_parameters/1x70_ligand.am1bcc.frcmod
 lig=loadmol2 ../002_parameters/1x70_ligand_antechamber.mol2

 ###create gase-phase complex
 gascomplex= combine {rec lig}

 ###write gas-phase pdb
 savepdb gascomplex 1x70.gas.complex.pdb
 ###write gase-phase toplogy and coord files for MMGBSA calc
 saveamberparm gascomplex 1x70.complex.parm7 1x70.gas.complex.rst7
 saveamberparm rec 1x70.gas.receptor.parm7 1x70.gas.receptor.rst7
 saveamberparm lig 1x70.gas.ligand.parm7 1x70.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 1x70.wet.complex.pdb

 ###check the system
 charge solvcomplex
 check solvcomplex

 ###write solvated toplogy and coordinate file
 saveamberparm solvcomplex 1x70.wet.complex.parm7 1x70.wet.complex.rst7
 quit

Note the commented out bond section in this code. This is needed for any disulfide bonds that are found in the receptor. For this system, the disulfide bonds were already made in the .pdb file. They have a tag of SSBOND in that file. In the event that your .pdb file does not contain the disulfide bond information, you will have to manually add those bonds in using the same format as above. With the input file copied, we can use the command

tleap -f leap.in

to create the system. Once this is done, there should be 10 parameter/coordinate files in your directory. These can be loaded into chimera with the following command

Tools -> MD/Ensemble Analysis -> MD Movie --> Select the parameter and coordinate file

It is important to check them in Chimera just to make sure that everything is running properly.

Equilibration

Before we can move into the primary production MD simulation, we have to equlibrate the system. Our equilibration process is 9 steps. Copy the following files into your 004_equil directory.

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-729@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-728@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
/

The restraint mask values for steps 8 and 9 are system-specific. The one for step 8 should include both the receptor and the ligand, and the one for step 9 should include only the receptor. The values for this can be found by looking at the wet.complex.pdb file that was made earlier. Search for the residue labelled LIG, and that will be the number for step 8 (in our case it was 729), and for step 9, it will be 1 number less (removes ligand from the mask). To run the equilibration process, create the following input file titled mdequilibration.sh

#!/bin/sh
#SBATCH --job-name=1x70_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/1x70.wet.complex.parm7"
coords="../003_leap/1x70.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' "

And submit it to the queue with the command

sbatch mdequilibration.sh

Production MD

With all of the set up and equilibration steps done, we can move onto the production MD for our system. Switch over to the 005_production directory and create the file 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-728@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 that the restraint mask should be the same as in step 9 of the equilibration. We can now create the slurm file mdproduction.sh with the following script

!/bin/sh
#SBATCH --job-name=1x70_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/1x70.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` "

Submit this to the queue with

sbatch mdproduction.sh

This will take a while to run, so try not to wait until the last minute to get this job in.

More Analysis

With the production MD finished, we can move on to some additional analysis. This will include RMSD, hydrogen bonding, and MMGBSA.

RMSD

The analysis will start with looking at the RMSD. This looks at how much the ligand and receptor moved during the MD simulation. To start, create a file named cpptraj_strip_wat.in to remove the water from the simulation. Paste the following code in

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

Note that like the previous steps, the reference mask is dependent on the specific system. It should include both the receptor and the ligand. Use the following command to execute the code.

cpptraj -i cpptraj_strip_wat.in

Now we will look at the RMSD for the ligand. Create a file named cpptraj_rmsd_lig.in and paste the following code

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

In this case, the mask should only include the ligand. Use the following command to execute the code.

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

Now we will do the same process but with the receptor instead. Make a file titled cpptraj_rmsd_rec.in and past in the following code.

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

This time, the mask should only include the receptor. And now use the following command to run this code.

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

Hydrogen bonding

Next, we will see if the ligand is forming any hydrogen bonds with the residues on the receptor. To start, make a file titled cpptraj_hbond.in and paste in the following code.

!/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-729 out 1HW9_calcipotriol_hbond.out avgout 1x70_calcipotriol_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
nointramol brid\
geout 1x70_calcipotriol_bridge-water.dat dist 3.0 angle 140

The line starting with hbond should contain the residue numbers for both the receptor and the ligand. Run this code with the following command.

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

MMGBSA

Lastly, we will be looking at the MMGBSA, which measures the binding affinity of the ligand and the receptor. To start, make a file titled mmgbsa.in and put in the following code.

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

This is a more involved script to run, so we will need to submit a job to run it on the cluster. Create a shell script titled mmgbsa.sh and put in the following commands.

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

Then just submit it to the queue with the following command

sbatch mmgbsa.sh