Generating and Running Lambda Windows for ABFE Calculations (Tyk2)

Zeke A. Piskulich¹, Patricio Barletta¹, Ryan Snyder¹, and Darrin M. York¹

¹Laboratory for Biomolecular Simulation Research, Institute

for Quantitative Biomedicine and Department of Chemistry and Chemical

Biology, Rutgers University, Piscataway, NJ 08854, USA

Learning objectives

Learn how to set up and generate lambda windows for absolute binding free energy (ABFE) calculations starting from an equilbrated end-state.
Learn how to run ABFE calculations using Amber

Relevant literature

Coming Soon!

Tutorial

In this Tutorial, you will learn how to set up and run absolute binding free energy calculations using Amberflow, a workflow engine developed in the York Lab for automating Alchemical Free Energy Calculations with Amber.

Setting Up and Running Annihilation Simulations 

For this tutorial, we will need two things to get started. The first is the equilibrated end-state files from the previous tutorial:ref:Equilibrating the End States of Tyk2 for ABFE Calculations <equilibrate_endstates>; and the second is the Boresch restraint files generated in the tutorial:ref:Generating Boresch Restraints for Tyk2 ABFE Calculations <boresch_restraints>. For your convenience, we have provided the necessary files below to get you started.

Attention

Update the download link to point to the finalized tutorial files once they are uploaded.

Choosing A Lambda Schedule 

The first step in setting up an ABFE calculation is to define a lambda schedule. A later tutorial will cover this topic in more detail with information on how to optimize lambda schedules for your specific system; however, for the purposes of this tutorial we will use a linear lambda schdule with 11 windows.

Warning

A linear lambda schedule is usually not optimizal for ABFE calculations, and ABFE calculations often need far more than 11 windows to converge. This tutorial is meant to demonstrate the basic steps of setting up and running an ABFE calculation, and is not intended to produce converged or realistic results.

The lambda schedule we will use is as follows:

Window  Lambda Value
---------------------
      0.0
      0.1
      0.2
      0.3
      0.4
      0.5
      0.6
      0.7
      0.8
      0.9
      1.0

Annihilation Simulation Setup 

With a lambda schedule in place, we can now slowly annihilate the ligand from the equilibrated complex and binder end-states. Roughly, in both cases this involves starting from the end-state where \(\lambda = 0\) and gradually turning off the ligand’s interactions with its environment until it is fully decoupled at \(\lambda = 1\).

Annhilate (Binder)

&cntrl
imin = 1, 
maxcyc = 5000,
ncyc = 100, 
ntx = 1, 
ntmin = 2,
ntwe = 0, 
ntwr = 0, 
ntpr = 1000,
ntb = 1, 
ntp = 0,
cut = 10.0, 
ntr = 0,
ioutfm=1, 
ntxo=1,

ifsc            = 1
icfe            = 1

clambda         = REPLACE_WITH_LAMBDA
timask1         = ':1'
timask2         = ''
scmask1         = ':1'
scmask2         = ''
scalpha         = 0.5
scbeta          = 1.0

gti_cut         = 1
gti_output      = 1
gti_add_sc      = 25
gti_scale_beta  = 1
gti_cut_sc_on   = 8
gti_cut_sc_off  = 10
gti_lam_sch     = 1
gti_ele_sc      = 1
gti_vdw_sc      = 1
gti_cut_sc      = 2
gti_ele_exp     = 2
gti_vdw_exp     = 2
gti_syn_mass    = 0
gremd_acyc      = 1

gti_bat_sc      = -1
/

Annhilate (Complex)

&cntrl
imin = 1, 
maxcyc = 5000,
ncyc = 100, 
ntx = 1, 
ntmin = 2,
ntwe = 0,
ntwr = 0,
ntpr = 1000,
ntb = 1, 
ntp = 0,
cut = 10.0, 
ntr = 0,
ioutfm=1, 
ntxo=1,

ifsc            = 1
icfe            = 1

clambda         = REPLACE_WITH_LAMBDA
timask1         = ':1'
timask2         = ''
scmask1         = ':1'
scmask2         = ''
scalpha         = 0.5
scbeta          = 1.0

gti_cut         = 1
gti_output      = 1
gti_add_sc      = 25
gti_scale_beta  = 1
gti_cut_sc_on   = 8
gti_cut_sc_off  = 10
gti_lam_sch     = 1
gti_ele_sc      = 1
gti_vdw_sc      = 1
gti_cut_sc      = 2
gti_ele_exp     = 2
gti_vdw_exp     = 2
gti_syn_mass    = 0
gremd_acyc      = 1

gti_bat_sc      = -1
nmropt          = 1
/
&wt type = 'END', /
DISANG=rest.in
/

Above, we have included the mdin files used to annihilate the ligand in both the binder and the complex environments. Note. These mdin files contain a placeholder for the specific lambda value for each window. This is REPLACE_WITH_LAMBDA.

