A predictive computational model reveals that GIV/girdin serves as a tunable valve for EGFR-stimulated cyclic AMP signals

INTRODUCTION

Cells constantly sense cues from their external environments and relay them to the interior. Sensing and relaying signals from cell-surface receptors involves second messengers such as cyclic nucleotides (Beavo and Brunton, 2002; Newton et al., 2016). Of the various cyclic nucleotides, the first to be identified was cyclic adenosine 3,5-monophosphate (cAMP), a universal second messenger. cAMP relays signals triggered by hormones, ion channels, and neurotransmitters (Sassone-Corsi, 2012) and also binds and regulates cAMP-binding proteins such as PKA and Epac1 (Sutherland and Rall, 1958).

Intracellular levels of cAMP are regulated by the antagonistic action of two classes of enzymes: adenylyl cyclases (ACs) and cyclic nucleotide phosphodiesterases (PDEs). ACs are membrane-bound enzymes that utilize ATP to generate cAMP. PDEs, on the other hand, are soluble enzymes that catalyze the degradation of the phosphodiester bond, resulting in the conversion of cAMP to AMP. Overall, the level of cellular cAMP in physiology is finely balanced between synthesis by AC, degradation by PDEs, and regulatory feedback loops from PKA to PDE (Tremblay et al., 1985; Sette and Conti, 1996; MacKenzie et al., 2002; Murthy et al., 2002; Sassone-Corsi, 2012) or others that act on ACs and PDEs (Yoshimasa et al., 1987; Bruce et al., 2003; Goraya and Cooper, 2005). Too much or too little cAMP is seen in many diseases. For example, high levels of cAMP have been shown to be generally protective in diverse cancers (e.g., breast, Tagliaferri et al., 1988; melanoma, Dumaz et al., 2006; pancreas, Boucher et al., 2001), whereas low cAMP levels fuel cancer progression (reviewed in Fajardo et al., 2014). This is because cAMP inhibits several harmful properties of tumor cells such as proliferation, invasion, stemness, and chemoresistance, while enhancing differentiation and apoptosis. Although drugs targeting the canonical GPCR/G-protein-cAMP signaling pathway have successfully translated to the clinic for tackling a wide range of ailments (Filmore, 2004), from hypertension to glaucoma, such strategies have largely failed to impact cancer care or outcomes. Thus, how tumor cells avoid high levels of cAMP appears to be incompletely understood, and therapeutic strategies to elevate cAMP remain unrealized.

Recently, the regulation of cAMP by noncanonical G-protein signaling that is initiated by growth factors (Ghosh, 2015, 2016; Ghosh et al., 2017; Aznar et al., 2016) has emerged as a new signaling paradigm. Growth factor signaling is a major form of signal transduction in eukaryotes, and dysregulated growth factor signaling (e.g., copy number variations or activating mutations in RTKs, increased growth factor production/concentration) is often encountered in advanced tumors and is frequently targeted with various degrees of success (Lowery and Yu, 2012). A body of work published by us and others has revealed that RTKs bind and activate trimeric G proteins via a family of proteins called guanine-nucleotide exchange modulators (GEMs; Ghosh et al., 2017). GEMs act within diverse signaling cascades and couple activation of these cascades to G-protein signaling via an evolutionarily conserved motif of ∼30 amino acids that directly binds and modulates Gα_i and Gα_s proteins. Most importantly, GIV-GEM serves as a GEF for Gα_i and as a GDI for Gα_s (Gupta et al., 2016). Despite this apparent paradox, both forms of modulation lead to suppression of cellular cAMP (Ghosh, 2016). By demonstrating how GIV, a prototypical member of a family of cytosolic guanine-nucleotide exchange modulators (GEMs; Gupta et al., 2016; Ghosh et al., 2017), uses a SH2-like module (Lin et al., 2014) to directly bind cytoplasmic tails of ligand-activated RTKs such as EGFR (Ghosh et al., 2010), we provided a definitive structural basis for several decades of observations made by researchers that G proteins can be coupled to and modulated by growth factors (reviewed in Marty and Richard, 2010). A series of studies since have revealed that growth factor–triggered noncanonical G protein→cAMP signaling through GIV has unique spatiotemporal properties and prolonged dynamics that are distinct from canonical GPCR-dependent signaling (reviewed in Aznar et al., 2016). In parallel, studies have also found that high levels of GIV expression fuel multiple ominous properties of cancer cells, such as invasiveness, chemoresistance, stemness, and survival, and are associated with poorer outcome in multiple cancers. Inhibition of GIV’s G protein–modulatory function has emerged as a plausible strategy to combat aggressive traits of cancers (reviewed in Ghosh, 2015; Ghosh et al., 2017; Ma et al., 2015). These findings provide us a unique opportunity to investigate, from a systems level, how modulation of trimeric GTPase Gα_i and Gα_s by GIV downstream of growth factors regulates cAMP and what impact such regulation might have on the aggressiveness of cancers.

In this study, we develop a mathematical model of cAMP signaling that is triggered by ligand stimulation of EGFR, and investigate how cAMP dynamics in cells is affected by the GIV-Gα_i/s and PDE axes. Further, we also seek to connect findings from the cell-based model to survival data from cancer-afflicted patients by identifying the most consequential variables within this signaling pathway. In doing so, this model not only interrogates the cross-talk between two of the most widely studied eukaryotic signaling hubs (RTKs and G proteins), but also reveals surprising insights into the workings of GIV-GEM and provides a mechanistic and predictive framework for experimental design and clinical outcomes.

RESULTS

Phenomenological model reveals that GIV-associated timescales modulate cyclic AMP dynamics

The emerging paradigm of noncanonical modulation of Gα_i/Gα_s proteins by growth factor RTKs comprises several temporally and spatially separated components (Figure 1A). We first developed a phenomenological model to identify the network topology of RTK–G protein–cAMP signaling (Figure 1B). This network captures the key events of the steps shown in Figure 1A. Briefly, receptor (R) stimulation is modeled using a time-dependent function S(t) to result in active receptor R^*. R^* then acts on GIV on two time scales—τ₁ for GIV-GEF activation and τ₂, for GIV-GDI activation. cAMP synthesis is directly proportional to the level of internalized endosomal receptor R* with a time scale of τ₃. GIV-GDI inhibits the internalized receptor and GIV-GEF inhibits cAMP production. Even though none of the components in this model reflects actual biochemical species, the simplified model has the advantage of capturing the key time scales of the events leading up to cAMP production from RTKs. Varying these time scales alters the dynamics of GIV-GEF, GIV-GDI, and cAMP (Figure 1, C and D, and Supplemental Figure S1E). Additionally, this model has the advantage of being parameterized by a small number of variables and parameters (see the Supplemental Material for details). Simulations from this model predict that GIV-GEF activation is rapid, whereas GIV-GDI activation is slow (Figure 1C). This temporal response is observed for a wide range of parameters with the internalization and degradation rates, k₄, k₅, as leading-order contributors across most time points (Supplemental Figure S1). τ delays are shown to make smaller contributions to cAMP response then the major rates, k₄, k₅ (Supplemental Figure S1). cAMP dynamics predicted from this simple model shows a delayed increase in cAMP corresponding to the time scale of GIV-GDI in Figure 1D. Furthermore, changing the receptor density shows that the presence of GIV suppresses the cAMP production; the role played by GIV in modulating cAMP is stronger at higher RTK concentrations because of the competing effects of Gα_s and Gα_i. This leads us to two conclusions: First, the network topology in the toy model, with GIV as the central regulator, is able to capture cAMP dynamics. Second, RTK copy number alone is an incomplete determinant of cAMP; RTK and GIV together determine cAMP concentrations.

Construction and experimental validation of a compartmental model for noncanonical G-protein signaling triggered by growth factors

Although the phenomenological model in Figure 1 allowed us to identify key features of RTK–G protein–cAMP signaling, it does not contain enough information to compare simulation output against experimental measurements. Therefore, the topology model was expanded to a larger biochemical reaction network so that the modules reflected the timescales τ₁, τ₂, and τ₃ within a larger network model (Figure 2A). The model consists of four modules—Module 1 focuses on the well-established dynamics of EGFR (Berkers et al., 1991; French et al., 1995; Schoeberl et al., 2002; Supplemental Table S4); Module 2 represents the dynamics of the formation of the EGFR·GIV Gα_i complex, representing τ₁ (Supplemental Table S5; within this complex GIV-GEM serves as a GEF for Gα_i); Module 3 represents the dynamics of the formation of the Gα_s·GIV-GDI complex, representing τ₂ for GEF-to-GDI conversion (Supplemental Table S6; within this complex GIV-GEM serves as a GDI for Gα_s); and Module 4 represents the dynamics of cAMP formation and represents τ₃ for Gα_s activation by internalized receptors (Supplemental Tables S7 and S8).

