Abstract
Background
Celltocell variability in protein expression can be large, and its propagation through signaling networks affects biological outcomes. Here, we apply deterministic and probabilistic models and biochemical measurements to study how network topologies and celltocell protein abundance variations interact to shape signaling responses.
Results
We observe bimodal distributions of extracellular signalregulated kinase (ERK) responses to epidermal growth factor (EGF) stimulation, which are generally thought to indicate bistable or ultrasensitive signaling behavior in single cells. Surprisingly, we find that a simple MAPK/ERKcascade model with negative feedback that displays graded, analog ERK responses at a single cell level can explain the experimentally observed bimodality at the cell population level. Model analysis suggests that a conversion of graded input–output responses in single cells to digital responses at the population level is caused by a broad distribution of ERK pathway activation thresholds brought about by celltocell variability in protein expression.
Conclusions
Our results show that bimodal signaling response distributions do not necessarily imply digital (ultrasensitive or bistable) single cell signaling, and the interplay between protein expression noise and network topologies can bring about digital population responses from analog single cell dose responses. Thus, cells can retain the benefits of robustness arising from negative feedback, while simultaneously generating populationlevel on/off responses that are thought to be critical for regulating cell fate decisions.
Background
Development, growth and homeostasis of multicellular organisms depend on the ability of individual cells to convert noisy, analog signals into clear, yesorno cell fate decisions, such as apoptosis, proliferation and differentiation. One way that cells make such decisions is through the use of signal transduction systems that sense the strength of an analog input signal, and then convert it into one of several distinct activity states, such as “on” or “off” output states of highly ultrasensitive or bistable systems [13]. For example, various mitogen concentrations can cause bistable activation of cyclindependent kinases to drive cell cycle transition decisions [46]. Theoretical studies have shown that signaling networks containing positive or double negative feedback loops [3], opposing modification enzymes exhibiting saturation kinetics [1] and multisite modification cycles [2,7] can exhibit digital (bistable or ultrasensitive) behavior. However, not all networks that contain such motifs will necessarily exhibit digital behavior; such behavior arises from the cell’s precise tuning of quantitative, spatiotemporal aspects of the network. Indeed, the signal transduction network connecting epidermal growth factor (EGF) to activation of extracellular signalregulated kinase 1/2 (ERK) contains many elements that potentially can lead to switchlike behavior. However, previous single cell studies in different mammalian cell lines have reported both graded [8,9] and “allornothing” [10] EGFinduced ERK activation responses. One determinant of whether signaling is graded or switchlike is the spatial localization of signal processing proteins [11].
Under idealized conditions of celltocell homogeneity, experimental techniques such as immunoblotting that measure average population responses may be able to detect allornone signaling responses, as long as the celltocell variability in response activation thresholds are negligible [12]. However, it is becoming clear that the fundamental processes of transcription and translation are inherently stochastic, and give rise to significant celltocell variability in protein levels [1320]. The primary stochastic factors are (i) the rate of transcription, which is burstlike due to the low number (two) of genes for a particular protein in a cell [21,22] and (ii) the number of proteins produced per mRNA, which is random due to competition between ribosomes and RNase for the mRNA [13,23,24]. Protein degradation also contributes to expression noise, but usually to a lesser extent, since protein copy numbers are typically large enough to dampen the comparatively small stochastic fluctuations in degradation rate. Thus, even genetically identical cells show substantial variations in protein and mRNA abundance, and as a result, may also show differences in their signaling responses [25]. Because of such heterogeneity in protein abundance, population average measurements are not sufficient for investigating “allornothing” responses; singlecell measurement techniques capable of capturing the dynamics of digital signal transduction are needed [12].
Here, we use flow cytometry to measure EGFinduced, singlecell ERK activation responses in a HEK293 cell population. We observe bimodal response distributions in cell populations that are usually thought to indicate switchlike behavior in single cells. Surprisingly, an ERK cascade signaling model incorporating negative feedback and a graded, analog single cell dose response is shown to be consistent with the observed population responses. Our model analysis suggests that such a conversion of analog responses in single cells to digital responses at the population level is due to protein abundance variability, which gives rise to a broad distribution of ERK pathway activation thresholds and RasGTP levels. Thus, bimodal response distributions do not necessarily imply digital single cell signaling; such distributions can arise from the interplay between protein expression noise and negative feedbackmediated, analog singlecell responses.
Results
Analyses of ERK responses to EGF in individual cells and populations
We used a flow cytometrybased phosphorylation assay (FCPA) [26] to determine the kinetics and dose response of ERK activation by EGF in HEK293 cells. We show that population averages obtained from FCPA results correspond well to traditional Western blot measurements of activated (dually phosphorylated) ppERK levels in cell populations (Additional file 1: Figure S1). However, the FCPA also reveals how individual cells contribute to this collective population response (Figure 1AD; Additional file 1: Figure S2). The increase in mean values of ppERK was dosedependent after two minutes of EGF stimulation, suggesting that analog signaling has occurred in individual cells. However, a fraction of cells contain ppERK levels similar to those of the basal state. We refer to this feature of the distribution as a shoulder. Although the height of this shoulder decreases with increasing EGF dose, its position remains unchanged, indicating a dosedependent fraction of cells failing to activate ERK. At five minutes after EGF stimulation, the ppERK distribution is unambiguously bimodal, implying digital “onoff” behavior. Higher EGF doses increase the fraction of cells with high ppERK (”ERKon”) at the expense of the “ERKoff” population. Thus, in a dosedependent manner, EGF increases the probability that a cell will have ERK turned on. At later time points, a bimodal distribution persists at some EGF doses, while data from other doses show “shouldering” patterns similar to the behavior at 2 minutes. Thus, the EGFinduced ERK response on the population level is complex consisting of both analog and digital elements.
Additional file 1. This additional file contains all the supplementary figures along with their legends.
Format: PPT Size: 1.4MB Download file
This file can be viewed with: Microsoft PowerPoint Viewer
Figure 1. Cell population dose and dynamic response of ppERK to EGF. AD. Each panel corresponds to a fixed time after EGF stimulations: 2 min (A), 5 min (B), 10 min (C) and 30 min (D). In each panel, the different colors correspond to different EGF doses as indicated by the visual legend. Each distribution is compiled from 10,000 individual HEK293 cell responses as measured by flow cytometry, and is representative of between three and six independent experiments. Events were gated based on forward and side scatter to exclude debris, dead cells, and cell clusters. The xaxis is the magnitude of activated ERK (ppERK) in arbitrary fluorescence units, and the yaxis is the frequency of observing a particular level of fluorescence in a cell. E. Total ERK abundance data. The bestfit gamma distribution curve is depicted by the green line, while the black line shows experimental data. Data are representative of five independent experiments, and were fit to a gamma distribution. Mean and standard deviation of the fit parameters are k = 5.4; θ = 2.7 × 10^{5} [AU]. FI. Cells were stimulated with either 0.1 nM (F,H) or 1 nM (G,I) EGF for five minutes, and then analyzed by flow cytometry to measure ppERK and total ERK levels simultaneously. In FG, black curves correspond to ppERK distributions, and green curves correspond to normalized distributions where ppERK levels in each cell were divided by the total ERK signal intensity in the same cell. To compare the green and black curves on the same axis, intensities for each distribution are divided by their respective mean. In the HI dot plots, ppERK levels are on the xaxis, whereas total ERK levels are on the yaxis.
Next, we investigated how celltocell variability in total ERK abundance affects the ppERK responses. Measurements of the total ERK distribution by flow cytometry, as expected, revealed substantial celltocell variability in total ERK levels (Figure 1E). The data are wellapproximated by a gamma distribution, which has been postulated by others to be a good representation of celltocell variability in protein levels (Figure 1Egreen line) [22,2731]. We then stimulated cells with 0.1 and 1 nM EGF for 5 minutes and measured both ppERK and ERK levels simultaneously (Figures 1FI). Normalizing the ppERK levels by the amount of total ERK in each individual cell does not change the variance of “ERKoff” population (Figures 1FG—compare green to black lines). This is most likely because measurement variability is dominant at these low ppERK levels, and normalizing by total ERK levels does not correct for measurement variability. Normalizing the ppERK levels by total ERK levels does reduce the variability of the “ERKon” population, but does not change the fraction of cells in the “ERKon” and “ERKoff” populations (Figures 1FG). This assertion is reinforced by the fact that ppERK levels in both the “ERKoff” and “ERKon” populations span the entire spectrum of total ERK levels (Figures 1HI). Moreover, there is significant positive correlation between total ERK and ppERK levels in both the ERKoff and ERKon populations (Figures 1HI). Thus, although celltocell variability in ERK abundance contributes to ppERK response variability, it does not control bimodality, raising the question of what other factors contribute to the observed bimodality.
Stochastic, dynamic modeling explanation of the data
EGF activates the small GTPase Ras, which activates ERK downstream of the Raf and MEK kinases. Although we were not able to measure GTPbound active Ras (RasGTP) by flow cytometry, the population average dose and dynamic responses were assayed via pulldown and Western blotting, and then quantified (Figure 2A). These population average data show a rapid rise and dosedependent peak in RasGTP levels after EGF stimulation, followed by a fast decline. Although the most direct interpretation of these RasGTP responses (where the population mean changes as a function of time and EGF dose, Figure 2A) is a unimodal RasGTP distribution, a recent study suggested that in T lymphocytes, a positive feedback between RasGTP and its activator guanine exchange factor Son of Sevenless (SOS) leads to bistability and hysteresis in Ras activation [3234]. If Ras activation was also bistable in HEK293 cells, then two distinct RasGTP populations would exist with high mean and low mean RasGTP levels (Figure 2B). Stimulation by EGF would only affect the relative fraction of cells in the two populations, but not their means. Since under basal conditions ppERK levels are negligible (Additional file 1: Figure S1A), the low mean RasGTP population would not contribute to ERK activation, implying that there is a threshold above which RasGTP levels cause ERK activation (Figure 2B). If we assume a simple sigmoidal dose–response relationship between RasGTP and ppERK levels (typical in MAPK cascades—reviewed in [35]), then a defined high mean RasGTP population would induce a defined high mean ppERK population with boundaries E_{onlow} and E_{onhigh} (Figure 2B). However, the flow cytometry data in Figure 1AD show that when clear bimodality is present, E_{onlow} and E_{onhigh} are different for various high mean ppERK populations. Thus in HEK293 cells, our single cell ppERK signaling data seem to be inconsistent with a bistable RasGTP model.
Figure 2. RasGTP Dynamics. A. HEK293 cells were stimulated with 0.1 nM (low, L), 1 nM (medium, M) or 10 nM (high, H) EGF for the indicated times and then cell lysates were assayed for RasGTP as described in “Materials and Methods” section. IP denotes the pulldown fraction, TCL denotes total cell lysate, and IB denotes immunoblot. Each data point corresponds to the average of three independent experiments, and error bars correspond to 90% confidence intervals. Data were normalized by dividing by basal (no EGF) RasGTP levels. B. A bimodal RasGTP distribution as would be obtained from a bistable RasGTP model for low, medium and high levels of EGF, and how it would map onto a dose–response relationship between RasGTP and ppERK.
If the RasGTP response to EGF is unimodal, then how might these mixed analogdigital responses emerge from salient features of the MAPK/ERK cascade? At the single cell level, dynamic responses are encoded by the pathway topology and reaction kinetics. Therefore, we examined different configurations of the MAPK/ERK cascade for their ability to reproduce the experimentally observed behavior. Specifically, we sought topologies where simulations showed that (i) distributions of active ERK display bimodal/shouldering behavior with increasing EGF dose, and (ii) the “ERKon” population mean increases with increasing EGF dose at early time points, but decreases with time at constant EGF dose. To explore this, we used a previously developed mechanistic model that relates active Ras to ppERK [36], and investigated in silico the ability of different network topologies to reproduce our experimental observations (Figure 3A). By changing the feedback strength parameter (F_{a}) in this model, we created three different topologies: positive feedback (PF; F_{a} = 5), ultrasensitive (US; F_{a} = 1), and negative feedback (NF; F_{a} = 0.5), all of which have been experimentally observed for MAPK cascades under various circumstances (reviewed in [35]).
Figure 3. Modeling and analysis of single cell characteristics of the ppERK dynamics and dose response. A. Schematic of a mechanistic model of ERK activation and its steadystate response properties. The positive feedback (PF), nofeedback, ultrasensitive (US), and negative feedback (NF) models have the feedback strength F_{a} set to different values (5, 1 and 0.5, respectively). The input is RasGTP, and the output is ppERK. BD. Steadystate, deterministic input/output response curves for the PF, US and NF models. EG. Steadystate, cell population input/output response curves for the PF, US and NF models. In BG, red denotes increasing input from low levels, while blue denotes decreasing input from high levels. When only one color is shown, there is no difference between the increasing and decreasing input curves. Dashed lines indicate the 95^{th} percentile of all simulations, and dotted lines indicate the 5^{th} percentile.
Steadystate analysis
First, we characterized the steadystate input–output behavior of these three models by changing the input (RasGTP) from zero to 100 nM at 1 nM increments and allowing the system reach a steadystate between each step change. Then, we reversed the stimulation, this time changing the input from 100 to zero nM. The PF model exhibits bistability/hysteresis, whereas the US and NF models do not (Figure 3BD). In fact, due to the inherent properties of a negative feedback loop coupled with a kinase amplifier module, the NF model exhibits a smooth, analog input–output relationship [3740]. However, the NF model also exhibits a threshold of ERK activation at low RasGTP levels as a result of the multitier, multisite phosphorylation structure of the MAPK/ERK cascade [2].
These deterministic simulations correspond to input–output curves for an average cell. To incorporate stochastic fluctuations in reaction rates, we applied the Gillespie algorithm to integrate the differential equation. However, these solutions did not appreciably change the steadystate dose responses (data not shown), indicating that under these conditions and model parameters, reaction rate fluctuations do not constitute a significant source of signaling variability. This is most likely due to the relatively high abundance of the MAPK/ERK cascade components.
We therefore explicitly included protein expression variability in the models. We first investigated whether the gamma distribution provides a generally valid model for the distribution of protein levels, as others have suggested [22,2731]. We found that there is good agreement between gamma distribution fits and both experimental and stochastic simulation data from the literature (Additional file 1: Figure S3AE) [22]. Next, we performed our own stochastic simulations using a simple protein expression model where a gene can be active or inactive, an active gene can produce mRNA, mRNA can produce protein, and both mRNA and protein can degrade, all with first order kinetics. We then analyzed the resulting distribution of steadystate protein abundance obtained from multiple independent simulations under 6400 different parameter conditions (see Additional file 1: Figure S3 legend). For most conditions, the steadystate protein abundance distribution is well represented by a gamma distribution (Additional file 1: Figure S3FG). Therefore, for the steadystate analysis we sampled total levels of Raf, MEK and ERK from a gamma distribution, and computed the dose response curves for 1000 cells, each cell having different, sampled levels of Raf, MEK and ERK (Figures 3EG). The means of these stochastic, steadystate response curves (solid lines) have the same qualitative features as the deterministic curves, and the PF model remains bistable. However, there is substantial celltocell variability in the dose responses. The RasGTP levels eliciting halfmaximal ppERK responses vary significantly, as do the maximum ppERK levels. According to these results, stochastic variability in protein expression is a major contributor to steadystate, celltocell signaling variability, inducing a wide distribution of ERK activation thresholds.
Analysis of transient responses
To simulate the dynamic behavior of ppERK, we first specified the RasGTP input kinetics, according to the unimodal RasGTP distribution hypothesis discussed above. Experimental data show that in EGFstimulated HEK293 cells, RasGTP levels peak between 1–5 minutes after EGF stimulation and then, approximately 10 minutes later, decay to a steadystate value that is slightly higher than basal RasGTP levels (Figure 2A and [41]). Moreover, increasing the EGF dose increases the peak magnitude of RasGTP levels, and shortens the rise time. We incorporated these experimentally observed trends into a simple mathematical model (see Methods), and obtained simulated RasGTP dynamics. We then used these simulated dynamics as input to the MAPK/ERK cascade model for determining the ppERK dynamic and dose responses (Figure 4A). To incorporate celltocell variability in Ras levels, we sampled the peak RasGTP values from a gamma distribution whose mean increases with increasing input magnitude (with fixed shape parameter k—see Methods).
Figure 4. Simulations of ppERK dynamics in cell populations. A. Simulated RasGTP dynamics for different EGF doses. Simulations were done as described in the “Materials and Methods” section. BD. Simulated dose and dynamic ppERK responses for the PF (B), US (C) and NF (D) models. To facilitate comparison of these simulations with the experimental data, normally distributed noise with mean and standard deviation of 10 nM was added to the raw simulation data. (EG) Dynamics and dose response of the ERKon population mean for the PF (E), US (F), and NF (G) models. Simulated population responses were parsed into ERKon and ERKoff populations based on a cutoff of 100 nM.
Using these RasGTP dynamics, we then investigated which models (NF, PF or US) reproduce the experimental observations described above. As expected, the PF and US models show bimodal population behavior because of their switchlike input–output responses (Figure 4BC). But surprisingly, so too does the NF model, despite exhibiting an analog input–output relationship (Figure 4D). This bimodality in the NF model is due to the wide range of ERK activation thresholds introduced by protein expression variability (Figure 3G), combined with variability in EGFinduced RasGTP levels. Thus, all three topologies exhibit time and dosedependent bimodality or “shouldering”. However, only the NF model simulations, and not those of the US or PF models, reproduce proper behavior of the ERKon population mean, namely that the mean increases as a function of dose at short times (Figure 4EG; Figure 1), and decreases as a function of time at a particular EGF dose (Figure 4EG; Additional file 1: Figure S2).
We conclude that for the realistic parameter values used here, the NF model with protein expression variability is most consistent with experimental data. To examine if this conclusion holds over a wide range of parameter values, we employed parameter sensitivity analysis (see Methods and Additional file 1: Figure S4). This analysis showed that models with negative feedback preferentially demonstrated the experimentally observed signaling characteristics over the examined parameter ranges (Additional file 1: Figure S4). Yet, we cannot rule out the possibility that positive feedback and ultrasensitive systems may also exhibit the experimentally observed behavior. Indeed, sensitivity analysis also showed that under some rare parameter conditions, the mean ppERK levels in the ERKon population increase as a function of dose at short times for the PF and US models (Additional file 1: Figure S4A,C). One mechanism that may lead to this PF and US model behavior is if the ppERK activation kinetics were slow, such that the behavior at 2 and 5 min post EGF stimulation were due to transient effects, rather than a pseudosteady state phenomenon. Yet, for PF models, simulated ppERK signaling remains high over the 30minute time course (Additional file 1: Figure S4B,D), rather than returning closer to basal levels as the experimental data show (Additional file 1: Figure S2). Thus, the ERK cascade model with negative feedback seems to be the most consistent with our experimental observations over a wide range of parameter values.
Test of the negative feedback prediction
Although the preceding analysis suggests that in our HEK293 cell system the most likely net feedback strength from ERK is negative, parameter sensitivity analysis showed that ultrasensitive or positive feedback systems might also account for such data, albeit in rare circumstances. If the feedback were negative, blocking ERK activity should increase the activation of upstream elements, such as RasGTP. Therefore, we measured the dynamic and dose response of RasGTP with and without the MEK inhibitor U0126, and found that blocking ERK activation increased RasGTP levels, confirming the presence of strong negative feedback (Figure 5A). Although positive feedback and ultrasensitivity have been observed in various MAPK cascades (reviewed in [35]), in HEK293 cells the major feedback regulation is negative, confirming the predictions of the modeling. Notably, this feedback is less significant at five minutes after EGF stimulation, when the RasGTP response is saturated and ppERK levels are at their peak, implying that either this feedback is slow (which may introduce instability and oscillations under certain conditions [42]), or perhaps that there are alternative negative feedback mechanisms.
Figure 5. Confirming the presence of negative feedback. HEK293 cells were pretreated with 5 μM U0126 or vehicle alone (DMSO) for 30 min prior to stimulation with 0.1 nM, 1 nM or 10 nM of EGF (A) or TGFα (B) for 5 or 30 minutes. Control cells were left unstimulated (−). Total cell lysates were assayed for activated RasGTP as described in “Materials and Methods” section. IP denotes the pulldown fraction, TCL denotes total cell lysate, and IB denotes immunoblot.
To investigate whether alternative negative feedback mechanisms may explain the weak feedback effects at 5 minutes poststimulation, we repeated the U0126 experiment with the EGF receptor ligand TGFα. Although both EGF and TGFα activate the EGF receptor and induce receptor endocytosis, EGF preferentially targets the receptor to multivesicular bodies and lysosomal degradation, while TGFα enhances receptor recycling and surface availability [43,44]. Thus, it is possible that EGFinduced receptor degradation or sequestration may be influencing our results. We found that the TGFαinduced RasGTP levels do not differ from those induced by EGF in the presence or absence of the MEK inhibitor U0126 over a 30minute time course (Figure 5B). Therefore we conclude that negative feedback from ERK seems to dominate traffickingmediated effects.
Discussion
We have studied EGFinduced signal transduction to ERK in single HEK293 cells, finding that the conversion of an analog signal at the single cell level to an apparent digital response at the population level can be mediated by a combination of celltocell variability in protein expression and a pathway design that incorporates negative feedback (Figure 6). A uniform step increase in EGF concentration causes a wide distribution of RasGTP levels due to celltocell heterogeneity in protein expression. Celltocell heterogeneity in protein expression also causes significant variability in the sigmoidal dose response relationship between RasGTP and ppERK, and in particular, in the ppERK activation threshold (Figure 3G and Figure 6). Because celltocell variability in RasGTP levels can span the range of ERK pathway activation thresholds, the pathway is activated to various degrees in individual cells. A distribution of ppERK levels ensues across the cell population. The mean of the ppERK distribution depends on EGF dosage and agrees with results obtained from Western blots. Despite the fact that the negative feedback smooths the RasGTP/ppERK dose–response relationship, a threshold for ppERK activation persists. This threshold element further enhances celltocell variability in ppERK levels, and results in bimodal responses at the population level. Thus, the resulting bimodal distribution relies on a combination of a threshold behavior and a linear ppERK increase followed by saturation behavior with increasing EGF dose (Figure 3G). Surprisingly and counterintuitively, bimodality does not require switchlike behavior at the singlecell level, but can arise from celltocell variability in protein expression and a pathway activation threshold. Thus, cells can retain the robustness benefits offered by negative feedback [3740], while generating on/off responses at the cell population level that are thought to be critical for cell fate decisions.
Figure 6. Conversion of Analog Inputs Into Bimodal Responses by a Negative Feedback System Combined with Protein Expression Noise. An analog EGF stimulus (blue, green, and red correspond respectively to small, medium, and large stimulation magnitudes) induces variable but dosedependent RasGTP responses in the cell population due to expression variability in the EGF pathway proteins. RasGTP responses are converted into ppERK responses in single cells according to a thresholdlinear response governed by negative feedback (NF model). However, variability in RasGTP levels coupled with variability in ERK activation thresholds creates bimodal active ERK distributions at the population level despite the analog input and linear dose response at the single cell level.
Our observations are unlikely to be caused by a fraction of cells simply not binding ligand. First, under our experimental conditions (~10^{6} cells/mL), at the lowest ligand dose (0.01 nM), the ratio of EGF molecules to cells is approximately 1000, making it very unlikely that a cell does not encounter a ligand molecule. Second, for nearly all EGF doses, a significant fraction of cells is in the “ERKon” population at some point in time, indicating that most cells have been activated and therefore had bound ligand.
How might cells still generate reliable signals despite protein expression noise? One possibility is that cells have a reliable foldchange response of ppERK from basal levels, and that downstream of ppERK cells employ systems that sense foldchanges rather than absolute levels. In fact this foldchange scenario has recently been shown to be the case. In cells stably expressing ERK2YFP from the endogenous promoter, EGF stimulation led to widely varying maximum nuclear ERK2YFP accumulation, with a coefficient of variation (CV) of approximately 0.3 [15]. However, normalizing the maximum nuclear ERK2YFP signal by the basal levels of ERK2YFP in the same cell, which yields foldchange responses, lowers the CV by approximately 3fold [15]. This is consistent with our observed effects of total ERK abundance variability on the total variance of ppERK in the ERKon population (Figure 1FG). To sense these foldchanges, rather than absolute levels, a cell may use a type1 incoherent feedforward loop (I1FFL), where an input X activates both an intermediate Y and the output Z, but Y represses Z [45]. Such a network structure may in principle be downstream of ppERK (X), which causes the immediateearly expression of multiple genes including cfos, which can mediate general transcriptional repression perhaps even of itself [46,47].
Although protein expression noise is certainly a hindrance to some biological functions, and evolution has selected for mechanisms such as the I1FFL that allow a cell to deal with this noise, there are potential benefits of and perhaps even essential functions for such noise. Tissue homeostasis may in fact require protein expression variability. Consider that there is no protein expression variability, and all cells that are involved with, for instance, hematopoiesis, respond identically to the various proliferation and differentiation cues. The body needs to produce, from the hematopoietic stem cells, a balance between the lymphoid and myeloid progenitors. If all the hematopoietic stem cells responded identically, then it would be nearly impossible for the body to maintain a finely tuned balance between the production of these two lineages. The same logic applies to the further differentiation of lymphoid and myeloid progenitors into various other downstream cell types, such as megakaryocytes, erythrocytes, B cells, T cells, and natural killer cells, where finely tuned control of differential cellfate decisions is even more critical. Thus, it is likely that without protein expression noiseinduced phenotypic variability, homeostasis of hematopoiesis, and probably other tissues, would not be possible. This logic argues for a conceptual model whereby growth factor concentration, in tissues, controls the probability a cell will choose a particular fate.
Conclusions
It is commonly thought that the existence of bimodal signaling behavior on the population level is indicative of socalled digital behavior (such as allornone switches) of the underlying signaling network in single cells. Our work demonstrates that this is not necessarily the case; protein expression noise coupled with nonlinear network dynamics can bring about digital population responses from analog single cell dose responses. In particular, we show that a network combining an activation threshold and strong negative feedback also robustly displays such bimodal population behavior due to celltocell variability in protein expression levels. This system retains the benefits of robustness arising from negative feedback, while simultaneously generating populationlevel on/off responses thought to be critical for cell fate decisions. Overall, the results extend our understanding of the amazing behavioral complexity that can be displayed by even small molecular networks [48].
Methods
Cell culture
Human Embryonic Kidney 293 (HEK293) cells were obtained from the American Type Culture Collection (Manassas, VA). Cells were maintained in a humidified 5% CO_{2} incubator at 37°C and cultured in Dulbecco's modified Eagle's medium/F12 supplemented with 10% fetal bovine serum (Life TechnologiesInvitrogen, Carlsbad, CA) and penicillinstreptomycin solution (100 μg/ml, Thermo Fisher Scientific).
Flow cytometry
HEK293 cells were serum starved for 16 hours before the experiment. The cells were then lifted (by scraping or trypsinization), washed twice with serumfree medium (containing soybean trypsin inhibitor in the case of tryptic lifting), allowed to equilibrate for 30 minutes, and stimulated with EGF (SigmaAldrich, St. Louis, MO). We verified that the bimodal ppERK behavior was not affected by cell detachment (Additional file 1: Figure S5). After EGF stimulation for the desired time interval, cells were fixed with 2% paraformaldehyde (SigmaAldrich) for 10 minutes at 37°C, and then cooled on ice. After centrifugation, the cells were permeabilized in icecold 90% methanol (SigmaAldrich) for 30 minutes. The cells were then washed by centrifugation and 5x10^{5} cells were resuspended in 90 μL incubation/blocking buffer (0.5% BSA in PBS) for 10 minutes. The cells were then incubated for 60 minutes in the dark at room temperature with phosphoERK1/2 (T202/Y204) mouse mAb (E10) Alexa 488 Conjugate for active ERK and ERK1/2 rabbit mAb (4695) detected by secondary staining with an antirabbit Alexa 647conjugate (Cell Signaling Technologies, Beverley, MA). The cells were washed by centrifugation with PBS and resuspended in 0.5 mL of PBS. The samples were then analyzed with a BectonDickinson FACSCalibur or on an Accuri C6. For each sample, 10,000 events (cells) were analyzed. Data were processed using FlowJo^{™} software (Tree Star, Inc.) and MATLAB^{™} (The Mathworks). Postgating by forward and side scatter was performed to remove events corresponding to dead cells, debris, and cell clusters (i.e. doublets). As controls we stained cells with nonspecific, isotypematched control antibodies (also obtained from Cell Signaling). We verified the specificity of the antibodies (Additional file 1: Figure S1).
Western blotting
The above procedure for cell preparation was followed, but instead of fixing cells in paraformaldhyde, cells were lysed and processed for Western blotting analysis as described previously [49,50]. RasGTP pulldowns were performed as described previously [49,50].
Mechanistic model simulations
MATLAB and the function ode15s was used to simulate a previously developed, ordinary differential equationbased ERK cascade model [36], which is described in detail in Tables 1 and 2. The function gamrnd was used to generate realizations of peak RasGTP, Raf, MEK, and ERK levels for individual “cells” in the stochastic simulations according to the gamma distribution
where N specifies a protein level, k is the shape parameter, and θ is the scale parameter. We specified the k (shape) parameter of each gamma distribution as 5.4, as was measured for total ERK (see Figure 1E), assuming roughly similar expression regulation. Since the mean of a gamma distribution is equal to kθ, the θ parameter of each gamma distribution was changed as needed to attain the desired distribution mean (see Table 1 for values of mean protein levels).
Table 1. Kinetic description of the ERK signaling cascade
Table 2. Ordinary differential equations for the ERK signaling cascade model
To estimate the parameters for the RasGTP dynamics, which are described by a simple exponential rise and decay model (see Table 2 for differential equations), we used least squares optimization to ensure that desired initial magnitude (I_{o}), peak magnitude (I_{max}), timetopeak (τ_{max}), timetoinflection (τ_{infl}), timetosteadystate (τ_{ss}), and steadystate magnitude (I_{ss}) of the RasGTP dynamics matches well to that which the model prescribes. Additional file 1: Figure S6 describes these RasGTP dynamics metrics graphically. As there are four unknown parameters in the RasGTP dynamics model (Table 2K_{1}, K_{2}, τ_{1}, τ_{2}), we need four equations, which we take as the following (their origin is described immediately below):
where w_{i} corresponds to a weight for optimization purposes (all w’s are 1 except for w_{2} which is 100). Eq. 2 specifies the proper steadystate magnitude; Eq. 3 specifies that the 1^{st} derivative at the timetopeak is zero; Eq. 4 specifies the proper magnitude at the timetosteady state (defined as 1% of the true steadystate value—see Additional file 1: Figure S6); and Eq. 5 specifies the proper peak magnitude. The following constraints are placed on this optimization problem:
Eq.6 specifies that there is a maximum at the timetopeak (2^{nd} derivative less than zero) and Eq. 7 specifies that the 1^{st} derivative is negative at the inflection point (RasGTP is decreasing towards the steadystate value). Mean peak RasGTP levels (I_{max}) were increased to simulate increasing input, and were linearly spaced between 10 nM and 200 nM using 6 points (10, 48, 86, 124, 162, and 200 nM), which correspond to EGF doses (in nM) of 0.01, 0.1, 0.5, 1, 5, and 10. Following the trends of the experimental data in Additional file 1: Figure 2A and [41], peak times for RasGTP (τ_{max}) were sampled linearly between 7 min and 2 min (7, 6, 5, 4, 3, 2), with 7 min corresponding to the lowest peak RasGTP level (EGF dose). Also, we took τ_{ss} as 10 min, I_{ss} as 15% of I_{max} realizations, I_{o} as 0, and τ_{infl} as (τ_{max} + τ_{ss})/2.
All code is available upon request.
Parameter sensitivity analysis
Five hundred different parameter sets were generated via latin hypercube sampling (MATLAB function lhsdesign) from a 23dimensional uniform distribution that spans +/− 1 order of magnitude around each nominal parameter value (taken from Table 1 with the exception of F_{a}—the feedback strength). For each of these parameter sets stochastic simulations were performed as described above. Briefly, total protein and RasGTP levels were sampled from a gamma distribution and 500 individual cell responses were simulated for each parameter set and feedback condition (negative, ultrasensitive with no feedback, positive). The results of these simulations were then analyzed for three features: the “analogicity” of the ERKon population, the “transience” of the ERKon population, and bimodality. The analogicity of a particular feedback/parameter set combination was calculated as follows, and is illustrated in Additional file 1: Figure S4A. First, the ERKon population was defined by those cells having ppERK levels over 200 nM. Then, the mean ppERK levels in the ERKon populations were calculated for those that contained greater than 10 cells. The analogicity of a given time point is defined as the maximum ERKon population mean minus the minimum (as compared across EGF doses). The analogicity of a feedback/parameter set combination is the sum of the 2 and 5 minute time point analogicities. The 10 and 30 minute time points are left out because these show very little analogicity in the experimental data (Figure 1 and Additional file 1: Figure S2). Parameter sets showing zero analogicity were discarded as inconsistent with experimental data. The transience of a particular feedback/parameter set combination is defined for a particular EGF dose as follows, and is pictorially illustrated in Additional file 1: Figure S4B. First, the ERKon population was defined as described above for analogicity, and any EGF dose where the ERKon population did not exist for all time points was not used for further transience calculations. The transience of an individual EGF dose is the mean of the ERKon population at 2 and 5 minutes minus that at 10 and 30 min. The transience of a feedback/parameter set combination is the sum over those from the individual EGF doses. Bimodality was evaluated via Hartigan’s Dip Test [51,52]. MATLAB code for this test was downloaded from http://www.nicprice.net/diptest/ webcite. The result is a pvalue associated with the hypothesis test that the empirical distribution of interest is unimodal as opposed to the alternative that it is not. We rejected the null hypothesis at the 0.05 level of significance. The bimodal fraction for a particular feedback/parameter set combination is defined as the number of nonunimodal distributions divided by the total number of dose/time point combinations. Parameter sets showing no bimodality were discarded as inconsistent with experimental data.
Competing interests
The authors declared that they have no competing interests.
Authors’ contributions
MB and JR performed and designed research and wrote the paper. MD, AK, and EA performed research. JH designed research. WK wrote the paper. BO and BK designed research and wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
We thank Claudio Gelmi and Erik Welf for helpful discussions. This work was supported by Science Foundation Ireland [06/CE/B1129]. MRB acknowledges a Marie Curie International Incoming Fellowship [236758] and an EMBO longterm fellowship [ALTF 815–2010].
References

