Boosting Convolutional Neural Networks' Protein Binding Site Prediction
  Capacity Using SE(3)-invariant transformers, Transfer Learning and
  Homology-based Augmentation
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$\unicode{x2013}$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$\unicode{x2013}$the model architecture, transfer learning and the
homology augmentation$\unicode{x2013}$would serve as useful components in
future works.
Daeseok Lee, Jeunghyun Byun, Bonggun Shin
Boosting convolutional neural networks' protein binding site prediction capacity using SE(3)-invariant transformers, transfer learning and homology-based augmentation

###### 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_ ([1], [2], [3]) and is required for _in-silico_ methods such as docking ([1], [2]). 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 ([1]). Also, preparing a suitable binding site is the first step in any virtual structure-based drug discovery pipeline ([1]).

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 ([2], [3], [4]). 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 ([2], [1]), and docking results can be prioritized according to whether the docked molecule has favored interactions with the key residues ([1], [2]).

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 ([1]), 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 ([1]).

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 ([12]). 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 ([10], [11]). 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 ([13]). Besides using an SE(3)-invariant attention mechanism ([12]) 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 [1]. 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 ([14], [15], [16], [17]), 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 ([15], [17]), 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.

### 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 ([11], [12], [13], [14]). 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. ([15], [16], [17])

One example is Frocket ([17]), 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 Frocket typically produces an excessive number of binding pockets, it has relatively good recall (96.4% on scPDB v.2017, according to [1]). Therefore, there are ML-based algorithms ([1], [1]) 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 ([18], [19], [17], [18]). 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 [17] 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. Putative binding pockets are identified in the query protein by relying on an external software ([10]).
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. The ligands in the putative templates are projected to the binding pocket.
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. The query protein sequence is aligned with the template sequences based on their position-specific scoring matrices (PSSM) profiles and secondary structure information.
2. Templates with the highest alignment _quality scores_ are chosen as the putative templates.
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.

### deep learning based BSP Methods

#### 2.2.1. DeepSite

