Evidential deep learning-based drug-target interaction prediction

Datasets

Three datasets DrugBank⁶⁸, Davis⁵⁹, and KIBA⁶⁰ were used for model evaluation. Supplementary Table 1 summarizes these datasets.

DrugBank: the data used in this study were published on January 3, 2020 (version 5.1.5). Drugs with inorganic compounds, very small molecule compounds [e.g., iron (DB01592), zinc (DB01593)], or SMILES strings that the RDKit Python package could not recognize were removed. Drugs with a SMILES length greater than 300 were also removed to comply with the maximum input of MG-BERT. Ultimately, 6148 drugs and 4085 targets, totaling 16,531 pairs of known drug-target interactions, were collected, and these known drug-target interactions were used as positive samples for training. We randomly disrupted these positive pairs to generate negative samples, specifically by randomly selecting a drug-target interaction pair and forming a negative drug-target interaction pair between the target in that pair and another drug with a known drug-target interaction pair. Negative samples with the same number of positive samples, totaling 16,531 negative sample pairs, were generated through the above steps. Thus, a total of 33,062 pairs of samples were collected in the DrugBank dataset, with a ratio of positive to negative samples of 1:1.

Davis and KIBA: The Davis and KIBA dataset records wet-lab values for determining the binding affinity between drugs and proteins. We refer to previous work and divide positive and negative samples using 5.0 and 12.1 as thresholds to construct a binary classification dataset.

In order to assess the ability of the model to predict novel compound interactions in real-world drug discovery, we implemented a cold-start scenario following the practice established by Wang et al.¹² Selected 10% of drugs from DrugBank dataset, and then all the DTI pairs associated with these drugs were used as a test set. Resulting in 3170 cold-start test samples VS. 29,894 training samples.

To perform the unseen data prediction (see “Uncertainty prediction accelerates drug discovery”), we collected all new drugs approved by the FDA in 2022 and their targets from the DrugBank dataset. In 2022, the FDA approved 37 new drugs for marketing, including 22 new molecular entities and 15 new biologic applications. Only molecular drugs were selected among these new approvals, while drugs already appearing in the training set were removed. Ultimately, we collected 13 drugs, 22 targets, and 24 pairs of interactions. We used the same method for constructing the negative samples as for the DrugBank training set, randomly generating as many negative samples as positive samples. Thus, we collected 48 pairs of samples for out-of-domain performance testing. In these samples, the ratio of positive to negative samples was 1:1.

Framework architecture

Protein feature encoder

The protein feature encoder consists of an initial encoding module and a LA module⁸¹.

In the initial coding module, we extracted sequence information using the widely used protein language model ProtTrans³⁰ (denoted as ProtT5), which is a self-supervised autocoder based on the transformer model. A 1024-dimensional sequence embedding for each residue was extracted using ProtT5. For each protein, we obtained a \(L\times 1024\)-dimensional feature matrix, where L is the length of the protein chain. We did not train ProtT5, but used its locked parameters to obtain the representations.

The protein feature matrix is encoded using light attention. The input to light attention is a protein embedding \(x\) of size \(L\times {d}_{{in}}\), where \(L\) is the length of the protein chain, i.e., the number of amino acids, and \({d}_{{in}}\) is the length of the initial representation of each amino acid, in this case, 1024. Two independent 1D convolutions transform the input. Convolution is used in the length dimension to generate attention coefficients and values \(e,v\).

where \(b\) is a learned bias, \({s}\) is filter sizes, \({W}^{(e)}\) is the weights of the 1D convolutions.

To obtain the attention weights, softmax normalization is performed on \({e}_{{ij}}\). The attention weight \({\alpha }_{{ij}}\) for the j-th residue and the i-th feature dimension is computed as:

The weight distribution of each feature dimension i is independent, and they can produce different attention patterns. The attention distribution calculates a weighted sum of the transformed residual embeddings \({v}_{{ij}}\). Thus, we obtain a representation \(x{\prime}\) independent of protein length.

Drug feature encoder

The drug feature encoder consists of a 2D topology encoding module and a 3D structure encoding module.

