In total, 32 network properties were calculated and adopted in subsequent analysis. These properties were popular for analyzing a complex biological network, which included: (1) Average Closeness Centrality: the average number of steps required to reach the studied node from any node in a network (Ma et al., 2016); (2) Average Shortest Path Length: the average length of shortest paths between the studied node and all other ones (Zhang et al., 2014); (3) Betweenness Centrality: the number of times the studied node serving as a linking bridge along shortest path between any two nodes (Zeidán-Chuliá et al., 2015); (4) Bridging Centrality: the product of the bridging coefficient and betweenness centrality (Hwang et al., 2008); (5) Bridging Coefficient: the extent of the studied node lying between any other densely connected nodes in the network (Paladugu et al., 2008); (6) Closeness Centrality Sum: the reciprocal of the sum of the shortest paths between the studied node and all other nodes in the network (Costenbader and ValenteFontanesi, 2003); (7) Clustering Coefficient: the number of the connected pairs between all neighbors of node (Watts and Strogatz, 1998); (8) Current Flow Betweenness: a centrality index measuring the level of information travels along all possible paths within network (Paladugu et al., 2008); (9) Current Flow Closeness: the variant of current flow betweenness (Zhang et al., 2017b); (10) Degree: the number of edges linked to a node (Braeuning, 2013); (11) Degree Centrality: the number of links incident upon a studied node (Batool and Niazi, 2014); (12) Deviation: the variation between sum of node distances and network unipolarity (Zhang et al., 2017a); (13) Distance Deviation: the absolute difference between nodes' distance sum and network's average distance (Rogelj et al., 2013); (14) Distance Sum: the sum of all shortest paths starting from the studied node (Bolser et al., 2003); (15) Eccentric: the absolute difference between nodes' eccentricities and network's average eccentricity (Zhang et al., 2017a); (16) Eccentricity: the maximum non-infinite shortest path length between the studied node and all other nodes in the network (Bolser et al., 2003); (17) Eccentricity Centrality: the largest geodesic distance between the node and any other node (Batool and Niazi, 2014); (18) Eigenvector Centrality: the sum of its neighbors' centrality values (Solá et al., 2013); (19) Harmonic Closeness Centrality: the sum of the reciprocals of the average shortest path lengths of each node in network (Zhang et al., 2017b); (20) Interconnectivity: a connectivity index indicating the quality of the studied nodes being connected together (Emig et al., 2013); (21) Load Centrality: the fraction of all the shortest paths that pass through the studied node (Kivimäki et al., 2016); (22) Neighborhood Connectivity: the average connectivity of all neighbors (Carson and Lu, 2015); (23) Normalized Betweenness: the fraction of network shortest paths that a given protein lies on (Paladugu et al., 2008); (24) Number of Self Loops: the number of edges starting and ending at the same node (Garlaschelli and Loffredo, 2004); (25) Number of Triangles: the number of triangles that include the studied node as a vertex (Rubinov and Sporns, 2010); (26) Page Rank Centrality: an adjustment of Katz by considering the diluted issue (Li et al., 2013); (27) Radiality: the level of reachability of a studied node via various shortest paths within the entire network (Koschützki and Schreiber, 2008); (28) Residual Closeness Centrality: the closeness measured by removing the studied node (Dangalchev, 2006); (29) Scaled Degree: the degree of a studied node relative to the most connected node within the same module (Sormani, 2012); (30) Stress: the number of shortest paths passing through a given node (Shannon et al., 2003); (31) Topological Coefficient: the extent to which a node in network shares interaction partners with other nodes (Zhu M. et al., 2009); (32) Z Score: a connectivity index based on degree distribution of a network (Rubinov and Sporns, 2010).

### Assessing the Biological System Profile for Each Studied Target

The biological system profile for each studied target included: (1) the number of target-affiliated and target immediate-downstream signaling pathways in KEGG database (Kanehisa et al., 2017). The target-affiliated pathways were determined by considering that (a) the pathways of the studied target should be lifeessential in both patients and healthy people and (b) the studied target should be in the pathway upstream with the capacity of regulating the biological function of the pathways. (2) The number of human tissues each target distributed in, assessed by the TissueDistributionDBs (Kogenaru et al., 2010) and Uniprot (UniProt Consortium, 2018) databases. A target was assumed to distribute in a given tissue if >5% of the total proteins are distributed in that tissue or the target concentration is higher than the average concentration of proteins in that tissue. (3) The number of human similarity proteins of a target outside the corresponding target family for probing off-target collateral effect (Zheng et al., 2006; Zhu F. et al., 2009). This was determined by BLAST similarity screening of human proteome in Uniprot database (UniProt Consortium, 2018) with a cutoff (E-value < 0.005; Song et al., 2006; Singh et al., 2007). (4) The differential expressions of the studied target in the diseasespecific tissue between patients and healthy individuals (Li et al., 2018). The relevant data were collected directly from TTD (Li et al., 2018) and calculated based on the human gene expression raw data of Affymetrix U133 Plus 2.0 platform in GEO (Barrett et al., 2013).

### Selecting the Differential Features Indicating NTI Drugs by Artificial Intelligence

The artificial intelligence (AI) has been recently proposed as a powerful technique for drug target discovery (Xu and Wang, 2014; Zhu et al., 2018), protein function prediction (Li et al., 2016a; Seo et al., 2018; Yu et al., 2018) and biomarker identification (Li B. et al., 2016; Li et al., 2017) through mimicking the human thinking procedures, learning processes and information extractions, which included the machine learning algorithm (Zhu et al., 2008a; Wang P. et al., 2015), the deep learning method (van der Burgh et al., 2017; Seo et al., 2018), and the cognitive-computing (Krittanawong et al., 2017). As one of the most popular machine learning algorithms, the Boruta algorithm based on wrapper method built around a random forest classifier (Kursa, 2014) was selected and adopted in this study. It is an extension to determine the relevance via comparing the relevance of the real features to that of the random probes (Pan et al., 2018). Since Boruta was constructed by an AI-based technique (machine learning), it was considered to be the most powerful approach with the stability in the variable selection, especially suitable for the low-dimensional dataset among other available strategies (Degenhardt et al., 2017). In this study, the differential features between NTI and NNTI drugs were therefore identified by R package Boruta (Shang et al., 2017). Particularly, human PPI network properties and biological system features of each target were first calculated, and the results of feature selection were then acquired using R package Boruta by setting the p-value < 0.05, maxRuns = 100, and doTrace = 2. In the meantime, the getImp was set to "getImpRfZ," and the mcAdj and holdHistory were set to "TRUE."

## RESULTS AND DISCUSSION

### Network Properties and Biological System Profile of NTI and NNTI Drugs

As reported, the human PPI network properties and biological system profile were key factors determining efficacy-safety balance (Zheng et al., 2006; Ragusa et al., 2010; Guo et al., 2018). Network properties were inherent feature of a target in the human PPI network, while biological system profile could reflect both the on-target and off-target pharmacology (Bender et al., 2007; Han et al., 2018; Zhu et al., 2018). Herein, 32 features of human PPI network together with 4 biological system properties were therefore adopted and calculated for further analyses. To the best of our knowledge, these were the most comprehensive sets of features ever applied for TI-related analysis. **Table 1** listed the calculated values of ten properties based on the connectivity and adjacency in human PPI network. These connectivity/adjacency-based network properties were designed to describe the level of connectivity among human proteins or the neighborhood features of the studied proteins (Chen et al., 2016). The properties included bridging coefficient, clustering coefficient, degree, degree centrality, interconnectivity, neighbor connectivity, number of triangles, scaled degree, topological coefficient, and Z-score (corresponding definitions were provided in section Materials and Methods). As shown in **Table 1**, 8 (80.0%) out of 10 properties were significantly different (p-value < 0.05, highlighted by bold font) between the targets of NTI and NNTI drugs, and half of those 10 properties were with the most significant differences (p-value < 0.01, highlighted by bold-underline).

Similar to the connectivity/adjacency-based network property, the calculated values of 16 properties based on the shortest path length in the human PPI network were provided in **Table 2** (corresponding definitions of these properties were provided in section Materials and Methods). As shown in TABLE 1 | The calculated values of 10 properties based on the connectivity and adjacency in the human PPI network.

