2023

These authors contributed equally to this work.

[1]Ye Wei These authors contributed equally to this work.

These authors contributed equally to this work.

[4]Ji Wu [1]Dan Zhao [1,5]Jianyang Zeng

1]Institute for Interdisciplinary Information Sciences, Tsinghua University,
China 2]Department of Intelligent Systems, KTH Stockholm, Sweden 3]Microsoft
Research Asia 4]Department of Electrical Engineering, Tsinghua University,
China 5]Silexon AI Technology Co., Ltd, Nanjing, China

# From Static to Dynamic Structures: Improving Binding Affinity Prediction
with a Graph-Based Deep Learning Model

Yaosen Min<EMAIL_ADDRESS>Peizhuo Wang  Xiaoting Wang  Han Li
Nian Wu  Stefan Bauer  Shuxin Zheng  Yu Shi  Yingheng Wang
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [

###### Abstract

Accurate prediction of the protein-ligand binding affinities is an essential
challenge in the structure-based drug design. Despite recent advance in data-
driven methods in affinity prediction, their accuracy is still limited,
partially because they only take advantage of static crystal structures while
the actual binding affinities are generally depicted by the thermodynamic
ensembles between proteins and ligands. One effective way to approximate such
a thermodynamic ensemble is to use molecular dynamics (MD) simulation. Here,
we curated an MD dataset containing 3,218 different protein-ligand complexes,
and further developed Dynaformer, which is a graph-based deep learning model.
Dynaformer was able to accurately predict the binding affinities by learning
the geometric characteristics of the protein-ligand interactions from the MD
trajectories. In silico experiments demonstrated that our model exhibits
state-of-the-art scoring and ranking power on the CASF-2016 benchmark dataset,
outperforming the methods hitherto reported. Moreover, we performed a virtual
screening on the heat shock protein 90 (HSP90) using Dynaformer that
identified 20 candidates and further experimentally validated their binding
affinities. We demonstrated that our approach is more efficient, which can
identify 12 hit compounds (two were in the submicromolar range), including
several newly discovered scaffolds. We anticipate this new synergy between
large-scale MD datasets and deep learning models will provide a new route
toward accelerating the early drug discovery process.

###### keywords:

Computational biology, Graph transformer, Molecular dynamics, Drug discovery

##  1 Introduction

Protein-ligand binding plays a key role in a wide range of biological
processes, such as enzyme catalysis, signaling pathways, and maintenance of
cell structure rube2022prediction ; mcfedries2013methods . Understanding the
binding mechanisms of ligand molecules is essential for the design and
development of novel drugs with high affinity and selectivity swinney2011were
; erlanson2016twenty . Therefore, developing accurate computational methods
for modeling protein-ligand binding and predicting affinities is one of the
most fundamental tasks in computational chemistry and drug discovery fields.

Recent advance in computational methods has enabled protein-ligand binding
affinity calculation through physics- and knowledge-based approaches
mobley2017predicting ; kumar2021binding . The molecular docking algorithms
utilize coarse-grained physical models to predict and evaluate the binding
conformation of small molecules. However, such methods usually incorporate
expert knowledge, hand-tuned features and parameters, thus often leading to
biased results fischer2021decision . All-atom molecular dynamics (MD) is
another popular computational technique to investigate the structural and
dynamical properties of biological systems. It is mainly based on Newtonian
mechanics to reveal the full atomic detail at a fine spatial-temporal
resolution, which is the basis of numerous binding affinity calculation
approaches, such as MM/PB(GB)SA and alchemical free energy
frenkel2001understanding ; valdes2021gmx_mmpbsa ; song2020evolution .
Nevertheless, molecular dynamics simulation usually consumes a large amount of
computational resources, which may thus impede its applications to high-
throughput virtual screening liu2016applying .

Data-driven approaches are another emerging trend for solving the binding
affinity prediction task lo2018machine . Primarily benefiting from the rapidly
growing amount of experimental protein-ligand complexes, data-driven
approaches have achieved a certain accuracy level through learning the binding
modes directly from three-dimensional structures. Nonetheless, current state-
of-the-art methods often require sophisticated feature engineering while
offering limited generalizability zhu2022assessment . More importantly,
various attempts such as docking-based verdonk2003improved ; friesner2004glide
; trott2010autodock ; korb2009empirical , feature-based dittrich2018converging
; mooij2005general ; su2018comparative ; wang2002further ;
ballester2010machine ; wang2017improving ; sanchez2021extended , voxel-based
stepniewska2018development ; jimenez2018k ; zheng2019onionnet  and graph-based
nguyen2021graphdta ; yang2019analyzing ; lim2019predicting ;
maziarka2020molecule ; klicpera2020directional ; song2020communicative ;
li2021structure  methods, mainly focus on tuning models on the static
structures, i.e., the crystallographic structures. In fact, both proteins and
ligands are not rigid structures and they constantly interact with each other
and the aqueous environment, following a specific conformation distribution.
Such a conformation distribution in an equilibrium state is called a
thermodynamic ensemble, and all the conformations collectively contribute to
the affinity of the protein-ligand binding. Unfortunately, most of the
experimental structures are resolved by X-ray crystallography, which only
catches a snapshot of the most populated states, and the conformation
distributions of the structures are inevitably lost. On the other hand, the
data-driven approach could fully make use of previously experimentally
resolved structures and measured binding affinities for better generalization.
Despite all the benefits, to the best of our knowledge, no large-scale MD
protein-ligand dataset has been hitherto reported, and few previous data-
driven approaches have taken such thermodynamic ensembles into account to
predict protein-ligand binding affinities.

To alleviate the aforementioned issues, we first curated a large-scale MD
trajectory dataset, containing 3,218 different protein-ligand complexes, based
on the PDBBind dataset wang2004pdbbind . Then, we introduced a graph
transformer framework, named Dynaformer, that was trained based on this MD
trajectory dataset. To capture the binding modes between proteins and ligands,
Dynaformer utilizes a roto-translation invariant feature encoding scheme,
taking various interaction characteristics into account, including interatomic
distances, angles between bonds, and various types of covalent or non-covalent
interactions. Compared to other existing methods, Dynaformer showed better
scoring and ranking power on the CASF-2016 benchmark dataset su2018comparative
. Through several case studies, we also discussed the underlying thermodynamic
binding mechanisms that may contribute to the performance improvement achieved
by Dynaformer. Furthermore, we experimentally validated several novel hits
identified by Dynaformer against heat shock protein 90 (HSP90), an attractive
drug target for cancer. We demonstrated that our model can be applied to a
real-world hit discovery task and serve as an effective tool to accelerate the
drug discovery process.

##  2 Results

###  2.1 The thermodynamic ensemble determines the binding affinity

A thermodynamic ensemble of a protein-ligand system represents a distribution
of complex conformations in equilibrium. Theoretically, these conformations,
with varying probabilities of occurrence, collectively depict the free energy
state of the system laio2008metadynamics ; christ2010basic ; wei2016protein .
Therefore, conformation snapshots sampled from the thermodynamic ensemble can
approximate the distribution of conformation space and further reflect the
free energy state of ligand binding. In other words, an MD trajectory can
provide more information about ligand binding than a single static structure
(Fig. 1A). More specifically, the binding affinity $K_{i}$ is primarily
determined by the Gibbs free energy change $\Delta G$ of the binding process,
which can be defined as:

$\begin{split}-RT\textrm{ln}K_{i}=\Delta G=\Delta G_{gas}+\Delta
G_{solv},\end{split}$ (1)

where $R$ represents the molar gas constant, $T$ represents the temperature,
$\Delta G_{gas}$ stands for the binding free energy between the pair of
protein and ligand in the gas phase, and $\Delta G_{solv}$ stands for the
solvation free energy difference between the protein-ligand complex and the
sum of stand-alone protein and ligand systems. Both $\Delta G_{gas}$ and
$\Delta G_{solv}$ can be represented by the ensemble-averaged energy terms
related to the atomic features derived from each snapshot (see Methods for
more details). Therefore, in principle, with sufficiently sampled snapshots
from an MD simulation, the binding affinity can be expressed as a function of
the atom features of the snapshots.

Based on the above principle, training deep learning models to model the
energy terms using the thermodynamic ensembles sampled from MD simulations can
enable one to predict the protein-ligand binding affinities. Here, we first
built a large-scale dataset containing MD simulations of 3,218 protein-ligand
complexes based on experimental structures derived from the Protein Data Bank
(PDB) berman2000protein . All of the complexes were further filtered based on
the PDBBind dataset wang2004pdbbind , consisting of protein and small molecule
complex structures with known binding affinities. In this study, we performed
a 10-nanosecond (ns) simulation for each complex and sampled 100 snapshots
from each simulation to characterize the conformational space. All snapshots
were collected to derive the MD trajectory dataset.

To gain comprehensive insight into the relationships between complex
stabilities and corresponding binding affinities, we also performed a deep
analysis of our MD trajectory dataset. More specifically, we first derived the
relative position movement of each ligand during the simulation by calculating
the average root mean squared deviation (RMSD) between its initial and current
positions of 100 snapshots at the binding site. Then, we examined the
association between the mean RMSDs of ligands and experimentally measured
binding affinities. Our analysis revealed that conformations from MD
simulations exhibited varying levels of stability, and the stabilities of
ligands were roughly inversely correlated with binding affinities (Fig. 1B).
More specifically, the upper bound for the mean RMSD of the ligand decreased
as the binding affinity increased. In total, 78.3% of trajectories in the
dataset exhibited high stability, with mean RMSDs of ligands $<3\mbox{\r{A}}$;
20.6% of trajectories displayed intermediate flexibility, where certain atoms
or functional groups in the ligands were flexible but did not leave the
original binding sites, with mean RMSDs of ligands between $3\mbox{\r{A}}$ and
$10\mbox{\r{A}}$; and 1.1$\%$ of complexes were unstable, with mean RMSDs of
ligands $>10\mbox{\r{A}}$, exhibiting a high level of ligand flexibility. In
the unstable cases, the binding affinities tended to be low, and ligands even
left the original binding pockets and entered the solvent. Several typical
examples are shown in Fig. 1C. In the unstable example (PDB ID: 4y3j) with the
lowest binding affinity, the ligand left the binding site shortly after the
simulation began, while in the intermediate (PDB ID: 3udh) and stable (PDB ID:
2yge) examples, the ligands remained in the binding sites, displaying
relatively higher binding affinities. This preliminary data analysis validated
that ligand dynamics can provide extra information with respect to the
corresponding binding affinity. Therefore, our MD trajectory dataset can
contain more information than the static crystal structure counterparts, and
thus should be able to contribute to binding affinity prediction.

###  2.2 Our virtual screening pipeline with Dynaformer

The overview of our virtual screening pipeline using Dynaformer is shown in
Fig. 2A, including the MD trajectory dataset preparation, model training,
performance evaluation, virtual screening and wet-lab validation. To feed data
into the deep learning model, snapshots from MD trajectories are converted
into graphs representing the ligand binding structures. Each snapshot is
converted into a graph containing nodes and edges that represent features from
atoms and their relationships, such as chemical bonds and non-covalent
interactions (Fig. 2B). A graph includes atoms and the inteatomic interactions
from both the ligand and the protein within a certain distance cutoff. After
preparing the MD trajectory dataset, we propose a novel graph transformer
architecture, named Dynaformer (Fig. 2C), to fully take advantage of the
characteristics of molecular dynamics from the dataset and improve the
predictions of protein-ligand binding affinities. Dynaformer is developed
based on the Graphormer, which has been proven to be effective for the
property prediction of small molecules ying2021transformers . Nevertheless,
the original Graphormer version was designed for encoding the features of
small molecules, that is, it can only encode the feature of isolated molecules
with covalent bonds, and failed to handle the non-covalent bonds between the
proteins and ligands, despite that these non-covalent interactions with
divergent directionality and saturation are actually the essential
intermolecular forces that “glue” the ligand to the binding site. To capture
these intricated non-covalent interaction features, Dynaformer introduces
additional attention bias terms to encode the structural information regarding
the interatomic distances and angles between these bonds. In addition, the
Gaussian basis function (GBF) is used for structural encodings, which is a
common form of energy function in computational chemistry to map structural
information to energy contributions (Eq. 5 in Methods) bartok2010gaussian . In
this way, from the graph representation of a protein-ligand complex, we are
able to model the protein-ligand interactions of a snapshot from an MD
trajectory. Moreover, to model longer-range interactions that may affect the
ligand binding beyond the distance cutoff of the graph representation, three
fingerprints from ECIF sanchez2021extended , RF-score ballester2010machine ,
and GB-score rayka_firouzi_2022  are employed and fused into the final linear
layer.

As previously mentioned, MD simulations are notoriously time-consuming, making
MD trajectory data rarely used in high-throughput virtual screening. To ensure
practical applicability and accurate prediction in realistic scenarios, we
employ a pretraining-finetuning strategy in our framework. More specifically,
we first pretrain Dynaformer on the graphs derived from the MD trajectory
dataset by predicting the corresponding binding affinities. Then, for the
benchmarking test and the virtual screening applications, the pretrained model
is further finetuned on the crystal structures from the same dataset, PDBBind.
In this way, the valuable patterns learned from the MD trajectories that may
be lost in the crystal structures can be effectively transferred to the
finetuning process. More details about the Dynaformer architecture and the
pretraining-finetuning strategy can be found in the Methods section. After the
pretraining and finetuning steps, the Dynaformer model can serve as a scoring
function for novel hit discovery in a real-world virtual screening
application. More specifically, the docking structures between the target of
interest and the molecules from a compound library are fed into the model to
predict the binding affinity. The top-ranked molecules are then selected as
the hit molecules and further investigated in the wet-lab experiments (Fig.
2A).

###  2.3 Molecular dynamics improves the binding affinity prediction

To evaluate the performance of Dynaformer on binding affinity prediction, we
tested it on the CASF-2016 benchmark dataset, which consisted of 285 protein-
ligand complexes covering 57 different target classes su2018comparative . More
specifically, we compared the scoring and ranking power between Dynaformer and
a range of baseline methods, including docking-based verdonk2003improved ;
friesner2004glide ; trott2010autodock ; korb2009empirical , feature-based
dittrich2018converging ; mooij2005general ; su2018comparative ;
wang2002further ; ballester2010machine ; wang2017improving ;
sanchez2021extended , voxel-based stepniewska2018development ; jimenez2018k ;
zheng2019onionnet  and graph-based nguyen2021graphdta ; yang2019analyzing ;
lim2019predicting ; maziarka2020molecule ; klicpera2020directional ;
song2020communicative ; li2021structure  methods, following the same
evaluation protocols as provided by CASF-2016. The scoring power of a model
refers to its ability to accurately predict the binding affinities, which can
be measured using metrics Pearson’s correlation coefficient (Pearson
$\bm{r}$), root-mean-square error (RMSE), and standard deviation (SD), while
the ranking power of a model indicates its ability to correctly rank the
relative order of binding ligands, which can be measured using metrics
Spearman’s coefficient (Spearman $\bm{\rho}$), Kendall’s coefficient (Kendall
$\bm{\tau}$), and predictive index (PI). The scoring and ranking performances
of different methods are shown in Fig. 3A and Fig. 3B, respectively. The test
results showed that, Dynaformer outperformed all the baseline methods. In
addition, the prediction scatter plot revealed a high correlation with Pearson
$r$ of 0.858 and a low prediction bias with RMSE of 1.114 between the
Dynaformer predicted and experimentally measured binding affinities Fig. 3C.

We next carried out comprehensive ablation studies on Dynaformer, to validate
the efficacy of its specific designs (Fig. 3D). We first introduced multiple
variations of the Dynaformer model, such as -F (i.e., Dynaformer without using
the fingerprints), -3D (i.e., without using structural encoding module), and
-M (i.e., without pretraining on the MD dataset). Then, we compared the
prediction performances of Dynaformer and its variations using Pearson
$\bm{r}$ and Spearman $\bm{\rho}$. Dynaformer achieved superior performances
in comparison with its variations, thus indicating the effectiveness of
individual modules employed in Dynaformer. More specifically, as expected, the
pretraining on the MD trajectory data played the most significant role in
enhancing the scoring and ranking power of Dynaformer. In addition, the 3D
structural encoding module, which encodes distance and angle information using
the GBF module, also played a crucial role in improving the prediction
performances. Moreover, incorporating the pre-calculated fingerprints also
yielded performance gain probably due to the increased receptive field. It is
noteworthy that the -F-3D-M model, a standard transformer model trained on
static crystal structures, performed comparably to the graph neural network
baseline models. Overall, this ablation study highlighted the relative
importance of individual feature encoding modules in Dynaformer.

To demonstrate that incorporating the pretraining on the MD dataset is
necessary, we further investigated whether Dynaformer can capture the
intrinstic binding-related patterns from the MD trajectory dataset, such as
the entropy feature, the enthalpy feature, and the subtle differences in
protein-ligand interactions. More specifically, we performed three case
studies from the test set of CASF-2016, where Dynaformer offered more accurate
predictions in comparison with baseline methods. We employed ECIF and SIGN as
baselines for comparison in these case studies, which were a fingerprint-based
and a graph-based binding affinity prediction model trained on static crystal
structures, respectively.

In the first case, we showed that Dynaformer can predict the binding affinity
more accurately through learning the ligand flexibility, i.e., the additional
entropic information derived from the MD trajectory. In the protein-ligand
complex structure shown in Fig. 4A (PDB ID: 2v7a), which is a ligand bound to
the T315I Abl kinase domain modugno2007crystal . As can be seen in the zoom-in
figure, the head of the ligand (i.e., the pyrrolo[3,4-c]pyrazole group) is
tightly buried in the binding pocket, whose contacts with the protein are
highly conserved during the MD simulation. In detail, the carbonyl oxygen of
residue E316 and the amide nitrogen of residue M318 are found to interact with
the two nitrogen atoms of the pyrrolopyrazole scaffold, and the nitrogen of
the amide group forms a hydrogen bond with the carbonyl oxygen of residue
M318. In addition, the benzyl group is also involved in the hydrophobic
interactions with residue L370. These interactions on the head of the ligand
lead to a considerably favorable enthalpy change, due to the formation of non-
covalent interactions. In addition, the tail of the ligand (i.e., the
N-methylpiperazine group) is exposed to the solvent and wiggles randomly,
which is also entropically favorable in ligand binding. This wiggling can be
confirmed from the MD trajectories, as shown in the root mean square
fluctuation (RMSF) plot in Fig. 4A, which measures the fluctuation of an atom
around its average position. Consequently, this ligand showed a strong binding
affinity, i.e., $\mbox{p}K_{i}=8.3$. In the prediction results, ECIF and SIGN
underestimated the binding affinity by a relatively large margin, while
Dynaformer predicted the binding affinity more accurately. In this case, both
the static crystal structure and MD trajectory data contain enthalpy
information, but the MD trajectory data includes additional entropic
information, i.e., the features of ligand flexibility, which actually plays a
crucial role in determining the binding free energy. Through incorporating
this extra information, Dynaformer was able to achieve better prediction
accuracy.

In the second case, we show that Dynaformer can more accurately predict the
binding affinity by modeling the enthalpy change, which is the energy term
caused by the ligand-protein interactions. As demonstrated in Fig. 4B, the
aspartyl protease $\beta$-secretase (BACE) is in complex with a compound in
the spiropyrrolidine scaffold (PDB ID: 3udh), which was screened from a
fragment library efremov2012discovery . Unlike the molecule shown in the first
case, the ligand here has a rigid structure, i.e., a single conformational
state and zero rotatable bonds. Therefore, its binding affinity depends mainly
on the enthalpy change caused by the interactions with the binding site. In
the determined structure, the pyrrolidine nitrogen is located between the two
catalytic acids on the side chain of residues D32 and D228, forming stable
hydrogen bonds. In addition, the oxindole and the phenyl ring next to the
pyrrolidine ring exhibit a strong complementarity in terms of the occupancy of
the binding site. However, the adjacent subpockets are still left unoccupied,
as it was the primary optimization direction as discussed in the original work
efremov2012discovery . As a result, the binding affinity of the ligand was as
low as $\mbox{p}K_{i}=2.8$. In our MD trajectories, the ligand gradually moved
away from its initial position after a few nanoseconds, as evidenced by an
average ligand RMSD of over $5\mbox{\r{A}}$. Such behavior is typically
associated with a weak binding affinity, despite some favorable protein-ligand
interactions. Our test validated that the Dynaformer can learn such knowledge
through training based on the MD trajectories, thus providing a more accurate
prediction ($\mbox{p}K_{i}=3.6$). Similar behaviors were also observed in
other protein-ligand complexes with low binding affinities, such as the
complexes with PDB IDs 4y3j and 3gv9.

With the third case shown in Fig. 4C, Dynaformer demonstrated its ability to
distinguish ligands with subtle differences in protein-ligand interactions,
thus allowing it to tackle the “activity cliff”, where similar ligands have
distinct binding affinities. Here, three ligands with a thiophene scaffold
bind to PTP1B (protein tyrosine phosphatase 1B) with varying potency, in which
the ligand of PDB ID 2qbp showed ten times stronger binding affinity compared
to that of PDB ID 2qbq, and the ligand of PDB ID 2qbq showed ten times
stronger binding affinity compared to that of PDB ID 2qbr wilson2007structure
. 2-phenylthiophene, as the structural core of these three compounds, adopts
an identical binding pose to occupy the active site, and the tail groups enter
another hydrophobic subpocket. The ligand RMSD plots in Fig. 4D indicated that
the ligand from 2qbp exhibited high conformational stability at the binding
site (mean ligand $\mbox{RMSD}<1$Å), which was likely due to the highly
conserved $\pi$-$\pi$ stacking with residue F182 and hydrogen bonds with
residues K120 and R221. On the other hand, even though the other two ligands
shared similar interactions as 2qbp, the ligand of PDB ID 2qbq partially left
the original binding site after 1 ns of simulation, and the structural core of
the ligand of PDB ID 2qbr finally left the binding site after around 9 ns of
simulation. Fig. 4E displayed the percentage of snapshots that key
interactions observed during the simulation. Such statistics information was
another piece of evidence showing how conserved the interactions were. The
occupancy of the hydrophobic subpocket by the different tail groups decreased
progressively with each modification in the ligand structure, from 2qbp to
2qbq and finally to 2qbr. Therefore, these structurally similar ligands
provided distinct binding strength, with $\mbox{p}K_{i}=8.4,7.4,6.3$ for 2qbp,
2qbq, and 2qbr, respectively. Through the rich interaction information learned
from the MD trajectories and its ability to model such interactions,
Dynaformer achieved a better ranking result compared to those of baselines
ECIF and SIGN (Fig. 4C), and thus may provide more reliable guidance toward
the downstream hit discovery task.

In summary, the above detailed analyses showed that the MD trajectories may
provide useful insights in understanding the protein-ligand binding
affinities, and our Dynaformer model was capable of capturing such rich
features for achieving better prediction results.

###  2.4 Hit discovery for the HSP90 target

To verify the applicability of Dynaformer in realistic scenarios such as hit
discovery in early drug discovery process, we applied Dynaformer to discover
potent hit compounds through scoring the docked protein-ligand structures. The
original docking score often contains many false positives among the top-
scoring candidates lyu2019ultra . Consequently, numerous downstream
unnecessary experiments are often required to identify the true hit compounds.
In this study, we demonstrate that Dynaformer can more efficiently deliver hit
compounds through scoring the docked poses. In particular, we selected the
heat shock protein 90 (HSP90) as our target protein, which is a vital
chaperone protein involved in the important biological pathways of many
refractory diseases, including cancer, neurodegenerative diseases, and viral
infection. To inhibit the HSP90 function, a common strategy is to design small
molecules targeting its ATP binding pocket to suppress the ATPase activity in
the N-terminal domain. Here, our training data from the PDBBind dataset only
included 18 protein-ligand structures associated with ATPase activity with
$\mbox{p}K_{i}$ values ranging from 3.84 to 8.32. All of these structures can
be categorized as known inhibitor scaffolds summarized in the literature
janin2010atpase ; li2019heat .

Next, we performed docking on HSP90 against the ChemBridge DIVERSet-EXP
library containing 50,000 small molecules that cover a diverse pharmacophore
space. The docked conformers were obtained using Autodock Vina, a widely used
open-source docking software eberhardt2021autodock ; trott2010autodock . The
reference protein-ligand structure (PDB ID: 2xdl) was chosen from the test
set, CASF-2016, to identify the binding pocket. Then, the docked poses were
then used to score and rank the molecules based on the Dynaformer predictions.
After that, the top-ranked compounds were filtered and visually inspected for
prioritization, following the approach and criteria described in the
literature bender2021practical ; fischer2021decision  (see Methods for
details). To ensure the fairness of prioritization, we created a pool
consisting of molecules among the top 10% from both Dynaformer and the
original scoring function in Autodock Vina. From the pool, 20 molecules were
chosen without knowing their rankings. Finally, the surface plasmon resonance
(SPR) experiments were conducted for each of the 20 molecules to measure the
inhibition constant $K_{i}$. In total, 12 of the 20 molecules had measurable
binding affinities.

Fig. 5A displays the chemical structures, experimentally measured $K_{i}$ and
SPR sensorgrams of the 12 molecules that showed measurable binding affinities.
The predicted and experimentally measured $K_{i}$ values are summarized in
Fig. 5B, for comparing the prediction power of Dynaformer and Autodock Vina.
The only difference here between Dynaformer and Autodock Vina was the scoring
function used for binding affinity prediction and the input dockeing
structures remained the same. The results demonstrate that Dynaformer is able
to predict highly potent molecules more effectively than the scoring function
used in Autodock Vina, and the Pearson $\bm{r}$ coefficient indicates that the
performance of Dynaformer (0.72) outperforms Autodock Vina (0.24) by a large
margin (Fig. 5C). More specifically, Dynaformer ranked the top three compounds
in the library at 0.78%, 5.10%, and 0.06%, which is much higher than Autodock
Vina’s ranking. In addition, the possible binding modes of the top three
compounds were consistent with the known binding patterns from the previous
studies of rational drug design (Fig. 5D). For example, previous
investigations revealed that 1) the hydrogen bond networks between inhibitors
and residues D93, T184, and K58 are crucial for stable binding; and 2)
enhanced affinity and selectivity can be achieved when the lipophilic bottom
of the pocket (including residues M98, F138, Y139, V150 and V186) is occupied
by the aromatic moiety of the ligand kung2008dihydroxylphenyl ;
murray2010fragment ; stebbins1997crystal ; janin2010atpase ; li2019heat ;
yoshimura2021thermodynamic .

