#  PharmacoNet: Accelerating Large-Scale
Virtual Screening by Deep Pharmacophore Modeling

Seonghwan Seo
Department of Chemistry, KAIST
<EMAIL_ADDRESS>
&Woo Youn Kim
Department of Chemistry, KAIST
AI Institute, KAIST
HITS Incorporation
<EMAIL_ADDRESS>

###### Abstract

As the size of accessible compound libraries expands to over 10 billion, the
need for more efficient structure-based virtual screening methods is emerging.
Different pre-screening methods have been developed to rapidly screen the
library, but the structure-based methods applicable to general proteins are
still lacking: the challenge is to predict the binding pose between proteins
and ligands and perform scoring in an extremely short time. We introduce
PharmacoNet, a deep learning framework that identifies the optimal 3D
pharmacophore arrangement which a ligand should have for stable binding from
the binding site. By coarse-grained graph matching between ligands and the
generated pharmacophore arrangement, we solve the expensive binding pose
sampling and scoring procedures of existing methods in a single step.
PharmacoNet is significantly faster than state-of-the-art structure-based
approaches, yet reasonably accurate with a simple scoring function.
Furthermore, we show the promising result that PharmacoNet effectively retains
hit candidates even under the high pre-screening filtration rates. Overall,
our study uncovers the hitherto untapped potential of a pharmacophore modeling
approach in deep learning-based drug discovery.

##  1 Introduction

In the realm of modern drug discovery, the search for novel and effective drug
molecules has necessitated the exploration of vast chemical spaces. Structure-
based virtual screening with ultra-large chemical libraries has attracted
great attention as a new paradigm for early drug discovery as make-on-demand
libraries of drug-like molecules have recently expanded into the tens of
billions [1, 2, 3, 4, 5]. Recent work reported that increasing the chemical
library size from 115 million to 11 billion significantly improved the hit
rate from 15% to 30% [6].

Molecular docking is a common approach for structure-based virtual screening.
Conventional docking methods involve two steps: generating a large number of
binding poses for each initial ligand conformation and ranking binding poses
using a well-defined scoring function [7, 8, 9]. As the scoring functions are
based on the physical energy equation, the highest score among the binding
poses – docking score – has been considered as a criterion for virtual
screening. Nevertheless, the pose space is large, continuous, and contains
numerous local optima, requiring binding pose search processes spanning from 1
to 100 seconds for each initial conformation [10, 11, 12]. Considering the
virtual screening uses over 100 initial conformations for accurate evaluation,
the computational demands of docking become a limiting factor in exploring a
wider screening space of over a billion molecules [1, 5]. As a result, the
expansion of chemical libraries has significantly increased the cost of
virtual screening, highlighting the need for more efficient screening
approaches.

Recent works [5, 6, 13] have focused on how to effectively select molecules
for docking instead of brute force docking. Sadybekov et al. [6] used a
structured library Enamine REAL [14] to perform hierarchical enumeration
screening, and Graff et al. [13] employed active learning to efficiently
explore molecular libraries. However, these approaches have the disadvantage
that they do not perform any evaluation for most of the molecules in the
library. To overcome this limitation, a two-step strategy involving pre-
screening&fine-screening can be employed in virtual screening procedures. Pre-
screening involves the evaluation of whole molecules in a library using a fast
scoring function, followed by fine-screening of selected molecules using
computationally intensive methods such as molecular docking or molecular
dynamics simulations.

To this end, we introduce PharmacoNet, the first deep learning framework for
structure-based pharmacophore modeling, – Figure 1. A pharmacophore111IUPAC
definition: ”an ensemble of steric and electronic features that is necessary
to ensure the optimal supramolecular interactions with a specific biological
target and to trigger (or block) its biological response” [15]222In this
paper, we use hydrogen-bond donor & acceptor, aromatic ring, cation, anion,
hydrophobic carbon, and halogen atom is an abstract description of molecular
features with the potential to form non-covalent interactions333e.g.
hydrophobic interaction, hydrogen bond, salt-bridge, $\pi$–$\pi$ stacking,
cation–$\pi$ interaction, halogen bond critical for stabilizing protein-ligand
complexes, such as aromatic rings and hydrogen bond (H-bond) donors &
acceptors. Pharmacophore modeling identifies the optimal 3D pharmacophore
arrangement that a ligand should have to ensure optimal supramolecular
interactions[15, 16, 17, 18]. Rather than energy-based scoring by modeling
pairwise interactions between ligand and protein atoms for a large number of
possible binding poses, our framework uses scoring based on how well a ligand
fits into the pharmacophore model, enabling extremely rapid yet reasonably
accurate scoring. To achieve this objective, PharmacoNet:

  1. 1.

identifies the desirable 3D pharmacophore arrangement (pharmacophore model)
that a ligand should have for stable binding to the binding site.

  2. 2.

performs coarse-grained graph matching to solve expensive binding pose
generation and scoring procedures into a single step: implicit binding pose
prediction and scoring.

  3. 3.

abstracts the ligand into a pharmacophore-level graph for reasonable accuracy
and high generalization ability without energy scoring.

Figure 1:  Overview of PharmacoNet. Pharmacophore modeling is performed only
once during the virtual screening process.

Our novel approach incorporates prior knowledge concerning pharmacological
properties and non-covalent molecular interactions as an inductive bias into
the deep learning model. We simplify the topology of complex structures by
coarse-grained modeling at the pharmacophore level and incorporating distance
constraints for non-covalent interactions. This effectively addresses the
challenge of model generalization to proteins not encountered during training,
a common issue in structure-based models due to the limited structure data
relative to the topological complexity of ligand-protein complexes. As a
result, our deep learning model performs well even when trained on a limited
dataset comprising 5,500 complex crystal structures, accounting for 30% of the
PDBBind v2020 dataset [19].

Figure 2:  Model architecture of PharmacoNet. For visualization, 3D feature
maps are represented as 2D maps.

##  2 PharmacoNet Architecture

The high-level protocol of our framework is described in Figure 1. The model
architecture used in PharmacoNet for pharmacophore modeling is shown in Figure
2, Section 2.1, and Appendix 5.2. The graph matching algorithm and scoring
function are explained in Section 2.2 and Appendix 5.4. The code is available
at https://github.com/SeonghwanSeo/PharmacoNet.

###  2.1 Instance Image Segmentation for Pharmacophore Modeling

Similar to our approach, Skalic et al [20] developed LigVoxel, the 3D CNN-
based model that performs inpainting the chemical functionality map from the
pocket image. However, LigVoxel generates only one map for each of three
functional types: aromatic, H-bond donor, and H-bond acceptor, so it does not
recognize individual pharmacophoric features or ensure supramolecular
interactions with the binding site.

In constructing a pharmacophore model, it is essential to include the nature
of each pharmacophore: the appropriate coordinates, counterpart protein
functional groups (or residues), and the supramolecular interactions between
them [15, 21]. To facilitate the detection of individual pharmacophores and
their inherent characteristics, we frame the pharmacophore modeling as an
instance image segmentation task – Figure 2. The details of the model
architecture are in Appendix 5.2.

Compared to the original image segmentation task, the main differences are:

  1. 1.

Instead of detecting instances from images, our model selects hot spots
(instances) among protein functional groups (tokens) in the binding site. If a
single functional group can form multiple types of non-covalent interactions,
we consider it as multiple tokens.

  2. 2.

Each instance contains a non-covalent interaction type (class) and
coordinates. Based on prior knowledge about the maximum length for each
interaction type, we can set the possible area of the corresponding
pharmacophore (bounding box) as constraints. Therefore, the model uses class
and bounding box information for binary mask classification instead of
predicting the class or the bounding box of the instance.

  3. 3.

A single voxel can be included in multiple instances.

  4. 4.

The model also addresses an image inpainting problem to segment an empty space
in the binding site. This means that there is no definitive answer for the
segmentation map, so we cannot use common evaluation metrics used in image
segmentation tasks, such as mIoU (mean Intersection over Union).

###  2.2 Graph Matching and Scoring

The pharmacophore model and the ligand pharmacophore arrangement can be
represented as 3D complete graphs of pharmacophores, and the similarity
evaluation between them can be framed as an inexact graph-matching task. Since
inexact graph matching allows for some tolerance between the correspondence
between nodes and edges in two graphs, there are many possible matches, unlike
exact graph matching. To select the best match, a scoring function that
reflects the number of matched nodes and the error in the length of edges is
required. Finally, inexact graph matching can be performed between a
pharmacophore model and the ligand pharmacophore arrangement based on the
scoring function and the score of the best match means the graph similarity.
This is conceptually the same as the principle of conventional molecular
docking. Details about the graph matching algorithm are in Appendix 5.4.

Recently, Shen et al. [22] reported the state-of-the-art scoring model based
on the pairwise atom distance likelihood between protein and ligand instead of
scoring in the unit of energy. This has the advantage of allowing relative
comparisons between ligands without the need for mapping to experimentally
measured affinity values. Similarly, we propose a distance likelihood-based
scoring function to evaluate the graph match between a pharmacophore model and
pharmacophore arrangements of the ligand.

