Learning Data-Driven Drug-Target-Disease Interaction via Neural Tensor Network

1 Introduction

Data-Driven Drug Discovery. Personalized medicine recommendation is one of the most promising assets to treat human disease [Li et al., 2015]. A critical step in personalized medicine is to understand drugs' mechanism of actions (MoAs) by exploring the biological interactions among drugs, targets, and diseases. In vitro experiments can be performed to identify potential associations of drug-targetdisease, but such systematic screening remains an expensive and time-consuming process. It takes more than 13 years and $2.87 billion to bring a new drug into the market [Hauser et al., 2017]. Researchers are thus resorting to machine learning to understand the drugs' MoAs by mining the emergence of large-scale chemical and genomic data [Chen and Li, 2018].

Two most prominent data-driven tasks among these recent developments are drug-disease and drug-target prediction (see [Ezzat et al., 2018; Li et al., 2015] for surveys). In these tasks, researchers have attempted to collect a variety of omics data, and predict new interactions of drug-disease or drugtarget through network inference [Chen and Li, 2017], multiview learning [Zheng et al., 2013], and deep learning [Tsubaki et al., 2018]. Beyond pairwise drug-disease or drug-target relationships, some recent studies have pointed at the importance of identifying triple-wise interactions of drug-targetdisease in human metabolic systems [Capuzzi et al., 2018; Chen and Li, 2019; Wang et al., 2018]. Among different methods, tensor factorization is a commonly used method to infer the missing entries of a drug-target-disease tensor [Chen and Li, 2019; Wang et al., 2018].

Tensor Factorization. Tensor factorizations aim to extract latent structure from high dimensional data [Kolda and Bader, 2009]. CP (CANDECOMP/PARAFAC) and Tucker are two popular tensor models with diverse variants being successfully applied in health data analysis [Wang et al., 2015; Chen and Li, 2019]. However, the CP and Tucker models (or their variants) suffer from two weaknesses.

First, their performance can be limited by linearity, which might not be expressive well for nonlinear data manifolds. Recently, a series of studies have shown that nonlinear tensor factorizations have superior performances over multilinear tensor models [Xu et al., 2012; Liu et al., 2019; Wu et al., 2019]. For example, InfTucker [Xu et al., 2012] was proposed to use a nonlinear Gaussian kernel. However, they rely on a prior Gaussian process over tensor data, which might be difficult to estimate in practice [Zhe et al., 2016].

The second drawback of CP or Tucker models is the data sparsity issue. To alleviate this issue, coupled tensor-matrix models are extended to jointly analyze tensor together with auxiliary information [Wang et al., 2015; Narita et al., 2012; Acar et al., 2011]. However, these methods are inherently limited by encoding feature matrices, which require tedious feature engineering [Shan et al., 2016]. As the number of features grows, designing and deploying them become challenging, especially for healthcare data.

Contributions. To tackle these challenges, we propose a Neural Tensor Network (NeurTN), which combines tensor algebra and deep neural networks to provide effective medicine recommendations. By replacing the multilinear multiplication with a neural network, NeurTN is able to characterize nonlinear dependencies among tensor data. Moreover, NeurTN incorporates various heterogeneous information to alleviate the data sparsity issue. Instead of constructing them as feature matrices, NeurTN uses geometric neural networks to learn the embeddings frommolecular graphs and target sequences, which allows to be trained end-to-end. Our datadriven model opens up opportunities to use large-scale omics data to discover drug' MoAs in pharmacological studies.

2 Preliminaries

2.1 Tensor Algebra

Tensors are multidimensional arrays that extend the concept of matrices [Kolda and Bader, 2009]. An N-order tensor is denoted as X ∈ RI1×I2×...×IN , where the (i1, . . . , iN)- th element of X is denoted as Xi1,...,iN . Matricization is the process of flattening a tensor into a matrix. The mode-n matricization of tensorX is represented as X(n) ∈ RIn×I1...In-1×In+1...IN and is arranging the mode-n fibers of the tensor as columns in X(n).

