Systems Perturbation Analysis of a Large-Scale Signal Transduction Model Reveals Potentially Influential Candidates for Cancer Therapeutics

Dysregulation in signal transduction pathways can lead to a variety of complex disorders, including cancer. Computational approaches such as network analysis are important tools to understand system dynamics as well as to identify critical components that could be further explored as therapeutic targets. Here, we performed perturbation analysis of a large-scale signal transduction model in extracellular environments that stimulate cell death, growth, motility, and quiescence. Each of the model’s components was perturbed under both loss-of-function and gain-of-function mutations. Using 1,300 simulations under both types of perturbations across various extracellular conditions, we identified the most and least influential components based on the magnitude of their influence on the rest of the system. Based on the premise that the most influential components might serve as better drug targets, we characterized them for biological functions, housekeeping genes, essential genes, and druggable proteins. The most influential components under all environmental conditions were enriched with several biological processes. The inositol pathway was found as most influential under inactivating perturbations, whereas the kinase and small lung cancer pathways were identified as the most influential under activating perturbations. The most influential components were enriched with essential genes and druggable proteins. Moreover, known cancer drug targets were also classified in influential components based on the affected components in the network. Additionally, the systemic perturbation analysis of the model revealed a network motif of most influential components which affect each other. Furthermore, our analysis predicted novel combinations of cancer drug targets with various effects on other most influential components. We found that the combinatorial perturbation consisting of PI3K inactivation and overactivation of IP3R1 can lead to increased activity levels of apoptosis-related components and tumor-suppressor genes, suggesting that this combinatorial perturbation may lead to a better target for decreasing cell proliferation and inducing apoptosis. Finally, our approach shows a potential to identify and prioritize therapeutic targets through systemic perturbation analysis of large-scale computational models of signal transduction. Although some components of the presented computational results have been validated against independent gene expression data sets, more laboratory experiments are warranted to more comprehensively validate the presented results.

