Introduction

The study of molecular interaction is the basis of drug or toxicity research. In vitro and in vivo experiments can be preceded by in silico filtering. These numerical simulations need huge computing resources, making them a relevant high-performance computing (HPC) use case. Moreover, the rise of artificial intelligence (AI) and specifically machine learning (ML) over the last decade has permitted better accuracy and reduced computational time. In that context, molecular docking is one of the most prominent modern techniques and is used to predict the likelihood that a small molecule, namely the ligand, will bind to a macromolecule, the receptor. In the case of that study, the former is a small chemical compound while the latter is a protein. Unlike ligand-based molecular docking which uses known ligand-receptor interactions to determine the score of the complex, structure-based molecular docking uses three-dimensional (3D) structures of both molecules to predict the complex conformation and its strength (the score). In the remainder of this article, all molecular dockings are structure-based. Molecular docking is a two-step process [1]. The first step, also called sampling, explores the ligand conformational space and in some cases the protein conformational space as well. Numerous exploration methods exist such as shape matching or systematic (iterative or incremental search), and both are based on simulation processes. However, the most popular methods are stochastic ones and rely on genetic algorithms. The second step consists in scoring each protein-ligand complex pose produced by the aforementioned step. That second step lists a great number of classic AI-free methods, which can be categorized into three classes. The first category is physics-based methods which use a weighted sum of energy terms to produce a score, the second is empirical methods where functions use a weighted sum of computationally-easy physicochemical terms, and the last is knowledge-based and uses experimental knowledge about docked protein-ligand complexes and their associated score. Finally, a mix of scoring functions from those three classes can act as an election quorum to produce a consensus score. In the case of stochastic exploration methods, that score is also used as a support guide for sampling.

Numerical simulations can also be used to reduce the exploration surface. Blind docking is used when no experimental knowledge or intuition-driven search is available, then the whole protein surface must be explored, and this tends to result in poorer models and less accurate scoring. Instead, targeted docking focuses only on promising sites, which is why binding site detection methods now exist to predict the putative binding site of a protein using its 3D structure.

Most recently, AI and ML have been investigated to speed up molecular docking. Indeed, the use of AI and ML has grown considerably in all scientific fields, and, in the case of molecular docking, ML started in the 2010s with random forest (RF) based methods such as RF-score [2]. The first use of neural network methods also dates to that period, but their use soared with convolutional neural networks (CNNs), introduced in 2015 with AtomNet [3]. The simplest way to represent data in these methods is the atom-based discretization of both molecules on a 3D grid of predefined cell size. Physicochemical properties are then added as features on each atom.

Lately, new variants of deep learning (DL) methods have emerged, using graph convolutional networks (GCNs) that generalize the well-known grid-based convolution to unstructured or graph-structured data. The graph structure allows for representing a molecule more reliably and accurately, allowing the proposal of more powerful scoring functions.

Hence, this paper proposes the heterogeneous graph score (HGScore) method, a protein-ligand complex scoring function based on a heterogeneous GCN (HGCN). The method is available on GitHub: github.com/KevinCrp/HGScore. After introducing existing AI methods for scoring, a dataset often used in the literature and the data representation are presented. Then, the proposed scoring method is described and this article is concluded by presenting a comparative study against existing methods.

Materials and methods

Dataset

To compare the performance of the proposed method and other works, the Protein Data Bank-bind (PDBbind) 2020 dataset [12] is used. It contains 19,443 ligand-protein complexes, each input of which is composed of a Protein Data Bank (PDB) file representing the protein, a PDB file for the protein binding pocket, and a MOL2 file for the ligand. In PDBbind, a binding pocket contains all residues whose at least one heavy atom is closer than 10 Å to one of the ligand atoms. Each ligand-protein complex has its associated binding affinity, express by –log(K_i) or –log(K_d) as a function of the complexes. Hereafter, only the protein binding pocket is used to describe the protein part. PDBbind authors have proposed splitting their database into several sets, as shown in Figure 1. The first is the general set containing all database inputs, while the second is the refined set created by selecting only the best-quality complexes from that first set. Both sets are updated each year. Finally, the third set consists in a core set that contains very few inputs but with high-quality [13].

