Equilibrating Endstates in ABFE Calculations (Tyk2)
Learning objectives
Learn how to equilibrate the endstates in preparation for running absolute binding free energy calculations (ABFE).
Relevant literature
Coming Soon!
Tutorial
In this Tutorial, you will learn how to equilibrate the endstates in preparation for running ABFE calculations for TYK2.
Setting Up and Running End-State Equilibration
The starting point for this tutorial is the topology (parm7) and coordinate (rst7) files for the complex and aqueous endstates. The files for this tutorial are located:
Once you have downloaded the files, extract them into a directory where you will be working on this tutorial.
tar -xvzf equil_tyk2.tar.gz
cd equil_tyk2
mkdir working_dir
cd working_dir
The files you should see in the top level directory are parm7 and rst7 files for both the complex and aqueous end states. The parm7 files contain the hmr tag indicating that hydrogen mass repartitioning has been used to increase the timestep of the simulations to 4 fs.
Note
While this tutorial focuses on equilibrating for ABFE calculations, these are generally good practices for equilibrating any protein system using Amber.
Equilibrating the Aqueous Endstates
To equilibrate the aqueous endstates, we will be using a number of steps to gradually relax the system, the MD input files are included below, and again can be found in the tutorial files you downloaded above.
&cntrl
imin = 1, maxcyc = 1000,
ncyc = 100, ntx = 1, ntmin = 1,
ntwe = 0, ntwr = 0, ntpr = 1000,
ntc = 2, ntf = 2, ntb = 1, ntp = 0,
cut = 10.0, ntr = 1, restraintmask = ':1&!@H=', restraint_wt = 50,
ioutfm=1, ntxo=1,
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 0
ntx = 1
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 250
ntwe = 0
ntwr = 250
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 1.0
ntb = 1
ntp = 0
barostat = 2
ntr=1
restraint_wt= 50
restraintmask=':1&!@H='
nmropt=1
/
&wt TYPE="TEMP0"
istep1=0
istep2=1000000
value1=298
value2=298
/
&wt TYPE="END"
/
&cntrl
imin = 0
nstlim = 20000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 2.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 50
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 20000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 2.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 50
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 25
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 10
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 5
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 2
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 1
restraintmask=':1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.002
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 250
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 5.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 0.5
restraintmask='@CA,N,C | :1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.004,
irest = 1
ntx = 5
ig = -1,
ntc = 2
ntf = 2
tol = 1e-05,
nscm = 0
ioutfm=1
ntxo=1,
ntwx = 250
ntwe = 0
ntwr = 500
ntpr = 1000,
cut = 10.0
iwrap = 0,
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 5.0
ntb = 2
ntp = 1,
barostat = 2,
/
There are a total of 11 steps to equilibrate the aqueous endstate; however, many of these steps can be grouped together into general goals.
These goals are:
Minimize the system to remove bad contacts (step 1)
Heat the system to the target temperature while restraining the solute (step 2)
Equilibrate the density of the system while restraining the solute (steps 3 and 4)
Gradually relax the solute restraints (steps 5-10)
Run a short unrestrained simulation to relax the system (step 11)
Below, we have included the commands to run each of these steps using amber.
pmemd.cuda -O -i ../binder_mdins/minimization.mdin -c ../binder_ejm_31.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_minimization.nc -r binder_ejm_31_minimization.rst7 -o binder_ejm_31_minimization.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/heat.mdin -c binder_ejm_31_minimization.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_heat.nc -r binder_ejm_31_heat.rst7 -o binder_ejm_31_heat.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/density1.mdin -c binder_ejm_31_heat.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_density1.nc -r binder_ejm_31_density1.rst7 -o binder_ejm_31_density1.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/density2.mdin -c binder_ejm_31_density1.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_density2.nc -r binder_ejm_31_density2.rst7 -o binder_ejm_31_density2.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/lower_restraints1.mdin -c binder_ejm_31_density2.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_lower_restraints1.nc -r binder_ejm_31_lower_restraints1.rst7 -o binder_ejm_31_lower_restraints1.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/lower_restraints2.mdin -c binder_ejm_31_lower_restraints1.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_lower_restraints2.nc -r binder_ejm_31_lower_restraints2.rst7 -o binder_ejm_31_lower_restraints2.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/lower_restraints3.mdin -c binder_ejm_31_lower_restraints2.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_lower_restraints3.nc -r binder_ejm_31_lower_restraints3.rst7 -o binder_ejm_31_lower_restraints3.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/lower_restraints4.mdin -c binder_ejm_31_lower_restraints3.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_lower_restraints4.nc -r binder_ejm_31_lower_restraints4.rst7 -o binder_ejm_31_lower_restraints4.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/lower_restraints5.mdin -c binder_ejm_31_lower_restraints4.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_lower_restraints5.nc -r binder_ejm_31_lower_restraints5.rst7 -o binder_ejm_31_lower_restraints5.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/backbone.mdin -c binder_ejm_31_lower_restraints5.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_backbone.nc -r binder_ejm_31_backbone.rst7 -o binder_ejm_31_backbone.mdout -ref ../binder_ejm_31.rst7
pmemd.cuda -O -i ../binder_mdins/free.mdin -c binder_ejm_31_backbone.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_unrestrained.nc -r binder_ejm_31_unrestrained.rst7 -o binder_ejm_31_unrestrained.mdout -ref ../binder_ejm_31.rst7
Note
On Running Amber commands
The above commands use the pmemd.cuda executable to run the simulations on a GPU. If you are running on a CPU, you can replace pmemd.cuda with pmemd.MPI or sander as appropriate.
Every command uses a number of common flags, which we will break down below for the first pmemd.cuda command.
pmemd.cuda -O -i ../binder_mdins/minimization.mdin -c ../binder_ejm_31.rst7 -p ../binder_ejm_31_hmr.parm7 -x binder_ejm_31_minimization.nc -r binder_ejm_31_minimization.rst7 -o binder_ejm_31_minimization.mdout -ref ../binder_ejm_31.rst7
The flags used here are:
-O: Overwrite output files if they already exist
-i: Input mdin file
-c: Input coordinate file (rst7)
-p: Input topology file (parm7)
-x: Output trajectory file (nc)
-r: Output restart file (rst7)
-o: Output mdout file (text)
-ref: Reference coordinate file for restraints
In the mdin files, the restraints are defined using the reference coordinate file provided with the -ref flag. We are currently using Root Mean Square (RMS) restrains on the ligand heavy atoms to keep them close to this reference structure during equilibration.
Equilibrating the Complex Endstates
To equilibrate the complex endstates, we will follow essentially the same steps as we did for the aqueous endstate, with the same general goals. The primary difference between these mdins and the previous is the protein restraint mask 1-289, which restrains the protein residues in addition to the ligand. The mdin files are included below, and again can be found in the tutorial files you downloaded above.
Warning
Note that the mdin restrains here have been tailored to the TYK2 system, and should be adjusted if you are working with a different system. The key change to make is to adjust the protein residue range in the restraint mask (e.g. 1-289) to match the protein you are working with, plus 1 to account for the ligand being the first residue.
&cntrl
imin = 1, maxcyc = 1000,
ncyc = 100, ntx = 1, ntmin = 1,
ntwe = 0, ntwr = 0, ntpr = 1000,
ntc = 2, ntf = 2, ntb = 1, ntp = 0,
cut = 10.0, ntr = 1, restraintmask = ':1-289&!@H=', restraint_wt = 50,
ioutfm=1, ntxo=1,
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 0
ntx = 1
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 250
ntwe = 0
ntwr = 250
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 1.0
ntb = 1
ntp = 0
barostat = 2
ntr=1
restraint_wt= 50
restraintmask=':1-289&!@H='
nmropt=1
/
&wt TYPE="TEMP0"
istep1=0
istep2=1000000
value1=298
value2=298
/
&wt TYPE="END"
/
&cntrl
imin = 0
nstlim = 20000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 2.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 50
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 20000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 2.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 50
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 25
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 10
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 5
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 2
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.001
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 0
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 4.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 1
restraintmask=':1-289&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.002
irest = 1
ntx = 5
ig = -1
ntc = 2
ntf = 2
tol = 1e-05
nscm = 0
ioutfm=1
ntxo=1
ntwx = 250
ntwe = 0
ntwr = 500
ntpr = 1000
cut = 10.0
iwrap = 0
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 5.0
ntb = 2
ntp = 1
barostat = 2
ntr=1
restraint_wt= 0.5
restraintmask='@CA,N,C | :1&!@H='
/
&cntrl
imin = 0
nstlim = 1000000
dt = 0.004,
irest = 1
ntx = 5
ig = -1,
ntc = 2
ntf = 2
tol = 1e-05,
nscm = 0
ioutfm=1
ntxo=1,
ntwx = 250
ntwe = 0
ntwr = 500
ntpr = 1000,
cut = 10.0
iwrap = 0,
ntt = 3
temp0 = 298
tempi = 298 ! tempi is the initial temperature, temp0 is the target
gamma_ln = 5.0
ntb = 2
ntp = 1,
barostat = 2,
/
The commands to run each of these steps using amber are included below.
pmemd.cuda -O -i ../complex_mdins/minimization.mdin -c ../target_ejm_31.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_minimization.nc -r complex_ejm_31_minimization.rst7 -o complex_ejm_31_minimization.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/heat.mdin -c complex_ejm_31_minimization.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_heat.nc -r complex_ejm_31_heat.rst7 -o complex_ejm_31_heat.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/density1.mdin -c complex_ejm_31_heat.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_density1.nc -r complex_ejm_31_density1.rst7 -o complex_ejm_31_density1.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/density2.mdin -c complex_ejm_31_density1.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_density2.nc -r complex_ejm_31_density2.rst7 -o complex_ejm_31_density2.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/lower_restraints1.mdin -c complex_ejm_31_density2.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_lower_restraints1.nc -r complex_ejm_31_lower_restraints1.rst7 -o complex_ejm_31_lower_restraints1.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/lower_restraints2.mdin -c complex_ejm_31_lower_restraints1.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_lower_restraints2.nc -r complex_ejm_31_lower_restraints2.rst7 -o complex_ejm_31_lower_restraints2.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/lower_restraints3.mdin -c complex_ejm_31_lower_restraints2.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_lower_restraints3.nc -r complex_ejm_31_lower_restraints3.rst7 -o complex_ejm_31_lower_restraints3.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/lower_restraints4.mdin -c complex_ejm_31_lower_restraints3.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_lower_restraints4.nc -r complex_ejm_31_lower_restraints4.rst7 -o complex_ejm_31_lower_restraints4.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/lower_restraints5.mdin -c complex_ejm_31_lower_restraints4.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_lower_restraints5.nc -r complex_ejm_31_lower_restraints5.rst7 -o complex_ejm_31_lower_restraints5.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/backbone.mdin -c complex_ejm_31_lower_restraints5.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_backbone.nc -r complex_ejm_31_backbone.rst7 -o complex_ejm_31_backbone.mdout -ref ../target_ejm_31.rst7
pmemd.cuda -O -i ../complex_mdins/free.mdin -c complex_ejm_31_backbone.rst7 -p ../target_ejm_31_hmr.parm7 -x complex_ejm_31_unrestrained.nc -r complex_ejm_31_unrestrained.rst7 -o complex_ejm_31_unrestrained.mdout -ref ../target_ejm_31.rst7
Again, run these, and you will be left with equilibrated endstate files for both the complex and aqueous endstates, which can be used to start the production ABFE calculations. You can repeat these steps for each of the remaining three ligands; the mdin files are included in the tutorial files you downloaded above.