##### Ligand pharmacophore arrangement.

The complete graph of ligand pharmacophore arrangement is denoted as
$\displaystyle{\mathcal{G}}^{L}$, where each node $\displaystyle
v^{l}_{i}\in{\mathcal{G}}^{L}$ represents the ligand pharmacophore and it
contains the set of its possible pharmacophore types $T^{l}_{i}$. For ligand
conformation $\displaystyle{\bm{C}}$, the position of node $v^{l}_{i}$ is
$\displaystyle{\bm{x}}^{l}_{i}({\bm{C}})$.

##### Pharamcophore model.

The complete graph of the pharmacophore model $\displaystyle{\mathcal{G}}^{P}$
uses pharmacophore as nodes, and each node $\displaystyle
v^{p}_{i}\in{\mathcal{G}}^{P}$ has a position
$\displaystyle{\bm{x}}^{p}_{i}\in\mathbb{R}^{3}$, a radius $\displaystyle
r^{p}_{i}\in\mathbb{R}$, a token score $s^{p}_{i}$ (see Appendix 5.2) and a
pharmacophore type $t^{p}_{i}$. Since one ligand pharmacophore can form
interactions with multiple protein functional groups, i.e. one $v^{l}$ can be
matched with multiple $v^{p}$, we perform clustering for the pharmacophore
model. (see Algorithm 1 in Appendix 5.4) The set of clusters is
$\displaystyle{\mathbb{C}}^{P}$ where each cluster $\displaystyle
C^{p}_{i}\in{\mathbb{C}}^{P}$ is a set of $v^{p}_{i^{\prime}}$. During the
graph matching, each ligand node $v^{l}$ corresponds to a cluster $C^{p}$
instead of a single pharmacophore model node $v^{p}$.

##### Graph match and scoring.

For the total set of possible graph matches $\displaystyle{\mathbb{M}}$, each
graph match $M_{i}\in\displaystyle{\mathbb{M}}$ is
$\\{(v^{l}_{1},I^{i}_{1}),(v^{l}_{2},I^{i}_{2}),...,(v^{l}_{|{\mathcal{G}}^{L}|},I^{i}_{|{\mathcal{G}}^{L}|})\\}$,
where the indexing set between $v^{l}_{i}$ and its corresponding cluster
$C^{p}_{j}$ is denoted as
$I_{i}=\\{i^{\prime}|i^{\prime}\in\\{1,2,...,|{\mathcal{G}}^{P}|\\},v^{p}_{i^{\prime}}\in
C^{p}_{j},t^{p}_{i^{\prime}}\in T^{l}_{i}\\}$. If $I^{i}_{j}$ is an empty set
$\varnothing$, we consider that the ligand node $v^{l}_{j}$ is not matched to
any cluster. The mixed gaussian distribution model is used to express the
probability density function of the distance between $v^{l}_{j}$ and
$v^{l}_{k}$ when both $I^{i}_{j}$ and $I^{i}_{k}$ are not empty set
$\varnothing$:

$\displaystyle
P(d|M_{i},v^{l}_{j},v^{l}_{k})=F^{i}_{jk}(d)=\frac{1}{N}\sum_{a\in
I^{i}_{j}}\sum_{b\in
I^{i}_{k}}\frac{{\bm{W}}[t^{p}_{a}]{\bm{W}}[t^{p}_{b}]s^{p}_{a}s^{p}_{b}}{\sqrt{2\pi({r^{p}_{a}}^{2}+{r^{p}_{b}}^{2})}}\exp\left[-\frac{(d-{\displaystyle||{\bm{x}}^{p}_{a}-{\bm{x}}^{p}_{b}||)^{2}}}{2({r^{p}_{a}}^{2}+{r^{p}_{b}}^{2})}\right]$
(1)

where ${\bm{W}}[t]$ is the weight for pharmacophore type $t$ and $N$ is the
normalized constant of the mixed gaussian distribution. The weights for each
pharmacophore type are in Appendix 5.4. Note that $M_{i}$, $v^{l}_{j}$,
$v^{l}_{k}$ is independent of ligand 3D conformation.

Then, the score of ligand conformation $\displaystyle{\bm{C}}$ for the graph
match $M_{i}$ can be represented as the weighted sum of pharmacophore pairwise
distance likelihood:

$\displaystyle\text{EdgeScore}({\bm{C}}|M_{i},v^{l}_{j},v^{l}_{k})=\begin{cases}w^{i}_{j}w^{i}_{k}F^{i}_{jk}(||{\bm{x}}^{l}_{j}({\bm{C}})-{\bm{x}}^{l}_{k}({\bm{C}})||)&\text{if}~{}I^{i}_{j}\neq\phi,I^{i}_{k}\neq\phi\\\
0&\text{else}\end{cases}$ (2)
$\displaystyle\text{Score}({\bm{C}},M_{i})=\sum_{j<k}^{|{\mathcal{G}}^{L}|}\text{EdgeScore}({\bm{C}}|M_{i},v^{l}_{j},v^{l}_{k})$
(3)

where $w^{i}_{j}=\sum_{a\in
I^{i}_{j}}{\displaystyle{\bm{W}}[t^{p}_{a}]s^{p}_{a}}/|I^{i}_{j}|$.

##  3 Experiments

###  3.1 Data

We used the PDBbind v2020 dataset [19], a collection of high-resolution
crystal structures, and measured binding affinities for 19,443 protein-ligand
complexes deposited in the Protein Data Bank (PDB) [23]. Following Shen et al.
[22], the dataset is divided into 17,658 training complexes and 1,500
validation complexes, excluding the 285 CASF-2016 [24] test complexes. We
excluded 8 complexes from the training set which contain multiple ligands
detected by Open Babel [25].

In this work, we used an image size of $64\times 64\times 64$ with a
resolution of $0.5$ Å to represent the binding site, resulting in a box side
length of $32$ Å. Note that the side length of the recommended maximum search
box of Autodock Vina is $30$ Å. Additionally, we applied random rotations and
translations to the binding site for equivariance [26, 27].

To establish the ground truth pharmacophore model of each complex for model
training, we considered the protein-ligand functional group pairs that form
non-covalent interactions as the pairs of the hot spot and the pharmacophore.
The identification of non-covalent interactions was performed using the
Protein-Ligand Interaction Profiler (PLIP) [28].

Typically, the length of non-covalent interactions is less than 6.0 Å.
Therefore, there is no need to consider regions distant from the pocket cavity
when predicting hot spots or pharmacophores. We delineate two cavity areas
based on the average length of interaction types: one within 5.0 Å from the
ligand atoms in crystal structure for short-range interactions (hydrophobic
interaction, H-bond, and halogen bond) and another within 7.0 Å for long-range
interactions ($\pi$–$\pi$ stacking, cation–$\pi$ interaction, salt bridge).

###  3.2 Benchmark Study

##### Virtual screening benchmark.

The widely-adopted virtual screening benchmarks are CASF-2016 [24], DUD-E
[29], and DEKOIS2.0 [30]. However, CASF-2016 only contains 285 molecules and
is not well-suited as a benchmark for the pre-screening task. In contrast,
both the DUD-E and DEKOIS2.0 benchmarks are more similar to real-world
scenarios. The DUD-E benchmark comprises 102 targets and offers several
hundred actives, along with 50 decoys for each active, all sharing similar
physicochemical properties. Likewise, DEKOIS2.0 contains 81 distinct targets,
with 40 actives and a substantial 1200 decoys for each target. Consequently,
our benchmark tests were conducted on both DUD-E and DEKOIS2.0 datasets.

In our evaluation, we compared PharmacoNet with state-of-the-art conventional
docking methods such as GLIDE (commercial) [7], AutoDock Vina [8], and SMINA
[9]. In addition, the comparison with the non-structure-based deep learning
model for pre-screening, DeepBindGCN [31], is in Appendix 5.1.1.

##### Runtime benchmark.

We followed the benchmark proposed by Stärk et al. [10], which measures the
average runtime on 363 complexes discovered since 2019 from PDBBind v2020.
Along with the methods in the conventional docking programs, we also include
structure-based deep learning models for binding pose prediction, GNINA [32],
EquiBind [10], TANKBind [11], and DiffDock [12].

##### Evaluation metric.

The screening power is assessed with the top $\alpha\%$ enrichment factor,
EFα%, which is defined as the ratio of the percentages of the active ligands
in the top $\alpha\%$ predictions to the percentages of the active ligands in
the total library, and AUROC, the average area under the receiver operating
characteristic. For runtime measurements, the processes that are performed
before virtual screening process are not considered, such as preprocessing of
protein or generation of initial conformers for ligand.

##### Inference details.

We generated ligand conformers with RDKit ETKDG [33, 34] from the SMILES and
used the top graph matching score among the conformers for screening. As the
default exhaustiveness of Autodock Vina and SMINA is 8, we compared the
screening performance of our model when using a single conformer and 8
conformers. We ran our model on the Intel Xeon Gold 6234.

##### Result.