Additionally, we found that some scaffolds of molecules with measurable
binding affinities shown in Fig. 5B have been validated by previous studies.
For instance, Compound 2, 7, 10 and 11 were inhibitors containing the
resorcinol-like scaffold discovered from natural products. Compound 5 and 9
shared a similar backbone belonging to the purine scaffold, while Compound 8
was similar to an aminotriazine compound previously reported in a fragment-
based screening li2019heat ; huth2007discovery . Notably, we found that
Compounds 1, 3, 4, 6, 12 have the potential to serve as promising novel
candidates for further hit to lead optimization. To the best of our knowledge,
there is no known precedent in the literature for these compounds, making them
particularly valuable for future lead optimization. In summary, our wet-lab
validation results demonstrated that Dynaformer can deliver effective hit
compounds with favorable binding affinities and novel scaffolds, thus
facilitating the drug discovery process.

##  3 Conclusion

Protein-ligand binding affinity prediction is a fundamental problem in early
drug discovery, as it can greatly help reduce the immense costs associated
with wet-lab experiments in lead compound design. In this study, we first
curated a comprehensive MD trajectory dataset consisting of 3,218 complexes
containing proteins and bound small molecules. Each trajectory, with
conformation snapshots from the thermodynamic ensemble, approximates the
conformation space distribution and represents the ligand binding free energy
state. We then developed Dynaformer, a graph transformer model that learns the
underlying physico-chemical patterns of ligand binding from this MD trajectory
dataset. Through extensive testing and ablation studies, we showed that the MD
trajectories can contain rich features related to enthalpy and entropy, and
thus provide more information than static crystal structures. Our results
demonstrated that Dynaformer significantly outperformed baseline approaches on
the CASF-2016 benchmark dataset, highlighting the importance of incorporating
MD trajectory data to enhance the scoring and ranking power of binding
affinity prediction. We further illustrated the efficacy of Dynaformer as a
scoring function in a real-world hit discovery process against HSP90. By
scoring the docked poses of a library containing 50,000 molecules, Dynaformer
can effectively identified hit compounds with favorable binding affinities and
novel scaffolds. Among the 20 experimentally tested compounds, 12 molecules
exhibited measurable binding affinities, including two compounds with
submicromolar $K_{i}$ values.

