Emergence of bimodal cell population responses from the interplay between analog single-cell signaling and protein expression noise

Abstract

Background

Cell-to-cell 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 cell-to-cell protein abundance variations interact to shape signaling responses.

Results

We observe bimodal distributions of extracellular signal-regulated 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/ERK-cascade 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 cell-to-cell 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 population-level on/off responses that are thought to be critical for regulating cell fate decisions.

Background

Development, growth and homeostasis of multi-cellular organisms depend on the ability of individual cells to convert noisy, analog signals into clear, yes-or-no 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 [1–3]. For example, various mitogen concentrations can cause bistable activation of cyclin-dependent kinases to drive cell cycle transition decisions [4–6]. Theoretical studies have shown that signaling networks containing positive or double negative feedback loops [3], opposing modification enzymes exhibiting saturation kinetics [1] and multi-site 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 signal-regulated kinase 1/2 (ERK) contains many elements that potentially can lead to switch-like behavior. However, previous single cell studies in different mammalian cell lines have reported both graded [8, 9] and “all-or-nothing” [10] EGF-induced ERK activation responses. One determinant of whether signaling is graded or switch-like is the spatial localization of signal processing proteins [11].

Under idealized conditions of cell-to-cell homogeneity, experimental techniques such as immunoblotting that measure average population responses may be able to detect all-or-none signaling responses, as long as the cell-to-cell 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 cell-to-cell variability in protein levels [13–20]. The primary stochastic factors are (i) the rate of transcription, which is burst-like 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 “all-or-nothing” responses; single-cell measurement techniques capable of capturing the dynamics of digital signal transduction are needed [12].

Here, we use flow cytometry to measure EGF-induced, single-cell ERK activation responses in a HEK293 cell population. We observe bimodal response distributions in cell populations that are usually thought to indicate switch-like 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 feedback-mediated, analog single-cell responses.

Results

Analyses of ERK responses to EGF in individual cells and populations

We used a flow cytometry-based 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 1A-D; Additional file 1: Figure S2). The increase in mean values of ppERK was dose-dependent 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 dose-dependent fraction of cells failing to activate ERK. At five minutes after EGF stimulation, the ppERK distribution is unambiguously bimodal, implying digital “on-off” behavior. Higher EGF doses increase the fraction of cells with high ppERK (”ERK-on”) at the expense of the “ERK-off” population. Thus, in a dose-dependent 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 EGF-induced ERK response on the population level is complex consisting of both analog and digital elements.

Figure 1

Cell population dose and dynamic response of ppERK to EGF. A-D. 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 x-axis is the magnitude of activated ERK (ppERK) in arbitrary fluorescence units, and the y-axis is the frequency of observing a particular level of fluorescence in a cell. E. Total ERK abundance data. The best-fit 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⁵ [AU]. F-I. 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 F-G, 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 H-I dot plots, ppERK levels are on the x-axis, whereas total ERK levels are on the y-axis.

