Interpretable Drug Target Prediction Using Deep Neural Representation

1 Introduction

The identification of drug-target interactions (DTI) is a key task in drug discovery, where drugs are chemical compounds and targets are proteins. The high profits on patented new drugs have motivated pharmaceutical labs to examine most potential interactions. Since experimental assays take time and are expensive, efficient computational methods for predicting and understanding potential DTIs are useful and urgently demanded.

Two major computational approaches have been investigated for DTI prediction: molecular docking and machine learning. Molecular docking using three-dimensional (3D) simulation is widely used for its reasonable accuracy and visual interpretability. However, it suffers from two significant shortcomings: first, getting 3D structures of proteins itself is a challenging task; second, large scale simulation can be time consuming.

Meanwhile, machine learning approaches are very promising as they enable large scale testing of candidates in a relatively short time. DTI prediction can be viewed as a binary classification problem, where the input is a pair of a drug candidate and a protein and the output label indicates whether there is an interaction between them. To the best of our knowledge, three main approaches have been investigated to represent DTI pairs: (1) chemical and biological descriptors such as molecule fingerprints or amino acid sequence descriptors [Faulon et al., 2007]; (2) an aggregation of expert designed similarity measures to derive candidate pairs from known DTI pairs [Ding et al., 2013; Yamanishi et al., 2010]; and (3) representations learned from descriptors from (1), for example, using techniques such as Restricted Boltzmann machine or Autoencoder [Wen et al., 2017; Chan et al., 2016]. Although relatively effective, most methods are usually black-boxes and less biologically interpretable.

In this paper, we propose an interpretable end-to-end neural network model that predicts DTIs directly from low level representations. Specifically, the inputs of the model are raw amino acids sequences and chemical structures, and it produces structure-level interpretations in addition to the DTI predictions. As shown in Figure 1, we use long short term memory recurrent neural networks and graph-based convolutional neural network to project proteins and drug structures into dense vector spaces. A two-way attention mechanism (shown as αpi and αdi) is used to calculate how the pair interacts and thus enables interpretability. Finally, the attentionbased vector representations are used by a classifier, a simple sigmoid function in this paper, to make a prediction. We have also shown that our model is extensible to incorporate highlevel information such as Gene Ontology annotations.

In the experiments, we pay special attention to the method of constructing testing examples. Our testing dataset is constructed in a way that simulates the practical situations, where, given a pair of drug and protein at testing time, the drug, the protein, or both of them may have not been observed at training time. Such experimental setting demands great generalization ability. Compared with previous methods, our model yields superior results, while using less feature engineering and domain expertise, specifically in the difficult cases that were not covered well by human-designed features and where neither the drug nor the protein from a testing pair is observed at training. At the end, we present a case study to demonstrate the visualizable interpretation and its usefulness during the drug discovery process.

2 RelatedWork

Predicting drug-target binding has been an interesting topic to the pharmaceutical industry and researchers. Molecular docking [Luo et al., 2016; Trott and Olson, 2010] is widely used to predict the binding mode and score given the 3D structure inputs of a drug molecule and a target. In order to make the binding prediction, the docking program looks for the optimal binding position of the drug molecule inside the binding pocket of the target and estimates the binding score according to predefined force fields. Though molecular docking is a popular tool for high throughput screening, it takes time for large scale predictions and it is limited by the availability of the 3D structures of the protein targets.

Machine learning methods have been implemented to predict DTIs. [Ding et al., 2013] reviewed similarity-based machine learning methods based on sequence, protein-protein interaction (PPI) network and Gene Ontology (GO) semantic information. Their machine learning approaches include nearest neighbor, bipartite local models, pairwise kernel method, kernelized Bayesian matrix factorization, networkbased inference and so on. [Faulon et al., 2007] developed support vector machine (SVM) models to predict DTIs and catalytic effect based on chemical structures and enzymatic reactions. [Yamanishi et al., 2010] integrated chemical, genomic and pharmacological data together to predict drugtarget binding via similarity-based methods.

With the increasing popularity of deep learning, researchers are adopting deep neural models to predict DTIs. [Wen et al., 2017] developed a Deep Belief Networks (DBN) model consisting of stacked Restricted Boltzmann machines (RBM). For a DTI pair, the inputs are Extended Connectivity Fingerprints of the drug and sequence composition descriptors of the protein. First the DBN is pre-trained in an unsupervised manner using only the training feature vectors, and then it is fine tuned with both feature vectors and labels from the training dataset. Instead of using RBM and DBN, [Chan et al., 2016; Wang et al., 2018] used stacked Autoencoder for representation learning and SVM or rotation forest for classification.

