Relative Binding Free Energy Calculations for Ligands with Diverse Scaffolds with the Alchemical Transfer Method

We review here the theory underpinning the Alchemical Transfer Method (ATM). The presentation below tries to be mathematically rigorous, however simplified notation is often employed to unclutter the equations. For example, in intermediate formulas, we often omit limits of integration and Jacobian factors for curvilinear coordinates when they do not affect the form and interpretation of the final result. To shorten the expressions of the partition functions, the degrees of freedom of the solvent are often not indicated and included by means of the solvent potential of mean force $\Psi$.

We start with the statistical mechanics expression for the bimolecular binding constant between a receptor molecule $R$ and a ligand $A$^{43, 44, 14}

$$K_{b}(A)=\frac{C^{\circ}}{8\pi^{2}}\frac{z_{RA}}{z_{R}z_{A}},$$

(1)

where

$$z_{A}=\int dx_{A}e^{-\beta\Psi_{A}(x_{A})}$$

(2)

and similarly for R, and

$$z_{RA}=\int dx_{R}dx_{A}d\zeta_{A}I(\zeta_{A})e^{-\beta\Psi_{RA}(x_{R},x_{A},\zeta_{A})}$$

(3)

where $x_{R}$ and $x_{A}$ are the internal coordinates of the receptor and ligand respectively, $\zeta_{A}=({\mathbf{c}}_{A},{\mathbf{\omega}}_{A})$ denotes, collectively, the three position and three orientation coordinates of ligand $A$ relative, in this case, to the reference frame of the receptor (Fig. 1), $\Psi_{A}$ and $\Psi_{RA}$ denote the effective potential energy functions in the solvent potential of mean force representation^{45, 14} of the ligand in the solvent and of the ligand bound to the solvated receptor, respectively, The function $I(\zeta_{A})$ in Eq. (3) is an indicator function that is set to $1$ if the position and orientation of the ligand are such that receptor and ligand are considered bound, and zero otherwise.^{43, 46, 44} The volume of configurational space encompassed by the indicator function is denoted here by $V_{\rm site}\Omega_{\rm site}$:

$$\int d\zeta I(\zeta)=V_{\rm site}\Omega_{\rm site}$$

(4)

Eq. (1) is turned into a statistical mechanics average by first considering the product $z_{R}z_{A}$ as the partition function of a system containing both $R$ and $A$, with $A$ located at a position $\zeta_{A}^{\ast}$ in the solvent bulk far away from $R$ so that it does not interact with it:

$$z_{R}z_{A}=\int dx_{R}dx_{A}e^{-\beta\Psi_{RA}(x_{R},x_{A},\zeta_{A}^{\ast})}$$

(5)

and then multiplying it and dividing it by

$$\int d\zeta_{A}^{\ast}I^{\ast}(\zeta_{A}^{\ast})=V_{\rm site}\Omega_{\rm site}$$

(6)

where

$$I^{\ast}({\mathbf{c}},{\mathbf{\omega}})=I({\mathbf{c}}-{\mathbf{d}},{\mathbf{\omega}})$$

(7)

is an indicator function identical to the one that defines the complex, but centered at a point ${\mathbf{d}}$ in the solvent bulk relative to the reference frame of the receptor.

After these manipulations, Eq. (1) becomes

$$K_{b}(A)=C^{\circ}V_{\rm site}\frac{\Omega_{\rm site}}{8\pi^{2}}\frac{\int dx_{R}dx_{A}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}I({\mathbf{c}}_{A},{\mathbf{\omega}}_{A})e^{-\beta\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A})}}{\int dx_{R}dx_{A}d{\mathbf{c}}^{\ast}_{A}d{\mathbf{\omega}}_{A}I^{\ast}({\mathbf{c}}^{\ast}_{A},{\mathbf{\omega}}_{A})e^{-\beta\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}^{\ast}_{A},{\mathbf{\omega}}_{A})}}$$

(8)

where the external coordinates $\zeta$ are expanded into their position and orientation components and, without loss of generality, the orientation dummy variables are denoted by ${\mathbf{\omega}}$ in both integrals. Next, we perform the change of variable ${\mathbf{c}}^{\ast}_{A}={\mathbf{c}}_{A}+{\mathbf{d}}$, noting that, using Eq. (7), $I^{\ast}({\mathbf{c}}^{\ast}_{A},{\mathbf{\omega}}_{A})=I^{\ast}({\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A})=I({\mathbf{c}}_{A},{\mathbf{\omega}}_{A})$ to obtain