Table 1 highlights the impressive speed of our framework in comparison to
existing energy-based docking methods and emerging deep learning-based models
for binding pose prediction. Our method outperforms the fastest conventional
docking programs by a factor of 10,000, and it is still 10 times faster than
the fastest deep learning model. It is worth noting that the computational
costs do not proportionally increase with the number of ligand conformers
used, as detailed in Appendix 5.1.2. These results underscore the significant
efficiency of our framework in virtual screening using the numerous ligand
conformations.

To achieve this significant speed-up, our framework uses a coarse-grained
graph-matching algorithm that predicts approximate binding poses implicitly
without considering atomic-level modeling. In addition, our scoring function
only contains less than 10 parameters, in contrast to physics-based or deep
learning-based scoring functions. Despite these simplifications, our framework
still shows reasonable accuracy through the pharmacophore-level abstraction,
as shown in Table 2.

Table 1:  Average runtime per prediction of conventional physics-based docking programs and several deep learning-based approaches. The top three lines contain energy-based binding pose optimization approaches and the meddle two lines contain data-driven approaches based on complex crystal structures. The last two lines show our models and the number after our model means the number of RDKit conformers used. "*" indicates that the method runs on a single CPU core.  |  | Average Runtime (s)
---|---|---
Methods | Type | 16-CPU | GPU
GLIDE(c.) | Energy-based binding pose prediction & scoring | 1405* | -
Autodock Vina |  | 205 | -
QVina-W |  | 49 | -
SMINA |  | 146 | -
GNINA |  | 247 | 146
EquiBind | Data-driven binding pose prediction | 0.16 | 0.04
TANKBind |  | 0.54 | 0.28
DiffDock |  | - | 10/40
PharmacoNet (8) | Coarse-grained graph matching & scoring | 0.0022 | -
PharmacoNet (1) |  | 0.0017 | -
Table 2:  Results of the virtual screening benchmarks. The top half contains conventional molecular docking methods, and the bottom half contains our framework.  | DEKOIS2.0 | DUD-E
---|---|---
Methods | EF0.5% | EF1% | EF5% | AUROC | EF1% | AUROC
GLIDE(c.) | 14.6 | 12.5 | 6.30 | 0.747 | 23.6 | 0.820
Autodock Vina | 5.46 | 4.51 | 2.82 | 0.633 | 9.93 | 0.725
SMINA | 4.62 | 4.04 | 2.72 | 0.627 | 9.57 | 0.722
PharmacoNet (8) | 3.75 | 3.42 | 2.41 | 0.624 | 5.73 | 0.645
PharmacoNet (1) | 3.43 | 3.23 | 2.55 | 0.625 | 5.36 | 0.652

###  3.3 Case Study: Pre-screening for Large-Scale Virtual Screening

To access the effectiveness of our model’s pre-screening capabilities, we
performed a large-scale virtual screening of 1 million ChEMBL molecules for
KRAS (G12C), a challenging drug target. One of the key objectives in the pre-
screening is to retain the hit candidates, the molecules with the best docking
scores. Note that only the top 100 to 1000 molecules in the entire library are
considered as hit candidates for experimental validation. In evaluating the
pre-screening performance of our framework, we used the top-$N$ enrichment
factor, denoted EFN, which is defined as the ratio of the percentages of hit
candidates in the retained molecules to the percentages of these candidates in
the entire library, when the hit candidates are considered as the top $N$
molecules with the highest docking score in the entire library. Furthermore,
we measured the average docking score in conjunction with the enrichment
factor. For the pre-screening stage, we used the top graph matching score
among the 10 RDKit ETKDG conformers for each ligand, and for the fine-
screening stage, we used SMINA with default settings.

Table 3:  Results of the pre-screening. ‘Filtration rate’ and ‘Number of retained’ indicates the filtration rates of pre-screening and the number of remaining molecules after pre-screening, respectively. ‘Total time’ is the time spent on the entire virtual screening including pre-screening (PharmacoNet) and fine-screening (SMINA). The last line represents the baseline scenario without any pre-screening.  Filtration | Number of | KRAS | Average | Total time (hours)
---|---|---|---|---
rate | retained | EF100 | EF500 | EF1000 | docking score | 16-CPU
90% | 100,000 | 5.46 | 4.86 | 4.76 | -8.35 | 150
95% | 50,000 | 8.53 | 7.74 | 7.50 | -8.43 | 75
98% | 20,000 | 9.92 | 9.33 | 9.58 | -8.50 | 31
99% | 10,000 | 15.9 | 12.5 | 12.0 | -8.55 | 15
- | 992,235 | - | - | - | -7.66 | 1300

(a) Pharmacophore model

(b) AMG-510, crystal structure

(c) CHEMBL183526

Figure 3:  Visualizations of the pharmacophore model for KRAS. The color of
pharmacophore and hot spot is as follows: orange for hydrophobic carbon,
purple for aromatic ring, red for anion, blue for cation, cyan for H-bond
donor, magenta for H-bond acceptor, and yellow for halogen atom and halogen
bond acceptor. (a) The generated pharmacophore model of the binding site (PDB
ID: 6OIM). (b) The crystal structure of the known inhibitor (AMG-510). (c) The
SMINA docking pose of the ligand with the highest score of pre-screening.

The pre-screening process was remarkably efficient, taking only 11 minutes for
1 million molecules on 16 CPU cores. The results in Table 3 show that our
model can significantly reduce the total virtual screening time while
effectively retaining hit candidate molecules. By filtering out 90% of the
molecules from the library, our model retained more than half of the top 100
hit candidates (0.01% of the library), and 16 of those with a 99% filtration
rate. Moreover, it’s worth noting that increasing the filtration rate of pre-
screening improved the average docking scores of the remaining molecules.
These observations suggest that there is a trend between our scores and
energy-based docking scores, even though energy is not considered in model
training.

Figure 3 shows the generated pharmacophore model for KRAS. Interestingly, in
the SMINA docking pose of CHEMBL183526, which achieved the highest matching
score, the chlorine atom (indicated by a red circle) aligns with the halogen
atom pharmacophores (represented by yellow spheres) in our model,
demonstrating the ability of our model to identify pharmacophore sites.
Notably, the corresponding halogen bond was absent from both the X-ray (PDBID:
6OIM) and cryo-EM (PDBID: 8G47) structures of AMG-510, a known inhibitor of
KRAS. This indicates the potential of our model to reveal unknown
pharmacophores contributing to energy stabilization, which are not present in
the complex structure of the active ligand.

###  3.4 Ablation Study

Table 4:  Ablation results. All methods are performed with 8 RDKit ETKDG conformers. "*" means the threshold for hot spot prediction is based on the training set instead of the validation set.  | DEKOIS2.0 | DUD-E
---|---|---
Methods | EF0.5% | EF1% | EF5% | AUROC | EF0.5% | EF1% | EF5% | AUROC
No weight | 1.57 | 1.73 | 2.04 | 0.603 | 4.85 | 4.19 | 2.72 | 0.624
Limited trainset | 3.53 | 3.07 | 2.54 | 0.622 | 5.61 | 4.78 | 3.28 | 0.643
Baseline | 3.75 | 3.42 | 2.41 | 0.624 | 7.23 | 5.73 | 3.23 | 0.645

To assess the impact of the size of the training dataset, we conducted
experiments using a restricted dataset of 5000 training complexes and 500
validation complexes randomly sampled from the original training and
validation sets, respectively. Table 4 shows that the model trained on this
limited dataset shows a slightly lower overall performance. However, it
remains competitive in terms of AUROC and outperforms the baseline in EF5%.
This result indicates that our model can learn overall pharmacological
features even with a limited training set.

Furthermore, we observed that the weights for pharmacophore types, which are
the only parameters of the scoring function, have a significant impact on
accuracy. Without using weights to pharmacophore types, the scoring function
does not reflect physics and chemistry, resulting in a degradation of
accuracy. In our study, we used straightforward parameters based on prior
knowledge such as the fact that salt bridges are stronger than H-bond, as
detailed in Appendix 5.4. Given that the current scoring function only used 7
parameters based on 6 interaction types, higher performance can be reached by
adding features such as atom type and parameter fitting.

###  3.5 Limitations

PharmacoNet relies on a scoring and graph-matching algorithm at the
pharmacophore level, which has reasonably accurate and strong generalization
capabilities without the need for a detailed parameter-fitting procedure.
However, this approach has a trade-off; it does not capture atomic-level
features. For example, it cannot distinguish between variations in the
strength of the same salt bridge due to differences in charge and atom type,
nor does it account for intramolecular energy.

To address these limitations, it may be useful to integrate force-field-based
approaches such as UFF [35] and MMFF [36] and atomic-level details into the
graph-matching and scoring algorithm, where a parameter-fitting process
becomes essential. In particular, a promising avenue for improvement is the
integration of a machine learning or deep learning-based approach into the
graph-matching algorithm. This integration will significantly enhance our
framework’s ability to account for fine-grained details and account for
intramolecular energy factors, ultimately leading to improved accuracy and
performance.

Furthermore, designing a balanced scoring model for both binding pose
prediction and screening has long been a challenge in the molecular docking
task [22, 37]. The success of the re-scoring models [22, 37, 38, 39] based on
docking pose indicates that our framework has the potential to improve
accuracy by integrating additional scoring functions tailored to each task,
rather than relying on a single function for both.