Goldbeter A, Koshland DE Jr: An amplified sensitivity arising from covalent modification in biological systems.
Proc Natl Acad Sci U S A 1981, 78(11):68406844. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Markevich NI, Hoek JB, Kholodenko BN: Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades.
J Cell Biol 2004, 164(3):353359. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ferrell JE Jr: Selfperpetuating states in signal transduction: positive feedback, doublenegative feedback and bistability.
Curr Opin Cell Biol 2002, 14(2):140148. PubMed Abstract  Publisher Full Text

Pomerening JR, Sontag ED, Ferrell JE Jr: Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2.
Nat Cell Biol 2003, 5(4):346351. PubMed Abstract  Publisher Full Text

Sha W, Moore J, Chen K, Lassaletta AD, Yi CS, Tyson JJ, Sible JC: Hysteresis drives cellcycle transitions in Xenopus laevis egg extracts.
Proc Natl Acad Sci U S A 2003, 100(3):975980. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tyson JJ, CsikaszNagy A, Novak B: The dynamics of cell cycle regulation.
Bioessays 2002, 24(12):10951109. PubMed Abstract  Publisher Full Text

Huang CY, Ferrell JE Jr: Ultrasensitivity in the mitogenactivated protein kinase cascade.
Proc Natl Acad Sci U S A 1996, 93(19):1007810083. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Santos SD, Verveer PJ, Bastiaens PI: Growth factorinduced MAPK network topology shapes Erk response determining PC12 cell fate.
Nat Cell Biol 2007, 9(3):324330. PubMed Abstract  Publisher Full Text

