Deep Learning Can Identify Explainable Reasoning Paths of Mechanism of Drug Action for Drug Repurposing from Multilayer Biological Network


 The discovery and repurposing of drugs require a deep understanding of the mechanism of drug action (MODA). Existing computational methods model these mechanisms with the help of protein-protein interaction (PPI) network. However, the molecular interactions of drugs in the human body are far beyond the PPI network, and the lack of interpretability of these models hinders their practical applications. In this study, we propose iDPath, an interpretable deep learning-based path-reasoning framework to identify potential drugs for the treatment of diseases by capturing the MODA based on simulating the paths from drugs to diseases in the human body. We create a multilayer biological network consisting of a gene regulatory layer, a PPI layer, a protein-chemical interaction layer, and a chemical-chemical interaction layer to comprehensively characterize the molecular interactions in the human body. Based on this multilayer biological network, iDPath utilizes a graph convolutional network module to capture the global connectivity information of the human molecular network; and a long short-term memory neural network module to capture the detailed mechanisms of drug action by modeling the sequential dependencies of the MODA-related biological paths between drugs and diseases. Moreover, iDPath applies two attention modules - node attention and path attention - to further enhance its interpretability. Experiments with drug screen data show that iDPath outperforms state-of-the-art machine learning methods on a general drug repurposing task featuring 1,993 drugs and 2,794 diseases. Further investigations demonstrate that iDPath can identify explicit critical paths that are consistent with clinical evidence. To demonstrate the practical value of iDPath, we apply it to identify potential drugs for the treatment of prostate cancer. Results show that iDPath can successfully discover new FDA-approved drugs. This research provides a novel interpretable artificial intelligence perspective on drug discovery.


