# CASTELO: Clustered Atom Subtypes aidEd Lead Optimization – a combined
machine learning and molecular modeling method

Leili Zhang1∗†, Giacomo Domeniconi2∗†, Chih-Chieh Yang1, Seung-gu Kang1,
Ruhong Zhou1, Guojing Cong1

1IBM Thomas J Watson Research Center,
1101 Kitchawan Rd, Yorktown Heights, NY 10598
2IBM Research Zürich,
Säumerstrasse 4, 8803 Rüschlikon, Switzerland

†Contributed equally to this work;
∗To whom correspondence should be addressed;
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Drug discovery is a multi-stage process that comprises two costly major steps:
pre-clinical research and clinical trials. Among its stages, lead optimization
easily consumes more than half of the pre-clinical budget. We propose a
combined machine learning and molecular modeling approach that automates lead
optimization workflow in silico. The initial data collection is achieved with
physics-based molecular dynamics (MD) simulation. Contact matrices are
calculated as the preliminary features extracted from the simulations. To take
advantage of the temporal information from the simulations, we enhanced
contact matrices data with temporal dynamism representation, which are then
modeled with unsupervised convolutional variational autoencoder (CVAE).
Finally, conventional clustering method and CVAE-based clustering method are
compared with metrics to rank the submolecular structures and propose
potential candidates for lead optimization. With no need for extensive
structure-activity relationship database, our method provides new hints for
drug modification hotspots which can be used to improve drug efficacy. Our
workflow can potentially reduce the lead optimization turnaround time from
months/years to days compared with the conventional labor-intensive process
and thus can potentially become a valuable tool for medical researchers.

_Keywords_ Lead optimization $\cdot$ Drug discovery $\cdot$ Molecular dynamics
simulation $\cdot$ Machine learning $\cdot$ Variational autoencoder

## Introduction

At a time of global health crisis, drug discovery is of utter importance to
bring the society back to its order. Beyond the crisis, drugs have helped to
improve life quality and increase the life expectancy. However, despite the
growing research and development expenditure every year [20, 40], the yearly
FDA-approval of drugs has mostly stalled since 1993 [29]. In fact, there were
a total of 3,437 FDA approved small-molecule and large-molecule drugs or
therapeutics in 2018 [9], with a yearly addition of only $\sim$1.2% (2014-2018
average). The so-called “Eroom’s law” overshadows the pharmaceutical industry,
with new molecular entity (NME) drugs per billion dollars spent decreasing
ever since 1950 [18]. While millions of people globally are suffering from
incurable diseases such as Parkinson’s disease [6], Huntington’s disease [21],
HIV/AIDS [11] and hepatitis B [26], technological revolutions are required to
find the cures and possibly avert Eroom’s law.

Computer-aided drug discovery (CADD) methods have been widely used in the
pharmaceutical industry before preclinical trials, from quantitative
structure-activity relationships (QSAR, [10]), pharmacophore modeling [43] to
drug-target docking algorithms [23]. Along with the conventional high-
throughput screening (HTS [27]) method and rising fragment-based drug
discovery (FBDD [30]) method, computational methods have helped to lower the
initial cost of drug discovery (target-to-hit and hit-to-lead, specifically).
The resulting lead products usually have a dissociation constant
($K_{\text{d}}$) in the micromolar range. The subsequent process would aim to
optimize the drug lead by lowering the $K_{\text{d}}$ value to the nanomolar
range while ensuring the safety of the drug (and therefore termed “lead
optimization”). However, the optimization processes highly rely on domain
knowledge and luck because of the astronomical possibilities in the chemical
space. Several computational methods have been developed to assist lead
optimization. For example, a semi-automated optimization method was suggested
by Lewis, combining 3D-QSAR model, random structure permutation methods and
human decisions [24]. Jain applied a combined QSAR-docking method to an
existing library of CDK2 inhibitors and showed encouraging agreement with the
experiments [17]. Tang et al. improved zinc endopeptidase inhibitor by
intuitively suggesting drug modifications from a binding mode suggested by the
molecular dynamics (MD) simulation [37]. Jorgensen et al. adopted Monte Carlo
(MC) simulations and free energy perturbation (FEP) simulations to virtually
screen chemical groups on a drug lead [19]. However, despite these
developments, lead optimization remains the most costly step before the
clinical trials, constituting $\sim$74% of the pre-clinical costs [31]. A
common issue of the previous computational methods lies in the generation of
new structures, which mostly relies on inefficient random guesses. Therefore,
a systematic way of identifying improvable molecular fragments is needed to
assist lead optimization decisions. Notably, computer-aided FBDD algorithms
[4] can suggest leads by combining the fragment hits (fragments with
relatively strong affinity, typically with $K_{\text{d}}$ $<$ 10 mM [34]). But
the use is limited once the lead has been identified.

Figure 1: The general pipeline for CASTELO. The starting point is the
generation of MD trajectories, with tools such as GROMACS. RMSD clustering can
be done with VMD software. In another route, we process MD trajectories with
python scripts to obtain contact matrices. Atom subtype information is used to
aggregate the calculated contact matrices. Following that, dynamism tensors
with temporal information is generated on top of the contact matrices using
python scripts. CVAE model is used to encode the dynamism data, before
clusters are calculated with tools such as HDBSCAN. Finally, we converge the
two routes by comparing clusters from conventional RMSD clustering and CVAE
clustering with proposed comparison metrics. The atom subtypes are ranked, as
the final output of CASTELO. With domain knowledge, we suggest modifications
for the lowest ranked atoms. Methods such as free energy perturbation
calculations can be used to verify CASTELO’s suggestions.

In this study, we propose a computational method, coined Clustered Atom
Subtypes aidEd Lead Optimization (CASTELO), that identifies modifiable
submolecular moieties in a lead molecule to narrow down the random guesses to
a few possibilities. The process is depicted in Figure 1. Briefly, we first
obtain the target-lead binding complexes from either crystal structures or
reasonable computational methods (such as homology modeling and MD
simulations). Using this structure, $\sim$100 ns of MD simulations are then
conducted. Subsequently, contact matrices, containing the relative distances
between atoms in the lead and atoms in the target protein, are extracted from
the simulations. The contact matrices are processed to contain temporal
information (coined as "dynamism tensor") and atom subtype information (see
details below). Compressed vectors are then created with convolutional
variational autoencoder (CVAE) using the processed contact matrices. CVAE
automatically reduces the high dimensionality of the contact matrix into a
latent space where states that share similar structural and energetic
characteristics [3] to each other. In order to consider the dynamic behavior
of the target and the lead, each time step is modeled with a dynamism tensor
(a 2-dimensional contact matrix) that contains target-lead interacting
information of the MD simulation. This tensor, closely resembling as a
2-channel image, is then fed to the CVAE. The latent space representation
generated by the CVAE allows the clustering of the time steps of the MD
simulation, grouping time steps with similar behavior. Stable snapshots will
be clustered in big clusters, while unstable ones will generate a number of
small clusters.

