Visualization of drug target interactions in the contexts of pathways and networks with ReactomeFIViz

Visualization of evidence of drug-target interactions

ReactomeFIViz provides features for visualization of drug-target interaction information for 171 FDA-approved cancer drugs from the Cancer Targetome²⁴ and for 2,102 world-wide approved drugs from DrugCentral^25,26, facilitating researchers in investigating supporting evidence for interactions between a drug and all its targets. These features address a major bottleneck in precision therapy research - the lack of rigorous and transparent supporting information on drugs and their targets. While it is well-established that many approved drugs are promiscuous and interact with many targets, it is more difficult to obtain and visualize reliable information about such additional targets^27,28. Resources such as the Cancer Targetome and DrugCentral provide thorough drug-target information, but do not provide visualizations for these drug-target interactions. By adding functionality for researchers to visualize drug-target interactions, ReactomeFIViz empowers investigation of drug-targets in a supporting evidence-based manner.

All drugs collected in the Cancer Targetome and DrugCentral can be accessed within ReactomeFIViz. For each drug, target interaction evidence can be filtered according to strength of the supporting assay values and displayed as either a table or as a histogram (Figure 2). This allows the user to assess drug-target relationships on the basis of supporting evidence. For example, the histogram in Figure 2 for target assay values for the drug sorafenib shows that there are many potential targets with assay values under 100 nM. Targets that have multiple assay values under 100 nM include FLT3, RET, KIT, RAF1, and BRAF. Additional targets with very low assay values include DDR1, DDR2, FLT4, and KDR, among others. Overall there is evidence across multiple assay types to support sorafenib interactions with a variety of targets. Sorafenib is often referred to as a “multi-kinase” inhibitor, and promiscuity (interactions with many targets) has been well noted in the literature^29,30. The ability to display and threshold the full range of supporting evidence is a key requirement for researchers engaged in drug discovery efforts. Some applications, such as the nomination of compounds for drug screens, will require very strong binding assay evidence and would benefit from being able to compare results across multiple assay types. Other applications, especially those that perform in silico modeling, may have less stringent requirements for binding assay values. For instance, researchers interested in exploring compounds for drug repurposing may want to look at more weakly binding compounds as a starting place before further optimization.

Figure 2. Visualizing drug-target interaction evidence for FDA-approved drug sorafenib via a histogram of drug-target assay values categorized by assay types (KD, EC50, IC50, and Ki).

Sorafenib interacts with many targets, even when restricting to target interactions supported by binding assay evidence < 100 nM.

Drug-centric perspective on drug-targeted pathways

Depending on the application, different perspectives on drug-target-pathway interaction data are necessary. If we are focused on investigating a particular drug or a small number of related drugs, we would want to explore targets and pathways. For instance, if we want to investigate off-target or toxic effects of a certain drug, we may want to consider all possible targets and pathways with which the drug interacts. In such a scenario, we can look up all the target interactions for a particular drug and map them to pathways. Furthermore, performing enrichment analysis identifies pathways with a significant number of targeted entities, suggesting pathways most perturbed by the drug. The top enriched pathways for sorafenib targets with supporting assay values <= 100 nM are shown in Table 1. Sorafenib is a receptor tyrosine kinase inhibitor, which is known experimentally to interact with a variety of targets. These targets are involved in several signaling pathways. For instance, we can see (Table 1) that many of sorafenib’s targets are involved in the RAF/MAP kinase cascade as well as other pathways involving VEGF signaling and PIP3/AKT signaling. Examining the full range of pathways targeted by a drug allows us to better understand the drug’s mechanism of action for both efficacy and side effects.

Pathway-centric perspective on drug-target interactions