Display full size

Figure 1. 

The PDBbind dataset and its components. The general set with 19,443 protein-ligand complexes, the refined set (included in the general) with 5,316 complexes, and both comparative assessment of scoring functions (CASF) benchmarks with 195 and 285 complexes respectively for the 2013 and the 2016 versions

In this work, the versions 2013 and 2016 of the core set are used as test sets (195 and 285 items, respectively), and 1,000 inputs from the refined set are randomly extracted to build the validation set and the remaining input as the training set (18,070 items). It is ensured that there are no duplicates between sets, except for the core sets which naturally have overlapping entries (elements from 2016 already present in the 2013 dataset). Also, in the pre-processing step, the protein files are cleaned by keeping only protein atoms, so that water and other hetero atoms are removed.

Data representation

A graph $\mathcal{G}\left(\mathcal{V},\mathrm{ }\mathcal{E}\right)$ is a set of V nodes (also called vertices) $\mathcal{V}$ linked by a set of E edges $\mathcal{E}$. Each node $v_i\mathrm{ }\in \mathrm{ }\mathcal{V}$ can have a set of features named $\overrightarrow{v_{\iota }}\in {𝕉}^n$. Idem for the edges with their respective $\overrightarrow{e_{\iota }}\in {𝕉}^m$.

A heterogeneous graph, as presented in Figure 2, is composed of several graphs but in that case, only a pair of graphs, [${\mathcal{G}}_1\left({\mathcal{V}}_1,\mathrm{ }{\mathcal{E}}_1\right)$ and ${\mathcal{G}}_2\left({\mathcal{V}}_2,\mathrm{ }{\mathcal{E}}_2\right)$] is used, each can have a different set of features: ${\mathcal{V}}_1\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{N_1 \times { n}_1}$, ${\mathcal{V}}_2\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{N_2 \times { n}_2}$, ${\mathcal{E}}_1\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{E_1 \times { m}_1}$ and ${\mathcal{E}}_2\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{E_2 \times m_2}$. In addition, two new sets of edges (${\mathcal{E}}_{1\to 2}$ and ${\mathcal{E}}_{2\to 1}$) allow connecting nodes from one graph to the second: ${\mathcal{E}}_{1\to 2}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{E_{12} \times l_{12}}$ and ${\mathcal{E}}_{2\to 1}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{E }_{21}\times l_{21}}$.

Display full size

Figure 2. 

Example of a heterogeneous graph. Blue nodes belong to ${\mathcal{G}}_1$ and green nodes belong to ${\mathcal{G}}_2$ . Blue edges belong to ${\mathcal{E}}_1$, green edges belong to ${\mathcal{E}}_2$, purple arrows correspond to ${\mathcal{E}}_{1\to 2}\mathrm{ }$and red arrows correspond to ${\mathcal{E}}_{2\to 1}\mathrm{ }$

Each molecule (protein and ligand) is represented with a graph of heavy atoms, thus, each graph node represents an atom, and the graph edges represent molecular bonds. A total of 23 physicochemical features are computed using Open Drug Discovery Toolkit (ODDT) [14] and are presented in Table 1. The 23 features are used as the second dimension for both the protein and ligand nodes. Then each ligand-protein atom pair closer than 4 Å is connected, also completed by a set of features enumerated in Table 2. Then after featurization is complete, each complex is represented by a protein graph ${\mathcal{G}}_{\mathrm{P}}\left({\mathcal{V}}_{\mathrm{P}},\mathrm{ }{ \mathcal{E}}_{\mathrm{P}}\right)$, ${\mathcal{V}}_{\mathrm{P}}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{\mathrm{N}}_{\mathrm{P}}\mathrm{ }\times \mathrm{ }23}$, ${\mathcal{E}}_{\mathrm{P}}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{\mathrm{E}}_{\mathrm{P}}\mathrm{ }\times \mathrm{ }6}$, a ligand graph ${\mathcal{G}}_{\mathrm{L}}\left({\mathcal{V}}_{\mathrm{L}},\mathrm{ }{ \mathcal{E}}_{\mathrm{L}}\right)$,$V_{\mathrm{L}}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{\mathrm{N}}_{\mathrm{L}}\mathrm{ }\times \mathrm{ }23}$ ,${\mathcal{E}}_{\mathrm{L}}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{\mathrm{E}}_{\mathrm{L}}\mathrm{ }\times \mathrm{ }6}$ and two inter-molecular edge sets ${\mathcal{E}}_{\mathrm{P}\to \mathrm{L}}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{\mathrm{E}}_{\mathrm{P}\mathrm{L}}\mathrm{ }\times \mathrm{ }6}$ and ${\mathcal{E}}_{\mathrm{L}\to \mathrm{P}}\mathrm{ }\in \mathrm{ }{\mathbb{R}}^{{\mathrm{E}}_{\mathrm{L}\mathrm{P}}\mathrm{ }\times \mathrm{ }6}$. A protein-ligand complex and its corresponding heterogeneous graph are represented in Figure 3. Heterogeneous graphs are the way to get the most accurate representation, allowing to represent atoms and bonds while differentiating them, in that case, according to their belonging to the protein or ligand.

