Abstract
Cellular levels of the versatile second messenger, cyclic-(c)AMP are regulated by the antagonistic actions of the canonical G protein→adenylyl cyclase pathway that is initiated by G-protein-coupled receptors (GPCRs) and by phosphodiesterases (PDEs); dysregulated cAMP signaling drives many diseases, including cancers. Recently, an alternative paradigm for cAMP signaling has emerged, in which growth factor-receptor tyrosine kinases (RTKs; e.g., EGFR) access and modulate G proteins via cytosolic guanine-nucleotide exchange modulator (GEM), GIV/Girdin; dysregulation of this pathway is frequently encountered in cancers. Here we present a comprehensive network-based compartmental model for the paradigm of GEM-dependent signaling that reveals unforeseen crosstalk and network dynamics between upstream events and the various feedback-loops that fine-tune the GEM action of GIV, and captures the experimentally determined dynamics of cAMP. The model also reveals that GIV acts a tunable control-valve within the RTK→cAMP pathway; hence, it modulates cAMP via mechanisms distinct from the two most-often targeted classes of cAMP modulators, GPCRs and PDEs.
- cAMP
- Cyclic adenosine monophosphate
- RTK
- Receptor Tyrosine Kinase
- EGFR
- Epidermal growth factor receptor
- GEM
- G protein exchange modulator
- GEF
- Guanine nucleotide exchange factor
- GDI
- Guanine nucleotide dissociation inhibitor
- SH2
- Src homology 2
- PM
- plasma membrane
- DAG
- Diacylglycerol
- EGF
- Epidermal growth factor
- AUC
- Area under the curve
- AC
- Adenylyl cyclase
- PDE
- Phosphodiesterase
- AMP
- adenosine monophosptate
- ATP
- adenosine triphosphate
- CDK5
- cyclin dependent kinase 5
- PKC-θ
- Protein kinase C θ
- PLC-γ
- Phospholipase C γ
- Epac1
- Exchange factor directly activated by cAMP 1
- RIA
- Radioimmunoassay
- IBMX
- 3-isobutyl-1-methylxant hene
1 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 [1, 2]. Of the various cyclic nucleotides, the first to be identified was cyclic adenosine 3,5-monophosphate (cAMP), a universal second messenger that is used by diverse forms of life, from the unicellular bacteria, to fungi, protozoans and mammals. cAMP relays signals triggered by hormones, ion channels, and neurotransmitters [3], and also binds and regulates other cAMP-binding proteins such as PKA and Epac1 [4].
The 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; they transmit signals from cell-surface receptors to second messengers. PDEs, on the other hand, are soluble, and catalyze the degradation of the phosphodiester bond resulting in the conversion of cAMP to AMP. PDEs are activated by protein kinase A (PKA), a downstream effector of cAMP, resulting in a negative feedback loop between cAMP and PDEs [5–8]. Thus, the level of cAMP in cells is a fine balance of synthesis by ACs, degradation by PDEs, and feedback through the PKA-PDE loop [3]. Both ACs and PDEs are also subject to positive and negative regulation by numerous other signaling pathways [9–11], which coordinately maintain cAMP levels in normal cells within a finite physiologic range. Dysregulated circuits that give rise to too much or too little cAMP can be unhealthy, and many diseased states in humans are characterized by signaling programs driven by abnormal levels of cellular cAMP [see legend, Figure 1A]. For example, in the context of cancers, multiple studies across different cancers (e.g., breast [12], melanoma [13], pancreas [14], etc.) agree that high levels of cAMP are generally protective, whereas low cAMP levels fuel cancer progression [reviewed in [15]]. High cAMP in-hibits several sinister phenotypes of tumor cells such as proliferation, invasion, stemness and chemoresistance, while enhancing differentiation and apoptosis [see legend, Figure 1A]. Therapies that target the canonical GPCR/G-protein-cAMP signaling pathway have been successfully translated to the clinics, and they account for 40% of currently marketed drugs that can treat a wide range of ailments [16], from hypertension to glaucoma. However, such strategies have largely failed to impact cancer care or outcome. 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 non-canonical G protein signaling that is initiated by growth factors [17–20] 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 also often encountered in advanced tumors and is frequently targeted with varying degrees of success [21]. A body of work published by us and others have revealed that RTKs bind and activate trimeric G proteins via a family of proteins called Guanine-nucleotide Exchange Modulators (GEMs; [17]). While the members of this family share very little sequence homology, and act within diverse signaling cascades, what unites them is the ability to couple activation of these cascades to G protein signaling via an evolutionarily conserved motif of approximately 30 amino acids that directly binds and modulates G proteins. By demonstrating how GIV, a prototypical member of a family of cytosolic guanine nucleotide exchange modulators (GEMs; [17, 22]), uses a SH2-like module [23] to directly bind cytoplasmic tails of ligand-activated RTKs such as EGFR [24] we provided a definitive structural basis for several decades of observations made by researchers that G-proteins can be coupled to and activated by growth factors (reviewed in [25]). A series of studies that have followed since have revealed that growth factor-triggered non-canonical G protein→cAMP signaling through GIV has unique spatiotemporal properties and prolonged dynamics that are distinct from canonical GPCR-dependent signaling [reviewed in [18]]. More importantly, by straddling two major eukaryotic signaling hubs [RTKs and G proteins] that are most frequently targeted for their therapeutic significance, GIV has its own growing list of pathophysiologic importance, as a therapeutic target in a variety of disease states, most prominently in cancers. High levels of GIV expression fuels multiple ominous properties of cancer cells, e.g., invasiveness, chemoresistance, stemness, survival, etc., and is associated with poorer outcome in multiple cancers, and inhibition of GIV’s G protein modulatory function has emerged as a plausible strategy to combat aggressive traits of cancers [26]; reviewed in [17, 20].
Despite these insights, the paradigm of RTK-dependent cAMP signaling remains nascent with many unknowns. Although the pathway may appear to be a linear connection between input (the growth factor RTKs) and output (G-proteins) elements, experimental data shows that non-canonical G protein→cAMP signaling via GIV-GEM integrates multiple input and output signals, with multiple feedback loops, and yet remains spatiotemporally segregated [initiated and terminated by specific phosphoevents, at specific times, on distinct membranes]. However, signaling through this pathway can be disrupted by disassembling any one of the key signaling interfaces (reviewed in [20]). Therefore, the behavior of such complex systems is hard to grasp by intuition.
Here, we use systems biology approaches, namely mathematical and computational modeling as the primary tools to generate a comprehensive model that can provide insight into some of these issues. This model, the first of its kind, not only connects two of 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 outcome.
Results and discussion
Construction and experimental validation of a compartmental model for non-canonical G-protein signaling triggered by growth factors
The emerging paradigm of non-canonical modulation of Gi/s proteins by growth factor RTKs is comprised of several temporally and spatially separated components (Figure 1B, C); each component is analyzed as distinct modules within a larger network model (Figure 1D). The model consists of four modules – Module 1 focuses on the dynamics of EGFR (Table S2); Module 2 represents the dynamics of the formation of the EGFR·GIV·Gαi complex (Table S3); Module 3 represents the dynamics of the formation of the Gαs·GIVGDI complex (Table S4); and Module 4 represents the dynamics of cAMP formation (Table S5,S6).
Within each module, the biochemical reaction network includes several known interactions. For example, binding of EGF to EGFR at the plasma membrane (PM) initiates a cascade of events, which includes receptor dimerization, cross-phosphorylation of the cytoplasmic tails, recruitment of signaling adaptors, and phosphoactivation of a plethora of enzymes to relay downstream signaling. Of relevance to this paradigm, GIV, a multi-domain signal transducer contains a SH2-like domain in its C-terminus, which enables its recruitment to autophosphorylated cytoplasmic tails of EGFR [23]. 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.
EGFR dynamics at the plasma membrane and the endosomal membrane
Module 1 of the reaction network models the dynamics of EGFR at the PM and the endosome (Figure 2A, Table S2). At the PM, EGFR is activated by ligand binding, receptor dimerization, and cross-phosphorylation; activated EGFR is internalized to the early endosome through endocytosis, from where it can be either recycled or degraded [27]. Active PM EGFR also forms a complex with GIV-GEF, and via GIV with Gαi, leading to the activation of Gαi [22]. On the other hand, while it remains unclear when and where EGF/EGFR activates Gαs, it is known that a pool of GIV that is on endosomes containing internalized EGFR binds and inactivates Gαs on the endosomal membrane. Once inactivated, Gαs-GDP enhances the degradation rate of internalized, endosomal EGFR, thereby limiting the pool of receptors available for recycling to the PM and serves to attenuate growth factor signaling [28].
Simulations from the model show that EGFR dynamics is governed by multiple time-scales when ligand stimulation triggers the redistribution of receptors from the PM to different pools (Figure 2B). The PM-pool of active receptors increases rapidly upon ligand stimulation (Figure 2B, red line) and subsequently recruits GIV, forming GIV-GEF·EGFR complexes (Figure 2B, purple line). The endosomal pool of active receptors increases at a slower time scale (Figure 2B, yellow line) than the PM-pool of active receptors. Recycling of the endosomal pool of receptors to the PM leads to a small second burst in the PM pool of receptors around 10 min (Figure 2B, yellow line). These findings are in agreement with Schoeberl et al. [27], indicating that our model accurately captures the EGFR dynamics. The total number of active receptors decreases over time because of receptor degradation (Figure 2B, blue line). The pool of receptors in the GIV-GEF-EGFR complex subsequently interact with Gαi at the PM to form the EGFR-GIV-Gαi complex. The effect of kinetic parameters of EGFR dynamics is shown in Figure S1 and we find that the balance of PM-pool versus internalized pool of EGFR is closely regulated by both the internalization rate and the Gαs-GDP dependent receptor degradation rate [28].
Dynamics of Gαi signaling: activation kinetics are shaped by both upstream EGFR dynamics and downstream PLC-γ→ DAG →PKC-θ signaling events
We next asked how EGFR dynamics affect the dynamics of Gαi signaling at the PM. Activation of EGFR at the PM triggers a series of downstream events, including the activation of CDK5 at the PM by its cofactor, p35 [29]. CDK5 phosphorylates GIV at Ser(S)1675 and enhances GIV’s ability to bind Gαi, i.e., CDK5 turns inactive GIV to into active GIV-GEF [30]. This allows GIV to couple Gαi to EGFR by assembling ternary EGFR-GIV-Gαi complexes at the PM [31] and activate Gαiin the vicinity of ligand-activated EGFR (Module 2 in the model, Figure 2C, Table S2, S3). EGFR also triggers the activation of the PLC-γ-DAG-PKC-θ pathway [32]; PKC-θ phosphorylates GIV at S1689 and terminates GIV GEF activity towards Gαi [33]. Such sequential phosphorylation has another function – it converts GIV that is a GEF for Gαi(GIV-GEF) into GIV that now serves as a GDI for Gαs(GIV-GDI); GIV-GDI binds and inhibits GDP exchange on Gαs [22].
We asked, how do the CDK5 and the PLC-γ pathways regulate dynamics of the EGFR·GIV·Gαi com-plex formation, which is the key precursor event essential for transactivation of Gαiby EGF/EGFR [23,31]. 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 normalized density of the EGFR·GIV·Gαi complex formation. The temporal dynamics of these normalized densities of the EGFR·GIV·Gαi complexes generated from simulations were in good agreement with experimental measurements, as determined by FRET imaging [26,31] and GST pulldown assays [22,30] carried out in HeLa cells responding to EGF (Figure 2D).
Sensitivity analyses showed that despite the substantial number of model parameters (Tables S10 and S13), the formation of the EGFR·GIV·Gαi complex is sensitive only to a few kinetic parameters and initial conditions over time (Tables S10 and S13, Figure S2). For example, a ten-fold variation of the forward rate for the binding of GIV-GEF to the activated receptor (kf in reaction 15, Table S3) affected the peak values of the complex formation (Figure S2C) but not the temporal features of the EGFR·GIV·Gαi complex formation (Figure S2D). Similarly, the activation of the GIV-GEF function by CDK5 (reaction 14, Table S3) affected the density of the complex (Figure S2E) but not the temporal dynamics (Figure S2F).
The dynamics of the EGFR·GIV·Gαi complexes, however, were sensitive to the initial concentrations of Gαi (expected), GIV (expected) and PLC-γ (unexpected) (Table S10). The sensitivity of EGFR·GIV·Gαi complex formation to PLC-γ likely stems from network cross-talk, because the PLC-γ→DAG→PKC-θ pathway terminates GIV-GEF, triggering the dissociation of GIV and Gαi, which triggers the disassembly of the EGFR·GIV·Gαi complexes (Figures 1C, 1D, 2A). Changes in PLC-γ impacted both the density and temporal dynamics of the EGFR·GIV·Gαi complexes. As expected, when the PLC-γ→DAG→PKC-θ pathway is inhibited, the lifetime of GIV-GEF is prolonged and vice versa (Figure S2G). This effect is evident when comparing the normalized densities against experiments (Figure S2H).
We conclude that early activation of GIV-GEF, and the observed dynamics of the assembly of EGFR·GIV·Gαi complexes are not only dependent on the upstream kinetics of EGFR activation, but also on the downstream conversion of GIV-GEF to GIV-GDI, mediated by the PLC-γ→DAG→PKC-θ pathway. Findings also indicate that the connections within the network effectively capture the dynamics of transactivation of Gαi by EGFR via GIV-GEF.
Dynamics of Gαs activation is most compatible with delayed activation triggered by internalized EGFR and inactivation by GIV-GEM on endosomes
Although GIV-GDI inhibits the activity of Gαs-GTP [22], the exact mechanism of Gαs activation by EGFR is currently unknown. Prior studies have shown that Gαs is located on early, sorting and recycling endosomes [34] and that upon EGF stimulation, its activation/inactivation on endosomes regulates endosome maturation and EGFR degradation [35]; in cells without Gαs, or in those expressing a constitutively active mutant Gαs, internalized EGFR stays longer in endosomes, thereby, prolonging signaling from that compartment [28]. We asked when and where Gαs is activated. Because compartmentalized EGFR signaling (PM versus endosomes) occurs are different time scales (Figure 2B), and Gαi and Gαs have different timescales of activation [5 min and 15 min respectively] [22], we reasoned that computationally predicted dynamics of all three possible scenarios of compartmentalized Gαs activation i.e. [1) exclusively at the PM (Figure 3A, blue box); 2) exclusively at the endosomes (Figure 3A, red box); and 3) both at the PM and then on the endosomes (Figure 3A, both)], can provide insight into which option might be in accordance with the actual observed time scales for the same.
In the first scenario, where ligand-activated EGFR triggers Gαs activation exclusively at the PM, activation is predicted to be rapid, with peak concentration at 15 s; this kinetic pattern mimics the dynamics of rapid EGFR activation at the PM [compare blue line in Figure 3C with red line in Figure 2B]. In the second scenario, where ligand-activated EGFR triggers Gαs activation exclusively on endosomal membranes, the time of peak activity is around 15 min (Figure 3C, red line), in accordance with the time scales of Gαs activation and cAMP production [22]. And finally, if we consider a scenario where ligand-activated EGFR triggers Gαs activation both at the PM and on endosomes, we observe a first peak of rapid activation at around 15 sec, followed by a second burst at around 15 min In all three scenarios, activation of Gαs [concentration of Gαs-GTP] was higher in the absence of GIV’s GDI activity (i.e., when the concentration of GIV is set to zero; Figure 3C). Based on the dynamics of EGFR at the PM [rapid, almost instantaneous] and on the endosome [approximately 10 min] (Figure 2B) and similar timescales for Gαs activation observed from the different modes of Gαs activation (Figure 3C), we predict that Gαs is likely activated on the endosomes. It is noteworthy that the spatiotemporal features of such non-canonical Gαs signaling does not affect the kinetics of Gαi activation, either in the presence (Figure 3B, solid line) or absence of GIV (Figure 3B, interrupted line), indicating that EGFR transmodulates Gαi and Gαs independently.
To validate model predictions, we used a Gαs conformational biosensor, nanobody Nb37-GFP that binds and helps detect the nucleotide-free intermediate during Gαs activation [22]. In control cells, no significant Gαs activity was detected, either before or after ligand stimulation, indicating that Gαs is either not activated after EGF stimulation or that its activity is efficiently suppressed by some modulator, presumably GIV, for sustained periods of time. In GIV-depleted cells [80-85% depletion of endogenous GIV by shRNA sequence targeting the 3’ UTR [22], Gαs activity was easily detected roughly 15 min after ligand stimulation and exclusively on vesicular structures, likely to be endosomes (Figure 3D,E); no such signal was noted at the PM, which is where canonical activation of Gαs by GPCR is initiated [22]. These results obtained in live cells using conformation sensitive antibodies are in agreement with our in vitro enzymology assays [22] in that they confirm a role of GIV’s GDI function in the inhibition of Gαs activity, and confirm a much delayed and compartmentalized pattern of non-canonical cyclical activation/inactivation of Gαs downstream of EGF. As for what activates Gαs downstream of EGF/EGFR, few studies have shown that EGFR binds Gαs [36,37] through its juxtamembrane region [38], and that this interaction triggers phosphoactivation of Gαs [36]. Such transactivation of Gαs by EGFR in cardiomyocytes is accompanied by augmented AC activation, elevation of cAMP, increased heart rate and contractility [36, 39]. Our model neither proves nor disproves this model for direct transactivation of Gαs by EGFR, but reveals that activation of Gαi is delayed and pinpoints endosomes as the site of such activation.
Finally, we evaluated the dynamics of formation of the Gαs·GIV-GDI complex [Module 2; Figure 3F], the precursor event that is essential for transinhibition of Gαs by EGF/EGFR [22]. Our model for the dynamics of assembly of Gαs·GIV-GDI complexes (Table S2, S4, S5) included the kinetics of receptor internalization, Gαs activation by internalized receptors, conversion of GIV-GEF to GIV-GDI by the PLC-γ→DAG→PKC-θ pathway, and the Gαs-GDP-dependent degradation of endosomal EGFR (Figure 2A). Simulations from this model showed a good qualitative agreement between normalized Gαs·GIV-GDI complex formation between model and cell-based experiments [22], particularly until 50 min after EGF stimulation (Figure 3G). However, our model was unable to capture the normalized concentrations of Gαs·GIVGDI complexes at 60 min, exposing limitations of network modeling, i.e., missing interactions in the network that may operate specifically at later time points. One plausible group of unknown proteins that are missing in our model are downstream phosphatases that presumably act on GIV-GDI on endosomes, and are responsible for the decline in the number of Gαs·GIV-GDI complexes at later time points. The role of kinetic parameters and initial conditions affecting the formation of the Gαs·GIV-GDI complex were explored in detail (Figure S3) and we found that the dynamics of Gαs-GDI complex formation is more sensitive to internalization and degradation of EGFR than to other kinetic parameters (Figure S3, Table S14).
Compartmentalized modulation of Gαi and Gαs governs EGF-triggered cAMP 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 (Figure 4A, Table S5). Because EGFR triggers Gαi activation predominantly at the PM, we assumed that the PM-pool of Gαi inhibits the PM AC→cAMP pathway at the PM. Similarly, because Gαsis activated predominantly on endosomes and endosomal ACs [eACs] can be stimulated at that location to synthesize cAMP locally [40], we assumed that the endosomal-pool of Gαs likely stimulates the eAC→cAMP pathway (Table S5). To capture the dynamics of cAMP in our model network, we included such compartmentalized G protein-AC interactions (Figure 4A).
Our model shows that the early inhibition of cAMP is due to the Gαi-mediated inhibition of AC (the green regime) (Figure 4B); cAMP production is increased later due to the activation of Gαs on the endosome (the blue regime) (Figure 4B). These dynamics are consistent with previously published GIVdependent cAMP dynamics, measured by FRET [22]. 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 1C); such temporally separated compartmentalized modulation of two Gα-proteins with opposing effects on AC ensures suppression of cAMP at both early and later times during EGF signaling [22]. Because GIV modulates both Gαi and Gαs in different compartments and at 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 4B, compare GIV = 0 μM with GIV = 10 μM). As expected, sensitivity analyses confirmed that cAMP is sensitive to the initial concentrations of AC, Gαi, and Gαsand the reaction rates associated with AC, PKA, and PDE (Table S11, S15).
To understand how the relative contributions of the two G protein modulatory functions of GIV [GEF versus GDI] on cAMP production, we investigated cAMP dynamics in three conditions (Figure 4B) – 1) GEF-deficient but GDI-proficient (mimicked experimentally by the GIV-S1764D/S1689D mutant, GIV-DD [22], 2) both GEF- and GDI-deficient (mimicked experimentally by GIV-F1685A mutant, GIV-FA [22,41], 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 4B, dashed cyan line) as observed previously in cells expressing the GIV-DD mutant [22]. In the second scenario, where both GEF and GDI functions were lost, increase in cAMP concentration occurred early and such elevation was sustained (Figure 4B, dashed dark green line), as observed previously in cells expressing the GIV-FA mutant [22]; this mirrored the profile observed in GIV-depleted cells (Figure 4B, 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 4B, dot-dashed blue line).
While the dynamics of cAMP production provide 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 from time scale variations, and averages the effect of fluctuations in concentrations [42]. AUCs for cAMP at different time points were calculated to investigate how the cumu-lative cAMP signal varies under different GIV conditions (Figure 4D). 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 GIV-FA [a mutant that is defective in both GDI and GEF functions] as well as in the absence of GIV, i.e., it increases progressively through 60 min (Figure 4D, 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].
We next investigated the effect of cAMP dynamics on PICA and CREB. Activation of PICA by cAMP was modeled as a Hill function [43] and the fitting was compared against the experimental data as shown in Figure S4. The temporal dynamics of cAMP are reflected in the downstream dynamics of PICA and PDE (Figure 4E, F). The delay in PICA and CREB activation at short time scales [5 min] corresponds to the GIV-GEF mediated decrease of cAMP for control conditions (Figure 4E,F). As cAMP increases, PICA activity and CREB phosphorylation increase. In cases where there is no GIV-GEF mediated decrease in cAMP at the short time scale, there is no delay in PICA and CREB activation, resulting in an immediate and increased accumulation of PICA and CREB.
We also investigated how GIV concentration affects cAMP dynamics [output signal] with varying EGF/EGFR numbers [input signals]. When GIV concentrations were set to 0 μM in the model (to simulate cells that don’t have GIV), the network showed sensitivity, in that, increased input signals triggered increased output signals Figure 4F); this effect was even more pronounced in the absence of PDE (Figure 4J). Such sensitivity was replaced by robustness when GIV concentrations were set to high levels (GIV = 5 μM Figure 4G), i.e., increased input signals failed to initiate proportional output signals; this effect was virtually unchanged and robustness was preserved despite the absence of PDE (Figure 4K). These effects are also evident by studying the AUC (Figure 4H, L). These predictions were tested by measuring cAMP as determined by a radioimmunoassay (RIA) in control and GIV-depleted HeLa cells responding to varying doses of EGF [expermental 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 IBMX (Figure 4I, M). In the presence of GIV, cAMP production is robustly suppressed in response to increasing EGF ligand (Figure 4I). 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 4M). 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.
Clinical predictions from the model – from math to man
Concurrent upregulation of both GIV and EGFR maximally reduces cAMP and carries poor prognosis in colorectal tumors
We next asked how cAMP levels that are triggered by growth factors and modulated via GIV may impact tumor aggressiveness and clinical outcome. First, we conducted simulations for a wide range of GIV (0-5 μM) and EGFR (120 to 2400 EGFR molecules/μm2) concentrations to identify how crosstalk between these two variable components regulates cAMP levels (Figure 5A-D). Within each category of EGFR concentration [low or high], cAMP levels are the highest when GIV levels are lowest, and vice versa. In the setting of high EGFR expression (Figure 5B), the impact of changing GIV was the highest, i.e., the range of cAMP response was the widest. By contrast, in the setting of low EGFR expression (Figure 5A) the impact of changing GIV on the cAMP levels was minimal. Because the EGFR copy number and GIV copy number variables in our model are a surrogate for EGFR or GIV signaling states, respectively, which cumulatively represents all perturbations that can impact the functions of both EGFR and GIV [i.e., activating mutations, phosphomodifications, gene copy number variation, or simply protein overexpression in tumors], we conclude that the Gαi/Gαs/cAMP modulatory function of GIV exerts its maximal impact in the setting of high EGFR signaling states. These findings indicate that the two modules [EGFR and GIV-dependent Gαs/Gαi/cAMP signaling] are intertwined via a functional crosstalk.
In order to quantify the extent of this crosstalk, we conducted simulations for a wider range of EGFR [48-2400 molecules/μm2] and GIV concentrations [0-10 μM] and calculated the AUC for the cAMP dynamics (Figure 5C, D). Varying GIV concentrations resulted in cAMP changes within a narrow range in low EGFR state, but showed larger variance in the high EGFR state (Figure 5C). Because cAMP has potent anti-cancer function, these findings point to the possibility that there are regimes of operation in the EGFR-GIV space that can be exploited for prognostication in cancers. To further dissect this space, we plotted the variations in AUC for EGFR and GIV variations (Figure 5D). The value of AUC corresponding to the control (GIV 1μM and EGFR 120 molecules/μm2, Table S9) is approximately 0.5 μM · min and is denoted by the yellow color and 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 for low GIV concentrations. But as GIV concentrations increased, even at high EGFR, the cAMP AUC decreased, indicating that the impact of increasing GIV on cAMP AUC was higher than the impact of increasing EGFR. Therefore, GIV levels, rather than EGFR copy number alone, can be thought of as a prognosticator for decreased cAMP levels.
To determine the impact of crosstalk 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 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; [44]) (Figure 5E). 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 compared to those with low GIV (Figure 5F). Among patients with low EGFR, expression of GIV at high or low levels did not impact survival (Figure 5G). Conversely, among patients with high levels of GIV, survival was significantly different between those with high versus low EGFR; no such trend was noted among patients with low GIV (Figure 5). Thus, the high GIV/high EGFR signature carried a poorer prognosis compared to all other patients. More importantly, patients with tumors expressing high EGFR did as well as those expressing low EGFR provided the levels of GIV in those tumors was low. Consistent with the fact that cAMP is a potent anti-tumor second messenger, these findings reveal that 1) high levels of EGFR signaling does not, by itself, fuel aggressive traits or carry a poor prognosis, but does 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, critical for dampening several aggressive tumor traits [Figure 1A].
The crosstalk between EGFR and GIV the we define here, and its impact on clinical outcome provide a plausible explanation for some long-standing 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 [21]. Although activating EGFR mutations, copy number variations, and levels of EGFR protein expression seem to be closely related to each other [45], the prognostic impact of EGFR expression in cancers has been ambiguous [46]. In some cancers, high EGFR copy numbers is associated with poor outcome [47, 48]; in others, high EGFR expression unexpectedly favors better overall and progression-free survival [49–51]. Thus, due to 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 [19]. Our findings that GIV levels in tumors with high EGFR is a key determinant of the levels of the anti-tumor second messenger cAMP, have provided a potential molecular basis for why elevated EGFR signaling in some tumors can be a beneficial in some, but a driver of metastatic progression in others. Because cAMP levels in tumor cells and GIV levels have been previously implicated in anti-apoptotic [52] and the development of chemoresistance [53], it is possible that the GIV-EGFR crosstalk 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 such is the case, remains to be evaluated.
Concurrent downregulation of both GIV and PDE activity maximally increases cAMP, carries good prognosis in colorectal tumors
Another factor that plays an important role in regulating cAMP levels is PDE. We next conducted simulations for different GIV (0 to 1 μM) and PDE (0.04 to 2 μM) concentrations to identify how the crosstalk between these two variable components regulates cAMP levels. Within each category of PDE concentration [low vs high PDE states; (Figure 6A-B)], cAMP levels are the highest when GIV levels are lowest, and vice versa. In the setting of low PDE activity (Figure 6B), the impact of changing GIV was the highest, i.e., the range of cAMP response was the widest. By contrast, in the setting of high PDE activity (Figure 6A) the impact of changing GIV on the cAMP levels was minimal. These effects can be seen also when comparing the AUCs for the low vs high PDE states, calculated over 1 hr (Figure 6C). While there is no significant change in the AUC with increasing GIV in a high-PDE state (red bars), increase in GIV leads to a decrease in cAMP in 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 6D) 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, 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 high PDE state virtually abolishes all effects of GIV-dependent inhibition of cAMP production, we also conclude that in this GIV-PDE crosstalk, PDE is a dominant node and GIV is the subordinate node. Because the ‘PDE copy number’ variable in our model is a surrogate for PDE activity/signaling, which cumulatively represents all perturbations that can impact the functions of PDEs in cells [i.e., phosphomodifications, mislocalization, or simply protein overexpression / hyperactivation in tumors; [54]], we conclude that the Gαi/Gαs modulatory function of GIV exerts its maximal impact on whether cAMP is high or low in the setting of low PDE signaling states.
To determine the impact of crosstalk between various PDE isoforms and GIV on clinical outcome, we carried out using the StepMiner algorithm, implemented within the Hegemon software on the same set of 466 patients with colorectal cancers as before, except patients were now stratified into low and high subgroups with regard to GIV (CCDC88A) and PDE gene-expression (Figure 6E). Among the 11 known PDE isoforms, we evaluated those that have previously been linked to colon cancer progression [PDE5A (Figure 6E-G), 4A and 10A (Figure S9)]. 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 high PDE state (Figure 6F) but highly significant in low PDE state (Figure 6G). Thus, the low GIV/low PDE signature carried a better prognosis compared to 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, critical for dampening several aggressive tumor traits (Figure 1A).
In the context of PDE, it has been demonstrated that overexpression of PDE isoforms in various cancer leads to impaired cAMP and/or cGMP generation [55]. PDE inhibitors in tumour models in vitro and in vivo have been shown to induce apoptosis and cell cycle arrest in a broad spectrum of tumour cells [56]. Despite the vast amount of preclinical evidence, there have been no PDE inhibitors that have successfully translated to the 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 carries 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.
Conclusions
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 non-canonical 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 over an hour, which is distinct from the canonical GPCR/G protein pathway; GPCRs initiate more rapid and finite cAMP dynamics in the order of msec to min (Figure 7A) [57]. 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 (Figure 7A) [58]. By contrast, in the case of RTK-mediated cAMP dynamics via GIV-GEM, the post-endocytic (i.e., endosomal signaling) component constitutes a dominant component of the overall cAMP dynamics that are triggered by RTKs (Figure 7A). 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 min to > 60 min, which coincides with other RTK signaling, trafficking events and transcriptional response, i.e., the major temporal domain of RTK activity, the so-called “window of activity” [59]. The 5 min to 1 hour 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 hours); these transcriptional targets not only generate feedback within the RTK-signaling cascade, but also set up crosstalk with other signaling pathways [59, 60]. In fact, GIVGEM has indeed been found to modulate myriad downstream signaling pathways from the activity of small GTPases, kinases and phosphatases, to transcription factors [reviewed in [18]]; how GIV-GEM has such a widespread and broad impact had remained a mystery. It is possible that such 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 our model and simulations were carried out using the prototype RTK, EGFR, fundamentals identified here are likely to be relevant also in the case of others RTKs, e.g., VEGFR, IGF-1R, InsulinR, PDGFR, etc. Each of these RTKs are known to engage GIV and rely on it’s G protein modulatory function [reviewed in [18]] and can signal both at the PM and within endosomes [61–65]. Future work will explore the impact of such multi-receptor integration by GIV-GEM.
Second, our network model has helped us identify key design principles of the action of GIV-GEM within signaling circuits by enabling the construction of a map to identify the relationship between input [EGFR]→value[GIV]→output[cAMP]→sink[PDE] relationship (Figure 7B). The crosstalk between the input, the valve, and the sink 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. The maps shows that the 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 7B-E). At low concentrations of GIV, such as those found in normal tissues [Figure S11, S12], cAMP levels are sensitive to increased signal input via EGF/EGFR, i.e., higher input elicits higher output. Such sensitivity is virtually abolished and replaced by robustness at higher GIV concentrations found in a variety of cancers [Figure 7 B, C; S11, S12], i.e., higher input fails to elicit higher output and instead, cAMP levels stay at low and relatively constant. Regulation of cellular cAMP concentration by the EGFR-GIV interplay appears to be dependent on PDE concentration (Figure 7D); when PDE is high, there is virtually little or no effect of changes in EGFR or GIV on cAMP concentrations. The PDE-GIV relationship for cAMP production responds to different EGFR inputs proportionally, with increasing EGFR copy numbers resulting in higher cAMP AUC (Figure 7E). This 3-way interplay between EGFR, GIV and PDE is obvious also in experimental data derived from HeLa cells (Figure 4, I, M). Heat maps derived from that data (Figure 7F-G) show that the EGFR→cAMP pathway is most sensitive [i.e., higher input (EGF) signal, produces higher output (cAMP)] when the activities of both GIV and PDE are at their lowest (shGIV, with PDE inhibition; Figure 7G). Conversely, the EGFR→cAMP pathway is most robust [i.e., output (cAMP) is maintained at low levels despite higher input (EGF) signal] when both GIV and PDE are high (shC, without PDE inhibition; Figure 7F). When GIV is low and PDE is high (shGIV, without PDE inhibition; Figure 7G), cAMP levels do not go up, likely because increased production is balanced by increased degradation. Why would a cell waste energy (ATP) in such a ‘futile cycle’? This situation is reminiscent of the maintenance of steady-state cGMP levels in the sub-μM range in thalamic neurons by concomitant guanylyl cyclase and PDE2 activities [66] and cAMP levels in pyramidal cortical neurons by concomitant AC and PDE4 activities [67]. Prior studies have suggested that such tonic cAMP production and PKA activity enable signal integration and crosstalk with other cascades [68]; unlike an on/offsystem gated exclusively by Gαsproteins, tonic activity allows both up- and downregulation by activation of Gαi or inhibition of Gαs (via GIV-GEM) and by PDEs. Our findings suggest that such up/down tunability is best achieved by changing the cellular concentrations of GIV. Because the flow of information in layers within signal transduction circuits in general [69–71], and more specifically for RTKs like EGFR [59, 72] is believed to conform to bow-tie microarchitecture, and cAMP is considered as one of the universal carrier molecules at the knot of such bow-ties which determines robustness [71], we conclude that GIV-GEM operates at the knot of the bow-tie as a tunable valve for controlling robustness within the circuit. Because layering of control of [information] flow is believed to conform to an hourglass architecture [70], 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 [18]. In an hourglass architecture, the lower and higher layers tend to see frequent evolutionary changes, while the carriers at the waist of the hourglass appear to be constant/invariant and sometimes, virtually ‘ossified’. Of relevance to our model, the importance of cAMP appears to be indeed ossified from unicellular organism to (wo)man alike, and GIV-GEM is expressed in ubiquitously in all tissues from fish to (wo)man and GIV-like GEMs have so far been identified as early as in C. elegans [73].
Third, our work also provides valuable clues into the impact of increased robustness at 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 like cancers; it is widely accepted that tumor cells hijack such robustness to gain growth and survival advantage and during the development of cancer [59, 74, 75]. Consistently, we found that GIV mRNA levels and DNA copy numbers are invariably higher across multiple cancers when compared to their respective normal tissue of origin (Figure S11, S12). Because GIV has been found to regulate several sinister properties of tumor cells across a variety of cancers (multiple studies, reviewed in [20]), 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. Such a phenomenon could be a 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.
Although our model captures experimentally observed time courses and generates testable hypotheses, it has 3 major limitations. First, 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 [76, 77]. Second, 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. Third, our model focuses exclusively on EGFR and does not account for the diverse classes of receptors [multiple RTKs, GPCRs, integrins, etc.] also use GIV to access and modulate G proteins. Despite these restrictions, we can identify some fundamental features of growth factor-triggered cAMP signaling for the first time using systems biology, including the role 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.
We conclude that GIV utilizes compartmental segregation to modulate the dynamics of RTK-G protein-cAMP signaling and confers robustness to these dynamics by functioning as a tunable control valve. Future systems efforts will build on this model to unravel further exciting features of GIV as a critical hub for signaling regulation at the knot of a bowtie [71] and elucidate the hidden complexity that arises from network architecture in non-canonical G protein signaling.
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 1D). We constructed the compartmental computational model in a modular manner, where each module represents the key events in the network. 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 function. The modules are as follows.
Module 1 consists of EGFR activation through EGF and internalization dynamics leading to the formation of EGFR·GIV·Gαi trimeric complex. This module includes the phenomenon that endosomal maturation and EGFR degradation in lysosomes requires the presence of inactive Gαs [GDP-bound state] [28]. 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 [28]. 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 [28] (Figure 2A).
Module 2 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 [33]. This phosphoevent does not impact GIV’s ability to inhibit Gαs [22]. Consequently, when it comes to G protein modulatory functions of GIV, phosphorylation by PKC-θ converts GIV-GEF into GIV-GDI (Figure 3F).
Module 3 contains the dynamics of the endosomal EGFR and how it activates Gαs, and subsequently adenylyl cyclase (AC) leading to the synthesis of cAMP (Figure 3A). Overall, the model contains 57 reactions. The complete set of reactions for each of the modules, their parameters and interactions, and the list of assumptions underlying network construction are provided as online supplementary materials (Tables S2 – S5).
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. Such assumption allowed us to generate a deterministic dynamical model. The model contains six different compartments: (i)PM, (ii) extracellular space, (iii) cytosol, (iv) endosomes, (v) endosomal membranes, and (vi) nucleus. 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 can be tracked. Each interaction was modeled as a chemical reaction either using mass-action kinetics for binding-unbinding reactions, and Michaelis-Menten kinetics for enzyme-catalyzed reactions, as is standard for models such as this [78, 79].
The network of interactions was constructed using the Virtual Cell modeling platform (http://www.nrcam.uchc.edu). We chose this platform because it is a user-friendly computational cell biology software, 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 high degree of numerical accuracy [80–83]. Also, the Virtual Cell platform has built-in capabilities to conduct dynamic sensitivity analysis, which is an important aspect of dynamical 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
In order to characterize the dynamics of the different protein activities, we use the area under the curve for the concentration versus time curve [42]. The area under the curve gives the total signal activated over the time of observation and for the ith species is given by AUCi in Eq. 1. This gives a measure of the total signal for different conditions.
Comparison with experimental data
The experimental data was extracted from Figure 1D of [22] using ImageJ for GEF-Gαi-EGFR complex and Gαs-GDI complex. The data was normalized such that the maximum value was 1. Parameter fitting using COPASI [84] was used to then match the normalized experimental data against the model output. Goodness of fit between experimental values and model output was determined using a root mean squared error (RMSE).
Dynamic parametric sensitivity analysis
Since 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 Ci, with respect to parameter kj is given by [85,86]
Since we are studying a dynamical system and not steady state behavior, we used the Virtual cell software to calculate the change in log sensitivity over time . The resulting time course gives us information about the time dependence of parametric sensitivity coefficients for the system. The variable of interest, Ci is said to be robust with respect to a parameter kj if the log sensitivity is of the order 1 [85]. We refer the reader to [85, 86] for a complete introduction to dynamical sensitivity analysis. We conducted dynamic sensitivity analysis for all the kinetic parameters for the reactions and initial concentrations of the different species in the model.
Measurement of cAMP
HeLa cells were serum starved (0.2 % FBS, 16 h) and incubated with isobutylmethylxanthine (IBMX, 200 μM, 20 min) followed by EGF (60 min). Reactions were terminated by aspiration of media and addition of 150 μl of ice-cold TCA 7.5% (w/v). cAMP content in TCA extracts was determined by radioimmunoassay (RIA) and normalized to protein [(determined using a dye binding protein assay (Bio-Rad)] [18, 87]. Data is 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 cohort of 466 patients where each tumor had been annotated with the disease-free survival (DFS) information of the corresponding patient. This cohort included gene expression data from four publicly available NCBI-GEO data-series (GSE14333, GSE17538, GSE31595, GSE37892) [88–91], 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, Florida, USA (n = 164); 2) the Vanderbilt Medical Center in Nashville, Tennessee, USA (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 datasets) 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 dataset has been published previously [44]. 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 dataset, we applied the Hegemon software tool [44]. The Hegemon software is an upgrade of the BooleanNet software [92], 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 logrank tests. Since all 466 samples contained in the dataset 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 NCBIGEO discovery dataset in 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), 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.
Acknowledgments
This work was supported by ARO W911NF-16-1-0411, AFOSR FA9550-15-1-0124, and NSF PHY-1505017 grants to P.R. and NIH grants CA100768, CA160911 and DK099226 (to P.G). M.G. was supported by the UCSD Frontiers of Innovation Scholars Program (FISP) G3020. L.S. was supported by T32DK0070202. The Virtual Cell is supported by NIH Grant Number P41 GM103313 from the National Institute for General Medical Sciences.
References
- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].
- [7].
- [8].↵
- [9].↵
- [10].
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].
- [51].↵
- [52].↵
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
- [59].↵
- [60].↵
- [61].↵
- [62].
- [63].
- [64].
- [65].↵
- [66].↵
- [67].↵
- [68].↵
- [69].↵
- [70].↵
- [71].↵
- [72].↵
- [73].↵
- [74].↵
- [75].↵
- [76].↵
- [77].↵
- [78].↵
- [79].↵
- [80].↵
- [81].
- [82].
- [83].↵
- [84].↵
- [85].↵
- [86].↵
- [87].↵
- [88].↵
- [89].
- [90].
- [91].↵
- [92].↵