# Boosting Convolutional Neural Networks’ Protein Binding Site Prediction
Capacity using SE(3)-invariant transformers, transfer learning and homology-
based augmentation

Daeseok Lee<EMAIL_ADDRESS>,  Jeunghyun Byun<EMAIL_ADDRESS>and
Bonggun Shin<EMAIL_ADDRESS>
###### Abstract.

Figuring out small molecule binding sites in target proteins, in the
resolution of either pocket or residue, is critical in many virtual and real
drug-discovery scenarios. Since it is not always easy to find such binding
sites based on domain knowledge or traditional methods, different deep
learning methods that predict binding sites out of protein structures have
been developed in recent years. Here we present a new such deep learning
algorithm, that significantly outperformed all state-of-the-art baselines in
terms of the both resolutions—pocket and residue. This good performance was
also demonstrated in a case study involving the protein human serum albumin
and its binding sites. Our algorithm included new ideas both in the model
architecture and in the training method. For the model architecture, it
incorporated SE(3)-invariant geometric self-attention layers that operate on
top of residue-level CNN outputs. This residue-level processing of the model
allowed a transfer learning between the two resolutions, which turned out to
significantly improve the binding pocket prediction. Moreover, we developed
novel augmentation method based on protein homology, which prevented our model
from over-fitting. Overall, we believe that our contribution to the literature
is twofold. First, we provided a new computational method for binding site
prediction that is relevant to real-world applications, as shown by the good
performance on different benchmarks and case study. Second, the novel ideas in
our method—the model architecture, transfer learning and the homology
augmentation—would serve as useful components in future works.

##  1\. Introduction

In structure-based drug discovery, the knowledge of ligand binding sites
(hereafter binding sites) on target proteins is crucial. It can aid rational
drug design ([Hal09], [NKC+16], [SSS+19]) and is required for in-silico
methods such as docking ([BGL+21], [NKC+16]). Such knowledge of binding sites
can be attained by analyses of experimental structures of the target protein
in complex with ligands. However, if no such structure is at one’s disposal,
it may be necessary to rely on computational means in order to identify the
binding sites.

In general, this computational task of Binding Site Prediction (BSP) can be
regarded as composition of two sub-tasks: (1) Binding Site Detection (BSD) and
(2) Binding Residue Identification (BRI).

Firstly, BSD aims to identify the binding sites in a coarse-grained manner and
score their druggability. A successful detection of highly druggable binding
sites can aid medical chemists in many ways when designing better drug
compounds. For example, the medical chemists can draw valuable insights in
improving drug compounds’ binding affinity or physical properties by examining
the receptor structure at the potential binding site ([Hal09]). Also,
preparing a suitable binding site is the first step in any virtual structure-
based drug discovery pipeline ([BGL+21]).

Secondly, BRI aims to identify residues in a given binding site that plays key
roles in interactions with ligands. Identification of the key residues has
been pursued in many previous research papers, due to its importance in
rational drug design ([NKC+16], [SSS+19], [BCM+10]). In particular, it has
several applications in virtual structure-based drug discovery. For example,
structural pharmacophore features can be selected based on the identified key
residues ([NKC+16], [BMR08]), and docking results can be prioritized according
to whether the docked molecule has favored interactions with the key residues
([BGL+21], [NKC+16]).

In this paper, we will focus in particular on the structure-based Deep
Learning methods to tackle the BSP problem. This choice reflects two recent
trends. Firstly, Deep Learning has been widely adopted for BSP ([ZCZ20]), and
has shown good performance. Secondly, it has become easier to identify protein
structures as a result of (1) the rapidly accumulating experimental data in
databases such as PDB, and (2) advances in Deep Learning methods e.g.
Alphafold ([JEP+21]).

Recent structure-based Deep Learning methods for BSP are predominantly based
solely on 3D CNN architectures that operate on grid-shaped inputs. This
results in several limitations, in terms of the performance of two sub-tasks
explained previously.

Firstly, the way the CNN-based methods aggregate the local (or short-range)
information to recognize the global patterns may be sub-optimal. They achieve
this through either clustering algorithms applied on top of the CNN outputs,
or through the CNN layers themselves. In the first case, the clustering
algorithm could be replaced by highly parameterized modules such as Neural
Networks. These can potentially outperform the clustering algorithm, since
many parameters can be optimized for the given task. In the second case, the
deep layers of convolution operations may suffer from the problem of “long-
term dependency”. This means, since a convolution layer’s operation is only
local, a deep hierarchy of convolutions must be applied in order to allow a
neuron to have a receptive field large enough to capture the global patterns.
This long-term dependency is known to impede the training ([HBF+01]). This
problem can be resolved if we adopt a neural network architecture whose
operation is not local in nature.

Secondly, these grid-based models do not directly operate on the protein
residues. Therefore, they need an ad-hoc conversion of outputs to obtain
predictions about the binding residues. For example, they posit that an atom
close to a predicted 3D point (e.g. a predicted voxel center) is part of a
binding residue ([JDMR+17], [SDZS20]). This may lead to sub-optimal BRI
performances, since the loss function used to train the model do not compare
the ground truth with the final output, but with the intermediate output
before the conversion. For a better BRI performance, having a parameterized
model that directly outputs the residue-level predictions would be preferable.

To resolve these problems, we devised a model that has geometric attention
layers that operate on top of residue-level CNN outputs. Our model is composed
of two modules dedicated to the sub-tasks BSD and BRI respectively. Both
modules (1) divide local surroundings around each protein residue into grids,
(2) process the grid features in parallel using a CNN model to obtain residue-
level local features, (3) use the geometric self-attention layers to update
the features, and (4) compute the final reductions to produce subtask-specific
outputs.

To further improve our model’s performance, we devised several additional
elements—transfer learning, SE(3)-invariance, and augmentation. All these
elements are intimately connected to our model’s architecture.

Firstly, we applied the transfer learning between the two modules. More
specifically, the BSD module was not trained from scratch but a portion of its
parameters was initialized from the trained parameters of the BRI module. This
was meant to facilitate the learning process of the BSD module, which may
suffer from a relative lack of data.

Secondly, in order to improve the robustness of our model’s predictions, we
adapted the model to be SE(3)-invariant. In other words, our model is
invariant to rotation and translation of the input structure ([FWFW20]).
Besides using an SE(3)-invariant attention mechanism ([JEP+21]) at the outset,
we made the grid featurization process SE(3)-invariant as well. This was
achieved by aligning the axes of the grids with a specific orientation. This
is similar to the grid alignment method used in [MAD21]. While the grids
aligned with that method is not completely deterministic (as one degree of
freedom remains), the grids aligned in our method do and hence achieve the
full SE(3)-invariance.

Lastly, we devised the data augmentation techniques that can be used along
with our model. This was necessary because existing augmentation methods based
on translations or rotations have no effect on training the SE(3)-invariant
model. In particular, we came up with a novel augmentation method based on
protein homology search and sequence alignment. Protein homology has been
utilized for the BSP problem ([SB09], [WKS10], [RZ12], [YRZ13]), but to the
best of our knowledge, we were the first to use it as an augmentation method
for the same problem.

The resulting method achieved significant performance gains over previous
methods when evaluated on various BSP datasets, both in terms of BSD and BRI.
A BSD metric increased by 3.8% on average, and a BRI metric increased by 16.9%
on average.

Through an ablation study, we showed that all the components of our model—the
model architecture, transfer learning, SE(3)-invariance, and augmentation—made
significant contributions, with a few exceptions for the BSD metric. At the
end of the experiment section, we provided potential reasons for such
exceptions.

Finally, we performed a case study on human serum albumin to show the
effectiveness of our model in real world applications. We based our case study
on two previous studies on the binding sites of human serum albumin ([WHR+13],
[Car03]), and examined how well our method could make predictions that are
compatible with the studies. The promising results showed the potential
usefulness of our method in in-silico drug discovery.

To summarize, our contributions in this paper are:

  * •

We developed a new SE(3)-invariant deep learning model for BRI and BSD that
combine CNN with the geometric self-attention layers.

  * •

We developed data augmentation methods, in particular homology augmentation,
that can be used when training our model.

  * •

We found that our new model, trained with the proposed data augmentation
methods, achieved significant performance gains over state-of-the-art deep
learning methods for BRI and BSD.

  * •

By an ablation study, we found that all elements of our method contributed
significantly to BRI performance.

  * •

By an ablation study, we examined which elements contributed unambiguously to
BSD performance. For those that did not, we provided possible explanations.

##  2\. Related Works

In this section, we discuss (1) existing BSP methods, focusing on the
traditional and the deep learning methods, and (2) a similar problem of
predicting ligand-specific binding sites.

###  2.1. Traditional BSP methods

####  2.1.1. Probe-based Methods

These methods use a fixed set of small molecules called ”probes” to determine
the binding sites in a query protein ([LJ05], [HGS09], [NHZ+12]).
Specifically, they place the probes at different positions on the surface of
the protein, and calculate the physical energy at the positions. The low-
energy positions are predicted to be the potential binding sites.

####  2.1.2. Geometry-based Methods

These methods rely on 3D geometric characterization of binding sites to detect
them. ([BNL03], [LGST09], [ZP11])

One example is Fpocket ([LGST09]), that tries to find concave regions of
appropriate sizes on the protein surface. It does so by approximating the
local curvatures by radii of alpha spheres, which are spheres with four heavy
atoms on them but no heavy atom inside. More specifically, it finds all alpha
spheres within a radius range, clusters them, and filters them according to
the number of constituent alpha spheres to produce a binding pocket.

Although Fpocket typically produces an excessive number of binding pockets, it
has relatively good recall (96.4% on scPDB v.2017, according to [DBRK15]).
Therefore, there are ML-based algorithms ([KH18], [AGC+21]) that make use of
it as a means to generate initial candidate binding sites. Our method employs
the same strategy.

####  2.1.3. Template-based Methods

These methods predict binding sites of a query protein based on the templates,
which are other similar proteins whose binding sites are already known
([SB09], [WKS10], [RZ12], [YRZ13]). A portion of the query protein is regarded
as a binding site if it resembles binding sites of templates either
sequentially or structurally.

For example, the authors of [YRZ13] suggested combining two template-based
approaches, one based on substructure comparison (TM-SITE) and the other on
sequence profile alignment (S-SITE). TM-SITE works as follows:

  1. (1)

Putative binding pockets are identified in the query protein by relying on an
external software ([CLT+09]).

  2. (2)

For each putative binding pocket, the template binding sites similar to it are
collected as putative templates. The similarity measure is based on both
structural and sequential comparisons.

  3. (3)