Our work underscores the potential of employing MD trajectories and deep
learning models to enhance binding affinity prediction. In both a benchmark
dataset and a real-world drug discovery scenario, we have illustrated
effectiveness of Dynaformer as a scoring function and its ability to discover
novel hit compounds. We anticipate that prediction performance can be further
improved by incorporating more high-quality data and carefully designed
training tasks.

##  4 Methods

###  4.1 Decomposition of the free energy related to the protein-ligand
binding affinity

In this section, we show that the free energy $\Delta G$ of protein-ligand
binding, which determines the binding affinity, can be decomposed into
multiple components. Each component can be represented as a function
associated with atomic features from the MD trajectories. This provides a
solid theoretical basis for learning structure-to-affinity relationships using
deep learning models. Generally, $\Delta G$ can be expressed as the sum of two
components, the gas-phase binding free energy and the solvation free energy,
that is,

$\begin{split}-RT\textrm{ln}K_{i}&=\Delta G=\Delta G_{gas}+\Delta
G_{solv},\end{split}$ (2)

where $\Delta G_{gas}$ represents the interaction energy between the given
pair of protein and ligand at their gas phase, and $\Delta G_{solv}$ stands
for the solvation free energy difference between the protein-ligand complex
and the sum of stand-alone protein and ligand systems. $R$ and $T$ stand for
the universal gas constant and temperature, respectively, and $K_{i}$ stands
for the thermodynamic equilibrium constant of the binding process.