visualization, comprehension, and interpretation of big data in biomedical research. These fields provide an array of methodologies including computer simulations that can be used to generate new hypotheses and identify which hypotheses might be more productive to undertake experimentally, and eliminate hypotheses with little chance of success (Kitano, 2002a;Ghosh et al., 2011). These methods can be effective in navigating complex network problems associated with diseases. Many diseases and pathologies can be characterized by the dysregulation or dysfunction of multiple molecular components that are connected within these highly intertwined biological and biochemical networks (Loscalzo and Barabasi, 2011). Biological networks, including biochemical signal transduction networks, consist of a large number of highly interconnected pathways that give rise to complex, non-linear dynamics governing various cellular functions (Helikar et al., 2008;Helikar and Rogers, 2009). Disruptions of these networks, such as mutations or disease states can have drastic effects upon the whole system. These effects are difficult to predict from static network diagrams.
However, understanding the hierarchy of these changes remains a paramount problem. Often the specific causal interactions of the disease state are hidden within the massive cell-wide alterations, making attempts to reverse a disease state more challenging. In addition, the specific causal interactions are difficult to predict making the development of a potential therapeutic target results in unforeseen side effects (Singh and Singh, 2012). The unwanted effects of these drugs are often drastic as seen with many cancer medications (Kayl and Meyers, 2006;Lotfi-Jam et al., 2008;Singh and Singh, 2012). Therefore, it is necessary to gain a systems level understanding of the components associated with the disease states.
In recent years, targeted therapy has been used for multiple diseases, e.g. cancer (Vanneman and Dranoff, 2012), and often involve the activation or inactivation of a specific component in a biological network by a small molecule or drug, for instance. Perturbation analyses allow one to interrogate the structure and dynamic footprint of the underlying biological system. Perturbation biology has been proposed as an approach to reduce the collateral damage caused by nonspecific drugs. Computational network perturbations and new methods to evaluate the robustness of a given network can help identify more effective network components to target in order to obtain desired outcomes with minimal disruption to the rest of the network (Molinelli et al., 2013).
In order to fully leverage the potential of computational network perturbation analyses large dynamical models are necessary. A wide spectrum of modeling approaches exists ranging from detailed (but less scalable) differential equation-based systems to large (but not dynamic) static networks. In the middle are approaches such as logical modeling that are relatively scalable while capable of capturing the dynamic nature of biological systems (Le Novère, 2015). Logical networks, namely Boolean networks, have been used to describe and simulate a wide spectrum of biological systems ranging in size as well as contextual application (Naldi et al., 2010;Helikar et al., 2012;Madrahimov et al., 2013;Rocha et al., 2013;Conroy et al., 2014).
Herein, we present results from a system-wide perturbation analysis of a large-scale Boolean model of a signal transduction network widely present in many types of cells. Specifically, the model (previously described in (Helikar et al., 2008)) represents signaling events within the integrated epidermal growth factor, G-protein coupled receptor, and integrin signaling network. The model consists of 139 components (mostly proteins) and 557 biochemical interactions. The simulation-based, system-wide perturbation analyses enabled us to identify the most and least influential components (ones with the most and least impact on the rest of the network). To explore the role and effects of these perturbations in the context of the complex extracellular environment, the simulations and analyses were conducted under four biologically relevant environmental conditions known to stimulate cell growth, cell death, motility, and quiescence (in addition to a set of randomly generated environmental stimuli). In order to investigate potential therapeutic targets, we performed functional annotation and analysis of the most influential signal transduction components under both inactivating (e.g. knock-out) and activating (e.g. over-expression) perturbations. The most influential components were found to be enriched with many biological processes and druggable targets. Also, the most influential components under activating perturbations were enriched with more essential genes than the least influential components. We identified a network of the most influential components consisting of drug targets considered in multiple cancer types. The highest ranked among the most influential components were already explored as drug targets against cancer. We used all the most influential components and their upstream regulators to identify novel interactions. We used this approach to identify novel drug targets in the signal transduction network. As a result of the systemic analysis, we identified one novel combinatorial target, PI3K-IP3R1, with consistent occurrence in all simulated environmental conditions. We simulated the effect of combinatorial perturbation and the results were correlated with the literature, further supporting our predictions.

Computational Model
The computational model analyzed in this work is a Boolean model of signal transduction in a generic cell type previously described in (Helikar et al., 2008). The signal transduction model was constructed manually from around 500 published papers. The model consists of several main signaling pathways, including the receptor tyrosine kinase (epidermal growth factor receptor), G-protein coupled receptors (G-alpha i, G-alpha q, G-alpha s, and G-alpha 12/13) and the integrin signaling pathways. Each component in the model corresponds to a signaling molecule (mainly protein). The model also contains nine external components that represent the extracellular environment (mostly composed of receptor ligands). It is fully annotated and freely available for simulations and/or download via the Cell Collective software Helikar et al., 2013) at www.thecellcollective.org.
Each model component can assume an active (1) or inactive (0) state at any time t. The activity state of the model's internal components is determined by the regulatory mechanisms of other directly interacting components. These regulatory mechanisms are described with Boolean functions (in the form of truth tables or Boolean expressions). To represent the milieu of stimuli in the extracellular environment, the model contains external components that represent various ligands. The activity level of these components is specified as a probability to simulate different levels of concentrations. This methodology was previously detailed and exemplified in (Helikar et al., 2008;Helikar and Rogers, 2009;Helikar et al., 2012;Todd and Helikar, 2012).

Model Simulations
The Cell Collective platform was used to perform all computational simulations of the model. Although the model is built by using discrete mathematics the output activity levels (AL) can be continuous (ranging from 0 to 100) as previously described in (Helikar et al., 2008;Helikar and Rogers, 2009). Each simulation is synchronous and consists of 800 steps, where the activity level of the measured output component is calculated as the fraction of ones (active states) over the last 300 iterations that describe the network's steady behavior (Helikar et al., 2008;Helikar and Rogers, 2009).