Mackeigan JP, Murphy LO, Dimitri CA, Blenis J: Graded mitogenactivated protein kinase activity precedes switchlike cFos induction in mammalian cells.
Mol Cell Biol 2005, 25(11):46764682. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Harding A, Tian T, Westbury E, Frische E, Hancock JF: Subcellular localization determines MAP kinase signal output.
Curr Biol 2005, 15(9):869873. PubMed Abstract  Publisher Full Text

Inder K, Harding A, Plowman SJ, Philips MR, Parton RG, Hancock JF: Activation of the MAPK module from different spatial locations generates distinct system outputs.
Mol Biol Cell 2008, 19(11):47764784. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ferrell JE Jr, Machleder EM: The biochemical basis of an allornone cell fate switch in Xenopus oocytes.
Science 1998, 280(5365):895898. PubMed Abstract  Publisher Full Text

McAdams HH, Arkin A: Stochastic mechanisms in gene expression.
Proc Natl Acad Sci U S A 1997, 94(3):814819. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

BarEven A, Paulsson J, Maheshri N, Carmi M, O'Shea E, Pilpel Y, Barkai N: Noise in protein expression scales with natural protein abundance.
Nat Genet 2006, 38(6):636643. PubMed Abstract  Publisher Full Text

CohenSaidon C, Cohen AA, Sigal A, Liron Y, Alon U: Dynamics and variability of ERK2 response to EGF in individual living cells.
Mol Cell 2009, 36(5):885893. PubMed Abstract  Publisher Full Text

