Atommap Corporation, New York, NY 10017 Atommap Corporation, New York, NY
10017 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

# Performance and Analysis of the Alchemical Transfer Method for Binding Free
Energy Predictions of Diverse Ligands

Lieyang Chen  Roivant Sciences, 151 W 42nd Street, 15th Floor, New York, NY
10036, United States  Yujie Wu  Roivant Sciences, 151 W 42nd Street, 15th
Floor, New York, NY 10036, United States  Chuanjie Wu  Roivant Sciences, 151 W
42nd Street, 15th Floor, New York, NY 10036, United States  Ana Silveira
Psivant Therapeutics, 451 D Street, Boston, MA 02210, United States  Woody
Sherman  Psivant Therapeutics, 451 D Street, Boston, MA 02210, United States
Huafeng Xu  Roivant Sciences, 151 W 42nd Street, 15th Floor, New York, NY
10036, United States  Emilio Gallicchio<EMAIL_ADDRESS>Department of Chemistry, Brooklyn College of the City University of New York,
New York, NY

###### Abstract

The Alchemical Transfer Method (ATM) is herein validated against the relative
binding free energies of a diverse set of protein-ligand complexes. We
employed a streamlined setup workflow, a bespoke force field, and the AToM-
OpenMM software to compute the relative binding free energies (RBFE) of the
benchmark set prepared by Schindler and collaborators at Merck KGaA. This
benchmark set includes examples of standard small R-group ligand modifications
as well as more challenging scenarios, such as large R-group changes, scaffold
hopping, formal charge changes, and charge-shifting transformations. The novel
coordinate perturbation scheme and a dual-topology approach of ATM address
some of the challenges of single-topology alchemical relative binding free
energy methods. Specifically, ATM eliminates the need for splitting
electrostatic and Lennard-Jones interactions, atom mapping, defining ligand
regions, and post-corrections for charge-changing perturbations. Thus, ATM is
simpler and more broadly applicable than conventional alchemical methods,
especially for scaffold-hopping and charge-changing transformations. Here, we
performed well over 500 relative binding free energy calculations for eight
protein targets and found that ATM achieves accuracy comparable to existing
state-of-the-art methods, albeit with larger statistical fluctuations. We
discuss insights into specific strengths and weaknesses of the ATM method that
will inform future deployments. This study confirms that ATM is applicable as
a production tool for relative binding free energy (RBFE) predictions across a
wide range of perturbation types within a unified, open-source framework.

##  1 Introduction

Alchemical binding free energy prediction tools are emerging as a best-in-
class standard for in silico prediction of binding free energies (i.e.,
protein-ligand binding affinity) in structure-based drug design.1, 2, 3, 4, 5,
6 While many challenges remain,7, 8 the increased utilization of RBFE
calculations has been fueled in part by the promising results of large-scale
validation campaigns against benchmark sets representative of actual drug
discovery projects.9, 10, 11, 12, 8, 13, 14, 15, 16, 17, 18, 19

Despite decades of progress, reliable prediction of binding affinities of
protein-ligand complexes remains a challenging problem with many unresolved
issues, especially related to large chemical modifications, core
transformations, and formal charge changes. In principle, the dissociation
constant of a complex can be measured via brute force molecular dynamics (MD)
simulations by sampling many binding/unbinding events.20 However, this
approach is generally not practical because the typical residence time of the
ligand in the binding site (from milliseconds to hours) is too computationally
expensive for practical applications in drug discovery. The alternative
physical pathway methods21, 22, 23, 24, 25, 26, 27, 28 consist of measuring
the reversible work of dragging the ligand from the solution to the binding
site (or vice versa) along a chosen route. While physically appealing and
computationally more efficient than brute force MD,29, 26, 30, 31 physical
pathway methods are rarely used in small molecule structure-based drug design
because they require overcoming transition states and they do not readily
apply to the common enclosed binding sites that lack a clear entryway.27
Furthermore, physical pathway methods do not apply to the direct estimation of
relative binding free energies useful in drug discovery applications.

In drug discovery applications, knowledge of the change in binding affinity
resulting from modifying a ligand into another is usually more relevant than
that of their absolute binding affinities due to the nature of the iterative
design-make-test (DMT) cycle that is almost always employed to advance a hit
to a development candidate.10, 2, 12 Furthermore, while sampling of the full
binding/unbinding pathway can yield information about binding kinetics, it is
unnecessary for the computation of binding thermodynamics. Accordingly,
Relative Binding Free Energy (RBFE) alchemical protocols can estimate the
difference in the binding free energies of a pair of related ligands more
directly than computing each of their absolute binding free energies.32, 19,
33, 34 That being said, absolute binding free energies still have great
potential value in the context of virtual screening of diverse molecules.35
This work focuses on relative binding free energies in the context of
optimizing screening hits to drugs.

In most RBFE software implementations, the transformation of one molecule to
another is accomplished by parameter interpolation approaches36, 33 that scale
parameters of the potential energy function to convert one ligand into another
through an alchemical transformation. However, parameter interpolation schemes
require complex and often non-transferable customization of the energy
subroutines of molecular dynamics engines and custom soft-core pair-potentials
used to correct singularities of the alchemical potential energy function.37,
38 Moreover, parameter interpolation implementations typically do not directly
connect the two ligands in their bound states. Rather, they rely on a
thermodynamic cycle and alchemical calculations in solution and receptor
environments separately, and often each step is further split into the
decoupling of electrostatic and non-electrostatic interactions to avoid
numerical instabilities.33, 8 Furthermore, transformations involving changes
in the net ligand charge require the additional calculation of correction
factors.39, 40, 41, 42, 43

Atom mapping procedures to find corresponding atom pairs for interpolation and
the creation or annihilation of atoms to dummy types add complexities to
current single-topology alchemical RBFE protocols.44, 11, 45, 46 Single-
topology transformations often encode non-standard molecular topology formats
that require custom system setup workflows that are incompatible with standard
molecular visualization and trajectory analysis tools. In many
implementations, RBFE alchemical schemes are limited to R-group
transformations involving ligand pairs sharing a common scaffold.8, 47 Except
for a few commercial products,48, 49, 50 scaffold-hopping RBFE calculations
involving cyclization, ring expansion, linking, or any other transformation
that necessitate the breaking or the formation of chemical bonds,51 are not
generally supported.

We have recently developed the Alchemical Transfer Method (ATM) that resolves
many of the aforementioned complexities of RBFE calculations. The key
innovations of the method are 1) the mapping of the potential energy functions
of the unbound and bound states by a coordinate transformation rather than a
variation of parameters and 2) expressing the alchemical potential energy
function in a dual-topology formulation as the combination of the energy
functions of the physical end states rather than the interpolation of their
parameters.52, 53, 28, 34

ATM is not as affected by the complexities of traditional alchemical methods.
It supports absolute and relative binding free energy calculations in a
unified way, computes free energies directly employing a single simulation box
with standard chemical topologies, and natively supports diverse perturbations
(standard R-groups, charge-changing, and scaffold-hopping transformations)
without correction factors and ancillary calculations. Furthermore, since ATM
employs a dual-topology formalism and does not use parameter interpolation or
custom soft-core pair potentials, it is easier to implement and more
straightforward to transfer across MD engines because it uses the unmodified
energy routines of the underlying molecular dynamics engine. For the same
reason, it applies to any molecular energy function, including the next
generation of polarizable,54, 55, 56, 57 quantum-mechanical,58, 59, 60, 61 and
machine-learning potentials62, 63 that are starting to be employed in
macromolecular simulations. The current open-source software release of ATM
employs the OpenMM molecular dynamics engine and has been successfully tested
on a series of drug discovery targets with the AMBER molecular mechanics force
field in academic and industrial settings.34, 18

The simplifications and greater range of applicability afforded by the
alchemical transfer approach can be particularly useful in drug-discovery
deployments to screen large and diverse ligand libraries in a more streamlined
fashion. In this work, we validate ATM against the community benchmark
prepared by Schindler et al.,12 which contains examples of standard peripheral
group transformations as well as more challenging scaffold-hopping and charge-
changing transformations representative of real-world drug-discovery
applications.

