Creating Boresch Restraints 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 determine Boresch restraints from end-state trajectories. .. end-learning-objectives Relevant literature ------------------- - Coming Soon! Tutorial -------- In this Tutorial, you will learn how to create Boresch restraints for use in Absolute Binding Free Energy (ABFE) calculations using the Tyk2 system as an example. .. contents:: :local: :depth: 4 .. start-tutorial Determining Boresch Restraints ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For this tutorial, we will be starting from the equilibrated end-state trajectory of a Tyk2-ligand complex in its bound state. These files were generated in the tutorial:ref:`Equilibrating the End States of Tyk2 for ABFE Calculations `; however, for 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. Theory of Boresch Restraints ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Boresch restraints are a type of positional restraint used in molecular dynamics simulations to maintain the relative position and orientation of a ligand within a binding site. They are particularly useful in absolute binding free energy calculations, where it is essential to keep the ligand in place while decoupling its interactions with the environment. The Boresch restraint method involves defining a set of six restraints based on three atoms from the ligand and three atoms from the protein. These restraints include: - One distance Restraint - Two Angle Restraints - Three Dihedral Restraints These restraints work together to fix the ligand's position and orientation relative to the protein, allowing for accurate free energy calculations. .. note:: This theory section should be expanded with a more detailed explanation of the mathematical formulation of Boresch restraints, including equations and diagrams to illustrate the concepts. .. note:: Much of the below theoretical discussion comes from Boresch et al., J. Phys. Chem. B, 2003, 107 (35), pp 9535–9551. DOI: `10.1021/jp0217839 `_. The goal of Boresch restraints is to define a known set of lambda-dependent restraints that can be applied to the ligand in the bound state to keep it in place while decoupling its interactions with the environment, while also allowing for an analytical correction to the free energy to account for the restraints. Roughly, the situation that needs to be corrected for is: .. math:: (P \cdot L)_{H_2O} \rightarrow^{\Delta A_r} (P)_{H_2O} + (L)_{H_2O}, Where here, (P \cdot L)_{H_2O} is the protein-ligand complex in water with restraints applied, (P)_{H_2O} is the protein in water, and (L)_{H_2O} is the ligand in water. $\Delta A_r$ is the free energy change associated with removing the restraints from the complex in water. Because the ligand is in the dummy state when decoupled, the Ligand does not have interactions with the protein or solvent and thus the restraints are the only interactions with external sources it feels. The potential energy function for the system is: .. math:: U = U_P(r_P) + U_P(r_L) + U_{r}(r_P, r_L). Here, there are terms coming from the protein (:math:`U_P`), the ligand (:math:`U_L`), and the restraints (:math:`U_r`). The contribution to the free energy that comes from the restraints can be written as: .. math:: \Delta A_r = -k_B T \ln \left( \frac{Z_P Z_L}{Z_{P \cdot L}} \right) Where :math:`Z_{P \cdot L}` is the partition function of the restrained complex, :math:`Z_P` is the partition function of the protein, and :math:`Z_L` is the partition function of the ligand. This calculation is feasible analytically when the restraints are chosen such that .. math:: Z_{P \cdot L} = Z_P \times Z_L \times Z_r, or in other words, choosing the restraints such that the partition function can be factored into independent contributions from the protein, ligand, and restraints. After a few assumptions (Rigid Rotor) the final expression for the partition function is given by: .. math:: Z_{P \cdot L} = r^2_{aA,0} \sin(\theta_{A,0}) \frac{(2 \pi k_B T)^3}{\sqrt{K_r K_{\theta^A} K_{\theta^B} K_{\phi^A} K_{\phi^B} K_{\phi^C}}}, Finally, the free energy contribution from the restraints can be expressed as: .. math:: \Delta A_r = -k_B T \ln \left[\frac{8\pi^2V\sqrt{K_rK_{\theta^A}K_{\theta^B}K_{\phi^A}K_{\phi^B}K_{\phi^C}}}{r^2_{aA,0} \sin(\theta_{A,0}) \sin(\theta_{B,0}) (2\pi k_B T)^3} \right] Where :math:`r_{aA,0}` and :math:`\theta_{A,0}` are the equilibrium distance and angle values for the restraints, and :math:`K_r`, :math:`K_{\theta^A}`, :math:`K_{\theta^B}`, :math:`K_{\phi^A}`, :math:`K_{\phi^B}`, and :math:`K_{\phi^C}` are the force constants for the respective restraints. How to Determine Boresch Restraints ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ To determine the Boresch restraints for your system, you will need to analyze an equilibrated trajectory for the protein-ligand complex. .. note:: If you have been following along with the Tyk2 ABFE Tutorials, it is very likely that you have already generated a trajectory. It very likely that your trajectory will be different from the one included here, but the steps to generate the boresch restraints will be the same. Don't be surprised if your final results differ from those presented here! For this tutorial, we will be using the following python code, which you should save as `generate_boresch_restraints.py`: .. toggle:: .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/boresch_files/generate_boresch_restraints.py This code uses the MDAnalysis package to analyze the trajectory and select appropriate atoms for the Boresch restraints. You can run the code with the following command: .. code-block:: bash python generate_boresch_restraints.py This will create a couple of image files, as well as a "rest.in" file that contains the necessary restraint information for use in an Amber simulation. There are two image files that are generated by the script. The first is `boresch_histograms.png`, which shows the histograms of the distance, angle, and dihedral values for the selected atoms over the course of the trajectory. This is useful for visualizing how well the selected atoms maintain their relative positions and orientations. .. image:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/boresch_files/boresch_histograms.png :alt: Boresch Histograms :align: center :width: 60% The second image is `boresch_timeseries.png`, which shows the time series of the distance, angle, and dihedral values for the selected atoms over the course of the trajectory. This is useful for visualizing how stable the selected atoms are over time. .. image:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/boresch_files/boresch_timeseries.png :alt: Boresch Time Series :align: center :width: 60% The final output is the `rest.in` file, which contains the necessary restraint information for use in an Amber simulation. This file includes the atom indices, equilibrium values, and force constants for the Boresch restraints. An example of the `rest.in` file generated is shown below: .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/boresch_files/rest.in :language: text .. tip:: The rest.in file should look familiar to you if you are familiar with Amber's NMR restraint format, as Boresch restraints are implemented using the same underlying machinery. Specifically, each &rst block in the rest.in file corresponds to a specific restraint type (distance, angle, dihedral) and contains the necessary parameters for that restraint. You can tell what type of restraint it is based on the number of atoms involved: - Distance restraints involve 2 atoms - Angle restraints involve 3 atoms - Dihedral restraints involve 4 atoms So the first line in rest.in is a distance restraint, since iat=1489,15 Now, the line after the &rst marker specifies the restraint values. Specifically, the fields always are: - r1: the lower bound where the restraint should reach infinity - r2: a mid-point of the restraint well where it should start taking effect (on its way to r1) - r3: another mid-point of the restraint well where it should start taking effect (on its way to r4) - r4: the upper bound where the restraint should reach infinity - rk2: the force constant for the left side of the well (between r1 and r2) - rk3: the force constant for the right side of the well (between r3 and r4) .. note:: In between r2 and r3, there is no force applied (though, here you'll note for all restraints r2=r3, meaning the well is harmonic in all cases.) With these boresch restraints generated, you are now ready to use them in your ABFE calculations! .. end-tutorial