∎

11institutetext:  S. Azimi 22institutetext: Department of Chemistry, Brooklyn
College of the City University of New York
PhD Program in Biochemistry, Graduate Center of the City University of New
York 33institutetext: J. Z. Wu 44institutetext: Department of Chemistry,
Brooklyn College of the City University of New York
PhD Program in Chemistry, Graduate Center of the City University of New York
55institutetext: S. Khuttan 66institutetext: Department of Chemistry, Brooklyn
College of the City University of New York
PhD Program in Biochemistry, Graduate Center of the City University of New
York 77institutetext: T. Kurtzman 88institutetext: Department of Chemistry,
Lehman College of the City University of New York
PhD Program in Chemistry, Graduate Center of the City University of New York
PhD Program in Biochemistry, Graduate Center of the City University of New
York 99institutetext: N. Deng 1010institutetext: Department of Chemistry and
Physical Sciences, Pace University, New York, New York 1111institutetext: E.
Gallicchio 1212institutetext: Department of Chemistry, Brooklyn College of the
City University of New York
PhD Program in Chemistry, Graduate Center of the City University of New York
PhD Program in Biochemistry, Graduate Center of the City University of New
York
1212email<EMAIL_ADDRESS>
# Application of the Alchemical Transfer and Potential of Mean Force Methods
to the SAMPL8 Host-Guest Blinded Challenge

Solmaz Azimi  Joe Z. Wu  Sheenam Khuttan  Tom Kurtzman  Nanjie Deng  Emilio
Gallicchio

###### Abstract

We report the results of our participation in the SAMPL8 GDCC Blind Challenge
for host-guest binding affinity predictions. Absolute binding affinity
prediction is of central importance to the biophysics of molecular association
and pharmaceutical discovery. The blinded SAMPL series have provided an
important forum for assessing the reliability of binding free energy methods
in an objective way. In this blinded challenge, we employed two binding free
energy methods, the newly developed alchemical transfer method (ATM) and the
well established potential of mean force (PMF) physical pathway method, using
the same setup and force field model. The calculated binding free energies
from the two methods are in excellent quantitative agreement. Importantly, the
results from the two methods were also found to agree well with the
experimental binding affinities released subsequently, with an $R^{2}$ of 0.89
(ATM) and 0.83 (PMF). Given that the two free energy methods are based on
entirely different thermodynamic pathways, the close agreement between the
results from the two methods and their general agreement with the experimental
binding free energies are a testament to the the high quality achieved by
theory and methods. The study provides further validation of the novel ATM
binding free energy estimation protocol and it paves the way to to further
extensions of the method to more complex systems.

††journal: jcamd

##  1 Introduction

The Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) series
of community challengesgeballe2010sampl2 ; mobley2014blind ; amezcua2021sampl7
have been organized to validate computational methods of molecular solvation
and binding in an unbiased way. SAMPL participants are asked to quantitatively
predict experimental measurements that are publicly disclosed only after the
predictions are submitted. The format of the challenges allows the robust
assessment of computational methods and have significantly contributed to
their advancement.mobley2017predicting  As computational models of small
molecule binding to protein receptors increasingly emerge as important
elements of structure-based drug discovery,Jorgensen2009 ; armacost2020novel
it is critical that the reliability of these models is independently assessed
and validated. We have contributed to several editions of the SAMPL challenges
to validate the ability of our computational models to accurately predict
host-guest and protein-ligand binding affinities.Gallicchio2012a ;
Gallicchio2014octacid ; GallicchioSAMPL4 ; deng2016large ; pal2016SAMPL5 .

In this work, we apply two conceptually orthogonal yet equivalent binding free
energy estimation methods, the Alchemical Transfer Method
(ATM)wu2021alchemical  and the Potential of Mean Force (PMF)deng2018comparing
method, to the SAMPL8 GDCC challenge
set111github.com/samplchallenges/SAMPL8/tree/master/host_guest/GDCC . The
modeled predictions are tested against each other, as well as with the blinded
experimental binding free energies measured by the Gibb
Group.suating2020proximal 222
github.com/samplchallenges/SAMPL8/blob/master/host_guest/Analysis/ExperimentalMeasurements/Final-
Data-Table-031621-SAMPL8.docx

In principle, computational models should yield equivalent binding free energy
predictions as long as they are based on the same chemical model and physical
description of inter-atomic interactions. By ensuring consistency between two
independent computational estimates, we can achieve an increased level of
confidence in the theoretical accuracy of the models and in the correctness of
their implementation. Furthermore, by comparing the computational predictions
to the experimental measurements in a blinded, unbiased fashion, we can assess
the predictive capability that can be expected of the models in actual
chemical applications.

While a variety of empirical methods are commonly used to model the binding
affinities of molecular complexes,sledz2018protein ; seidel2020applications
here we are concerned with methods based on physical models of inter-atomic
interactions and a rigorous statistical mechanics theory of the free energy of
molecular binding.Gilson:Given:Bush:McCammon:97 ; Gallicchio2011adv ;
cournia2020rigorous  Binding free energy methods are classified as physical or
alchemical depending on the nature of the thermodynamic path employed to
connect the unbound to the bound states of the molecular complex for computing
the reversible work of binding.Gallicchio2021binding  Physical pathway methods
define a physical path in coordinate space in 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 artificial
intermediate states in which the ligand is progressively decoupled from the
solution environment and coupled to the receptor.

In this work, we compare the results of the PMF method,deng2018comparing  a
physical pathway method, to that of the ATM alchemical methodwu2021alchemical
on identically prepared molecular systems. Because free energy is a
thermodynamic state function, binding free energy estimates should be
independent of the specific path employed, whether physical or alchemical.
Obtaining statistically equivalent estimates of the binding free energies
using these two very different thermodynamic paths constitutes a robust
validation of both methods. The very recently developed ATM, in particular,
benefits from the backing of the more established PMF method in this
application.

This paper is organized as follows. We first review the PMF and ATM methods,
describe the host-guest systems included in the SAMPL8 GDCC challenge, and
provide the system setup and simulation details of our free energy
calculations. We then present the binding free energy estimates we obtained
with the PMF and ATM approaches and compare them to each other and with the
experimental measurements that were disclosed only after the predictions were
submitted to the SAMPL8 organizers. Overall, the work shows that the ATM and
PMF methods provide consistent binding free energy estimates that, in
conjunction with the force field model employed here, are in statistical
agreement with experimental observations.

##  2 Theory and Methods

###  2.1 The Potential of Mean Force Method

The Potential of Mean Force method, hereon PMF, employed in this work is a
physical binding pathway approach fully described in reference 13. Here, we
briefly summarize the statistical mechanics basis of the method.
Implementation details specific to this work are described in the
Computational Details section.

The PMF method estimates the standard free energy of binding as the sum of the
free energy changes of the following processes:

  1. 1.

The transfer of one ligand molecule from an ideal solution at the standard
concentration $C^{\circ}=1M$ to a region in the solvent bulk of volume equal
to the volume of the receptor binding site, followed by the imposition of
harmonic restraints that keep the ligand in a chosen reference binding
orientation. The free energy term corresponding to this process, denoted as
$\Delta G^{\rm bulk}_{\rm restr}$, is evaluated analytically.

  2. 2.

The transfer of the ligand molecule from the solvent bulk to the receptor
binding site along a suitable physical pathway (see Computational Details).
The free energy change along this pathway is described by a potential of mean
force parameterized by the distance between two reference atoms of the ligand
and the receptor (Figure 1). The free energy change for this process, denoted
by $w(r_{\rm min})-w(r^{\ast})$, is given by the value at the minimum of the
potential of mean force relative to the value in the bulk.

  3. 3.

$\Delta G_{\rm vibr}$ is related to the ratio of the configurational partition
functions of the ligand within the binding site of the receptor vs. when it is
harmonically restrained at the bulk location $r^{\ast}$.

  4. 4.

The release of the harmonic restraints while the ligand is bound to the
receptor. The free energy change for this process, denoted by $-\Delta G_{\rm
restr}^{\rm bound}$, is evaluated by Bennett’s Acceptance Ratio method (BAR).

Hence, the PMF estimate of the free energy of binding is given by

$\Delta G^{\circ}_{b}=\Delta G^{\rm bulk}_{\rm restr}+[w(r_{\rm
min})-w(r^{\ast})]+\Delta G_{\rm vibr}-\Delta G_{\rm restr}^{\rm bound}$ (1)

Additional computational details and parameters used in this work to implement
the PMF calculations are described in the Computational Details section.

Figure 1: Schematic of Potential of Mean Force (PMF) method. From left to
right, the figure represents the physical pathway that the ligand undergoes
from the bound to unbound state. Shown above is a sequence of 3 snapshots
representing 3 of the 20 umbrella windows, where the ligand gets pulled at
varying distances along the physical pathway away from the host (through the
use of reference atoms assigned to both the ligand and host). The red dots
represent the oxygen atoms of water molecules. The big bulky molecule
represents the TEMOA host, while the small molecule represents the G1 guest.

###  2.2 The Alchemical Transfer Method

The Alchemical Transfer Method, hereon ATM, is a recently-developed method to
compute the absolute binding free energy of molecular complexes. The method is
fully described in reference 12. Here, we give only a brief overview of ATM,
particularly focusing on the aspects specific to this work. Further
implementation details are described in the Computational Details section.

Given the standard free energy of binding $\Delta G^{\circ}_{b}$, defined as
the difference in free energy between the bound complex and the unbound
components, $\Delta G^{\circ}_{\rm b}=\Delta G^{\circ}_{\rm site}+\Delta
G^{\ast}_{b}$. ATM computes the excess component of the binding free energy,
$\Delta G^{\ast}_{b}$, defined as the reversible work for transferring the
ligand from a region of volume $V_{\rm site}$ in the solvent bulk to a region
of the same volume in the receptor binding site.Gallicchio2011adv  The
standard free energy of binding is given by the excess component plus the
ideal component, $\Delta G^{\circ}_{\rm site}=-k_{B}T\ln C^{\circ}V_{\rm
site}$, which corresponds to the free energy change of transferring one ligand
molecule from an ideal solution at the standard concentration $C^{\circ}=1M$
to a region in the solvent bulk of volume that is equal to the volume of the
receptor binding site, $V_{\rm site}$.Gilson:Given:Bush:McCammon:97  The
concentration-dependent ideal term is computed analytically and the excess
component is computed by ATM using numerical molecular simulations described
in Computational Details and below.

In ATM, the transfer of the ligand from the solvent bulk to the receptor
binding site is carried out in two alchemical steps that connect the bound and
unbound end states to one alchemical intermediate (Figure 2), in which the
ligand molecule interacts equally with both the receptor and the solvent bulk
at half strength. The potential energy function of the alchemical intermediate
is defined as

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

