Determining the Boresch Restraint Contribution to the Free Energy ================================================================= | 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 the contribution of Boresch restraints to the free energy in ABFE calculations. .. end-learning-objectives Activities ---------- In this Activity, you will learn how to determine the contribution of Boresch restraints to the free energy in ABFE calculations. .. contents:: :local: :depth: 4 .. start-tutorial Determining the Boresch Restraint Contribution to the Free Energy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ To properly account for the contribution of Boresch restraints to the free energy in ABFE calculations, we need to calculate the free energy associated with applying these restraints. This can be done using the following formula: .. math:: \Delta A_{rest,BOR2003} = - k_B T \ln \left\{ \frac{8 \pi^2 V_0}{r_{a,A,0}^2 \sin \theta_{A,0} \sin \theta_{B,0}} \frac{\sqrt{K_r K_{\theta_A} K_{\theta_B} K_{\phi_A} K_{\phi_B} K_{\phi_C}}}{(2\pi k_B T)^3} \right\} or .. math:: \Delta A_{rest,BOR2003} = - k_B T \ln \left\{ \frac{8 \pi^2 V_0}{r_{a,A,0}^2 \sin \theta_{A,0} \sin \theta_{B,0}} \frac{\sqrt{c K_r c K_{\theta_A} c K_{\theta_B} c K_{\phi_A} c K_{\phi_B} c K_{\phi_C}}}{(\pi k_B T)^3} \right\} Here, :math:`V_0` is the standard state volume (e.g. 1660 Å^3 for a 1 M standard state), :math:`r_{a,A,0}` is the equilibrium distance between the anchor atoms in the Boresch restraint, :math:`\theta_{A,0}` and :math:`\theta_{B,0}` are the equilibrium angles defined by the anchor atoms and the ligand atom in the Boresch restraint, 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 distance, angle, and dihedral components of the Boresch restraint. The first formula is used when the potential is defined as K*(x-x0)^2, and the second formula is used when the potential is defined as 0.5*K*(x-x0)^2 (i.e. when there is a 0.5 factor included in the potential). .. warning:: Amber defines the Boresch restraint potential as K*(x-x0)^2, which does not include the 0.5 factor. Therefore, when calculating the contribution of Boresch restraints to the free energy for Amber simulations, you should use the first formula and not include the 0.5 factor in the force constants. In this activity, you will learn how to calculate this free energy contribution using the parameters from an example Boresch restraint for the Tyk2 system. First, however, we need code for calculating the Boresch restraint contribution to the free energy. If you are using FE-Toolkit, there is a script that is used to automatically calculate this value (edgembar/src/python/lib/edgembar/VBA.py); however, we have included a standalone script below that you can use to calculate the contribution by hand. Let's start with a simple example from the literature. In the recent paper by Boresch, "On the Analytical Corrections for Restraints in Absolute Binding Free Energy Calculations", :footcite:t:`Boresch_JChemInfModel_2024_v64_p3605` Boresch defines a Boresch restraint with the following parameters: - :math:`K_r` = 5.0 kcal/mol/Å^2 - :math:`K_{\theta_A}` = 5.0 kcal/mol/rad^2 - :math:`K_{\phi_A}` = 2.5 kcal/mol/rad^2 - :math:`r_0` = 5.1 Å - :math:`\theta_{A,0}` = 67.5 degrees - :math:`\theta_{B,0}` = 84.5 degrees Below is an example python script you can download and try. First try entering the parameters from the Boresch paper to see if you can reproduce the value of :math:`\Delta A_{rest,BOR2003}` reported in the paper. The output should be around :math:`-6.27 kcal/mol` (which matches the table). Note - when asked about whether the :math:`0.5` factor is included in the potential, answer "n" or "no" since the potential used in the Boresch paper does not include the :math:`0.5` factor. .. code-block:: python import numpy as np # Constants kt = 300.0 * 8.314 / 4184.0 # kcal/mol V0 = 1660.0 # Å^3 # User option use_half = input("Does your restraint potential include 0.5*K*(x-x0)^2? [y/n]: ").strip().lower() if use_half in ["y", "yes"]: c = 0.5 elif use_half in ["n", "no"]: c = 1.0 else: raise ValueError("Please enter 'y' or 'n'.") # Inputs r0 = float(input("r0 (Å): ")) thetaA0_deg = float(input("thetaA0 (degrees): ")) thetaB0_deg = float(input("thetaB0 (degrees): ")) k_r = float(input("k_r (kcal/mol/Å^2): ")) k_thetaA = float(input("k_thetaA (kcal/mol/rad^2): ")) k_thetaB = float(input("k_thetaB (kcal/mol/rad^2): ")) k_phiA = float(input("k_phiA (kcal/mol/rad^2): ")) k_phiB = float(input("k_phiB (kcal/mol/rad^2): ")) k_phiC = float(input("k_phiC (kcal/mol/rad^2): ")) # Convert angles to radians thetaA0 = np.deg2rad(thetaA0_deg) thetaB0 = np.deg2rad(thetaB0_deg) # Effective force constants in the Boltzmann factor: # exp[-beta * c*K*(x-x0)^2] k_r_eff = c * k_r k_thetaA_eff = c * k_thetaA k_thetaB_eff = c * k_thetaB k_phiA_eff = c * k_phiA k_phiB_eff = c * k_phiB k_phiC_eff = c * k_phiC # Boresch RRHO correction prefactor = (8 * np.pi**2 * V0) / ( r0**2 * np.sin(thetaA0) * np.sin(thetaB0) ) force_term = np.sqrt( k_r_eff * k_thetaA_eff * k_thetaB_eff * k_phiA_eff * k_phiB_eff * k_phiC_eff ) / ((np.pi * kt) ** 3) delta_A = -kt * np.log(prefactor * force_term) print() print(f"kt = {kt:.8f} kcal/mol") print(f"0.5 factor included: {use_half in ['y', 'yes']}") print(f"Delta A_rest = {delta_A:.6f} kcal/mol") Below, we have included a sample input file for a Boresch restraint that was generated in :doc:`the Boresch restraints tutorial `. .. literalinclude:: /_static/files/ModularTutorials/Alchemical/ABFE_of_Tyk2/abfe_tyk2_tutorial/boresch_files/rest.in :language: text Try entering the parameters from this file into the script above to calculate the contribution of the Boresch restraint to the free energy for your Tyk2 ligand. Here, the parameters should also use "n" for the :math:`0.5` factor since in Amber the potential is defined as :math:`K*(x-x0)^2`. The value you get should be around :math:`-10.68` kcal/mol, which matches the value calculated by FE-Toolkit for this restraint. .. jupyter-execute:: :hide-code: import numpy as np def boresch_restraint_free_energy( r0, thetaA0_deg, thetaB0_deg, k_r, k_thetaA, k_thetaB, k_phiA, k_phiB, k_phiC, use_half=False, temperature=300.0, V0=1660.0, ): """Return the analytical Boresch restraint free energy in kcal/mol. Parameters ---------- r0 : float Equilibrium distance in Å. thetaA0_deg, thetaB0_deg : float Equilibrium angles in degrees. k_r, k_thetaA, k_thetaB, k_phiA, k_phiB, k_phiC : float Force constants in Amber-style units. use_half : bool, default=False True if the restraint potential is defined as 0.5*K*(x-x0)^2. False if it is defined as K*(x-x0)^2, as in Amber. temperature : float, default=300.0 Temperature in K. V0 : float, default=1660.0 Standard-state volume in Å^3. """ kt = temperature * 8.314 / 4184.0 # kcal/mol c = 0.5 if use_half else 1.0 thetaA0 = np.deg2rad(thetaA0_deg) thetaB0 = np.deg2rad(thetaB0_deg) k_r_eff = c * k_r k_thetaA_eff = c * k_thetaA k_thetaB_eff = c * k_thetaB k_phiA_eff = c * k_phiA k_phiB_eff = c * k_phiB k_phiC_eff = c * k_phiC prefactor = (8 * np.pi**2 * V0) / ( r0**2 * np.sin(thetaA0) * np.sin(thetaB0) ) force_term = np.sqrt( k_r_eff * k_thetaA_eff * k_thetaB_eff * k_phiA_eff * k_phiB_eff * k_phiC_eff ) / ((np.pi * kt) ** 3) delta_A = -kt * np.log(prefactor * force_term) return delta_A # Example: values from the literature example delta_A = boresch_restraint_free_energy( r0=4.44575, thetaA0_deg=83.997, thetaB0_deg=66.38, k_r=15.06, k_thetaA=53.21, k_thetaB=40.55, k_phiA=29.83, k_phiB=62.83, k_phiC=49.20, use_half=False, ) print(f"Delta A_rest = {delta_A:.6f} kcal/mol") .. end-tutorial Relevant Literature ------------------- .. footbibliography::