The goal of our method is to discover specific submolecular moieties of the
lead molecule that harm the target-lead interactions (thus termed as
"malicious atoms"). One way to group the atoms in a molecule to "submolecular
moieties" is to categorize them in subtypes based on their physical properties
(see Method for details). To differentiate the contributions of the subtypes,
CVAE and clustering are not only applied on the whole dynamism tensors, but
also to the tensors of each subtype. Clustering information are thus generated
for each subtype and compared with conventional clustering information (such
as from root mean sqaure displacement (RMSD) clustering) of the whole
molecule. Finally, using comparison metrics such as cosine similarity or
average difference (see details below), we rank the subtypes from malicious
atoms to beneficial atoms. If the overall simulation renders a stable binding
structure for the lead (which is often the case), the atom subtypes with lower
values of comparison metrics are labeled as malicious atoms, i.e. modifiable
atoms for lead optimization, because of their “deviation” from the stable
binding state of the rest of the molecule. To verify if these atoms are indeed
improvable, we modify the suggested atoms and calculate binding free energy
change using FEP calculations.

As an example, we applied this process to the sweetness taste receptor T1R2
with five well-known sweeteners: sucrose (the reference for relative
sweetness, RS = 1 [1]),
(1R,2R,3R,4R,5R)-4-Chloro-1-[(2R,3S,4S,5R)-3,4-dihydroxy-2,5-bis(hydroxymethyl)oxolan-2-yl]oxy-6-(hydroxymethyl)oxane-3,4-diol
(4R-Cl-sucrose thereafter, RS = 5 [14]), sucralose (RS = 600 [1]), dulcin (RS
= 250 [1]), and an isovanillyl sweetener (isovanillyl thereafter, RS = 400
[2]). The sweetness taste receptor structure is taken from Perez-Aguilar et
al. [32], shown in Figure 2A. The sweetener molecules are shown in Figure 2C.
The search of the binding modes for all five sweeteners follows the the
combined docking-MD method reported previously [46]. With CASTELO, we are able
to suggest modifications that improve molecular sweetness of the simulated
molecules by identifying malicious atoms. Importantly, no structure-activity
relationship data are needed for this process. The minimal requirement is a
working lead molecule and a target protein structure.

Figure 2: The structure of sucrose bound T1R2 is shown in (A). The detailed
interactions between sucrose and its surrounding residues are depicted in (B).
The structures of the five sweeteners tested in this study are illustrated in
(C).

## Methods

### Molecular docking

The protein structure of T1R2 was taken from Perez-Aguilar et al. [32]. The
structures of sucrose, 4R-Cl-sucrose, sucralose, dulcin and isovanillyl were
constructed with Jmol. [12] We refered to previous studies [44] and assumed
that the "flytrap" domain (referring to the boxed region in Figure 2A due to
the resemblance between T1R2 and the venus flytrap) of T1R2 should be the
binding domain for sweeteners. Initial binding structures were then
constructed by selecting the "flytrap" domain in T1R2 using Autodock Vina.
[38]

### Molecular dynamics simulation (dataset)

The initial binding structures constructed from Autodock Vina were solvated
with water in a simulation box of 10${\times}$10${\times}$10 nm3. 0.1 M of
NaCl was then used to ionize and neutralize the simulation box. The resulting
simulation systems contain roughly 105,000 atoms. CHARMM36 force field [15]
was used for the T1R2 protein, water and ions. CGenFF force field was
generated with the online server [41] for the five sweeteners. NAMD2.9 [33]
was used to run all simulations on Blue Gene supercomputer. Particle mesh
Ewald (PME) was used for electrostatic interaction calculations with a grid
size of 1 Å. A switching function was used for van der Waals (VDW)
interactions calculations, where switching distance was 10 Å, cutoff distance
was 12 Å and VDW pairlist distance was 13.5 Å. Before production runs, 20,000
steps of energy minimization and 250 ps of equilibration using 0.5 fs time
step were performed. For the production runs we used 2 fs as the time step.
Depending on the stability of the structures, 50 to 100 ns of simulations were
performed for each of the docked structures. Snapshots of the trajectories
were saved every 20 ps for the data collection.

### Free energy perturbation (FEP)

Following previous studies [46, 45], we used Zwanzig equation [47] and
thermodynamic cycle shown in Figure S1 to calculate the relative binding free
energy of the five sweeteners with "muted ethane" (refered to as ethane
thereafter for simplicity) as the reference point. Muted ethane had the
structure of ethane but was chargeless. We performed two sets of simulations
with NAMD2.9 for each sweeteners (the thermodynamic cycle is shown in Figure
S1): the "bind" simulations where T1R2 was present and sweeteners were mutated
to ethane; and "free" simulations where T1R2 was absent and sweeteners were
mutated to ethane. The resulting free energy change was denoted as
${\Delta}F_{bind}$ and ${\Delta}F_{free}$. The relative binding free energy of
the sweeteners can then be calculated as:
${\Delta\Delta}F_{SWT}={\Delta}F_{SWT}-{\Delta}F_{ethane}={\Delta}F_{free}-{\Delta}F_{bind}$.
If we assume that the relative sweetness of the sweeteners was a direct
measurement of the dissociation constant $K_{d}$, then we can compute the
relative sweetness as follows:

$CRS(SWT)=e^{-({\Delta\Delta}F_{SWT}-{\Delta\Delta}F_{sucrose})/RT}$ (1)

### Grouping submolecular moieties using atoms subtypes

We categorized atoms using the VDW parameters of the atom types (such as
CG331) from the CHARMM36 force field. Atom types that fell in the range of 10%
of $\sigma$ (particle size) and $\epsilon$ (dispersion energy) with each other
were classified into the same "atom subtype". Such classification resulted in
46 atom subtypes from the approximate 300 atom types from the CHARMM36
protein, carbohydrate and cgenff force field. The 46 atom subtypes were used
to differentiate "beneficial atoms" (atom subtypes that likely strengthened
sweetener-T1R2 interactions) from "malicious atoms" (atom subtypes that likely
weakened sweetener-T1R2 interactions). Note that we chose these the atom
subtypes based on our experience with the current dataset. Other subtyping
method was not excluded in our patent application and future studies.

### AI-enhanced Clustering

Figure 3: CVAE pipeline. Left part of the figure shows the input preparation:
the contact matrix of a time step $t$ is grouped over the $y$-axis by drug
atom subtypes (yellow cells are contacts). Each contact matrix of a time step
$t$ is paired to its temporal difference with a previous one $t-\delta$ to
form a dynamism tensor (violet/yellow cells are new/disappearing contacts,
green cells are stable). Each atom subtype is selected in turn and fed to the
CVAE. The representation in the latent space is used to cluster the time steps
by similar behaviors.

Autoencoders are unsupervised neural networks widely used for dimensionality
reduction and pattern recognition. Their architecture is composed of an
encoder part $g_{\phi}$ that compresses an input into a latent space and a
decoder part $f_{\theta}$ that reconstructs the original input using the low-
dimensional features. Variational autoencoders (VAEs) introduce the
optimization constraint to the latent space to be normally distributed, this
coerces the network to distribute the information more evenly into the latent
space [8]. In MD simulation analysis, a common input is a contact matrix from
a single MD time step; hence, convolutional layers (CVAE) rather than regular
feedforward are typically used. This leads to filter maps that can better
recognize local patterns independently of the position.

