# Quantum Secrecy in Thermal States II

Elizabeth Newton  Anne Ghesquière  Freya L. Wilson  Benjamin T. H. Varcoe
Quantum Experimental Group, School of Physics and Astronomy, University of
Leeds, Leeds LS2 9JT, United Kingdom<EMAIL_ADDRESS>Martin Moseley
Airbus Defense & Space

###### Abstract

In this paper we consider a scheme for cryptographic key distribution based on
a variation of continuous variable quantum key distribution called central
broadcast. In the continuous variable central broadcast scheme, security
arises from discord present in the Hanbury Brown and Twiss effect from a
thermal source. The benefit of this scheme is that it expands the range of
frequencies into the microwave regime. Longer wavelengths, where the thermal
photon number is higher and correlations remain robust over long distances,
may even be preferable to optical wavelengths. Assming that Alice controls the
source but not the distribution of the light (eg satellite broadcasts), then
we demonstrate that the central broadcast scheme is robust to an entangling
cloner attack. We establish the security of the protocol both experimentally
and theoretically.

Quantum key distribution (QKD) is rapidly gaining widespread acceptance
Swiss:2007  as a method of secure key exchange and several high bandwidth
devices have been demonstrated.However, having distributed information across
a network, there remains a limitation of key exchange at the user access
point. For the end user, wireless access is the ideal use model. The user
access system must be both inexpensive and accessible without compromising
security and maintaining the ability to work on scales of the order of metres
to tens of metres.

Recently, the potential of thermal states for QKD has been established
sakuya:2017_2 ; Qi:2018 . Although thermal states have sometimes been
described as too noisy weedbrookprl:2010 ; weedbrookpra:2012 , they exhibit
Hanbury Brown and Twiss correlations which have been found to exhibit positive
discord ragy:2013 , a necessary condition for QKDpirandola:2014 .

Consider a central broadcast protocol in which the radiation is split between
twol parties, who now have correlated signals from which they can build a key.
Another advantage to using thermal states is that they are easy and low-cost
to produce. Whereas large-scale implementations of QKD such as those described
above require specific infrastructure, thermal states central broadcasting
protocols can be implemented over short distances, with low-power devices.

In the scheme proposed in sakuya:2017_2 : a thermal source is incident on a
beamsplitter, with one output port connected to Alice and the other to Bob. We
assumed that Alice controls the source, the channel leading to the
beamsplitter and the beamsplitter itself. They also control the channel
separating them from the beamsplitter. The only part opened to Eve resides on
the branch between the beamsplitter and Bob. We found that there is both a
positive key rate and positive discord between the legal parties, both at
optical frequencies (experimental result) and microwave frequencies
(theoretical analysis). In Qi:2018 , the authors shine a thermal source on a
beamsplitter to prepare the states used by the legal parties in QKD, and find
that such a source average photon number of 100 allows for efficient passive
QKD.

In this paper, we relax our security by surrendering control of the channel
between the source and the beamsplitter, leaving it open to attack. Alice
retains control of the source, the beamsplitter, their channel and detector.
Eve can then attack the channel going to Bob (sakuya:2017_2 ) or the channel
going to the beamsplitter. This aims to provide Eve greater knowledge of the
states making up the thermal radiation.

In the following, we describe the protocol and its modelling in more details.
Finally we present theoretical and experimental results that demonstrate the
security of this scheme.

##  I Protocol

####

This protocol is illustrated on Figure 1. A source (for instance a trusted
satellite) emits thermal radiation which is picked by the legal parties and
the eavesdropper. We can consider that Eve can access quite a large portion of
the signal, intercepting much of what should go to Alice and Bob. We model
this by giving Eve an entangling cloner, so she can divert as much of the
signal to her as convenient. However, we consider that the source is trusted;
this means that the eavesdropper does not use the satellite to relay her own
signal.

Figure 1: In this situation, a satellite beams down a signal, which is
received by Alice, Bob and Eve. Eve can have a very large portion of the
signal, but she does not control the signal being emitted.

We express the protocol formally as follows :

  * •

Alice creates a beam from a trusted thermal source.

  * •

On the way to their trusted beamsplitter with transmittance $\eta_{2}$, the
signal is interfered with by Eve, via an entangling cloner denoted $\eta_{1}$.

  * •

Alice uses $\eta_{2}$ to divert part of the signal to her detector and send
the rest on to Bob.

  * •

Similarly to sakuya:2017_2 , the bunched nature of the pairs coming out of
$\eta_{1}$ means that fluctuations present at Alice’s detector are correlated
with those at Bob’s detector.

  * •

To derive their data, Alice and Bob slice these fluctuations as convenient; as
an example, a fluctuation above the mean could be a 1 and a fluctuation below
the mean, a 0.

  * •

Like any QKD scheme, our protocol requires quantum correlations. To confirm
that the signal from Alice and Bob are correlated is done through verifying
the thermal nature of their signal. Thus, Alice sends Bob small chunks of data
for him to perform a $g^{(2)}$ calculation. A $g^{(2)}>1$ means that the
signal is thermal.

  * •

Alice and Bob now have a stream of independent and randomly correlated bits
from which they can derive a key, the security of which they can improve with
Cascade and Advantage Distillation, as per any QKD scheme.

This scheme was implemented as shown on Figure 2. In order to simulate high
levels of noise, we consider two attenuator channels between $\eta_{2}$ and
the legal parties, equivalent to adding a beamsplitter of transmittance
$\eta_{3}$ between $\eta_{2}$ and Alice ($\eta_{4}$ for Bob) , with a input
state of variance $N_{3}$ (and $N_{4}$) at the second input arm.

Figure 2: This is the schematic for this iteration of the protocol. The change
in slight : effectively, Eve and Alice switch places, hence the deceiving
resemblance between the schematics as shown in sakuya:2017_2  and this one.

Once again, this is not a prepare-and-send scheme. Alice controls the source,
but the process of splitting pairs happening at the beamsplitter is
stochastic, therefore unpredictable. Eve has no access to the channels between
$\eta_{2}$ and either Alice or Bob, nor any control over their detectors.

###  I.1 Theoretical modelling

####

Let us recall that thermal states can be modelled using Gaussian statistics,
which makes them easily defined and manipulated through their first and second
moments Eisert:2003 ; GarciaPatron:2007 . The former are contained in the
displacement vector $\left\langle\hat{r}\right\rangle$, where $\hat{r}$ is the
system’s operator, and $\rho$ the state’s density operator. The second moments
are contained in the covariance matrix $\gamma$ defined as

$\gamma_{ij}={\rm
Tr}\left[\rho\left\\{(\hat{r}_{i}-\left\langle\hat{r}_{i}\right\rangle),(\hat{r}_{j}-\left\langle\hat{r}_{j}\right\rangle)\right\\}\rho\right]\,,$

where we write the anti-commutator using $\left\\{\right\\}$. A thermal state
has covariance matrix $\gamma_{in}=2(\bar{n}+1)\boldsymbol{I}$, where
$\bar{n}$ is the average photon number and $\boldsymbol{I}$ the identity
matrix, and null displacement. We use the Bose Einstein distribution

$\bar{n}=\frac{1}{e^{\nicefrac{{\hbar\omega}}{{k_{B}T}}}-1}\,,$ (1)

and consider detectors measuring radiation at $30GHz$ and $T=300K$, so that
$\bar{n}=1309$.

The beamsplitters are modelled as

$\boldsymbol{V}_{i}=\left(\begin{array}[]{cc}\sqrt{\eta_{i}}\boldsymbol{I}&\mu_{i}\boldsymbol{I}\\\
-\mu_{i}\boldsymbol{I}&\sqrt{\eta_{i}}\boldsymbol{I}\end{array}\right)\,,$

where $\eta_{i}$ is the transmittance and $\mu_{i}=\sqrt{1-\eta_{i}}$
represents the noise. They act on the state as
$\gamma_{out}=\boldsymbol{V}\gamma_{in}\boldsymbol{V}^{T}$.

The input state at the first beamsplitter contains the thermal source and
Eve’s source; it has covariance matrix and displacement vector

$\gamma_{in}=\left(\begin{array}[]{cccc}V_{s}^{x}&0&0&0\\\ 0&V_{s}^{p}&0&0\\\
0&0&V_{e}^{x}&0\\\
0&0&0&V_{e}^{p}\end{array}\right)\qquad\textbf{r}_{in}=(x_{s};p_{s};x_{e};p_{e})^{T}\,.$

We note the structure of the covariance matrix as
$\gamma_{in}=\gamma_{source}\bigoplus\gamma_{Eve}$. The two empty sub-matrices
would represent potential pre-existing correlations between the source and
Eve, which in our set-up, is unrealistic.

The output of the second beamsplitter is

$\gamma_{out}=\left(\begin{array}[]{ccc}\tilde{\gamma_{b}}&\tilde{\gamma_{ab}}&\tilde{\gamma_{eb}}\\\
\tilde{\gamma_{ab}}&\tilde{\gamma_{a}}&\tilde{\gamma_{ea}}\\\
\tilde{\gamma_{eb}}&\tilde{\gamma_{ea}}&\tilde{\gamma_{e}}\end{array}\right)\,.$

We make the channel between $\eta_{2}$ and Alice, and between $\eta_{2}$ and
Bob thermal noise channels by inputting states of variance $N_{3}$ on Alice’s
branch and $N_{4}$ on Bob’s as

$N_{i}=\frac{\eta_{i}\chi_{i}}{1-\eta_{i}}\,,\qquad\text{with}\qquad\chi_{i}=\frac{1-\eta_{i}}{\eta_{i}}+\epsilon_{i}\,,$

and $\epsilon_{i}$ the channel excess noise GarciaPatron:2007 . The input
state at $\eta_{3}$ and $\eta_{4}$ is then

$\gamma_{int}=\left(\begin{array}[]{cc}N_{3}&0\\\
0&N_{3}\end{array}\right)\bigoplus\gamma_{out}\bigoplus\left(\begin{array}[]{cc}N_{4}&0\\\
0&N_{4}\end{array}\right)\,,$

where $\gamma_{out}$ is the state at the output of $\eta_{2}$, the first block
sub-matrix is the input state at $\eta_{3}$ and the last sub-matrix, the input
state at $\eta_{4}$.

The output covariance matrix is

$\Gamma_{out}=\left(\begin{array}[]{ccccc}\tilde{\Gamma_{v}}&\tilde{\Gamma_{va}}&\tilde{\Gamma_{ve}}&\tilde{\Gamma_{vb}}&\tilde{\Gamma_{vv^{\prime}}}\\\
\tilde{\Gamma_{va}}&\tilde{\Gamma_{a}}&\tilde{\Gamma_{ea}}&\tilde{\Gamma_{ab}}&\tilde{\Gamma_{av^{\prime}}}\\\
\tilde{\Gamma_{ve}}&\tilde{\Gamma_{ea}}&\tilde{\Gamma_{e}}&\tilde{\Gamma_{eb}}&\tilde{\Gamma_{ev^{\prime}}}\\\
\tilde{\Gamma_{vb}}&\tilde{\Gamma_{ab}}&\tilde{\Gamma_{eb}}&\tilde{\Gamma_{b}}&\tilde{\Gamma_{bv^{\prime}}}\\\
\tilde{\Gamma_{vv^{\prime}}}&\tilde{\Gamma_{av^{\prime}}}&\tilde{\Gamma_{ev^{\prime}}}&\tilde{\Gamma_{bv^{\prime}}}&\tilde{\Gamma_{v^{\prime}}}\end{array}\right)$