Although aforementioned (deep) machine learning methods have proved to be able to make relatively effective predictions, the lack of interpretability limits their practicality from biological perspectives. Recently, differentiable representation learning methods that can be directly applied on low-level representations enable the potential of interpretable DTI predictions. For example, [Altae-Tran et al., 2017; Duvenaud et al., 2015] explored using graph convolutional network to model chemical structures, and while it is intuitive to apply recurrent neural network (RNN) on protein sequences [Pollastri et al., 2002], [Schwaller et al., 2017] also used RNN to model SMILES strings, which are sequential encoding of chemical structures.

To the best of our knowledge, this paper presents the first end-to-end deep machine learning work that produces interpretable DTI predictions directly from low level representations. By learning directly from molecular structures and protein sequences, our approach saves the effort of designing biochemical descriptors or similarity measures, both of which can be expensive processes of feature engineering.

3 Model

3.1 ProblemFormulation

On the one hand, a protein sequence consists of a list of amino acids p = (a1, . . . , an), where ai is one of the 23 types of amino acids (20 standard, 2 additional, and 1 for unknown). Additionally, each protein has a set of gene ontology (GO) annotations [Ashburner et al., 2000] GOp = {g1, . . . , gm} that provides high level information (e.g., protein function).

On the other hand, a drug is represented by a SMILES [Weininger, 1988] sequence, which essentially encodes a chemical structure graph d = {V,E}, where V is a set of atoms and E is a set of chemical bonds that bind two atoms as undirected edges. We use RDKit 1 to transform SMILES string to chemical structure graph.

The goal of DTI prediction is to learn a model that takes a pair (p, d) as input and outputs y ∈ {0, 1}, where y = 1 indicates an interaction between p and d.

3.2 Recurrent Neural Network

In the situation where proteins are represented by amino acid sequences and drugs are represented by SMILES strings, we use recurrent neural network (RNN) to project sequential input to dense vector representations. Specifically, because, in reality, protein sequences fold in 3 dimensional space and because SMILES strings are contextual by design, both of which can be viewed as long-distance dependencies, we use long short term memory (LSTM) RNNs [Hochreiter and Schmidhuber, 1997] for their ability to memorize long term information. At each time step t, the LSTM unit takes the t-th input token embedding xt ∈ RM , the hidden states from the previous time step h(t-1) ∈ RH , and the cell states from previous step c(t-1) ∈ RH . Then it produces new hidden and cell state ht, ct:ht, ct = LSTM(xt, ht-1, ct-1)

Here, M and H are two hyper parameters for the dimension of the embedding space and the dimension of the hidden space respectively. We initialize h0 = c0 = 0H as a vector of zeros. Suppose the input tokens belong to a vocabulary V = {t1, · · · , t|V |}, the input embeddings are obtained as xi = W T v Ii where Wv ∈ R|V |×M is a learnable parameter and Ii ∈ R|V |×1 is a vector whose i-th value is one and all others zero.

3.3 Convolutional Neural Network on Graph

When drugs are represented by chemical structural graphs, we use convolutional neural network (CNN) to project chemical structural graphs to dense vector representations. It is more intuitive than using RNN to model them because it eliminates the step of linearizing the graph structures. [Duvenaud et al., 2015] recently demonstrated that, as a differentiable generalization of circular fingerprint, a CNN based neural fingerprint provides better descriptive drug modeling in a datadriven manner.

1 http://www.rdkit.org/

Algorithm 1: Pseudocode of graph CNN.

Input: Molecule graph G = (V, E), radius R, hidden weights H1 1 . . .H 5R Output: A vector ra for each atom a Initialize: for each node a ∈ V do ra ← g(a) // Look up initial feature vectors end 1 for L = 1 to R do 2 for each node a ∈ V do 3 N = neighbors(a) 4 v ← ra + ∑u∈N ru 5 ra ← σ(vH|N| L ) 6 end 7 end