The energy term $\Delta G_{gas}$ is the binding free energy between protein
and ligand in the gas phase, and it can be further split into two parts, i.e.,
the enthalpy term $\Delta H_{gas}$ and the entropy term $\Delta S_{gas}$.
Following the calculation scheme of the interaction entropy method
duan2016interaction ; sun2017interaction , both $\Delta H_{gas}$ and $\Delta
S_{gas}$ can also be expressed with the interaction energies $E^{int}$, that
is,

$\begin{split}\Delta G_{gas}&=\Delta H_{gas}-T\Delta S_{gas},\\\ &=\langle
E^{int}\rangle+RT\textrm{ln}\langle e^{\beta\Delta
E^{int}}\rangle,\end{split}$ (3)

where the enthalpy $\Delta H_{gas}$ and entropy $-T\Delta S_{gas}$ terms
correspond to $\langle E_{pl}^{bind}\rangle$ and $RT\textrm{ln}\langle
e^{\beta\Delta E^{int}}\rangle$, respectively, $\beta$ stands for a constant
$1/RT$, $\langle E^{int}\rangle$ stands for the trajectory-averaged
interaction energy of protein and ligand contributed by different interaction
types, such as hydrogen and hydrophobic interactions, Van der Waals
interactions, $\pi$-$\pi$ stackings and other long-range interactions. More
specifically, $E^{int}$ could be expressed as $\sum\sum u(i,j)$, which is the
sum of the interaction energy between each pair of atoms $i$ and $j$ from the
protein and ligand, respectively. Here, $\Delta E^{int}=E^{int}-\langle
E^{int}\rangle$ represents the fluctuation of the interaction energy around
the trajectory-averaged value.

The $\Delta G_{solv}$ term is the solvation free energy difference between the
protein-ligand complex and the sum of stand-alone protein and ligand systems.
It can be further split into two parts, i.e., the non-polar solvation free
energy $\Delta G_{np}$ and the polar solvation free energy $\Delta G_{pol}$,
that is,

$\begin{split}\Delta G_{solv}&=\Delta G_{np}+\Delta G_{pol},\\\
&=\langle\sum\left(\gamma\mbox{SASA}(i)+\delta\right)\rangle+\langle\sum
q_{i}\Phi(i)\rangle,\end{split}$ (4)

where the non-polar solvation free energy $\Delta G_{np}$ term, can be
obtained using an empirical solvent-accessible surface area (SASA) formula
parameterized by $\gamma$ and $\delta$, and $\Delta\mbox{SASA}$ stands for the
change of solvent accessible surface area (SASA) upon binding huang2018fast ;
pratt1977theory . The polar solvation free energy $\Delta G_{pol}$ term is
controlled by the electrostatic potential $\Phi(i)$, which can be obtained
through solving the Poisson-Boltzmann (PB) equation or the Generalized Born
(GB) equation massova2000combined . The $\Delta G_{pol}$ can then be
calculated through summing over the electrostatic potential $\Phi(i)$ at the
position of each atom $i$ weighted by the partial charge $q_{i}$ of atom $i$.

In summary, each of the above components of the binding free energy change can
be expressed as a function of the features of the atoms or atom pairs involved
in the ligand binding. Therefore, they can be naturally represented and
approximated through training our Dynaformer on the large-scale MD trajectory
data. Since each term is actually averaged over the MD trajectory, the results
can be more robust compared to that obtained from a single static crystal
structure.

###  4.2 Preparation of the MD trajectory dataset

We used the protein-ligand complexes from the PDBBind dataset wang2004pdbbind
as the initial conformations for our molecular dynamics simulations. PDBBind
contains a comprehensive collection of high-quality crystal complex structures
of proteins and ligands curated from the Protein Data Bank (PDB)
berman2000protein , together with the corresponding experimentally measured
binding affinity data. More specifically, the PDBBind dataset consists of
three parts, i.e., general, refined, and core sets. The general set includes
21,382 complexes, while the refined set is a subset of the general set
containing 4,852 structures, selected by a set of rules based on the quality
of crystal structures and binding data liu2017forging . The core set, which is
separately organized as CASF (Comparative Assessment of Scoring Functions), is
a stand-alone test set consisting of 285 high-quality complexes covering 57
classes of targets for evaluating different computational models. The general
and refined sets are updated annually, while the CASF set is not. Here, we
evaluated our model and all baseline models on the latest version of CASF,
i.e., CASF-2016.

Subsequently, we used the web-based CHARMM-GUI to clean up the data and
prepare their corresponding MD simulation input files jo2008charmm ;
lee2016charmm . The MD simulations were performed based on the CHARMM36
forcefield lee2016charmm  and using NAMD engine, which is a parallel molecular
dynamics program designed for high-performance simulation of large
biomolecular systems phillips2020scalable . For each MD simulation, the
complex was solvated in a truncated periodic TIP3P water box, and the minimum
distance from the surfaces of the box to the complex atoms was set to 10 Å.
Counter ions were added to neutralize systems, and the initial configuration
was decided using a short Monte Carlo simulation. The ligand topology and
parameter file were generated using the ParamChem service ghosh2011molecular .
The missing residue was modeled using GalaxyFill coutsias2004kinematic . The
simulation temperature was maintained at 303.15 K. For a stable protein-ligand
complex, several nanoseconds are usually sufficient for free energy
calculation, as in MM/PB(GB)SA genheden2010obtain . However, severe structural
stability issues, such as ligand drifting and protein unfolding, become more
obvious if the MD simulation duration is longer. Therefore, to tradeoff the
simulation time and computational complexity, we finally chose a simulation
time of 10 ns. 0.5 ns NVT (constant volume and temperature) was set for the
equilibrating stage before 10 ns NPT (constant pressure and temperature)
simulation was performed.

In practice, we excluded any complex associated with the following conditions,
which led to failed MD simulations: 1) covalent ligands that had a covalent
bond with the protein; 2) complicated or multiple ligands that led to failed
forcefield parametrization; 3) the protein structure contained too many
missing residues and failed to repair; and 4) specific types of proteins, such
as membrane proteins, which were not suitable for simulation in a water box.
After removing the above cases, 3,218 protein-ligand complexes were selected
from the PDBBind refined set. For each complex, 10 ns simulation was
performed, using the previously mentioned protocol, and from each simulation,
100 snapshots were sampled to form a trajectory.

###  4.3 Data preprocessing for Dynaformer

In our study, we converted the structure of a protein-ligand complex, obtained
from either an MD trajectory or a crystal structure, into a graph
representation, which was subsequently fed into Dynaformer as input data. More
specifically, we denote the ligand, the binding pocket, and their interatomic
relationships as a graph, $\mathcal{G}=(\mathcal{V},\mathcal{E})$, as depicted
in Fig. 2B. Here, the node set
$\mathcal{V}=\left\\{\mathbf{v}_{0},\mathbf{v}_{1},\mathbf{v}_{2},\dots,\mathbf{v}_{n}\right\\}$
represents the set of input atom feature vectors, where $n$ stands for the
number of atoms and $\mathbf{v}_{0}$ refers to a virtual node employed as a
global representation of the graph $\mathcal{G}$. The graph $\mathcal{G}$
comprised of atom features derived from the ligand and the binding pocket. The
binding pocket here is characterized by atoms that belong to the protein and
are situated within a specific distance ($d_{p}$) from any atom of the ligand.
The edge set
$\mathcal{E}=\left\\{\mathbf{e}_{i,j}~{}|~{}i,j\in\mathcal{V}\right\\}$ stands
for the set of covalent or non-covalent bond feature vectors, with
$\mathbf{e}_{i,j}$ indicating such bond feature vector between atoms $i$ and
$j$. Here, the covalent interactions occur between atoms possessing chemical
bonds, whereas non-covalent interactions transpire between any two atoms
within a distance threshold ($d_{s}$).

We adopt the same features as used in Open Graph Benchmark (OGB) hu2020open ,
and the comprehensive definitions of node and edge features can be found in
Supplementary Table 1. We extract the snapshots from MD trajectories using the
MDAnalysis library michaud2011mdanalysis  and feature extraction is performed
with the RDKit library landrum2013rdkit . In our experiments, we set
$d_{p}=d_{s}=5$Å.

###  4.4 The model architecture of Dynaformer

Dynaformer is a transformer-based model that learns the protein-ligand binding
features from molecular dynamics data. Building based on Graphormer
ying2021transformers , we propose a roto-translation invariant encoding scheme
to encode the atomic features. At a high level, the Dynaformer predicts the
binding affinity given the collective interaction patterns of a ligand at the
active binding site.

First, we highlight how Dynaformer encodes the distance and angle features
between atoms. As previously stated, binding free energy depends on the atomic
features and interactions of atom pairs in the complex. To represent these
scalar values for distances and angles, we utilize a Gaussian basis function
(GBF). The GBF is advantageous in modeling non-linear relationships between
atom pairs, such as Van der Waals and electrostatic potentials
cisneros2012application . Moreover, the GBF is a smooth function that can be
seamlessly integrated into neural networks. The specific calculations for
structural encodings are performed using the subsequent formulas:

$\displaystyle\mbox{GBF}(x_{k})$
$\displaystyle=\exp\left[-\frac{(x_{k}-\mu_{k})^{2}}{2\sigma_{k}^{2}}\right],k=1,\dots,K,$
(5) $\displaystyle d(i,j)$
$\displaystyle=\mbox{GBF}\left(\left[\mathbf{v}_{i}|r_{ij}|\mathbf{v}_{j}\right]\mathbf{W}_{d}\right),$
$\displaystyle a(i,j)$
$\displaystyle=\mbox{GBF}\left(\left[\mathbf{v}_{i}|\sum_{k}\angle{ijk}|\mathbf{v}_{j}\right]\mathbf{W}_{a}\right),$