Table 1.

Nodes features

Feature	Size	Description
Atom type	8	One-hot encoded (B, C, N, O, F, P, S, and others)
Hybridization	7	One-hot encoded (other, sp, sp2, sp3, sq planar, trigonal bipyramidal, octahedral)
Partial charge	1	Float (Gasteiger charge model)
Is hydrophobic	1	Boolean
Is in an aromatic cycle	1	Boolean
Is H-bond acceptor	1	Boolean
Is H-bond donor	1	Boolean
Is H-bond donor hydrogen	1	Boolean
Is negatively charged/chargeable	1	Boolean
Is positively charged/chargeable	1	Boolean

Display full size

Table 2.

Intra- and Inter-molecular edge features

Feature	Size	Intra	Inter	Description
Length/distance	1	X	X	Float
Belongs to an aromatic cycle	1	X	-	Boolean
Belongs to a ring	1	X	-	Boolean
Type	3	X	-	One-hot (simple, double, or triple)
Is a hydrogen bond	1	-	X	Boolean
Is a salt bridge	1	-	X	Boolean
Is a hydrophobic contact	1	-	X	Boolean
Has π-stacking	1	-	X	Boolean
Has π-cation	1	-	X	Boolean

Display full size

π-stacking, and π-cation are only available for inter-molecular edges whose terminal atom from the protein side belongs to a ring. X: applicable; -: not applicable

Display full size

Figure 3. 

A protein-ligand complex and its associated heterogeneous graph. On the left, the 29G11 protein (in red) complexed with the ligand HEP in green (PDB: 1a0q). On the right, the heterogeneous graph represents the same complex, red nodes represent the protein binding site atoms, and the green nodes are the ligand atoms, both without hydrogen. Red and green edges represent intramolecular bonds, red edges belong to the protein while green edges belong to the ligand. Finally, the cyan edges are the inter-molecular links representing the protein-ligand interaction

Method

Graph convolution

Graph neural networks (GNNs) are a flexible type of neural objects that exist in various flavors, as classified by Wu et al. [15]. While recurrent GCNs (RGCNs) were the first attempt at modern graph neural modeling, nowadays GCNs are the most popular kind. GCNs are created by stacking multiple graph convolutional layers and updating the state of the graph at each step.

The implementation is based on PyTorch Geometric [16], a framework relying on the concept of message-passing neural network (MPNN, Equation 1) that includes libraries for sparse computation. Hence, the principle of spatial convolution as presented in Figure 4 is used.

(1)$$h_i^{(l+1)}=U_l\left(h_i^{(l)},{\sum }_{j\in 𝒩(i)}{M}_l(h_i^{(l)},h_j^{(l)},e_{j\to i})\right)$$

Where $h_i^{(l)}$ is the hidden features vector of node i at layer l, $h_l^{(0)}=\overrightarrow{v_{\iota }}$, U_l and M_l are functions with learnable parameters and N(i) is the neighborhood of node i.

Display full size

Figure 4. 

The spatial convolution process. In spatial convolution, the current node, here identified by a star, aggregates information coming from its direct neighbors here represented in green, purple, and blue. The current node feature vector is then updated as shown in the graph on the right