##  4 Conclusion

In this work, we have demonstrated for the first time the ability and
potential of a deep learning approach based on pharmacophore modeling. Our
novel model, PharmacoNet, a 3D image instance segmentation model for
generating pharmacophore models, opens up fresh avenues of the structure-based
deep learning approach by coarse-grained modeling of proteins and ligands at
the pharmacophore level. With a coarse-grained graph-matching algorithm and a
simple distance likelihood-based scoring function, our framework achieves
significant speed-up with reasonable accuracy and high generalization ability.
Moreover, our framework has revealed its potential as a pre-screening method
in the realm of large-scale virtual screening. We believe that the
incorporation of pharmacophore modeling and pharmacophoric knowledge into deep
learning-based methods paves the way for innovative breakthroughs in drug
discovery.

## Acknowledgments and Disclosure of Funding

This work was supported by Basic Science Research Programs through the
National Research Foundation of Korea (NRF), grant funded by the Ministry of
Science and ICT (NRF-2018R1A5A1025208, NRF-2023R1A2C2004376,
RS-2023-00257479).

##  5 Supplementary Material

###  5.1 Additional Result

####  5.1.1 Comparison with Non-structure-based Scoring Model for Pre-
screening

Zhang et al. [31] developed a deep learning model DeepBindGCN for pre-
screening that works for multiple targets by integrating the 1D molecular
vector of a pocket and a ligand. The model uses only the information of
protein pocket and ligand 2D graph without the information of complex, so it
can perform scoring by combining the pre-calculated vector of ligands and
proteins.

However, spatial and topological information is lost during global pooling to
represent them as a 1D vector. Since the spatial coordination of protein and
ligand is crucial in determining their interactions, the loss of spatial
information leads to low performance and generalization ability for unseen
target [40, 38]. Table 5 shows that DeepBindGCN is faster than our model in
the GPU environment, but our model is faster in the CPU environment and shows
better accuracy.

Table 5:  Comparison with non-structure-based deep learning model for pre-screening, DeepBindGCN-RG. EF0.5%, EF1%, and EF5% are average top 0.5%, 1.0%, 5.0% enrichment factor, respectively. PDBBind is the runtime benchmark where ’16-CPU’ and ’GPU’ mean the average runtime for each device. The number in parenthesis after our model means the number of RDKit conformers used.  | DEKOIS2.0 | DUD-E | PDBBind
---|---|---|---
Methods | EF0.5% | EF1% | EF5% | EF1% | 16-CPU (s) | GPU (s)
DeepBindGCN-RG | 0.89 | 0.80 | 0.82 | 2.11 | 0.003 | 0.0005
PharmacoNet (8) | 3.75 | 3.42 | 2.41 | 5.73 | 0.0022 | -
PharmacoNet (1) | 3.43 | 3.23 | 2.55 | 5.36 | 0.0017 | -

####  5.1.2 Runtime according to the number of Ligand Conformers

Figure 4:  Average runtime according to the number of RDKit ETKDG
conformations.

We have measured execution times using the PDBbind test set while varying the
number of conformers. Notably, our coarse-grained graph-matching algorithm can
process multiple conformations simultaneously, resulting in improved
efficiency as the number of conformers increases. Specifically, the execution
time for a single conformation was found to be 1.7 ms, while for 20 conformers
it was 2.8 ms. All measurements were made in a 16-CPU environment with the
same machine, Intel Xeon Gold 6234.

###  5.2 Detailed Model Architecture

##### Atomic feature.

For each atom in the binding site, we use the residue type (20 common amino
acids and unknown for others), atom type (C, N, O, S, and unknown for others),
and functional group type (hydrophobic carbon, H-bond donor & acceptor,
halogen bond acceptor, aromatic ring, cation, and anion) as an atomic feature.

##### Voxelization of binding site.

We represent the protein binding site as a 3D voxelized image. For the protein
binding site $\displaystyle{\mathcal{B}}$, the atomic feature, and coordinates
of the binding site atom are denoted ${\bm{h}}^{b}_{i}$ and
${\bm{x}}^{b}_{i}$, respectively. Then, the 3D input image
$\displaystyle{\bm{\mathsfit{I}}}\in\mathbb{R}^{D\times H\times W\times C}$ is
represented as follows:

$\displaystyle{\bm{\mathsfit{I}}}_{d,h,w,:}=\sum_{i}^{N^{b}}{f(||T^{{\bm{\mathsfit{I}}}\rightarrow{\mathcal{B}}}(d,h,w)-{\bm{x}}^{b}_{i}||)\times{\bm{h}}^{b}_{i}}$
(4)

where $N^{b}$ is the number of atoms in binding site, $D$, $H$, $W$ are the
spatial dimensions, and $C$ is the number of atomic features. To better
represent the position of the atoms, we use a transform function
$T^{A\rightarrow B}$ that converts coordinates from space $A$ to space $B$ and
a kernel function $\displaystyle f(x)=\bm{1}_{x\leq 1.5}e^{-2x^{2}}$ where
$\bm{1}_{x\leq 1.5}$ is $1$ when $x\leq 1.5$ and $0$ otherwise.

##### Voxelization of ground truth pharmacophore model.

The ground truth (see Section 3.1) is voxelized to 3D binary masks
corresponding to each hot spot. Voxels within 1.0 A of each pharmacophore are
marked as 1, otherwise 0. The shape of the ground truth 3D binary masks are
$(D,H,W)$.

##### Hot spot prediction.

To detect hot spots, PharmacoNet leverages the Feature Pyramid Network [41]
with a 3D extension of the Swin Transformer V2 Encoder [42], which is
typically used in object detection tasks, to obtain multi-scale 3D feature
maps denoted as ${\bm{\mathsfit{F}}}^{s}\in\mathbb{R}^{sD\times sH\times
sW\times F_{s}}$, where $s$ represents the scale, and $F_{s}$ is the
corresponding hidden dimension of each feature map.

$\displaystyle{\mathbb{F}}=\\{{\bm{\mathsfit{F}}}^{1},{\bm{\mathsfit{F}}}^{1/2},...\\}=\phi^{\text{backbone}}({\bm{\mathsfit{I}}})$

Each protein functional group (token), denoted as $t_{i}^{t}$, can be
represented as a combination of two components: the hidden feature obtained
from the closest voxel to its coordinates ${\bm{x}}_{i}^{t}$, and an embedding
vector ${\mathcal{E}}_{c_{i}^{t}}$ corresponding to its non-covalent
interaction type (class), $c_{i}^{t}$. Since some functional group has the
potential to form multiple types of interactions, PharmacoNet distinguish
between them with their embedding vectors.

$\displaystyle
d_{i},h_{i},w_{i}=T^{{\mathcal{B}}\rightarrow{\bm{\mathsfit{F}}}^{1}}({\bm{x}}^{t}_{i})$
$\displaystyle{\bm{f}}^{t}_{i}=\phi^{\text{token}}({\bm{\mathsfit{F}}}^{1}_{d_{i},h_{i},w_{i},:},{\mathcal{E}}_{c^{t}_{i}})$
$\displaystyle s(t^{t}_{i})=\sigma({\bm{f}}^{t}_{i})$ (5)

Given that the number of hot spots is a small fraction of all tokens, we use
the score’s relative position in the score distribution of the validation set
as token score instead of directly using the sigmoid score. The top 40% for
$\pi-\pi$ stacking, cation-$\pi$ interactions, and salt bridge, and the top
20% for hydrophobic interaction, H-bond, and halogen bond of the score
distribution were selected as the threshold of hot spot prediction. (see
Appendix 5.5)

##### Cavity prediction.

Typically, the length of the non-covalent interaction is shorter than 6.0 Å,
and there is no need to consider regions far from the pocket cavity when
predicting hot spots or pharmacophores. Our model predicts the two cavity
areas, one for short-range interaction, such as hydrophobic interaction,
H-bond, halogen bond, and another for long-range interaction, such as
$\pi-\pi$ stacking, cation-$\pi$ interactions, and salt bridge. The hot spots
and pharmacophore prediction are not performed outside of the corresponding
cavity.

$\displaystyle\displaystyle{\bm{\mathsfit{C}}}^{\text{long}}$
$\displaystyle=\sigma(\phi^{\text{cavity}}_{\text{long}}({\bm{\mathsfit{F}}}^{1}))$
(6) $\displaystyle{\bm{\mathsfit{C}}}^{\text{short}}$
$\displaystyle=\sigma(\phi^{\text{cavity}}_{\text{short}}({\bm{\mathsfit{F}}}^{1}))$
(7)

where
${\bm{\mathsfit{C}}}^{\text{long}},{\bm{\mathsfit{C}}}^{\text{short}}\in\\{0,1\\}^{D\times
H\times W}$ represent binarized outputs with a threshold of 0.5 for the cavity
areas corresponding to long-range and short-range interactions, respectively.

##### Pharmacophore prediction.

For the selected hot spots, our model predicts the optimal pharmacophore
positions, ensuring the non-covalent interactions in an instance segmentation
manner.

