2023 AMBER tutorial 2 with PDBID 3WZE

From Rizzo_Lab
Jump to: navigation, search

Introduction

AMBER is employed to compute biomolecular simulations between a receptor and a ligand. AMBER allows us to look at changes in interaction over time. We will be using the protein and receptor, 3WZE.

Directory Setup

As always, we set up folders to keep us organized as we move generate files:

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

3WZE Structures

We will generate fresh mol2 and pdb files of both the receptor and the ligand.

Receptor

To avoid any unintended errors, we start with a fresh PDB file and generate the receptor in Chimera. First, open the PDB in Chimera. Then select the protein

Select -> Structure -> protein

Then

Select -> Invert (all models)

to select everything but the protein. Next we delete everything but the protein by

Actions -> Atoms/Bonds -> delete

Note: we do NOT add hydrogens or add charge. This will be addressed in TLEaP below. And now we save the file as BOTH a mol2 and a PDB (3wze_rec.mol2; 3wze_rec.pdb). Sometimes mol2 files are tricky to load in AMBER so it is best to have both just in case.

Ligand

To generate only the ligand, open 3wze.pdb and repeat:

Select -> Structure -> protein

This time we want to delete the protein so we delete the selection by

Actions -> Atoms/Bonds -> delete

Next we want to remove the detergent, acetone, and waters by

Select -> Residue -> HOH -> Actions -> Atoms/Bonds -> delete
Select -> Residue -> DTT -> Actions -> Atoms/Bonds -> delete
Select -> Residue -> ACT -> Actions -> Atoms/Bonds -> delete

We will once again add hydrogens to the ligand by selecting the ligand and then

Tools > Structure editing > Add H

And we will add charge by

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

Now we save this as a mol2 file, 3wze_lig.mol2.

Go to the directory:

cd 001_structure

and transfer/save all files there.

Amber Simulation Parameters

Next, change directory to

cd 002_parameters

To generate parameters for the simulation, we must implement the following:

antechamber -i ../001_structure/3wze_lig_wH.mol2 -fi mol2 -o 3wze_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 0

nc = 0 because the charge on the ligand is 0. If your ligand is non-zero, enter the appropriate charge at the end of this line. It may be helpful to check the protonation state of the ligand at pH 7. Once 3wze_ligand_antechamber.mol2 output file is generated, run parmch2:

parmchk2 -i 3wze_ligand_antechamber.mol2 -f mol2 -o 3wze_ligand.am1bcc.frcmod

TLEaP

Next we will generate the AMBER topology file and coordinate files. Switch to the directory:

 cd 003_leap

Two types of files will be generated, parm7 (topology) and rst7 (coordinates). Create the input file:

vi leap.in

And then input the following (make sure you change the username):

   #!/USERNAME/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 ../000_structure/3wze_rec.pdb
   ##@make disulfide bond
   ###load ligand frcmod/mol2
   loadamberparams ../000_structure/3wze_ligand.am1bcc.frcmod  
   lig=loadmol2 ../000_structure/3wze_ligand_antechamber.mol2
   ###create gase-phase complex
   gascomplex= combine {rec lig}
   ###write gas-phase pdb
   savepdb gascomplex 3wze.gas.complex.pdb
   ###write gas-phase toplogy and coord files for MMGBSA calc
   saveamberparm gascomplex 3wze.complex.parm7 3wze.gas.complex.rst7
   saveamberparm rec 3wze.gas.receptor.parm7 3wze.gas.receptor.rst7
   saveamberparm lig 3wze.gas.ligand.parm7 3wze.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 3wze.wet.complex.pdb
   ###check the system
   charge solvcomplex 
   check solvcomplex
   ###write solvated toplogy and coordinate file
   saveamberparm solvcomplex 3wze.wet.complex.parm7 3wze.wet.complex.rst7
   quit

Once the files are generated, transfer the parm7 and rst7 files to the local environment. You can run them in Chimera to check the build. Open the protein, then open TOOLS--MD/ENSEMBLE ANALYSIS--MD MOVIE. Open the parm7 in prmtop box and then add the rst7 as a trajectory. Then click OK.

Equilibration

Next, we will minimize and equilibrate the system over nine successive steps. For information about the different parameters see the AMBER16 manual.

Change directory to:

cd 004_equil

The first step is to create the new file:

vi 01.min.md

And insert the following: 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
/

You'll want to pay close attention to the restraint mask - it lifts in later steps. Next:

vi 02.equil.mdin

Insert:

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
/

