# Potential Distribution Theory of Alchemical Transfer

Solmaz Azimi  Department of Chemistry and Biochemistry, Brooklyn College 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  Emilio
Gallicchio  Ph.D. Program in Biochemistry, 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
<EMAIL_ADDRESS>Department of Chemistry and Biochemistry,
Brooklyn College of the City University of New York, New York, NY

###### Abstract

We present an analytical description of the Alchemical Transfer Method (ATM)
for molecular binding using the Potential Distribution Theory (PDT) formalism.
ATM models the binding free energy by mapping the bound and unbound states of
the complex by translating the ligand coordinates. PDT relates the free energy
and the probability densities of the perturbation energy along the alchemical
path to the probability density at the initial state, which is the unbound
state of the complex in the case of a binding process. Hence, the ATM
probability density of the transfer energy at the unbound state is first
related by a convolution operation to the probability densities for coupling
the ligand to the solvent and coupling it to the solvated receptor, for which
analytical descriptions are available with parameters obtained from maximum
likelihood analysis of data from double-decoupling alchemical calculations.
PDT is then used to extend this analytical description along the alchemical
transfer pathway. We tested the theory on the alchemical binding of five
guests to the TEMOA host from the SAMPL8 benchmark set. In each case, the
probability densities of the perturbation energy for transfer along the
alchemical transfer pathway obtained from numerical calculations match those
predicted from the theory and double-decoupling simulations. The work provides
a solid theoretical foundation for alchemical transfer, offers physical
insights on the form of the probability densities observed in alchemical
transfer calculations, and confirms the conceptual and numerical equivalence
between the alchemical transfer and double-decoupling processes.

##  I Introduction

The modeling of free energies is critical to the characterization of
materials, chemical processes, molecular recognition, and many other areas of
investigation where molecular-level insights are pursued.[1] Bridging the gap
between theoretical models and experimental observations is perhaps the most
important role of free energy calculations in modern chemical research. By
comparing calculated free energies to measured values, scientists can refine
their models and gain deeper insights into chemical composition, chemical
interactions, and the dynamical behavior of complex molecular systems.
Alchemical free energy models exploit the state function property of free
energy to compute free energy differences between physical states by
transforming the system’s Hamiltonian along a sequence of non-physical states
that cannot be realized in the laboratory.[2]

The study of protein-ligand binding in drug discovery is one of the most
popular applications of alchemical free energy models.[3, 4, 5, 6, 7, 8] In
the common double-decoupling approach, the binding free energy of a receptor-
ligand complex in solution is obtained as the difference between the free
energy of alchemically turning on the ligand-solvent interactions in solution
and the free energy of turning off the interactions between the bound ligand
and the solvated receptor, each obtained from separate calculations.[9, 10]
The term Double-Decoupling Method (DDM) that refers to this approach stems
from considering each step as “decoupling” the ligand from its environment
(the solvent or the receptor in the solvent) to reach a vacuum state.[11] The
same process can be equivalently described in terms of the reverse steps of
“coupling” the ligand to either the solution or the solvated receptor from
vacuum.

Conceptually, each coupling step of DDM is closely related to the particle
insertion method originally introduced to model solvation free energies.[12,
13, 14, 15] The particle insertion method obtains the solvation free energy by
performing an exponential average of the solute-solvent interaction energies
resulting from random insertions of the solute into an ensemble of pure
solvent configurations. The Potential Distribution Theory (PDT)[12, 16]
applies to particle insertion when considering the solute-solvent energy as
the perturbation energy $u$. The main statement of PDT is that the solvation
free energy and the probability density distribution of the solute-solvent
energy $p_{1}(u)$ in the coupled state are determined by the probability
density distribution $p_{0}(u)$ of the perturbation energy in the uncoupled
ensemble.[16] The PDT relationship extends to the sequence of alchemical
intermediate states corresponding to the progressive introduction of the
solute-solvent interaction by means of a progress parameter $\lambda$,[17]
implying that knowledge of $p_{0}(u)$ determines the perturbation energy
distributions $p_{\lambda}(u)$ and the free energy profile $\Delta G(\lambda)$
along the entire alchemical coupling process.[16, 18]

We used these PDT results to model alchemical binding processes from a
statistical perspective,[19] and to optimize alchemical potential energy
functions.[20] In particular, Kilburg and Gallicchio[21] developed an
analytical model of $p_{0}(u)$ with parameters learned from the distribution
of perturbation energies collected from alchemical molecular simulations of
molecular coupling processes. Using the PDT formula, the model reproduced the
perturbation energies and the free energy profile along the alchemical binding
pathway of host-guest complexes with implicit solvation.[21] The model was
later employed to study alchemically-induced phase transitions, their role in
the rate of convergence of free energy estimates, and in devising optimized
alchemical potential energy functions applicable with explicit solvation
models.[20, 22]

Recently, a direct alchemical transfer route to the free energy of molecular
binding in explicit solution that bypasses the ligand’s vacuum state has been
developed. In the resulting Alchemical Transfer Method (ATM),[23], the
alchemical transformation is encoded in a coordinate transformation that
directly translates the ligand from a position in the solvent bulk into the
receptor binding site. ATM has been extensively validated against host-guest
systems[23, 24, 25] and relative binding free energy protein-ligand
benchmarks,[26, 27, 28] and is considered a viable alchemical approach in
applied research, especially with advanced many-body potential models not yet
supported by standard alchemical models.[29, 30]

This work extends the Potential Distribution Theory (PDT) formalism to
describe alchemical transfer processes. Unlike coupling processes, where
interatomic interactions are created from the uncoupled state, alchemical
transfer involves the gain of ligand-receptor interactions accompanied by the
simultaneous loss of ligand-solvent interactions. Following the PDT formalism,
we approach the problem by seeking the probability density of the perturbation
energy for alchemical transfer at the initial state, where the ligand is in
solution. Because the perturbation energy for alchemical transfer is the
difference between the receptor-ligand and the solvent-ligand interaction
energies, we model its probability distribution as the convolution of the
probability distributions of the two components that we obtain by double-
decoupling alchemical simulations. We show that the PDT applied to the
convolution function successfully reproduces the perturbation energy
probability densities throughout alchemical transfer processes.

We illustrate the PDT theory developed here by applying it to the alchemical
transfer binding of a series of guests to a molecular host. The coupling
processes of the guest to the solution and the receptor from vacuum are
simulated and analyzed in terms of the analytical model of alchemical binding
of Kilburg and Gallicchio.[21] This procedure provides optimized parameters
for the analytical expressions of the $p_{0}(u)$ functions for each coupling
process. We show that the convolution of the analytical models learned from
alchemical coupling simulations matches the probability densities of the
alchemical transfer perturbation energy at the solvated state obtained from
alchemical transfer simulations of binding.

These results demonstrate that the PDT is applicable to alchemical transfer
and that double-decoupling and transfer processes are statistically equivalent
because the probability densities of the perturbation energy of the second can
be determined from the first. More generally, we illustrate that the PDT,
which is traditionally applied to coupling processes, is also suitable to
describe more complex alchemical processes such as alchemical transfer. The
work also illustrates the fundamental concept underlying the PDT that the
alchemical pathways connecting the same endpoints are interrelated because
they originate from the same probability density kernel $p_{0}(u)$. Hence, a
$p_{0}(u)$ model learned from one pathway yields information about all other
alchemical pathways.

##  II Theory

###  II.1 Alchemical Transfer and Double-Decoupling for Modeling Molecular
Binding Equilibria

Consider the standard free energy, $\Delta G^{\circ}_{b}$, of the non-covalent
association equilibrium between receptor R and ligand L to form the receptor-
ligand complex RL

$R_{(aq)}+L_{(aq)}\rightleftharpoons RL_{(aq)}\,,$ (1)

that is related to the binding constant $K_{b}$ through

$\Delta G^{\circ}_{b}=-k_{B}T\ln K_{b},$ (2)

where $k_{B}$ is Boltzmann’s constant and $T$ is the temperature. A
statistical mechanics expression for $K_{b}$ is[11, 31, 18, 32]

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

where $\beta=1/(k_{B}T)$,

$u(x)=U_{1}(x)-U_{0}(x)$ (4)

is the binding energy of the configuration $x$ of the complex defined as the
potential energy of the bound complex, $U_{1}(x)$, relative to the potential
energy, $U_{0}(x)$, of the unbound configuration when receptor and ligand are
uncoupled at large separation, $\langle\ldots\rangle_{0}$ represents the
ensemble average in the uncoupled state, $V_{\rm site}$ is the volume of the
receptor binding site, and $C^{\circ}$ is the one molar standard
concentration. The free energy $\Delta G^{\circ}_{\rm id}=-k_{B}T\ln
C^{\circ}V_{\rm site}/(8\pi^{2})$, is the ideal component of the standard
binding free energy (the value of the standard binding free energy that would
be observed if the ligand did not interact with the receptor) while the term

$\Delta G_{b}=-k_{B}T\ln\langle e^{-\beta u}\rangle_{0}$ (5)

is the excess component of the binding free energy.

The objective of alchemical computational binding free models is to estimate
the quantity $\Delta G_{b}$ as accurately and rapidly as possible, using a
series of non-physical potential energy functions $U_{\lambda}(x)$, where
$0\leq\lambda\leq 1$ is the alchemical progress parameter, interpolating
between the potential energy functions, $U_{0}(x)$ and $U_{1}(x)$, that
describe the unbound and bound states of the complex.[33, 34, 35, 36, 37]
$\Delta G_{b}$ for the association process in Eq. (1) can be modeled
alchemically in explicit solvent by either alchemical transfer[23] or double-
decoupling.[11, 31, 38, 10] In the alchemical transfer method (ATM), the
receptor and ligand are simulated together in a solvent box using a
$\lambda$-dependent hybrid potential of $U_{0}(x)$ and $U_{1}(x)$, where
$U_{1}(x)$ is sampled by translating the ligand into the receptor binding site
from an arbitrary position in the solvent.

The double-decoupling method (DDM), which is more widely used, models the
binding process in two steps: first, L is alchemically transferred from vacuum
to solution (the solvent coupling step), and then, in a separate simulation, L
is alchemically transferred from vacuum to the binding site of R (the receptor
coupling step). The excess binding free energy of the association between L
and R is the difference between the free energies of the solvent and receptor
coupling steps (Figure 1).

Figure 1: The binding free energy of ligand L to receptor R (top panel) in
solution (top panel) is estimated directly by alchemical transfer of the
ligand into the receptor (horizontal arrow), or by double-decoupling, whereby
the ligand is alchemically decoupled from the solvent to vacuum (blue arrow)
followed by alchemical coupling of the ligand to the solvated receptor from
vacuum (pink arrow). The probability densities of the ligand-receptor
interaction energy, $u_{\rm h}$, along the coupling process are derived by the
probability density, $p_{0{\rm h}}(u_{\rm h})$, at the initial decoupled
state. Similarly, the decoupling process is described by the probability
density $\tilde{p}_{1{\rm s}}(u_{\rm s})$ of the loss of ligand-solvent
interaction energy at the solvent-coupled state. The probability density of
the perturbation energy for alchemical transfer, $u_{\rm t}$, is given by the
convolution of these two functions (horizontal purple arrow). The statistic of
the alchemical dissociation process (bottom panel) is obtained analogously
from the probability densities of the loss of ligand-receptor interactions and
the gain of ligand-solvent interactions.

The two coupling steps of DDM can be modeled analytically (see below), taking
advantage of the fact that the binding energy $u$ corresponds, in these cases,
to the interaction energy between the ligand and the environment–the solvent
or the receptor in the solvent.[21] In alchemical transfer,[23, 26] the
perturbation energy is the sum of the loss of ligand-solvent interactions and
gain of ligand-receptor interactions (Figure 1), and the corresponding
perturbation energy distributions are obtained by the convolution of the
distributions of the decoupling and coupling models (see below).

For either coupling or transfer, the ensemble average in Eq. (5) can be
expressed in the form[18]

$\langle e^{-\beta u}\rangle_{0}=\int_{-\infty}^{+\infty}du\ e^{-\beta u}\
p_{0}(u)$ (6)

where $p_{0}(u)$ is the probability density of the binding energy at the
initial state of the alchemical process, in which receptor and ligand are
uncoupled, either when the ligand is placed in the solvent far away from the
receptor as in ATM, or when the ligand is in vacuum as in the two coupling
steps of DDM. As further discussed below, the function $p_{0}(u)$ plays a
central role in this work.

The general expression of the $\lambda$-dependent alchemical potential energy
function employed in this work is

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

where the perturbation energy $u(x)$ is defined by Eq. (4), $W_{\lambda}(u)$
is the alchemical perturbation energy function with the property that
$W_{0}(u)=0$, and $W_{1}(u)=u$ such that $U_{\lambda}(x)$ in Eq. (7) yields
$U_{0}(x)$ and $U_{1}(x)$ at the endpoints of the alchemical transformation.
The standard linear form $W_{\lambda}(u)=\lambda u$ satisfies this criteria,
but non-linear versions can be more efficient in numerical applications.[20]
The specific expression and parameterization of $W_{\lambda}(u)$ employed in
the calculations presented in this work are given in the Computational Details
section.

According to the Potential Distribution Theorem (PDT),[16] the probability
density $p_{0}(u)$ of the perturbation energy at the initial state encodes all
the information of the alchemical process, including the behavior of the
intermediate $\lambda$ states. In particular, the PDT states that the
probability distributions of the perturbation energy at the intermediate
$\lambda$-states are given by[18, 20]

$p_{\lambda}(u)=\frac{e^{-\beta W_{\lambda}(u)}\ p_{0}(u)}{K(\lambda)}\,,$ (8)

where

$K(\lambda)=\int_{-\infty}^{+\infty}e^{-\beta W_{\lambda}(u)}\ p_{0}(u)\,du$
(9)