The mean values (together with standard deviation) and median values of these properties between the targets of NTI and NNTI drugs were provided, and the statistical difference (p-value) for each property between targets of NTI and NNTI drugs were also calculated (p-values <0.05 and <0.01 were highlighted by bold and bold-underline, respectively).

TABLE 2 | The calculated values of 16 properties based on the shortest path length in human PPI network.

Mean values (together with standard deviation) and median values of these properties between the targets of NTI and NNTI drugs were provided, and the statistical difference (p-value) for each property between targets of NTI and NNTI drugs were also calculated (p-values <0.05 and <0.01 were highlighted by bold and bold-underline, respectively).

**Table 2**, all properties were found to be significantly different (p-values < 0.05, in bold font) between the targets of NTI and NNTI drug, and 14 (87.5%) of the 16 properties were with the most significant difference (p-value < 0.01, bold-underline). Moreover, the calculated values of 4 human biological system properties were shown in **Table 3** (definition of these properties was given in section Materials and Methods). As reported, these properties were frequently adopted to analyze the druggability of therapeutic targets for not only approved drugs but also the drugs in clinical trial development or withdrawn from market (Li et al., 2018). Herein, two properties were identified as significantly different (p-value < 0.01, bold-underline) between targets of NTI and NNTI drugs, which included the number of pathways affiliated by the targets of the studied drugs and the number of similarity proteins outside target's functional family. One thing needed to be emphasized was that the standard deviation of many properties was even larger than their mean value (such as bridging coefficient, clustering coefficient, and Z-score). These deviations indicated that the corresponding p-value may not be enough to measure the difference between the targets of NTI and NNTI drug. Moreover, any of the individual feature (p-value < 0.05 shown in **Tables 1**–**3**) could not be used to satisfactorily differentiate the targets of NTI drugs from that of the NNTI ones. Thus, this finding inspired us to discover the differential features using more advanced computational algorithm and collectively considering multiple properties.

TABLE 3 | The calculated values of four human biological system properties.

The mean values (together with standard deviation) and median values of these properties between the targets of NTI and NNTI drugs were provided, and the statistical difference (p-value) for each property between targets of NTI and NNTI drugs were also calculated (p-values <0.05 and <0.01 were highlighted by bold and bold-underline, respectively).

TABLE 4 | 19 substantially overlapped network properties grouped into 5 property groups based on their innate mutual dependence.

### Discovering the Key Features of NTI Drug Targets by Artificial Intelligence

Based on the in-depth investigation of 36 properties in **Tables 1**–**3**, several properties were found to be not fully independent or even duplicate in their descriptions (like degree vs. scaled degree). In this study, all 36 properties were systematically reviewed, and 19 of these 36 were identified to be substantially overlapped with some other properties (**Table 4**). Since there was significant dependence among the 19 properties, the use of all 36 properties for statistical feature selection may introduce strong biases. Thus, the 19 properties were grouped based on their innate mutual dependence. As shown in **Table 4**,

five property groups were generated by considering equation and description of these 19 properties, and each group was named by the first property (ordered alphabetically) in the corresponding group. As a result, these five groups included: the average closeness centrality, average shortest path length, betweenness centrality, degree, eccentricity. To minimize the possible bias induced by the innate mutual dependence among properties, only these five properties were considered in subsequent feature selection analysis, instead of investigating all 19 properties. Taking the remaining 17 relatively independent properties into consideration, 22 properties in total of each target were selected for subsequent feature selection.

As one of the most popular feature selection strategies based on AI, the Boruta algorithm based on a wrapper method built around a random forest classifier (Kursa, 2014) was adopted in this study. Boruta was considered the most powerful method with the stability in variable selection, especially suitable for the low-dimensional dataset among other reported strategies (Degenhardt et al., 2017). In this study, the key differential features were thus selected from 22 properties using R package Boruta by setting the p-value < 0.05. As a result, eight properties were selected as able to collectively reflect the target's mechanism underlying NTI drugs. As illustrated in **Figure 1**, the boxplots colored in red and green referred to the targets of NTI and NNTI drugs, respectively. Some key features increased from the targets of NTI drug to that of NNTI one (such as average shortest path length), while others demonstrated a decrease (such as average closeness centrality). Based on the comprehensive literature review, some of those 8 key features had been reported to be indirectly relevant to drugs' efficacy-safety balances. For example, the lower value of average closeness centrality of target was reported to demonstrate a less lethality risk (Chen et al., 2011), which was consistent with the findings of this study (a much higher average closeness centrality of the targets of NTI drugs was observed compared with that of NNTI ones, shown in **Figure 1**). Moreover, the higher level (lower value) of interconnectivity was frequently observed in lethal diseases such as cardiovascular disorder and cancer (Muhammd et al., 2018).

Oncological and cardiovascular disorder had been recognized as life-threatening diseases, and the majority of their drugs were reported to be NTI ones (Muller and Milton, 2012; Yu et al., 2015). Thus, the result of interconnectivity in **Figure 1** was consistent with these previous reports, which further validated the effectiveness of applied algorithm in identifying key target features underlying NTI drugs.

Moreover, there were four groups of targets as defined in section Materials and Methods: (a) targets of NTI drugs, (b) targets of both NTI and NNTI drugs, (c) targets of drugs without reported TI, and (d) targets of NNTI drugs. Apart from the target groups (a) and (d), the remaining groups provided more complicated and informative data for illustrating the mechanism underlying NTI drugs. On one hand, the targets in group (b) were affected by both NTI and NNTI drugs, which might reflect properties from both sides, but might also be significantly affected by the properties of confirmed NTI drugs. On the other hand, no TI data of the group (c) targets was reported based on literature review. It was possible that some NTI drugs were not discovered for those targets. But considering the large number of group (c) targets (339 in total), it was highly possible that most of those group (c) targets were only aimed by NNTI drugs, and just a small fraction of which could find new NTI drug in the future. The value of 8 properties of those 4 target groups were illustrated in **Figure 1**. It was interesting that all properties followed a clear descending/ascending trend from the targets of group (a) to (d), which was in accordance with the analyses provided above. Thus, these findings could be another line of evidence that validated the effectiveness of the feature identification algorithm applied in this study.

### Target Mechanism Underlying NTI Drugs Collectively Determined by Multiple Profiles

By collectively considering **Figure 1** and **Tables 1**–**3**, seven out of those eight selected key features showed significant difference (p-value < 0.05), but it was clear that these significant differences did not guarantee the corresponding feature as the key differential one (57.7% of the features with significant difference (p-value < 0.05) were not selected as key differential ones). Moreover, significant difference was not observed for the selected key feature bridging coefficient (p-value = 0.22). This finding indicated that those eight features collectively determined the target mechanism of NTI drugs, and the TI-related mechanism might be the result of the synergistical effects among those features. Moreover, the majority of these eight key features were identified for the first time by this study, and this work was also the first analysis on the collective effects of both PPI network properties and biological system profile on the drug efficacysafety balance.

Further analysis on these eight identified key features (shown in **Figure 1**) revealed that these key features were found to belong to three feature groups. These feature groups were connectivity and centrality of targets in human PPI network together with human biological system features. By combining the data in **Figure 1**, the key features within the same feature group (illustrated in **Figure 2**) followed the same ascending/descending trends, which were colored by the same background. As shown in **Figure 2**, the targets of NTI drugs were highly centralized and connected, and the number of similarity proteins and the number of affiliated pathways were substantially higher than those of NNTI drug. Since the number of similarity proteins and affiliated pathways was reported to be good indicator of target druggability (Zhu F. et al., 2009; Li et al., 2018), the NTI profile identified in this study was in accordance with that of reported target druggability.

### CONCLUSION

This work is the first study conducting comprehensive review on the TI data of all FDA approved drugs (**Supplementary Table S1**) and revealing the collective effects of both human PPI network properties and biological system profiles on drug efficacy-safety balance. Eight key features were identified here as collectively differentiating the target mechanisms between NTI and NNTI drugs. These features revealed that the targets of NTI drugs were highly centralized and connected in human PPI network, and the numbers of similarity proteins and target-affiliated pathways were both much higher than those of NNTI drugs. These

### REFERENCES

findings together with the newly discovered features/feature groups clarified the key factors indicating drug's narrow TI and could therefore provide a novel direction for determining the delicate drug efficacy-safety balance.