Within each module, the biochemical reaction network includes several known interactions curated from the literature along with kinetic parameters (see Supplemental Tables S2–S6 for references). However, to our knowledge, no mathematical models of GIV-GEM interactions within the RTK→G protein→cAMP pathway exist. Therefore, we had to estimate kinetic parameters for certain reactions in each module. Of the 76 kinetic parameters in the model, 57 parameters were from the literature while 19 were fitted to experimental data. We fitted the model to the experimentally measured dynamics of the EGFR·GIV Gα_i complex (Figure 2B) and the Gα_s·GIV-GDI complex (Figure 2C). Because the actual concentration of this complex in cells is not known, and is likely to vary from cell to cell, we analyzed peak times and fold change of these complexes. The temporal dynamics of the normalized densities of both these complexes generated from simulations were in good agreement with experimental measurements, as determined by PLA and GST pull-down assays (Bhandari et al., 2015; Gupta et al., 2016) carried out in HeLa cells responding to EGF. We provide a detailed discussion of parameter estimation, goodness of fit, and uncertainty quantification in later sections. Our choice of experimental assays for validating the model and fitting parameters were carefully chosen in consideration of the strengths and weaknesses of each assay (see Materials and Methods). We used this modular network to investigate the role played by GIV in regulating the dynamics of EGFR, EGFR·GIV Gα_i complex, Gα_s·GIV-GDI complex, and cAMP.

Biological prediction: compartmentalized modulation of Gα_i and Gα_s by GIV-GEM governs EGF-triggered cyclic AMP dynamics

Because EGF/EGFR triggers activation of Gα_i at the PM first, followed by activation of Gα_s on the endosomes later, production of cAMP must be a balance between the antagonistic actions of these two G proteins on membrane-bound ACs (Supplemental Table S7). We assumed that the PM pool of Gα_i inhibits the AC→cAMP pathway at the PM. Similarly, because Gα_s is activated predominantly on endosomes and endosomal ACs (eACs) can be stimulated to synthesize cAMP locally (Vilardaga et al., 2014), we assumed that the endosomal pool of Gα_s likely stimulates the eAC→cAMP pathway (Supplemental Table S7). To capture the dynamics of cAMP in our model network, we included such compartmentalized G protein–AC interactions.

Our model predicted that the early inhibition of cAMP is due to the Gα_i-mediated inhibition of AC (the green regime; Figure 2D); cAMP production is increased later due to the activation of Gα_s on the endosome (the blue regime; Figure 2D). This dynamics is consistent with previously published GIV-dependent cAMP dynamics, measured by FRET (Gupta et al., 2016). While activation of GIV-GEF occurs earlier (within 5 min) at the PM, conversion of GIV-GEF to GIV-GDI occurs later (15–30 min) when EGFR is already compartmentalized in endosomes (Figure 1A); such temporally separated compartmentalized modulation of two G_α-proteins with opposing effects on AC ensures suppression of cAMP both early and later during EGF signaling (Gupta et al., 2016). Because GIV modulates both Gα_i and Gα_s in different compartments and on different time scales, the model predicts that increasing GIV concentration should dampen overall cAMP response to EGF and that decreasing GIV concentration should do the opposite (Figure 2D; compare GIV = 0.01 μM with GIV = 10 μM). These predictions were validated experimentally by measuring cAMP in control (shControl) and GIV-depleted (shGIV; >95% depletion by band densitometry, see Supplemental Figure S2E) cells at various time points after EGF stimulation (Figure 2E) by a radioimmunoassay (RIA). We found that in comparison with control cells, cellular levels of cAMP were always higher after EGF stimulation, at both early and late time points. In fact, when superimposed, the model and experiment showed good agreement throughout 60 min (Supplemental Figures S3, S7, and S8). As expected, sensitivity analyses confirmed that cAMP is sensitive to the initial concentrations of PDE and AC, and the reaction rates are associated with AC, internalization, PKA, and PDE (Supplemental Tables S14 and S17).

To distinguish what might be the relative contributions of the two G protein–modulatory functions of GIV (GEF versus GDI) on cAMP production, we investigated cAMP dynamics under three conditions (Figure 2D)—1) GEF-deficient but GDI-proficient (mimicked experimentally by the GIV-S1764D/S1689D mutant GIV-DD; Gupta et al., 2016); 2) both GEF- and GDI-deficient (mimicked experimentally by the GIV-F1685A mutant GIV-FA; Garcia-Marcos et al., 2009; Gupta et al., 2016); and 3) GEF-proficient but GDI-deficient (an in silico mutant, because there is no known mutant yet that can mimic this situation in experiments). In the first scenario, where GIV’s GEF function is selectively lost, but GDI function is preserved, increase in cAMP concentration occurred early (Figure 2D, dashed cyan line), as observed previously in cells expressing the GIV-DD mutant (Gupta et al., 2016). In the second scenario, where both GEF and GDI functions were lost, increase in cAMP concentration occurred early and such elevation was sustained (Figure 2D, dashed dark green line), as observed previously in cells expressing the GIV-FA mutant (Gupta et al., 2016); this mirrored the profile observed in GIV-depleted cells (Figure 2D, solid green line). Finally, in the third scenario, selective blocking of GIV’s GDI function using an in silico mutant resulted in an early decrease followed by a prolonged increase in cAMP concentration (Figure 2D, dot–dashed blue line).

While the dynamics of cAMP production provides insight into how different conditions lead to changes in concentration, the area under the curve (AUC) for cAMP concentration provides information critical for decision making, buffers time-scale variations, and averages the effect of fluctuations in concentrations (Atay and Skotheim, 2017). AUCs for cAMP at different time points were calculated to investigate how the cumulative cAMP signal varies under different GIV conditions (Figure 2F). For the control bars (in orange), we observe that at the 5-min time point, the AUC is negative. This represents the initial decrease in cAMP concentration. The AUC becomes positive and increases by 15 min, signifying a net accumulation of cAMP. The AUCs look similar in the GIV-FA mutant (defective in both GDI and GEF functions) as well as in the absence of GIV; that is, they increase progressively for 60 min (Figure 2F; compare the light green and dark green bars). If GIV levels are increased (10 μM, red bars), the AUCs remain negative throughout, showing the sustained nature of the dampening effect of GIV on cAMP. This dampening effect on cAMP is achieved primarily via activation of Gα_i in the short term (GEF regime) and via inhibition of Gα_s in the long term (GDI regime). These findings are in keeping with the previously suspected role of GIV in reducing cellular cAMP, but the model reveals that the relative contributions of GIV’s GEF and GDI functions are separated in time and space at a resolution that is experimentally unachievable.

Biological prediction: GIV dampens EGF/EGFR-triggered cyclic AMP production

Common wisdom from canonical signaling suggests that an increase in stimulus through receptor copy number leads to a proportional increase in cAMP concentration. Therefore, we would expect that an increase in EGFR density would lead to an increase in cAMP concentration. But the phenomenological model indicated that cAMP concentration depends on both the receptor copy number and the GIV concentration (Figure 1D). We investigated how GIV concentration affects cAMP dynamics with various EGF/EGFR numbers. When GIV concentrations were set to 0.01 μM in the model (to simulate cells that do not have GIV), increased input signals triggered increased output signals (Figure 3A; Supplemental Figure S10); this effect was even more pronounced in the absence of PDE (Figure 3E). This proportional response was lost when GIV concentrations were set to high levels (GIV = 10 μM; Figure 3B); that is, increased input signals failed to initiate proportional output signals. This effect was virtually unchanged and robustness was preserved despite the absence of PDE (Figure 3F). These effects are also evident from comparing the AUCs under the simulated conditions (Figure 3, C and G).

To test these predictions, we measured cAMP by RIA in control and GIV-depleted HeLa cells as in Figure 2E, except that in this instance we measured only at 60 min, but with various doses of EGF (experimental equivalent of variable input in simulations). To recapitulate simulations in the presence or absence of PDE, assays were carried out in parallel in the presence or absence of 3-isobutyl-1-methylxanthine (IBMX), an inhibitor of PDE (Figure 3, D and H). In the presence of GIV, cAMP production is robustly suppressed in response to increasing EGF ligand (Figure 3G). In the absence of GIV, cAMP production is sensitive to increased EGF, an effect that is further accentuated when PDEs are inhibited with IBMX (Figure 3K). Taken together, these results indicate that GIV primarily serves as a dampener of cellular cAMP that is triggered downstream of EGF. Unlike PDE, which reduces cellular cAMP by degrading it, GIV does so by fine tuning its production by G proteins and membrane ACs. Thus, our model identified that GIV is a critical determinant of cAMP concentrations in response to EGFR signaling.