is the $\lambda$-dependent excess binding constant. In turn, the excess
binding free energy profile is given by $\Delta G_{b}(\lambda)=-k_{B}T\ln
K(\lambda)$, which, at $\lambda=1$, yields the excess binding free energy (Eq.
(5)).[18, 21, 20] Hence, knowledge of $p_{0}(u)$ determines the free energy
profile and the perturbation energy distributions at all intermediate states
along any alchemical pathway joining two given states. This work aims to
construct a model for $p_{0}(u)$ applicable to direct alchemical transfer.

Note that the PDT results summarized by Eqs. (8) and (9) apply to alchemical
models based on energy interpolation[39] whose perturbation energy functions
depend on only one or a few collective variables for which it is meaningful to
consider probability densities as a function of $\lambda$. For example, it
applies to alchemical transfer because its alchemical perturbation energy
function $W_{\lambda}(u)$ in Eq. (7) depends on atomic coordinates only
through the perturbation energy $u(x)=U_{1}(x)-U_{0}(x)$.[23] However, PDT
does not apply to parameter interpolation alchemical models[36, 10, 9] or
models based on $\lambda$-dependent soft-core pair potentials[40] whose
potential energy functions depend directly on atomic coordinates in complex
ways.

###  II.2 Analytical Theory of Alchemical Coupling

Consider $p_{0}(u)$, the probability density of the interaction energy, $u$,
of two molecular entities, such as a solute with a solvent, in the decoupled
state where the two molecular entities are not interacting. Here and
elsewhere, the subscript "$0$" in $p_{0}(u)$ refers to the $\lambda=0$ state
of the alchemical transformation. For an alchemical coupling process, the
initial state is the decoupled state of the system. The coupling energy of a
configuration $x$ of the system is the perturbation energy [Eq. (4)] where
$U_{1}(x)$ is the potential energy of the system when the two molecular
entities are interacting (_coupled state_), and $U_{0}(x)$ is the potential
energy when their interactions are turned off (_decoupled state_). In this
work, $U_{1}(x)$ is obtained from $U_{0}(x)$ by rigidly translating the solute
into the solvent or the receptor from an arbitrary position in vacuum.[19]

Because it does not interact with the environment, a solute explores many
positions and orientations in the decoupled state. Hence, multiple atomic
collisions and overlaps will likely be found when a solute configuration
generated in vacuum is transferred into the solvent. As a result of these
collisions, perturbation energies $u$ corresponding to these configurations
are likely to be large and positive.[14] In addition to these short-ranged
repulsive interactions, solute-solvent interactions are characterized by long-
ranged, slowly varying, and mostly favorable electrostatic and dispersion
interactions.[41, 42]

Kilburg and Gallicchio[21] exploited the distinct nature of collisional and
long-range interactions to develop an analytical model for $p_{0}(u)$. They
expressed the total interaction energy $u$ as the sum of (i) a collisional
interaction energy ($u_{c})$, representing unfavorable short-ranged, repulsive
interactions, and (ii) a background interaction energy ($u_{b}$), representing
mostly favorable, long-ranged, attractive interactions. They reasoned that the
background interaction energy should follow linear response and central limit
statistics because many individual interatomic interactions contribute to
it.[43, 44, 45, 46, 47, 48, 49, 50, 51] Conversely, the collisional
interaction energy is dominated by the closest, most repulsive pairwise atomic
interaction and is thus expected to follow max statistics.[52] Starting with a
Lennard-Jones pair-potential description of collisions,[20] Kilburg and
Gallicchio developed an analytical statistical model of the collision energy
and expressed the probability density in the decoupled ensemble, $p_{0}(u)$,
of the total interaction energy, $u=u_{c}+u_{b}$, as the convolution of the
collisional and background statistical models.

Specifically, $p_{0}(u)$ is written as[21, 20]

$p_{0}(u)=b\mathcal{N}(u)+(1-b)C(u)$ (10)

where

$\mathcal{N}(u_{b})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\
e^{\frac{-(u_{b}-\bar{u}_{0})^{2}}{2\sigma^{2}}}.$ (11)

is the normal distribution with mean $\bar{u}_{0}$ and standard deviation
$\sigma$, $b$ is the probability that no collisions occur in the decoupled
ensemble, and

$C(u)=\int_{-\infty}^{+\infty}du^{\prime}\
\mathcal{N}(u^{\prime})\mathcal{F}(u-u^{\prime})=(\mathcal{N}\circledast\mathcal{F})(u).$
(12)

is the convolution of the probability density, ${\mathcal{F}}(u_{c})$, of the
collisional interaction energy with the normal distribution of Eq. (11) that
represents the background interaction energy. Here and elsewhere in this work,
the convolution operation ($\circledast$) arises whenever we interrogate the
statistical behavior of a random energy variable $u$, which is the sum of two
random variables whose statistics are known or assumed.

The analytical expression of ${\mathcal{F}}(u_{c})$ is[20]

${\mathcal{F}}(u_{c})=n_{l}\left[1-\frac{(1+x_{c})^{1/2}}{(1+x)^{1/2}}\right]^{n_{l}-1}\frac{H(u)}{4\epsilon}\frac{(1+x_{c})^{1/2}}{x(1+x)^{3/2}}$
(13)

in which, $x=\sqrt{1+u_{c}/\epsilon+{\tilde{u}}/{\epsilon}}$,
$x_{c}=\sqrt{1+{\tilde{u}}/{\epsilon}}$, and $H(\cdot)$ is the Heaviside step
function. The parameters of the collisional model have the following physical
interpretations. The parameter $\tilde{u}$ represents the interaction energy
above which the solute-solvent interaction energy follows max statistics.
$n_{l}$, which scales as the solute size, describes the number of
statistically independent atom groups of the solute. Finally, $\epsilon$ is an
effective Lennard-Jones potential energy prefactor that describes the rate of
increase of the collisional energy as two atoms approach each other. The
collisional parameters, together with the linear response parameters
$\bar{u}_{0}$ and $\sigma$, specify the analytical model of alchemical
coupling of Eq. (10). As described below, the parameters of the analytical
model for $p_{0}(u)$ are obtained by maximum likelihood analysis of simulation
data.

To model the alchemical coupling of flexible polyatomic ligands that can adopt
more than one conformation, in this work, we express $p_{0}(u)$ as the
weighted average of modes described by the model above:[22]

$p_{0}(u)=\sum^{m}_{i}w_{i}p_{0,i}(u)$ (14)

where $p_{0,i}(u)$, with parameters $p_{b,i}$, $\bar{u}_{0,i}$, etc., is the
$p_{0}(u)$ model specific for mode $i$, and the weight parameters $w_{i}$
represents the population of each mode in the decoupled ensemble.

###  II.3 Analytical Model of Alchemical Transfer

In this section we derive an analytical model of alchemical transfer using the
alchemical coupling formalism described above. Consider the thermodynamic
scheme in Fig. 1, in which a guest is transferred from the solution to the
binding site of a molecular host in order to alchemically estimate the binding
free energy of the host-guest complex. Even though an analytical model of
direct alchemical transfer is not available, the transfer process can be
represented by sequential decoupling and coupling processes covered by the
analytical formalism. The change in potential energy, $u_{\rm t}$, for
transferring the ligand from the solvent to the host, is the sum of the loss
of ligand-solvent interaction energy, $u_{\rm s}$, and the gain, $u_{\rm h}$,
of the interaction energy between the ligand and the host (including the
surrounding solvent). Hence, the probability density of $u_{\rm t}$ can be
expressed as the convolution of the probability densities of $u_{\rm s}$ and
$u_{\rm h}$ collected in the initial ensemble where receptor and ligand are
dissociated in solution.

Denoting $p_{1\rm s}(u_{\rm s})$ as the probability density of the ligand-
solvent interaction energy of the ligand in solution (the end-state of the
solvation alchemical process), and $p_{0\rm h}(u_{\rm h})$ as that of the
ligand-receptor coupling energy in the initial state, the probability density
of the perturbation energy for solvent to receptor transfer, $u_{\rm t}=u_{\rm
h}-u_{\rm s}$, is given by the convolution

$p^{+}_{0\rm t}(u_{\rm t})=(p_{0\rm h}\circledast{\tilde{p}}_{1\rm s})(u_{\rm
t})$ (15)

where ${\tilde{p}}_{1\rm s}(u_{\rm s})=p_{1\rm s}(-u_{\rm s})$, is the
probability density of the solute-solvent interaction energy loss in the
solvated coupled ensemble (the initial state of the alchemical process, which
is denoted by the blue arrow in the bottom panel of Fig. 1). (The meaning of
the superscripted plus symbol in Eq. (15) is specified below.)

Because in the fully solvated state the ligand does not experience collisions
with the solvent molecules, the solute-solvent interaction energy is expected
to follow linear response and the corresponding probability density can be
described by a normal distribution

$\tilde{p}_{1\rm s}(u_{\rm s})=\tilde{\mathcal{N}}_{1\rm s}(u_{\rm s})$ (16)

with mean $-\bar{u}_{\rm s}$ and standard deviation $\sigma_{\rm s}$, where
$\bar{u}_{\rm s}$ and $\sigma_{\rm s}$ are the mean and standard deviation of
the ligand-solvent interaction energy in the solvated coupled ensemble,
respectively. Generally, as in Eq. (14), $\tilde{p}_{1\rm s}(u_{\rm s})$ is
represented by a weighted sum of normal distributions. However, we assume one
solvation mode for now to keep the notation simple. Furthermore, the
assumption of linear response is expected to break down when, during the
alchemical decoupling process, the interactions of the ligand with the solvent
are weakened to the point that atomic clashes occur with significant
probability. Nevertheless, following the Alchemical Transfer Method (ATM)
protocol,[23] and as further discussed below, this model is applied only up to
the alchemical intermediate state at $\lambda=1/2$ before atomic clashes are
observed.

The $p_{0\rm h}(u_{\rm h})$ function of the second alchemical process
corresponding to the gain of interactions between the solvated ligand and the
receptor is modeled in this work as that of the coupling processes of the
guest in vacuum to the host (see Fig. 1 top, pink arrow). This approximation
is justified by the small influence of the solvent on the distribution of the
internal degrees of freedom of the rigid guests considered here. Hence
$p_{0\rm h}(u_{\rm h})$ is represented by the analytical model of alchemical
coupling described in the previous section. Here too, we generally consider
multiple binding modes in practice, each of them described by the form in Eq.
(10)

$p_{0\rm h}(u_{\rm h})=b\mathcal{N}_{0\rm h}(u_{\rm h})+(1-b)C_{0\rm h}(u_{\rm
h})$ (17)

where $\mathcal{N}_{0\rm h}(u_{\rm h})$ is a normal distribution with mean
$\bar{u}_{\rm h}$ and standard deviation $\sigma_{h}$, and

$C_{0\rm h}(u_{\rm h})=(\mathcal{N}_{0\rm h}\circledast{\mathcal{F}}_{\rm
h})(u_{\rm h})$ (18)

is the convolution of the linear response and collisional models for coupling
the guest to the host.

By inserting Eqs. (16), (17), and (18) in Eq. (15) and using the linearity of
the convolution operator, we finally obtain the following model for the
probability density of alchemical transfer

$p^{+}_{0\rm t}(u_{\rm t})=b\mathcal{N}^{+}_{0\rm t}(u_{\rm
t})+(1-b)C^{+}_{0\rm t}(u_{\rm t})$ (19)

where

$\mathcal{N}^{+}_{0\rm t}(u_{\rm t})=(\tilde{\mathcal{N}}_{1\rm
s}\circledast\mathcal{N}_{0\rm h})(u_{\rm t})$ (20)

is the convolution of the normal distributions for desolvation and receptor
coupling, which is itself a normal distribution with mean $\bar{u}_{0\rm
t}=\bar{u}_{0\rm h}-\bar{u}_{1\rm s}$ and standard deviation $\sigma_{0\rm
t}=\sqrt{\sigma^{2}_{0\rm h}+\sigma^{2}_{1\rm s}}$, and

$C^{+}_{0\rm t}(u_{\rm t})=(\mathcal{N}^{+}_{0\rm
t}\circledast{\mathcal{F}}_{\rm h})(u_{\rm t})$ (21)

is the convolution of the collisional distribution for the coupling to the
host with the normal distribution in Eq. (20).

Eqs. (19)–(21) establish that, under the present assumptions, the $p_{0}$
distribution for the transfer process has the same form as that of a coupling
process with parameters determined by specific combinations of those of the
coupling process. In particular, the mean perturbation energy linear response
parameter of the transfer model ($\bar{u}_{0\rm t}$) is the difference between
the corresponding parameter of the receptor coupling model and the average
solute-solvent interaction energy in the fully solvated state. The variance
linear response parameter ($\sigma^{2}_{0\rm t}$) is the sum of the variances
of the desolvation and receptor coupling processes. The parameters of the
collisional model ($b$, $\epsilon$, $\tilde{u}$, and $n_{l}$) are inherited
directly from the receptor coupling collisional model.

The expressions above have been derived for the simplest case of one solvation
coupling mode and one receptor coupling mode. In general, each possible pair
of modes with weights $w_{i\rm h}$ and $w_{j\rm s}$ of $m_{\rm h}$ receptor
coupling modes and $m_{\rm s}$ solvation coupling modes, respectively, combine
in the manner above to yield $m_{h}m_{s}$ transfer modes each with weight
$w_{ij\rm t}=w_{i\rm h}w_{j\rm s}$.

The formalism so far describes the transfer process in the binding direction
illustrated by the upper panel of Fig. (1) and denoted by a ‘+’ superscript in
the expressions above. A similar prescription applies to the transfer
unbinding process, where the bound ligand decouples from the receptor and
couples to the solution. The corresponding analytical model for the
$p^{-}_{0\rm t}(u_{\rm t})$ unbinding distribution has the same form (19) as
for the transfer model for binding but with collisional parameters obtained
from the model of coupling to the solvent and linear response parameters
obtained from combining those of decoupling from the receptor and coupling to
the solvent.