The ligands in the putative templates are projected to the binding pocket.

  4. (4)

A consensus voting by the projected ligands determines whether the residues in
the binding pocket are in the binding site or not.

On the other hand, S-SITE works as follows:

  1. (1)

The query protein sequence is aligned with the template sequences based on
their position-specific scoring matrices (PSSM) profiles and secondary
structure information.

  2. (2)

Templates with the highest alignment quality scores are chosen as the putative
templates.

  3. (3)

A consensus voting by the templates determines whether the residues in the
query sequence are in the binding site.

Our homology-based augmentation algorithm is inspired by these methods. The
overall flow of it resembles that of TM-SITE, and the use of global sequential
alignment is shared by S-SITE. However, our aim in applying the algorithm is
not to make a final prediction of the model, but rather to augment the
training dataset.

###  2.2. deep learning based BSP Methods

####  2.2.1. DeepSite

DeepSite([JDMR+17]) predicts the binding sites by using a 3D CNN model and a
clustering algorithm. The inference steps of DeepSite are as follows: (1) it
generates points spanning the entire protein-occupied 3D space, (2) predicts
the ligandability of the points using the CNN model computed on a 3D grid
centered at the points, and (3) clusters the ligandable points to produce
binding sites.

####  2.2.2. DeepSurf

DeepSurf([MAD21]) has the same overall procedure with DeepSite([JDMR+17]), but
it uses more sophisticated approaches in several aspects. More specifically,
it tries to improve (1) the generation of initial points, (2) the formation of
input grids, and (3) the architecture of the 3D CNN model. The initial points
are sampled on the Solvent Accessible Surface (SAS) of the protein, rather
than the entire span of the protein. Then, the axes of the grids formed at
those points are not arbitrarily oriented, but one axis is set to be the
normal vector of the SAS. Finally, rather than using a plain CNN model, they
used 3D equivalents of ResNet and Bottleneck ResNet ([HZRS16]).

####  2.2.3. Kalasanty

Kalasanty ([SDZS20]) tries to solve BSP by viewing it as a 3D image
segmentation problem. Therefore, it uses a 3D equivalent of the U-net model
([RFB15]), which was originally developed for 2D images. It applies the U-net
model to large grids that cover most of the query proteins, and outputs the
connected components consisting of positively predicted voxels as binding
sites.

####  2.2.4. Deeppocket

Similar to our proposed method, Deeppocket ([AGC+21]) relies on the binding
site candidates generated by Fpocket. It has separate detection and
segmentation models, where the former is a plain 3D CNN model, and the latter
is a U-net model similar to the one used in Kalasanty. The detection model is
used to rank the binding site candidates generated by Fpocket, and the
segmentation model is used to segment the 3D voxels centered at the top-ranked
sites.

###  2.3. Predicting ligand-specific binding sites

Recently, deep learning models that predict protein-ligand complex structures,
given a protein-ligand pair, have been developed ([LYH+19], [MLAdRCW21],
[SGP+22]). In principle, these models can be used to find ligand-specific
binding sites. Therefore, one may argue that the BSP models are strictly less
useful than these models, since their predictions on binding sites do not take
into account the partner ligands. However, we argue that they are still useful
in their own right. Firstly, for many applications, predicting ligand-agnostic
binding sites is enough, and even desirable. For example, a typical docking
experiment requires a binding pocket location as a prerequisite, and docks all
molecules in a virtual library to the pocket. To predict binding pockets using
the models that do consider the partner ligands, preparing appropriate ligands
may add additional complexity to the problem. This is similar to the problem
of preparing appropriate probe molecules in the previously mentioned prob-
based methods. Secondly, the performances of the methods that predict protein-
ligand complex structures are not satisfactory at this point. For example,
EquiBind ([SGP+22]) scored median and mean ligand rmsd of $6.2\AA$ and
$8.2\AA$, which suggests that the model is not accurate enough to correctly
identify binding residues. Therefore, we might want to focus on the easier and
well-studied BSP problem.

##  3\. Problem Definition

Thus far, the BSP task has not been addressed explicitly under a common
definition across different literatures, even though they used different model
compositions. For example, while Deeppocket is comprised of separate
“detection” and “segmentation” models, Kalasanty only uses single segmentation
model, whose output is post-processed by a clustering algorithm.

In order to fairly assess different BSP models, it is necessary to envisage an
unified definition of the BSP task. To be more specific, we will formally
establish standards on the input and output. All baseline models can be
regarded as following the standard, which will be explained in the experiment
section.

Moreover, we will also explain the decomposition of the task into sub-tasks
(including BSD and BRI), which is employed in our method and Deeppocket.

###  3.1. The BSP task

BSP is the task of identifying the ligand binding sites in a given protein. In
the task, we are given as input: a protein structure $P$ and the number of
binding sites $n$. We assume that there are known structures of ligands
$l_{i}$ ($i=1,\cdots,n$) that correspond to the binding sites.

The goal of the task is to predict an unordered set of $n$ binding sites of
$P$ where the ligands $l_{1},\cdots,l_{n}$ bind.

A predicted binding site is of the form $(\hat{c}_{i},\hat{R}_{i})$, where
$\hat{c}_{i}\in\mathbb{R}^{3}$ is the binding site center and
$\hat{R}_{i}\subset\left\\{1,\cdots,size(P)\right\\}$ is the set of indices of
binding residues. For example, an ideal prediction
$\left\\{(\hat{c}_{1},\hat{R}_{1}),\cdots,(\hat{c}_{n},\hat{R}_{n})\right\\}$
is such that

  * •

$\hat{c}_{i}$ is close (e.g. within the radius threshold $4\AA$) to $l_{i}$

  * •

$\hat{R}_{i}$ is the set of indices of residues close (e.g. within the radius
threshold $4\AA$) to $l_{i}$

The methods that we used to evaluate such predictions will be explained in
5.4.

###  3.2. Decomposition into sub-tasks

Our method divides the BSP task into sub-tasks (1) candidate generation, (2)
Binding Site Detection (BSD) and (3) Binding Residue Identification (BRI),
each corresponding to a dedicated module (this is similar to TM-SITE ([YRZ13])
and Deeppocket ([AGC+21])). To be more specific, let $(P,n)$ be an input to
perform BSP on. First, the candidate generation module takes the protein
structure $P$ as an input then generates the candidate binding site centers
$\hat{c}^{\prime}_{1},\hat{c}^{\prime}_{2},...,\hat{c}^{\prime}_{m}\in\mathbb{R}^{3}$,
where typically $m\gg n$. Next, the BSD module takes $(P,\hat{c}_{i})$ ($1\leq
i\leq m$) as the inputs and outputs the predicted druggability of
$\hat{c}^{\prime}_{i}$ in $P$. The druggability scores are then used to rank
the candidate centers, the top $n$ of which form a filtered list
$\hat{c}_{1},\cdots,\hat{c}_{n}$ of candidate centers. Lastly, for each $1\leq
i\leq n$, the BRI module takes as input $(P,\hat{c}_{i})$, and outputs
$\hat{R}_{i}$, that is the set of binding residues within the binding site.
The resulting set
$\left\\{(\hat{c}_{1},\hat{R}_{1}),\cdots,(\hat{c}_{n},\hat{R}_{n})\right\\}$
becomes the final output of the model.

##  4\. Key Components of Our Method

Figure 1.  (a) Our BRI module. (b) Our BSD module. (c) Our BSD module is
trained in two stages, where in the second stage, the parameters of the shared
parts are initialized from the result of training a BRI module in the first
stage. (d) Ground truth generation. A candidate site is a true binding site if
the center is within $4\AA$ from a ligand atom (i.e. $DCA<4\AA$). A residue is
a binding site residue if it is within $4\AA$ from a ligand atom.  Figure 2.
The figure illustrates, given hidden vectors
$\left\\{x_{i}\right\\}_{i=1}^{n}$ and geometric information
$\left\\{(R_{i},t_{i})\right\\}_{i=1}^{n}$, how our attention layer produces
updated hidden vectors $\left\\{x_{i}^{\prime}\right\\}_{i=1}^{n}$. (a)
Obtaining query, key and value vectors required for the computation. (a-1)
Standard ($q_{i}^{h}$) and geometric ($\textbf{q}_{i}^{hp}$) query vectors.
(a-2) Standard ($k_{i}^{h}$) and geometric ($\textbf{k}_{i}^{hp}$) key
vectors. (a-3) Standard ($v_{i}^{h}$) and geometric ($\textbf{v}_{i}^{hp}$)
value vectors. (b) Obtaining attention weights. (b-1) The attention weights
are based on the standard terms ($w_{ij}^{standard}$) and the geometric terms
$w_{ij}^{geometric}$. (b-2) The standard terms are based on the match
(measured by inner product) between the standard query and key vectors. (b-3)
The geometric terms are based on the match (measured by the distance between
their global coordinates) between the geometric query and key vectors. (c-1)
The updated hidden vectors are calculated based on the aggregated standard
value vectors ($o_{i}^{h}$) and the aggregated global coordinates of the value
vectors ($\textbf{o}_{i}^{hp}$). (c-2) Aggregating the standard value vectors.
(c-3) Aggregating the global coordinates of the geometric value vectors.

This section briefly illustrates the key components of our method, which will
be explained in more details in Section 7. These include the details of the
modules (candidate generation, BSD and BRI) as well as other aspects
independent to the model architecture. In particular, the latter includes the
transfer learning and the homology-based augmentation method.

###  4.1. The candidate generation module

To generate the binding site candidate centers, we use an external software
Fpocket ([LGST09]). Given a protein structure, Fpocket finds sets of heavy
atoms $\hat{S}_{1},\cdots,\hat{S}_{m}$, each corresponding to a region
geometrically likely to be a binding pocket. Then, we find the candidate
centers $\hat{c}^{\prime}_{i}$ ($i=1,\cdots,m$) by taking the center of the
mass of the atoms in $\hat{S}_{i}$.

We chose Fpocket as the candidate generation method because it achieves a
sufficiently high recall rate (96.4% on scPDB v.2017, according to [AGC+21]).
This means that, for a given protein and its binding site, it is likely that
at least one of the generated candidates corresponds to the binding site.
Then, provided that the BSD module ranks the candidates properly, the top-$n$
candidates may approximate the true binding site centers with a high accuracy.

###  4.2. The BSD module

The BSD module takes as input the protein structure and a candidate binding
site center $\hat{c}^{\prime}$, and outputs the predicted druggability at
$\hat{c}^{\prime}$.