Introduction
Artificial intelligence (AI) has recently shown the huge potential to subvert the typical drug discovery process 1 . Scientists are now using deep learning technologies to discover candidate drugs for the treatment of COVID-19 2-4 , Alzheimer's disease 5 , cancers 6 , and so on. Among various applications, drug repurposing, the identification of the new use of approved or investigational drugs that are outside of the original medical indication, can shorten the time of drug development while ensuring safety, and thus attracts special attention from drug discovery communities and industries 7,8 . Existing computational approaches aim to precisely model the effects of drugs in the human body. A common approach is to model the outcomes of drugs in a clinical environment, such as using electronic health records to discover the efficacy of drugs on a specific population 9 , or emulating clinical trials on real-world patient data 10 .
With the development of high-throughput omics technologies, the detailed characterization of the molecular interactions of drugs in the human body becomes possible 11 . The protein-protein interaction (PPI) network serves as a "skeleton" for body's signaling circuitry 12 and shows tremendous power in guiding drug discovery [13][14][15][16] . A series of studies explored the potential of mining the network properties of drugs in the PPI network in synergistic drug combination identification 13 and drug repurposing 14 . A recent study introduced advanced graph-based deep learning approaches to the identification of anti-cancer drug combinations by learning the graphic representations of the PPI network 17 . However, the molecular network in the human body is not limited to the PPI network, the gene regulatory mechanisms 18 , the binding work of the proteins and chemicals 19 , and the interactions of the chemicals 20, 21 also play a role in the mechanism of drug action (MODA). These processes rely on the drug's interactions with various proteins and chemicals in the human body 22,23 . Usually, the MODA is described by biological pathways, a series of biochemical and molecular steps to achieve a specific function or to produce a certain 2/26 product 24 . Such biological pathways can be naturally denoted as a series of paths in the biological network.
Furthermore, instead of targeting specific proteins, some drugs need to take further chemical reactions to be effective 25 . For example, cytarabine, an important drug in the treatment of acute myeloid leukemia 26 , must be phosphorylated intracellularly to a nucleotide (cytarabine 5'-triphosphate, Ara-CTP) before it can exert its cytotoxic effect 27 .
Recent studies have introduced novel machine learning models to drug discovery based on human biological networks [28][29][30][31] . Most of these networks are combinations of PPI network and other networking information such as the drug interactions and drug similarities 30,31 . The rich patterns hidden in the molecular interactions (such as chemical-chemical interactions) in the human body have not been fully used in drug discovery and repurposing tasks. Whether the addition of these molecular interactions can facilitate the accurate modeling of drug effects is under-researched. Furthermore, nearly all of these machine learning models are black-box. Users cannot explain why the model generates a certain prediction.
The lack of interpretability hinders machine learning models' potential in practical drug discovery tasks.
The need for explainable machine learning models led to the development of a series of novel neural network architectures, such as attribution methods 32,33 and knowledge-graph-based models 34,35 . These models have been further applied in healthcare [36][37][38][39] , such as using biological-informed neural networks to identify anti-cancer drug combinations 38 and predict disease prognostic based on comorbidity network 39 .
However, whether these explainable modules can accelerate drug discovery and further enhance the knowledge of the MODA is unknown.
To fill the aforementioned research gaps, based on our previous research on interpretable machine learning models 4,17,39,40 , we propose the interpretable Deep learning-based Path-reasoning framework for drug repurposing (iDPath), which captures the mechanisms of drug actions by identifying the paths from drugs to diseases in the human body. To accurately characterize the MODA, we build a multilayer 3/26 biological network instead of using the PPI network alone. The multilayer biological network is the integration of a gene regulatory layer, a protein-protein interaction layer, a protein-chemical interaction layer, and a chemical-chemical interaction layer. Starting from this multilayer biological network, iDPath utilizes a graph convolutional network module to capture the global connectivity information of the human molecular network; and a long-short term memory neural network module to capture the detailed mechanisms of drug action based on the shortest paths between drugs and diseases. Furthermore, iDPath introduces two attention modules, namely the path attention and the node attention, to enhance model interpretability. Experiments with drug screen data demonstrate the superior performance of iDPath in a general drug repurposing task featuring 1,993 drugs and 2,794 diseases. Further investigations demonstrate that iDPath can identify explicit critical paths that are consistent with clinical evidence. To demonstrate the practical value of iDPath, we apply it to identify potential drugs for the treatment of prostate cancer.
Results show that iDPath can successfully discover new FDA-approved drugs. These results indicate that iDPath can facilitate drug discovery and repurposing, and has the potential to address other computational chemistry and biology tasks involving the understanding of the molecular interactions in the human body.

Methods
In this section, we describe the datasets and the proposed iDPath model, as well as baseline models for drug repurposing, including DeepWalk, Graph Convolutional Networks (GCN), Long Short-Term Memory (LSTM) networks, and Knowledge-aware Path Recurrent Network (KPRN).

Multilayer Biological Network
Gene regulatory network (GRN) layer. The gene regulatory network is adopted from RegNetwork 41 , which collects the experimentally validated and the predicted regulations based on the transcription factor 4/26 binding sites. The edges in RegNetwork start from transcription factor (TF) and microRNA (miRNA) and target the regulated genes. In total, RegNetwork provides us with 369,277 gene regulations between 1,456 TFs, 1,904 miRNAs and 19,719 genes.
Protein-protein interaction (PPI) layer. The PPI network consists of information from two sources.
The first dataset, STRING dataset 42 , is the most comprehensive database of known and predicted PPIs till now, with more than 1,380 million PPIs among over 9 million proteins. We only keep the PPIs in the human body (homo sapiens) and at high confidence or better (confidence level > 0.7). Another PPI dataset is the human interactome built by Cheng et al. 14 . This dataset is harnessed from multiple databases with experimental evidence. After pre-processing, our PPI network contains 614,970 interactions connected by 13,758 proteins.

Protein-chemical interaction (PCI) layer.
We obtain a PCI network by curating from STITCH database 43 , which is the most comprehensive database of known and predicted interactions between chemicals and proteins till now. We select PCIs in the human body (homo sapiens) at high confidence or better (confidence level > 0.7). The processed PCI network consists of 203,551 interactions among 9,393 proteins and 73,199 chemicals.
Chemical-chemical interaction (CCI) layer. CCI network is curated from STITCH 43 database and further processed by selecting CCIs in the human body (homo sapiens) at high confidence or better (confidence level > 0.7). The processed CCI network has 396,284 interactions among 107,055 chemicals.
Constructing the multilayer biological network. We construct a multilayer biological network by mapping all the entities in GRN, PPI, PCI and CCI to same nomenclature. The proteins are named by their encoded genes and the miRNAs are mapped to their corresponding genes by BioMart 44 . All the genes are encoded to their Entrez IDs 45 . All the chemicals are denoted by their PubChem CIDs (Compound ID 5/26 number) 46 .

Therapeutic drug-disease pairs
For the drug repurposing task, we collect therapeutic drug-disease pairs from Therapeutic Target Database (TTD) 47 , which provides the known and explored therapeutic protein and nucleic acid targets, the targeted disease, and the pathway information of tens of thousands of drugs. We only keep the FDA-approved drugs in TTD and map them to PubChem CID in order to be consistent with the chemicals in our multilayer biological network. The diseases in TTD are in the ICD-11 coding system and are mapped to their corresponding ICD-10 codes. The cleaned dataset of therapeutic drug-disease pairs includes 1,993 drugs and 2,794 diseases and constitutes 19,500 pairs.

Drug-protein associations and drug-chemical associations
We collect drug-protein associations from 4 datasets: the PCI network from STITCH 43

Disease-gene associations and disease-miRNA associations
The disease-gene associations include genes and variants associated with human diseases, curated from DisGeNET 49 by selecting expert-curated repositories. The miRNAs associated with human diseases come from the Human microRNA Disease Database (HMDD) 50 , which is a database about curated experiment-

6/26
supported evidence for human microRNA (miRNA) and disease associations. All the genes, variants, and miRNAs are mapped to Entrez IDs, and diseases are mapped to ICD-10 codes. After processing, we have 230,837 associations among 7,559 genes, 6,830 variants, 705 miRNAs with 5,602 diseases.

Overall architecture of iDPath
The iDPath framework for drug repurposing is presented in Figure 1

GCN to capture the global connectivity information of the multilayer biological network
With the uniform nomenclature of the nodes, we aggregate the drug-related information (drug-protein associations, drug-chemical associations), the disease-related information (disease-gene associations and disease-miRNA associations), and the multilayer biological network to the network G = (V, E), where V and E are the node set and edge set, respectively, iDPath introduces a three-layer GCN, following a spatial-based GCN architecture 51 , to encode the global connectivity topological information of the multilayer biological network. Suppose there are N nodes in total, the initial embeddings of these nodes are E ∈ R N×d (d is the dimension of the embedding), and the adjacency matrix of network G is A, the  Figure 1. The framework of iDPath on drug repurposing tasks. a The multilayer biological network consists of four layers: gene regulatory network layer (one-way red arrows), protein-protein interaction layer (two-way black arrows), protein-chemical interaction layer (two-way purple arrows), and chemical-chemical interaction layer (two-way orange arrows). The blue two-way dashed arrows represent that the two corresponding nodes in different layers are identical. The nodes associated with the drugs and diseases are marked by green dashed lines. b The MODA-related biological paths are identified by the shortest paths between drug and disease generated in the multilayer biological network. c The schematic representation of the algorithm: the multilayer biological network is fed into three-layer GCN to learn the embeddings of all nodes. The GCN embeddings of nodes along the shortest path between one drug and one disease are fed into an LSTM module to learn the sequential dependencies. Node attention and path attention modules are introduced to aggregate the embeddings of the nodes and paths. The final prediction is the probability that one drug is effective for one disease.
computation formula of layer l of GCN is shown as below: where D is the diagonal node degree matrix, I is the identity matrix.

MODA-related biological paths
The MODA is dependent on the interactions of drugs with molecules in the human body, which can be represented as a series of paths in our multilayer biological network 24 . To accurately model the effects of drugs, we need to identify informative paths to represent the MODA in an efficient way.

8/26
Among various approaches to generate these paths, the shortest paths are the most informative ones, which have already been applied in many drug-related studies 13,14,52 . Considering the huge amount of nodes and edges of our multilayer biological network, finding the shortest paths between all drugdisease pairs would be time-consuming. To overcome this challenge, we introduced GPU-accelerated graph algorithms implemented by cuGraph 53 . For a drug and a disease, the shortest path between them is connected by their associated nodes in the multilayer biological network. Given a drug and a disease, the shortest paths between this pair form a set PAT H = {path 1 , path 2 , ..., path L }, where .. → node disease }, node m1 and node m2 denote the middle nodes of one path, and L is the number of shortest paths. Since L differs among drug-disease pairs, we choose a fixed value for L by randomly sampling from the shortest paths set PAT H.