N-mode product. The n-mode product of a tensor X ∈ RI1×I2×...×IN with a matrix A ∈ RJ×In is denoted as X ×nA with size RI1···×In-1×J×In+1···×IN . Also, we have:

CP tensor factorization. The CP is the most common tensor model, which assumes a compact hidden structure in the data [Kolda and Bader, 2009]. Formally, the CP model approximates the element of a third-order X as:

where U,V,W are the factor matrices; r denotes the tensor rank. The CP model has its ease of interpretation: each rankone component serves as a latent concept or clusters for the data. However, the CP model adopts multilinear assumption, which may be insufficient to capture more complex and nonlinear feature interactions. Our work here addresses this issue by designing the predictor f(·) using neural networks.

2.2 ProblemDefinition

The medicine recommendation task can be formulated as a tensor completion problem [Chen and Li, 2019; Wang et al., 2018]. To be specific, the input can be organized as a drugtarget-disease tensor X ∈ RM×N×L, where M, N, and L denote the number of drugs, targets, and diseases, respectively. An entry Xijk = 1 if an interaction among drug i, target j, and disease k is observed; Xijk = 0, otherwise. Our aim is to estimate the scores of unobserved elements ˆXijk, which can be used to infer novel interactions of drug-target-disease.

2.3 Input and Embedding Layer

Our goal is to estimate a score between drug i, target j and disease k, i.e., ˆXijk. We first describe the biological features with respect to drugs, targets, and diseases.

One-hot Embeddings

Given a drug i, a target j, and a disease k, their one-hot encodings ai ∈ RM , bj ∈ RN , and ck ∈ RL can be obtained based on their identities. We can obtain their dense embeddings via three lookup tables:

where ̂ui ∈ Rd1 , ̂vj ∈ Rd2 , and ̂wk ∈ Rd3 are new embeddings for drug i, target j, and disease k, respectively.

Here we also incorporate medical knowledge of drugs and targets. We leave the exploration of diseases as a future work. Instead of constructing these auxiliary information as feature matrices [Zheng et al., 2013; Chen and Li, 2019], we apply geometric neural networks to learn the structural information of both drugs' chemical structures and targets' sequences.

Graph Neural Network for Molecular Graph

The chemical structure of a drug determines its physicochemical properties, further affects its pharmacological activity. In this study, we extract features of chemical structures by using Graph Neural Networks (GNN) [Duvenaud et al., 2015; Gilmer et al., 2017]. A chemical molecule can be represented as a graph G = (V, E), where nodes denote the atoms and edges denote the chemical bonds between atoms. The key idea behind GNN is to iteratively aggregate information from atoms' local neighbors with convolutional operator. Such information can be further propagated across the graph by stacking multiple graph convolutional layers, resulting in a learned featurization for each molecule.

Formally, given a graph Gi for drug i, the GNN converts the graph Gi into a vector representation riG by aggregating the representation rv of each atom node v in the graph [Duvenaud et al., 2015]. GNN first obtains the representation r(t) v by the following transition function:

where r(t) v is the representation of node v in the t-th step and r(0) v can be initialized by its atom features [Duvenaud et al., 2015]. N(v) is the neighbors of v and |N(v)| denotes the number of neighbors around v. H|N(v)| (t) is a weight and σ(·) is the sigmoid. A global pooling function then acquires the graph representation by combining features from all the atoms as:

here W(t) gnn is a weight and riG ∈ Rg is the dense vector representations for the graph Gi.

Convolutional Neural Network for Target Sequence

Similarly, we can obtain vector representations of target sequences by using a convolutional neural network (CNN). Here we adopt a recent developed CNN model to map target sequences to a real-valued vectors by applying hierarchical filter functions [Tsubaki et al., 2018]. Given that a target sequence can be treated as a string of amino acids, it can be transformed into a list of words by using n-gram feature extraction technique. Let Sj = x1, x2, . . . , x|Sj| be a sequence for target j, where xi is the i-th word after transformation and |Sj| is the length of sequence. All words can be then translated into randomly word embeddings as:

