Equilibrating Endstates in 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 equilibrate the endstates in preparation for running absolute binding free energy calculations (ABFE).

.. end-learning-objectives

Relevant literature
-------------------

- Coming Soon!


Tutorial
--------

In this Tutorial, you will learn how to equilibrate the endstates in preparation for running ABFE calculations for TYK2. 

.. contents::
   :local:
   :depth: 4

.. start-tutorial

.. mermaid::

   flowchart LR

      %% Aqueous endstate boxed path
      subgraph AQ[Aqueous Endstate Equilibration]
         direction LR
         A0[Aqueous Start<br>Inputs: binder parm7 and rst7<br>Outputs: Initial aqueous rst7] --> A1
         A1[Minimization<br>Inputs: binder rst7<br>Subprograms: pmemd minimization.mdin<br>Outputs: binder minimization rst7] --> A2
         A2[Heating<br>Inputs: minimized aqueous rst7<br>Subprograms: pmemd heat.mdin<br>Outputs: heated aqueous rst7] --> A3
         A3[Density 1<br>Inputs: heated aqueous rst7<br>Subprograms: pmemd density1.mdin<br>Outputs: density1 aqueous rst7] --> A4
         A4[Density 2<br>Inputs: density1 aqueous rst7<br>Subprograms: pmemd density2.mdin<br>Outputs: density2 aqueous rst7] --> A5
         A5[Lower Restraints 1<br>Inputs: density2 aqueous rst7<br>Subprograms: pmemd lower_restraints1.mdin<br>Outputs: restrained aqueous rst7] --> A6
         A6[Lower Restraints 2<br>Inputs: restrained aqueous rst7<br>Subprograms: pmemd lower_restraints2.mdin<br>Outputs: restrained aqueous rst7] --> A7
         A7[Lower Restraints 3<br>Inputs: restrained aqueous rst7<br>Subprograms: pmemd lower_restraints3.mdin<br>Outputs: restrained aqueous rst7] --> A8
         A8[Lower Restraints 4<br>Inputs: restrained aqueous rst7<br>Subprograms: pmemd lower_restraints4.mdin<br>Outputs: restrained aqueous rst7] --> A9
         A9[Lower Restraints 5<br>Inputs: restrained aqueous rst7<br>Subprograms: pmemd lower_restraints5.mdin<br>Outputs: lightly restrained aqueous rst7] --> A10
         A10[Backbone Binder Restraints<br>Inputs: lightly restrained aqueous rst7<br>Subprograms: pmemd backbone.mdin<br>Outputs: backbone restrained aqueous rst7] --> A11
         A11[Unrestrained MD<br>Inputs: backbone restrained aqueous rst7<br>Subprograms: pmemd free.mdin<br>Outputs: equilibrated aqueous rst7]
      end

      %% Complex endstate boxed path
      subgraph CX[Complex Endstate Equilibration]
         direction LR
         B0[Complex Start<br>Inputs: target parm7 and rst7<br>Outputs: Initial complex rst7] --> B1
         B1[Minimization<br>Inputs: complex rst7<br>Subprograms: pmemd minimization.mdin<br>Outputs: complex minimization rst7] --> B2
         B2[Heating<br>Inputs: minimized complex rst7<br>Subprograms: pmemd heat.mdin with protein and ligand restraints<br>Outputs: heated complex rst7] --> B3
         B3[Density 1<br>Inputs: heated complex rst7<br>Subprograms: pmemd density1.mdin<br>Outputs: density1 complex rst7] --> B4
         B4[Density 2<br>Inputs: density1 complex rst7<br>Subprograms: pmemd density2.mdin<br>Outputs: density2 complex rst7] --> B5
         B5[Lower Restraints 1<br>Inputs: density2 complex rst7<br>Subprograms: pmemd lower_restraints1.mdin<br>Outputs: restrained complex rst7] --> B6
         B6[Lower Restraints 2<br>Inputs: restrained complex rst7<br>Subprograms: pmemd lower_restraints2.mdin<br>Outputs: restrained complex rst7] --> B7
         B7[Lower Restraints 3<br>Inputs: restrained complex rst7<br>Subprograms: pmemd lower_restraints3.mdin<br>Outputs: restrained complex rst7] --> B8
         B8[Lower Restraints 4<br>Inputs: restrained complex rst7<br>Subprograms: pmemd lower_restraints4.mdin<br>Outputs: restrained complex rst7] --> B9
         B9[Lower Restraints 5<br>Inputs: restrained complex rst7<br>Subprograms: pmemd lower_restraints5.mdin<br>Outputs: lightly restrained complex rst7] --> B10
         B10[Backbone Binder Restraints<br>Inputs: lightly restrained complex rst7<br>Subprograms: pmemd backbone.mdin<br>Outputs: backbone restrained complex rst7] --> B11
         B11[Unrestrained MD<br>Inputs: backbone restrained complex rst7<br>Subprograms: pmemd free.mdin<br>Outputs: equilibrated complex rst7]
      end

      OUT['Final Equilibrated Systems']

      A11 --> OUT
      B11 --> OUT


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:


- :download:`Download equil_files.tar.gz from <https://drive.google.com/drive/folders/1RmAqaQhzcij8bMIafnlaBmFub26k0vCP?usp=drive_link>`

Once you have downloaded the files, extract them into a directory where you will be working on this tutorial.

.. code-block:: bash

    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.

.. tab-set::

    .. tab-item:: Minimization

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/minimization.mdin

    .. tab-item:: Heating

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/heat.mdin

    .. tab-item:: Density 1

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/density1.mdin
    
    .. tab-item:: Density 2

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/density2.mdin

    .. tab-item:: Lower Protein Restraints 1

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/lower_restraints1.mdin

    .. tab-item:: Lower Protein Restraints 2

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/lower_restraints2.mdin

    .. tab-item:: Lower Protein Restraints 3

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/lower_restraints3.mdin

    .. tab-item:: Lower Protein Restraints 4

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/lower_restraints4.mdin

    .. tab-item:: Lower Protein Restraints 5

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/lower_restraints5.mdin

    .. tab-item:: Backbone Binder
        
        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/backbone.mdin

    .. tab-item:: Unrestrained

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/binder_mdins/free.mdin

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.

.. code-block:: bash

    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.

    .. code-block:: bash

        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.

.. tab-set::

    .. tab-item:: Minimization

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/minimization.mdin

    .. tab-item:: Heating

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/heat.mdin

    .. tab-item:: Density 1

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/density1.mdin
    
    .. tab-item:: Density 2

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/density2.mdin

    .. tab-item:: Lower Protein Restraints 1

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/lower_restraints1.mdin

    .. tab-item:: Lower Protein Restraints 2

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/lower_restraints2.mdin

    .. tab-item:: Lower Protein Restraints 3

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/lower_restraints3.mdin

    .. tab-item:: Lower Protein Restraints 4

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/lower_restraints4.mdin

    .. tab-item:: Lower Protein Restraints 5

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/lower_restraints5.mdin

    .. tab-item:: Backbone Binder
        
        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/backbone.mdin

    .. tab-item:: Unrestrained

        .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/equil_files/complex_mdins/free.mdin

The commands to run each of these steps using amber are included below.

.. code-block:: bash

    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.

.. end-tutorial