Abstract

Cellular signal transduction is governed by multiple feedback mechanisms to elicit robust cellular decisions. The specific contributions of individual feedback regulators, however, remain unclear. Based on extensive time‐resolved data sets in primary erythroid progenitor cells, we established a dynamic pathway model to dissect the roles of the two transcriptional negative feedback regulators of the suppressor of cytokine signaling (SOCS) family, CIS and SOCS3, in JAK2/STAT5 signaling. Facilitated by the model, we calculated the STAT5 response for experimentally unobservable Epo concentrations and provide a quantitative link between cell survival and the integrated response of STAT5 in the nucleus. Model predictions show that the two feedbacks CIS and SOCS3 are most effective at different ligand concentration ranges due to their distinct inhibitory mechanisms. This divided function of dual feedback regulation enables control of STAT5 responses for Epo concentrations that can vary 1000‐fold in vivo. Our modeling approach reveals dose‐dependent feedback control as key property to regulate STAT5‐mediated survival decisions over a broad range of ligand concentrations.

Visual Overview

Synopsis

Cells interpret information encoded by extracellular stimuli through the activation of intracellular signaling networks and translate this information into cellular decisions. A prime example for a system that is exposed to extremely variable ligand concentrations is the erythroid lineage. The key regulator Erythropoietin (Epo) facilitates continuous renewal of erythrocytes at low basal levels but also secures compensation in case of, e.g., blood loss through an up to 1000‐fold increase in hormone concentration. The Epo receptor (EpoR) is expressed on erythroid progenitor cells at the colony forming unit erythroid (CFU‐E) stage. Stimulation of these cells with Epo leads to rapid but transient activation of receptor and JAK2 phosphorylation followed by phosphorylation of the latent transcription factor STAT5. Although STAT5 is known to be an essential regulator of survival and differentiation of erythroid progenitor cells, a quantitative link between the dynamic properties of STAT5 signaling and survival decisions remained unknown. STAT5‐mediated responses in CFU‐E cells are modulated by multiple attenuation mechanisms that operate on different time scales. Fast‐acting mechanisms such as depletion of Epo by rapid receptor turnover and recruitment of the phosphatase SHP‐1 control the initial signal amplitude at the receptor level. Transcriptional feedback regulators such as suppressor of cytokine signaling (SOCS) family members CIS and SOCS3 operate at a slower time scale. Despite the ample knowledge of the individual components involved, only little is known about the specific contributions of these regulators in controlling dynamic properties of STAT5 in response to a broad range of input signals. Therefore, dynamic pathway modeling is required to understand the complex regulatory network of feedback regulators.

To address these questions, we established a dual negative feedback model of JAK2/STAT5 signaling in primary erythroid progenitor cells isolated from mouse fetal livers. We provide a large data set of JAK2/STAT5 signaling dynamics employing quantitative immunoblotting, mass spectrometry and quantitative RT–PCR measured under different perturbation conditions to calibrate our model (Figure 3). The structure of our model was constructed to comprise the minimal number of parameters necessary to explain the data. Thereby, we aimed at a model with fully identifiable parameters that are essential to obtain high predictive power. Parameter identifiability was analyzed by the profile likelihood approach. Applying this method, we could establish a dual negative feedback model of JAK2‐STAT5 signaling with structurally and in most cases practically identifiable parameters.

A major bottle‐neck in combining signal transduction events with cellular phenotypes is the discrepancy in the time scale and stimuli concentrations that are applied in the different experiments. The sensitivity of biochemical assays to determine phosphorylation events within minutes or hours after stimulation is usually lower than the threshold of sensitivity in assays to determine the physiological response after one or more days. Facilitated by the model, we were able to compute the integrated response of JAK2/STAT5 signaling components for experimentally unaddressable Epo concentrations. Our results demonstrate that the integrated response of pSTAT5 in the nucleus accurately correlates with the experimentally determined survival of CFU‐E cells. This provides a quantitative link of the dependency of primary CFU‐E cells on pSTAT5 activation dynamics. By correlation analysis, we could identify the early signaling phase (⩽1 h) of STAT5 to be the most predictive for the fraction of surviving cells, which was determined ∼24 h later. Thus, we hypothesize that as a general principle in apoptotic decisions, ligand concentrations translated into kinetic‐encoded information of early signaling events downstream of receptors can be predictive for survival decisions 24 h later.

After the first hour of stimulation, it is important to constrain signaling to a residual steady‐state level. Constitutive phosphorylation of the JAK2/STAT5 pathway has a crucial role in the onset of polycythemia vera (PV), a disease associated with Epo‐independent erythroid differentiation. The two identified transcriptional feedback proteins, CIS and SOCS3, are responsible for adjusting the phosphorylation level of STAT5 after 1 h of stimulation. Since the Epo input signal can vary over a broad range of ligand concentrations, we asked how CIS and SOCS3 can facilitate control of STAT5 long‐term phosphorylation levels over the entire physiological relevant hormone concentrations. By using model simulations, we revealed that the two feedbacks are most effective at different Epo concentration ranges. Predicted by our mathematical model, the major role of CIS in modulating STAT5 phosphorylation levels is at low, basal Epo concentrations, whereas SOCS3 is essential to control the STAT5 phosphorylation levels at high Epo doses (Figure 6). As a potential molecular mechanism of this dose‐dependent inhibitory effect, we could identify the quantity of pJAK2 relative to pEpoR that increases with higher Epo concentrations. Since SOCS3 can inhibit JAK2 directly via its KIR domain to attenuate downstream STAT5 activation, SOCS3 becomes more effective with the relative increase of JAK2 activation. Hence, CIS and SOCS3 act in a concerted manner to ensure tight regulation of STAT5 responses over the broad physiological range of Epo concentrations.

In summary, our mathematical approach provided new insights into the specific function of feedback regulation in STAT5‐mediated life or death decisions of primary erythroid cells. We dissected the roles of the transcriptionally induced proteins CIS and SOCS3 that operate as dual feedback with divided function thereby facilitating the control of STAT5 activation levels over the entire range of physiological Epo concentrations. The detailed understanding of the molecular processes and control distribution of Epo‐induced JAK/STAT signaling can be further applied to gain insights into alterations promoting malignant hematopoietic diseases.

We show that the amount of nuclear phosphorylated STAT5 integrated for 60 min post Epo stimulation directly correlates with the fraction of surviving cells 24 h later.

CIS and SOCS3 were identified as the most relevant transcriptional feedback regulators of JAK2/STAT5 signaling in primary erythroid progenitor cells. Applying the model, we revealed that CIS‐mediated inhibitory effects are most important at low ligand concentrations, whereas SOCS3 inhibition is more effective at high ligand doses.

The distinct modes of inhibition of CIS and SOCS3 at various Epo concentrations provide a strategy for achieving control of JAK2/STAT5 signaling over the entire range of physiological Epo concentrations.

Introduction