We used CVAE to model each MD trajectory in a low-dimensional space. Similar
to previous works [39, 3], we model the input as a contact matrix between
atoms of the drug molecule and atoms of the protein in each structure of the
simulation. Specifically, the contact matrix has size of $N\times M$, where N
is the number of drug atoms and M is the number of protein atoms. Our goal is
to find binding patterns during the simulation, hence the temporal information
over time is important to locate stable and unstable states. We enriched each
time step $t$ input with a contacts’ dynamism representation: the difference
of contacts between the current time $t$ and a previous one $t-\delta$, where
$\delta$ is a parameter of our framework. To maintain the the dynamism matrix
binary as well as independent to the direction of the movement, we use the
absolute value of the difference. In this matrix, 0 cells represents a stable
situation, while 1s show a dynamic behavior (either a new contact or the
disappearance of a previous contact). Each time step of the simulation is then
modeled as a tensor of $2\times N\times M$ (termed as the dynamism tensor),
with the first two dimensions representing the contact matrix in $t$ and the
temporal difference of the contacts in $t$ and $t-\delta$. We used
$\delta=500$ (i.e. 10 ns) in our experimental setting.

CVAE represents each contact tensor in a $d$-dimensional latent space, i.e. a
vector of $d$ elements. This intrinsically brings time steps with similar
contact matrices to closer positions in the latent space, allowing clustering
techniques to group sets of time steps with similar behavioral patterns. In
other words, applying clustering over the latent space vectors allows the
discrimination of big clusters with several time steps of similar behavior,
i.e. stable situation, or small clusters, i.e. unstable. In CASTELO, we used
HDBSCAN [7] as our clustering method. This algorithm overcomes the limitation
of knowing in advance the number of clusters, that is the typical drawback of
partitioning techniques such as K-means, and also the density threshold,
typically required by the standard DBSCAN.

The procedure described above provides general binding information of a MD
simulation, i.e. states of the whole drug molecule. CASTELO provides a fine-
grained perspective of the drug molecule behavior: a specific binding view of
each atom, or group of atoms. At this purpose, we group the drug atoms
according to their physical properties (as described in previous Section) and
we apply the CVAE and clustering pipeline described above only focusing on the
contacts of each individual atom subtype. In this case, each time step is
sampled by a tensor of size $2\times 1\times M$, and squeezing the 1-dimension
becomes a 2-channel image of size $2\times M$. Figure 3 shows the pipeline
with the input preparation of one targeted atom subtype. We repeat the CVAE
and clustering procedure for each atom subtype of the drug. This approach
allows the comparison of the overall binding behavior with that of each type
(detailed in next section).

CASTELO pipeline is completely formed by unsupervised techniques, that are
hard to tune and evaluate. In order to provide reliable suggestions, we train
slightly different CVAE architectures varying the hyper-parameters. In our
experiments we varied the latent dimensions $d$ and the number of
convolutional filters $f$. This choice was made because we noticed few cases
of all non-clustered samples in cases with latent space dimension $d$ too big.
We also tested the results by varying the time $\delta$ of the contacts’
dynamism, but the contribution was not pronounced.

### Experimental settings

As descrived above, we trained a number of CVAE models for each atom subtype.
Each input sample, i.e. a time step $t$, is represented by a $2\times M$
tensor selecting the contact and dynamism vectors related to the targeted
subtype. In our experimental settings the CVAE encoder part is formed by four
convolutional layers with $f$ filters of size 1x7, we decided to keep the
convolutions between contacts and dynamism separated, thus the 1-sized filter
in the first dimension. We trained models varying the filter numbers
$f\in\\{32,64\\}$. Considering the highly narrowed shape of the input, we used
a stride of 2 in the second dimension, i.e. on the protein atoms. We did not
used pooling or padding. The decoder part is specular to the encoder. We used
three different values for the latent space size: $d\in\\{3,5,10\\}$. For each
atom subtype we thus train $A=6$ different CVAE architectures, varying $f$ and
$d$. We combined BCE and KLE as loss function [22], using RMSProp as optimizer
and a learning rate of $0.005$. We trained each CVAE model for a maximum of
600 epochs with an early stopping mechanism after 10 epochs without
improvement in the loss. We set HDBSCAN with a minimum cluster size of 50, all
the other parameters are set to the default values.

### Submolecular moiety suggestion for lead optimization

The suggestion of malicious atoms to increase the binding affinity of the
molecule is made by a comparison of the clustering result of each atom subtype
with respect to the clusters obtained while considering the whole molecule. In
this comparison, a clustering that considers the whole molecule is needed as
reference. A CVAE model that has the contact matrices with all the drug atoms
as input may be used as reference. Otherwise, in order to avoid possible
similarity biases due to the same architecture used by the compared clusters,
any other traditional clustering methodology can be used as reference. In our
experiments, we decided to use the default RMSD clustering provided by
VMD1.9.3. [16]

As described in the previous section, each atom subtype has a number of
clustering results, one for each trained CVAE architecture. The final
comparison is made with the averaged values of the comparison of each
architecture with the whole molecule. This comparison also provides a standard
deviation value that may be used as agreement score of the different CVAE
models of the same atom subtype. I.e. averaged comparison values with small
standard deviation indicate a well grounded suggestion.

More specifically, two types of comparison metrics were explored in this
study, with the first one being the cosine similarity (CosSim):

$CosSim=\frac{\sum\limits_{t\in T}C_{t,A}\cdot
C_{t,S}}{\sqrt{\sum\limits_{t\in T}C_{t,A}^{2}}\cdot\sqrt{\sum\limits_{t\in
T}C_{t,S}^{2}}}$ (2)

where $t$ is a time frame that belongs to the trajectory ($T$), $C_{t,A}$ is
the cluster size of the cluster at time $t$ for atom subtype A, calculated
with the AI-enhanced clusering, $C_{t,S}$ is the cluster size of the cluster
at time $t$ for the whole molecule (S), calculated with conventional methods
such as RMSD clustering. Note that CosSim has a range of $[0,1]$, with lower
values indicating deviation between the two arrays and higher values
indicating similarity.

The second comparison metric examined was average difference (AvgDiff):

$AvgDiff=\frac{\sum\limits_{t\in T}(C_{t,A}-C_{t,S})}{\text{count}(T)}$ (3)

where $\text{count}(T)$ is the total number of time frames for the trajectory.
The values of AvgDiff have a wider range than CosSim, depending on the cluster
sizes. Generally speaking, a negative AvgDiff indicates that subtype A is less
stable than the whole molecule (S), while a positive AvgDiff indicates the
opposite.

After assigning a comparison metric value, each of the atom subtypes is given
a final ranking, with the lowest ranked subtype as the suggestion for lead
optimization. As a follow-up in this study, we used domain knowledge to
suggest some modifications and tested the modifications with FEP calculations
(see results for more details).

## Results and discussion