$$K_{b}(A)=C^{\circ}V_{\rm site}\frac{\Omega_{\rm site}}{8\pi^{2}}\frac{\int dx_{R}dx_{A}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}I({\mathbf{c}}_{A},{\mathbf{\omega}}_{A})e^{-\beta\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A})}}{\int dx_{R}dx_{A}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}I({\mathbf{c}}_{A},{\mathbf{\omega}}_{A})e^{-\beta\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A})}}$$

(9)

Finally, by multiplying and dividing the integrand in the numerator by
$\exp[-\beta\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A})]$, Eq. (9) becomes

$$K_{b}(A)=C^{\circ}V_{\rm site}\frac{\Omega_{\rm site}}{8\pi^{2}}\langle e^{-\beta u}\rangle_{0}$$

(10)

where $u$ is the perturbation function

$$u=\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A})-\Psi_{RA}(x_{R},x_{A},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A})$$

(11)

corresponding to the change in potential energy for rigidly translating the ligand from a position in the bulk ${\mathbf{c}}_{A}+{\mathbf{d}}$ to a position in the binding site ${\mathbf{c}}_{A}$, and the ensemble average $\langle\ldots\rangle_{0}$ is conducted when the ligand is in the bulk.

Eqs. (10) and (11) form the basis of the Alchemical Transfer Method (ATM) for the estimation of Absolute Binding Free Energies (ABFE).²⁷

Consider now the ratio, $K_{b}({\rm B})/K_{b}({\rm A})$, of the equilibrium binding constants of two ligands A and B to the same receptor R in the same binding mode described by the indicator function $I(\zeta)$. To simplify the notation, in the following we will also assume that $I(\zeta)$ is independent of the orientation coordinates. That is we assume that $I(\zeta)=I({\mathbf{c}})$. Under this assumption $\Omega_{\rm site}=8\pi^{2}$ which cancels the same factor in the denominator of Eq. (10).

From Eq. (1) and canceling the common factor $z_{\rm R}$, we have

$$\frac{K_{b}({\rm B})}{K_{b}({\rm A})}=\frac{z_{{\rm RB}}z_{\rm A}}{z_{{\rm RA}}z_{\rm B}}\,,$$

(12)

where $z_{\rm RB}$ is the intramolecular configurational partition function of the complex between R and B, $z_{\rm B}$ is the intramolecular partition function of ligand B, and similarly for $z_{\rm RA}$ and $z_{\rm A}$. Analogously to manipulations that lead to Eq. (8), the products $z_{{\rm RB}}z_{\rm A}$ and $z_{{\rm RA}}z_{\rm B}$ are expressed as partition function integrals of systems containing the receptor and the two ligands such that one of the ligands is in the binding site and the other is centered on a position of the bulk at distance ${\mathbf{d}}$ relative to the receptor coordinate frame. Specifically we multiply the numerator of Eq. (12) by Eq. (6) for A and the denominator by the corresponding integral for B. The result is

$$\frac{K_{b}(B)}{K_{b}(A)}=\frac{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}^{\ast}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}_{B}d{\mathbf{\omega}}_{B}I^{\ast}({\mathbf{c}}^{\ast}_{A})I({\mathbf{c}}_{B})e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}^{\ast}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}_{B},{\mathbf{\omega}}_{B})}}{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}^{\ast}_{B}d{\mathbf{\omega}}_{B}I({\mathbf{c}}_{A})I^{\ast}({\mathbf{c}}^{\ast}_{B})e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B},{\mathbf{\omega}}_{B})}}$$

(13)

Next, we perform the changes of variables ${\mathbf{c}}^{\ast}_{A}={\mathbf{c}}_{A}+{\mathbf{d}}$ and ${\mathbf{c}}_{B}={\mathbf{c}}^{\ast}_{B}-{\mathbf{d}}$ for the integral in the numerator to obtain