Given the large number of calculations involved, aspects of the ATM workflow
have been automated, which was facilitated by the nature of the ATM approach
that avoids custom alchemical topologies and atom mapping typical of
conventional RBFE workflows.8, 36, 4, 64, 65 Knowing that the quality of the
potential energy functions can have a substantial effect on the free energy
prediction accuracy,66, 67 we applied a bespoke force field parameter
generation protocol for each of the ligands. While the force field generation
engine (FFEngine) is not publicly available, the parameters for each of the
ligands in this work have been included in the Supporting Information. Thus,
the work presented here is fully reproducible using the open-source version of
ATM and the published force field parameters.

The paper is organized as follows: We first introduce the theory of the
Alchemical Transfer Method and present some of the key implementation details.
We then describe the benchmarks sets, system setup, and alchemical simulation
details. The results are presented and analyzed next. The paper concludes with
a discussion of examples illustrating the strengths and weaknesses of the
method and identifying areas of improvement for the application of ATM to drug
discovery projects.

##  2 Theory and Methods

###  2.1 The Alchemical Transfer Method (ATM)

The alchemical Transfer Method (ATM) models the free energy difference between
two chemical states related by a coordinate transformation. One such example
is the molecular binding processes investigated here, represented as the
translation of the ligand from the solvent to the receptor binding site. The
method details are fully described in previously published works.53, 34, so an
abridged account is provided here. Briefly, consider the potential energy,
$U_{0}(x)$, of the unbound state ($R+A$) of the complex between a receptor R
and a ligand A when the ligand is in solution and where
$x=(x_{R},x_{A},x_{S})$ represents the coordinates of the receptor, ligand,
and solvent, respectively. ATM expresses the potential energy function,
$U_{1}(x)$, that describes the state $RA$ when the ligand is bound to the
receptor as

$U_{1}(x)=U_{0}(x_{R},x_{A}+h,x_{S})\,,$ (1)

where $h$ is a fixed displacement vector that brings the ligand from its
position in the solvent bulk to the receptor binding site (Figure 1). The free
energy difference between the bound and unbound states is then calculated by
Free Energy Perturbation (FEP), similar to standard binding free energy
methods. To this end, we define the perturbation energy function as

$u(x)=U_{1}(x)-U_{0}(x)=U_{0}(x_{R},x_{A}+h,x_{S})-U_{0}(x_{R},x_{A},x_{S})$
(2)

and introduce a $\lambda$-dependent alchemical potential energy function

$U_{\lambda}(x)=U_{0}(x)+W_{\lambda}[u(x)]\,,$ (3)

where $W_{\lambda}(u)$ is an alchemical perturbation function with the
properties $W_{0}(u)=0$ and $W_{1}(u)=u$, ensuring that Eq. (3) interpolates
from the initial unbound state $U_{0}(x)$ at $\lambda=0$ and the final bound
state $U_{1}(x)$ at $\lambda=1$. The linear alchemical perturbation function
$W_{\lambda}(u)=\lambda u$ is a common choice. In this work, we adopt a non-
linear expression described in Computational Details that yields faster
convergence than the default linear version.68

As elaborated elsewhere,53, 34 the alchemical path between the unbound and
bound endpoints is divided into two legs: one starting from the unbound state
at $\lambda=0$ using the alchemical potential in Eq. (3), and a second leg
starting from the bound state $U_{1}(x)$ morphing in the other direction
towards the unbound state using the alchemical potential
$U_{\lambda}(x)=U_{1}(x)+W_{1-\lambda}[-u(x)]$. Both legs terminate at
$\lambda=1/2$ at the ATM alchemical intermediate with the potential energy
function $U_{1/2}(x)=[U_{0}(x)+U_{1}(x)]/2$ that is an equally weighted
average of the unbound and bound states.

Figure 1:  Illustration of the Relative Binding Free Energy calculation with
the Alchemical Transfer Method. The first ligand is placed into the protein
binding pocket, and the second in the solvent. During the alchemical
transformation, the two ligands switch places. The resulting free energy
corresponds to the difference in their binding free energies.

The ATM formulation above is for an Absolute Binding Free Energy (ABFE)
calculation. Here, we are concerned with Relative Binding Free Energy (RBFE)
prediction, in which the binding of a ligand $B$ occurs simultaneously as
another ligand, $A$, leaves the receptor binding site. The free energy change
of this process is the difference between the binding free energies of the two
ligands. More specifically, an RBFE ATM calculation computes the free energy
change from the state $RA+B$ with ligand $A$ bound to the receptor to the
state $RB+A$ where ligand $B$ is bound to the receptor. In ATM, this process
is described by a coordinate transformation that translates ligand $B$ from a
position in the solvent to the binding site and simultaneously translates
ligand $A$ from the receptor binding site to the solvent. Using the notation
introduced above, the initial state $RA+B$ is described by the potential
energy function $U_{0}(x)$, $x=(x_{R},x_{A},x_{B},x_{S})$, and the potential
energy function $U_{1}(x)$ of the final state $RB+A$, where the positions of
the two ligands are switched, is written as

$U_{1}(x)=U_{0}(x_{R},x_{A}-h,x_{B}+h,x_{S})\,,$ (4)

where the displacement vector $h$ translates ligand $B$ into the binding site
while the opposite displacement $-h$ translates ligand $A$ in the reversed
direction (Figure 1). The nature of the coordinate transformation (a
translation of one ligand vs. the translation of two ligands in opposite
directions) is the only fundamental difference between ABFE and RBFE protocols
in ATM. The definitions of the perturbation energy, the alchemical potential
energy function [Eqs. (2)–(3)] and the alchemical intermediate are the same.

The main advantage of ATM for relative versus absolute binding free energy
calculations is computational efficiency. For absolute binding free energy
calculations it is necessary to converge the unbound state of the receptor,
which can be a slow process given the possibility of high-barrier
conformational rearrangements, differences in binding site solvent structure,
changes to the protonation/tautomerization states of binding site residues,
and other differences that may exist between the bound and unbound states. ATM
RBFE calculations always have a ligand in the binding site, thus minimizing
these effects. The actual free energy of the unbound receptor is
inconsequential in the context of RBFE calculations because it is a constant
for each ligand.

The alchemical intermediates in RBFE calculations are unphysical states in
which the ligands are present simultaneously in the receptor binding site and
solution, each at a strength proportional to the coupling parameter $\lambda$.
The more dissimilar the ligands, the more strained the conformations of the
system required to accommodate the alchemical intermediate states. This
characteristic is reflected in the high free energy barrier encountered at the
alchemical midpoint for some ligand pairs. The height of the free energy
barrier is used below as one of the proxies to judge the quality of RBFE
calculations and the confidence level of their estimates. The specific
settings of the RBFE protocol used in this work are described in Computational
Details.

####  2.1.1 FFEngine Ligand Force Field Assignment

FFEgine is a force field toolkit built with ParmEd69, RDKit70, GAFF271,
GFN2-xTB72, GPU-powered QM package Terachem73, 74 that provides high-quality
parameters for drug-like molecules based on a bespoke workflow where quantum
mechanical calculations are performed on each molecule to obtain the potential
surface for force field fitting. The atom types in FFEngine are defined
following the hierarchical structure described by Jin Z. et al.75 In total,
FFEngine utilizes approximately 200 atom types, with the objective of covering
the chemical space of drug-like molecules.