Sweetener | log(RS)* |  ${\Delta\Delta}F$** | log(CRS)***
---|---|---|---
Sucrose | 0[1] | -6.9 ${\pm}$ 0.7 | 0 ${\pm}$ 0
4R-Cl-sucrose | 0.70[14] | -10.2 ${\pm}$ 0.8 | 2.33 ${\pm}$ 0.76
Sucralose | 2.78[1] | -11.7 ${\pm}$ 0.8 | 3.38 ${\pm}$ 0.76
Dulcin | 2.40[1] | -10.6 ${\pm}$ 0.4 | 2.61 ${\pm}$ 0.59
Isovanillyl | 2.60[2] | -11.1 ${\pm}$ 0.8 | 2.96 ${\pm}$ 0.74
Table 1: Free energy. *Relative sweetness. **Binding free energy relative to
ethane (kcal/mol). ***Computed relative sweetness is calculated with Equation
1 in reference to sucrose.

### RMSD clustering identified stable binding mode and binding motifs

MD simulations were used to identify stable binding modes of the sweeteners in
the binding pocket of T1R2 flytrap domain (shortened for T1R2 thereafter) of
human sweetness taste receptor [44]. We prepared the initial guesses of the
binding structures using molecular docking software Autodock Vina [38]. Around
five highest scored docking result(s) were taken as the starting structures in
molecular dynamics simulations. Upon the completion of 50-100 ns of
simulations, binding mode clusters were classified using RMSD clustering
method provided in VMD [16] using 2.0 Å as the cutoff. RMSD values were
calculated on the sweetener molecules only, while T1R2 structures were
aligned. The identified clusters (binding modes) were considered as "stable"
if the cluster size persisted longer than 50 ns, following the practices in
previous studies [46, 45]. The persistent time of the clusters for all five
sweeteners were plotted in Figure S2 (only the trajectories with the largest
clusters were shown), indicating that stable binding modes were found for each
of the sweeteners.

In Figure 2A, we illustrated the stable binding mode extracted from the
T1R2-sucrose simulation. The sucrose molecule resided well within the binding
pocket popularly referred to as the "venus flytrap". A closer look (Figure 2B)
revealed that the crucial binding residues are mostly hydrophilic, including
the hydrogen bond receptors D142 D278, E302, among other residues such as
Y164, S165, P185, Y215, S303, A305, V309, L310, T326, I327, R383, V384 and
S387. If the two loop domains on the right of the boxed region in Figure 2A
were considered as the "lips" of the flytrap, we noticed that sucrose bound
deeply in the "mouth". One indication of the sweetener being closer to the
"lips" would be lower residue numbers among the interacting residues (see
Figure S3A for a reference of residue numbering). These findings agreed well
with the previous binding structure reported by Zhang et al. [44], where the
binding pocket facilitated hydrogen bonds between the 8 hydroxyl groups in
sucrose and hydrophilic residues in T1R2. Such hydroxyl group-facilitated
sweetness was suggested more than half a century ago and has remained a
heavily investigated topic [25, 36, 5].

The stable binding structures of the rest of the five sweeteners were
illustrated in Figure S3. We found that the crucial binding residues for
4R-Cl-sucrose included D142, E302, Y103, N143, Y164, S165, A166, T184, Y215,
S303, W304, I306, V309, L310, I325, T326, I327, R383, and S387, highly similar
to those of sucrose. This indicated that the binding mechanism mainly remained
the same after changing one hydroxyl group of the sucrose to the chlorine atom
in 4R-Cl-sucrose. In contrast, the crucial binding residues for sucralose
(D142, E302, L41, I67, L71, S144, Y164, S165, A166, I167, T184, H190, S303,
A305, V309, and V384) included more residues on the slightly more hydrophobic
lips region of the flytrap. This shift was expected because 3 relatively more
hydrophilic hydroxyl groups of the sucrose were modified to 3 chlorine atoms
in sucralose. Nevertheless, D142 and E302 remained crucial as hydrogen bond
receptors for sucralose.

The stable binding structure of dulcin interacted with mostly hydrophilic
residues L41, Y103, D142, Y164, S165, T184, P185, H190, E302, S303, A305,
I306, T326, I327, R383, V384, S387, and L448, similar to those of sucrose and
4R-Cl-sucrose. The main hydrogen bonds appeared to be between the urea group
of dulcin and D142. Finally, the stable binding structure of isovanillyl
interacted with F39, S40, L41, V64, I67, L71, Y103, D142, Y164, S165, Y215,
E302, S303, V309, V384, V385, S387, V388, indicating a binding domain that was
more hydrophobic and similar to sucralose. Only occasional (weak) hydrogen
bonds were found between isovanillyl and E302, along with backbone carbonyl
groups of several residues.

### Binding free energy of the sweeteners agreed well with experimental
relative sweetness

Starting from the stable binding structures of the five sweeteners in this
study, we performed FEP calculations for each by constructing "dual
topologies" between muted ethane and the corresponding sweetener (see Methods
for more details). The computed binding free energy ${\Delta\Delta}F$ was
therefore the relative binding free energy between the sweetener and muted
ethane in T1R2 binding pocket. We assumed that muted ethane had a flat binding
energy surface in the T1R2 binding site due to its chargeless nature. Then the
calculated binding free energy should be relatively accurate among the five
sweeteners. Note that similar approaches have been verified in a number of
previous publications [46, 45, 28, 35, 42].

We noticed that the general trend of ${\Delta\Delta}F$ values and experimental
values of $RS$ agreed well with each other (Table 1). The binding mechanisms
of the five sweeteners were, however, vastly different as suggested by the
interacting residues listed above. To bridge between our speculations on the
structural information and quantitation, we reported the electrostatic and VDW
contributions in ${\Delta\Delta}F$ by performing energy decomposition for all
five FEP calculations (see Table S1). We found that sucrose, 4R-Cl-sucrose,
and dulcin mostly utilized polar interactions. In contrast, isovanillyl mostly
utilized non-polar interactions, agreeing well with the structural
speculations. Surprisingly, despite being in a similar binding domain to
isovanillyl, sucralose still utilized mostly polar interactions. Comparing the
binding modes between sucrose, 4R-Cl-sucrose and sucralose, it was likely that
the replacement of hydroxyl groups to chlorine lowered the penalty of the
remaining hydrophilic hydroxyl groups in hydrophobic environment, rather than
directly boosting the hydrophobic interactions.

We further computed the relative sweetness ($CRS$), listed in Table 1. The
computed results agreed generally well with the experimental $RS$ values. Such
agreement showcased that FEP method could be used to calculate the $RS$ of
artificial sweeteners. In later sections, we used FEP calculations to validate
that "lead optimization" (meaning to increase the $RS$ value) of the sweetener
was possible using CASTELO.

### AI-enhanced clustering with all atoms agree well with the RMSD clustering

Figure 4: T-SNE projection of the CVAE latent representations. The top row
uses a color map based on the simulation time. The middle and bottom rows use
colors to distinguish clustering results respectively of our AI-enhanced model
and of the RMSD. Here we depict the encoded representations generated by the
CVAE model with $f=32$ filters and $d=5$ latent dimensions.