where the block sub-matrices are given in the appendix.

####

Maurer and Wolf Maurer:1999  have proved a theorem providing conditions to be
satisfied for a scheme such as ours to be secure. The theorem reads as follows
:

[quote]

######  Theorem 1.

In Scenario 1, the following conditions are equivalent :

  1. 1.

$I(A:B|E)>0$

  2. 2.

$K(A:B\parallel E)>0$

  3. 3.

$I(A:B\downarrow E)>0$

[end quote] 111even though this is a direct quote, we have adapted the
notation to our scheme.

where $K(A:B\parallel E)$ is the secret key rate. The third condition is
actually the most restrictive. $I(A:B\downarrow E)$ is the intrinsic
conditional mutual information; it determines the unreducible amount of
conditional mutual information between Alice and Bob, regardless of any
attemps by Eve at acquiring more information through local operations; in
other words, it is information inaccessible to Eve. Furthermore, it satisfies

$I(A:B\downarrow E)<I(A:B|E)\,,$

which makes it a tighter condition on the secret key rate.

We can see its relation to the quantum discord if we recall that the latter,
$D(B|A)$, is defined as the difference between the mutual information $I(A:B)$
and the classical mutual information $J(B|A)$ (or $J(A|B)$). $I(A:B)$
quantifies all possible correlations between Alice and Bob, but $J(B|A)$
quantifies those measured by local operations at Alice’s and Bob’s sites.
Therefore, it can be understood as the intrinsic conditional mutual
information as described previously. Let us therefore, rewrite the theorem as
:

######  Theorem 2.

In our central broadcast scheme, the following conditions are equivalent:

  1. 1.

$I(A:B|E)>0$

  2. 2.

$K(A:B\parallel E)>0$

  3. 3.

$D(B|A)>0$

It is therefore enough in principle, to demonstrate that either condition is
satisfied. We shall however, prove two, namely the positivity of the
conditional mutual information and that of the discord. The latter will allow
us to demonstrate the quantum nature of the secrecy.

####

The mutual information $I(A:B)$ is given by

$I(A:B)=S(\Gamma_{a})+S(\Gamma_{b})-S(\Gamma_{ab})\,,$

where $S(x)$ is the Von Neumann entropy and $\Gamma_{i}$ the covariance
matrices of A, B and AB respectively. The Von Neumann entropy is given by

$S(x)=\sum_{i=1}^{N}\left(\frac{x_{i}+1}{2}\right)\log\left(\frac{x_{i}+1}{2}\right)-\left(\frac{x_{i}-1}{2}\right)\log\left(\frac{x_{i}-1}{2}\right)$

where $x_{i}$ are the symplectic eigenvalues of $\Gamma$. The discord is
defined explicitely as

$D(B|A)=S(\Gamma_{a})-S(\Gamma_{ab})+\min_{\Gamma_{0}}S(\Gamma_{b|x_{A}})\,$

where $\Gamma_{b|x_{A}}$ is the covariance matrix of B conditionned by a
homodyne measurement on A Weedbrookrmp:2012

$\Gamma_{b|x_{A}}=\,\Gamma_{b}-\Gamma_{ab}(X\Gamma_{a}X)^{-1}\Gamma_{ab}^{T}\,,$

with $X=\left(\begin{array}[]{cc}1&0\\\ 0&0\end{array}\right)$ and $()^{-1}$
the pseudo-inverse. The conditional mutual information is

$I(A:B|E)=S(\Gamma_{ae})+S(\Gamma_{be})-S(\Gamma_{e})-S(\Gamma_{abe})\,.$

##  II Results and discussion

####

The protocol was realised experimentally. The thermal source is provided by a
superluminescent diode coupled to an external cavity, making it a tuneable
laser, run without any added modulation. The laser can be run separately in
coherent or in thermal mode, and the thermality of the source was established
in sakuya:2017_2 . The source bandwidth was measured at $\Delta\lambda=0.4$nm
spread around a centre wavelength of $\lambda_{0}=780.09$nm; this give a
coherence time of $\tau_{c}=4.8$ps. The detectors are ThorLabs Det36A
photodiodes, coupled to a LeCroy Waverunner 44xi oscilloscope; the combined
integration time is 14ns and the oscilloscope samples at 5GSps.

The conditional mutual information is calculated from the sliced data strings
using Shannon entropies $H(x)=-\sum p(x)\log(p(x))$ in terms of the measured
frequencies $p(x)$.

Figure 3 shows that the scheme works experimentally as predicted. $I(A:B|E)$
is best as $\eta_{1}$ tends to 1, and at $\eta_{2}=0.5$, so when Alice and Bob
gets equal shares of most of the thermal source signal. This corresponds to a
situation where the eavesdropper is absent, and where there is minimal loss.
As long as the $\eta_{1}>0.5$, the eavesdropper gets little of the signal and
the advantage is to the legal parties. However, no matter how much signal Eve
receives, the conditional mutual information is always positive, and never
exhibits a sharp fall-off, typical of point-to-point schemes over the $3dB$
limit. This means that it is always possible to build key, albeit slowly.

Figure 3: Conditional mutual information for thermal states (left) versus
coherent states (right). We can see that when $\eta_{1}\rightarrow 1$, so when
there is no amount of Eve coming between the legal parties, the conditional
mutual information peaks.

Figure 3 allows us also to illustrate that this scheme cannot work in the
coherent regime. As mentioned before in sakuya:2017_2 , coherent radiation is
not bunched; therefore, it holds none of the intrinsic correlations contained
in bunched pairs. There is no splitting of pairs occuring at the
beamsplitters, because there are no such pairs; single photons travel through
uncorrelated to Alice and Bob, who as a result can build no key from them.
This is shown on the right-hand graph of the figure. $I(A:B|E)$ remains
constant, no matter how much Eve lets through, no matter the split between
Alice and Bob.

####

Let us now compare these results to those obtained through our theoretical
modelling.

(a)

(b)

Figure 4: We plot the conditional mutual information $I(A:B|E)$ (left) and the
discord $D(B|A)$ (right)against $\eta_{1}$, with $\eta_{2}=0.5$,
$\eta_{3}=\eta_{4}=0.2$ and $\epsilon_{3}=\epsilon_{4}=10^{-2}$. At
$\eta_{2}=0.5$, Alice and Bob share equal part of the signal.

Figure 4a shows the behaviour of the conditional mutual information as Eve
lets more and more of the signal through. The plots match our experimental
results. The higher $\eta_{1}$, the higher $I(A:B:|E)$. Also, since $I(A:B|E)$
is always positive, we conclude there always is secrecy in our scheme.

We can also explore how the initial state of Eve influences the secrecy
between Alice and Bob. For that, we vary $V_{e}$ and see that as it increases,
$I(A:B|E)$ increases also. The reason for this, we have mentioned before and
will detail further in the lines below.

Figure 4b illustrates the positivity of the discord, regardless of $\eta_{1}$.
This means that there always are quantum correlations between Alice and Bob.
This satisfies the third of the conditions from our theorem, and we can affirm
quantum secrecy.

####

What is remarkable is the value of the discord when $\eta_{1}$ is null, so
before Eve begins to let the source signal through. In this case, what is
actually measured is the amount of quantum correlations within Eve’s state. We
have seen that the higher $V_{e}$ is, the higher $I(A:B|E)$, but here we see
that the discord follows a similar trend. This is particularly evident when
$V_{e}=250SNU$.

This is a result of the physics of thermal states. To understand this, let us
step back and consider a single beamsplitter (input arms labelled 1 and 2,
output arms labelled 3 and 4) with a thermal state at one input. Since it is
bunched, there will be correlated photon pairs travelling into the
beamsplitter. If both photons travel into the same input (say arm 1), we can
expect three outputs Loudon:2000  :

  * •

both photons are travelling through on arm 3 $P(2_{3},0_{4})$,

  * •

both photons travel onto arm 4 $P(0_{3},2_{4})$ or

  * •

one photon for each arm $P(1_{3},1_{4})$.

This corresponds to Eve inputting a vacuum or a coherent state at arm 2 and
why we can in fact equate her to any loss in the channel.

On the other hand, if Eve inputs a thermal state as well, there is now a
correlated pair of photons travelling into each input arm. This will gives us
the following outputs: $P(4_{3},0_{4})$, $P(0_{3},4_{4})$, $P(2_{3},2_{4})$,
$P(3_{3},1_{4})$, and $P(1_{3},3_{4})$. The third case $P(2_{3},2_{4})$, is
three-degenerate; either both pairs get to the other side in one piece (which
accounts for two degeneracies) or both pairs are split (the remaining
degeneracy). This means that accounting for all possible outcomes, there are
only two cases where there will not be at least one correlated pair travelling
into $\eta_{2}$ to Alice and Bob: either one pair is split at $\eta_{1}$ and
Eve gets three photons $P(3_{3},1_{4})$ (mitigated by the fact that Eve would
choose to let most of the signal through at $\eta_{1}$ in order not to be
noticed) or both pairs are split at $\eta_{1}$, which is one of the
$P(2_{3},2_{4})$ degeneracies.

If $V_{e}=1SNU$, then Eve inputs a vacuum state, and Alice and Bob build key
solely from the pairs produced at the source. As a result, the discord is
minimal at $\eta_{1}\rightarrow 0$. If $V_{e}>1SNU$, Eve’s state can be
regarded as thermal; in this case, she contributes pairs to those coming from
the source. In fact, if the eavesdropper’s input is too significant, the legal
parties can build a quantum secure key, regardless of how much signal is
coming from the source. As in any QKD, we expect that the eavesdropper will
try to minimise her input, if only to escape detection. At best, she can hope
to merely “listen” in, in which case, her input is $V_{e}=1SNU$. Yet, as soon
as signal begins going through ($\eta_{1}>0.1$), the legal parties can build a
quantum secure key, albeit slowly.

Let us point out that these plots have been obtained for very high level of
noise on Alice’s and Bob’s branches. Indeed $\eta_{3}$ and $\eta_{4}$ are such
that $80\%$ of their signal is lost. Yet, even in this case, the legal parties
are able to construct a quantum secure key.

##  III Concluding remarks

####

In our previous protocol, the security arose from the quantum correlations
within a pair which would split between Alice and Bob. Since Eve placed
herself on the arm going to Bob, she would interfere with/capture photons on
their way to him, but she could not build a three way correlation sufficient
to attack key exchange. In this paper, she places herself prior to the
splitting of the pairs between Alice and Bob, interfering with the pairs
directly from the source. Unlike the situation in sakuya:2017_2 , she is not
limited by the Heisenberg uncertainty principle and can intercept and resend
bunched pairs at her leisure. Therefore, the legal parties cannot distinguish
the pairs coming from her to those coming from the source. This is not as bad
as it sounds, however, because since Eve has no interaction with the output of
$\eta_{2}$, any correlations within the pairs split at $\eta_{2}$ are
completely safe from tampering. So actually, as long as Alice and Bob have a
$g^{(2)}(0)>1$, their pairs, regardless of their origins, are correlated, and
quantum secrecy is possible sakuya:2017_2 .

