2022 AMBER tutorial 3 with PDBID 1X70
Contents
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