LSTM layer
Given a drug-disease pair, the embeddings E GCN generated by the GCN, and the shortest path set between this pair PAT H, we employ LSTM 54 to encode both long-term and short-term dependencies in a MODA-related biological path. Such sequential dependencies are crucial to the model intelligibility.
Meanwhile, we introduce the type of these nodes to strengthen the model's capability of identifying different nodes. Here, we consider four types: protein (gene), chemical, drug, and disease. Their embeddings E TY PE ∈ R 4×d are randomly initialized. Therefore, given one node node j of one path path p , the input to LSTM is the concatenation of its GCN embedding e j = E GCN [node j ] and type embedding e ′ j = E TY PE [node j ], that is: where ⊕ denotes the concatenation operation. The hidden state h j−1 generated by the previous node node j−1 in the same path and x j are used to learn the hidden state of the input of the next node node j+1 ,

9/26
which is defined as follows: where i j , f j , g j , and o j are the input, forget, cell, and output gates, respectively; c j ∈ R d ′ and h j ∈ R d ′ are the cell state and hidden state at path step j, and d ′ is the dimension of the output; and W h are learnable weights, and b i , b f , b g , and b o are bias; σ denotes the sigmoid function and ⊙ is the Hadamard product. The hidden states h j for each path step are aggregated to the attention modules for the representation of the whole path and final prediction. Since the shortest paths are not of equal length, we borrow the padding method as follows 55 . Suppose the maximum length of one path is set to l max , for the paths that are shorter than l max , we use a padding value pad (such as 0) to fill the path, and the following processing will ignore these padding positions not to affect the performance.