This is the strength of this scheme. Even if Eve succeeds in hiding in the
noise, if her input is not either vacuum or a perfect coherent state, she will
contribute correlations to the pool which Alice and Bob can build key from,
but she cannot know when or if these injected states have contributed to the
key. Another option for Eve is to actually become the source; we explore this
in a forthcoming publication.

This experiment was carried out at optical frequencies using a pseudo thermal
source, however, the theoretical modelling was performed at values of
$\bar{n}$ consistent with the microwave regime. Interferometeres used in radio
astronomy rely on the presence of thermal correlations being preserved over
astronmical distances, and as the results in this paper suggest that the
results are highly portable to the microwave regime. Hence, this method of key
exchange appears to be a viable option for long distance key exchange.

####

The authors are grateful to network collaborators J. Rarity, S. Pirandola, C.
Ottaviani, T. Spiller, N. Luktenhaus and W. Munro for very fruitful
discussions. This work was supported by funding through the EPSRC Quantum
Communications Hub EP/M013472/1 and additional funding for F.W. from Airbus
Defense & Space.

####

Data that support the findings of this study are available from the Research
Data Leeds Repository with the identifier https://doi.org/10.5518/587
Data:QSTS2 .

## References

  * [1] Paul Marks.  Quantum cryptography to protect swiss election.  News Scientist, pages https://institutions.newscientist.com/article/dn12786–quantum–cryptography–to–protect–swiss–election/, 2007\.
  * [2] Elizabeth Newton, Anne Ghesquière, Freya L. Wilson, Benjamin T. H. Varcoe, and Martin Moseley.  Quantum secrecy in thermal states.  Journal Physics B, 52(12):125501, 2019.
  * [3] Bing Qi, Philip G. Evans, and Warren P. Grice.  Passive state preparation in the gaussian-modulated coherent-states quantum key distribution.  Physical Review A, 97:012317, 2018.
  * [4] Christian Weedbrook, Stefano Pirandola, Seth Lloyd, and Timothy C. Ralph.  Physical Review Letters, 105:110501, 2010.
  * [5] Christian Weedbrook, Stefano Pirandola, and Timothy C. Ralph.  Physical Review A, 86(2):022318(12), 2012.
  * [6] Sammy Ragy and Gerardo Adesso.  Physica Scripta, 2013(T153):014052, 2013.
  * [7] Stefano Pirandola.  Scientific Reports, 4:6956, 2014.
  * [8] Jens Eisert and Martin Plenio.  Int. J. Quant. Inf., 1:479, 2003.
  * [9] Raúl García-Patrón Sanchez.  PhD thesis, Université Libre de Bruxelles, 2007.
  * [10] Ueli M. Maurer and Stefan Wolf.  IEEE Transactions on Information Theory, 45(2):499–514, 1999.
  * [11] Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J. Cerf, Timothy C. Ralph, Jeffrey H. Shapiro, and Seth Lloyd.  Review of Modern Physics, 84(2):621(49), 2012.
  * [12] Rodney Loudon.  The Quantum Theory of Light.  Oxford University Press, 3rd edition, 2000.
  * [13] Elizabeth Newton.  Data for Thermal State QKD, https://doi.org/10.5518/587.  University of Leeds, 2019.

##  Appendix A Protocol 2

###  A.1 After $\eta_{2}$

The submatrices are as follows

$\displaystyle\tilde{\gamma_{b}}=$
$\displaystyle\left(\begin{array}[]{cc}\eta_{2}+\mu_{2}^{2}(\eta_{1}V_{s}^{x}+\mu_{1}^{2}V_{e}^{x})&0\\\
0&\eta_{2}+\mu_{2}^{2}(\eta_{1}V_{s}^{p}+\mu_{1}^{2}V_{e}^{p})\end{array}\right)$
(4) $\displaystyle\tilde{\gamma_{a}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{2}^{2}+\eta_{2}(\eta_{1}V_{s}^{x}+\mu_{1}^{2}V_{e}^{x})&0\\\
0&\mu_{2}^{2}+\eta_{2}(\eta_{1}V_{s}^{p}+\mu_{1}^{2}V_{e}^{p})\end{array}\right)$
(7) $\displaystyle\tilde{\gamma_{e}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{1}^{2}V_{s}^{x}+\eta_{1}V_{e}^{x}&0\\\
0&\mu_{1}^{2}V_{s}^{p}+\eta_{1}V_{e}^{p}\end{array}\right)$ (10)
$\displaystyle\tilde{\gamma_{ea}}=$
$\displaystyle\left(\begin{array}[]{cc}-\mu_{1}\sqrt{\eta_{1}}\sqrt{\eta_{2}}(V_{s}^{x}-V_{e}^{x})&0\\\
0&-\mu_{1}\sqrt{\eta_{1}}\sqrt{\eta_{2}}(V_{s}^{p}-V_{e}^{p})\end{array}\right)$
(13) $\displaystyle\tilde{\gamma_{eb}}=$
$\displaystyle\left(\begin{array}[]{cc}-\mu_{1}\sqrt{\eta_{1}}\mu_{2}(V_{s}^{x}-V_{e}^{x})&0\\\
0&-\mu_{1}\sqrt{\eta_{1}}\mu_{2}(V_{s}^{p}-V_{e}^{p})\end{array}\right)$ (16)
$\displaystyle\tilde{\gamma_{ab}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{2}\sqrt{\eta_{2}}(\eta_{1}V_{s}^{x}+\mu_{1}^{2}V_{e}^{x}-1)&0\\\
0&\mu_{2}\sqrt{\eta_{2}}(\eta_{1}V_{s}^{p}+\mu_{1}^{2}V_{e}^{p}-1)\end{array}\right)$
(19)

###  A.2 Ater $\eta_{3}$ and $\eta_{4}$

The submatrices are as follows

$\displaystyle\tilde{\Gamma_{e}}=$
$\displaystyle\left(\begin{array}[]{cc}\left\langle\tilde{X_{e}}^{2}\right\rangle&0\\\
0&\left\langle\tilde{P_{e}}^{2}\right\rangle\end{array}\right)\,,\quad\tilde{\Gamma_{a}}=\left(\begin{array}[]{cc}\mu_{3}^{2}N_{3}+\eta_{3}\left\langle\tilde{X_{a}}^{2}\right\rangle&0\\\
0&\mu_{3}^{2}N_{3}+\eta_{3}\left\langle\tilde{P_{a}}^{2}\right\rangle\end{array}\right)$
(24) $\displaystyle\tilde{\Gamma_{b}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{4}^{2}N_{4}+\eta_{4}\left\langle\tilde{X_{b}}^{2}\right\rangle&0\\\
0&\mu_{4}^{2}N_{4}+\eta_{4}\left\langle\tilde{P_{b}}^{2}\right\rangle\end{array}\right)\,,\quad\tilde{\Gamma_{v}}=\left(\begin{array}[]{cc}\eta_{3}N_{3}+\mu_{3}^{2}\left\langle\tilde{X_{a}}^{2}\right\rangle&0\\\
0&\eta_{3}N_{3}+\mu_{3}^{2}\left\langle\tilde{P_{a}}^{2}\right\rangle\end{array}\right)$
(29) $\displaystyle\tilde{\Gamma_{v^{\prime}}}=$
$\displaystyle\left(\begin{array}[]{cc}\eta_{4}N_{4}+\mu_{4}^{2}\left\langle\tilde{X_{b}}^{2}\right\rangle&0\\\
0&\eta_{4}N_{4}+\mu_{4}^{2}\left\langle\tilde{P_{b}}^{2}\right\rangle\end{array}\right)$
(32) $\displaystyle\tilde{\Gamma_{ea}}=$
$\displaystyle\left(\begin{array}[]{cc}\sqrt{\eta_{3}}\left\langle\tilde{X_{a}}\tilde{X_{e}}\right\rangle&0\\\
0&\sqrt{\eta_{3}}\left\langle\tilde{P_{a}}\tilde{P_{e}}\right\rangle\end{array}\right)\,\quad\tilde{\Gamma_{eb}}=\left(\begin{array}[]{cc}\sqrt{\eta_{4}}\left\langle\tilde{X_{b}}\tilde{X_{e}}\right\rangle&0\\\
0&\sqrt{\eta_{4}}\left\langle\tilde{P_{b}}\tilde{P_{e}}\right\rangle\end{array}\right)$
(37) $\displaystyle\tilde{\Gamma_{ab}}=$
$\displaystyle\left(\begin{array}[]{cc}\sqrt{\eta_{3}}\sqrt{\eta_{4}}\left\langle\tilde{X_{a}}\tilde{X_{b}}\right\rangle&0\\\
0&\sqrt{\eta_{3}}\sqrt{\eta_{4}}\left\langle\tilde{P_{a}}\tilde{P_{b}}\right\rangle\end{array}\right)\,,\quad\tilde{\Gamma_{vv^{\prime}}}=\left(\begin{array}[]{cc}-\mu_{3}\mu_{4}\left\langle\tilde{X_{a}}\tilde{X_{b}}\right\rangle&0\\\
&-\mu_{3}\mu_{4}\left\langle\tilde{P_{a}}\tilde{P_{b}}\right\rangle\end{array}\right)$
(42) $\displaystyle\tilde{\Gamma_{vb}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{3}\sqrt{\eta_{4}}\left\langle\tilde{X_{a}}\tilde{X_{b}}\right\rangle&0\\\
0&\mu_{3}\sqrt{\eta_{4}}\left\langle\tilde{P_{a}}\tilde{P_{b}}\right\rangle\end{array}\right)\,,\quad\tilde{\Gamma_{av^{\prime}}}=\left(\begin{array}[]{cc}-\sqrt{\eta_{3}}\mu_{4}\left\langle\tilde{X_{a}}\tilde{X_{b}}\right\rangle&0\\\
0&-\sqrt{\eta_{3}}\mu_{4}\left\langle\tilde{P_{a}}\tilde{P_{b}}\right\rangle\end{array}\right)$
(47) $\displaystyle\tilde{\Gamma_{va}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{3}\sqrt{\eta_{3}}\left(\left\langle\tilde{X_{a}}^{2}\right\rangle-
N_{3}\right)&0\\\
0&\mu_{3}\sqrt{\eta_{3}}\left(\left\langle\tilde{P_{a}}^{2}\right\rangle-
N_{3}\right)\end{array}\right)$ (50)
$\displaystyle\tilde{\Gamma_{bv^{\prime}}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{4}\sqrt{\eta_{4}}\left(N_{4}-\left\langle\tilde{X_{b}}^{2}\right\rangle\right)&0\\\
0&\mu_{4}\sqrt{\eta_{4}}\left(N_{4}-\left\langle\tilde{P_{b}}^{2}\right\rangle\right)\end{array}\right)$
(53) $\displaystyle\tilde{\Gamma_{ve}}=$
$\displaystyle\left(\begin{array}[]{cc}\mu_{3}\left\langle\tilde{X_{a}}\tilde{X_{e}}\right\rangle&0\\\
0&\mu_{3}\left\langle\tilde{P_{a}}\tilde{P_{e}}\right\rangle\end{array}\right)\,,\quad\tilde{\Gamma_{ev^{\prime}}}=\left(\begin{array}[]{cc}-\mu_{4}\left\langle\tilde{X_{b}}\tilde{X_{e}}\right\rangle&0\\\
0&-\mu_{4}\left\langle\tilde{P_{b}}\tilde{P_{e}}\right\rangle\end{array}\right)$
(58)


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