$$\frac{K_{b}(B)}{K_{b}(A)}=\frac{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}^{\ast}_{B}d{\mathbf{\omega}}_{B}I({\mathbf{c}}_{A})I^{\ast}({\mathbf{c}}^{\ast}_{B})e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B}-{\mathbf{d}},{\mathbf{\omega}}_{B})}}{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}^{\ast}_{B}d{\mathbf{\omega}}_{B}I({\mathbf{c}}_{A})I^{\ast}({\mathbf{c}}^{\ast}_{B})e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B},{\mathbf{\omega}}_{B})}}$$

(14)

Again, by multiplying and dividing the integrand in the numerator by the Boltzmann factor in the denominator, Eq. (14) is expressed by an ensemble average

$$\frac{K_{b}(B)}{K_{b}(A)}=\langle e^{-\beta u}\rangle_{0}$$

(15)

with the perturbation energy

$$u=\Psi(x_{R},x_{A},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B}-{\mathbf{d}},{\mathbf{\omega}}_{B})-\Psi(x_{R},x_{A},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B},{\mathbf{\omega}}_{B})$$

(16)

corresponding to the change in potential energy for rigidly translating ligand B from a position in the bulk ${\mathbf{c}}^{\ast}_{B}$ to a position in the binding site ${\mathbf{c}}_{B}-{\mathbf{d}}$ while at the same time ligand A is translated from its position in the binding site ${\mathbf{c}}_{A}$ to a position in the bulk ${\mathbf{c}}_{A}+{\mathbf{d}}$, and the ensemble average $\langle\ldots\rangle_{0}$ is conducted while ligand A is the binding site and ligand B is in the bulk.

Eqs. (15) and (16) are the basis of the Relative Binding Free Energy (RBFE) estimation protocol of the Alchemical Transfer Method (ATM) presented in this work. In analogy with the double-topology free energy perturbation framework,⁷ we call this approach dual-ligand to emphasize the fact that the calculations is performed with a solvent box containing the two ligands.

While, Eq. (15) is formally correct, it converges slowly near the end-points because placing ligand B anywhere in the binding site and placing ligand A anywhere in the solvent bulk is likely to produce severe clashes with either receptor atoms or solvent molecules. For ligands of similar shape, the occurrences of severe clashes with nearby atoms would be reduced if the two ligands land in a similar location and orientation when their positions are exchanged. In order to improve convergence, we introduce a position restraint potential $U_{r}(c_{AB})$, as well as an orientation restraint potential $U_{\Omega}(\omega_{AB})$ to keep the two ligands at approximately the same distance and approximately aligned with each other. Here, $c_{AB}$ represents the distance between the reference frames of the two ligands and $\omega_{AB}$ represents, schematically, their mutual orientation (Figure 1). The specific form of the alignment restraint potentials used in this work is given below. For now, we simply assume that they are symmetric with respect to the exchange of the A and B ligand labels and we show that, under this assumption, the ATM RBFE formulation of Eqs. (15) and (16) is independent of the specific alignment restraint potential.

To see why this is the case, note that to obtain Eq. (13), we multiplied and divided Eq. (12) by Eq. (6). By instead multiplying the numerator by the integrals

$$V_{r}=\int d{\mathbf{c}}^{\ast}_{A}I^{\ast}({\mathbf{c}}^{\ast}_{A})e^{-\beta U_{r}(c_{AB})}$$

(17)

and

$$\Omega_{r}=\int d{\mathbf{\omega}}_{A}e^{-\beta U_{\Omega}(\omega_{AB})}$$

(18)

and by multiplying the denominator by the corresponding integrals for ligand B that have the same values as those for ligand A, we obtain (notice that the integrals $V_{r}$ and $\Omega_{r}$ above and the corresponding ones for ligand B, are constants independent of the external degrees of freedom of the ligand that is not integrated over)

$$\frac{K_{b}(B)}{K_{b}(A)}=\frac{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}^{\ast}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}_{B}d{\mathbf{\omega}}_{B}I^{\ast}({\mathbf{c}}^{\ast}_{A})I({\mathbf{c}}_{B})e^{-\beta[U_{r}(c_{AB})+U_{\Omega}(\omega_{AB})]}e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}^{\ast}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}_{B},{\mathbf{\omega}}_{B})}}{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}^{\ast}_{B}d{\mathbf{\omega}}_{B}I({\mathbf{c}}_{A})I^{\ast}({\mathbf{c}}^{\ast}_{B})e^{-\beta[U_{r}(c_{AB})+U_{\Omega}(\omega_{AB})]}e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B},{\mathbf{\omega}}_{B})}}$$