where $\mu_{k}$ and $\sigma_{k}$ stand for the learnable parameters, $K$
stands for the number of encoding heads, $x_{k}$ stands for a distance or
angle value for encoding, $\mathbf{v}_{i}$ stands for the atomic features of
atom $i$, $~{}|$ stands for the concatenation operation between vectors,
$r_{ij}$ represents the Euclidean distance scalar in Ås between atoms $i$ and
$j$, $\angle{ijk}$ indicates the angle scalar in degrees between bonds formed
by atoms $i$, $j$, $k$, and $\mathbf{W}_{d}$ and $\mathbf{W}_{a}$ correspond
to the weight matrices for distance and angle features, respectively.

Then, we briefly introduce the Dynaformer architecture and emphasize the self-
attention module. Dynaformer has two main components in each encoder layer: a
multi-head self-attention module (MHA) and a feed-forward network (FFN). The
encodings of structural features as described in Fig. 2C are fused into an MHA
module, which then updates the feature representations of each atom and their
contributions to the final binding affinity. In each head of MHA, the
attention is calculated as follows:

$\displaystyle\mbox{Attention}(\mathbf{Q},\mathbf{K},\mathbf{V})=\mbox{softmax}\left(\mathbf{A}\right)\cdot\mathbf{V},$
(6) $\displaystyle
A_{ij}=\frac{\left(\mathbf{h}_{i}\mathbf{W_{Q}}\right)\left(\mathbf{h}_{j}\mathbf{W_{K}}\right)^{T}}{\sqrt{d_{K}}}+d(i,j)+a(i,j)+\frac{1}{P}\sum_{p=1}^{P}\mathbf{e}_{p}\mathbf{w}_{p}^{T},$
(7)

where
$\mathbf{Q}=\mathbf{H}\mathbf{W_{Q}},\mathbf{K}=\mathbf{HW_{K}},\mathbf{V}=\mathbf{HW_{V}}$
stands for the query, key and value matrices in transformer model,
$\mathbf{H}=\\{\mathbf{h}_{i}\\}$ stands for hidden representations of nodes,
matrix $\mathbf{A}$ stands for the self-attention matrix, $A_{ij}$ is the
unnormalized attention score between atoms $i$ and $j$, and $d_{K}$ stands for
the dimensionality of the query, key, and value, which serves as a
normalization term. In addition, the term
$\frac{1}{P}\sum_{p=1}^{P}\mathbf{e}_{p}\mathbf{w}_{p}^{T}$ encodes the
interatomic interactions on the shortest path between atom $i$ and $j$, where
$\mathbf{w}_{p}$ denotes the weight that projects edge feature vector
$\mathbf{e}_{p}$ to a scalar and $P$ stands for the maximum length of the
shortest path to keep the calculation feasible. This term is an extension of
the adjacency matrix that encodes the topological relationships between any
pairs of atoms.

In this MHA module, the feature representations of nodes are updated through
the attention matrix $\mathbf{A}$, which is a weighted summation of all node
representations. In every layer of Dynaformer, the atomic representations
gradually exchange information and aggregate to the representation of the
virtual node. After the last encoder layer, the feature representation of the
virtual node is fused with the three pre-calculated fingerprints from ECIF
sanchez2021extended , RF-score ballester2010machine  and GB-score
rayka_firouzi_2022 , for predicting the binding affinities. Such a scheme
allows the Dynaformer to learn the effective structural patterns contributing
to the final binding affinities. A detailed explanation of the model
architecture is stated in Supplementary Information. In summary, the
Dynaformer model offers a novel approach to predict protein-ligand binding
affinities by utilizing a roto-translation invariant encoding scheme for
atomic and structural features.

###  4.5 Training and hyperparameter settings

The training of Dynaformer consisted of two stages, i.e., pretraining and
finetuning. The model configuration and the pretraining hyperparameters are
summarized in Supplementary Table 2. We labeled the binding affinities of MD
trajectories according to their PDB IDs, which means that all snapshots of one
trajectory share the same label. The mean squared error (MSE) loss function
was used to pretrain the model. Different from the vanilla transformer, we
chose 4 layers with 512 hidden dimensions through grid searching. We also used
an adversarial graph data augmentation method, named FLAG, to mitigate the
potential over-fitting problem in the pretraining stage kong2020flag . For the
finetuning stage, Dynaformer was finetuned on the PDBBind general set. In this
stage, the mean absolute error (MAE) was used as a loss function as it
resulted in larger loss values. In addition, the peak learning rate was
adjusted to $10^{-5}$ and the batch size was set to 16. At both stages, the
datasets were split into a training and a validation subsets with a ratio of
9:1. The validation set was used to guide the training or finetuning process.
NVIDIA V100-32G GPU cards were used to train our model, and the training
process took about 24 hours.

###  4.6 Hit discovery details against HSP90

We prepared the HSP90 protein structure and the small molecules from a
compound library for docking. The protein structure and the corresponding
pocket position are determined from a reference structure (PDB ID: 2xdl) from
the RCSB PDB. We used PyMOL PyMOL  (version 2.5.0, open-source) to remove the
solvent molecules, ions and the original ligand. The average coordinate of the
original ligand, $x=65.4,y=35.0,z=25.5$, was utilized as the center of docking
box. Then we added hydrogens and charges to the protein, and then converted
the conformation representation into the PDBQT format, using ADFRsuite
ravindranath2015autodockfr  (version 1.0, build 5). For small molecules, we
used the ETKDG algorithm from the RDKit library landrum2013rdkit  (version
2022.3.3) to generate the initial conformations of molecules from the
Chembridge DIVERSet-EXP library. Up to 10 conformers were generated for each
molecule with a minimum heavy atom $\mbox{RMSD}=0.5$Å between conformers.
Here, redocking can help stride over the energy barrier and thus reach the
stable docking poses. All initial conformers were processed using the Meeko
Python library, including salt removal, hydrogen addition and partial charge
assignment. Finally, all initial conformers were converted to the PDBQT format
before docking. Then, Autodock Vina eberhardt2021autodock  (version 1.2.3) was
utilized for docking the protein and small molecule conformers. The docking
box size was chosen to be 25Å and the exhaustiveness of the global search was
set to 32. Autodock Vina was run in parallel on 3200 CPUs and all docking jobs
were finished within 20 hours. The best docking pose and the corresponding
docking energy of each molecule were kept. The pipeline of receptor
preparation, ligand preparation and docking is shown in Supplementary Figure
1.

Subsequently, the docked structures undergo binding affinity prediction using
Dynaformer within one hour. The top 10% of molecules from both Dynaformer
predictions and Autodock Vina docking energy were combined together. Following
the literature lyu2019ultra ; bender2021practical , we conducted a filtering
and visual inspection procedure. In particular, a filtering process according
to Veber’s rule and Lipinski’s rule of five lipinski1997experimental ;
veber2002molecular  was used. More specifically, we removed molecules that
have more than 5 hydrogen bond donors, more than 10 hydrogen bond acceptors, a
calculated octanol-water partition coefficient (Clog P) exceeding 5, more than
10 rotatable bonds, or a polar surface area larger than 140
$\mbox{\r{A}}^{2}$. Additionally, molecules with high intrinsic energy, such
as those containing 3-, 4-, 7-, or more than 12-membered rings, were excluded.
Then, visual inspection was performed based on the HSP90 ATP binding pocket
following the same principle as in stebbins1997crystal ; huth2007discovery ;
kung2008dihydroxylphenyl ; murray2010fragment ; janin2010atpase ;
yoshimura2021thermodynamic , and promising molecules were selected according
to their structures, that is, whether they possess a hydrophobic head and a
hydrophilic tail. Up to 20 molecules were selected for the subsequent wet lab
validation experiments (Supplementary Figure 2).

Finally, we experimentally validated the binding affinities of the above 20
selected molecules. Surface plasmon resonance (SPR) measurements were
performed by Biacore S200 instrument (GE Healthcare Life Sciences). HSP90
(9-236) was cross-linked on the surface of the Series S Sensor Chip CM5 (GE
Healthcare) at pH=4.0 by an amine coupling kit with
1-ethyl-(3-dimethylaminopropyl)-carbodiimide hydrochloride (EDC) and NHS.
Experiments were performed in a running buffer (1$\times$ PBS with 5% DMSO).
Each compound was serially diluted by a factor of two in the running buffer.
The compound solution flowed through the surface of the chip at a flow rate of
30 µL/min. Sensorgrams were analyzed and generated by Biacore Evaluation
Software. The inhibition constants $K_{i}$ of 8 out of 20 molecules could not
be determined, probably due to the poor aqueous solubility or non-specific
binding.

###  4.7 Evaluation Metrics

RMSF. The root mean squared fluctuation (RMSF) is the fluctuation around an
average, per atom or residue, over a sequence of structures such as a
trajectory. More specifically, it can be expressed as:

$\mbox{RMSF}=\sqrt{\langle(x_{i}-\overline{x_{i}})^{2}\rangle},$ (8)

where $x_{i}$ stands for the coordinates of the particle $i$, and
$\overline{x_{i}}$ stands for the trajectory-averaged position of $i$. While
the RMSD quantifies how much a structure diverges from a reference over time,
the RMSF can measure the mobility of parts of the system. Thus parts of the
structure with high RMSF values diverges from the average, indicating high
mobility.

Pearson $\bm{r}$. The Pearson correlation coefficient, also known as Pearson
$\bm{r}$, measures the linear relationship between two sequences of data
points $\mathbf{x}=\\{x_{i}\\}$ and $\mathbf{y}=\\{y_{i}\\}$. More
specifically, it can be defined as:

$\bm{r}(\mathbf{x},\mathbf{y})=\frac{\sum_{i=1}^{N}(x_{i}-\overline{x})(y_{i}-\overline{y})}{\sqrt{\sum_{i=1}^{N}(x_{i}-\overline{x})^{2}\sum_{i=1}^{N}(y_{i}-\overline{y})^{2}}},$
(9)