In the 2D topology encoding module, the small molecule pre-training model MG-BERT⁵⁸ is used to extract drugs’ 2D topological graph information. MG-BERT is a pre-trained model based on a self-supervised learning approach for learning atomic-level representations of drug molecules. The model uses the transformer architecture and employs a mask-based self-supervised learning approach. Again, we did not train MG-BERT but used its frozen parameters to obtain an initial representation of the drug. In this work, the representation obtained by MG-BERT is extracted as an initial representation of the drug molecule, which is used as input to 1DCNN to obtain a final representation of the 2D information of the drug.

In the 3D structure encoding module, given a 3D structure graph, an atom-bond graph and a bond-angle graph are constructed, under which the representation vectors of atoms and bonds are learned iteratively⁸². In the atom-bond graph, the nodes of the graph are atoms, and the edges of the graph are covalent bonds of atomic bonds. In the bond-angle graph, the nodes of the graph are the bonds, and the edges are the bond angles between the two bonds. The scalar atomic distances and angles are mapped into high-dimensional vectors using the Gaussian kernel function. In a more concrete sense, the representation vectors of atom \(u\) and bond \((u,v)\) for the kth iteration are denoted as \({h}_{u}\) and \({h}_{{uv}}\), respectively. The initialization is set as \({h}_{u}^{(0)}={x}_{u}\) and \({h}_{{uv}}^{(0)}={x}_{{uv}}\).

Given the bond \((u,v)\), its representation vector \({h}_{{uv}}^{(k)}\) at the kth iteration is formulated by

\(N(u)\) and \(N(v)\) represent the adjacent atoms of u and \(\,v\),\(\{(u,w):w\in N(u)\}\cup \{(v,w):w\in N(v)\}\) are the neighboring bonds of \((u,v)\).

The function “Aggregate” aggregates messages, while “Combine” updates the bond-angle graph.

The information from adjacent bonds and corresponding bond angles is aggregated into \({a}_{{uv}}^{\left(k\right)}\), and subsequently, the bond representation \(\left(u,v\right)\) is updated based on the aggregated information.

For an atom \(u\), at the kth iteration, its representation vector \({{h}}_{u}^{\left(k\right)}\) can be expressed as:

The messages from the bonds are learned from the bond-angle graph. The aggregated messages then update the representation vector of atom u.

The atom’s representation vectors at the final iteration are integrated to obtain the molecular representation vector \({h}_{G}\) through the READOUT function.

Evidential layer

The learned drug and protein representations are concatenated and serve as inputs to the evidential layer. The evidence layer is a fully connected neural network with N layers. Unlike classical neural networks, the softplus layer is used to ascertain non-negative output. In the deep evidential classification model, the joint evidence distribution is presented as a Dirichlet distribution, which is used to represent the interaction prediction and the degree of evidence associated with that prediction. The outputs of the evidence layer are used as parameters of the predicted Dirichlet distribution. For the DTI prediction task, the network has two outputs.

Evidential learning for classification

In subjective logic (SL), for each class \(k=1,\,.\,..,K\), SL will provide a belief mass b_k and an overall uncertainty mass u to compose a framework for evaluating the K classification. These mass values are all non-negative and sum to 1, i.e.,

where \({{\rm{u}}}\ge \,0\) and\(\,{{{\rm{b}}}}_{k}\,\ge \,0\) for \({k}=1,\,.\,..,K\).

The b_k value is computed from the amount of evidence for the k^th class. Let \({{{\rm{e}}}}_{k}\,\ge \,0\) be the evidence for the kth class, then the belief b_k and uncertainty u are calculated as

where\(\,S=\,{\sum }_{i=1}^{K}{e}_{i}-1\).

Here, the evidence is defined as a measure of the support from the data collection to a certain class to support the sample classification. A belief mass assignment, i.e., subjective opinion, corresponds to a Dirichlet distribution with parameters \({\alpha }_{k}={e}_{k}+1\).

That is:

Where \({{\rm{S}}}=\,{\sum }_{i=1}^{K}{\alpha }_{i}\)

The expected probability of the K^th class is the mean of the corresponding Dirichlet distribution, calculated as

