Quantitative Systems Pharmacological Analysis of Drugs of Abuse Reveals the Pleiotropy of Their Targets and the Effector Role of mTORC1

Existing treatments against drug addiction are often ineffective due to the complexity of the networks of protein-drug and protein-protein interactions (PPIs) that mediate the development of drug addiction and related neurobiological disorders. There is an urgent need for understanding the molecular mechanisms that underlie drug addiction toward designing novel preventive or therapeutic strategies. The rapidly accumulating data on addictive drugs and their targets as well as advances in machine learning methods and computing technology now present an opportunity to systematically mine existing data and draw inferences on potential new strategies. To this aim, we carried out a comprehensive analysis of cellular pathways implicated in a diverse set of 50 drugs of abuse using quantitative systems pharmacology methods. The analysis of the drug/ligand-target interactions compiled in DrugBank and STITCH databases revealed 142 known and 48 newly predicted targets, which have been further analyzed to identify the KEGG pathways enriched at different stages of drug addiction cycle, as well as those implicated in cell signaling and regulation events associated with drug abuse. Apart from synaptic neurotransmission pathways detected as upstream signaling modules that “sense” the early effects of drugs of abuse, pathways involved in neuroplasticity are distinguished as determinants of neuronal morphological changes. Notably, many signaling pathways converge on important targets such as mTORC1. The latter emerges as a universal effector of the persistent restructuring of neurons in response to continued use of drugs of abuse.


INTRODUCTION
Drug addiction is a chronic relapsing disorder characterized by compulsive, excessive, and self-damaging use of drugs of abuse. It is a debilitating condition that potentially leads to serious physiological injury, mental disorder and death, resulting in major health, and social economic impacts worldwide (Nestler, 2013;Koob and Volkow, 2016). Substances with diverse chemical structures and mechanisms of action are known to cause addiction. Except for alcohol and tobacco, substances of abuse are commonly classified into six groups based on their primary targets or effects: cannabinoids (e.g., cannabis), opioids (e.g., morphine, heroin, fentanyl), central nervous system (CNS) depressants (e.g., pentobarbital, diazepam), CNS stimulants (e.g., cocaine, amphetamine), hallucinogens (e.g., ketamine, lysergic acid diethylamide), and anabolic steroids (e.g., nandrolone, oxymetholone).
The primary actions of drugs of abuse have been well studied. In spite of the pleiotropy and heterogeneity of drugs of abuse, they share similar phenotypes: from acute intoxication to chronic dependence (Taylor et al., 2013), the reinforcement shift from positive to negative through a three-stage cycle involving binge/intoxication, withdrawal/negative effect, and preoccupation/anticipation . Notably, virtually all drugs of abuse augment dopaminergic transmission in the reward system (Wise, 1996). However, the detailed cellular pathways of addiction processes are still far from known. For example, cocaine acts primarily as an inhibitor of dopamine (DA) transporter (DAT) and results in DA accumulation in the synapses of DA neurons (Shimada et al., 1991;Volkow et al., 1997). However, it has been shown that DA accumulation per se is not sufficient to account for the rewarding process associated with cocaine addiction; serotonin (5-HT) and noradrenaline (or norepinephrine, NE) also play important roles (Rocha et al., 1998;Sora et al., 1998). Another example is ketamine, a non-selective antagonist for N-methyl-d-aspartate (NMDA) receptor (NMDAR), notably most effective in the amygdala and hippocampal regions of neurons (Collingridge et al., 1983). In addition to its primary action, ketamine affects a number of other neurotransmitter receptors, including sigma-1 (Mendelsohn et al., 1985), substance P (Okamoto et al., 2003), opioid (Hustveit et al., 1995), muscarinic acetylcholine (mACh) (Hirota et al., 2002), nicotinic acetylcholine (nACh) (Coates and Flood, 2001), serotonin (Kapur and Seeman, 2002), and γ-aminobutyric acid (GABA) receptors (Hevers et al., 2008). The promiscuity of drugs of abuse brings an additional layer of complexity, which prevents the development of efficient treatment against drug addiction.
In recent years, there has been significant progress in the characterization of drug/target/pathway relations driven by the accumulation of drug-target interactions and pathways data, as well as the development of machine learning, in silico genomics, chemogenomics, and quantitative systems pharmacology (QSP) tools. Several innovative studies started to provide valuable information on substance abuse targets and pathways. For example, Li et al. curated 396 drug abuse related genes from the literature and identified five common pathways underlying the reward and addiction actions of cocaine, alcohol, opioids, and nicotine (Li et al., 2008). Hu et al. analyzed the genes related to nicotine addiction via a pathway and network-based approach (Hu et al., 2018). Biernacka et al. performed genome-wide analysis on 1,165 alcohol-dependence cases and identified two pathways associated with alcohol dependence (Biernacka et al., 2013). Xie et al. generated chemogenomics knowledgebases focused on G-protein coupled receptors (GPCRs) related to drugs of abuse in general (Xie et al., 2014), and cannabinoids in particular (Xie et al., 2016). Notably, these studies have shed light on selected categories or subgroups of drugs. There is a need to understand the intricate couplings between multiple pathways implicated in the cellular response to drugs of abuse, identify mechanisms common to various categories of drugs while distinguishing those unique to selected categories.
We undertake here such a systems-level approach using a dataset composed of six different categories of drugs of abuse. Following a QSP approach proposed earlier (Stern et al., 2016), we provide a comprehensive, unbiased glimpse of the complex mechanisms implicated in addiction. Specifically, as shown in Figure 1, a set of 50 drugs of abuse with a diversity of chemical structures (Supplementary Figure 1) and pharmacological actions were collected as probes, and the known targets of these drugs as well as the targets predicted using our probabilistic matrix factorization (PMF) method (Cobanoglu et al., 2013) were analyzed to infer biological pathways associated with drug addiction. Our analysis yielded 142 known and 48 predicted targets and 173 pathways permitting us to identify both generic mechanisms regulating the responses to drug abuse as well as specific mechanisms associated with selected categories, which could facilitate the development of auxiliary agents for treatment of addiction.
A key step in our approach is to identify the targets for drugs of abuse. There exists various drug-target interaction databases (DBs), web servers and computational models, as summarized recently (Chen et al., 2016). The DBs utilized in this work are the drug-target database DrugBank (Wishart et al., 2018) and the protein-chemical database STITCH (Szklarczyk et al., 2016). DrugBank is a bioinformatics and cheminformatics resource that combines drug data with comprehensive target information. It is frequently updated, with the current version containing 10,562 drugs, 4,493 targets and corresponding 16,959 interactions. Since most of drugs of abuse are approved or withdrawn drugs, DrugBank is a good source for obtaining information on their interactions. STITCH, on the other hand, is much more extensive. It integrates chemical-protein interactions from experiments, other DBs, literature and predictions, resulting in data on 430,000 chemicals and 9,643,763 proteins across 2,031 genomes. We have used in the present analysis the subset of human protein-chemicals data supported by experimental evidence. The method of approach adopted here is an important advance over our original PMF-based machine learning methodology for predicting drug-target interactions (Cobanoglu et al., 2013). First, the approach originally developed for mining DrugBank has been extended to analyzing the STITCH DB, the content of which is 2-3 orders of magnitude larger than DrugBank (based on the respective numbers of interactions). Second, the information on predicted drug-target associations is complemented by pathway data on humans inferred from the KEGG pathway DB (December 2017 version; Kanehisa et al., 2017) upon pathway enrichment analysis of known and predicted targets. Third, the outputs are subjected to extensive analyses to detect recurrent patterns and formulate new hypotheses for preventive or therapeutic strategies against drug abuse. FIGURE 1 | Workflow of the quantitative systems pharmacological analysis. (A) 50 drugs of abuse with a diversity of chemical structures and pharmacological actions were collected as probes. (B) 142 known targets of these drugs were identified through drug-target interaction database DrugBank and chemical-protein interaction database STITCH. (C) 48 predicted targets were predicted using our probabilistic matrix factorization (PMF) method (Cobanoglu et al., 2013). (D) 173 human pathways were inferred from the KEGG pathways database by mapping the known and predicted targets. (E,F) The pathways were grouped into 5 clusters. The functioning of identified targets and pathways and their involvement in drug addiction were comprehensively examined.