We firstly tested the reliability of our AI-enhanced clustering with a direct
comparison to a reference clustering using RMSD in VMD1.9.3. In order to get a
real comparison, we trained the CVAE models using all the sweetener atoms,
i.e. with the same contact matrices used by RMSD as input. To have a
2-dimensional tensor as input of the CVAE, we concatenated the contacts and
dynamism tensors in the first dimension, forming a $2N\times M$ input. We
noticed that the obtained clustering results agreed well with those obtained
by RMSD. Figure 4 shows the t-SNE [13] visualization of the latent
representations by the CVAE models of each sweetener, each data point is a
time step of the simulation. In the top row we used a color map based on the
simulation time. The middle and bottom rows used colors to distinguish
clusters found respectively by our AI-enhanced model and by the RMSD. Firstly,
we noted that the HDBSCAN behaved well in clustering time steps projected into
the latent space by the CVAE model, this was noticeable by the agreement of
colors mapping in the middle row and actual clusters of points. It was also
noticeable that, generally, clusters are formed by time frames close in the
simulation. This supported the use of the dynamism tensor, that can help the
model to distinguish time steps with similar contact matrices but different
dynamism. The general agreement between the middle and bottom rows proved the
reliability of our AI-enhance clustering methodology, allowing us to proceed
in a fine-grained analysis of the drug molecule behavior.

### AI-enhanced clustering varies with individual atom subtypes, hinting
malicious atom subtypes

(a)  (b)  (c)  (d)  (e)

Figure 5: Clustering results for all atom subtypes in the five sweeteners. The
first bar shows the RMSD clustering with VMD1.9.3. The second bar is
calculated with CVAE and HDBSCAN when all atoms are considered. The following
bars represent the clustering results for each atom subtype of the
corresponding sweetener. Each plot depicts the cluster size at time $t$
against simulation time $t$. For CVAE and HDSCAN results, each cluster size
$t$ is the average of the sizes given by the $A$ models. Bigger clusters are
colored with dark red while small clusters are colored with dark blue.

We further analyzed the behavior of each atom subtype and compared it to the
RMSD reference. With the first bar as reference from RMSD, Figure 5 offers a
visual comparison of the differences in clustering. In particular, each bar
was created with a color map that shows the cluster size of each time step $t$
of the simulation (with bigger clusters in dark red and smaller in dark blue).
The second bar shows the clusters obtained by CVAE and HDBSCAN when using all
the sweetener atoms. The remaining bars show the clustering results for each
individual subtype. For all the AI-enhanced bars, the averaged cluster size
among the clustering results given by the different model trained is showed.
All trajectories have a stable binding mode according to RMSD clustering, that
showed big red clusters along almost all the simulations. Only sucralose
showed a less stable behavior, especially in the final part of the trajectory.
As a confirmation of the agreement between RMSD and AI-enhanced clustering,
the first two bars of each sweetener agreed well. Only for isovanillyl this
agreement was less marked (as also showed in Figure 4). The different binding
behaviors of the individual subtypes were observable by comparing the subtypes
cluster sizes with the RMSD ones. Reddish bars identified stable subtypes, for
instance subtype 11 and 27 showed a stable result for sucrose. Inversely, bars
with a different coloring with respect to the reference, hinted less stable
subtypes that could be marked as "malicious atoms", because they might weaken
the protein-sweetener interactions. In order to quantitatively identify the
malicious atoms, we adopted comparison metrics, as reported below.

### Comparison metrics such as cosine similarity pointed out the malicious
atoms.

If the sweeteners remained stable in T1R2, the atom subtypes that deviated
from the overall stability should be the "malicious atoms". CASTELO was
designed to process MD simulations with stable binding modes to find these
malicious atoms that weakened T1R2-sweetener binding affinity. We used
Equations 2 and 3 to compare the stability between the whole molecule
(calculated with RMSD clustering) and each of the atom subtypes (calculated
with AI-enhanced clustering). The resulting CosSim and AvgDiff metrics were
plotted in Figure 6A from all simulations. It was not surprising that these
two metrics were positively correlated (R2 = 0.31).

Figure 6: Comparison metrics are calculated by comparing the clusters with
conventional RMSD clustering and AI-enhanced clusering with Equations 2 and 3.
In (A), we show that the two comparison metrics CosSim and AvgDiff are
positively correlated. In (B) and (C), we calculated the RMSF of the atom
subtypes for all T1R2-sweetener trajectories. Clearly, CosSim and AvgDiff are
not correlated with RMSF.

It is important to mention that the deviation between atom subtypes and the
whole molecule should not be comprehended as a time-averaged measurement such
as root mean square fluctuation (RMSF). If we consider an atom subtype that
only fluctuates drastically when the whole molecule is unstable, the RMSF of
this atom subtype might be high but it does not deviate from the overall
behavior of the molecule. Therefore, each of the timeframes must be compared
separately between the whole molecule and its atom subtypes. With this in
mind, Equations 2 and 3 are used as the comparison metrics. Another advantage
of AI-enhanced clusters is that they contain rotational and vibrational
information, compared to the conventional RMSF calculations that mainly
consider translational fluctuations. The additional rotational information
from AI-enhanced clusters might be particularly useful to monitor atom groups
such as -CH3 and -CH2\- that would have been easily ignored otherwise. To
clearly show that our metrics differ from conventional RMSF calculations, we
compared the CosSim and AvgDiff to RMSF values from all simulations in Figure
6B and 6C. Both plots yielded R2 = 0.0, indicating that CosSim and AvgDiff
contained independent information from RMSF.

To find the malicious atoms, we focused on the trajectories that comprised
largest clusters, such as the ones shown in Figure 5. CosSim and AvgDiff
metrics were plotted in Figure 7 for these trajectories. A ranking system was
used to pick out the atom subtype most prone to undermine the overall
T1R2-sweetener binding stability. For example, in Figure 7A, subtype 12 of
dulcin was clearly the malicious atom subtype, because it ranked the lowest in
both CosSim and AvgDiff metrics. In Figure 7B, subtype 7 and 8 were mostly
likely the malicious atom subtypes, because of their low values in both CosSim
and AvgDiff. The selected atom subtypes were chosen as the candidate for lead
optimization (see results in the next section).

Figure 7: We rank the atom subtypes from the five sweeteners with the
calculated comparison metrics CosSim and AvgDiff. Negative values mean that
the corresponding atoms are "malicious". Positive values mean that the
corresponding atoms are "beneficial". The results for dulcin and isovanillyl
are plotted separately in (A) and (B). We group sucrose, 4R-Cl sucrose and
sucralose due to their structural similarity. The CosSim and AvgDiff values of
the 8 atom subtypes shared among these three sweeteners are plotted in (C) and
(D), respectively. The relative binding free energy (E, with ethane as the
reference point) is plotted in juxtaposition to the comparison metrics.

In Figure 7C and 7D, we compared the comparison metrics of sucrose, 4R-Cl-
sucrose and sucralose, which could be considered as a lead optimization chain
due to their structural similarity. We noticed that the added chlorine atoms
(subtype 4) in 4R-Cl-sucrose and sucralose were increasingly beneficial to the
overall structural stability (both CosSim and AvgDiff metrics were
increasingly positive). The suggested malicious atom subtypes were
carbohydrate carbons that connect with hydroxyl groups (-CHOH-, subtype 10).
Upon modifying several of the hydroxyl groups to chlorine atoms, both CosSim
and AvgDiff metrics of subtype 10 were increasingly more positive, indicating
that the replacement of chlorine atoms affected the stability of its
connecting carbons. As a comparison, we plotted the relative binding free
energy of sucrose, 4R-Cl-sucrose and sucralose in Figure 7E, which clearly
correlated with CosSim and AvgDiff metrics for subtype 4 and subtype 10.
Finally, we provided a possibility to further improve sucralose by suggesting
that the malicious atoms in sucralose were probably from subtype 27, which
were the aliphatic hydrogens in the -CH2OH group.