(19)

As before, we set ${\mathbf{c}}^{\ast}_{A}={\mathbf{c}}_{A}+{\mathbf{d}}$ and ${\mathbf{c}}_{B}={\mathbf{c}}^{\ast}_{B}-{\mathbf{d}}$ in the numerator, then use Eq. (7) and the symmetry property of the restraint potential functions to obtain:

$$\frac{K_{b}(B)}{K_{b}(A)}=\frac{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}^{\ast}_{B}d{\mathbf{\omega}}_{B}I({\mathbf{c}}_{A})I^{\ast}({\mathbf{c}}^{\ast}_{B})e^{-\beta[U_{r}(c_{AB})+U_{\Omega}(\omega_{AB})]}e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B}-{\mathbf{d}},{\mathbf{\omega}}_{B})}}{\int dx_{R}dx_{A}dx_{B}d{\mathbf{c}}_{A}d{\mathbf{\omega}}_{A}d{\mathbf{c}}^{\ast}_{B}d{\mathbf{\omega}}_{B}I({\mathbf{c}}_{A})I({\mathbf{c}}^{\ast}_{B})e^{-\beta[U_{r}(c_{AB})+U_{\Omega}(\omega_{AB})]}e^{-\beta\Psi(x_{R},x_{A},x_{B},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B},{\mathbf{\omega}}_{B})}}$$

(20)

which is the same as Eq. (14) with the addition of the alignment restraint potential terms. The corresponding ensemble average expression is the same as Eq. (15) with Eq. (16), except that A and B remain aligned. Eq. (20) confirms that the ATM RBFE formula Eq. (15) yields an unbiased estimate of the relative binding free energy with or without alignment restraints.

The specific forms of the restraining potentials used in this work are as follows. The positional restraint is a spherical harmonic potential between the reference atoms $1$ and $1^{\prime}$ of the two ligands:

$$U_{r}(c_{AB})=\frac{1}{2}k_{r}|{\mathbf{r}}^{\ast}_{1A}-{\mathbf{d}}-{\mathbf{r}}_{1B}|^{2}$$

(21)

where ${\mathbf{r}}^{\ast}_{1A}$ is the position of the first reference atom of ligand A in the bulk and ${\mathbf{r}}_{1B}$ is the position of the first reference atom of ligand B in the binding site. Given the relationship ${\mathbf{r}}^{\ast}={\mathbf{r}}+{\mathbf{d}}$, this potential is symmetric with respect to the exchange of A and B. The orientation restraints are implemented by two potentials. One, named $U_{\theta}$, is based on the angle $\theta$ between the $z$-axis of the reference frame of A relative to that of B (Fig. 1), which is naturally symmetric with respect to the exchange of the two frames. The second restraint is based on the dihedral angle $\psi$ in Figure 1. Because the $\psi$ angle is not symmetric with respect to the exchange of A and B, we use a symmetrized restraining potential of the form

$$\frac{1}{2}[U_{\psi}(\psi_{BA})+U_{\psi}(\psi_{AB})]$$

(22)

where $U_{\psi}$ is the restraining potential below, $\psi_{BA}$ and $\psi_{AB}$ are the $\psi$ angles of the reference frame of ligand B with respect to A and of A relative to B, respectively. Each orientation restraint potential is based on the cosine of the angle. For example

$$U_{\theta}(\theta)=\frac{1}{2}k_{\theta}(1-\cos\theta)$$

(23)

and similarly for the $\psi$ angle, where $k_{\theta}$ (and $k_{\psi}$ in the case of the $\psi$ angle) is a dimensionless factor that controls the strength of the restraint.

For convenience of notation, Eqs. (15) and (16), were derived in the solvent potential of mean force formalism to avoid carrying the degrees of freedom of the solvent. When the definition of the solvent potential of mean force is substituted into Eq. (14) to show explicitly the degrees of freedom of the solvent,⁴⁴ the RBFE formula in Eq. (15) is unchanged and the formula for the perturbation energy Eq. (16) becomes

$$u=U(x_{R},x_{A},{\mathbf{c}}_{A}+{\mathbf{d}},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B}-{\mathbf{d}},{\mathbf{\omega}}_{B},x_{S})-U(x_{R},x_{A},{\mathbf{c}}_{A},{\mathbf{\omega}}_{A},{\mathbf{c}}^{\ast}_{B},{\mathbf{\omega}}_{B},x_{S})$$

