These authors contributed equally to this work Ph.D. Program in Biochemistry,
The Graduate Center of the City University of New York, New York, NY These
authors contributed equally to this work Ph.D. Program in Biochemistry, The
Graduate Center of the City University of New York, New York, NY These authors
contributed equally to this work Ph.D. Program in Chemistry, The Graduate
Center of the City University of New York, New York, NY Ph.D. Program in
Chemistry, The Graduate Center of the City University of New York, New York,
NY Ph.D. Program in Biochemistry, The Graduate Center of the City University
of New York, New York, NY

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

Solmaz Azimi  Department of Chemistry, Brooklyn College of the City University
of New York, New York, NY  Sheenam Khuttan  Department of Chemistry, Brooklyn
College of the City University of New York, New York, NY  Joe Z. Wu
Department of Chemistry, Brooklyn College of the City University of New York,
New York, NY  Rajat K. Pal  Roivant Sciences, Boston, MA  Emilio Gallicchio
<EMAIL_ADDRESS>Department of Chemistry, Brooklyn College of the
City University of New York, New York, NY

###### Abstract

We present an extension of Alchemical Transfer Method (ATM) for the estimation
of relative binding free energies of molecular complexes applicable to
conventional as well as scaffold-hopping alchemical transformations. The
method, named ATM-RBFE, implemented in the free and open-source OpenMM
molecular simulation package, aims to provide a simpler and more generally
applicable route to the calculation of relative binding free energies than is
currently available. The method is based on sound statistical mechanics theory
and a novel coordinate perturbation scheme designed to swap the positions of a
pair of ligands such that one is transferred from the bulk solvent to the
receptor binding site while the other moves simultaneously in the opposite
direction. The calculation is conducted directly using a single solvent box
prepared using conventional setup tools, without splitting of electrostatic
and non-electrostatic transformations, and without pairwise soft-core
potentials. ATM-RBFE is validated here against the absolute binding free
energies of the SAMPL8 GDCC host-guest benchmark set and against a benchmark
set of estrogen receptor $\alpha$ complexes. In each case the method yields
self-consistent and converged relative binding free energy estimates in
agreement with absolute binding free energies, reference literature values as
well as experimental measurements.

##  1 Introduction

Relative Binding Free Energy (RBFE) calculations are a critical component of
modern structure-based drug discovery efforts.1, 2, 3 Generally, the goal of
these calculations is to estimate the difference of the standard binding free
energies of two molecular ligands to the same receptor more directly than
taking the difference of the corresponding Absolute Binding Free Energies
(ABFE). While a variety of approaches to obtain RBFEs have been proposed, most
large-scale deployments to date employ a thermodynamic cycle in which one
ligand is alchemically transformed into another ligand while both are bound to
the receptor and, separately, while both are in the solvent bulk. The
difference of the free energy changes of these two transformations yields the
RBFE of the ligand pair.4, 5, 6

In RBFE methods, the transformations from one ligand into another are termed
alchemical as they are based on progressively modifying the Hamiltonian of the
system through a series of artificial chemical states that interpolate between
the physical end states corresponding to the two complexes. In RBFE
calculations, alchemical transformations are commonly carried out by parameter
interpolation, in which the parameters of the energy function (partial
charges, Lennard-Jones parameters, force constants, etc.) are progressively
changed by means of an alchemical progress parameter $\lambda$.7, 8 Generally,
the implementation of parameter interpolation methods require specific
customization of the energy routines of the molecular dynamics engine,
especially when soft-core pairwise potentials are used to treat end-point
singularities.9, 10

A variety of implementations of alchemical transformations are in common use.
In single-topology implementations, a selected set of atoms of the first
ligand are mapped to corresponding atoms of the second ligand. The first set
of atoms is then transformed into the other by interpolating the force field
parameters assigned to these atoms.11 In the case of two ligands with
different atom counts, the atoms of one ligand that do not correspond in the
other ligand are converted to and from dummy atoms.12 In dual-topology
implementations by the means of congeneric pairing, a whole functional group
of one ligand is typically converted to that of the second ligand by
converting the first group to a dummy group and converting the second from a
dummy group to the target group simultaneously.7

Single- and dual-topology approaches have distinct advantages and limitations
that make them more or less preferable depending on the specific application.
The correct treatment of dummy atoms in single-topology methods can be
particularly subtle.12 In dual-topology approaches, the mode of attachment of
the dummy group to the ligand scaffold can also require care in order to avoid
artifacts in the resulting free energies.13, 14 These and other implementation
issues can be quite complex to address especially when considered in the
context of hybrid approaches involving both single- and dual-topology
transformations15, as well as the more specialized separated-topology
approaches.16 In all of these methods, specialized system setup tools are
employed to deal with the non-chemical molecular topologies to treat, for
example, dummy atoms and provide ways to make sure that the atoms of the
groups being transformed in dual topologies do not interact with each other.
To avoid numerical instabilities, the alchemical transformations are also
typically split into distinct simulation legs to separately treat van der
Waals and electrostatic interactions.5, 8