Cells interpret information encoded by extracellular stimuli through the activation of intracellular signaling networks and link this to cellular decisions (Kholodenko et al, 2010). Feed‐forward signaling events are counteracted by negative feedback mechanisms facilitating adaptation, stability and desensitization of input signals. Transcriptional feedback regulators that are induced as early genes have been identified as major constituent of signaling pathways such as DUSP family members for the MAPK cascade, Smad7 and Ski family members for the TGFβ/Smad pathway and suppressor of cytokine signaling (SOCS) proteins for JAK/STAT signaling (Legewie et al, 2008). Importantly, these transcriptional feedback regulators are downregulated in multiple carcinomas, implying a crucial role in constraining growth factor and cytokine‐induced signaling (Freeman, 2000; Amit et al, 2007). Despite the ample knowledge of the individual components involved, only little is known about the specific contributions of these regulators in controlling dose‐dependent dynamic properties of signaling networks in response to a broad range of input signals.

A prime example for a system that is exposed to extremely variable ligand concentrations is the erythroid lineage. Key regulator of erythropoiesis is the hormone erythropoietin (Epo) that facilitates continuous renewal of short‐lived erythrocytes at low basal levels of ∼15 mU/ml Epo, but also secures compensation of blood loss through an up to 1000‐fold increase in hormone concentration in acute hypoxic conditions (Jelkmann, 2004). Epo binds to the hematopoietic cytokine receptor, the Epo receptor (EpoR), that is primarily present on erythroid progenitor cells, but recently also has been uncovered on tumor cells (Lappin et al, 2002; Lai et al, 2005; Liang et al, 2010). In the erythroid lineage signaling through the EpoR is essential for survival, proliferation and differentiation of erythroid progenitor cells at the colony forming unit erythroid (CFU‐E) stage. Stimulation of these cells with Epo leads to rapid but transient activation of receptor and JAK2 phosphorylation followed by phosphorylation of the latent transcription factor STAT5. The role of STAT5a/b in erythropoiesis has been controversially discussed, but evidence based on the analysis of STAT5ΔN/ΔN mice that express partially functional STAT5a/b proteins and mice completely deficient of STAT5a/b (Socolovsky et al, 1999; Cui et al, 2004; Yao et al, 2006; Kerenyi et al, 2008) suggests that STAT5 is a crucial regulator of survival and differentiation in erythroid progenitor cells. Initial studies showed that erythroid cells expressing STAT5ΔN/ΔN, the hypomorphic STAT5a/b proteins, exhibit elevated rates of apoptosis due to a failure of Bcl‐xL upregulation (Socolovsky et al, 2001). Subsequently, generated STAT5a/b−/− mice displayed fatal defects in all lymphoid lineages, suffered from microcytic anemia due to enhanced apoptosis and died prenatally (Zhu et al, 2008). Recently, it was shown that STAT5 can enable erythropoiesis even in the absence of EpoR and JAK2 underscoring the essential role of STAT5 in regulating apoptosis and differentiation in erythropoiesis (Grebien et al, 2008). However, a quantitative link between dynamic properties of STAT5 signaling and survival decisions remained to be established.

The STAT5‐mediated responses in CFU‐E cells are modulated by multiple attenuation mechanisms that operate on different time scales. Fast‐acting mechanisms such as depletion of Epo by rapid receptor turnover (Becker et al, 2010) and recruitment of the phosphatase SHP‐1 (Klingmüller et al, 1995) control the initial signal amplitude at the receptor level. Transcriptional feedback regulators such as SOCS family members operate at a slower time scale and the precise impact on the kinetics of STAT5 signaling has not yet been determined. Among the eight known SOCS family members, four have been proposed in studies primarily performed with transformed (erythro‐) leukemic cell lines as putative transcriptional negative regulators involved in attenuation of EpoR‐JAK2/STAT5 signaling. They include cytokine‐induced SH2‐domain containing protein (CIS), SOCS1, SOCS2 and SOCS3 (Yoshimura et al, 1995; Jegalian and Wu, 2002; Sarna et al, 2003). Furthermore, overexpression of CIS in erythroid progenitor cells has been shown to inhibit STAT5‐mediated survival signals (Ketteler et al, 2003). Although the inhibitory mechanisms of these proteins are well characterized and are partially distinct (Endo et al, 1997; Matsumoto et al, 1997; Sasaki et al, 2000), a major open question was whether they have redundant or specific functions in restraining STAT5 signaling in erythroid progenitor cells.

Recently, several kinetic models of JAK/STAT signaling pathways have been established. The majority of these studies rely on previously published data and simulations focusing on control mechanisms of the IFN activated JAK1/JAK2/STAT1 pathway and feedback by SOCS1 (Zi et al, 2005; Soebiyanto et al, 2007) or robustness in the JAK/STAT1 pathway (Shudo et al, 2007). Up to now only few studies combine quantitative data with mathematical models. Examples are a core model of JAK2/STAT5 signaling investigating nucleo‐cytoplasmic shuttling (Swameye et al, 2003) and an IFNα‐induced model investigating a positive feedback loop (Maiwald et al, 2010). Understanding how transcriptional negative feedback proteins specifically regulate JAK/STAT pathway dynamics and the extent of cellular responses requires mathematical models that combine extensive time‐resolved biochemical and phenotypic data.

In this study, we establish a transcriptional feedback model of JAK2/STAT5 signaling in primary erythroid progenitor cells and demonstrate that the STAT5 response integrated over time is quantitatively linked to survival decisions. Applying the model, we dissect the divided function of CIS and SOCS3 that enables to confine STAT5 phosphorylation levels at distinct ligand concentration levels. Our results imply that dual feedback regulation ensures stringent and fine‐tuned control of cell survival decisions over the entire broad spectrum of physiologically relevant Epo concentrations.

Results

To examine the impact of transcriptional feedback regulation on Epo‐induced JAK2‐STAT5 signaling, we used the generic inhibitor of transcription actinomycin D and monitored the time‐dependent activation of the signaling network in primary erythroid progenitors at the CFU‐E stage. To confirm the effect of the actinomycin D treatment, the expression of the known Epo‐induced negative regulator CIS was monitored exemplarily by quantitative immunoblotting (Figure 1A and B). Stimulation with Epo during 3 h of observation leads to transient phosphorylation of EpoR, JAK2 and STAT5 with rapid kinetics, followed by a decrease to a residual steady‐state level. A comparison of time course data in the presence and absence of inhibitor demonstrated that the steady‐state phosphorylation level of STAT5 is elevated in actinomycin D‐treated cells, indicating that attenuation of this species involves de novo gene transcription. The phosphorylation kinetics of the upstream signaling components, EpoR and JAK2, however, showed no major effects compared with the profile of the untreated control. These results imply that transcriptional feedback regulators induced as immediate early genes have a central role in modulating the phosphorylation kinetics of STAT5.