Model reveals that GIV-GEM may serve as a tunable valve for cAMP flux: high GIV implies low flux, whereas low GIV implies high flux

To quantify the extent of cross-talk between EGFR and GIV, we conducted simulations for a wider range of EGFR [36–1800 molecules/μm²] and GIV concentrations [0.01–10 μM] and calculated the AUC for the cAMP dynamics (Figure 4, A and B; Supplemental Figure S11). In a low-EGFR state, varying GIV concentrations resulted in cAMP changes only within a narrow range; however, in a high-EGFR state, varying GIV concentrations achieved a larger variance in cAMP (Figure 4A). To further dissect this space, we plotted the variations in AUC for EGFR and GIV variations (Figure 4B). The value of AUC corresponding to the control (GIV 1μM and EGFR 240 molecules/μm², Supplemental Table S11) is ∼0.45 μM min and is denoted by the yellow color and is marked as a black solid line for different EGFR and GIV concentrations in the heat map; elevated cAMP level is denoted by green and reduced cAMP by red. We observed that increasing EGFR increased cAMP AUC in the setting of low GIV concentrations. But when GIV concentrations are high, cAMP AUC remained low regardless of increasing levels of EGFR, indicating that the impact of increasing GIV on cAMP AUC was higher than the impact of increasing EGFR. Therefore, GIV levels, in conjunction with EGFR levels, can be thought of as key determinants, and high GIV in the setting of high EGFR may facilitate tonic suppression of cAMP levels regardless of pathway stimulation.

Another factor that plays an important role in regulating cAMP levels is PDE. We next conducted simulations for different GIV (0.01–1 μM) and PDE (0.04–2.5 μM) concentrations to identify how the cross-talk between these two variable components regulates cAMP levels. Within each category of PDE concentration (low vs. high PDE states; Figure 4C; Supplemental Figure S12), cAMP levels are highest when GIV levels are lowest, and vice versa. In the setting of low PDE activity, the impact of changing GIV was highest, that is, the range of cAMP response was widest. In contrast, in the setting of high PDE activity, the impact of changing GIV on the cAMP levels was minimal. These effects can be seen when the AUCs for the low and high PDE states, calculated over 1 h, are compared (Figure 4C). While there is no significant change in the AUC with increasing GIV in a high-PDE state (red bars), an increase in GIV leads to a decrease in cAMP in the low-PDE state (green bars). That is, for a given GIV concentration, the effect of PDE is always stronger. Furthermore, a heat map of cAMP AUCs (Figure 4D) shows the interplay between PDE and GIV concentrations over a wide range. For low PDE concentration, increasing GIV decreases cAMP AUC, but the cAMP AUC is well above the yellow value (marked as control). However, an increase in PDE concentration leads to a dramatic decline in cAMP AUC even when GIV levels are low; this condition is likely to result in futile cycling (high cAMP production due to low GIV and high cAMP clearance due to high PDE signaling). Together, these findings indicate that the effect of GIV concentration on cAMP levels in cells is discernible only when PDE activity is low. Because the high-PDE state virtually abolishes all effects of GIV-dependent inhibition of cAMP production, we also conclude that in this GIV-PDE cross-talk, PDE is a dominant node and GIV is the subordinate node.

Our network model has helped us identify key design principles of the action of GIV-GEM within the EGF/EGFR signaling circuit by enabling construction of a map to identify the relationship between the key components—input (EGFR)→valve (GIV)→output (cAMP)→sink (PDE; Figure 4E), validate the impact of this relationship by experimental assessment of cellular cAMP (Figures 2 and 3), and interpret the role of each component in the context of network architecture. That there is complex, nonlinear, and nonintuitive cross-talk between EGFR, GIV, and PDE in regulating cAMP levels is evident from the fact that the isoplanes, which capture the same cAMP AUC, are not flat but are bent surfaces in this space. It appears that variation of cellular concentration of functionally active GIV-GEM molecules serves as the most tunable component that regulates the flow of signal from EGF/EGFR (input) to cAMP (output; Figure 4E). At low concentrations of GIV, such as those found in normal tissues, cAMP levels are sensitive to increased signal input via EGF/EGFR; that is, higher input elicits higher output. Such sensitivity is virtually abolished and replaced by robustness at higher GIV concentrations found in a variety of cancers; that is, higher input fails to elicit higher output, and instead, cAMP levels stay low and relatively constant. This three-way interplay between EGFR, GIV, and PDE is obvious also in experimental data derived from HeLa cells (Figure 3, D and H) suggesting that GIV acts a tunable valve for the input–output relationships that govern RTK–G protein–cAMP signaling.

Clinical predictions from the model—from math to man

The signaling network model built using dynamic protein–protein interactions and enzymatic reactions during signal transduction predicts that GIV dampens cAMP signaling in response to EGF/EGFR stimulation over a 60-min time course (Ghosh, 2015, 2016; Ghosh et al., 2017; Aznar et al., 2016). We chose to stick to 60 min because this is the time period in which almost all the tyr-phosphorylated EGFR pool disappears, which coincides with the entry of EGFR into late endosomes, where 60% total EGFR remains while the remaining 40% is degraded in the endosomes (Newton et al., 2016). In fact, modeling studies using a plethora of experimentally determined parameters (Sassone-Corsi, 2012) have concluded that transient responses exhibit pronounced maxima, reached within 15–30 s of EGF stimulation and followed by a decline to relatively low (quasi-steady-state) levels in all parameters tested at 60, 90, and 120 min. Findings emphasized that 60 min is the earliest time that is reflective of steady state, whether it is immediate postreceptor events or gene change response or cellular phenotypes. We next asked whether we could relate this prediction to patient clinical data for survival that reflect a steady state that evolves over years, not just weeks or months. Numerous prior studies (Tagliaferri et al., 1988; Boucher et al., 2001; Dumaz et al., 2006; Fajardo et al., 2014) have shown that low cAMP facilitates several ominous tumor cell traits, and hence is permissive to cancer progression and worse outcomes. Consistently, numerous studies have confirmed that high GIV levels (which our model predicts will lead to a tonic suppressed cAMP state regardless of the degree of stimulus) are associated with aggressive tumor cell traits and poorer clinical outcomes.

We begin by redefining the input and output for our system. The inputs for these analyses are EGFR mRNA, GIV mRNA, and PDE mRNA, which we use as a surrogate measure of copy number of proteins because others have shown that mRNA can indeed predict protein copy numbers per cell (Edfors et al., 2016) and that mRNA abundance positively correlates with protein levels in healthy and cancer tissues (Kosti et al., 2016). The output is patient survival probability over time, an outcome that reflects tumor aggressiveness, which we use as a surrogate measure of low or suppressed cAMP states. The prediction from the signaling model can be recast as the following hypotheses—EGFR gene-expression levels or PDE gene-expression levels, which by themselves are yet to emerge as clinically useful prognosticators of survival, should become important determinants when analyzed within context of GIV. We formulate and test two hypotheses in the following sections.

Concurrent up-regulation of both GIV and EGFR maximally reduces cAMP and carries poor prognosis in colorectal tumors

To determine the impact of cross-talk between EGFR and GIV on clinical outcome, we compared the mRNA expression levels to disease-free survival (DFS) in a data set of 466 patients with colorectal cancers (see Materials and Methods). Patients were stratified into negative (low) and positive (high) subgroups with regard to GIV (CCDC88A) and EGFR gene-expression levels with the use of the StepMiner algorithm, implemented within the Hegemon software (hierarchical exploration of gene-expression microarrays online; Dalerba et al., 2011; Figure 5A). Kaplan–Meier analyses of DFS over time showed that among patients with high EGFR, expression of GIV at high levels carried a significantly poorer prognosis than among those with low GIV (Figure 5B; Supplemental Figure S13). Among patients with low EGFR, expression of GIV at high or low levels did not impact survival (Figure 5C; Supplemental Figure S13). Among patients with low EGFR, expression of GIV at high or low levels did not impact survival (Figure 5C). That the impact of EGFR-GIV interplay on patient survival is significant in a rigorous Kaplan–Meier analysis of a sufficiently large cohort of patients, despite numerous independent variables, indicates that the interplay between EGFR and the G protein modulator GIV is an important determinant of cancer progression. More importantly, patients with tumors expressing high EGFR did as well as those expressing low EGFR provided the levels of GIV in those tumors were low. These findings reveal that 1) high levels of EGFR signaling do not, by themselves, fuel aggressive traits or carry a poor prognosis, but do so when GIV levels are concurrently elevated; 2) in tumors with low GIV, the high EGFR signaling state may be critical for maintaining high cAMP levels and therefore for dampening several aggressive tumor traits.