$\displaystyle{\bm{\mathsfit{D}}}^{i}=\sigma(\phi^{\text{mask}}({\mathbb{F}},{\bm{f}}^{t}_{i},{\bm{x}}^{t}_{i}))\circ{\bm{\mathsfit{C}}}^{\text{short}}$
(8)

where the $\displaystyle{\bm{\mathsfit{D}}}^{i}\in[0,1]^{D\times H\times W}$
is the model output, and ${\bm{\mathsfit{M}}}^{i}\in\\{0,1\\}^{D\times H\times
W}$ is the binary mask of ${\bm{\mathsfit{D}}}^{i}$ (threshold of 0.5). When
there are separated segments in the binary mask, we consider each segment as a
distinct pharmacophore since one protein functional group can form the same
type of interaction with multiple pharmacophores.

Each pharmacophore contains the voxels, and the center and radius are defined
as follows:

center
$\displaystyle=\frac{1}{\sum_{i}^{N}{{\bm{\mathsfit{D}}}_{d_{i},h_{i},w_{i}}}}\sum_{i}^{N}{{\bm{\mathsfit{D}}}_{d_{i},h_{i},w_{i}}T^{{\bm{\mathsfit{M}}}\rightarrow{\mathcal{B}}}\left([d_{i},h_{i},w_{i}]\right)}$
(9) radius $\displaystyle=\sqrt[3]{N/(4\pi/3)}\times\text{resolution}$ (10)

where $N$ is the number of voxels, $d_{i},h_{i},w_{i}$ is the coordinates of
each voxel, and ${\bm{\mathsfit{D}}}$ and ${\bm{\mathsfit{M}}}$ are the model
output and binary mask corresponding to the pharmacophore.

###  5.3 Training & Inference Detail

##### Training details.

Unlike the Encoder based on Swin Transformer V2, predicting the Binary Mask
for each Pharmacophore is done with 3D-CNN. This requires a very large amount
of computation. Therefore, we train by sampling up to 6 pharmacophores in each
iteration of training.

Also, the number of hotspots (instances) is very small compared to the total
number of functional groups (tokens) contained in the binding site.
Furthermore, distant functional groups have less predictive importance.
Therefore, only a maximum of 200 tokens are utilized for training, prioritized
in the order of hotspots, tokens close to the ligand, and tokens distant from
the ligand.

In contrast to the Encoder based on Swin Transformer V2, the prediction of
Binary Masks for each Pharmacophore is carried out using a 3D-CNN which is
computationally expensive. As a result, we train the model by sampling up to 6
pharmacophores for each complex in each training iteration.

Furthermore, the number of hotspots (instances) is considerably smaller
compared to the total number of functional groups (tokens) present within the
binding site. Consequently, we sample up to 200 tokens for each complex, with
priority given to hotspots, tokens close to the ligand, and tokens far from
the ligand, to reduce the data imbalance.

We use AdamW optimizer [43] for 200,000 iterations with a batch size of 8, the
learning rate of $1\times 10^{-4}$, weight decay of 0.05, and
$(\beta_{1},\beta_{2})$ of $(0.9,0.999)$. The model training takes about 48
hours on a 4 NVIDIA RTX 3090 GPUs. The hidden size of the model is based on
the SwinV2-T, and the total number of parameters is 31M parameters.

##### Inference details.

The voxel outside the maximum distance of the non-covalent interaction (based
on the PLIP rule) from the corresponding hot spot is masked with 0. In
addition, the voxels within 1.0 Å of all protein atoms are masked. Finally,
the model performed gaussian smoothing with a kernel size of 5 and a sigma of
0.5 to reduce the noise in the segmentation map.

###  5.4 Graph-Matching Algorithm.

All notations are defined in Section 2.2. In the context of inexact graph
matching, the number of all possible graph matches is
$\displaystyle|{\mathbb{C}}^{P}|^{|{\mathcal{G}}^{L}|}$. This means that the
time complexity of the graph matching is
$\displaystyle{\mathcal{O}}(|{\mathbb{C}}^{P}|^{|{\mathcal{G}}^{L}|})$, and
the computational requirements increase exponentially with the complexity of
the ligand or pharmacophore model. Therefore, the demand for an efficient
graph-matching algorithm becomes mandatory.

##### Parameters for scoring function.

The only parameters of the scoring function are the weights for pharmacophore
types. Hydrogen bonds and halogen bonds exhibit comparable characteristics and
strengths, both significantly outweighing hydrophobic interactions.
Additionally, salt bridges are stronger than hydrogen bonds or halogen bonds.
$\pi-\pi$ stackings are the major driving force that stabilizes the complex in
terms of entropy. Based on this prior knowledge, we set the weight as follows:

$\displaystyle{\bm{W}}[t]=\begin{cases}1&\text{if }t\text{ is Hydrophobic
carbon}\\\ 4&\text{if }t\text{ in }\\{\text{H-bond donor, H-bond acceptor,
Halogen atom}\\}\\\ 8&\text{if }t\text{ in }\\{\text{Cation, Anion, Aromatic
ring}\\}\end{cases}$

##### Ligand pharmacophore arrangement clustering.

In some cases, a ligand’s functional group may comprise numerous
pharmacophores. For example, benzene contains 7 pharmacophores: an aromatic
ring and six hydrophobic carbons. To improve the efficiency of the graph
matching, we perform clustering according to the following criteria:

  * •

The cations, anions, H-bond acceptors, and H-bond donors in the same
functional group are grouped.

  * •

The aromatic ring and hydrophobic carbons in the same functional group are
grouped.

  * •

The hydrophobic carbons in the carbon chain are grouped.

##### Distance-constraint graph matching.

Indeed, most graph match is physically infeasible, for example, it is
impossible for two ligand pharmacophores that are only 1Å apart to correspond
to two pharmacophore model clusters that are 10Å apart. To account for this,
we use the following distance constraint between $(v^{l}_{j},I^{i}_{j})$ and
$(v^{l}_{k},I^{i}_{k})$:

$\displaystyle
f(d_{1},d_{2},r_{i},r_{j})=\begin{cases}1&\text{if}~{}|d_{1}-d_{2}|\leq
2\sqrt{r_{i}^{2}+r_{j}^{2}}\\\ 0&\text{else}\\\ \end{cases}$
$\displaystyle\sum_{a\in I^{i}_{j}}\sum_{b\in
I^{i}_{k}}f(||{\bm{x}}^{l}_{j}-{\bm{x}}^{l}_{k}||,||{\bm{x}}^{p}_{a}-{\bm{x}}^{p}_{b}||,r^{p}_{a},r^{p}_{b})\geq
0.5|I^{i}_{j}||I^{i}_{k}|$ (11)

##### Pharmacophore model clustering.

Algorithm 1 describes the clustering algorithm for pharmacophore model.

Algorithm 1 Pharmacophore Model Graph Clustering Algorithm

1:UsedNodes $\leftarrow\\{\\}$

2:$\displaystyle{\mathbb{C}}^{P}\leftarrow\\{\\}$

3:for $v^{p}_{i}$ in $\displaystyle{\mathcal{G}}^{P}$ do

4: if $v^{p}_{i}$ in UsedNodes then

5: continue

6: end if

7: $C^{p}\leftarrow\\{\\}$

8: if $t^{p}_{i}$ in {Cation, Anion, Aromatic} then

9: for $v^{p}_{j}\neq v^{p}_{i}$ in $\displaystyle{\mathcal{G}}^{P}$ do

10: if $t^{p}_{i}=t^{p}_{j}$ and
$\displaystyle||{\bm{x}}^{p}_{i}-{\bm{x}}^{p}_{j}||<1.5$Å then

11: add $\\{v^{p}_{j}\\}$ to $C^{p}$

12: else if $\displaystyle||{\bm{x}}^{p}_{i}-{\bm{x}}^{p}_{j}||<3.0$Å then

13: if $t^{p}_{i}$ in {Cation, Anion} and $t^{p}_{j}$ in {H-bondDonor,
H-bondAcceptor} then

14: add $v^{p}_{j}$ to $C^{p}$

15: else if $t^{p}_{i}$ is Aromatic and $t^{p}_{j}$ is HydrophobicCarbon then

16: add $v^{p}_{j}$ to $C^{p}$

17: end if

18: end if

19: end for

20: end if

21: add $C^{p}$ to $\displaystyle{\mathbb{C}}^{P}$

22: UsedNodes $\leftarrow$ UsedNodes $\cup~{}C^{p}$

23:end for

24:for $v^{p}_{i}$ in $\displaystyle{\mathcal{G}}^{P}$ do

25: if $v^{p}_{i}$ in UsedNodes then

26: continue

27: end if

28: $C^{p}\leftarrow\\{\\}$

29: if $t^{p}_{i}$ in {HydrophobicCarbon, HalogenAtom} then

30: for $v^{p}_{j}\neq v^{p}_{j}$ in $\displaystyle{\mathcal{G}}^{P}$ do

31: if $t^{p}_{i}=t^{p}_{j}$ and
$\displaystyle||{\bm{x}}^{p}_{i}-{\bm{x}}^{p}_{j}||<3.0$Å then