∑ ∑
The model was simulated and analyzed under four biologically-relevant environmental conditions: cell growth, cell death, quiescence, motility (and randomly generated conditions). Each environmental condition was defined by different combinations of the activity levels of external components (ligands). The activity level ranges of the environmental conditions were determined by an optimization method whereby 2,000 simulations were run with all external stimuli ranging from 0-100 (except for IL1_TNF and Stress that were limited to low activity levels). Subsequently, specific environmental conditions for cell growth, cell death, motility, and quiescence were determined by selecting for environmental conditions that met the activity levels of biologically relevant model components: Akt, Erk, Cdc42, Rac (Helikar et al., 2008). From the initial 2,000 random simulations, environmental conditions that yielded the appropriate biological response were averaged and the inclusion of all environmental conditions within one standard deviation created an activity range for each of the environmental components (Table 1).
A wild type (WT) experiment was conducted under each environmental condition without any perturbations. Subsequently, systematic perturbation experiments were conducted under each condition, whereby each component of the model was constitutively activated (gain-of-function/overexpression) or inactivated (loss-of-function/knock-out). Each experiment consisted of randomly selecting 100 combinations of activity levels of the external stimuli from each condition activity range. (The only exception was the random environmental condition, which was simulated 2,000 times.) Each of the 100 combinations were simulated 30 times (i.e., 30 replicates) to ensure consistency of the dynamics in response to a specific combination of stimuli. These replicates were subjected to a Fligner Killeen test of homogeneity of variables which confirmed that the measured activity levels of the network components were homologous for identical combinations of activity levels of the environmental stimuli.

Model Analysis
The Kolmogorov-Smirnov (KS) test (Wang et al., 2003) was used to compare the WT dynamics (under each environmental condition) with the dynamics of each perturbation experiment. If the KS test resulted in a p-value less than 0.05, then it has a difference value (DV) equal to the test statistic; otherwise, the DV for a component is zero. In order to avoid skewed results from the perturbation itself, its difference value from the WT is set to zero.

Most and Least Influential Components
The most influential components are defined as components that induce the largest changes in the network under a given perturbations. The ranking of the perturbations is derived by calculating an influence score (IS), which is found by the summing of each DV for all the components per perturbation experiment. The top ten percent are considered most influential, and the bottom components with IS value 0 were considered the least influential.

Most Affected Components to a Specific Perturbation
For each perturbation induced, the components that are most sensitive to that perturbation are ranked in decreasing order to be able to characterize downstream effects of the perturbation on the network.

Annotation and Biological Relevance of Signal Transduction Components
All model components were first annotated using the appropriate NCBI gene IDs (Pruitt et al., 2007) for associated genes and UniProt IDs (Consortium, 2011) for protein products of the genes. All components were then further characterized using online resources such as DrugBank (Wishart et al., 2006).
Essentiality data were obtained from the Online GEne Essentiality (OGEE) database and mapped on the most and least influential components (Chen et al., 2012). DrugBank data were used to obtain druggability information for each component in the network. Data on cancer associated genes were obtained from The Cancer Genome Atlas (TCGA) (Weinstein et al., 2013) and mapped on the most influential components to identify cancer associated most influential components. The enrichment of essential genes and druggable proteins was computed based on the number of genes mapped on most or least influential components out of the total number of most and least influential components.