Niepel M, Spencer SL, Sorger PK: Nongenetic celltocell variability and the consEquationuences for pharmacology.
Curr Opin Chem Biol 2009, 13(5–6):556561. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK: Nongenetic origins of celltocell variability in TRAILinduced apoptosis.
Nature 2009, 459(7245):428432. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S: Transcriptomewide noise controls lineage choice in mammalian progenitor cells.
Nature 2008, 453(7194):544547. PubMed Abstract  Publisher Full Text

Raser JM, O'Shea EK: Noise in gene expression: origins, consEquationuences, and control.
Science 2005, 309(5743):20102013. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems.
Nat Rev Genet 2009, 10(2):122133. PubMed Abstract  Publisher Full Text

Pedraza JM, Paulsson J: Effects of molecular memory and bursting on fluctuations in gene expression.
Science 2008, 319(5861):339343. PubMed Abstract  Publisher Full Text

Raj A, Peskin CS, Tranchina D, Vargas DY, Tyagi S: Stochastic mRNA synthesis in mammalian cells.
PLoS Biol 2006, 4(10):e309. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yu J, Xiao J, Ren X, Lao K, Xie XS: Probing gene expression in live cells, one protein molecule at a time.
Science 2006, 311(5767):16001603. PubMed Abstract  Publisher Full Text