Next, we investigated how cell-to-cell variability in total ERK abundance affects the ppERK responses. Measurements of the total ERK distribution by flow cytometry, as expected, revealed substantial cell-to-cell variability in total ERK levels (Figure 1E). The data are well-approximated by a gamma distribution, which has been postulated by others to be a good representation of cell-to-cell variability in protein levels (Figure 1E-green line) [22, 27–31]. We then stimulated cells with 0.1 and 1 nM EGF for 5 minutes and measured both ppERK and ERK levels simultaneously (Figures 1F-I). Normalizing the ppERK levels by the amount of total ERK in each individual cell does not change the variance of “ERK-off” population (Figures 1F-G—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 “ERK-on” population, but does not change the fraction of cells in the “ERK-on” and “ERK-off” populations (Figures 1F-G). This assertion is reinforced by the fact that ppERK levels in both the “ERK-off” and “ERK-on” populations span the entire spectrum of total ERK levels (Figures 1H-I). Moreover, there is significant positive correlation between total ERK and ppERK levels in both the ERK-off and ERK-on populations (Figures 1H-I). Thus, although cell-to-cell 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 GTP-bound active Ras (RasGTP) by flow cytometry, the population average dose and dynamic responses were assayed via pull-down and Western blotting, and then quantified (Figure 2A). These population average data show a rapid rise and dose-dependent 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 [32–34]. 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 _on-low and E _on-high (Figure 2B). However, the flow cytometry data in Figure 1A-D show that when clear bimodality is present, E _on-low and E _on-high 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 pull-down 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 analog-digital 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 “ERK-on” 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 steady-state response properties. The positive feedback (PF), no-feedback, 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. B-D. Steady-state, deterministic input/output response curves for the PF, US and NF models. E-G. Steady-state, cell population input/output response curves for the PF, US and NF models. In B-G, 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.

Steady-state analysis

First, we characterized the steady-state 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 steady-state 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 3B-D). 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 [37–40]. However, the NF model also exhibits a threshold of ERK activation at low RasGTP levels as a result of the multi-tier, multi-site 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 steady-state 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, 27–31]. 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 S3A-E) [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 steady-state protein abundance obtained from multiple independent simulations under 6400 different parameter conditions (see Additional file 1: Figure S3 legend). For most conditions, the steady-state protein abundance distribution is well represented by a gamma distribution (Additional file 1: Figure S3F-G). Therefore, for the steady-state 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 3E-G). The means of these stochastic, steady-state response curves (solid lines) have the same qualitative features as the deterministic curves, and the PF model remains bistable. However, there is substantial cell-to-cell variability in the dose responses. The RasGTP levels eliciting half-maximal ppERK responses vary significantly, as do the maximum ppERK levels. According to these results, stochastic variability in protein expression is a major contributor to steady-state, cell-to-cell 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 EGF-stimulated HEK293 cells, RasGTP levels peak between 1–5 minutes after EGF stimulation and then, approximately 10 minutes later, decay to a steady-state 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 cell-to-cell 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. B-D. 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. (E-G) Dynamics and dose response of the ERK-on population mean for the PF (E), US (F), and NF (G) models. Simulated population responses were parsed into ERK-on and ERK-off 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 switch-like input–output responses (Figure 4B-C). 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 EGF-induced RasGTP levels. Thus, all three topologies exhibit time and dose-dependent bimodality or “shouldering”. However, only the NF model simulations, and not those of the US or PF models, reproduce proper behavior of the ERK-on population mean, namely that the mean increases as a function of dose at short times (Figure 4E-G; Figure 1), and decreases as a function of time at a particular EGF dose (Figure 4E-G; 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 ERK-on 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 pseudo-steady state phenomenon. Yet, for PF models, simulated ppERK signaling remains high over the 30-minute 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 pull-down 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 post-stimulation, 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 multi-vesicular bodies and lysosomal degradation, while TGFα enhances receptor recycling and surface availability [43, 44]. Thus, it is possible that EGF-induced 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 30-minute time course (Figure 5B). Therefore we conclude that negative feedback from ERK seems to dominate trafficking-mediated effects.

Discussion

We have studied EGF-induced 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 cell-to-cell 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 cell-to-cell heterogeneity in protein expression. Cell-to-cell 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 cell-to-cell 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 cell-to-cell 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 switch-like behavior at the single-cell level, but can arise from cell-to-cell variability in protein expression and a pathway activation threshold. Thus, cells can retain the robustness benefits offered by negative feedback [37–40], 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 dose-dependent 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 threshold-linear 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⁶ 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 “ERK-on” 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 fold-change response of ppERK from basal levels, and that downstream of ppERK cells employ systems that sense fold-changes rather than absolute levels. In fact this fold-change scenario has recently been shown to be the case. In cells stably expressing ERK2-YFP from the endogenous promoter, EGF stimulation led to widely varying maximum nuclear ERK2-YFP accumulation, with a coefficient of variation (CV) of approximately 0.3 [15]. However, normalizing the maximum nuclear ERK2-YFP signal by the basal levels of ERK2-YFP in the same cell, which yields fold-change responses, lowers the CV by approximately 3-fold [15]. This is consistent with our observed effects of total ERK abundance variability on the total variance of ppERK in the ERK-on population (Figure 1F-G). To sense these fold-changes, rather than absolute levels, a cell may use a type-1 incoherent feedforward loop (I1-FFL), 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 immediate-early expression of multiple genes including c-fos, 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 I1-FFL 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 cell-fate decisions is even more critical. Thus, it is likely that without protein expression noise-induced 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 so-called digital behavior (such as all-or-none 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 cell-to-cell variability in protein expression levels. This system retains the benefits of robustness arising from negative feedback, while simultaneously generating population-level 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₂ incubator at 37°C and cultured in Dulbecco's modified Eagle's medium/F-12 supplemented with 10% fetal bovine serum (Life Technologies-Invitrogen, Carlsbad, CA) and penicillin-streptomycin 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 serum-free medium (containing soybean trypsin inhibitor in the case of tryptic lifting), allowed to equilibrate for 30 minutes, and stimulated with EGF (Sigma-Aldrich, 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 (Sigma-Aldrich) for 10 minutes at 37°C, and then cooled on ice. After centrifugation, the cells were permeabilized in ice-cold 90% methanol (Sigma-Aldrich) for 30 minutes. The cells were then washed by centrifugation and 5x10⁵ 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 phospho-ERK1/2 (T202/Y204) mouse mAb (E10) Alexa 488 Conjugate for active ERK and ERK1/2 rabbit mAb (4695) detected by secondary staining with an anti-rabbit Alexa 647-conjugate (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 Becton-Dickinson 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). Post-gating 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 non-specific, isotype-matched 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 pull-downs were performed as described previously [49, 50].

Mechanistic model simulations

MATLAB and the function ode15s was used to simulate a previously developed, ordinary differential equation-based 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

$$f\left(N\right)=\frac{N^{\left(k-1\right)}{e}^{-\frac{N}{\theta }}}{{\theta }^k\Gamma \left(k\right)}$$

(1)

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 ), time-to-peak (τ _max ), time-to-inflection (τ _infl ), time-to-steady-state (τ _ss ), and steady-state 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 2-K₁, K₂, τ₁, τ₂), we need four equations, which we take as the following (their origin is described immediately below):

$$0=\left(I_o+\left(K_1+K_2\right)-I_{\mathit{ss}}\right)w_1$$

(2)

$$0=\left(\frac{K_1}{{\tau }_1}exp\left(-\frac{{\tau }_{max}}{{\tau }_1}\right)+\frac{K_2}{{\tau }_2}exp\left(-\frac{{\tau }_{max}}{{\tau }_2}\right)\right)w_2$$

(3)

$$\begin{array}{ll}0=& \left(I_o+K_1\left(1-exp\left(-\frac{{\tau }_{\mathit{ss}}}{{\tau }_1}\right)\right)\right.\\ \left.+,K_2,\left(1-exp\left(-\frac{{\tau }_{\mathit{ss}}}{{\tau }_2}\right)\right),-,1.01,·,I_{\mathit{ss}}\right)w_3\end{array}$$

(4)

$$\begin{array}{ll}0=& \left(I_o+K_1\left(1-exp\left(-\frac{{\tau }_{max}}{{\tau }_1}\right)\right)\right.\\ \left.+,K_2,\left(1-exp\left(-\frac{{\tau }_{max}}{{\tau }_2}\right)\right),-,I_{max}\right)w_4\end{array}$$

(5)

where w _i corresponds to a weight for optimization purposes (all w’s are 1 except for w ₂ which is 100). Eq. 2 specifies the proper steady-state magnitude; Eq. 3 specifies that the 1^st derivative at the time-to-peak is zero; Eq. 4 specifies the proper magnitude at the time-to-steady state (defined as 1% of the true steady-state value—see Additional file 1: Figure S6); and Eq. 5 specifies the proper peak magnitude. The following constraints are placed on this optimization problem:

$$0>-\frac{K_1}{{{\tau }_1}^2}exp\left(-\frac{{\tau }_{max}}{{\tau }_1}\right)-\frac{K_2}{{{\tau }_2}^2}exp\left(-\frac{{\tau }_{max}}{{\tau }_2}\right)$$

(6)

$$0>\frac{K_1}{{\tau }_1}exp\left(-\frac{{\tau }_{\mathit{infl}}}{{\tau }_1}\right)+\frac{K_2}{{\tau }_2}exp\left(-\frac{{\tau }_{\mathit{infl}}}{{\tau }_2}\right)$$

(7)

Eq.6 specifies that there is a maximum at the time-to-peak (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 steady-state 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 23-dimensional 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 ERK-on population, the “transience” of the ERK-on 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 ERK-on population was defined by those cells having ppERK levels over 200 nM. Then, the mean ppERK levels in the ERK-on populations were calculated for those that contained greater than 10 cells. The analogicity of a given time point is defined as the maximum ERK-on 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 ERK-on population was defined as described above for analogicity, and any EGF dose where the ERK-on 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 ERK-on 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/. The result is a p-value 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 non-unimodal distributions divided by the total number of dose/time point combinations. Parameter sets showing no bimodality were discarded as inconsistent with experimental data.