##  III Methods

###  III.1 The Alchemical Transfer Method

The double-decoupling and alchemical transfer binding free energy calculations
reported in this work have been conducted using the Alchemical Transfer Method
(ATM).[23, 24, 26, 25, 27, 28, 30] Unlike alchemical approaches that modify
the parameters of the energy function,[36, 5] ATM relates the potential energy
function of the final state ($U_{1}(x)$) to that of the initial state
($U_{0}(x)$) by a coordinate transformation.

Specifically, for the case of the solvation process of a ligand L from vacuum,
denoted by $U_{0}(x_{\rm S},x_{\rm L})$ the potential energy of the system
when the ligand’s coordinates $x_{\rm L}$ are such that the ligand is placed
in vacuum far away from the solvent, whose molecules have coordinates $x_{\rm
S}$. The potential energy of the system when the ligand is placed in the
solvent is expressed in terms of $U_{0}(x_{\rm S},\ x_{\rm L})$ as
$U_{1}(x_{\rm S},\ x_{\rm L})=U_{0}(x_{\rm S},\ x_{\rm L}+h)$, where $h$ is a
displacement vector that brings the ligand from its position in vacuum to the
corresponding position in the solvent. The binding process between a receptor
R and a ligand L is described in a similar way by using a displacement vector
that transfers the ligand from a position into the solvent to the binding site
of the receptor. This formalism represents the unbound and bound states of the
system by a single set of degrees of freedom and to define the perturbation
energy as the difference in the system’s potential energy before and after the
application of the ligand displacement. For example, for the binding process
we define the perturbation energy as

$u^{+}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})=U_{0}(x_{\rm R},\ x_{\rm S},\ x_{\rm
L}+h)-U_{0}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})$ (22)

and the corresponding alchemical potential energy function as

$U^{+}_{\lambda}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})=U_{0}(x_{\rm R},\ x_{\rm
S},\ x_{\rm L})+W_{\lambda}[u^{+}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})].$ (23)

The alchemical potential energy function (23) cannot cover the entire
alchemical binding pathway when the solvent is represented explicitly.[23]
Instead, the alchemical process is decomposed into two legs. In the first leg,
the system is taken from the unbound state to an alchemical intermediate state
(at $\lambda=1/2$, typically) using the potential (23). The second leg
proceeds in the unbinding direction starting from the bound state at
$\lambda=1$ until it reaches the same alchemical intermediate, using the
alchemical potential energy function

$U^{-}_{1-\lambda}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})=U_{0}(x_{\rm R},\ x_{\rm
S},\ x_{\rm L}+h)+W_{\lambda}[u^{-}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})]$ (24)

where

$u^{-}(x_{\rm R},\ x_{\rm S},\ x_{\rm L})=U_{0}(x_{\rm R},\ x_{\rm S},\ x_{\rm
L})-U_{0}(x_{\rm R},\ x_{\rm S},\ x_{\rm L}+h)$ (25)

is the perturbation energy corresponding to the second leg.

To compute the binding free energy, the reversible work values ($\Delta G^{+}$
and $\Delta G^{-}$) along each alchemical leg is calculated by thermodynamic
reweighting[53, 54] and the excess binding free energy is estimated from their
difference:

$\Delta G_{b}=\Delta G^{+}-\Delta G^{-}\,.$ (26)

##  IV Computational Details

###  IV.1 Benchmark Systems

We tested the proposed PDT theory for alchemical transfer on a subset of the
SAMPL8 host-guest benchmark (Figure 2). Specifically, we considered the
binding of the five small guest compounds, named G1 to G5, plus water to the
tetramethyl octa-acid (TEMOA) host.[55] Experimentally, the SAMPL8
measurements at a pH where, with the exception of G2, the guests are expected
to be deprotonated in solution and when bound to the host.[24] However, due to
the difficulties of applying the double-decoupling process to ionized
species,[56] in this work, we opted to carry out the numerical test of our
theory on the neutral forms of the guests.

(a)

(b)

Figure 2: The SAMPL8 benchmark set considered in this work. Top: tetramethyl
octa acid host (TEMOA) shown in licorice representation. Gray corresponds to
carbon atoms, red to oxygen, and white to hydrogen. Middle: the five guests of
the SAMPL8 benchmark set bound to TEMOA. Bottom: the chemical structures of
the five guests, G1 to G5, shown in ball-and-stick (CPK) representation. Gray
corresponds to carbon atoms, red to oxygen, brown to bromine, and white to
hydrogen. The color of the labels corresponds to the color of the guests in
the middle panel.

###  IV.2 System Setup and Simulation Settings

In Eq. (7), we employ the softplus alchemical potential energy function[20,
22]