Concurrent down-regulation of both GIV and PDE activity maximally increases cAMP and carries a good prognosis in colorectal tumors

Next, to determine the impact of cross-talk between various PDE isoforms and GIV on clinical outcome, we used the StepMiner algorithm, implemented within the Hegemon software on the same set of 466 patients with colorectal cancers as before, except that patients were now stratified into low and high subgroups with regard to GIV (CCDC88A) and PDE gene expression (Figure 5, D–F). Among the 11 known PDE isoforms, we evaluated those that have previously been linked to colon cancer progression (PDE5A, Figures 5, D–F, 4A, and 10A and Supplemental Figure S13). Kaplan–Meier analyses of DFS over time showed that although expression of GIV at high levels was associated with disease progression and poorer survival in both low- and high-PDE groups, the risk of progression was not statistically significant in the high-PDE state (Figure 5E) but was highly significant in the low-PDE state (Figure 5F). Thus, the low GIV/low PDE signature carried a better prognosis than for all other patients. Consistent with the fact that cAMP is a potent anti-tumor second messenger, these findings reveal that 1) high levels of PDE signaling may not be a bad thing, especially when GIV levels are low; 2) in tumors with low PDE signaling, the low GIV signaling state may serve as a key synergy for driving up cAMP levels and therefore be critical for dampening several aggressive tumor traits.

DISCUSSION

Systems biology aims to understand and control the properties of biological networks; experimental data collected using top-down approaches are used to construct in silico bottom-up models, with the ultimate goal of generating experimentally testable predictions. In this work, we used a systems biology approach to construct the first-ever compartmental-network model of growth factor–triggered cAMP signaling and identified two key features of noncanonical G-protein signaling via GIV-GEM.

First, we identified that compartmentalized RTK signaling at the PM and on the endosomes directly imparts a delayed and prolonged cAMP dynamics lasting more than 1 h, which is distinct from the canonical GPCR/G protein pathway; GPCRs initiate more rapid and finite cAMP dynamics on the order of ms to min (Supplemental Figure S5; Violin et al., 2008). In the case of GPCRs, the PM-based signals are believed to be the dominant component of the overall cAMP dynamics with signal attenuation during endocytosis (Supplemental Figure S5; Irannejad et al., 2013). In contrast, in the case of RTK-mediated cAMP dynamics via GIV-GEM, the postendocytic (i.e., endosomal signaling) component constitutes a dominant component of the overall cAMP dynamics that is triggered by RTKs (Supplemental Figure S5). What may be the impact of these distinct temporal features on RTK signaling? It is noteworthy that RTK-triggered cAMP dynamics that is modulated by GIV-GEM spans 5 to >60 min, which coincides with other RTK signaling, trafficking events, and transcriptional response, that is, the major temporal domain of RTK activity, the so-called “window of activity” (Amit et al., 2007b). The 5 min-to-1 h time scale encompasses the time of peak mRNA expression of many immediate-early genes (which peak at 20 min) and delayed-early genes (which peak between 40 min and 2 h); these transcriptional targets not only generate feedback within the RTK-signaling cascade, but also set up cross-talk with other signaling pathways (Amit et al., 2007a,b). In fact, GIV-GEM has been found to modulate myriad downstream signaling pathways from the activity of small GTPases, kinases, and phosphatases to transcription factors (reviewed in Aznar et al., 2016); how GIV-GEM has such a widespread and broad impact had remained a mystery. It is possible that this broad impact could stem from GIV’s ability to modulate the cellular levels of the versatile second messenger cAMP in a sustained manner throughout the window of RTK activity, although other mechanisms might also be in play.

Second, our model identifies that GIV-GEM acts a tunable valve for cAMP by operating at the knot of a bow-tie architecture. Because layering of control of (information) flow is believed to conform to an hourglass architecture (Doyle and Csete, 2011), in which diverse functions and diverse components are intertwined via universal carriers, GIV’s ability to control the universal carrier, cAMP, could explain why GIV has been found to be important for diverse cellular functions and impact diverse components (Aznar et al., 2016).

Third, our work also provides valuable clues into the impact of increased robustness in high-GIV states in cancers. Robustness in signaling is an organizing principle in biology, not only for the maintenance of homeostasis but also in the development and progression of chronic debilitating diseases such as cancers; it is widely accepted that tumor cells hijack such robustness to gain growth and survival advantages during the development of cancer (Kitano, 2004; Amit et al., 2007b; Iadevaia et al., 2014). Consistently, we found that GIV mRNA levels and DNA copy numbers are higher across multiple cancers than in their respective normal tissue of origin (Supplemental Figures S14 and S15). Because GIV has been found to regulate several harmful properties of tumor cells across a variety of cancers (multiple studies, reviewed in Ghosh, 2015), it is possible that the high GIV–driven robustness maintains cAMP at low constant levels despite increasing input signals as a tumor evolves when targeted by biologicals or chemotherapy agents. This phenomenon could be part of a higher-order organizing principle in most aggressive cancers, and therefore justify GIV as a potential target for network-based anti-cancer therapy.

Furthermore, the cross-talk between EGFR and GIV that we define here and its impact on clinical outcomes provide a plausible explanation for some longstanding conundrums in the field of oncology. Deregulated growth factor signaling (e.g., copy number variations or activating mutations in EGFR, increased growth factor production/concentration) is often encountered and targeted for therapy in advanced cancers (Lowery and Yu, 2012). Although activating EGFR mutations, copy number variations, and levels of EGFR protein expression seem to be closely related to each other (Liang et al., 2010), the prognostic impact of EGFR expression in cancers has been ambiguous (Nicholson et al., 2001). In some cancers, high EGFR copy numbers are associated with poor outcomes (Selvaggi et al., 2004; Park et al., 2014); in others, high EGFR expression unexpectedly favors better overall and progression-free survival (Jiang et al., 2013; Llovet et al., 2015; Park et al., 2016). Thus, for reasons that are unclear, not all tumors with high EGF/EGFR signaling have an aggressive clinical course. Dysregulated GIV expression, on the other hand, is consistently associated with poorer outcome across a variety of cancers (Ghosh, 2016). Our finding that GIV levels in tumors with high EGFR are a key determinant of the levels of the anti-tumor second messenger cAMP has provided a potential molecular basis for why elevated EGFR signaling can be beneficial in some tumors, but a driver of metastatic progression in others. Because cAMP levels in tumor cells and GIV levels have been previously implicated in anti-apoptotic signaling (McEwan et al., 2007) and the development of chemoresistance (Waugh, 2012), it is possible that the GIV-EGFR cross-talk we define here also determines how well patients may respond to anti-EGFR therapies, and who may be at highest risk for developing drug resistance. Whether this is the case remains to be evaluated. Similarly, in the context of PDE, it has been demonstrated that overexpression of PDE isoforms in various cancers leads to impaired cAMP and/or cGMP generation (Bender and Beavo, 2006). PDE inhibitors in tumor models in vitro and in vivo have been shown to induce apoptosis and cell cycle arrest in a broad spectrum of tumor cells (Savai et al., 2010). Despite the vast amount of preclinical evidence, there have been no PDE inhibitors that have successfully translated to cancer clinics. For example, based on the role of cAMP in apoptosis and drug resistance, our model predicts that those with low GIV/high EGFR (high cAMP state) are likely to respond well to anti-EGFR therapy inducing tumor cell apoptosis, whereas those with high GIV/high EGFR (low cAMP state) may be at highest risk for developing drug resistance. Similarly, our finding that low PDE levels in the setting of high GIV carry a poor prognosis predicts that the benefits of PDE inhibitors may be limited to patients who have low GIV expression in their tumors. Whether such predictions hold true remains to be investigated.

Model limitations