System-wide Perturbation Analysis Reveals Core Components of the Signal Transduction Network
A critical objective of biomedical research is the identification and prioritization of novel therapeutic targets. In this context, we performed systematic perturbation analysis in a generic signal transduction model. The workflow used in this work is illustrated in Figure 1.
The activating/inactivating perturbation experiments for each component in the model were carried out across four environmental conditions (as described in Methods). Additional randomly generated extracellular conditions were used to check the robustness of the model and results. Perturbation analysis enabled us to identify and rank components of the signaling network that are most and least influential (Supplementary Table 1). The heatmaps for all the environmental conditions (Supplementary Figures S1-S10) indicate that a few components had high influence on rest of the system. Therefore, we considered the top 10% of the components from each condition as the most influential. In contrast, the components that had no influence on the system were considered as the least influential (KS=0).
Also, the most influential components correspond to network components that, when perturbed, affect the largest part of the network in terms of the number of affected components and the magnitude of the effect. The most influential components were found for both inactivating (Figure 2 A) and activating (Figure 2 B) perturbations under the different environmental conditions. It is interesting to note that many of the most influential components overlap across all environmental conditions. However, the most influential components do not overlap between two types of perturbations (inactivating or activating).

Housekeeping Genes are Enriched in the Most Influential Components Common in Different Environments
Next, we investigated if the overlapped most influential components among different environmental conditions, have constitutive expression. Under inactivating perturbations, out of the seven components common among the different environmental conditions, PI4K, PI5K, ARF and PI3K were associated with housekeeping genes (Eisenberg and Levanon, 2013). Under activating perturbations, Trafs, Erk, Mek and SHP2 (out of nine common components), were associated with housekeeping genes. Housekeeping genes associated with the common components are displayed in Table 2. This observation suggests that the most influential components that are common among different environmental conditions are likely to function as housekeeping genes.

Unique components associated with each environmental condition are found to be condition specific
Under both types of perturbations, certain environmental conditions had several uniquely associated components ( Figure 2, Table 3). Under inactivating perturbations, components uniquely associated with the cell death condition are Calmodulin (CaM), RGS, and Palpha_iR. Out of these, CaM and RGS have been previously associated with cell death and apoptosis (Fisher, 2009;Berchtold and Villalobo, 2014). In fact, CaM plays a central role in the regulation of several cellular functions, including programmed cell death (Berchtold and Villalobo, 2014). It is also known that RGS protein can regulate cell death, cell cycle and cell division (Fisher, 2009). Under activating perturbations, the most influential components associated with the cell death-inducing condition include Gbg_i and Alpha_iR. On the other hand, PP2A was found to be most influential under the growth condition, Ras and Sos under motility condition, and PAK under quiescence condition. These results are also further supported by published studies that reported Gbg_i (GNB) to be involved in mTOR-mediated antiapoptotic pathways; Gbg_i was also functionally annotated with apoptosis and cell death (Wazir et al., 2013). PP2A was reported as a highly regulated Ser/Thr phosphatase involved in cell growth and signaling (Janssens and Goris, 2001). In pancreatic cancer cell lines, the knock down of KRAS has been found to lead to the decrease in cell motility and proliferation (Rachagani et al., 2011;Birkeland et al., 2012). Furthermore, the Grb2-Sos1 complex has been found to most likely promote cell motility, and tumerogenesis (Qu et al., 2014). These observations suggest that the proteins which were uniquely associated with simulated environmental conditions are most likely to have the association with that condition. Finally, the literature evidence obtained for housekeeping, or condition associated genes, further supports our simulation results.

Key Biological Processes are Enriched in the Most Influential Components
Next, we assessed the enrichment of biological processes or pathways in the most influential components. The most influential components across all four conditions under both types of perturbation showed significant enrichment with key biological processes. The counts and fold differences of enriched biological terms in all the conditions are shown in Figure 3 and 4. In the case of inactivating perturbations, inositol phosphate metabolism was enriched under all environmental conditions (Figure 3). In the case of activating perturbations, the significantly enriched biological processes include phosphate metabolic processes, kinase activity, apoptosis and interestingly, the nonsmall lung cancer pathway (Figure 4). These results illustrate that the group of proteins with similar biological functions appear as the influential components under each type of perturbation.