32: add $\\{v^{p}_{j}\\}$ to $C^{p}$

33: end if

34: end for

35: else if $t^{p}_{i}$ in {H-bondDonor, H-bondAcceptor} then

36: for $v^{p}_{j}\neq v^{p}_{j}$ in $\displaystyle{\mathcal{G}}^{P}$ do

37: if $t^{p}_{j}$ in {H-bondDonor, H-bondAcceptor} and
$\displaystyle||{\bm{x}}^{p}_{i}-{\bm{x}}^{p}_{j}||<3.0$Å then

38: add $\\{v^{p}_{j}\\}$ to $C^{p}$

39: end if

40: end for

41: end if

42: add $C^{p}$ to $\displaystyle{\mathbb{C}}^{P}$

43: UsedNodes $\leftarrow$ UsedNodes $\cup~{}C^{p}$

44:end for

###  5.5 Hot Spot Detection

Figure 5:  The score distribution of hot spots in the validation set. The red
line is the threshold for hot spot prediction.

###  5.6 Related Work

##### Structure-based pharmacophore modeling.

Besides molecular docking, structure-based pharmacophore modeling is the major
strategy of virtual screening. The pharmacophore modeling approaches are
divided into complex-based [44, 45, 46] and protein-based (receptor-based)
[47, 48, 49, 50, 51, 52]. Since the active ligand for the target binding site
is often unknown in real-world applications, our work focuses on protein-based
pharmacophore modeling.

Protein-based pharmacophore search has been shown in many case studies to
perform competitively with molecular docking on virtual screening task [53,
54]. However, protein-based approaches often identify too many pharmacophoric
features compared to complex-based approaches using active ligand information,
making it difficult to determine the optimal supramolecular interaction
pattern [55]. Therefore, many pharmacophore modeling programs still rely on
expert intuition to select the appropriate number of pharmacophores, which may
require experimentally/computationally expensive procedures such as molecular
dynamics, molecular docking, or fragment crystal structures [16].

##### Protein-ligand binding pose prediction.

Traditional binding pose prediction method is based on physical energy-based
scoring function and optimization algorithm [7, 8, 9]. Recent works improve
the performance by incorporating deep learning-based scoring function [32,
56]. Another promising approach is to frame the pose prediction task as a
regression problem [10, 11] or a generative modeling problem [12] with deep
learning. EquiBind [10] and TANKBind [11] directly predict the atomic
positions of the binding pose, so they are much faster than the traditional
docking methods. However, their performance has not reached that of
traditional docking, and they often cause steric crashes, leading to incorrect
scores with physics-based and deep learning-based scoring functions.[12, 57]
DiffDock [12] succeeds in generating high-quality conformations through
diffusion generative modeling, but it is not fast enough to be used for pre-
screening due to multiple denoising steps.

##### Rigid exhaustive docking.

The flexibility of the molecules is one of the biggest challenges to
predicting conformations. It makes optimization difficult by not only
increasing the pose space but also creating numerous local optima. To avoid
this problem, one of the main docking strategies is rigid-body docking with an
ensemble of ligand conformers [7, 58, 59, 60, 61, 62, 63]. Rigid exhaustive
docking is a practical docking strategy with sets of numerous core
conformations of ligands [10].

##### Machine learning-based scoring function for virtual screening.

Recently, machine learning-based and deep learning-based scoring functions
have been developed to boost virtual screening. The scoring strategies are
usually divided into non-structure-based [31, 64] and structure-based [38, 65,
66, 67, 37, 39, 22] whether they take into account the conformations of
protein-ligand complexes. Non-structure-based scoring methods have the
advantage of not requiring the use of complex structures However, since the
spatial relationship between protein and ligand is essential in determining
their interaction, structure-based scoring models show better performance and
generalization ability by recognizing key interactions from the complex
conformation [68]. Therefore, non-structure-based models are usually used in
ligand-based virtual screening with known active ligands of the target protein
[69, 70], and structure-based models are used in virtual screening strategies
for general targets.

##### Model generalization to unseen protein.

The generalization ability to unseen proteins is one of the biggest challenges
for structure-based scoring models. Several studies have pointed out that deep
learning models are often trained on the biases inherent in the data instead
of the physical interaction between proteins and ligands [22, 71]. To improve
the generalization ability, incorporating pharmacological properties as an
inductive bias is a promising means in the various tasks of deep learning-
based drug discovery [38, 37, 72, 73].

##### Instance image segmentation.

Image segmentation is the task of partitioning an image into multiple
segments. These segments may be delineated based on either semantic attributes
(e.g. person, dog, building) [74] or instance-level criteria [75]. Semantic
segmentation was tackled as a pixel classification problem, where each pixel
is classified into one of the categories. Therefore, different objects can be
grouped into one segment in semantic segmentation. In contrast, instance
segmentation recognizes the segment and category for each object together.
Instance segmentation is divided into a two-stage approach that detects
objects and then determines the mask of each object [75], and a one-stage
approach that performs segmentation on the entire image and then assigns a
category to each segment [76, 77].