Alternatively, we may start with a pathway of interest and investigate all drugs that target components of the pathway within a selected affinity range. This perspective allows a broad view for assessing how confident we may be able to find a drug targeted to a particular pathway, as many pathways have only weakly or no interacting drugs²⁴. Figure 3 highlights drugs targeting KIT in the pathway “SCF-KIT signaling” (https://reactome.org/content/detail/R-HSA-1433557), which are retrieved by using the “Fetch Cancer Drugs” feature in ReactomeFIViz to visualize drug-target interactions collected in the Cancer Targetome. As displayed in Figure 3, KIT is targeted by ten different drugs, which are all kinase inhibitors: erlotinib, imatinib, bosutinib, sunitinib, pazopanib, dasatinib, sorafenib, axitinib, vandetanib, and nilotinib^24,29, and supported by assay values less than 100 nM. This view suggests a list of potential drug candidates for researchers to perturb the activity of the SCF-KIT signaling pathway.

Figure 3. Drugs targeting KIT in the “SCF-KIT Signaling” pathway.

Drug-target interactions were fetched from the Cancer Targetome and supported with multiple assay types having values <= 100 nM.

Constrained fuzzy logic modeling of drug action on pathways

Overlaying drug-target interactions onto a pathway reveals potential drug action targets in the pathway. However, the topological structure of a pathway may be complicated and contain one or more feedback loops. This makes it difficult to predict the dynamic effect of a drug on pathway activity by inspection only. To assist researchers in performing computational modeling of drug-induced pathway perturbation, we developed an automated approach to convert Reactome pathways into Boolean networks (BNs). These BNs were then subject to constrained fuzzy logic simulation^31,32 to model a drug’s perturbation on the activity of a particular pathway. We conducted two constrained fuzzy logic simulations: a reference simulation without considering drug action, and a perturbation simulation considering drug action. To model drug action within the simulation, we mapped drug-target interaction affinity values as “activation” or “inhibition” strength according to the drug’s action mechanism and its assay values collected in the drug data sources. We then used this information to modify the transfer functions related to the drug's target(s) in the fuzzy logic model. Finally, we calculated the relative impact scores for individual entities in the pathway by comparing the model that includes drug effects to the one that does not. For details, see Methods.

To illustrate this capability, we take the example of modeling the perturbation induced by the kinase inhibitor sorafenib on the “Signaling by SCF-KIT” pathway. Stem cell factor (SCF), together with its receptor, c-KIT, a tyrosine kinase receptor, regulates pathways related to proliferation, migration, survival, and differentiation of multiple cell types³³. These regulated pathways include RAF/MAP kinase signaling, PIP3-activated AKT signaling, and the JAK/STAT signaling pathway. The Reactome-annotated pathway “Signaling by SCF-KIT” includes a set of linked reactions that produce entities activating these regulated pathways.

Drugs targeting proteins involved in SCF-KIT signaling have been under active study for many years³⁴. Evidence shows that sorafenib binds strongly to c-KIT (Figure 2 and Figure 3) with a minimum binding constant of 16 nM (KD) according to the Cancer Targetome. We performed four constrained fuzzy logic simulations to study the perturbation of sorafenib on the activity of a complex called “p-STAT dimers”, which is one of the major outputs of the SCF-KIT pathway. “p-STAT dimers” activates expression of genes regulated by STAT proteins and is produced by the reaction “Disassociation and translocation of STATs to the nucleus”. This reaction, together with three other reactions, “Recruitment of STATs”, “Phosphorylation of STATs”, and “Dimerization of STATs”, forms a loop, providing a way to produce activated STAT dimer by recycling the “p-JAK2:SFKs:p-KIT" complex (Figure 4). The four simulations were:

Figure 4. Partial diagram of the “Signaling by SCF-KIT” pathway showing the feedback loop that generates the “p-STAT dimers” complex annotated in Reactome.

The entry point of the loop and “p-STAT dimers”, together with four reactions forming the loop, are highlighted in light blue. The reactions are labeled with names in the light-grey boxes. The target of drug sorafenib, KIT protein, is in the upstream of this loop and not shown here.

In Figure 5, we show the activity (labeled as Logic Fuzzy Value in the y-axis) of “p-STAT dimers” over time steps in the simulation. It is well known that pathways having feedback loops usually show periodic cycle patterns or attractors for activities of entities involved in the loops³⁵. As expected, we saw these attractors in all four simulations. In the first drug perturbation modeling (reference_1 and sorafenib_1, Figure 5A) with the initial value of PRKCA equal to 1.0 as fully activated, the application of sorafenib did not notably perturb the activity of “p-STAT dimers” (relative impact score = 8.63E-4), most likely because the loop used to produce “p-STAT dimers” was not affected by sorafenib under this simulation’s condition. However, by reducing the initial value of PRKCA from 1.0 to 0.5, the cycle attractor of “p-STAT dimers” changed from (0, 1, 0, 0, 1, 1) to (0.5, 1, 0.5, 0.5, 1, 1) in the absence of sorafenib (simulation reference_2, Figure 5B), increasing the lower activity in the attractor. Reduction of the initial value of PRKCA brought down the activity of “p-KIT (S741,746):PKC alpha”, the inhibitor of the reaction, “Autophosphorylation of KIT”, resulting in an increase of the activity of “p-JAK2:SFKs:p-KIT complex” via several intermediate reactions. “p-JAK2:SFKs:p-KIT complex” is the entry point of the loop that produces “p-STAT dimers” (Figure 4). Increasing the activity level of this complex levels up the p-STAT dimers, thereby increasing the lower value of its cycle attractor. Intriguingly, the application of sorafenib abolished the effect of the reduction of PRKCA’s initial value, bringing down the lower value of the cycle attractor from 0.5 to 0.0024, the same value when PRKCA’s initial value was set at 1.0. Hence, sorafenib had a much stronger perturbation when PRKCA’s initial activity was lower (relative impact score = -0.17). For the detailed pathway annotation, see Figure 6 and also refer to the Reactome pathway diagram: https://reactome.org/PathwayBrowser/#/R-HSA-1433557.

Figure 5. Constrained fuzzy logic simulation results for complex “pSTAT-dimers” in pathway “Signaling by SCF-KIT” from 4 simulations.

Panel A: Simulations conducted with the default setup provided by ReactomeFIVIz. Logic fuzzy values for two simulations, reference_1 and sorafenib_1, are almost the same except bottom values starting from step 18: 0.0 for reference_1 and 0.0024 for sorafenib_1 (see the insert for an example). Panel B: Same as in Panel A except the initial value of PRKCA was reduced from 1.0 to 0.5. Logic fuzzy values prior to time step 18 are the same for both reference and sorafenib simulations. The results were exported from ReactomeFIViz and plotted in Microsoft Excel.

The above comparison analysis of the simulation results demonstrates the importance of mathematical modeling, even without experimentally measured parameters, by predicting that alterations in the PRKCA gene may impact the effects of the drug sorafenib on SCF-KIT signaling. For instance, a loss-of-function mutation in the PRKCA gene could influence how sorafenib affects the SCF-KIT pathway. Alternatively, another drug that inhibits PRKCA could be used in combination with sorafenib to potentially impact the pathway in a different manner than with sorafenib alone. This level of mechanistic and dynamic investigation is not possible by merely overlaying drug-target interactions onto the pathway.

We also performed analysis for another output of this pathway, “PI(3,4,5)P3” (https://reactome.org/PathwayBrowser/#/R-HSA-1433557&SEL=R-ALL-179838) (Extended data: Figure S1), which was annotated as an activator for pathway “PIP3 activates AKT signaling”. We observed similar behavior as in the case of “p-STAT dimers”: reducing the initial value of PRKCA increased the activity of “PI(3,4,5)P3”. The application of sorafenib modestly increased PI(3,4,5)P3's activity using default PKCA activity level (1.0), but reduced its activity when PRKCA’s activity was reduced. However, the relative impact scores for PI(3,4,5)P3 are much larger compared to the case of the p-STAT dimers described above, most likely because this entity has a converged single stable activity, in contrast to p-STAT dimer’s cycle attractor.

We’d like to emphasize that the above simulations were performed using a set of initial values and transfer function parameters that we developed to ensure simulations with converged solutions. In reality, the actual tumor cell conditions will likely differ from the parameters we used. ReactomeFIViz provides a set of intuitive user interfaces for users to try different parameters in simulation. Future work for this modeling framework will include approaches for estimating and learning these parameters based on large scale omics datasets. We hope to address this daunting issue soon.

To assist researchers in investigating drug impact on pathway activities in a systematic way, ReactomeFIViz provides a feature to highlight entities based on relative impact scores. Figure 6 shows SCF-KIT Signaling pathway highlighted according to the relative impact scores with PRKCA’s initial value equal to 1.0 (Figure 6A) and 0.5 (Figure 6B) for sorafenib.

Figure 6. The “Signaling by SCF-KIT” pathway with entities highlighted according to relative impact scores calculated from constrained fuzzy logic simulations.

A: Relative impact scores for sorafenib with PRKCA’s initial value equal to 1.0 as fully activated. B: Same as A except that the initial value of PRKCA was reduced from 1.0 to 0.5.

Network perspective on drug-target interactions

A network perspective unites multiple pathways together, allowing a bird’s eye view of drug-target relationships. This perspective is especially important when considering drugs with targets that span multiple pathways or considering potential crosstalk between pathways (e.g. Table 1). Projecting drug-target interactions onto a network helps prioritization of targets, targeted pathways, and targeted therapies. Here, we use the TCGA ovarian cancer data³⁶ as an example to show the utility of this network perspective for prioritizing targeted therapies for ovarian cancer.

The Reactome functional interaction (FI) network was constructed by extracting functional relationships from manually curated pathways in Reactome and several other pathway databases, and then combining them with predicted interactions based on a machine learning approach, Naïve Bayes Classifier²³. Using genes mutated in TCGA ovarian cancer samples, we constructed a FI network, where each node represents a gene mutated in at least 5 patients in the TCGA cohort, and each edge represents a functional interaction between two genes.

Combining the TCGA ovarian cancer FI network with the Cancer Targetome reveals that there are 18 different cancer drugs in the Cancer Targetome that potentially interact with target proteins in the ovarian cancer set, at affinity values <= 100 nM. One strategy for prioritizing drugs interacting with members of this network would be to prioritize drugs targeting proteins that are also members of multiple significantly enriched pathways for the ovarian cancer cohort. For instance, epidermal growth factor receptor (EGFR) may be prioritized as a target for ovarian cancer. It has roles in multiple pathways that are significantly enriched for mutated genes, including “L1CAM interactions”, “Focal adhesions”, “Calcium signaling pathway”, “PI3K-AKT signaling pathway”, and “HIF-1 signaling pathway” (Extended data: Table S1). EGFR is targeted by both monoclonal antibody drugs and kinase inhibitors³⁷. 14 of the 18 drugs that interact with targets in the network bind to EGFR (Figure 7). EGFR is a member of the well-studied ERBB protein family and is a target of interest for ovarian cancer³⁷. Many studies have shown encouraging results for the effect of EGFR and other ERBB-family inhibitors on ovarian cancer cell lines³⁸. However, despite promising preclinical evidence, inhibitors targeting the ERBB signaling have up to this point shown little efficacy in patient clinical trials³⁸. Current efforts in this area are focused on the use of ERBB-family inhibitors in combination with cytotoxic or targeted therapies^39,40.

Figure 7. Functional interaction network for frequently mutated genes in TCGA ovarian cancer samples.

Only part of the FI network was shown here for drugs targeting protein products of genes in the network. Network nodes are genes mutated in at least 5 TCGA samples, where size of the node indicates the number of samples with the mutated gene. Black edges are FIs between genes, blue edges drug-target interactions with <= 100 nM supporting assay values, dashed edges predicted FIs, solid edges extracted from annotated FIs, -> activation or catalysis FIs, -| inhibition FI, - FIs extracted from complexes or inputs of reactions. EGFR, TNIK, PAK3 and sunitinib malate are highlighted in yellow.

An alternative strategy for prioritizing drugs interacting with this FI network would be to select drugs interacting with targets in multiple, disjoint pathways to potentially achieve a greater coverage of key genes in the network. For instance, the drug sunitinib malate targets both “TRAF2 And NCK Interacting Kinase” (TNIK) and “P21 (RAC1) Activated Kinase 3” (PAK3) (Figure 7), which are individually involved in different pathways enriched by mutated genes (Extended data: Table S1).

Cancer Targetome

The Cancer Targetome²⁴ aggregates drug-target interactions from four public databases: DrugBank⁴⁶, Therapeutic Targets Database⁴⁷, the IUPHAR/BPS Guide to Pharmacology⁴⁸, and BindingDB⁴⁹. This resource provides interaction information across 171 FDA-approved cancer drugs, 880 protein targets, and over 6800 interactions between these drugs and targets. Each drug interaction is accompanied by supporting evidence across three tiers²⁴: parent database, PubMed IDs for literature references, and experimental assay values. Experimental assay values include IC50, EC50, KD, and Ki values. Each drug-target interaction may be supported by multiple types of assay values. The companion supporting evidence allows for assessment of putative drug-target interactions in an empirical data-driven manner, in particular, for consideration and inclusion of secondary or “off” targets of drugs that may play an important role in a drug’s effects on a particular pathway. For more information, see 24. The version currently supported in ReactomeFIViz was developed in 2016.

DrugCentral

The DrugCentral²⁵ database collects extensive drug information, including but not limited to compound structure, active ingredients, and mechanism of action for approved drugs from the US Food and Drug Administration (FDA), European Medicines Agency (EMA), and the Japan Pharmaceuticals and Medical Devices Agency (PMDA). Target information in DrugCentral is mined from bioactivity resources such as ChEMBL and IUPHAR, and may also be supplemented with manually curated target information from drug approval labels. Overall, DrugCentral (version released in August, 2017) includes information for 2102 compounds, 1550 targets, and 13,407 drug-target interactions. The drug-target interaction file supported by ReactomeFIViz, drug.target.interaction.08292017.tsv, was downloaded from the DrugCentral web site (http://drugcentral.org) in 2017.

TCGA ovarian cancer dataset

For the functional interaction network in Figure 7, we used the Gene/Mutation Set Analysis feature in ReactomeFIViz to create a network with genes mutated in at least 5 samples in the TCGA Ovarian Cancer dataset³⁶. Pathway enrichment p-values in ReactomeFIViz were calculated using the binomial test and were false discovery rate (FDR)-corrected based on the Benjamini-Hochberg approach⁵⁰.

Construction of Boolean networks from Reactome pathways

In Reactome, the biochemical reaction is the basic unit for pathway annotation. Pathways are annotated as a set of biochemical reactions that are linked together. We developed a scheme to automatically convert Reactome pathways into Boolean networks (Figure 8) by adopting our previous approach used to convert Reactome pathways into probabilistic graphical models²⁰ based on the PARADIGM approach¹⁵. A typical reaction in Reactome has one or more inputs, one catalyst, one or more activators, one or more inhibitors, and one or more outputs. A logical “AND” relationship was created among inputs, the catalyst, and activators. Inhibitors were negated and then added into this “AND” relationship. To make the relationship simple, if multiple outputs, activators, or inhibitors are annotated for a reaction, we introduced accessory nodes for them. Therefore, a typical reaction in Figure 8A was converted into a set of Boolean relationships (· for AND, + for OR, ! for NOT):

Figure 8. An automatic scheme to convert Reactome pathways into Boolean networks by handling Reaction, Complex and EntitySet instances.

A: Convert a typical Reactome reaction into a set of Boolean relationships. B: Convert a complex into an AND relationship between complex subunits and the complex. C: Convert an EntitySet into an OR relationship between set members and the set.

Reactome has annotated many complex assembly reactions. However, some complexes are annotated directly for reactions without explicit description of how subunits are assembled together. We expanded a complex like this into an “AND” relationship among its subunits and generated the following Boolean relationship for the complex as illustrated in Figure 8B:

EntitySet is a class in the Reactome data model⁵¹ used to collect a set of entities with the same function in a biochemical reaction and is routinely used to collect protein isoforms. We expanded an EntitySet instance into an “OR” relationship among annotated members:

Thus, by converting biochemical reactions, complexes, and EntitySets into a set of Boolean relationships, we converted a Reactome pathway into a Boolean network.

Constrained fuzzy logic simulation of drug perturbation on pathways

We adopted the constrained fuzzy logic simulation approach^31,32 to simulate the dynamic behavior of Boolean networks that are automatically converted from Reactome pathways by extending two-state Boolean variables to continuous variables with values between 0 and 1, inclusively. This allows our modeling approach to be applied to continuous ‘omics data types and drug-target assay values without the need for discretization. Based on the minimum supporting assay values for drug-target interactions loaded from the drug data sources, we calculated inhibition or activation strength by using the following hard-coded piecewise function:

Function f(input) may be a Hill function^31,32 or an Identity function. For the four simulations conducted with the “Signaling by SCF-KIT” pathway, the Identity function was used.

To measure the perturbation impact caused by drug action, we developed a relative impact score based on the time step plots for individual entities in the pathway. We conducted the simulation twice: one with drugs as the perturbation, and another without drugs as the reference. The relative impact score was then calculated as follows (AUC for area under curve of fuzzy logic variable value vs. time step):

Our test results indicated that the above defined relative impact score, ranging between -1.0 and 1.0, inclusively, converged within a limited number of time steps in most cases.

Software implementation

We adopted the three-tier software architecture: The data tier provides drug-target interactions from the Cancer Targetome and DrugCentral, contents from the Reactome database, and the Reactome FI network; The server tier fetches the content from the data tier and then send it to the front-end app via a RESTful API (application programming interface). The API also provides methods for Boolean network construction and fuzzy logic simulation. The server was implemented upon the Spring framework (https://spring.io) along with the Hibernate OR (object/relational) mapping (http://hibernate.org); the front end provides all user interfaces.

Our implementation of Boolean network in Java referred to the Java implementation of open source CellNOptR⁵², a comprehensive Boolean network modeling toolkit for 'omics data, hosted at https://github.com/saezlab/cytocopter/tree/master. The implementation of constrained fuzzy logic simulation referred to the code (version 1.2) hosted at https://github.com/saezlab/CNORfuzzy by Saez Lab. The simulation ran until either the maximum difference between two consecutive iterations was less than 1.0 × 10^-6, or the total number of iterations reached the larger of 100 or 1.2 × total number of variables in the Boolean network. To calculate the relative impact scores, we expanded the plots of fuzzy variable values vs. time steps based on reached attractors by 20 time steps each time. This expansion stopped when the maximum difference between two consecutive expansions was either less than 0.01, or the expansion reached a total of 1000 time steps. As the default setup for simulation, 1.0 was assigned as the initial values for those entities annotated as inputs of reactions forming loops or inputs that are not annotated as outputs in the pathway.

We used Java 8 (https://java.com) and Eclipse (https://www.eclipse.org) as the integrated development environment (IDE) for programming. ReactomeFIViz is released through the Cytoscape App Store (http://apps.cytoscape.org/apps/reactomefiplugin). The code is open source, hosted at GitHub (https://github.com/reactome-fi/CytoscapePlugIn).

Underlying data

The TCGA ovarian cancer mutation data file: http://cpws.reactome.org/caBigR3WebApp2018/ov.maf.txt.zip

Extended data

Zenodo: Reactome_FIViz_Drug_Visualization_Supp, https://doi.org/10.5281/zenodo.3239441⁵⁴

This project contains the following extended data:

Table S1. Significant pathways for the TCGA ovarian cancer FI network. Pathways enriched by the mutated genes in the TCGA ovarian cancer FI network with FDR <= 0.05 are listed. Highlighted pathways contain EGFR (in yellow except Focal adhesion(K)), PAK3 (in brown), or TNIK (in green), as discussed in the main text. (R) indicates pathways collected from Reactome, (K) from KEGG (1), (N) from NCI Pathway Interaction Database (2), (P) from Panther Pathways.

1. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017 Jan 4;45(D1):D353–61.

2. Schaefer CF, Anthony K, Krupa S, Buchoff J, Day M, Hannay T, et al. PID: the Pathway Interaction Database. Nucleic Acids Res. 2009 Jan;37(Database issue):D674–9.

Figure S1. Constrained fuzzy logic simulation results for entity PI(3,4,5)P3 [plasma membrane] (https://reactome.org/PathwayBrowser/#/R-HSA-1433557&SEL=R-ALL-179838&PATH=R-HSA-162582,R-HSA-9006934) in the “Signaling by SCF-KIT” pathway from 4 simulations. A: Simulations conducted with the default setup provided by ReactomeFIVIz. B: Same as A except the initial value of PRKCA was reduced from 1.0 to 0.5. Logic fuzzy values prior to time step 11 are the same for all four simulations, which is 0.0.

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).