Cai L, Friedman N, Xie XS: Stochastic protein expression in individual cells at the single molecule level.
Nature 2006, 440(7082):358362. PubMed Abstract  Publisher Full Text

ColmanLerner A, Gordon A, Serra E, Chin T, Resnekov O, Endy D, Pesce CG, Brent R: Regulated celltocell variation in a cellfate decision system.
Nature 2005, 437(7059):699706. PubMed Abstract  Publisher Full Text

Perez OD, Nolan GP: Phosphoproteomic immune analysis by flow cytometry: from mechanism to translational medicine at the singlecell level.
Immunol Rev 2006, 210:208228. PubMed Abstract  Publisher Full Text

Taniguchi Y, Choi PJ, Li GW, Chen H, Babu M, Hearn J, Emili A, Xie XS: Quantifying E. coli proteome and transcriptome with singlemolecule sensitivity in single cells.
Science 2010, 329(5991):533538. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Shahrezaei V, Swain PS: Analytical distributions for stochastic gene expression.
Proc Natl Acad Sci U S A 2008, 105(45):1725617261. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cohen AA, Kalisky T, Mayo A, GevaZatorsky N, Danon T, Issaeva I, Kopito RB, Perzov N, Milo R, Sigal A, et al.: Protein dynamics in individual human cells: experiment and theory.
PLoS One 2009, 4(4):e4901. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Paulsson J, Berg OG, Ehrenberg M: Stochastic focusing: fluctuationenhanced sensitivity of intracellular regulation.
Proc Natl Acad Sci U S A 2000, 97(13):71487153. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Friedman N, Cai L, Xie XS: Linking stochastic dynamics to population distribution: an analytical framework of gene expression.
Phys Rev Lett 2006, 97(16):168302. PubMed Abstract  Publisher Full Text