The EviDTI is trained end-to-end using a multi-objective loss function \({{\boldsymbol{L}}}(\varTheta )\) to achieve the maximum fit and minimum incorrectly predicted evidence for the model.

Where \({\lambda }_{t}\) is the annealing coefficient, t is the index of the current training epoch, \({D}({p}_{i}\left|\left\langle 1,…,1\right\rangle \right.)\) is the uniform Dirichlet distribution, and \({\widetilde{\alpha }}_{i}={y}_{i}+(1-{y}_{i})\odot {\alpha }_{i}\) is the Dirichlet parameters after removal of the non-misleading evidence from predicted parameters \({\alpha }_{i}\) for sample i. More details about the Evidential Deep Learning to Quantify Classification Uncertainty.

Comparison methods

For the DrugBank, Davis and KIBA dataset, we retrained and tuned all models based on the parameters provided in their respective GitHub repositories. All methods were conducted with five repeated experiments, and the results are reported as the mean and variance.

DeepConv-DTI¹⁹

DeepConv-DTI utilizes CNNs to extract features from proteins and fully connected networks (FCNs) to extract features from drugs, followed by making predictions.

GraphDTA⁶⁴

GraphDTA employs GNN and CNN to represent drugs and proteins, respectively. As Nguyen et al. reported, the GAT_GCN combination is chosen as the feature extractor in subsequent experimental comparisons.

For a fair comparison, GraphDTA were adapted according to the following steps:(1) A Sigmoid activation function was added after the last layer to convert the model output to binary probabilities. (2) A binary cross-entropy was used as the loss function instead of the original mean square error (MSE) loss. (3) A grid search was performed to optimise key hyperparameters such as learning rate and batch size, and a validation set was used to determine the best configuration.

MolTrans²⁰

MolTrans employs the transformer model to encode protein and drug sequences into feature embeddings. These embeddings are used to construct an interaction matrix, which is subsequently processed by CNNs and FCNs to make predictions.

TransformerCPI⁶⁵

TransformerCPI utilizes the Transformer architecture, treating drugs and proteins as distinct sequences. By employing CNN and GCN, the model generates representations for protein sequences and atom structures, respectively. Subsequently, TransformerCPI extracts interaction features through the transformer decoder and employs linear layers to produce the final output, representing the interaction probability.

HyperAttentionDTI¹⁷

HyperAttention DTI is an end-to-end bio-inspired model built upon the CNN and attention mechanism. Deep CNNs are employed to learn feature matrices for both drugs and proteins adeptly. The attention mechanism was leveraged on these feature matrices to capture intricate non-covalent intermolecular interactions among atoms and amino acids, assigning an attention vector to each atom or amino acid.

AIGO-DTI⁶⁶

AIGO-DTI proposes an Adaptive Iterative Graph Optimization (AIGO)-DTI prediction framework. This framework integrates atomic cluster information and enhances molecular features through the design of functional group prompts and graph encoders, thereby optimizing the construction of DTI association networks. Furthermore, the optimization of graph structure is transformed into a node similarity learning problem, utilizing multi-head similarity metric functions to iteratively update the network structure to improve the quality of DTI information.

GraphormerDTI²⁴

GraphormerDTI integrates the Graph Transformer neural network with a 1D-CNN to extract the representations of drug and target, respectively. Specifically, GraphormerDTI embeds molecular graphs into vector representations through iterative Transformer-based message passing. The structural representations of the molecules are encoded through node centrality encoding, node spatial encoding and edge encoding. The attentional mechanisms are used to model interactions.

DLM-DTI²⁷

DLM-DTI was comprised of three primary components: the drug encoder, the target encoder, and the interaction prediction head. The target encoder incorporates both the teacher and student models of language models for protein sequences. The teacher model employed for target sequence encoding was the ProtBERT model. The student model was designed to be consistent with the original teacher model with fewer layers to conduct knowledge distillation. The function of the drug encoder is to convert SMILES sequences to molecular representations, and this is achieved by employing the ChemBERTa encoder. The interaction prediction header predicts interaction probabilities by concatenating drug and target representations.

