Determining the Boresch Restraint Contribution to the Free Energy
Learning objectives
Learn how to determine the contribution of Boresch restraints to the free energy in ABFE calculations.
Activities
In this Activity, you will learn how to determine the contribution of Boresch restraints to the free energy in ABFE calculations.
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:
or
Here, \(V_0\) is the standard state volume (e.g. 1660 Å^3 for a 1 M standard state), \(r_{a,A,0}\) is the equilibrium distance between the anchor atoms in the Boresch restraint, \(\theta_{A,0}\) and \(\theta_{B,0}\) are the equilibrium angles defined by the anchor atoms and the ligand atom in the Boresch restraint, and \(K_r\), \(K_{\theta_A}\), \(K_{\theta_B}\), \(K_{\phi_A}\), \(K_{\phi_B}\), and \(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”, Boresch[1] Boresch defines a Boresch restraint with the following parameters:
\(K_r\) = 5.0 kcal/mol/Å^2
\(K_{\theta_A}\) = 5.0 kcal/mol/rad^2
\(K_{\phi_A}\) = 2.5 kcal/mol/rad^2
\(r_0\) = 5.1 Å
\(\theta_{A,0}\) = 67.5 degrees
\(\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 \(\Delta A_{rest,BOR2003}\) reported in the paper. The output should be around \(-6.27 kcal/mol\) (which matches the table). Note - when asked about whether the \(0.5\) factor is included in the potential, answer “n” or “no” since the potential used in the Boresch paper does not include the \(0.5\) factor.
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 the Boresch restraints tutorial.
&rst iat=1489,15
r1=0.00000,r2=4.44575,r3=4.44575,r4=999.000,rk2=15.06, rk3=15.06/
&rst iat=1477,1489,15
r1=-180.00000,r2=83.99734,r3=83.99734,r4=180.000,rk2=53.21, rk3=53.21/
&rst iat=1489,15,14
r1=-180.00000,r2=66.38281,r3=66.38281,r4=180.000,rk2=40.55, rk3=40.55/
&rst iat=1490,1477,1489,15
r1=-180.00000,r2=11.03062,r3=11.03062,r4=180.000,rk2=29.83, rk3=29.83/
&rst iat=1477,1489,15,14
r1=-180.00000,r2=30.50819,r3=30.50819,r4=180.000,rk2=62.83, rk3=62.83/
&rst iat=1489,15,14,16
r1=-180.00000,r2=-137.04995,r3=-137.04995,r4=180.000,rk2=49.20, rk3=49.20/
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 \(0.5\) factor since in Amber the potential is defined as \(K*(x-x0)^2\). The value you get should be around \(-10.68\) kcal/mol, which matches the value calculated by FE-Toolkit for this restraint.
Delta A_rest = -10.682981 kcal/mol