where $x_{S}$ denotes the coordinates of the atoms of the receptor and of the
solvent, $x_{L}$ denotes the coordinates of the atoms of the ligand while in
the receptor binding site, and $h$ is the constant displacement vector that
brings the atoms of the ligand from the receptor site to the solvent bulk
site. In this scheme, $U(x_{S},x_{L})$ is the potential energy of the system
when the ligand is in the binding site, $U(x_{S},x_{L}+h)$ is the potential
energy after translating the ligand rigidly into the solvent bulk, and
$U_{1/2}(x_{S},x_{L})$ is the hybrid alchemical potential given by the average
of the two. In the alchemical intermediate state, receptor atoms and solvent
molecules interact with the ligand at half strength but at both ligand
locations. Similarly, the force that ligand atoms interact with receptor atoms
and solvent molecules at the intermediate state is an average of the forces
exerted by the ligand at the two distinct locations. As discussed in reference
12, the ATM alchemical intermediate has an analogous role as the vacuum
intermediate state in the conventional double-decoupling
method,Gilson:Given:Bush:McCammon:97  but without fully dehydrating the
ligand.

Figure 2: The Alchemical Transfer Method (ATM) involves two simulation legs,
which, in total, transfer the ligand from the solvent bulk to the binding site
of the receptor. The two legs connect the bound and unbound end states through
an alchemical intermediate that involves the ligand molecule interacting
equally with both the receptor and the solvent bulk at half strength. Here,
the receptor is the TEMOA host and the ligand is the G4 guest. The green box
represents the solvent box with water molecules designated in blue. In the
TEMOA structure, carbon atoms are represented in cyan and oxygen atoms in red.

The bound and unbound states of the complex are connected to the common
intermediate by means of alchemical potentials of the form

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

where $U_{0}(x)$ denotes the potential energy function of the initial state,
which is either $U(x_{S},x_{L})$, corresponding to the bound complex in Leg 1
(Figure 2), or $U(x_{S},x_{L}+h)$, corresponding to Leg 2 (Figure 2),
$\lambda$ is a progress parameter that goes from $0$ to $1/2$,

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

is the binding energy function.Gallicchio2010  In Equation 4, $U_{1}(x)$
denotes the potential energy function of the end state which is either
$U(x_{S},x_{L}+h)$, corresponding to the unbound complex in Leg 1 of Figure 2,
or $U(x_{S},x_{L})$, corresponding to the bound complex in Leg 2 (Figure 2).
Finally,

$u_{\rm sc}(u)=u;\quad u\leq u_{c}$ (5) $u_{\rm sc}(u)=(u_{\rm
max}-u_{c})f_{\rm sc}\left[\frac{u-u_{c}}{u_{\rm
max}-u_{c}}\right]+u_{c};\quad u>u_{c}$ (6)

with

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

and

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

is a soft-core perturbation function that avoids singularities near the
initial states of each leg ($\lambda=0$). The parameters of the soft-core
function, $u_{\rm max}$, $u_{c}$, and $a$ used in this work are listed in
Computational Details.

The free energy change for each leg is obtained by multi-state thermodynamic
reweightingTan2012  using the perturbation energies $u_{\rm sc}[u(x)]$
collected during the molecular dynamics runs at various values of $\lambda$.
As illustrated by the thermodynamic cycle in Figure 2, the excess component of
the binding free energy is obtained by the difference of the free energies of
the two legs:

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

Because the end states of ATM are similar to that of the PMF method summarized
above, the two methods compute the same free energy of binding. However, each
employs a different thermodynamic path. The PMF method progressively displaces
the ligand from the binding site to the bulk along a physical path, whereas
ATM employs an unphysical alchemical path, in which the ligand is displaced
directly from the binding site to the solvent bulk.

###  2.3 SAMPL8 Systems

The chemical structures of the two hosts and 5 guests molecules are shown in
Fig. 3. Both the hosts TEETOA and TEMOA are octaacids that carry a net charge
of -8 at the pH value of 11.5 used in the experiment. The five guests, with
the exception of the protonated G2 (namely G2P), are carboxylate derivatives
that are also negatively charged at the same pH. The computational
calculations employed the initial host and guest structure files provided in
the SAMPL8 dataset found at
https://github.com/samplchallenges/SAMPL8/tree/master/host_guest/GDCC.

Figure 3: Superimposed benchmark systems in this study. The two hosts,
tetramethyl octa acid (TEMOA) and tetraethyl octa acid (TEETOA), are shown in
licorice representation, with light gray corresponding to TEETOA and dark gray
to TEMOA. Both light and dark gray represent carbon atoms and red, oxygen
atoms. The six guests that are bound to the hosts are shown in ball-and-stick
(CPK) representation, for which the color of the structure corresponds to the
label of the guest. G2D designates deprotonated G2 and G2P, protonated G2.
Note that ball-and-stick representation undermines the aromaticity of the six-
membered ring. For the guests, green corresponds to carbon atoms, red oxygen
atoms, and white hydrogen atoms.

###  2.4 Computational Details

The guests were manually docked to each host using Maestro (Schrödinger, Inc.)
to render a set of host-guest molecular complexes that were then used to
derive forcefield parameters with AmberTools. The complexes were assigned
GAFF2/AM1-BCC parameters and solvated in a water box with a 12 Angstrom
solvent buffer and sodium counterions to balance the negative charge. The
position and orientation of the host for each complex were restrained near the
center of the box and along the diagonal with a flat-bottom harmonic potential
of force constant 25.0 kcal/(mol Å2) and a tolerance of 1.5 Å was set on the
heavy atoms at the lower cup of the molecule (the first 40 atoms of the host
as listed in the provided files). The systems were energy minimized and
thermalized at 300 K prior to proceeding with the ATM and PMF calculations.

####  2.4.1 PMF setup

The computation of the standard binding free energies using the PMF method
involves the following steps:deng2018comparing  (1) applying a harmonic
restraint on the three Euler angles of the guest in the bound state to
restrain guest orientation; (2) applying a harmonic restraint on the polar and
azimuthal angles in spherical coordinates to restrain the guest center along a
fixed axis when it binds/unbinds; (3) reversibly extracting the guest from the
binding pocket along the chosen axis until it reaches the bulk region; (4)
release the restraints on the guest center and guest orientation, which allows
the guest to occupy the standard volume and rotate freely in the bulk solvent.
The standard binding free energy is then obtained by summing up the reversible
work associated with each of the above steps using Eq. (1).

The position and orientation of the guest relative to the host was controlled
using coordinate systems which consisted of 3 reference atoms of the host (P1,
P2, and P3) and 3 reference atoms of the guest (L1, L2, and
L3).Boresch:Karplus:2003  For all the hosts, P1 was chosen to be the center of
the bottom ring of each host and L1 the center of each guest molecule which
lies approximately 4 Angstroms away from P1. The PMF was calculated along the
P1-L1 distance using umbrella sampling with biasing potentials having a force
constant of 1000 kJ/(mol nm2). The three Euler angles and two polar and
azimuthal angles were restrained using harmonic potentials with a force
constant of 1,000 kJ/(mol rad2) centered on the angles of the thermalized
structures such that the guest is pulled straight out of the pocket of the
host while minimizing collisions with the sidechains of the rim of the host.
It is important to note that an unobstructed path is necessary for the guest’s
pull axis for the PMF method.

Equilibration (1.2 ns) and production (20 ns) umbrella sampling was then
initiated over 20 umbrella windows to cover a distance of 4.0 to 18.0
Angstroms, i.e. from within the binding region to the bulk along the P1-L1
axis. WHAM analysis was used to generate the PMF and the corresponding
uncertainties by bootstrapping. The free energy of releasing the angular
restraints in the bulk and in the bound state were computed using BAR as
implemented in GROMACS.pronk2013gromacs

####  2.4.2 ATM Setup

Each of the Cartesian components of the translation vector $h$ were set to
approximately half of the longest diagonal of the simulation box to place the
ligand near the corner of the solvent box farthest away from the host and its
periodic images (Fig. 2). Beginning at the bound state at $\lambda=0$, the
systems were then progressively annealed to the symmetric alchemical
intermediate at $\lambda=1/2$ during a $250$ ps run using the ATM alchemical
potential energy function for Leg 1 [Eq. (2)]. This step yields a suitable
initial configuration of the system without severe unfavorable repulsive
interactions at either end states of the alchemical path so that molecular
dynamics replica exchange alchemical simulation can be conducted for each leg
as described below.

In order to prevent large attractive interactions between opposite charges at
small distances in nearly uncoupled states, polar hydrogen atoms with zero
Lennard-Jones parameters were modified to $\sigma_{\rm LJ}=0.1$ Å and
$\epsilon_{\rm LJ}=10^{-4}$ kcal/mol. khuttan2021single  We established that
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. Alchemical MD
calculations were conducted with the OpenMM 7.3eastman2017openmm  MD engine
and the SDM integrator plugin (github.com/Gallicchio-
Lab/openmm_sdm_plugin.git) using the OpenCL platform. In order to maintain the
temperature at 300 K, a Langevin thermostat with a time constant of 2 ps was
implemented. For each ATM leg, Hamiltonian Replica Exchange in $\lambda$ space
was conducted every 5 ps with the ASyncRE software gallicchio2015asynchronous
that is customized for OpenMM and SDM (github.com/Gallicchio-Lab/async_re-
openmm.git). Each leg employed 11 $\lambda$ states uniformily distributed
between $\lambda=0$ and $\lambda=1/2$. All ATM calculations employed the soft-
core perturbation energy with parameters $u_{\rm max}=300$ kcal/mol,
$u_{c}=100$ kcal/mol, and $a=1/16$. A flat-bottom harmonic potential between
the centers of mass of the host and the guest with a force constant of 25
kcal/mol $\AA^{2}$ was applied for a distance greater than $4.5\AA$ to define
the binding site region ($V_{\rm site}$). The concentration-dependent term,
$\Delta G^{\circ}_{\rm site}=-k_{B}T\ln C^{\circ}V_{\rm site}=0.87$, which
corresponds to 300 K temperature and the volume $V_{\rm site}$ of a sphere
with a radius of $4.5\AA$, was added to yield the final free energy estimate.
Perturbation energy samples and trajectory frames were collected every 5 ps.
Every replica was simulated for a minimum of 10 ns. For ATM, UWHAM was used to
compute binding free energies and the corresponding uncertainties with the
first 5 ns of the trajectory discarded.

####  2.4.3 Free Energy of Binding for Ligands in Multiple Protonation States

When multiple chemical species contribute to binding, we use the free energy
combination formulaGallicchio2011adv

$\Delta G_{b}^{\circ}=-kT\ln\sum_{i}P_{0}(i)e^{-\beta\Delta
G_{b}^{\circ}(i)},$ (10)

where $\Delta G_{b}^{\circ}(i)$ is the standard binding free energy for
species $i$ and $P_{0}(i)$ is the population of that species in the unbound
state. In the case of an acid/base equilibrium with acidity constant

$K_{a}=\frac{[A^{-}][H^{+}]}{[HA]}=\frac{[A^{-}]}{[HA]}10^{-pH}=\alpha
10^{-pH},$ (11)

where $[\ldots]$ are concentration in molar units,

$\alpha=10^{pH-pKa},$ (12)

is the concentration ratio of the deprotonated and protonated forms, the
population fraction of the deprotonated species is

$P_{0}(A^{-})=\frac{[A^{-}]}{[HA]+[A^{-}]}=\frac{\alpha}{1+\alpha}$ (13)

and the population fraction of the protonated species is

$P_{0}(HA)=\frac{[HA]}{[HA]+[A^{-}]}=1-P_{0}(A^{-})=\frac{1}{1+\alpha}.$ (14)

