Deregulation of ErbB signaling plays a key role in the progression of multiple human cancers. To help understand ErbB signaling quantitatively, in this work we combine traditional experiments with computational modeling, building a model that describes how stimulation of all four ErbB receptors with epidermal growth factor (EGF) and heregulin (HRG) leads to activation of two critical downstream proteins, extracellular‐signal‐regulated kinase (ERK) and Akt. Model analysis and experimental validation show that (i) ErbB2 overexpression, which occurs in approximately 25% of all breast cancers, transforms transient EGF‐induced signaling into sustained signaling, (ii) HRG‐induced ERK activity is much more robust to the ERK cascade inhibitor U0126 than EGF‐induced ERK activity, and (iii) phosphoinositol‐3 kinase is a major regulator of post‐peak but not pre‐peak EGF‐induced ERK activity. Sensitivity analysis leads to the hypothesis that ERK activation is robust to parameter perturbation at high ligand doses, while Akt activation is not.
The ErbB receptors, ErbB1 (EGFR), ErbB2 (HER2/NEU), ErbB3, and ErbB4, are widely expressed throughout the human body and control cell fate decisions, such as proliferation, apoptosis, and differentiation (Yarden and Sliwkowski, 2001). Given the wide physiological expression of ErbB receptors and their importance in controlling cell fate decisions, it is not surprising that deregulation of ErbB signaling is implicated in the progression of multiple human cancers (Yarden and Sliwkowski, 2001). Because the ErbB signaling network is a complex, dynamic system with multiple feedbacks, it is difficult to understand its overall signaling behavior through traditional, qualitative methods. In the current work, we combine traditional experimental methods with quantitative modeling to understand ErbB signaling in MCF‐7 breast cancer cells. We specifically focus on the short‐term (⩽30 min) dynamic behavior of all four ErbB receptors and the key downstream intermediates ERK and Akt in response to stimulation with the ligands epidermal growth factor (EGF) and heregulin (HRG) (Figure 1). Based on the schematic in Figure 1, we build an ordinary differential equation model consisting of 117 species, 235 parameters, and 96 net reactions to represent the dynamics of ligand‐dependent ErbB network activation. We use this model to generate hypotheses about how the cell dynamically controls ERK and Akt activation in a ligand‐dependent fashion, and experimentally corroborate the predictions that EGF‐induced ERK activity is much more sensitive than HRG‐induced ERK activity to the MEK inhibitor U0126 (Figure 7A and B), and that the phosphoinositol‐3 kinase (PI‐3K) inhibitor wortmannin has a large affect on post‐ but not pre‐peak EGF‐induced ERK signaling.
Gaining insight into the behavior of such a large model for the purposes of hypothesis generation is a difficult task, and in this work, we relied mainly on parametric sensitivity analysis to accomplish this. Naturally, the sensitivity analysis results revealed many obvious modes of regulation, such as control of ERK activation by MEK, and control of Akt activation by PI‐3K. However, these results also revealed non‐obvious, ligand‐dependent modes of regulation: (i) ERK activation is robust to parametric perturbations at high ligand doses, while Akt activation is not. This difference in robustness to parameter variation may be due to the involvement of ERK in multiple negative feedback loops, and the absence of such feedback from Akt. It is well known from control theory that negative feedback endows a system with robustness to disturbances (Ogunnaike and Ray, 1994; Freeman, 2000). (ii) PI‐3K abundance is a dominant factor for the post‐peak (declining portion of the ERK response curve) but not pre‐peak (rising portion of the ERK response curve) EGF‐induced ERK activation. This phenomenon is due to degradation of active ErbB1 homodimers. At short times, there are still many 1‐1 homodimers that can signal to ERK; as time progresses, however, these homodimers are degraded. Therefore, at longer times, signaling to ERK must rely on alternative mechanisms—in this case, the PI‐3K‐Grb2‐associated binder 1 (Gab1) pathway.
The results of the sensitivity analysis allowed us to gain general understanding of the ligand‐dependent control of ERK and Akt activity. While this general understanding is valuable by itself, understanding how known, cancer‐correlated network abnormalities affect signaling is of primary and immediate interest. We used the model to investigate how such a network abnormality, ErbB2 overexpression, which occurs in approximately 25% of breast cancer, would affect the dynamics of ERK and Akt activation. Simulation results suggested that it should have a large effect on signaling at high EGF doses by transforming normally transient signals into more sustained signals. Significantly, this model prediction quantitatively agrees with the data of Wolf‐Yadlin et al (2006) and Kumar et al (2007), which shows that at 10 and 30 min after EGF stimulation, ERK activation is between 1.15‐ and 2‐fold higher for human mammary epithelial cells overexpressing ErbB2. This effect of ErbB2 overexpression is a direct consequence of ligand‐induced ErbB1 homodimer internalization and degradation. EGF transmits signals only through ErbB1 homodimers and ErbB1‐ErbB2 heterodimers, and ErbB1 homodimers undergo preferential ligand‐induced internalization and degradation relative to other ErbB dimers (Baulida et al, 1996). As ErbB2 overexpression shifts the ErbB dimer distribution toward more ErbB1‐ErbB2 heterodimers rather than ErbB1 homodimers, it therefore sustains signaling due to decreased internalization and degradation.
Our computational analyses of the ErbB signaling network generated a number of hypotheses regarding how MCF‐7 cells control ERK and Akt activities in a dynamic and ligand‐dependent fashion. These hypotheses can and should be tested against the experiments, and their results will help us and others refine and improve the model. In this study, we performed two experiments to test our hypotheses. To test the hypothesis that HRG‐induced ERK activation should be robust to the MEK inhibitor U0126, while EGF‐induced ERK activation should not, we performed a U0126 titration on 5‐min ERK activity (Figure 7A and B). These results confirmed the prediction that HRG‐stimulated ERK activation was more robust than EGF‐stimulated ERK activation to U0126. We experimentally investigated the effect of the PI‐3K inhibitor wortmannin on ERK signaling, and the results corroborated our prediction that PI‐3K is a dominant factor for post‐peak but not pre‐peak EGF‐induced ERK signaling. However, from the wortmannin experiments we also discovered limitations of our current model, which suggested directions for future model improvements, such as incorporation of the tyrosine phosphatase SHP2.
Building a quantitative model of the ErbB signaling network that incorporates stimulation of all four ErbB receptors with two ligands simultaneously was a difficult endeavor. This endeavor was worthwhile, however, as we generated many hypotheses that were only realizable using model‐based analysis. Experimental investigation of all these hypotheses is beyond the scope of a single study; however, over time the combined results of many studies can be used to test these hypotheses. When conflicts arise, our current understanding and model of the ErbB signaling network will be revised accordingly.
We developed a quantitative kinetic model that relates EGF and HRG stimulation of all four ErbB receptors to ERK and Akt activation in MCF‐7 breast cancer cells.
Experiments and model simulations show that EGF‐induced responses are transient while HRG‐induced responses are sustained. Significantly, model analysis shows that ErbB2 overexpression transforms transient EGF‐induced signaling into sustained signaling, which is corroborated by the experimental data.
Model simulations predict and experiments confirm that HRG‐induced ERK activity is more robust to the ERK cascade inhibitor U0126 than EGF‐induced ERK activity.
We hypothesize through model analysis that PI‐3K is a major regulator of post‐peak but not pre‐peak EGF‐induced ERK activity, and confirm this hypothesis experimentally using the PI‐3K inhibitor wortmannin.
The ErbB signaling network, a major player in tumorigenesis, comprises the following: multiple extracellular ligands; the four trans‐membrane receptors ErbB1 (EGFR), ErbB2 (HER2/NEU), ErbB3, and ErbB4; cytoplasmic adapters, scaffolds, enzymes, and small molecules. Signaling is initiated when ligand binds to a receptor and causes the receptors to homo‐ or heterodimerize. Receptor dimerization activates the receptor's tyrosine kinase domain, which leads to autophosphorylation of tyrosine residues on receptor cytoplasmic tails. Multiple cytoplasmic adapter, scaffold, and enzymatic proteins are then recruited to the plasma membrane by binding to receptor phosphotyrosines. A complex network of interactions between the activated receptors, recruited proteins, and plasma membrane molecules eventually culminates in the activation of multiple downstream effectors, including extracellular‐signal‐regulated kinase (ERK) and protein kinase B/Akt, which are implicated in the control of proliferation and survival.
Abnormalities within the ErbB signaling network correlate with the development of several cancer types, and multiple drugs that target these defects have been used to treat cancer successfully (Yarden and Sliwkowski, 2001). For example, in the approximately 25% of patients whose metastatic breast cancer overexpresses the ErbB2 receptor, median patient survival time is improved by administering the ErbB2‐targeted, monoclonal antibody trastazumab (Hortobagyi, 2001). Although knowledge of some specific ErbB signaling network defects associated with tumorigenesis has led to the development of successful cancer treatments, these targeted treatments are rarely ‘magic bullets’. Furthermore, there are instances where potentially drug‐sensitive cancers either do not respond and/or eventually become resistant to treatment (Robinson et al, 2006). Improving the efficacy of these targeted treatments requires a more detailed understanding of the mechanisms by which cancer‐correlated network properties cause deregulation of the entire ErbB signaling network.
Spatio‐temporal signaling aspects play a key role in the ErbB network's control of cell fate, as different inputs stimulate different activation kinetics, which ultimately lead to different cell fates. For example, in AU565 cells, stimulation with the ligand heregulin (HRG) causes sustained network activation and leads to differentiation (Lessor et al, 1998); in PC12 cells, stimulation with the ligand epidermal growth factor (EGF) causes transient network activation and leads to proliferation (Marshall, 1995). In MCF‐7 cells, Thottassery et al (2004) showed that EGF and HRG cause transient and sustained network activation, respectively. Although it is clear that (i) different ErbB ligands can stimulate different network activation dynamics, and (ii) that there is a connection between ligand‐dependent activation kinetics and cell fate, to understand how the ErbB signaling network controls cell fate, we must first elucidate the mechanisms that control ligand‐dependent activation kinetics. Similarly, understanding ligand‐dependent signaling mechanisms is a key step in understanding how the ErbB network's deregulation contributes to tumorigenesis.
Because the ErbB signaling system is a highly interconnected, dynamic network containing multiple feedback loops, it is difficult to predict the response of the network solely by qualitative means. It is becoming increasingly clear that quantitative methods are required to understand the mechanisms by which signaling networks function. Therefore, in this work, we take a combined experimental and computational model‐based approach to understand the ErbB network that was pioneered by Kholodenko et al (1999), and expanded upon by Schoeberl et al (2002), Hatakeyama et al (2003), Hendriks et al (2003), Resat et al (2003), Blinov et al (2006), Shankaran et al (2006), and many others. This approach employs a combination of mechanistic, ordinary differential equation (ODE) modeling (for simulation) with quantitative immunoblotting (for experimental measurements of signaling dynamics).
Current methods for dynamic modeling of the interactions between proteins that contain multiple phosphorylation sites and binding domains requires dealing with a combinatorial explosion of potential species, significantly complicating the development and simulation of signaling network models. For example, a mechanistic description of the ErbB1 receptor that simultaneously accounts for the ligand‐binding domain, the dimerization site, the kinase domain, and 10 phosphorylation sites requires more than 106 differential equations. This phenomenon, referred to as ‘combinatorial complexity’, is a fundamental problem in developing mechanistic, differential equation models of signal transduction networks (Goldstein et al, 2004; Blinov et al, 2006). Due to combinatorial complexity, previous models of ErbB signaling were either limited to a single ErbB receptor and ligand (Kholodenko et al, 1999; Schoeberl et al, 2002; Hatakeyama et al, 2003) or only accounted for complexity at the level of ligand binding and receptors (Hendriks et al, 2003, Shankaran et al, 2006). Here, we develop a novel technique for handling the problem of combinatorial complexity, allowing us to incorporate a greater number of ErbB network entities than in previous models.
Using a combined quantitative modeling and experimental approach, in the current work we investigate the short‐term (⩽30 min) response of the key ErbB signaling network intermediates ERK and Akt to stimulation with the EGF and HRG ligands, in MCF‐7 breast cancer cells. We build a model to represent the measured dynamics of ligand‐dependent ErbB network activation, use the model to generate hypotheses about the control of ERK and Akt activation, and experimentally corroborate the prediction that EGF‐induced ERK activity is much more sensitive than HRG‐induced ERK activity to the ERK cascade inhibitor U0126. Through model analysis, we find that control of ERK and Akt activation is distributed among many cellular mechanisms, and the responsible mechanisms change as a function of time and ligand dose combinations.
To help understand ErbB network signaling, we developed a computational model that relates EGF and HRG stimulation of the ErbB receptors to activation of ERK and Akt. It is an ODE model that consists of two compartments (extracellular and cytoplasmic), 117 species, 235 parameters, and 96 (net) reactions (details of the rate equations, parameter values, and parameter estimation procedures can be found in the Supplementary information online). The model, although reasonably large, is not a combinatorially complex, in silico replica of all potential distinct biochemical species and processes. Such a microscopically comprehensive model would be impractical to develop, both computationally and experimentally. The goals for this model are to reflect the experimental data measured in this study to help provide insight into mechanisms that drive the observed phenomena. In this regard, our goals are similar to the goals of those who developed previous models of ErbB signalling. A simplified schematic representation of the model structure is shown in Figure 1, the reaction network is shown in Figure 2, and the model is described as follows.
Ligand binding and dimerization.
EGF has high affinity for ErbB1, HRG has high affinity for both ErbB3 and ErbB4, and no natural ligand is known for ErbB2. Ligand‐bound ErbB1, ErbB3, and ErbB4 can dimerize with other ligand‐bound ErbB1, ErbB3, or ErbB4, whereas ErbB2 is constitutively dimerization prone. Because ErbB2 is constitutively dimerization competent, it typically is referred to as the preferred dimerization partner in the ErbB family and tends to form heterodimers with other ErbB family members (Graus‐Porta et al, 1997). Therefore, ligand‐bound ErbB1, ErbB3, and ErbB4 and free ErbB2 can form homo‐ and/or heterodimers that consist of 1‐1, 1‐2, 2‐3, 2‐4, 3‐4, 4‐4, 1‐3, and 1‐4 receptors. We neglect 3‐3 dimers, because Berger et al (2004) showed that these dimers do not form, and further, ErbB3 receptor is kinase dead (Citri et al, 2003). To our knowledge no studies have shown that ErbB2 homodimers will initiate significant levels of signaling in the absence of other ErbB receptors. However, wild‐type ErbB2 can form a small population of homodimers, and the E664V mutation found in some ErbB2‐dependent cancers can cause a large amount of ErbB2 homodimers to form (Burke et al, 1997). In cells bearing this mutation, ErbB2 homodimers may be important to consider, as ErbB2 homodimerization might represent a sequestration mechanism that combats ErbB2 overexpression. However, MCF‐7 cells are estrogen receptor positive (indicative of early stage breast cancers), do not depend on or overexpress ErbB2, and we are unaware of any experimental evidence to indicate that the E664V mutation associated with later‐stage, ErbB2‐dependent breast cancers is present in MCF‐7 cells. As Burke et al (1997) showed that only approximately 5% of all wild‐type ErbB2 dimers exist in oligomeric form, sequestration of ErbB2 through homodimerization should have minimal impact on signaling in MCF‐7 cells, and therefore we neglect 2‐2 homodimers.
Receptor dimer autophosphorylation and the ‘virtual phosphorylation site’.
Once a receptor dimer is formed, it gains tyrosine kinase activity and can autophosphorylate on several tyrosine residues. Simultaneously accounting for all these phosphorylation sites results in a combinatorial explosion of potential species, thus, we represent all autophosphorylation sites as a single ‘virtual phosphorylation site’ as similar to previous models of ErbB signaling (e.g. Kholodenko et al, 1999). To reduce the model in this manner, we must assume that that either all or no sites are phosphorylated at a particular time. This assumption may be close to reality, as the phosphorylation and dephosphorylation steps are fast relative to downstream events. Phosphorylation is fast because the receptor cytoplasmic tails (substrate) are tethered to the dimer (kinase), concentrating substrate and enzyme in an extremely small volume. Because the tyrosine phosphatase PTB‐1B is recruited to ErbB receptor dimers and activated in a ligand‐dependent fashion (Liu and Chernoff, 1997), it is also likely that dephosphorylation is fast. Also supporting this assumption, autophosphorylation sites on ErbB1 and ErbB2 show very similar phosphorylation dynamics in human mammary epithelial cells (HMECs) in response to EGF and HRG (Wolf‐Yadlin et al, 2006), and autophosphorylation sites on ErbB1 show very similar phosphorylation kinetics in HeLa cells in response to EGF (Schulze et al, 2005).
Membrane recruitment and ligand‐induced degradation.
Membrane recruitment of signaling proteins to the plasma membrane is an early step of signal propagation through the ErbB network, and is critical for signal propagation, because it colocalizes key network enzymes (Kholodenko et al, 2000). Following receptor phosphorylation, we consider that the proteins Shc, Grb2, phosphoinositol‐3 kinase (PI‐3K), PTP‐1B, and RasGAP are recruited by binding directly to tyrosine‐phosphorylated receptors. We assume that PI‐3K cannot bind directly to ErbB1 or ErbB2 (Citri and Yarden, 2006). The proteins SOS and Gab1 are secondary recruits, and have proline‐rich domains (PRD), which bind to the N‐terminal or C‐terminal Src homology 3 (SH3) domains of Grb2, respectively (Figure 3). Gab1 can also be recruited to the membrane via its pleckstrin homology (PH) domain binding to PIP3 (Figure 3). Following membrane recruitment, Shc, RasGAP, and Gab1 are tyrosine phosphorylated by receptor dimers, and then Shc and Gab1 can recruit additional downstream proteins. We represent the multiple Gab1 tyrosine phosphorylation sites as a virtual phosphorylation site.
ErbB1 homodimers undergo rapid EGF‐induced internalization and subsequent degradation through a multistep process, whereas other dimers stay much longer at the plasma membrane (Sorkin et al, 1993; Baulida et al, 1996). For this model, we approximate the rate of ligand‐induced 1‐1 dimer degradation as a first order process, which is reasonable as the initial internalization and recycling steps are fast relative to the subsequent degradation step (Wiley, 2003) (see Supplementary information online for derivation).
Reducing the complexity of membrane recruitment.
The receptor cytoplasmic tails, which contain the tyrosine phosphorylation sites, are of similar size to their downstream binding proteins, and the tyrosine phosphorylation sites considered on Gab1 are in close proximity. It is therefore likely that steric hindrance will interfere with the simultaneous binding of more than one adapter protein to a virtual phosphorylation site. Thus, we represent binding to virtual sites as fully competitive (mutually exclusive) between all downstream proteins. The complement of binding sites that comprise each virtual phosphorylation site can be found in the Supplementary information online. If a virtual phosphorylation site contains more than one actual site that is capable of recognizing a downstream protein (e.g. ErbB3 has multiple PI‐3K binding sites), the binding on‐rate of the downstream protein to the virtual site should be multiplied by the number of sites it can potentially bind, a factor we term the on‐rate multiplier (Supplementary information online).
We assume that tyrosine‐phosphorylated Shc, RasGAP, and Gab1 are concentrated in a relatively thin layer near the membrane, where their ErbB kinases are localized (Brown and Kholodenko, 1999). Once these phosphorylated proteins dissociate from the membrane, they are rapidly dephosphorylated; therefore, we neglect the concentration of these phosphorylated forms in the cytoplasm. We also do not consider the double membrane recruitment of Grb2 or Gab1. For example, we neglect that Gab1 membrane‐bound through its PH domain can simultaneously bind via its C‐terminal SH3 domain to membrane‐bound Grb2 (although we do consider that this same Gab1 can recruit free Grb2 through its PRD). As only one of these binding events is sufficient for membrane recruitment, accounting for both membrane localizing binding events simultaneously should have a small effect on downstream signaling.
Membrane recruitment is mediated by specific binding reactions between multiple network entities, and it results in the formation of several large multi‐protein complexes. Accounting for the microscopic details of all these large membrane complexes creates a combinatorial explosion of potential species. To reduce this complexity, we assume that for propagating signals, the route of membrane recruitment does not affect the action of the recruited protein. For example, we assume that it is not important whether Grb2 is bound to different receptor dimers or to tyrosine‐phosphorylated Shc (the route of membrane recruitment), all that matters is that Grb2 is at the membrane such that it can recruit downstream proteins. Mathematically, we assume (in this example) the rate law describing how downstream proteins bind to Grb2 is identical regardless of how Grb2 is recruited to the membrane. This assumption allows us to reduce model complexity by defining lumped, membrane‐localized states (denoted by Σ) that are sums over plasma membrane‐localized proteins that have identical downstream signaling action. The definitions of membrane‐localized states can be found in Table I, and lumping is depicted by double solid‐head arrows in Figure 2.
To use this Σ‐approach, we must modify traditional methods of modeling. First, ‘seeding’ species that undergo subsequent binding or modification (e.g. EijG) represent a mosaic of different membrane‐bound states. Consequently, these seeding species cannot be used in moiety conservation calculations for adapter and scaffold proteins. Moiety conservation is determined by adding all membrane‐localized Σ‐states and free states (Supplementary information online). Second, we must change the driving force for dissociation reactions of seeding species. For example, dissociation of E11G cannot be driven by the concentration of E11G, because other things could have bound to the G. Therefore, we multiply the dissociation driving force concentration by a fraction that represents the amount of the particular membrane‐bound species that has not been subsequently modified (Supplementary information online).
Downstream signaling and feedback.
The main effects of membrane localization of recruited proteins are dephosphorylation of phosphotyrosines by PTP‐1B, production of PIP3 by PI‐3K, which leads to Akt activation, and conversion of Ras‐GDP to Ras‐GTP by SOS, which leads to Raf, MEK, and ERK activation. Ras‐GTP is converted back to Ras‐GDP by both RasGAP and phosphorylated RasGAP, and PIP3 is converted back to PIP2 by PTEN. To reduce model complexity, when appropriate, we describe complex yet sequential multistep processes, such as the activation of Raf by Ras and Akt by PIP3, as a single, semimechanistic step. Because the condensed multistep processes are sequential, these simplifications allow the reduced model to retain the original network topology. However, these simplifications necessarily require us to describe the simplified steps empirically.
Whether a reduced, one‐step description is appropriate or not depends on if the activation mechanism is processive or distributed. In processive mechanisms, the substrate does not dissociate and rebind the enzyme between intermediate activation steps, and therefore a reduced, one‐step mechanism using a Michaelis–Menten rate law is adequate in these cases. In distributed mechanisms, the substrate dissociates and rebinds the enzyme between catalysis steps, and this can lead to potential ultrasensitivity and bistability (Markevich et al, 2004). As it is known that the activation of ERK by MEK follows a distributed mechanism (Zhao and Zhang, 2001), we use a full, mass action description for ERK activation. In the absence of evidence for a distributed mechanism, we use a simplified, one‐step description.
Although it is clear that internalized 1‐1 homodimers can initiate signaling (Wiley, 2003), we do not explicitly distinguish between signaling from the plasma membrane and internalized 1‐1 dimers. Thus, the 1‐1 homodimer states in the model represent a sum over internalized and plasma membrane receptors. However, there is evidence that PI‐3K cannot signal from endosomes due to lack of PIP2 (Haugh and Meyer, 2002). To account for this phenomenon yet keep from significantly increasing model complexity, we define an empirical, time‐dependent parameter that represents the endosomal fraction of ErbB1 receptors that have indirectly recruited PI‐3K. This parameter is represented by a general first order system, and the details are in the Supplementary information online. We then exclude this fraction of PI‐3K from converting PIP2 to PIP3.
Activated ERK mediates multiple modes of negative feedback. Activated ERK can inactivate SOS (Dong et al, 1996) and Gab1 (Lehr et al, 2004) by direct phosphorylation. Additionally, many literature studies show that activated ERK phosphorylates the ErbB1 receptor by direct phosphorylation of T669 (e.g. Northwood et al, 1991; Sanghera et al, 1992), and that this can lead to decreased receptor internalization and substrate specificity (Heisermann et al, 1990) and can mediate decreased kinase activity (Takishima et al, 1988). Therefore, we allow for the possibility that ERK can phosphorylate ErbB1, causing protection from degradation and decreased kinase activity. Interestingly, this threonine site and the motif recognized by ERK (PLTP) is conserved on ErbB2 and ErbB4, but not on ErbB3. Therefore, we also allow for the possibility of ERK phosphorylation of ErbB2 and ErbB4.
Dynamic and dose response of the ErbB signaling network to EGF and HRG
Investigation of any dynamic system begins with input perturbations followed by monitoring of system variables. Therefore, we made different magnitude step changes in the concentrations of EGF and HRG and observed the dynamic response of multiple components of the ErbB signaling network. We constrain our study to the first 30 min of signaling, as past this point immediate early gene responses begin to become important in determining the network activation state. Until such responses are mechanistically well understood, it will be difficult to build a model that can describe signaling past 30 min. However, even if the immediate early gene responses are well understood, it will still be a technically difficult task to model this next layer of complexity.
Experimental results along with model simulations are shown in Figure 4. Overall, the model simulations show excellent quantitative and qualitative agreement with the experimental data, although the model slightly overpredicts 30‐min ERK activity. We think this may be due to not incorporating the immediate early gene expression response of MKP1, which is an ERK phosphatase known to be upregulated in response to EGF within 30 min (Brondello et al, 1997). Adding such a mechanism is out of the scope of this paper and we are presently working on such a model.
Figure 4A and B show results from simultaneous stimulation with both EGF and HRG. At constant HRG stimulation, increasing EGF slightly increases the response of both ERK and Akt. In contrast, at constant EGF stimulation, increases in HRG cause a large increase in ERK and Akt activation, and notably, make the signaling more sustained. The effect that HRG has on EGF signaling is more pronounced at low EGF doses. These data collectively show that HRG acts as a dominant ligand, that is, when both HRG and EGF are present, the downstream signaling resembles that of HRG.
Figure 4C shows the results when the EGF and HRG responses are compared directly. We see that globally HRG stimulates higher peak and more sustained activation of Shc, MEK, ERK and Akt. However, both model simulations and experimental data show that EGF‐induced 5‐min ERK activation is approximately 80% of the HRG‐induced value, whereas EGF‐induced MEK activation is only approximately 50% of the HRG‐induced value. This suggests that HRG stimulates excess MEK activation above what is needed for peak ERK activation, and leads us to hypothesize that HRG‐induced ERK activity should be less sensitive to a MEK inhibitor than EGF‐induced ERK activity. We test this hypothesis experimentally, and the results are presented in a subsequent section of the current study.
Model parameter sensitivity analysis generates novel hypotheses as to the dominant mechanisms for ERK and Akt activation
To begin to understand the behavior of this complex system, we performed linear sensitivity analysis by making a 1% change in each model parameter and looking at the fractional effect on each observable. (This analysis is analogous to calculating the control coefficients and response coefficients in metabolic control analysis; e.g. Hornberg et al, 2005). The results from this analysis contain a wealth of information about how the ErbB signaling network might be functioning in MCF‐7 cells, and the full parameter sensitivity matrix can be found in the Supplementary information online. Here, we focus our efforts on investigating the mechanisms that may be controlling the peak (5‐min) and long‐time (30‐min) activation of ERK and Akt at different dose combinations, and these results are shown in Table II. Although Table II contains many results that are to be expected (e.g. parameters directly affecting PIP2/PIP3 conversion control Akt activity and parameters directly affecting Ras, Raf, and MEK activation control ERK activity), we find that there are several non‐obvious results, which we elaborate on below.
ERK activation is robust to parameter perturbation at high ligand doses, whereas Akt activation is not.
This statement is evident by observing the magnitude of parameter perturbation effects at 10 nM EGF, 10 nM HRG, and 10 nM EGF with 10 nM HRG (Table II). For ERK, at both short and long time points, the majority of the top positive and negative sensitivity coefficients are well below 0.01 (10 nM EGF‐stimulated 30 min ERK activation to a lesser extent than the others). This indicates that the parameters with the greatest effect on ERK activation can be changed by significant amounts without having a large impact. On the other hand, for Akt, the largest sensitivities are well above 0.01, indicating that small changes in these parameters can cause large changes in Akt activity. This difference in robustness to parameter variation may be due to the involvement of ERK in multiple negative feedback loops, and the absence of such feedback from Akt. It is well known from control theory that negative feedback endows a system with robustness to disturbances (Ogunnaike and Ray, 1994; Freeman, 2000).
An alternative way to interpret robustness properties of a particular quantity from sensitivity analysis is to evaluate statistics of the sensitivity coefficient ‘population’ (Table III). If the mean of the population is close to or essentially zero (in this case ∣〈S〉∣≪0.01), then a small standard deviation is indicative of robustness, as most sensitivity coefficients will then be close to zero. Again if the mean of the population is essentially zero, then a high kurtosis (the 4th moment; a measure of the ‘peakedness’ of the distribution around the mean) is also indicative of robustness. High kurtosis indicates that the majority of the sensitivities are clustered tightly around the mean, with more of the variance due to large, infrequent deviations rather than moderate, frequent deviations. Interestingly, a higher kurtosis gives a higher probability that a sensitivity with extreme deviation from the mean will occur, indicating that although a particular quantity may be robust to many perturbations, it may be very fragile to a few. This interpretation of kurtosis is also in line with the current thinking that there is an inevitable tradeoff between robustness and fragility (Kitano, 2004). We see from Table III that these statistical interpretations corroborate the arguments presented above that ERK activation is robust at high ligand doses, whereas Akt is not, and additionally represent a simple and quick way to evaluate robustness.
PI‐3K is a dominant factor for post‐ but not pre‐peak EGF‐induced ERK activation.
For both low and high doses of EGF, Table II shows that PI‐3K is a dominant‐positive factor for 30 min ERK activation but not for 5 min ERK activation. This implies that EGF‐induced ERK activity is controlled by different cellular mechanisms at different times; the descending portion (post‐peak) but not the rising portion (pre‐peak) of the response curve is highly dependent on PI‐3K activity. This phenomenon is due to degradation of active 1‐1 homodimers. At short times, there are still many 1‐1 homodimers that can signal to ERK, however, as time progresses these homodimers are degraded, and therefore signaling to ERK must rely on alternative mechanisms at longer times. In this case, an alternative mechanism is the PI‐3K‐PIP3‐Gab1 pathway. We note that this conclusion is not sensitive to the time‐dependent, non‐mechanistic parameter that characterizes the amount of ErbB1 in endosomes incapable of signaling through PI‐3K.
Modeling analysis of the multiple negative feedback loops
There are multiple negative feedback loops incorporated into the model, and with the exception of ligand‐induced degradation, it is not obvious what distinct roles, if any, these feedback loops may play in regulating ligand‐dependent signaling. To investigate the potential roles of these various negative feedback loops, we inhibited these feedbacks in silico and observed the predicted ERK and Akt activation at different ligand doses (Figure 5). As negative feedback loops are being inhibited, we expected that ERK and Akt activity should always increase. However, Figure 5 shows that this is not always the case. Most notably, ERK negative feedback to receptors (Figure 5B) positively affects EGF‐induced peak ERK and Akt activity. Further simulations suggested that this is because ERK inhibits ErbB2 less than ErbB1, manifested as decreased RasGAP membrane recruitment mediated by a shift toward more 1‐2 heterodimers rather than 1‐1 homodimers.
While simulations predict that PTP1‐B is the dominant‐negative feedback for Akt signaling (Figure 5A) with ERK‐Gab feedback (Figure 5D) and ligand‐induced degradation (Figure 5F) playing a smaller role, there is no dominant feedback controlling ERK signaling; it differs based on time and dose combination. Long‐term, EGF‐induced ERK signaling is mainly controlled by ERK receptor feedback (Figure 5B) and ligand‐induced degradation (Figure 5F), short‐ and long‐term HRG‐induced signaling is mainly controlled by RasGAP phosphorylation (Figure 5E) and to a lesser extent the ERK‐SOS feedback (Figure 5C) and ERK‐Gab feedback (Figure 5D). Such distributed control of ERK activity makes it a very versatile signaling entity for the cell, being finely tunable in many different ways, perhaps explaining, in part, why ERK is important in so many different signaling systems.
Model predictions, validation, and limitations
ErbB2 overexpression sustains EGF signaling.
The above model analyses allowed us to gain general understanding of the ligand‐dependent control of ERK and Akt activity. While this general understanding is valuable by itself, understanding how known, cancer‐correlated network abnormalities affect signaling is of primary and immediate interest. We used the model to investigate how such a network abnormality, ErbB2 overexpression, which occurs in approximately 25% of breast cancer, would affect the dynamics of ERK and Akt activation, and the results are shown in Figure 6A. ErbB2 overexpression is not predicted to cause significant changes in HRG signaling, whereas ErbB2 overexpression leads to more sustained EGF signaling, while leaving peak ERK activation unaffected and actually decreasing peak Akt activation. This model prediction qualitatively and quantitatively agrees with the data of Wolf‐Yadlin et al (2006) and Kumar et al (2007), which shows that at 10 and 30 min after EGF stimulation, ERK activation is between 1.15‐ and 2‐fold higher for HMECs overexpressing ErbB2. The agreement between our model prediction and an independent data set gives a strong point of model validation.
As ErbB2 overexpression shifts the receptor dimer distribution toward more 1‐2 heterodimers rather than 1‐1 homodimers and only 1‐1 homodimers undergo ligand‐induced degradation, the more sustained EGF‐signaling is expected. However, the predicted decrease in peak Akt signaling is not expected, as any increase in receptor abundance should lead to increased signaling. Looking at the kinetic scheme and model parameters, we hypothesized that this phenomenon might be due to increased recruitment of PTP1‐B by 1‐2 heterodimers, as 1‐2 heterodimers had a slightly higher predicted affinity for PTP1‐B than 1–1 homodimers. We therefore performed simulations where the rate constants for reaction 74 (PTP1‐B binding to 1‐2 heterodimers) were set equal to the rate constants for reaction 73 (PTP1‐B binding to 1‐1 homodimers), and the results are shown in Figure 6B. We found that under these conditions ErbB2 overexpression then increases peak Akt activation.
HRG‐stimulated ERK activation is more robust than EGF‐stimulated ERK activation to the MEK inhibitor U0126.
The experimental and simulation data in Figure 4C show that both EGF and HRG induce similar peak ERK activation, but the magnitude of peak MEK activation is much smaller for EGF. This suggests that HRG may stimulate more MEK activity than is necessary to induce the observed peak ERK activation, and led us to hypothesize that HRG‐stimulated ERK activation would be less sensitive to a MEK inhibitor (U0126) than would EGF‐stimulated ERK activation. We investigated this hypothesis both computationally and experimentally, and the results are shown in Figures 7A and B. As can be seen, both the experimental data and model simulations show that EGF‐stimulated ERK activation is indeed more sensitive to U0126 than HRG‐stimulated ERK activation. Although the model slightly overpredicts the effect of U0126 on HRG‐stimulated ERK activation, the model simulations have excellent quantitative agreement with the EGF data.
Ligand‐dependent effects of the PI‐3K inhibitor wortmannin.
To test further the predictive powers of our model, we experimentally measured the effect of a PI‐3K inhibitor, wortmannin, on EGF‐ and HRG‐induced ERK activation over a 30‐min time course (Figure 8A), and we compared these data to the model predictions (Figure 8B). The model correctly predicts several features of the experimental data. First, we observe that wortmannin has significant effect on post‐ but not pre‐peak EGF‐induced ERK activity. Importantly, this corroborates our previous prediction based on sensitivity analysis. Second, for all time points less than 30 min, both the model and experiments show that wortmannin has a small effect on HRG‐induced ERK activity. Third, experiments and model simulations show that wortmannin affects the magnitude of EGF‐induced ERK signaling at times greater than 5 min, but does not change the qualitative shape of the response.
While there are several points of corroboration between model predictions and experimental data in Figure 8, there are also points of disagreement. The model underestimates the effect of wortmannin on EGF‐induced signaling at long times. Also, the experimental data show that wortmannin causes HRG‐induced signaling to become slightly transient at 30 min, but the model does not show this. As mentioned and explained above in ‘Dynamic and Dose Response of the ErbB Signaling Network to EGF and HRG’, the model over estimates the magnitude of EGF‐induced signaling past 10 min.
While these discrepancies show that the model cannot predict all the effects of wortmannin on ERK activation, they give us direction as to how the model may be improved in future studies. An attractive candidate for addition to the model, which may help to explain these wortmannin data, is the tyrosine phosphatase SHP2. SHP2 is recruited to Gab1 in response to ligand stimulation (Yart et al, 2001), but is actually a positive regulator of ERK signaling due to its selective dephosphorylation of RasGAP binding sites on receptors (Agazie and Hayman, 2003) and Gab1 (Montagner et al, 2005). Thus SHP2 represents an additional PI‐3K‐dependent pathway that positively affects ERK signaling.
In the present work, we developed a quantitative kinetic model that relates EGF and HRG stimulation of the ErbB receptors to ERK and Akt activation in MCF‐7 breast cancer cells. To our knowledge, this is first model to take into account all four ErbB receptors, simultaneous stimulation with two ligands, and both the ERK and Akt pathways. In the context of this study, we used this model to help gain mechanistic insight into ligand‐dependent responses of the ErbB signaling network.
When a mechanistic modeling approach is used to represent the ErbB signaling network, combinatorial complexity makes the modeling task infeasible. Rule‐based modeling tools that facilitate building combinatorially complex models have been recently developed and promise to be valuable for dealing with the complexity of signal transduction modeling (e.g. Blinov et al, 2006). However, until current parameter identification and bifurcation analysis techniques are capable of handling extremely large models generated by rule‐based algorithms, we should take steps to reduce combinatorially complex models and make them mathematically tractable. As similar to previously published models on ErbB signaling (e.g. Kholodenko et al, 1999), we represented the entire complement of phosphorylated cytoplasmic receptor tyrosines as a ‘virtual phosphorylation site’. We showed a subtle yet important consequence of this assumption: when a protein binds to multiple locations on the receptor tail, the on‐rate should be multiplied by the number of locations. We introduce the concept of membrane‐localized ‘Σ‐states’, where we assume that the action of membrane recruited proteins can be decoupled from its path to the membrane. Although all these assumptions impose approximations, they also make the modeling exercise computationally tractable. As the model reflects the measured system dynamics reasonably well and gives insight into the differences in ligand‐dependent signaling, the assumptions seem reasonable for the analysis carried out.
To begin to understand the behavior of this complex system, and specifically how peak and long‐time ERK and Akt activities are regulated, we employed sensitivity analysis and in silico inhibition studies. In general, such modeling analyses are indispensable for generating hypotheses as to how quantities of interest are controlled in complex dynamic systems. Naturally, these analyses revealed many obvious modes of regulation, for instance, control of Akt activation by PIP2/PIP3 reaction parameters and control of ERK activation by Ras, Raf, and MEK reaction parameters. However, these modeling analyses also revealed non‐obvious modes of potential regulation: (i) ERK activation is robust to parameter perturbation at high ligand doses, whereas Akt activation is not and (ii) PI‐3K is a dominant factor for long‐term but not short‐term EGF‐induced ERK activation. Analysis of the multiple negative feedback loops revealed specialized roles for each feedback in regulating ERK and Akt activity at different time points and dose combinations. By combining logical reasoning with simulation analysis, these observations led to many new insights as to how the ErbB signaling network functions in MCF‐7 cells, such as the signaling specificity of different ErbB dimers and robustness of signaling responses. Investigating all of these hypotheses experimentally is outside the scope of a single study. However, over time the combined results of many studies can be used to test these hypotheses, and when conflicts arise, revise our model and understanding of the ErbB signaling network.
We also used the model to generate hypotheses about a specific topic of known relevance: the effect of ErbB2 overexpression on ErbB network signaling. Our analysis revealed that while ErbB2 overexpression should have a small effect on HRG‐induced signaling, it has causes EGF‐induced signaling to become more sustained. The prediction that ErbB2 overexpression will sustain EGF‐induced ERK signaling is corroborated by studies of ErbB2 overexpression in HMECs, which showed that at 10 and 30 min after EGF stimulation, ERK activation is between 1.15‐ and 2‐fold higher for cells overexpressing ErbB2 (Wolf‐Yadlin et al, 2006; Kumar et al, 2007). The current hypothesis to explain the sustainment of EGF‐signaling is that ErbB2 overexpression causes a higher proportion of 1‐2 heterodimers rather than 1‐1 homodimers to be activated upon ligand stimulation, and as only 1‐1 homodimers undergo ligand‐induced internalization, this overexpression sustains signaling. While Sorkin et al (1993) and Baulida et al (1996) were among the first to show the preferential retention of ErbB2‐containing heterodimers on the plasma membrane, recent experimental studies also support this hypothesis. Haslekas et al (2005) showed that ErbB2 inhibits internalization of EGF bound to ErbB1 by driving heterodimerization. Lynch et al (2004) reported that internalization defective ErbB1 mutants commonly found in ErbB1‐kinase‐inhibitor‐responsive non‐small cell lung cancers display sustained signaling in response to EGF. Hendriks et al (2005) showed that 184A1 HMECs with high ErbB2 expression had more sustained EGF‐induced ERK activity compared to ERK responses from the same cell line with low ErbB2 expression. Furthermore, their modeling results imply that this difference may be due to ErbB2‐induced alterations in ErbB1 receptor trafficking. Overall, our model predictions and the corresponding literature lend insight into potential reasons why ErbB2‐overexpressing or ErbB1 internalization‐defective mutants are found in cancer.
We also investigated a hypothesis resulting from preliminary examination of the experimental data and model simulation results: the differential sensitivity of EGF‐ and HRG‐stimulated ERK activity to the MEK inhibitor U0126. The model predicted that HRG‐stimulated ERK activity would be less sensitive to U0126 than the EGF‐stimulated activity, and this result was confirmed experimentally. This result is also consistent with the experimental studies of Thottassery et al (2004), who found that HRG but not EGF reduced the negative effect of U0126 on proliferation and on ERK activation in MCF‐7 cells. Additionally, the model had excellent quantitative agreement for EGF signaling, but slightly overpredicted the effect of U0126 on HRG signaling. It may be possible that mechanisms not incorporated into the model, such as the ERK cascade scaffolds KSR and MP1 or B‐Raf, which most likely decrease the sensitivity of the ERK cascade to inhibition, may be able to account for the HRG discrepencies.
While our model is able to correctly predict many features of ligand‐dependent ErbB signaling, we also found limitations to predict all the effects of wortmannin on ERK signaling. The current model correctly predicted that post‐peak but not pre‐peak EGF‐induced ERK signaling should be sensitive to a PI‐3K inhibitor, but underestimated the post‐peak effects of wortmannin on EGF‐ and HRG‐induced ERK signaling. This lack of agreement between the model predictions and experimental data is not entirely unexpected, as all models have limitations. These limitations suggest where the model needs to be refined and improved in future studies. The incorporation of the tyrosine phosphatase SHP2 is a logical next step to help explain these wortmannin effects, as SHP2 is a PI‐3K‐dependent positive regulator of ERK signaling. Incorporation of SHP2 is a significant modeling undertaking, however, as SHP2 selectively dephosphorylates RasGAP‐binding sites on Gab1 and ErbB receptors and its addition will require description of individual phosphorylation sites, as opposed to virtual phosphorylation sites. Such a description will require further development of novel reduction methods to keep the model tractable.
Studies similar to the current work hold the potential to complement the field of targeted cancer treatment. Currently, many ErbB‐targeted pharmaceuticals are clinically used to treat cancer, including the small molecule ErbB1 kinase inhibitors erlotinib and gefitinib (Harari and Huang, 2006), and a monoclonal ErbB2 antibody, trastuzumab (Adams and Weiner, 2005). Although correlation analyses have provided insight as to when these pharmaceuticals will be effective, there is still much work to be done to determine when to use these pharmaceuticals, and in what combinations. A computational model of the ErbB signaling systems in the cancer type of interest can help predict factors that could guide the choice of when to use a particular targeted pharmaceutical, and in what combinations. As a simple example, we saw in the current analysis that for an ErbB2‐overexpressing cancer, if ErbB1 homodimers do not undergo ligand‐induced degradation, then inhibiting ErbB2 with a drug such as trastuzumab may have minimal effect on EGF‐induced signaling. More generally, the results of sensitivity analysis can be used to conjecture the critical determinants of signaling behavior, and thus give potential drug targets. A multivariate, global sensitivity analysis could give insight into the drug combinations that may work synergistically.
Materials and methods
Cell culture and western blots
The MCF‐7 human breast cancer cell line was obtained from American Type Culture Collection (ATCC) and maintained in DMEM (Gibco BRL, Githersburg, MD) supplemented with 10% fetal bovine serum. Before growth hormone treatment, the cells were synchronized by serum‐starvation for 16–24 h, and then EGF (PeproTech House, London, England) or HRG‐β176‐246 (R&D Systems Inc., Minneapolis, MN) was added. In experiments that involved U0126 (Calbiochem, San Diego, CA) treatment, the inhibitor was added 20 min before growth hormone treatment. Cells were incubated with the growth hormone in the presence or absence of kinase inhibitors for the indicated period of time, washed three times with phosphate‐buffered saline (PBS), and lysed with Bio‐Plex lysis buffer (Bio‐Rad laboratories, Hercules, CA).
Cell lysate was cleared by centrifugation, and the total protein concentration of the supernatant was determined using a protein assay reagent. Lysates were normalized to the same total protein concentration by dilution. Gel electrophoresis and western blotting were performed as described previously (Hatakeyama et al, 2003). Antibodies against doubly phosphorylated p44/42 ERK, ERK, phospho‐Akt (Ser473), Akt, phospho‐MEK1/2 (Ser217/221), and MEK were purchased from Cell Signaling Technology Inc. (Beverly, MA), anti‐phospho‐Shc (Tyr317) and anti‐Shc antibodies were purchased from Upstate Biotechnology (Lake Placid, NY), and all anti‐ErbB receptor antibodies and anti‐phosphotyrosine antibody (PY20) were purchased from Santa Cruz Biotechnology (Santa Cruz, CA). To detect ErbB receptor phosphorylation, the total cell lysate was immunoprecipitated with the corresponding ErbB antibodies, and immunoblotted with anti‐phosphotyrosine antibody (PY20, Santa Cruz Biotechnology). Protein band intensities were quantified using a densitometer (Fuji Film Corp., Japan). To analyze a large set of samples at a time, as required with the dual‐ligand stimulation experiments, we used the Bio‐Plex Suspension Assay System (Bio‐Rad Laboratories) to detect phosphorylated and non‐phosphorylated forms of ERK and Akt. We previously confirmed that the data obtained from the Bio‐Plex assay are reproducible and correlate well with western blot data in our cell system (unpublished data).
From the kinetic scheme (Figure 2), which describes the connectivity of the reaction network, a deterministic ODE model was derived using the law of mass action and saturating rate laws to describe the corresponding reactions. The rate of change of a species concentration with time (time derivative) is calculated by summing all reaction rates that produce this species and subtracting all reaction rates that consume this species. When a reaction involves species whose concentrations are defined in compartments different from where the reaction is taking place, these species concentrations are rescaled to reflect the reaction compartment volume, and the reaction rate is calculated based on the rescaled species concentrations. The species concentrations are rescaled back to the reference compartment volume for calculating the rate of change. The cytoplasm is chosen as the reference compartment in this study. Details about the rate laws and parameter values can be found in the Supplementary information online.
Model fitting & simulation
Fitting & simulation was carried out with MATLAB 7.0 Serivce Pack 2 software (The Mathworks; Natick, MA) on AMD 64 Dual Core 2.0 GHz processor computers running CENTOS 4.4 linux (www.centos.org). Differential equations were integrated using the function ode15s, which is a variable order solver, based on the numerical differentiation formulas (NDFs) and is designed for stiff systems. Local minimization was carried out with the function lsqnonlin, which is a subspace trust region method and is based on the interior‐reflective Newton method. Details of the parameter estimation can be found in the Supplementary information online.
Normalization of experimental data for direct comparison with model simulations
Several normalization steps were taken to make quantitative immunoblot data (relative concentrations) compatible with data from model simulations (absolute concentrations). First, all experimental data were divided by their respective control intensity. Second, all data collected from the same experiment/blot were divided by a normalization point, chosen to be the largest intensity point. This normalization point is chosen in favor of the zero point, because there is less uncertainty associated with more intense signals (higher signal to noise ratio). Lastly, the zero point from each blot is subtracted from every data point collected from that blot.
We thank the NIH for funding (Grants GM59570, AA07463, and R33HL088283).
Rate Equations, Species, and Other Additional Model Information
Sensitivity Analysis Results
Virtual Phosphorylation Site Compositions
On‐rate Multiplier Derivation
First Order Internalization and Degradation Derivation
Parameter Estimation Procedure
SBML File of the ErbB Signaling Model
This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.
- Copyright © 2007 EMBO and Nature Publishing Group