The most Influential Components under Activating Perturbations are Enriched with Essential Genes
Mutations in an essential gene can be lethal. Based on the hypothesis that the influential components might serve as essential for the survival of the cell, we performed essentiality analysis. Under activating perturbations, more essential genes were found within the most influential components than the least influential components (Figure 5 A). Under the cell death environmental condition, a total of 69% of the most influential components were essential; this is in contrast to the least influential components that contained 31% essential genes. Under other environmental conditions-growth, motility and quiescence, the difference of essential genes between the most influential and the least influential components are 23%, 15%, and 32% respectively.
On the other hand, under inactivating perturbations we found either an equal or larger number of essential genes in the least influential components (Figure 5 B). The most significant differences were observed under the cell death condition: the least influential components have 66% of essential genes in contrast to the 46% essential genes in the most influential. Also, under the growth condition 68% and 53% of essential genes were contained within the least and the most influential components, respectively. Under the motility and quiescence conditions, there were 3% and 9% more essential genes within the least influential components than the most influential components, respectively. We found that under inactivating perturbations, the number of essential genes among the least influential components were slightly larger than the activating perturbation (Figure 5 C and Figure 5 D). On the other hand, under activating perturbations, the more essential genes mapped within the most influential components than the least influential components.
Thus, the most influential components are essential under activating perturbations, suggesting an environmental condition-specific essentiality.

The Most Influential Components are enriched with Druggable Proteins
To further investigate the importance of the most influential components, we evaluated the distribution of known druggable targets. We obtained druggability data from the DrugBank database (Wishart et al., 2006) and mapped them on the most and least influential components. A total of 51 components in the whole network were enriched with druggable proteins. We found that under both types of perturbations and across all environmental conditions more druggable proteins were found among the most influential than the least influential components ( Figure 6). The enrichment for druggable proteins within the most influential components implicates these as critical network components.

Ranked most influential components based on downstream components
We identified the most affected components of the most influential components under both types of perturbations. We combined all environmental conditions to construct networks of the most influential components with their downstream targets. We subsequently mapped druggable proteins and cancer associated genes on these networks. Under inactivating perturbations, we obtained a network consisting of the most influential components: PI3K, EGFR, PP2A, GRK and CaM (Figure 7 A). Under activating perturbations, we obtained a network composed of influential components: EGFR, IL1_TNFR, ERK, SHP2, RKIP, Ras, Gbg_i, Fak, Integrins, and PP2A (Figure 7 B).
Total number of downstream targets for each of the most influential druggable component under both inactivating and activating perturbations are listed in Table 4. We observed that EGFR, a validated cancer drug target (Mendelsohn, 2001), affects the largest number of components under activating and inactivating perturbations.

The Most influential components mainly affect other most influential components
Here, we identified all components that directly affect the activity of each most influential component (KS=1). Interestingly, most of these direct upstream components were also ranked as the most influential in at least one environmental condition (Figure 8). Under inactivating perturbations, 22 components were directly upstream of the most influential components. Of these, 19 were the most influential under at least one environmental condition. On the other hand, under activating perturbations, out of 45 upstream components, 19 were also ranked as most influential. Additionally, under inactivating perturbations, nine (CaM, EGFR, Gbg_i, GRK, IP3R1, PP2A, PI3K, Ras, and Src) out of total 22 upstream components are druggable. Out of these 22 components, six components (CaM, EGFR, Gbg_i, GRK, IP3R1, and PP2A) were upstream to the most influential druggable components. Under activating perturbations, 21 (CaM, Cdc42, EGFR, Erk, Fak, Gbg_i, Grb2, GRK, IL1_TNFR, Integrins, IP3R1, PDK1, PI3K, PKA, PP2A, Rac, Raf, Ras, RKIP, SHP2, and Src ) out of 45 upstream to the most influential components are associated with druggable proteins. Out of these 21, ten components were also the most influential. Under both types of perturbations, a total of 18 (alpha_iR, ARF, B_Arrestin, Ca, CaM, EGFR, Gbg_i, GRK, IP3R1, Palpha_iR, PI5K, PIP2_45, PIP3_345, PP2A, RGS, PI3K, Ras, Src) upstream components were common. Nine of these components (CaM, EGFR, Gbg_i, GRK, IP3R1, PP2A, PI3K, Ras, and Src) were druggable or these were used as the drug targets. The important drug targets such as EGFR, PI3K, Ras, Raf also appeared as influential upstream components. Together, these results suggest that under inactivating perturbations the activity of the most influential components are likely to be modulated by the other most influential components.