# Alchemical Transfer Approach to Absolute Binding Free Energy Estimation

Joe Z. Wu  Department of Chemistry, Brooklyn College of the City University of
New York, New York, NY  Solmaz Azimi  Department of Chemistry, Brooklyn
College of the City University of New York, New York, NY  Sheenam Khuttan
Department of Chemistry, Brooklyn College of the City University of New York,
New York, NY  Nanjie Deng  Department of Chemistry and Physical Sciences, Pace
University, New York, New York, 10038  Emilio Gallicchio
<EMAIL_ADDRESS>Department of Chemistry, Brooklyn College of the
City University of New York, New York, NY

###### Abstract

The Alchemical Transfer Method (ATM) for the calculation of standard binding
free energies of non-covalent molecular complexes is presented. The method is
based on a coordinate displacement perturbation of the ligand between the
receptor binding site and the explicit solvent bulk, and a thermodynamic cycle
connected by a symmetric intermediate in which the ligand interacts with the
receptor and solvent environments with equal strength. While the approach is
alchemical, the implementation of ATM is as straightforward as for physical
pathway methods of binding. The method is applicable in principle with any
force field, it does not require splitting the alchemical transformations into
electrostatic and non-electrostatic steps, and it does not require soft-core
pair potentials. We have implemented ATM as a freely available and open-source
plugin of the OpenMM molecular dynamics library. The method and its
implementation are validated on the SAMPL6 SAMPLing host-guest benchmark set.
The work paves the way to streamlined alchemical relative and absolute binding
free energy implementations on many molecular simulation packages and with
arbitrary energy functions including polarizable, quantum-mechanical, and
artificial neural network potentials.

##  1 Introduction

The binding free energy of a molecular complex is a rigorous thermodynamic
measure of the degree of affinity of two molecules for each other.
Measurements of binding free energies are useful in a wide range of chemical
and medicinal applications ranging from drug discovery to chemical detection
and toxicology. The ability to estimate binding free energies by computational
modeling adds an important dimension to this probe.1 Relative binding free
energy models are widely used, for example, in drug lead optimization.2

This work is concerned in particular with atomistic models of the absolute
binding free energy. These methods can be divided into two classes that are
based on the nature of the thermodynamic path used to connect the bound and
unbound states of the molecular complex.3 Physical pathway methods define a
spatial coordinate along which the reversible work for bringing the two
molecules together is calculated. Conversely, alchemical methods connect the
bound and unbound states by a series of unphysical intermediate states. The
Single-Decoupling Method (SDM),4, 5 for example, is based on progressively
turning on the effective interaction between the ligand and the receptor with
an implicit representation of the solvent. The related Double-Decoupling
Method (DDM),6 which is widely used to estimate the absolute binding free
energy with an explicit representation of the solvent, estimates the binding
free energy as the difference between the free energies of coupling the ligand
to a pure solvent system and to the solvated receptor.

Physical and alchemical binding free energy methods are characterized by
distinct challenges and limitations. While it is broadly employed, the
alchemical DDM approach is known to suffer from poor convergence and strong
bias, especially for large and charged ligands, and these consequences are
characterized by large and compensating decoupling free energies that result
in the amplification of statistical and systematic errors in the binding free
energy estimate.7 Another significant shortcoming of DDM is the large
perturbation of the environment of the ligand in going from the solvated
states to the vacuum intermediate state. Unless the ligand is conformationally
restrained, the transition from a solvated state to vacuum can trigger
substantial intramolecular conformational changes that relax slowly to the
bound or unbound configurations. Other technical limitations of DDM stem from
differences in the composition and size of the molecular systems used for the
two decoupling legs,8 and from the inconsistent treatment of long-range
interactions.9, 10 Even for relatively simple systems, these and other
alchemical transformations have to be conducted with care to avoid
singularities and slow convergence. It is recommended for example, to couple
electrostatic interactions separately from other interactions,11 and to employ
customized soft-core pair potentials12 to avoid end-point singularities.

Because they are often based on modifying the parameters and the form of the
energy function, software implementations of alchemical binding free energy
methods also tend to be significantly complex and require in-depth knowledge
of the data structures of the target molecular simulation package. For
example, core energy routines are usually customized to implement the specific
modified pair potentials that represent the interaction of the ligand with the
rest of the system. These modified pair potentials depend on the alchemical
progress parameter (generally denoted by $\lambda$) which, together with other
alchemical variables, become additional system parameters. Methods such as
Thermodynamic Integration (TI)13, 14 require additional routines to implement
the calculations of the gradients of the energy function with respect to the
parameters that are alchemically transformed. The implementation of alchemical
transformations involving many-body potential terms, where the $\lambda$
dependence affects more than individual pair interaction energies, is
particularly challenging. These include alchemical applications with
polarizable potentials,15 Ewald long-range electrostatic treatments,16
implicit solvent models,17 as well as some conventional intramolecular
potential terms.18

Physical pathway methods19, 20, 7 address some the limitations of DDM by
physically moving the ligand from the solvent bulk to the binding site.
Physical pathway binding free energy calculations are seen as preferable over
DDM for large and charged ligands because they are typically performed as one
continuous path in a single solvent box without transferring the ligand to a
different phase. In addition, software implementations of physical pathway
methods do not require as much customization of the underlying molecular
simulation package as alchemical methods do. The primary limitation of
physical pathway methods is the high computational cost due to the requirement
of equilibrating the complex at many intermediate receptor-ligand separations
that might not be of interest. Because they require a physical exit and entry
channel, physical pathway methods are also not generally applicable to
occluded binding sites.21

Building upon on our Single Decoupling Method (SDM) for absolute binding free
energy estimation22, 23 implemented in OpenMM,24 we have been investigating
ways to streamline alchemical calculations with explicit solvation. SDM, which
has been designed for binding free energy calculations with implicit
solvation, computes the alchemical perturbation energy by translating the
ligand from the solvent medium to the receptor binding site, rather than
attempting to selectively turn-on and turn-off individual ligand-receptor
interactions. The approach treats all interactions in one concerted step and
employs the standard molecular mechanics force field without soft-core pair
potentials. End-point singularities are addressed by a suitable non-linear
alchemical energy function.23 We have recently shown that the same approach is
applicable to the estimation of the concerted hydration free energies of drug-
sized solutes in water droplets with explicit solvation.25

In this work, we extend this concerted alchemical scheme to the calculations
of absolute binding free energies with explicit solvation. In the resulting
alchemical scheme, called the Alchemical Transfer Method (ATM) (Figure 1B),
the unbound and bound states of the molecular complex are related by a
translation vector that brings the ligand from the solvent bulk to the
receptor binding site in a single solvent box. We show that the proposed
method addresses some of the aforementioned challenges of binding free energy
calculations by exploiting the best characteristics of the alchemical and
physical methods. Like DDM, the method is based on alchemical transformations
that aim conformational sampling only in the solvent bulk and the receptor
binding site, and like tphysical methods, ATM is based on moving the ligand in
physical space, in a single simulation box, and without transferring the
ligand to vacuum. The ATM method, implemented as a freely available plugin of
the OpenMM molecular simulation package, did not require any modifications of
the OpenMM core energy routines. We validate the ATM approach on a rigorous
benchmarking dataset developed by Rizzi et al.26

##  2 Theory and Methods

###  2.1 Alchemical Transformations

Free energy changes are estimated using alchemical transformations based on a
$\lambda$-dependent potential energy functions $U_{\lambda}(x)$ that brings
the system from an initial state at $\lambda=0$, described by the potential
function $U_{0}(x)$, to a final state at $\lambda=1$, corresponding to the
potential function $U_{1}(x)$.

For each transformation, the alchemical potential energy is expressed as

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

where $x$ represents the set of atomic coordinates of the system,

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

is the perturbation energy, and $W_{\lambda}(u)$ is the generalized softplus
alchemical perturbation function

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

The parameters $\lambda_{2}$, $\lambda_{1}$, $\alpha$, $u_{0}$, and $w_{0}$
are functions of $\lambda$ (see Computational Details), and

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

with

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

and

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

is the soft-core perturbation energy. The soft-core function is monotonic map
that avoids singularities near the initial state of the alchemical
transformation at $\lambda=0$ without affecting the distribution of
perturbation energies at the final state at $\lambda=1$. As expressed in Eq.
(4), the soft-core perturbation energy is designed to cap the perturbation
energy $u(x)$ to a maximum positive value $u_{\rm max}$ and to be equal to the
perturbation energy when this is below a cutoff value $u_{0}$. The $u_{0}$
cutoff is selected to be sufficiently large so that $u_{\rm sc}(u)=u$ for all
observed samples collected at the end state. The specific values of $u_{0}$,
$u_{\rm max}$ and of the scaling parameter $a$ used in this work are listed in
Computational Details.

In order to reproduce the desired end points, it is necessary that the
alchemical perturbation function is defined such that $W_{0}(u)=0$ and
$W_{1}(u)=u$ at $\lambda=0$ and $\lambda=1$, respectively. This requirement is
satisfied by the linear function, $W_{\lambda}(u)=\lambda u_{\rm sc}(u)$,
which is special case of the softplus function in Eq. (3) for which
$\lambda_{1}=\lambda_{2}=\lambda$. The linear function is the standard choice
for the alchemical perturbation function. As it can be verified from Eq. (3),
in general the end-point requirement is satisfied whenever
$\lambda_{1}=\lambda_{2}=0$ at $\lambda=0$ and $\lambda_{1}=\lambda_{2}=1$ at
$\lambda=1$.

###  2.2 The Alchemical Transfer Method for Binding Free Energy Estimation

Consider the non-covalent binding process between a receptor R and a ligand L.
The standard free energy of binding, $\Delta G^{\circ}_{b}$, is defined as the
difference in free energy between the bound complex and the unbound components
at the standard concentration of $C^{\circ}=1$ M,

$\Delta G^{\circ}_{b}=\Delta G^{\circ}_{\rm site}+\Delta G^{\ast}_{b}.$ (7)

where $\Delta G^{\ast}_{b}$ is the excess component, defined as the reversible
work for transferring the ligand to the binding site region of the receptor of
volume $V_{\rm site}$ from a region of the same volume in the solvent bulk
(Figure 1) plus a concentration-dependent term

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