### FEP calculations verified that malicious atom types could be modified to
strengthen T1R2-sweetener binding affinity, providing an opportunity for lead
optimization.

Based on the ranked comparison metric values in Figure 7A and 7B, we
visualized the malicious atoms of dulcin and isovanillyl in Figure 8A and 8B,
as a way to guide chemists in lead optimization. For example, subtype 12 of
dulcin could be easily changed to a carbonyl (-CO-) or a chloronated carbon
(-CHCl-). Subtype 8 of isovanillyl could be easily changed to an aliphatic
carbon (-CH2-). In addition to these malicious atoms, we also selected some
neutral atoms as control groups (like subtype 40 of dulcin and subtype 10 of
isovanillyl). To verify if the identified atoms were indeed "malicious", FEP
calculations were adopted to compare the binding free energy of the modified
lead molecules to the lead molecules (exemplified in Figure 8; a full list can
be found in Table S2 and Table S3). For example, a dual topology was made
between dulcin and modified dulcin where group X (also subtype 12, see Table
S2) was changed from CH2 to CHCl. Two sets of simulations were performed: one
in T1R2 and one solvated with 0.1 M NaCl solution. The calculated binding free
energy (${\Delta}{\Delta}F$) was thus a direct comparison between dulcin and
the modified dulcin.

Figure 8: Attempted lead optimization for dulcin (A) and isovanillyl (B). The
molecules are illustrated with two representations: solid stick and ball model
for atom identification (cyan for carbon, white for hydrogen, blue for
nitrogen, red for oxygen and green for chlorine); bubbles for CosSim values
(blue for positive/beneficial, green for neutral and red for
negative/malicious). Two examples are illustrated, one for optimization and
one as control, for both dulcin and isovanillyl. The resulting binding free
energy changes are listed next to the examples, with negative values meaning
strengthened binding affinity and positive values meaning weakened binding
affinity. A complete list of the lead optimization can be found in Table S2
and Table S3.

We tried 4 combinations for dulcin optimization varying neutral moiety R and
malicious moiety X, shown in Table S2. Modifications on moiety R always
weakened T1R2-dulcin binding affinity. However, modifications on moiety X
resulted in one favorable modifications. -2.4 $\pm$ 1.0 kcal/mol was obtained
when moiety X was changed from -CH2\- to -CHCl-. Interestingly, when moiety X
was changed to -CO-, the binding affinity was weakened, probably because its
surroundings was slightly hydrophobic (I327, seen in Figure S3E).

Similarly, we attempted 2 combinations for isovanillyl optimization, varying
neutral moiety R and malicious moiety X, shown in Table S3. Modifications on
moiety R neither strengthened nor weakened T1R2-isovanillyl binding affinity
(${\Delta}{\Delta}F$ was not significant compared to its standard deviation).
In contrast, modifications on moiety X significantly strengthened
T1R2-isovanilly binding affinity by -3.8 $\pm$ 1.0 kcal/mol. It was counter-
intuitive that changing moiety X to a hydrophobic group would significantly
benefit the binding affinity, given that its surroundings included D142 and
Y164 (Figure S3F). On the other hand, the surroundings of moiety R included
V388, F39 and L71, which would likely favor a hydrophobic modification.
Further analysis pointed out that by modifying moiety X, its nearby phenol-
hydroxyl group was freed from forming intramolecular hydrogen bond and
established a hydrogen bond with D142. This transition lowered the
electrostatic penalty of charged D142 residue around the hydrophobic
isovanillyl (evidenced in Table S1). Our findings blindly agreed with a
previous isovanillyl optimization which suggested the optimal moiety X to be
-CH2-, [2] reinforcing that CASTELO could potentially assist in lead
optimization.

## Conclusions

In this work, we introduced CASTELO, a ML-MD pipeline that processed MD
simulation data with drug targets and their known leads and suggested
modifiable submolecular moieties for lead optimization. We generated dynamism
tensors by including temporal information in the conventional contact
matrices. CVAE method was adopted to compress the dynamism tensors into latent
space before the data clustering with HDBSCAN. The lead molecule was grouped
into atom subtypes to pin down submolecular contributions to the target-lead
binding affinity. The resulting cluster information was compared to a
traditional RMSD-based clustering method for the whole molecule with the
proposed comparison metrics CosSim and AvgDiff. Finally, we ranked the
submolecular moieties to find clues for lead optimization. With T1R2-sweetener
as a model system, we proved that our pipeline nicely explained the
improvement of the sweetness from sucrose to 4R-Cl-sucrose and sucralose. Most
notably, we suggested two brand new molecules based on the CASTELO pipeline
using T1R2-dulcin and T1R2-isovanillyl simulations. With free energy
calculations, we verified that the newly improved dulcin was $\sim$57 times
sweeter than dulcin, $\sim$14,000 times sweeter than sucrose. The newly
improved isovanillyl was computed to be $\sim$600 times sweeter than
isovanillyl, 240,000 times sweeter than sucrose. Our physics-based model
should be transferable to other systems, given that newest development of the
force field could be trusted. We plan to use more target-lead systems to test
CASTELO’s scalability and interpretability. For example, a similar approach
could be adapted for major histocompatibility complex (MHC) and epitope
complexes (paper submitted). The identification of any moiety may indicate
destabilizing motions of its surrounding atoms rather than moiety itself.
Generative models could be developed on top of the malicious atom
identifications in this study. Although CASTELO’s dependency on experimental
data is low, we foresee that our tool could be further enhanced with more use
cases, especially with more experimental data to reference to.

## Supplementary material

Table S1 provides the binding free energy decomposition on the T1R2-sweetener
systems. Table S2 is provided for the in silico lead optimization of dulcin.
Table S3 is provided for the in silico lead optimization of isovanillyl.
Figure S1 is attached to show the thermodynamic cycle used for FEP
calculations. Figure S2 plots the stable clusters identified by RMSD
clustering. Figure S3 illustrates the stable binding structures for the five
sweeteners in the T1R2 flytrap domain.

## Author’s contributions

L.Z. and G.D. conceived of the presented work. L.Z. and S.K. performed the
MD/FEP simulations. L.Z. and G.D. extracted and processed contact data from
the simulations. G.D. and C.Y. trained the CVAE models. L.Z. and G.D.
conceived of the comparison metrics for the final ranking of submolecular
moieties. R.Z. and G.C. supervised the findings of this work. All authors
discussed the results and contributed to the final manuscript.

## Acknowledgements

We thank Josef Klucik and Paul Winget for valuable discussions on the subject
of the sweeteners. R.Z. and G.C. gratefully acknowledge the financial support
from the IBM Bluegene Science Program (W125859, W1464125 and W1464164),
Computing Cloud Clusters and Witherspoon supercomputer in IBM.