The Most influential Components as Drug Targets, and Drug Resistance
The top most influential components such as EGFR, PI3K, ERK, and Ras etc. have been previously explored as drug targets in multiple cancer types. However, it is also evident from literature that several most influential components have been associated with drug resistance. For example, in non-small cell lung cancer, mutation within the kinase domain of EGFR, and epithelial-mesenchymal transition are responsible for the development of resistance to gefitinib (Holohan et al., 2013). In colorectal, and head and neck cancers, KRAS mutation, EGFR-S492R mutation, and increased ErBb signaling are responsible for resistance against Cetuximab (Dienstmann et al., 2012;Holohan et al., 2013). Furthermore, PI3K showed drug resistance in breast cancer against rapamycin through the expression of RSK3 and RSK4 (Rodon et al., 2013). Mutations in ERK1 or ERK2 have shown resistance against ERK inhibitors or RAF/MEK inhibitors (Wagle et al., 2014). Tumors with mutation in BRAF V600E can adapt to the RAF inhibitors (Lito et al., 2013;Perna et al., 2015). As such, the identification and prediction of drug targets alone is not sufficient to identify completely useful drug targets.

Regulatory Interactions between the Most Influential components and their Upstream Components
To develop a better strategy that can account for drug resistance of the most important drug targets, we sought to investigate novel regulatory interactions. We analyzed the previously described interactions between the most influential components and their direct upstream components. We found that some interactions consistently occur in more than one environmental condition. For example, the inactivation of IP3R1 increases the activity of PI3K under all four environmental conditions. However, the maximal effect was observed under the death environmental condition. Additionally, the inactivation of IP3R1 leads to inactive RGS under three conditions: cell growth, motility, and quiescence. These finding also correlate with published studies that found that RGS positively regulates apoptosis (Fisher, 2009). Other examples of consistently occurred interactions include: the activation of Grb2 leads to increased activity levels of Ras under all four environmental conditions, and increased Sos activity under death and quiescence conditions. The activation of Rac increases the activation of PAK under cell death and growth conditions. Overall, we found three types of interactions: inactivation of one component leads to the increase of activity of another component (PI3K-IP3R1, IP3R1-PI3K, RGS-IP3R1), inactivation of a component leads to decreased activity of another component (IP3R1-RGS), and activation of a component leads to increased activity of another component (Grb2-Ras, Grb2-Sos, Rac-PAK).
The fold differences of all these interactions are displayed in the Table 5. Under the cell death condition, the inactivation of IP3R1 results in PI3K activity increase by 2.38 fold. Similarly, PI3K inactivation leads to a 5.42 fold increase in IP3R1 activity. In the case of other interactions, the inactivation of IP3R1 leads to inactive RGS under the cell growth, motility and quiescence conditions. Under the motility and quiescence conditions, the inactivation of Gbg_i leads to inactive CaM. The activation of Grb2 increases the activity of Ras 7.40 fold under the cell death condition, and 2.13 fold under the quiescence condition.Grb2 activation also affects Sos 7.8 fold under the cell death condition and 2.18 fold under the quiescence condition. An activating perturbation of Rac increases the activity of PAK more than 18 fold under the cell death condition, and 5.59 fold under the growth condition.
These results suggest different types of regulatory effects of activating and inactivating perturbations of direct upstream components of most influential components.