We will start by generating the lambda windows for both the binder and complex annihilation simulations. To do this, we will start with our simple lambda schedule defined above.

Start by creating a working_directory for these steps within the abfe_files directory. We will make separate directories for the binder and complex annihilation simulations within this folder.

mkdir working_dir
cd working_dir
mkdir binder complex

Next, we will use a bash for loop to generate the mdin files for each lambda window for both the binder and complex annihilation simulations.

# Binder Annhilation Windows
lambda_values=(0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0)
for i in ${!lambda_values[@]}
do
   lambda_value=${lambda_values[$i]}
   sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../binder_mdins/annhilate.mdin > binder/annihilate_lambda_${i}.mdin
done

# Complex Annhilation Windows
for i in ${!lambda_values[@]}
do
   lambda_value=${lambda_values[$i]}
   sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../complex_mdins/annhilate.mdin > complex/annihilate_lambda_${i}.mdin
done

Now, you should have 11 mdin files for both the binder and complex annihilation simulations, each corresponding to a different lambda window.

Now, we can proceed to run these lambda windows iteratively starting from the \(\lambda = 0\) window and moving to the next window using the output of the previous window as input.

We will start with the binder simulations. Here, we will use the equilibrated end-state files (provided with the tutorial files) as input for the first window (\(\lambda = 0\)), and then use the output of each window as input for the next window.

cd binder
# Run Lambda Windows Iteratively
for i in {0..10}
do
   if [ $i -eq 0 ]
   then
      # First window uses equilibrated end-state files
      pmemd.cuda -O -i annihilate_lambda_${i}.mdin -p  ../../topology/binder_ejm_31_hmr.parm7 -c ../../topology/binder_ejm_31_unrestrained.rst7 -o annihilate_binder_${i}.mdout -r annihilate_binder_${i}.rst7 -x annihilate_binder_${i}_traj.nc -ref ../../topology/binder_ejm_31_unrestrained.rst7
   else
      # Subsequent windows use output of previous window
      pmemd.cuda -O -i annihilate_lambda_${i}.mdin -p ../../topology/binder_ejm_31_hmr.parm7 -c annihilate_binder_$((i-1)).rst7 -o annihilate_binder_${i}.mdout -r annihilate_binder_${i}.rst7 -x annihilate_binder_${i}_traj.nc -ref annihilate_binder_$((i-1)).rst7
   fi
done

These commands will run the binder annilation simulations for each lambda window iteratively. These are minimizations, and thus are relatively quick to run.

For the complex simulations, we will follow the same procedure, except we also need to include the Boresch restraint rest.in file, as well as the lambda.sch file that defines how the restraints are scaled during the annihilation process.

This lambda.sch file controls how the Boresch restraints are scaled during the annihilation process to ensure that the ligand remains properly restrained as its interactions with the environment are turned off.

TypeRestBA, smooth_step2, symmetric, 1.0, 0.0

Warning

The lambda.sch is essential for the proper scaling of Boresch restraints during the annihilation process. If not included, the restraints will be applied to the real state and then slowly turned off, which is exactly the opposite of what we want to achieve. This is a very common mistake when setting up ABFE calculations with Boresch restraints, and will lead to incorrect results.

Thus, we will copy these files into the complex working directory before running the simulations.

cd ../complex
cp ../../boresch_restraints/rest.in .
cp ../../boresch_restraints/lambda.sch .
# Run Lambda Windows Iteratively
for i in {0..10}
do
   if [ $i -eq 0 ]
   then
      # First window uses equilibrated end-state files
      pmemd.cuda -O -i annihilate_lambda_${i}.mdin -p ../../topology/target_ejm_31_hmr.parm7 -c ../../topology/complex_ejm_31_unrestrained.rst7 -o annihilate_complex_${i}_out.log -r annihilate_complex_${i}.rst7 -x annihilate_complex_${i}.nc -ref ../../topology/complex_ejm_31_unrestrained.rst7
   else
      # Subsequent windows use output of previous window
      pmemd.cuda -O -i annihilate_lambda_${i}.mdin -p ../../topology/target_ejm_31_hmr.parm7 -c annihilate_complex_$((i-1)).rst7 -o annihilate_complex_${i}.mdout -r annihilate_complex_${i}.rst7 -x annihilate_complex_${i}.nc -ref annihilate_complex_$((i-1)).rst7
   fi
done

These calculations will run relatively quickly, as they are minimizations. After this, we can move on to equilibrating the lambda windows.

Generating and Running Lambda Windows for ABFE Calculations (Tyk2)

Learning objectives

Relevant literature

Tutorial

Setting Up and Running Annihilation Simulations 

Choosing A Lambda Schedule 

Annihilation Simulation Setup 

Equilibrating the Lambda Windows 

Binder Equilibration 

Complex Equilibration 

Running Production Simulations 

Binder Production 

Complex Production 