The populations of each protonation state of the ligands and the corresponding
standard binding free energies $\Delta G_{b}^{\circ}(A^{-})$ and $\Delta
G_{b}^{\circ}(HA)$ are combined using Eq. (10) to obtain an estimate of the
observed free energy of binding.

This strategy was employed for the guest G2, 4-bromophenol, which exists in
two protonation states. A pH of 11.5, as indicated in the SAMPL8 GitHub site,
and a pKa of 9.17 (pubchem.ncbi.nlm.nih.gov/compound/4-bromophenol) was used
to calculate the concentrations of the protonation states and combine them
with the calculated binding free energies to yield a binding free energy
estimate for G2 (see Table 5).

##  3 Results

The results are presented as follows. Table 1 summarizes the absolute binding
free energy predictions from ATM and PMF submitted to the SAMPL8 organizers,
compared to the experimental values which were disclosed to us only after
submission. The results of the constituent calculations for each method that
led to the binding free energy predictions are listed in Tables 3 and 4 for
the ATM and PMF methods, respectively. These tables report the values of the
free energy changes for each leg of the ATM calculations and the components of
the PMF estimates, including those of the vibrational free energy and the
restraint free energy that contribute to the overall PMF process. The free
energy analysis for the protonated and deprotonated species implicated in the
complexes of the G2 guest is illustrated in Table 5.

###  3.1 Absolute Binding Free Energy Estimates by ATM and PMF

Table 1: PMF and ATM standard binding free energy predictions compared to the
experimental values.

Complex | Experimenta | ATMa | PMFa
---|---|---|---
TEMOA-G1 | $-6.96\pm 0.2$ | $-6.71\pm 0.3$ | $-6.43\pm 0.4$
TEMOA-G2 | $-8.41\pm 0.1$ | $-9.90\pm 0.8$ | $-9.37\pm 0.8$
TEMOA-G3 | $-5.78\pm 0.1$ | $-8.26\pm 0.3$ | $-8.71\pm 0.4$
TEMOA-G4 | $-7.72\pm 0.1$ | $-8.63\pm 0.3$ | $-8.79\pm 0.6$
TEMOA-G5 | $-6.67\pm 0.1$ | $-7.70\pm 0.3$ | $-8.15\pm 0.8$
TEETOA-G1 | $-4.49\pm 0.2$ | $-1.07\pm 0.3$ | $-1.38\pm 0.8$
TEETOA-G2 | $-5.16\pm 0.1$ | $-4.76\pm 0.3$ | $-6.22\pm 1.8$
TEETOA-G3 | NB | $-1.65\pm 0.3$ | $-1.42\pm 0.8$
TEETOA-G4 | $-4.47\pm 0.2$ | $-2.51\pm 0.3$ | $-2.25\pm 0.8$
TEETOA-G5 | $-3.32\pm 0.1$ | $-2.82\pm 0.3$ | $-3.36\pm 1.9$

a In kcal/mol.

Table 2: Agreement metrics (root mean square error, RMSE, correlation
coefficient of determination, $R^{2}$, slope of the linear regression, $m$,
and Kendall rank order correlation coefficient, $\tau$) between the computed
binding free energies and the experimental measurements.

| RMSE | $R^{2}$ | m | $\tau$
---|---|---|---|---
ATM/PMF | 0.60 | 0.99 | 1.05 | 1.00
Exp./ATM | 1.71 | 0.89 | 1.65 | 0.69a
Exp./PMF | 1.79 | 0.83 | 1.50 | 0.69a

a TEETOA-G3, a non-binder experimentally, was included in the $\tau$
calculation as the weakest complex.

The binding free energy estimates obtained from the two complementary
computational methods, ATM and PMF, are in very good agreement with an $R^{2}$
value of 0.965 and an RMSE value of 0.989(?) kcal/mol. In addition, the
ranking of the binding free energies of the complexes between the ATM and PMF
datasets is in perfect agreement. Both methods consistently estimated the
complex with the most favorable binding free energy to be TEMOA-G2, with a
free energy value of -9.90 kcal/mol predicted by ATM and -9.37 kcal/mol by
PMF. The least favorable binding free energy was predicted for the complex
TEETOA-G1 by both methods, -1.07 kcal/mol by ATM and -1.38 kcal/mol by PMF.
Both methods predict that all of the guests bind TEMOA more favorably than
TEETOA.

All of the carboxylic acid guests were modeled as ionic. We modeled both
protonation states of the G2 guest (Tables 3 and 4) and combined the
corresponding binding free energies using the experimental pKa of the guest
(Table 5). With a discrepancy of 2.77 kcal/mol, the deprotonated G2 molecule
(hereon G2D) yielded the most divergent binding free energy estimate between
the ATM and PMF datasets. Nevertheless, since this protonation state is found
to contribute little to binding (Table 5), the observed discrepancy did not
affect significantly the correspondence between the two sets of SAMPL8 binding
free energy predictions.

The molecular dynamics trajectories consistently yielded the expected binding
mode of the guests to the TEMOA and TEETOA hosts. The polar/ionic end of the
guests is oriented towards the water solvent while the more non-polar end of
the molecule is inserted into the binding cavity of the hosts (Figure 3). In
the complexes, the ethyl sidechains of the TEETOA host point outward extending
further the host binding cavity and the surface of contact between the guests
and the hosts. In the apo state, however, the ethyl sidechains are observed to
be mostly folded into the TEETOA cavity (not shown). We hypothesize that the
conformational reorganization of TEETOA, the lack of favorable water
expulsion, and the poorer hydration of the bound guests are responsible for
the weaker binding capacity of TEETOA relative to TEMOA. We intend to
investigate further these aspects of the binding mechanism in future work.

ATM and PMF both predict that G2D is one of the weakest binders for TEMOA and
TEETOA (Tables 3 and 4). G2D is expected to be frustrated in the bound state
because the bromine atom prefers to be in the cavity of the host, whereas the
oxide group strongly prefers to remain hydrated (Figure 3). The side chains of
both hosts prevent the hydration of the negative oxygen atom. This steric
hindrance is especially evident in TEETOA, which possesses four ethyl groups
on its outer ring. Due to its poor binding affinity, the deprotonated G2D is
not predicted to contribute significantly to binding despite its higher
concentration in solution at the experimental pH. Conversely, due to its
smaller desolvation penalty, both the PMF and ATM methods indicate that
protonated G2 (hereon G2P) is the strongest binder in the set for both TEMOA
and TEETOA (Tables 3 and 4). G2P is in fact predicted to be the dominant
species for binding even after factoring in the protonation penalty at the
experimental pH of 11.5.

The ATM free energy components $\Delta G_{1}$ and $\Delta G_{2}$ for each leg
of the ionic hosts (Table 3), being in the 40 to 50 kcal/mol range, are
significantly larger in magnitude than the resulting binding free energies.
These free energies correspond to the reversible work to reach the alchemical
intermediate state in which the guest interacts with both the receptor and the
solvent bulk intermediates. The high free energy of the alchemical
intermediate relative to the bound and solvated states suggests that the ionic
group can not be properly accommodated to simultaneously interact effectively
with both environments. This hypothesis is confirmed by the much smaller ATM
leg free energies for the neutral G2P guest. While large, the ATM leg free
energies of the ionic guests are expected to be significantly smaller than
those that would have obtained in a double-decoupling
calculationdeng2018comparing  that would involve displacing the guests into
vacuum where hydration interactions are completely removed. The statistical
uncertainties of the ATM binding free energy estimates, generally around $1/3$
of a kcal/mol, are relatively small.

While still moderate, the PMF binding free energy estimates (Table 4) come
with somewhat larger uncertainties than ATM. The source of uncertainties is
approximately equally split between the reversible work of releasing the
restraints (2nd column) and work of ligand extraction (3rd column). However,
in some cases (TEETOA-G2 and TEETOA-G5) the uncertainty of the work of
extraction is particularly large and probably indicative of sampling
bottlenecks at intermediate stages of the extraction process for this host.

Table 3: ATM absolute binding free energy estimates for the TEMOA and TEETOA
complexes.

Complex |  $\Delta G_{1}$a |  $\Delta G_{2}$a |  $\Delta G^{\circ}_{\rm site}$a |  $\Delta G^{\circ}_{b}$a
---|---|---|---|---
TEMOA-G1 | $53.27\pm 0.21$ | $45.69\pm 0.21$ | $0.87$ | $-6.71\pm 0.30$
TEMOA-G2D | $42.37\pm 0.18$ | $35.48\pm 0.21$ | $0.87$ | $-6.02\pm 0.28$
TEMOA-G2P | $22.57\pm 0.27$ | $8.60\pm 0.78$ | $0.87$ | $-13.10\pm 0.83$
TEMOA-G3 | $56.42\pm 0.18$ | $47.29\pm 0.18$ | $0.87$ | $-8.26\pm 0.25$
TEMOA-G4 | $53.13\pm 0.24$ | $43.63\pm 0.18$ | $0.87$ | $-8.63\pm 0.30$
TEMOA-G5 | $53.49\pm 0.24$ | $44.92\pm 0.18$ | $0.87$ | $-7.70\pm 0.30$
TEETOA-G1 | $51.65\pm 0.27$ | $49.71\pm 0.21$ | $0.87$ | $-1.07\pm 0.34$
TEETOA-G2D | $42.26\pm 0.24$ | $39.83\pm 0.27$ | $0.87$ | $-1.57\pm 0.36$
TEETOA-G2P | $22.31\pm 0.24$ | $13.48\pm 0.15$ | $0.87$ | $-7.95\pm 0.28$
TEETOA-G3 | $55.31\pm 0.24$ | $52.79\pm 0.18$ | $0.87$ | $-1.65\pm 0.30$
TEETOA-G4 | $52.28\pm 0.24$ | $48.90\pm 0.18$ | $0.87$ | $-2.51\pm 0.30$
TEETOA-G5 | $53.58\pm 0.21$ | $49.89\pm 0.18$ | $0.87$ | $-2.82\pm 0.28$

a In kcal/mol.

Table 4: PMF absolute free energy estimates for TEMOA and TEETOA complexes.

Complex |  $-\Delta G_{\rm restr}^{\rm bound}$a |  $[w(r_{\rm min})-w(r^{\ast})]$a |  $\Delta G_{\rm vibr}$a |  $\Delta G_{\rm restr}^{\rm bulk}$a |  $\Delta G^{\circ}_{b}$a
---|---|---|---|---|---
TEMOA-G1 | $-4.09\pm 0.23$ | $-12.27\pm 0.36$ | $0.24$ | $9.69$ | $-6.43\pm 0.43$
TEMOA-G2D | $-2.05\pm 0.33$ | $-11.01\pm 0.18$ | $0.12$ | $9.69$ | $-3.25\pm 0.38$
TEMOA-G2P | $-5.31\pm 0.78$ | $-17.12\pm 0.21$ | $0.17$ | $9.69$ | $-12.57\pm 0.81$
TEMOA-G3 | $-5.61\pm 0.30$ | $-12.83\pm 0.30$ | $0.04$ | $9.69$ | $-8.71\pm 0.42$
TEMOA-G4b | $-5.00\pm 0.47$ | $-13.72\pm 0.36$ | $0.24$ | $9.69$ | $-8.79\pm 0.59$
TEMOA-G5 | $-5.36\pm 0.81$ | $-12.74\pm 0.15$ | $0.26$ | $9.69$ | $-8.15\pm 0.82$
TEETOA-G1 | $-3.76\pm 0.60$ | $-7.60\pm 0.54$ | $0.28$ | $9.69$ | $-1.38\pm 0.81$
TEETOA-G2D | $-5.50\pm 0.84$ | $-5.25\pm 2.73$ | $0.20$ | $9.69$ | $-0.86\pm 2.86$
TEETOA-G2P | $-4.85\pm 0.57$ | $-14.51\pm 1.68$ | $0.25$ | $9.69$ | $-9.42\pm 1.77$
TEETOA-G3 | $-3.70\pm 0.24$ | $-7.36\pm 0.81$ | $-0.05$ | $9.69$ | $-1.42\pm 0.84$
TEETOA-G4 | $-3.77\pm 0.12$ | $-8.39\pm 0.75$ | $0.22$ | $9.69$ | $-2.25\pm 0.76$
TEETOA-G5 | $-4.47\pm 0.06$ | $-8.81\pm 1.89$ | $0.23$ | $9.69$ | $-3.36\pm 1.89$