where xi ∈ Rd is the embedding of the i-th word. Let input c(0) i = xi:i+w-1 refer to [xi, . . . , xi+w-1], the concatenation of w contiguous embeddings. The CNN uses one filter as:

where c(1) i ∈ Rd is a hidden vector, Wcnn ∈ Rd×dw is the weight and bcnn is the bias. Note that the filter projects a dw-dimensional input vector into a d-dimensional hidden vector. As such, one can obtain a set of hidden vector: C = {c(t) 1 , c(t) 2 , . . . , c(t) C } by applying the filters hierarchically [Tsubaki et al., 2018]. Finally, an average function is used to aggregate features from hidden vectors as:

here rjS ∈ Rd is the representation for the target sequence Sj.

The Unified Embeddings

We fuse all available features for more sophisticated embeddings by using fully connected neural layers. Given drug's features (̂ui, riG), target's features (̂vj, rjS ), and disease's feature ̂wk, we can obtain new embeddings ui, vj , and wk 1 via:

where ⊕ is the concatenation. In addition, by choosing the parameters Θu and Θv properly, we can reshape all embeddings with the same size, i.e., {ui, vj,wk} ∈ Rr . As such, a nonlinear tensor model can be built upon these embeddings.

3 The ProposedModel

The proposed NeurTN contains three components: MLP, GTF, and CTN as shown in Figure 1.

3.1 Multi-Layer Perceptron

It is straightforward to feed all the features into a MultiLayer Perceptron (MLP) [He et al., 2017; Wu et al., 2019], in which each hidden layer can learn nonlinear feature interactions from ui, vj , and wk:

where WL and bL denote the weight and bias. We choose the ReLU(·) as activation function in hidden layers and the sigmoid function σ(·) as predictor function.

The MLP is capable of learning nonlinearity of the concatenated features. Nevertheless, the concatenated features, i.e., ui ⊕ vj ⊕wk, may lose some information in the original embeddings that are useful for later interaction learning. To avoid such information loss, we further propose two triplewise layers for learning multi-aspect features.

1wk only contains one-hot feature since there is no auxiliary information for diseases in this study.

3.2 Generalized Tensor Factorization

Inspired by the success of tensor CP model in Eq. (1), we propose a Generalized Tensor Factorization (GTF) to learn nonlinear triple-wise feature interactions. Specifically, given the drug embedding ui ∈ Rr , target embedding vj ∈ Rr , and disease embedding wk ∈ Rr , we define a novel triple-wise interaction layer as:

here ⊙ denotes the element-wise product. Then, we can project the hidden vector φ(ui, vj,wk) into an output layer:

here fout(·) and h denote the activation function and weights. Relation to tensor CP model: Intuitively, if we use an identity function for fout (e.g., fout(x) = x) and enforce the weights h to be a uniform vector of 1, we can rewrite the Eq. (8) as ˆXijk = ∑rt=1 uitvjtwkt, which exactly recovers the tensor CP model as Eq. (1). In other words, GTF extends the CP model into a nonlinear tensor machine. As neural networks have strong ability to fit the data, GTF is thus more expressive than the linear CP model. In this work, we implement the GTF layer by using the sigmoid function fout(x) = 1/(1 + e-x ) and learning h from the data.

3.3 Compressed Tensor Network

To better capture the triple-wise feature interactions among drug-target-disease, we further consider the outer product on the embeddings ui, vj , and wk:

where E ∈ Rr×r×r is a tensor feature map and ⊗ denotes the outer product. Our motivation for E is straightforward. The map E captures more signals than elementwise product in Eq. (7) since it encodes any triple-wise feature interactions. Such strategy has been widely used to boost system performance in deep learning [He et al., 2018; Cohen et al., 2016]. Moreover, the feature map E lies in a 3D tensor format, which allows us to design an efficient feedforward neural network by using tensor algebra latter.

