Generating and Running Lambda Windows for ABFE Calculations (Tyk2) ================================================================== | Zeke A. Piskulich\ :sup:`1`, Patricio Barletta\ :sup:`1`, Ryan Snyder\ :sup:`1`, and Darrin M. York\ :sup:`1` | :sup:`1`\ Laboratory for Biomolecular Simulation Research, Institute | for Quantitative Biomedicine and Department of Chemistry and Chemical | Biology, Rutgers University, Piscataway, NJ 08854, USA Learning objectives ------------------- .. start-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 .. end-learning-objectives 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. .. contents:: :local: :depth: 4 .. start-tutorial 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 `; and the second is the Boresch restraint files generated in the tutorial:ref:`Generating Boresch Restraints for Tyk2 ABFE Calculations `. 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: .. code-block:: text Window Lambda Value --------------------- 0 0.0 1 0.1 2 0.2 3 0.3 4 0.4 5 0.5 6 0.6 7 0.7 8 0.8 9 0.9 10 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 :math:`\lambda = 0` and gradually turning off the ligand's interactions with its environment until it is fully decoupled at :math:`\lambda = 1`. .. tab-set:: .. tab-item:: Annhilate (Binder) .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/binder_mdins/annhilate.mdin .. tab-item:: Annhilate (Complex) .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/complex_mdins/annhilate.mdin 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. .. code-block:: bash 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. .. code-block:: bash # 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 :math:`\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 (:math:`\lambda = 0`), and then use the output of each window as input for the next window. .. code-block:: bash 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. .. code-block:: text 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. .. code-block:: bash 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. Equilibrating the Lambda Windows ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We will now equilibrate each of the lambda windows generated in the previous step. Unlike the annihilation simulations, we will run the equilibration simulations for each lambda window in parallel, as they are independent of each other. Currently, these parallel runs will not be using replica exchange - replica exchange won't turn on until the production phase. Binder Equilibration ^^^^^^^^^^^^^^^^^^^^ We will again start with the binder, equilibrating each lambda window using the corresponding mdin files provided below. . .. tab-set:: .. tab-item:: Heat .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/binder_mdins/heat.mdin .. tab-item:: Equil1 .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/binder_mdins/equil1.mdin .. tab-item:: Equil2 .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/binder_mdins/equil2.mdin .. tip:: These equilibration mdin files contain a placeholder for the specific lambda value for each window. This is REPLACE_WITH_LAMBDA. .. warning:: Remember that mdin files REQUIRE a blank line at the end of the file to function properly with Amber. If you get an error that there was an unexpected end of file while reading the mdin, check to make sure there is a blank line at the end of your mdin file. If you use the provided mdin files, they already have this blank line. Just like before, we will generate the mdins programmatically for each lambda window using a bash for loop. .. code-block:: bash cd ../binder # Generate Equilibration MDINs for Each Lambda Window 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 {0..10} do lambda_value=${lambda_values[$i]} sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../binder_mdins/heat.mdin > heat_lambda_${i}.mdin sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../binder_mdins/equil1.mdin > equil1_lambda_${i}.mdin sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../binder_mdins/equil2.mdin > equil2_lambda_${i}.mdin done Now, we also need a groupfile for the replica exchange simulations. This file will look a lot like the pmemd.CUDA command calls that we made before, except the groupfile will contain all of the lambda windows to be run in parallel. The script below will generate the groupfile for the binder equilibration simulations. .. code-block:: bash rm binder_heat_groupfile.in binder_equil1_groupfile.in binder_equil2_groupfile.in for i in {0..10} do echo "-O -i heat_${i}.mdin -p ../../topology/binder_ejm_31_hmr.parm7 -c annihilate_binder_${i}.rst7 -o binder_heat_${i}.mdout -r binder_heat_${i}.rst7 -x binder_heat_${i}.nc -ref ../../topology/binder_ejm_31_unrestrained.rst7" >> binder_heat_groupfile.in echo "-O -i equil1_${i}.mdin -p ../../topology/binder_ejm_31_hmr.parm7 -c binder_heat_${i}.rst7 -o binder_equil1_${i}.mdout -r binder_equil1_${i}.rst7 -x binder_equil1_${i}.nc -ref ../../topology/binder_ejm_31_unrestrained.rst7" >> binder_equil1_groupfile.in echo "-O -i equil2_${i}.mdin -p ../../topology/binder_ejm_31_hmr.parm7 -c binder_equil1_${i}.rst7 -o binder_equil2_${i}.mdout -r binder_equil2_${i}.rst7 -x binder_equil2_${i}.nc -ref ../../topology/binder_ejm_31_unrestrained.rst7" >> binder_equil2_groupfile.in done Now we can run these calculations in parallel using pmemd.cuda.MPI. .. code-block:: bash mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile binder_heat_groupfile.in mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile binder_equil1_groupfile.in mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile binder_equil2_groupfile.in .. hint:: Here, note that the usual pmemd commands are included in the group file. On the actual command line, our calls are much simpler now. Take this for example: .. code-block:: bash mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile binder_heat_groupfile.in - mpirun -np 11: This tells mpirun to use 11 processes (one for each lambda window). - pmemd.cuda.MPI: This is the MPI-enabled version of pmemd.cuda, which allows for parallel execution across multiple processes. - -ng 11: This flag specifies the number of groups (or simulations) to run in parallel, which is 11 in this case. - -groupfile binder_heat_groupfile.in: This specifies the group file that contains the individual pmemd commands for each lambda window. We generated this file in the previous step. While these won't take too long to run, they will take longer than the annihilation simulations, as these are full MD simulations (on the order of an hour). Complex Equilibration ^^^^^^^^^^^^^^^^^^^^^ .. tab-set:: .. tab-item:: Heat .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/complex_mdins/heat.mdin .. tab-item:: Equil1 .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/complex_mdins/equil1.mdin .. tab-item:: Equil2 .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/complex_mdins/equil2.mdin As with the binder, we will generate the mdins programmatically for each lambda window using a bash for loop. .. code-block:: bash cd ../complex # Generate Equilibration MDINs for Each Lambda Window 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 {0..10} do lambda_value=${lambda_values[$i]} sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../complex_mdins/heat.mdin > heat_${i}.mdin sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../complex_mdins/equil1.mdin > equil1_${i}.mdin sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../complex_mdins/equil2.mdin > equil2_${i}.mdin done Now, we also need a groupfile for the replica exchange simulations. This file will look a lot like the pmemd.CUDA command calls that we made before, except the groupfile will contain all of the lambda windows to be run in parallel. The script below will generate the groupfile for the binder equilibration simulations. .. code-block:: bash rm complex_heat_groupfile.in complex_equil1_groupfile.in complex_equil2_groupfile.in for i in {0..10} do echo "-O -i heat_${i}.mdin -p ../../topology/complex_ejm_31_hmr.parm7 -c annihilate_complex_${i}.rst7 -o complex_heat_${i}.mdout -r complex_heat_${i}.rst7 -x complex_heat_${i}.nc -ref ../../topology/complex_ejm_31_unrestrained.rst7" >> complex_heat_groupfile.in echo "-O -i equil1_${i}.mdin -p ../../topology/complex_ejm_31_hmr.parm7 -c complex_heat_${i}.rst7 -o complex_equil1_${i}.mdout -r complex_equil1_${i}.rst7 -x complex_equil1_${i}.nc -ref ../../topology/complex_ejm_31_unrestrained.rst7" >> complex_equil1_groupfile.in echo "-O -i equil2_${i}.mdin -p ../../topology/complex_ejm_31_hmr.parm7 -c complex_equil1_${i}.rst7 -o complex_equil2_${i}.mdout -r complex_equil2_${i}.rst7 -x complex_equil2_${i}.nc -ref ../../topology/complex_ejm_31_unrestrained.rst7" >> complex_equil2_groupfile.in done Now we can run these calculations in parallel using pmemd.cuda.MPI. .. code-block:: bash mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile complex_heat_groupfile.in mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile complex_equil1_groupfile.in mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile complex_equil2_groupfile.in These will take much longer than the binder equilibration simulations, as the complex simulations involve a much larger system size. Expect these to take several hours to complete. .. hint:: If you want to run these in the background, you can use `nohup` or `screen` to keep the processes running after you log out of your session. .. code-block:: bash nohup mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile complex_heat_groupfile.in & Note - if you do this, you need to run these commands for each equilibration stage (heat, equil1, equil2) separately, as each command needs to finish before starting the next one. Running Production Simulations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In the final step, we will run production simulations for each lambda window for both the binder and complex environments. As with the equilibration simulations, we will run a quick equilibration without replica exchange first, followed by the production simulations with replica exchange enabled. .. note:: *On Generating Independent Trials* In practice, ABFE calculations are often run with multiple independent trials to improve convergence and sampling. This can be achieved by starting each trial from different initial conditions, such as different random seeds or different starting structures. In this tutorial, we will demonstrate how to set up and run a single trial for simplicity. However, in a real-world scenario, you would repeat the equilibration and production steps multiple times with different starting conditions to generate independent trials. This approach helps to ensure that the results are robust and not dependent on a single set of initial conditions. There are a few ways to achieve this, such as modifying the random seed in the mdin files or using different starting structures for each trial. The key is to ensure that each trial is independent and samples different regions of conformational space. Binder Production ^^^^^^^^^^^^^^^^^ Like before, we will start with the binder, running each lambda window using the corresponding mdin files provided below. .. tab-set:: .. tab-item:: Trial Equilibration .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/binder_mdins/trial_equil.mdin .. tab-item:: Trial Production .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/binder_mdins/trial_production.mdin We will generate the mdins programmatically for each lambda window using a bash for loop. .. code-block:: bash cd ../binder # Generate Production MDINs for Each Lambda Window 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 {0..10} do lambda_value=${lambda_values[$i]} sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../binder_mdins/trial_equil.mdin > trial_equil_lambda_${i}.mdin sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../binder_mdins/trial_production.mdin > trial_production_lambda_${i}.mdin done We also need a groupfile for the replica exchange simulations, similar to what we did for the equilibration simulations. The script below will generate the groupfile for the binder production simulations. .. code-block:: bash rm binder_trial_equil_groupfile.in binder_trial_production_groupfile.in for i in {0..10} do echo "-O -i trial_equil_${i}.mdin -p ../../topology/binder_ejm_31_hmr.parm7 -c binder_equil2_${i}.rst7 -o binder_trial_equil_${i}.mdout -r binder_trial_equil_${i}.rst7 -x binder_trial_equil_${i}.nc -ref ../../topology/binder_ejm_31_unrestrained.rst7" >> binder_trial_equil_groupfile.in echo "-O -i trial_production_${i}.mdin -p ../../topology/binder_ejm_31_hmr.parm7 -c binder_trial_equil_${i}.rst7 -o binder_trial_production_${i}.mdout -r binder_trial_production_${i}.rst7 -x binder_trial_production_${i}.nc -ref ../../topology/binder_ejm_31_unrestrained.rst7" >> binder_trial_production_groupfile.in done With that, you are now set up to run production simulations for the binder environment for each lambda window. You can run these simulations in parallel using pmemd.cuda.MPI. .. code-block:: bash mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile binder_trial_equil_groupfile.in mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile binder_trial_production_groupfile.in -rem 3 -remlog remd_binder_ejm31.log .. note:: Note the second command includes the flags `-rem 3 -remlog remd_binder_ejm31.log`. These flags enable replica exchange during the production simulations, allowing for enhanced sampling across the lambda windows. The `-rem 3` flag specifies that Hamiltonian molecular dynamics will be used for the replica exchange, while the `-remlog` flag specifies the name of the log file to record the exchange attempts and outcomes. With the binder phase calculated, we can now move on to the complex phase. Complex Production ^^^^^^^^^^^^^^^^^^ .. tab-set:: .. tab-item:: Trial Equilibration .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/complex_mdins/trial_equil.mdin .. tab-item:: Trial Production .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/abfe_files/complex_mdins/trial_production.mdin We will generate the mdins programmatically for each lambda window using a bash for loop like we did before. .. code-block:: bash cd ../complex # Generate Production MDINs for Each Lambda Window 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 {0..10} do lambda_value=${lambda_values[$i]} sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../complex_mdins/trial_equil.mdin > trial_equil_lambda_${i}.mdin sed "s/REPLACE_WITH_LAMBDA/$lambda_value/g" ../../complex_mdins/trial_production.mdin > trial_production_lambda_${i}.mdin done We also need a groupfile for the replica exchange simulations, similar to what we did for the equilibration simulations. The script below will generate the groupfile for the complex production simulations. .. code-block:: bash rm complex_trial_equil_groupfile.in complex_trial_production_groupfile.in for i in {0..10} do echo "-O -i trial_equil_${i}.mdin -p ../../topology/target_ejm_31_hmr.parm7 -c complex_equil2_${i}.rst7 -o complex_trial_equil_${i}.mdout -r complex_trial_equil_${i}.rst7 -x complex_trial_equil_${i}.nc -ref ../../topology/complex_ejm_31_unrestrained.rst7" >> complex_trial_equil_groupfile.in echo "-O -i trial_production_${i}.mdin -p ../../topology/target_ejm_31_hmr.parm7 -c complex_trial_equil_${i}.rst7 -o complex_trial_production_${i}.mdout -r complex_trial_production_${i}.rst7 -x complex_trial_production_${i}.nc -ref ../../topology/complex_ejm_31_unrestrained.rst7" >> complex_trial_production_groupfile.in done With that, you are now set up to run production simulations for both the binder and complex environments for each lambda window. You can run these simulations in parallel using pmemd.cuda.MPI. .. code-block:: bash mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile complex_trial_equil_groupfile.in mpirun -np 11 pmemd.cuda.MPI -ng 11 -groupfile complex_trial_production_groupfile.in -rem 3 -remlog remd_complex_ejm31.log Following the completetion of these steps, you will have successfully set up and run absolute binding free energy calculations for the Tyk2 system. You can then proceed to analyze the results and calculate the binding free energies using appropriate analysis tools. .. end-tutorial