Boykevisch S, Zhao C, Sondermann H, Philippidou P, Halegoua S, Kuriyan J, BarSagi D: Regulation of ras signaling dynamics by Sosmediated positive feedback.
Curr Biol 2006, 16(21):21732179. PubMed Abstract  Publisher Full Text

Das J, Ho M, Zikherman J, Govern C, Yang M, Weiss A, Chakraborty AK, Roose JP: Digital signaling and hysteresis characterize ras activation in lymphoid cells.
Cell 2009, 136(2):337351. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Prasad A, Zikherman J, Das J, Roose JP, Weiss A, Chakraborty AK: Origin of the sharp boundary that discriminates positive and negative selection of thymocytes.
Proc Natl Acad Sci U S A 2009, 106(2):528533. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kholodenko BN, Birtwistle MR: Fourdimensional dynamics of MAPK information processing systems.
Wiley Interdiscip Rev Syst Biol Med 2009, 1(1):2844. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Markevich NI, Tsyganov MA, Hoek JB, Kholodenko BN: Longrange signaling by phosphoprotein waves arising from bistability in protein kinase cascades.
Mol Syst Biol 2006, 2:61. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sauro HM, Kholodenko BN: Quantitative analysis of signaling networks.
Prog Biophys Mol Biol 2004, 86(1):543. PubMed Abstract  Publisher Full Text