## References

  * [1] Jessica Ahmed, Saskia Preissner, Mathias Dunkel, Catherine L Worth, Andreas Eckert, and Robert Preissner.  Supersweet—a resource on natural and artificial sweetening agents.  Nucleic acids research, 39(suppl_1):D377–D382, 2010.
  * [2] Angela Bassoli, Lucio Merlini, and Gabriella Morini.  Isovanillyl sweeteners. from molecules to receptors.  Pure and applied chemistry, 74(7):1181–1187, 2002.
  * [3] Debsindhu Bhowmik, Shang Gao, Michael T Young, and Arvind Ramanathan.  Deep clustering of protein folding simulations.  BMC bioinformatics, 19(18):47–58, 2018.
  * [4] Yuemin Bian and Xiang-Qun Sean Xie.  Computational fragment-based drug design: current trends, strategies, and applications.  The AAPS journal, 20(3):59, 2018.
  * [5] F Bruni, C Di Mino, S Imberti, SE McLain, NH Rhys, and MA Ricci.  Hydrogen bond length as a key to understanding sweetness.  The journal of physical chemistry letters, 9(13):3667–3672, 2018\.
  * [6] Verónica Cabreira, Patrício Soares-da Silva, and João Massano.  Contemporary options for the management of motor complications in parkinson’s disease: Updated clinical review.  Drugs, 79(6):593–608, 2019.
  * [7] Ricardo JGB Campello, Davoud Moulavi, and Jörg Sander.  Density-based clustering based on hierarchical density estimates.  In Pacific-Asia conference on knowledge discovery and data mining, pages 160–172. Springer, 2013.
  * [8] Carl Doersch.  Tutorial on variational autoencoders.  arXiv, page 1606.05908, 2016.
  * [9] Tayebeh Farhadi and Seyed MohammadReza Hashemian.  Computer-aided design of amino acid-based therapeutics: A review.  Drug design, development and therapy, 12:1239, 2018.
  * [10] Paola Gramatica.  Principles of qsar models validation: internal and external.  QSAR & combinatorial science, 26(5):694–701, 2007.
  * [11] Roy M Gulick and Charles Flexner.  Long-acting hiv drugs for treatment and prevention.  Annual review of medicine, 70:137–150, 2019.
  * [12] Robert M Hanson.  Jmol–a paradigm shift in crystallographic visualization.  Journal of Applied Crystallography, 43(5):1250–1260, 2010.
  * [13] Geoffrey E Hinton and Sam Roweis.  Stochastic neighbor embedding.  Advances in neural information processing systems, 15:857–864, 2002\.
  * [14] Les Hough and Riaz Khan.  Intensification of sweetness.  Trends in Biochemical Sciences, 3(1):61–63, 1978.
  * [15] Jing Huang and Alexander D MacKerell Jr.  Charmm36 all-atom additive protein force field: Validation based on comparison to nmr data.  Journal of computational chemistry, 34(25):2135–2145, 2013.
  * [16] William Humphrey, Andrew Dalke, Klaus Schulten, et al.  Vmd: visual molecular dynamics.  Journal of molecular graphics, 14(1):33–38, 1996.
  * [17] Ajay N Jain.  Qmod: physically meaningful qsar.  Journal of computer-aided molecular design, 24(10):865–878, 2010\.
  * [18] Richard Jones and JR Wilsdon.  The biomedical bubble: Why uk research and innovation needs a greater diversity of priorities, politics, places and people.  2018\.
  * [19] William L Jorgensen, Juliana Ruiz-Caro, Julian Tirado-Rives, Aravind Basavapathruni, Karen S Anderson, and Andrew D Hamilton.  Computer-aided design of non-nucleoside inhibitors of hiv-1 reverse transcriptase.  Bioorganic & medicinal chemistry letters, 16(3):663–667, 2006\.
  * [20] Kenneth I Kaitin.  Deconstructing the drug development process: the new face of innovation.  Clinical Pharmacology & Therapeutics, 87(3):356–361, 2010.
  * [21] Hongsuk Kang, Francisco X Vázquez, Leili Zhang, Payel Das, Leticia Toledo-Sherman, Binquan Luan, Michael Levitt, and Ruhong Zhou.  Emerging $\beta$-sheet rich conformations in supercompact huntingtin exon-1 mutant structures.  Journal of the American Chemical Society, 139(26):8820–8827, 2017\.
  * [22] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [23] Douglas B Kitchen, Hélène Decornez, John R Furr, and Jürgen Bajorath.  Docking and scoring in virtual screening for drug discovery: methods and applications.  Nature reviews Drug discovery, 3(11):935–949, 2004.
  * [24] Richard A Lewis.  A general method for exploiting qsar models in lead optimization.  Journal of medicinal chemistry, 48(5):1638–1648, 2005.
  * [25] Michael G Lindley, Gordon G Birch, and Riaz Khan.  Sweetness of sucrose and xylitol. structural considerations.  Journal of the Science of Food and Agriculture, 27(2):140–144, 1976\.
  * [26] Uri Lopatin.  Drugs in the pipeline for hbv.  Clinics in liver disease, 23(3):535–555, 2019.
  * [27] Ricardo Macarron, Martyn N Banks, Dejan Bojanic, David J Burns, Dragan A Cirovic, Tina Garyantes, Darren VS Green, Robert P Hertzberg, William P Janzen, Jeff W Paslay, et al.  Impact of high-throughput screening in biomedical research.  Nature reviews Drug discovery, 10(3):188–195, 2011.
  * [28] Shuichi Miyamoto and Peter A Kollman.  Absolute and relative binding free energy calculations of the interaction of biotin and its analogs with streptavidin using molecular dynamics/free energy perturbation approaches.  Proteins: Structure, Function, and Bioinformatics, 16(3):226–245, 1993.
  * [29] Asher Mullard.  2019 fda drug approvals.  Nature reviews. Drug discovery, 19(2):79, 2020.
  * [30] Christopher W Murray and David C Rees.  The rise of fragment-based drug discovery.  Nature chemistry, 1(3):187–192, 2009.
  * [31] Steven M Paul, Daniel S Mytelka, Christopher T Dunwiddie, Charles C Persinger, Bernard H Munos, Stacy R Lindborg, and Aaron L Schacht.  How to improve r&d productivity: the pharmaceutical industry’s grand challenge.  Nature reviews Drug discovery, 9(3):203–214, 2010.
  * [32] Jose Manuel Perez-Aguilar, Seung-gu Kang, Leili Zhang, and Ruhong Zhou.  Modeling and structural characterization of the sweet taste receptor heterodimer.  ACS Chemical Neuroscience, 10(11):4579–4592, 2019.
  * [33] James C Phillips, Rosemary Braun, Wei Wang, James Gumbart, Emad Tajkhorshid, Elizabeth Villa, Christophe Chipot, Robert D Skeel, Laxmikant Kale, and Klaus Schulten.  Scalable molecular dynamics with namd.  Journal of computational chemistry, 26(16):1781–1802, 2005.
  * [34] Amanda J Price, Steven Howard, and Benjamin D Cons.  Fragment-based drug discovery and its application to challenging drug targets.  Essays in biochemistry, 61(5):475–484, 2017.
  * [35] Devleena Shivakumar, Joshua Williams, Yujie Wu, Wolfgang Damm, John Shelley, and Woody Sherman.  Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the opls force field.  Journal of chemical theory and computation, 6(5):1509–1519, 2010\.
  * [36] Tetsuo Suami and Leslie Hough.  Molecular mechanisms of sweet taste 2: Glucopyranose, fructopyranose and sucrose.  Journal of carbohydrate chemistry, 11(8):953–967, 1992.
  * [37] Jing Tang, Jewn Giew Park, Charles B Millard, James J Schmidt, and Yuan-Ping Pang.  Computer-aided lead optimization: improved small-molecule inhibitor of the zinc endopeptidase of botulinum neurotoxin serotype a.  PLoS One, 2(8):e761, 2007.
  * [38] Oleg Trott and Arthur J Olson.  Autodock vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading.  Journal of computational chemistry, 31(2):455–461, 2010.
  * [39] Yuko Tsuchiya, Kei Taneishi, and Yasushige Yonezawa.  Autoencoder-based detection of dynamic allostery triggered by ligand binding based on molecular dynamics.  Journal of chemical information and modeling, 59(9):4043–4051, 2019\.
  * [40] ÖNER Tulum.  Innovation and financialization in the us biopharmaceutical industry.  Unpublished doctoral dissertation, Faculty of Economics, University of Ljubljana, 2018.
  * [41] Kenno Vanommeslaeghe and Alexander D MacKerell Jr.  Automation of the charmm general force field (cgenff) i: bond perception and atom typing.  Journal of chemical information and modeling, 52(12):3144–3154, 2012.
  * [42] Billy J Williams-Noonan, Elizabeth Yuriev, and David K Chalmers.  Free energy methods in drug design: Prospects of “alchemical perturbation” in medicinal chemistry: miniperspective.  Journal of medicinal chemistry, 61(3):638–649, 2018.
  * [43] Sheng-Yong Yang.  Pharmacophore modeling and applications in drug discovery: challenges and recent advances.  Drug discovery today, 15(11-12):444–450, 2010.
  * [44] Feng Zhang, Boris Klebansky, Richard M Fine, Haitian Liu, Hong Xu, Guy Servant, Mark Zoller, Catherine Tachdjian, and Xiaodong Li.  Molecular mechanism of the sweet taste enhancers.  Proceedings of the National Academy of Sciences, 107(10):4752–4757, 2010.
  * [45] Leili Zhang and Ruhong Zhou.  Structural basis of potential binding mechanism of remdesivir to sars-cov-2 rna dependent rna polymerase.  The Journal of Physical Chemistry B, 2020.
  * [46] Shulin Zhuang, Leili Zhang, Tingjie Zhan, Liping Lu, Lu Zhao, Haifei Wang, Joseph A Morrone, Weiping Liu, and Ruhong Zhou.  Binding specificity determines the cytochrome p450 3a4 mediated enantioselective metabolism of metconazole.  The Journal of Physical Chemistry B, 122(3):1176–1184, 2018.
  * [47] Robert W Zwanzig.  High-temperature equation of state by a perturbation method. i. nonpolar gases.  The Journal of Chemical Physics, 22(8):1420–1426, 1954.