FFEngine generates ligand parameters for the AMBER force field functional
form76, for use in AMBER, GROMACS, OpenMM, or other software packages. The
charge assignment uses GFN2-xTB/BCC, which is similar to the AM1/BCC model
from AMBER. Atomic partial charges are assigned based on atom types and the
bond charge correction (BCC) parameters fitted to the HF/6-31G* electrostatic
potential (ESP) from 50,000 drug-like compounds and their conformations. In
GFN2-xTB/BCC, the pre-charges are assigned with the semiempirical method
GFN2-XTB. The molecular structures are relaxed with the machine learning force
field GFN-FF to prepare the structure before partial charges assignment.77 The
vdW, bond, angle, and torsion parameters are assigned based on a more refined
set of atom types. GAFF2 was used as the fallback parameters for bond, angle,
and torsional degrees of freedom.

####  2.1.2 DiffNet Analysis

ATM yields estimates of the differences between the binding free energies of
pairs of ligands. We employed DiffNet64, 78 to estimate the binding free
energies of the ligands in the set. Diffnet finds the set of binding free
energies most consistent with the network of their differences (referred to as
edges) and their uncertainties based on the maximum likelihood principle. The
DiffNet solution is known up to an arbitrary constant that we set to match the
average of the experimental binding free energies.

The uncertainties associated with the edges of the network of free energy
transformations are an important element of the DiffNet protocol. Edges with
small uncertainties weigh on the final solution more than edges with larger
uncertainties. To define confidence levels in the ATM RBFE predictions, we
used two measures of the quality of the alchemical calculations: 1) gaps in
the perturbation energy distributions and 2) height of the free energy barrier
along the alchemical path.

To assess the presence of gaps in the perturbation energy distributions, the
binding energy samples of each leg were histogrammed with bins of size $\Delta
u=k_{B}T/\Delta\lambda$, where $\Delta\lambda=1/N_{\rm states}$ is the size of
the subdivision of the alchemical path among $N_{\rm states}$ alchemical
states (approximately $13$ kcal/mol in this application). A case with a
sequence of two or more bins with zero counts was flagged as unlikely to be
converged, and its uncertainty was increased by $\Delta u$. A case with a
sequence of two or more bins with less than 10% of the expected number of
samples based on a uniform distribution was flagged as possibly unconverged,
and its uncertainty was increased by a tenth of $\Delta u$. Similarly, to
capture the lower confidence of predictions of large transformations with high
free energy intermediates, we increased the uncertainty of predictions where
the free energy of the alchemical intermediate, $\Delta\Delta G(1/2)$, exceeds
$30$ kcal/mol relative to either end state. Specifically, in these cases, we
increased the uncertainty linearly by $a(\Delta\Delta G(1/2)-g)$ where $a=0.1$
kcal/mol and $g=30$ kcal/mol. While this assessment successfully flagged
possibly problematic calculations, the specific parameters used here to assign
confidence levels were set empirically and would benefit from validation in
future studies.

####  2.1.3 Simulation Setup Workflow

To conduct this work, we wrote a Python-based ATM RBFE setup and analysis
workflow. The workflow performs water placement, force field generation,
reference atom selection, displacement vector searching, AToM-OpenMM
simulation, relative binding free energy calculation, and DiffNet analysis. As
input, the workflow requires fully prepared proteins and docked ligand poses.

In order to get better solvation structures for the receptors, 3D-RISM79 from
AmberTools80 was first applied to estimate the implicit water distribution on
the receptor-primary ligand complex surface, followed by Placevent81 to place
the explicit water molecules on the complex surface. Any water molecules that
clashed with other ligands in the set were removed. The workflow then assigned
ligand force field parameters using FFEngine as described above.

An ATM RBFE calculation requires the choice of three corresponding reference
atoms for the alignment of the ligand pair.34 The alignment reference atoms
cannot be colinear and are best chosen among rigid core atoms of the two
ligands. This task involves manual selection the reference atoms of only one
representative ligand placed into the receptor binding site. All other ligands
are then automatically aligned to the first by a minimum distance search based
on their initial poses.

In the simulation box for an ATM calculation, one ligand is placed in the
binding pocket, and the other one is translated into the solvent phase by a
displacement vector $h$. While the choice of the displacement vector could be
arbitrary, a reasonable choice should ensure that any atom of the second
ligand is at least 10 Å away from any of the receptor atoms. On the other
hand, a large displacement would unnecessarily increase the box size,
affecting performance. We used TLeap80 to generate a preliminary rectangular
solvation box for the receptor with a 10 Å water buffer to select a good
displacement vector automatically. The ATM displacement vector $h$ was then
obtained by displacing the ligand from the binding pocket towards the center
of the box face with the smallest surface area until all atoms of the ligand
were found outside the box.

The Amber input chemical topologies with the assigned force field parameters
generated by our automated workflow are listed in the simulation input files
available on the GitHub repository at
https://github.com/EricChen521/ATM_MerckSet.

###  2.2 Benchmark Systems

The benchmark sets prepared by Schindler et al.12 include eight receptor
targets with 24 to 44 ligands for each target (Table 1). The benchmark also
specifies the ligand pairs forming the edges of the graph of RBFE
calculations.82 In total, the benchmark set includes 264 protein-ligand
complexes and 550 RBFE edges. Schindler et al. reported free energy estimates
for 525 of the 550 transformations using Schrödinger’s FEP+ package.10 In this
work, we considered all of the protein-ligand complexes and the corresponding
edges except the two involving compound 28 of the CDK8 set, which we suspect
to have been misidentified. This compound lacks a measured binding free energy
and is classified as a non-binder, even though close analogs are strong
binders.83

The benchmark sets provided by Schindler et al.12 are considered challenging
because, in addition to standard small R-group modifications, they include a
significant number of more challenging transformations. The set includes 43
large R-group transformations (more than 10 added/removed atoms), 66 charge
transformations (either changing the formal charge or moving the location of
the formal charge), and 60 scaffold-hopping transformations. Large R-group
transformations are considered challenging because they often induce changes
in the conformation and hydration pattern of the complex to accommodate the
new groups of atoms. Charge-changing RBFEs, which involve ligands with
different net charges, have traditionally required specialized strategies or
correction terms20, 84, 43 unnecessary in our alchemical transfer approach.34
Nevertheless, charge-changing transformations and the related charge-shifting
transformations remain challenging because of the receptor and solvent
reorganization that they often induce. Finally, scaffold-hopping
transformations that include cyclization, ring-breaking, and ring
expansion/reduction transformations that involve the breaking or forming of
chemical bonds, which traditionally require specialized strategies,48, 50 are
generally straightforward and are treated here with ATM as any other
transformation.34

###  2.3 Computational Details

###  2.4 Molecular Systems Setup

We employed the structures of the proteins and bound ligands as provided by
Schindler et al.12 The simulation systems were prepared for the RBFE
calculations using the automated workflow described above. The AMBER FF14SB
force field85 was used for the protein and the TIP3P model86 for water. K+ and
Cl- ions were added to the system if needed for neutralization. The force
constants of the ligand alignment restraints were set to $k_{r}=2.5$ kcal/(mol
Å2) for the position restraint, and $k_{\theta}=k_{\psi}=50$ kcal/mol for the
orientational restraints.34 The C$\alpha$ atoms’ positions of the receptors
were restrained to their initial values using flat-bottom harmonic restraints
with a tolerance of $1.5$ Å and a force constant of $25$ kcal/(mol Å2).

The solvated systems were energy minimized, thermalized, and relaxed for 400
ps at 298 K and 1 bar constant pressure and annealed to the $\lambda=1/2$
alchemical intermediate in 500 ps while restraining the receptor and the
ligands’ atoms. The system was then equilibrated in the NVT ensemble at 298 K
for 300 ps at the alchemical intermediate after releasing the restraints on
the ligands and the receptor (except for the C$\alpha$ atoms positional
restraints described above). The resulting structures were used as starting
configurations for the alchemical replica exchange simulations described next.

###  2.5 Alchemical Transfer Relative Binding Free Energy Protocol