Selection of Drugs of Abuse and Their Known Targets
We selected as input 50 drugs commonly known as drugs of abuse using two basic criteria: (i) diversity in terms of structure and mode of action, and (ii) availability of information on at least one human target protein in DrugBank v5 (Wishart et al., 2018) or STITCH v5 (Szklarczyk et al., 2016). The selected drugs represent six different categories: CNS stimulants, CNS depressants, opioids, cannabinoids, anabolic steroids, and hallucinogens (see Supplementary Table 1 and Supplementary Figure 1).
A dataset of 142 known targets, listed in Supplementary Table 2, were retrieved from DrugBank and STITCH DBs for these 50 drugs. The list includes all targets reported for these drugs in DrugBank, and those with high confidence score, based on experiments, reported in STITCH. Each chemical-target interaction is annotated with five confidence scores in STITCH: experimental, DB, text-mining, prediction, and a combination score of the previous four, each ranging from 0 to 1. We selected the human protein targets with experimental confidence scores of 0.4 or higher. Supplementary Table 2 summarizes the 142 targets we identified as well as the associated 445 drug-target interactions.
Structure-based and interaction-pattern-based similarities between pairs of drugs were evaluated using two different criteria. The former was based on structure-based distance calculated as the Tanimoto distance between their 2D structure fingerprints. Tanimoto distances were evaluated using Python RDKit suite (RDKit: Open-Source Cheminformatics Software. https://www.rdkit.org/). Similarities based on their interactions patterns with known targets were evaluated by evaluating target-based distances. To this aim, we represented each drug i by a 142-dimensional "target vector" d i , the entries of which represent the known targets and are assigned values of 0 or 1, depending on the existence/observation of an interaction between the corresponding target and drug i. Interaction-pattern similarities between drug pairs i and j were evaluated by calculating the correlation cosine cos(d i . d j ) = (d i . d j )/(|d i | |d j |) between these vectors, and the corresponding cosine distance is [1-cos(d i . d j )]. Likewise, ligand-based distances between target pairs i and j were evaluated as the cosine distance between the 50-dimensional vectors t i and t j corresponding to the two targets, the entries of which are 0 or 1 depending on absence or existence of an interaction between the target and the corresponding drug of abuse.

Probabilistic Matrix Factorization (PMF) Based Drug-Target Interaction Prediction
Novel targets for each drug were predicted using our probabilistic matrix factorization (PMF) based machine learning approach (Cobanoglu et al., 2013(Cobanoglu et al., , 2015. Briefly, we start with a sparse matrix R representing the known interactions between N drugs and M targets. Using the PMF algorithm, we decomposed R into a drug matrix U and a target matrix V, by learning the optimal D latent variables to represent each drug and each target. The product of U T and V assigns values to the unknown (experimentally not characterized) entries of the reconstructed R, each value representing the confidence score for a novel drug-target interaction prediction

R N×M = U T N×D V D×M
Using this method, we trained two PMF models, one based on 11,681 drug-target interactions between 6,640 drugs and 2,255 targets from DrugBank v5, and the other based on 8,579,843 chemical-target interactions for 311,507 chemicals and 9,457 targets from STITCH v5 human experimentally confirmed subset, respectively. We evaluated the confidence scores in the range [0, 1] for each predicted drug-target interaction, in both cases. We selected the interactions with confidence scores higher than 0.7 within the top 10 predicted targets for each input drug. This led to 161 novel interactions identified between 27 out of the 50 input drugs and 89 targets (composed of 41 known and 48 novel targets; Supplementary Table 3).

Pathway Enrichment Analysis
We mapped the 50 drugs with 142 known and 48 predicted targets to the KEGG pathways (version December 2017, homo sapiens) (Kanehisa et al., 2017). 114 and 173 pathways were mapped by 142 known targets and all targets (both known and predicted) respectively (see Supplementary Table 4). In order to prioritize enriched pathways, we calculated the hypergeometric p-values based on the targets as the enrichment score as follows. Given a list of targets, the enrichment p-value for pathway A (P A ) is the probability of randomly drawing k 0 or more targets that belong to pathway A: where M is the total number of human proteins in the KEGG Pathway, m is the total number of proteins/targets we identified, and K is the number of proteins that belong to pathway A, while k 0 is the number of targets we identified that belong to pathway A. The obtained p-values are adjusted by a False Discovery Rate (FDR) correction to account for multiple testing, using the widely used Benjamini-Hochberg method (Benjamini and Yekutieli, 2001). The cutoff of the adjusted p-values gives us an upper bound of the false discovery rate. The false discovery rate is the fraction of false significant pathways maximally expected from the significant pathways identified in our case. We sort pvalues from smallest to largest, with m being the total number of pathways. The adjusted p-value, p * i , corresponding to the ith pathway is:   Figure 2D). However, we also note outliers, such as cocaine lying among opioids, as opposed to its categorization as a CNS stimulant, or promethazine, a CNS depressant, lying among hallucinogens (shown by arrows). The peculiar behavior of cocaine is consistent with its high promiscuity (see Figure 3A for the number of targets associated with each examined drug). This type of promiscuity becomes even more apparent when the drugs are organized based on their structure (or 2D fingerprints; see section Materials and Methods) as may be seen in Figure 2A.
For example, opioids (cyan labels/arc; clustered together in Figures 2B,D based on their interactions) are now distributed in two or more branches of the structure-based dendrogram in Figure 2C; likewise, CNS depressants (blue) and cannabinoids (light brown), grouped each as a single cluster in target-based dendrograms in Figure 2D, are now distributed into two or more clusters in Figure 2C.
Overall these results suggest that the functional categorization of the drugs does not necessarily comply with their structural characteristics. The similar functionality presumably originates from targeting similar pathways, but the difference in the structure suggests that either their targets, or the binding sites on the same target, are different; or the binding is not selective enough such that multiple drugs can bind the same site. Consequently, a diversity of pathways or a multiplicity of cellular responses are triggered by the use and abuse of these drugs.

The Selected Drugs and Identified Targets Are Highly Diverse and Promiscuous
We evaluated the similarities between proteins targeted by drugs of abuse, based on their interaction patterns with the studied drugs of abuse. Figures 2E,F display the respective target-target distances, and corresponding dendrogram. Supplementary Table 2 lists the full names of these targets, organized in the same order as the Figure 2E axes. We discern several groups of targets clustered together in consistency with their biological functions. For example, practically all GABA receptor subtypes (brown) are clustered together. This large cluster also includes the riboflavin transporter 2A (SLC52A2), which may be required for GABA release (Tritsch et al., 2012). On the other hand, the different subtypes of serotonin (or 5-hydroxytryptamine, 5-HT) receptors (5HTRs) participate in distinct clusters pointing to the specificity of different subtypes vis-à-vis different drugs of abuse (labeled in Figure 2F).
The large majority of neurotransmitter transporters, such as Na + /Cl − -dependent GABA transporters (SLC6A1) and glycine transporter (SLC6A9) are in the same cluster (pink, labeled). Acetylcholine receptors also lie close to (or are even interspersed among) Na + /Cl − -dependent neurotransmitter transporters, presumably due to shared drugs such as cocaine. However, the three transporters playing a crucial role in developing drug addiction, DAT, NE transporter (NET) and serotonin transporter (SERT) (labeled SLC6A2: NET, SLC6A3: DAT, SLC6A4: SERT) are distinguished by from all other neurotransmitter transporters as a completely disjoint group. The corresponding branch of the dendrogram (highlighted by the yellow circle) also includes vesicular amino acid transporters and trace amine-associated receptor 1 (TAAR1) known to interact with these transporters (Miller, 2011). We also note in the same branch two seemingly unrelated targets: flavin monoamine oxidase which draws attention to the role of oxidative events; and α2-adrenergic receptor subtypes A-C, which uses NE as a chemical messenger for mediating stimulant effects such as sensitization and reinstatement of drug seeking, and adenylate cyclase as another messenger to regulate cAMP levels (Sofuoglu and Sewell, 2009).
Supplementary Table 2 summarizes the 445 known interactions between these 50 drugs and 142 targets. We observe an average of 8.9 interactions per drug and 3.1 interactions per target. There are 23 promiscuous drugs that target at least 10 proteins as shown in Figure 3A. Cocaine, the most promiscuous psychostimulant, interacts with 45 known, and 3 predicted targets. It is known that cocaine binds DAT to lock it in the outward-facing state (OFS) and block the reuptake of DA. It similarly antagonizes SERT and NET (Heikkila et al., 1975;Sora et al., 1998), and also affects muscarinic acetylcholine receptors (mAChRs) M1 and M2 (Williams and Adinoff, 2008). Our PMF model also predicted a potential interaction between cocaine and M5. While this interaction is not listed in current DBs, there is experimental evidence suggesting that muscarinic AChR M5 plays an important role in reinforcing the effects of cocaine (Fink-Jensen et al., 2003), in support of the PMF model prediction.
The PMF model enables us to predict novel targets. For example, anabolic steroid nandrolone has only two known interactions, and cannabinoid cannabichromene has one. However, 10 new targets were predicted with high confidence scores for each of them (Supplementary Table 3 and Supplementary Figure 2A). This is due to the data available in STITCH DB, which offers a large training dataset that enhances the performance of our machine learning approach. Overall, 89 new interactions were predicted for known targets, and 42 novel targets were predicted with 72 interactions. Figure 3C displays the distribution of all targets among different protein families. As will be further elaborated below, among the newly identified drug-target pairs, nandrolone-MAPK14 (mitogen-activated protein kinase 14, also known as p38α) and canabichromene-IKBKB (inhibitor of NFκ-B kinase subunit β) play a role in regulating mTORC1 signaling, which will be shown to be a potential effector of drug addiction.
Turning to targets, three opioid receptors (OPRM1, OPRD1, and OPRL1) exhibit the highest level of promiscuity (Supplementary Figure 2B). The µ-type opioid receptor (OPRM1) interacts with 14 known drugs including all opioids as well as ketamine and dextromethorphan. We also predicted a novel interaction between OPRM1 and the CNS stimulant methylphenidate. This is consistent with experimental observations that methylphenidate upregulates OPRM1's activity in the reward circuitry in a mouse model (Zhu et al., 2011). Furthermore, tissue-based transcriptome analysis (Uhlén et al., 2015) shows that 69% of our 190 targets are expressed in the brain, and 49 of them show elevated expression levels in the brain compared to other tissue types (Supplementary Table 5). Among all the targets, NMDA receptor 1 (GRIN1) shows the highest elevated expression. It is also one of the top 5 enriched genes overall in the brain (Uhlén et al., 2015).
Taken together, the 50 selected drugs of abuse and the 142 known and 48 novel targets we identified cover a diversity of biological functions, are involved in many cellular pathways, and are generally promiscuous. In order to reveal the common mechanisms that underlie the development and escalation of drug addiction and also distinguish the effects specific to selected drugs, we proceed now to a detailed pathway analysis, presented next.

Pathway Enrichment Analysis Reveals the Major Pathways Implicated in Various Stages of Addiction Development
Our QSP analysis yielded a total of 173 pathways, including 114 associated with the known targets of the examined dataset of drugs of abuse, and 59 associated with the predicted targets. The detailed pathway enrichment results can be found in Supplementary Table 4. These pathways can be grouped in five categories (Figure 4; Supplementary Figures 3, 4, and Supplementary Table 4):

Synaptic Neurotransmission (NT)
Six significantly enriched (with adjusted p-value < 0.05) pathways are associated with synaptic neurotransmission: dopaminergic, serotonergic, glutamatergic, synaptic vesicle cycle, cholinergic, and GABAergic synapses pathways. Sixty-eight known targets and 7 predicted targets are involved in these pathways. This is consistent with the fact that neurotransmission plays a dominant role in the rewarding system and is key to drug addiction (Volkow and Morales, 2015).

Signal Transduction (SG)
Forty-six intracellular signaling pathways were mapped by 92 targets comprised of 66 known and 25 predicted targets. Notably, many of these pathways have been reported to play a role in mediating the effects of drugs of abuse. These include the top five [calcium signaling (Li et al., 2008), retrograde endocannabinoid signaling (Mechoulam and Parker, 2013), cGMP-PKG signaling (Shen et al., 2016), cAMP signaling (Philibin et al., 2011), and Rap1 signaling (Cahill et al., 2016)] as well as some pathways with relatively low enrichment score (i.e., 0.2 < adjusted p-value), such as TNF signaling (Zhu et al., 2018), MAPK signaling , PI3K-Akt signaling (Neasta et al., 2011), NF-κB signaling (Nennig and Schank, 2017), and mTOR signaling (Neasta et al., 2014). We note that many receptors targeted by drugs of abuse take part in the KEGG neuroactive ligand-receptor interaction pathway. In the interest of focusing on intracellular signaling effects, we have not included these in the SG category; they are listed in the "Other Pathways" in Supplementary Table 4.

Autonomic Nervous System (ANS)-Innervation (ANS)
We also identified 10 pathways regulating ANS-innervated systems such as endocrine secretion, taste transduction, and circadian entrainment. Recent evidences suggested drugs of abuse such as morphine (Al-Hasani and Bruchas, 2011) and cocaine (Moeller et al., 1997;Prosser et al., 2014) can influence ANS-innervated systems and may contribute to the withdrawn symptoms associated with drug addiction. Thirty-seven known and 9 predicted targets take part in these pathways.
Neuroplasticity (NP). Eight enriched pathways with potential to alter the morphology of neurons, were found to be related to drug addiction. Among them, long-term potentiation (LTP) and long-term depression (LTD) are key to reward-related learning and addiction by modifying the fine tuning of dopaminergic firing (Jones and Bonci, 2005). Axon guidance pathway regulates the growth direction of neuron cells (Bahi and Dreyer, 2005). Regulation of actin cytoskeleton plays important role in morphological development and structural changes of neurons (Luo, 2002). Gap junctions connect neighboring neurons via intercellular channels that allow direct electrical communication (Belousov and Fontes, 2013) and regulate the efficiency of communication between electrical synapses (Belousov and Fontes, 2013). Nineteen known targets and 5 predicted targets are involved in these pathways. Insulin-like growth factor 1 receptor (IGF1R) is predicted as a target of drug triazolam (Supplementary Table 4). IGF1R is involved in LTP, adherens junction and focal adhesion pathways. It functions via canonical signaling pathways noted above in the SG category, such as the PI3K-Akt-mTOR and Ras-Raf-MAPK pathways (Lee et al., 2016) and it plays important role in neuroplasticity (Lee et al., 2016). We note that the NP group involves many pathways directly relevant to drug addiction (Bahi and Dreyer, 2005;Kalivas and Volkow, 2011;Moradi et al., 2013;Rothenfluh and Cowan, 2013). There is no target unique to this particular group of pathways ( Figure 4B). However, the fact that the targets belonging to the NP group are also shared by other groups consolidates the significance of these targets.

Disease-Associated Pathways (DS)
Fifty enriched pathways mapped by 51 known and 17 predicted targets are associated with diverse diseases in different organs such as brain, liver, and lung. They also cover various drug addiction mechanisms including: nicotine addiction, morphine addiction, cocaine addiction, amphetamine addiction, and alcoholism. Additionally, there are "other pathways" such as those involved in cell migration, differentiation, immune responses, and metabolic events, which can be seen in Supplementary Table 4.
Taken together, the enrichment analysis reveals five major categories of pathways that regulate the three stages of drug addiction cycle: (1) binge and intoxication, (2) withdrawal and negative affect, and (3) preoccupation and anticipation (or craving) (Koob and Volkow, 2010). Drugs of abuse directly affect neurotransmission pathways: they increase the accumulation of DA and other neurotransmitters in the synaptic and extrasynaptic regions, which in turn results in the hedonic feeling (stage 1) and triggers the DA reward system. Dysregulation of ANS-innervation pathways may cause negative effects and feelings (stage 2) and feedback to the CNS. Addictive drugs impair executive processes by disrupting the reward system (neurotransmission pathways) and imparting morphological changes via neuroplasticity pathways (e.g., LTD and LTP), which then result in craving (stage 3). Below, we present an in-depth analysis of the role of these pathways or their shared targets in drug addiction.

Selected Targets Shared by Dominant Pathways Emerge as Common Mediators of Drug Addiction
We next analyzed the overlapping targets between the pathways in different functional categories.
AMPA receptor plays a crucial role in LTP and LTD, which are vital to neuroplasticity, memory and learning . Serotonin receptors, expressed in both the CNS and the peripheral nervous system (e.g., gastrointestinal tract), are responsible for anxiety, impulsivity, memory, mood, sleep, thermoregulation, blood pressure, gastrointestinal motility, and nausea (Pytliak et al., 2011). They have been proposed to be therapeutic targets for treating cocaine use disorder (Howell and Cunningham, 2015). RAC1 is involved in five neuroplasticity pathways, including axon guidance, adherens junction and tight junction pathways (Supplementary Table 4), and 13 intracellular signal transduction pathways. It regulates neuroplasticity, as well as apoptosis and autophagy (Natsvlishvili et al., 2015). DA receptor D 2 is a target of 28 drugs of abuse (out of 50 examined here) and is involved in cAMP signaling, and gap junction pathways, in addition to dopaminergic signaling. It is implicated in reward mechanisms in the brain (Blum et al., 1996) and the regulation of drug-seeking behaviors (Edwards et al., 2006). Finally, PI3K turns out to be the most pleiotropic target among those targeted by drugs of abuse, being involved in 61 pathways identified here, including neuroplasticity pathways such as axon guidance, and several downstream signaling pathways such as PI3K-Akt, mTOR, Ras and Jak-STAT pathways.
Overall, the above listed 23 proteins shared by at least four different groups of pathways are distinguished here as highly pleiotropic proteins involved in the large majority of pathway categories implicated in drug abuse. Most of them are ligandor voltage-gated ion channels or neurotransmitter receptors, mainly AMPAR, NMDAR, Cav2.1, mAChR, and serotonin and DA receptors. However, it is interesting to note the targets PI3K and p38α, not currently reported in DrugBank and STITCH, emerge as highly pleiotropic targets of the drugs of abuse. These are suggested by the current analysis to directly or indirectly affect addiction development and await future experimental validation. Finally, a number of proteins take part in specific drug-abuse-related pathways and might serve as targets for selective treatments. Supplementary Table 6 provides a list of such targets uniquely implicated in distinctive pathways.

Pathway Enrichment Highlights the Interference of Drugs of Abuse With Synaptic Neurotransmission
It is broadly known that neurotransmitters such as DA, 5-HT, NE, endogenous opioids, ACh, endogenous cannabinoids, Glu, and GABA are implicated in drug addiction (Tomkins and Sellers, 2001;Everitt and Robbins, 2005;Parolaro and Rubino, 2008;Benarroch, 2012). Our analysis also showed that the serotonergic synapse (adjusted p-value p * i = 2.01E-18), GABAergic synapse (p * i = 1.19E-17), cholinergic synapse (p * i = 2.36E-07), dopaminergic synapse (p * i = 1.66E-06) and glutamatergic synapse (p * i = 1.86E-03) pathways were significantly enriched (Supplementary Table 4). A total number of 34 drugs (across six different groups) target at least one of these pathways. However, the identification of a pathway does not necessarily mean that the drug directly affects that particular neurotransmitter transport/signaling. There may be indirect effects due to the crosstalks between synaptic signaling pathways. For example, the ionotropic glutamate receptors NMDAR and AMPAR are also the downstream mediators in the dopaminergic synapse pathway. Likewise, GABARs are downstream mediators in the serotonergic synapse pathway.
In Figure 5, we highlight five major neurotransmission events that directly mediate addiction, and illustrate how eight drugs of abuse interfere with them. Despite the promiscuity of the drugs of abuse, some selectively map onto a single synaptic neurotransmission pathway. For example, psilocin [a hallucinogen whose structure is similar to 5HT (Diaz, 1997)] interacts with several types of 5HTRs, regulating serotonergic synapse exclusively (see Figure 5 and Supplementary Table 4). In contract, loperamide (not shown) affects all neurotransmission pathways by interacting with the voltage-dependent P/Q-type calcium channel (VGCC), regulating calcium flux on synapses. Cocaine targets four of these synaptic neurotransmission events (serotonergic, GABAergic, cholinergic, and dopaminergic synapses), through its interactions with 5-HT3R, sodiumand chloride-dependent GABA transporter (GAT), muscarinic (M1 and M2) and nicotinic AChRs, and DAT, respectively. Methadone affects three synaptic neurotransmissions, including serotonergic synapse, dopaminergic synapse, and glutamatergic synapse through the interactions with SERT, DAT, and glutamate receptors (NMDAR), respectively.
It is worth noting that the current analysis helps us generate new hypotheses, yet to be experimentally validated, on the ways drugs of abuse affect neurotransmission. In addition to the new role of the muscarinic AChR M5 suggested by the current analysis in section the selected drugs and identified targets are highly diverse and promiscuous, our PMF model suggested that cannabichromene, a cannabinoid whose primary target is the transient receptor (TRPA1), could interact with DAT and thus regulate dopaminergic transmission, which will require further examination.
The above synaptic neurotransmission events act as upstream signaling modules that "sense" the early effects of drug abuse. In the next section, we focus on the downstream signaling events elicited by drug abuse.

mTORC1 Emerges as a Potential Downstream-Effector Activated by Drugs Abuse
The calcium-, cAMP-, Rap1-, Ras-, AMPK-, ErbB-, MAPK-, and PI3K-Akt-signaling pathways in the SG category (Supplementary Table 4) crosstalk with each other and form a unified signaling network. As shown in Figure 6, FIGURE 5 | The impact of drugs of abuse on synaptic neurotransmission. Five major neurotransmission events are highlighted, mediated by (counterclockwise, starting from top): GABA receptors and transporters, ionotropic glutamate receptors (NMDAR and AMPAR) and cation channels, serotonin (5HT) receptors (5-HTR) and transporters (SERT), muscarinic or nicotinic AChRs, and dopamine (DA) receptors and transporters. Vesicular monoamine transporters (VMAT) that translocate DA are also shown. Drugs affecting the different pathways are listed, color coded with their categories, as presented in Figure 2. Solid red arrows indicate a known drug-target interaction, dashed red arrows indicate predicted drug-target interactions. Other molecules shown in the diagram are: KA, kainate receptor; MAO, monoamine oxidase; HVA, homovanillate; 3-MT, 3-methoxytyramine; MOR, mu-type opioid receptor; AChE, acetylcholinesterase; and 5-H1AA, 5-hydroxyindoleacetate.
ligand-binding to GPCRs modulates the production of cAMP, which leads to the activation of Rap1. Activated Rap1 modules the Ca 2+ signaling by inducing the production of inositol triphosphate (IP 3 ) and also activates the PI3K-Akt signaling cascade. Stimulations of ErbB family of receptor tyrosine kinases (related to epidermal growth factor receptor EGFR) as well as insulin-like growth factor receptor IGF1R trigger both PI3K-Akt and MAPK signaling cascades (proteins colored blue in Figure 6). Notably all these pathways merge and regulate a group of downstream proteins (shown in dark yellow in Figure 6); and at the center of this cluster lies the mammalian target of rapamycin (mTOR) complex 1 (mTORC1) which is likely to be synergistically regulated by all these merging pathways. mTORC1 is not only a master regulator of autophagy (Rabanal-Ruiz et al., 2017), but also controls protein synthesis and transcription (Ma and Blenis, 2009). It has been reported to promote neuroadaptation following exposure to drugs of abuse including cocaine, alcohol, morphine and 9tetrahydrocannabinol (THC) (Neasta et al., 2014). Our results lead to the hypothesis that mTORC1 may act as a universal effector of the cellular response to drug abuse at an advanced FIGURE 6 | A unified signaling network mediates the effects of drugs of abuse. Black arrows represent the activation, inhibition, and translocation events during signal transduction. Solid gray arrows represent the known drug-target interactions. Dashed gray arrows represent predicted drug-target interactions. The diagram illustrates the targets of several drugs of abuse belonging to different categories: loperamide, fentanyl, heroin, morphine, and methadone from opioids; midomafetamine, ketamine, dextromethorphan, LSD, and psilocin from hallucinogens; triazolam, diazepam, alprazolam, pentobarbital, eszopiclone, flunitrazepam, and zaleplon from CNS depressants; cannabichromene, 2-AG, cannabinol, and dronabinol from cannabinoids; methamphetamine, cocaine, AMPH, and phendimetrazine from CNS stimulants; and nandrolone from anabolic steroids. mTORC1 emerges as a hub where the effects on several targets of addictive drugs appear to be consolidated to lead to cell death and/or protein synthesis in the CNS, and in particular, to AMPAR/PSD95 synthesis that induces morphological changes in the dendrites.
(preoccupation and anticipation, or craving) stage, controlling the synthesis of selected proteins and ensuing cell growth, which may result in persistent alterations in the dendritic morphology and neuronal circuitry.
In Figure 6, selected interactions between drugs from different substance groups and their targets are highlighted using gray arrows. The figure illustrates that not only many known drug-target interactions, but also predicted ones involved in the unified signaling network. For example, our PMF model predicted that diazepam would interact with PI3K to influence mTORC1 signaling (dashed gray arrows denote predictions). It has been reported that Ro5-4864, a benzodiazepine derivative of diazepam suppresses activation of PI3K (Yousefi et al., 2013), which corroborates our prediction. We further predicted that cannabichromene may interact with IκB kinase β (IKKβ) to regulate mTORC1 by inhibiting TSC1/2. Interestingly, another cannabinoid, arachidonoyl ethanolamine, is known to directly inhibits IKKβ (Sancho et al., 2003). Taken together, our results suggest a unified network that underlies the development of drugs addiction, in which mTORC1 appears to play a key effector role.

DISCUSSION
In the present study we focused on the targets and pathways affected by drugs of abuse, toward gaining a systems-level understanding of key players and dominant interactions that control the response to drug abuse and the development of drug addiction. Using machine learning methods, we focused on 50 drugs of abuse that form a chemically and functionally diverse set, and analyzed their 142 targets as well as the corresponding cellular pathways and their crosstalk. Our analysis identified: (i) 48 additional proteins targeted by drugs of abuse, including PIK3CA, IKBKB, EGFR, and IGF1R, are shown to be key mediators of downstream effects of drug abuse.
(ii) 161 new interactions between the drugs of abuse and the known and predicted targets, including those between cocaine and M5, methylphenidate and OPRM1, and diazepam and PI3K, not reported in existing DBs, but supported by prior experiments, and others (e.g., the interactions of cannabichromene with IKBKB and DAT) that await experimental validation. (iii) A dataset of 70 pathways, composed of 6 neurotransmission pathways, 46 signal transduction pathways, 8 neuroplasticity pathways and 10 autonomic nervous system innervation pathways which are proposed to govern different stages of the molecular, cellular and tissue level responses to drug abuse and in addiction development.
Overall, our comprehensive analysis led to new hypotheses on drug-target interactions and signaling and regulation mechanism elicited by drugs of abuse in general, along with those on selected targets and pathways for specific drugs. Below we elaborate on the biological and biomedical implications of these findings.

Persistent Restructuring in Neuronal Systems as a Feature Underlying Drug Addiction
Enriched pathways in the neuroplasticity category include gap junction, LTP, LDP, adherens junction, regulation of actin cytoskeleton, focal adhesion, axon guidance, and tight junction (Supplementary Table 4). These are responsible for the changes in the morphology of dendrites. For instance, DA regulates excitatory synaptic plasticity by modulating the strength and size of synapses through LTP and LTD (De Roo et al., 2008;Volkow and Morales, 2015). The restructuring of dendritic spines involves the rearrangements of cytoskeleton and actin-myosin (Volkow and Morales, 2015). The axon guidance molecules guide the direction of neuronal growth. Drugs of abuse can induce the changes in CNS through these pathways. For example, chronic exposure to cocaine increases dendritic spine density in medium spiny neurons (Russo et al., 2010). The disruption in axon guidance pathway and alteration in synaptic geometry can result in drugrelated plasticity (Bahi and Dreyer, 2005). The persistent restructuring in the CNS caused by drugs of abuse is responsible for long-term behavioral plasticity driving addiction (Volkow et al., 2003;Russo et al., 2010;Volkow and Morales, 2015). As will be further discussed below, mTORC1 plays a central role in the synthesis of new proteins (e.g., AMPARs) and thereby neuronal (dendrites) growth, alteration of the synaptic geometry and therefore rewiring of the neuronal circuitry.

ANS May Mediate the Negative-Reinforcement of Drug Addiction
The current study further points to pathways regulating the ANS-innervated systems. As the NP pathways influence the neuroplasticity in the ANS, we hypothesize that drugs of abuse might induce a persistent restructuring in the ANS as well. The drug-related plasticity in ANS may lead to the dysregulation of ANS-innervated systems and cause negative effects and feelings during the second stage of drug addiction. Drug addiction is well known as a brain disease (Volkow and Morales, 2015). However, many drugs of abuse can disrupt the activity of ANS and cause disorders in ANS-innervated systems (Al-Hasani and Bruchas, 2011;Huang, 2017). For example, opioids (e.g., morphine) alter neuronal excitability and neurotransmission in the ANS (Wood and Galligan, 2004), and induce disorders in gastrointestinal system, smooth muscle, skin, cardiovascular, and immune system (Al-Hasani and Bruchas, 2011). Cannabinoids (e.g., THC) modulate the exocytotic NE release in ANSinnervated organs through presynaptic cannabinoid receptors (Ishac et al., 1996).
The pathways we identified in the ANS category regulate insulin secretion, gastric acid secretion, vascular smooth muscle contraction, pancreatic secretion, salivary secretion, and renin secretion (Supplementary Table 4). Their dysfunction may be associated with the autonomic withdrawal syndrome, such as thermoregulatory disorder (chills and sweats) and gastrointestinal upset (abdominal cramps and diarrhea), which has been observed in drug/substance users (Wise and Koob, 2014). In addition, the stress and depression caused by these negative effects may be part of the negative reinforcement of drug addiction (Self and Nestler, 1995;Koob and Le Moal, 2001). In other words, the drug induced ANS disorders can feedback to CNS and mediate the negative reinforcement. Compared to the structural changes in CNS, the disorder and persistent restructuring in ANS is less studied and it could be a future direction in the study of development of drug addiction and related diseases.

mTORC1 Appears as a Key Mediator of Cellular Morphological Changes Elicited in Response to Continued Drug Abuse
The functioning and regulation of mTOR signaling has been elucidated over the past two decades. It became clear that mTORC1 plays a crucial role in regulating diverse cellular processes including protein synthesis, autophagy, lipid metabolism, and mitochondrial biogenesis (Saxton and Sabatini, 2017). In the brain, mTORC1 coordinates neural development, circuit formation, synaptic plasticity, and long-term memory (Lipton and Sahin, 2014). The dysregulation of mTORC1 pathway is associated with many neurodevelopmental and neurodegenerative diseases such as Parkinson's disease and Alzheimer's disease. mTORC1 has been noted to be an important mediator of the development of drug addiction and relapse vulnerability (Dayas et al., 2012). Accumulating evidences show that pharmacological inhibition of mTORC1 (often through rapamycin treatment) can prevent sensitization of methamphetamine-induced place preference (Narita et al., 2005), reduce craving in heroin addicts (Shi et al., 2009), attenuate the expression of alcohol-induced locomotor sensitization (Neasta et al., 2010), suppress the expression of cocaine-induced place preference (Bailey et al., 2012), protect against the expression of drug-seeking and relapse by reducing AMPAR (GluA1) and CaMKII levels (James et al., 2014), and inhibit reconsolidation of morphine-associated memories (Lin et al., 2014).
Our unbiased computational analysis based on a diverse set of 50 drugs of abuse supports the hypothesis that mTORC1 may act as a universal effector or controller of neuroadaptations induced by drugs of abuse (Neasta et al., 2014). The major signal transduction pathways we identified that involve targets of drugs of abuse interconnect and converge to the mTORC1 signaling cascade (Figure 6). Most drugs of abuse in our list target upstream regulators of mTORC1, including membrane receptors (e.g., GPCRs, RTKs and NMDAR), kinases (e.g., PI3K, p38α, and IKKβ), and ion channels (e.g., Ca V 2.1 and TRPV2). Notably, the impact of some of these known or predicted targets has been experimentally confirmed. For example, blockade of the known target NMDAR using MK801 reduces the amnesic-like effects of cannabinoid THC (Puighermanal et al., 2009). Likewise, inhibition of PI3K (a predicted target) by LY294002 suppresses morphine-induced place preference in rats (Cui et al., 2010) and the expression of cocaine-sensitization (Izzo et al., 2002). Our results thus provide a pool of candidate targets implicated in cellular responses to addictive drugs, which await to be consolidated by further tests.
The downstream effectors of mTORC1, which specifically mediate drug behavioral plasticity is far from known. mTORC1 can mediate the activation of S6Ks and 4E-BPs, which leads to increased production of proteins required for synaptic plasticity including AMPAR and PSD-95 (Dayas et al., 2012). EM reconstruction of hippocampal neuropil showed the variability in the size and shape of dendrites depending on synaptic activity (Bartol Jr et al., 2015), which in turn correlates with information storage. Recently studies have revealed that Atg5-and Atg7-dependent autophagy in dopaminergic neurons regulates cellular and behavioral responses to morphine (Su et al., 2017). Cocaine exposure results in ER stress-induced and mTORC1-dependent autophagy (Guo et al., 2015). Fentanyl induces autophagy via activation of ROS/MAPK pathway (Yao et al., 2016). Methamphetamine induces autophagy through the κ-opioid receptor (Ma et al., 2014). These observations are consistent with the currently inferred role of mTORC1 as a downstream effector of cellular responses to drug addiction.

Drug Repurposing Opportunities for Combatting Drug Addiction
Autophagy modulating drugs have been shown to have therapeutic effects against liver and lung diseases. The signaling network presented in Figure 6 involves many targets of such drugs. For instance, carbamazepine affects IP 3 production and enhances autophagy via calcium-AMPK-mTORC1 pathway (Hidvegi et al., 2010). It has been identified as a potential drug for treating α1-antitrypsin deficiency, hepatic fibrosis, and lung proteinopathy (Hidvegi et al., 2010(Hidvegi et al., , 2015. Rapamycin is a potential drug for lung disease such as fibrosis (Abdulrahman et al., 2011;Patel et al., 2012). Other liver and lung drugs which facilitate the removal of aggregates by promoting autophagy may also affect drug-related neurodegenerative disorders. Supplementary Table 7 summarizes 15 autophagymodulating drugs for liver and lung diseases. Target identification and pathway analysis of this subset of drugs using the same protocol as those adopted for the 50 drugs of abuse indeed confirmed that drugs of abuse and liver/lung drugs share many common pathways (Supplementary Figure 5). Notably, among those pathways, neuroactive ligand-receptor interactions, calcium signaling, and serotonergic synapse pathways are among the top 10 enriched pathways of both drugs of abuse and liver/lung drugs. Amphetamine addiction and alcoholism are also enriched by targets of liver/lung drugs. Thus, an interesting future direction is to examine whether autophagy modulating drugs for liver and lung diseases could be repurposed, if necessary by suitable refinements to increase their selectivity, for treating drug addiction.
In summary, our results invite attention to new targets of addictive drugs and pathways implicated in the development of addiction, as well as new therapeutic opportunities. Recent studies support the utility of such computationally-driven QSP predictions. The validation of these predictions requires comprehensive wet-lab bioactivity assays (Pahikkala et al., 2015). In particular, the establishment of the proposed role of mTORC1 would require in vitro and in vivo longitudinal studies given that our current study points to the involvement of mTORC1 at later stages of drug addiction. In a recent study, we identified the role of protein kinase A (PKA) pathway in Huntington's disease using a QSP approach and verified experimentally (Pei et al., 2017). A similar combined computational-experimental framework could be adopted to extend the current study and establish new strategies. Though these experiments are beyond the scope of the current paper, our unbiased computational study provides insights into the pleiotropy of the targets of addictive drugs as well as the common signaling platforms that may serve as mediators of drug addiction.
Knowledge of pathways implicated in drug addiction may be used, as a next step, to construct kinetic models to quantitatively assess the orchestration of signals induced by pathway crosstalks. Our previous studies on Toll-like receptors  and cell fate decision processes  have demonstrated the utility of identifying such crosstalks for detecting synergistic response mechanisms and designing polypharmacological strategies. Therefore, the computational data presented here presents a milestone toward developing new therapies against drug addiction by identifying new targets beyond those usually investigated by focused studies. Finally, our analysis framework is generic and could be adopted for characterizing the targets and pathways of other complex disorders by suitable redefinition of the input set of drugs of interest.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript and/or the supplementary files.

AUTHOR CONTRIBUTIONS
FP, HL, and IB conceived and designed the research. FP and HL performed the research. FP, HL, BL, and IB analyzed the results and wrote the manuscript.

FUNDING
This work was supported by the National Institutes of Health awards P30DA035778, and P41GM103712.

ACKNOWLEDGMENTS
FP wish to thank Dr. D. Lansing Taylor for his mentorship.