(24)

in which the actual potential energy function $U$ replaces the effective potential energy function $\Psi$, and, in addition to the coordinates of the receptor and the ligands, it includes the coordinates, $x_{S}$, of the solvent.

While Eqs. (15) and (16) formally yield the correct relative binding free energy, their direct implementation in terms of the exponential average of the perturbation energy is numerically unstable.¹⁴ Instead, the free energy change is computed by well-established protocols based on stratification using a sequence of alchemical states defined by a $\lambda$-dependent Hamiltonian⁴⁷ and multi-state thermodynamic reweighting methods.^{40, 41} The specific alchemical protocol we employ is described in detail in previous publications.^{48, 49, 27} Briefly, in this work we employ a linear $\lambda$-dependent alchemical potential energy function of the form

$$U_{\lambda}(x)=U_{0}(x)+W_{\lambda}(u)$$

(25)

where $x$ represents the set of atomic coordinates of the system (the internal coordinates of the ligands and the receptor, the external coordinates of the ligands, and the coordinates of the solvent), the perturbation energy $u(x)$ is defined by Eq. (24) and $W_{\lambda}(u)$ is the generalized softplus alchemical perturbation function

$$W_{\lambda}(u)=\frac{\lambda_{2}-\lambda_{1}}{\alpha}\ln\left[1+e^{-\alpha(u_{\rm sc}(u)-u_{0})}\right]+\lambda_{2}u_{\rm sc}(u)+w_{0}.$$

(26)

The parameters $\lambda_{2}$, $\lambda_{1}$, $\alpha$, $u_{0}$, and $w_{0}$ are functions of $\lambda$ which are adjusted to enhance convergence,⁴⁹ and

$$u_{\rm sc}(u)=\begin{cases}u&u\leq u_{c}\\ (u_{\rm max}-u_{c})f_{\rm sc}\left[\frac{u-u_{c}}{u_{\rm max}-u_{c}}\right]+u_{c}&u>u_{c}\end{cases}$$

(27)

with

$$f_{\text{sc}}(y)=\frac{z(y)^{a}-1}{z(y)^{a}+1}\,,$$

(28)

and

$$z(y)=1+2y/a+2(y/a)^{2}$$

(29)

is the soft-core perturbation energy function to avoid singularities near the initial state of the alchemical transformation.⁴⁸ The parameters $u_{\rm max}$, $u_{c}$, and $a$ are set to cap the perturbation energy $u(x)$ to a maximum positive value without affecting it away from the singularity. The specific values of $u_{0}$, $u_{\rm max}$, and of the scaling parameter $a$ used in this work are listed in the Computational Details.

The standard linear alchemical interpolation potential

$$U_{\lambda}(x)=U_{0}(x)+\lambda u_{sc}[u(x)]$$

(30)

is a special case of Eq. (26) for $\lambda_{1}=\lambda_{2}=\lambda$. With one exception (see Calculation Details) all of the RBFE calculations reported in this work employed the linear alchemical potential.

Similar to the ABFE protocol,²⁷ the RBFE alchemical protocol involves two free energy legs. The initial state of the first leg corresponds to the physical state in which ligand 1 is bound to the receptor and ligand 2 is in the solvent bulk. The initial state of the second leg is the physical state in which the positions of the two ligands are swapped relative to the first leg. The initial states of the two legs are connected to the same alchemical intermediate defined by the potential function (25) at $\lambda=1/2$ corresponding to a non-physical state in which both ligands are present simultaneously in the binding site and in the solvent bulk, while not interacting with each other but each only interacting at 50% strength with the respective environments. The difference of the free energy changes, $\Delta G_{1}$ and $\Delta G_{2}$, of each leg gives the relative standard binding free energy of the two ligands

$$\Delta\Delta G^{\circ}_{b}(1,2)=\Delta G^{\circ}_{b}(2)-\Delta G^{\circ}_{b}(1)=\Delta G_{1}-\Delta G_{2}$$

(31)

