On Analytical Corrections for Restraints in Absolute Binding Free Energy Calculations

Double decoupling absolute binding free energy simulations require an intermediate state at which the ligand is held solely by restraints in a position and orientation resembling the bound state. One possible choice consists of one distance, two angle, and three dihedral angle restraints. Here, I demonstrate that in practically all cases the analytical correction derived under the rigid rotator harmonic oscillator approximation is sufficient to account for the free energy of the restraints.

Eq. 1 (= Eq. 4 of Ref. 1) can be obtained with the single Mathematica command listed as Eq. 1b below.The corresponding rigid rotator harmonic oscillator approximation is Angle restraint: Here, Chen et al. also use an approximation, 1 i.e., changing the limits of integration from (0, π) to ±∞.
The following Mathematica command leads immediately to the result given in Eq. 2. A full derivation is given in the Section Exact solution for the angle configurational integral below.
Dihedral angle restraint Here, since there is no Jacobian factor, the only approximation needed to obtain Z RR φ is extending the limits of integration to ±∞.
For completeness, the Mathematica expression to obtain the exact expression for Z φ :

Hints on obtaining the numerical results
The numerical examples given in the main text are based on some results/tests provided in Ref.
Eq. 4 was derived using U (x) = K /2(x − x 0 ) 2 as the expression for the harmonic restraint potential function.If one does not include the factor 1/2, i.e., U (x) = K(x − x 0 ) 2 , this leads to a slightly modified expression for the restraint contribution Eq. 5 is more in line with the internal workings of CHARMM, 4 which uses U (x) = K(x−x 0 ) 2 for harmonic potentials, i.e., the factor 1/2 is absent.The force constants reported in Ref. 3 are the values for the CHARMM internal functional form without the factor 1/2.Therefore, upon inserting numbers in Eq. 4 (Eq.32 of Ref. 3), one has to double each force constant since K = 2K.For example, when using Eq. 4, one has to employ force constants of 40 kcal/(mol Å 2 ) [40 kcal/(mol rad 2 )], i.e., 2×20 kcal/(mol Å 2 ) [2×20 kcal/(mol rad 2 )], to reproduce the CL1 results in Tables 3 and 4 of Ref. 3 and in Table 1 of the main text.By contrast, when using the RR approximations given in Eqs.1-3 above, which lead to Eq. 5, one can use the force constants of 20 kcal/(mol Å 2 ) [20 kcal/(mol rad 2 )] "as is" to obtain the CL1 result.
( To use the harmonic oscillator approximation, the first order Taylor approximation was used, i.e., Note that because of the Taylor formula, the factor 1/2 appears in the harmonic approximation for the dihedral restraint.In these cases, the situation is reversed.
Using Eq. 4, the multiplication by 2 must not be applied to the dihedral force constants.
Consider, e.g., restraint L1A-F, Table 5 of Ref. 3, for which we also report results in Table 1 of the main text.Here, different force constants for the individual terms were used; that is, the distance restraint had a force constant of 4 kcal/(mol Å 2 ), the two angle restraints had a force constant of 8 kcal/(mol rad 2 ), and the three dihedral angle restraints had a force constant of 10 kcal/(mol rad 2 ).When calculating ∆A RR rest using Eq. 4 (Eq.32 of Ref 3), the distance and angle force constants need to be doubled, whereas the dihedral force constants can be used as is.By contrast, when utilizing the expressions derived here (i.e., Eq. 5 above), the force constants for the distance and angle restraints can be used as is, but one has to scale the dihedral force constants by 1/2 (k φ = 10 kcal/(mol rad 2 ) corresponds to k φ = k φ /2 = 5 kcal/(mol rad 2 )).

Exact solution for the angle configurational integral Z θ
For completeness, we derive here the exact expression for the configurational integral of the angle restraint.Modern CAS will provide closed solutions for both the antiderivative, and for the full integral from 0 to π, but the resulting expressions are difficult to understand.
Furthermore, the result contains error functions of complex arguments, suggesting that the configurational integral might be complex.Breaking the problem into a series of simpler integrals and exploiting properties of the error function, one not only obtains a more readable expression but can show that the result is purely real.
We start with the substitution t = θ − θ 0 , i.e., and now focus on the antiderivative Using a CAS, e.g., Mathematica 2 with the rule-based integration system 5 loaded, 1 one finds so if the limits of integration are extended from 0 to π to ±∞, the imaginary part of the error function vanishes and the real part gives +1 for +∞, −1 for −∞.Hence, we obtain immediately which is Eq. 5 of Ref. 1.This approximation can be derived in various ways more easily; any modern CAS gives the result directly when using ±∞ as limits of integration.Nevertheless, it is a consistency check that the result can also be obtained in a straightforward manner from the exact antiderivative.
Two additional observations follow from Eq. 7 in combination with Eqs.11 and 12.
First, if the equilibrium angle θ 0 is near 90 • (π/2), then cos θ 0 is close to 0 and most of the contributions to Z θ will come from the sin θ 0 × I2 term in Eq. 7. Second, compared to the exact expression Eq. 6, both the RR approximation and Schrödinger's approximation (cf. Eq. 2) are incorrect for θ 0 = 0 and θ 0 = π.Both approximations give zero in these cases, whereas the exact expression leads to a non-zero value.In practice, however, target angle values near 0 or π need to be avoided anyway to avoid singularities in the forces, 1,7 so in practice this limitation of both approximate expressions can be safely ignored.
3. When operating with the expressions provided by Boresch et al. (specifically, (1) We reproduce here Eq. 32 of Ref. 3 with minor adjustments of the nomenclature