Article pubs.acs.org/crt
Drug-Induced Rhabdomyolysis: From Systems Pharmacology Analysis to Biochemical Flux Junguk Hur,†,‡,⊥ Zhichao Liu,§,⊥ Weida Tong,§ Reijo Laaksonen,∥ and Jane P. F. Bai*,‡ †
Department of Neurology, University of Michigan, Ann Arbor, Michigan 48109, United States Office of Clinical Pharmacology, Center for Drug Evaluation and Research, U.S. Food and Drug Administration, Silver Spring, Maryland 20993, United States § Divisions of Bioinformatics and Biostatistics, National Center for Toxicological Research, U.S. Food and Drug Administration, Jefferson, Arkansas 72079, United States ∥ The School of Medicine, University of Tampere, Tampere 33520, Finland ‡
S Supporting Information *
ABSTRACT: The goal of this study was to integrate systems pharmacology and biochemical flux to delineate drug-induced rhabdomyolysis by leveraging prior knowledge and publicly accessible data. A list of 211 rhabdomyolysis-inducing drugs (RIDs) was compiled and curated from multiple sources. Extended pharmacological network analysis revealed that the intermediators directly interacting with the pharmacological targets of RIDs were significantly enriched with functions such as regulation of cell cycle, apoptosis, and ubiquitin-mediated proteolysis. A total of 78 intermediators were shown to be significantly connected to at least five RIDs, including estrogen receptor 1 (ESR1), synuclein gamma (SNCG), and janus kinase 2 (JAK2). Transcriptomic analysis of RIDs profiled in Connectivity Map on the global scale revealed that multiple pathways are perturbed by RIDs, including ErbB signaling and lipid metabolism pathways, and that carnitine palmitoyl transferase 2 (CPT2) was in the top 1 percent of the most differentially perturbed genes. CPT2 was downregulated by nine drugs that perturbed the genes significantly enriched in oxidative phosphorylation and energy-metabolism pathways. With statins as the use case, biochemical pathway analysis on the local scale implicated a role for CPT2 in statin-induced perturbation of energy homeostasis, which is in agreement with reports of statin−CPT2 interaction. Considering the complexity of human biology, an integrative multiple-approach analysis composed of a biochemical flux network, pharmacological on- and off-target networks, and transcriptomic signature is important for understanding drug safety and for providing insight into clinical gene−drug interactions.
■
INTRODUCTION
cascade effect of disrupting intracellular physiological homeostasis and leading to apoptosis or necrosis. Considering the complex cross-talk among intracellular pathways and the accumulated knowledge bases as well as databases of systems biology, the extended pharmacological network of rhabdomyolysis-induced drugs (RIDs) can provide insight into the pathways related to drug-induced rhabdomyolysis. The extended pharmacological network is defined as consisting of drugs, their targets, and proteins or genes that interact, under a prespecified statistical criterion, with the drugs via their targets. To understand serious ADR retrospectively, individual laboratories in the scientific community conducted studies and published results, typically with a focus on a specific component of the systems of human biology and physiology.
Drug-induced rhabdomyolysis, although rare, is a serious adverse reaction and can be fatal. Cerivastatin was withdrawn from the market because of treatment-related rhabdomyolysis associated with its use.1 Myopathy and myalgia induced by drugs may be related to drug-induced rhabdomyolysis, although they are not as serious as rhabdomyolysis and may result in a dose reduction or switching of medications. Prospective prediction of adverse drug reactions (ADRs) at the patient level has been very challenging. A recent study by Vilar et al.2 compiled 196 drugs across a wide range of pharmacologic classes. In addition to the structural similarities investigated by Vilar et al., the biological pathways perturbed by these 196 drugs and their pharmacologic networks could overlap to varying degrees with parts of the perturbed biological networks; however, they did so completely independent of one another. In addition to engaging specific pharmacological target(s), each of these drugs can interact with off-targets, resulting in a © 2014 American Chemical Society
Special Issue: Systems Toxicology Received: November 5, 2013 Published: January 14, 2014 421
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
■
For example, there are reports that statins reduce in vitro ATP production3 and perturb several functional pathways in human muscles and change the lipid profiles in circulation;4 there are also case reports of carnitine palmitoyl transferase 2 (CPT2) deficiency carriers developing treatment-related rhabdomyolysis while receiving statins.5 Clearly, there is a need to change the current paradigm of scientific investigation of ADRs to achieve the goal of predicting drug safety. There is also a critical need to leverage existing knowledge and curated data for an integrative analysis for the purpose of developing new drugs with a better effectiveness/risk index. It is our goal to develop predictive models for drug-induced adverse reactions, including drug-induced rhabdomyolysis. We hypothesize that some gene−drug interactions that underpin specific ADRs may not be detectable by genome-wide association studies but could be understood and predicted using systems pharmacology analysis. To establish a predictive framework and to lay the groundwork for achieving future success in developing predictive models for ADRs, we propose that both the extended pharmacological network of RIDs, including their on- and off-target effects, and the genes perturbed by RIDs, as shown by transcriptomic analysis, could be linked to the biochemical pathways that are essential for intracellular energy homeostasis. We first constructed the extended pharmacological base network by utilizing the existing knowledge of drug targets (such as http://www.drugbank.ca/) and then expanded it by leveraging our knowledge of protein− protein and genetic interactions (http://thebiogrid.org/). We propose to use this computational exercise to define, as broadly as possible, the scope of the biological plausibility for druginduced rhabdomyolysis. Then we utilized the published transcripomic database of Connectivity Map6 to delineate the functional and biochemical pathways that are perturbed by RIDs. We further compared these results from cancer cell lines and extended pharmacological network analysis with the transcriptomic profiles of human muscles that were obtained from subjects receiving statins4 to determine and understand how these different approaches interact among one another in the context of understanding the toxicity of drugs at the molecular level. On the basis of our preliminary transcriptomic analysis of RIDs, oxidative phosphorylation was identified as a key pathway perturbed. Because statins are widely used for preventive care of cardiovascular diseases, their interactions with CPT2 via biochemical pathways were used as a use case for gaining insight into drug-induced rhabdomyolysis. The SLCO1B1 gene, encoding the organic anion-transporting polypeptide (OATP1B1) that is involved in hepatic uptake of stains, is associated with statin-induced myopathy.7 In addition to the genes involved in absorption, distribution, metabolism, and elimination (ADME), there are non-ADME genes reportedly associated with drug-induced rhabdomyolysis8 and statin myopathy.5 For example, CPT2 deficiency resulting from mutations in the CPT2 gene9 and gene−drug interactions between CPT2 mutation and drug action8 cause genetic rhabdomyolysis and drug-induced rhabdomyolysis, respectively. To demonstrate further the impact of a specific perturbed pathway on cellular physiological homeostasis, we adopted a qualitative balance analysis of ATP flux that is produced through the biochemical pathway/network involving CPT2 and the target of statins. This study demonstrates the advantage of utilizing three different approaches for an integrative systemspharmacology understanding of drug-induced rhabdomyolysis.
Article
EXPERIMENTAL PROCEDURES
Compilation and Curation of the Drug List. The list of RIDs was compiled and curated in multiple steps, as illustrated in Figure 1.
Figure 1. Overall workflow of our study, consisting of construction of pharmacological networks, transcriptomic analysis, and biochemical flux balance analysis. To be consistent across different sets, we limited our search term specifically to rhabdomyolysis. A text-mining approach was developed to screen automatically the large volume of drug labels and packaging. Briefly, the most recent drug labels (in PDF file format) of each approved drug were retrieved from the Drugs@FDA database (http:// www.accessdata.fda.gov/scripts/cder/drugsatfda/). The PDF files of drug labels were converted to text files using Xpdf (http://www. foolabs.com/xpdf/). Text-converted drug labels from DrugBank and XML-formatted drug information from DailyMed (http://dailymed. nlm.nih.gov/) were analyzed to identify any drug names and active ingredient information for the list of drugs for which the term rhabdomyolysis was included in the adverse event section. The flanking text surrounding the rhabdomyolysis term was manually checked to make sure that rhabdomyolysis was indeed a treatmentrelated adverse event. Another list was compiled from the SIDER2 database (http://sideeffects.embl.de/), a computer-readable drug and human side-effects relationship data source. The SIDER2 database consists of information from different version of drug labels and health-related public documents. We also collected RIDs from a published paper by Vilar et al.2 and a curated list from the FDA (http://www.fda.gov/). To reduce redundancy and to ensure consistency, we normalized drug names. The official drug names used in DrugBank were adopted. Any drug with no corresponding record in the DrugBank database was processed by a Perl script to identify automatically the active pharmaceutical ingredient rather than generic drug names. For example, Olmesartan Medoxomil was normalized to Olmesartan. The list of drugs obtained from aforementioned resources were normalized, combined, and carefully reviewed. Of note, when compiling the list of RIDs, we did not differentiate the drugs with incidence rates reported during phase 3 clinical trials from those with postmarketing reports only. Only the drug products of a single active pharmaceutical ingredient were compiled. To curate the final list, we further referred to the clinical criteria of rhabdomyolysis, which is defined as muscle symptoms plus significant creatinine kinase elevation (usually >40 times the upper limit of the normal range (ULN)) and creatinine elevation (with brown urine/urine myoglobin).10,11 Extended Pharmacological Network (On Targets Plus Their Interacting Proteins) Analysis. A base pharmacological network, consisting of RIDs and their known drug targets, was constructed by referencing the drug databases, including DrugBank, Therapeutic Target Database (http://bidd.nus.edu.sg/group/cjttd/) (TTD), and PharmGKB (http://www.pharmgkb.org/). Drugs with any known targets were included in the base network. To identify common but 422
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
Article
Biochemical Flux Analysis. Qualitative biochemical flux balance was analyzed to complement the pharmacological targets- and transcriptomics-based network analyses for the use case of statininduced rhabdomyolysis. CPT2 deficiency or mutation has been reported in statin-induced myopathy or rhabdomyolysis.5 Statins were used to illustrate gene−drug interaction, with the focus on the CPT2 gene encoding carnitine palmitoyl transferase 2 (CPT2). In light of the biological complexity, we focused on the linear paths from fatty acid oxidation to ATP production involving CPT2 by ignoring any feedback loops as well as on the biosynthesis of prenylated Rab and coenzyme Q10 (CoQ10) (also known as ubiquinone),17 which are reportedly related to statin-induced myopathy or rhabdomyolysis.18 We did this by focusing on ATP production from fatty acid oxidation (ATP is produced from the tricarboxylic acid cycle) and from the biological processes involving CoQ1017 as well as prenylated Rab1.18 For the purpose of flux analysis, the immediate paths forward from fatty acid oxidation are designated as CPT2 and non-CPT2 involved (Figure 2). Throughout the text, the flux denoted as J on the left and
potentially crucial regulator or effector in the pathogenesis of druginduced rhabdomyolysis, we extended the base network by adding genes that directly interact with the drug targets via either protein− protein interactions (PPIs) or genetic interactions (GIs). These additional genes, referred to as intermediator in the current study, are not a direct target of these drugs; however, they may play important roles in the pathogenesis of drug-induced rhabdomyolysis by interacting with the targets of RIDs. The PPIs and GIs were derived from the Biological General Repository for Interaction Data sets (BioGRID; release 3.2.97). Ten extended networks, denoted as Min01−Min10, were generated using the number of interacting drugs from 1−10 as the criterion for intermediators to be included in the extended networks. The Min01 network included any intermediators with known interactions with the targets of any RIDs, whereas Min10 included only the intermediators interacting with the targets of at least 10 RIDs. The enriched biological functions of drug targets and intermediators were analyzed by the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/).12,13 Gene ontology terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms were used as biological functions in the current study. Terms with a Benjamini−Hochberg (BH) corrected p value < 0.05 were deemed significant. A heat map with −log(BH p value) as a color index was generated for the top 10 most significantly over-represented terms (by p value) of the sets of drug targets and intermediators to visualize the different contents of enriched biological functions across different gene sets. To identify the intermediators that are significantly associated with the RIDs, we generated 1000 random networks, with each network consisting of the same number of drugs as the RIDs randomly selected from the 6825 drugs available in DrugBank. Z-test was used to evaluate the significance of each intermediator’s drug-degree (the number of interacting drugs) in the network of RIDs against those from the 1000 random networks. Significant intermediators, defined by a false discovery rate (FDR) less than 0.05, which is the expected proportion of false discoveries,14 and the number of interacting drugs ≥5 were deemed as specifically associated with the on-target pharmacological network of RIDs and were subject to further functional-enrichment analysis as described above. Transcriptomics-Based Toxicological Network Analysis. Transcriptomic profiles of the RIDs were obtained from the Connectivity Map from the Broad Institute of MIT and Harvard (http://www.broadinstitute.org/cmap/).6 The Connectivity Map contains a collection of genome-wide transcriptional expression data from cultured human cancer cells treated with bioactive small molecules profiled using Affymetrix Human Genome U133A 2.0 arrays. Seventy five of the curated 211 RIDs were found to be included in the Connectivity Map build 02. To create a representative gene signature for each drug, the 250 most upregulated genes, 250 most downregulated genes, gene-expression profiles from different cell types, dosages, and treatment duration were combined using Kru−Bor merging method.15 A simple frequency-based approach was used to identify the most frequently perturbed genes by the RIDs. The genes in the signatures of 75 drugs were ranked by frequency; the top 1% most frequently perturbed genes were further examined. Drugs with CPT2 as one of the signature genes were subject to a pathwayenrichment analysis to identify the commonly perturbed pathways by these drugs. Integration of the Extended Pharmacological and Transcriptional Networks. Two whole-genome gene-expression profiles of RIDs were collected and compared against the pharmacological networks of RIDs. The first microarray data set included the 191 significantly modulated genes in the muscles of horses suffering from recurrent exertion rhabdomyolysis.16 These modulated genes included human, mouse, and horse genes, which were all mapped to human genes. The second profiling data included 105 genes reported to be differentially expressed in human muscle tissues by statin treatment.4 The modulated genes by disease and drug treatment were compared with the components in the extended pharmacological network of RIDs.
Figure 2. Scheme of biochemical flux leading to ATP production through oxidation of fatty acids. Oxidation of fatty acids takes place in mitochondria with their entry through CPT2 and non-CPT2 routes. Their common intermediate product of acetyl coA will then be used to produce ATP via the TCA cycle or via the 3-hydroxyl-3methylbutanyl-CoA (HMG-CoA) route. For the use case of stastins, detailed branches leading to ATP production from HMS-CoA are provided in this scheme. The enzymes are highlighted in orange, ATP produced, in red, comments of statins, in navy blue, and the rest, in black. Italicized x represented the fraction x. right of a mathematical sign or relationship always represents consumption and production, respectively. The biochemical items following the letter J, therefore, represent substrate or product corresponding to consumption and production, respectively. On the basis of Figure 2, we can express the flow of flux in the following, with x representing the physiological fraction.
JacetylcoA = xJATP(TCA) + (1 − x)Jmevalonate
(1)
Jmevalonate = Jfarnesylpyrophosphate =
∑ (Jgeranylgeranylpp + Jcholesterol + JcoQ10 + Jothers ) (2)
Jgeranylgeranylpp → JRabprenylation → JATP(rab) Jmevalonate → JATP(coQ10) + JATP(rab) CPT2 is upstream of the acetyl coA production pathway and is several reactions away from the step of 3-hydroxyl-3-methylglutanyl-CoA to mevalonic acid by HMG-CoA reductase that statins (HMG-CoA reductase inhibitors) inhibit. Because CPT2 is involved only in 423
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
Article
Figure 3. Base on-target pharmacological network of drug-inducing rhabdomyolysis drugs. A pharmacological network was generated including 167 rhabdomyolysis-inducing drugs with at least one known target, 237 human drug targets, and 35 nonhuman drug targets. The node color represents different type of entities: cyan, drug; red, human drug target; and magenta, nonhuman drug-target. (A) Base pharmacological network. (B) Biggest subnetwork, consisting mostly of psycholeptic or neurological drugs. (C) Subnetwork consisting of antiseizure drugs. (D) Subnetwork consisting of statins. (E) Subnetwork consisting of anticancer agents.
■
RESULTS Construction of On-Target and Extended Pharmacological Networks. We curated and compiled a list of 211 RIDs from multiple resources (Table S1 in the Supporting Information). The currently known drug targets were collected from three databases (DrugBank, TTD, and PharmGKB) for the 211 RIDs, among which 167 drugs had at least one known protein target. The base pharmacological network of RIDs was generated using these 167 drugs and 272 known targets (237 human proteins and 35 nonhuman proteins) (Figure 3). To identify those proteins that are not on-targets of RIDs but are potentially critical in the molecular network of RIDs, we extended this base network by adding intermediators. The criterion for an intermediator to be considered is that it has either genetic or protein−protein interactions with any of the drug targets in the network, where interaction is derived from the BioGRID database. Five extended networks were generated using from 1 to 10 minimum drug degrees (with interacting drugs in the range of 1−10, respectively), and their statistics are summarized in Table S2 in the Supporting Information. A heat map of statistical significance of over-represented biological functions among the drug targets and intermediators was generated (Figure 4). This heat map illustrates the 10 most significantly over-represented biological functions in each set, where the color gradient corresponds to the level of
oxidation of a portion of medium-chain fatty acids and all long-chain fatty acids, for understanding the additive interaction of CPT2 deficiency and statins, flux analysis could be focused on CPT2 only. If considering the fraction involving CPT2
JacetylcoA(CPT2) = xJATP(TCA)(CPT2) + (1 − x)Jmevalonate(CPT2)
(3)
Jfao(CPT2) = Jmfao(CPT2) + Jlfao(CPT2) = JacetylcoA(CPT2)
(4)
where fao represents fatty acid oxidation, m is medium chain, and is l long chain. Jmevalonate(CPT2) = Jfarnesylpp(CPT2) =
∑ Jgeranylgeranylpp(CPT2) + Jcholesterol(CPT2) + JcoQ10(CPT2) + Jothers(CPT2)
(5)
The fraction not involving CPT2 would be as follows. Jmevalonate(non ‐ CPT2) = Jfarnesylpp(non ‐ CPT2) =
∑ Jgeranylgeranylpp(non ‐ CPT2) + Jcholesterol(non ‐ CPT2) + JcoQ10(non ‐ CPT2) + Jothers(non ‐ CPT2)
(6) 424
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
Article
Figure 4. Over-represented biological function. A heat map was generated using the top 10 most over-represented biological functions in each gene set (the human target set and four sets of intermediators, with the number of interacting drugs ranging from 2 to 5). Benjamini−Hochberg (BH)corrected p values of the top 10 most significant functional terms are represented in a heat map with −log(BH-corrected p value) as color index and number. A black cell indicates that the corresponding biological function in the gene set was not statistically significant.
significance, defined as −log(BH p value). The lower the BH p value is, the darker the red color is, suggesting a higher level of significance. The complete details of the functional enrichment results of the drug targets and Min05 networks are in Tables S3 and S4 (Supporting Information). Because the numbers of significantly over-represented functional terms (BH p value
Drug-Induced Rhabdomyolysis: From Systems Pharmacology Analysis to Biochemical Flux Junguk Hur,†,‡,⊥ Zhichao Liu,§,⊥ Weida Tong,§ Reijo Laaksonen,∥ and Jane P. F. Bai*,‡ †
Department of Neurology, University of Michigan, Ann Arbor, Michigan 48109, United States Office of Clinical Pharmacology, Center for Drug Evaluation and Research, U.S. Food and Drug Administration, Silver Spring, Maryland 20993, United States § Divisions of Bioinformatics and Biostatistics, National Center for Toxicological Research, U.S. Food and Drug Administration, Jefferson, Arkansas 72079, United States ∥ The School of Medicine, University of Tampere, Tampere 33520, Finland ‡
S Supporting Information *
ABSTRACT: The goal of this study was to integrate systems pharmacology and biochemical flux to delineate drug-induced rhabdomyolysis by leveraging prior knowledge and publicly accessible data. A list of 211 rhabdomyolysis-inducing drugs (RIDs) was compiled and curated from multiple sources. Extended pharmacological network analysis revealed that the intermediators directly interacting with the pharmacological targets of RIDs were significantly enriched with functions such as regulation of cell cycle, apoptosis, and ubiquitin-mediated proteolysis. A total of 78 intermediators were shown to be significantly connected to at least five RIDs, including estrogen receptor 1 (ESR1), synuclein gamma (SNCG), and janus kinase 2 (JAK2). Transcriptomic analysis of RIDs profiled in Connectivity Map on the global scale revealed that multiple pathways are perturbed by RIDs, including ErbB signaling and lipid metabolism pathways, and that carnitine palmitoyl transferase 2 (CPT2) was in the top 1 percent of the most differentially perturbed genes. CPT2 was downregulated by nine drugs that perturbed the genes significantly enriched in oxidative phosphorylation and energy-metabolism pathways. With statins as the use case, biochemical pathway analysis on the local scale implicated a role for CPT2 in statin-induced perturbation of energy homeostasis, which is in agreement with reports of statin−CPT2 interaction. Considering the complexity of human biology, an integrative multiple-approach analysis composed of a biochemical flux network, pharmacological on- and off-target networks, and transcriptomic signature is important for understanding drug safety and for providing insight into clinical gene−drug interactions.
■
INTRODUCTION
cascade effect of disrupting intracellular physiological homeostasis and leading to apoptosis or necrosis. Considering the complex cross-talk among intracellular pathways and the accumulated knowledge bases as well as databases of systems biology, the extended pharmacological network of rhabdomyolysis-induced drugs (RIDs) can provide insight into the pathways related to drug-induced rhabdomyolysis. The extended pharmacological network is defined as consisting of drugs, their targets, and proteins or genes that interact, under a prespecified statistical criterion, with the drugs via their targets. To understand serious ADR retrospectively, individual laboratories in the scientific community conducted studies and published results, typically with a focus on a specific component of the systems of human biology and physiology.
Drug-induced rhabdomyolysis, although rare, is a serious adverse reaction and can be fatal. Cerivastatin was withdrawn from the market because of treatment-related rhabdomyolysis associated with its use.1 Myopathy and myalgia induced by drugs may be related to drug-induced rhabdomyolysis, although they are not as serious as rhabdomyolysis and may result in a dose reduction or switching of medications. Prospective prediction of adverse drug reactions (ADRs) at the patient level has been very challenging. A recent study by Vilar et al.2 compiled 196 drugs across a wide range of pharmacologic classes. In addition to the structural similarities investigated by Vilar et al., the biological pathways perturbed by these 196 drugs and their pharmacologic networks could overlap to varying degrees with parts of the perturbed biological networks; however, they did so completely independent of one another. In addition to engaging specific pharmacological target(s), each of these drugs can interact with off-targets, resulting in a © 2014 American Chemical Society
Special Issue: Systems Toxicology Received: November 5, 2013 Published: January 14, 2014 421
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
■
For example, there are reports that statins reduce in vitro ATP production3 and perturb several functional pathways in human muscles and change the lipid profiles in circulation;4 there are also case reports of carnitine palmitoyl transferase 2 (CPT2) deficiency carriers developing treatment-related rhabdomyolysis while receiving statins.5 Clearly, there is a need to change the current paradigm of scientific investigation of ADRs to achieve the goal of predicting drug safety. There is also a critical need to leverage existing knowledge and curated data for an integrative analysis for the purpose of developing new drugs with a better effectiveness/risk index. It is our goal to develop predictive models for drug-induced adverse reactions, including drug-induced rhabdomyolysis. We hypothesize that some gene−drug interactions that underpin specific ADRs may not be detectable by genome-wide association studies but could be understood and predicted using systems pharmacology analysis. To establish a predictive framework and to lay the groundwork for achieving future success in developing predictive models for ADRs, we propose that both the extended pharmacological network of RIDs, including their on- and off-target effects, and the genes perturbed by RIDs, as shown by transcriptomic analysis, could be linked to the biochemical pathways that are essential for intracellular energy homeostasis. We first constructed the extended pharmacological base network by utilizing the existing knowledge of drug targets (such as http://www.drugbank.ca/) and then expanded it by leveraging our knowledge of protein− protein and genetic interactions (http://thebiogrid.org/). We propose to use this computational exercise to define, as broadly as possible, the scope of the biological plausibility for druginduced rhabdomyolysis. Then we utilized the published transcripomic database of Connectivity Map6 to delineate the functional and biochemical pathways that are perturbed by RIDs. We further compared these results from cancer cell lines and extended pharmacological network analysis with the transcriptomic profiles of human muscles that were obtained from subjects receiving statins4 to determine and understand how these different approaches interact among one another in the context of understanding the toxicity of drugs at the molecular level. On the basis of our preliminary transcriptomic analysis of RIDs, oxidative phosphorylation was identified as a key pathway perturbed. Because statins are widely used for preventive care of cardiovascular diseases, their interactions with CPT2 via biochemical pathways were used as a use case for gaining insight into drug-induced rhabdomyolysis. The SLCO1B1 gene, encoding the organic anion-transporting polypeptide (OATP1B1) that is involved in hepatic uptake of stains, is associated with statin-induced myopathy.7 In addition to the genes involved in absorption, distribution, metabolism, and elimination (ADME), there are non-ADME genes reportedly associated with drug-induced rhabdomyolysis8 and statin myopathy.5 For example, CPT2 deficiency resulting from mutations in the CPT2 gene9 and gene−drug interactions between CPT2 mutation and drug action8 cause genetic rhabdomyolysis and drug-induced rhabdomyolysis, respectively. To demonstrate further the impact of a specific perturbed pathway on cellular physiological homeostasis, we adopted a qualitative balance analysis of ATP flux that is produced through the biochemical pathway/network involving CPT2 and the target of statins. This study demonstrates the advantage of utilizing three different approaches for an integrative systemspharmacology understanding of drug-induced rhabdomyolysis.
Article
EXPERIMENTAL PROCEDURES
Compilation and Curation of the Drug List. The list of RIDs was compiled and curated in multiple steps, as illustrated in Figure 1.
Figure 1. Overall workflow of our study, consisting of construction of pharmacological networks, transcriptomic analysis, and biochemical flux balance analysis. To be consistent across different sets, we limited our search term specifically to rhabdomyolysis. A text-mining approach was developed to screen automatically the large volume of drug labels and packaging. Briefly, the most recent drug labels (in PDF file format) of each approved drug were retrieved from the Drugs@FDA database (http:// www.accessdata.fda.gov/scripts/cder/drugsatfda/). The PDF files of drug labels were converted to text files using Xpdf (http://www. foolabs.com/xpdf/). Text-converted drug labels from DrugBank and XML-formatted drug information from DailyMed (http://dailymed. nlm.nih.gov/) were analyzed to identify any drug names and active ingredient information for the list of drugs for which the term rhabdomyolysis was included in the adverse event section. The flanking text surrounding the rhabdomyolysis term was manually checked to make sure that rhabdomyolysis was indeed a treatmentrelated adverse event. Another list was compiled from the SIDER2 database (http://sideeffects.embl.de/), a computer-readable drug and human side-effects relationship data source. The SIDER2 database consists of information from different version of drug labels and health-related public documents. We also collected RIDs from a published paper by Vilar et al.2 and a curated list from the FDA (http://www.fda.gov/). To reduce redundancy and to ensure consistency, we normalized drug names. The official drug names used in DrugBank were adopted. Any drug with no corresponding record in the DrugBank database was processed by a Perl script to identify automatically the active pharmaceutical ingredient rather than generic drug names. For example, Olmesartan Medoxomil was normalized to Olmesartan. The list of drugs obtained from aforementioned resources were normalized, combined, and carefully reviewed. Of note, when compiling the list of RIDs, we did not differentiate the drugs with incidence rates reported during phase 3 clinical trials from those with postmarketing reports only. Only the drug products of a single active pharmaceutical ingredient were compiled. To curate the final list, we further referred to the clinical criteria of rhabdomyolysis, which is defined as muscle symptoms plus significant creatinine kinase elevation (usually >40 times the upper limit of the normal range (ULN)) and creatinine elevation (with brown urine/urine myoglobin).10,11 Extended Pharmacological Network (On Targets Plus Their Interacting Proteins) Analysis. A base pharmacological network, consisting of RIDs and their known drug targets, was constructed by referencing the drug databases, including DrugBank, Therapeutic Target Database (http://bidd.nus.edu.sg/group/cjttd/) (TTD), and PharmGKB (http://www.pharmgkb.org/). Drugs with any known targets were included in the base network. To identify common but 422
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
Article
Biochemical Flux Analysis. Qualitative biochemical flux balance was analyzed to complement the pharmacological targets- and transcriptomics-based network analyses for the use case of statininduced rhabdomyolysis. CPT2 deficiency or mutation has been reported in statin-induced myopathy or rhabdomyolysis.5 Statins were used to illustrate gene−drug interaction, with the focus on the CPT2 gene encoding carnitine palmitoyl transferase 2 (CPT2). In light of the biological complexity, we focused on the linear paths from fatty acid oxidation to ATP production involving CPT2 by ignoring any feedback loops as well as on the biosynthesis of prenylated Rab and coenzyme Q10 (CoQ10) (also known as ubiquinone),17 which are reportedly related to statin-induced myopathy or rhabdomyolysis.18 We did this by focusing on ATP production from fatty acid oxidation (ATP is produced from the tricarboxylic acid cycle) and from the biological processes involving CoQ1017 as well as prenylated Rab1.18 For the purpose of flux analysis, the immediate paths forward from fatty acid oxidation are designated as CPT2 and non-CPT2 involved (Figure 2). Throughout the text, the flux denoted as J on the left and
potentially crucial regulator or effector in the pathogenesis of druginduced rhabdomyolysis, we extended the base network by adding genes that directly interact with the drug targets via either protein− protein interactions (PPIs) or genetic interactions (GIs). These additional genes, referred to as intermediator in the current study, are not a direct target of these drugs; however, they may play important roles in the pathogenesis of drug-induced rhabdomyolysis by interacting with the targets of RIDs. The PPIs and GIs were derived from the Biological General Repository for Interaction Data sets (BioGRID; release 3.2.97). Ten extended networks, denoted as Min01−Min10, were generated using the number of interacting drugs from 1−10 as the criterion for intermediators to be included in the extended networks. The Min01 network included any intermediators with known interactions with the targets of any RIDs, whereas Min10 included only the intermediators interacting with the targets of at least 10 RIDs. The enriched biological functions of drug targets and intermediators were analyzed by the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/).12,13 Gene ontology terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms were used as biological functions in the current study. Terms with a Benjamini−Hochberg (BH) corrected p value < 0.05 were deemed significant. A heat map with −log(BH p value) as a color index was generated for the top 10 most significantly over-represented terms (by p value) of the sets of drug targets and intermediators to visualize the different contents of enriched biological functions across different gene sets. To identify the intermediators that are significantly associated with the RIDs, we generated 1000 random networks, with each network consisting of the same number of drugs as the RIDs randomly selected from the 6825 drugs available in DrugBank. Z-test was used to evaluate the significance of each intermediator’s drug-degree (the number of interacting drugs) in the network of RIDs against those from the 1000 random networks. Significant intermediators, defined by a false discovery rate (FDR) less than 0.05, which is the expected proportion of false discoveries,14 and the number of interacting drugs ≥5 were deemed as specifically associated with the on-target pharmacological network of RIDs and were subject to further functional-enrichment analysis as described above. Transcriptomics-Based Toxicological Network Analysis. Transcriptomic profiles of the RIDs were obtained from the Connectivity Map from the Broad Institute of MIT and Harvard (http://www.broadinstitute.org/cmap/).6 The Connectivity Map contains a collection of genome-wide transcriptional expression data from cultured human cancer cells treated with bioactive small molecules profiled using Affymetrix Human Genome U133A 2.0 arrays. Seventy five of the curated 211 RIDs were found to be included in the Connectivity Map build 02. To create a representative gene signature for each drug, the 250 most upregulated genes, 250 most downregulated genes, gene-expression profiles from different cell types, dosages, and treatment duration were combined using Kru−Bor merging method.15 A simple frequency-based approach was used to identify the most frequently perturbed genes by the RIDs. The genes in the signatures of 75 drugs were ranked by frequency; the top 1% most frequently perturbed genes were further examined. Drugs with CPT2 as one of the signature genes were subject to a pathwayenrichment analysis to identify the commonly perturbed pathways by these drugs. Integration of the Extended Pharmacological and Transcriptional Networks. Two whole-genome gene-expression profiles of RIDs were collected and compared against the pharmacological networks of RIDs. The first microarray data set included the 191 significantly modulated genes in the muscles of horses suffering from recurrent exertion rhabdomyolysis.16 These modulated genes included human, mouse, and horse genes, which were all mapped to human genes. The second profiling data included 105 genes reported to be differentially expressed in human muscle tissues by statin treatment.4 The modulated genes by disease and drug treatment were compared with the components in the extended pharmacological network of RIDs.
Figure 2. Scheme of biochemical flux leading to ATP production through oxidation of fatty acids. Oxidation of fatty acids takes place in mitochondria with their entry through CPT2 and non-CPT2 routes. Their common intermediate product of acetyl coA will then be used to produce ATP via the TCA cycle or via the 3-hydroxyl-3methylbutanyl-CoA (HMG-CoA) route. For the use case of stastins, detailed branches leading to ATP production from HMS-CoA are provided in this scheme. The enzymes are highlighted in orange, ATP produced, in red, comments of statins, in navy blue, and the rest, in black. Italicized x represented the fraction x. right of a mathematical sign or relationship always represents consumption and production, respectively. The biochemical items following the letter J, therefore, represent substrate or product corresponding to consumption and production, respectively. On the basis of Figure 2, we can express the flow of flux in the following, with x representing the physiological fraction.
JacetylcoA = xJATP(TCA) + (1 − x)Jmevalonate
(1)
Jmevalonate = Jfarnesylpyrophosphate =
∑ (Jgeranylgeranylpp + Jcholesterol + JcoQ10 + Jothers ) (2)
Jgeranylgeranylpp → JRabprenylation → JATP(rab) Jmevalonate → JATP(coQ10) + JATP(rab) CPT2 is upstream of the acetyl coA production pathway and is several reactions away from the step of 3-hydroxyl-3-methylglutanyl-CoA to mevalonic acid by HMG-CoA reductase that statins (HMG-CoA reductase inhibitors) inhibit. Because CPT2 is involved only in 423
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
Article
Figure 3. Base on-target pharmacological network of drug-inducing rhabdomyolysis drugs. A pharmacological network was generated including 167 rhabdomyolysis-inducing drugs with at least one known target, 237 human drug targets, and 35 nonhuman drug targets. The node color represents different type of entities: cyan, drug; red, human drug target; and magenta, nonhuman drug-target. (A) Base pharmacological network. (B) Biggest subnetwork, consisting mostly of psycholeptic or neurological drugs. (C) Subnetwork consisting of antiseizure drugs. (D) Subnetwork consisting of statins. (E) Subnetwork consisting of anticancer agents.
■
RESULTS Construction of On-Target and Extended Pharmacological Networks. We curated and compiled a list of 211 RIDs from multiple resources (Table S1 in the Supporting Information). The currently known drug targets were collected from three databases (DrugBank, TTD, and PharmGKB) for the 211 RIDs, among which 167 drugs had at least one known protein target. The base pharmacological network of RIDs was generated using these 167 drugs and 272 known targets (237 human proteins and 35 nonhuman proteins) (Figure 3). To identify those proteins that are not on-targets of RIDs but are potentially critical in the molecular network of RIDs, we extended this base network by adding intermediators. The criterion for an intermediator to be considered is that it has either genetic or protein−protein interactions with any of the drug targets in the network, where interaction is derived from the BioGRID database. Five extended networks were generated using from 1 to 10 minimum drug degrees (with interacting drugs in the range of 1−10, respectively), and their statistics are summarized in Table S2 in the Supporting Information. A heat map of statistical significance of over-represented biological functions among the drug targets and intermediators was generated (Figure 4). This heat map illustrates the 10 most significantly over-represented biological functions in each set, where the color gradient corresponds to the level of
oxidation of a portion of medium-chain fatty acids and all long-chain fatty acids, for understanding the additive interaction of CPT2 deficiency and statins, flux analysis could be focused on CPT2 only. If considering the fraction involving CPT2
JacetylcoA(CPT2) = xJATP(TCA)(CPT2) + (1 − x)Jmevalonate(CPT2)
(3)
Jfao(CPT2) = Jmfao(CPT2) + Jlfao(CPT2) = JacetylcoA(CPT2)
(4)
where fao represents fatty acid oxidation, m is medium chain, and is l long chain. Jmevalonate(CPT2) = Jfarnesylpp(CPT2) =
∑ Jgeranylgeranylpp(CPT2) + Jcholesterol(CPT2) + JcoQ10(CPT2) + Jothers(CPT2)
(5)
The fraction not involving CPT2 would be as follows. Jmevalonate(non ‐ CPT2) = Jfarnesylpp(non ‐ CPT2) =
∑ Jgeranylgeranylpp(non ‐ CPT2) + Jcholesterol(non ‐ CPT2) + JcoQ10(non ‐ CPT2) + Jothers(non ‐ CPT2)
(6) 424
dx.doi.org/10.1021/tx400409c | Chem. Res. Toxicol. 2014, 27, 421−432
Chemical Research in Toxicology
Article
Figure 4. Over-represented biological function. A heat map was generated using the top 10 most over-represented biological functions in each gene set (the human target set and four sets of intermediators, with the number of interacting drugs ranging from 2 to 5). Benjamini−Hochberg (BH)corrected p values of the top 10 most significant functional terms are represented in a heat map with −log(BH-corrected p value) as color index and number. A black cell indicates that the corresponding biological function in the gene set was not statistically significant.
significance, defined as −log(BH p value). The lower the BH p value is, the darker the red color is, suggesting a higher level of significance. The complete details of the functional enrichment results of the drug targets and Min05 networks are in Tables S3 and S4 (Supporting Information). Because the numbers of significantly over-represented functional terms (BH p value