The RBFE alchemical calculation for leg 1, for example, consists of a series of molecular dynamics simulations with the potential function (25) distributed at increasing values of $\lambda$ varying from $0$ to $1/2$ corresponding to ligand 1 bound and ligand 2 unbound, and to the alchemical intermediate, respectively. At each MD step, the perturbation energy (24) is obtained by displacing ligand 1 into the solvent bulk by translating its position along the vector $h$ while simultaneously translating ligand 2 in the opposite direction. The values of the perturbation energies are saved at regular intervals and processed by conventional multi-state reweighting analysis⁴¹ to compute the free energy change. The simulations for leg 2 are performed in the same way but starting from placing ligand 2 is in the binding site and ligand 1 is in the solvent bulk.

We have tested the ATM RBFE protocol on the SAMPL8 GDCC host-guest set⁵⁰¹¹1https://github.com/samplchallenges/SAMPL8/tree/master/host_guest/GDCC and on a set of complexes of the estrogen receptor $\alpha$ (ER$\alpha$).^{51, 21}

The SAMPL8 GDCC host-guest set consists of five small molecule guests binding to two host receptors (Figure 2). The hosts and the guests are known to form complexes with 1:1 stoichiometry under experimental conditions. The GDCC hosts (named TEMOA and TEETOA) are cup-shaped with a deep hydrophobic cavity occupied by the non-polar end of the bound guests. In the host-guest complexes, the charged/polar groups of the guests tend to be oriented towards the solvent. The hosts differ in the nature of the sidechains that surround the rim of the cavity (methyl groups in TEMOA and longer ethyl groups in TEETOA). In this work, all pairwise relative binding free energies for this set were compared to the differences of the ATM absolute binding free energy estimates that were obtained previously.⁵⁰ The ABFE estimates for this set are considered as a particularly trustworthy reference because they closely agree with independent potential of mean force estimates and with the experimentally measured binding free energies.⁵⁰ As shown in Figure 2, the majority of the pairwise RBFE calculations in this set involve scaffold hopping transformations. For example, all of the transformations to or from the G3 guest, among others, involve the deletion of a ring or the conversion of an aromatic 6-membered ring to an aliphatic 5-memberd ring. All of the guests except G2 are modeled in their anionic deprotonated forms. The G2 guest is modeled as neutral with the phenol oxygen protonated, as in our previous work, the protonated form of G2 demonstrated to contribute more to binding to both TEMOA and TEETOA hosts.⁵⁰ Hence, because all of the transformations involving G2 involve a change of the net charge of the ligand, this set tests the ability of the ATM RBFE protocol to handle changes in net charge in addition to changes of ligand topology.

The ER$\alpha$ benchmark consists of four molecular complexes of the ER$\alpha$ receptor with a series of benzopyran agonists.⁵¹ The four ligands (named 2e, 2d, 3a, and 3b) bind in a similar fashion to a large solvent-occluded cavity of the receptor (Figure 3). The ligands share a similar topology consisting of a common scaffold along the long molecular axis ending in two polar phenol groups, flanked by a variable side group with either a five- or six-membered ring. Alchemical RBFE transformations in this set can be cleanly classified as either straightforward terminal R-group modifications (such as 3a to 3b) or more challenging scaffold-hopping transformations involving ring expansion (such as 2e to 3a). The common long molecular axis and the conserved portion of the side group provide a natural definition of the reference frame to the four ligands.3