Birtwistle MR, Kolch W: Biology using engineering tools: the negative feedback amplifier.
Cell Cycle 2011, 10(13):20692076. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sturm OE, Orton R, Grindlay J, Birtwistle M, Vyshemirsky V, Gilbert D, Calder M, Pitt A, Kholodenko B, Kolch W: The mammalian MAPK/ERK pathway exhibits properties of a negative feedback amplifier.
Sci Signal 2010, 3(153):ra90. PubMed Abstract  Publisher Full Text

FritscheGuenther R, Witzel F, Sieber A, Herr R, Schmidt N, Braun S, Brummer T, Sers C, Bluthgen N: Strong negative feedback from Erk to Raf confers robustness to MAPK signalling.
Mol Syst Biol 2011, 7:489. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Borisov N, Aksamitiene E, Kiyatkin A, Legewie S, Berkhout J, Maiwald T, Kaimachnikov NP, Timmer J, Hoek JB, Kholodenko BN: Systemslevel interactions between insulinEGF networks amplify mitogenic signaling.
Mol Syst Biol 2009, 5:256. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kholodenko BN: Negative feedback and ultrasensitivity can bring about oscillations in the mitogenactivated protein kinase cascades.
Eur J Biochem 2000, 267(6):15831588. PubMed Abstract  Publisher Full Text