### AUTHOR CONTRIBUTIONS

FZ conceived the idea and supervised the work. XL, JY, and JT performed the research. XL, JY, JT, YL, QY, ZX, RZ, YW, JH, LT, and WX prepared and analyzed the data. FZ wrote the manuscript. All authors have read and approved this manuscript.

### FUNDING

This work was funded by National Natural Science Foundation of China (81872798); Innovation Project on Industrial Generic Key Technologies of Chongqing (cstc2015zdcy-ztzx120003); and Fundamental Research Funds for Central Universities (10611CDJXZ238826, CDJZR14468801, CDJKXB14011, 2015CDJXY).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar. 2018.01245/full#supplementary-material

of enzymes and cascade interactions of neighboring biological processes to identify drug-targets. Mol. Biosyst. 7, 1033–1041. doi: 10.1039/c0mb00249f

modulators in clinical trials by molecular dynamics simulations. ACS Chem. Neurosci. 9, 1492–1502. doi: 10.1021/acschemneuro.8b00059

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Li, Yin, Tang, Li, Yang, Xiao, Zhang, Wang, Hong, Tao, Xue and Zhu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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

Fen Pei † , Hongchun Li † , Bing Liu\* and Ivet Bahar\*

*Department of Computational and Systems Biology, School of Medicine, University of Pittsburgh, Pittsburgh, PA, United States*

#### Edited by:

*Zhi-Liang Ji, Xiamen University, China*

#### Reviewed by:

*Amit P. Bhavsar, University of Alberta, Canada Huanmei Wu, Indiana University, Purdue University Indianapolis, United States*

#### \*Correspondence:

*Bing Liu<EMAIL_ADDRESS>Ivet Bahar<EMAIL_ADDRESS>
*†These authors have contributed equally to this work*

#### Specialty section:

*This article was submitted to Translational Pharmacology, a section of the journal Frontiers in Pharmacology*

Received: *08 November 2018* Accepted: *14 February 2019* Published: *08 March 2019*

#### Citation:

*Pei F, Li H, Liu B and Bahar I (2019) Quantitative Systems Pharmacological Analysis of Drugs of Abuse Reveals the Pleiotropy of Their Targets and the Effector Role of mTORC1. Front. Pharmacol. 10:191. doi: 10.3389/fphar.2019.00191* 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.

Keywords: drug abuse, quantitative systems pharmacology, pleiotropic proteins, mTOR complex 1, drug-target interactions, neurotransmission, machine learning, cellular pathways

#### 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 (Koob and Volkow, 2016). 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.

## MATERIALS AND METHODS

#### 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**<sup>i</sup> , 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<sup>i</sup> . dj**) = (**d<sup>i</sup> . dj**)/(|**d<sup>i</sup>** | |**d<sup>j</sup>** |) between these vectors, and the corresponding cosine distance is [1–cos(**d<sup>i</sup> . dj**)]. Likewise, ligand-based distances between target pairs i and j were evaluated as the cosine distance between the 50-dimensional vectors **t<sup>i</sup>** and **t<sup>j</sup>** 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, 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