Node attention and path attention
Attention mechanism 56 is widely used in various deep learning tasks to enhance model intelligibility 17  infinity for the following Softmax transformation. Then, we applied a linear layer with the Softmax activation to aggregate the embeddings to one numeric value denoting the importance (weight). That is where W n ∈ R d ′ ×1 is the learnable parameter, Ω p ∈ R l max ×1 denotes the weights of each node in path path p , where Ω p = {ω 1 , ω 2 , ..., ω l max }. For one node j in path p , its weight is computed as follows: Then, we aggregate the hidden states of these nodes weighted byŵ i to get the embedding of path p : Path attention. After the aggregation operation, for one drug-disease pair, the embeddings of their shortest paths are: E PAT H = {e path 1 , e path 2 , ...e path L }, the path attention layer is similar to the node attention: where σ p is the sigmoid function. The final predictionŷ is the aggregation of the weighted embedding of each path, indicating the probability that the drug is effective in treating the disease.

Objective function
In this study, we treat the training of iDPath as a binary classification task following common practice.
That is, besides all the therapeutic drug-disease pairs (marked as 1), we introduce the negative sampling to get an equal number of non-therapeutic drug-disease pairs (marked as 0). The objective function used by iDPath is the binary-cross-entropy loss with l 2 regularization: where Θ is the set of parameters to be learnt in iDPath, λ is the l 2 regularizer to prevent over-fitting.