ATM RBFE calculations were conducted using the AToM-OpenMM package version
3.2.3,87 the ATM MetaForce OpenMM plugin version 0.3.1,88 and the OpenMM MD
engine version 7.7.89. The ATM alchemical schedule is comprised of two
legs.53, 34 Using the notation introduced in Theory and Methods, the first leg
corresponds to the transformation from the $RA+B$ state described by the
potential energy function $U_{0}(x)$ to the alchemical intermediate state
$R(AB)_{1/2}+(AB)_{1/2}$ described by the potential energy function
$U_{0}(x)/2+U_{1}(x)/2$. The second leg reaches the same alchemical
intermediate state but starts from the final $RB+A$ state. The difference
between the free energy changes of the first leg and second legs yields the
relative binding free energy $\Delta\Delta G_{b}^{\circ}(B,A)$ of ligand B
with respect to ligand A. The alchemical potential energy function for the
first leg is given in Eq. (3) while that for the second leg is

$U_{1-\lambda}(x)=U_{1}(x)+W_{\lambda}[-u(x)]\,,$ (5)

where $0\leq\lambda\leq 1/2$, the perturbation energy is given in Eq. (2),
and, for both legs, the alchemical perturbation function is the softplus
function68

$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).$ (6)

The parameters $\lambda_{2}$, $\lambda_{1}$, $\alpha$, and $u_{0}$ are
functions of $\lambda$ (see below).52. The function

$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}$
(7)

with

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

and

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

with, in this work, $u_{\rm max}=200$ kcal/mol, $u_{c}=100$ kcal/mol, and
$a=1/16$, is the soft-core perturbation energy function designed to avoid
singularities near the initial state of the alchemical transformation.52

A schedule of 11 equispaced $\lambda$-states from $\lambda=0$ to $\lambda=1/2$
was employed for each of the two legs. The corresponding schedules of softplus
alchemical parameters in Eq. (6) were:
$\lambda_{1}=0,0,0,0,0,0,0.1,0.2,0.3,0.4,0.5$,
$\lambda_{2}=0,0.1,0.2,0.3,0.4,0.5,0.5,0.5,0.5,0.5,0.5$, and $\alpha=0.1$
(kcal/mol)-1 and $u_{0}=110$ kcal/mol for all $\lambda$-states. We employed
this alchemical schedule for all the ligand pairs of the benchmark sets in
this work without further optimization.

For each run, asynchronous Hamiltonian replica exchange90 molecular dynamics
conformational sampling was performed with a 2 fs timestep and a replica
running time of 40 ps on 4 GPUs using the AToM-OpenMM software.87 Exchanges
between the two equivalent alchemical intermediate states at $\lambda=1/2$
allow replicas to transition from one alchemical leg to the other.
Perturbation energy samples were collected every 40 ps. Relative binding free
energies were computed from replica trajectories at least 5 ns long
(approximately 110 ns of MD in aggregate per ligand pair), discarding the
first third of the samples for equilibration. UWHAM multi-state analysis91 was
used for free energy and statistical error estimation.

The molecular system files, AToM-OpenMM input files, and UWHAM analysis code
are available on GitHub at https://github.com/EricChen521/ATM_MerckSet.

##  3 Results

The relative binding free energy performance of the automated ATM protocol
compared with the experimental measurements of the benchmark sets is
summarized in Figure 2, where each point represents a ligand pair. For each
target set, we report the Average Unsigned Error (AUE), the Pearson’s
correlation coefficient $r$, and the fraction of concordant predictions
$f_{\rm concd.}$ relative to the experiments. A concordant prediction is a
case in which the direction of the change in binding affinity of more than
$0.5$ kcal/mol is correctly predicted (see Methods).

Figure 2: Scatter plots of the calculated RBFEs vs. the corresponding
differences of experimental binding free energies. The 1:1 line represents
perfect agreement. The dark and light-shaded areas correspond to predictions
within 1 kcal/mol and 2 kcal/mol of the experiments, respectively. For each
set the Average Unsigned Error (AUE), Pearson’s correlation coefficient ($r$),
and the fraction of concordant predictions ($f_{\rm concd.}$) are reported.

According to these statistics, the c-Met, PFKFB3, CDK8, and Hif-2$\alpha$
predictions are better than the other targets, with correlation coefficients
of 0.6 or higher with over 80% of concordant pairs. The performance for the
c-Met set with $r=0.84$ and $f_{\rm concd.}=86$% is particularly encouraging.
However, due to a relatively small number of outliers, the AUEs for these sets
are quite high (from $0.98$ to $1.5$ kcal/mol). The RBFE prediction
performance for the remaining four sets (Eg5, SHP2, Syk, and TNKS2) is not as
good. The Syk and TNKS2 set have fewer major outliers and display relatively
small AUEs. However, the calculated RBFEs are poorly correlated with the
experiments with correlation coefficients of 0.52 and 0.40, respectively. The
SHP2 and, particularly, the Eg5 sets have both a significant number of
outliers and poor correlation with the experiments.

The confidence levels of the RBFE predictions represented by the error bars in
Figure 2 vary significantly from set to set. As described in the Methods,
confidence levels reflect the difficulty of the alchemical calculation. They
are assigned based on statistical fluctuations, connectedness, and the height
of the free energy barrier of the alchemical pathway. Surprisingly, the
calculations for c-Met, which are the closest to the experiments, generally
have higher uncertainty than the other sets. The RBFE calculations for the Eg5
and Syk sets are also considered challenging. With a few exceptions, the RBFEs
for the other sets are assigned with high confidence. The types of RBFE
transformations and their challenges are further explored in the Discussion.

The accuracy of the absolute binding free energies (ABFE) predictions
calculated from the RBFEs using the Diffnet protocol (see Methods) is
summarized in Table 1 and illustrated in Figure 3. The Pearson’s correlation
coefficients, Kendall’s rank order correlation coefficients, and average
unsigned errors relative to the experimental binding affinities averaged over
the eight systems are $r=0.60$, $\tau=0.46$, and AUE $=1.11$ kcal/mol,
respectively. These overall statistical measures are encouraging, given the
unsupervised nature of the ATM calculation workflow and the challenges of the
benchmarks. Results are slightly inferior to those obtained for the same
benchmark set with the FEP+ software package ($r=0.62$, $\tau=0.50$, and AUE
$=0.98$ kcal/mol).12

Interestingly, the triangulation of the data afforded by DiffNet results in
ABFEs significantly tighter agreement with the experiments than the RBFEs
estimates. This is particularly noticeable for the c-Met set, which, despite
the major outliers and large uncertainties of the RBFEs (Figure 2), yields
high-confidence predictions with low AUEs (Figure 3).

The ATM ABFE predictions for three of the eight systems (c-Met, PFKBF3, and
CDK8) have high correlation coefficients and generally low AUEs relative to
the experiments. FEP+, based on a different alchemical approach, performed
similarly well for these systems, suggesting that they are tractable for
diverse free energy methods. Although affected by a few major outliers, the
performance on the Hif-2$\alpha$ set is also generally good. Conversely, ATM
performed noticeably worse for the SHP2 and Syk sets, where it resulted in an
unusually poor AUE relative to the experiments. ATM and FEP+ performed equally
poorly for the Syk and Eg5 systems indicating a common set of challenges. The
TNKS2 set is less informative due to the small range of experimental
affinities (Figure 3).

Table 1:  Statistical measures of the agreement between experimental and
computed binding free energy values obtained with the ATM method compared with
literature results using FEP+.12

Target |  |  | ATM |  | FEP+12
---|---|---|---|---|---
|  $N$a |  |  $r$b |  $\tau$c | AUEd |  |  $r$b |  $\tau$c | AUEd
CDK8 | 32 |  | 0.79 | 0.66 | 1.35 |  | 0.62 | 0.57 | 1.20
c-Met | 24 |  | 0.94 | 0.82 | 0.73 |  | 0.90 | 0.73 | 0.82
Eg5 | 28 |  | 0.49 | 0.39 | 1.06 |  | 0.47 | 0.54 | 1.09
HIF-2$\alpha$ | 42 |  | 0.49 | 0.39 | 0.95 |  | 0.61 | 0.45 | 0.84
PFKFB3 | 40 |  | 0.72 | 0.50 | 1.07 |  | 0.79 | 0.60 | 1.09
SHP-2 | 26 |  | 0.58 | 0.39 | 1.32 |  | 0.71 | 0.61 | 0.74
Syk | 44 |  | 0.44 | 0.24 | 1.34 |  | 0.50 | 0.29 | 0.85
TNKS2 | 27 |  | 0.36 | 0.27 | 1.02 |  | 0.40 | 0.29 | 1.23
Cumulative | 263 |  | 0.60 | 0.46 | 1.11 |  | 0.62 | 0.50 | 0.98