where $N$ stands for the number of samples, and $\overline{x},\overline{y}$
stand for the mean values of $\mathbf{x}=\\{x_{i}\\}$ and
$\mathbf{y}=\\{y_{i}\\}$, respectively. Pearson $\bm{r}$ varies between
$[-1,1]$, $0$ indicates no correlation and $-1$ or $1$ indicate the exact
negative or positive linear relationship. Note that Pearson $\bm{r}$ is
invariant to the scale or location of the two variables, which means that it
does not suggest the error between prediction and groundtruth.

Spearman $\bm{\rho}$. The Spearman correlation coefficient, denoted as
Spearman $\bm{\rho}$ is a nonparametric measure of the ranking order between
two sequences of data points. In other words, Spearman $\bm{\rho}$ has a
similar form as Pearson $\bm{r}$ but calculates the ranking values to assess
the monotonic relationships between variables. More specifically, it can be
defined as:

$\bm{\tau}(\mathbf{x},\mathbf{y})=\bm{r}(\mathbf{rx},\mathbf{ry})=\frac{\sum_{i=1}^{N}(rx_{i}-\overline{rx})(ry_{i}-\overline{ry})}{\sqrt{\sum_{i=1}^{N}(rx_{i}-\overline{rx})^{2}\sum_{i=1}^{N}(ry_{i}-\overline{ry})^{2}}},$
(10)

where $N$ stands for the number of samples, $\mathbf{rx}=\\{rx_{i}\\}$ and
$\mathbf{ry}=\\{ry_{i}\\}$ stand for the ranking orders of complexes between
the predicted and experimental binding affinities, respectively, and
$\overline{rx},\overline{ry}$ stand for the mean rank orders, respectively.
Spearman $\bm{\rho}$ is invariant to the specific values, but it is sensitive
to the order of variables.

Kendall $\bm{\tau}$. The Kendall correlation coefficient, denoted as Kendall
$\bm{\tau}$, is a measure of the ordinal association between two sequences of
ranking data points. More specifically, it can be defined as:

$\bm{\tau}(\mathbf{x},\mathbf{y})=\frac{N_{concord}-N_{discord}}{\sqrt{(N_{concord}+N_{discord}+\mathcal{T})(N_{concord}+N_{discord}+\mathcal{U})}},$
(11)

where $N_{concord}$ and $N_{discord}$ stand for the numbers of concordant and
discordant pairs, respectively. Here, given two pairs of prediction and
groundtruth $\\{(x_{1},y_{1}),(x_{2},y_{2})\\}$, if
$\mbox{sign}(x_{2}-x_{1})=\mbox{sign}(y_{2}-y_{1})$, this pair is called a
concordant pair. In contrast, if
$\mbox{sign}(x_{2}-x_{1})=-\mbox{sign}(y_{2}-y_{1})$, this pair is called a
discordant pair. A pair is tied if $x_{1}=x_{2}$ or $y_{1}=y_{2}$.
$\mathcal{T}$ and $\mathcal{U}$ are the numbers of tied pairs in $\mathbf{x}$
and $\mathbf{y}$, respectively. Kendall $\bm{\tau}$ is invariant to the
specific values, but it is sensitive to the order of variables.

PI. The predictive index (PI) is an evaluation metric mainly proposed for
virtual screening pearlman2001free . In particular, PI measures how well a
scoring function can correctly rank different ligands for a specific target by
taking the differences of binding affinities into account. More specifically,
it can be defined as:

$\mbox{PI}(\mathbf{x},\mathbf{y})=\frac{\sum_{j>i}^{N}\sum_{i=1}^{N}\omega_{ij}S_{ij}}{\sum_{j>i}^{N}\sum_{i=1}^{N}\omega_{ij}},$
(12)

where $\omega_{ij}$ stands for the absolute difference of binding affinities
between ligands $i$ and $j$, i.e.,
$\omega_{ij}=|\mbox{p}K_{i}^{(i)}-\mbox{p}K_{i}^{(j)}|=|y_{i}-y_{j}|$,
$S_{ij}\in\\{1,0\\}$ indicates if the ranking orders of binding affinities
$x_{i}$ and $x_{j}$ are consistent between prediction and groundtruth.

##  5 Data availability

All MD trajectories, scripts for MD simulation and preprocessed data for
Dynaformer are available at
https://1drv.ms/f/s!Ah9r82oejjV8piQHq_qAieio_86z?e=B1E53d.

##  6 Code availability

The Dynoformer is implemented in Python and the source code, trained models as
well as evaluation results are available at
https://github.com/Minys233/Dynaformer.

Supplementary information The supplementary material is available in the
attachment.

Acknowledgements We thank Hao Zhang and Tristan Bereau for the insightful
discussions.

Author contributions Y.W., D.Z. and J.Z conceived the idea; J.Z. and J.W.
planned the dry lab and wet lab experiments; Y.W., P.W. and Y.M. performed MD
simulations and quality assessment; Y.M., Y.W., P.W., J.Z., S.Z and J.W.
designed the Dynaformer model; Y.M., Y.S. and Y.H.W trained and tested the
Dynaformer model; Y.M. conducted virtual screening; J.Z., D.Z and X.W.
directed SPR assays; Y.M., Y.W., J.Z., D.Z. and J.W. analysised the data;
Y.M., P.W., Y.W., H.L. and X.W. wrote the original drafts and all authors
commented on the manuscript.

Competing interests Authors declare that they have no competing interests.

Figure 1: Illustrative overview of the thermodynamic ensemble and the MD
trajectory dataset. A The conformations in a thermodynamic ensemble generally
follow a specific distribution. In principle, the binding affinity is
determined by the properties of the thermodynamic ensemble rather than a
single conformation. Based on this fact, we built a dataset containing various
thermodynamic ensembles sampled by MD simulations. More details can be found
in the main text. B Relationship between the mean ligand RMSD and the binding
affinity in the MD trajectory dataset. The mean ligand RMSD is calculated as
the average distance deviation of the ligand from its original position in the
crystal structure during the simulation. The heatmap and the mean RMSD range
in the inset indicate that the ligand dynamics are related to binding
affinities, i.e., the higher the mean ligand RMSD, the lower the binding
affinity. C Representative examples of ligand RMSD and their binding
affinities, including unstable (PDB ID: 4y3j), intermediate (PDB ID: 3udh),
and stable (PDB ID: 2yge) MD trajectories. For the unstable case, the ligand
RMSD increased sharply at 60th snapshot, indicating that the ligand left the
binding site. The ligand with intermediate stability remained relatively
flexible but stayed near the binding site. In the stable case, the ligand
stuck tightly to the binding pocket, yielding low ligand RMSDs. Figure 2:
Illustration of our virtual screening pipeline with Dynaformer. A Hit
discovery through our deep learning model Dynaformer based on an MD trajectory
dataset. First, crystal structures of protein-ligand complexes are collected
to form an MD trajectory dataset. Then, the simulated trajectories are
utilized for training the binding affinity prediction model. Next, the trained
model is evaluated against a benchmark dataset. The model with the best
performance is then used as a scoring function for the virtual screen of the
compound candidates. Finally, the wet-lab experiment validates the hit
molecules for further lead optimization. B Example snapshot of protein-ligand
from the MD trajectory dataset and its graph representation. In order to
create such graph representation to feed into the model, atoms from both the
ligand and the protein are chosen based on their proximity to the ligand,
within a specified distance cutoff. The nodes and the edges represent the
atoms and their covalent or non-covalent interactions, respectively. C The
architecture of Dynaformer. The node features (i.e., atomic features) are
first fed into a multi-head attention module. The structural encodings,
including distance, angle, and edge features, are encoded and added as the
attention bias. After the final layer, the feature representations of the
graph are globally pooled and fused with the pre-calculated fingerprints,
which are knowledge-based features capturing structural and chemical
properties. The fused representations are then used for the final prediction
of the binding affinity. Abbreviations in the figure: Linear, linear layer;
Feed-forward, feed-forward neural network; MatMul, matrix multiplication.
Figure 3: Performance evaluation of binding affinity prediction. A Scoring
power evaluation on the 285 CASF-2016 complexes. The scoring power illustrates
the inter-target performances of how well the predictions correlate with the
experimentally measured binding affinities. Pearson coefficient (Pearson
$\bm{r}$), root mean squared error (RMSE), and standard deviation (SD) are
used to as evaluation metrics. B Ranking power evaluation on the 57 classes of
target proteins from CASF-2016, in which each target bound to five different
ligands. The goal of evaluating the ranking power is to assess the intra-
target performances of scoring functions in terms of how accurately they
predict the relative order of molecules. Evaluation metrics are measured in
terms of the average values of Spearman’s coefficient (Spearman $\rho$),
Kendall’s coefficient (Kendall $\tau$), and predictive index (PI) on each
target. C Detailed scatter plot of the Dynaformer predicted $\mbox{p}K_{i}$
and experimentally measured $\mbox{p}K_{i}$ values of the 285 protein-ligand
complexes from CASF-2016, demonstrating the effectiveness of Dynaformer in
predicting binding affinity with high correlation and low prediction bias. D
Ablation studies of Dynaformer. Here, -F indicates a modified version of
Dynaformer without any pre-calculated fingerprints; -3D stands for a modified
version of Dynaformer without the structural encoding module; -M stands for a
modified version of Dynaformer that is trained without MD trajectories, i.e.,
only crystal structures are used for training; The results showed the
consistent improvements in both scoring and ranking power by incorporating MD
data. Figure 4: Illustrative examples of Dynaformer trained with MD
trajectories improves binding affinity prediction. A An example of improved
affinity prediction through incorporating additional entropic information. The
head of the ligand is tightly bound in the pocket, resulting in a significant
favorable enthalpy change. The random movement of the solvent-exposed ligand
tail contributes to a favorable entropy change during ligand binding. B An
example of improved affinity prediction through better protein-ligand
interaction modeling. With a rigid structure and no rotatable bonds, the
binding affinity relies primarily on the enthalpy changes from binding site
interactions. Despite a few favorable interactions, subpockets around the
ligand remain unoccupied. Therefore, the ligand moves away from its initial
position during the simulation, leading to low binding affinity. The ability
of Dynaformer to model these interactions results in better prediction. C An
example of ranking ligands with subtle interaction differences. Dynaformer
addresses activity cliffs by differentiating between ligands sharing the same
scaffold but with varying binding affinities. D Conformational stability of
the three ligands shown in C. The 2qbp ligand exhibited high stability during
the simulation, while the 2qbq ligand partially left the original binding site
after 1 ns, and the 2qbr ligand left the binding site after 9 ns. The varying
binding affinities of these ligands can be attributed to the hydrophobic
interactions with the second subpocket on the right side shown in C. E The
interaction fraction analysis of the three ligands shown in C. The percentages
of simulation snapshots with different protein-ligand interactions are shown.
Higher binding affinity corresponded to more stable interactions observed
during the simulation. With such interaction features and better modeling
capability, Dynaformer can provide more accurate ligand rankings. Figure 5:
Virtual screening against the target HSP90 and the corresponding experimental
validation. A Structures, experimentally measured $K_{i}$ values and the SPR
sensorgrams of the top 12 molecules with measurable binding affinities from
the virtual screening. B Detailed $K_{i}$ values of the selected 12 molecules
and their ranking statistics by Dynaformer and Autodock Vina. To best of our
knowledge, Compounds 1, 3, 4, 6 and 12 show promise as novel candidates for
hit-to-lead optimization. C The correlations of the predictions by Dynaformer
and Autodock Vina predictions versus experimental $K_{i}$ values of the 12
compounds. D The possible binding modes of the top three compounds, in which
the interacting residues are consistent with the previous studies. Overall,
the wet-lab validations demonstrated that Dynaformer can effectively identify
hit compounds with favorable binding affinities and novel scaffolds.