## References

  * Lyu et al. [2019] Jiankun Lyu, Sheng Wang, Trent E Balius, Isha Singh, Anat Levit, Yurii S Moroz, Matthew J O’Meara, Tao Che, Enkhjargal Algaa, Kateryna Tolmachova, et al.  Ultra-large library docking for discovering new chemotypes.  _Nature_ , 566(7743):224–229, 2019.
  * Stein et al. [2020] Reed M Stein, Hye Jin Kang, John D McCorvy, Grant C Glatfelter, Anthony J Jones, Tao Che, Samuel Slocum, Xi-Ping Huang, Olena Savych, Yurii S Moroz, et al.  Virtual discovery of melatonin receptor ligands to modulate circadian rhythms.  _Nature_ , 579(7800):609–614, 2020.
  * Gorgulla et al. [2020] Christoph Gorgulla, Andras Boeszoermenyi, Zi-Fu Wang, Patrick D Fischer, Paul W Coote, Krishna M Padmanabha Das, Yehor S Malets, Dmytro S Radchenko, Yurii S Moroz, David A Scott, et al.  An open-source drug discovery platform enables ultra-large virtual screens.  _Nature_ , 580(7805):663–668, 2020.
  * Luttens et al. [2022] Andreas Luttens, Hjalmar Gullberg, Eldar Abdurakhmanov, Duy Duc Vo, Dario Akaberi, Vladimir O Talibov, Natalia Nekhotiaeva, Laura Vangeel, Steven De Jonghe, Dirk Jochmans, et al.  Ultralarge virtual screening identifies sars-cov-2 main protease inhibitors with broad-spectrum activity against coronaviruses.  _Journal of the American Chemical Society_ , 144(7):2905–2920, 2022.
  * Gentile et al. [2022] Francesco Gentile, Jean Charle Yaacoub, James Gleave, Michael Fernandez, Anh-Tien Ton, Fuqiang Ban, Abraham Stern, and Artem Cherkasov.  Artificial intelligence–enabled virtual screening of ultra-large chemical libraries with deep docking.  _Nature Protocols_ , 17(3):672–697, 2022.
  * Sadybekov et al. [2022] Arman A Sadybekov, Anastasiia V Sadybekov, Yongfeng Liu, Christos Iliopoulos-Tsoutsouvas, Xi-Ping Huang, Julie Pickett, Blake Houser, Nilkanth Patel, Ngan K Tran, Fei Tong, et al.  Synthon-based ligand discovery in virtual libraries of over 11 billion compounds.  _Nature_ , 601(7893):452–459, 2022.
  * Friesner et al. [2004] Richard A Friesner, Jay L Banks, Robert B Murphy, Thomas A Halgren, Jasna J Klicic, Daniel T Mainz, Matthew P Repasky, Eric H Knoll, Mee Shelley, Jason K Perry, 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.
  * Trott and Olson [2010] 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.
  * Koes et al. [2013] David Ryan Koes, Matthew P Baumgartner, and Carlos J Camacho.  Lessons learned in empirical scoring with smina from the csar 2011 benchmarking exercise.  _Journal of chemical information and modeling_ , 53(8):1893–1904, 2013.
  * Stärk et al. [2022] Hannes Stärk, Octavian Ganea, Lagnajit Pattanaik, Regina Barzilay, and Tommi Jaakkola.  Equibind: Geometric deep learning for drug binding structure prediction.  In _International conference on machine learning_ , pages 20503–20521. PMLR, 2022.
  * Lu et al. [2022] Wei Lu, Qifeng Wu, Jixian Zhang, Jiahua Rao, Chengtao Li, and Shuangjia Zheng.  Tankbind: Trigonometry-aware neural networks for drug-protein binding structure prediction.  _Advances in neural information processing systems_ , 35:7236–7249, 2022.
  * Corso et al. [2022] Gabriele Corso, Hannes Stärk, Bowen Jing, Regina Barzilay, and Tommi Jaakkola.  Diffdock: Diffusion steps, twists, and turns for molecular docking.  _arXiv preprint arXiv:2210.01776_ , 2022.
  * Graff et al. [2021] David E Graff, Eugene I Shakhnovich, and Connor W Coley.  Accelerating high-throughput virtual screening through molecular pool-based active learning.  _Chemical science_ , 12(22):7866–7881, 2021.
  * Grygorenko et al. [2020] Oleksandr O Grygorenko, Dmytro S Radchenko, Igor Dziuba, Alexander Chuprina, Kateryna E Gubina, and Yurii S Moroz.  Generating multibillion chemical space of readily accessible screening compounds.  _Iscience_ , 23(11), 2020.
  * Wermuth et al. [1998] Camille-Georges Wermuth, CR Ganellin, Per Lindberg, and LA Mitscher.  Glossary of terms used in medicinal chemistry (iupac recommendations 1998).  _Pure and applied Chemistry_ , 70(5):1129–1143, 1998.
  * Yang [2010] Sheng-Yong Yang.  Pharmacophore modeling and applications in drug discovery: challenges and recent advances.  _Drug discovery today_ , 15(11-12):444–450, 2010.
  * Schaller et al. [2020] David Schaller, Dora Šribar, Theresa Noonan, Lihua Deng, Trung Ngoc Nguyen, Szymon Pach, David Machalz, Marcel Bermudez, and Gerhard Wolber.  Next generation 3d pharmacophore modeling.  _Wiley Interdisciplinary Reviews: Computational Molecular Science_ , 10(4):e1468, 2020.
  * Borkar et al. [2023] Maheshkumar Borkar, Arati Prabhu, Abhishek Kanugo, and Rupesh Kumar Gautam.  Pharmacophore modeling.  In _Computational Approaches in Drug Discovery, Development and Systems Pharmacology_ , pages 159–182. Elsevier, 2023.
  * Liu et al. [2015] Zhihai Liu, Yan Li, Li Han, Jie Li, Jie Liu, Zhixiong Zhao, Wei Nie, Yuchen Liu, and Renxiao Wang.  Pdb-wide collection of binding data: current status of the pdbbind database.  _Bioinformatics_ , 31(3):405–412, 2015.
  * Skalic et al. [2019] Miha Skalic, Alejandro Varela-Rial, José Jiménez, Gerard Martínez-Rosell, and Gianni De Fabritiis.  Ligvoxel: inpainting binding pockets using 3d-convolutional neural networks.  _Bioinformatics_ , 35(2):243–250, 2019.
  * Giordano et al. [2022] Deborah Giordano, Carmen Biancaniello, Maria Antonia Argenio, and Angelo Facchiano.  Drug design by pharmacophore and virtual screening approach.  _Pharmaceuticals_ , 15(5):646, 2022.
  * Shen et al. [2023] Chao Shen, Xujun Zhang, Chang-Yu Hsieh, Yafeng Deng, Dong Wang, Lei Xu, Jian Wu, Dan Li, Yu Kang, Tingjun Hou, et al.  A generalized protein–ligand scoring framework with balanced scoring, docking, ranking and screening powers.  _Chemical Science_ , 14(30):8129–8146, 2023.
  * Berman et al. [2000] Helen M Berman, John Westbrook, Zukang Feng, Gary Gilliland, Talapady N Bhat, Helge Weissig, Ilya N Shindyalov, and Philip E Bourne.  The protein data bank.  _Nucleic acids research_ , 28(1):235–242, 2000.
  * Su et al. [2018] Minyi Su, Qifan Yang, Yu Du, Guoqin Feng, Zhihai Liu, Yan Li, and Renxiao Wang.  Comparative assessment of scoring functions: the casf-2016 update.  _Journal of chemical information and modeling_ , 59(2):895–913, 2018.
  * O’Boyle et al. [2011] Noel M O’Boyle, Michael Banck, Craig A James, Chris Morley, Tim Vandermeersch, and Geoffrey R Hutchison.  Open babel: An open chemical toolbox.  _Journal of cheminformatics_ , 3(1):1–14, 2011.
  * Ragoza et al. [2022] Matthew Ragoza, Tomohide Masuda, and David Ryan Koes.  Generating 3d molecules conditional on receptor binding sites with deep generative models.  _Chemical science_ , 13(9):2701–2713, 2022.
  * Long et al. [2022] Siyu Long, Yi Zhou, Xinyu Dai, and Hao Zhou.  Zero-shot 3d drug design by sketching and generating.  _Advances in Neural Information Processing Systems_ , 35:23894–23907, 2022.
  * Salentin et al. [2015] Sebastian Salentin, Sven Schreiber, V Joachim Haupt, Melissa F Adasme, and Michael Schroeder.  Plip: fully automated protein–ligand interaction profiler.  _Nucleic acids research_ , 43(W1):W443–W447, 2015.
  * Mysinger et al. [2012] Michael M Mysinger, Michael Carchia, John J Irwin, and Brian K Shoichet.  Directory of useful decoys, enhanced (dud-e): better ligands and decoys for better benchmarking.  _Journal of medicinal chemistry_ , 55(14):6582–6594, 2012.
  * Bauer et al. [2013] Matthias R Bauer, Tamer M Ibrahim, Simon M Vogel, and Frank M Boeckler.  Evaluation and optimization of virtual screening workflows with dekois 2.0–a public library of challenging docking benchmark sets.  _Journal of chemical information and modeling_ , 53(6):1447–1462, 2013.
  * Zhang et al. [2023] Haiping Zhang, Konda Mani Saravanan, and John ZH Zhang.  Deepbindgcn: Integrating molecular vector representation with graph convolutional neural networks for protein–ligand interaction prediction.  _Molecules_ , 28(12):4691, 2023.
  * McNutt et al. [2021] Andrew T McNutt, Paul Francoeur, Rishal Aggarwal, Tomohide Masuda, Rocco Meli, Matthew Ragoza, Jocelyn Sunseri, and David Ryan Koes.  Gnina 1.0: molecular docking with deep learning.  _Journal of cheminformatics_ , 13(1):1–20, 2021.
  * Landrum et al. [2006] Greg Landrum et al.  Rdkit: Open-source cheminformatics, 2006.
  * Riniker and Landrum [2015] Sereina Riniker and Gregory A Landrum.  Better informed distance geometry: using what we know to improve conformation generation.  _Journal of chemical information and modeling_ , 55(12):2562–2574, 2015.
  * Rappé et al. [1992] Anthony K Rappé, Carla J Casewit, KS Colwell, William A Goddard III, and W Mason Skiff.  Uff, a full periodic table force field for molecular mechanics and molecular dynamics simulations.  _Journal of the American chemical society_ , 114(25):10024–10035, 1992.
  * Halgren [1996] Thomas A Halgren.  Merck molecular force field. i. basis, form, scope, parameterization, and performance of mmff94.  _Journal of computational chemistry_ , 17(5-6):490–519, 1996.
  * Moon et al. [2023] Seokhyun Moon, Sang-Yeon Hwang, Jaechang Lim, and Woo Youn Kim.  A versatile deep learning-based protein-ligand interaction prediction model for accurate binding affinity scoring and virtual screening.  _arXiv preprint arXiv:2307.01066_ , 2023.
  * Moon et al. [2022] Seokhyun Moon, Wonho Zhung, Soojung Yang, Jaechang Lim, and Woo Youn Kim.  Pignet: a physics-informed deep learning model toward generalized drug–target interaction predictions.  _Chemical Science_ , 13(13):3661–3673, 2022.
  * Shen et al. [2022] Chao Shen, Xujun Zhang, Yafeng Deng, Junbo Gao, Dong Wang, Lei Xu, Peichen Pan, Tingjun Hou, and Yu Kang.  Boosting protein–ligand binding pose prediction and virtual screening based on residue–atom distance likelihood potential and graph transformer.  _Journal of Medicinal Chemistry_ , 65(15):10691–10706, 2022.
  * Moovarkumudalvan et al. [2022] Balasubramanian Moovarkumudalvan, Anupriya Madhukumar Geethakumari, Ramya Ramadoss, Kabir H Biswas, and Borbala Mifsud.  Structure-based virtual screening and functional validation of potential hit molecules targeting the sars-cov-2 main protease.  _Biomolecules_ , 12(12):1754, 2022.
  * Lin et al. [2017] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2117–2125, 2017.
  * Liu et al. [2022] Ze Liu, Han Hu, Yutong Lin, Zhuliang Yao, Zhenda Xie, Yixuan Wei, Jia Ning, Yue Cao, Zheng Zhang, Li Dong, et al.  Swin transformer v2: Scaling up capacity and resolution.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 12009–12019, 2022.
  * Loshchilov and Hutter [2017] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Wolber and Langer [2005] Gerhard Wolber and Thierry Langer.  Ligandscout: 3-d pharmacophores derived from protein-bound ligands and their use as virtual screening filters.  _Journal of chemical information and modeling_ , 45(1):160–169, 2005.
  * Chen and Lai [2006] Jing Chen and Luhua Lai.  Pocket v. 2: further developments on receptor-based pharmacophore modeling.  _Journal of chemical information and modeling_ , 46(6):2684–2691, 2006.
  * Ortuso et al. [2006] Francesco Ortuso, Thierry Langer, and Stefano Alcaro.  Gbpm: Grid-based pharmacophore model: concept and application studies to protein–protein recognition.  _Bioinformatics_ , 22(12):1449–1455, 2006.
  * Böhm [1992] Hans-Joachim Böhm.  The computer program ludi: a new method for the de novo design of enzyme inhibitors.  _Journal of computer-aided molecular design_ , 6:61–78, 1992.
  * Kirchhoff et al. [2001] Paul D Kirchhoff, Rob Brown, Scott Kahn, Marvin Waldman, and CM Venkatachalam.  Application of structure-based focusing to the estrogen receptor.  _Journal of Computational Chemistry_ , 22(10):993–1003, 2001.
  * Schuller et al. [2006] Andreas Schuller, Uli Fechner, Steffen Renner, Lutz Franke, Lutz Weber, and Gisbert Schneider.  A pseudo-ligand approach to virtual screening.  _Combinatorial chemistry & high throughput screening_, 9(5):359–364, 2006.
  * Roland et al. [2015] Wibke SU Roland, Marijn PA Sanders, Leo van Buren, Robin J Gouka, Harry Gruppen, Jean-Paul Vincken, and Tina Ritschel.  Snooker structure-based pharmacophore model explains differences in agonist and blocker binding to bitter receptor htas2r39.  _PLoS One_ , 10(3):e0118200, 2015.
  * Johnson and Karanicolas [2016] David K Johnson and John Karanicolas.  Ultra-high-throughput structure-based virtual screening for small-molecule inhibitors of protein–protein interactions.  _Journal of chemical information and modeling_ , 56(2):399–411, 2016.
  * Lee et al. [2020] Ji Young Lee, James M Krieger, Hongchun Li, and Ivet Bahar.  Pharmmaker: Pharmacophore modeling and hit identification based on druggability simulations.  _Protein Science_ , 29(1):76–86, 2020.
  * Barillari et al. [2008] Caterina Barillari, Gilles Marcou, and Didier Rognan.  Hot-spots-guided receptor-based pharmacophores (hs-pharm): a knowledge-based approach to identify ligand-anchoring atoms in protein cavities and prioritize structure-based pharmacophores.  _Journal of chemical information and modeling_ , 48(7):1396–1410, 2008.
  * Tran-Nguyen et al. [2018] Viet-Khoa Tran-Nguyen, Franck Da Silva, Guillaume Bret, and Didier Rognan.  All in one: Cavity detection, druggability estimate, cavity-based pharmacophore perception, and virtual screening.  _Journal of chemical information and modeling_ , 59(1):573–585, 2018.
  * Qing et al. [2014] Xiaoyu Qing, Xiao Yin Lee, Joren De Raeymaeker, Jeremy RH Tame, Kam YJ Zhang, Marc De Maeyer, and Arnout RD Voet.  Pharmacophore modeling: advances, limitations, and current utility in drug discovery.  _Journal of Receptor, Ligand and Channel Research_ , pages 81–92, 2014.
  * Méndez-Lucio et al. [2021] Oscar Méndez-Lucio, Mazen Ahmad, Ehecatl Antonio del Rio-Chanona, and Jörg Kurt Wegner.  A geometric deep learning approach to predict binding conformations of bioactive molecules.  _Nature Machine Intelligence_ , 3(12):1033–1039, 2021.
  * Yu et al. [2023] Yuejiang Yu, Shuqi Lu, Zhifeng Gao, Hang Zheng, and Guolin Ke.  Do deep learning models really outperform traditional approaches in molecular docking?  _arXiv preprint arXiv:2302.07134_ , 2023.
  * Lape et al. [2010] Michael Lape, Christopher Elam, and Stefan Paula.  Comparison of current docking tools for the simulation of inhibitor binding by the transmembrane domain of the sarco/endoplasmic reticulum calcium atpase.  _Biophysical chemistry_ , 150(1-3):88–97, 2010.
  * Ganea et al. [2021] Octavian-Eugen Ganea, Xinyuan Huang, Charlotte Bunne, Yatao Bian, Regina Barzilay, Tommi Jaakkola, and Andreas Krause.  Independent se (3)-equivariant models for end-to-end rigid protein docking.  _arXiv preprint arXiv:2111.07786_ , 2021.
  * Kuntz et al. [1982] Irwin D Kuntz, Jeffrey M Blaney, Stuart J Oatley, Robert Langridge, and Thomas E Ferrin.  A geometric approach to macromolecule-ligand interactions.  _Journal of molecular biology_ , 161(2):269–288, 1982.
  * Mcgann et al. [2003] Mark R Mcgann, Harold R Almond, Anthony Nicholls, J Andrew Grant, and Frank K Brown.  Gaussian docking functions.  _Biopolymers: Original Research on Biomolecules_ , 68(1):76–90, 2003.
  * McGaughey et al. [2007] Georgia B McGaughey, Robert P Sheridan, Christopher I Bayly, J Chris Culberson, Constantine Kreatsoulas, Stacey Lindsley, Vladimir Maiorov, Jean-Francois Truchon, and Wendy D Cornell.  Comparison of topological, shape, and docking methods in virtual screening.  _Journal of chemical information and modeling_ , 47(4):1504–1519, 2007.
  * Vigers and Rizzi [2004] Guy PA Vigers and James P Rizzi.  Multiple active site corrections for docking and virtual screening.  _Journal of medicinal chemistry_ , 47(1):80–89, 2004.
  * Gentile et al. [2020] Francesco Gentile, Vibudh Agrawal, Michael Hsing, Anh-Tien Ton, Fuqiang Ban, Ulf Norinder, Martin E Gleave, and Artem Cherkasov.  Deep docking: a deep learning platform for augmentation of structure based drug discovery.  _ACS central science_ , 6(6):939–949, 2020.
  * Zheng et al. [2022] Liangzhen Zheng, Jintao Meng, Kai Jiang, Haidong Lan, Zechen Wang, Mingzhi Lin, Weifeng Li, Hongwei Guo, Yanjie Wei, and Yuguang Mu.  Improving protein–ligand docking and screening accuracies by incorporating a scoring function correction term.  _Briefings in Bioinformatics_ , 23(3):bbac051, 2022.
  * Torng and Altman [2019] Wen Torng and Russ B Altman.  Graph convolutional neural networks for predicting drug-target interactions.  _Journal of chemical information and modeling_ , 59(10):4131–4149, 2019.
  * Lim et al. [2019] Jaechang Lim, Seongok Ryu, Kyubyong Park, Yo Joong Choe, Jiyeon Ham, and Woo Youn Kim.  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.
  * Panday and Ghosh [2019] Shailesh Kumar Panday and Indira Ghosh.  In silico structure-based prediction of receptor–ligand binding affinity: current progress and challenges.  _Structural Bioinformatics: Applications in Preclinical Drug Discovery Process_ , pages 109–175, 2019.
  * Bahi and Batouche [2018] Meriem Bahi and Mohamed Batouche.  Deep learning for ligand-based virtual screening in drug discovery.  In _2018 3rd international conference on pattern analysis and intelligent systems (PAIS)_ , pages 1–5. IEEE, 2018.
  * Sharma et al. [2022] Tanuj Sharma, Venu Venkatarame Gowda Saralamma, Duk Chul Lee, Mohammad Azhar Imran, Jaehyuk Choi, Mohammad Hassan Baig, and Jae-June Dong.  Combining structure-based pharmacophore modeling and machine learning for the identification of novel btk inhibitors.  _International Journal of Biological Macromolecules_ , 222:239–250, 2022.
  * Chen et al. [2019] Lieyang Chen, Anthony Cruz, Steven Ramsey, Callum J Dickson, Jose S Duca, Viktor Hornak, David R Koes, and Tom Kurtzman.  Hidden bias in the dud-e dataset leads to misleading performance of deep learning in structure-based virtual screening.  _PloS one_ , 14(8):e0220113, 2019.
  * Zhung et al. [2023] Wonho Zhung, Hyeongwoo Kim, and Woo Youn Kim.  A protein-ligand interaction-focused 3d molecular generative framework for generalizable structure-based drug design.  _chemrxiv_ , 2023.
  * Imrie et al. [2021] Fergus Imrie, Thomas E Hadfield, Anthony R Bradley, and Charlotte M Deane.  Deep generative design with 3d pharmacophoric constraints.  _Chemical science_ , 12(43):14577–14589, 2021.
  * Long et al. [2015] Jonathan Long, Evan Shelhamer, and Trevor Darrell.  Fully convolutional networks for semantic segmentation.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 3431–3440, 2015.
  * He et al. [2017] Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross Girshick.  Mask r-cnn.  In _Proceedings of the IEEE international conference on computer vision_ , pages 2961–2969, 2017.
  * Cheng et al. [2022] Bowen Cheng, Ishan Misra, Alexander G Schwing, Alexander Kirillov, and Rohit Girdhar.  Masked-attention mask transformer for universal image segmentation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 1290–1299, 2022.
  * Jain et al. [2023] Jitesh Jain, Jiachen Li, Mang Tik Chiu, Ali Hassani, Nikita Orlov, and Humphrey Shi.  Oneformer: One transformer to rule universal image segmentation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2989–2998, 2023.