The setup of RBFE calculations based on two simulation systems—the receptor
with the two ligands in the solvent and, separately, the ligands in the
solvent—is further complicated by the preparation and simulation of multiple
systems. RBFE implementations based on separate hydration and receptor-binding
systems, for example, are not naturally suited for ligand pairs with different
net charges.17 In these cases, the change of system charge during the
alchemical transformation causes issues with the treatment of long-range
electrostatic interactions and introduces artifacts in the free energy
estimates unless complex correction factors are introduced.18 To overcome some
of these challenges, coupled transformations have also been used in which the
change in the net charge of the ligands is counterbalanced by the opposite
change of charge of an ion to a solvent molecule or to another ion.19, 20

Scaffold hopping transformations that, unlike the simpler terminal R-group
modifications, include structural changes in the chemical topologies of the
two ligands such as breaking and forming rings, reducing the size of rings,
and modifying linker groups are particularly challenging with traditional RBFE
approaches. Wang. et al.21 introduced a customized soft-bond stretch potential
to address the numerical instabilities encountered when forming and breaking
covalent bonds. More recently, Zou et al.22 presented a method that uses a
series of restraining potentials to address the same problem.

As a result of the above discussed issues and other complexities associated
with RBFE methods23 and despite the strong need in structure-based drug
discovery, software tools in this arena remain of limited availability to many
practitioners in the field. This deficiency is particularly felt in academic
settings where commercial solutions might be prohibitively expensive or not
sufficiently customizable, and the substantial amount of human resources and
the specialized expertise required to develop and deploy RBFE solutions anew
are not available. As an example, there are only a handful of open-source
academic codes, at various stages of development, that are claimed to support
scaffold-hopping RBFE transformations.24, 25, 26

In an attempt to address some of the aforementioned challenges, in this work
we present a streamlined RBFE protocol based on the Alchemical Transfer Method
(ATM).27 The method aims to provide a simpler and more generally applicable
route to relative binding free energy calculations than what is currently
available. Like the parent ATM Absolute Binding Free Energy (ABFE) approach
and unlike conventional parameter interpolation-based methods, the ATM-RBFE
protocol is based on a straightforward coordinate transformation perturbation.
As described here in this method, and similar to the double-system/single-box
method28, 29 for ABFE calculations, the receptor and the two ligands are
initially placed in a single solvated simulation box in such a way that one
ligand is bound to the receptor and the other is placed in an arbitrary
position in the solvent bulk. Molecular dynamics simulations are then
conducted with a $\lambda$-dependent perturbation alchemical potential energy
function based on swapping the positions of the two ligands.

The ATM RBFE protocol proposed here has the following notable features.

  * •

It is applicable to simple terminal R-group transformations, as well as to
scaffold-hopping transformations to connect ligand pairs that do not share the
same topology.

  * •

It does not require splitting the alchemical transformations into
electrostatic and non-electrostatic steps, and it does not require the
implementation of soft-core pair potentials.

  * •

It does not require modifications of the core energy routines of the molecular
dynamics engine. It uses only the energies and forces returned by the
unmodified existing routines used for molecular dynamics time propagation.

  * •

By design, it is applicable with any potential energy function, including
conventional fixed-charge molecular mechanics force fields with and without
long-range electrostatics, as well as many-body potentials, such as
polarizable,30, 31 quantum-mechanical,32, 33, 34 machine learning
potentials,35, 36 implicit solvation models,37 and coarse-grained
potentials,38 which are generally poorly supported by free energy alchemical
protocols.

  * •

The molecular simulation system contains only chemically valid topologies
without dummy atoms and is prepared with conventional tools.

  * •

Because perturbation energies are $\lambda$-independent and can be evaluated
algebraically rather than by rescoring trajectories by calling the energy
routines of the MD engine, the method is easily implemented in conjunction
with advanced conformational sampling algorithms, such as Hamiltonian replica
exchange,39 and with multi-state free energy analysis tools.40, 41

  * •

It is amenable to simple, compact and self-contained software implementations.
The method presented here is implemented into a freely available and open-
source plugin of the popular OpenMM library for molecular simulations.42

The work is organized as follows. We first derive the formulation of the ABFE
and RBFE ATM protocols from a well established statistical mechanics theory of
non-covalent bimolecular binding. Here we also introduce ligand alignment
restraints used to enhance the convergence of the calculations for similar
ligand pairs and we show mathematically that they yield unbiased binding free
energy estimates. We then present the application of the ATM RBFE method to a
set of host-guest complexes and a set of protein-ligand complexes.
Specifically, we compute all pairwise relative binding free energies for the
SAMPL8 GDCC host-guest sets and for four proposed inhibitors of the estrogen
receptor $\alpha$ nuclear receptor. Despite the fact that a majority of the
ligand pairs considered in this work are related by challenging scaffold-
hopping transformations, we show that the ATM RBFE protocol yields self-
consistent relative binding free estimates that are in good agreement with not
only absolute binding free energy estimates, but also consistent with
literature reference values and with experimental measurements. This work
paves the way to streamlined and more readily available RBFE estimation tools
in structure-based drug discovery and other fields.

##  2 Theory and Methods

###  2.1 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$.

###  2.2 Absolute Binding Free Energy Formalism

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 representation45, 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).27

###  2.3 Relative Binding Free Energy Formalism

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).