Layers

The proposed method is based on Attentive FP [17], which was adapted to work with molecular complexes. While the neural model was originally designed to produce a molecular fingerprint allowing the prediction of chemical properties such as toxicity or solubility, it was adapted to be applied to a molecular complex (ligand-protein) to produce a binding affinity score. That model is composed of graph attentional convolution (GATConv) and gated recurrent unit cell (GRUCell) layers, explained below.

The GATConv layer proposed by Veličković et al. [18] allows using an attention mechanism on GCN. In this work, the second version of the GATConv (GATv2Conv) layer proposed by Brody et al. [19] that fixes the static attention problem of GATConv is used. Thus, the network can learn which neighbor nodes are more relevant thanks to implicit weights. Indeed, in addition to the matrix of weights W, the GATConv operators learn an attention coefficient a_ij representing the magnitude of the edge j → i or in other words the importance of node features of j for explaining node i. All neighbors j of node i have a coefficient representing their importance for i. Equations 2 and 3 are used to compute the attention coefficient, and Equation 4 updates the current node features according to this coefficient.

(2)$$e(h_i,h_j)=a^{⊺}\text{LeakyReLU}(W{h}_i\parallel W{h}_j)$$

Where the function e allows scoring the edge from j to i nodes, a and W are learned, the ‖ is the vector concatenation operator, ⊺ the transposition operator, and LeakyReLU the leaky rectified linear unit.

(3)$${\alpha }_{ij}={\text{softmax}}_j(e(h_i,h_j))=\frac{exp(e(h_i,h_j))}{{\sum }_{j^{\prime }\in {𝒩}_i}\text{exp}(e(h_i,h_{j^{\prime }}))}$$

As explained by Veličković et al. [18], all attention coefficients are normalized by a soft-max function over all node neighbors.

(4)$$h_i^{\prime }=\sigma ({\sum }_{j\in 𝒩(i)}{\alpha }_{ij}W{h}_j)$$

Where σ is a nonlinearity function.

The GRUCell layer proposed by Cho et al. [20] is a recurrent neural network layer that keeps the memory of items over long sequences of data. A GRUCell layer is described by four equations to compute the layer update. The first two equations are the update gate (Equation 5) and the reset gate (Equation 6), respectively allowing relevant information to flow from the past (h_t–1) or to be forgotten. The current memory content follows Equation 7 and stores the relevant information from the past. Finally, the layer is updated with Equation 8.

(5)$$z_t=\sigma \left(W^{(z)}{x}_t+U^{(z)}{h}_{t-1}\right)$$

(6)$$r_t=\sigma \left(W^{(r)}{x}_t+U^{(r)}{h}_{t-1}\right)$$

(7)$$h_t^{\prime }=\text{tanh}(W{x}_t+r_t\odot U{h}_{t-1})$$

(8)$$h_t=z_t\odot h_{t-1}+(1-z_t)\odot h_t^{\prime }$$

Where ⊙ is the Hadamard product.

HGScore

As mentioned above, the heterogeneous graph contains two sub-graphs (protein and ligand) in which nodes are connected. Thus, 4 types of edges: 𝒢_p ↔ 𝒢_p, 𝒢_L ↔ 𝒢_L, 𝒢_p ↔ 𝒢_L, and 𝒢_L ↔ 𝒢_p are used, resulting in 6 independent GCNs (one for each edge type in atomic part and one for each molecule in molecular part), and ending with a multi-layer perceptron (MLP) to compute the final score as presented in Figures 5–7. Different ways to construct the network including one that does not use the molecular embedding part and feeds the atomic embedding pooling to the MLP have been tested. The proposed method, used in the remainder of this paper, is the strategy that showed the best results.

Display full size

Figure 5. 

The HGScore atomic embedding part. The atomic embedding part extends Attentive FP atomic embedding to work on the heterogeneous graph to produce an atomic embedding of both molecules. The aggregation (Aggr) module allows a combination of previous layer outputs using mean, add, or max Aggr. ReLu: rectified linear unit

Display full size

Figure 6. 