that corresponds to the free energy of transfer of a ligand molecule from an
ideal solution at concentration $C^{\circ}$ to a region of volume $V_{\rm
site}$ in the solvent.3 In the remainder, we will focus on the calculation of
the excess free energy component $\Delta G^{\ast}_{b}$ by alchemical molecular
simulations.

(A) (B)

Figure 1: Schematic illustrations of (A) the double decoupling (DDM) and (B)
the Alchemical Transfer (ATM) methods for the calculation of the binding free
energy between a molecular receptor (orange doughnut) and a ligand (black
circle). The dashed circle within the receptor represents the binding site
region. The blue boxes represent the solvent. The unbound and bound end states
for the two methods are considered thermodynamically equivalent. In both
methods, the end states are transformed to a common intermediate state and the
excess binding free energy is the difference of the free energy changes of the
two legs, $\Delta G^{\ast}_{b}=\Delta G_{2}-\Delta G_{1}$. Double decoupling
defines an intermediate state in which the ligand is in vacuum (white). The
alchemical transformations in the Alchemical Transfer Method are instead
performed in the same solvent box and the ligand does not leave the solvated
system.

The Double Decoupling Method (DDM)6 has emerged as one of the gold standards
for the alchemical calculation of absolute binding free energies in condensed
phases. DDM is based on the thermodynamic cycle illustrated in Fig. 1A,
whereby the bound and unbound states of the molecular complex are
thermodynamically linked by an intermediate state in which the ligand is
placed in vacuum. The excess binding free energy is expressed as the
difference of the free energies of alchemically decoupling the ligand from the
unbound (leg 1) and bound (leg 2) states into the intermediate vacuum state,
as

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

The Alchemical Transfer Method (ATM) proposed here avoids the vacuum
intermediate and requires only one molecular system. As illustrated in Fig.
1B, ATM is based on an alchemical intermediate in which the ligand interacts
simultaneously with the solvent bulk and the receptor. We assume, without loss
of generality, that there is a suitable coordinate frame attached to the
receptor and that the binding site region (represented by the dashed circle in
Fig. 1) is fixed relative to this coordinate frame. Under these assumptions,
every point in the binding site region maps to a unique point into a region of
the same shape in the solvent bulk by means of a constant translation vector
$h$, also at rest relative to the receptor coordinate frame (Figure 2). The
bound state of the system is defined as any configuration in which the ligand
is in the binding site region. Conversely, the unbound state of the system is
defined as any configuration of the system in which the ligand is placed into
the bulk solvent region (Fig. 1B).

Under these assumptions, any configuration of the bound state maps to an
unique configuration of the unbound state by a rigid translation of the ligand
atoms by the vector $h$ (Figure 2). The reverse is also true. Any
configuration of the unbound system maps uniquely to a configuration of the
bound state by translation of the vector $-h$. Hence, the translation vector
$h$ can be used as a perturbation parameter to connect, in a statistical
thermodynamic sense, the bound and unbound states.

Figure 2: Illustration of the unbound state of the complex between the CB8
host (center) and the G3 guest (lower left). The small red spheres represent
the oxygen atoms of the water molecules. The displacement vector that
translates the ligand from the binding site to the bulk solvent position is
indicated.

For example, consider leg 1 of the ATM cycle in Fig. 1B. Denoting the
potential energy function of the system as $U(x)$, with $x=(x_{S},x_{L})$
being the coordinates of the bound system, where $x_{L}$ are the coordinates
of the ligand in the receptor binding site and $x_{S}$ are the coordinates of
the receptor and the solvent (the surroundings), the potential energy
functions $U_{0}$ and $U_{1}$ of the initial and final states of the leg (the
bound and unbound states in this case) are respectively

$U_{0}(x_{S},x_{L})=U(x_{S},x_{L})$ (10)

and

$U_{1}(x_{S},x_{L})=U(x_{S},x_{L}+h)\,.$ (11)

To connect the bound and unbound states, we consider the hybrid alchemical
potential of Eq. (1) with the perturbation energy

$u_{1}(x_{S},x_{L})=U(x_{S},x_{L}+h)-U(x_{S},x_{L}),$ (12)

which is defined as the change in the potential energy of the system for
rigidly translating the ligand atoms from the binding site region to the bulk
solvent region while all other degrees of freedom of the system remain
unchanged.

With these definitions, the alchemical hybrid potential used for the first leg
of the ATM cycle is:

${\rm Leg\ 1:}\quad
U_{\lambda}(x_{S},x_{L})=U(x_{S},x_{L})+W_{\lambda}(u_{1})\,,\quad
0<=\lambda<=1/2$ (13)

where the alchemical perturbation function $W_{\lambda}(u)$ is defined in Eq.
(3) and the perturbation energy $u_{1}$ is defined in Eq. (12). As indicated
in Eq. (13), the alchemical pathway for the first leg is terminated at
$\lambda=1/2$, where the ligand interacts with half strength with both the
receptor environment and the solvent bulk. This ensures that severe steric
clashes are not likely to occur at $\lambda=1/2$. The $u_{0}$ parameter of the
soft-core perturbation potential is set to a large enough value so that the
perturbation potential $u$ does not exceed $u_{0}$ at $\lambda=1/2$. Under
these conditions, the original and the soft-core perturbation potentials are
the same and it follows from Eq. (13) that the alchemical potential energy at
$\lambda=1/2$ is

$U_{1/2}(x_{S},x_{L})=\frac{1}{2}U(x_{S},x_{L})+\frac{1}{2}U(x_{S},x_{L}+h).$
(14)

Eq. (14) defines the potential energy function of the alchemical intermediate
of the ATM cycle in Fig. 1B. Accordingly, in the alchemical intermediate
ensemble the ligand interacts symmetrically with the receptor and bulk solvent
environments. Conversely, in this state the receptor atoms and the solvent
molecules interact with the ligand at half strength. The alchemical
calculation that corresponds to the alchemical potential Eq. (13) yields the
free energy change $\Delta G_{1}$ in going from the bound state to the ATM
alchemical intermediate.

To formulate the second leg of the ATM cycle connecting the unbound state to
the intermediate state (Figure 1), the role of the end states is reversed
relative to the first leg. As before, $x_{L}$ describes the coordinates of the
ligand in the receptor binding site and $x_{L}+h$ describes the coordinates of
the ligand in the solvent bulk and the initial and final states are now
defined by the potential energy functions of, respectively, the unbound and
bound states:

$U_{0}(x_{S},x_{L})=U(x_{S},x_{L}+h)$ (15)

and

$U_{1}(x_{S},x_{L})=U(x_{S},x_{L})\,.$ (16)

The alchemical potential energy function for the second leg of the ATM cycle
is

${\rm Leg\ 2:}\quad
U_{\lambda}(x_{S},x_{L})=U(x_{S},x_{L}+h)+W_{\lambda}(u_{2})\,,\quad
0<=\lambda<=1/2$ (17)

with the perturbation energy

$u_{2}(x_{S},x_{L})=U(x_{S},x_{L})-U(x_{S},x_{L}+h)$ (18)

that corresponds to the potential energy change of transferring the ligand
from the bulk to the receptor binding site. Note that under the same
assumptions that led to Eq. (14), Eq. (17) reaches $\lambda=1/2$ at the same
alchemical intermediate as Eq. (13). Thus, Eq. (9) holds for the ATM
thermodynamic cycle in Fig. 1B.

We end the presentation of the Alchemical Transfer Method by discussing the
requirement of splitting the alchemical path into two legs. Although the
symmetric hybrid potential Eq. (14) can formally cover the direct path from
the unbound to the bound states by extending the $\lambda$ range from $0$ to
$1$, in practice it suffers severe end-point singularities at both of the end
states. The potential energy, $U(x_{S},x_{L}+h)$, when the ligand is placed in
the solvent bulk region, is ill-defined near $\lambda=0$ due to the atomic
clashes that occur when ligand-solvent interactions are turned off.
Conversely, near $\lambda=1$, the potential energy $U(x_{S},x_{L})$, when the
ligand is placed in the receptor site, is ill-defined due to clashes between
ligand atoms and receptor atoms. The soft-core protocol we employ is based on
an asymmetric definition of the perturbation potential [see Eqs. (12) and
(18)] that can address the singularity at one end-point or the other but not
both simultaneously with only one continuous alchemical perturbation potential
energy function.

###  2.3 Software Implementation

The method is implemented as an integrator plugin
(github.com/rajatkrpal/openmm_sdm_plugin)23 of the OpenMM library.24 The
integrator is based on the Langevin thermostat and high-level routines that
displace the ligand, issue calls to energy and forces calculation routines,
and compute the alchemical potential energy [Eq. (3)] and its gradients by
combining the returned system energies and forces. For example for leg 1, at
each MD-step the plugin first computes and saves in temporary buffers the
potential energy and the forces when the ligand is in the binding site. Then
it displaces the ligand in the bulk by the displacement vector $h$ and
recalculates the energy and forces (except for the binding site restraining
potential, see below). The perturbation energy $u_{1}$ [Eq. (12)] and its
gradients are obtained by taking the corresponding differences before and
after the displacement. The gradients of the alchemical perturbation energy
(13) are derived from those of $u_{1}$ and of the undisplaced potential
$U(x_{S},x_{L})$ by application of the chain rule. The resulting forces are
then used to propagate atomic coordinates by one MD step. The same process is
used for leg 2, except that the ligand is initially placed in the solvent bulk
and it is translated into the binding site by reversing the displacement. In
each case the binding restraint potential (see below) is applied when the
ligand is in the binding site. No modifications of the core OpenMM energy
routines are applied.

###  2.4 Benchmark Systems

The benchmark systems here were drawn from the host-guest systems presented in
the SAMPL6 SAMPLing challenge26. The octa-acid (OA) and cucurbit[8]uril (CB8)
hosts are well-studied supramolecular systems that have been featured in
previous host-guest binding SAMPL challenges and the three guests resemble
conventional druglike small molecules and fragments. In total, the benchmark
systems presented here include 5-hexenoic acid (OA-G3) and 4-methylpentanoic
acid (OA-G6) for the OA host, and quinine (CB3-G3) for the CB3 host. Despite
the name, the guest G3 in CB8-G3 is distinct from the G3 guest in OA-G3. The
parametrized systems, including their solvent descriptions, partial charges,
and initial geometries, are provided at the github SAMPL6 site:
github.com/samplchallenges/SAMPL6/tree/master/host_guest/SAMPLing.

Host-guest systems are practical alternatives to protein-ligand systems
because of their minimal atom count and the improbability of undergoing major
conformational reorganization. These systems thus enable the investigation of
novel simulation techniques, such as longer timescales and appropriate
sampling of multiple equivalent binding modes. Each host-guest complex has
five conformations that differ in the position within the binding site, as
well as in torsion angles. The procedure of obtaining five replicate free
energy calculations allows for assessing the statistical uncertainty and
reproducibility of the methodologies. The benchmark systems employed here
create an accessible platform to improve both the predictive accuracy and
computational efficiency of free energy calculations.26

###  2.5 Computational Details