Figure 1:  Diagrams illustrating the external degrees of freedom
$\zeta=({\mathbf{c}},{\mathbf{\omega}})$, with
${\mathbf{\omega}}=(\theta,\varphi,\psi)$, of a target molecule (red) relative
to a reference molecule (black). The reference frame of each molecule is
defined by three non co-linear reference atoms denoted by $1,2,3$ and
$1^{\prime},2^{\prime},3^{\prime}$ for the reference and target molecules,
respectively. $\mathbf{c}$ (top) is defined by the vector distance between
$1^{\prime}$ and $1$. The atoms $1,2,$ and $3$ define a Cartesian coordinate
system (bottom) with $1$ at the origin, and where the $z$ axis is along the
$1$ to $2$ direction and the $x$ axis is oriented such that $3$ is in the $xz$
plane. When the target molecule is translated by $-\mathbf{c}$ so that
$1^{\prime}$ coincides to $1$, $\theta$ and $\varphi$ are the spherical polar
angles of $2^{\prime}$. The angle $\psi$ is defined by the dihedral angle
formed by atoms $3$, $1$, $2^{\prime}$, and $3^{\prime}$ in this sequence. The
alignment of the two molecules is maintained by restraining potentials such
that $2^{\prime}$ is approximately co-linear with $2$ ($\theta\simeq 0$)–in
which case the $\varphi$ degree of freedom can be ignored–and $3^{\prime}$ is
approximately in the $xz$ plane ($\psi\simeq 0$).

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,7
we call this approach _dual-ligand_ to emphasize the fact that the
calculations is performed with a solvent box containing the two ligands.

###  2.4 Alignment Restraints

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.

###  2.5 Calculation Protocol

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,44 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.14 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 Hamiltonian47
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,49 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.48 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,27 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 analysis41 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.

###  2.6 Molecular Systems

We have tested the ATM RBFE protocol on the SAMPL8 GDCC host-guest
set50111https://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.50 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.50 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.50 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.

Figure 2: The SAMPL8 GDCC molecular systems

The ER$\alpha$ benchmark consists of four molecular complexes of the
ER$\alpha$ receptor with a series of benzopyran agonists.51 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

Figure 3: The ERalpha complexes

###  2.7 Computational Details

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.21
Force-field parameter assignments and solvation of the molecular systems were
performed using Ambertools 18. The GAFF1.8/AM1-BCC force field parameters52,
53 were assigned for the hosts and the ligands and the Amber ff14SB
parameters54, 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 Å2) 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.

Table 1:  Alchemical schedule of the softplus perturbation function for the
two legs of the TEMOA-G1 to TEMOA-G4 transformation.

$\lambda$ | $\lambda_{1}$ | $\lambda_{2}$ |  $\alpha$a |  $u_{0}$b |  $w_{0}$b
---|---|---|---|---|---
0.00 | 0.00 | 0.00 | 0.10 | 150 | 0
0.05 | 0.00 | 0.05 | 0.10 | 135 | 0
0.10 | 0.00 | 0.10 | 0.10 | 120 | 0
0.15 | 0.00 | 0.15 | 0.10 | 105 | 0
0.20 | 0.00 | 0.20 | 0.10 | 90 | 0
0.25 | 0.00 | 0.25 | 0.10 | 75 | 0
0.30 | 0.10 | 0.30 | 0.10 | 60 | 0
0.35 | 0.20 | 0.35 | 0.10 | 40 | 0
0.40 | 0.30 | 0.40 | 0.10 | 40 | 0
0.45 | 0.40 | 0.45 | 0.10 | 40 | 0
0.50 | 0.50 | 0.50 | 0.10 | 40 | 0

aIn (kcal/mol)-1. bIn kcal/mol.

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 Å2 applied for separation greater than $4.5$
Å. Perturbation energy samples and trajectory frames were saved every 10 ps.
Asynchronous Hamiltonian replica exchanges39 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 UWHAM41, 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/Å2,
$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/Å2, $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.342 MD engine and the SDM integrator plugin (github.com/Gallicchio-
Lab/openmm_sdm_plugin.git) using the OpenCL platform. The ASyncRE software,39
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.

##  3 Results

###  3.1 SAMPL8 Host-Guest Systems

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.50 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.

Table 2: ATM absolute binding free energy estimates for the host-guest complexes from reference 50. Complex | $\Delta G^{\circ}_{b}$
---|---
TEMOA-G1 | $-6.71\pm 0.30$
TEMOA-G2P | $-13.10\pm 0.83$
TEMOA-G3 | $-8.26\pm 0.25$
TEMOA-G4 | $-8.63\pm 0.30$
TEMOA-G5 | $-7.70\pm 0.30$
TEETOA-G1 | $-1.07\pm 0.34$
TEETOA-G2P | $-7.95\pm 0.28$
TEETOA-G3 | $-1.65\pm 0.30$
TEETOA-G4 | $-2.51\pm 0.30$
TEETOA-G5 | $-2.82\pm 0.28$
Table 3: Relative binding free energies of the TEMOA and TEETOA host-guest
complexes.