In doing so, it featurizes the surroundings of $\hat{c}^{\prime}$ into a set
of per-residue 3D grids, and processes the grids through a neural network to
produce the output. Here, each grid in the set corresponds to a residue close
enough to $\hat{c}^{\prime}$ (distance threshold $17\AA$), and encodes the
local environment of the residue.

The neural network of the BSD module is composed of (1) a residue-local
feature extracting unit which runs in parallel for each grid (2) an
aggregation unit which globally aggregates the local features (3) the
reduction unit which maps the aggregated feature to single scalar quantity.
The feature extracting unit is a 3D CNN model, and the aggregation unit is
composed of several geometric self-attention layers. The reduction part is
composed of a point-wise feed-forward layer and a mean-reduction operation.

###  4.3. The BRI module

The BRI module takes in the protein structure and a putative binding site
$\hat{c}$ as inputs and outputs the set of predicted binding residue indices.

The BRI module shares the residue-local feature extraction and global
aggregation units with the BSD module. To be more specific, the BRI module
shares the following units with the BSD module: (1) the CNN feature extractor
and (2) the stack of geometric attention layers up to the penultimate one in
the BSD module. However, the remaining part of the BRI module is only
comprised of a point-wise feed-forward layer without a mean-reduction
operation. Hence, the outputs of the last layer are used to determine (with a
threshold value) whether the corresponding residues are binding site residues
or not.

###  4.4. Transfer learning

Transfer learning can be applied between the BSD and BRI sub-tasks thanks to
the shared architectures between the BSD and BRI modules. More specifically,
we initialize the weights of the BSD module’s shared parameters with the
weights obtained from BRI module’s training. The rationale behind this
procedure is the following intuition: the protein’s binding site can be
determined based on the patterns of the binding residues. Under this
rationale, we hypothesize that a well-performing BRI module will learn useful
features that can transfer well to the BSD task.

In addition, the transfer learning allows BSD module to leverage relatively
more abundant labels present in BRI dataset. While there is one label per
binding site for the BSD task, there are a multiple number of binding residues
per binding site. Thus, it is desirable to exploit such abundance in labels
from the BRI task for BSD task via the transfer learning.

###  4.5. The geometric self-attention layers

As mentioned previously, we adopted the geometric self-attention layers to
globally aggregate local features obtained by the CNN. First, let
$\left\\{x_{i}\right\\}{i=1}^{n}$ denote a sequence of hidden vectors that is
either the initial local features computed by the CNN or the output of the
preceding attention layer. Then a geometric self-attention layer $f_{att}$
transforms $\left\\{x_{i}\right\\}{i=1}^{n}$ into another sequence of hidden
vectors $\left\\{x_{i}^{\prime}\right\\}{i=1}^{n}$ based on the protein
structure. In our method, we represent the protein structure as local frames
([JEP+21])
$\left\\{T_{i}\right\\}_{i=1}^{n}=\left\\{(R_{i},t_{i})\right\\}_{i=1}^{n}$
associated to the residues, where $R_{i}$ is the residue orientation and
$t_{i}$ is the residue center. Then the geometric self-attention layer
$f_{att}$ takes the following form:

(4.1)
$\left\\{x_{i}^{\prime}\right\\}_{i=1}^{n}=f_{att}\left(\left\\{x_{i}\right\\}_{i=1}^{n},\left\\{T_{i}\right\\}_{i=1}^{n}\right)$

The geometric attention mechanism has several advantages over the modules from
previous works used for the feature aggregation—the traditional clustering
algorithms and the CNN models. First, unlike the Neural Network based methods
(such as our), the clustering algorithms are less flexible in terms of the
number of adjustable parameters and they do not fit their parameters based on
the gradients. As for the CNN models, the attention layers are arguably more
effective in terms of modeling the long-distance dependency. While an
attention layer can emulate arbitrarily distant interaction by a single step
of operation, convolution layers require several steps to do so due to their
local nature. This long-term dependency may lead to sub-optimal learning
outcomes ([HBF+01]).

For the attention mechanism, we used a modified version of Invariant Point
Attention from Alphafold ([JEP+21]). The essential ideas of its computation,
compared to the standard attention mechanism ([VSP+17]), are that (1) it uses
not only the standard query, key and value vectors but also geometric ones,
(2) it calculates the attention weights based on not only the inner products
of the standard query and key vectors but also the distances between the
geometric query and key vectors, and (3) its output is determined by not only
the aggregated standard value vectors but also the aggregated geometric value
vectors. More detailed descriptions are provided in Figure 2 and the Methods
section. The effect of using the geometric vectors in the attention mechanism
is discussed as a part of our ablation study.

###  4.6. Satisfaction of SE(3)-invariance

A function is SE(3)-invariant if its output remains unchanged when SE(3)
transformations (translations, rotations or compositions of them) are applied
to the input. In our context, when $\left\\{v_{i}\right\\}_{i=1}^{N}$ are the
coordinates of the protein atoms, $\left\\{f_{i}\right\\}_{i=1}^{N}$ are the
feature vectors of the protein atoms, and $\left\\{T_{i}\right\\}_{i=1}^{n}$
are the local frames related to the protein residues, a function
$f(\left\\{v_{i}\right\\}_{i=1}^{N},\left\\{f_{i}\right\\}_{i=1}^{N}\left\\{T_{i}\right\\}_{i=1}^{n})$
is SE(3)-invariant if, for any SE(3) transformation $T$, we have

$f(\left\\{Tv_{i}\right\\}_{i=1}^{N},\left\\{f_{i}\right\\}_{i=1}^{N},\left\\{TT_{i}\right\\}_{i=1}^{n})=f(\left\\{v_{i}\right\\}_{i=1}^{N},\left\\{f_{i}\right\\}_{i=1}^{N},\left\\{T_{i}\right\\}_{i=1}^{n})$

The SE(3)-invariance is a desired property for the structure-based BSP models.
This is because the binding site information should remain unchanged
regardless of the reference frame. To reflect this, we are injecting this
inductive bias via incorporating the SE(3)-invariance property into our model
and hence achieve robustness in our model’s prediction ([FWFW20]).

Two strategies were employed to achieve SE(3)-invariance in our model. First,
we adopted SE(3)-invariant attention mechanism at the outset. To be more
specific, we required the attention layer in (4.1) to be invariant to any
SE(3) transformation $T=(R,t)$, that is:

(4.2)
$f_{att}(\left\\{x_{i}\right\\}_{i=1}^{n},\left\\{TT_{i}\right\\}_{i=1}^{n})=f_{att}(\left\\{x_{i}\right\\}_{i=1}^{n},\left\\{T_{i}\right\\}_{i=1}^{n}).$

The details of the SE(3)-invariant attention mechanism will be explained in
the Methods section. Then, the other strategy employed to meet the
SE(3)-invariance was altering the grid-based residue featurization process.
When constructing the grids, we did not use an arbitrary (xyz-) axes, but
rather used axes aligned with respect to the orientations of the residues.
This will be also explained in more detail in the Methods section.

###  4.7. Augmentation strategies

Data augmentation is one of the most essential elements comprising the modern
Deep Learning. It plays a crucial role in (1) alleviating the problem of
overfitting which occurs commonly due to the lack of data and (2) improving
the trained model’s generalization capability ([SK19]). In essence, the data
augmentation enlarges the effective size of the training set by applying
various transformations to the inputs.

Previous Deep Learning methods for BSP ([JDMR+17], [AGC+21], [SDZS20])
employed various augmentations. These augmentation methods were mainly based
on applying transformations in the SE(3) class (rotations, translations and
compositions of them) to the inputs of the CNN models.

In the absence of data augmentation, our model suffered from a clear pattern
of over-fitting (Figure 4). We believe that this is due to the lack of
diversity in the training dataset as compared to the large complexity of the
model.

However, the existing augmentation methods based on SE(3) transformations have
no effect on the training dynamics when used in conjunction with our
SE(3)-invariant model. From this aspect, we devised two novel augmentation
methods: (1) a method based on a class of geometric transformations and (2) a
method based on the homologous proteins.

Firstly, the augmentation based on the geometric transformation introduces the
random perturbations to the residue orientations. This differs to the usual
“random rotation” in that it permits only rotations of small magnitude. We
applied independent random perturbations on the residue orientations once then
used the modified residue orientations as an input to every attention layer.
The rationale behind applying such random perturbation is to promote diversity
in geometric information while not completely forgetting the original residue
orientation.

Next, we employed an augmentation scheme based on the protein homology. This
fundamentally differs from the usual augmentation methods, in that it relies
on an unlabelled external database of unbound protein structures to pre-
compute the transformations. Essentially, we first align the protein sequences
in the training set (referred to as a internal protein sequence) with the
protein sequences from the external database (referred to as a external
protein sequence). Then, we assign the ground-truth binding site and binding
residue labels of the internal protein sequence to the aligned external
protein sequence as its labels. Finally, we augment our training dataset with
these aligned external protein sequences and their labels assigned. The
extended dataset may contain some label noises since the assignment of labels
based on the homology relations is subject to inaccuracies to some extent. In
spite of this problem, our augmentation scheme enhances diversity of input
data compared to the previous augmentation schemes. This is because the
extended dataset may include protein species absent in the training dataset.

##  5\. Experiments

###  5.1. Datasets

We used scPDB v.2017 ([DBRK15]) as a main dataset for the training and the
validation. In addition, we used three other datasets for the tests: COACH420
([KH18]), HOLO4K ([KH18]) and CHEN ([CMGK11]). To be more specific, we used
the training subset of scPDB dataset provided by [MAD21] for 5-fold cross-
validation. This subset excludes proteins which have sequence similarity
higher than 90% to the proteins in one of the external test datasets. We used
the remaining part of the scPDB as a test dataset. Thus, the test datasets
were comprised of the scPDB test set and the external test dataests —
COACH420, HOLO4k and CHEN. The CHEN dataset had holo and apo subsets. Thus in
the tests using the apo subset, we obtained the ground-truth binding sites
from the structural alignments with the corresponding holo structures. More
specifically, the ligands of the holo structures were superimposed onto the
apo structures according to the structural alignments. The characteristics of
each test dataset and the details of the structural alignments are described
in the Supplementary Information.

###  5.2. Baseline Methods