Then:

vi 03.min.mdin

and insert the following text:

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
/

Next:

vi 04.min.mdin

input:

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
/

Then:

vi 05.min.mdin

and insert the following text:

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
/
Next:
vi 06.equil.mdin

And:

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
/

Then:

vi 07.equil.mdin

Insert:

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
/

Next:

vi 08.equil.mdin

Here, you will need to specify the exact residues of your COMPLEX to enter at the restraintmask line. This means the total protein residues plus the ligand. Residue numbers which can be found in the wet.complex.pdb file generated in the 003_leap directory. So here, 3WZE is 298 residues plus the ligand is 299.

Insert the following:
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-299@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
/

Next:

vi 09.equil.mdin

Insert the following paying attention to the restraintmask line once more:

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

All nine files will be run simultaneously. First we create the file:

vi mdequilibration.sh

And insert the script:

#!/bin/bash
#
#SBATCH --job-name=3wze_equilibration
#SBATCH --output=equilibration_output.txt
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=1
#SBATCH --time=48:00:00
#SBATCH -p long-24core

module load amber/16

do_parallel="mpirun -np 80 pmemd.MPI"

prmtop="../003.tleap/3wze.wet.complex.prmtop"
coords="../003.tleap/3wze.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
  • Keep in mind that the 24 core is likely discontinued as of May 2023.

Submit to the queue:

sbatch mdequilibration.sh

Production

To generate the actual MD simulation, change to the following directory:

cd 005_production

And open the following file:

vi 10.prod.mdin

Note the restraint mask should include ONLY protein residues this time. Insert the following text:

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

This will take a very long time to run on CPUs so here it may be advisable to run on the GPU. Create the file:

vi mdproduction.sh

And insert the following text:

!/bin/sh
#SBATCH --job-name=3wze_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/3wze.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 to the queue.

sbatch mdproduction.sh

MD Analysis

Analysis and inspection of the RMSD, hydrogen bonds, and free energies of binding can now be evaluated.

RMSD

First, we create a file to remove the waters and undo changes from the production trajectory.

vi cpptraj_strip_wat.in

And insert the following text:

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

Reference mask once again contains protein residues + ligand. Now, run the following command:

cpptraj -i cpptraj_strip_wat.in 

Next, we create the file to generate the ligand RMSD:

vi cpptraj_rmsd_lig.in

And input the following text:

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

Run the command:

cpptraj -p ../003.tleap/3wze.complex.parm7 -i cpp_rmsd_lig.in

Next, we will calculate RMSD of the receptor. Create the file:

vi cpptraj_rmsd_rec.in

And input the following:

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

Restraint mask here only contains protein residues. Now run the command:

cpptraj -p ../003.tleap/3wze.complex.parm7 -i cpp_rmsd_receptor.in

Hydrogen Bonding

Next, we can calculate the number of hydrogen bonds. Create the file:

vi cpptraj_hbond.in

And insert the text:

!/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-299 out 3wze_hbond.out avgout 3wze_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
nointramol brid\
geout 3wze_bridge-water.dat dist 3.0 angle 140

Once again, include both protein and ligand residues (here, 299). Run the command:

cpptraj -p ../003.tleap/3wze.wet.complex.prmtop -i cpp_hbonds.in

MM-GBSA

And lastly, we calculate the free energy of binding. Create the file:

vi mmgbsa.in

And input the following text:

mmgbsa 3wze 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 long computation and so we create a script:

vi mmgbsa.sh

With the text:

#bin/bash
#SBATCH --time=2-00:00:00
#SBATCH --nodes=2
#SBATCH --ntasks=28
#SBATCH --job-name=3wze_MMGBSA
#SBATCH --output=3wze_MMGBSA.log
#SBATCH -p extended-28core
#Define topology files
solv_parm="../003.tleap/3wze.wet.complex.prmtop"
complex_parm="../003.tleap/3wze.complex.parm7"
receptor_parm="../003.tleap/3wze.gas.receptor.parm7"
lig_parm="../003.tleap/3wze.gas.ligand.parm7"
trajectory="../005.production/10_prod.trj"
MMPBSA.py -O -i mmgbsa.in \
     -o  3wze_mmgbsa_results.dat \ 
     -eo 3wze_mmgbsa_perframe.dat \
     -sp ${solv_parm} \
     -cp ${complex_parm} \
     -rp ${receptor_parm} \
     -lp ${lig_parm} \
      -y ${trajectory}

And execute with:

sbatch mmgbsa.sh