Phosphorylation profile of Epo‐induced EpoR‐JAK2‐STAT5 signaling in the presence and absence of actinomycin D. (A) Time course experiment of STAT5 phosphorylation and CIS expression. Primary CFU‐E cells were starved and stimulated with 5 U/ml Epo after 10 min pre‐treatment with actinomycin D or vehicle (0.1% dimethylsulfoxide) alone. Cellular lysates were subjected to immunoprecipitation and immunoblotting with STAT5, CIS and anti‐phosphotyrosine antibodies. (B) Cellular lysates from primary CFU‐E cells treated as described were subjected to immunoprecipitation and quantitative immunoblotting with EpoR, JAK2, STAT5, CIS and anti‐phosphotyrosine antibodies. Splines (dashed lines) are displayed for visual guidance. Enhanced steady‐state phosphorylation level of STAT5 is highlighted (gray). (For immunoblots, see Supplementary Figure S2.) (C) Expression profiling of Epo‐induced JAK‐STAT signaling regulators. Primary CFU‐E cells were starved and stimulated with 0.5 U/ml Epo. RNA was subjected to microarray analysis and profiles of known SOCS genes were compiled. Log2‐fold change of mRNA levels was calculated relative to mean gene expression at time point t=0 h. Error bars were estimated from the gene expression variation of the biological duplicates at 0 h. Only CIS and SOCS3 show a significant regulation over time with respect to control (P‐values <0.0001, Student's t‐test). (For controls, see Supplementary Figure S1.) Source data is available for this figure at www.nature.com/msb.

To identify the relevant transcriptional feedback regulators that are involved in attenuation of Epo‐induced JAK2‐STAT5 signaling in addition to CIS, we performed genome‐wide expression profiling. Murine CFU‐E cells were stimulated with Epo and the time‐resolved mRNA expression profiles were analyzed to identify significantly upregulated negative feedback regulators of the JAK‐STAT signaling pathway (Figure 1C; Supplementary Figure S1). The results revealed that among the eight known SOCS proteins exclusively CIS and SOCS3 showed rapid and permanent upregulation at the mRNA level in particular within the first hour, which is consistent with the presence of STAT5 consensus binding sites in the promoters of these two genes (Matsumoto et al, 1997; Emanuelli et al, 2000). In contrast to studies in other cells that reported upregulation of additional SOCS transcripts upon Epo stimulation (Jegalian and Wu, 2002), neither SOCS1 nor SOCS2 mRNA induction was detected in the first hours after stimulation in CFU‐E cells. This suggests that other SOCS members are either primarily expressed in transformed cell lines or as hypothesized by Sarna et al (2003) may be induced at other stages during erythroid maturation. The slight upregulation of SOCS2 mRNA after ∼20 h supports the latter hypothesis.

Additional negative regulators of Epo‐induced JAK‐STAT signaling such as PIAS proteins (PIAS1‐4) and hematopoietic protein tyrosine phosphatases (PTPs) usually do not show stimulation‐dependent transcription but are constitutively expressed. To confirm this in primary erythroid progenitor cells, we compiled expression profiles of the PIAS1‐4 proteins and the hematopoietic phosphatases SHP‐1, SHP‐2, PTP1B and PTPRC (CD45) (Supplementary Figure S1). In line with previous reports, stimulation‐dependent induction of these genes was not detected. Hence, CIS and SOCS3 are the two prime candidates that can act as dual transcriptional feedback regulators of STAT5 phosphorylation levels in erythroid progenitor cells.

To dissect the specific roles of CIS and SOCS3 in controlling the JAK2‐STAT5 signaling behavior, we implemented an ordinary differential equation (ODE) model of the Epo‐induced JAK2/STAT5 pathway. Our model is divided into three submodules describing events at the plasma membrane, in the cytoplasm and in the nucleus (Figure 2). As model input, we assumed a constant Epo concentration, because the depletion of Epo in primary CFU‐E cells is marginal at the high concentrations used in this study. To explain the events at the plasma membrane, EpoR and the tyrosine kinase JAK2 were modeled as a complex that is present in different activation states corresponding to the different phosphorylated species (see Supplementary Material Section 2.2). The EpoRpJAK2 variable also describes phosphorylated JAK2 that is dissociated from the EpoR. Ligand‐dependent attenuation, for instance by receptor internalization and dephosphorylation, was summarized by including deactivation of JAK2 and EpoR by the tyrosine phosphatase SHP‐1 as described in a previous model of Epo‐induced MAPK signaling in CFU‐E cells (Schilling et al, 2009). Briefly, the tyrosine phosphatase SHP‐1, which is constitutively expressed (Supplementary Figure S2B) has previously been shown to require binding to the specific phosphorylated tyrosine residue Tyr429 on the EpoR to activate phosphatase activity (Pei et al, 1994; Klingmüller et al, 1995). Thus, we included an activated form of SHP‐1 in our model that in turn catalyzes dephosphorylation of JAK2 and EpoR. Phosphorylation of STAT5 by JAK2 is mediated by recruiting STAT5 to phosphotyrosine 343 and 401 of the EpoR cytoplasmic domain that functions as a scaffold (Gobert et al, 1996; Klingmüller et al, 1996; Barber et al, 2001). Consistent with the observation that phosphorylation of STAT5 is impaired in the absence of docking sites (Gobert et al, 1996), we considered that the STAT5 phosphorylation by the complex pEpoRpJAK2 must be faster than by the complex EpoRpJAK. After dimerization, the phosphorylated STAT5 translocates into the nucleus, binds to consensus sequences on the DNA and migrates back into the cytoplasm to enter a new cycle of activation (Swameye et al, 2003). These different forms of STAT5 were integrated in the model by three distinct variables, unphosphorylated STAT5 in the cytosol (STAT5), phosphorylated STAT5 dimers in the cytosol (pSTAT5) and phosphorylated STAT5 dimers in the nucleus (npSTAT5). To enter a new cycle of activation, STAT5 is dephosphorylated and exported from the nucleus into the cytoplasm in our model by a single reaction. Since various phosphatases are known to facilitate dephosphorylation of STAT5 in the nucleus and cytoplasm (Aoki and Matsuda, 2000; Yu et al, 2000; Hoyt et al, 2007), we assumed constitutive inactivation of STAT5 during traffic from nucleus to cytoplasm (npSTAT5 → STAT5).

Mathematical model of dual negative feedback regulation of JAK2‐STAT5 signaling. The model is represented as process diagram and reactions are modulated by enzyme catalysis (circle‐headed lines) or inhibition (bar‐headed lines). Dashed‐dotted lines indicate delayed reactions used for RNA transcription. Prefix ‘p’ represents phosphorylated species and ‘n’ represents nuclear species. Binding of the ligand Epo to its cognate receptor results in phosphorylation of JAK2, which becomes activated. In a second step, EpoR is phosphorylated in the complex EpoRpJAK2. The pEpoRpJAK2 complex was modeled by considering different subclasses of phosphorylated tyrosines at the receptor (see Supplementary Material Section 2.2). STAT5 is recruited by pEpoR and phosphorylated by pJAK2, dimerizes and translocates to the nucleus. As target genes of STAT5, the negative feedback proteins CIS and SOCS3 are expressed, which inhibit the pathway. For an extended model description, see Supplementary information.

Based on our experimental evidence shown in Figure 1, the STAT5 target genes CIS and SOCS3 were included in the model as transcriptional feedback regulators. The production of SOCS3 and CIS transcripts was modeled using a delay that is based on the linear chain trick (MacDonald, 1976). Based on their previously identified molecular mechanisms of inhibition, the effect of CIS and SOCS3 was included in the model. CIS was considered as inhibitor of STAT5 activation mediated by the pEpoRpJAK2 complex, because it competes with STAT5 for the same pTyr binding site on the receptor. In contrast, SOCS3 is known to bind directly to JAK2 and inhibit JAK2 activation and downstream events, i.e., phosphorylation of EpoR and STAT5. Additionally, it can bind to the phosphorylated EpoR and inhibit the recruitment of STAT5 to the specific pTyr binding sites on the EpoR as described before for CIS. Therefore, SOCS3 in the model attenuates STAT5 phosphorylation mediated by both EpoRpJAK2 and pEpoRpJAK2. (A detailed description with references is provided in the Supplementary Material Section 2.2.) Although the effect of JAK2 inhibition by transcriptional feedback regulators was not significant as shown in Figure 1, we included the inhibition of JAK2 by SOCS3 since the effect could be obscured by other attenuation mechanisms, i.e., deactivation of JAK2 by SHP‐1. To reduce the model complexity to the requirements of our biological question, which is important for achieving structurally and practically identifiable model parameters (Raue et al, 2009; Chen et al, 2010), we systematically eliminated short‐lived protein–protein complexes and fast reactions. Finally, the model comprised 29 parameters describing 36 reactions and 86 nuisance parameters, such as scaling and offset parameters of the experimental data. Considering the total of 541 experimental data points the estimation uncertainties of the model parameters are small enough to allow for accurate model predictions (see Supplementary Material Section 2.5).

Calibration of the dual transcriptional feedback model

To estimate the model parameters multiple data sets of Epo‐induced JAK2/STAT5 signaling were acquired in primary CFU‐E cells using three different experimental techniques for quantitative data generation. By quantitative immunoblotting, the time‐resolved activation kinetics of the experimental observables phosphorylated and total EpoR and JAK2, phosphorylated and total STAT5 as well as induction of CIS and SOCS3 expression were monitored (Figure 3A). Furthermore, the initial values of all proteins were determined (Supplementary Figure S3; Supplementary Table S1) and dose–response experiments were performed by determining the phosphorylation levels at the time point of maximal activation for different Epo concentrations (Supplementary Figures S20–S23). Since for phosphorylated proteins only relative values were obtained by immunoblotting, we used quantitative mass spectrometry (MS) with isotope‐labeled standard peptides to determine site‐specific phosphorylation degrees of STAT5 (Figure 3B; Supplementary Figure S4). Both STAT5a and STAT5b are activated in response to Epo stimulation and share functional roles in erythroid cells. Therefore, we performed MS analysis on STAT5b as representative for both isoforms. This approach revealed that after 10 min of Epo stimulation 73% of STAT5 molecules are phosphorylated on Tyr694. In combination with the absolute values of the total proteins that were determined by quantitative immunoblotting, total concentrations of phosphorylated STAT5 were calculated with the model. Additionally to the observations at the protein level, expression kinetics of CIS and SOCS3 were determined at mRNA levels using quantitative RT–PCR to validate the data from the microarray experiments and to improve the temporal resolution at early time points for parameter estimation (Figure 3C).

Model calibration with experimental data of JAK2‐STAT5 signaling obtained by different experimental techniques. For all experiments, primary CFU‐E cells were starved and stimulated with 5 U/ml Epo. At the indicated time points, samples were subjected to (A) quantitative immunoblotting, (B) mass spectrometry analysis or (C) qRT–PCR. Experimental data (black circles) with estimated standard errors and trajectories of the best fit (solid lines) are represented. Mass spectrometry data represent replicates of four independent experiments. (For additional experimental data used for the model calibration, see Figure 4 and Supplementary Figures S11–S23.) In total, 531 data points representing 18 different experimental conditions were used for model calibration. Source data is available for this figure at www.nature.com/msb.

To provide the required information content necessary to disentangle the role of the dual negative feedback regulation, we also performed experiments applying specific perturbations. Time‐resolved activation of the pathway was monitored in primary CFU‐E cells that overexpressed SHP‐1, CIS or SOCS3 (Figure 4). We focused on overexpression since RNAi‐mediated knockdown was not addressable with the currently available techniques in these cells (see Supplementary Material Section 1.2). The number of time points that could be examined was limited since in primary cells the retroviral transduction efficiency is quite low (∼20%) and only positively transduced cells were used for the experiments. Since the level of the overexpressed proteins is high at the start of each experiment, already a reduction of the initial peak is observed in addition to the impact on the steady‐state phosphorylation level. SOCS3 overexpression leads to a reduction in the phosphorylation of JAK2, whereas CIS overexpression does not have a major effect on JAK2 phosphorylation. This is in line with the previously proposed inhibitory mechanisms of both proteins, i.e., SOCS3 is blocking not only the binding of STAT5 to the receptor but can also inhibit JAK2 via its kinase inhibitory region (KIR) domain. CIS instead does not contain a KIR domain but is known to compete with STAT5 for binding to the EpoR phosphotyrosine 401 (Yoshimura et al, 1995; Ketteler et al, 2003).

Experimental data of JAK2‐STAT5 signaling under perturbed conditions used for model calibration. CFU‐E cells were retrovirally transduced with SHP‐1 (SHP1oe), SOCS3 (SOCS3oe) or CIS (CISoe). Positively transduced cells were starved and stimulated with 5 U/ml Epo for 60 or 120 min, respectively. Cellular lysates were subjected to immunoprecipitation and immunoblotting to determine activation profiles of JAK2‐STAT5 pathway components for overexpression conditions compared with control cells. Trajectories of the best fit are indicated (solid lines) with experimental data (circles) and estimated standard errors. For pEpoR in the CIS overexpression experiment, two data sets were combined. (For further information and additional experimental data used for the model calibration, see Figure 3 and Supplementary Figures S11–S23.) Source data is available for this figure at www.nature.com/msb.

Based on these data sets, the parameters were simultaneously estimated by maximizing the likelihood with the MATLAB implementation of the trust‐region method with user supplied derivatives (Coleman and Li, 1996). In order to allow for normally distributed measurement noise the likelihood was evaluated in logarithmic concentration space (Kreutz et al, 2007). To account for uncertainties in the parameter estimates, the identifiability was investigated and 95% confidence intervals were calculated by applying the profile likelihood approach (Raue et al, 2009). The model parameters were structurally and, in most cases, also practically identifiable (see Supplementary Model Results 2.5). In summary, our dual feedback regulation model was able to adequately represent the experimental data sets for normal and perturbed conditions, dose–response behavior of various signaling components as well as mRNA expression and absolute phosphorylation amounts determined based on MS analysis.

Next, we applied the calibrated dual feedback regulation model to examine the link of the behavior of STAT5 and the two negative feedback regulators with regard to physiological responses of cells. Since STAT5 is known to induce anti‐apoptotic target gene expression in CFU‐E cells, we asked if a direct relationship exists between the amount of phosphorylated STAT5 in the nucleus (npSTAT5) and the extent of survival in CFU‐E cells. We focused on the integral response, which is the area under the curve, because this entity captures both the kinetics and magnitude of the signal. It has been demonstrated that processes with slow kinetics (i.e., genetic networks and cell fate decisions) downstream of processes with fast kinetics (i.e., phosphorylation‐based signaling networks) act as integrators capable of measuring how long the upstream signal has been on (Behar et al, 2007). Moreover, besides cytokine responses, in MAPK signaling the integral response of a transcription factor was used previously to link ERK activity and DNA synthesis (Asthagiri et al, 2000). Using the dual feedback regulation model we could quantitatively assess the integral response of STAT5 even for very low Epo concentrations, which are intractable in experimental approaches but are of important physiological significance. For instance, the STAT5 phosphorylation levels elicited by Epo concentrations lower than 0.1 U/ml are below the detection limit of immunoblotting experiments, but sufficient to promote cell survival. We predicted the trajectory of npSTAT5 response for wild‐type cells and additionally for cells overexpressing either SOCS3 or CIS at different Epo concentrations and calculated the area under the curve. As an example of this procedure, the model prediction of the integral STAT5 response at a single Epo concentration is depicted (Figure 5A). As shown in the overlays of the different experimental conditions, SOCS3 overexpression is predicted to reduce the integral of pSTAT5 in the nucleus much more than overexpression of CIS.

Linking the integral response of phosphorylated STAT5 in the nucleus to the survival rate of CFU‐E cells. (A) The calibrated model was used to simulate the integral response of phosphorylated STAT5 in the nucleus (npSTAT5) over the broad physiological range of Epo concentrations in wild‐type, CIS and SOCS3 overexpressing (oe) cells. A representative example at Epo=10−6.78 U/cell is depicted. (B) TUNEL assay to determine the fraction of apoptotic cells. Wild‐type CFU‐E cells and cells overexpressing SHP‐1, CIS or SOCS3 were cultured 24 h in various Epo concentrations. Histograms show the representative result of a TUNEL assay with Epo=10−6.78 U/cell. Cells for positive control were treated with DNaseI for 10 min. (C) Overlay of scaled integral response of npSTAT5 (60 min) including 95% confidence bands (shades) and experimentally determined survival rates of CFU‐E cells for wild‐type cells, CIS and SOCS3 overexpressing cells (circles). Source data is available for this figure at www.nature.com/msb.

To evaluate whether the predicted integral STAT5 response correlates with the extent of survival, we experimentally determined the percentage of apoptotic and viable primary CFU‐E cells using the terminal deoxynucleotidyl transferase‐mediated dUTP nick end‐labeling (TUNEL) assay in combination with flow cytometry at various Epo concentrations for wild‐type and CIS or SOCS3 overexpressing cells. Here again, only four different Epo concentrations could be determined in CIS and SOCS3 overexpressing cells due to the low transduction efficiency that results in limited primary cell material. A representative example at Epo=10−6.78 U/cell is displayed (Figure 5B), which demonstrates that SOCS3 overexpression affects survival to a larger extent than CIS considering similar overexpression levels (Supplementary Figure S5A).

To systematically test if npSTAT5 is the major factor that contributes to survival decisions, we additionally parameterized the contributions to the survival signal of the different pathway species (pEpoR, pJAK2, npSTAT5, CIS and SOCS3) in an additive way (see Supplementary Material Section 1.1). By comparing the contributions of various signaling components to the anti‐apoptosis response, we detect that the npSTAT5 signal contributed >99%.

As apoptosis is an all‐or‐none decision, we performed simulations describing the effect of intrinsic (Supplementary Figure S8) and extrinsic stochasticity (Supplementary Figure S9) on npSTAT5. In both cases, the trajectory of npSTAT5 is rather negligibly affected by noise and does not show bi‐stability. To further analyze the direct correlation between the integral response of phosphorylated STAT5 in the nucleus and the survival decision of erythroid progenitor cells (Figure 6C), we investigated the induction of anti‐apoptotic target genes. First, global gene expression data were used to identify a panel of Epo‐regulated anti‐apoptotic genes. Among the preselected candidates (Bcl‐xL, BIM, Bcl‐2 and Pim‐1) that were previously associated with Epo‐dependent regulation, exclusively Pim‐1 showed rapid and significant induction compared with untreated control cells (Supplementary Figure S7A). The dose‐dependent behavior of Pim‐1 (Supplementary Figure S7B) mirrored the behavior of the integral npSTAT5 response (Figure 6C) and continuously increased over different Epo doses in an averaged cell population. Thus, the integral of nuclear phosphorylated STAT5 correlates with Pim‐1 induction, while the individual cell fate may be influenced by cell‐to‐cell variability (Spencer et al, 2009).

Dual negative feedback with divided function in JAK2‐STAT5 signaling. (A) The steady‐state level of phosphorylated STAT5 in the nucleus was simulated in the presence of only one transcriptional negative regulator, CIS or SOCS3, and in the absence of both. The increase of pSTAT5 steady‐state levels was calculated relative to wild‐type cells (black line) at t=360 min. Dashed lines indicate upper and lower 95% confidence bands for the prediction. For two exemplary Epo concentrations, (I) Epo=10−9 U/cell and (II) Epo=10−6 U/cell, the time profiles of npSTAT5 are shown. (B) Model prediction and experimental data for the increase of phosphorylated JAK2 relative to phosphorylated Epo receptor with rising Epo concentrations at t=30 min. (C) Mechanistic scheme explaining dual negative feedback with divided function in JAK2‐STAT5 signaling for different Epo concentration ranges. The concentration dependency of the inhibitory effects of CIS and SOCS3 is caused by the increasing fraction of phosphorylated JAK2 compared with phosphorylated EpoR. At low Epo levels, CIS mostly impacts pEpoR‐dependent STAT5 activation by preventing binding of STAT5 via its SH2 domain to the specific phosphotyrosine sites on the EpoR receptor. At high Epo concentrations, SOCS3 mostly impacts pJAK2‐dependent STAT5 phosphorylation by inhibiting the activation of JAK2 via its kinase inhibitory region (KIR).

As the duration of nuclear activated STAT5 that is necessary to induce a survival response is unknown, we tested which integration time showed the strongest correlation with the experimentally determined survival. Interestingly, an integration time of ∼60 min correlated best with the experimental results as demonstrated by the log‐likelihood profile (Supplementary Figure S5B). The overlay of the simulated npSTAT5 response integrated for 60 min and the experimental survival data of CFU‐E cells at various Epo concentrations demonstrate a close relationship of both quantities (Figure 5C). Since the model was calibrated with biochemical experiments, in which different cell densities were applied than in the apoptosis assays, Epo concentrations U/ml were converted to U/cell in order to achieve comparability of conditions. The different effects of SOCS3 and CIS overexpression were adequately reflected by the model considering the obtained confidence intervals for the parameter estimates and model trajectories (see Supplementary Model Results 2.5).

We concluded that the dual feedback regulation model captured the inhibitory effects of the negative feedback regulators over a broad Epo concentration range.

Since CIS and SOCS3 inhibit the pathway by different mechanisms, we asked whether they might have distinct roles in controlling STAT5 phosphorylation over the entire range of physiological Epo concentrations that can vary by 1000‐fold (Jelkmann, 2004). To investigate the impact of the transcriptional feedback regulators CIS and SOCS3 on the STAT5 steady‐state phosphorylation level, we calculated the fold change of STAT5 steady‐state phosphorylation at 360 min post Epo stimulation for in silico cells lacking either both negative feedback regulators (CIS‐SOCS3 K.o.) or solely CIS (CIS K.o.) or SOCS3 (SOCS3 K.o.) relative to phosphorylation level in wild‐type cells (Figure 6A). The increase of the pSTAT5 steady‐state level was simulated over a broad Epo concentration range per cell encompassing basal (I) as well as acute (II) conditions. Exemplarily, the impact of the absence of both feedback regulators or of both individually on the time course of STAT5 phosphorylation is depicted for basal and acute Epo concentrations (Figure 6A, lower panels). As expected, we observed an increase in the pSTAT5 steady‐state level for the double knockout simulation over the entire Epo range. Interestingly, the single knockout simulations revealed that CIS and SOCS3 have a distinct impact on pSTAT5 phosphorylation at different Epo concentrations. This was unexpected since both proteins are expressed to similar extents and time frames after stimulation with Epo doses applied in the experiments. At low Epo concentrations, the absence of CIS resulted in an increase in the amplitude and the steady‐state level of STAT5 phosphorylation while the absence of SOCS3 had no major impact. Conversely, at high Epo concentrations the steady‐state level of STAT5 phosphorylation is increased in the absence of SOCS3 whereas CIS has no major effect. This synergistic action of the two feedbacks allows adjusting STAT5 phosphorylation levels over the entire concentration range of physiological Epo concentrations from basal levels up to extremely high values that occur under acute hypoxic stress.

To gain insights into the underlying molecular mechanism of this dose‐dependent inhibitory effect of CIS and SOCS3, we investigated the model structure and dynamics of each species of the pathway. By analyzing the components in dependence of ligand doses, we identified the observable pJAK2 relative to pEpoR as a potential explanation for the different inhibitory effects dependent on the Epo dose. With rising Epo concentrations the ratio of pJAK2 to pEpoR changes, because more JAK2 is phosphorylated relative to the EpoR (Supplementary Figure S6). To provide experimental evidence for this model hypothesis, we compared dose responses of pJAK2 relative to the response of pEpoR determined by immunoblotting with the model predictions (Figure 6B). Interestingly, we observed that the model simulation of increasing pJAK2 relative to pEpoR is in good agreement with the experimentally determined ratio of both quantities at 30 min with rising Epo concentrations. Since SOCS3 is inhibiting the activation of STAT5 via both pJAK2 and pEpoR whereas CIS only blocks binding of STAT5 to the EpoR, the impact of SOCS3 increases with increasing ratio of pJAK2 versus pEpoR. Taken together, these data and the model suggest that SOCS3 inhibition becomes more important with the relative increase of JAK2 phosphorylation and thus at higher Epo doses (Figure 6C). CIS instead inhibits the binding of STAT5 to EpoR; and thus, its inhibitory effect dominates for lower Epo doses, when the ratio of pJAK2 to pEpoR is not enhanced.

Discussion

Based on extensive time‐resolved data sets in primary erythroid progenitor cells and mathematical modeling, we dissected the different contributions of the two feedback regulators CIS and SOCS3 and identified dual transcriptional feedback effective at discrete ligand concentrations as key property of the EpoR‐STAT5 signaling system. This result underlines the general principle of negative feedback regulation to increase control and stability of transcriptional responses over a broad range of input values (Becskei and Serrano, 2000; Freeman, 2000). Our observation that multiple feedbacks are necessary to orchestrate tight regulation of transcription factor activity from very low to very high ligand doses suggests a new strategy of SOCS‐mediated feedback control.

Most of the previous JAK/STAT models have a complex structure comprising a large number of variables and literature‐based parameters (Yamada et al, 2003; Zi et al, 2005; Soebiyanto et al, 2007). In this work, as an alternative modeling strategy, a bottom‐up approach was employed, which aims at fully identifiable parameters that are essential to obtain models with high predictive power (Bruggeman et al, 2002; Aldridge et al, 2006). We therefore employed a large set of data measured under different perturbation conditions in combination with a model structure comprising the minimal number of parameters necessary to explain the data. Parameter identifiability was analyzed by the profile likelihood approach (Raue et al, 2009). Applying this method, we could establish a dual negative feedback model of JAK2‐STAT5 signaling with structurally and in most cases also practically identifiable parameters.

Qualitatively, the important role of STAT5 as crucial regulator of survival and differentiation of erythroid progenitor cells has been established previously (Socolovsky et al, 2001; Yao et al, 2006; Grebien et al, 2008; Zhu et al, 2008). We demonstrate that the calculated integrated response of pSTAT5 in the nucleus accurately correlates with the experimentally determined survival of CFU‐E cells, thereby providing a quantitative link of the dependency of primary CFU‐E cells on pSTAT5 activation dynamics. In line with this, a recent study that quantitatively determined pSTAT5 activity levels in IL‐2 induced single T cells reported on a causal link between enhanced pSTAT5 levels and decreased cell death (Feinerman et al, 2010). We concluded that fine‐tuned changes in STAT5 activation levels are a critical determinant for cell fate decisions. This observation is also true during an earlier differentiation stage of the erythroid lineage. A previous study showed that partial depletion of STAT5 from CD34+ cells by a lentiviral RNAi approach in the presence of thrombopoietin and stem cell factor resulted in a decrease of erythroid progenitors (BFU‐E). Conversely, overexpression of an activated STAT5a mutant strongly induced erythroid differentiation (Olthof et al, 2008).

Our finding that the integral of npSTAT5 is a good predictor for survival does not imply that survival is exclusively depending on STAT5. The integral of STAT5 activity in the nucleus transfers quantitative information about extracellular ligand concentrations to downstream signals, i.e., expression of anti‐apoptotic target genes such as Pim‐1 that contribute to the ultimate cellular response. Additional pro‐survival factors such as the PI3K/AKT pathway that have been shown previously to be involved in prevention of apoptosis in CFU‐E cells (Bouscary et al, 2003) may also contribute. However, overexpression of constitutively active AKT could not substitute for the apoptosis‐suppressing function of the EpoR‐STAT5 pathway in JAK2−/− erythroid cells (Ghaffari et al, 2006). Hence, though we cannot rule out the involvement of other pathways, our study underlines the direct relationship of the integral STAT5 response and survival decisions of primary erythroid progenitor cells.

A major bottle‐neck in combining signal transduction events with cellular phenotypes is the discrepancy in the time scale and stimuli concentrations that are applied in the different experiments. The sensitivity of biochemical assays to determine phosphorylation events within minutes or hours after stimulation is usually lower than the threshold of sensitivity in assays to determine the physiological response after one or more days. By employing our mathematical model, we were able to compute the integrated response of STAT5 at Epo concentrations that are beyond the threshold of biochemical experiments. This allowed us to directly link the activation status of the transcription factor with a cellular response. By correlation analysis, we could identify the early signaling phase (⩽1 h) of STAT5 to be most predictive for survival decisions, which was determined ∼24 h later. Interestingly, an earlier study that investigated a compendium of signaling pathways and responses in TNF‐, EGF‐ and insulin‐treated HT‐29 cells revealed similar results in the predictability of apoptosis‐survival cell fate decisions by early signaling events (Gaudet et al, 2005). The authors report that early signaling events between 5 and 90 min immediately downstream of the receptors were more predictive than 2–8 h, presumably because information is forwarded to downstream transcriptional events. Thus, we hypothesize that as a general principle in apoptotic decisions, ligand concentrations translated into kinetic‐encoded information of early signaling events downstream of receptors (⩽1 h) can be predictive for survival decisions 24 h later.

We could show by experimental data and mathematical modeling that SOCS3 and CIS reduce the steady‐state STAT5 phosphorylation level after ∼1 h of stimulation. Our results at a single Epo concentration are consistent with other studies that investigated the regulation of STAT phosphorylation dynamics by SOCS proteins. As one example, SOCS3 that is induced upon IL‐6 stimulation controls the late phosphorylation profile of STAT1 and STAT3 in the liver. This was shown in an experimental study using a conditional SOCS3 knockout, which resulted in enhanced phosphorylation of STATs, although the amplitude was unaffected (Croker et al, 2003). The importance for adjusting the steady‐state phase of transcription factor activity by transcriptional feedback regulators was shown also for other signaling pathways. EGF‐induced ERK signaling in HeLa cells is increased 1 h after stimulation when translation is blocked by cycloheximide (Amit et al, 2007). This temporal regulation pattern, including short‐term deactivation by constitutively expressed phosphatases and late‐phase modulation by slow transcriptional feedbacks may evolve as a general paradigm for tightly controlled cytokine‐induced signaling pathways (Legewie et al, 2008). The existence of multiple overlapping feedback regulation mechanisms, each with distinct temporal characteristics, ensures the effective control of the signaling dynamics to appropriately fine‐tune cellular responses.

The fact that one or two SOCS family members are induced upon stimulation appears to be a general principle in cytokine‐induced signaling, although there is no simple relationship between a particular cytokine and the pattern of respectively induced SOCS mRNAs (Krebs and Hilton, 2001). Induced in a cell type‐specific and cytokine‐dependent manner, SOCS proteins can distinctively coordinate cytokine‐induced responses. A prime example for the evidence that two SOCS proteins uniquely attenuate STAT dynamics is the study of Wormald et al (2006) showing that SOCS1 and SOCS3 regulate STAT1 and STAT3 phosphorylation dynamics with different kinetics. Yet, the impact of different ligand concentrations has not been considered so far. By using model simulations to analyze the inhibitory effects of CIS and SOCS3 on the STAT5 phosphorylation level at previously unobserved Epo concentrations, our results revealed that the two feedbacks are most effective at different Epo concentration ranges. Conventional experimental techniques, however, cannot resolve the entire dynamic range of ligand concentrations. According to our mathematical model, the major role of CIS is modulating STAT5 phosphorylation levels at low, basal Epo concentrations, whereas SOCS3 is essential to control the STAT5 phosphorylation levels at high Epo doses.

As potential molecular mechanism of this dose‐dependent inhibitory effect, we could identify the quantity of pJAK2 relative to pEpoR that increases with higher Epo concentrations. Since SOCS3 can inhibit JAK2 directly via its KIR domain to attenuate downstream STAT5 activation, SOCS3 becomes more effective with the relative increase of JAK2 activation. These observations are further supported by the CIS and SOCS3 overexpression experiments that showed the strong inhibition of JAK2 phosphorylation only by SOCS3 in primary CFU‐E cells. Moreover, in cells overexpressing SOCS3 the STAT5 phosphorylation level was more extensively reduced and survival was more decreased than in CIS overexpressing cells considering similar overexpression levels.

Our findings raise the question why it is important for the cell to tightly control the long‐term steady‐state signaling level of STAT5 by transcriptional feedback regulators over the entire range of high and low Epo doses although the first hour of STAT5 activation is predictive for the survival decision. We hypothesize that when the decision for survival has occurred, it is essential to constrain signaling to a residual steady‐state level in order to prevent aberrant events that could lead to uncontrolled erythroid progenitor growth. Constitutive phosphorylation of the JAK2/STAT5 pathway caused by activating JAK2 mutations has a crucial role in the onset of polycythemia vera (PV), a disease that is characterized by the formation of endogenous colonies with Epo‐independent differentiation (Prchal and Axelrad, 1974; Weinberg et al, 1989; Kota et al, 2008). Moreover, human erythroid progenitor cells transduced with a constitutive phosphorylated form of STAT5 were reported to survive, proliferate and differentiate in the absence of Epo and in this way mimic the PV phenotype. The essential requirement for progenitor cells to tightly constrict the Epo input signal after 1 h of stimulation over the broad range of physiological Epo concentrations that can vary over 1000‐fold is already apparent at the upstream receptor level. Different studies have shown that activation of EpoR and JAK2 is rapidly terminated within the first hour by dephosphorylation as well as internalization and degradation of Epo (Gross and Lodish, 2006; Becker et al, 2010).

We propose that the transcriptional feedback proteins CIS and SOCS3 are required to tightly adjust the phosphorylation level of STAT5 after 1 h of stimulation. This hypothesis is supported by reports demonstrating the crucial role of SOCS3 in embryonic development. Mice lacking the SOCS3 gene exhibit embryonic lethality at days E12–16. Marine et al (1999) showed that these mice display erythrocytosis with dramatic expansion of erythropoiesis within the fetal liver as well as throughout the embryo. Roberts et al (2001) demonstrated that the death is associated with abnormalities in the placenta. Vice versa, enforced expression of SOCS3 in vivo specifically suppressed fetal liver erythropoiesis (Marine et al, 1999). Moreover, a loss‐of‐function mutation of SOCS3 has been proposed to contribute to the onset of myeloproliferative disease in PV patients (Suessmuth et al, 2009). In cancer, high input doses are mimicked by aberrant activation of JAKs that are frequently mutated in many tumor cells. Interestingly, there are numerous recent studies showing that in several malignant tumors JAK activating mutations are complemented by gene silencing of SOCS3 and SOCS1, the two SOCS members which contain a KIR domain (He et al, 2003; Chim et al, 2004; Johan et al, 2005; Jost et al, 2007). This underlines the essential role of SOCS3 to control high input doses. Thus, the absence of SOCS3 severely impacts the growth and survival of erythroid progenitor cells and is essential to abrogate signaling downstream of the EpoR at long‐term steady‐state phosphorylation levels upon high upstream input conditions.

In contrast to SOCS3, our model predicts that CIS acts most efficiently at a single point of the network at low ligand concentrations. This is line with other reports that demonstrated the major role of CIS as a specific competitive binding inhibitor of STAT5 at the pY401 position of the EpoR (Matsumoto et al, 1997; Verdier et al, 1998; Ketteler et al, 2002). The specific inhibition of STAT5‐mediated responses by CIS is also supported by transgenic mice that overexpress CIS. These animals display diminished expression of STAT5‐mediated responses in growth hormone and prolactin signaling, similar to STAT5ΔN/ΔN knockout mice (Matsumoto et al, 1999). In contrast to SOCS3 knockout mice, CIS knockout mice are viable, but show an increase in hematopoietic progenitor cells (Kubo et al, 2003), which is in line with our model prediction of CIS as modulator of fine‐tuned pSTAT5 responses at basal level Epo input.

In summary, our mathematical approach provided new insights into the specific function of feedback regulation in STAT5‐mediated life or death decisions of primary erythroid cells. We dissected the roles of the transcriptionally induced proteins CIS and SOCS3 that operate as dual feedback with divided function thereby facilitating the control of STAT5 activation levels over the entire range of physiological Epo concentrations. The detailed understanding of the molecular processes and control distribution of Epo‐induced JAK/STAT signaling can be further applied to gain insights into alterations promoting malignant hematopoietic diseases.

Normalization was performed in the R environment together with the Bioconductor toolbox (http://www.bioconductor.org/) using the Robust Multichip Average (RMA) for background adjustment, quantile normalization and summarization (Irizarry et al, 2003). Subsequent probe annotation was handled with the Affymetrix mouse4302 annotation package in Version 2.4.5. If multiple probes map to the same Gene ID, the one with the largest test inter‐quartile range (IQR) among all time points was selected. The expression data were deposited in the GEO database under accession number GSE26151. All expression levels were log2 transformed and gene fold expression was calculated with respect to the mean gene expression at 0 h. The error in fold expression ±0.43 has been estimated from the sample IQR of the gene expression differences of the biological duplicates at 0 h.

Quantitative RT–PCR

Total RNA from 3 × 106 primary erythroid progenitor cells stimulated with 0.5 U/ml Epo was isolated using the RNeasy Mini Plus Kit (Qiagen) according to the manufacturer's instructions. For quantitative RT—PCR, cDNA was generated with the QuantiTect Reverse Transcription Kit (Qiagen) and analyzed using the LightCycler 480 with the hydrolysis‐based Universal ProbeLibrary platform (Roche Diagnostics). Crossing point values were calculated using the Second Derivative Maximum method of the LightCycler 480 Basic Software (Roche Diagnostics). PCR efficiency correction was performed for each PCR set‐up individually. Relative concentrations were normalized using HPRT as a reference gene. UPL Probes and primer sequences were selected with the Universal ProbeLibrary Assay Design Center (Roche Diagnostics). According to the manufacturer's instruction, the selected probes are intron spanning, thereby ensuring the amplification of fully matured mRNA. UPL Probes and primer sequences were CIS_for 5′‐gacatggtcctttgcgtaca‐3′, CIS_rev 5′‐atgccccagtgggtaagg‐3′, Probe#1 5′‐cctggagc‐3′ and SOCS3_for 5′‐atttcgcttcgggactagc‐3′, SOCS3_rev 5′‐aacttgctgtgggtgaccat‐3′, Probe#83 5′‐cagccacc‐3′.

Detection of apoptosis by TUNEL assay

Labeling of DNA strand breaks using the TUNEL assay and analysis by flow cytometry allows for quantitative analysis of apoptosis at single cell level. Briefly, cells (3 × 105 cells/ml) were starved 1 h and stimulated with different Epo concentrations over 20 h. 2 × 106 CFU‐E cells were fixed in 2% paraformaldehyde and permeabilized with 0.1% Triton X‐100, 0.1% sodium citrate for 2 min. The positive control was treated with DNase I (Roche Diagnostics) for 15 min. After washing with PBS/0.3% BSA, cells were resuspended in TUNEL reaction mixture (Roche Diagnostics) to label free 3′‐OH groups of DNA. Fluorescein‐dUTP incorporated in nucleotide polymers was detected and quantified by flow cytometry using the FACSCalibur (Becton Dickinson).

MS analysis

MS analyses were performed to measure site‐specific phosphorylation degrees of STAT5b in starved unstimulated (0 min) and Epo‐stimulated (10 min) murine CFU‐E cells. Following immunoprecipitations of endogenous STAT5, proteins were separated on 1D‐PAGE and stained with Coomassie SimplyBlue™ SafeStain (Invitrogen). After broad cutting, STAT5 bands were in‐gel digested with trypsin as described (Seidler et al, 2009). MS analyses were performed using a nanoAcquity UPLC (Waters) coupled to a LTQ‐Orbitrap XL (Thermo) as described (Seidler et al, 2010). Data were acquired using a Top3 Data‐dependent acquisition: MS/MS spectra of the three most intense precursors with charge ⩾2 were recorded. Site‐specific phosphorylation degrees of STAT5b were determined by using a homologous and labeled one‐source peptide/phosphopeptide standard (Hahn et al, 2010). For STAT5b, it was observed that the homologous standard pair used delivers virtually identical quantitative data compared with an isotopically labeled standard pair. MS data were analyzed using Xcalibur 2.0.6 and MASCOT 2.2.2.

Model simulation and parameter estimation

The ODE system consists of 25 dynamic variables and was solved by applying the CVODES solver (Hindmarsh et al, 2005) with an absolute and relative accuracy of 10−8. The solver was compiled as C‐executable for MATLAB together with ODE equations. The ODE system was evaluated for 24 different experimental conditions. For numerical efficiency, these 24 independent calculations of the ODE system were parallelized. Including all experimental settings the model consisted of 115 unknown parameters. Kinetic parameters depending on unknown concentration scales have been disentangled from the respective scaling parameters by suitable transformations (as described in the Supplement). This facilitates a more efficient estimation of the unknown scaling factors. In the case of mRNA measurements of CIS and SOCS3, no experimental information was available about the absolute concentrations. Taking into consideration the independence of the kinetic parameters from the scaling parameters due to the applied parameter transformations, the internal scale of these two mRNA components was fixed to one. Consequently, the model cannot provide predictions of absolute mRNA concentration. For the absolute concentration of the EpoR_JAK complex, a literature value was included as prior information for concentration scale of the receptor complex.

The model parameters were estimated by maximum likelihood estimation taking into account a total of 541 data points. To allow for global parameter estimates, we scanned the parameter space, beginning from reasonable initial guesses by taking into account the time scale of the measurements, over orders of magnitude along the profile likelihood (Raue et al, 2009). For each step in the calculation of the profile likelihood, we applied the MATLAB implementation of the trust‐region method (lsqnonlin) with user supplied derivatives. For the calculation of the derivatives, the sensitivity equations (Leis and Kramer, 1988) were simultaneously solved together with the original ODE systems, which is featured by the CVODES solver efficiently. To improve convergence of the estimation and due to the positive definite nature of the parameters such as rate constant, initial concentrations, concentration offsets and scales, and noise magnitudes, the parameters were estimated in logarithmic parameter space. In order to provide for normally distributed measurement noise the likelihood was evaluated in logarithmic concentration space (Kreutz et al, 2007) transforming the mainly multiplicative nature of the noise to an additive one. The magnitude of the additive measurement noise was estimated simultaneously. For each experimental technique and, in the case of immunoblotting, for each detection antibody utilized, an individual noise parameter is included in the model. The calculations of the profile likelihood suggest that the global optimum has been found.

Conflict of Interest

The authors declare that they have no conflict of interest.

Supplementary Information

Acknowledgements

We thank Sandra Manthey, Susen Lattermann and Maria Saile for excellent technical assistance; Nao Iwamoto and Verena Becker for helpful discussions. We also thank Andrea C Pfeifer for providing information on the CFU‐E cell volume. We acknowledge funding by German Federal Ministry of Education and Research (BMBF)–funded MedSys‐Network LungSys (JB, AR, DK, MEB, HB, WL, JT and UK), the SBCancer network in the Helmholtz Alliance on Systems Biology (MS, WL, JT and UK), the EU FP6 project ‘Computational Systems Biology of Cell Signalling’ (LSHG‐CT‐2004‐512060) (JB, JT and UK) and the Excellence Initiative of the German Federal and State Governments (HB, CK and JT).

Author contributions: JB performed experiments and wrote the manuscript; AR performed mathematical modeling and wrote the manuscript; MS performed experiments (Figure 6B; Supplementary Figures S16 and S17) and designed research; MEB and WDL carried out MS data generation and analysis (Figure 3B); CK and DK performed statistical analysis; NG carried out microarray data generation; HB performed microarray data analysis (Figure 1C); JT and UK designed research and wrote the manuscript.

This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.