$$\mathcal{R}\_{\mathbf{N}\times\mathbf{M}} = U\_{\mathbf{N}\times\mathbf{D}}^T V\_{\mathbf{D}\times\mathbf{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<sup>0</sup> or more targets that belong to pathway A:

$$P^{A} = \sum\_{k\_0 \le k \le m} \frac{\binom{K}{k} \binom{M-K}{m-k}}{\binom{M}{m}}$$

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<sup>0</sup> 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:

#### pi <sup>∗</sup> = mink=i...m{min(pkm/i, 1)}

**Supplementary Table 4** lists these p-values for pathway enrichments based on both known and predicted targets.

The source code used for generating the results reported in this study is available at https://github.com/Fengithub/DA.

#### RESULTS

### Functional Similarity of Drugs of Abuse Does Not Imply Structural Similarity, Consistent With the Multiplicity of Their Actions

**Figure 2** presents a quantitative analysis of the functional and structural diversity of the examined n = 50 drugs of abuse, and the similarities among the m = 142 known targets of these addictive drugs. The n × n maps in **Figures 2A,B** display the drug-drug pairwise distances/dissimilarities based on their 2D fingerprints (**Figure 2A**), and their interaction patterns with their targets. **Figures 2C**,**D** display the corresponding dendrograms. The drugs are indexed and color-coded as in **Supplementary Table 1** and **Supplementary Figure 1**. As expected, drugs belonging to the same functional category (same color) exhibit more similar interaction patterns (**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).

labels in (C,D) are color-coded based on their categories: CNS stimulants (*green*), CNS depressants (*blue*), opioids (*cyan*), cannabinoids (*light brown*), anabolic steroids (*black*) and hallucinogens (*magenta*). Note that the drugs of abuse in the same category do not necessarily show structural similarities nor similar interaction pattern with targets. (E) Pairwise distance map for the 142 known targets based on their interaction patterns with the 50 drugs. The indices in (E) follows the same order as those listed clockwise in the dendrogram (F). The tree maps in (C,D,F) are generated based on the respective distances values in the (A,B,E).

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

FIGURE 3 | Promiscuity of drugs of abuse and their targets, and major families of proteins targeted by drugs of abuse. Number of known (*gray*) and predicted (*white*) interactions are shown by bars for (A) drugs of abuse and (B) their targets. The examined set consists of 50 drugs of abuse and a total of 142 known and 48 predicted targets, involved in 445 (known) and 161 (predicted) interactions. (A) Displays the number of interactions known or predicted for all 50 drugs. (B) Displays the results for the targets that interact with at least 4 known drugs (36 targets). The colors used for names of drugs and targets are same as those used in Figure 2. (C) Displays the distribution of families of proteins targeted by drugs of abuse.

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.  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

categories. See the complete list of pathways and targets in Supplementary Table 4.

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.

First, we note that eight pleiotropic proteins are shared by all five categories (at the intersection of the five Venn diagrams in **Figure 4B):** AMPA receptor (subtype GluA2; GRIA2), NMDA receptors 1 and 2A-D (designated as GRIN1, GRIN2A, GRIN2B, GRIN2C, and GRIN2D) and voltagedependent calcium channel Cav2.1 (or CACNA1A) as well as the predicted target phosphatidylinositol 3-kinase class 1A catalytic subunit α (PIK3CA) (**Supplementary Table 4**).

Second, 15 proteins are distinguished as targets of four of these major pathways: Serotonin receptors 5HTR2-A, -B and - C), GABA<sup>A</sup> receptors 1-6 (GABRA1- GABRA6), β-1 adrenergic receptor 1 (ADRB1), Ras-related C3 botulinum toxin substrate 1 (RAC1; member of Rho family of GTPases), mAChR M<sup>3</sup> (CHRM3) and DA receptor D<sup>2</sup> (DRD2), and two predicted targets - p38α (MAPK14) and DA receptor D<sup>1</sup> (DRD1).

AMPA receptor plays a crucial role in LTP and LTD, which are vital to neuroplasticity, memory and learning (Volkow et al., 2016). 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<sup>2</sup> 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 ∗ <sup>i</sup> = 2.01E-18), GABAergic synapse (p ∗ <sup>i</sup> = 1.19E-17), cholinergic synapse (p ∗ <sup>i</sup> = 2.36E-07), dopaminergic synapse (p ∗ <sup>i</sup> = 1.66E-06) and glutamatergic synapse (p ∗ <sup>i</sup> = 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**,

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 Ca2<sup>+</sup> signaling by inducing the production of inositol triphosphate (IP3) 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 1<sup>9</sup> tetrahydrocannabinol (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

(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.

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., CaV2.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<sup>3</sup> 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, 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 (Liu et al., 2016) and cell fate decision processes (Liu et al., 2014, 2017) 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.

#### REFERENCES

#### ACKNOWLEDGMENTS

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

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar. 2019.00191/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Pei, Li, Liu and Bahar. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Ontological and Non-Ontological Resources for Associating Medical Dictionary for Regulatory Activities Terms to SNOMED Clinical Terms With Semantic Properties

*Cédric Bousquet1,2\*, Julien Souvignet1,2, Éric Sadou1, Marie-Christine Jaulent1 and Gunnar Declerck3*

*Edited by: Lixia Yao, Mayo Clinic, United States*

#### *Reviewed by:*

*Yongqun Oliver He, University of Michigan Health System, United States Zhe He, Florida State University, United States*

*\*Correspondence:*

*Cédric Bousquet<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Translational Pharmacology, a section of the journal Frontiers in Pharmacology*

*Received: 18 July 2018 Accepted: 31 July 2019 Published: 10 September 2019*

#### *Citation:*

*Bousquet C, Souvignet J, Sadou É, Jaulent M-C and Declerck G (2019) Ontological and Non-Ontological Resources for Associating Medical Dictionary for Regulatory Activities Terms to SNOMED Clinical Terms With Semantic Properties. Front. Pharmacol. 10:975. doi: 10.3389/fphar.2019.00975*

*1 Laboratoire d'Informatique Médicale et d'Ingénierie des Connaissances en e-Santé, LIMICS, Sorbonne Université, Inserm, Université Paris 13, Paris, France, 2 Unit of Public Health and Medical Informatics, University of Saint Etienne, Saint Etienne, France, 3 EA 2223 Costech (Connaissance, Organisation et Systèmes Techniques), Centre de Recherche, Sorbonne Universités, Université de technologie de Compiègne, Compiègne, France*

Background: Formal definitions allow selecting terms (e.g., identifying all terms related to "Infectious disease" using the query "has causative agent organism") and terminological reasoning (e.g., "hepatitis B" is a "hepatitis" and is an "infectious disease"). However, the standard international terminology Medical Dictionary for Regulatory Activities (MedDRA) used for coding adverse drug reactions in pharmacovigilance databases does not beneficiate from such formal definitions. Our objective was to evaluate the potential of reuse of ontological and non-ontological resources for generating such definitions for MedDRA.

Methods: We developed several methods that collectively allow a semiautomatic semantic enrichment of MedDRA: 1) using MedDRA-to-SNOMED Clinical Terms (SNOMED CT) mappings (available in the Unified Medical Language System metathesaurus or other mapping resources, e.g., the MedDRA preferred term "hepatitis B" is associated to the SNOMED CT concept "type B viral hepatitis") to extract term definitions (e.g., "hepatitis B" is associated with the following properties: has finding site liver structure, has associated morphology inflammation morphology, and has causative agent hepatitis B virus); 2) using MedDRA labels and lexical/syntactic methods for automatic decomposition of complex MedDRA terms (e.g., the MedDRA systems organ class "blood and lymphatic system disorders" is decomposed in blood system disorders and lymphatic system disorders) or automatic suggestions of properties (e.g., the string "cyclic" in preferred term "cyclic neutropenia" leads to the property has clinical course cyclic).

Results: The Unified Medical Language System metathesaurus was the main ontological resource reusable for generating formal definitions for MedDRA terms. The non-ontological resources (another mapping resource provided by Nadkarni and Darer in 2010 and MedDRA labels) allowed defining few additional preferred terms. While the Ci4SeR tool helped the curator to define 1,935 terms by suggesting potential supplemental relations based on the parents' and siblings' semantic definition, defining manually all MedDRA terms remains expensive in time.

Discussion: Several ontological and non-ontological resources are available for associating MedDRA terms to SNOMED CT concepts with semantic properties, but providing manual definitions is still necessary. The ontology of adverse events is a possible alternative but does not cover all MedDRA terms either. Perspectives are to implement more efficient techniques to find more logical relations between SNOMED CT and MedDRA in an automated way.

Keywords: adverse drug reaction, Medical Dictionary for Regulatory Activities, SNOMED Clinical Terms, ontology, clinical terminology, pharmacovigilance

#### INTRODUCTION

Formal representation of semantics as provided by computational ontologies and associated semantic Web techniques have been extensively used in medical data integration systems in the last decade (Sheth et al., 2005), and they now tend to be acknowledged as a powerful means to improve the quality of the processing chain of medical data, process automatic extraction of information and knowledge from large databases or ensure semantic interoperability between disparate data processing systems (Park and Hardiker, 2009; Schriml et al., 2012; Schulz and Jansen, 2013).

In the medical domain, classic terminologies are gradually giving way to clinical terminologies, in which terms are defined using knowledge representation languages (Rossi Mori et al., 1998). An example is SNOMED Clinical Terms (SNOMED CT), a general clinical terminology whose objective is to represent all possible terms required for coding the patient record and other applications for representation of biomedical information (Khorrami et al., 2018). SNOMED CT presents several advantages compared with classic terminologies, especially the ability to apply techniques of semantic reasoning in order to build new groups of terms, whereas classic terminologies are limited to default groupings (generally made manually by experts) that are already specified as part of the terminology (Bousquet et al., 2005).

Medical Dictionary for Regulatory Activities (MedDRA) is a classic terminology used by regulatory authorities and pharmaceutical companies for coding adverse drug reactions (ADR) in pharmacovigilance databases (Brown et al., 1999). MedDRA terms are not formally defined and search is therefore limited to existing categories (Bousquet et al., 2005). It is frequently difficult to identify the exact MedDRA category that represents a given medical condition under investigation in a sufficiently specific and exhaustive way, for example, during a pharmacovigilance database search (Brown, 2003).

Since several years, we have performed studies that showed that a knowledge-based approach is efficient for building new groups of ADR terms with World Health Organization Adverse Reaction Terminology (WHO ART) (Alecu et al., 2006) (Iavindrasana et al., 2006) and with MedDRA (Henegar et al., 2006; Declerck et al., 2012; Asfari et al., 2016; Souvignet et al., 2016a) in an automated way. This means that starting from a resource containing formal definitions of ADR terms, it is possible to make queries that correspond to a case definition in order to retrieve the related set of terms. This strategy was applied in Pharmacovigilance Adverse Reaction Terminology Server (Alecu et al., 2007) where building a knowledge base for all WHO ART terms was a challenge. Indeed, all definitions were to be set manually, and we therefore focused on automated ways to enrich WHO ART (Iavindrasana et al., 2006). We found that mapping of WHO ART with SNOMED CT by means of the Unified Medical Language System (UMLS) metathesaurus proved to be a very efficient method to build formal definitions of WHO ART terms in an automated way (Alecu et al., 2008).

Difficulties we encountered for enriching WHO ART now appear at a larger scale in MedDRA due to a growing number of terms and a more complex organization of MedDRA. Indeed, only about 50% of MedDRA terms [excluding lowest level term (LLT)] were associated with a SNOMED CT concept in UMLS (Bodenreider, 2009). Therefore, the mapping method we applied to WHO ART was a fair starting point but proved to be insufficient for obtaining an exhaustive enrichment of MedDRA.

Our objective was to evaluate the potential of reuse of ontological and non-ontological resources for defining and/or enriching definitions of MedDRA terms. We present in this article

**Abbreviations:** ADR, Adverse drug reactions; AERS, Adverse Events Reporting System; Ci4SeR, Curation Interface for Semantic Resources; HLGT, Higher Level Group Term; HLT, High Level Term; ICD-10, International classification of diseases, 10th edition; LALR, Lexically assign logically refine; LLT, Lowest Level Term; LOINC, Logical Observation Identifiers Names & Codes; MedDRA, Medical dictionary for drug regulatory activities; NEC, Not elsewhere classified; NOS, Not otherwise specified; OAE, Ontology of Adverse Events; OWL, Web Ontology Language; PT, Preferred Term; SMQ, Standardized MedDRA Queries; SNOMED CT, SNOMED Clinical Terms; SOC, System Organ Class; UMLS, Unified Medical Language System; WHO ART, World Health Organization-Adverse Reaction Terminology.

several complementary methods that may benefit different levels of automation and could also be reused in order to semantically enrich other terminologies. These include methods such as i) extracting SNOMED CT definitions based on MedDRA-to-SNOMED CT mappings available in the UMLS metathesaurus or other mapping resources and ii) developing lexical and syntactic methods using MedDRA term label information. The ability to reuse the selected ontological and non-ontological resources was measured by comparing the number of MedDRA terms associated with a formal definition after processing of these resources with the number of MedDRA terms to define. Additionally, we manually curated some term definitions using expert knowledge that allowed us to evaluate the time necessary and validated a sample of the formal definitions provided by the previous methods. We stored formal definitions of MedDRA terms in a semantic resource named OntoADR (Bousquet et al., 2014).