Algorithm 1 shows the pseudo-code of our variant of the neural fingerprint algorithm that produces a dense vector representation for each atom in the molecule. We omit the steps required to calculate a single vector fingerprint, because instead we will exploit atom vectors with attention mechanism (Section 3.4) to enable interpretability on the drug side. At the initialization phase, as described in Section 4.2 of [Duvenaud et al., 2015], the atom features are initialized as a 62 dimension sparse vector that indicates both chemical and topological properties of the atom. Then the algorithm iteratively applies convolutional operation on the graph R times and updates the atom vectors. The radius parameter R controls how many hops information can be propagated to. In this paper, R = 3.

While CNNs are usually applied on tensors, e.g. images, this algorithm is convolutional in the sense that it applies filters to each atom and its neighborhood to capture local signals, and then the aggregated local signals are pooled to get the final vector representation. Different from images where each pixel always has 8 neighbor pixels, an atom can have from one to five neighbor atoms. So instead of using one convolutional filter, the algorithm uses 5 types of linear filters H1 . . .H 5 for atoms with corresponding number of neighbors.

3.4 Attentive Pooling Network

Neural networks with attention mechanism have been effectively applied to vision tasks (e.g., image captioning) and natural language processing tasks (e.g., machine translation), where the output components selectively focus on subsets of the input based on attention weights. Extending the one-way attention for pairwise inference, attentive pooling network [dos Santos et al., 2016] provides a two-way attention mechanism that enables the input pairs to be aware of each other.

Suppose P ∈ RHp×Lp is the context matrix of a given protein, where Hp, Lp are the dimension of the protein hidden space and the number of inputs, it can be formed in 3 ways as proteins have two input sources: (1) columns of the matrix P are the LSTM hidden vectors with the amino acids sequence as input so that Lp equals the number of amino acids in the sequence, (2) columns of P are GO annotations embeddings so that Lp equals the number of GO terms for the protein, and (3) the concatenation of both (1) and (2). Note that since there is no order between GO terms, for situation (2) and (3), the embeddings of GO terms are fed to the attention module directly without going through RNN.

Similarly, suppose D ∈ RHd×Ld is the context matrix of a given drug, Hd, Ld being the dimension of the drug hidden space and the number of inputs. The columns of D can be (1) the LSTM hidden vectors with SMILES string as input so that Ld equals the number of tokens in the SMILES string or (2) the atom vectors obtained from graph CNN so that Ld equals the number of atoms in the molecule.

A soft alignment matrix A ∈ RLp×Ld is calculated as

where U ∈ RHp×Hd is a trainable parameter. For an intuitive example, when proteins are represented by amino acid sequences and drug chemical structure graphs, A presents the interaction between each amino acid and each atom.

Next, the attention weights αp ∈ RLp , αd ∈ RLd , which can be interpreted as importance scores on the input units, are calculated by applying row-wise and column-wise maxpooling operations to A.

Finally, αp and αd are normalized by Softmax function, and the results of which are used as weights to weighted sums of the context vectors:

where the softmax function is defined as

3.5 InferenceWith Siamese Network

A Siamese network [Bromley et al., 1994] has two input multi-layer networks and one output whose value corresponds to the similarity, possibility of interaction in the case of this work, between the input pair. As shown in Figure 1, two networks with 3 linear layers and 2 rectifier layers are used. To reduce the hyper-parameter space, we require all the linear layers to have the same input and output dimension Hs except the first one, whose input dimension corresponds to previous outputs.

The attention based vector representations rp and rd are fed separately into the two networks. Then we take the inner product of the outputs and use a sigmoid function to predict the probability that an interaction exists between a pair of protein and drugvp = fp(rp) and vd = fd(rd) (4)

where fp, fd are the transformations of the siamese networks for protein and drugs respectively.

In a classification scenario, a hyper-parameter threshold δ is selected as classification boundary

3.6 Training

Given a dataset D = {(pi, di, yi)}, i = 1 · · · n, the model can be trained by maximizing the likelihood of observing the training data, which is equivalent to minimizing the cross entropy loss function

(7)

where Θ is the set of neural network parameters.

However, although, in this paper, we use a dataset with both positive and negative pairs as described in section 4, negative pairs (yi = 0) are usually not available for similar tasks especially when a dataset is from a knowledge graph that stores only existing triples. Therefore, we employ a pairwise ranking loss [Bordes et al., 2011] that, for each given protein p, maximizes the margin between interacting drugs and non-interacting drugs, i.e. ranking positive drugs higher than negative drugs as much as possible.