The 3D tensor format of E naturally suggests a CNN model to use shared parameters to perform "convolution". However, the feature E from the outer product space does not have spatial locality. Thus, we cannot employ the CNN to E by enforcing an unexpected local connectivity pattern. Inspired by the tensor algebra in n-mode product (Sec 2.1), we propose a simple but efficient Compressed Tensor Network (CTN) to perform feedforward computations. Basically, the tensor map E is successively compressed by a sequence of three weight matrices along each mode with nonlinear transformation:

where Hl, {A(l) ,B(l) ,C(l) }, and B(l) denote the hidden tensor, the weights, and bias in the l-th layer, respectively. Reshape(·) flattens HL into a vector. The output layer is a fully connected layer with the sigmoid function as predictor.

Dropout: Dropout [Srivastava et al., 2014] is a regularization technique for neural networks to prevent overfitting. The idea of dropout is to randomly turn off neurons with probability ρ during training, and use all neurons in the testing. Here we also apply a dropout layer on the feature E, i.e., randomly dropping ρ percent of its elements.

3.4 The Overall Model

Joint Training. We have developed three instantiations of nonlinear tensor models: MLP, GTF, and CTN. We present our unified model NeurTN by joint learning these three modules. Formally, let zL, φ, and HL denote the outputs of the last hidden layers of MLP, GTF, and CTN, respectively. Then, we have:

here we choose the sigmoid function σ(·) as final predictor.

Relation to wide&deep learning: The proposed NeurTN is similar to the well-known Wide&Deep Learning [Cheng et al., 2016; He et al., 2017], which has the benefits of memorization and generalization. Our GTF can be regarded as a wide component whereas the MLP and CTN can be viewed as deep components.

Model Optimization. We adopt pairwise learning methods to optimize model parameters [Bordes et al., 2013]. The idea behind pairwise learning is that an observed triplet should be predicted with a higher score than an unobserved one. This can be achieved by minimizing:

here f (·) and Θ denote our predictive function and model parameters in Eq. (11). D+ denotes the set of positive triplets

(e.g., Xijk = 1), and D- denotes the set of negative triplets corresponding to D+ by sampling from unobserved elements. Following [Bordes et al., 2013], for each positive training triplet (i, j, k), we randomly sample one negative training triplet (i′ , j ′ , k′ ) in the training step.

It is worth mentioning that our model is highly instructive and allows the use of deep learning techniques, such as dropout regularization, which is more effective than the kernel-based nonlinear tensor models (see Sec 4.2).

4 Experiments

4.1 Datasets and Baselines

Datasets. We obtain data from three public databases [Chen and Li, 2019; Wang et al., 2018]: CTD 2 , DrugBank 3 , and UniProt 4 . We only focus on drugs that have DrugBank identifier for later collecting auxiliary information. As such, we obtain 436, 322 triplets of drug-target-disease, involving 1, 901 drugs, 2, 514 targets, and 2, 923 diseases. For drugs, their SMILES, a string encoding of chemical structures, are downloaded from DrugBank. These SMILES strings can be converted to molecular graphs using RDKit tool 5 , which can be then fed into the GNN module. For targets, their amino acid sequences are collected from UniProt and can be used by the CNN module without any pre-processing.

Baselines. We mainly compare with tensor-based models that can learn the drug-target-disease data. (1) CP and Tucker [Kolda and Bader, 2009]: both adopt multilinear assumption. (2) nTucker [Zhe et al., 2016]: a nonlinear Tucker based on Gaussian process. (3) CoSTCo [Liu et al., 2019]: a recent CNN-based tensor model. (4) CMTF [Acar et al., 2011]: a tensor-matrix model regarding auxiliary information as feature matrices. (5) TFAI [Narita et al., 2012]: a tensor model with mode regularization. (6) AirCP [Ge et al., 2016]: a tensor-matrix model with graph regularization. (7) NTF [Wu et al., 2019]: a neural network with MLP. (8) Rubik [Wang et al., 2015]: a tensor model with non-negativity and sparsity constraints. (9) DTD [Chen and Li, 2019]: a recent tensor-matrix model in drug discovery.