a Number of complexes. b Pearson’s correlation coefficient. c Kendall’s rank
order correlation coefficient. d Averaged unsigned error in kcal/mol.

Figure 3: Scatter plots of the Diffnet ABFE estimates from the ATM RBFE values
vs. the corresponding experimental binding free energies. The 1:1 line
represents perfect agreement. The horizontal directed lines denote ligands for
which the experimental affinity is known only within a range. The dark and
light-shaded areas correspond to predictions within 1 kcal/mol and 2 kcal/mol
of the experiments, respectively. For each set, the Average Unsigned Error
(AUE), Pearson’s correlation coefficient ($r$), and the Kendall’s $\tau$ rank
order correlation coefficient are reported.

##  4 Discussion

The results obtained in this validation study confirm the applicability of the
ATM approach to large-scale binding free energy estimation campaigns for
challenging and diverse ligand libraries. In this section, we present some
examples illustrating the strengths and weaknesses of the method that we
observed during this work. We plan to incorporate these insights in deploying
this approach in future work.

###  4.1 Small Charge-Preserving R-group RBFE Estimates are Generally Well-
Converged

Over 70% of the ligand pairs of the benchmark (384 out of 548) are classified
as conventional R-group alchemical transformations involving one or more small
peripheral substitutes of the same net charge with 10 or fewer heavy atoms.
The calculations of these pairs are expected to converge rapidly, and the
corresponding relative binding free energy estimates to be more reliable than
the other more challenging cases. As illustrated by the example in Figure 4,
this expectation is largely confirmed by the consistently good quality
measures in terms of the overlaps between perturbation energy distributions
and the moderate height of the free energy profile of R-group transformations.
The RBFE estimates of R-group transformations are also generally in closer
agreement and have fewer outliers relative to the experiments than other
transformation types. For example, the RMSE of the R-group transformations for
the CDK8, c-Met, and HIF-2$\alpha$ sets are $1.68$, $1.21$, and $1.53$
kcal/mol, respectively, compared to $1.87$, $1.38$, and $1.81$ kcal/mol for
the entire sets.

Figure 4: The results of the RBFE calculation between the CHEMBL3402758 and
CHEMBL3402760 ligands of the c-Met receptor illustrating a reliable small
R-group transformation. (A) Chemical structures of the ligands and calculated
and experimental RBFE. (B) Perturbation energy distributions for the
alchemical states from $\lambda=0$ to $\lambda=1/2$ (from right to left in
color sequence) showing good overlaps between nearby distributions. (C) Free
energy profile along the alchemical transformation displaying a moderate free
energy barrier.

Conversely, as illustrated by the Syk pair in Figure 5, RBFE estimates between
very dissimilar ligands, even if involving only peripheral groups without
charge variations, are often deemed unconverged and unlikely to be predictive
of the actual difference in binding affinity. We generally observe this
behavior when R-group modifications involve more than 10 heavy atoms. These
cases (43 out of 548) often present gaps in the sequence of perturbation
energy distributions along the alchemical path and display a free energy
barrier of the free energy profile of 40 kcal/mol or more (Figure 5). Many of
the transformations in this category are outliers compared to the experiments.
For example, nearly half of the pairs of the Syk set in this category deviate
from the experiments by more than 2 kcal/mol. Transformations between very
dissimilar ligands such as this are probably best implemented as sequences of
smaller and more manageable alchemical steps.

Figure 5: The results of the RBFE calculation between the CHEMBL3265035 and
CHEMBL3264996 ligands of the Syk receptor illustrating an unreliable large
R-group transformation. (A) Chemical structures of the ligands and calculated
and experimental RBFE. (B) Perturbation energy distributions for the
alchemical states from $\lambda=0$ to $\lambda=1/2$ (from right to left in
color sequence) showing a large gap between high-energy and low-energy
distributions. (C) Free energy profile along the alchemical transformation
displaying a large free energy barrier.

###  4.2 Charge-changing and charge-shifting Transformations are Challenging

In contrast to small R-group transformations, RBFE calculations involving a
change of net charge of the ligand or a shift of a charged group from one
position to another often display large uncertainties. The c-Met pair in
Figure 6 where a positively charged alkylamine moiety is added to a neutral
carbamate substituent, is an example of this class. While the free energy of
the alchemical intermediate is high (above 40 kcal/mol), the sequence of
alchemical states is relatively well connected. The Eg5 case in Figure 7
includes a shift by three bonds of a charged ammonium group. This alchemical
transformation is assigned a large uncertainty because it is characterized by
both a disconnected alchemical path and a high free energy intermediate.

Overall, the 66 RBFE predictions classified as either charge-changing or
charge-shifting transformations have larger uncertainties and poorer agreement
with the experiments than average. This is especially so for the
transformations in this class for the SHP2 set that are found to have an RMSE
relative to the experiments in excess of 3 kcal/mol. Nearly all edges with
large uncertainties for the c-Met and Eg5 sets (Figure 2) correspond to
charge-changing or charge-shifting transformations. Achieving better
convergence for these transformations would significantly improve ATM’s
promising overall prediction accuracy for these sets.

Figure 6: The results of the RBFE calculation between the CHEMBL3402745 and
CHEMBL3402743 ligands of the c-Met receptor illustrating a possibly reliable
charge-changing transformation. (A) Chemical structures of the ligands and
calculated and experimental RBFE. (B) Perturbation energy distributions for
the alchemical states from $\lambda=0$ to $\lambda=1/2$ (from right to left in
color sequence) showing mostly good overlaps but an undersampled region near
$u_{\rm sc}\simeq 160$ kcal/mol. (C) Free energy profile along the alchemical
transformation displaying a large free energy barrier. Figure 7: The results
of the RBFE calculation between the CHEMBL1084431 and CHEMBL1085666 ligands of
the Eg5 receptor illustrating an unreliable charge-shifting transformation.
(A) Chemical structures of the ligands and calculated and experimental RBFE.
(B) Perturbation energy distributions for the alchemical states from
$\lambda=0$ to $\lambda=1/2$ (from right to left in color sequence) showing a
gap between high-energy and low-energy distributions and a bimodal
distribution (pink). (C) Free energy profile along the alchemical
transformation displaying a large free energy barrier. Figure 8: The
structures of the complexes of CHEMBL1084431 (cyan carbon atoms, ribbon) and
CHEMBL1085666 (green carbon atoms, ribbon) with the Eg5 receptor of Figure 7.
A conformational reorganization occurs whereby the salt bridge between the
ligand and the receptor (yellow bonds) switches from Glu102 to the Glu104
residue going from the first complex to the second. The structures are taken
at the end of the trajectory of the alchemical end-states.

Because the system remains neutral, the poorer outcomes observed here for
alchemical transfer transformations involving charge variations are not
obviously due to systematic boundary conditions and finite system biases
present in some double-decoupling protocols.84, 43 The convergence issues
observed with charged groups are more likely related to the conformational
reorganization of the complex that occurs when interactions are varied by
placing or relocating a charged group of the ligand. This is illustrated by
the case in Figure 7 where an amide hydrogen bonding acceptor group is
replaced by a strong ammonium hydrogen bonding donor originally placed three
bonds away. In this case, we observe that the two ligands shift position and
form salt bridges with two different glutamate residues of the receptor
(Figure 8). The equilibration between these two conformational states of the
complex hinders convergence because it occurs slowly relative to the
timescales of the alchemical simulations. In addition, because the receptor
needs to accommodate both charged groups of the ligand simultaneously,
conformational frustration at the alchemical intermediate leads to high free
energies and gaps and lack of state overlaps along the alchemical pathway
(Figure 7).