Co-targeting IP3R1 with PI3K
As discussed earlier, although PI3K was identified as one of the most influential components, it has been also associated with drug resistance. Based on the interactions of upstream regulators of most influential components discussed above, we further investigated the interactions involving PI3K and IP3R1 with the objective of identifying a secondary drug target that could be potentially used to address the issue of PI3K-associated drug resistance. In contrast to PI3K/Akt signaling, IP3R1 positively regulates apoptosis. We hypothesized that the rate of apoptosis will increase when IP3R1 is overactivated (activating perturbation) and PI3K is inactivated (inactivation perturbation). Despite the strong dynamical relationship between IP3R1 and PI3K, these two components are only connected indirectly through a sub-network. In this sub-network, Gbg_i is upstream of and directly activates both components. IP3R1 regulates PI3K through a Ca->EGFR route, whereas PI3K regulates IP3R1 via a PTEN -> PIP2_45 -> IP3 route (Figure 9).
The inactivating perturbation of PI3K resulted in the inactivation of 29 components across all four environmental conditions. Under PI3K inactivation, the average activity of IP3R1 increased from 71.9% in wild type to 85.18%. This perturbation also led to down-regulation of positive regulators of apoptosis phospholipase A2 (PLA2) and arachidonic acid (AA). AA released by PLA2 triggers Ca2+ dependent apoptosis through mitochondrial pathways (Penzo et al., 2004). The elevation in Ca2+ is thought to be involved in apoptosis (Pinton et al., 2008). It was shown that blocking calcium channels can directly lead to tumor promotion (Mason, 1999). Thus, inactivation of PI3K can block cell proliferation; simultaneously it can lower the rate of apoptosis.
Under growth condition, the activating perturbation of IP3R1 increased the activity of apoptosisassociated components: Ca, CaM, CaMK, CaMKK and RGS in the range of +1.41 Fold to +2.09 fold when compared to wild type.
To simulate the cell death effect under growth condition, we carried out a double perturbation of IP3R1 and PI3K whereby we constitutively activated IP3R1 and inactivated PI3K1. Under this combinatorial perturbation, we found 27 proteins including proto-oncogenes such as Akt (which suppresses apoptosis) and Raf to be down-regulated. Here, we found eight proteins with more than 19% increased activity than in the case of a single inactivating perturbation of PI3K. These proteins include Rap1 (+ 1.19 Fold), Ca (+ 1.21 Fold), CaM (+1.21 Fold), CaMKK (+ 1.21 Fold), Myosin (+ 1.22 Fold), CaMK (+ 1.65 Fold), PLA2 (+ 1.98 Fold), and AA (+1.98 Fold) ( Table 6; full list of all affected components is given in Supplementary Table 2). It is noteworthy that these components have been found to be associated with apoptosis or cell death. Therefore, under combinatorial perturbations, components involved in cell proliferation were downregulated through the inactivation of PI3K, and the activity of tumor-suppressor genes (PLA2) with arachidonic acid (AA) and other components, including Ca, CaM, and CaMK, was increased as a result of the IP3R1 overactivation.
Together, these results suggest a novel regulatory interaction between PI3K and IP3R1, and that cotargeting both of these components may serve as therapeutic strategy rather than targeting PI3K alone.