a In kcal/mol.

###  3.2 Calculated Free Energy Estimates Relative to Experimental
Measurements

The two computational methods employed in this work reproduced the
experimental binding free energy estimates relatively well, particularly more
so for the TEMOA host than for the TEETOA host (Table 1). Both methods
correctly predict TEMOA-G2 as the highest affinity complex in the set with
good quantitative accuracy in the binding free energy predictions ($-8.41$
kcal/mol experimentally compared to calculated $-9.90$ and $-9.37$ kcal/mol
from ATM and PMF, respectively). Concomitantly, both methods correctly predict
relatively weak absolute binding free energies of -1.65 kcal/mol and -1.42
kcal/mol, respectively, for TEETOA-G3 which is an experimental non-binder.
Excluding TEETOA-G3, the least favorable binding affinity measurement was
obtained for TEETOA-G5, which is correctly scored as one of the weakest
complex by both computational methods. Overall, despite the the narrow range
of the moderate binding free energies, the computational rankings based on the
binding free energies are in good agreement with the experimental rankings
with a Kendall rank-order correlation coefficient of 0.69. (Table 2)

As illustrated in Figure 4 the calculated binding free energies are highly
correlated to the experimental values with Pearson $R^{2}$ correlation
coefficients of 89% and 83% for ATM and PMF, respectively (Table 2). The
calculations are also in reasonable quantitative agreement with the
experimental measurements with RMSE deviations of $1.71$ kcal/mol for ATM and
$1.79$ kcal/mol for PMF. Interestingly, the computational models tend to
overestimate the binding affinity of the TEMOA complexes and to underestimate
those of the complexes with TEETOA. The largest deviation occurs for TEETOA-G1
which has a moderate observed binding free energy of $-4.47$ kcal/mol, which
is underestimated by the computational predictions by around $-1$ kcal/mol. A
large deviation, but in the opposite direction, is also observed for TEMOA-G3
($-5.78$ kcal/mol experimentally compared to $-8.26$ and $-8.71$ kcal/mol
computationally) (Table v1). A poor prediction for this complex was expected
based on previous efforts with the GAFF/AM1-BCC force field with TIP3P
solvation used here.rizzi2018overview

In summary, the blinded predictions reported here were scored as among the
best of the SAMPL8 GDCC challenge and second only to those obtained with the
more accurate AMOEBA force fieldshi2021amoeba
(github.com/samplchallenges/SAMPL8/blob/master/host_guest/Analysis/Ranked_Accuracy).

Table 5: Binding free energy contributions of the protonated and deprotonated
G2 complexes to the ATM and PMF binding free estimates.

| TEMOA-G2/ATM | TEMOA-G2/PMF | TEETOA-G2/ATM | TEETOA-G2/PMF
---|---|---|---|---
$\Delta G_{b}^{\circ}$(HA)a | $-13.10\pm 0.83$ | $-12.57\pm 0.81$ | $-7.95\pm 0.28$ | $-9.42\pm 1.77$
$P_{0}({\rm HA})$ | $4.66\times 10^{-3}$ | $4.66\times 10^{-3}$ | $4.66\times 10^{-3}$ | $4.66\times 10^{-3}$
$\displaystyle P_{0}(HA)e^{-\beta\Delta G_{b}^{\circ}(HA)}$ | $1.65\times 10^{7}$ | $6.77\times 10^{7}$ | $2.92\times 10^{3}$ | $3.42\times 10^{4}$
$\Delta G_{b}^{\circ}({\rm A}^{-})$a | $-6.02\pm 0.28$ | $-3.25\pm 0.38$ | $-1.57\pm 0.36$ | $-0.86\pm 2.86$
$P_{0}({\rm A}^{-})$ | $0.995$ | $0.995$ | $0.995$ | $0.995$
$\displaystyle P_{0}(A^{-})e^{-\beta\Delta G_{b}^{\circ}(A^{-})}$ | $2.43\times 10^{4}$ | $232$ | $13.6$ | $4.22$
$\Delta G_{b}^{\circ}$a | $-9.90\pm 0.83$ | $-9.37\pm 0.81$ | $-4.76\pm 0.28$ | $-6.22\pm 1.8$

a In kcal/mol.

Figure 4: Linear regression of combined TEMOA and TEETOA predictions with ATM
and PMF.

##  4 Discussion and Conclusions

In this study, we employed two independent binding free energy approaches, the
newly developed alchemical transfer method (ATM)khuttan2021single ;
wu2021alchemical  and the well established PMF physical pathway
methoddeng2018comparing  to blindly predict the absolute binding affinities of
the host-guest systems as part of the SAMPL8 GDCC blind challenge. The SAMPL
series of community challenges has consistently yielded high-quality datasets
to test computational models of binding,geballe2010sampl2 ; mobley2014blind ;
amezcua2021sampl7 ; GallicchioSAMPL4 ; deng2016large ; pal2016SAMPL5  and we
decided to use it here to stringently validate the ATM and PMF methods in an
unbiased fashion.

Despite their radical differences in spirit and in practice, we find that the
calculated binding affinities from the two methods are in remarkable
quantitative agreement with an RMSE of only 0.6 kcal/mol and an $R^{2}$ of
$99$%. This level of agreement, well within statistical fluctuations, gives
high confidence in the theoretical foundations and in the correctness of
implementation of each approach. The level of consistency of the computational
methods also adds confidence that their predictions are unbiased and primarily
reflective of the force field model.

We find that the standard GAFF/AM1-BCC/TIP3P model employed here tends to
overestimate the binding free energies of strongly bound complexes while it
tends to understimate those of more weakly bound complexes, as also indicated
by the larger than one slope of the linear regressions (Tables 1, 2). While it
may be a result, in this case, of specific aspects of the TEMOA and TEETOA
hosts, this trend has been generally observed with this force field
combination.rizzi2018overview  The more accurate AMOEBA force
fieldshi2021amoeba  appears to correctly predict these trends
(github.com/samplchallenges/SAMPL8/blob/master/host_guest/Analysis/Ranked_Accuracy).

The stringent blinded test conducted in this work is a further validation of
the ATM binding free energy method that we have recently
proposed.wu2021alchemical  ATM, implemented on top of the versatile OpenMM
molecular dynamics engine,eastman2017openmm  promises to provide an accurate
and streamlined route to absolutewu2021alchemical  and relative binding free
calculations.Azimi2021RBFE  While alchemical, ATM, similar to the PMF pathway
method,deng2018comparing  makes use of a single simulation system, and it
avoids problematic vacuum intermediates and the splitting of the alchemical
path into electrostatic and non-electrostatic transformations. ATM also does
not require soft-core pair potentials and modifications of energy routines,
and can be easily implemented as a controlling routine on top of existing
force routines of MD engines.

In summary, this work provides a rare blinded and stringent test of binding
free energy models. It shows that the application of sound statistical
mechanics theories of binding and careful modeling of chemical systems can
lead to reliable predictions limited only by the quality of the force field
model.

##  5 Acknowledgements

We acknowledge support from the National Science Foundation (NSF CAREER
1750511 to E.G.). Molecular simulations were conducted on the Comet and
Expanse GPU clusters at the San Diego Supercomputing Center supported by NSF
XSEDE award TG-MCB150001. We appreciate the National Institutes of Health for
its support of the SAMPL project via R01GM124270 to David L. Mobley.