This and the previous examples underscore the various ways in which
differences between ligand pairs prevent successful RBFE predictions. Two
ligands can have very different sizes and shapes, as in the example of Figure
5, or, as in the present example of Figure 7, they can be structurally similar
but differ radically in the way that corresponding groups interact with their
environments. Cases with very dissimilar ligands such as these can probably be
addressed by breaking the alchemical transformation into smaller steps by
inserting suitable chemical intermediates. For example, the transformation in
Figure 7 that attempts to replace one charged group with another
simultaneously causing a large conformational reorganization (Figure 8) could
be implemented by first removing or neutralizing one ammonium group and then
inserting the other in a second step. The case in Figure 6 and similar others
we observed indicate that individual insertions of charged groups are more
manageable than replacing one charged group with another placed elsewhere.

###  4.3 Scaffold-hopping Transformations are as Straightforward as R-group
Transformations

The Alchemical Transfer Method is based on a dual-topology representation,
making it easier to set up scaffold-hopping transformations between ligands
that do not share the same core topology. Indeed, the ATM setup procedure
described here for scaffold-hopping transformations is the same as any other
transformation. Sixty of the 548 RBFE calculations conducted in this work were
classified as scaffold hopping transformations.

As illustrated by the examples in Figures 9 and 10, the results of this
validation campaign have also shown that scaffold-hopping RBFE calculations
converge just as efficiently as small R-group transformations of the same
kind. In the HIF-2$\alpha$ case in Figure 9, for example, a 5-membered ring is
formed by the cyclization of two substituents of the central aromatic group.
In the case of Figure 10, a 5-membered ring is expanded to a 6-membered ring.
In both of these cases and many similar ones we observed, there are good
overlaps between perturbation energy distributions, and the free energy of the
alchemical intermediate is moderate, just as in small R-group transformations
(Figure 4)

Figure 9: The results of the RBFE calculation between the ligands 43 and 235
of the HIF-2$\alpha$ receptor illustrating a reliable small scaffold-hopping
cyclization transformation. (A) Chemical structures of the ligands and
calculated and experimental RBFE. (B) Perturbation energy distributions for
the alchemical states from $\lambda=0$ to $\lambda=1/2$ (from right to left in
color sequence) showing good overlaps between nearby distributions. (C) Free
energy profile along the alchemical transformation displaying a small free
energy barrier.  Figure 10: The results of the RBFE calculation between the
ligands Example 14 and 6 of the SHP2 receptor illustrating a reliable small
scaffold-hopping ring expansion transformation. (A) Chemical structures of the
ligands and calculated and experimental RBFE. (B) Perturbation energy
distributions for the alchemical states from $\lambda=0$ to $\lambda=1/2$
(from right to left in color sequence) showing good overlaps between nearby
distributions. (C) Free energy profile along the alchemical transformation
displaying a small free energy barrier.

##  5 Conclusions

We tested the Alchemical Transfer Method (ATM) as implemented in the open
source AToM-OpenMM software package with a bespoke ligand force field on the
large and challenging relative binding free energy benchmark sets developed by
Schindler et al.12 In this work, 548 RBFE calculations were set up and
performed largely unsupervised using our Python workflow facilitated by the
minimal ATM customization requirements. The raw RBFE estimates from ATM showed
higher statistical fluctuations relative to FEP+. Nevertheless, the absolute
binding free energy estimates produced by the DiffNet algorithm on this data
provided an accuracy relative to the experiments comparable on average to
FEP+. The results indicate that standard small R-group transformations
converge rapidly, as expected. RBFE estimates for large R-group
transformations and charge-changing and charge-shifting transformations are
less reliable due to the slow convergence rate of conformational
reorganization effects induced by the significant changes in molecular size
and ligand-receptor interactions.

Based on the insights from this work, we plan to break up difficult RBFE
transformations into multiple, more manageable steps. For example, large
R-group transformations will be implemented as a series of smaller
transformations involving suitable chemical intermediates. Similarly, to
minimize the impact of conformational reorganization effects, we will consider
breaking up transformations involving the replacement of an exclusive hydrogen
bond donor with an exclusive hydrogen bond acceptor using chemical
intermediates containing groups, such as hydroxyls, with promiscuous hydrogen
bond behavior. This study further confirms that scaffold-hopping
transformations involving the formation or breaking of chemical bonds do not
provide additional challenges with this method. In conclusion, this study
confirms that ATM is a promising production tool for lead optimization in
structure-based drug discovery.

##  6 Acknowledgments

We acknowledge support from the National Science Foundation (NSF CAREER
1750511 to E.G.). We are grateful for the computer time provided on Roivant’s
Neo computational cluster.