French AR, Tadaki DK, Niyogi SK, Lauffenburger DA: Intracellular trafficking of epidermal growth factor family ligands is directly influenced by the pH sensitivity of the receptor/ligand interaction.
J Biol Chem 1995, 270(9):43344340. PubMed Abstract  Publisher Full Text

Roepstorff K, Grandal MV, Henriksen L, Knudsen SL, Lerdrup M, Grovdal L, Willumsen BM, van Deurs B: Differential effects of EGFR ligands on endocytic sorting of the receptor.
Traffic 2009, 10(8):11151127. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Goentoro L, Shoval O, Kirschner MW, Alon U: The incoherent feedforward loop can provide foldchange detection in gene regulation.
Mol Cell 2009, 36(5):894899. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Nakakuki T, Birtwistle MR, Saeki Y, Yumoto N, Ide K, Nagashima T, Brusch L, Ogunnaike BA, OkadaHatakeyama M, Kholodenko BN: Ligandspecific cFos expression emerges from the spatiotemporal control of ErbB network dynamics.
Cell 2010, 141(5):884896. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gius D, Cao XM, Rauscher FJ 3rd, Cohen DR, Curran T, Sukhatme VP: Transcriptional activation and repression by Fos are independent functions: the C terminus represses immediateearly gene expression via CArG elements.
Mol Cell Biol 1990, 10(8):42434255. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kholodenko BN: Cellsignalling dynamics in time and space.
Nat Rev Mol Cell Biol 2006, 7:165176.
PMID: 16482094
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
Kiyatkin A, Aksamitiene E, Markevich NI, Borisov NM, Hoek JB, Kholodenko BN: Scaffolding protein Grb2associated binder 1 sustains epidermal growth factorinduced mitogenic and survival signaling by multiple positive feedback loops.
J Biol Chem 2006, 281(29):1992519938. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Aksamitiene E, Achanta S, Kolch W, Kholodenko BN, Hoek JB, Kiyatkin A: Prolactinstimulated activation of ERK1/2 mitogenactivated protein kinases is controlled by PI3kinase/Rac/PAK signaling pathway in breast cancer cells.
Cell Signal 2011, 23(11):17941805. PubMed Abstract  Publisher Full Text

Hartigan JA, Hartigan PM: The Dip Test of Unimodality.
Ann Stat 1985, 13(1):7084. Publisher Full Text

Hartigan PM: Computation of the Dip Statistic to Test for Unimodality.