## References

  * (1) Matthew T Geballe, A Geoffrey Skillman, Anthony Nicholls, J Peter Guthrie, and Peter J Taylor.  The SAMPL2 blind prediction challenge: introduction and overview.  J. Comp. Aided Mol. Des., 24(4):259–279, 2010.
  * (2) David L Mobley, Shuai Liu, Nathan M Lim, Karisa L Wymer, Alexander L Perryman, Stefano Forli, Nanjie Deng, Justin Su, Kim Branson, and Arthur J Olson.  Blind prediction of hiv integrase binding from the SAMPL4 challenge.  J. Comp. Aided Mol. Des., pages 1–19, 2014.
  * (3) Martin Amezcua, Léa El Khoury, and David L Mobley.  SAMPL7 host–guest challenge overview: assessing the reliability of polarizable and non-polarizable methods for binding free energy calculations.  J. Comp.-Aid. Mol. Des., 35(1):1–35, 2021.
  * (4) David L Mobley and Michael K Gilson.  Predicting binding free energies: frontiers and benchmarks.  Ann. Rev. Bioph., 46:531–558, 2017.
  * (5) William L Jorgensen.  Efficient drug lead discovery and optimization.  Acc Chem Res, 42:724–733, 2009.
  * (6) Kira A Armacost, Sereina Riniker, and Zoe Cournia.  Novel directions in free energy methods and applications, 2020.
  * (7) E. Gallicchio and R. M. Levy.  Prediction of SAMPL3 host-guest affinities with the binding energy distribution analysis method (BEDAM).  J. Comp. Aided Mol. Design., 25:505–516, 2012.
  * (8) Emilio Gallicchio, Haoyuan Chen, He Chen, Michael Fitzgerald, Yang Gao, Peng He, Malathi Kalyanikar, Chuan Kao, Beidi Lu, Yijie Niu, Manasi Pethe, Jie Zhu, and Ronald M Levy.  BEDAM binding free energy predictions for the SAMPL4 octa-acid host challenge.  J. Comp. Aided Mol. Des., 29:315–325, 2015.
  * (9) Emilio Gallicchio, Nanjie Deng, Peng He, Alexander L. Perryman, Daniel N. Santiago, Stefano Forli, Arthur J. Olson, and Ronald M. Levy.  Virtual screening of integrase inhibitors by large scale binding free energy calculations: the SAMPL4 challenge.  J. Comp. Aided Mol. Des., 28:475–490, 2014.
  * (10) Nanjie Deng, William F Flynn, Junchao Xia, RSK Vijayan, Baofeng Zhang, Peng He, Ahmet Mentes, Emilio Gallicchio, and Ronald M Levy.  Large scale free energy calculations for blind predictions of protein–ligand binding: the d3r grand challenge 2015.  J. Comp.-Aided Mol. Des., 30(9):743–751, 2016.
  * (11) Rajat Kumar Pal, Kamran Haider, Divya Kaur, William Flynn, Junchao Xia, Ronald M. Levy, Tetiana Taran, Lauren Wickstrom, Tom Kurtzman, and Emilio Gallicchio.  A combined treatment of hydration and dynamical effects for the modeling of host-guest binding thermodynamics: The SAMPL5 blinded challenge.  J. Comp. Aided Mol. Des., 31:29–44, 2016.
  * (12) Joe Z Wu, Solmaz Azimi, Sheenam Khuttan, Nanjie Deng, and Emilio Gallicchio.  Alchemical transfer approach to absolute binding free energy estimation.  J. Chem. Theory Comput., 17:3309, 2021.
  * (13) Nanjie Deng, Di Cui, Bin W Zhang, Junchao Xia, Jeffrey Cruz, and Ronald Levy.  Comparing alchemical and physical pathway methods for computing the absolute binding free energy of charged ligands.  Phys. Chem. Chem. Phys., 20(25):17081–17092, 2018.
  * (14) Paolo Suating, Thong T Nguyen, Nicholas E Ernst, Yang Wang, Jacobs H Jordan, Corinne LD Gibb, Henry S Ashbaugh, and Bruce C Gibb.  Proximal charge effects on guest binding to a non-polar pocket.  Chemical Science, 11(14):3656–3663, 2020.
  * (15) Paweł Śledź and Amedeo Caflisch.  Protein structure-based drug design: from docking to molecular dynamics.  Curr. Op. Struct. Biol., 48:93–102, 2018.
  * (16) Thomas Seidel, Oliver Wieder, Arthur Garon, and Thierry Langer.  Applications of the pharmacophore concept in natural product inspired drug design.  Molecular Informatics, 39(11):2000059, 2020.
  * (17) 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.
  * (18) Emilio Gallicchio and Ronald M Levy.  Recent theoretical and computational advances for modeling protein-ligand binding affinities.  Adv. Prot. Chem. Struct. Biol., 85:27–80, 2011.
  * (19) Zoe Cournia, Bryce K Allen, Thijs Beuming, David A Pearlman, Brian K Radak, and Woody Sherman.  Rigorous free energy simulations in virtual screening.  Journal of Chemical Information and Modeling, 2020.
  * (20) Emilio Gallicchio.  Free energy-based computational methods for the study of protein-peptide binding equilibria.  In Thomas Simonson, editor, Computational Peptide Science: Methods and Protocols, Methods in Molecular Biology. Springer Nature, 2021.
  * (21) Emilio Gallicchio, Mauro Lapelosa, and Ronald M. Levy.  Binding energy distribution analysis method (BEDAM) for estimation of protein-ligand binding affinities.  J. Chem. Theory Comput., 6:2961–2977, 2010.
  * (22) Zhiqiang Tan, Emilio Gallicchio, Mauro Lapelosa, and Ronald M. Levy.  Theory of binless multi-state free energy estimation with applications to protein-ligand binding.  J. Chem. Phys., 136:144102, 2012.
  * (23) S Boresch, F Tettinger, M Leitgeb, and M Karplus.  Absolute binding free energies: A quantitative approach for their calculation.  J. Phys. Chem. B, 107:9535–9551, 2003.
  * (24) Sander Pronk, Szilárd Páll, Roland Schulz, Per Larsson, Pär Bjelkmar, Rossen Apostolov, Michael R Shirts, Jeremy C Smith, Peter M Kasson, David van der Spoel, Berk Hess, and Erik Lindahl.  Gromacs 4.5: a high-throughput and highly parallel open source molecular simulation toolkit.  Bioinformatics, 29:845–854, 2013.
  * (25) S Khuttan, Solmaz Azimi, Joe Z Wu, and E Gallicchio.  Alchemical transformations for concerted hydration free energy estimation with explicit solvation.  J. Chem. Phys, 154:054103, 2021.
  * (26) Peter Eastman, Jason Swails, John D Chodera, Robert T McGibbon, Yutong Zhao, Kyle A Beauchamp, Lee-Ping Wang, Andrew C Simmonett, Matthew P Harrigan, Chaya D Stern, et al.  Openmm 7: Rapid development of high performance algorithms for molecular dynamics.  PLoS Comp. Bio., 13(7):e1005659, 2017.
  * (27) Emilio Gallicchio, Junchao Xia, William F Flynn, Baofeng Zhang, Sade Samlalsingh, Ahmet Mentes, and Ronald M Levy.  Asynchronous replica exchange software for grid and heterogeneous computing.  Computer Physics Communications, 196:236–246, 2015.
  * (28) Andrea Rizzi, Steven Murkli, John N McNeill, Wei Yao, Matthew Sullivan, Michael K Gilson, Michael W Chiu, Lyle Isaacs, Bruce C Gibb, David L Mobley, et al.  Overview of the sampl6 host–guest binding affinity prediction challenge.  J. Comp.-Aid. Mol. Des., 32(10):937–963, 2018.
  * (29) Yuanjun Shi, Marie L Laury, Zhi Wang, and Jay W Ponder.  Amoeba binding free energies for the sampl7 trimertrip host–guest challenge.  J. Comp.-Aid. Mol. Des., 35(1):79–93, 2021.
  * (30) Solmaz Azimi, Sheenam Khuttan, Joe Z. Wu, Rajat Pal, and Emilio Gallicchio.  Relative binding free energy calculations for ligands with diverse scaffolds with the alchemical transfer method.  ArXiv Preprint, XXX:XXX–XXX, 2021.


# Learning to Explicitate Connectives with Seq2Seq Network
for Implicit Discourse Relation Classification