Complex Pair |  $\Delta G_{1}$a,b |  $\Delta G_{2}$a,b |  $\Delta\Delta G^{\circ}_{\rm b}$(RBFE)a,b |  $\Delta\Delta G^{\circ}_{b}$(ABFE)a,c
---|---|---|---|---
TEMOA-G1-G2p | $41.68\pm 0.21$ | $47.34\pm 0.30$ | $-5.66\pm 0.37$ | $-6.39\pm 0.88$
TEMOA-G1-G3 | $20.67\pm 0.18$ | $23.01\pm 0.24$ | $-2.34\pm 0.30$ | $-1.55\pm 0.39$
TEMOA-G1-G4 | $15.98\pm 0.30$ | $18.45\pm 0.39$ | $-2.47\pm 0.49$ | $-1.92\pm 0.42$
TEMOA-G1-G5 | $14.53\pm 0.24$ | $16.39\pm 0.39$ | $-1.86\pm 0.46$ | $-0.99\pm 0.42$
TEMOA-G2p-G3 | $49.52\pm 0.30$ | $45.59\pm 0.36$ | $3.93\pm 0.47$ | $4.84\pm 0.87$
TEMOA-G2p-G4 | $45.19\pm 0.27$ | $42.05\pm 0.27$ | $3.14\pm 0.38$ | $4.47\pm 0.88$
TEMOA-G2p-G5 | $44.07\pm 0.24$ | $40.22\pm 0.27$ | $3.85\pm 0.36$ | $5.40\pm 0.88$
TEMOA-G3-G4 | $16.62\pm 0.39$ | $16.95\pm 0.39$ | $-0.33\pm 0.55$ | $-0.37\pm 0.39$
TEMOA-G3-G5 | $16.82\pm 0.48$ | $16.80\pm 0.45$ | $0.02\pm 0.66$ | $0.56\pm 0.39$
TEMOA-G4-G5 | $10.80\pm 0.27$ | $10.35\pm 0.24$ | $0.44\pm 0.36$ | $0.93\pm 0.42$
TEETOA-G1-G2p | $40.97\pm 0.27$ | $48.52\pm 0.27$ | $-7.55\pm 0.38$ | $-6.88\pm 0.67$
TEETOA-G1-G3 | $21.64\pm 0.24$ | $23.24\pm 0.27$ | $-1.60\pm 0.36$ | $-0.58\pm 0.66$
TEETOA-G1-G4 | $16.20\pm 0.30$ | $17.70\pm 0.39$ | $-1.50\pm 0.49$ | $-1.44\pm 0.66$
TEETOA-G1-G5 | $14.51\pm 0.24$ | $15.70\pm 0.39$ | $-1.19\pm 0.46$ | $-1.75\pm 0.65$
TEETOA-G2p-G3 | $51.89\pm 0.27$ | $44.94\pm 0.27$ | $6.95\pm 0.38$ | $6.30\pm 0.61$
TEETOA-G2p-G4 | $46.33\pm 0.30$ | $40.95\pm 0.27$ | $5.38\pm 0.40$ | $5.44\pm 0.61$
TEETOA-G2p-G5 | $44.46\pm 0.24$ | $39.08\pm 0.27$ | $5.38\pm 0.36$ | $5.13\pm 0.60$
TEETOA-G3-G4 | $17.90\pm 0.30$ | $18.82\pm 0.24$ | $-0.92\pm 0.38$ | $-0.86\pm 0.62$
TEETOA-G3-G5 | $16.69\pm 0.45$ | $17.09\pm 0.48$ | $-0.40\pm 0.66$ | $-1.17\pm 0.61$
TEETOA-G4-G5 | $10.43\pm 0.30$ | $10.47\pm 0.24$ | $-0.04\pm 0.38$ | $-0.31\pm 0.61$

a In kcal/mol, errors are reported as 3 times the standard deviation. b This
work. c from the differences of absolute binding free energies from Table 2.50

Figure 4: Validation of RBFE against the ABFE values for the SAMPL8 benchmark
set. The line represents perfect agreement between the RBFE values and the
differences of the corresponding ABFE results with the root mean square
deviation of 0.8 kcal/mol, a value well within statistical uncertainties .

###  3.2 Test of the Alignment Restraints

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,48 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.

Table 4: Calculated binding free energy of TEMOA-G4 relative to TEMOA-G1 with
various alignment restraints

Protocol | $\Delta G_{1}$ | $\Delta G2$ |  $\Delta\Delta G^{\circ}_{b}$(RBFE)
---|---|---|---
ABFE differencea |  |  | $-1.92\pm 0.42$
RBFE full restraintsa | $15.98\pm 0.30$ | $18.45\pm 0.39$ | $-2.47\pm 0.48$
RBFE pos. restraints | $15.26\pm 0.20$ | $17.54\pm 0.27$ | $-2.28\pm 0.34$
RBFE no restraints | $20.94\pm 0.26$ | $23.39\pm 0.30$ | $-2.45\pm 0.40$

a from Table 3.

###  3.3 ER$\alpha$ Complexes

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.21 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.21 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.

Table 5: Relative binding free energies of the ER$\alpha$ complexes.