where γ > 0 is a hyper-parameter that specifies the width of the margin, and N +(p) and N -(p) give the set of drugs that interact with p and those do not interact with p respectively. In this setting, negative examples can be generated by sampling pseudo-negative drugs with heuristic criteria if a dataset does not have any.

4 Experiments

4.1 Dataset

BindingDB [Gilson et al., 2016] is a public, web-accessible database of measured binding affinities, focusing chiefly on the interactions of small molecules (drugs/drug candidates) and proteins (targets/target candidates). We took a snapshot of BindingDB that contains 1.3 million data records, each of which contains information such as the identifiers of involved entities, the observed experiment results, etc. By the following criteria we construct a binary classification dataset 2 with 39,747 positive examples and 31,218 negative examples.

1. Record has chemical identifier (PubChemCID), and the small molecule has chemical structure represented by SMILES 3 .
2. Record has protein identifier (Uniprot ID), and the protein has both sequence representation and Gene Ontology annotations [Ashburner et al., 2000].
3. Record has IC50 value, a primary measure of binding effectiveness.
4. The chemical molecule weight is less than 1,000Da, due to our focus on small molecule drugs.
5. By following the activity threshold discussion in [Wang et al., 2016], record is positive if its IC50 is less than 100nm, negative if IC50 greater than 10,000nm.

Suppose Ptrain and Dtrain are the sets of proteins and drugs that are observed in the training dataset. [Pahikkala et al., 2014] suggested that there are 4 experimental settings under which the model can be learned and applied to predict the label between a drug (candidate) and a protein target d and p:

1. Both p and d are observed in the training dataset: p ∈ Ptrain and d ∈ Dtrain;
2. The protein p is observed in the training dataset but the drug d is not: p ∈ Ptrain and d /∈ Dtrain;
3. The drug d is observed in the training dataset but the gene p is not: p /∈ Ptrain and d ∈ Dtrain;

2 https://github.com/IBM/InterpretableDTIP

3 Simplified molecular-input line-entry system [Weininger, 1988]

1. Neither p nor d is observed in the training dataset: p /∈ Ptrain and d /∈ Dtrain.

We split proteins and drugs into those that should be observed in training and those that should not with four experimental settings; we then allocate DTI pairs into training, development, and testing datasets. Statistics of the datasets are shown in Table 1.

4.2 Training Details

During training, the parameters are initialized randomly from an uniform distribution Θ ∼ (-0.08, 0.08). In each step, with batch size equals 32, a batch of proteins or drugs is randomly selected from the training data. Then we retrieve positive and negative interactions that involve them as input to the model. Dropout is applied on the output of embedding lookups and between Siamese network layers. According to Equation 5 and Equation 8, we use Adam gradient descent optimization with initial learning rate equals to 0.001 to train the parameters. We train the model for 30 epochs, where each epoch consists of 100 steps. The model is evaluated using the validation dataset after each epoch, and the one with the best ROC score is reported.

The validation set is also used to select model configurations, hyperparameters, and the classification boundary δ. As mentioned in Section 3.4, there are 3 ways to represent proteins and 2 ways to represent drugs. We found that the best model uses both amino acid sequences and GO annotations for proteins, and it uses atom vectors from graph CNN for drugs. We use gradient-boosted-tree search from scikit-optimize 4 to find effective hyperparameters. All the space dimensions are selected from powers of two that are between 8 and 64, embedding dropout probability from 0 to 0.2, Siamese dropout probability from 0 to 0.5, and training margin γ from 0.0001 to 0.5. All ranges are inclusive. The

4 https://scikit-optimize.github.io/

values of hyperparameters of the best model are shown in Table 2, and the best classification boundary is δ = 0.4995.

4.3 Baselines

We compare our system with 3 baselines: matrix factorization, a similarity-based machine learning method, and a previous deep learning based method.

Matrix Factorization Instead of viewingDTI prediction as a classification task, it can also be treated as a collaborative filtering task. Given the IC50 values of the known pairs in our dataset, the task is to predict the IC50 values of other pairs. From this perspective, matrix factorization has been widely used and proved effective [Koren et al., 2009]. We use the implementation from LIBMF 5 as the first baseline.