We compared our method to the previous state-of-the-art Deep Learning methods,
which are based on CNN: Deeppocket ([AGC+21]), Kalasanty ([SDZS20]) and
DeepSurf ([MAD21]). All these methods are briefly explained in 2.2. For
Deeppocket and DeepSurf, we have trained their parameters from scratch
according to our dataset splits. However for Kalasanty, we used the parameters
released by the authors due to the high computational costs to train. It is
important to note that the parameters of Kalasanty were trained on the entire
scPDB v.2017 dataset. To be specific, the training data of Kalasanty may have
included data whose protein sequences are similar (similarity above 90%) with
those in the test dataset. Thus, Kalasanty method have an advantage in terms
of the coverage of the training dataset compared to the other methods when
they are evaluated on the external test datasets.

###  5.3. Ablated Methods

We conducted an ablation study to assess the effectiveness of each component
of our proposed method. We considered the omission of the following
components:

  * •

The geometric information used in the self-attention layers

  * •

The use of local features extracted by the CNN

  * •

The grid alignment process adopted to achieve SE(3)-invariance in our model

  * •

The transfer learning from the BRI to the BSD module

  * •

The data augmentation methods as a whole

  * •

Each of the data augmentation methods

We used the most basic form of transformers, that is BERT ([VSP+17]), as a
geometric information omitted version of self-attention layers to conduct the
ablation study on ‘the geometric information used in self-attention layers’.

In the ablation study of ‘the use of local features extracted by the CNN’, we
removed the CNN component from our model. To be more specific, the hidden
vectors for the attention layers are directly obtained from the One-hot-
encoding layer for the amino acid types followed by the token embedding layer.
To compensate for the loss of model complexity, we added two more attention
layers to the default configuration of BRI and BSD model architecture.

In the ablation study on the grid alignment process, we removed the alignment
process but used the standard random rotation augmentation for a fair
comparison. By doing so, we attempted to argue the efficacy of our method as
compared to the common practice, rather than the case where no additional
technique is applied.

###  5.4. Evaluation

In 3, we formalized the inputs and outputs of the BSP methods. According to
the formalism, a BSP method (1) takes as inputs a protein structure and the
number of binding sites ($n$), and (2) outputs $n$ predictions
$(\hat{c}_{1},\hat{R}_{1}),\cdots,(\hat{c}_{n},\hat{R}_{n})$ with
$\hat{c}_{i}$ and $\hat{R}_{i}$ being the binding site centers and the binding
residue indices, respectively. Now that we have formalized the problem, we
need to establish an evaluation scheme. This requires (1) interpreting each
method’s input and output as per the problem definition and (2) defining the
evaluation metrics in terms of the input and output formats specified in the
problem definition.

####  5.4.1. How the baseline methods fit into our problem definition

The definitions of $\hat{c}_{1},\cdots,\hat{c}_{n}$ for the baseline methods
are mostly natural, and derived directly from the original papers. All methods
produce a ranked list of predictions, thus we limit them to produce only the
top-$n$ outputs. Also, they compute the centers of predicted binding sites in
their evaluation, so we can compute $\hat{c}_{i}$ as they prescribed.

However, not all baseline models output the predicted binding sites at the
residue level. Thus, it is necessary to map their outputs to sets of residues
$\hat{R}_{1},\cdots,\hat{R}_{n}$. For example in the Deeppocket ([AGC+21]),
the authors used the distance threshold $2.5\AA$ (performed best in their
validation set) to determine the binding residues from the segmented voxels;
therefore, we followed the same procedure. For Kalasanty ([SDZS20]) and
DeepSurf ([MAD21]), the authors introduced a method to convert their
predictions to the atom-level predictions (which was implemented in their
code); therefore, we regarded the residues having at least one such predicted
binding atoms as the predicted binding residues.

####  5.4.2. The evaluation metrics

We use three evaluation metrics: (1) the success rate for detection (success
rate) (2) the average IOU of binding residues with respect to the closest
ligands (IOU) (3) the average IOU of the binding residues with respect to the
successfully detected ligands (conditional IOU). The metrics evaluate
different combinations of BSD and BRI performances. The success rate metric
evaluates the BSD performance, the IOU metric evaluates the BRI performance
but it is also influenced by the BSD performance, and the conditional I metric
aims to evaluate the BRI performance alone. We give additional details with
regards to how these evaluation metrics compare to their counterpart metrics
introduced in the previous literatures in Supplementary Information.

In order to provide a formal definition of each metric, we shall adopt
following notations:

  * •

$n^{(i)}$ is the number of ground-truth ligands bound to the $i$-th protein.

  * •

$\left\\{l^{(i)}_{1},\cdots,l^{(i)}_{n^{(i)}}\right\\}$ is the set of ground-
truth ligands bound to the $i$-th protein.

  * •

$\left\\{(c^{(i)}_{1},BR^{(i)}_{1})\cdots,(c^{(i)}_{n^{(i)}},BR^{(i)}_{n^{(i)}})\right\\}$
is the set of predictions of the method to evaluate.

  * •

$\left\\{TBR^{(i)}_{1},\cdots,TBR^{(i)}_{n^{(i)}}\right\\}$ is the set of true
binding residue indices, where $TBR^{(i)}_{j}$ is defined to be the set of
residues in the $i$-th protein that is within $4\AA$ from $l^{(i)}_{j}$

The success rate metric measures the correspondence between the predicted
binding site centers $\left\\{c^{(i)}_{1},\cdots,c^{(i)}_{n^{(i)}}\right\\}$
and the positions of the ground-truth ligands
$\left\\{l^{(i)}_{1},\cdots,l^{(i)}_{n^{(i)}}\right\\}$. For each $i$, we
compute the F1 score (the harmonic mean of precision and recall) based on the
definition of detection. Specifically, we define that $c^{(i)}_{j}$ is a
correct detection of $l^{(i)}_{k}$ when $c^{(i)}_{j}$ is within $4\AA$ (a
threshold commonly used in the literature e.g. [AGC+21] and [MAD21]) from any
ligand of $l^{(i)}_{k}$. In other words, we define that detection is correctly
performed when the Distance from Center to Atom (DCA) is $<4\AA$. Then, the F1
scores are weighted-averaged (weighted by $n^{(i)}$) over the proteins. In
summary, we obtain this metric as

$\left(\sum_{i}n^{(i)}\cdot\frac{2}{\frac{1}{P^{(i)}}+\frac{1}{R^{(i)}}}\right)\bigg{/}\left(\sum_{i}n^{(i)}\right)$

, where $P^{(i)}$ is the precision defined as follows:

$P^{(i)}=\frac{\\#\left\\{1\leq j\leq n^{(i)}:c^{(i)}_{j}\text{ detects one of
}l^{(i)}_{1},\cdots l^{(i)}_{n^{(i)}}\right\\}}{n^{(i)}}$

and $R^{(i)}$ is the recall defined as follows:

$R^{(i)}=\frac{\\#\left\\{1\leq k\leq n^{(i)}:l^{(i)}_{k}\text{ is detected by
one of }c^{(i)}_{1},\cdots,c^{(i)}_{n^{(i)}}\right\\}}{n^{(i)}}$

This is a BSD metric since this involves only the predicted binding site
centers, not the predicted binding residues.

The IOU metric compares the predicted binding residues $BR^{(i)}_{j}$ with the
true binding residues $TBR^{(i)}_{\phi^{(i)}(j)}$ of the ligand
$l^{(i)}_{\phi^{(i)}(j)}$ closest to the predicted binding site center
$c^{(i)}_{j}$. Here, the index $\phi^{(i)}(j)$ of the closest ligand is
defined as

$\phi^{(i)}(j)=\operatorname*{argmin}_{k=1}^{n^{(i)}}DCA(c^{(i)}_{j},l^{(i)}_{k})$

The comparison is performed in terms of Intersection Over Union (IOU), and the
quantity is averaged over all pairs of $(i,j)$. In summary, we obtain the
second metric as

$\left(\sum_{i}\sum_{j=1}^{n^{(i)}}\frac{\\#(BR^{(i)}_{j}\cap
TBR^{(i)}_{\phi^{(i)}(j)})}{\\#(BR^{(i)}\cup
TBR^{(i)}_{\phi^{(i)}(j)})}\right)\bigg{/}\left(\sum_{i}n_{i}\right)$

Although this is essentially a BRI metric, it also depends on the BSD
performance due to the definition of $\phi^{(i)}(j)$. In particular, if the
predicted center $c^{(i)}(j)$ is far from any ligand, the set of predicted
binding residues $\hat{R}_{j}^{(i)}$ does not contribute to the metric.

The conditional IOU metric is almost the same as the IOU metric, but it aims
to eliminate the previously mentioned problem that the BRI performance is
bound by the BSD performance. It does so by focusing on the case that the
predicted binding sites is close to at least one ligand. In summary, we obtain
the metric as