Complex Pair |  $\Delta G_{1}$a,b |  $\Delta G_{2}$a,b |  $\Delta\Delta G^{\circ}_{\rm b}$(RBFE)a,b |  $\Delta\Delta G^{\circ}_{b}$(Wang et al.)a,c
---|---|---|---|---
2d-2e | $12.67\pm 0.28$ | $15.74\pm 0.39$ | $-3.06\pm 0.48$ | $-1.20\pm 0.12$
2d-3a | $17.01\pm 0.39$ | $14.83\pm 0.43$ | $2.18\pm 0.58$ | $3.77\pm 0.30$
2d-3b | $16.17\pm 0.37$ | $16.30\pm 0.40$ | $-0.13\pm 0.54$ | $1.37\pm 0.36$
2e-3a | $17.66\pm 0.33$ | $13.35\pm 0.34$ | $4.31\pm 0.48$ | $5.26\pm 0.30$
2e-3b | $13.87\pm 0.27$ | $11.33\pm 0.35$ | $2.53\pm 0.45$ | $2.86\pm 0.33$
3a-3b | $15.10\pm 0.26$ | $17.70\pm 0.32$ | $-2.60\pm 0.41$ | $-2.36\pm 0.33$

a In kcal/mol, errors are reported as 3 times the standard deviation. b This
work. c from reference 21.

##  4 Discussion

Binding free energy estimation remains one of the most challenging problem in
modern computational chemistry. This is particularly so for the
conformationally variable macromolecular targets relevant to structure-based
drug discovery. Relative binding free energy (RBFE) methods based on
alchemical transformations address some of these challenges by comparing
directly pairs of ligands that share a similar binding mode thereby bypassing,
for example, the modeling of the transition from the apo to the holo state of
the receptor necessary for absolute binding free energy estimation (ABFE).
However, alchemical RBFE methods based on conventional single- and dual-
topology setups have been traditionally limited to simple R-group
modifications, which considerably limits their range of applicability.

Recently, advanced alchemical scaffold-hopping methods have been introduced
which allow for more complex changes in molecular topologies involving
breaking and formation of covalent bonds. Wang et al.21 introduced a
customized soft bond stretch potential to smoothly form and break bonds. More
recently, Zou et al.22 introduced a scheme based on auxiliary restraining
potentials to achieve the same goal. Multi-site $\lambda$-dynamics has been
presented as an alternative to scaffold hopping transformations involving ring
opening and ring expansion/reduction.55 These advances have required a
considerable amount of development and testing effort and are currently only
available on commercial platforms–the Schrödinger FEP+, the XTalPi XFEP, the
BIOVIA Discovery Studio platforms respectively. As academic efforts have
lagged behind in this arena,24, 25, 26 open-source software packages have not
yet reached a comparable level of robustness for production applications.

Scaffold-hopping transformations are only one of the many challenges that
prevent a widespread implementation and use of alchemical RBFE tools. Among
others, single- and dual-topology RBFE implementations generally require the
customization of energy and force routines, the design of pairwise soft-core
potentials, and the treatment of dummy groups and their attachment to the core
of the ligand. In addition, RBFE calculations generally require the mapping of
the corresponding atoms of the ligands, the setup of separate simulations for
the bound and unbound states, and the treatment of changes in net charge of
the ligands.7, 8 The substantial intellectual and development efforts involved
in the correct design, testing, and applications of alchemical RBFE software
tools can be discouraging to many practitioners.

In this work, we have presented a streamlined RBFE protocol based on the
Alchemical Transfer Method (ATM)27 that promises to remove many of these
obstacles while retaining a similar level of computational efficiency as the
most advanced RBFE methods. The proposed ATM-RBFE protocol makes use of a
single simulation system and does not require soft-core pair potentials and
non-physical chemical topologies with dummy atoms. The method also does not
require modifications of energy routines and can be easily implemented as a
controlling routine on top of existing force routines of MD engines. The
method has been implemented as open-source and freely available plugin of the
OpenMM molecular simulation library.42

The method has its roots in the development of a statistical mechanics
formulation of alchemical binding56 and the discovery of coordinate-based free
energy perturbation schemes and novel alchemical potential functions capable
of removing sampling bottlenecks along the alchemical thermodynamic path.48
The resulting ATM method has been recently extended to the modeling of
hydration49 and molecular binding27, 50 with explicit solvation. In the
present work, we combine ATM with a ligand alignment restraints protocol to
model relative binding free energies. We tested the method on two stringent
benchmarks that include scaffold-hopping and net charge modifications
obtaining good agreement with ABFE calculations, literature values, and
experimental observations.

##  5 Conclusions

We presented the theory and the implementation of the ATM-RBFE protocol for
the streamlined calculation of relative binding free energies. We successfully
tested the method on the SAMPL8 GDCC host-guest set and on a benchmark set of
estrogen receptor $\alpha$ complexes including challenging scaffold-hopping
transformations.

##  6 Acknowledgements

We acknowledge support from the National Science Foundation (NSF CAREER
1750511 to E.G.). Molecular simulations were conducted on the Comet and
Expanse GPU clusters at the San Diego Supercomputing Center supported by NSF
XSEDE award TG-MCB150001.