Note that although this baseline does not use content-based information from amino acids sequences and molecules, it uses extra information from the IC50 values, including those of neither positive nor negative examples, from BindingDB. The raw IC50 values are truncated at an upper threshold 100,000 and fed into LIBMF in their logarithmic values.

Similarity Based Method We also compare our approach against a similarity based prediction approach. In particular, we use as our second baseline a state-of-the-art similarity based prediction system adapted to DTI predictions: Tiresias [Fokoue et al., 2016]. Tiresias system takes as input a pair of drug and protein, builds a feature vector using the statistics of a large number of similarity measures against known DTI pairs, and outputs the prediction using a logistic regression classification model.

Deep Learning Based Method At last, we compare our approach with [Wen et al., 2017] as mentioned in Section 2. Its hyperparameters, including the number of units of RBM layers, the number of pre-train and fine tuning epochs, and the learning rates are obtained by the same gradient-boostedtree search as for our model. We refer this method as DBN in the rest of this paper.

5 Results

5.1 Effectiveness

Before diving into the results, recall that different amounts and levels of information are exposed to the systems. Since matrix factorization is not content-based, it does not learn from any low-level representation and cannot handle unobserved instance, and so it is only evaluated in the first setting. However, it uses extra information from the raw IC50 values and uses more training examples whose IC50 values (larger than 100nm and smaller than 10,000nm) fall outside of our classification range. Tiresias, as the state-of-the-art similarity-based system, uses a set of expert designed similarity measures based on protein and drug properties. For deep learning approaches, DBN uses middle-level features from expert designed molecular fingerprints and protein descriptors, whereas our best model (E2E) learns from molecular chemical structures, amino acids, and GO annotations in an end-to-end manner. For a fair comparison with DBN, we also include a version of our model without using GO annotations (E2E/GO).

5 https://www.csie.ntu.edu.tw/ cjlin/libmf/

Figure 2 presents the performance comparison of our model and the baselines. As mentioned in Section 4.1, the systems are evaluated in four different testing datasets, to which the systems have different visibility of entities in the training data. In terms of overall performance, our best system (E2E) is the only one that consistently performs well across all datasets and all metrics (all greater than 0.8), whereas the performance of the baselines drops dramatically when the proteins are not observed in the training data. If we aggregated the 4 cases, our method outperforms other baselines in all metrics, e.g. the average AUPR of E2E is 0.91 while that of DBN is 0.81.

The baselines are more effective than our model when the tested proteins are observed. However, when the tested proteins are not observed, Tiresias fails as its similarity measures for proteins are not effective, DBN loses to its counterpart E2E/GO in most cases except for its accuracy in the third setting, and, finally, our best model E2E introduces more improvement on top of E2E/GO by exploiting its flexibility of incorporating high level information.

5.2 Case Study for Interpretability

The key advantage of our model over all baselines is its interpretability. Before conducting costly lab tests on potential interactions, it is helpful for researchers to gain insights by telling them where to look at.

To demonstrate that, we conducted a case study of a top predicted interaction between chemical (PubChem ID: 117793281) and protein (UniProt ID: Q9HBH9) by comparing our interpretation with that obtained by molecular docking, which is currently more biologically interpretable. We use molecular docking to generate their structural interaction pattern and the result was analyzed in Discovery Studio Visualizer 2017 R2 (Figure 3). We found that, according to the attention mechanism, the top seven contributing protein residues are either within or surrounding the protein binding pocket, and two of them (GLU 160 and PHE 150 as labeled in Figure 3) may closely interact with the drug. From the drug's perspective, among the most contributing atoms (marked in saturated red), three are nitrogen atoms, which have potential to formulate hydrogen bonds, and one is a carbon atom that has a relatively big solvent accessible surface, indicating potential solvent interaction can be formulated. Thus, our model gives reasonable cues on the factors for the binding, which may have broad pharmaceutical applications.

6 Conclusion

We have presented an interpretable end-to-end deep learning architecture to predict drug-target interactions from low level representations. Experimental evaluation shows that this approach overall outperforms all baselines, and it is the only one capable of generalizing well to new proteins (i.e., not seen in the training data), which is critical for drug discovery given that only a small fraction of proteins are known to be targets of chemical compounds. Furthermore, we have illustrated, via a case study, the ability of our approach to provide biological insights to understand the nature of predicted interactions.