## References

  *   * (1) Rube, H.T., Rastogi, C., Feng, S., Kribelbauer, J.F., Li, A., Becerra, B., Melo, L.A., Do, B.V., Li, X., Adam, H.H., et al.: Prediction of protein–ligand binding affinity from sequencing data with interpretable machine learning. Nature Biotechnology 40(10), 1520–1527 (2022)
  * (2) McFedries, A., Schwaid, A., Saghatelian, A.: Methods for the elucidation of protein-small molecule interactions. Chemistry & biology 20(5), 667–673 (2013)
  * (3) Swinney, D.C., Anthony, J.: How were new medicines discovered? Nature reviews Drug discovery 10(7), 507–519 (2011)
  * (4) Erlanson, D.A., Fesik, S.W., Hubbard, R.E., Jahnke, W., Jhoti, H.: Twenty years on: the impact of fragments on drug discovery. Nature reviews Drug discovery 15(9), 605–619 (2016)
  * (5) Mobley, D.L., Gilson, M.K.: Predicting binding free energies: frontiers and benchmarks. Annual review of biophysics 46, 531–558 (2017)
  * (6) Kumar, V.G., Agrawal, S., Suresh Kumar, T.K., Moradi, M.: Binding affinity estimation from restrained umbrella sampling simulations. bioRxiv, 2021–10 (2021)
  * (7) Fischer, A., Smiesko, M., Sellner, M., Lill, M.A.: Decision making in structure-based drug discovery: visual inspection of docking results. Journal of Medicinal Chemistry 64(5), 2489–2500 (2021)
  * (8) Frenkel, D., Smit, B.: Understanding Molecular Simulation: from Algorithms to Applications, Second edition edn., pp. 167–200. Academic Press, San Diego (2002)
  * (9) Valdés-Tresanco, M.S., Valdés-Tresanco, M.E., Valiente, P.A., Moreno, E.: gmx_mmpbsa: a new tool to perform end-state free energy calculations with gromacs. Journal of chemical theory and computation 17(10), 6281–6291 (2021)
  * (10) Song, L.F., Merz Jr, K.M.: Evolution of alchemical free energy methods in drug discovery. Journal of Chemical Information and Modeling 60(11), 5308–5318 (2020)
  * (11) Liu, T., Lu, D., Zhang, H., Zheng, M., Yang, H., Xu, Y., Luo, C., Zhu, W., Yu, K., Jiang, H.: Applying high-performance computing in drug discovery and molecular simulation. National science review 3(1), 49–63 (2016)
  * (12) Lo, Y.-C., Rensi, S.E., Torng, W., Altman, R.B.: Machine learning in chemoinformatics and drug discovery. Drug discovery today 23(8), 1538–1546 (2018)
  * (13) Zhu, H., Yang, J., Huang, N.: Assessment of the generalization abilities of machine-learning scoring functions for structure-based virtual screening. Journal of Chemical Information and Modeling 62(22), 5485–5502 (2022)
  * (14) Verdonk, M.L., Cole, J.C., Hartshorn, M.J., Murray, C.W., Taylor, R.D.: Improved protein–ligand docking using gold. Proteins: Structure, Function, and Bioinformatics 52(4), 609–623 (2003)
  * (15) Friesner, R.A., Banks, J.L., Murphy, R.B., Halgren, T.A., Klicic, J.J., Mainz, D.T., Repasky, M.P., Knoll, E.H., Shelley, M., Perry, J.K., et al.: Glide: a new approach for rapid, accurate docking and scoring. 1. method and assessment of docking accuracy. Journal of medicinal chemistry 47(7), 1739–1749 (2004)
  * (16) Trott, O., Olson, A.J.: 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)
  * (17) Korb, O., Stutzle, T., Exner, T.E.: Empirical scoring functions for advanced protein-ligand docking with plants. Journal of chemical information and modeling 49(1), 84–96 (2009)
  * (18) Dittrich, J., Schmidt, D., Pfleger, C., Gohlke, H.: Converging a knowledge-based scoring function: Drugscore2018. Journal of chemical information and modeling 59(1), 509–521 (2018)
  * (19) Mooij, W.T., Verdonk, M.L.: General and targeted statistical potentials for protein–ligand interactions. Proteins: Structure, Function, and Bioinformatics 61(2), 272–287 (2005)
  * (20) Su, M., Yang, Q., Du, Y., Feng, G., Liu, Z., Li, Y., Wang, R.: Comparative assessment of scoring functions: the casf-2016 update. Journal of chemical information and modeling 59(2), 895–913 (2018)
  * (21) Wang, R., Lai, L., Wang, S.: Further development and validation of empirical scoring functions for structure-based binding affinity prediction. Journal of computer-aided molecular design 16(1), 11–26 (2002)
  * (22) Ballester, P.J., Mitchell, J.B.: A machine learning approach to predicting protein–ligand binding affinity with applications to molecular docking. Bioinformatics 26(9), 1169–1175 (2010)
  * (23) Wang, C., Zhang, Y.: Improving scoring-docking-screening powers of protein–ligand scoring functions using random forest. Journal of computational chemistry 38(3), 169–177 (2017)
  * (24) Sánchez-Cruz, N., Medina-Franco, J.L., Mestres, J., Barril, X.: Extended connectivity interaction features: improving binding affinity prediction through chemical description. Bioinformatics 37(10), 1376–1382 (2021)
  * (25) Stepniewska-Dziubinska, M.M., Zielenkiewicz, P., Siedlecki, P.: Development and evaluation of a deep learning model for protein–ligand binding affinity prediction. Bioinformatics 34(21), 3666–3674 (2018)
  * (26) Jiménez, J., Skalic, M., Martinez-Rosell, G., De Fabritiis, G.: Kdeep: protein-ligand absolute binding affinity prediction via 3d-convolutional neural networks. Journal of chemical information and modeling 58(2), 287–296 (2018)
  * (27) Zheng, L., Fan, J., Mu, Y.: Onionnet: a multiple-layer intermolecular-contact-based convolutional neural network for protein–ligand binding affinity prediction. ACS omega 4(14), 15956–15965 (2019)
  * (28) Nguyen, T., Le, H., Quinn, T.P., Nguyen, T., Le, T.D., Venkatesh, S.: Graphdta: Predicting drug–target binding affinity with graph neural networks. Bioinformatics 37(8), 1140–1147 (2021)
  * (29) Yang, K., Swanson, K., Jin, W., Coley, C., Eiden, P., Gao, H., Guzman-Perez, A., Hopper, T., Kelley, B., Mathea, M., et al.: Analyzing learned molecular representations for property prediction. Journal of chemical information and modeling 59(8), 3370–3388 (2019)
  * (30) Lim, J., Ryu, S., Park, K., Choe, Y.J., Ham, J., Kim, W.Y.: Predicting drug–target interaction using a novel graph neural network with 3d structure-embedded graph representation. Journal of chemical information and modeling 59(9), 3981–3988 (2019)
  * (31) Maziarka, L., Danel, T., Mucha, S., Rataj, K., Tabor, J., Jastrzebski, S.: Molecule attention transformer. arXiv preprint arXiv:2002.08264 (2020)
  * (32) Klicpera, J., Groß, J., Günnemann, S.: Directional message passing for molecular graphs. arXiv preprint arXiv:2003.03123 (2020)
  * (33) Song, Y., Zheng, S., Niu, Z., Fu, Z.-H., Lu, Y., Yang, Y.: Communicative representation learning on attributed molecular graphs. In: IJCAI, vol. 2020, pp. 2831–2838 (2020)
  * (34) Li, S., Zhou, J., Xu, T., Huang, L., Wang, F., Xiong, H., Huang, W., Dou, D., Xiong, H.: Structure-aware interactive graph neural networks for the prediction of protein-ligand binding affinity. In: Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 975–985 (2021)
  * (35) Wang, R., Fang, X., Lu, Y., Wang, S.: The pdbbind database: Collection of binding affinities for protein- ligand complexes with known three-dimensional structures. Journal of medicinal chemistry 47(12), 2977–2980 (2004)
  * (36) Laio, A., Gervasio, F.L.: Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Reports on Progress in Physics 71(12), 126601 (2008)
  * (37) Christ, C.D., Mark, A.E., Van Gunsteren, W.F.: Basic ingredients of free energy calculations: a review. Journal of computational chemistry 31(8), 1569–1582 (2010)
  * (38) Wei, G., Xi, W., Nussinov, R., Ma, B.: Protein ensembles: how does nature harness thermodynamic fluctuations for life? the diverse functional roles of conformational ensembles in the cell. Chemical reviews 116(11), 6516–6551 (2016)
  * (39) Berman, H.M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T.N., Weissig, H., Shindyalov, I.N., Bourne, P.E.: The protein data bank. Nucleic acids research 28(1), 235–242 (2000)
  * (40) Ying, C., Cai, T., Luo, S., Zheng, S., Ke, G., He, D., Shen, Y., Liu, T.-Y.: Do transformers really perform badly for graph representation? Advances in Neural Information Processing Systems 34, 28877–28888 (2021)
  * (41) Bartók, A.P., Payne, M.C., Kondor, R., Csányi, G.: Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. Physical review letters 104(13), 136403 (2010)
  * (42) Rayka, M., Firouzi, R.: Gb-score: Minimally designed machine learning scoring function based on distance-weighted interatomic contact features. Molecular Informatics (2023)
  * (43) Modugno, M., Casale, E., Soncini, C., Rosettani, P., Colombo, R., Lupi, R., Rusconi, L., Fancelli, D., Carpinelli, P., Cameron, A.D., et al.: Crystal structure of the t315i abl mutant in complex with the aurora kinases inhibitor pha-739358. Cancer research 67(17), 7987–7990 (2007)
  * (44) Efremov, I.V., Vajdos, F.F., Borzilleri, K.A., Capetta, S., Chen, H., Dorff, P.H., Dutra, J.K., Goldstein, S.W., Mansour, M., McColl, A., et al.: Discovery and optimization of a novel spiropyrrolidine inhibitor of $\beta$-secretase (bace1) through fragment-based drug design. Journal of medicinal chemistry 55(21), 9069–9088 (2012)
  * (45) Wilson, D.P., Wan, Z.-K., Xu, W.-X., Kirincich, S.J., Follows, B.C., Joseph-McCarthy, D., Foreman, K., Moretto, A., Wu, J., Zhu, M., et al.: Structure-based optimization of protein tyrosine phosphatase 1b inhibitors: from the active site to the second phosphotyrosine binding site. Journal of medicinal chemistry 50(19), 4681–4698 (2007)
  * (46) Lyu, J., Wang, S., Balius, T.E., Singh, I., Levit, A., Moroz, Y.S., O’Meara, M.J., Che, T., Algaa, E., Tolmachova, K., et al.: Ultra-large library docking for discovering new chemotypes. Nature 566(7743), 224–229 (2019)
  * (47) Janin, Y.L.: Atpase inhibitors of heat-shock protein 90, second season. Drug discovery today 15(9-10), 342–353 (2010)
  * (48) Li, L., Wang, L., You, Q.-D., Xu, X.-L.: Heat shock protein 90 inhibitors: an update on achievements, challenges, and future directions. Journal of medicinal chemistry 63(5), 1798–1822 (2019)
  * (49) Eberhardt, J., Santos-Martins, D., Tillack, A.F., Forli, S.: Autodock vina 1.2. 0: New docking methods, expanded force field, and python bindings. Journal of Chemical Information and Modeling 61(8), 3891–3898 (2021)
  * (50) Bender, B.J., Gahbauer, S., Luttens, A., Lyu, J., Webb, C.M., Stein, R.M., Fink, E.A., Balius, T.E., Carlsson, J., Irwin, J.J., et al.: A practical guide to large-scale docking. Nature protocols 16(10), 4799–4832 (2021)
  * (51) Kung, P.-P., Funk, L., Meng, J., Collins, M., Zhou, J.Z., Johnson, M.C., Ekker, A., Wang, J., Mehta, P., Yin, M.-J., et al.: Dihydroxylphenyl amides as inhibitors of the hsp90 molecular chaperone. Bioorganic & medicinal chemistry letters 18(23), 6273–6278 (2008)
  * (52) Murray, C.W., Carr, M.G., Callaghan, O., Chessari, G., Congreve, M., Cowan, S., Coyle, J.E., Downham, R., Figueroa, E., Frederickson, M., et al.: Fragment-based drug discovery applied to hsp90. discovery of two lead series with high ligand efficiency. Journal of medicinal chemistry 53(16), 5942–5955 (2010)
  * (53) Stebbins, C.E., Russo, A.A., Schneider, C., Rosen, N., Hartl, F.U., Pavletich, N.P.: Crystal structure of an hsp90–geldanamycin complex: targeting of a protein chaperone by an antitumor agent. Cell 89(2), 239–250 (1997)
  * (54) Yoshimura, C., Nagatoishi, S., Kuroda, D., Kodama, Y., Uno, T., Kitade, M., Chong-Takata, K., Oshiumi, H., Muraoka, H., Yamashita, S., et al.: Thermodynamic dissection of potency and selectivity of cytosolic hsp90 inhibitors. Journal of Medicinal Chemistry 64(5), 2669–2677 (2021)
  * (55) Huth, J.R., Park, C., Petros, A.M., Kunzer, A.R., Wendt, M.D., Wang, X., Lynch, C.L., Mack, J.C., Swift, K.M., Judge, R.A., et al.: Discovery and design of novel hsp90 inhibitors using multiple fragment-based design strategies. Chemical biology & drug design 70(1), 1–12 (2007)
  * (56) Duan, L., Liu, X., Zhang, J.Z.: Interaction entropy: A new paradigm for highly efficient and reliable computation of protein–ligand binding free energy. Journal of the American Chemical Society 138(17), 5722–5728 (2016)
  * (57) Sun, Z., Yan, Y.N., Yang, M., Zhang, J.Z.: Interaction entropy for protein-protein binding. The Journal of chemical physics 146(12), 124124 (2017)
  * (58) Huang, H., Simmerling, C.: Fast pairwise approximation of solvent accessible surface area for implicit solvent simulations of proteins on cpus and gpus. Journal of chemical theory and computation 14(11), 5797–5814 (2018)
  * (59) Pratt, L.R., Chandler, D.: Theory of the hydrophobic effect. The Journal of chemical physics 67(8), 3683–3704 (1977)
  * (60) Massova, I., Kollman, P.A.: Combined molecular mechanical and continuum solvent approach (mm-pbsa/gbsa) to predict ligand binding. Perspectives in drug discovery and design 18, 113–135 (2000)
  * (61) Liu, Z., Su, M., Han, L., Liu, J., Yang, Q., Li, Y., Wang, R.: Forging the basis for developing protein–ligand interaction scoring functions. Accounts of chemical research 50(2), 302–309 (2017)
  * (62) Jo, S., Kim, T., Iyer, V.G., Im, W.: Charmm-gui: a web-based graphical user interface for charmm. Journal of computational chemistry 29(11), 1859–1865 (2008)
  * (63) Lee, J., Cheng, X., Swails, J.M., Yeom, M.S., Eastman, P.K., Lemkul, J.A., Wei, S., Buckner, J., Jeong, J.C., Qi, Y., et al.: Charmm-gui input generator for namd, gromacs, amber, openmm, and charmm/openmm simulations using the charmm36 additive force field. Journal of chemical theory and computation 12(1), 405–413 (2016)
  * (64) Phillips, J.C., Hardy, D.J., Maia, J.D., Stone, J.E., Ribeiro, J.V., Bernardi, R.C., Buch, R., Fiorin, G., Hénin, J., Jiang, W., et al.: Scalable molecular dynamics on cpu and gpu architectures with namd. The Journal of chemical physics 153(4), 044130 (2020)
  * (65) Ghosh, J., Marru, S., Singh, N., Vanomesslaeghe, K., Fan, Y., Pamidighantam, S.: Molecular parameter optimization gateway (paramchem) workflow management through teragrid asta. In: Proceedings of the 2011 TeraGrid Conference: Extreme Digital Discovery, pp. 1–8 (2011)
  * (66) Coutsias, E.A., Seok, C., Jacobson, M.P., Dill, K.A.: A kinematic view of loop closure. Journal of computational chemistry 25(4), 510–528 (2004)
  * (67) Genheden, S., Ryde, U.: How to obtain statistically converged mm/gbsa results. Journal of computational chemistry 31(4), 837–846 (2010)
  * (68) Hu, W., Fey, M., Zitnik, M., Dong, Y., Ren, H., Liu, B., Catasta, M., Leskovec, J.: Open graph benchmark: Datasets for machine learning on graphs. Advances in neural information processing systems 33, 22118–22133 (2020)
  * (69) Michaud-Agrawal, N., Denning, E.J., Woolf, T.B., Beckstein, O.: Mdanalysis: a toolkit for the analysis of molecular dynamics simulations. Journal of computational chemistry 32(10), 2319–2327 (2011)
  * (70) Landrum, G., et al.: Rdkit: A software suite for cheminformatics, computational chemistry, and predictive modeling. Greg Landrum (2013)
  * (71) Cisneros, G.A.: Application of gaussian electrostatic model (gem) distributed multipoles in the amoeba force field. Journal of Chemical Theory and Computation 8(12), 5072–5080 (2012)
  * (72) Kong, K., Li, G., Ding, M., Wu, Z., Zhu, C., Ghanem, B., Taylor, G., Goldstein, T.: Flag: Adversarial data augmentation for graph neural networks. arXiv preprint arXiv:2010.09891 (2020)
  * (73) Schrödinger, LLC: The PyMOL Molecular Graphics System, Version 1.8 (2023)
  * (74) Ravindranath, P.A., Forli, S., Goodsell, D.S., Olson, A.J., Sanner, M.F.: Autodockfr: advances in protein-ligand docking with explicitly specified binding site flexibility. PLoS computational biology 11(12), 1004586 (2015)
  * (75) Lipinski, C.A., Lombardo, F., Dominy, B.W., Feeney, P.J.: Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Advanced drug delivery reviews 23(1-3), 3–25 (1997)
  * (76) Veber, D.F., Johnson, S.R., Cheng, H.-Y., Smith, B.R., Ward, K.W., Kopple, K.D.: Molecular properties that influence the oral bioavailability of drug candidates. Journal of medicinal chemistry 45(12), 2615–2623 (2002)
  * (77) Pearlman, D.A., Charifson, P.S.: Are free energy calculations useful in practice? a comparison with rapid scoring functions for the p38 map kinase protein system. Journal of medicinal chemistry 44(21), 3417–3423 (2001)