## References

  * Jorgensen 2004 Jorgensen, W. L. The many roles of computation in drug discovery. _Science_ 2004, _303_ , 1813–1818
  * Abel et al. 2017 Abel, R.; Wang, L.; Harder, E. D.; Berne, B.; Friesner, R. A. Advancing drug discovery through enhanced free energy calculations. _Acc. Chem. Res._ 2017, _50_ , 1625–1632
  * Armacost et al. 2020 Armacost, K. A.; Riniker, S.; Cournia, Z. Novel directions in free energy methods and applications. 2020
  * Wang et al. 2015 Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Mass, C.; Knight, L. J.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcho, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field. _J. Am. Chem. Soc._ 2015, _137_ , 2695–2703
  * Cournia et al. 2017 Cournia, Z.; Allen, B.; Sherman, W. Relative binding free energy calculations in drug discovery: recent advances and practical considerations. _J. Chem. Inf. Model._ 2017, _57_ , 2911–2937
  * Cournia et al. 2020 Cournia, Z.; Allen, B. K.; Beuming, T.; Pearlman, D. A.; Radak, B. K.; Sherman, W. Rigorous free energy simulations in virtual screening. _Journal of Chemical Information and Modeling_ 2020,
  * Mey et al. 2020 Mey, A. S. J. S.; Allen, B. K.; Macdonald, H. E. B.; Chodera, J. D.; Hahn, D. F.; Kuhn, M.; Michel, J.; Mobley, D. L.; Naden, L. N.; Prasad, S.; Rizzi, A.; Scheen, J.; Shirts, M. R.; Tresadern, G.; Xu, H. Best Practices for Alchemical Free Energy Calculations [Article v1.0]. _Living Journal of Computational Molecular Science_ 2020, _2_ , 18378
  * Lee et al. 2020 Lee, T.-S.; Allen, B. K.; Giese, T. J.; Guo, Z.; Li, P.; Lin, C.; McGee Jr, T. D.; Pearlman, D. A.; Radak, B. K.; Tao, Y.; Tsai, H.-C.; Xu, H.; Sherman, W.; York, D. M. Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery. _J. Chem. Inf. Model._ 2020,
  * Steinbrecher et al. 2011 Steinbrecher, T.; Joung, I.; Case, D. A. Soft-core potentials in thermodynamic integration: Comparing one- and two-step transformations. _J. Comput. Chem._ 2011, _32_ , 3253–3263
  * Lee et al. 2020 Lee, T.-S.; Lin, Z.; Allen, B. K.; Lin, C.; Radak, B. K.; Tao, Y.; Tsai, H.-C.; Sherman, W.; York, D. M. Improved Alchemical Free Energy Calculations with Optimized Smoothstep Softcore Potentials. _J. Chem. Theory Comput._ 2020,
  * Liu et al. 2013 Liu, S.; Wu, Y.; Lin, T.; Abel, R.; Redmann, J. P.; Summa, C. M.; Jaber, V. R.; Lim, N. M.; Mobley, D. L. Lead optimization mapper: automating free energy calculations for lead optimization. _J. Comp. Aided Mol. Des._ 2013, _27_ , 755–770
  * Fleck et al. 2021 Fleck, M.; Wieder, M.; Boresch, S. Dummy Atoms in Alchemical Free Energy Calculations. _J. Chem. Theory Comput._ 2021, _XXX_ , XXX–XXX
  * Zou et al. 2019 Zou, J.; Tian, C.; Simmerling, C. Blinded prediction of protein–ligand binding affinity using Amber thermodynamic integration for the 2018 D3R grand challenge 4. _J. Comp. Aid. Mol. Des._ 2019, _33_ , 1021–1029
  * Gallicchio 2021 Gallicchio, E. In _Computational Peptide Science: Methods and Protocols_ ; Simonson, T., Ed.; Methods in Molecular Biology; Springer Nature, 2021
  * Jiang et al. 2019 Jiang, W.; Chipot, C.; Roux, B. Computing relative binding affinity of ligands to receptor: An effective hybrid single-dual-topology free-energy perturbation approach in NAMD. _J. Chem. Inf. Model._ 2019, _59_ , 3794–3802
  * Rocklin et al. 2013 Rocklin, G. J.; Mobley, D. L.; Dill, K. A. Separated topologies—A method for relative binding free energy calculations using orientational restraints. _The Journal of chemical physics_ 2013, _138_ , 02B614
  * Chen et al. 2018 Chen, W.; Deng, Y.; Russell, E.; Wu, Y.; Abel, R.; Wang, L. Accurate calculation of relative binding free energies between ligands with different net charges. _J. Chem. Theory Comput._ 2018, _14_ , 6346–6358
  * Rocklin et al. 2013 Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects. _The Journal of chemical physics_ 2013, _139_ , 11B606_1
  * Dixit and Chipot 2001 Dixit, S. B.; Chipot, C. Can absolute free energies of association be estimated from molecular mechanical simulations? The biotin-streptavidin system revisited. _J. Phys. Chem. A_ 2001, _105_ , 9795–9799
  * Chen et al. 2013 Chen, W.; Wallace, J. A.; Yue, Z.; Shen, J. K. Introducing titratable water to all-atom molecular dynamics at constant pH. _Biophys. J._ 2013, _105_ , L15–L17
  * Wang et al. 2017 Wang, L.; Deng, Y.; Wu, Y.; Kim, B.; LeBard, D. N.; Wandschneider, D.; Beachy, M.; Friesner, R. A.; Abel, R. Accurate modeling of scaffold hopping transformations in drug discovery. _J. Chem. Theory Comput._ 2017, _13_ , 42–54
  * Zou et al. 2021 Zou, J.; Li, Z.; Liu, S.; Peng, C.; Fang, D.; Wan, X.; Lin, Z.; Lee, T.-S.; Raleigh, D. P.; Yang, M.; Simmerling, C. Scaffold Hopping Transformations Using Auxiliary Restraints for Calculating Accurate Relative Binding Free Energies. _J. Chem. Theory Comput._ 2021,
  * Loeffler et al. 2018 Loeffler, H. H.; Bosisio, S.; Duarte Ramos Matos, G.; Suh, D.; Roux, B.; Mobley, D. L.; Michel, J. Reproducibility of free energy calculations across different molecular simulation software packages. _J. Chem. Theory Comput._ 2018, _14_ , 5567–5582
  * Jespers et al. 2019 Jespers, W.; Esguerra, M.; Åqvist, J.; Gutiérrez-de Terán, H. QligFEP: an automated workflow for small molecule free energy calculations in Q. _J. Cheminf._ 2019, _11_ , 26
  * Vilseck et al. 2019 Vilseck, J. Z.; Sohail, N.; Hayes, R. L.; Brooks III, C. L. Overcoming challenging substituent perturbations with multisite $\lambda$-dynamics: a case study targeting $\beta$-secretase 1. _J. Phys. Chem. Lett._ 2019, _10_ , 4875–4880
  * Ries et al. 2021 Ries, B.; Normak, K.; Weiss, R. G.; Rieder, S.; Champion, C.; König, G.; Riniker, S. Relative binding free energies with scaffold-hopping type transformations using RE-EDS. _Submitted for publication_ 2021, _XXX_ , XXX
  * Wu et al. 2021 Wu, J. Z.; Azimi, S.; Khuttan, S.; Deng, N.; Gallicchio, E. Alchemical Transfer Approach to Absolute Binding Free Energy Estimation. _J. Chem. Theory Comput._ 2021, _17_ , 3309
  * Gapsys et al. 2015 Gapsys, V.; Michielssens, S.; Peters, J. H.; de Groot, B. L.; Leonov, H. _Molecular Modeling of Proteins_ ; Springer, 2015; pp 173–209
  * Macchiagodena et al. 2020 Macchiagodena, M.; Pagliai, M.; Karrenbrock, M.; Guarnieri, G.; Iannone, F.; Procacci, P. Virtual Double-System Single-Box: A Nonequilibrium Alchemical Technique for Absolute Binding Free Energy Calculations: Application to Ligands of the SARS-CoV-2 Main Protease. _J. Chem. Theory Comput._ 2020, _16_ , 7160–7172
  * Harger et al. 2017 Harger, M.; Li, D.; Wang, Z.; Dalby, K.; Lagardère, L.; Piquemal, J.-P.; Ponder, J.; Ren, P. Tinker-OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs. _J. Comp. Chem._ 2017, _38_ , 2047–2055
  * Panel et al. 2018 Panel, N.; Villa, F.; Fuentes, E. J.; Simonson, T. Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations. _Biophys. J._ 2018, _114_ , 1091–1102
  * Beierlein et al. 2011 Beierlein, F. R.; Michel, J.; Essex, J. W. A simple QM/MM approach for capturing polarization effects in protein- ligand binding free energy calculations. _J. Phys. Chem. B_ 2011, _115_ , 4911–4926
  * Lodola and De Vivo 2012 Lodola, A.; De Vivo, M. _Adv. Prot. Chem. Struct. Biol._ ; Elsevier, 2012; Vol. 87; pp 337–362
  * Hudson et al. 2019 Hudson, P. S.; Woodcock, H. L.; Boresch, S. Use of interaction energies in QM/MM free energy simulations. _J. Chem. Theory Comput._ 2019, _15_ , 4632–4645
  * Smith et al. 2019 Smith, J. S.; Nebgen, B. T.; Zubatyuk, R.; Lubbers, N.; Devereux, C.; Barros, K.; Tretiak, S.; Isayev, O.; Roitberg, A. E. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. _Nature Comm._ 2019, _10_ , 1–8
  * Rufa et al. 2020 Rufa, D. A.; Macdonald, H. E. B.; Fass, J.; Wieder, M.; Grinaway, P. B.; Roitberg, A. E.; Isayev, O.; Chodera, J. D. Towards chemical accuracy for alchemical free energy calculations with hybrid physics-based machine learning/molecular mechanics potentials. _BioRxiv_ 2020,
  * Zhang et al. 2017 Zhang, B.; Kilburg, D.; Eastman, P.; Pande, V. S.; Gallicchio, E. Efficient Gaussian Density Formulation of Volume and Surface Areas of Macromolecules on Graphical Processing Units. _J. Comp. Chem._ 2017, _38_ , 740–752
  * Spiriti et al. 2019 Spiriti, J.; Subramanian, S. R.; Palli, R.; Wu, M.; Zuckerman, D. M. Middle-way flexible docking: Pose prediction using mixed-resolution Monte Carlo in estrogen receptor $\alpha$. _PloS one_ 2019, _14_ , e0215694
  * Gallicchio et al. 2015 Gallicchio, E.; Xia, J.; Flynn, W. F.; Zhang, B.; Samlalsingh, S.; Mentes, A.; Levy, R. M. Asynchronous replica exchange software for grid and heterogeneous computing. _Computer Physics Communications_ 2015, _196_ , 236–246
  * Shirts and Chodera 2008 Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. _J. Chem. Phys._ 2008, _129_ , 124105
  * Tan et al. 2012 Tan, Z.; Gallicchio, E.; Lapelosa, M.; Levy, R. M. Theory of binless multi-state free energy estimation with applications to protein-ligand binding. _J. Chem. Phys._ 2012, _136_ , 144102
  * Eastman et al. 2017 Eastman, P.; Swails, J.; Chodera, J. D.; McGibbon, R. T.; Zhao, Y.; Beauchamp, K. A.; Wang, L.-P.; Simmonett, A. C.; Harrigan, M. P.; Stern, C. D., et al. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. _PLoS Comp. Bio._ 2017, _13_ , e1005659
  * Gilson et al. 1997 Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The Statistical-Thermodynamic Basis for Computation of Binding Affinities: A Critical Review. _Biophys. J._ 1997, _72_ , 1047–1069
  * Gallicchio and Levy 2011 Gallicchio, E.; Levy, R. M. Recent Theoretical and Computational Advances for Modeling Protein-Ligand Binding Affinities. _Adv. Prot. Chem. Struct. Biol._ 2011, _85_ , 27–80
  * Roux and Simonson 1999 Roux, B.; Simonson, T. Implicit Solvent Models. _Biophys. Chem._ 1999, _78_ , 1–20
  * Boresch et al. 2003 Boresch, S.; Tettinger, F.; Leitgeb, M.; Karplus, M. Absolute binding free energies: A quantitative approach for their calculation. _J. Phys. Chem. B_ 2003, _107_ , 9535–9551
  * Chipot and Pohorille (2007) Eds. Chipot and Pohorille (Eds.), _Free Energy Calculations. Theory and Applications in Chemistry and Biology_ ; Springer Series in Chemical Physics; Springer, Berlin Heidelberg: Berlin Heidelberg, 2007
  * Pal and Gallicchio 2019 Pal, R. K.; Gallicchio, E. Perturbation potentials to overcome order/disorder transitions in alchemical binding free energy calculations. _J. Chem. Phys._ 2019, _151_ , 124116
  * Khuttan et al. 2021 Khuttan, S.; Azimi, S.; Wu, J. Z.; Gallicchio, E. Alchemical Transformations for Concerted Hydration Free Energy Estimation with Explicit Solvation. _J. Chem. Phys_ 2021, _154_ , 054103
  * Azimi et al. 2021 Azimi, S.; Khuttan, S.; Wu, J. Z.; Deng, N.; Gallicchio, E. Application of the Alchemical Transfer and Potential of Mean Force Methods to the SAMPL8 Cavitand Host-Guest Blinded Challenge. _J. Comp.-Aid. Mol. Des._ 2021, _XXX_ , XXX–XXX
  * Norman et al. 2006 Norman, B. H.; Dodge, J. A.; Richardson, T. I.; Borromeo, P. S.; Lugar, C. W.; Jones, S. A.; Chen, K.; Wang, Y.; Durst, G. L.; Barr, R. J.; Montrose-Rafizadeh, C.; Osborne, H. E.; Amos, R. M.; Guo, S.; Boodhoo, A.; Krishnan, V. Benzopyrans are selective estrogen receptor $\beta$ agonists with novel activity in models of benign prostatic hyperplasia. _J. Med. Chem._ 2006, _49_ , 6155–6157
  * Wang et al. 2006 Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. _J. Mol. Graph. Model._ 2006, _25_ , 247–260
  * He et al. 2020 He, X.; Liu, S.; Lee, T.-S.; Ji, B.; Man, V. H.; York, D. M.; Wang, J. Fast, Accurate, and Reliable Protocols for Routine Calculations of Protein–Ligand Binding Affinities in Drug Design Projects Using AMBER GPU-TI with ff14SB/GAFF. _ACS Omega_ 2020, _5_ , 4611–4619
  * Maier et al. 2015 Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. _J. Chem. Theory Comput._ 2015, _11_ , 3696–3713
  * Raman et al. 2020 Raman, E. P.; Paul, T. J.; Hayes, R. L.; Brooks III, C. L. Automated, accurate, and scalable relative protein–ligand binding free-energy calculations using lambda dynamics. _Journal of Chemical Theory and Computation_ 2020, _16_ , 7895–7914
  * Kilburg and Gallicchio 2018 Kilburg, D.; Gallicchio, E. Analytical Model of the Free Energy of Alchemical Molecular Binding. _J. Chem. Theory Comput._ 2018, _14_ , 6183–6196