Statistical analysis of binding sites

To verify the correlation between residues with higher attention values in LA and drug-target binding sites, we obtained the real three-dimensional structures of the DTI interaction pairs in the DrugBank dataset from the PDB database for statistical analysis of binding sites.

First, drug-target pairs lacking the true structure of the protein were excluded. Next, for each protein, only structures containing small molecules that interact with them in the PDB file were retained. The binding sites of the drugs were then calculated, defining a binding residue as any atom of that residue within a 4 Å Euclidean distance of the drug molecule. Finally, to ensure accurate mapping of attention values to residues, the amino acid sequence indexes in the PDB were aligned with those in UniProt, which were used during training.

In the Light Attention network, for a given amino acid sequence of length L, the network will output an \(L\times {d}_{{in}}\) dimensional attention vector, where \({d}_{{in}}\) denotes the dimension of the attention vector. We get the attention value for each amino acid by taking the average of this matrix in the\(\,{d}_{{in}}\) dimension. A larger attention value means that the amino acid contributes more to DTI prediction.

Experimental settings

EviDTI is implemented in Python 3.8 and PyTorch 1.12.0, as well as functions in PYG 2.2.0, Scikit-learn 1.2.2, Numpy 1.24.3, Pandas 1.4.2, and RDKit 2020.09.1. The batch size was set to 32, and the Adam optimizer was used with a learning rate of 5e-4 and a weight decay of 1e-5. We allowed the model to run for a maximum of 500 epochs for all datasets, and if the model’s accuracy on the validation set did not decrease within 30, training would stop. The parameters of all neural networks are listed in the Supplementary Table 11 and Supplementary Fig. 9.

Compounds collection from Targetmol-Tyrosine kinase modulator library (L2200)

There were 746 potential tyrosine kinase modulators collected from the Targetmol-Tyrosine kinase modulator library (L2200). To ensure the credibility of the experiments, we excluded kinase modulators with drug similarity greater than 0.7 to the DrugBank dataset and ultimately retained 52 potential modulators. Drug similarity was calculated based on the Tanimoto coefficient.

In vitro kinase activity assays

The FAK (Cat#V1971) and FLT3 (Cat#V4064) were purchased from Promega. Compounds Tyrphostin 9 (Cat#T2479), GSK180736A(Cat#3513), Vodobatinib (Cat# T8882), BIX02188 (Cat# T1744), PKG drug G1 (Cat# T4661), ST271 (Cat# T4511), Flumatinib mesylate (Cat# T7861), Mubritinib (Cat# T6124), PF-562271 (Cat# T6124) and Sorafenib (Cat# T0093L) were purchased from Targetmol, USA. All compounds were tested from 1000 nM, with 3-fold dilution for 10 points.

In vitro kinase activity assays were conducted through ADP-Glo assays (#V9101) provided by Promega. The protocol for the FAK and Flt3 assay is described as follows: Enzyme, substrate, ATP, and compounds were diluted in 1X reaction buffer composed of 40 mM Tris (pH 7.5), 2 mM MnCl, 100 µM sodium vanadate. In a 384-well low-volume plate, l μl of the compound at indicated doses or 5% dimethyl sulfoxide (DMSO), 2 μl of FAK or FLT3 enzyme (15 ng/well), and 2 ul of substrate/ATP mix (final concentration: 20 μM ATP) were added to each well. The plate was then incubated at 25 °C for 60 min to allow for kinase activity. Following the enzymatic reaction, 5 μl of ADP-GloT Reagent was added to each well, and the plate was incubated at 25 °C for an additional 40 min. Subsequently, 10 μl of Kinase Detection Reagent was added to convert ADP to ATP and introduce luciferase and luciferin to detect ATP, and the plate was incubated for a final 30 min at 25 °C. Luminescence was recorded with an integration time of 0.5 s using the SpectraMax iD3 instrument (Molecular Devices).

The 50% effective concentration was calculated using Prism 8 by fitting the following equation:

where X is a log of concentration, Y is a response, and top and bottom are the responses of controls, each assay was repeated at least three times, and we computed the mean and standard deviation for the values.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.