Although our model captures experimentally observed time courses and generates testable hypotheses, it has a few major limitations. From a model development standpoint, the compartmental well-mixed model we used does not account for the spatial location and geometries of the different compartments and cell shape, many of which can affect the dynamics of cell signaling (Rangamani and Iyengar, 2007; Rangamani et al., 2013). Furthermore, a major concern is the estimation of kinetic parameters for the different reactions. Of the 76 kinetic parameters in our model, a large majority (57) were from models published before. Nineteen parameters, all of which are related to GIV interactions with internalization, Gα_s, and Gα_i, were estimated from experimental data. The uncertainty in some of these parameters was quite large (Supplemental Figure S6). While sloppy parameter space is a problem common to many signaling networks (Gutenkunst et al., 2007), in this case, it is exacerbated by the fact that the model we have constructed is the first of its kind for this pathway. The issue with kinetic parameters also reflects the fact that the field of RTK-G protein regulation is relatively young and makes the case for more quantitative investigations of GIV-GEM modulated signaling. By themselves, these facts can lower confidence in the exact temporal dynamics predicted by the model. However, our confidence in the model is bolstered by the time scales predicted by the phenomenological model and implications for patient survival data.

Additionally, our model focuses exclusively on cAMP as output signal and does not account for other EGF/EGFR-driven signaling pathways that are known to regulate cellular responses such as the Ras-Raf-MEK-ERK pathway. This pathway is known to modulate cAMP and be modulated by GIV via the ability of the latter to affect adaptor protein recruitment to the cytosolic tail of EGFR (Ghosh et al., 2010). However, the vast parameter space associated with model building indicates that this pathway would require its own study. In addition, the dynamics of ERK1/2 activation on the endosomal structure would have to be explored due to GIV-GEM’s interaction with endosomal maturation through Gα_s. The model expansion would lead to an even larger parameter space, with less available data.

Moreover, our model focuses exclusively on EGFR and does not account for the diverse classes of receptors (multiple RTKs, GPCRs, integrins, etc.) that also use GIV to access and modulate G proteins. Model validations used HeLa cells not only because this is the most common model cell line, which has been used exhaustively by us and others to study both GIV and EGFR biology, but also because it has a modest level of expression for both EGFR and GIV. Despite these restrictions, we can identify some fundamental features of growth factor–triggered cAMP signaling for the first time using systems biology, including the roles of compartmentalization, cross-talk between EGFR and GIV, GIV-dependent robustness within the RTK-cAMP signaling axis, and cross-talk between PDE and GIV in controlling cAMP concentration. Although there exist many different signaling pathways downstream of EGFR, GIV-GEM is the only direct link between EGFR and G proteins to date. Owing to the high congruency between the model and validation this leads to two possibilities: either GIV–GEM interaction operates on a standalone basis, or there exist additional feedback loops where GIV–GEM acts as a leading-order contributor (Supplemental Figure S9).

Conclusions

We conclude that GIV utilizes compartmental segregation to modulate the dynamics of RTK→G protein→cAMP signaling and confers robustness on this dynamics by functioning as a tunable control valve. Future systems efforts will build upon this model to unravel further exciting features of GIV as a critical hub for signaling regulation at the knot of a bowtie (Friedlander et al., 2015) and elucidate the hidden complexity that arises from network architecture in noncanonical G-protein signaling.

MATERIALS AND METHODS

Modular construction of the reaction network

A biochemical network model was constructed to capture the main events in the signal transduction cascade from EGF to cAMP through GIV (Figure 2A). We constructed the compartmental computational model in a modular manner, where each module represented key events within the network. The model was trained using key data sets published over the past decade on GIV-GEM, most notably those that defined the spatiotemporal kinetics of EGFR·GIV, EGFR·GIV Gα_i interactions (Ghosh et al., 2010; Ma et al., 2015; Midde et al., 2015), dynamics of phosphoregulation of GIV-GEM (López-Sánchez et al., 2013; Bhandari et al., 2015), and most importantly, the dual modulation of Gα_i/s by GIV-GEM that is brought about by temporally and spatially separated phosphorylation events (Gupta et al., 2016). Last, but not least, we also used the published role of Gα_s in the feedback regulation of endocytic down-regulation of EGF/EGFR signaling (Beas et al., 2012).

The model contains 76 kinetic parameters. Nineteen parameters were fitted with various confidence values (Supplemental Figure S6). Each kinetic parameter used in this model originated from peer-­reviewed publications of computational models or from experimental measurement. The EGF/EGFR and trimeric GTPase–related kinetic parameters were taken from work done by multiple independent groups, often cross-validated across groups engaged in studying each of these paradigms/pathways. For the GIV-related parameters that originate at the interface between the two pathways (EGFR and G proteins), we have used Gupta et al. (2016).

We note here that while there are many more biochemical components involved in signaling from EGF to cAMP, our choice of components was based on experimentally measured temporal dynamics of GIV-GEF and GIV-GDI functions. The modules are as follows.

Module 1 consists of EGFR activation through EGF and internalization dynamics leading to degradation due to the Gα_s·GIV-GDI complex. This module includes the phenomenon that endosomal maturation and EGFR degradation in lysosomes requires the presence of inactive Gα_s [GDP-bound state] (Beas et al., 2012). The presence of Gα_s in the inactive state promotes maturation of endosomes, shuts down the mitogenic MAPK-ERK1/2 signals from endosomes, and suppresses cell proliferation (Beas et al., 2012). In the absence of Gα_s or in cells expressing a constitutively active mutant Gα_s, EGFR stays longer in endosomes, MAPK-ERK1/2 signals are enhanced, and cells proliferate (Beas et al., 2012; Supplemental Figure S5).

Module 2 contains EGFR-mediated activation of p35 and downstream activation of Gα_i though GIV-GEM; the former activated CDK5, which phosphorylates GIV at S1764, allowing the ability to activate Gα_i (López-Sánchez et al., 2013). This phosphoevent does not impact GIV’s ability to inhibit Gα_s (Gupta et al., 2016; Supplemental Figure S2B).

Module 3 contains EGFR-mediated activation of PLC-γ and downstream activation of PKC-θ; the latter phosphorylates GIV at S1689 and terminates its ability to activate Gα_i (López-Sánchez et al., 2013). This phosphoevent does not impact GIV’s ability to inhibit Gα_s (Gupta et al., 2016). Consequently, when it comes to G-protein modulatory functions of GIV, phosphorylation by PKC-θ converts GIV-GEF into GIV-GDI (Supplemental Figure S2D).

Module 4 contains the dynamics of the AC and how it synthesizes cAMP, leading to downstream effectors and controllers (Supplemental Figure S2C). We currently only consider PDE feedback for cAMP reduction. Overall, the model contains 56 reactions. The complete set of reactions for each of the modules, their para­meters and interactions, and the list of assumptions underlying network construction are provided as online supplemental materials (Supplemental Tables S4–S8; Supplemental Figure S8).

We assumed that the signaling components were present in large enough quantities, and different concentrations of each component were computed to explore how varying expression levels in different tissues/cell types impact the signaling pathway. These assumptions allowed us to generate a deterministic dynamical model. The model contains five different compartments: 1) PM, 2) extracellular space, 3) cytosol, 4) endosomes, and 5) endosomal membranes. It was assumed that each compartment is well mixed, and fluxes were used to depict transport across the different compartments so that the dynamic changes in the concentrations of the different components could be tracked. Each interaction was modeled as a chemical reaction using mass-action kinetics for binding–unbinding reactions and Michaelis–Menten kinetics for enzyme-catalyzed reactions, as is standard for models such as this (Bhalla and Iyengar, 1999; Rangamani and Iyengar, 2008).