Baselines
In this section, we describe a set of baseline models to compare with, including DeepWalk, GCN, LSTM, and Knowledge aware Path Recurrent Network (KPRN).
DeepWalk. DeepWalk 57 is a widely used graph embedding approach by modeling a stream of short random walks and has already been introduced to several drug-related tasks, such as drug-target identification 58 .
GCN. GCN 51 is widely applied to many drug-related tasks, such as drug discovery using the drug's Smiles features 59 , and anti-cancer drug combination identification 17 . As a component of iDPath, we apply GCN individually to test its performance on the drug repurposing task. Specifically, besides the basic graph convolutional layer, we introduce a fully connected layer to combine the embeddings of drug and disease for the final prediction.
LSTM. LSTM has been applied to drug discovery 60,61 . We apply a vanilla LSTM network and use the last hidden states of each path as its representation. A two-layer fully connected layer followed by the

12/26
LSTM layer is employed to generate the final prediction.
KPRN. Knowledge aware Path Recurrent Network (KPRN) 35 is an advanced path-based model for a reasonable recommendation based on a knowledge graph. We apply KPRN to the drug repurposing task by feeding it the same input as iDPath.

Performance evaluation
All models are trained on a binary classification task, where we utilize commonly used metrics to evaluate and fine-tune the models. The metrics used in the binary classification task include accuracy, recall, the area under the receiver operating characteristic curve (AUROC), and the area under the precision-recall curve (AUPRC). The drug repurposing task can be viewed as a recommendation task, that is, for each disease in our dataset, we use the model to calculate the probability that one drug is effective in treating one disease and then rank all the drugs based on the probability. We introduce two commonly used metrics in the recommendation system, NDCG@K and Hit@K. Details of these metrics and the experiment setup can be found in the supplementary information.

Results and Discussions iDPath consistently outperforms baselines in drug repurposing.
In general, baseline models can be classified into two categories: graph-embedding-based models (GCN and DeepWalk), and path-based models (LSTM and KPRN). iDPath presents a modeling framework that combines graph-embedding-based and path-based approaches. We compare the performance of iDPath with baselines in the drug repurposing problem. As shown in Figure 2a and  not been correctly classified (Figure 2b). The utilization of the shortest paths can significantly improve the performance, as demonstrated by the superior performance of iDPath and path-based models over graph-embedding-based models (Figure 2a and c). These results indicate that the extracted MODA-related biological paths have pharmacological relevance.

Incorporating multiple biological network layers improves the prediction performance.
We further investigate the performance of iDPath with different biological network layers. Existing studies mainly make predictions using the PPI network alone 13,14 . Here, we evaluate the performance of iDPath with only one layer (PPI, GRN or PCI). Note that CCI cannot be directly linked to diseases, so we do not evaluate the model with only the CCI layer. As shown in Figure 2d, the full multilayer biological network can improve the performance of iDPath. Comparing individual networks, PPI performs the best, followed by GRN, and finally PCI. Note that the iDPath variants combining two or three layers cannot compete with the model with all four layers. Then, we investigate the proportion of nodes and interactions at each network layer in the identified MODA-related biological paths (Figure 3), to examine the impact of different network layers on iDPath performance. We find that nodes and interactions in the PPI layer are not the most prevalent in the identified paths. Instead, GRN nodes and PCI interactions are more common.
Combining these results with the dominating role of the PPI network in prediction performance, we find that although the connectivity in the PPI network can capture the key relationships between drugs and diseases, it requires much more information at the GRN, PCI, and CCI layers to further reveal the hidden biological paths related to MODA. By revealing these hidden paths, iDPath achieves higher prediction accuracy. shortest paths between abiraterone and prostate cancer, iDPath prefers the paths traversing through the gene targeted by both abiraterone and prostate cancer, such as abiraterone → CYP3A4 → prostate cancer and abiraterone → AR → prostate cancer. Previous studies have shown that abiraterone is a moderate inhibitor of CYP3A4 62 , and CYP3A4 is associated with oxidative deactivation of testosterone, which is the etiology of prostate cancer 63 . Androgen receptor (AR) is highly relevant to the growth and differentiation of prostate cancer 64 , and abiraterone inhibits androgen biosynthesis to control the progression of prostate cancer 65 . Abiraterone is found to be an inhibitor of CYP17A1, 66 which has also been identified by iDPath.