Parameter Settings. For tensor-matrix models CMTF, TEAI, AirCP, Rubik, and DTD, the feature matrices for auxiliary information are constructed using feature engineering as [Zheng et al., 2013; Chen and Li, 2017], The parameter settings for all the baselines are carefully tuned to achieve optimal performance. For NeurTN, the embedding size r in Eq. (5) is searched within [16, 32, 64, 128]. For MLP and CTN, we both employ three hidden layers with dropout ratio ρ = 0.3 and each layer sequentially decreases the half size of inputs. Our models are built upon PyTorch 6 with Adam optimizer [Kingma and Ba, 2015]. We search the batch size and the learning rate within {128, 256, 512, 1024} and {0.001, 0.005, 0.01, 0.05, 0.1}, respectively. We use gridbased search to find the best parameter settings. We tune model parameters using validation set and terminate training if the performance does not improve for 100 epochs.

2 http://ctdbase.org/downloads/

3 https://www.drugbank.ca/

4 https://www.uniprot.org/

5 http://rdkit.org/

6 https://pytorch.org/

Evaluation Protocols. We randomly split the dataset into 80% training, 10% validation, and 10% test sets. The validation set is used for tuning hyper-parameters and the final results are evaluated on the test set. To better construct negative test triplets, we adopt similar procedures as [Ge et al., 2016; Chen and Li, 2019; Bordes et al., 2013]: 1) Scenario 1 (random sample): for each positive test triplet (i, j, k), we randomly sample a negative triplet (i′ , j ′ , k′ ) such that (i′ , j ′ , k′ ) is unobserved from data, i.e., Xi′j′k′ = 0; 2) Scenario 2 (sample drug mode): we corrupt the drug mode by replacing the drug i with a new drug i′ so that the (i′ , j, k) is unobserved; 3) Scenario 3 (sample target mode): we corrupt the triplet (i, j, k) by (i, j ′ , k) in the target mode; 4) Scenario 4 (sample disease mode): the triplet (i, j, k) is replaced with (i, j, k′ ). Such scenarios are meaningful in clinics. For example, given a pair of drug-target (i, j), we can find their potential indications for new disease k′ [Hauser et al., 2017]. In addition, we apply filtered settings [Bordes et al., 2013] such that those test negative samples will not appear in the training step.

We use two common top-n metrics, Hit@n and NDCG@n, to evaluate recommendation performance [He et al., 2017; Bordes et al., 2013]. The metric of Hit@n measures if a test triplet is among the top-n ranked list, while NDCG@n is a position-aware metric which assigns larger weights on higher positions. Moreover, we use the strategy in [Bordes et al., 2013] to avoid heavy computation on all triplets. For example, in Scenario 1, we randomly generate 100 negative samples (i′ , j ′ , k′ ) for each test (i, j, k). Based on the ranking of these triplets, Hit@n and NDCG@n can be computed.

4.2 Experimental Results

Pure Tensor Completion. We have proposed four variants of nonlinear tensor models: MLP (Eq. (6)), GTF (Eq. (8)), CTN (Eq. (10)), and the unified NeurTN (Eq. (11)). In this section, We first compare them with the baselines CP, Tucker, nTucker, and CoSTCo, which are pure tensor machines without any auxiliary data of drugs or targets. For fair comparison, we only use the one-hot features in Eq. (2) as inputs for the proposed models. Due to page limitation, we only show the top-n performances for Scenario 1, and similar trends can be observed under different scenarios. Table 1 shows the results in terms of Hit@n and NDCG@n.

It is observed that nonlinear tensor models consistently outperform the multilinear models. For example, nonlinear nTucker performs better than Tucker; GTF outperforms the CP model and gains average improvements of 7.39% on Hit@n and 11.63% on NDCG@n. These improvements are statistically significant and mainly stem from the powerful representation of nonlinear models. Second, the performance of CoSTCo is limited for drug-target-disease data that do not have local connectivity. As such, CoSTCo may be insufficient to capture nonlinear patterns for biological data using a CNN. Third, the performances of CTN are better than GTF, implying a good representations of a deeper neural network. Finally, NeurTN achieves the best performances, presumably this owes to the the benefits of its joint learning.