DeepSite([1]) 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([1] has the same overall procedure with DeepSite([1]), 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 ([18]).
#### 2.2.3. Kalasanty

Kalasanty ([16]) tries to solve BSP by viewing it as a 3D _image segmentation_ problem. Therefore, it uses a 3D equivalent of the U-net model ([12]), 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 ([1]) 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.

### Predicting ligand-specific binding sites

Recently, deep learning models that predict protein-ligand complex structures, given a protein-ligand pair, have been developed ([1], [10], [15]). 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 ([15]) 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.

### 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\{1,\cdots,size(P)\}\) 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.

### 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 ([13]) and Deeppocket ([1])). 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 _drugability_ 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

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.

### The candidate generation module

To generate the binding site candidate centers, we use an external software Fpocket ([1]). 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 \(\tilde{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 [1]). 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.

### The BSD module

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

In doing so, it featurizes the surroundings of \(\tilde{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 \(\tilde{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.

### 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.

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.
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.

### 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.

### 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 \(\{x_{i}\}\,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 \(\{x_{i}\}\,i=1^{n}\) into another sequence of hidden vectors
\(\{x^{\prime}_{i}\}\,i=1^{n}\) based on the protein structure. In our method, we represent the protein structure as _local frames_ ([1]) \(\left\{T_{i}\right\}_{i=1}^{n}=\left\{\left(R_{i},t_{i}\right)\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:

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

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 ([1]).

For the attention mechanism, we used a modified version of _Invariant Point Attention_ from Alphafold ([1]). The essential ideas of its computation, compared to the standard attention mechanism ([13]), 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.

### 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\{Tu_{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 ([13]).

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:

\[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}). \tag{4.2}\]

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.

### 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 ([15]). 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 ([1], [1], [1]) 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

### Datasets

We used scPDB v.2017 ([1]) as a main dataset for the training and the validation. In addition, we used three other datasets for the tests: COACH420 ([11]), HOLO4K ([1]) and CHEN ([12]). To be more specific, we used the training subset of scPDB dataset provided by [13] 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.

### Baseline Methods

We compared our method to the previous state-of-the-art Deep Learning methods, which are based on CNN: Deeppocket ([1]), Kalasanty ([1]) and DeepSurf ([1]). 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.

### 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 ([10]), 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.

### 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_{1}^{(i)},\cdots,l_{n^{(i)}}^{(i)}\right\}\) is the set of ground-truth ligands bound to the \(i\)-th protein.
* \(\left\{(c_{1}^{(i)},BR_{1}^{(i)})\cdots,(c_{n^{(i)}}^{(i)},BR_{n^{(i)}}^{(i)})\right\}\) is the set of predictions of the method to evaluate.
* \(\left\{TBR_{1}^{(i)},\cdots,TBR_{n^{(i)}}^{(i)}\right\}\) is the set of _true binding residue indices_, where \(TBR_{j}^{(i)}\) is defined to be the set of residues in the \(i\)-th protein that is within \(4\AA\) from \(l_{j}^{(i)}\)

The success rate metric measures the correspondence between the predicted binding site centers \(\left\{c_{1}^{(i)},\cdots,c_{n^{(i)}}^{(i)}\right\}\) and the positions of the ground-truth ligands \(\left\{l_{1}^{(i)},\cdots,l_{n^{(i)}}^{(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_{j}^{(i)}\) is a _correct detection_ of \(l_{k}^{(i)}\) when \(c_{j}^{(i)}\) is within \(4\AA\) (a threshold commonly used in the literature e.g. [1] and [1]) from any ligand of \(l_{k}^{(i)}\). 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_{j}^{(i)}\text{ detects one of }l_{1}^{(i)},\cdots l_{n^{(i)}}^{(i)}\right\}}{n^{(i)}}\]

and \(R^{(i)}\) is the _recall_ defined as follows:

\[R^{(i)}=\frac{\#\left\{1\leq k\leq n^{(i)}:l_{k}^{(i)}\text{ is detected by one of }c_{1}^{(i)},\cdots,c_{n^{(i)}}^{(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_{j}^{(i)}\) with the true binding residues \(TBR_{\phi^{(i)}(j)}^{(i)}\) of the ligand \(l_{\phi^{(i)}(j)}^{(i)}\) closest to the predicted binding site center \(c_{j}^{(i)}\). Here, the index \(\phi^{(i)}(j)\) of the closest ligand is defined as

\[\phi^{(i)}(j)=\operatorname*{arg\,min}_{k=1}^{n^{(i)}}DCA(c_{j}^{(i)},l_{k}^{( i)})\]
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}^{(i)}_{j}\) 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},t^{(i)}_{\phi^{(i)}(j)})<4 \hat{A}\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.

### 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.

### 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.

### 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".

### Experiment Results

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
\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline  & s\(\epsilon\)PDB(field-out) & COACH420 & HOLO4K & CHE=holo & CHE=apo \\ \hline 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\) \\ Kalasuardy & \(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\) \\ DeepPosed & \(67.9\pm 0.4\) & \(55.7\pm 0.7\) & \(72.9\pm 0.2\) & \(\mathbf{42.4\pm 0.3}\) & \(34.5\pm 1.2\) \\ \hline Ours & \(\mathbf{70.1\pm 0.4}\) & \(59.1\pm 0.3\) & \(\mathbf{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\) & \(\mathbf{36.7\pm 1.0}\) \\ Ours(no alignment) & \(69.0\pm 0.6\) & \(\mathbf{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\) & \(58.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\) \\ \hline \end{tabular}
\end{table}
Table 1. (BSD metric) F1 Success rate for detection. The mean and standard deviation are calculated based on the metric values of 5 different loss-validation folds.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline  & s\(\epsilon\)PDB(field-out) & COACH420 & HOLO4K & CHE=holo & CHE=apo \\ \hline 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\) \\ Kalasuardy & \(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\) \\ DeepPosed & \(0.095\pm 0.002\) & \(0.506\pm 0.005\) & \(0.371\pm 0.004\) & \(0.339\pm 0.008\) & \(0.382\pm 0.000\) \\ Ours & \(\mathbf{0.643\pm 0.004}\) & \(\mathbf{0.685\pm 0.007}\) & \(\mathbf{0.415\pm 0.002}\) & \(\mathbf{0.495\pm 0.010}\) & \(\mathbf{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\) \\ \hline \end{tabular}
\end{table}
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.

Figure 4. The effect of the augmentation methods on the training

Figure 3. The effect of transfer learning on the BSD training
(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.
### The structure and Binding Sites 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 ([4], [10]). The authors of [10] 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 [12] provided more detailed analysis on the binding sites at the IB subdomain. They did so by inspecting the crystal structures 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.

### 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 ([10]. 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.

### Binding Site Detection

#### 6.3.1. experiment procedure

First, we obtained our model and Deeppocket's BSD module's predictions on 15 binding sites identified by [10]. 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 [10]. 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 [10].

#### 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.

### Binding Residue Identification

Figure 5. The subdomains of HSA
Figure 6. A comparison between our BSD module’s predictions and the _frequencies_ from [1]

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 [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 [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

### 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 [1]. More precisely, when \(\mathbf{x}_{1}\), \(\mathbf{x}_{2}\) and \(\mathbf{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=(\mathbf{e}_{1}\ \mathbf{e}_{2}\ \mathbf{e}_{3})\) that we refer to as "orientation" is obtained as follows:

\[\mathbf{v}_{1} =\mathbf{x}_{3}-\mathbf{x}_{2}\] \[\mathbf{v}_{2} =\mathbf{x}_{1}-\mathbf{x}_{2}\] \[\mathbf{e}_{1} =\mathbf{v}_{1}/\left\|\mathbf{v}_{1}\right\|\] \[\mathbf{u}_{2} =\mathbf{v}_{2}-(\mathbf{e}_{1}\cdot\mathbf{v}_{2})\mathbf{e}_{1}\] \[\mathbf{e}_{2} =\mathbf{u}_{2}/\left\|\mathbf{u}_{2}\right\|\] \[\mathbf{e}_{3} =\mathbf{e}_{1}\times\mathbf{e}_{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. Collect coordinates \(\mathbf{x}_{i}\in\mathbb{R}^{3}\) (\(1\leq i\leq n\)) of the heavy (non-hydrogen) atoms in the protein.
2. Collect feature vectors \(f_{i}\in\mathbb{R}^{d_{feature}}\) (\(1\leq i\leq n\)) of the heavy atoms.
3. Given a choice of _grid axes_\((\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{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})\mathbf{e}_{i}\quad((a_{1},a_{2},a_{3}) \in\{0,\cdots,L-1\}^{3})\], where \(t\) is the grid center (the alpha carbon coordinate), \(r\) is the _grid resolution_ and \(L\) is the _grid size_.
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 [1], [1] and [25], 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. ([25], [1].)

### 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).
### 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 [1]. The model is adapted from the bottleneck ResNet model introduced in [15] for the image classification. The bottleneck architecture reduces the number of parameters thus enables the use of a deeper network. [1] showed that the 3D BottleNeck ResNet model, while being light, performed competitively well compared to its non-bottleneck counterpart.

### 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}^{\neg}u\) plays the reverse role.

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

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. The geometric query and key vectors \(\mathbf{q}_{i}^{hp}\) and \(\mathbf{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. 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}\mathbf{q} _{i}^{hp}-T_{j}\mathbf{k}_{j}^{hp}\right\|\] by applying a softmax operation. More precisely, the attention weight becomes \[w_{ij}^{h}=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. 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. The geometric value vectors \(\mathbf{v}_{j}^{h}\) are computed by a linear map from \(x_{j}\), and aggregated as (7.2) \[\mathbf{o}_{i}^{hp}=T_{i}^{-1}(\sum_{j}w_{ij}^{h}T_{j}\mathbf{v}_{j}^{hp})\]
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},\mathbf{o}_{i}^{hp},\left\| \mathbf{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.

### 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 \(\{x_{i}\}_{i=1}^{n}\) and \(\{T_{i}\}_{i=1}^{n}\), satisfy the same equality. For any \(T\in SE(3)\), the (7.1) satisfies the equality because \(\left\|TT_{i}\mathbf{q}_{i}^{hp}-TT_{j}\mathbf{k}_{j}^{hp}\right\|=\left\|T_{ i}\mathbf{q}_{i}^{hp}-T_{j}\mathbf{k}_{j}^{hp}\right\|=\left\|\text{(since $T$ preserves the norm)}\), and (7.2) satisfies the equality because \(\left(TT_{i}\right)^{-1}(\sum_{j}w_{ij}^{h}(TT_{j})\mathbf{v}_{j}^{hp})=T_{i}^ {-1}T^{-1}T(\sum_{j}w_{ij}^{h}T_{j}\mathbf{v}_{j}^{hp})=T_{i}^{-1}(\sum_{j}w_{ ij}^{h}T_{j}\mathbf{v}_{j}^{hp})\). Note that this derivation is almost identical to that presented in [JEP\({}^{+}\)21](Supplementary Information, page 28).
### The Homology-based Augmentation

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

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 protein. The bright green region must comprise at least 50% of the entire green region in order for the binding site to count.

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.
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. a query chain \(x\) in \(\mathcal{S}\) 2. a ligand \(l\) associated to \(x\) found in step 1 of the procedure and 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. 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. 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.