Discussion
We have presented a systemic perturbation analysis of a signal transduction network model to identify and characterize functionally important components. We used these components to explore novel therapeutic strategies against cancer. Specifically, we used a logical modeling approach to analyze the dynamics of a large-scale signal transduction model. Logical modeling approaches have been used, for example, to understand the dynamics of signal transduction and gene regulation networks to identify drug synergies in gastric cancers, and to identify potential drug combinations (Flobak et al., 2015). In biochemical networks, combined effect of topology and dynamical features have been shown to have the most significant impact on the dynamics of the network (Kochi et al., 2014). Computational approaches have become indispensable tools to understand biological pathways and disease phenotypes. Examples include computational methods such as molecular modeling, text mining, and network modeling to identify drug targets in a vast array of diseases from pathogens to complex disorders (Flórez et al., 2010;Yao et al., 2010;Folger et al., 2011;Madrahimov et al., 2013;Puniya et al., 2013).
In the present work, the identified most influential components were characterized for biological functions. The relevance of identified influential components was established with pathway analysis, mapping of housekeeping genes, essential proteins, and association with druggable proteins. Interestingly, we found enrichment of housekeeping genes in the most influential components that were independent of the extracellular environments. A notable agreement is obtained from literature surveys for the most influential components, which were unique to specific environmental conditions. Because essential components are important from a disease perspective, the identified most influential components may serve as potential candidates and essential proteins under specific conditions. Under activating perturbations, we found that essential genes were enriched more within the most influential components than within the least influential components. The high association of dysregulated signal transduction proteins with different subtypes of cancers suggests that these components may be important candidates for drug targets. Notably, the most influential components are enriched with several already known drug targets. However, many of these drug targets (EGFR, ERK, Ras, PI3K etc.) have been associated with drug resistance (West et al., 2002;Kobayashi et al., 2005;Linardou et al., 2008;Wheeler et al., 2010;Dienstmann et al., 2012). The mechanism of drug resistance includes mutation in the targeted protein or expression of other genes (altered expression) to bypass the effect caused by perturbation, and deregulation in apoptosis, etc. (Gottesman, 2002;Holohan et al., 2013).
Thus, to identify novel regulatory interactions, we explored components that are upstream to the most influential components associated with drug resistance. Interestingly, several upstream components (more than 90% in the case of inactivating perturbations) to the most influential components were also identified as most influential. Thus, the most influential components form a tightly connected subnetwork of proteins interacting with each other. In yeast, it has previously shown that the essential proteins are hubs in the network and have more interconnections than non-essential proteins, and form a module or sub-network (Song and Singh, 2013).
The interaction between IP3R1 and PI3K was observed under all environmental conditions. IP3R1 activation, when combined with PI3K inactivation, increases the activities of PLA2 and AA, which are decreased with a single PI3K knockdown. It was already shown that AA released by PLA2 helps to initiate apoptosis (Penzo et al., 2004). In a Dictyostelium discoideum chemotaxis experiment, it was also shown that cells with PI3K deficiency were more sensitive to PLA2 inhibition (Chen et al., 2007), which supports our predicted interaction between PI3K and PLA2. To this end, we hypothesized that the PI3K inactivation could be combined with the over-activation of IP3R1 to increase the activity of proteins involved in apoptosis. IP3R1 inactivation can lead to the downregulation of RGS, and reversibly, the overexpression of IP3R1 can lead to increased activity of RGS. Similar to IP3R1, RGS subtype RGS3T has been found to be involved in inducing cell death (Fisher, 2009), and it has also been found that RGS can suppress the PI3K activity downstream of the receptor (Liang et al., 2009). Therefore, the constitutive activation of IP3R1 might also negatively regulate the activity of PI3K. Systemic analysis of the most influential components and their upstream components has led us to identify novel combinations of drug targets. In various studies, combinatorial therapies have shown a decrease in drug resistance in pathogens. In combinatorial therapy, a protein associated with drug resistance can be targeted in combination with different protein of either the same or different pathway (Fischbach, 2011). Clinical trials have also suggested that the efficiency of cytotoxic drugs increases when given in combinations (Al-Lazikani et al., 2012).
In conclusion, by combining IP3R1 (activation) and PI3K (inactivation), we were able to stimulate cell death under the cell growth condition. Based on this, one can hypothesize that it might be possible that the decrease in cell proliferation with increased apoptosis as a result of this combinatorial intervention could subsequently increase the rate of clearance of tumor cells, and serve as a novel strategy for important targets associated with drug resistance. However, more laboratory validations will be required to test this hypothesis.

Acknowledgments
We would like to thank Resa Helikar for providing feedback on the manuscript.