## References

  * 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
  * Zhang et al. 2021 Zhang, C.-H.; Stone, E. A.; Deshmukh, M.; Ippolito, J. A.; Ghahremanpour, M. M.; Tirado-Rives, J.; Spasov, K. A.; Zhang, S.; Takeo, Y.; Kudalkar, S. N., et al. Potent noncovalent inhibitors of the main protease of SARS-CoV-2 from molecular sculpting of the drug perampanel guided by free energy perturbation calculations. _ACS Central Science_ 2021, _7_ , 467–475
  * Allen et al. 2022 Allen, B. K.; Kulkarni, M. M.; Chamberlain, B.; Dwight, T.; Koh, C.; Samant, R.; Jernigan, F.; Rice, J.; Tan, D.; Li, S.; Marino, K.; Huang, H.; Chiswick, E.; Tesar, B.; Sparks, S.; Lin, Z.; McGee, T. D.; Kolossváry, I.; Lin, C.; Shechter, S.; Soutter, H.; Bastos, C.; Taimi, M.; Lai, S.; Petrin, A.; Kane, T.; Swann, S.; Gardner, H.; Winter, C.; Sherman, W. Design of a systemic small molecule clinical STING agonist using physics-based simulations and artificial intelligence. _bioRxiv_ 2022,
  * Ganguly et al. 2022 Ganguly, A.; Tsai, H.-C.; Fernández-Pendás, M.; Lee, T.-S.; Giese, T. J.; York, D. M. AMBER Drug Discovery Boost Tools: Automated Workflow for Production Free-Energy Simulation Setup and Analysis (ProFESSA). _J. Chem. Inf. Model._ 2022, _62_ , 6069–6083
  * Xu 2022 Xu, H. The slow but steady rise of binding free energy calculations in drug discovery. _J. Comp.-Aid. Mol. Des._ 2022, 1–8
  * Mobley 2012 Mobley, D. L. Lets get honest about sampling. _J. Comp. Aided Mol. Des._ 2012, _26_ , 93–95
  * 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, _60_ , 5595–5623
  * Gallicchio et al. 2014 Gallicchio, E.; Deng, N.; He, P.; Perryman, A. L.; Santiago, D. N.; Forli, S.; Olson, A. J.; Levy, R. M. Virtual Screening of Integrase Inhibitors by Large Scale Binding Free Energy Calculations: the SAMPL4 Challenge. _J. Comput.-Aided Mol. Des._ 2014, _28_ , 475–490
  * 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
  * 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. Comput.-Aided Mol. Des._ 2019, _33_ , 1021–1029
  * Schindler et al. 2020 Schindler, C. E.; Baumann, H.; Blum, A.; Böse, D.; Buchstaller, H.-P.; Burgdorf, L.; Cappel, D.; Chekler, E.; Czodrowski, P.; Dorsch, D.; I, E. M. K.; Follows, B.; Fuchs, T.; Grädler, U.; Gunera, J.; Johnson, T.; Lebrun, C. J.; Karra, S.; Klein, M.; Knehans, T.; Koetzner, L.; Krier, M.; Leiendecker, M.; Leuthner, B.; Li, L.; Mochalkin, I.; Musil, D.; Neagu, C.; Rippmann, F.; Schiemann, K.; Schulz, R.; Steinbrecher, T.; Tanzer, E.-M.; Unzue Lopez, A.; Viacava Follis, A.; Wegener, A.; Kuhn, D. Large-scale assessment of binding free energy calculations in active drug discovery projects. _J. Chem. Inf. Model._ 2020, _60_ , 5457–5474
  * Kuhn et al. 2020 Kuhn, M.; Firth-Clark, S.; Tosco, P.; Mey, A. S.; Mackey, M.; Michel, J. Assessment of binding affinity via alchemical free-energy calculations. _J. Chem. Inf. Model._ 2020, _60_ , 3120–3130
  * Gapsys et al. 2020 Gapsys, V.; Pérez-Benito, L.; Aldeghi, M.; Seeliger, D.; Van Vlijmen, H.; Tresadern, G.; De Groot, B. L. Large scale relative protein ligand binding affinities using non-equilibrium alchemy. _Chemical Science_ 2020, _11_ , 1140–1152
  * Bieniek et al. 2021 Bieniek, M. K.; Bhati, A. P.; Wan, S.; Coveney, P. V. TIES 20: Relative binding free energy with a flexible superimposition algorithm and partial ring morphing. _J. Chem. Theory Comput._ 2021, _17_ , 1250–1265
  * Hahn et al. 2022 Hahn, D.; Bayly, C.; Boby, M. L.; Bruce Macdonald, H.; Chodera, J.; Gapsys, V.; Mey, A.; Mobley, D.; Perez Benito, L.; Schindler, C.; Tresadern, G.; Warren, G. Best Practices for Constructing, Preparing, and Evaluating Protein-Ligand Binding Affinity Benchmarks [Article v1.0]. _Living Journal of Computational Molecular Science_ 2022, _4_ , 1497
  * Gapsys et al. 2022 Gapsys, V.; Hahn, D. F.; Tresadern, G.; Mobley, D. L.; Rampp, M.; de Groot, B. L. Pre-exascale computing of protein–ligand binding free energies with open source software for drug design. _J. Chem. Inf. Model._ 2022, _62_ , 1172–1177
  * Sabanés Zariquiey et al. 2023 Sabanés Zariquiey, F.; Pérez, A.; Majewski, M.; Gallicchio, E.; De Fabritiis, G. Validation of the Alchemical Transfer Method for the Estimation of Relative Binding Affinities of Molecular Series. _J. Chem. Inf. Model._ 2023, _63_ , 2438–2444
  * 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
  * Pan et al. 2017 Pan, A. C.; Xu, H.; Palpant, T.; Shaw, D. E. Quantitative characterization of the binding and unbinding of millimolar drug fragments with molecular dynamics simulations. _J. Chem. Theory Comput._ 2017, _13_ , 3372–3377
  * Deng and Roux 2009 Deng, Y.; Roux, B. Computations of standard binding free energies with molecular dynamics simulations. _J. Phys. Chem. B_ 2009, _113_ , 2234–2246
  * Gumbart et al. 2013 Gumbart, J. C.; Roux, B.; Chipot, C. Efficient determination of protein–protein standard binding free energies from first principles. _J. Chem. Theory Comput._ 2013, _9_ , 3789–3798
  * Velez-Vega and Gilson 2013 Velez-Vega, C.; Gilson, M. K. Overcoming dissipation in the calculation of standard binding free energies by ligand extraction. _J. Comp. Chem._ 2013, _34_ , 2360–2371
  * Limongelli et al. 2013 Limongelli, V.; Bonomi, M.; Parrinello, M. Funnel metadynamics as accurate binding free-energy method. _Proc. Natl. Acad. Sci._ 2013, _110_ , 6358–6363
  * Lapelosa 2017 Lapelosa, M. Free energy of binding and mechanism of interaction for the MEEVD-TPR2A peptide–protein complex. _J. Chem. Theory Comput._ 2017, _13_ , 4514–4523
  * Deng et al. 2018 Deng, N.; Cui, D.; Zhang, B. W.; Xia, J.; Cruz, J.; Levy, R. Comparing alchemical and physical pathway methods for computing the absolute binding free energy of charged ligands. _Phys. Chem. Chem. Phys._ 2018, _20_ , 17081–17092
  * Cruz et al. 2020 Cruz, J.; Wickstrom, L.; Yang, D.; Gallicchio, E.; Deng, N. Combining Alchemical Transformation with a Physical Pathway to Accelerate Absolute Binding Free Energy Calculations of Charged Ligands to Enclosed Binding Sites. _J. Chem. Theory Comput._ 2020, _16_ , 2803–2813
  * Azimi et al. 2022 Azimi, S.; Wu, J. Z.; Khuttan, S.; Kurtzman, T.; Deng, N.; Gallicchio, E. Application of the alchemical transfer and potential of mean force methods to the SAMPL8 host-guest blinded challenge. _J. Comp. Aid. Mol. Des._ 2022, _36_ , 63–76
  * Comer et al. 2015 Comer, J.; Gumbart, J. C.; Hénin, J.; Lelièvre, T.; Pohorille, A.; Chipot, C. The adaptive biasing force method: Everything you always wanted to know but were afraid to ask. _The Journal of Physical Chemistry B_ 2015, _119_ , 1129–1151
  * Mahinthichaichan et al. 2021 Mahinthichaichan, P.; Vo, Q. N.; Ellis, C. R.; Shen, J. Kinetics and Mechanism of Fentanyl Dissociation from the $\mu$-Opioid Receptor. _JACS Au_ 2021, _1_ , 2208–2215
  * Tse et al. 2020 Tse, C.; Wickstrom, L.; Kvaratskhelia, M.; Gallicchio, E.; Levy, R.; Deng, N. Exploring the free-energy landscape and thermodynamics of protein-protein association. _Biophys. J._ 2020, _119_ , 1226–1238
  * Jorgensen 2004 Jorgensen, W. L. The many roles of computation in drug discovery. _Science_ 2004, _303_ , 1813–1818
  * 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
  * Azimi et al. 2022 Azimi, S.; Khuttan, S.; Wu, J. Z.; Pal, R. K.; Gallicchio, E. Relative binding free energy calculations for ligands with diverse scaffolds with the alchemical transfer method. _J. Chem. Inf. Model._ 2022, _62_ , 309–323
  * 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. _J. Chem. Inf. Model._ 2020, _60_ , 4153–4169
  * 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. Comput.-Aided Mol. Des._ 2013, _27_ , 755–770
  * 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, _16_ , 5512–5525
  * 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
  * Wallace and Shen 2012 Wallace, J. A.; Shen, J. K. Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH. _J. Chem. Phys._ 2012, _137_ , 184105
  * 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
  * 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. _J. Chem. Phys._ 2013, _139_ , 11B606_1
  * Fleck et al. 2021 Fleck, M.; Wieder, M.; Boresch, S. Dummy atoms in alchemical free energy calculations. _J. Chem. Theory Comput._ 2021, _17_ , 4403–4419
  * 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
  * Gallicchio 2022 Gallicchio, E. In _Computational Peptide Science: Methods and Protocols_ ; Simonson, T., Ed.; Methods in Molecular Biology; Springer Nature, 2022; pp 303–334
  * Zhang et al. 2021 Zhang, H.; Kim, S.; Giese, T. J.; Lee, T.-S.; Lee, J.; York, D. M.; Im, W. CHARMM-GUI Free Energy Calculator for Practical Ligand Binding Free Energy Simulations with AMBER. _J. Chem. Inf. Model._ 2021, _61_ , 4145–4151
  * 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
  * 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. _J. Chem. Theory Comput._ 2020, _16_ , 7895–7914
  * 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,
  * Liu et al. 2015 Liu, S.; Wang, L.; Mobley, D. L. Is ring breaking feasible in relative binding free energy calculations? _J. Chem. Inf. Model._ 2015, _55_ , 727–735
  * 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
  * Wu et al. 2021 Wu, J. Z.; Azimi, S.; Khuttan, S.; Deng, N.; Gallicchio, E. Alchemical transfer approach to absolute binding free energy estimation. _Journal of Chemical Theory and Computation_ 2021, _17_ , 3309–3319
  * 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
  * Huang et al. 2018 Huang, J.; Lemkul, J. A.; Eastman, P. K.; MacKerell Jr., A. D. Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks. _J. Comput. Chem._ 2018, _39_ , 1682–1689
  * Das et al. 2022 Das, A. K.; Liu, M.; Head-Gordon, T. Development of a Many-Body Force Field for Aqueous Alkali Metal and Halogen Ions: An Energy Decomposition Analysis Guided Approach. _J. Chem. Theory and Comput._ 2022, _18_ , 953–967
  * 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. Protein 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
  * Casalino et al. 2020 Casalino, L.; Nierzwicki, Ł.; Jinek, M.; Palermo, G. Catalytic mechanism of non-target DNA cleavage in CRISPR-Cas9 revealed by ab initio molecular dynamics. _ACS catalysis_ 2020, _10_ , 13596–13605
  * 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 Commun._ 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,
  * Xu 2019 Xu, H. Optimal measurement network of pairwise differences. _J. Chem. Inf. Model._ 2019, _59_ , 4720–4728
  * Li et al. 2021 Li, P.; Li, Z.; Wang, Y.; Dou, H.; Radak, B. K.; Allen, B. K.; Sherman, W.; Xu, H. Precise binding free energy calculations for multiple molecules using an optimal measurement network of pairwise differences. _J. Chem. Theory Comput._ 2021, _18_ , 650–663
  * Lu et al. 2021 Lu, C.; Wu, C.; Ghoreishi, D.; Chen, W.; Wang, L.; Damm, W.; Ross, G. A.; Dahlgren, M. K.; Russell, E.; Von Bargen, C. D., et al. OPLS4: Improving force field accuracy on challenging regimes of chemical space. _Journal of chemical theory and computation_ 2021, _17_ , 4291–4300
  * Chung et al. 2023 Chung, M. K.; Miller, R. J.; Novak, B.; Wang, Z.; Ponder, J. W. Accurate Host–Guest Binding Free Energies Using the AMOEBA Polarizable Force Field. _Journal of Chemical Information and Modeling_ 2023, _63_ , 2769–2782
  * 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
  * Shirts et al. 2017 Shirts, M. R.; Klein, C.; Swails, J. M.; Yin, J.; Gilson, M. K.; Mobley, D. L.; Case, D. A.; Zhong, E. D. Lessons learned from comparing molecular dynamics engines on the SAMPL5 dataset. _J Comput Aided Mol Des_ 2017, _31_ , 147–161
  * 70 Landrum, G. RDKit: Open-source cheminformatics.
  * He et al. 2020 He, X.; Man, V. H.; Yang, W.; Lee, T.-S.; Wang, U. A fast and high-quality charge model for the next generation general AMBER force field. _J. Chem. Phys._ 2020, _153_ , 114502
  * Bannwarth et al. 2019 Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-xTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions. _J. Chem. Theory Comput._ 2019, _15_ , 1652––1671
  * Ufimtsev and Martínez 2009 Ufimtsev, I.; Martínez, T. J. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients and First Principles Molecular Dynamics. _J. Chem. Theo. Comp._ 2009, _5_ , 2619
  * Titov et al. 2013 Titov, A. V.; Ufimtsev, I. S.; Luehr, N.; Martínez, T. J. Generating Efficient Quantum Chemistry Codes for Novel Architectures. _J. Chem. Theo. Comp._ 2013, _9_ , 213
  * Jin et al. 2016 Jin, Z.; Yang, C.; Cao, F.; Li, F.; Jing, Z.; Chen, L.; Shen, Z.; Xin, L.; Tong, S.; Sun, H. Hierarchical atom type definitions and extensible all-atom force fields. _J Comput Chem._ 2016, _37_ , 653–664
  * 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
  * Spicher and Grimme 2020 Spicher, S.; Grimme, S. Robust Atomistic Modeling of Materials, Organometallic, and Biochemical Systems. _Angewandte Chemie_ 2020, _59_ , 15665–15673
  * Xu 2019 Xu, H. DiffNet. 2019; https://github.com/forcefield/DiffNet
  * Stumpe et al. 2011 Stumpe, M. C.; Blinov, N.; Wishart, D.; Kovalenko, A.; Pande, V. S. Calculation of local water densities in biological systems: a comparison of molecular dynamics simulations and the 3D-RISM-KH molecular theory of solvation. _J. Phys. Chem. B_ 2011, _115_ , 319–328
  * Case et al. 2019 Case, D.; Aktulga, H.; Belfon, K.; Ben-Shalom, I.; Berryman, J.; Brozell, S.; Cerutti, D.; Cheatham, T.; III,; Cisneros, G.; Cruzeiro, V.; Darden, T.; Duke, R.; Giambasu, G.; Gilson, M.; Gohlke, H.; Goetz, A.; Harris, R.; Izadi, S.; Izmailov, S.; Kasavajhala, K.; Kaymak, M.; King, E.; Kovalenko, A.; Kurtzman, T.; Lee, T.; LeGrand, S.; Li, P.; Lin, C.; Liu, J.; Luchko, T.; Luo, R.; Machado, M.; Man, V.; Manathunga, M.; Merz, K.; Miao, Y.; Mikhailovskii, O.; Monard, G.; Nguyen, H.; OHearn, K.; Onufriev, A.; Pan, F.; Pantano, S.; Qi, R.; Rahnamoun, A.; Roe, D.; Roitberg, A.; Sagui, C.; Schott-Verdugo, S.; Shajan, A.; Shen, J.; Simmerling, C.; Skrynnikov, N.; Smith, J.; Swails, J.; Walker, R.; Wang, J.; Wang, J.; Wei, H.; Wolf, R.; Wu, X.; Xiong, Y.; Xue, Y.; York, D.; Zhao, S.; ; Kollman, P. Amber 2019. 2019; https://ambermd.org/
  * Sindhikara et al. 2012 Sindhikara, D.; Yoshida, N.; Hirata, F. Placevent: an algorithm for prediction of explicit solvent atom distribution-application to HIV-1 protease and F-ATP synthase. _J. Comp. Chem._ 2012, _33_
  * Schindler and Kuhn 2019 Schindler, C.; Kuhn, D. fep-benchmark. 2019; https://github.com/MCompChem/fep-benchmark
  * Schiemann et al. 2016 Schiemann, K.; Mallinger, A.; Wienke, D.; Esdar, C.; Poeschke, O.; Busch, M.; Rohdich, F.; Eccles, S. A.; Schneider, R.; Raynaud, F. I., et al. Discovery of potent and selective CDK8 inhibitors from an HSP90 pharmacophore. _Bioorg. Med. Chem. Lett._ 2016, _26_ , 1443–1451
  * Öhlknecht et al. 2020 Öhlknecht, C.; Lier, B.; Petrov, D.; Fuchs, J.; Oostenbrink, C. Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies. _J. Comput. Chem._ 2020, _41_ , 986–999
  * 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
  * Jorgensen et al. 1983 Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. _J. Chem. Phys._ 1983, _79_ , 926–935
  * ATo 2022 AToM-OpenMM. https://github.com/Gallicchio-Lab/AToM-OpenMM, 2022
  * 88 The ATM Meta Force Plugin for OpenMM. https://github.com/Gallicchio-Lab/openmm-atmmetaforce-plugin, Version 0.3.1
  * 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.; Wiewiora, R. P.; Brooks, B. R.; Pande, V. S. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. _PLoS Comput. Biol._ 2017, _13_ , e1005659
  * 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. _Comput. Phys. Commun._ 2015, _196_ , 236–246
  * 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