The network of interactions was constructed using the Virtual Cell modeling platform and was later transferred to COPASI (version 4.24, build 197; available at www.nrcam.uchc.edu, http://copasi.org/). We chose this platform because it is a user-friendly computational cell biology program, which allows us to generate the system of differential equations based on the input reactions and has been used successfully to model signaling networks of various sizes with a high degree of numerical accuracy (Loew and Schaff, 2001; Slepchenko et al., 2003; Ditlev et al., 2012; Falkenberg et al., 2013). The model was later exported into COPASI to leverage the built-in fitting techniques. Also, the Virtual Cell and COPASI platform has built-in capabilities to conduct dynamic sensitivity analysis, which is an important aspect of dynamic systems modeling. As we discuss in later sections, we use this capability to identify sources of system robustness and sloppiness.

Characteristics of the signaling cascade

To characterize the dynamics of the different protein activities, we use the area under the concentration vs. time curve (Atay and Skotheim, 2017). The area under the curve gives the total signal activated over the time of observation and for the ith species:

This gives a measure of the total signal for different conditions.

Comparison with experimental data

Raw data corresponding to Figure 1, C and D, of Gupta et al. (2016) were used for model fitting. The data were normalized so that the initial value was 1. Parameter fitting using COPASI (Mendes et al., 2009) was used to match the normalized experimental data against the model output, with corresponding expected initial values based on the experimental method. Goodness of fit between experimental values and model output was determined using root-mean-squared error (RMSE).

Choice of assays for model validation

We chose to validate our temporal–spatial model for the dynamic assembly of the EGFR·GIV Gα_i and Gα_s·GIV-GDI complexes using previously published protein–protein interaction assays carried out in cells responding to EGF (Gupta et al., 2016). For EGFR·GIV Gα_i complexes that are formed within 5 min after ligand stimulation at the plasma membrane, we modeled the GST pull-down assays carried out using GST-GIV-CT that is expressed in cells (and hence phosphomodified in response to EGF stimulation) and endogenous Gα_i. The findings of these pull-down assays mirrored observations by FRET-based (Midde et al., 2015) and co-IP assays (Ghosh et al., 2010; Lin et al., 2014). For Gα_s·GIV-GDI complexes that are formed later and on endosomes, we modeled the proximity-ligation assays (PLA) on endogenous GIV and Gα_s proteins, which provide a crude estimate of complexes on endomembranes.

Dynamic parametric sensitivity analysis

Because a continuing challenge in building computational models of signaling networks is the choice of kinetic parameters, we conducted a dynamic parametric sensitivity analysis. This sensitivity analysis of the model was performed with the goal of identifying the set of parameters and initial concentrations that the model response is most sensitive to. The log sensitivity coefficient of the concentration of the ith species, C_i, with respect to parameter k_j is given (Ingalls and Sauro, 2003; Varma et al., 2005) by

Because we are studying a dynamical system and not steady state behavior, we used COPASI to calculate the local log sensitivity at the 5-, 15-, 30-, and 60-min points (Supplemental Figure S4). The resulting values give information about the time dependence of parametric sensitivity coefficients for the system at those points. The variable of interest, C_i, is said to be robust with respect to a parameter k_j if the log sensitivity is of order 1 (Varma et al., 2005). We refer the reader to (Ingalls and Sauro, 2003; Varma et al., 2005) for a complete introduction to dynamic sensitivity analysis. We conducted dynamic sensitivity analyses for all the kinetic parameters, initial concentrations of the different species, and compartment sizes in the model (Supplemental Figure S4). The variation (delta factor) used was 0.001 with a delta minimum of 1 × 10^–1; for a value X, S∈ [0.999X, 1.001X]. Sensitive parameters or corresponding outputs of interest (cAMP,vGα_s·GIV-GDI, EGFR·GIVGα_i) are reported in Supplemental Tables S12–S18.

Measurement of cAMP

HeLa cells were serum-starved (0.2% fetal bovine serum, 16 h) and incubated with isobutylmethylxanthine (IBMX, 200 μM, 20 min) followed by EGF. Stimulation was carried out using either fixed EGF concentrations followed by assessment of cAMP at various time points (as in Figure 3B) or various EGF concentrations followed by an assessment of cAMP at 60 min (as in Figure 3, F–J). Reactions were terminated by aspiration of media and addition of 150 μl of ice-cold TCA 7.5% (wt/vol). cAMP content in TCA extracts was determined by RIA and normalized to protein (determined using a dye-binding protein assay [Bio-Rad]; Lopez-Sanchez et al., 2014; Aznar et al., 2016). Data are expressed as fmol cAMP/μg total protein.

Stratification of colon cancer patients in distinct gene-expression subgroups and comparative analysis of their survival outcomes

The association between the levels of GIV (CCDC88A) and either EGFR or PDE mRNA expression and patient survival was tested in a cohort of 466 patients where each tumor had been annotated with the DFS information of the corresponding patient. This cohort included gene expression data from four publicly available National Center for Biotechnology Information (NCBI)-Gene Expression Omnibus (GEO) data series (GSE14333, GSE17538, GSE31595, GSE37892; Jorissen et al., 2009; Smith et al., 2010; Laibe et al., 2012; Thorsteinsson et al., 2012) and contained information on 466 unique primary colon carcinoma samples, collected from patients at various clinical stages (AJCC Stage I–IV/Duke’s Stage A–D) by five independent institutions: 1) the H. Lee Moffit Cancer Center in Tampa, FL (n = 164); 2) the Vanderbilt Medical Center in Nashville, TN (n = 55); 3) the Royal Melbourne Hospital in Melbourne, Australia (n = 80); 4) the Institut PaoliCalmette in Marseille, France (n = 130); 5) the Roskilde Hospital in Copenhagen, Denmark (n = 37). To avoid redundancies (i.e., identical samples replicated two or more times across multiple NCBI-GEO data sets) all 466 samples contained in this subset were cross-checked to exclude the presence of duplicates. A complete list of all GSMIDs of the experiments contained within the NCBI-GEO discovery data set has been published previously (Dalerba et al., 2011). To investigate the relationship between the mRNA expression levels of selected genes (i.e., CCDCDDC, Wnt5a, EGFR, and FZD7) and the clinical outcomes of the 466 colon cancer patients represented within the NCBI-GEO discovery data set, we applied the Hegemon software tool (Dalerba et al., 2011). The Hegemon software is an upgrade of the BooleanNet software (Sahoo et al., 2008), where individual gene-expression arrays, after having been plotted on a two-axis chart based on the expression levels of any two given genes, can be stratified using the StepMiner algorithm and automatically compared for survival outcomes using Kaplan–Meier curves and log-rank tests. Because all 466 samples contained in the data set had been analyzed using the Affymetrix HG-U133 Plus 2.0 platform (GPL570), the threshold gene-expression levels for GIV/CCDC88A, PDE, and EGFR were calculated using the StepMiner algorithm based on the expression distribution of the 25,955 experiments performed on the Affymetrix HG-U133 Plus 2.0 platform. We stratified the patient population of the NCBI-GEO discovery data set into different gene-expression subgroups, based on either the mRNA expression levels of GIV/CCDC88A alone (i.e., CCDC88A neg vs. pos), PDE alone (i.e., PDE neg vs. pos), or EGFR alone (i.e., EGFR neg vs. pos) or a combination of GIV and either EGFR or PDE. Once grouped based on their gene-expression levels, patient subsets were compared for survival outcomes using both Kaplan–Meier survival curves and multivariate analysis based on the Cox proportional-hazards method.

FOOTNOTES

ACKNOWLEDGMENTS

REFERENCES