16/26
prostate cancer contributes to the prediction the most among all the CYP17A1 related paths, which is also consistent with previous biological studies 67 .
To further validate that the paths with higher weights are more relevant to the progression of prostate cancer, we perform the KEGG Pathway enrichment analysis 68 on the proteins of the top 50 paths and bottom 50 paths for the abiraterone-prostate cancer pair. As shown in Figure 4b and 4c, the paths with higher weights focus on the P13K-Akt signaling pathway 69 , regulation of actin cytoskeleton 70 , prostate cancer pathway, and so on, which are highly related to the progression of prostate cancer. For example, the activation of P13K-Akt signaling pathway appears to be characteristic of many aggressive prostate cancers and is more frequently observed as prostate cancer progresses toward a resistant and metastatic disease 69 . In contrast, the paths with lower weights (Figure 4c) are more enriched in the pathways related to other cancers or more general cancer progression, not specific to prostate cancer.

Drug repurposing for prostate cancer.
To demonstrate iDPath's utility in the real-world setting, we apply it to the discovery of potential drugs for treating prostate cancer among 1,572 FDA-approved drugs that have not been labeled as therapeutic drugs for prostate cancer in our dataset. The ten drugs with the highest score, together with their top-3

18/26
critical paths and top-3 KEGG pathways, are shown in Table 1. Among the ten drugs identified for prostate cancer, six drugs have already been proved effective in previous studies, including dutasteride 73 , aspirin 74 , erlotinib 75 , midostaurin 76 , apalutamide 77 , and atorvastatin 78 . Among them, apalutamide has recently been approved for the treatment of prostate cancer 79 , but has not been labeled in our dataset. Specifically, iDPath finds that the most relevant paths for the efficacy of Apalutamide are through Enzalutamide or Abiraterone (both are FDA-approved drugs for the treatment of prostate cancer and labeled as therapeutic in our dataset), where the combination with Abiraterone has already been proved synergistic in a recent study 80 .
For other drugs identified as therapeutic but not officially approved, the KEGG pathway enrichment analysis shows that the proteins that existed in their critical paths enriched in prostate cancer-related pathways, such as PI3K-Akt signaling pathway 81 and prostate cancer pathway.

Conclusion
In this study, we propose iDPath, an advanced deep learning framework to identify explainable biological paths to characterize the mechanism of drug actions and predict the drugs that can be repurposed for treating certain diseases. iDPath is built on a multilayer biological network consisting of GRN, PPI, PCI, and CCI networks. The proposed model achieves superior prediction performance compared to state-of-the-art models on a general drug repurposing task. Furthermore, we find that extending the PPI network to a multilayer biological network of the human body can significantly improve the prediction performance in drug repurposing. We investigate the identified critical paths of drugs to prostate cancer and find that the critical paths are consistent with the known mechanism of the drug action. At last, we apply iDPath to the challenging problem of identifying potential drugs for the treatment of prostate cancer.
Results show that iDPath can effectively identify the newly approved drugs not recorded by the database.
We believe iDPath can bring revelation to the explainable deep learning technologies to drug discovery. Chemical carcinogenesis -reactive oxygen species Table 1. The top-3 critical paths and top-3 KEGG pathways of the potential drugs for the treatment of prostate cancer. These drugs are ranked top 10 by iDPath among all the FDA-Approved drugs used in this study. The head (drug) and tail (prostate cancer) of these critical paths are ignored due to the limit of space. The top-3 critical paths are determined by the weights generated by the path attention module. The KEGG pathways are identified by KEGG enrichment analysis on the proteins existed in the top-50 critical paths and ranked by p-adjust values.

19/26
As a deep learning approach, iDPath is limited to in silico study, which can be extended by in virto and in vivo experiments to further validate its practical value and consistency with clinical evidence in future studies. In addition, the identified paths may contain rich biological knowledge beyond this study, such as some popular paths may be associated with common mechanisms of action in a class of diseases, which is worth further study.