The host-guest systems were prepared from the provided MOL2 files provided by the SAMPL8 organizers (https://github.com/samplchallenges/SAMPL8/tree/master/host_guest/GDCC). The guests were manually docked to the hosts and aligned relative to each other using Maestro (Schrödinger Inc.). The structures of the ER$\alpha$ protein complexes were also prepared in Maestro starting from the provided files.²¹ Force-field parameter assignments and solvation of the molecular systems were performed using Ambertools 18. The GAFF1.8/AM1-BCC force field parameters^{52, 53} were assigned for the hosts and the ligands and the Amber ff14SB parameters^{54, 13} were assigned to the ER$\alpha$ protein receptor. The complexes were assembled using the tleap program in Ambertools 18. Prior to adding the solvent, the receptor and a pair of aligned ligands were loaded into tleap where one of the ligands was translated to a position corresponding to the solvent bulk (see below). Solvent molecules were added with a 10 Å buffer around the complex and the unbound ligand of the ligand pair. The resulting solvent box was neutralized by adding sodium or chloride ions as needed.

The complexes were energy minimized and thermalized at 300 K. Beginning at the bound state at $\lambda=0$, the systems were then annealed to the symmetric alchemical intermediate at $\lambda=1/2$ for $250$ ps using the ATM alchemical potential energy function for Leg 1 [Eq. (30)]. This step provides a suitable initial configuration of the system without severe unfavorable repulsion interactions at alchemical intermediate state to seed the molecular dynamics simulations. The host molecules were loosely restrained with a flat-bottom harmonic potential of force constant 25.0 kcal/(mol Ų) and a tolerance of 1.5 Å was set on the heavy atoms at the lower cup of the molecule (the first 40 atoms of the host as listed in the provided files). The same positional restraints were applied to the C$\alpha$ atoms of the ER$\alpha$ protein receptor.

With one exception (see below), the linear alchemical perturbation potential, $W_{\lambda}(u)=\lambda u_{\rm sc}(u)$ was used for all alchemical calculations with $11$ $\lambda$-states uniformly distributed between $\lambda=0$ and $1/2$ for each of the two ATM legs. The calculation for the more challenging test without alignment restraints from TEMOA-G4 to TEMOA-G1 employed the softplus perturbation potential (25) with the parameters listed in Table 1.

The soft-core perturbation energy Eq. (27) was used for all calculations with $u_{\rm max}=200$ kcal/mol, $u_{c}=100$ kcal/mol, and $a=1/16$. The ligands was sequestered within the binding site by means of a flat-bottom harmonic potential between the centers of mass of the receptor and the ligand with a force constant of $25$ kcal/mol Ų applied for separation greater than $4.5$ Å. Perturbation energy samples and trajectory frames were saved every 10 ps. Asynchronous Hamiltonian replica exchanges³⁹ in $\lambda$-space were performed every 10 ps. The Langevin thermostat with a time constant of 2 ps was used to maintain the temperature at 300 K. Each replica was simulated for a minimum of 20 ns for the host-guest systems and for 14 ns for the protein-ligand systems. Binding free energies and the corresponding uncertainties were computed from the perturbation energy samples using UWHAM⁴¹, discarding the first half of the trajectory.

Alignment restraints were imposed using three reference atoms of the ligands. For the SAMPL8 guests, alignment settings were customized for each pair to maintain an approximate alignment between the the polar groups oriented towards the solvent and the body of the guest oriented towards the binding cavity of the host. The values of the force constants for the restraining potentials in Eqs. (21) and (23) were adjusted to allow for more flexibility between more dissimilar guests. $k_{r}$ varied from $25$ to $0.5$ kcal/mol/Ų, $k_{\theta}$ varied from $50$ to $10$ kcal/mol, and $k_{\psi}$ varied from $50$ to $0$ kcal/mol. The reference atoms for the ER$\alpha$ ligands are illustrated in Figure 3 and the corresponding force constants were set to $k_{r}=2.5$ kcal/mol/Ų, $k_{\theta}=10$ kcal/mol, and $k_{\psi}=10$ kcal/mol. The specific alignment settings for each system can be found in the simulation input files posted at github.com/Gallicchio-Lab/ATM-relative-binding-free-energy-paper.

The alchemical molecular dynamics simulations were performed with the OpenMM 7.3⁴² MD engine and the SDM integrator plugin (github.com/Gallicchio-Lab/openmm_sdm_plugin.git) using the OpenCL platform. The ASyncRE software,³⁹ customized for OpenMM and SDM (github.com/Gallicchio-Lab/async_re-openmm.git), was used for the Hamiltonian Replica Exchange in $\lambda$ space for each ATM leg. Simulation input files for this work are freely available at github.com/Gallicchio-Lab/ATM-relative-binding-free-energy-paper. The replica exchange simulations were performed on the XSEDE Comet and Expanse GPU HPC cluster at the San Diego Supercomputing Center, each using four GPUs per node.

Table 3 and Figure 4 report the calculated relative binding free energies for the SAMPL8 host-guest systems compared to the corresponding estimates from the absolute binding free energies reported in reference 50 and reproduced here in Table 2. The results of the tests of the alignment restraints on the TEMOA-G1-G4 system are shown in Table 4.

As shown in Table 3, the calculated RBFEs for the twenty pairs of complexes (5th column), ten for each SAMPL8 host, are in good agreement with the corresponding estimates from the ABFEs values (4th column). The root mean square deviation between the two sets of free energy differences is only 0.8 kcal/mol, a value well within statistical uncertainties. The statistical uncertainties of the RBFEs estimates are found to be systematically smaller than those of the differences of ABFEs estimates. For example, the uncertainty for the large G1 to G3 scaffold-hopping transformation is 0.38 kcal/mol compared to an uncertainty of 0.66 kcal/mol from the differences of the corresponding ABFE estimates. The RBFE estimates are also found to have a high level of self-consistency with an average 3-points cycle closure error of only 0.25 kcal/mol. The latter value is computed by evaluating the free energy change along all of the triangular cycles in the network of transformations in Figure 2.

As illustrated in Table 3 the RBFE estimates are obtained from the differences between the free energy changes for the two legs ($\Delta G_{1}$ and $\Delta G_{2}$ in Table 3). These tend to be of moderate magnitude for the transformations between guests with the same net charge–between 10 to 20 kcal/mol range in magnitude, compared to the much larger free energy changes for the two legs of the ABFE calculations, which range between 40 to 50 kcal/mol.⁵⁰ This trend, which is remarkably maintained even for large scaffold-hopping transformations such as G1 to G3, can be justified by the large desolvation penalty occurring for binding a negatively-charged guest to the host that cancels out in the RBFE calculations. Conversely, the trend is reversed for the RBFE transformations involving the neutral protonated G2 guest (G2P in Tables 3 and 2), in which the ABFE leg free energies tend to be significantly smaller than for the RBFE transformations. Evidently, in these cases, the desolvation penalty is much smaller whereas the RBFE transformations are contending with a change of net charge. Nevertheless, even in these cases the RBFE calculations yield consistent free energy estimates with reasonably small statistical uncertainties.

As illustrated in Theory and Methods, the restraining potentials to keep the positions and orientations of the ligands in approximate alignment are designed to improve the convergence of the calculations without biasing the results. We tested this fact numerically on the RBFE calculation from G1 to G4 for the TEMOA host. The results, shown in Table 4, largely confirm the theoretical expectations. The RBFE estimate obtained without alignment restraints (that is with only binding site restraints) agrees within statistical uncertainty with the estimate with full alignment restraints and with only positional restraints (that is by turning off orientation restraints).

In addition to confirming the theory, this data shows that, despite the similarities between the G1 and G4 scaffolds (Figure 2, a strict alignment between two ligands is not necessarily beneficial for the RBFE calculation. For example, the estimate obtained without angular restraints shows slightly smaller free energy changes with smaller statistical uncertainties in the two legs compared to the full restraints. This suggests that allowing for more flexibility can lower the free energy of the intermediate state in which the two ligands occupy simultaneously the binding cavity.

The RBFE calculation without restraints, while also agreeing with the other estimates and with moderate uncertainty, suffered large entropic bottlenecks at the initial states of each leg (not shown) corresponding to the loss of translational and rotational degrees of freedom upon binding,⁴⁸ and required the use of the optimized softplus alchemical potential to reach convergence (see Computational Details). It can be thus concluded that some level of alignment restraints is beneficial to reach rapid convergence of the RBFE estimates with a general-purpose linear alchemical potential.

The calculated RBFEs for all six pairs of the four ER$\alpha$ complexes (Figure 3) are shown in Table 5 compared to the corresponding values obtained by Wang et al.²¹ Given that were obtained with different force field models, the relative binding free energies calculated here for these systems are in good agreement with the literature values, which were themselves in reasonable agreement with experimental observations.²¹ The values obtained here have similar uncertainties as the literature values and also display a good level of self-consistency with an average 3-point cycle closure error of $0.6$ kcal/mol, which is within statistical uncertainty.

One interesting and promising feature of the protein-ligand RBFEs obtained here is that the magnitudes of the free energy changes in each leg and their level of convergence do not appear to depend strongly on the difficulty of the alchemical transformation. For example, the large scaffold-hopping transformation from 2d to 3a which involves a ring expansion as well as the replacement and the displacement of a non-polar R-group with a polar one (Figure 3), presents leg free energies $\Delta G_{1}$ and $\Delta G_{2}$ which are only slightly larger in magnitude and with similar uncertainties that those of the 2d to 2e transformation which involves only a straightforward mutation of the three terminal atoms.