The organization of OntoADR and results of semantic queries on OntoADR have been already published (Bousquet et al., 2014; Souvignet et al., 2016b). This article presents methods we implemented for reusing ontological and non-ontological resources to enable the formalization of the semantics and results obtained with each of these methods, how they automate the development of formal representations of MedDRA terms, the limits related to these methods, and additional developments that would be required for a more complete semantic enrichment of MedDRA.

#### BACKGROUND

#### Hierarchical Organization of Medical Dictionary for Regulatory Activities

The MedDRA hierarchy consists of five levels (from broad to narrow), among which four are depicted in **Figure 1**: System Organ Class (SOC), e.g., hepatobiliary disorders; higher level group term (HLGT), e.g., hepatic and hepatobiliary disorders; high level term (HLT), e.g., hepatic viral infections; preferred term (PT), e.g., hepatitis B; and LLT not shown on the figure. The PT level is preferred for data analysis and retrieval. MedDRA was defined as multi-axial because one PT may be present in one primary SOC and also in several secondary SOC. However, one PT may exist only within one single HLT within a SOC. As HLT within a SOC constitutes disjoint classes, it is seldom reliable to consider only one HLT or higher level category when searching for MedDRA terms related to a pharmacovigilance safety topic (Bousquet et al., 2005; Asfari et al., 2016).

Moreover, it was recognized that HLT are not always sufficient to represent clinical conditions involving several organs (e.g., anaphylactic shock involving the kidney, liver, cardiovascular, and respiratory systems) because they only group together terms belonging to the same SOC. When searching for signals associated with a drug, MedDRA terms representing the suspected ADR must thus be identified prior to the running of signal detection algorithms Souvignet et al. (2012). For instance, if one suspects a given drug to cause acute renal failure, using the MedDRA term "renal failure acute" is generally not sufficient for the algorithms to extract a signal because the acute renal failure condition can be coded with several related MedDRA terms by health professionals (e.g., "renal impairment," "blood creatinine abnormal," or "dialysis"). Identifying clinically related terms in MedDRA is not an easy task, as those terms might exist in different locations of the MedDRA hierarchy.

Since several years, the Maintenance and Support Services Organization, which is responsible for MedDRA maintenance and diffusion, builds standardized MedDRA queries (SMQ) to address these issues (Mozzicato, 2007). SMQs consist of sets of PT from different branches of MedDRA that allow describing a particular medical condition and are intended to aid in case identification. SMQs are a way to describe safety topics relevant

for pharmacovigilance that are not covered by the HLT and HLGT present in the MedDRA hierarchy. However, the achievement of SMQ raises important difficulties.

SMQs are currently developed manually by experts from the Maintenance and Support Services Organization, which is time consuming. Furthermore, once defined, they should not be modified or customized, and the existing SMQs do not cover all issues possible with drugs. Because experts have (even slightly) different understandings of the medical condition targeted by an SMQ, the kind of terms and the rationale for their selection may differ from an SMQ to another. This means that from one group of experts to another, for the same safety topic, the list of MedDRA terms selected could be different. Because SMQs are manually implemented, they could also miss important MedDRA terms. For those different reasons, the development of methods for automated selection of MedDRA terms on the basis of semantic information is desirable. Indeed, an automation of the process of PT selection in SMQ, even partial, could increase the quality and reproducibility of the SMQ and allow an important saving of time.

#### Difficulties for Searching Terms in Medical Dictionary for Regulatory Activities

The performance of pharmacovigilance systems based on spontaneous reporting is dependent on the information systems in which case reports are stored. In particular, these systems are subordinated to the ability of users to retrieve and exploit case reports in order to 1) reinforce existing knowledge on drug

safety, 2) make assumptions about the existence of a causal relationship between a drug and an adverse event, and 3) evaluate the available information to implement regulatory measures to secure drug therapies. The search for pharmacovigilance case reports is difficult because it is necessary to identify the medical terms indicating the safety topic that one wishes to evaluate. In general, a term is not sufficient to designate this safety topic, and it is preferable to look for all case reports in relation to a set of terms (Hauben et al., 2006; Hansen et al., 2007). According to the MedDRA® Data Retrieval and Presentation: Points To Consider (ICH Working Group, 2018), "clinically related PTs might be overlooked or not recognized as belonging together because they might be in different groupings within a single SOC or they may be located in more than one SOC."

**Figure 2** shows the problem of finding terms in MedDRA. Terms associated with a green tick are related to valvulopathy, while terms marked with a red cross do not correspond to valvulopathy. It is observed that the search terms are located in different branches of the terminology, which requires the pharmacovigilant specialist more time and effort to carry out his query. In addition, several HLT or HLGT must be combined to arrive at the final result, and irrelevant terms are present in these groups, which means that a search based on HLT or HLGT groupings will be associated with a large number of irrelevant PT. Another method for searching MedDRA terms is the textual query, but the terminology seems complex and does not reveal discriminating strings in the search for valvulopathies. For example "stenos \* aort \*" gives as results "stenosis of the aortic valve," "congenital aortic stenosis of the

valve," and "stenosis of the mitral valve and insufficiency of the aortic valve," but it is necessary to perform additional text searches corresponding to findings for valvular involvement such as calcification or insufficiency.

#### Interface, Aggregation, and Reference Terminologies

Schulz et al. (2017) recently evaluated how interface, aggregation, and reference terminologies may interact in the context of the new 11th version of the International Classification of Diseases and its relations with SNOMED CT. Interface terminology was defined by Rosenbloom et al. (2006) as a "systematic collection of health care–related phrases (terms) that supports clinicians' entry of patient-related information into computer programs [ … ]." According to Spackman (1997), the "main purpose of a reference terminology [ … ] is the retrieval and analysis of data." The term "aggregation terminology" was first introduced by Rogers (2005) to designate a classification systems in which its main purpose is to enable "statistical aggregation," further defined by Schulz et al. (2017) as consisting of single hierarchies and disjoint classes.

We consider that such clarification would be useful in our case study, to better explain what we intend to do with MedDRA and SNOMED CT. Within our approach, MedDRA is the aggregation terminology, and SNOMED CT is the reference terminology. As our purpose is to improve retrieval and not coding of MedDRA terms, we did not work on building an interface terminology. Such interface terminology would be desirable to facilitate data entry in pharmacovigilance databases but is outside of our scope. In the following paragraph, we show how a graphical user interface implementing SNOMED CT as a reference terminology could help users' experience in selecting MedDRA terms and potentially improving search in pharmacovigilance databases.

#### Rationale for Supplementing Medical Dictionary for Regulatory Activities With Formal Definitions

We consider it is possible to overcome the limitations associated with the organization of MedDRA terminology and the difficulty to identify related MedDRA terms, by proposing an alternative method for the grouping of PTs based on their medical meaning rather than their position in the hierarchy. This new method is based on PT modeling in a form that allows logical inferences by a computer. From a technical point of view, the implementation of this method is based on knowledge engineering, a branch of artificial intelligence in which it is possible to describe MedDRA terms using a formal language (Bousquet et al., 2014). In the field of knowledge engineering, we define "ontology" as the set of objects of a domain and relations between these objects.

While McKnight (1999) recognized the need for user-directed composition of controlled health terminologies, and the required improvement of the user interface in the context of data entry, we believe that such user interface is also of great importance to enable composition for data retrieval. In a previous work, we implemented OntoADR query tools, a graphical user interface that relies on OntoADR, and compared the performances of eight users in selecting MedDRA PTs with the MedDRA web browser in a pilot study on five medical conditions (Souvignet et al., 2019). Although the number of medical conditions was low, we observed a statistically significant improvement by using OntoADR query tools compared with the MedDRA web browser for selecting MedDRA PTs (+27% precision and +34% recall). Similar to Maedche and Staab (2001), we consider that the target application may serve as a measure for validating the implemented ontology and believe that such criteria is more important than criteria relying only on the evaluation of the ontology without taking into account the context where it is used. This pilot study confirmed the validity of our approach and justifies that we continue the implementation of our mappings between MedDRA and SNOMED CT.

