2016 AMBER tutorial with Beta Trypsin
For additional Rizzo Lab tutorials see AMBER Tutorials.
In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).
Contents
I. Introduction
Yaping
II. Structural Preparation
Antechamber, Parmchk, tLeap
Omar, Katie
Antechamber, Parmchk, tLeap
Before beginning the Molecular Dynamics protocol using AMBER, you must first set up your files. In your 001.chimera folder, you will add the following files:
1BJU.lig.mol2 1BJU.rec.noH.pdb 1BJU.rec.noH.ZIN.pdb
To prepare the first two files, please see the 2025 DOCK tutorial at the following link: [1]
Note: delete any headers before the atoms/helix information.
In order to prepare 4TKG.rec.noH.ZIN.pdb, first open 4TKG.rec.noH.pdb and use the command "SHIFT + g" to reach the bottom of the pdb. The last lines of your files should look like this:
ATOM 1648 OE1 GLU B 205 23.839 -23.190 57.747 1.00 0.00 O TER 1649 GLU B 205 HETATM 1650 ZN ZN B 206 28.130 -3.467 55.482 1.00 0.00 Zn END
To make 4TKG.rec.noH.ZIN.pdb, you will need to change the "ZN" atom ID to "ZIN" so that AMBER can read the atom type.
ATOM 1648 OE1 GLU B 205 23.839 -23.190 57.747 1.00 0.00 O TER 1649 GLU B 205 ATOM 1650 ZIN ZIN B 206 28.130 -3.467 55.482 1.00 0.00 Zn END
To run a tleap, first generate a 002.tleap file, with the command:
mkdir 002.tleap
cp ~/../../tleap.lig.in
- tleap.lig.in
source leaprc.ff14SB #Load a force field source leaprc.gaff LIG = loadmol2 4TKG.lig.mol2 #Conformation file saveamberparm LIG 4TKG.lig.gas.leap.prm7 4TKG.lig.gas.leap.rst7 #Save the ligand gas phase AMBER topology and coordinate file solvatebox LIG TIP3PBOX 10.0 #Solvate the lig using TIP3P, solvent box radii 10 angstroms saveamberparm LIG 4TKG.lig.wat.leap.prm7 4TKG.lig.wat.leap.rst7 #Save the ligand water phase AMBER topology and coordinate files quit
tleap –s –f tleap.lig.in > tleap.lig.out
The atom types which used in '4TKG.lig.mol2' file are not recognized.
So the name of the atoms has to be redefined based on the gaff force field, antechamber can be used to solve the error.
antechamber –i 4TKG.lig.mol2 –fi mol2 –o 4TKG.lig.gaff.mol2 –fo mol2
Change the content of tleap.lig.in to (gaff file)
Re-run tleap by:
tleap –s –f tleap.lig.in > tleap.lig.out
So now run parmchk to fix the missing parameters
parmchk –i 4TKG.lig.gaff.mol2 –f mol2 –o 4TKG.lig.ante.frcmod
Change the content of tleap.lig.in to (load amber params) run tleap.lig.in, and all the errors are fixed.
cp ~/../../tleap.rec.in
- tleap.rec.in
source leaprc.ff14SB #Load a force field REC = loadpdb ../001.chimera/4TKG.rec.noH.ZIN.pdb #Conformation file for receptor saveamberparm REC 4TKG.rec.gas.leap.prm7 4TKG.rec.gas.leap.rst7 #Save the receptor gas phase AMBER topology and coordinate file solvateBox REC TIP3PBOX 10.0 #Solvate the lig using TIP3P, solvent box radii 10 angstroms saveamberparm REC 4TKG.rec.wat.leap.prm7 4TKG.rec.wat.leap.rst7 #Save the receptor water phase AMBER topology and coordinate file quit
tleap –s –f tleap.rec.in > tleap.rec.out
For the missing parameters, cp the two ionparam from ../000.amberfiles/
cp ../000.amberfiles/ions.lib cp ../000.amberfiles/ions.frcmod
Change the content of the tleap.rec.in to (load param)
run tleap.rec.in, and all the errors are fixed. (explain)
cp ~/../../tleap.com.in tleap –s –f tleap.com.in > tleap.com.out
This whole process can be run on a cluster by using the following script:
III. Simulation using pmemd
PMEMD
Agatha, Beilei
MD input file1 01mi.in:
mkdir 01mi.in
01mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 5.0,
Copy 01md.in and rename it to 02md.in:
cp 01mi.in 02md.in
02md.in: equilibration &cntrl imin = 0, ntx = 1, irest = 0, nstlim = 50000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.001, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-224 & !@H=', restraint_wt = 5.0,
Using the same way to make 03mi.in, 04mi.in, 05mi.in, 06md.in, 07md.in, 08md.in, 09md.in, 10md.in, 11md.in:
03mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 2.0,
04mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 0.1,
05mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 0.05,
06md.in: equilibration &cntrl imin = 0, ntx = 1, irest = 0, nstlim = 50000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.001, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-224 & !@H=', restraint_wt = 1.0,
07md.in: equilibration &cntrl imin = 0, ntx = 5, irest = 1, nstlim = 50000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.001, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-224 & !@H=', restraint_wt = 0.5,
08md.in: equilibration &cntrl imin = 0, ntx = 5, irest = 1, nstlim = 50000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.001, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,
09md.in: equilibration &cntrl imin = 0, ntx = 5, irest = 1, nstlim = 50000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.001, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,
10md.in: production &cntrl imin = 0, ntx = 5, irest = 1, nstlim = 500000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.002, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,
11md.in: production &cntrl imin = 0, ntx = 5, irest = 1, nstlim = 500000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.002, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,
IV. Simulation Analysis
Ptraj
Lauren, Haoyue
ptraj.strip.wat.in:
trajin ../003.PMEMD/10md.trj 1 1000 1 trajin ../003.PMEMD/11md.trj 1 1000 1 trajout 1BJU.trj.strip nobox strip :WAT
ptraj.rec.rmsd.in:
trajin 1BJU.trj.strip trajout 1BJU.com.trj.stripfit reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7 rms reference out 1BJU.rmsd.CA.dat :1-223@CA
ptraj.lig.rmsd.in:
trajin 1BJU.com.trj.stripfit reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7 rms reference out 1BJU.lig.rmsd.dat :224@C*,N*,O*,S* nofit
ptraj.strip.rec.in:
trajin 1BJU.com.trj.stripfit 1 2000 1 trajout 1BJU.lig.trj.stripfit strip :1-223
ptraj.strip.lig.in:
trajin 1BJU.com.trj.stripfit 1 2000 1 trajout 1BJU.rec.trj.stripfit strip :224
ptraj.dist.in:
trajin 1BJU.com.trj.stripfit distance 171OD2_224N1 :171@OD2 :224@N1 out 171OD2_224N1.dat distance 171OD1_224N2 :171@OD1 :224@N2 out 171OD1_224N2.dat
MM-GBSA Energy Calculation
begin by creating a new directory
mkdir 005.MMGBSA/
as well as an input file
vim.gbrescore.in
the following parameters will provide a good estimate of the energy of binding based on solvation enthalpy and binding enthalpy
Single point GB energy calc &cntrl ntf = 1, ntb = 0, ntc = 2, idecomp= 0, igb = 5, saltcon= 0.00, gbsa = 2, surften= 1.0, offset = 0.09, extdiel= 78.5, cut = 99999.0, nsnb = 99999, imin = 5, maxcyc = 1, ncyc = 0, /
We can now submit a job to the cluster we are working on. This script, inside a file named run.sander.rescore.csh , is designed for LIRED, however it can be modified for whatever supercomputer you are working on
#! /bin/tcsh #PBS -l nodes=1:ppn=1 #PBS -l walltime=48:00:00 #PBS -o zzz.qsub.out #PBS -e zzz.qsub.err #PBS -V #PBS -N mmgbsa set workdir = /nfs/guest46/amber_tutorial/005.mmgbsa cd $workdir sander -O -i gb.rescore.in \ -o gb.rescore.out.com \ -p ../002.tleap/4TKG.com.gas.leap.prm7 \ -c ../002.tleap/4TKG.com.gas.leap.rst7 \ -y ../004.ptraj/4TKG.com.trj.stripfit \ -r restrt.com \ -ref ../002.tleap/4TKG.com.gas.leap.rst7 \ -x mdcrd.com \ -inf mdinfo.com
sander -O -i gb.rescore.in \ -o gb.rescore.out.lig \ -p ../002.tleap/4TKG.lig.gas.leap.prm7 \ -c ../002.tleap/4TKG.lig.gas.leap.rst7 \ -y ../004.ptraj/4TKG.lig.trj.stripfit \ -r restrt.lig \ -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \ -x mdcrd.lig \ -inf mdinfo.lig sander -O -i gb.rescore.in \ -o gb.rescore.out.test.rec \ -p ../002.tleap/4TKG.rec.gas.leap.prm7 \ -c ../002.tleap/4TKG.rec.gas.leap.rst7 \ -y ../004.ptraj/4TKG.rec.trj.stripfit \ -r restrt.rec \ -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \ -x mdcrd.rec \ -inf mdinfo.rec exit
Then submit the job as a .csh
qsub run.sander.rescore.csh
Your ligand, receptor, and receptor-ligand complex will all have .out files. The format should look
FINAL RESULTS
NSTEP ENERGY RMS GMAX NAME NUMBER
1 5.9132E+03 2.0005E+01 1.2640E+02 C 159
BOND = 661.8980 ANGLE = 1751.7992 DIHED = 2581.7692
VDWAALS = -1696.6585 EEL = -13958.9335 EGB = -3125.9524
1-4 VDW = 747.0185 1-4 EEL = 7750.8118 RESTRAINT = 0.0000
ESURF = 11201.4791
minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02
We can now use a system of equations to elucidate various energies related to our system. VDWAALS = ΔGvdw E(ELS) = ΔGcoulombic E(GB) = ΔGpolar SASA = ESURF
ΔGnonpolar can be estimated by:
ΔGnonpolar = (Solvent Accessible Surface Area)*0.00542 + 0.92
This equation is based on empirical observations relating a molecules surface area and solvation energy
ΔGnonpolar and ΔGmmgbsa are related by anther equation:
ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar
all the information needed is in our output files.
ΔΔGbind can now be solved via the equation:
ΔΔGbind = ΔG(mmgbsa complex) – (ΔG(mmgbsa ligand) + ΔG(mmgbsa receptor) (3)
ΔΔGbind is a function of ligand position, so it is useful to plot it. For publications, it is standard protocol to include the mean and standard deviation for your ΔΔGbind.
Lets make a script named get.mmgbsa.csh to extract this sophisticated energetic estimation from the three initial output files of the calculations obtained solve for ΔΔGbind:
bash get.mmgbsa.csh
The output file will be .dat and we can easily use the program xmgrace to plot the information visually
xmgrace MMGBSA_vs_time.dat