Supplementary Material

for

CASTELO: Clustered Atom Subtypes aidEd Lead Optimization – a combined machine
learning and molecular modeling method

Sweetener | ${\Delta\Delta}F$ | ${\Delta\Delta}F_{elec}$ | ${\Delta\Delta}F_{vdw}$
---|---|---|---
Sucrose | -6.9 ${\pm}$ 0.7 | -8.0 ${\pm}$ 0.5 | 0.2 ${\pm}$ 0.5
4R-Cl-sucrose | -10.2 ${\pm}$ 0.8 | -9.3 ${\pm}$ 0.5 | -2.3 ${\pm}$ 0.6
Sucralose | -11.7 ${\pm}$ 0.8 | -11.5 ${\pm}$ 0.6 | -1.9 ${\pm}$ 0.6
Dulcin | -10.6 ${\pm}$ 0.4 | -7.1 ${\pm}$ 0.3 | -3.2 ${\pm}$ 0.6
Isovanillyl | -11.1 ${\pm}$ 0.8 | 0.8 ${\pm}$ 0.5 | -12.8 ${\pm}$ 0.6
Table S1: Binding free energy decomposition. ${\Delta\Delta}F_{elec}$ or
${\Delta\Delta}F_{vdw}$ is calculated with Zwanzig equation when only
electrostatic energy or VDW energy is considered. Note that due to the
construction of Zwanzig equation, ${\Delta\Delta}F$ is not a simple addition
of ${\Delta\Delta}F_{elec}$ or ${\Delta\Delta}F_{vdw}$.

R X ${\Delta\Delta}F$ NH2 CH2 0 NH2 CO +2.5 ${\pm}$ 1.1 NH2 CHCl -2.4 ${\pm}$
1.0 CH3 CH2 +2.4 ${\pm}$ 1.3 NHCH2CH3 CH2 +4.9 ${\pm}$ 1.3

Table S2: Lead optimization for dulcin. Relative free energy ${\Delta\Delta}F$
is calculated using the original dulcin molecule as the reference (where R is
NH2 and X is CH2). Negative values of ${\Delta\Delta}F$ indicate that the
modifications strengthen the binding affinity. The unit of ${\Delta\Delta}F$
is kcal/mol.

R X ${\Delta\Delta}F$ CO CO 0 CO CH2 -3.8 ${\pm}$ 1.0 CH2 CO -1.0 ${\pm}$ 1.4

Table S3: Lead optimization for isovanillyl. Relative free energy
${\Delta\Delta}F$ is calculated using the original isovanillyl molecule as the
reference (where R and X are both CO). Negative values of ${\Delta\Delta}F$
indicate that the modifications strengthen the binding affinity. The unit of
${\Delta\Delta}F$ is kcal/mol. Figure S1: The thermodynamic cycle for FEP
calculations. See method for details. Figure S2: From the MD simulations, we
identify the clusters of sweetener binding states with RMSD calculations on
the sweetener alone, while the structure of T1R2 is aligned. The cutoff for
the RMSD clustering method is 2 Å. The clusters are plotted using their
"persistent time" (defined as the total time each cluster persists in the
corresponding simulation). From top to bottom, we plot the persistent time of
sucrose, 4R-Cl-sucrose, sucralose, dulcin and isovanillyl, respectively.

Figure S3: (A) We illustrate T1R2 protein structure with N-terminus colored by
red and C-terminus colored by blue. Stable binding structures of the five
sweeteners in T1R2 are identified from the MD simulations and drawn in (B)
with sucrose, (C) with 4R-Cl-sucrose, (D) with sucralose, (E) with dulcin and
(F) with isovanillyl, respectively. Protein structures are drawn with
transparent "newcartoon" representation in VMD. Sweeteners and their
interacting residues are drawn with sticks models. The colors for the sticks
models are selected as follows: cyan for carbon in the sweeteners, grey for
carbon in T1R2, white for hydrogen, red for oxygen, blue for nitrogen, and
green for chlorine.