In order to address a safety issue, it is necessary to identify case reports in pharmacovigilance databases relative to this issue. A safety issue may concern the causality assessment of drug D in the occurrence of medical condition C. **Figure 3** shows the example of three use cases where one is evaluating the causal role of suspected drug D in three medical conditions: a) upper gastrointestinal hemorrhage, b) medical conditions with symptom of erythema, and c) fungal infectious disorders. A single MedDRA term is usually not sufficient to characterize a given medical condition. OntoADR is intended to support selection of MedDRA terms according to different criteria. **Figure 4** depicts the query performed in OntoADR to retrieve MedDRA terms associated to these medical conditions, e.g., "finding site: upper gastrointestinal tract," and "associated morphology: hemorrhage," and the 10 first MedDRA terms retrieved by this query, e.g., anastomotic ulcer hemorrhage, aorto-esophageal fistula, chronic gastrointestinal bleeding, etc.

**Table 1** shows parts of formal definitions for 10 MedDRA terms associated with upper gastrointestinal hemorrhage among 27, in particular, SNOMED CT concepts that are filler of relations "finding site" and "associated morphology." Fillers that are not relevant to upper gastrointestinal hemorrhage are in italic, e.g., *Stenosis* for the MedDRA PT "anastomotic ulcer hemorrhage." In case the filler of the "associated morphology" relation is "hemorrhage," the query is immediately satisfied for this condition, e.g., three MedDRA terms (aorto-esophageal fistula, gastric antral vascular ectasia, and gastric hemorrhage) are defined as having "hemorrhage" as their associated morphology. When "hemorrhage" is not the filler of the "associated morphology" relation, other relevant fillers may be retrieved thanks to the subsumption mechanism that establishes hierarchical relations between a parent concept and its children concept. For example, "hemorrhage" subsumes "acute bleeding ulcer," 'bleeding varices," "chronic hemorrhage," and "hemorrhagic inflammation." "Upper gastro intestinal tract structure" subsumes "duodenal structure," "esophageal structure," "gastrojejunal junction structure," "pyloric antrum structure," and "stomach structure."

### Semantic Enrichment

The traditional process of domain ontology construction is based on expert intervention (Bedini and Nguyen, 2007). Although this manual procedure guarantees a fair quality of the generated

upper gastrointestinal hemorrhage, (B) MedDRA terms describing medical conditions associated with erythema, (C) MedDRA terms describing infectious diseases induced by fungi.

resource, it suffers from several difficulties; among those are the cold start problem (starting from scratch) and the lack of availability of domain experts (Qawasmeh et al., 2018). In fact, the high cost of experts' interventions is the major bottleneck identified early in the state of the art of ontology construction (Cullen and Bryman, 1988; Simperl et al., 2006; Balakrishna et al., 2010). This bottleneck justifies reusing and linking existing resources, when available, to create new ontologies (Alani, 2006). Reuse is not always possible because ontologies may not exist in the field of interest. For example, Mazo et al. (2017) describe a histological ontology of the human cardiovascular system and report that they "did not find in the State-of-the-Art an ontology of TABLE 1 | Finding site and associated morphology of 10 MedDRA terms describing upper gastrointestinal hemorrhage among 27.

histology neither a similar organization of hierarchies of histology terms that [they] may be able to reuse." At the opposite, when all the ontologies that are needed are available, it is sufficient to reuse and assemble them, such as in the example of the development of the orthology ontology (Fernández-Breis et al., 2016).

Ontology enrichment is the task of extending an existing ontology with additional concepts and semantic relations and placing them at the correct position in the ontology (Petasis et al., 2011). Automatic ontological construction is often based on learning (Maedche and Staab, 2001; Buitelaar et al., 2005). Such approach can be based on unstructured texts (Asim, 2018; Cimiano et al., 2006; Emani et al., 2015; Costa et al., 2016; Dasgupta et al., 2018), informal ontologies (Astrakhantsev and Turdakov, 2013), or linked data (Gavankar et al., 2012; Tiddi et al., 2012; Riga et al., 2017). A particular case of unstructured data corresponds to the labels of ontology identifiers that may be very dense with information. Such information described as "hidden semantics" by Third (2012) received little attention, at the exception of the gene ontology identifiers (e.g., Quesada-Martínez et al., 2015). SNOMED CT identifiers may also benefit from such approach, but this was limited to taking into account the "acute" and "chronic" qualifiers for evolution of diseases (Rector and Iannone, 2012) and the occurrence of congenital diseases (van Damme et al., 2018). Such "hidden semantics" were also detected by Nadkarni and Darer (2010) to build correspondences between MedDRA and SNOMED CT, but these correspondences were not associated with relations, which makes this work interesting for reuse but requires reengineering to transform the correspondences into semantic relations.

One of the major difficulties of these approaches is the extraction of non-taxonomic relationships (Dahab et al., 2008; Sánchez and Moreno, 2008; Villaverde et al., 2009; Petasis et al., 2011; Serra et al., 2014). Furthermore, several automatic approaches for ontology reusing and engineering still require domain experts and knowledge engineers (Bobed et al., 2012). Thus, semiautomatic approaches could be a good alternative (Balakrishna et al., 2010).

Semiautomatic approaches employ intelligent methods to significantly reduce, without completely replacing, human efforts (Huang et al., 2014). In such approaches, the role of experts could be limited to validating final automatic learning results (Wächter and Schroeder, 2010) or suggesting improvements at the end of ontology life cycle (Alobaidi, 2018). Expert intervention can be achieved with the help of graphical user interface (Wächter and Schroeder, 2010), spreadsheets (Blfgeh et al., 2017; Judkins et al., 2018), or specified pipelines such as the eXtensible ontology development (He et al., 2018).

#### NeON Methodology

After comparing several methods for ontology development, we selected the NeON methodology (Suárez-Figueroa et al., 2012) because it was the most appropriate to illustrate the strategy that we followed for designing OntoADR. While other methodologies may also be considered for knowledge engineering and may be more relevant in other contexts, we considered that dimensions of reuse were the most important features when selecting NeON.

While other approaches for ontology engineering provide methodological guidance, the NeON Methodology does not prescribe a rigid workflow. It suggests a variety of pathways based on nine flexible scenarios that address common issues, such as reusing, reengineering, and merging ontological resources. These ontological resources also comprise ontology design patterns (Aranguren, 2008; Gangemi, 2005; Blomqvist, 2008; Presutti and Gangemi, 2008), which are generic templates or abstract descriptions proposed to enforce best practices in ontology implementation. One particularity of NeON matching well with our specific approach is that it also takes into account reusing and reengineering of non-ontological resources, which is not the case of other methodologies such as METHONTOLOGY (Fernández-López et al., 1997) and On-To-Knowledge (Sure et al., 2004). These non-ontological resources may consist of structured data such as terminologies (Jimeno-Yepes et al., 2009) or databases, unstructured data (e.g., articles), or semi-structured data (e.g., XML, JSON) (Qawasmeh et al., 2018). In addition to these nine scenarios, NeON also integrates support activities such as knowledge acquisition, documentation, and evaluation that should be carried out during the whole ontology development cycle.

### MATERIAL AND METHODS

#### Summary of the Method Application of the NeON methodology

We applied the following scenarios of the NeON methodology for implementing OntoADR (**Figure 5**).

• Scenario 1. From Specification to Implementation: this includes four steps (specification, conceptualization, formalization, and implementation). We previously presented these steps in previous work [(Bousquet et al., 2014)] and limit here the

scope of this presentation to scenarios that emphasize the reuse of ontological and non-ontological resources.

In addition to these different scenarios, we also implemented "ontology support activities" for "knowledge acquisition" that comprises activities for (1) capturing knowledge from the MedDRA labels and work by a domain expert for adding formal definition using the Ci4SeR tool [(Souvignet et al., 2014)] and (2) "ontology validation" that consists in checking that the meaning of the ontology definitions are compliant with the definitions we intended the MedDRA terms to convey.

#### Flow chart of the method

**Figure 6** depicts a flow chart representing an overall representation of the several steps and tasks proposed in the article to get an overview of the algorithm at a glance. While this diagram could make readers believe that all these steps were conducted in parallel, it is proposed only as a convenient way to apprehend the method as a whole. The previous paragraph where we applied the NeON methodology shows a different perspective where different scenarios were applied at different time. In the flow chart, each MedDRA term is considered one after the other and can go through several parallel paths according to different conditions.