The ATM calculations employed the host-guest molecular complexes parametrized
using the GAFF1.8/AM1-BCC force-field and the TIP3P water model as provided by
Rizzi et al.26 from
(https://github.com/samplchallenges/SAMPL6/tree/master/host_guest/SAMPLing).
The Cartesian components of the displacement vector $h$ were set to
approximately half the dimensions of the simulation box in order to ensure
that the ligand is placed in the corner of the solvent box (Figure 2). This
position is the farthest from the host, which the solvent box is centered
around, and its periodic images. The complexes were energy minimized and
thermalized at 300 K. Then, using the ATM alchemical potential energy function
for leg 1 [Eq. (13)] and starting at the bound state at $\lambda=0$, the
systems were annealed to the symmetric intermediate $\lambda=1/2$ for $250$
ps. The purpose of this step is to obtain a suitable initial configuration
without severe unfavorable repulsion interactions at either end of the
alchemical paths in order to start the molecular dynamics replica exchange
alchemical calculations for each leg (see below). To limit the fluctuations of
the position of the bulk solvent region, which would impact convergence, the
position and orientation of the receptor was loosely restrained with a flat-
bottom harmonic potential of force constant 25.0 kcal/(mol Å2) and a tolerance
of 1.5 Å on all of the heavy atoms of the lower portion of the receptor (the
first 40 atoms of the host as listed in the provided files).

Polar hydrogen atoms with zero Lennard-Jones parameters were modified to
$\sigma_{\rm LJ}=0.1$ Å and $\epsilon_{\rm LJ}=10^{-4}$ kcal/mol to avoid
large attractive interactions between opposite charges at small distances in
nearly uncoupled states. The change in potential energy of the system in the
unbound, bound, and symmetric intermediate states due to this modification of
the Lennard-Jones parameters is below single floating point precision. Single
Decoupling alchemical calculations were prepared using the SDM workflow
(github.com/egallicc/openmm_sdm_workflow.git). MD calculations employed the
OpenMM24 MD engine and the SDM integrator plugins
(github.com/rajatkrpal/openmm_sdm_plugin.git) using the OpenCL platform. The
ASyncRE software,27 customized for OpenMM and SDM
(github.com/egallicc/async_re-openmm.git), was used for the Hamiltonian
Replica Exchange in $\lambda$ space for each ATM leg.

The linear alchemical perturbation potential, $W_{\lambda}(u)=\lambda u_{\rm
sc}(u)$, corresponding to Eq. (3) with $\lambda_{1}=\lambda_{2}=\lambda$, was
used for the octa-acid (OA) systems with $11$ $\lambda$-states uniformly
distributed between $\lambda=0$ and $1/2$ for each of the two ATM legs. The
ATM calculations for the more challenging CB8-G3 complex employed the softplus
perturbation potential Eq. (3) with $24$ $\lambda$-states and the parameters
listed in Tables 1 and 2 for each leg. The parameters of the softplus
perturbation potential were optimized using the scheme described previously23,
25 which involved running trial calculations with the linear potential,
fitting the analytical model of binding28 to each transformation, and
adjusting the parameters of the softplus potential to the resulting
$\lambda$-function23 to obtain a smooth alchemical transition.

The soft-core perturbation energy Eq. (4) was used for all calculations with
$u_{\rm max}=300$ kcal/mol, $u_{0}=100$ kcal/mol. The ligand was sequestered
within the binding site by means of a flat-bottom harmonic potential between
the centers of mass of the host and the ligand with a force constant of $25$
kcal/mol Å2 applied for separation greater than $4.5$ Å. Perturbation energy
samples and trajectory frames were saved every 5 ps. Hamiltonian replica
exchanges in $\lambda$-space were performed every 5 ps. The Langevin
thermostat with a time constant of 2 ps was used to maintain the temperature
at 300 K. Each replica was simulated for a minimum of 10 ns. Binding free
energies and the corresponding uncertainties were computed from the
perturbation energy samples using UWHAM29, discarding the first 5 ns of
trajectory, followed by the addition of the concentration-dependent term
$\Delta G^{\circ}_{\rm site}=-k_{B}T\ln C^{\circ}V_{\rm site}=0.87$ kcal/mol
that corresponds to $300$ K temperature and the volume $V_{\rm site}$ of a
sphere of radius $4.5$ Å. The replica exchange simulations were run on the
XSEDE Comet GPU HPC cluster at the San Diego Supercomputing Center each using
four GPUs per node.

Table 1:  Alchemical schedule of the softplus perturbation function for leg 1
for the CB8/G3 complex.

$\lambda$ | $\lambda_{1}$ | $\lambda_{2}$ |  $\alpha$a |  $u_{0}$b |  $w_{0}$b
---|---|---|---|---|---
0.000 | 0.000 | 0.000 | 0.250 | 230.0 | 0.000
0.022 | 0.000 | 0.050 | 0.250 | 230.0 | 0.000
0.043 | 0.000 | 0.100 | 0.250 | 220.0 | 0.000
0.065 | 0.000 | 0.150 | 0.250 | 215.0 | 0.000
0.087 | 0.000 | 0.200 | 0.250 | 210.0 | 0.000
0.109 | 0.000 | 0.250 | 0.250 | 205.0 | 0.000
0.130 | 0.000 | 0.300 | 0.250 | 200.0 | 0.000
0.152 | 0.000 | 0.350 | 0.250 | 195.0 | 0.000
0.174 | 0.000 | 0.400 | 0.250 | 188.0 | 0.000
0.196 | 0.000 | 0.400 | 0.250 | 178.0 | 0.000
0.217 | 0.000 | 0.400 | 0.200 | 170.0 | 0.000
0.239 | 0.000 | 0.400 | 0.150 | 160.0 | 0.000
0.261 | 0.000 | 0.400 | 0.150 | 152.0 | 0.000
0.283 | 0.000 | 0.425 | 0.140 | 145.0 | 0.000
0.304 | 0.000 | 0.450 | 0.130 | 140.0 | 0.000
0.326 | 0.000 | 0.475 | 0.140 | 135.0 | 0.000
0.348 | 0.000 | 0.500 | 0.150 | 128.0 | 0.000
0.369 | 0.000 | 0.500 | 0.200 | 120.0 | 0.000
0.391 | 0.000 | 0.500 | 0.200 | 112.0 | 0.000
0.413 | 0.100 | 0.500 | 0.200 | 110.0 | 0.000
0.435 | 0.200 | 0.500 | 0.150 | 110.0 | 0.000
0.457 | 0.300 | 0.500 | 0.150 | 110.0 | 0.000
0.478 | 0.400 | 0.500 | 0.100 | 110.0 | 0.000
0.500 | 0.500 | 0.500 | 0.100 | 110.0 | 0.000

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

Table 2:  Alchemical schedule of the softplus perturbation function for leg 2
for the CB8/G3 complex.

$\lambda$ | $\lambda_{1}$ | $\lambda_{2}$ |  $\alpha$a |  $u_{0}$b |  $w_{0}$b
---|---|---|---|---|---
0.000 | 0.000 | 0.000 | 0.250 | 175 | 0.000
0.022 | 0.000 | 0.084 | 0.250 | 172.5 | 0.000
0.043 | 0.000 | 0.167 | 0.250 | 170 | 0.000
0.065 | 0.000 | 0.233 | 0.225 | 165 | 0.000
0.087 | 0.000 | 0.333 | 0.200 | 160 | 0.000
0.109 | 0.000 | 0.417 | 0.175 | 155 | 0.000
0.130 | 0.000 | 0.500 | 0.150 | 150 | 0.000
0.152 | 0.000 | 0.550 | 0.125 | 145 | 0.000
0.174 | 0.000 | 0.600 | 0.100 | 140 | 0.000
0.196 | 0.000 | 0.650 | 0.095 | 135 | 0.000
0.217 | 0.000 | 0.700 | 0.090 | 130 | 0.000
0.239 | 0.000 | 0.700 | 0.085 | 125 | 0.000
0.261 | 0.000 | 0.700 | 0.080 | 120 | 0.000
0.283 | 0.000 | 0.700 | 0.075 | 112 | 0.000
0.304 | 0.050 | 0.700 | 0.070 | 110 | 0.000
0.326 | 0.100 | 0.700 | 0.0675 | 102.5 | 0.000
0.348 | 0.150 | 0.700 | 0.065 | 105 | 0.000
0.369 | 0.200 | 0.700 | 0.060 | 100 | 0.000
0.391 | 0.250 | 0.650 | 0.055 | 100 | 0.000
0.413 | 0.300 | 0.625 | 0.050 | 100 | 0.000
0.435 | 0.350 | 0.600 | 0.045 | 100 | 0.000
0.457 | 0.400 | 0.575 | 0.040 | 100 | 0.000
0.478 | 0.450 | 0.550 | 0.035 | 100 | 0.000
0.500 | 0.500 | 0.500 | 0.030 | 100 | 0.000

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

##  3 Results

The standard binding free energy estimates, $\Delta G_{b}^{\circ}$, for the
host-guest systems obtained using the ATM method are listed in Table 3
compared to the corresponding experimental measurements, $\Delta
G^{\circ}_{b}$(exp), and the reference computational estimates obtained
through the attach-pull-release (APR) methodology.26 Table 3 also lists the
number of energy and forces evaluations per replicate (a proxy for the
computational cost) for the reference APR calculations, $n_{\rm eval}$(ref),
and for the ATM calculations reported here ($n_{\rm eval}$). The APR method
was selected as representative for this comparison because APR, similar to ATM
and unlike DDM, displaces the ligand in a position in the solvent bulk at a
finite distance from the host.

The ATM binding free energy estimates reported in Table 3 are obtained as the
average of five replicates started from different initial conformations as
reported in Tables 4, 5, and 6. These tables also report the calculated free
energies, $\Delta G_{1}$ and $\Delta G_{2}$, of the two alchemical legs for
each replicate. The binding free energy of each replicate is the difference
between those of the two legs (the excess component) plus the standard state
term $\Delta G^{\circ}_{\rm site}$, which in this case measures out to be
approximately $0.87$ (see Computational Details). The statistical
uncertainties of the averages of each term reported in Tables 4, 5, and 6 are
expressed as 95% confidence interval of the mean based on the t-test
distribution with four degrees of freedom.26 The statistical fluctuations of
the binding free energies among the five conformations for each system were
consistently smaller than the those of each of the legs, suggesting some level
of systematic error cancellation.

The ATM results obtained for the OA-G3 and OA-G6 complexes are in good
agreement (within $0.5$ kcal/mol) with the reference values and well within
the range of estimates obtained with other methods.26 The ATM binding free
energy estimate for the more challenging CB-G3 complex deviates more
substantial from the APR reference ($2$ kcal/mol less favorable) and from
those of the other methods applied to this benchmark system.26 The origin of
this discrepancy is not obvious. However all ATM estimates appear to generally
underestimate binding affinities relative to the other methods bringing them,
perhaps coincidentally, closer to the experimental measurements. The
deviations between the experimental standard binding free energies and the ATM
estimates are $0.71$ kcal/mol, $1.35$ kcal/mol, and $2.08$ kcal/mol for,
respectively, the OA-G3, OA-G6, and CB8-G3 complexes, compared to $1.12$
kcal/mol, $1.83$ kcal/mol, and $4.05$ kcal/mol with APR.

The range of the spread between ATM replicates obtained here for the octaacid
systems is generally larger than with APR and other methods26 albeit at a
generally higher computational cost. For CB8-G3, ATM yields a spread similar
to APR and the other methods with significantly less computational cost ($480$
vs. $2,135$ million energy evaluations as compared to APR, Table 3).

In overall, ATM yields standard binding free energy estimates within the
general range displayed by the established methods tested on the SAMPL6
SAMPLing benchmark set at a similar computational expense.26 These initial
results confirm the validity of the ATM approach.

Table 3:  Standard binding free energy estimates and corresponding
computational effort for the three host-guest complexes with the Alchemical
Transfer Method compared to experimental and reference computed values.

Complex |  $\Delta G^{\circ}_{b}$(exp)a,b |  $\Delta G^{\circ}_{b}$(ref)a,c,d |  $n_{\rm eval}$(ref)c,e |  $\Delta G^{\circ}_{b}$a,e |  $n_{\rm eval}$f
---|---|---|---|---|---
OA-G3 | $-5.18\pm 0.02$ | $-6.3\pm 0.1$ | $458\times 10^{6}$ | $-5.89\pm 0.33$ | $220\times 10^{6}$
OA-G6 | $-4.97\pm 0.02$ | $-6.8\pm 0.1$ | $305\times 10^{6}$ | $-6.32\pm 0.21$ | $220\times 10^{6}$
CB8-G3 | $-6.45\pm 0.06$ | $-10.5\pm 0.6$ | $2135\times 10^{6}$ | $-8.53\pm 0.64$ | $480\times 10^{6}$

aIn kcal/mol. bFrom references 30 and 31. cFrom reference 26. dAPR method.
eThis work, from Tables 4, 5, and 6. fThis work.

Table 4:  Free energy estimates for the two legs of the Alchemical Transfer
Method and corresponding standard binding free energy estimates for the OA-G3
complex starting with each of the the five initial SAMPL6 SAMPLing
conformations.

Conformation |  $\Delta G_{1}$a |  $\Delta G_{2}$a |  $\Delta G^{\circ}_{\rm site}$a |  $\Delta G^{\circ}_{b}$a
---|---|---|---|---
OA-G3-0 | $57.00$ | $50.47$ | $0.87$ | $-5.66$
OA-G3-1 | $57.79$ | $50.97$ | $0.87$ | $-5.95$
OA-G3-2 | $57.74$ | $50.79$ | $0.87$ | $-6.08$
OA-G3-3 | $57.65$ | $51.21$ | $0.87$ | $-5.57$
OA-G3-4 | $57.25$ | $50.19$ | $0.87$ | $-6.18$
Averageb | $57.49\pm 0.42$ | $50.73\pm 0.50$ | $0.87$ | $-5.89\pm 0.33$

aIn kcal/mol. bWith t-test 95% confidence intervals with 4 degrees of freedom
based on the standard deviation of the mean.

Table 5:  Free energy estimates for the two legs of the Alchemical Transfer
Method and corresponding standard binding free energy estimates for the OA-G6
complex starting with each of the the five initial SAMPL6 SAMPLing
conformations.

Conformation |  $\Delta G_{1}$a |  $\Delta G_{2}$a |  $\Delta G^{\circ}_{\rm site}$a |  $\Delta G^{\circ}_{b}$a
---|---|---|---|---
OA-G6-0 | $57.80$ | $50.64$ | $0.87$ | $-6.29$
OA-G6-1 | $58.19$ | $50.78$ | $0.87$ | $-6.54$
OA-G6-2 | $58.27$ | $50.97$ | $0.87$ | $-6.43$
OA-G6-3 | $58.74$ | $51.70$ | $0.87$ | $-6.17$
OA-G6-4 | $58.27$ | $51.25$ | $0.87$ | $-6.14$
Averageb | $58.25\pm 0.41$ | $51.07\pm 0.52$ | $0.87$ | $-6.32\pm 0.21$

aIn kcal/mol. bWith t-test 95% confidence intervals with 4 degrees of freedom
based on the standard deviation of the mean.

Table 6:  Free energy estimates for the two legs of the Alchemical Transfer
Method and corresponding standard binding free energy estimates for the CB8-G3
complex starting with each of the the five initial SAMPL6 SAMPLing
conformations.

Conformation |  $\Delta G_{1}$a |  $\Delta G_{2}$a |  $\Delta G^{\circ}_{\rm site}$a |  $\Delta G^{\circ}_{b}$a
---|---|---|---|---
CB8-G3-0 | $70.50$ | $61.06$ | $0.87$ | $-8.57$
CB8-G3-1 | $70.65$ | $60.81$ | $0.87$ | $-8.97$
CB8-G3-2 | $71.27$ | $61.40$ | $0.87$ | $-9.00$
CB8-G3-3 | $71.74$ | $62.96$ | $0.87$ | $-7.91$
CB8-G3-4 | $69.09$ | $60.02$ | $0.87$ | $-8.20$
Average | $70.65\pm 1.24$ | $61.25\pm 1.34$ | $0.87$ | $-8.53\pm 0.64$

aIn kcal/mol. bWith t-test 95% confidence intervals with 4 degrees of freedom
based on the standard deviation of the mean.

##  4 Discussion

The Alchemical Transfer Method (ATM) presented here implements a perturbation
potential based on rigidly displacing the coordinates of the ligand atoms from
a region in the solvent bulk to the receptor binding site or viceversa. Like
in smart-darting Monte Carlo32, this is accomplished using a displacement
vector that can be thought as connecting a unique pair of points of two
conformational macrostates. At each MD time-step the ligand disappears from
one place and appears in another in a way that is physically not achievable or
”alchemical”. The change in potential energy of the system due to the ligand’s
displacement is the perturbation energy of the $\lambda$-dependent alchemical
potential energy function that is used for conformational sampling and free
energy estimation. While not presented here, the method is applicable to the
calculation of the relative free energy of binding between two ligands13 by
swapping their positions in the bulk and in the receptor site. This work is
ongoing and will be reported in a forthcoming publication.

Similarly to physical pathway methods,19, 33, 21 ATM is relatively easy to
implement in molecular simulation packages because it does not require
modifications of system parameters nor customized single- and dual-topologies
setups that characterize conventional alchemical binding free energy methods.
The method is illustrated here using a plugin implementation on top of the
core OpenMM library.24 The method is agnostic of the underlying energy
function. It has been validated here with explicit solvation and Particle Mesh
Ewald (PME) long range electrostatics. It is conceivably applicable without
approximations to any kind of many-body potential function, including
polarizable,15, 34 quantum-mechanical,35, 36, 37 and artificial neural
network38, 39 potentials.

Unlike alchemical approaches such as double-decoupling,6 which requires two
systems, and dual-system single box alchemical methods,40, 41, 42 which
require dual topologies, ATM works with a single standard model of the
receptor-ligand complex solvated in a solvent box as in conventional molecular
dynamics applications. In addition, ATM does not require soft-core pair
potentials nor the splitting of the alchemical transformations into separate
electrostatic and steric/dispersion steps.23, 25

Similar to single-box alchemical approaches40, 41, 42 ATM avoids alchemical
transformations that place the ligand in vacuum,6, 43 which are particularly
problematic for large and charged ligands.21 The desolvation step of double-
decoupling, for example, removes all of the ligand-water interactions, even
though those of the solvent-exposed region of the ligand are likely to form
again in the subsequent coupling step with the receptor. With ATM, instead,
existing hydration interactions in the bulk are more likely to be replaced by
similar interactions with the ligand displaced into the binding site.
Moreover, unless the ligand is properly restrained, the vacuum intermediate is
likely to introduce hard to converge free energy terms related to the
reorganization of the ligand conformational ensemble from vacuum to the
solvated environment.

ATM has some drawbacks, some of which are technical in nature and are likely
to be addressed in the future. Because it calculates the system energy and
forces twice for each MD time-step,4 once with the ligand in the bulk and
again with the ligand in the receptor pocket, the method is a factor of two
slower per step than standard molecular dynamics. The two energy evaluations
are however independent and can be conceivably run in parallel on two attached
computational devices for added performance. Currently the method also
requires the recalculation of the non-bonded neighbor list after each ligand
displacement resulting in an additional 10 to 15% slow-down per step with
OpenMM for these systems. As we observed here for the CB8-G3 system, the
binding of bulky ligands requires optimized softplus alchemical perturbation
functions trained on trial calculations with the linear alchemical
potential.28, 23, 25 In future work, we plan to implement an adaptive
algorithm to refine the parameters of the alchemical potential function on the
fly.

Here we have validated ATM on the rigorous SAMPL6 SAMPLing dataset.26 The
dataset includes well-studied systems prepared with a single set of force
field parameters and in different initial conformations to probe both
systematic and statistical errors. The binding free energies of the systems
have been computed and validated with a diverse collection of approaches,
including alchemical and physical pathway methods.26 The ATM results for the
octacid systems obtained here are well within the range of estimates reported
in reference 26, thereby adding confidence that the ATM approach is sound and
that it has been implemented correctly. We observed in particular good
agreement with the Attach Pull Release (APR) method,44 a physical pathway
approach in which the guest is progressively displaced into the solvent bulk
to a comparable distance from the host as in this work. ATM yields a
statistically significant less favorable binding free energy estimate than the
other methods for the more challenging CB8-G3 system. The source of the
deviation is unclear, however, ATM appears to generally yield binding free
energies of smaller magnitude and closer to the experimental measurements than
the other methods. Taking into account the relative computational expense, the
statistical uncertainties obtained here indicate that ATM estimates have a
comparable level of reproducibility and computational efficiency as the
methods tested in reference 26.

##  5 Conclusions

We have presented the Alchemical Transfer Method (ATM) for the calculation of
standard binding free energies of non-covalent molecular complexes. The method
is based on a coordinate displacement perturbation of the ligand between the
receptor binding site and the bulk solvent and a thermodynamic cycle connected
by a symmetric intermediate in which the ligand interacts with the receptor
and solvent environments equally. While the approach is alchemical, ATM’s
implementation is as straightforward as physical pathway methods of binding.
ATM does not require splitting the alchemical transformations into
electrostatic and non-electrostatic steps and it does not employ soft-core
pair potentials. We have implemented ATM as a freely available and open-source
plugin of the OpenMM molecular dynamics library. The method and its
implementation have been validated on the SAMPL6 SAMPLing host-guest benchmark
set.

##  6 Acknowledgments

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

## References

  * Simonson 2016 Simonson, T. The physical basis of ligand binding. _In Silico Drug Discovery and Design_ 2016, 3–43
  * Jorgensen 2009 Jorgensen, W. L. Efficient drug lead discovery and optimization. _Acc Chem Res_ 2009, _42_ , 724–733
  * Gallicchio and Levy 2011 Gallicchio, E.; Levy, R. M. Recent Theoretical and Computational Advances for Modeling Protein-Ligand Binding Affinities. _Adv. Prot. Chem. Struct. Biol._ 2011, _85_ , 27–80
  * Gallicchio et al. 2010 Gallicchio, E.; Lapelosa, M.; Levy, R. M. Binding Energy Distribution Analysis Method (BEDAM) for Estimation of Protein-Ligand Binding Affinities. _J. Chem. Theory Comput._ 2010, _6_ , 2961–2977
  * Pal et al. 2019 Pal, R. K.; Ramsey, S.; Gadhiya, S.; Cordone, P.; Wickstrom, L.; Harding, W. W.; Kurtzman, T.; Gallicchio, E. Inclusion of Enclosed Hydration Effects in the Binding Free Energy Estimation of Dopamine D3 Receptor Complexes. _PLoS One_ 2019, _14_ , e0222092
  * Gilson et al. 1997 Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The Statistical-Thermodynamic Basis for Computation of Binding Affinities: A Critical Review. _Biophys. J._ 1997, _72_ , 1047–1069
  * Deng et al. 2018 Deng, N.; Cui, D.; Zhang, B. W.; Xia, J.; Cruz, J.; Levy, R. Comparing alchemical and physical pathway methods for computing the absolute binding free energy of charged ligands. _Phys. Chem. Chem. Phys._ 2018, _20_ , 17081–17092
  * Rocklin et al. 2013 Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects. _The Journal of chemical physics_ 2013, _139_ , 11B606_1
  * Pan et al. 2017 Pan, A. C.; Xu, H.; Palpant, T.; Shaw, D. E. Quantitative characterization of the binding and unbinding of millimolar drug fragments with molecular dynamics simulations. _J. Chem. Theory Comput._ 2017, _13_ , 3372–3377
  * Ohlknecht et al. 2020 Ohlknecht, C.; Perthold, J. W.; Lier, B.; Oostenbrink, C. Charge-Changing Perturbations and Path Sampling via Classical Molecular Dynamic Simulations of Simple Guest–Host Systems. _J. Chem. Theory Comput._ 2020, _Article ASAP_
  * Lee et al. 2020 Lee, T.-S.; Lin, Z.; Allen, B. K.; Lin, C.; Radak, B. K.; Tao, Y.; Tsai, H.-C.; Sherman, W.; York, D. M. Improved Alchemical Free Energy Calculations with Optimized Smoothstep Softcore Potentials. _J. Chem. Theory Comput._ 2020,
  * Steinbrecher et al. 2007 Steinbrecher, T.; Mobley, D. L.; Case, D. A. Nonlinear scaling schemes for Lennard-Jones interactions in free energy calculations. _J Chem Phys_ 2007, _127_ , 214108
  * Mey et al. 2020 Mey, A. S. J. S.; Allen, B. K.; Macdonald, H. E. B.; Chodera, J. D.; Hahn, D. F.; Kuhn, M.; Michel, J.; Mobley, D. L.; Naden, L. N.; Prasad, S.; Rizzi, A.; Scheen, J.; Shirts, M. R.; Tresadern, G.; Xu, H. Best Practices for Alchemical Free Energy Calculations [Article v1.0]. _Living Journal of Computational Molecular Science_ 2020, _2_ , 18378
  * Lee et al. 2020 Lee, T.-S.; Allen, B. K.; Giese, T. J.; Guo, Z.; Li, P.; Lin, C.; McGee Jr, T. D.; Pearlman, D. A.; Radak, B. K.; Tao, Y.; Tsai, H.-C.; Xu, H.; Sherman, W.; York, D. M. Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery. _J. Chem. Inf. Model._ 2020,
  * Harger et al. 2017 Harger, M.; Li, D.; Wang, Z.; Dalby, K.; Lagardère, L.; Piquemal, J.-P.; Ponder, J.; Ren, P. Tinker-OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs. _J. Comp. Chem._ 2017, _38_ , 2047–2055
  * Darden et al. 1993 Darden, T. A.; York, D. M.; Pedersen, L. G. _J. Chem. Phys._ 1993, _98_ , 10089–10092
  * Gallicchio et al. 2009 Gallicchio, E.; Paris, K.; Levy, R. M. The AGBNP2 Implicit Solvation Model. _J. Chem. Theory Comput._ 2009, _5_ , 2544–2564
  * Zou et al. 2019 Zou, J.; Tian, C.; Simmerling, C. Blinded prediction of protein–ligand binding affinity using Amber thermodynamic integration for the 2018 D3R grand challenge 4. _J. Comp. Aid. Mol. Des._ 2019, _33_ , 1021–1029
  * Woo and Roux 2005 Woo, H.-J.; Roux, B. Calculation of absolute protein-ligand binding free energy from computer simulations. _Proc. Natl. Acad. Sci. USA_ 2005, _102_ , 6825–6830
  * Limongelli et al. 2013 Limongelli, V.; Bonomi, M.; Parrinello, M. Funnel metadynamics as accurate binding free-energy method. _Proc. Natl. Acad. Sci._ 2013, _110_ , 6358–6363
  * Cruz et al. 2020 Cruz, J.; Wickstrom, L.; Yang, D.; Gallicchio, E.; Deng, N. Combining Alchemical Transformation with a Physical Pathway to Accelerate Absolute Binding Free Energy Calculations of Charged Ligands to Enclosed Binding Sites. _J. Chem. Theory Comput._ 2020, _16_ , 2803–2813
  * Kilburg and Gallicchio 2018 Kilburg, D.; Gallicchio, E. Assessment of a Single Decoupling Alchemical Approach for the Calculation of the Absolute Binding Free Energies of Protein-Peptide Complexes. _Frontiers in Molecular Biosciences_ 2018, _5_ , 22
  * Pal and Gallicchio 2019 Pal, R. K.; Gallicchio, E. Perturbation potentials to overcome order/disorder transitions in alchemical binding free energy calculations. _J. Chem. Phys._ 2019, _151_ , 124116
  * Eastman et al. 2017 others,, et al. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. _PLoS Comp. Bio._ 2017, _13_ , e1005659
  * Khuttan et al. 2021 Khuttan, S.; Azimi, S.; Wu, J. Z.; Gallicchio, E. Alchemical Transformations for Concerted Hydration Free Energy Estimation with Explicit Solvation. _J. Chem. Phys_ 2021, In press.
  * Rizzi et al. 2020 Rizzi, A. et al. The SAMPL6 SAMPLing challenge: Assessing the reliability and efficiency of binding free energy calculations. _J. Comp. Aid. Mol. Des._ 2020, 1–33
  * Gallicchio et al. 2015 Gallicchio, E.; Xia, J.; Flynn, W. F.; Zhang, B.; Samlalsingh, S.; Mentes, A.; Levy, R. M. Asynchronous replica exchange software for grid and heterogeneous computing. _Computer Physics Communications_ 2015, _196_ , 236–246
  * Kilburg and Gallicchio 2018 Kilburg, D.; Gallicchio, E. Analytical Model of the Free Energy of Alchemical Molecular Binding. _J. Chem. Theory Comput._ 2018, _14_ , 6183–6196
  * Tan et al. 2012 Tan, Z.; Gallicchio, E.; Lapelosa, M.; Levy, R. M. Theory of binless multi-state free energy estimation with applications to protein-ligand binding. _J. Chem. Phys._ 2012, _136_ , 144102
  * Murkli et al. 2019 Murkli, S.; McNeill, J. N.; Isaacs, L. Cucurbit[8]uril guest complexes: blinded dataset for the SAMPL6 challenge. _Supramolecular Chemistry_ 2019, _31_ , 150–158
  * Sullivan et al. 2019 Sullivan, M. R.; Yao, W.; Gibb, B. C. The thermodynamics of guest complexation to octa-acid and tetra-endo-methyl octa-acid: reference data for the sixth statistical assessment of modeling of proteins and ligands (SAMPL6). _Supramolecular chemistry_ 2019, _31_ , 184–189
  * Andricioaei et al. 2001 Andricioaei, I.; Straub, J. E.; Voter, A. F. Smart darting monte carlo. _J. Chem. Phys._ 2001, _114_ , 6994–7000
  * Henriksen et al. 2015 Henriksen, N. M.; Fenley, A. T.; Gilson, M. K. Computational calorimetry: high-precision calculation of host–guest binding thermodynamics. _J. Chem. Theory Comput._ 2015, _11_ , 4377–4394
  * Panel et al. 2018 Panel, N.; Villa, F.; Fuentes, E. J.; Simonson, T. Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations. _Biophys. J._ 2018, _114_ , 1091–1102
  * Beierlein et al. 2011 Beierlein, F. R.; Michel, J.; Essex, J. W. A simple QM/MM approach for capturing polarization effects in protein- ligand binding free energy calculations. _J. Phys. Chem. B_ 2011, _115_ , 4911–4926
  * Lodola and De Vivo 2012 Lodola, A.; De Vivo, M. _Adv. Prot. Chem. Struct. Biol._ ; Elsevier, 2012; Vol. 87; pp 337–362
  * Hudson et al. 2019 Hudson, P. S.; Woodcock, H. L.; Boresch, S. Use of interaction energies in QM/MM free energy simulations. _J. Chem. Theory Comput._ 2019, _15_ , 4632–4645
  * Smith et al. 2019 Smith, J. S.; Nebgen, B. T.; Zubatyuk, R.; Lubbers, N.; Devereux, C.; Barros, K.; Tretiak, S.; Isayev, O.; Roitberg, A. E. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. _Nature Comm._ 2019, _10_ , 1–8
  * Rufa et al. 2020 Rufa, D. A.; Macdonald, H. E. B.; Fass, J.; Wieder, M.; Grinaway, P. B.; Roitberg, A. E.; Isayev, O.; Chodera, J. D. Towards chemical accuracy for alchemical free energy calculations with hybrid physics-based machine learning/molecular mechanics potentials. _BioRxiv_ 2020,
  * Gapsys et al. 2015 Gapsys, V.; Michielssens, S.; Peters, J. H.; de Groot, B. L.; Leonov, H. _Molecular Modeling of Proteins_ ; Springer, 2015; pp 173–209
  * Ekimoto et al. 2018 Ekimoto, T.; Yamane, T.; Ikeguchi, M. Elimination of finite-size effects on binding free energies via the warp-drive method. _Journal of chemical theory and computation_ 2018, _14_ , 6544–6559
  * Macchiagodena et al. 2020 Macchiagodena, M.; Pagliai, M.; Karrenbrock, M.; Guarnieri, G.; Iannone, F.; Procacci, P. Virtual Double-System Single-Box: A Nonequilibrium Alchemical Technique for Absolute Binding Free Energy Calculations: Application to Ligands of the SARS-CoV-2 Main Protease. _J. Chem. Theory Comput._ 2020,
  * Mobley and Gilson 2017 Mobley, D. L.; Gilson, M. K. Predicting binding free energies: frontiers and benchmarks. _Ann. Rev. Bioph._ 2017, _46_ , 531–558
  * Velez-Vega and Gilson 2013 Velez-Vega, C.; Gilson, M. K. Overcoming dissipation in the calculation of standard binding free energies by ligand extraction. _J. Comp. Chem._ 2013, _34_ , 2360–2371