$\left(\sum_{i}\sum_{j\in S^{(i)}}\frac{\\#(BR^{(i)}_{j}\cap
TBR^{(i)}_{\phi^{(i)}(j)})}{\\#(BR^{(i)}\cup
TBR^{(i)}_{\phi^{(i)}(j)})}\right)\bigg{/}\left(\sum_{i}\\#S^{(i)}\right)$

, where

$S^{(i)}=\left\\{j=1,\cdots,n^{(i)}:DCA(c^{(i)}_{j},l^{(i)}_{\phi^{(i)}(j)})<4\AA\right\\}$

This is similar to metrics used in Deeppocket ([AGC+21]) and DeepSurf
([MAD21]) to evaluate models’ binding site segmentation capability conditional
on the successful location of the binding sites.

###  5.5. Training the BSD Module

To implement the transfer learning described in 4.4, the BSD training consists
of two stages. The first stage is pre-training the part of the BSD module’s
architecture shared by the BRI module, as depicted in Figure 1. In doing so,
we attach the unshared part of the BRI module architecture on top of the
shared part in the BSD module then train the combined model for the BRI task.
The second stage is fine-tuning the entire original BSD module for the BSD
task. In the second stage, to promote smoother transfer learning, we freeze
the parameters of the parts trained in the first stage up to certain steps of
gradient descent.

In both stages, we use balanced sampling of binding site candidates of
positive and negative labels. This is because, among the binding site
candidates predicted by Fpocket (On average 33 per protein in scPDB), only few
of them are actual binding sites (typically only one). Training without such
balanced sampling may render the trained model biased toward the out-numbered
label. ([JK19])

In addition, in the first stage, we resolve the similar problem of unbalanced
residue labels by using a weighted loss function. This loss function consists
of a weighted sum of terms from different residues, where the binding and non-
binding residues attain the following weights:

$w_{pos}=\frac{1}{2n_{pos}},\quad w_{neg}=\frac{1}{2n_{neg}}$

, where $n_{pos}$ and $n_{neg}$ are the number of binding and non-binding
residues respectively.

###  5.6. Training the BRI module

To train the BRI module, we use only the positive binding site candidates. (We
do the same to train Deeppocket) This is because, in our intended usage, the
BRI module operates on the binding sites detected by the BSD module. Note that
this intention is reflected in the evaluation metric average IOU of binding
residues against the closest ligands as well. All the settings of the first
stage of BSD training were maintained, except the balanced sampling of the
binding site candidates.

###  5.7. Data Augmentations

We apply the data augmentation methods described in 4.7 in every training
scenario (pre-training the BSD module, fine-tuning the BSD module and training
the BRI module), except when they are omitted as a part of the ablation study.
Applying the random perturbation means transforming the inputs to the
geometric attention layers once for each data sample. Applying the homology-
based “augmentation” means adding to the original loss (originated from the
original dataset) an auxiliary loss calculated in the same way but originated
from the “extended dataset”.

###  5.8. Experiment Results

Table 1. (BSD metric) F1 Success rate for detection. The mean and standard deviation are calculated based on the metric values for 5 different cross-validation folds. | scPDB(held-out) | COACH420 | HOLO4k | CHEN-holo | CHEN-apo
---|---|---|---|---|---
DeepSurf | $62.4\pm 1.3$ | $43.6\pm 1.3$ | $59.7\pm 1.4$ | $24.5\pm 1.4$ | $22.3\pm 1.2$
Kalasanty | $70.0\pm 0.0$ | $50.8\pm 0.0$ | $44.9\pm 0.0$ | $28.5\pm 0.0$ | $27.1\pm 0.0$
DeepPocket | $67.9\pm 0.4$ | $55.7\pm 0.7$ | $72.2\pm 0.2$ | $\bm{42.4\pm 0.3}$ | $34.5\pm 1.2$
Ours | $\bm{70.1\pm 0.4}$ | $59.1\pm 0.3$ | $\bm{77.0\pm 0.6}$ | $41.2\pm 1.2$ | $36.5\pm 0.2$
Ours(BERT) | $65.6\pm 1.2$ | $54.8\pm 1.2$ | $70.8\pm 1.1$ | $38.2\pm 1.2$ | $32.8\pm 0.7$
Ours(no CNN) | $69.2\pm 0.7$ | $57.7\pm 0.3$ | $75.1\pm 0.2$ | $41.7\pm 1.0$ | $\bm{36.7\pm 1.0}$
Ours(no alignment) | $69.0\pm 0.6$ | $\bm{59.4\pm 0.7}$ | $75.6\pm 0.7$ | $41.0\pm 0.6$ | $35.4\pm 0.6$
Ours(no transfer) | $61.9\pm 0.3$ | $53.9\pm 0.3$ | $70.0\pm 0.5$ | $39.8\pm 0.5$ | $35.1\pm 0.6$
Ours(no augmentation) | $64.2\pm 2.0$ | $55.8\pm 1.4$ | $71.3\pm 1.5$ | $39.1\pm 1.7$ | $34.3\pm 0.5$
Ours(no homology) | $68.3\pm 1.1$ | $57.3\pm 0.9$ | $75.1\pm 1.1$ | $39.9\pm 0.5$ | $34.6\pm 0.6$
Ours(no perturbation) | $70.0\pm 0.5$ | $58.9\pm 0.7$ | $76.2\pm 0.3$ | $41.8\pm 0.7$ | $36.3\pm 1.1$
Table 2. (BSD + BRI metric) Average IOU of binding residues against the closest ligands. The mean and standard deviation are calculated based on the metric values for 5 different cross-validation folds. | scPDB(held-out) | COACH420 | HOLO4k | CHEN-holo | CHEN-apo
---|---|---|---|---|---
DeepSurf | $0.288\pm 0.007$ | $0.194\pm 0.006$ | $0.207\pm 0.005$ | $0.104\pm 0.003$ | $0.085\pm 0.005$
Kalasanty | $0.260\pm 0.000$ | $0.183\pm 0.000$ | $0.146\pm 0.000$ | $0.101\pm 0.000$ | $0.092\pm 0.000$
Deeppocket | $0.440\pm 0.002$ | $0.313\pm 0.003$ | $0.277\pm 0.03$ | $0.190\pm 0.005$ | $0.186\pm 0.003$
Ours | $\bm{0.490\pm 0.003}$ | $\bm{0.398\pm 0.004}$ | $\bm{0.346\pm 0.002}$ | $\bm{0.287\pm 0.004}$ | $\bm{0.264\pm 0.004}$
Ours(BERT) | $0.430\pm 0.008$ | $0.321\pm 0.010$ | $0.294\pm 0.005$ | $0.228\pm 0.007$ | $0.206\pm 0.005$
Ours(no CNN) | $0.467\pm 0.003$ | $0.357\pm 0.005$ | $0.315\pm 0.002$ | $0.247\pm 0.005$ | $0.227\pm 0.005$
Ours(no alignment) | $0.478\pm 0.005$ | $0.387\pm 0.004$ | $0.335\pm 0.003$ | $0.271\pm 0.005$ | $0.242\pm 0.004$
Ours(no augmentation) | $0.381\pm 0.003$ | $0.317\pm 0.006$ | $0.275\pm 0.004$ | $0.229\pm 0.007$ | $0.211\pm 0.003$
Ours(no homology) | $0.473\pm 0.003$ | $0.387\pm 0.007$ | $0.341\pm 0.003$ | $0.285\pm 0.010$ | $0.257\pm 0.004$
Ours(no perturbation) | $0.452\pm 0.002$ | $0.370\pm 0.007$ | $0.320\pm 0.005$ | $0.269\pm 0.005$ | $0.240\pm 0.006$
Table 3. (BRI metric) Average IOU of binding residues against the detected ligands. The mean and standard deviation are calculated based on the metric values for 5 different cross-validation folds. | scPDB(held-out) | COACH420 | HOLO4k | CHEN-holo | CHEN-apo
---|---|---|---|---|---
DeepSurf | $0.402\pm 0.010$ | $0.419\pm 0.013$ | $0.330\pm 0.007$ | $0.372\pm 0.019$ | $0.336\pm 0.020$
Kalasanty | $0.356\pm 0.000$ | $0.362\pm 0.000$ | $0.344\pm 0.000$ | $0.333\pm 0.000$ | $0.323\pm 0.000$
Deeppocket | $0.595\pm 0.002$ | $0.506\pm 0.005$ | $0.371\pm 0.004$ | $0.395\pm 0.008$ | $0.382\pm 0.006$
Ours | $\bm{0.643\pm 0.004}$ | $\bm{0.585\pm 0.007}$ | $\bm{0.415\pm 0.002}$ | $\bm{0.495\pm 0.010}$ | $\bm{0.473\pm 0.004}$
Ours(BERT) | $0.567\pm 0.012$ | $0.459\pm 0.013$ | $0.356\pm 0.007$ | $0.353\pm 0.016$ | $0.328\pm 0.010$
Ours(no CNN) | $0.624\pm 0.002$ | $0.548\pm 0.004$ | $0.395\pm 0.002$ | $0.450\pm 0.014$ | $0.419\pm 0.006$
Ours(no alignment) | $0.637\pm 0.005$ | $0.572\pm 0.001$ | $0.413\pm 0.002$ | $0.479\pm 0.005$ | $0.450\pm 0.006$
Ours(no augmentation) | $0.522\pm 0.001$ | $0.475\pm 0.007$ | $0.347\pm 0.004$ | $0.397\pm 0.010$ | $0.380\pm 0.007$
Ours(no homology) | $0.628\pm 0.003$ | $0.567\pm 0.009$ | $0.412\pm 0.002$ | $0.481\pm 0.007$ | $0.449\pm 0.010$
Ours(no perturbation) | $0.596\pm 0.006$ | $0.547\pm 0.008$ | $0.391\pm 0.004$ | $0.468\pm 0.008$ | $0.445\pm 0.008$
Figure 3. The effect of transfer learning on the BSD training  Figure 4. The
effect of the augmentation methods on the training

The experiment results are summarized in Table 1, Table 2, Table 3, Figure 3
and Figure 4.

####  5.8.1. BSD performance

In general, Table 1 shows that our method significantly outperformed the
baseline methods in the BSD task. Although there is an exception that
Deeppocket outperformed our method on CHEN-holo, the gap in performance
($1.2\%p$) is relatively insignificant compared to the average gain in
performance ($3.1\%p$) across other datasets.

The ablation results show that all components had positive effects on the
performance in general. In particular, (1) the use of geometric information in
the attention layers, (2) the transfer learning and (3) augmentation (in
particular the homology augmentation) showed positive contributions
consistently and significantly. Omitting each component resulted in the
decrease in the model’s performance as follows:

  * •

The geometric information: $4.3\%p$

  * •

The transfer learning: $4.64\%p$

  * •

The augmentation as a whole: $3.84\%p$

  * •

The homology augmentation: $1.74\%p$

  * •

The random perturbation augmentation: $0.14\%p$

However, we observed that not all components of our proposed model showed
consistent contributions to the model’s performance across different datasets.
This was the case for the use of CNN, the grid alignment and the random
perturbation. The questionable effects of the first two components may be
accounted to the insufficient size of the training dataset for the BSD task.
We believe that the model required a larger training dataset to compensate for
the increased model complexity from the introduction of the CNN component.
Similarly, aligning the grids rather than applying a random rotation to the
grids may have decreased the diversity of the training dataset. Lastly, the
random perturbation’s questionable contribution may be attributed to the
offset in the effects from the homology augmentation. Indeed, the random
perturbation resulted in a significant improvement in performance in the
absence of the homology augmentation, which can be confirmed by comparing “no
augmentation” with “no homology” in Table 1 and Figure 4.

####  5.8.2. The effect of transfer learning

Figure 3 shows the effect of transfer learning in the BSD training. According
to the plot, the transfer learning had two effects on the training process.
Firstly, it significantly accelerated the convergence speed. In fact,
validation loss almost dropped to the convergence level in 2000 steps. Note
that until 4000 steps, we updated only the weights of the un-pretrained parts.
Secondly, the transfer learning also significantly improved the validation
loss of the converged state. This finding is compatible with the ablation
result of Table 1.

####  5.8.3. BRI performance

In terms of the BRI performance, Table 2 and Table 3 show that our model
outperformed the baselines by a significant margin on all test datasets. In
particular, while the strongest baseline model Deeppocket performed poorly on
the external test datasets, our model did not. This shows our model
generalizes well to the proteins such that no similar proteins are encountered
during training.

The ablation results show that all key components of our method explained in 4
contributed significantly in the good performance. These components are (1)
the architectural aspects (the CNN and attention layers), (2) the grid
alignment process and (3) the augmentation methods.

####  5.8.4. The effect of augmentation

Figure 4 shows that the augmentation contributed significantly in alleviating
over-fitting in all training scenarios (pre-training the BSD module, fine-
tuning the BSD module and training the BRI module). While the augmentation as
a whole was shown to dramatically reduce the over-fitting, individual
augmentation methods (random perturbation and homology augmentation) were also
shown to be effective. This finding is compatible with the ablation results of
Table 1, Table 2 and Table 3.

###  5.9. Conclusions

In general, our model significantly outperformed the baseline methods in terms
of both BSD and BRI. All key components in 4 contributed to the performance,
with some exceptions with respect to BSD. Additional analyses confirmed that
the transfer learning and the augmentation worked as intended during training.
The transfer learning accelerated and improved the convergence, and the
augmentation methods helped overcome the over-fitting problem.

##  6\. A case study

We conducted a case study to demonstrate our model’s applicability in real-
world drug discovery scenarios. We chose Human Serum Albumin (HSA) as the
target protein, due to its relevance to the drug discovery and the fact that
it binds to various molecules at different sites ([PJ95], [Car03], [PKP+21]).
To examine our model’s performance on HSA, we drew on two prior studies on the
binding sites of HSA ([Car03], [WHR+13]), and analyzed how compatible our BSD
and BRI modules’ predictions were with the findings of the studies.

The remaining parts of this section are organized as follows:

  * •

In 6.1, we briefly explain the structure of HSA, and elucidate what the
reference papers ([Car03], [WHR+13]) revealed about the binding sites of it.

  * •

In 6.2, we explain basic settings of our case study experiments.

  * •

In 6.3, we analyze our BSD module’s performance on HSA based on the findings
of [Car03].

  * •

In 6.4, we analyze our BRI module’s performance on HSA based on the findings
of [WHR+13].

###  6.1. The structure and Binding Sites of HSA

Figure 5.  The subdomains of HSA

HSA is composed of three homologous domains (I, II, III), each composed of two
subdomains A and B, as depicted in Figure 5 ([PJ95], [Car03]). The authors of
[Car03] performed a large-scale survey of HSA binding sites based on 142
crystal structures involving HSA. By analyzing the complexes, they identified
14 different binding sites, alongside their frequencies. They noted that
binding sites “IB”, “IIA” and “IIIA” dominated clearly. Inspired from this
result, the authors of [WHR+13] provided more detailed analysis on the binding
sites at the IB subdomain. They did so by inspecting the crystal sructures of
6 oncology drugs (9-amino-camptothecin, camptothecin, idarubicin, teniposide,
etoposide and bicalutamide) in complex with HSA. In particular, for each
structure, they identified all key residues of HSA and their interaction types
with the drug molecule.

###  6.2. Basic settings

We designed the analyses such that they can faithfully assess our model’s
real-world applicability.

Firstly, before the analyses, we trained our model with a new dataset split
(different from the ones used in our main experiments) to prevent data
leakage. Specifically, we removed all 42 “albumin” structures from the scPDB
v.2017 dataset, randomly sampled a validation set of size 1000 from the
remaining, and took all the other proteins as the training set. Moreover, we
ensured that there is no leakage coming from the homology-based augmentation,
by using a new augmentation dataset generated from the new training set.

Secondly, all our model’s predictions are based on taking as input the HSA
structure provided by the Alphafold database ([TAW+21]. Therefore, a good
performance under this setting would imply that our model can be used to
predict the binding sites of a protein without any known experimental
structures. In particular, one can make use of the publicly available
Alphafold database in the predictions.

###  6.3. Binding Site Detection

Figure 6.  A comparison between our BSD module’s predictions and the
frequencies from [Car03]

####  6.3.1. experiment procedure

First, we obtained our model and Deeppocket’s BSD module’s predictions on 15
binding sites identified by [Car03]. This means that the predictions were made
based on 15 different inputs, where we set the “binding site center” to be the
mean alpha carbon coordinates of residues comprising one of the binding sites.
The indices of residues comprising each binding site were provided by [Car03].
Note that we re-trained the Deeppocket model with a new dataset split to avoid
data leakage, just as we did for our model.

Then, we assessed each model’s predictions by comparing them with the ranks of
the frequencies of the binding sites as recorded in [Car03].

####  6.3.2. results and analysis

The results are summarized in Figure 6. Note that our model successfully
assigned high (more than $70\%$) druggability scores on the second to sixth
most frequent binding sites. Moreover, these five were exactly the binding
sites that scored the highest. The probability of the latter condition being
met by a random ordering is only $0.2\%$, which shows the statistical
significance of our model’s ability in replicating binding sites’ frequency
ranks. On the other hand, although Deeppocket assigned high druggability
scores on the second and third most frequent binding sites, it failed to do so
on the fourth to sixth most frequent ones.

###  6.4. Binding Residue Identification

(a) All 6 drug molecules examined in [WHR+13]

(b) Our model’s prediction

(c) Deeppocket’s prediction

(d) Kalasanty’s prediction

(e) DeepSurf’s prediction

Figure 7.  The purple, red and blue residues indicate the true-positive,
false-positive and false-negative binding site residues respectively.
Therefore, the larger the purple region is compared to the regions with the
other colors, the better the prediction is.

####  6.4.1. experiment procedure

First, we predicted the binding site residues in the IB subdomain of HSA using
our model and all the baseline methods — Deeppocket, Kalsanty and DeepSurf.
When obtaining our model’s (resp. Deeppocket’s) predictions, we set the
“binding site center” to be the mean alpha carbon coordinate of the IB
subdomain residues, and run the BRI module (resp. the segmentation model).
Note that as we did for our model, we re-trained the Deeppocket and DeepSurf
models with a new dataset split to avoid data leakeage.

Then, we visually inspected the binding residues predicted by different
methods, comparing them with the ground truth provided by [WHR+13]. As the
ground truth, we took the union of all the “key residues” (in the IB
subdomain) of the 6 drug molecules as determined by [WHR+13].

####  6.4.2. results and analysis

The results are visualized in Figure 7. From the figures, it is clear that our
model’s prediction best matched with the ground truth. In particular, while
all the baseline models resulted in false-positives within the left-bottom
helix, our model did not. Also, while all the baseline models showed limited
precision in identifying the key residues in the other parts of the IB
subdomain, our model showed mostly successful identifications.

##  7\. Methods

###  7.1. The residue orientation

Our model uses a concept of “residue orientation” in the grid alignment
process and the geometric attention layers. In the grid alignment process, the
grid axes are aligned with respect to the orientations. In the geometric
attention layers, the orientations are used as a part of the input.

We define the residue orientation in terms of the relative positions of atoms
surrounding the alpha carbon, as in [JEP+21]. More precisely, when
$\textbf{x}_{1}$, $\textbf{x}_{2}$ and $\textbf{x}_{3}$ are the coordinates of
$N$, $C\alpha$ and $C$ (the carbon that is not $C\beta$ and adjacent to
$C\alpha$), the rotation matrix
$R=(\textbf{e}_{1}\enspace\textbf{e}_{2}\enspace\textbf{e}_{3})$ that we refer
to as “orientation” is obtained as follows:

$\displaystyle\textbf{v}_{1}=\textbf{x}_{3}-\textbf{x}_{2}$
$\displaystyle\textbf{v}_{2}=\textbf{x}_{1}-\textbf{x}_{2}$
$\displaystyle\textbf{e}_{1}=\textbf{v}_{1}/\left\|\textbf{v}_{1}\right\|$
$\displaystyle\textbf{u}_{2}=\textbf{v}_{2}-(\textbf{e}_{1}\cdot\textbf{v}_{2})\textbf{e}_{1}$
$\displaystyle\textbf{e}_{2}=\textbf{u}_{2}/\left\|\textbf{u}_{2}\right\|$
$\displaystyle\textbf{e}_{3}=\textbf{e}_{1}\times\textbf{e}_{2}$

###  7.2. The grid featurization

As explained in 4.2 and 4.3, our BSD and BRI modules take as input a sequence
of 3D voxelized images, where each image represents the local environment of a
residue. The process of making such voxelized images, the “grid
featurization”, proceeds as follows:

  1. (1)

Collect coordinates $\textbf{x}_{i}\in\mathbb{R}^{3}$ ($1\leq i\leq n$) of the
heavy (non-hydrogen) atoms in the protein.

  2. (2)

Collect feature vectors $f_{i}\in\mathbb{R}^{d_{feature}}$ ($1\leq i\leq n$)
of the heavy atoms.

  3. (3)

Given a choice of grid axes $(\textbf{e}_{1},\textbf{e}_{2},\textbf{e}_{3})$,
a cubical grid is laid such that the centers of the voxels in the grid become

$t+\sum_{i=1}^{3}r(a_{i}-\frac{L-1}{2})\textbf{e}_{i}\quad\text{($(a_{1},a_{2},a_{3})\in\left\\{0,\cdots,L-1\right\\}^{3}$)}$

, where $t$ is the grid center (the alpha carbon coordinate), $r$ is the grid
resolution and $L$ is the grid size.

  4. (4)

Compute feature vectors corresponding to each voxel in the grid, by summing up
those of nearby heavy atoms.

The same process was used in [AGC+21], [MAD21] and [SDZS20], although the
grids in those methods are not laid on the protein residues.

The specific details of our grid-featurization process are:

  * •

The grid resolution is $1\AA$ and the grid size is $16$.

  * •

The atom features are of dimension 18 and include atom types, hybridization,
degree, partial charge and aromaticity. ([SDZS20], [AGC+21].)

###  7.3. The grid alignment

As explained in 4.6, we apply the “grid alignment” during the grid
featurization to promote SE(3)-invariance. This means that we choose the grid
axes to be the vectors comprising the residue orientation (See 7.1).

###  7.4. The CNN model

Our model incorporates a 3D CNN to encode the residue-level grid features. For
the CNN architecture, it uses a 3D BottleNeck ResNet model introduced in
[MAD21]. The model is adapted from the bottleneck ResNet model introduced in
[HZRS16] for the image classification. The bottleneck architecture reduces the
number of parameters thus enables the use of a deeper network. [MAD21] showed
that the 3D BottleNeck ResNet model, while being light, performed
competitively well compared to its non-bottleneck counterpart.

###  7.5. The geometric attention layers

The geometric attention is an integral part of our model’s architecture. The
attention mechanism was adopted from Alphafold([JEP+21])’s IPA(Invariant Point
Attention) with an adjustment necessary to adapt it to our forms of inputs.

The inputs of the attention layers are composed of the following:

  * •

$x_{i}\in\mathbb{R}^{d_{hidden}}$ ($i=1,\cdots,n$), hidden vectors associated
to the residues.

  * •

$T_{i}=(R_{i},t_{i})\in SO(3)\times\mathbb{R}^{3}$ ($i=1,\cdots,n$), the local
frames associated to the residues, where $t_{i}$ is the position of the alpha
carbon and $R_{i}$ is the rotation matrix that represents the residue
orientation (See 7.1). Note that the operation $v\mapsto T_{i}v$ maps the
local coordinates (with respect to the local frame) to the corresponding
global coordinates, and the operation $u\mapsto T_{i}^{-1}u$ plays the reverse
role.

Then, the computation is carried out in the following steps (See also Figure
2):

  1. (1)

The standard query and key vectors $q_{i}^{h}$ and $k_{i}^{h}$ are computed by
the linear mappings from $x_{i}$. Here, $h$ stands for a “head”.

  2. (2)

The geometric query and key vectors $\textbf{q}_{i}^{hp}$ and
$\textbf{k}_{i}^{hp}$ in $\mathbb{R}^{3}$ are computed by the linear mappings
from $x_{i}$. Here, $h$ stands for a “head” and $p$ stands for a “point” of
attention.

  3. (3)

The attention weight from the $i$-th token to the $j$-th token is computed
from a linear combination of the standard attention weight

$w_{ij}^{h,standard}=\frac{1}{\sqrt{d_{hidden}}}q_{i}^{h}\cdot k_{j}^{h}$

and the geometric attention weight

(7.1)
$w_{ij}^{h,geometric}=\frac{1}{\sqrt{N_{points}}}\sum_{p}\left\|T_{i}\textbf{q}_{i}^{hp}-T_{j}\textbf{k}_{j}^{hp}\right\|$

by applying a softmax operation. More precisely, the attention weight becomes

$w^{h}_{ij}=softmax_{j}(\frac{1}{\sqrt{2}}(w_{ij}^{h,standard}-\log(1+\gamma^{h})w_{ij}^{h,geometric}))$

where $\gamma^{h}$ is a learnable parameter.

  4. (4)

The standard value vectors $v_{j}^{h}$ are computed by a linear map from
$x_{j}$, and aggregated as

$o_{i}^{h}=\sum_{j}w_{ij}^{h}v_{j}^{h}$

  5. (5)

The geometric value vectors $\textbf{v}_{j}^{h}$ are computed by a linear map
from $x_{j}$, and aggregated as

(7.2)
$\textbf{o}_{i}^{hp}=T_{i}^{-1}(\sum_{j}w_{ij}^{h}T_{j}\textbf{v}_{j}^{hp})$

  6. (6)

The aggregated vectors as well as their sizes are concatenated and linearly
mapped via $f_{final}$ to produce the output of the attention layer

$x_{i}^{\prime}=f_{final}(concat_{h,p}(o_{i}^{h},\textbf{o}_{i}^{hp},\left\|\textbf{o}_{i}^{hp}\right\|))$

The adjustment made to the original IPA is the omission of the “attention
bias” term. In the original paper, this term was based on the “pair
representation” computed at the earlier stages of the Alphafold architecture
using the evolutionary information of the protein. Since our model does not
involve this representation, the omission is necessary to use the IPA in our
model.

###  7.6. SE(3)-invariance of the geometric attention layers

As noted in 4.6, our model’s SE(3)-invariance relies on the attention layer’s
SE(3)-invariance. In order to prove that the layer is SE(3)-invariant, one has
to examine the equality (4.2). Essentially, the equality holds for our
attention layer because the quantities (7.1) and (7.2), as functions of
$\left\\{x_{i}\right\\}_{i=1}^{n}$ and $\left\\{T_{i}\right\\}_{i=1}^{n}$,
satisfy the same equality. For any $T\in SE(3)$, the (7.1) satisfies the
equality because
$\left\|TT_{i}\textbf{q}_{i}^{hp}-TT_{j}\textbf{k}_{j}^{hp}\right\|=\left\|T(T_{i}\textbf{q}_{i}^{hp}-T_{j}\textbf{k}_{j}^{hp})\right\|=\left\|T_{i}\textbf{q}_{i}^{hp}-T_{j}\textbf{k}_{j}^{hp}\right\|$
(since $T$ preserves the norm), and (7.2) satisfies the equality because
$(TT_{i})^{-1}(\sum_{j}w_{ij}^{h}(TT_{j})\textbf{v}_{j}^{hp})=T_{i}^{-1}T^{-1}T(\sum_{j}w_{ij}^{h}T_{j}\textbf{v}_{j}^{hp})=T_{i}^{-1}(\sum_{j}w_{ij}^{h}T_{j}\textbf{v}_{j}^{hp})$.
Note that this derivation is almost identical to that presented in
[JEP+21](Supplementary Information, page 28).

###  7.7. The Homology-based Augmentation

Figure 8.  These figures illustrate the way our homology-based augmentation
determines the positive and negative binding site candidates in augmented
proteins. Figure (A) depicts an augmented protein, where there are two
positive (the blue points) and one negative (the red point) binding site
candidate centers. Out of binding site candidates proposed by Fpocket, they
are labeled based on the distances to the proxy centers (the purple points) of
binding sites inferred from homology relations. Figures (B) and (C) depict the
homologous proteins in the original database that gave rise to the inferred
binding sites in the augmented protein. $X$ and $Y$ are their ligands. The
bright and dark green regions of the chains indicate the residues in close
proximity to the ligands, while only the bright green region has evolutionary
correspondence to residues in the augmented protien. The bright green region
must comprise at least 50% of the entire green region in order for the binding
site to count.

(a) structure of UniProt Q9VC32 from the Alphafold database

(b) structure of PDB 4G34 from scPDB

(c) structure of PDB 4BID from scPDB

Figure 9.  These figures illustrate the way our homology-based augmentation
assigns residue labels in augmented proteins with respect to a positive
binding site candidate. Figure (A) illustrates UniProt protein Q9VC32 as an
augmented protein. The red and purple residues correspond to the red residues
in figure (B) of the PDB protein 4G34, which are the ligand-binding residues.
Similarly, the blue and purple residues correspond to the blue residues in
figure (C) of the PDB protein 4BID, which are the ligand-binding residues. The
purple residues, the intersection, attain labels 1.0, while the other colored
residues attain labels 0.5. This means that our augmentation method regards
the purple residues as the most likely ligand-binding ones.

The “homology-based augmentation” is one of our new augmentation methods used
to overcome the problem of over-fitting. It is distinguished from the usual
augmentation methods in that it is not based on transformations applied to the
samples on the fly during the training. Instead, it pre-computes appropriate
“augmented samples” out of an external database of unlabelled protein
structures, and use the augmented dataset consisting of the augmented samples
during training. Essentially, the augmented samples are selected based on the
sequence alignments computed with respect to the proteins in the original
training set. In this section, we describe the augmentation method in more
detail, clarifying its inputs, outputs, the procedures and the underlying
rationale.

The augmentation method requires a seed database $\mathcal{S}^{*}$ of multi-
chain protein-ligand complexes and a target database $\mathcal{T}$ of single-
chain protein structures. In our instantiation, $\mathcal{S}^{*}$ was the
portion of the scPDB dataset used to train the current fold of the cross-
validation. For $\mathcal{T}$, we used the entire Alphafold Protein Structure
Database (as of April. 2022) that contained 992,316 protein structures from
the proteome of human and 47 other key organisms as well as the Swiss-Prot
entries.

The augmentation procedure outputs two types of information, which together
form the “augmented dataset” and are used during the training as described in
5.7. The first information denotes the centers of the binding site candidates
in proteins in a selected subset of $\mathcal{T}$, labeled either positive or
negative. This is used to augment the BSD train dataset. The second
information denotes, for each previous positive binding site candidates, the
likelihood of each nearby protein residue to be a ligand-binding residue. This
is used to augment the BRI train dataset.

In the following, we describe the steps of the procedure. The italized words
are general terms whose specification may vary depending on own’s needs.
Whenever there is an italized word, we provide our specification at the end of
the step.

  * (1)

In each holo structure of $\mathcal{S}^{*}$, find ligands associated to
exactly one chain. As a result, obtain a database $\mathcal{S}$ of protein
chains associated to at least one such single-chain ligands. (A chain can be
associated to multiple single-chain ligands) We define that a chain and a
ligand are associated to each other if they have heavy atoms within $4\AA$ to
each other.

  * (2)

Run a homology search algorithm with $\mathcal{S}$ as the query database and
$\mathcal{T}$ (the database of single-chain protein structures) as the target
database. Based on the results, obtain a MSA for each chain in $\mathcal{S}$.
For the homology search algorithm, we use the software HHBlits with its
default setting.

  * (3)

For each triplet $(x,l,y)$, composed of:

    1. (1)

a query chain $x$ in $\mathcal{S}$

    2. (2)

a ligand $l$ associated to $x$ found in step 1 of the procedure and

    3. (3)

a target chain $y$ aligned with $x$ in the MSA,

determine whether the ligand $l$’s binding site in $x$ is preserved in $y$.
The triplets for which the previous determination was affirmative will be
called preserving. We define that a triplet $(x,l,y)$ is preserving if at
least half of the residues of $x$ that are in close contact with $l$ (heavy
atoms within $4\AA$) are aligned with a residue of $y$ in the MSA.

  * (4)

For each preserving triplet $(x,l,y)$, find a proxy center of the binding site
in $y$ that corresponds to the ligand $l$’s binding site in $x$. We define the
proxy center to be the mean of the alpha carbon coordinates of the residues of
$y$ aligned in the MSA with a residue of $x$ in close contact with $l$.

  * (5-1)

On each chain $y$ in $\mathcal{T}$ that is involved in at least one preserving
triplet, run Fpocket to get an initial list of binding site center candidates.
Label a candidate center “positive” if it is within a lower threshold from a
proxy center obtained in the previous step. Label it “negative” if it is
farther than a higher threshold from any such proxy center. If a candidate
center does not fall into these categories, ignore it and exclude it from the
dataset. We define the lower threshold to be $7.5\AA$ and the upper threshold
to be $30\AA$. Figure 8 illustrates this step using schematic figures.

  * (5-2)

For each positively labeled binding site candidate from the previous step,
label residues of $y$ with the estimated likelihood of comprising the binding
site. The estimate is obtained as a result of “voting” of the homologous
chains in $\mathcal{S}$ that gave rise to the binding site. More specifically,
among the preserving triplets $(x,l,y)$ whose proxy center gave rise to the
binding site (in the sense of the step (5-1)), the proportion of such triplets
for which the residue at hand corresponds (in MSA) to a residue in the binding
site of $l$ is computed. Figure 9 illustrates this step using an actual
example.

The assignments of different labels in the previous procedure are based on the
following hypotheses:

  * •

The positive binding site label: if a pocket-like site (discovered by a
geometry-based BSP method) is surrounded by the sequence fragments that are
homologous to the binding site sequence fragments of other proteins, it is
likely to be a binding site.

  * •

The negative binding site label: even if a site is pocket-like, if it is far
from any sequence fragments that are homologous to the binding site sequence
fragments of other proteins (in a given seed database), it is unlikely to be a
binding site.

  * •

The residue labels: whether a residue nearby a binding site is a part of the
binding site or not can be determined by whether the same is true for
corresponding residues of the homologous binding sites.

These hypotheses except the second one have been the bases of the template-
based BSP methods introduced in 2.1.3.

However, the second hypothesis is, to our knowledge, has not been employed in
the previous BSP methods. The result of ablating the homology augmentation in
Table 1 provides a partial evidence that this hypothesis is valid at least to
some degree. Otherwise if the negative labels did not provide valuable
learning signals, the homology augmentation would have had only adversary
effects on the BSD module’s performance.

## References

  * [AGC+21] Rishal Aggarwal, Akash Gupta, Vineeth Chelur, CV Jawahar, and U Deva Priyakumar.  Deeppocket: ligand binding site detection and segmentation using 3d convolutional neural networks.  Journal of Chemical Information and Modeling, 2021.
  * [BCM+10] Natasja Brooijmans, Yu-Wei Chang, Dominick Mobilio, Rajiah A Denny, and Christine Humblet.  An enriched structural kinase database to enable kinome-wide structure-based analyses and drug discovery.  Protein Science, 19(4):763–774, 2010.
  * [BGL+21] Brian J Bender, Stefan Gahbauer, Andreas Luttens, Jiankun Lyu, Chase M Webb, Reed M Stein, Elissa A Fink, Trent E Balius, Jens Carlsson, John J Irwin, et al.  A practical guide to large-scale docking.  Nature protocols, 16(10):4799–4832, 2021.
  * [BMR08] 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\.
  * [BNL03] T Andrew Binkowski, Shapor Naghibzadeh, and Jie Liang.  Castp: computed atlas of surface topography of proteins.  Nucleic acids research, 31(13):3352–3355, 2003.
  * [Car03] Daniel C Carter.  Crystallographic survey of albumin drug interaction and preliminary applications in cancer chemotherapy.  Burger’s Medicinal Chemistry and Drug Discovery, pages 437–468, 2003.
  * [CLT+09] John A Capra, Roman A Laskowski, Janet M Thornton, Mona Singh, and Thomas A Funkhouser.  Predicting protein ligand binding sites by combining evolutionary sequence conservation and 3d structure.  PLoS computational biology, 5(12):e1000585, 2009.
  * [CMGK11] Ke Chen, Marcin J Mizianty, Jianzhao Gao, and Lukasz Kurgan.  A critical comparative assessment of predictions of protein-binding sites for biologically relevant organic compounds.  Structure, 19(5):613–621, 2011.
  * [DBRK15] Jérémy Desaphy, Guillaume Bret, Didier Rognan, and Esther Kellenberger.  sc-pdb: a 3d-database of ligandable binding sites—10 years on.  Nucleic acids research, 43(D1):D399–D404, 2015.
  * [FWFW20] Fabian Fuchs, Daniel Worrall, Volker Fischer, and Max Welling.  Se (3)-transformers: 3d roto-translation equivariant attention networks.  Advances in Neural Information Processing Systems, 33:1970–1981, 2020.
  * [Hal09] Thomas A Halgren.  Identifying and characterizing binding sites and assessing druggability.  Journal of chemical information and modeling, 49(2):377–389, 2009\.
  * [HBF+01] Sepp Hochreiter, Yoshua Bengio, Paolo Frasconi, Jürgen Schmidhuber, et al.  Gradient flow in recurrent nets: the difficulty of learning long-term dependencies, 2001.
  * [HGS09] Marylens Hernandez, Dario Ghersi, and Roberto Sanchez.  Sitehound-web: a server for ligand binding site identification in protein structures.  Nucleic acids research, 37(suppl_2):W413–W416, 2009.
  * [HZRS16] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [JDMR+17] José Jiménez, Stefan Doerr, Gerard Martínez-Rosell, Alexander S Rose, and Gianni De Fabritiis.  Deepsite: protein-binding site predictor using 3d-convolutional neural networks.  Bioinformatics, 33(19):3036–3042, 2017.
  * [JEP+21] John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Žídek, Anna Potapenko, et al.  Highly accurate protein structure prediction with alphafold.  Nature, 596(7873):583–589, 2021.
  * [JK19] Justin M Johnson and Taghi M Khoshgoftaar.  Survey on deep learning with class imbalance.  Journal of Big Data, 6(1):1–54, 2019.
  * [KH18] Radoslav Krivák and David Hoksza.  P2rank: machine learning based tool for rapid and accurate prediction of ligand binding sites from protein structure.  Journal of cheminformatics, 10(1):1–12, 2018.
  * [LGST09] Vincent Le Guilloux, Peter Schmidtke, and Pierre Tuffery.  Fpocket: an open source platform for ligand pocket detection.  BMC bioinformatics, 10(1):1–11, 2009.
  * [LJ05] Alasdair TR Laurie and Richard M Jackson.  Q-sitefinder: an energy-based method for the prediction of protein–ligand binding sites.  Bioinformatics, 21(9):1908–1916, 2005.
  * [LYH+19] Zhirui Liao, Ronghui You, Xiaodi Huang, Xiaojun Yao, Tao Huang, and Shanfeng Zhu.  Deepdock: enhancing ligand-protein interaction prediction by a combination of ligand and structure information.  In 2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pages 311–317. IEEE, 2019.
  * [MAD21] Stelios K Mylonas, Apostolos Axenopoulos, and Petros Daras.  Deepsurf: a surface-based deep learning approach for the prediction of ligand binding sites on proteins.  Bioinformatics, 37(12):1681–1690, 2021.
  * [MLAdRCW21] 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.
  * [NHZ+12] Chi-Ho Ngan, David R Hall, Brandon Zerbe, Laurie E Grove, Dima Kozakov, and Sandor Vajda.  Ftsite: high accuracy detection of ligand binding sites on unbound protein structures.  Bioinformatics, 28(2):286–287, 2012.
  * [NKC+16] Tony Ngo, Irina Kufareva, James LJ Coleman, Robert M Graham, Ruben Abagyan, and Nicola J Smith.  Identifying ligands at orphan gpcrs: current status using structure-based approaches.  British journal of pharmacology, 173(20):2934–2951, 2016.
  * [PJ95] Theodore Peters Jr.  All about albumin: biochemistry, genetics, and medical applications.  Academic press, 1995.
  * [PKP+21] Jimin Park, Mi-Sun Kim, Taeseong Park, Young Hwan Kim, and Dong Hae Shin.  Crystal structure of pharmaceutical-grade human serum albumin.  International Journal of Biological Macromolecules, 166:221–228, 2021.
  * [RFB15] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [RZ12] Ambrish Roy and Yang Zhang.  Recognizing protein-ligand binding sites by global structural alignment and local geometry refinement.  Structure, 20(6):987–997, 2012.
  * [SB09] Jeffrey Skolnick and Michal Brylinski.  Findsite: a combined evolution/structure-based approach to protein function prediction.  Briefings in bioinformatics, 10(4):378–391, 2009.
  * [SDZS20] Marta M Stepniewska-Dziubinska, Piotr Zielenkiewicz, and Pawel Siedlecki.  Improving detection of protein-ligand binding sites with 3d segmentation.  Scientific reports, 10(1):1–9, 2020.
  * [SGP+22] 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.
  * [SK19] Connor Shorten and Taghi M Khoshgoftaar.  A survey on image data augmentation for deep learning.  Journal of big data, 6(1):1–48, 2019.
  * [SSS+19] Jiyong Su, Chenyang Song, Yunlong Si, Linlin Cui, Tong Yang, Yuying Li, Hao Wang, Guihua Tai, and Yifa Zhou.  Identification of key amino acid residues determining ligand binding specificity, homodimerization and cellular distribution of human galectin-10.  Glycobiology, 29(1):85–93, 2019.
  * [TAW+21] Kathryn Tunyasuvunakool, Jonas Adler, Zachary Wu, Tim Green, Michal Zielinski, Augustin Žídek, Alex Bridgland, Andrew Cowie, Clemens Meyer, Agata Laydon, et al.  Highly accurate protein structure prediction for the human proteome.  Nature, 596(7873):590–596, 2021.
  * [VSP+17] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Advances in neural information processing systems, pages 5998–6008, 2017.
  * [WHR+13] Zhong-min Wang, Joseph X Ho, John R Ruble, John Rose, Florian Rüker, Melanie Ellenburg, Robert Murphy, James Click, Elizabeth Soistman, Leslie Wilkerson, et al.  Structural studies of several clinically important oncology drugs in complex with human serum albumin.  Biochimica et Biophysica Acta (BBA)-General Subjects, 1830(12):5356–5374, 2013.
  * [WKS10] Mark N Wass, Lawrence A Kelley, and Michael JE Sternberg.  3dligandsite: predicting ligand-binding sites using similar structures.  Nucleic acids research, 38(suppl_2):W469–W473, 2010.
  * [YRZ13] Jianyi Yang, Ambrish Roy, and Yang Zhang.  Protein–ligand binding site recognition using complementary binding-specific substructure comparison and sequence profile alignment.  Bioinformatics, 29(20):2588–2595, 2013.
  * [ZCZ20] Jingtian Zhao, Yang Cao, and Le Zhang.  Exploring the computational methods for protein-ligand binding site prediction.  Computational and structural biotechnology journal, 18:417–426, 2020.
  * [ZP11] Hongbo Zhu and M Teresa Pisabarro.  Mspocket: an orientation-independent algorithm for the detection of ligand binding pockets.  Bioinformatics, 27(3):351–358, 2011.