Some manual definitions can optionally be added (see the section *NeON Methodology*). All partial definitions acquired with the different algorithms are then automatically combined into a merged definition of the MedDRA term. We used MedDRA version 17 that consists of 26 SOC, 334 HLGT, 1,720 HLT, 20.559 PT, and 72,637 LLT, SNOMED CT version March 2015 and UMLS version 2014AB. SNOMED CT concepts were extracted from the Concepts\_Core\_INT file (Release Format 1) and the hierarchy and semantic properties from the Relationships\_Core\_INT file. This version of MedDRA was applied to the following paths: "using UMLS metathesaurus mappings," "automatic enrichment methods," and "manual definition of concepts." MedDRA 13 that consists of 26 SOC, 335 HLGT, 1,709 HLT, 18,786 PT, and 68,258 LLT was applied to the following paths: "Using other mapping resources" and "Using a decomposition algorithm and Metamap software to map complex MedDRA terms."

#### Problems With Mapping Other Layers Than the PT Level

We have tried to map other layers (SOC, HLT, and HLGT). For instance, the cardiac disorders SOC concept has for formal definition hasFindingSite some "Heart Structure." While

some SOC, HLGT, and HLT were accurately defined, essentially thanks to mappings in UMLS, it was decided not to present the results in this article. We used formal definitions associated with the three MedDRA higher levels only in the Ci4SeR tool (Souvignet et al., 2014), where the curator could use them if she considered them as relevant. We decided not to map the LLT because PT is the preferred level for case report analysis and search in pharmacovigilance databases.

We explained in previous work why the MedDRA hierarchy cannot be converted into a subsumption tree because this sometimes causes semantic inconsistencies (Bousquet et al., 2014). The reason is that most high level categories in MedDRA are intended to reflect the domain actors' practices (i.e., following the different medical specialties) and are not necessarily organized according to different semantic criteria as one would expect in a well-formed ontology. A first example concerns groups of symptoms (HLGT or HLT) that are placed under the general categories of disorders that they are the symptom of (SOC or HLGT). Such hierarchical organization would not be authorized in an ontology as the relation being-a-symptom-of does not imply an is-a relation. For instance, the PT "dyspnoea" and "dizziness or syncope" belong to the HLGT "cardiac disorder signs and symptoms" that is under the SOC cardiac disorders. While dyspnea for instance refers to conditions that may be associated to cardiac disorders, such symptom cannot be considered as a cardiac disorder. A second example is the MedDRA PT sudden death that belongs to the following hierarchies: 1) ventricular arrhythmias and cardiac arrest (HLT)/cardiac arrhythmias (HLGT), and 2) death and sudden death (HLT)/fatal outcomes (HLGT). While "cardiac arrhythmias" is defined in OntoADR with the hasFindingSite some "Heart Structure" property, the sudden death PT should not inherit from such property because sudden death could be the consequence of death that is not of cardiac origin.

#### Using Medical Dictionary for Regulatory Activities-to-SNOMED Clinical Terms Mappings From UMLS Metathesaurus

The UMLS may be used as a source of knowledge for adding formal definitions to medical terminologies (Schulz and Hahn, 2001). Based on our initial experience (Alecu et al., 2008), we assume that SNOMED CT is currently the best candidate for providing formal definitions to MedDRA. SNOMED CT terms are defined using description logic (DL) formalism, and a fair number of alignments between MedDRA and SNOMED CT are present in the UMLS metathesaurus. Therefore, most of the formal definitions attributed to MedDRA terms in OntoADR are based on semantic information extracted from the SNOMED CT clinical terminology. When a MedDRA term is mapped to a SNOMED CT concept, we reused the semantic information within SNOMED CT in order to build the formal definition of the MedDRA term. Identifying reliable MedDRA-to-SNOMED CT mappings is thus an essential step in our methodology to define MedDRA term semantics.

The UMLS (Lindberg et al., 1993) consists of a semantic network and a metathesaurus developed by the US National Library of Medicine to link terms from more than a hundred controlled vocabularies, including SNOMED CT and MedDRA. Terms from the different vocabularies are linked together by association to a unique UMLS concept defined by a concept unique identifier, e.g., "C0019163" that is mapped to both MedDRA term "hepatitis B" and SNOMED CT concept "type B viral hepatitis."

In OntoADR, MedDRA term "hepatitis B" has the formal definition: hasFindingSite some "Liver Structure," hasAssociatedMorphology some "Inflammation Morphology," and hasCausativeAgent some "Hepatitis B Virus" where has FindingSite, hasAssociatedMorphology, and hasCausativeAgent are OntoADR semantic relations inspired from SNOMED CT, and "liver structure," "inflammation morphology," and "hepatitis B virus" are SNOMED CT concepts we imported in OntoADR (**Figure 1**).

In order to map MedDRA and SNOMED CT terms, we developed an algorithm following these steps: i) search for the MedDRA PT in the UMLS using the MedDRA identifier; ii) for MedDRA PT without SNOMED CT mappings in the UMLS, if the PT has one or more related LLT considered as synonymous, then the LLT identifier is used for a new UMLS search. iii) If this second search is unsuccessful, the algorithm performs a last UMLS search using PT and LLT labels, seeking to pair these labels with SNOMED CT concepts by string matching.

All mapping propositions selected from the UMLS metathesaurus were validated, modified, or completed by knowledge engineers and pharmacovigilance experts of our team. i) All one-to-one mappings we decided to use were first validated by checking manually the correspondence between the meanings of terms.

Several SNOMED CT concepts may be proposed as synonyms (Fung et al., 2005) of the same MedDRA concept, although they have different meanings. Each MedDRA concept in OntoADR can have only one equivalent SNOMED CT concept. When several SNOMED CT concepts are proposed in UMLS as synonym of a MedDRA term, only one was selected by an expert for building the formal definition. Such selection should be based on synonymy between a SNOMED CT concept and a MedDRA term. According to Fung, synonymy between term X and Y may be defined according to linguistic criteria (Fung et al., 2005) such as enforcing that it is possible to replace X by Y in any sentence without modifying the meaning. Examples of such synonyms are "celiac disease" and "gluten enteropathy," or "kidney stone" and "renal calculus." Selection of a SNOMED CT concept was performed first by comparing its label with the MedDRA term label. In case both were identical, which occurred most of the time, it was obvious to select this mapping, but in other cases, we took into account the medical relevance of the mapping and had to rely on expert evaluation. For example, three SNOMED CT concepts are mapped to the MedDRA term "Spondylitis": "inflammatory spondylopathy," "undifferentiated spondylitis," and "spondylitis," and the later appears as a perfect match according to label comparison.

Such a validation process was necessary because UMLS mapping propositions are not always semantically valid. MedDRA terms and SNOMED CT concepts mapped together in UMLS can refer to different medical entities, even if they are homonyms. For instance, the MedDRA term "vascular disorders" and its SNOMED CT homonym "vascular disorder" are mapped together in UMLS; however, the former refers to disorders of blood and lymphatic systems and the later only to disorders of blood vessels (lymphatic system disorders are caught by the concept "disorder of lymphatic system" in SNOMED CT). ii) In case of one-to-n mappings, a manual expert choice was made to select the SNOMED CT concept whose definition best fitted the meaning of the correspondent MedDRA concept. iii) When no SNOMED CT concept among the ones suggested by UMLS was satisfactory, the definition of the correspondent MedDRA concept was made manually. iv) Mapping a MedDRA term with a SNOMED CT concept does not ensure that the former gets a complete (or even a satisfying) formal definition of its semantics: the formal definition can be incomplete or even be literally absent: it is common to find SNOMED CT concepts, for instance psychiatric concepts, that have no definitional properties. When necessary, the semantic properties from SNOMED CT attributed to MedDRA concepts through the mapping procedure were thus completed manually by additional assertions.

#### Using Another Medical Dictionary for Regulatory Activities-to-SNOMED Clinical Terms Mapping Resources

To complete the mappings selected from UMLS, we also made use of Nadkarni and Darer's propositions of mappings (Nadkarni and Darer, 2010). Using one year of data (recorded between July 1, 2008 and April 30, 2009) from the US Food and Drug Administration Adverse Events Reporting System (AERS) pharmacovigilance database, the authors identified 3,705 MedDRA PT that collectively accounted for 95% of case reports. The 3,705 selected MedDRA terms correspond to high-frequency terms in the US Food and Drug Administration database and potentially have a great added value. After eliminating terms already mapped to SNOMED CT concepts in UMLS, they attempted to map manually the remaining terms (786 in total) with software assistance. Most of those terms (733) could be mapped by Nadkarni and Darer with SNOMED CT concepts *via*  one-to-one or one-to-n mappings.