The HGScore molecular embedding part. The two atomic part embeddings are pooled by a sum to produce a molecular embedding which is used in the classical Attentive FP molecular embedding part to produce protein and ligand fingerprints

Display full size

Figure 7. 

The HGScore prediction part. The two molecular part fingerprints are concatenated into an interaction fingerprint allowing an MLP to predict an affinity score

The atomic-GCN embedding part of the model, presented in Figure 5, encodes atomic features in a latent space. For both molecules, nodes are updated by two GCNs, one for each edge type. This resulted in four GCNs, all Attentive FP-inspired and using GATv2Conv and GRUCell. The main idea is to dissociate the embedding for each edge type so that each model will mix an intertwined pair of input signals enriched by their specific context (protein-protein, protein-ligand, ligand-protein, ligand-ligand). The resulting signal propagates information from context to context several times over by repeating the same operation.

Once the atomic-GCN embedding is completed, the same algorithm (Figure 6) proposed by Attentive FP (the molecular embedding) is applied to each molecule fingerprint, which consists of a sum global pooling, i.e., all node vectors are summed to produce a unique molecular representation. This vector is then processed by a molecular-GCN composed of GATv2Conv and GRUCell, also built from Attentive FP.

Finally, the concatenation of the two out-coming molecular fingerprints is used to produce the interaction fingerprint, which is fed to a classic feed-forward MLP (Figure 7), resulting in an affinity score of the complex.

Results

The method’s implementation is based on PyTorch, Lightning, and PyTorch Geometric [16] to build the graph data, the datasets, and the model, and perform training and evaluation. Optuna was used to perform hyper-optimization to find the best architecture and perform the training with the same hyperparameter set several times over because the model is sensitive to weight initialization [23]. As a result, the model is composed of 3 atomic embedding layers with mean Aggr whose sizes are {48, 57, 68} for the protein part and {48, 76, 121} for the ligand part, followed by 4 molecular embedding layers. All atomic parts have 2 heads for their respective GATv2Convs. The interaction fingerprint is a vector of size 189 fed to an MLP with channels {189, 94, 1}.

During training, a learning rate of 10^–3 associated with the Adam optimizer and a weight decay equal to 10⁻⁴ are used. The neural model is trained on 300 epochs with a dropout set at 0.13 (applied to all layers), but unpromising, stagnant candidates are stopped (with early stopping). The model was trained on an Nvidia DGX [a server with two sockets of 20 cores each, 512 gigabytes (GB) of random access memory (RAM)] using one graphics processing unit [GPU, an Nvidia V100 with 16 GB of video RAM (VRAM)], 300 epochs taking about 3 h00. The mean inference time is 0.10 s using one central processing unit (CPU), however, the necessary time to construct the graph from the complex file must be considered. That processing takes nearly 0.05 s (on the CASF 2016 complexes), but that time does not include the time for extracting the pocket if necessary, which is approximately 0.30 s per complex. The different times were not deeply studied because the other articles do not provide a basis for comparison. The results for the scoring power are presented in Figures 8 and 9 and Table 3. The scoring power results show a good R_p, especially for the CASF 2016, and both error metrics (MAE and RMSE) are pretty low, thus, the model is validated with respect to the scoring power.

Display full size

Figure 8. 

Results on CASF 2013. The blue points represent the predicted scores as a function of the targeted ones, the blue line is the blue points linear regression, and the orange line is the ground truth linear regression

Display full size

Figure 9. 

Results on CASF 2016. The blue points represent the predicted scores as a function of the targeted ones, the blue line is the blue points linear regression, and the orange line is the ground truth linear regression

Ranking power results are presented in Table 4, which shows that the model has good performance for the ranking power, especially on the CASF 2016 dataset.

Docking power results are presented in Tables 5 and 6. As shown in Table 5, the method ranks perfectly 1/3 of the native ligands (or a closer one) among the decoys. Moreover, the method ranks more than half of the native ligands (or a closer one) in the top 2 and 3. However, Table 6 does not show a good ranking performance; that lack of performance may be caused by the large number of used decoys and the small range of all scores (the average score range is 3.0) so the score range being small the model has more chance of not respecting the RMSD ranking. It should be noted that the fact that the produced scores are in a small range is because all the experimental scores of the training set are themselves in a small range.