• Amit I, Citri A, Shay T, Lu Y, Katz M, Zhang F, Tarcic G, Siwak D, Lahad J, Jacob-Hirsch J, et al. (2007a). A module of negative feedback regulators defines growth factor signaling. Nat Genet 39, 503–512. Crossref, Medline, Google Scholar
• Amit I, Wides R, Yarden Y (2007b). Evolvable signaling networks of receptor tyrosine kinases: relevance of robustness to malignancy and to cancer therapy. Mol Syst Biol 3, 151. Crossref, Medline, Google Scholar
• Atay O, Skotheim JM (2017). Spatial and temporal signal processing and decision making by MAPK pathways. J Cell Biol 216, 317–330. Crossref, Medline, Google Scholar
• Aznar N, Kalogriopoulos N, Midde KK, Ghosh P (2016). Heterotrimeric G protein signaling via GIV/girdin: breaking the rules of engagement, space, and time. Bioessays 38, 379–393. Crossref, Medline, Google Scholar
• Beas AO, Taupin V, Teodorof C, Nguyen LT, Garcia-Marcos M, Farquhar MG (2012). Gas promotes EEA1 endosome maturation and shuts down proliferative signaling through interaction with GIV (Girdin). Mol Biol Cell 23, 4623–4634. Link, Google Scholar
• Beavo JA, Brunton LL (2002). Cyclic nucleotide research—still expanding after half a century. Nat Rev Mol Cell Bio 3, 710–718. Crossref, Medline, Google Scholar
• Bender AT, Beavo JA (2006). Cyclic nucleotide phosphodiesterases: molecular regulation to clinical use. Pharmacol Rev 58, 488–520. Crossref, Medline, Google Scholar
• Berkers JA, van Bergen en Henegouwen PM, Boonstra J (1991). Three classes of epidermal growth factor receptors on HeLa cells. J Biol Chem 266, 922–927. Medline, Google Scholar
• Bhalla US, Iyengar R (1999). Emergent properties of networks of biological signaling pathways. Science 283, 381–387. Crossref, Medline, Google Scholar
• Bhandari D, Lopez-Sanchez I, To A, Lo IC, Aznar N, Leyme A, Gupta V, Niesman I, Maddox AL, Garcia-Marcos M, et al. (2015). Cyclin-dependent kinase 5 activates guanine nucleotide exchange factor GIV/girdin to orchestrate migrationproliferation dichotomy. Proc Natl Acad Sci USA 112, E4874–E4883. Crossref, Medline, Google Scholar
• Boucher M-J, Duchesne C, Lainé J, Morisset J, Rivard N (2001). CAMP protection of pancreatic cancer cells against apoptosis induced by erk inhibition. Biochem Bioph Res Co 285, 207–216. Crossref, Medline, Google Scholar
• Bruce JI, Straub SV, Yule DI (2003). Crosstalk between cAMP and Ca2+ signaling in non-excitable cells. Cell Calcium 34, 431–444. Crossref, Medline, Google Scholar
• Dalerba P, Kalisky T, Sahoo D, Rajendran PS, Rothenberg ME, Leyrat AA, Sim S, Okamoto J, Johnston DM, Qian D, et al. (2011). Single-cell dissection of transcriptional heterogeneity in human colon tumors. Nat Biotechnol 29, 1120–1127. Crossref, Medline, Google Scholar
• Ditlev JA, Michalski PJ, Huber G, Rivera GM, Mohler WA, Loew LM, Mayer BJ (2012). Stoichiometry of Nck-dependent actin polymerization in living cells. J Cell Biol 197, 643–658. Crossref, Medline, Google Scholar
• Doyle JC, Csete M (2011). Architecture, constraints, and behavior. Proc Natl Acad Sci USA 108, 15624–15630. Crossref, Medline, Google Scholar
• Dumaz N, Hayward R, Martin J, Ogilvie L, Hedley D, Curtin JA, Bastian BC, Springer C, Marais R (2006). In melanoma, RAS mutations are accompanied by switching signaling from BRAF to CRAF and disrupted cyclic amp signaling. Cancer Res 66, 9483–9491. Crossref, Medline, Google Scholar
• Edfors F, Danielsson F, Hallström BM, Käll L, Lundberg E, Pontén F, Forsström B, Uhlén M (2016). Gene-specific correlation of RNA and protein levels in human cells and tissues. Mol Syst Biol 12, 883. Crossref, Medline, Google Scholar
• Fajardo AM, Piazza GA, Tinsley HN (2014). The role of cyclic nucleotide signaling pathways in cancer: targets for prevention and treatment. Cancers (Basel) 6, 436–458. Crossref, Medline, Google Scholar
• Falkenberg CV, Blinov ML, Loew LM (2013). Pleomorphic ensembles: formation of large clusters composed of weakly interacting multivalent molecules. Biophys J 105, 2451–2460. Crossref, Medline, Google Scholar
• Filmore D (2004). It’s a GPCR world. Modern Drug Discovery 7, 24–28. Google Scholar
• French AR, Tadaki DK, Niyogi SK, Lauffenburger DA (1995). Intracellular trafficking of epidermal growth factor family ligands is directly influenced by the pH sensitivity of the receptor/ligand interaction. J Biol Chem 270, 4334–4340. Crossref, Medline, Google Scholar
• Friedlander T, Mayo AE, Tlusty T, Alon U (2015). Evolution of bow-tie architectures in biology. PLoS Comput Biol 11, e1004055. Crossref, Medline, Google Scholar
• Garcia-Marcos M, Ghosh P, Farquhar MG (2009). GIV is a nonreceptor gef for gαi with a unique motif that regulates akt signaling. Proc Natl Acad Sci USA 106, 3178–3183. Crossref, Medline, Google Scholar
• Ghosh P (2015). Heterotrimeric G proteins as emerging targets for network based therapy in cancer: end of a long futile campaign striking heads of a Hydra. Aging (Albany NY) 7, 469–474. Crossref, Medline, Google Scholar
• Ghosh P (2016). The untapped potential of tyrosine-based G protein signaling. Pharmacol Res 105, 99–107. Crossref, Medline, Google Scholar
• Ghosh P, Beas AO, Bornheimer SJ, Garcia-Marcos M, Forry EP, Johannson C, Ear J, Jung BH, Cabrera B, Carethers JM, et al. (2010). A Galphai-GIV molecular complex binds epidermal growth factor receptor and determines whether cells migrate or proliferate. Mol Biol Cell 21, 2338–2354. Link, Google Scholar
• Ghosh P, Rangamani P, Kufareva I (2017). The GAPs, GEFs, GDIs and … now, GEMs: new kids on the heterotrimeric G protein signaling block. Cell Cycle 0, 1–6. Google Scholar
• Goraya TA, Cooper DM (2005). Ca2+-calmodulin-dependent phosphodiesterase (PDE1): current perspectives. Cell Signal 17, 789–797. Crossref, Medline, Google Scholar
• Gupta V, Bhandari D, Leyme A, Aznar N, Midde KK, Lo IC, Ear J, Niesman I, López-Sánchez I, Blanco-Canosa JB, et al. (2016). GIV/girdin activates Gα_i and inhibits Gα_s via the same motif. Proc Natl Acad Sci USA 113, E5721–E5730. Crossref, Medline, Google Scholar
• Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP (2007). Universally sloppy parameter sensitivities in systems biology models. PLOS Computational Biology 3, e189. Crossref, Google Scholar
• Iadevaia S, Nakhleh LK, Azencott R, Ram PT (2014). Mapping network motif tunability and robustness in the design of synthetic signaling circuits. PLoS One 9, e91743. Crossref, Medline, Google Scholar
• Ingalls BP, Sauro HM (2003). Sensitivity analysis of stoichiometric networks: an extension of metabolic control analysis to non-steady state trajectories. J Theor Biol 222, 23–36. Crossref, Medline, Google Scholar
• Irannejad R, Tomshine JC, Tomshine JR, Chevalier M, Mahoney JP, Steyaert J, Rasmussen SGF, Sunahara RK, El-Samad H, Huang B, et al. (2013). Conformational biosensors reveal GPCR signalling from endosomes. Nature 495, 534–538. Crossref, Medline, Google Scholar
• Jiang Z, Li C, Li F, Wang X (2013). EGFR gene copy number as a prognostic marker in colorectal cancer patients treated with cetuximab or panitumumab: a systematic review and meta analysis. PLoS One 8, e56205. Crossref, Medline, Google Scholar
• Jorissen RN, Gibbs P, Christie M, Prakash S, Lipton L, Desai J, Kerr D, Aaltonen LA, Arango D, Kruhøffer M, et al. (2009). Metastasis-associated gene expression changes predict poor outcomes in patients with Dukes stage B and C colorectal cancer. Clin Cancer Res 15, 7642–7651. Crossref, Medline, Google Scholar
• Kitano H (2004). Cancer as a robust system: implications for anticancer therapy. Nat Rev Cancer 4, 227–235. Crossref, Medline, Google Scholar
• Kosti I, Jain N, Aran D, Butte AJ, Sirota M (2016). Cross-tissue analysis of gene and protein expression in normal and cancer tissues. Sci Rep 6, 24799. Crossref, Medline, Google Scholar
• Laibe S, Lagarde A, Ferrari A, Monges G, Birnbaum D, Olschwang S (2012). A seven-gene signature aggregates a subgroup of stage II colon cancers with stage III. OMICS 16, 560–565. Crossref, Medline, Google Scholar
• Liang Z, Zhang J, Zeng X, Gao J, Wu S, Liu T (2010). Relationship between EGFR expression, copy number and mutation in lung adenocarcinomas. BMC Cancer 10, 376. Crossref, Medline, Google Scholar
• Lin C, Ear J, Midde K, Lopez-Sanchez I, Aznar N, Garcia-Marcos M, Kufareva I, Abagyan R, Ghosh P (2014). Structural basis for activation of trimeric Gi proteins by multiple growth factor receptors via GIV/girdin. Mol Biol Cell 25, 3654–3671. Link, Google Scholar
• Llovet P, Sastre J, Ortega JS, Bando I, Ferrer M, García-Alfonso P, Donnay O, Carrato A, Jiménez A, Aranda E, et al. (2015). Prognostic value of BRAF, PI3K, PTEN, EGFR copy number, amphiregulin and epiregulin status in patients with KRAS codon 12 wild-type metastatic colorectal cancer receiving first-line chemotherapy with anti-EGFR therapy. Mol Diagn Ther 19, 397–408. Crossref, Medline, Google Scholar
• Loew LM, Schaff JC (2001). The Virtual Cell: a software environment for computational cell biology. Trends Biotechnol 19, 401–406. Crossref, Medline, Google Scholar
• Lopez-Sanchez I, Dunkel Y, Roh YS, Mittal Y, Minicis SD, Muranyi A, Singh S, Shanmugam K, Aroonsakool N, Murray F, et al. (2014). GIV/girdin is a central hub for profibrogenic signalling networks during liver fibrosis. Nat Commun 5, 4451. Crossref, Medline, Google Scholar
• López-Sánchez I, Garcia-Marcos M, Mittal Y, Aznar N, Farquhar MG, Ghosh P (2013). Protein kinase C-theta (PKCθ) phosphorylates and inhibits the guanine exchange factor, GIV/Girdin. Proc Natl Acad Sci USA 110, 5510–5515. Crossref, Medline, Google Scholar
• Lowery FJ, Yu D (2012). Growth factor signaling in metastasis: current understanding and future opportunities. Cancer Metast Rev 31, 479–491. Crossref, Medline, Google Scholar
• Ma GS, Aznar N, Kalogriopoulos N, Midde KK, Lopez-Sanchez I, Sato E, Dunkel Y, Gallo RL, Ghosh P (2015). Therapeutic effects of cell-permeant peptides that activate G proteins downstream of growth factors. Proc Natl Acad Sci USA 112, E2602–E2610. Crossref, Medline, Google Scholar
• MacKenzie SJ, Baillie GS, McPhee I, MacKenzie C, Seamons R, McSorley T, Millen J, Beard MB, van Heeke G, Houslay MD (2002). Long PDE4 cAMP specific phosphodiesterases are activated by protein kinase A-mediated phosphorylation of a single serine residue in Upstream Conserved Region 1 (UCR1). Br J Pharmacol 136, 421–433. Crossref, Medline, Google Scholar
• Marty C, Richard DY (2010). Heterotrimeric g protein signaling outside the realm of seven transmembrane domain receptors. Mol Pharmacol 78, 12–18. Crossref, Medline, Google Scholar
• McEwan DG, Brunton VG, Baillie GS, Leslie NR, Houslay MD, Frame MC (2007). Chemoresistant KM12C colon cancer cells are addicted to low cyclic AMP levels in a phosphodiesterase 4-regulated compartment via effects on phosphoinositide 3-kinase. Cancer Res 67, 5248–5257. Crossref, Medline, Google Scholar
• Mendes P, Hoops S, Sahle S, Gauges R, Dada J, Kummer U (2009). Computational modeling of biochemical networks using COPASI. Methods Mol Biol 500, 17–59. Crossref, Medline, Google Scholar
• Midde KK, Aznar N, Laederich MB, Ma GS, Kunkel MT, Newton AC, Ghosh P (2015). Multimodular biosensors reveal a novel platform for activation of G proteins by growth factor receptors. Proc Natl Acad Sci USA 112, E937–E946. Crossref, Medline, Google Scholar
• Murthy KS, Zhou H, Makhlouf GM (2002). PKA-dependent activation of PDE3A and PDE4 and inhibition of adenylyl cyclase V/VI in smooth muscle. Am J Physiol-Cell Ph 282, C508–C517. Crossref, Medline, Google Scholar
• Newton AC, Bootman MD, Scott JD (2016). Second messengers. Cold Spring Harb Perspect Biol 8, a005926. Crossref, Medline, Google Scholar
• Nicholson RI, Gee JM, Harper ME (2001). EGFR and cancer prognosis. Eur J Cancer 37 (Suppl 4), 9–15. Crossref, Medline, Google Scholar
• Park HS, Jang MH, Kim EJ, Kim HJ, Lee HJ, Kim YJ, Kim JH, Kang E, Kim SW, Kim IA, et al. (2014). High EGFR gene copy number predicts poor outcome in triple-negative breast cancer. Mod Pathol 27, 1212–1222. Crossref, Medline, Google Scholar
• Park JS, Kim HS, Bae YS, Cheong JH, Rha SY, Noh SH, Kim H (2016). Prognostic significance and frequency of EGFR expression and amplification in surgically resected advanced gastric cancer. Jpn J Clin Oncol 46, 507–516. Crossref, Medline, Google Scholar
• Rangamani P, Iyengar R (2007). Modelling spatio-temporal interactions within the cell. J Biosci 32, 157–167. Crossref, Medline, Google Scholar
• Rangamani P, Iyengar R (2008). Modelling cellular signalling systems. Essays Biochem 45, 83–94. Crossref, Medline, Google Scholar
• Rangamani P, Lipshtat A, Azeloglu EU, Calizo RC, Hu M, Ghassemi S, Hone J, Scarlata S, Neves SR, Iyengar R (2013). Decoding information in cell shape. Cell 154, 1356–1369. Crossref, Medline, Google Scholar
• Sahoo D, Dill DL, Gentles AJ, Tibshirani R, Plevritis SK (2008). Boolean implication networks derived from large scale, whole genome microarray datasets. Genome Biol 9, R157. Crossref, Medline, Google Scholar
• Sassone-Corsi P (2012). The cyclic AMP pathway. Cold Spring Harb Perspect Biol 4, a011148. Crossref, Medline, Google Scholar
• Savai R, Pullamsetti SS, Banat GA, Weissmann N, Ghofrani HA, Grimminger F, Schermuly RT (2010). Targeting cancer with phosphodiesterase inhibitors. Expert Opin Investig Drugs 19, 117–131. Crossref, Medline, Google Scholar
• Schoeberl B, Eichler-Jonsson C, Gilles ED, Müller G (2002). Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol 20, 370–375. Crossref, Medline, Google Scholar
• Selvaggi G, Novello S, Torri V, Leonardo E, De Giuli P, Borasio P, Mossetti C, Ardissone F, Lausi P, Scagliotti GV (2004). Epidermal growth factor receptor overexpression correlates with a poor prognosis in completely resected non-small-cell lung cancer. Ann Oncol 15, 28–32. Crossref, Medline, Google Scholar
• Sette C, Conti M (1996). Phosphorylation and activation of a cAMP-specific phosphodiesterase by the cAMP-dependent protein kinase. Involvement of serine 54 in the enzyme activation. J Biol Chem 271, 16526–16534. Crossref, Medline, Google Scholar
• Slepchenko BM, Schaff JC, Macara I, Loew LM (2003). Quantitative cell biology with the Virtual Cell. Trends Cell Biol 13, 570–576. Crossref, Medline, Google Scholar
• Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, Lu P, Johnson JC, Schmidt C, Bailey CE, et al. (2010). Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. Gastroenterology 138, 958–968. Crossref, Medline, Google Scholar
• Sutherland EW, Rall TW (1958). Fractionation and characterization of a cyclic adenine ribonucleotide formed by tissue particles. J Biol Chem 232, 1077–1092. Medline, Google Scholar
• Tagliaferri P, Katsaros D, Clair T, Ally S, Tortora G, Neckers L, Rubalcava B, Parandoosh Z, Chang Y, Revankar GR, et al. (1988). Synergistic inhibition of growth of breast and colon human cancer cell lines by site-selective cyclic amp analogues. Cancer Res 48, 1642–1650. Medline, Google Scholar
• Thorsteinsson M, Kirkeby LT, Hansen R, Lund LR, Sørensen, LT, Gerds TA, Jess P, Olsen J (2012). Gene expression profiles in stages II and III colon cancers: application of a 128-gene signature. Int J Colorectal Dis 27, 1579–1586. Crossref, Medline, Google Scholar
• Tremblay J, Lachance B, Hamet P (1985). Activation of cyclic GMP-binding and cyclic AMP-specific phosphodiesterases of rat platelets by a mechanism involving cyclic AMP-dependent phosphorylation. J Cyclic Nucl Prot 10, 397–411. Medline, Google Scholar
• Varma A, Morbidelli M, Wu H (2005). Parametric Sensitivity in Chemical Systems, New York: Cambridge University Press. Google Scholar
• Vilardaga J-P, Jean-Alphonse FG, Gardella TJ (2014). Endosomal generation of cAMP in GPCR signaling. Nat Chem Biol 10, 700–706. Crossref, Medline, Google Scholar
• Violin JD, DiPilato LM, Yildirim N, Elston TC, Zhang J, Lefkowitz RJ (2008). β2-adrenergic receptor signaling and desensitization elucidated by quantitative modeling of real time cAMP dynamics. J Biol Chem 283, 2949–2961. Crossref, Medline, Google Scholar
• Waugh MG (2012). Phosphatidylinositol 4-kinases, phosphatidylinositol 4-phosphate and cancer. Cancer Lett 325, 125–131. Crossref, Medline, Google Scholar
• Yoshimasa T, Sibley DR, Bouvier M, Lefkowitz RJ, Caron MG (1987). Cross-talk between cellular signalling pathways suggested by phorbol-ester-induced adenylate cyclase phosphorylation. Nature 327, 67–70. Crossref, Medline, Google Scholar