Tensor CompletionwithAuxiliary Information. Nowwe compare the overall performance of NeurTN with the baselines with auxiliary information. The results of CP, Tucker, nTucker, and CoSTCo are omitted due to their inferior performances without auxiliary data. Figure 2 shows the performance of Top-n performances, where n = [5, 10, 15].

From the figures, we can observe that our proposed NeurTN outperforms coupled tensor-matrix factorization methods in all scenarios. The superior performance of NeurTN mainly benefits from its deep neural networks. It is intuitive that neural networks would have stronger ability to fit the data, while the multilinear assumptions in coupled tensor-matrix factorizations do not. More importantly, NeurTN utilizes geometric neural networks GNN and CNN, which can learn features from molecular graphs and protein sequences in the training process. Such end-to-end representation learning can potentially obtain more interpretable data-driven features instead of predefining hand-crafted feature matrices in coupled tensor-matrix factorizations. Given these encouraging results, we use our model to predict novel interactions in the diabetes-specific pathway.

Case Study on Diabetes. In clinics, given a disease, it is critical to know which drugs can treat this disease as well as the targets involving in the disease-specific pathway. The pairs of (drug, target) related to a special disease are important in personalized treatments. Here we use NeurTN to predict the novel (drug, target) pairs related to diabetes's metabolic pathway. The top-10 novel pairs of (drug, target) corresponding to diabetes are listed in Table 2. The predictive results can be then evaluated by domain experts to see whether such recommendations are clinically meaningful. They can be further validated in vivo experiments, such as animal studies and clinical trials. In summary, the proposed model provides a systematic approach to narrow down the search space for further clinical trials in drug discovery.

4.3 Ablation Study

We perform some ablation studies to investigate the impact of each module in NeurTN. Table 3 shows the performance of our default method and its variants in the case of Scenario 1.

Our results are summarized as follows: (1) Remove MLP: Without MLP layers, we find that the performance is slightly worse. Although the black-box nature of MLP, its hierarchical structure is still helpful to learn more complex interactions; (2) Remove GTF: Not surprisingly, the results are worse than the default method, which suggests that the GTF can capture triple-wise feature interactions in a nonlinear fashion; (3) Remove CTN: This variant substantially decreases the overall performance with a large margin, verifying the effectiveness the CTN in capturing the useful feature interactions from outer product feature space; (4) Remove GNN: The chemical structure of a drug determines its pharmacological activity. Removing GNN thus decreases the overall performance; (5) Remove CNN: Similar, amino acid sequences determine the therapeutic function of targets. Deleting the CNN module thus hurts the final performance.

Embedding size of NeurTN: The embedding size r in Eq. (5) affects the representation ability of NeurTN. We vary r within [16, 32, 64, 128]. As shown in Figure 3(a), NeurTN benefits from a large embedding size in Scenario 1. Results on other scenarios have similar trends and are omitted.

Layers of NeurTN: We also conduct experiments to investigate the impact of the number of layers in NeurTN. We vary the number of layers in the NeurTN within L = [2, 3, 4, 5].

As shown in Figure 3(b), stacking more layers generally improves the system performance, which demonstrates the usage of more neural layers is able to model complex drugtarget-disease interactions.

Dropout regularization: Figure 3(c) shows the performances of NeurTN w.r.t. dropout ratio. Our results show that dropout offers better performance. Specifically, using a dropout ratio ρ ≈ 0.3 achieves an optimal accuracy.

5 Conclusion

A critical step in personalized medicine is to understand drugs' MoAs by exploring the biological interactions among drugs, targets, and diseases in human metabolic systems. In this study, we present a novel NeurTN, which seamlessly combines the tensor algebra and deep neural network to capture the nonlinear relationships among health data. In the future, we aim to incorporate the auxiliary information of diseases under the NeurTN.

Acknowledgments

This work has been supported in part by NSF CCF1815139.