Discussion

In addition to the intrinsic characterization of the new scoring method, the HGScore results are compared to some other DL scoring methods in Tables 7 and 8 respectively for CASF 2013 and 2016. The other methods compared are Pafnucy [4], DeepBindRG [5], OnionNet [7], IGN [10], and ECIF6::LD-GBT [11]. Some classical scoring function performances on CASF 2013 [24] and CASF 2016 [13] are presented to provide a no AI-centric comparison.

HGScore is among the best results on the CASF datasets for the scoring power metrics (Table 8), and the ability of the method to split the learning following the bond type may be an important part of the improvement. However, compared to existing scoring methods, the proposed new score does not bring significant improvement in the case of the CASF 2013 dataset (Table 7); it may be because the model of the present study has been trained on PDBbind version 2020. Indeed, the CASF datasets are the most representative complexes of PDBbind respectively in versions 2013 and 2016. Although the IGN method has a R_p of 0.83, the CASF 2013 test set is reduced to only 85 items which is far from the whole set, thus it is difficult to conclude about this R_p value. Concerning the ECIF6::LD-GBT method, the CASF 2016 is used for both the validation and test sets, whereas for HGScore the test set is unknown when the model is tested on it, which is the canonical way of proceeding. On the other hand, some authors have trained their method on the PDBbind refined set, which is composed of the highest quality complexes, but that choice may imply a lower prediction quality on the poorest quality complexes, which are the most frequent ones. It also should be noted that some methods do not provide any other metric than R_p, involving a greater difficulty in comparison.

The IGN method is the one that most resembles HGScore, indeed, the HGScore heterogeneous graphs can be considered as the assembly of the three independent graphs of the IGN method.

The present model has been trained on the 2020 version of PDBbind. The 2016 version of CASF is closer in terms of composition, which explains why HGScore improves the performance on the latter but not on the 2013 version. Nevertheless, it should be noted that on the CASF 2013, HGScore outperforms the others since it leads to lower error metrics values: the MAE and the RMSE.

Since the training set impacts the results of the scoring, it is of utmost importance to compare the performance of the HGScore with the performance of the score of other methods trained, validated, and tested identically. That is why the Pafnucy method [4] was also trained from scratch on the same HGScore dataset split. Tested on CASF 2016 (Table 8), the HGScore method still outperforms Pafnucy trained in a completely identical way (starred line in Table 8). Because of the difficulty to get the source codes and then to execute them easily, in a strictly identical way to HGScore, only the Pafnucy method has been retrained.

In conclusion, the new proposed DL method allows scoring protein-ligand complexes. Based on the Attentive FP method, the presented model brings the originality to split the learning process regarding the interaction class, ligand to ligand, protein to protein, ligand to protein, and protein to ligand. That general idea allows the model to learn more precisely the effect of the molecules themselves and the interaction between them. The method does not propose a new type of layer but rather a new way to represent ligand-protein complexes, using a heterogeneous graph. Trained on the PDBbind dataset, like almost all counterpart methods, HGScore reaches the top of the methods assessed on the CASF 2016 dataset. Moreover, the method largely outperforms the classical scoring methods presented in a previous article [1]. Although HGScore does not bring a huge improvement on the main test metric (Rp on CASF 2016), it outperforms all methods blindly tested on CASF. Moreover, the quality already achieved by the metrics implies that their improvement cannot be very large anymore.

As mentioned previously the network is split according to the interaction category. This method will allow to easily use transfer learning and thus reuse the feature extractor and high-level abstract representation based on the type of interaction between the two molecular groups. Of course, it will be necessary to ensure that the roles of the receptor (here the protein) and the ligand are preserved. Therefore, the use of transfer learning to adapt this new model to cases of protein-glycosaminoglycan (GAG) or protein-proteoglycan docking is very promising.

The next objective will be to integrate that scoring function into a complete docking process and analyze a) its impact on docking results and understand the root cause for better performance over other docking programs such as AutoDock and b) its impact on computation time. For now, HGScore can be used as an external scoring function allowing to reassess the score of a set of docked complexes to refine existing ranking.