Wei Shi† and Vera Demberg†,‡
†Dept. of Language Science and Technology
‡Dept. of Mathematics and Computer Science, Saarland University
Saarland Informatics Campus, 66123 Saarbrücken, Germany
{w.shi<EMAIL_ADDRESS>
###### Abstract

Implicit discourse relation classification is one of the most difficult steps
in discourse parsing. The difficulty stems from the fact that the coherence
relation must be inferred based on the content of the discourse relational
arguments. Therefore, an effective encoding of the relational arguments is of
crucial importance. We here propose a new model for implicit discourse
relation classification, which consists of a classifier, and a sequence-to-
sequence model which is trained to generate a representation of the discourse
relational arguments by trying to predict the relational arguments including a
suitable implicit connective. Training is possible because such implicit
connectives have been annotated as part of the PDTB corpus. Along with a
memory network, our model could generate more refined representations for the
task. And on the now standard 11-way classification, our method outperforms
the previous state of the art systems on the PDTB benchmark on multiple
settings including cross validation.

##  1 Introduction

Discourse relations describe the logical relation between two
sentences/clauses. When understanding a text, humans infer discourse relation
between text segmentations. They reveal the structural organization of text,
and allow for additional inferences. Many natural language processing tasks,
such as machine translation, question-answering, automatic summarization,
sentiment analysis, and sentence embedding learning, can also profit from
having access to discourse relation information. Recent years have seen more
and more works on this topic, including two CoNNL shared tasks (Xue et al.,
2015, 2016).

Penn Discourse Tree Bank (Prasad et al., 2008, PDTB) provides lexically-
grounded annotations of discourse relations and their two discourse relational
arguments (i.e., two text spans). Discourse relations are sometimes signaled
by explicit discourse markers (e.g., because, but). Example 1 shows an
explicit discourse relation marked by “because”; the presence of the
connective makes it possible to classify the discourse relation with high
reliability: Miltsakaki et al. (2005) reported an accuracy of 93.09% for 4-way
classification of explicits.

Discourse relations are however not always marked by an explicit connective.
In fact, implicit discourse relations (i.e. relations not marked by an
explicit discourse cue) outnumber explicit discourse relations in naturally
occurring text. Readers can still infer these implicit relations, but
automatic classification becomes a lot more difficult in these cases, and
represents the main bottleneck in discourse parsing today. Example 2 shows an
implicit contrastive relation which can be inferred from the two text spans
that have been marked Arg1 and Arg2. When annotating implicit relations in the
PDTB, annotators were asked to first insert a connective which expresses the
relation, and then annotate the relation label. This procedure was introduced
to achieve higher inter-annotator agreement for implicit relations between
human annotators. In the approach taken in this paper, our model mimics this
procedure by being trained to explicitate the discouse relation, i.e. to
insert a connective as a secondary task.

  1. 1.

_[I refused to pay the cobbler the full $95] Arg1 because [He did poor
work.]Arg2
_

_— Explicit, Contingency.Cause_

  2. 2.

_[In the energy mix of the future, bio-energy will also have a key role to
play in boosting rural employment and the rural economy in Europe .] Arg1
(Implicit = However) [At the same time , the promotion of bio-energy must not
lead to distortions of competition.]Arg2
_

_— Implicit, Comparison.Contrast_

The key in implicit discourse relation classification lies in extracting
relevant information for the relation label from (the combination of) the
discourse relational arguments. Informative signals can consist of surface
cues, as well as the semantics of the relational arguments. Statistical
approaches have typically relied on linguistically informed features which
capture both of these aspects, like temporal markers, polarity tags, Levin
verb classes and sentiment lexicons, as well as the Cartesian products of the
word tokens in the two arguments (Lin et al., 2009). More recent efforts use
distributed representations with neural network architectures (Qin et al.,
2016a).

The main question in designing neural networks for discourse relation
classification is how to get the neural networks to effectively encode the
discourse relational arguments such that all of the aspects relevant to the
classification of the relation are represented, in particular in the face of
very limited amounts of annotated training data, see e.g. Rutherford et al.
(2017). The crucial intuition in the present paper is to make use of the
annotated implicit connectives in the PDTB: in addition to the typical
relation label classification task, we also train the model to encode and
decode the discourse relational arguments, and at the same time predict the
implicit connective. This novel secondary task forces the internal
representation to more completely encode the semantics of the relational
arguments (in order to allow the model to decode later), and to make a more
fine-grained classification (predicting the implicit connective) than is
necessary for the overall task. This more fine-grained task thus aims to force
the model to represent the discourse relational arguments in a way that allows
the model to also predict a suitable connective. Our overall discourse
relation classifier combines representations from the relational arguments as
well as the hidden representations generated as part of the encoder-decoder
architecture to predict relation labels. What’s more, with an explicit memory
network, the network also has access to history representations and acquire
more explicit context knowledge. We show that our method outperforms previous
approaches on the 11-way classification on the PDTB 2.0 benchmark.

The remaining of the paper is organized as follows: Section 2 discusses
related work; Section 3 describes our proposed method; Section 4 gives the
training details and experimental results, which is followed by conclusion and
future work in section 5.

##  2 Related Work

###  2.1 Implicit Discourse Relation Classification

Implicit discourse relation recognition is one of the most important
components in discourse parsing. With the release of PDTB (Prasad et al.,
2008), the largest available corpus which annotates implicit examples with
discourse relation labels and implicit connectives, a lot of previous works
focused on typical statistical machine learning solutions with manually
crafted sparse features (Rutherford and Xue, 2014).

Recently, neural networks have shown an advantage of dealing with data
sparsity problem, and many deep learning methods have been proposed for
discourse parsing, including convolutional (Zhang et al., 2015), recurrent (Ji
et al., 2016), character-based (Qin et al., 2016a), adversarial (Qin et al.,
2017) neural networks, and pair-aware neural sentence modeling (Cai and Zhao,
2017). Multi-task learning has also been shown to be beneficial on this task
(Lan et al., 2017).

However, most neural based methods suffer from insufficient annotated data.Wu
et al. (2016) extracted bilingual-constrained synthetic implicit data from a
sentence-aligned English-Chinese corpus. Shi et al. (2017, 2018) proposed to
acquire additional training data by exploiting explicitation of connectives
during translation. Explicitation refers to the fact that translators
sometimes add connectives into the text in the target language which were not
originally present in the source language. They used explicitated connectives
as a source of weak supervision to obtain additional labeled instances, and
showed that this extension of the training data leads to substantial
performance improvements.

The huge gap between explicit and implicit relation recognition (namely, 50%
vs. 90% in 4-way classification) also motivates to incorporate connective
information to guide the reasoning process. Zhou et al. (2010) used a language
model to automatically insert discourse connectives and leverage the
information of these predicted connectives. The approach which is most similar
in spirit to ours, Qin et al. (2017), proposed a neural method that
incorporates implicit connectives in an adversarial framework to make the
representation as similar as connective-augmented one and showed that the
inclusion of implicit connectives could help to improve classifier
performance.

###  2.2 Sequence-to-sequence Neural Networks

Sequence to sequence model is a general end-to-end approach to sequence
learning that makes minimal assumptions on the sequence structure, and firstly
proposed by Sutskever et al. (2014). It uses multi-layered Long Short-Term
Memory (LSTM) or Gated Recurrent Units (GRU) to map the input sequence to a
vector with a fixed dimensionality, and then decode the target sequence from
the vector with another LSTM / GRU layer.

Sequence to sequence models allow for flexible input/output dynamics and have
enjoyed great success in machine translation and have been broadly used in
variety of sequence related tasks such as Question Answering, named entity
recognition (NER) / part of speech (POS) tagging and so on.

If the source and target of a sequence-to-sequence model are exactly the same,
it is also called Auto-encoder, Dai and Le (2015) used a sequence auto-encoder
to better represent sentence in an unsupervised way and showed impressive
performances on different tasks. The main difference between our model and
this one is that we have different input and output (the output contains a
connective while the input doesn’t). In this way, the model is forced to
explicitate implicit relation and try to learn the latent pattern and
discourse relation between implicit arguments and connectives and then
generate more discriminative representations.

Figure 1: The Architecture of Proposed Model.

##  3 Methodology

Our model is based on the sequence-to-sequence model used for machine
translation (Luong et al., 2015), an adaptation of an LSTM (Hochreiter and
Schmidhuber, 1997) that encodes a variable length input as a fix-length
vector, then decodes it into a variable length of outputs. As illustrated in
Figure 1, our model consists of three components: Encoder, Decoder and
Discourse Relation Classifier. We here use different LSTMs for the encoding
and decoding tasks to help keep the independence between those two parts.

The task of implicit discourse relation recognition is to recognize the senses
of the implicit relations, given the two arguments. For each discourse
relation instance, The Penn Discourse Tree Bank (PDTB) provides two arguments
($Arg_{1}$, $Arg_{2}$) along with the discourse relation (Rel) and manually
inserted implicit discourse connective ($Conn_{i}$). Here is an implicit
example from section 0 in PDTB:

  1. 3.

$\mathbf{Arg_{1}}$: This is an old story.
$\mathbf{Arg_{2}}$: We’re talking about years ago before anyone heard of
asbestos having any questionable properties.
$\mathbf{Conn_{i}}$: in fact
$\mathbf{Rel}$: Expansion.Restatement

During training, the input and target sentences for sequence-to-sequence
neural network are $\left[\textit{$Arg_{1}$};\textit{$Arg_{2}$}\right]$ and
$\left[\textit{$Arg_{1}$};\textit{$Conn_{i}$};\textit{$Arg_{2}$}\right]$
respectively, where “;” denotes concatenation.

###  3.1 Model Architecture

####  3.1.1 Encoder

Given a sequence of words, an encoder computes a joint representation of the
whole sequence.

After mapping tokens to Word2Vec embedding vectors (Mikolov et al., 2013), a
LSTM recurrent neural network processes a variable-length sequence
$x=(x_{1},x_{2},...,x_{n})$. At time step $t$, the state of memory cell
$c_{t}$ and hidden $h_{t}$ are calculated with the Equations 1:

$\small\begin{gathered}\left[\begin{array}[]{c}i_{t}\\\ f_{t}\\\ o_{t}\\\
\hat{c_{t}}\end{array}\right]=\left[\begin{array}[]{c}\sigma\\\ \sigma\\\
\sigma\\\ \tanh\end{array}\right]W\cdot[h_{t-1},x_{t}]\\\ c_{t}=f_{t}\odot
c_{t-1}+i_{t}\odot\hat{c_{t}}\\\ h_{t}=o_{t}\odot\tanh(c_{t})\\\
\end{gathered}$ (1)

where $x_{t}$ is the input at time step $t$, $i$, $f$ and $o$ are the input,
forget and output gate activation respectively. $\hat{c_{t}}$ denotes the
current cell state, $\sigma$ is the logistic sigmoid function and $\odot$
denotes element-wise multiplication. The LSTM separates the memory $c$ from
the hidden state $h$, which allows for more flexibility incombining new inputs
and previous context.

For the sequence modeling tasks, it is beneficial to have access to the past
context as well as the future context. Therefore, we chose a bidirectional
LSTM as the encoder and the output of the word at time-step $t$ is shown in
the Equation 2. Here, element-wise sum is used to combine the forward and
backward pass outputs.

$\small h_{t}=\left[\overrightarrow{h_{t}}\oplus\overleftarrow{h_{t}}\right]$
(2)

Thus we get the output of encoder:

$\small h_{e}=\left[h^{e}_{1},h^{e}_{2},...,h^{e}_{n}\right]$ (3)

####  3.1.2 Decoder

With the representation from the encoder, the decoder tries to map it back to
the targets space and predicts the next words.

Here we used a separate LSTM recurrent network to predict the target words.
During training, target words are fed into the LSTM incrementally and we get
the outputs from decoder LSTM:

$\small h_{d}=\left[h^{d}_{1},h^{d}_{2},...,h^{d}_{n}\right]$ (4)

#### Global Attention

In each time-step in decoding, it’s better to consider all the hidden states
of the encoder to give the decoder a full view of the source context. So we
adopted the global attention mechanism proposed in Luong et al. (2015). For
time step $t$ in decoding, context vector $c_{t}$ is the weighted average of
$h_{e}$, the weights for each time-step are calculated with $h_{t}^{d}$ and
$h_{e}$ as illustrated below:

$\small\alpha_{t}=\frac{\exp({h_{t}^{d}}^{\top}\mathbf{W_{\alpha}}h_{e})}{\sum\limits_{t=1}^{n}\exp({h_{t}^{d}}^{\top}\mathbf{W_{\alpha}}h_{e})}$
(5) $\small c_{t}=\alpha h_{e}$ (6)

#### Word Prediction

Context vector $c_{t}$ captured the relevant source side information to help
predict the current target word $y_{t}$. We employ a concatenate layer with
activation function $\tanh$ to combine context vector $c_{t}$ and hidden state
of decoder $h^{d}_{t}$ at time-step t as follows:

$\small\hat{h^{d}_{t}}=\tanh(\mathbf{W_{c}}\left[c_{t};h^{d}_{t}\right])$ (7)

Then the predictive vector is fed into the softmax layer to get the predicted
distribution $\hat{p}(y_{t}|s)$ of the current target word.

$\small\begin{gathered}\hat{p}(y_{t}|s)=softmax(\mathbf{W}_{s}\hat{h_{d}}+\mathbf{b}_{s})\\\
\hat{y_{t}}=\arg\max_{y}\hat{p}(y_{t}|s)\end{gathered}$ (8)

After decoding, we obtain the predictive vectors for the whole target sequence
$\hat{h_{d}}=\left[h^{d}_{1},h^{d}_{2},...,h^{d}_{n}\right]$. Ideally, it
contains the information of exposed implicit connectives.

#### Gated Interaction

In order to predict the coherent discourse relation of the input sequence, we
take both the $h_{encoder}$ and the predictive word vectors $h_{d}$ into
account. K-max pooling can “draw together” features that are most
discriminative and among many positions apart in the sentences, especially on
both the two relational arguments in our task here; this method has been
proved to be effective in choosing active features in sentence modeling
(Kalchbrenner et al., 2014). We employ an average k-max pooling layer which
takes average of the top k-max values among the whole time-steps as in
Equation 9 and 10:

$\small\bar{h}_{e}=\frac{1}{k}\sum\limits^{k}_{i=1}topk(h_{e})$ (9)
$\small\bar{h}_{d}=\frac{1}{k}\sum\limits^{k}_{i=1}topk(\hat{h^{d}})$ (10)

$\bar{h}_{e}$ and $\bar{h}_{d}$ are then combined using a linear layer (Lan et
al., 2017). As illustrated in Equation 11, the linear layer acts as a gate to
determine how much information from the sequence-to-sequence network should be
mixed into the original sentence’s representations from the encoder. Compared
with bilinear layer, it also has less parameters and allows us to use high
dimensional word vectors.

$\small
h^{*}=\bar{h}_{e}\oplus\sigma(\mathbf{W}_{i}\bar{h}_{d}+\mathbf{b}_{i})$ (11)

#### Explicit Context Knowledge

To further capture common knowledge in contexts, we here employ a memory
network proposed in Liu et al. (2018), to get explicit context representations
of contexts training examples. We use a memory matrix $M\in R^{K\times N}$,
where $K,N$ denote hidden size and number of training instances respectively.
During training, the memory matrix remembers the information of training
examples and then retrieves them when predicting labels.

Given a representation $h^{*}$ from the interaction layer, we generate a
knowledge vector by weighted memory reading:

$\small k=Msoftmax(M^{T}h^{*})$ (12)

We here use dot product attention, which is faster and space-efficient than
additive attention, to calculate the scores for each training instances. The
scores are normalized with a softmax layer and the final knowledge vector is a
weighted sum of the columns in memory matrix $M$.

Afterwards, the model predicts the discourse relation using a softmax layer.

$\small\begin{gathered}\hat{p}(r|s)=softmax(\mathbf{W}_{r}[k;h^{*}]+\mathbf{b}_{r})\\\
\hat{r}=\arg\max_{y}\hat{p}(r|s)\end{gathered}$ (13)

###  3.2 Multi-objectives

In our model, the decoder and the discourse relation classifier have different
objectives. For the decoder, the objective consists of predicting the target
word at each time-step. The loss function is calculated with masked cross
entropy with $\mathtt{L2}$ regularization, as follows:

$\small\mathit{Loss_{de}}=-\frac{1}{n}\sum\limits^{n}_{t=1}y_{t}\log(\hat{p_{y}})+\frac{\lambda}{2}\parallel\theta_{de}\parallel^{2}_{2}$
(14)

where $y_{t}$ is one-hot represented ground truth of target words,
$\hat{p_{y}}$ is the estimated probabilities for each words in vocabulary by
softmax layer, $n$ denotes the length of target sentence. $\lambda$ is a
hyper-parameter of $L2$ regularization and $\theta$ is the parameter set.

The objective of the discourse relation classifier consists of predicting the
discourse relations. A reasonable training objective for multiple classes is
the categorical cross-entropy loss. The loss is formulated as:

$\small\mathit{Loss_{cl}}=-\frac{1}{m}\sum\limits^{m}_{i=1}r_{i}\log(\hat{p_{r}})+\frac{\lambda}{2}\parallel\theta_{cl}\parallel^{2}_{2}$
(15)

where $r_{i}$ is one-hot represented ground truth of discourse relation
labels, $\hat{p_{r}}$ denotes the predicted probabilities for each relation
class by softmax layer, $m$ is the number of target classes. Just like above,
$\lambda$ is a hyper-parameter of $L2$ regularization.

For the overall loss of the whole model, we set another hyper-parameter $w$ to
give these two objective functions different weights. Larger $w$ means that
more importance is placed on the decoder task.

$\small\mathit{Loss}=\mathit{w}\cdot\mathit{Loss_{de}}+\mathit{(1-w)}\cdot\mathit{Loss_{cl}}$
(16)

###  3.3 Model Training

To train our model, the training objective is defined by the loss function we
introduced above. We use Adaptive Moment Estimation (Adam) (Kingma and Ba,
2014) with different learning rate for different parts of the model as our
optimizer. Dropout layers are applied after the embedding layer and also on
the top feature vector before the softmax layer in the classifier. We also
employ $L_{2}$ regularization with small $\lambda$ in our objective functions
for preventing over-fitting. The values of the hyper-parameters, are provided
in Table 2. The model is trained firstly to minimize the loss in Equation 14
until convergence, we use scheduled sampling (Bengio et al., 2015) during
training to avoid “teacher-forcing problem”. And then to minimize the joint
loss in Equation 16 to train the implicit discourse relation classifier.

##  4 Experiments and Results

###  4.1 Experimental Setup

We evaluate our model on the PDTB. While early work only evaluated
classification performance between the four main PDTB relation classes, more
recent work including the CoNLL 2015 and 2016 shared tasks on Shallow
Discourse Parsing (Xue et al., 2015, 2016) have set the standard to second-
level classification. The second-level classification is more useful for most
downstream tasks. Following other works we directly compare to in our
evaluation, we here use the setting where AltLex, EntRel and NoRel tags are
ignored. About 2.2% of the implicit relation instances in PDTB have been
annotated with two relations, these are considered as two training instances.

To allow for full comparability to earlier work, we here report results for
three different settings. The first one is denoted as PDTB-Lin (Lin et al.,
2009); it uses sections 2-21 for training, 22 as dev and section 23 as test
set. The second one is labeled PDTB-Ji (Ji and Eisenstein, 2015), and uses
sections 2-20 for training, 0-1 as dev and evaluates on sections 21-22. Our
third setting follows the recommendations of Shi and Demberg (2017), and
performs 10-fold cross validation on the whole corpus (sections 0-23). Table 1
shows the number of instances in train, development and test set in different
settings.

Settings | Train | Dev | Test
---|---|---|---
PDTB-Lin | 13351 | 515 | 766
PDTB-Ji | 12826 | 1165 | 1039
Cross valid. per fold avg. | 12085 | 1486 |  1486111Cross-validation allows us to test on all 15057 instances.
Table 1: Numbers of train, development and test set on different settings for
11-way classification task. Instances annotated with two labels are double-
counted and some relations with few instances have been removed.

The advantage of the cross validation approach is that it addresses problems
related to the small corpus size, as it reports model performance across all
folds. This is important, because the most frequently used test set (PDTB-Lin)
contains less than 800 instances; taken together with a lack in the community
to report mean and standard deviations from multiple runs of neural networks
(Reimers and Gurevych, 2018), the small size of the test set makes reported
results potentially unreliable.

#### Preprocessing

We first convert tokens in PDTB to lowercase and normalize strings, which
removes special characters. The word embeddings used for initializing the word
representations are trained with the CBOW architecture in
Word2Vec222https://code.google.com/archive/p/word2vec/ (Mikolov et al., 2013)
on PDTB training set. All the weights in the model are initialized with
uniform random.

To better locate the connective positions in the target side, we use two
position indicators ($\langle conn\rangle$, $\langle/conn\rangle$) which
specify the starting and ending of the connectives (Zhou et al., 2016), which
also indicate the spans of discourse arguments.

Since our main task here is not generating arguments, it is better to have
representations generated by correct words rather than by wrongly predicted
ones. So at test time, instead of using the predicted word from previous time
step as current input, we use the source sentence as the decoder’s input and
target. As the implicit connective is not available at test time, we use a
random vector, which we used as “impl_conn” in Figure 2, as a placeholder to
inform the sequence that the upcoming word should be a connective.

#### Hyper-parameters

There are several hyper-parameters in our model, including dimension of word
vectors $d$, two dropout rates after embedding layer $q_{1}$ and before
softmax layer $q_{2}$, two learning rates for encoder-decoder $lr_{1}$ and for
classifier $lr_{2}$, top $k$ for k-max pooling layer, different weights $w$
for losses in Equation (16) and $\lambda$ denotes the coefficient of
regularizer, which controls the importance of the regularization term, as
shown in Table 2.

$d$ | $q_{1}$ | $q_{2}$ | ${lr}_{1}$ | ${lr}_{2}$ | $k$ | $w$ | $\lambda$
---|---|---|---|---|---|---|---
100 | 0.5 | 0.2 | $2.5e^{-3}$ | $5e^{-3}$ | 5 | 0.2 | $5e^{-4}$
Table 2: Hyper-parameter settings.

###  4.2 Experimental Results

We compare our models with six previous methods, as shown in Table 3. The
baselines contain feature-based methods (Lin et al., 2009), state-of-the-art
neural networks (Qin et al., 2016a; Cai and Zhao, 2017), including the
adversarial neural network that also exploits the annotated implicit
connectives (Qin et al., 2017), as well as the data extension method based on
using explicitated connectives from translation to other languages (Shi et
al., 2017).

Additionally, we ablate our model by taking out the prediction of the implicit
connective in the sequence to sequence model. The resulting model is labeled
Auto-Encoder in Table 3. And seq2seq network without knowledge memory, which
means we use the output of gated interaction layer to predict the label
directly, as denoted as Seq2Seq w/o Mem Net.

Methods | PDTB-Lin | PDTB-Ji | Cross Validation
---|---|---|---
Majority class | 26.11 | 26.18 | 25.59
Lin et al. (2009) | 40.20 | - | -
Qin et al. (2016a) | 43.81 | 45.04 | -
Cai and Zhao (2017) | - | 45.81 | -
Qin et al. (2017) | 44.65 | 46.23 | -
Shi et al. (2017) (with extra data) | 45.50 | - | 37.84
Encoder only (Bi-LSTM) (Shi et al., 2017) | 34.32 | - | 30.01
Auto-Encoder | 43.86 | 45.43 | 39.50
Seq2Seq w/o Mem Net | 45.75 | 47.05 | 40.29
Proposed Method | 45.82 | 47.83 | 41.29
Table 3: Accuracy (%) of implicit discourse relations on PDTB-Lin, PDTB-Ji and
Cross Validation Settings for multi-class classification.

Our proposed model outperforms the other models in each of the settings.
Compared with performances in Qin et al. (2017), although we share the similar
idea of extracting highly discriminative features by generating connective-
augmented representations for implicit discourse relations, our method
improves about 1.2% on setting PDTB-Lin and 1.6% on the PDTB-Ji setting. The
importance of the implicit connective is also illustrated by the fact that the
“Auto-Encoder” model, which is identical to our model except it does not
predict the implicit connective, performs worse than the model which does.
This confirms our initial hypothesis that training with implicit connectives
helps to expose the latent discriminative features in the relational
arguments, and generates more refined semantic representation. It also means
that, to some extent, purely increasing the size of tunable parameters is not
always helpful in this task and trying to predict implicit connectives in the
decoder does indeed help the model extract more discriminative features for
this task. What’s more, we can also see that without the memory network, the
performances are also worse, it shows that with the concatenation of knowledge
vector, the training instance may be capable of finding related instances to
get common knowledge for predicting implicit relations. As Shi and Demberg
(2017) argued that it is risky to conclude with testing on such small test
set, we also run cross-validation on the whole PDTB. From Table 3, we have the
same conclusion with the effectiveness of our method, which outperformed the
baseline (Bi-LSTM) with more than 11% points and 3% compared with Shi et al.
(2017) even though they have used a very large extra corpus.

For the sake of obtaining a better intuition on how the global attention works
in our model, Figure 2 demonstrates the weights of different time-steps in
attention layer from the decoder. The weights show how much importance the
word attached to the source words while predicting target words. We can see
that without the connective in the target side of test, the word filler still
works as a connective to help predict the upcoming words. For instance, the
true discourse relation for the right-hand example is Expansion.Alternative,
at the word filler’s time-step, it attached more importance on the negation
“don’t” and “tastefully appointed”. It means the current representation could
grasp the key information and try to focus on the important words to help with
the task. Here we see plenty room for adapting this model to discourse
connective prediction task, we would like to leave this to the future work.

We also try to figure out which instances’ representations have been chosen
from the memory matrix while predicting. Table 4 shows two examples and their
context instances with top 2 memory attentions among the whole training set.
We can see that both examples show that the memory attention attached more
importance on the same relations. This means that with the Context Memory, the
model could facilitate the discourse relation prediction by choosing examples
that share similar semantic representation and discourse relation during
prediction.

Figure 2: Visualization of attention weights during predicting target sentence
in train and test, x-axis denotes the source sentence and the y-axis is the
targets. First two figures are examples from training set with implicit
connectives inside, while the following one, in which the implicit connective
has been replaced by the word filler “impl_conn”, is from test. In recent
years, U.S. steelmakers have supplied about 80% of the 100 million tons of
steel used annually by the nation. (in addition,) Of the remaining 20% needed,
the steel-quota negotiations allocate about 15% to foreign suppliers.
---
— Expansion.Conjunction
1\. The average debt of medical school graduates who borrowed to pay for their
education jumped 10% to $42,374 this year from $38,489 in 1988, says the
Association of American Medical Colleges. (furthermore) that’s 115% more than
in 1981
— Expansion.Conjunction
2\. … he rigged up an alarm system, including a portable beeper, to alert him
when Sventek came on the line. (and) Some nights he slept under his desk.
— Expansion.Conjunction
Prices for capital equipment rose a hefty 1.1% in September, while prices for
home electronic equipment fell 1.1%. (Meanwhile,) food prices declined 0.6%,
after climbing 0.3% in August.
— Comparison.Contrast
1\. Lloyd’s overblown bureaucracy also hampers efforts to update marketing
strategies. (Although) some underwriters have been pressing for years to tap
the low-margin business by selling some policies directly to consumers.
— Comparison.Contrast 2\. Valley National ”isn’t out of the woods yet.
(Specifically), the key will be whether Arizona real estate turns around or at
least stabilizes
— Expansion.Restatement
Table 4: Example of attention in Context Knowledge Memory. The sentences in italic are from PDTB test set and following 2 instances are the ones with top 2 attention weights from training set. Relation | Train | Dev | Test
---|---|---|---
Comparison | 1855 | 189 | 145
Contingency | 3235 | 281 | 273
Expansion | 6673 | 638 | 538
Temporal | 582 | 48 | 55
Total | 12345 | 1156 | 1011
Table 5: Distribution of top-level implicit discourse relations in the PDTB. Methods | Four-ways |  One-Versus-all Binary ($F_{1}$)
---|---|---
$F_{1}$ | Acc. | Comp. | Cont. | Expa. | Temp.
Rutherford and Xue (2014) | 38.40 | 55.50 | 39.70 | 54.42 | 70.23 | 28.69
Qin et al. (2016b) | - | - | 41.55 | 57.32 | 71.50 | 35.43
Liu et al. (2016) | 44.98 | 57.27 | 37.91 | 55.88 | 69.97 | 37.17
Ji et al. (2016) | 42.30 | 59.50 | - | - | - | -
Liu and Li (2016) | 46.29 | 57.17 | 36.70 | 54.48 | 70.43 | 38.84
Qin et al. (2017) | - | - | 40.87 | 54.46 | 72.38 | 36.20
Lan et al. (2017) | 47.80 | 57.39 | 40.73 | 58.96 | 72.47 | 38.50
Our method | 46.40 | 61.42 | 41.83 | 62.07 | 69.58 | 35.72
Table 6: Comparison of $F_{1}$ scores (%) and Accuracy (%) with the State-of-
the-art Approaches for four-ways and one-versus-all binary classification on
PDTB. Comp., Cont., Expa. and Temp. stand for Comparison, Contingency,
Expansion and Temporal respectively.

####  4.2.1 Top-level Binary and 4-way Classification

A lot of the recent works in PDTB relation recognition have focused on first
level relations, both on binary and 4-ways classification. We also report the
performance on level-one relation classification for more comparison to prior
works. As described above, we followed the conventional experimental settings
(Rutherford and Xue, 2015; Liu and Li, 2016) as closely as possible. Table 5
shows the distribution of top-level implicit discourse relation in PDTB, it’s
worth noticing that there are only 55 instances for Temporal Relation in the
test set.

To make the results comparable with previous work, we report the $F_{1}$ score
for four binary classifications and both $F_{1}$ and Accuracy for 4-way
classification, which can be found in Table 6. We can see that our method
outperforms all alternatives on Comparison and Contingency, and obtain
comparable scores with the state-of-the-art in others. For 4-way
classification, we got the best accuracy and second-best $F_{1}$ with around
2% better than in Ji et al. (2016).

##  5 Conclusion and Future Work

We present in this paper a novel neural method trying to integrate implicit
connectives into the representation of implicit discourse relations with a
joint learning framework of sequence-to-sequence network. We conduct
experiments with different settings on PDTB benchmark, the results show that
our proposed method can achieve state-of-the-art performance on recognizing
the implicit discourse relations and the improvements are not only brought by
the increasing number of parameters. The model also has great potential
abilities in implicit connective prediction in the future.

Our proposed method shares similar spirit with previous work in Zhou et al.
(2010), who also tried to leverage implicit connectives to help extract
discriminative features from implicit discourse instances. Comparing with the
adversarial method proposed by Qin et al. (2017), our proposed model more
closely mimics humans’ annotation process of implicit discourse relations and
is trained to directly explicitate the implicit relations before
classification. With the representation of the original implicit sentence and
the explicitated one from decoder, and the help of the explicit knowledge
vector from memory network, the implicit relation could be classified with
higher accuracy.

Although our method has not been trained as a generative model in our
experiments, we can see potential for applying it to generative tasks. With
more annotated data, minor modification and fine-tuned training, we believe
our proposed method could also be applied to tasks like implicit discourse
connective prediction, or argument generation in the future.

##  6 Acknowledgments

This work was supported by German Research Foundation (DFG) as part of SFB
1102 “Information Density and Linguistic Encoding”. We would like to thank the
anonymous reviewers for their careful reading and insightful comments.

## References

  * Bengio et al. (2015) Bengio, S., O. Vinyals, N. Jaitly, and N. Shazeer (2015).  Scheduled sampling for sequence prediction with recurrent neural networks.  In Proceeding of NIPS, pp. 1171–1179.
  * Cai and Zhao (2017) Cai, D. and H. Zhao (2017).  Pair-aware neural sentence modeling for implicit discourse relation classification.  In IEA-AIE, pp. 458–466. Springer.
  * Dai and Le (2015) Dai, A. M. and Q. V. Le (2015).  Semi-supervised sequence learning.  In Proceedings of NIPS, pp. 3079–3087.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S. and J. Schmidhuber (1997).  Long short-term memory.  Neural computation 9(8), 1735–1780.
  * Ji and Eisenstein (2015) Ji, Y. and J. Eisenstein (2015).  One vector is not enough: Entity-augmented distributional semantics for discourse relations.  TACL 3, 329–344.
  * Ji et al. (2016) Ji, Y., G. Haffari, and J. Eisenstein (2016).  A latent variable recurrent neural network for discourse relation language models.  In Proceedings of NAACL, pp. 332–342.
  * Kalchbrenner et al. (2014) Kalchbrenner, N., E. Grefenstette, and P. Blunsom (2014).  A convolutional neural network for modelling sentences.  In Proceedings of ACL.
  * Kingma and Ba (2014) Kingma, D. and J. Ba (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Lan et al. (2017) Lan, M., J. Wang, Y. Wu, Z.-Y. Niu, and H. Wang (2017).  Multi-task attention-based neural networks for implicit discourse relationship representation and identification.  In Proceedings of EMNLP, pp. 1299–1308.
  * Lin et al. (2009) Lin, Z., M.-Y. Kan, and H. T. Ng (2009).  Recognizing implicit discourse relations in the penn discourse treebank.  In Proceedings of EMNLP, pp. 343–351.
  * Liu et al. (2018) Liu, Q., Y. Zhang, and J. Liu (2018).  Learning domain representation for multi-domain sentiment classification.  In Proceedings of NAACL, pp. 541–550.
  * Liu and Li (2016) Liu, Y. and S. Li (2016).  Recognizing implicit discourse relations via repeated reading: Neural networks with multi-level attention.  In Proceedings of EMNLP, pp. 1224–1233.
  * Liu et al. (2016) Liu, Y., S. Li, X. Zhang, and Z. Sui (2016).  Implicit discourse relation classification via multi-task neural networks.  In Proceedings of AAAI, pp. 2750–2756.
  * Luong et al. (2015) Luong, M.-T., H. Pham, and C. D. Manning (2015).  Effective approaches to attention-based neural machine translation.  In Proceedings of EMNLP, pp. 1412–1421.
  * Mikolov et al. (2013) Mikolov, T., I. Sutskever, K. Chen, G. S. Corrado, and J. Dean (2013).  Distributed representations of words and phrases and their compositionality.  In Proceedings of NIPS, pp. 3111–3119.
  * Miltsakaki et al. (2005) Miltsakaki, E., N. Dinesh, R. Prasad, A. Joshi, and B. Webber (2005).  Experiments on sense annotations and sense disambiguation of discourse connectives.  In Proceedings of the Fourth Workshop TLT-2005.
  * Prasad et al. (2008) Prasad, R., N. Dinesh, A. Lee, E. Miltsakaki, L. Robaldo, A. Joshi, and B. Webber (2008).  The penn discourse treebank 2.0.  In Proceedings of LREC.
  * Qin et al. (2016a) Qin, L., Z. Zhang, and H. Zhao (2016a).  Implicit discourse relation recognition with context-aware character-enhanced embeddings.  In Proceedings of COLING.
  * Qin et al. (2016b) Qin, L., Z. Zhang, and H. Zhao (2016b).  A stacking gated neural architecture for implicit discourse relation classification.  In Proceedings of EMNLP, pp. 2263–2270.
  * Qin et al. (2017) Qin, L., Z. Zhang, H. Zhao, Z. Hu, and E. Xing (2017).  Adversarial connective-exploiting networks for implicit discourse relation classification.  In Proceedings of ACL, pp. 1006–1017.
  * Reimers and Gurevych (2018) Reimers, N. and I. Gurevych (2018).  Why comparing single performance scores does not allow to draw conclusions about machine learning approaches.  arXiv preprint arXiv:1803.09578.
  * Rutherford et al. (2017) Rutherford, A., V. Demberg, and N. Xue (2017).  A systematic study of neural discourse models for implicit discourse relation.  In Proceedings of EACL, pp. 281–291.
  * Rutherford and Xue (2014) Rutherford, A. and N. Xue (2014).  Discovering implicit discourse relations through brown cluster pair representation and coreference patterns.  In Proceedings of EACL, pp. 645–654.
  * Rutherford and Xue (2015) Rutherford, A. and N. Xue (2015).  Improving the inference of implicit discourse relations via classifying explicit discourse connectives.  In Proceedings of NAACL, pp. 799–808.
  * Shi and Demberg (2017) Shi, W. and V. Demberg (2017).  On the need of cross validation for discourse relation classification.  In Proceedings of EACL, pp. 150–156.
  * Shi et al. (2018) Shi, W., F. Yung, and V. Demberg (2018).  Acquiring annotated data with cross-lingual explicitation for implicit discourse relation classification.  arXiv preprint arXiv:1808.10290.
  * Shi et al. (2017) Shi, W., F. Yung, R. Rubino, and V. Demberg (2017).  Using explicit discourse connectives in translation for implicit discourse relation classification.  In Proceedings of IJCNLP, pp. 484–495.
  * Sutskever et al. (2014) Sutskever, I., O. Vinyals, and Q. V. Le (2014).  Sequence to sequence learning with neural networks.  In Proceedings of NIPS, pp. 3104–3112.
  * Wu et al. (2016) Wu, C., X. Shi, Y. Chen, Y. Huang, and J. Su (2016).  Bilingually-constrained synthetic data for implicit discourse relation recognition.  In Proceedings of EMNLP, pp. 2306–2312.
  * Xue et al. (2015) Xue, N., H. T. Ng, S. Pradhan, R. Prasad, C. Bryant, and A. Rutherford (2015).  The conll-2015 shared task on shallow discourse parsing.  In Proceedings of CoNLL-15 Shared Task, pp. 1–16.
  * Xue et al. (2016) Xue, N., H. T. Ng, A. Rutherford, B. Webber, C. Wang, and H. Wang (2016).  Conll 2016 shared task on multilingual shallow discourse parsing.  In Proceedings of CoNLL-16 shared task, pp. 1–19.
  * Zhang et al. (2015) Zhang, B., J. Su, D. Xiong, Y. Lu, H. Duan, and J. Yao (2015).  Shallow convolutional neural network for implicit discourse relation recognition.  In Proceedings of EMNLP, pp. 2230–2235.
  * Zhou et al. (2016) Zhou, P., W. Shi, J. Tian, Z. Qi, B. Li, H. Hao, and B. Xu (2016).  Attention-based bidirectional long short-term memory networks for relation classification.  In Proceedings of ACL, pp. 207–212.
  * Zhou et al. (2010) Zhou, Z.-M., Y. Xu, Z.-Y. Niu, M. Lan, J. Su, and C. L. Tan (2010).  Predicting discourse connectives for implicit discourse relation recognition.  In Proceedings of COLING, pp. 1507–1514.