$W_{\lambda}[u_{\rm
sc}(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)$ (27)

where the soft-core perturbation energy function $u_{\rm sc}(u)$ is defined
as:

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

with

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

and

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

The soft-core perturbation energy function $u_{\rm sc}(u)$ is a monotonically
increasing function of $u(x)=U_{1}(x)-U_{0}(x)$, designed to smoothly cap
large values of the perturbation energy encountered along the alchemical
transformation to the maximum value without affecting the end states.[22, 23,
26] In this work, we set $u_{c}=0$ kcal/mol and $u_{\rm max}=50$ kcal/mol for
the coupling and decoupling calculations. The transfer calculations for water-
in-water and TEMOA-H2O employed the parameters $u_{c}=0$ kcal/mol and $u_{\rm
max}=50$ kcal/mol, whereas the transfer calculations for the guests to the
host employed $u_{c}=100$ kcal/mol and $u_{\rm max}=200$ kcal/mol. The $a$
parameter of the soft-core function was set to $0.0625$ in all cases.

The parameters ($\lambda_{1}$, $\lambda_{2}$, $\alpha$, and $u_{0}$) of the
softplus alchemical perturbation energy function [Eq. (27)] are functions of
$\lambda$ and vary along the alchemical transformation according to a set
schedule (see the Appendix). The softplus alchemical perturbation function
reduces to the standard linear form $\lambda u_{\rm sc}$ when the schedule is
such that $\lambda_{1}=\lambda_{2}=\lambda$. The TEMOA-H2O coupling and
decoupling calculations employed a linear alchemical schedule. All other
alchemical transformations employed non-linear schedules to accelerate
conformational mixing. The analytical models’ parameters for $p_{0}(u)$ are
independent of the alchemical schedules.

The host-guest systems were prepared from the original MOL2 files provided by
the SAMPL8 organizers at
https://github.com/samplchallenges/SAMPL8/tree/master/host_guest/GDCC. All
five guests were protonated and manually placed into the inner cavity of the
TEMOA host with their polar ends directing out of the cavity (Figure 2) using
Maestro (Schrödinger, Inc.). Force-field parameter assignments with the
GAFF1.8/AM1-BCC force field and solvation of the systems were performed using
AmberTools 19 and the LEaP program.

The simulations were conducted in a water slab (Figure 3) of approximate
dimensions $40\times 60\times 42$ Å3 embedded in a $40\times 60\times 142$ Å3
periodic simulation box. The resulting system contains layers of water slabs
of 42 Å thickness separated by 100 Å-thick vacuum regions along the $z$
direction. The evaporation of water molecules from the slab was prevented by
imposing a flat-bottom harmonic restraint to the oxygen atoms of the water
molecules along the $z$-direction with a force constant of 1.9 kcal/mol Å2 and
a tolerance of 21 Å from the center of the slab. The water solvent in the slab
was minimized and thermalized at 300 K. In the alchemical decoupling
calculations, the guest was transferred from the water slab or the host
binding site to a position in the vacuum region displaced by 70 Å along the
$z$-direction. The alchemical transfer calculations employed a 30 Å
displacement vector along the $x$-direction parallel to the slab to bring the
guest from the solvent to the binding site of the host.

Figure 3: The water slab system that all calculations were conducted in. The
system shown here is TEMOA and G2. Structures are styled in licrorice, in
which cyan are carbon atoms, red oxygen, white hydrogen, and light pink
bromine.

The heavy atoms of the lower cup of the TEMOA host, which were the first 40
atoms of the host as listed in the SAMPL8-provided files, were restrained by a
flat-bottom harmonic potential with a force constant of 25 kcal/mol Å2 and a
tolerance of 0.5 Å.[24] A flat-bottom harmonic potential between the centers
of mass of the host and the guest with a force constant of 25 kcal/mol Å2 was
applied for a distance greater than $4.5$ Å to define the binding site region
($V_{\rm site}$). The corresponding ideal binding free energy term[18]

$\Delta G^{\circ}_{\rm id}=-k_{B}T\ln C^{\circ}V_{\rm site}$ (31)

is equal in this case to $0.87$ kcal/mol.

Alchemical transfer calculations started from a configuration at the
alchemical intermediate obtained by slow annealing to $\lambda=1/2$ in 250 ps
starting from the initial state at $\lambda=0$ prepared by conventional energy
minimization, thermalization, and relaxation steps. The annealing step
establishes a suitable initial configuration of the system at $\lambda=1/2$
without severe unfavorable repulsive interactions that would otherwise hamper
molecular dynamics simulations.

Asynchronous Hamiltonian molecular dynamics replica exchange[57, 58] with a 2
fs time-step and 22 replicas in $\lambda$-space was employed for production
conformational sampling using the AToM-OpenMM software package.[59] Replicas
were cycled in and out of the GPU devices every 40 ps. Perturbation energy
samples and trajectory frames were saved with the same frequency. Each replica
was simulated for 60 ns. Free energies and their corresponding uncertainties
were estimated using UWHAM thermodynamic reweighting[60] after discarding 1/3
of the initial trajectory.

###  IV.3 Parameter Optimization Protocol

We obtained optimized parameters of the analytical model for $p_{0}(u)$ by
means of the maximum-likelihood analysis of the distributions of perturbation
energy samples of all $\lambda$-states from the corresponding alchemical
molecular dynamics calculations.[21] Specifically, because the calculations
provide samples of the soft-core perturbation energies $u_{\rm sc}$, we
consider the probability density of the soft-core perturbation energy given by

$p_{0}(u_{\rm sc})=p_{0}(u)/u_{\rm sc}^{\prime}(u)$ (32)

where $u=u(u_{\rm sc})$ is the value of the inverse of the soft-core function
at $u_{\rm sc}$ [Eq. (28)] and $u_{\rm sc}^{\prime}(u)\geq 0$ is the
derivative of the soft-core function. Even though they are different
functions, to simplify the notation, here we use the same symbol, $p_{0}$, for
the probability density functions of the original, $p_{0}(u)$ and soft-core,
$p_{0}(u_{\rm sc})$, perturbation energies and use their arguments to
distinguish them.

The cost function is then expressed in terms of the likelihood function
$\mathcal{L}(\theta)$ as

$-\log\mathcal{L}(\theta)=-\sum_{i}\log p_{\lambda_{i}}(u_{{\rm
sc},i}|\theta)$ (33)

where $\theta$ represents the collection of the parameters of the model that
we seek to optimize and, $p_{\lambda_{i}}(u_{{\rm sc}}|\theta)$ is the
analytical expression of the probability density of $u_{\rm sc}$ at the
alchemical state at $\lambda=\lambda_{i}$. The latter is obtained from Eq. (8)
using the perturbation energy function $W_{\lambda}[u_{\rm sc}(u)]$. The sum
in Eq. (33) runs over the samples from the molecular simulations at all
$\lambda$-states, where $u_{{\rm sc},i}$ denotes the soft-core perturbation
energy of the sample and $\lambda_{i}$ the value of $\lambda$ of the
alchemical states from where the sample was collected.[21]

Initial guesses for the parameters of the analytical model for $p_{0}(u)$ were
derived from the shapes of kernel density estimates of the $\log p_{0}(u_{\rm
sc})$ function and of the corresponding $\lambda$-function[20, 22]

$\lambda_{0}(u_{\rm sc})=k_{B}T\frac{d\log p_{0}(u_{\rm sc})}{du_{\rm sc}}$
(34)

(see Figure 4). Briefly, the UWHAM statistical inference analysis provides a
statistical weight $W_{0,i}$ to each sample $i$ that represents the
probability of observing it at $\lambda=0$, even though it might have been
collected at some other $\lambda$-state. We obtained $p_{0}(u_{\rm sc})$ using
the weighted Gaussian kernel estimate

$p_{0}(u_{\rm sc})=\sum_{i}W_{0,i}\mathcal{N}(u_{\rm sc}|u_{{\rm
sc},i},\sigma)$ (35)

where $\mathcal{N}(u|{\bar{u}},\sigma)$ is the normal distribution with mean
${\bar{u}}$ and standard deviation $\sigma$. The kernel estimate for the
$\lambda$-function was obtained similarly using the derivative of the normal
distribution as a kernel function. In this work, we performed Gaussian kernel
estimates with $\sigma=1$ kcal/mol.

The $\log p_{0}(u_{\rm sc})$ and $\lambda_{0}(u_{\rm sc})$ functions have
characteristic shapes (see Figure 4) that provide information about the number
of modes and their parameters. For example, following linear response, $\log
p_{0}(u_{\rm sc})$ tends to vary quadratically at low energies when the system
is nearly coupled. Similarly, the $\lambda$-function often varies linearly in
this regime, and the onset and slope of the curve provide estimates for the
analytical model’s mean and standard deviation parameters. Conversely,
deviations from quadratic and linear behaviors can be ascribed to
contributions from multiple modes and can provide information about their
parameters and relative statistical weight. The relationship between the
collisional parameters ($b$, $\epsilon$, $\tilde{u}$, and $n_{l}$) and the
shapes of the $\log p_{0}(u_{\rm sc})$ and $\lambda_{0}(u_{\rm sc})$ is less
obvious. Nevertheless, it was helpful to find initial guesses by studying the
effect of varying the parameters on the agreement between the calculated and
analytical curves. We used the Mathematica program (Wolfram, Inc) for this
purpose.

After finding initial guesses for the number of modes and their parameters,
the parameters of the analytical model were refined by minimizing the cost
function (33) using a protocol implemented in TensorFlow available at
https://github.com/Gallicchio-Lab/femodel-tf-optimizer.[21] The procedure
involves the numerical integration steps to evaluate, for example, Eq. (12).
These were performed by Gauss-Hermite quadrature using 19 nodes. Parameter
optimization was performed on the same set of perturbation energy samples used
for the estimation of free energies.

(a) Coupling

(b) Solvation

Figure 4: The $\log p_{0{\rm h}}(u_{\rm h,sc})$ and $\log p_{0{\rm s}}(u_{\rm
s,sc})$ functions (left) and the $\lambda_{0{\rm h}}(u_{\rm h,sc})$ and
$\lambda_{0{\rm s}}(u_{\rm s,sc})$ $\lambda$-functions (right) of the soft-
core perturbation energies for (a) the coupling of water to TEMOA and (b) the
solvation of water in water. The yellow curves are Gaussian kernel estimates
from the samples of the alchemical molecular simulations. The red and blue
curves are from the corresponding optimized analytical models.

##  V Results

Here, we present a series of results to illustrate that the statistical
behavior of alchemical transfer processes can be predicted by analyzing
double-decoupling processes. The results are arranged as follows. First, a
simple proof-of-principle case study is demonstrated for the transfer of a
water molecule from one location in the solvent to another. We then present
the results of double-decoupling alchemical calculations for the SAMPL8 series
of TEMOA host-guest complexes and the analytical models of the alchemical
probability densities. In each case, we show that the free energies and
analytical models for the corresponding alchemical transfer processes obtained
by the convolution of the decoupling processes agree with the alchemical
transfer simulations’ results.

###  V.1 Transfer of Water in Water

The transfer free energy of a water molecule can be computed by first
transferring the water molecule from solution to vacuum and then coupling it
back to solvent (double-decoupling process). The transfer free energy can also
be obtained by directly moving the water molecule from one location in the
solvent to another (alchemical transfer process). Because the solvent is
uniform, the initial and final states of the transfer process are
thermodynamically equivalent, and the corresponding free energy is zero. Due
to this symmetry, the decoupling legs from the initial and final states and
the alchemical transfer’s forward and reverse legs are also equivalent. Hence,
we will present only one example of each.

This case study serves as a proof-of-principle test for investigating the
hypothesis posed in this work. The alchemical coupling simulation yields the
parameters of the analytical model for $p_{0{\rm s}}(u_{\rm s}$), the
probability density of the interaction energy between the distinguished water
molecule and the rest of the solvent. The optimized parameters of the
analytical model for water coupling are listed in Table 1 under the “H2O
coupling” header. They indicate that the coupling of water to water is well
described by one mode, that the mean and standard deviations of the background
interaction energy in the decoupled ensemble are $2.41$ and $3.46$ kcal/mol,
respectively (the $\bar{u}_{b}$ and $\sigma_{b}$ parameters), that the
probability of finding a configuration free of clashes is $5.77\times 10^{-3}$
(the $b$ parameter), that the effective Lennard-Jones $\epsilon$ parameter for
collisions is $3.9$ kcal/mol, that the minimum collision energy is $3.9$
kcal/mol (the $\tilde{u}$ parameter), and that the effective average number of
colliding atoms is $2.50$ (the $n_{l}$ parameter).

When it is transformed according to the Potential Distribution Theorem
prescription [Eq. (8)], the $p_{0{\rm s}}(u_{\rm s}$) model yields the
probability densities $p_{\lambda{\rm s}}(u_{\rm s}$) of the water-solvent
interaction energy $u_{s}$ as a function of the coupling parameter $\lambda$
(Figure 6b). The weakly coupled states near $\lambda=0$ are characterized by
wide distributions with long tails at high interaction energies characteristic
of frequent and severe atomic clashes. The interaction energies plotted in
Figure 6b and elsewhere are damped down by the soft-core function [Eq. (28)].
The raw interaction energies of these states can be orders of magnitude
greater and the corresponding distributions stretch towards large values. As
the coupling increases, the distributions shift to lower interaction energies.
The features at interaction energies just greater than zero are artifacts due
to the soft-core function that terminates there. For $\lambda$ values above
approximately $1/2$, the probability densities at $\lambda$ values above about
$1/2$ assume a Gaussian shape and shift towards lower energies proportionally
to $\lambda$ at constant width as expected from linear response.[50, 43] until
they reach the fully coupled state at $\lambda=1$ described by $p_{1{\rm
s}}(u_{\rm s}$).

Next, we take the convolution of $\tilde{p}_{{1}\rm s}(u_{\rm s})=p_{1{\rm
s}}(-u_{\rm s})$ and $p_{0{\rm s}}(u_{\rm s})$ that yields, according to Eq.
(19), a model for the probability density $p_{0{\rm t}}(u_{\rm t})$ for the
alchemical transfer of the water molecule from one place in the solvent to
another. The parameters of the analytical model of $p_{0{\rm t}}(u_{\rm t})$
are listed in Table 1 under the “H2O transfer” header. Furthermore,
application of Eq. (8) to $p_{0{\rm t}}(u_{\rm t})$ yields analytical
predictions of the probability densities of the transfer perturbation energies
at all $\lambda$ values along the alchemical path. As shown in Figure 5, there
is an excellent agreement between the analytical predictions of the
$p_{\lambda{\rm t}}(u_{\rm t})$ alchemical transfer probability densities and
the results of molecular simulation.

Figure 5: The probability densities $p_{\lambda{\rm t}}(u_{\rm t,sc})$ of the sot-core perturbation energy of alchemical transfer collected from simulations of water transfer in water (dots) compared with the analytical probability densities (lines) predicted from the model of alchemical coupling $p_{0{\rm s}}(u_{\rm s})$. $\lambda=0.1$ (yellow), $\lambda=0.2$ (green), $\lambda=0.3$ (blue), $\lambda=0.4$ (purple), and $\lambda=0.5$ (pink). Table 1: The parameters of the analytical model of alchemical coupling and alchemical transfer for a water molecule in explicit solvent. | $w_{i}$ | $b$ |  $\bar{u}_{b}$111kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
H2O coupling222from Table 4
mode 1 | $1.0$ | $5.77\times 10^{-3}$ | $2.41$ | $3.46$ | $3.9$ | $3.9$ | $2.50$
H2O transfer333Leg 1 and Leg 2
mode 1 | $1.0$ | $5.82\times 10^{-3}$ | $20.1$ | $4.89$ | $3.9$ | $3.9$ | $2.50$

###  V.2 Binding Free Energy Estimates of the Host-Guest Complexes

The results of the double-decoupling (DDM) and alchemical transfer (ATM)
binding free energy calculations for the complexes of TEMOA with water and the
five SAMPL8 guests are presented in Table 2. In general, the free energy
estimates between the two methods are within statistical uncertainty. In
particular, the two methods clearly agree for TEMOA-H2O, TEMOA-G2, TEMOA-G4,
and TEMOA-G5. The difference in the estimates for TEMOA-G1 and TEMOA-G3 is
about 1 kcal/mol, which is just outside the level of confidence and could
indicate slow convergence or a small bias. The statistical uncertainties of
the DDM and ATM estimates are similar, considering that the ATM calculations
are only half as expensive as the combinations of the two coupling steps of
the DDM calculations.

G2 and G4 are the most favorable binders of the five guests, a rank
consistently predicted by both DDM and ATM. The least favorable binder in the
DDM calculations is G3, whereas ATM designates G1 as the weakest binder. Of
the five guests, G1 is the bulkiest, containing 14 heavy atoms, and G3 is the
only non-planar and non-aromatic guest. Evidently, binding either of these
molecules to the cavity of TEMOA is a more challenging transformation than
binding G2, G4, and G5 guests, all of which contain a single aromatic ring.
The double-decoupling free energies indicate that the higher affinities of the
G2 and G4 guests are due to their stronger interactions with the host (the
“Host Coupling” free energies in Table 2) relative to the other guests, while
the variation of the hydration free energies (the “Solvent Coupling” free
energies in Table 2), which oppose binding, are comparatively smaller. The
computed binding free energies of the protonated guests presented here cannot
be compared directly with the modeled and experimental binding affinities of
the ionized forms of the guests reported earlier.[24, 61, 62] However, G2 and
G4 are the strong binders in both cases, probably reflecting small variations
of the guests’ pKa’s and the resulting ionization penalties.[25]

Water’s excess binding free energy, $\Delta
G^{\ast}_{b,w}=-k_{B}T\ln\rho/\rho_{0}$, to the TEMOA cavity measures the
average water density within the binding site volume relative to bulk.[63]
Taking into account the volume of the spherical binding site region (see
Computational Details), the transfer free energy calculated for water of
approximately $0.8$ kcal/mol (Table 2) is consistent with about three water
molecules present within TEMOA’s cavity displaced upon binding of the guest
molecule. This number is consistent with direct observations.[24]

Table 2: Excess free energy estimates of the case studies in this work obtained by two methods: coupling/decoupling (DDM) and direct transfer (ATM). The difference between the host coupling and the solvent coupling free energies yields the DDM binding free energy of the designated ligand to host TEMOA. Uncertainties reported are twice the standard deviation. | Host Coupling444kcal/mol | Solvent Couplinga | DDMa | ATMa
---|---|---|---|---
TEMOA-H2O | $-3.87\pm 0.07$ | $-4.50\pm 0.1$ | $0.628\pm 0.1$ | $0.845\pm 0.3$
TEMOA-G$1$ | $-14.2\pm 0.1$ | $-2.51\pm 0.1$ | $-11.7\pm 0.2$ | $-10.5\pm 0.2$
TEMOA-G$2$ | $-17.0\pm 0.09$ | $-4.32\pm 0.1$ | $-12.7\pm 0.1$ | $-12.7\pm 0.2$
TEMOA-G$3$ | $-14.3\pm 0.1$ | $-3.50\pm 0.1$ | $-10.8\pm 0.2$ | $-11.6\pm 0.3$
TEMOA-G$4$ | $-17.6\pm 0.1$ | $-3.54\pm 0.1$ | $-14.1\pm 0.2$ | $-13.9\pm 0.3$
TEMOA-G$5$ | $-15.1\pm 0.1$ | $-3.19\pm 0.1$ | $-11.9\pm 0.2$ | $-11.7\pm 0.3$

###  V.3 Analytical Models for Coupling to the Host

The binding of one water molecule to TEMOA from vacuum can be described
analytically by a probability density function $p_{0h}(u_{h})$ with two modes
(Table 1) of nearly equal statistical weight. The second mode corresponds to
configurations that are more likely to clash with the host’s atoms or water
molecules in the cavity (the smaller $b$ parameter in Table 3). The two modes
probably reflect the position of the uncoupled water molecule either at the
center or at the rim of the binding site volume–where it is more likely to
find the host’s atoms–or configurations that, by chance, have fewer or more
water molecules bound to the apo form of the host.

The coupling models for the guests to TEMOA are more complex, reflecting
higher conformational heterogeneity. This is particularly evident for the
bound state complexes of TEMOA with G1 and G3 whose perturbation energy
probability densities at $\lambda=1$ have two modes (Figures 6 and 7) each
corresponding to a conformational state with a significant population that
contributes to the binding. In general, we found that the coupling models of
the guests to TEMOA are described by at least three modes: a binding-competent
mode characterized by a small statistical weight (the $w$ parameter in Table
3) and a relatively small probability of clashes (large $b$), another mode
extremely unfavorable to binding dominated by clashes, and a moderate mode
that is not as unfavorable (Table 3). Consistently with their more favorable
coupling free energies (Table 2) and higher overall binding affinities, the
binding-competent modes of the TEMOA-G2 and TEMOA-G4 complexes (mode 1) tend
to have larger $b$ values (smaller chance of atomic clashes) and smaller
$\bar{u}_{b}$ values (stronger interactions with the host) than the other
complexes.

The modes unfavorable to binding correspond to positions and orientations of
the guest that cause frequent and severe atomic clashes with TEMOA, as
indicated by the large $n_{l}$ parameter, which is related to the number of
atoms of the ligand that experience very repulsive energies, and a small $b$
value, which reflects the probability of finding a configuration in which the
ligand binds to TEMOA without atomic collisions. As expected based on their
relative sizes, water has a greater probability of binding without collisions
than the larger guests as reflected by the smaller $b$ values of the latter.
We found that various combinations of collisional parameters ($\epsilon$,
$\tilde{u}$, and $n_{l}$) fit the simulation data equally well; the values
reported should be considered order-of-magnitude estimates. Nevertheless, the
optimized $\epsilon$, $\tilde{u}$, and $n_{l}$ values tend to be larger in
magnitude than those of the coupling model of TEMOA-H2O (first column of Table
3 and Fig. 6 top), reflecting the larger role of atomic collisions for the
larger molecular guests, and generally increase in step with how favorable a
mode is towards binding.

Table 3: Parameters for coupling models of the various molecules to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$555kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
TEMOA–H2O
mode 1 | $4.46\times 10^{-1}$ | $1.40\times 10^{-2}$ | $-0.51$ | $2.66$ | $1.0$ | $1.0$ | $3.0$
mode 2 | $5.54\times 10^{-1}$ | $6.90\times 10^{-4}$ | $0.00$ | $3.12$ | $1.0$ | $2.5$ | $5.4$
TEMOA–G1
mode 1 | $1.23\times 10^{-5}$ | $2.71\times 10^{-9}$ | $-19.14$ | $3.72$ | $1.0$ | $1.0$ | $15.9$
mode 2 | $7.10\times 10^{-3}$ | $2.12\times 10^{-18}$ | $-22.23$ | $4.69$ | $1.0$ | $1.0$ | $49.3$
mode 3 | $9.93\times 10^{-1}$ | $9.57\times 10^{-13}$ | $-8.91$ | $4.78$ | $20.0$ | $200.0$ | $60.0$
TEMOA–G2
mode 1 | $2.27\times 10^{-2}$ | $1.43\times 10^{-8}$ | $-23.85$ | $2.58$ | $2.1$ | $2.1$ | $7.4$
mode 2 | $1.99\times 10^{-1}$ | $1.49\times 10^{-6}$ | $-15.95$ | $3.17$ | $5.2$ | $22.4$ | $17.3$
mode 3 | $7.79\times 10^{-1}$ | $1.35\times 10^{-6}$ | $-9.48$ | $3.83$ | $9.0$ | $89.8$ | $46.3$
TEMOA–G3
mode 1 | $1.37\times 10^{-4}$ | $8.56\times 10^{-7}$ | $-15.38$ | $3.72$ | $1.0$ | $1.0$ | $11.7$
mode 2 | $3.48\times 10^{-2}$ | $1.52\times 10^{-8}$ | $-6.58$ | $4.51$ | $1.0$ | $1.0$ | $48.7$
mode 3 | $9.65\times 10^{-1}$ | $2.62\times 10^{-14}$ | $-11.87$ | $4.80$ | $18.0$ | $179.8$ | $60.0$
TEMOA–G4
mode 1 | $3.90\times 10^{-5}$ | $1.55\times 10^{-9}$ | $-27.20$ | $3.00$ | $1.0$ | $1.0$ | $14.3$
mode 2 | $3.07\times 10^{-2}$ | $3.29\times 10^{-10}$ | $-20.79$ | $3.54$ | $1.0$ | $1.0$ | $50.5$
mode 3 | $9.69\times 10^{-1}$ | $8.32\times 10^{-11}$ | $-11.57$ | $4.49$ | $18.8$ | $188.0$ | $60.0$
TEMOA–G5
mode 1 | $2.82\times 10^{-5}$ | $7.28\times 10^{-9}$ | $-24.56$ | $2.85$ | $1.0$ | $1.0$ | $11.4$
mode 2 | $3.00\times 10^{-2}$ | $4.75\times 10^{-9}$ | $-15.31$ | $3.80$ | $1.0$ | $1.0$ | $44.1$
mode 3 | $9.70\times 10^{-1}$ | $5.04\times 10^{-28}$ | $1.99$ | $5.23$ | $12.9$ | $123.8$ | $60.0$

(a) TEMOA-H2O coupling (b) H2O hydration

(c) TEMOA-G1 coupling (d) G1 hydration

(e) TEMOA-G2 coupling (f) G2 hydration

Figure 6: Probability densities $p_{\lambda{\rm h}}(u_{\rm h,sc})$ and
$p_{\lambda{\rm s}}(u_{\rm s,sc})$ of the soft-core interaction energies for
the coupling of a guest to TEMOA (left) and the coupling of a guest to water
(right) collected from simulations (dots) and predicted from the analytical
model’s descriptions of the probability densities at the decoupled state
$p_{0{\rm h}}(u_{\rm h})$ and $p_{0{\rm s}}(u_{\rm s})$ (lines). Each color
corresponds to an alchemical state: pink $\lambda=0$, purple $\lambda=0.1$,
blue $\lambda=0.25$, green $\lambda=0.4$, yellow $\lambda=0.5$, orange
$\lambda=0.7$, and red $\lambda=1$. The x-axis has units of kcal/mol.

(a) TEMOA-G3 coupling (b) G3 hydration

(c) TEMOA-G4 coupling (d) G4 hydration

(e) TEMOA-G5 coupling (f) G5 hydration

Figure 7: Probability densities $p_{\lambda{\rm h}}(u_{\rm h,sc})$ and
$p_{\lambda{\rm s}}(u_{\rm s,sc})$ of the soft-core interaction energies for
the coupling of a guest to TEMOA (left) and the coupling of a guest to water
(right) collected from simulations (dots) and predicted from the analytical
model’s descriptions of the probability densities at the decoupled state
$p_{0{\rm h}}(u_{\rm h})$ and $p_{0{\rm s}}(u_{\rm s})$ (lines). Each color
corresponds to an alchemical state: pink $\lambda=0$, purple $\lambda=0.1$,
blue $\lambda=0.25$, green $\lambda=0.4$, yellow $\lambda=0.5$, orange
$\lambda=0.7$, and red $\lambda=1$. The x-axis has units of kcal/mol.

###  V.4 Hydration Models

Table 4 lists the optimized parameters for the hydration models of the guests.
The model for coupling a water molecule to the water solvent is represented by
one mode (Table 4). With the exception of G3, which is the most flexible guest
(and G5, although its modes 2 and 3 are nearly equivalent), the $p_{0}(u)$
models for the guests are described by two modes. The interplay between the
modes is particularly evident for the hydration of G1 and G3 that present
bimodal distributions in the $\lambda=1$ hydrated state (Figures 6 and 7).
Similarly to the models for the coupling to the solvated host, the first mode
is more favorable toward hydration with stronger interactions with the solvent
(the $\bar{u}_{b}$ parameter) and a smaller chance of atomic clashes (the $b$
parameter).

The $b$ parameter, which measures the probability of finding configurations
free of clashes in the uncoupled ensemble,[63] has special significance in the
theory of solvation and hydrophobicity,[41, 64, 65, 66] since it is related to
the free energy cost of forming a cavity in the solvent of the size and shape
of the solute, $\Delta G_{\rm cav}=-k_{B}T\log b$. Specifically, the smaller
the $b$ parameter, the higher the free energy cost of cavity formation.
Indeed, the optimized $b$ parameters we obtained of the SAMPL8 guests track
the size of the guests (Table 4), with the smaller guests (G2 and G5) having
the largest $b$ for mode 1 while G1, the largest guest, has the smallest $b$
and the highest free energy of cavity formation. The model predicts that the
probability of finding a cavity the size of a water molecule is approximately
$5.77\times 10^{-3}$ (Table 4) corresponding to a cavity formation free energy
cost of about $3$ kcal/mol. Estimating free energies of cavity formation
through this approach is significant because, similar to early information-
theory models,[67] our analytical model does not assume a specific form of the
solute-solvent repulsive interaction potential.[41, 66] Rather, it exploits
the distinct statistical signatures of atomic collisions, which, arguably, are
the fundamental defining characteristic of the solute cavity.

Table 4: Parameters for hydration models of the SAMPL8 host-guest complexes. | $w_{i}$ | $b$ |  $\bar{u}_{b}$666kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
H2O
mode 1 | $1.00$ | $5.77\times 10^{-3}$ | $2.41$ | $3.46$ | $3.9$ | $3.9$ | $2.5$
G1
mode 1 | $2.79\times 10^{-5}$ | $2.20\times 10^{-11}$ | $-7.51$ | $4.27$ | $1.0$ | $1.0$ | $26.2$
mode 2 | $9.99\times 10^{-1}$ | $3.20\times 10^{-14}$ | $13.42$ | $6.27$ | $2.3$ | $2.3$ | $60.0$
G2
mode 1 | $1.27\times 10^{-3}$ | $4.60\times 10^{-7}$ | $-6.87$ | $3.37$ | $1.0$ | $1.0$ | $13.4$
mode 2 | $9.99\times 10^{-1}$ | $2.41\times 10^{-9}$ | $3.90$ | $4.45$ | $1.0$ | $1.0$ | $36.2$
G3
mode 1 | $2.98\times 10^{-6}$ | $1.17\times 10^{-9}$ | $-11.28$ | $3.12$ | $1.0$ | $1.0$ | $15.1$
mode 2 | $5.57\times 10^{-3}$ | $2.98\times 10^{-16}$ | $-10.69$ | $4.44$ | $1.0$ | $1.0$ | $33.5$
mode 3 | $9.94\times 10^{-1}$ | $1.57\times 10^{-13}$ | $-4.38$ | $4.05$ | $1.5$ | $1.5$ | $60.0$
G4
mode 1 | $8.64\times 10^{-5}$ | $1.23\times 10^{-9}$ | $-6.83$ | $4.11$ | $1.0$ | $1.0$ | $20.2$
mode 2 | $9.99\times 10^{-1}$ | $1.70\times 10^{-12}$ | $9.74$ | $5.68$ | $1.0$ | $1.0$ | $57.9$
G5
mode 1 | $7.94\times 10^{-8}$ | $3.68\times 10^{-6}$ | $-5.86$ | $3.99$ | $1.0$ | $1.0$ | $9.2$
mode 2 | $1.60\times 10^{-3}$ | $2.30\times 10^{-16}$ | $-11.04$ | $4.49$ | $1.0$ | $1.0$ | $26.3$
mode 3 | $9.98\times 10^{-1}$ | $0.00$ | $-10.40$ | $4.47$ | $1.0$ | $1.0$ | $56.9$

###  V.5 Analytical Models of Alchemical Transfer

As discussed in the Theory section and illustrated in Figure 1, the models for
the perturbation energy distributions of the alchemical transfer of a guest
from solution to the host and back are constructed from the convolution of the
coupling and hydration models described above [Eqs. (19)–(21)]. Alchemical
transfer is described by two thermodynamic legs, the first starting from the
host and the guest dissociated in solution and the other from the guest bound
to the host. Both legs terminate at the same alchemical intermediate state at
$\lambda=1/2$. We present the models for each leg individually (Figures 8 and
9). Unlike the coupling transformations, the initial states at $\lambda=0$ do
not represent uncoupled states but rather states where host and guest are
coupled to the bulk solvent or to each other.

Each pair of coupling modes, one for coupling to the host and one for coupling
to the solvent, combine to produce a transfer mode for the first leg and one
for the second leg. Hence, for example, the $p_{0}^{+}(u_{t})$ model for leg 1
of the TEMOA-G1 complex is composed of six modes (Table 6) one for each
combination of the three coupling modes of G1 to the host (Table 3) and the
two coupling modes to the solvent (4). One of the more intricate systems in
the set is the alchemical transfer for TEMOA-G3 which is modeled by nine modes
from all possible combinations of the three modes for hydration and host’s
coupling. Depending on their resulting statistical weights, some of these
modes are not apparent in the distributions in Figures 8 and 9, and others
contribute significantly only at some $\lambda$-values and in a specific range
of perturbation energies.

The parameters for each transfer mode are calculated using straightforward
relations described below Eqs. (19)–(21): the statistical weight is the
product of the statistical weights of the coupling modes, the $\bar{u}_{b}$
parameter is given by the difference between the corresponding parameter of
the host coupling model and the average solute-solvent interaction energy in
the fully solvated state, the $\sigma_{b}$ parameter is the geometric average
of the corresponding coupling parameters, and the $b$, $\epsilon$,
$\tilde{u}$, and $n_{l}$ parameters are inherited directly from the coupling
mode to the host. The parameters for the $p_{0}^{-}(u_{t})$ model for leg 2
are derived similarly, except that the initial and final states are reversed.
For example, the $\bar{u}_{b}$ parameter for a mode of leg 2 is given by the
difference between the corresponding parameter of the solvent coupling model
and the average host-guest interaction energy in the bound state.

The parameters of the alchemical transfer models are presented in Tables 6 to
10. Figure 9 shows the corresponding perturbation energy probability densities
(continuous curves). Despite their complexities, the analytical models’
predictions agree well with the perturbation energy distributions obtained and
numerical Alchemical Transfer (ATM) simulations (dots, in Figures 8 and 9).
This result confirms the validity of the theory and the assumptions on which
it is based, and it provides a physical interpretation of the distributions
observed from alchemical transfer calculations.

Similarly to the host coupling and hydration processes (Figures 6 and 7), near
$\lambda=0$, when the guest interacts only with the solvent, the distributions
for transfer towards binding (leg 1 in Figures 8 and 9) lie at very
unfavorable perturbation energies because they are dominated by clashes
between the guest and the atoms of the host (and bound water molecules) when
the guest is transferred from the solvent bulk. Analogously, clashes between
the guest and the solvent’s atoms dominate the probability distributions of
the perturbation energy for transfer in leg 2. The shapes of the distributions
in this regime are determined by the collisional parameters of the
corresponding coupling processes. The distributions progressively shift to
lower perturbation energies As $\lambda$ increases and the interactions at the
transfer position are turned on. The complex multimodal behavior near
$u_{sc}=100$ kcal/mol is an artifact of the soft-core function that terminates
there (see Computational Details). As expected, near $\lambda=1/2$, the
distribution display linear response behavior. The $\lambda=1/2$ states of leg
1 and leg 2 are equivalent and the corresponding distributions differ only in
the sign of the perturbation energy [Eqs. 22 and 25]. Hence, as evidenced in
Figures 8 and 9, they are related by mirror symmetry about zero.

Table 5: Parameters for transfer model of one water molecule binding to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$777kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
Leg 1
mode 1 | $4.62\times 10^{-1}$ | $5.82\times 10^{-3}$ | $14.8$ | $4.36$ | $3.65$ | $4.34$ | $2.50$
mode 2 | $5.38\times 10^{-1}$ | $5.82\times 10^{-3}$ | $18.8$ | $4.66$ | $3.65$ | $4.34$ | $2.50$
Leg 2
mode 1 | $4.18\times 10^{-1}$ | $1.49\times 10^{-2}$ | $17.1$ | $4.36$ | $1.00$ | $1.00$ | $2.95$
mode 2 | $5.82\times 10^{-1}$ | $6.50\times 10^{-4}$ | $17.7$ | $4.66$ | $1.00$ | $0.00$ | $6.51$

(a) TEMOA-H2O transfer leg 1 (b) TEMOA-H2O transfer leg 2

(c) TEMOA-G1 transfer leg 1 (d) TEMOA-G1 transfer leg 2

(e) TEMOA-G2 transfer leg 1 (f) TEMOA-G2 transfer leg 2

Figure 8:  Probability densities $p_{\lambda{\rm t}}(u_{\rm t,sc})$ of the
soft-core perturbation energy of alchemical transfer collected from
simulations (dots) and predicted from the analytical model’s description of
$p_{0{\rm t}}(u_{\rm t})$ at the initial state (lines). Each color represents
a distinct alchemical state: pink $\lambda=0.5$, purple $\lambda=0.4$, blue
$\lambda=0.3$, green $\lambda=0.2$, yellow $\lambda=0.1$, and orange
$\lambda=0$. Left: leg 2, which describes unbinding. Right: leg 1, which
describes binding. The x-axis is in units of kcal/mol.

(a) TEMOA-G3 transfer leg 1 (b) TEMOA-G3 transfer leg 2

(c) TEMOA-G4 transfer leg 1 (d) TEMOA-G4 transfer leg 2

(e) TEMOA-G5 transfer leg 1 (f) TEMOA-G5 transfer leg 2

Figure 9:  Probability densities $p_{\lambda{\rm t}}(u_{\rm t,sc})$ of the
soft-core perturbation energy of alchemical transfer collected from
simulations (dots) and predicted from the analytical model’s description of
$p_{0{\rm t}}(u_{\rm t})$ at the initial state (lines). Each color represents
a distinct alchemical state: pink $\lambda=0.5$, purple $\lambda=0.4$, blue
$\lambda=0.3$, green $\lambda=0.2$, yellow $\lambda=0.1$, and orange
$\lambda=0$. Left: leg 2, which describes unbinding. Right: leg 1, which
describes binding. The x-axis is in units of kcal/mol.

##  VI Discussion

The Potential Distribution Theorem (PDT)[16] offers a useful formalism to
describe alchemical transformations.[21, 20] It leads naturally to a
representation of alchemical coupling in terms of the probability density
distribution $p_{0}(u)$ of the interaction energy $u$ between a ligand and its
environment (a solvent or a receptor) collected in the uncoupled ensemble.[19,
18] The PDT formulas [Eqs. (8) and (9)] relate $p_{0}(u)$ to the free energy
profile and the sequence of perturbation energy distributions along the
alchemical pathway.

Because $p_{0}(u)$ does not depend on the alchemical $\lambda$-dependent
potential energy function, it describes any alchemical process that connects
the uncoupled and coupled states of the system. The central role of $p_{0}(u)$
in the statistical mechanics of alchemical transformations is analogous to
that of the density of states $\Omega(E)$ in standard statistical
mechanics.[68] Note, for instance, the parallel between the PDT relationship
[Eq. (8)], which gives the perturbation energy distributions as a function of
the alchemical progress parameter $\lambda$, and the well known canonical
ensemble relationship $p_{\beta}(E)\propto\Omega(E)\exp(-\beta E)$ relating
the energy distribution of a system’s energy as a function of temperature to
the density of states.[20] In both cases, knowledge of a single function (the
density of states or $p_{0}(u)$) uniquely determined the statistical behavior
of the system for all values of a system’s parameter (the temperature in the
case of the canonical ensemble or $\lambda$ in the case of an alchemical
process). Hence, $p_{0}(u)$ can be considered a master function to describe
the thermodynamics of alchemical states in the same way that the density of
states describes the thermodynamics of physical systems.

Taking advantage of that the perturbation energy for transfer is the sum of
the coupling energy of the ligand to the solvated receptor and the uncoupling
energy from the solvent, in this work, we developed an analytical PDT
description of alchemical transfer for binding by expressing the $p_{0}(u)$
function for transfer as the convolution of the corresponding functions for
uncoupling the ligand to the solvent and to the solvated receptor. The
$p_{0}(u)$ function of each coupling process is expressed in terms of the
parameters of the analytical model of alchemical coupling of Kilburg and
Gallicchio[21, 20] obtained by maximum likelihood analysis of double-
decoupling alchemical simulations. We showed that the resulting analytical
model for transfer reproduces the perturbation energy distributions observed
in alchemical transfer simulations of host-guest complexes in explicit
solvent.

Alchemical Transfer (ATM)[23] is a method developed recently to compute the
absolute and relative binding free energies of molecular complexes, including
those of protein-ligand complexes relevant to computer-aided drug
discovery.[69, 26, 27, 28] Because it is based on a direct coordinate
transformation and a simple dual-topology implementation,[70, 39] ATM is
easily applied to complex scaffold-hopping and charge-changing
transformations.[27] ATM is also applicable with any force field, including
many-body and machine-learned potentials that started to be deployed in drug
discovery projects.[29, 30] This work builds a solid theoretical foundation
for alchemical transfer and provides physical insights on the origin of the
complex perturbation energy distributions often observed in alchemical
transfer simulations. Furthermore the results of this work verify the
theoretical and numerical consistency betweeen alchemical transfer and the
more established double-decoupling alchemical descriptions of binding.[11, 10]

More generally, the work reinforces the benefits of a view of alchemical
processes in terms of the progressive modifications of the statistical
distributions of the system. For example, we adopted this approach to develop
a graphical scheme to optimize the form and the parameters of alchemical
potential energy functions to enhance convergence by avoiding alchemically
induced phase transitions.[20] In this work, we use a similar approach based
on the PDT to illustrate how all of the alchemical pathways originating from
one state and ending in another are interdependent. Specifically, we showed
that distributions, particularly those of intermediate $\lambda$ states that
are unphysical, learned from double-decoupling simulations of binding provide
information to reproduce those observed in alchemical transfer calculations.

Unfortunately, the insights obtained here do not immediately transfer to some
of the popular alchemical models in current use. The PDT formalism assumes an
alchemical potential energy function that depends on one or, at most, a few
collective variables with interpretable statistical distributions. Hence, the
PDT does not apply to the double-decoupling alchemical models based on
parameter interpolation and soft-core pair potentials implemented in some MD
engines,[71, 3, 9, 72] whose perturbation energy depends in complex ways on
atomic coordinates directly. However, the results obtained here should be
applicable to most dual-topology alchemical models based on energy
interpolation.[73, 74]

##  VII Conclusions

We presented a Potential Distribution Theory (PDT)[16, 18] description of the
Alchemical Transfer Method (ATM)[22, 23, 24] for molecular binding. The
probability density $p_{0}(u)$ of the perturbation energy for transfer at the
dissociated state of the solvated complex, which is the central quantity for
PDT, is expressed as the convolution of the probability densities for
decoupling the ligand from the solvent and coupling it to the solvated
receptor obtained from double-decoupling alchemical calculations. We tested
the theory on the alchemical binding of five guests to the TEMOA host from the
SAMPL8 benchmark set. In each case, the probability densities of the
perturbation energy for transfer along the alchemical transfer pathway
obtained from numerical calculations match those predicted from the double-
decoupling distributions represented using the analytical model of alchemical
coupling of Kilburg and Gallicchio.[21, 20] The results of the work provide a
solid theoretical foundation for alchemical transfer, provide physical
insights on the form of the probability densities observed in alchemical
transfer calculations, and confirm the conceptual and numerical equivalence
between the alchemical transfer and double-decoupling processes.

##  VIII Appendix

###  VIII.1 Other Representations of $p_{0}(u)$

The predicted $\log p_{0}(u_{\rm sc})$ functions and the corresponding
$\lambda$-functions $\lambda_{0}(u_{\rm sc})$ with respect to the soft-core
perturbation energies of the molecular complexes studied in this work. The
analytical models for the alchemical coupling to the solvated host are in
pink, and those for alchemical hydration are in blue. The yellow lines
represent Gaussian kernel estimates from samples collected from numerical
alchemical simulations.

(a)

(b)

Figure 10: TEMOA-G1.

(a)

(b)

Figure 11: TEMOA-G2.

(a)

(b)

Figure 12: TEMOA-G3.

(a)

(b)

Figure 13: TEMOA-G4.

(a)

(b)

Figure 14: TEMOA-G5

###  VIII.2 Optimized Parameters of the Analytical Model of Alchemical
Transfer

Table 6: Parameters for the transfer model of guest G1 binding to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$888kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
Leg 1
mode 1 | $8.91\times 10^{-7}$ | $6.40\times 10^{-11}$ | $34.4$ | $5.59$ | $13.8$ | $386$ | $4.89$
mode 2 | $1.07\times 10^{-1}$ | $1.39\times 10^{-14}$ | $52.6$ | $7.05$ | $13.6$ | $0.00$ | $41.0$
mode 3 | $1.23\times 10^{-7}$ | $6.40\times 10^{-11}$ | $36.1$ | $6.05$ | $13.8$ | $386$ | $4.89$
mode 4 | $5.38\times 10^{-1}$ | $1.39\times 10^{-14}$ | $54.3$ | $7.42$ | $13.6$ | $0.00$ | $41.0$
mode 5 | $5.38\times 10^{-1}$ | $6.40\times 10^{-11}$ | $51.1$ | $6.41$ | $13.8$ | $386$ | $4.89$
mode 6 | $5.38\times 10^{-1}$ | $1.39\times 10^{-14}$ | $69.2$ | $7.72$ | $13.6$ | $0.00$ | $41.0$
Leg 2
mode 1 | $6.69\times 10^{-6}$ | $1.63\times 10^{-9}$ | $18.1$ | $5.59$ | $13.5$ | $237$ | $3.83$
mode 2 | $7.71\times 10^{-1}$ | $2.26\times 10^{-13}$ | $25.5$ | $6.05$ | $24.4$ | $1220$ | $50.0$
mode 3 | $9.73\times 10^{-3}$ | $7.18\times 10^{-18}$ | $18.0$ | $6.41$ | $18.6$ | $0.00$ | $30.0$
mode 4 | $1.88\times 10^{-6}$ | $1.63\times 10^{-9}$ | $30.9$ | $7.05$ | $13.5$ | $237$ | $3.83$
mode 5 | $2.16\times 10^{-1}$ | $2.26\times 10^{-13}$ | $38.2$ | $7.42$ | $24.4$ | $1220$ | $50.0$
mode 6 | $2.73\times 10^{-3}$ | $7.18\times 10^{-18}$ | $30.8$ | $7.72$ | $18.6$ | $0.00$ | $30.0$
Table 7: Parameters for the transfer model of guest G2 binding to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$999kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
Leg 1
mode 1 | $1.03\times 10^{-2}$ | $2.60\times 10^{-6}$ | $7.23$ | $4.58$ | $3.50$ | $3.52$ | $5.29$
mode 2 | $9.02\times 10^{-1}$ | $4.08\times 10^{-8}$ | $18.9$ | $3.91$ | $21.9$ | $305$ | $31.0$
mode 3 | $8.05\times 10^{-2}$ | $1.44\times 10^{-5}$ | $13.3$ | $4.93$ | $8.75$ | $63.3$ | $7.32$
mode 4 | $7.33\times 10^{-5}$ | $2.60\times 10^{-6}$ | $9.86$ | $5.33$ | $3.50$ | $3.52$ | $5.29$
mode 5 | $6.42\times 10^{-3}$ | $4.08\times 10^{-8}$ | $21.6$ | $4.77$ | $21.9$ | $305$ | $31.0$
mode 6 | $5.73\times 10^{-4}$ | $1.44\times 10^{-5}$ | $15.9$ | $5.64$ | $8.75$ | $62.3$ | $7.32$
Leg 2
mode 1 | $6.57\times 10^{-4}$ | $7.54\times 10^{-7}$ | $28.0$ | $4.58$ | $13.9$ | $41.9$ | $4.17$
mode 2 | $8.50\times 10^{-1}$ | $1.85\times 10^{-9}$ | $37.9$ | $5.33$ | $7.77$ | $0.00$ | $22.8$
mode 3 | $8.76\times 10^{-16}$ | $7.54\times 10^{-7}$ | $6.76$ | $3.91$ | $13.9$ | $41.9$ | $4.17$
mode 4 | $1.13\times 10^{-12}$ | $1.85\times 10^{-9}$ | $16.7$ | $4.77$ | $7.77$ | $0.00$ | $22.8$
mode 5 | $1.15\times 10^{-4}$ | $7.54\times 10^{-7}$ | $27.5$ | $4.93$ | $13.9$ | $41.9$ | $4.17$
mode 6 | $1.49\times 10^{-1}$ | $1.85\times 10^{-9}$ | $37.4$ | $5.64$ | $7.77$ | $0.00$ | $22.8$
Table 8: Parameters for the transfer model of guest G3 binding to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$101010kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
Leg 1
mode 1 | $3.04\times 10^{-6}$ | $8.68\times 10^{-10}$ | $27.3$ | $4.89$ | $7.62$ | $58.9$ | $5.32$
mode 2 | $7.68\times 10^{-1}$ | $1.38\times 10^{-13}$ | $51.5$ | $5.39$ | $4.36$ | $39.7$ | $34.4$
mode 3 | $6.42\times 10^{-3}$ | $3.69\times 10^{-14}$ | $37.1$ | $6.10$ | $17.2$ | $47.5$ | $13.5$
mode 4 | $8.34\times 10^{-7}$ | $8.68\times 10^{-10}$ | $53.6$ | $6.09$ | $7.62$ | $58.9$ | $5.32$
mode 5 | $2.11\times 10^{-1}$ | $1.38\times 10^{-13}$ | $35.5$ | $6.50$ | $4.36$ | $39.8$ | $34.4$
mode 6 | $1.76\times 10^{-3}$ | $3.69\times 10^{-14}$ | $52.0$ | $7.10$ | $17.2$ | $47.5$ | $13.5$
mode 7 | $4.95\times 10^{-8}$ | $8.68\times 10^{-10}$ | $52.0$ | $5.50$ | $7.62$ | $58.9$ | $5.32$
mode 8 | $1.25\times 10^{-2}$ | $1.38\times 10^{-13}$ | $52.0$ | $5.95$ | $4.36$ | $39.7$ | $34.4$
mode 9 | $1.05\times 10^{-4}$ | $3.69\times 10^{-14}$ | $52.0$ | $6.61$ | $17.2$ | $47.5$ | $13.5$
Leg 2
mode 1 | $8.40\times 10^{-7}$ | $9.38\times 10^{-7}$ | $12.8$ | $4.89$ | $3.91$ | $0.00$ | $8.55$
mode 2 | $6.51\times 10^{-3}$ | $7.91\times 10^{-12}$ | $23.0$ | $6.09$ | $26.8$ | $396$ | $49.6$
mode 3 | $1.82\times 10^{-4}$ | $1.91\times 10^{-8}$ | $21.5$ | $5.50$ | $7.00$ | $70.0$ | $19.7$
mode 4 | $1.01\times 10^{-6}$ | $9.38\times 10^{-7}$ | $15.0$ | $5.39$ | $3.91$ | $0.00$ | $8.55$
mode 5 | $7.80\times 10^{-3}$ | $7.91\times 10^{-12}$ | $25.3$ | $6.50$ | $26.8$ | $396$ | $49.6$
mode 6 | $2.18\times 10^{-4}$ | $1.91\times 10^{-8}$ | $23.8$ | $5.95$ | $7.00$ | $70.0$ | $19.7$
mode 7 | $1.24\times 10^{-4}$ | $9.38\times 10^{-7}$ | $28.5$ | $6.10$ | $3.91$ | $0.00$ | $8.55$
mode 8 | $9.58\times 10^{-1}$ | $7.91\times 10^{-12}$ | $38.7$ | $7.10$ | $26.8$ | $396$ | $49.6$
mode 9 | $2.68\times 10^{-2}$ | $1.91\times 10^{-8}$ | $37.2$ | $6.61$ | $7.00$ | $70.0$ | $19.7$

Table 9: Parameters for the transfer model of guest G4 binding to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$111111kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
Leg 1
mode 1 | $4.53\times 10^{-7}$ | $2.99\times 10^{-9}$ | $35.0$ | $5.63$ | $9.94$ | $49.8$ | $6.69$
mode 2 | $1.05\times 10^{-2}$ | $1.99\times 10^{-12}$ | $51.5$ | $6.85$ | $12.3$ | $17.9$ | $27.4$
mode 3 | $9.61\times 10^{-7}$ | $2.99\times 10^{-9}$ | $37.1$ | $5.93$ | $9.94$ | $49.7$ | $6.69$
mode 4 | $2.24\times 10^{-2}$ | $1.99\times 10^{-12}$ | $53.6$ | $7.10$ | $12.3$ | $17.9$ | $27.4$
mode 5 | $4.15\times 10^{-5}$ | $2.99\times 10^{-9}$ | $35.5$ | $5.35$ | $9.94$ | $49.7$ | $6.69$
mode 6 | $9.67\times 10^{-1}$ | $1.99\times 10^{-12}$ | $52.0$ | $6.63$ | $12.3$ | $17.9$ | $27.4$
Leg 2
mode 1 | $5.19\times 10^{-3}$ | $1.65\times 10^{-9}$ | $18.2$ | $5.63$ | $10.8$ | $70.2$ | $13.1$
mode 2 | $9.62\times 10^{-1}$ | $7.08\times 10^{-11}$ | $21.8$ | $5.93$ | $28.6$ | $587$ | $38.8$
mode 3 | $7.79\times 10^{-6}$ | $5.62\times 10^{-8}$ | $12.5$ | $5.35$ | $7.82$ | $56.3$ | $3.92$
mode 4 | $1.76\times 10^{-4}$ | $1.65\times 10^{-9}$ | $27.4$ | $6.86$ | $10.8$ | $70.2$ | $13.1$
mode 5 | $3.26\times 10^{-2}$ | $7.08\times 10^{-11}$ | $31.0$ | $7.10$ | $28.6$ | $587$ | $38.8$
mode 6 | $2.64\times 10^{-7}$ | $5.62\times 10^{-7}$ | $21.7$ | $6.63$ | $7.82$ | $56.3$ | $3.92$
Table 10: Parameters for the transfer model of guest G5 binding to host TEMOA. | $w_{i}$ | $b$ |  $\bar{u}_{b}$121212kcal/mol | $\sigma_{b}^{\rm a}$ | $\epsilon^{\rm a}$ | $\tilde{u}^{\rm a}$ | $n_{l}$
---|---|---|---|---|---|---|---
Leg 1
mode 1 | $5.22\times 10^{-8}$ | $5.52\times 10^{-6}$ | $33.1$ | $4.95$ | $1.99$ | $19.9$ | $4.44$
mode 2 | $1.03\times 10^{-3}$ | $2.65\times 10^{-14}$ | $32.9$ | $5.63$ | $3.79$ | $22.9$ | $12.6$
mode 3 | $9.96\times 10^{-1}$ | $8.61\times 10^{-18}$ | $39.9$ | $5.99$ | $24.3$ | $0.00$ | $26.7$
mode 4 | $9.96\times 10^{-11}$ | $5.52\times 10^{-6}$ | $30.7$ | $5.63$ | $1.99$ | $19.9$ | $4.44$
mode 5 | $3.37\times 10^{-6}$ | $2.65\times 10^{-14}$ | $30.4$ | $6.23$ | $3.79$ | $22.9$ | $12.6$
mode 6 | $3.26\times 10^{-3}$ | $8.61\times 10^{-18}$ | $37.4$ | $6.560$ | $24.3$ | $0.00$ | $26.7$
mode 7 | $4.10\times 10^{-34}$ | $5.52\times 10^{-6}$ | $30.8$ | $6.42$ | $1.99$ | $19.9$ | $4.44$
mode 8 | $8.08\times 10^{-30}$ | $2.65\times 10^{-14}$ | $30.6$ | $6.95$ | $3.79$ | $22.9$ | $12.6$
mode 9 | $7.82\times 10^{-27}$ | $8.61\times 10^{-18}$ | $37.6$ | $7.25$ | $24.3$ | $0.00$ | $26.7$
Leg 2
mode 1 | $4.86\times 10^{-6}$ | $2.92\times 10^{-9}$ | $7.99$ | $4.95$ | $10.9$ | $152$ | $3.05$
mode 2 | $5.01\times 10^{-2}$ | $1.40\times 10^{-9}$ | $22.5$ | $5.63$ | $22.0$ | $217$ | $16.9$
mode 3 | $1.73\times 10^{-1}$ | $5.04\times 10^{-38}$ | $38.3$ | $6.42$ | $50.5$ | $646$ | $43.1$
mode 4 | $1.69\times 10^{-5}$ | $2.92\times 10^{-9}$ | $20.2$ | $5.63$ | $10.9$ | $152$ | $3.05$
mode 5 | $1.74\times 10^{-1}$ | $1.40\times 10^{-9}$ | $34.8$ | $6.23$ | $22.0$ | $217$ | $16.9$
mode 6 | $6.02\times 10^{-1}$ | $5.04\times 10^{-38}$ | $50.6$ | $6.95$ | $50.5$ | $646$ | $43.1$
mode 7 | $1.51\times 10^{-8}$ | $2.92\times 10^{-9}$ | $20.3$ | $5.99$ | $10.9$ | $152$ | $3.05$
mode 8 | $1.56\times 10^{-4}$ | $1.40\times 10^{-9}$ | $34.8$ | $6.56$ | $22.0$ | $217$ | $16.9$
mode 9 | $5.39\times 10^{-4}$ | $5.04\times 10^{-38}$ | $50.6$ | $7.25$ | $50.5$ | $646$ | $43.1$

## References

  * 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).
  * Tuckerman [2023] M. E. Tuckerman,  _Statistical mechanics: theory and molecular simulation_ (Oxford university press, 2023).
  * Abel  _et al._ [2017] R. Abel, L. Wang, E. D. Harder, B. Berne, and R. A. Friesner, “Advancing drug discovery through enhanced free energy calculations,” Acc. Chem. Res. 50, 1625–1632 (2017).
  * Schindler  _et al._ [2020] C. E. Schindler, H. Baumann, A. Blum, D. Böse, H.-P. Buchstaller, L. Burgdorf, D. Cappel, E. Chekler, P. Czodrowski, D. Dorsch, E. M. K. I, B. Follows, T. Fuchs, U. Grädler, J. Gunera, T. Johnson, C. J. Lebrun, S. Karra, M. Klein, T. Knehans, L. Koetzner, M. Krier, M. Leiendecker, B. Leuthner, L. Li, I. Mochalkin, D. Musil, C. Neagu, F. Rippmann, K. Schiemann, R. Schulz, T. Steinbrecher, E.-M. Tanzer, A. Unzue Lopez, A. Viacava Follis, A. Wegener, and D. Kuhn, “Large-scale assessment of binding free energy calculations in active drug discovery projects,” J. Chem. Inf. Model. 60, 5457–5474 (2020).
  * Ganguly  _et al._ [2022] A. Ganguly, H.-C. Tsai, M. Fernández-Pendás, T.-S. Lee, T. J. Giese, and D. M. York, “Amber drug discovery boost tools: Automated workflow for production free-energy simulation setup and analysis (professa),” J. Chem. Inf. Model. 62, 6069–6083 (2022).
  * Procacci [2022] P. Procacci, “Relative binding free energy between chemically distant compounds using a bidirectional nonequilibrium approach,” J. Chem. Theory Comput. 18, 4014–4026 (2022).
  * Baumann  _et al._ [2023] H. M. Baumann, E. Dybeck, C. L. McClendon, F. C. Pickard IV, V. Gapsys, L. Pérez-Benito, D. F. Hahn, G. Tresadern, A. M. Mathiowetz, and D. L. Mobley, “Broadening the scope of binding free energy calculations using a separated topologies approach,” J. Chem. Theory Comput. 19, 5058–5076 (2023).
  * Allen  _et al._ [2022] B. K. Allen, M. M. Kulkarni, B. Chamberlain, T. Dwight, C. Koh, R. Samant, F. Jernigan, J. Rice, D. Tan, S. Li, K. Marino, H. Huang, E. Chiswick, B. Tesar, S. Sparks, Z. Lin, T. D. McGee, I. Kolossváry, C. Lin, S. Shechter, H. Soutter, C. Bastos, M. Taimi, S. Lai, A. Petrin, T. Kane, S. Swann, H. Gardner, C. Winter, and W. Sherman, “Design of a systemic small molecule clinical sting agonist using physics-based simulations and artificial intelligence,” bioRxiv (2022), 10.1101/2022.05.23.493001.
  * Lee  _et al._ [2020a] T.-S. Lee, B. K. Allen, T. J. Giese, Z. Guo, P. Li, C. Lin, T. D. McGee Jr, D. A. Pearlman, B. K. Radak, Y. Tao, H.-C. Tsai, H. Xu, W. Sherman, and D. M. York, “Alchemical binding free energy calculations in AMBER20: Advances and best practices for drug discovery,” J. Chem. Inf. Model. 60, 5595–5623 (2020a).
  * Mey  _et al._ [2020] A. S. J. S. Mey, B. K. Allen, H. E. B. Macdonald, J. D. Chodera, D. F. Hahn, M. Kuhn, J. Michel, D. L. Mobley, L. N. Naden, S. Prasad, A. Rizzi, J. Scheen, M. R. Shirts, G. Tresadern, and H. Xu, “Best practices for alchemical free energy calculations [article v1.0],” Living Journal of Computational Molecular Science 2, 18378 (2020).
  * Gilson  _et al._ [1997] M. K. Gilson, J. A. Given, B. L. Bush, and J. A. McCammon, “The statistical-thermodynamic basis for computation of binding affinities: A critical review,” Biophys. J. 72, 1047–1069 (1997).
  * Widom [1982a] B. J. Widom, Chem. Phys. 86, 869 (1982a).
  * Valleau and Torrie [1977] J. Valleau and G. Torrie, “A guide to monte carlo for statistical mechanics: 2. byways,” in  _Statistical Mechanics: Part A: Equilibrium Techniques_ (Springer, 1977) pp. 169–194.
  * Simonson [1993] T. Simonson, “Free energy of particle insertion: an exact analysis of the origin singularity for simple liquids,” Molecular Physics 80, 441–447 (1993).
  * Stamatopoulou and Ben-Amotz [1998] A. Stamatopoulou and D. Ben-Amotz, “Cavity formation free energies for rigid chains in hard sphere fluids,” J. Chem. Phys. 108, 7294–7300 (1998).
  * Beck, Paulaitis, and Pratt [2006] T. L. Beck, M. E. Paulaitis, and L. R. Pratt,  _The Potential Distribution Theorem and Models of Molecular Solutions_ (Cambridge University Press, New York, 2006).
  * Swope and Andersen [1984] W. C. Swope and H. C. Andersen, “A molecular dynamics method for calculating the solubility of gases in liquids and the hydrophobic hydration of inert-gas atoms in aqueous solution,” The Journal of Physical Chemistry 88, 6548–6556 (1984).
  * Gallicchio and Levy [2011] E. Gallicchio and R. M. Levy, “Recent theoretical and computational advances for modeling protein-ligand binding affinities,” Adv. Prot. Chem. Struct. Biol. 85, 27–80 (2011).
  * Gallicchio, Lapelosa, and Levy [2010] E. Gallicchio, M. Lapelosa, and R. M. Levy, “Binding energy distribution analysis method (BEDAM) for estimation of protein-ligand binding affinities,” J. Chem. Theory Comput. 6, 2961–2977 (2010).
  * Pal and Gallicchio [2019] R. K. Pal and E. Gallicchio, “Perturbation potentials to overcome order/disorder transitions in alchemical binding free energy calculations,” J. Chem. Phys. 151, 124116 (2019).
  * Kilburg and Gallicchio [2018] D. Kilburg and E. Gallicchio, “Analytical model of the free energy of alchemical molecular binding,” J. Chem. Theory Comput. 14, 6183–6196 (2018).
  * Khuttan  _et al._ [2021] S. Khuttan, S. Azimi, J. Z. Wu, and E. Gallicchio, “Alchemical transformations for concerted hydration free energy estimation with explicit solvation,” J. Chem. Phys. 154, 054103 (2021).
  * Wu  _et al._ [2021] J. Z. Wu, S. Azimi, S. Khuttan, N. Deng, and E. Gallicchio, “Alchemical transfer approach to absolute binding free energy estimation,” J. Chem. Theory Comput. 17, 3309–3319 (2021).
  * Azimi  _et al._ [2022a] S. Azimi, J. Z. Wu, S. Khuttan, T. Kurtzman, N. Deng, and E. Gallicchio, “Application of the alchemical transfer and potential of mean force methods to the sampl8 host-guest blinded challenge,” J. Comp. Aid. Mol. Des. 36, 63–76 (2022a).
  * Khuttan  _et al._ [2023] S. Khuttan, S. Azimi, J. Z. Wu, S. Dick, C. Wu, H. Xu, and E. Gallicchio, “Taming multiple binding poses in alchemical binding free energy prediction: the $\beta$-cyclodextrin host–guest sampl9 blinded challenge,” Phys. Chem. Chem. Phys. 25, 24364–24376 (2023).
  * Azimi  _et al._ [2022b] S. Azimi, S. Khuttan, J. Z. Wu, R. K. Pal, and E. Gallicchio, “Relative binding free energy calculations for ligands with diverse scaffolds with the alchemical transfer method,” J. Chem. Inf. Model. 62, 309–323 (2022b).
  * Chen  _et al._ [2023] L. Chen, Y. Wu, C. Wu, A. Silveira, W. Sherman, H. Xu, and E. Gallicchio, “Performance and analysis of the alchemical transfer method for binding free energy predictions of diverse ligands,” arXiv , 2308.08671 (2023), arXiv:2308.08671 [physics.chem-ph] .
  * Sabanés Zariquiey  _et al._ [2023] F. Sabanés Zariquiey, A. Pérez, M. Majewski, E. Gallicchio, and G. De Fabritiis, “Validation of the alchemical transfer method for the estimation of relative binding affinities of molecular series,” J. Chem. Inf. Model. 63, 2438–2444 (2023).
  * Eastman  _et al._ [2023] P. Eastman, R. Galvelis, R. P. Peláez, C. R. Abreu, S. E. Farr, E. Gallicchio, A. Gorenko, M. M. Henry, F. Hu, J. Huang, A. Andreas Krämer, J. Michel, J. A. Mitchell, V. S. Pande, J. a. P. Rodrigues, J. Rodriguez-Guerra, A. C. Simonett, J. Swails, I. Zhang, J. D. Chodera, G. De Fabritiis, and T. E. Markland, “OpenMM 8: Molecular dynamics simulation with machine learning potentials,” J. Phys. Chem. B 128, 109–116 (2023).
  * Sabanés Zariquiey  _et al._ [2024] F. Sabanés Zariquiey, R. Galvelis, E. Gallicchio, J. D. Chodera, T. E. Markland, and G. De Fabritiis, “Enhancing protein–ligand binding affinity predictions using neural network potentials,” J. Chem. Inf. Model. (2024).
  * Gilson and Zhou [2007] M. K. Gilson and H.-X. Zhou, “Calculation of protein-ligand binding affinities.” Annu. Rev. Biophys. Biomol. Struct. 36, 21–42 (2007).
  * Gallicchio [2021] E. Gallicchio, “Methods in molecular biology, computational peptide science: Methods and protocols,” (Springer Nature, 2021) Chap. Free Energy-Based Computational Methods for the Study of Protein-Peptide Binding Equilibria, pp. 303–334.
  * McCammon and Straatsma [1992] J. A. McCammon and T. P. Straatsma, “Computational alchemy,” Annu. Rev. Phys. Chem. 43, 407 (1992).
  * Jorgensen [2010] W. L. Jorgensen, “Drug discovery: Pulled from a protein’s embrace,” Nature 466, 42–43 (2010).
  * Cournia, Allen, and Sherman [2017] Z. Cournia, B. Allen, and W. Sherman, “Relative binding free energy calculations in drug discovery: recent advances and practical considerations,” J. Chem. Inf. Model. 57, 2911–2937 (2017).
  * Cournia  _et al._ [2020] Z. Cournia, B. K. Allen, T. Beuming, D. A. Pearlman, B. K. Radak, and W. Sherman, “Rigorous free energy simulations in virtual screening,” J. Chem. Inf. Model. 60, 4153–4169 (2020).
  * Rizzi  _et al._ [2020] A. Rizzi, T. Jensen, D. R. Slochower, M. Aldeghi, V. Gapsys, D. Ntekoumes, S. Bosisio, M. Papadourakis, N. M. Henriksen, B. L. De Groot, Z. Cournia, A. Dickson, J. Michel, M. K. Gilson, M. R. Shirts, D. L. Mobley, and J. D. Chodera, “The SAMPL6 SAMPLing challenge: Assessing the reliability and efficiency of binding free energy calculations,” J. Comp. Aid. Mol. Des. , 1–33 (2020).
  * Deng and Roux [2008] Y. Deng and B. Roux, “Computation of binding free energy with molecular dynamics and grand canonical monte carlo simulations.” J Chem Phys 128, 115103 (2008).
  * König  _et al._ [2020] G. König, N. Glaser, B. Schroeder, A. Kubincová, P. H. Hünenberger, and S. Riniker, “An alternative to conventional $\lambda$-intermediate states in alchemical free energy calculations: $\lambda$-enveloping distribution sampling,” J. Chem. Inf. Model. (2020).
  * Lee  _et al._ [2020b] T.-S. Lee, Z. Lin, B. K. Allen, C. Lin, B. K. Radak, Y. Tao, H.-C. Tsai, W. Sherman, and D. M. York, “Improved alchemical free energy calculations with optimized smoothstep softcore potentials,” J. Chem. Theory Comput. 16, 5512–5525 (2020b).
  * Chandler, Weeks, and Andersen [1983] D. Chandler, J. D. Weeks, and H. C. Andersen, “Van der waals picture of liquids, solids, and phase transformations,” Science 220, 787–794 (1983).
  * Alper and Levy [1990] H. Alper and R. M. Levy, “Field strength dependence of dielectric saturation in liquid water,” J. Phys. Chem. 94, 8401–8403 (1990).
  * Alper and Levy [1993] H. Alper and R. M. Levy, “Dielectric and thermodynamic response of a generalized reaction field model for liquid state simulations,” J. Chem. Phys. 99, 9847–9852 (1993).
  * Levy, Belhadj, and Kitchen [1991] R. M. Levy, M. Belhadj, and D. B. Kitchen, “Gaussian fluctuation formula for electrostatic free-energy changes in solution,” The Journal of Chemical Physics 95, 3627–3633 (1991).
  * Åqvist, Medina, and Samuelsson [1994] J. Åqvist, C. Medina, and J.-E. Samuelsson, “A new method for predicting binding affinity in computer-aided drug design,” Prot. Engineering 7, 385–391 (1994).
  * Åqvist and Hansson [1996] J. Åqvist and T. Hansson, “On the validity of electrostatic linear response in polar solvents,” J. Phys. Chem. 100, 9512–9521 (1996).
  * Carlson and Jorgensen [1995] H. Carlson and W. Jorgensen, “An extended linear response method for determining free energies of hydration,” J. Phys. Chem. 99, 10667–10673 (1995).
  * Levy and Gallicchio [1998] R. M. Levy and E. Gallicchio, “Computer simulations with explicit solvent: Recent progress in the thermodynamic decomposition of free energies and in modeling electrostatic effects,” Annu. Rev. Phys. Chem. 49, 531–67 (1998).
  * Jones-Hertzog and Jorgensen [1997] D. Jones-Hertzog and W. Jorgensen, “Binding affinities for sulfonamide inhibitors with human thrombin using Monte Carlo simulations with a linear response method,” J. Med. Chem. 40, 1539–1549 (1997).
  * Simonson [2002] T. Simonson, “Gaussian fluctuations and linear response in an electron transfer protein,” Proc. Natl. Acad. Sci. 99, 6544–6549 (2002).
  * Su  _et al._ [2007] Y. Su, E. Gallicchio, K. Das, E. Arnold, and R. Levy, “Linear interaction energy (lie) models for ligand binding in implicit solvent: Theory and application to the binding of nnrtis to hiv-1 reverse transcriptase,” J. Chem. Theory Comput. 3, 256–277 (2007).
  * Gumbel [2012] E. J. Gumbel,  _Statistics of Extremes_ (Dover Publications, New York, 2012).
  * Tan  _et al._ [2012a] Z. Tan, E. Gallicchio, M. Lapelosa, and R. M. Levy, “Theory of binless multi-state free energy estimation with applications to protein-ligand binding,” J. Chem. Phys. 136, 144102 (2012a).
  * Shirts and Chodera [2008] M. R. Shirts and J. D. Chodera, “Statistically optimal analysis of samples from multiple equilibrium states,” J. Chem. Phys. 129, 124105 (2008).
  * Suating  _et al._ [2020] P. Suating, T. T. Nguyen, N. E. Ernst, Y. Wang, J. H. Jordan, C. L. Gibb, H. S. Ashbaugh, and B. C. Gibb, “Proximal charge effects on guest binding to a non-polar pocket,” Chem. Sci. 11, 3656–3663 (2020).
  * Hahn  _et al._ [2022] D. Hahn, C. Bayly, M. L. Boby, H. Bruce Macdonald, J. Chodera, V. Gapsys, A. Mey, D. Mobley, L. Perez Benito, C. Schindler, G. Tresadern, and G. Warren, “Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [article v1.0],” Living Journal of Computational Molecular Science 4, 1497 (2022).
  * Gallicchio, Levy, and Parashar [2008] E. Gallicchio, R. M. Levy, and M. Parashar, “Asynchronous replica exchange for molecular simulations,” J. Comp. Chem. 29, 788–794 (2008).
  * Gallicchio  _et al._ [2015] E. Gallicchio, J. Xia, W. F. Flynn, B. Zhang, S. Samlalsingh, A. Mentes, and R. M. Levy, “Asynchronous replica exchange software for grid and heterogeneous computing,” Comput. Phys. Commun. 196, 236–246 (2015).
  * ATo [2022] “AToM-OpenMM,” https://github.com/Gallicchio-Lab/AToM-OpenMM (2022).
  * Tan  _et al._ [2012b] Z. Tan, G. Emilio, L. Mauro, and L. Ronald, “Theory of binless multi-state free energy estimation with applications to protein-ligand binding,” J. Chem. Phys 136 (2012b).
  * Amezcua  _et al._ [2022] M. Amezcua, J. Setiadi, Y. Ge, and D. L. Mobley, “An overview of the SAMPL8 host–guest binding challenge,” J. Comp.-Aided Mol. Des. 36, 707–734 (2022).
  * Suating  _et al._ [2022] P. Suating, N. E. Ernst, B. D. Alagbe, H. A. Skinner, J. T. Mague, H. S. Ashbaugh, and B. C. Gibb, “On the nature of guest complexation in water: Triggered wetting–water-mediated binding,” J. Phys. Chem. B 126, 3150–3160 (2022).
  * Widom [1982b] B. Widom, “Potential-distribution theory and the statistical mechanics of fluids,” J. Phys. Chem. 86, 869–872 (1982b).
  * Berne [1996] B. J. Berne, “Inferring the hydrophobic interaction from the properties of neat water,” Proc. Natl. Acad. Sci. USA 93, 8800–8803 (1996).
  * Huang and Chandler [2000] D. M. Huang and D. Chandler, “Temperature and length sclae dependence of hydrophobic effects and their possible implications for protein folding,” Proc. Natl. Acad. Sci. USA 97, 8324–8327 (2000).
  * Gallicchio, Kubo, and Levy [2000] E. Gallicchio, M. M. Kubo, and R. M. Levy, “Enthalpy-entropy and cavity decomposition of alkane hydration free energies: Numerical results and implications for theories of hydrophobic solvation,” J. Phys. Chem. B 104, 6271–6285 (2000).
  * Hummer  _et al._ [1996] G. Hummer, S. Garde, A. E. García, A. Pohorille, and L. R. Pratt, “An information theory model of hydrophobic interactions,” Proc. Natl. Acad. Sci. USA 93, 8951–8955 (1996).
  * Hill [1986] T. L. Hill,  _An Introduction to Statistical Thermodynamics_ (Dover, New York, 1986).
  * Armacost, Riniker, and Cournia [2020] K. A. Armacost, S. Riniker, and Z. Cournia, “Novel directions in free energy methods and applications,” (2020).
  * Rocklin, Mobley, and Dill [2013] G. J. Rocklin, D. L. Mobley, and K. A. Dill, “Separated topologies—a method for relative binding free energy calculations using orientational restraints,” J. Chem. Phys. 138, 02B614 (2013).
  * Farhi and Singh [2017] A. Farhi and B. Singh, “A novel method for calculating relative free energy of similar molecules in two environments,” Comp. Phys. Comm. 212, 132–145 (2017).
  * Gapsys  _et al._ [2020] V. Gapsys, L. Pérez-Benito, M. Aldeghi, D. Seeliger, H. Van Vlijmen, G. Tresadern, and B. L. De Groot, “Large scale relative protein ligand binding affinities using non-equilibrium alchemy,” Chemical Science 11, 1140–1152 (2020).
  * Raman  _et al._ [2020] E. P. Raman, T. J. Paul, R. L. Hayes, and C. L. Brooks III, “Automated, accurate, and scalable relative protein–ligand binding free-energy calculations using lambda dynamics,” J. Chem. Theory Comput. 16, 7895–7914 (2020).
  * König  _et al._ [2021] G. König, B. Ries, P. H. Hünenberger, and S. Riniker, “Efficient alchemical intermediate states in free energy calculations using $\lambda$-enveloping distribution sampling,” Journal of Chemical Theory and Computation 17, 5805–5815 (2021).