Several problems have been encountered when trying to reuse Nadkarni and Darer's propositions of mappings (Nadkarni and Darer, 2010). i) First, in the case of one-to-n mappings, the authors broke down a MedDRA term in such a way to associate it to several SNOMED CT concepts but did not specify which semantic relation was relevant. For example, they mapped the MedDRA concept "tongue discoloration" with SNOMED CT concepts "abnormal color" and "entire tongue" but did not specify which semantic relation interconnected the first to the others. Obviously, it cannot be here an equivalence or synonymy (Fung et al., 2005) ("same as") relation as in the case of one-to-one mapping. When SNOMED CT concepts belong to branches such as *body structure* or *morphologic abnormality*, the relationship to use is easy to deduce, and its creation can be automated: it will be in the first case the hasFindingSite relationship and, in the second, the hasAssociatedMorphology relationship. However, when it comes to SNOMED CT concepts from branches *finding*, *qualifier value*, or *disorder*, the relationship to use is not obvious, and only a human expert can decide. A major part of our recovery work was to specify these relationships by making use of the set of relationships available in OntoADR.

ii) We have also occasionally been forced to revise the proposed Nadkarni and Darer's mappings, partly for reasons of pure semantic accuracy and partly because of the purpose of OntoADR, which we illustrate here using three examples:

(combined site)" is set in an isolated portion of the SNOMED CT branch "body structure" (branch called "group of anatomical entities"). Nothing connects it to the concepts of the digestive system structures (e.g., no relationship part-of). Due to the SEP decomposition (structure, entire, part) of the anatomical branch of SNOMED CT, the concept "anal structure" has no relation to the concept "anus and rectum." It would have been impossible to use this localization by semantic reasoning, for example, to identify concepts located on part of the anorectal system (principle of subsumption reasoning: concepts that have a relationship of location on parts of the anorectal structure are considered by inference as siblings of the concept of diseases that are located on the whole anorectal structure). We therefore preferred to use the SNOMED CT concept "anorectal structure" to define the relationship hasFindingSite of the MedDRA concept "anorectal discomfort" in OntoADR. This SNOMED CT concept allows the semantic reasoning operation described previously. Moreover, we can assume that the MedDRA term "anorectal discomfort" is sometimes used to encode ADRs in a non-specified way that may be anal or rectal, and not both, as is implied by the use of the SNOMED CT concept "anus and rectum (combined site)." It is therefore important to locate by subsumption concepts that are located within a substructure of the whole anorectal structure.

#### Using a Syntactic Decomposition Algorithm on Complex Medical Dictionary for Regulatory Activities Terms

Among MedDRA terms that are not mapped with SNOMED CT terms in UMLS, there are many complex terms, i.e., corresponding to composed expressions. MedDRA complex terms are of several kinds: a) expressions composed with an AND logical operator or commas (e.g., "acute and chronic thyroiditis" or "pregnancy, labour, delivery, and postpartum conditions"); b) expressions composed with "NEC" (not elsewhere classified), "unspecified," or with a text between brackets, usually to specify exclusion clauses [e.g., "autoimmune disorders NEC," "laryngeal neoplasms malignancy unspecified," or "ocular neoplasms malignant (excl. melanomas)"]; c) they can also combine these different kinds of complexity [e.g., "gastrointestinal and abdominal pains (excl. oral and throat)" or "ocular structural change, deposit, and degeneration NEC"]. These terms are usually terms of level HLT, HLGT, and SOC in the MedDRA hierarchy. However, their definitions have a great added value because some terms they subsume may inherit their properties. Indeed, defining one high level term with a morphology property may amount to defining all child terms with this property within the limits of what we have indicated in the section *Problems With Mapping Other Layers Than the PT Level*.

The complex MedDRA terms present two kinds of difficulties: i) the difficulty of mapping with SNOMED CT that tends to favor simple concepts probably because most complex concepts correspond to pure classifying artifacts, e.g., "not elsewhere classified," without real counterpart in the phenomena that are part of medicine; ii) difficulties for formalization of meaning: representing in OWL (Web Ontology Language) the meaning of a compound concept containing exclusions with logical operators is constrained by the expressiveness of the DL language used. In OntoADR, it is not possible to describe the exact same MedDRA semantics due to computability constraints. To date, we have developed a technical solution for the first point, but no satisfactory solution of conceptualization (especially in terms of human cost modeling) has yet been developed to meet the second point. It should be noted that this issue is regarding mainly terms of high levels and does not affect the progress of definitions for PT terms in OntoADR.

In order to map complex MedDRA terms, we developed an algorithm for syntactic decomposition. It consists of three routines: 1) a routine for "cleaning" terms; 2) a routine for identification of complex expressions; and 3) a routine for decomposition of an expression from a set of formal rules. Routine 1 begins by suppressing from the MedDRA labels unnecessary characters or characters that cannot be supported by the decomposition routine [stop words, content between brackets, terms as "unspecified," "NOS" (not otherwise specified), etc.]. Routine 2 identifies decomposable expressions: it searches for keywords that indicate a probable composition of the expression ("AND," "OR," "WITH," ",", etc.). Finally, routine 3 decomposes the complex expression in a set of simpler expressions (cf. **Table 2**), by applying different rules, for example:

$$\begin{aligned} \text{(A AND B).q} &\rightarrow \text{ A.q} + \text{B.q} \\ \text{q.(A AND B)} &\rightarrow \text{ q.A} + \text{q.B} \end{aligned}$$

We then used the MetaMap software (Aronson, 2001) to map all new decomposed concepts to existing SNOMED CT concepts.

#### Automatic Lexical Enrichment Methods

We have used a rule-based algorithm for automatic suggestion of properties from the MedDRA label to enrich the formal definition of concepts. Two key procedures have been implemented in the algorithm:

1. When the algorithm detects a given string Sx in a MedDRA label, it automatically adds a corresponding property Px in the OWL concept definition. For example: if the string "pain" or "algia" is found in a MedDRA concept's label, the semantic property hasDefinitionalManifestation some Pain is automatically added to the concept's definition. Similarly, if the string "perforation" is found, the formal definition hasAssociatedMorphology some Perforation is suggested. All created properties are then validated by an expert. Illegitimate properties are rejected. For example, the algorithm proposed to add the formal property hasAssociatedMorphology some Hernia to the MedDRA concept "hernia repair," as the string "hernia" was found in the label. Semantically, this assignment is obviously illegitimate: the "hernia repair" is not a type of hernia and cannot be defined by this morphological property. The property has therefore been rejected. The expert, however, took advantage of this suggestion to correct it in: OccursAfter TABLE 2 | Examples of parsing of complex MedDRA terms using different rules.

some Hernia. This validation step is also necessary due to the occasionally polysemic expressions used for automatic generation of properties. For example, the automatic generation of the hasClinicalCourse some Cyclic property, when the algorithm detects the "cyclic" string in a MedDRA label is valid for concepts such as "cyclic neutropenia" or "cyclic vomiting syndrome," where the term "cyclic" indicates the clinical course of the disease. However, it is not valid for the concept "cyclic AMP," which refers to a clinical test (a measure of the presence or amount of cyclic adenosine monophosphate, e.g., in urine). We could have improved the automatic processing in order to detect these problematic cases, but the formalization of these exceptions would have taken longer time than using a manual approval process.

A restriction is applied to prevent the property Px to be duplicated when it already exists in a MedDRA term definition, e.g., the detection of the "perforation" string in the label of a MedDRA concept CMed only results in the creation of the property hasAssociatedMorphology some Perforation if CMed does not already own a property hasAssociatedMorphology some <Morphology>. If it is the case, we assume that the relation Rx (in this case hasAssociatedMorphology) has already been filled in correctly.

2. A second procedure is implemented by the algorithm to automatically generate properties. Based on the same principles, but working with more complex patterns of recognition, it was designed to complete definitions of MedDRA concepts referring to investigations and their results (SOC « Investigations »).