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

doi:10.1186/1752-0509-6-109

BMC Syst Biol. 2012; 6: 109.

Published online 2012 August 24. doi: 10.1186/1752-0509-6-109

PMCID: PMC3484110

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

Marc R Birtwistle,^1,^2,^3,⁴ Jens Rauch,² Anatoly Kiyatkin,^#³ Edita Aksamitiene,^#³ Maciej Dobrzyński,^#² Jan B Hoek,³ Walter Kolch,² Babatunde A Ogunnaike,⁴ and Boris N Kholodenko^2,³

¹Dept. of Pharmacology and Systems Therapeutics, Mount Sinai School of Medicine, New York, NY, 10029, USA

²Systems Biology Ireland, University College Dublin, Belfield, Dublin 4, Ireland

³Dept. of Anatomy, Pathology, and Cell Biology, Thomas Jefferson University, Philadelphia, PA, 19107, USA

⁴Dept. of Chemical and Biomolecular Engineering, University of Delaware, Newark, DE, 19716, USA

Corresponding author.

^#Contributed equally.

Marc R Birtwistle: marc.birtwistle/at/gmail.com; Jens Rauch: jens.rauch/at/ucd.ie; Anatoly Kiyatkin: anatoly.kiyatkin/at/jefferson.edu; Edita Aksamitiene: exa140/at/jefferson.edu; Maciej Dobrzyński: maciej.dobrzynski/at/ucd.ie; Jan B Hoek: jan.hoek/at/jefferson.edu; Walter Kolch: walter.kolch/at/ucd.ie; Babatunde A Ogunnaike: ogunnaike/at/udel.edu; Boris N Kholodenko: boris.kholodenko/at/ucd.ie

Received May 7, 2012; Accepted July 27, 2012.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

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 (Figure1A-D;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 (more ...)

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 (Figure1E).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 (Figure1E-green1E-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 (Figures1F-I).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 (Figures1F-G—compare1F-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 (Figures1F-G).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 (Figures1H-I).1H-I). Moreover, there is significant positive correlation between total ERK and ppERK levels in both the ERK-off and ERK-on populations (Figures (Figures1H-I).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 (Figure2A).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 Figure2A)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 (Figure2B).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 (Figure2B).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 (Figure2B).2B). However, the flow cytometry data in Figure Figure1A-D1A-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 (more ...)

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 (Figure3A).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), (more ...)

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 (Figure3B-D).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 (Figures3E-G).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 (Figure2A2A 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 (Figure4A).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 (more ...)

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 (Figure4B-C).4B-C). But surprisingly, so too does the NF model, despite exhibiting an analog input–output relationship (Figure (Figure4D).4D). This bimodality in the NF model is due to the wide range of ERK activation thresholds introduced by protein expression variability (Figure (Figure3G),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 (Figure4E-G;4E-G; Figure Figure1),1), and decreases as a function of time at a particular EGF dose (Figure (Figure4E-G;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 (Figure5A).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 (more ...)

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 (Figure5B).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 (Figure6).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 (Figure3G3G and Figure Figure6).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 (Figure3G).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 (more ...)

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 (Figure1F-G).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 Tables11 and and2.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

Table 1

Kinetic description of the ERK signaling cascade

Table 2

Ordinary differential equations for the ERK signaling cascade model

(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 Figure1E),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 Table11 for values of mean protein levels).

To estimate the parameters for the RasGTP dynamics, which are described by a simple exponential rise and decay model (see Table Table22 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 (Table22-K₁, K₂, τ₁, τ₂), we need four equations, which we take as the following (their origin is described immediately below):

(2)

(3)

(4)

(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:

(6)

(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 Figure2A2A 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 Table11 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 (Figure11 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.

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.

Supplementary Material

Additional file 1

This additional file contains all the supplementary figures along with their legends.

Click here for file^{(1.4M, ppt)}

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 long-term 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):6840–6844. doi: 10.1073/pnas.78.11.6840. [PubMed] [Cross Ref]
Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004;164(3):353–359. doi: 10.1083/jcb.200308060. [PMC free article] [PubMed] [Cross Ref]
Ferrell JE Jr. Self-perpetuating states in signal transduction: positive feedback, double-negative feedback and bistability. Curr Opin Cell Biol. 2002;14(2):140–148. doi: 10.1016/S0955-0674(02)00314-9. [PubMed] [Cross Ref]
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):346–351. doi: 10.1038/ncb954. [PubMed] [Cross Ref]
Sha W, Moore J, Chen K, Lassaletta AD, Yi CS, Tyson JJ, Sible JC. Hysteresis drives cell-cycle transitions in Xenopus laevis egg extracts. Proc Natl Acad Sci U S A. 2003;100(3):975–980. doi: 10.1073/pnas.0235349100. [PubMed] [Cross Ref]
Tyson JJ, Csikasz-Nagy A, Novak B. The dynamics of cell cycle regulation. Bioessays. 2002;24(12):1095–1109. doi: 10.1002/bies.10191. [PubMed] [Cross Ref]
Huang CY, Ferrell JE Jr. Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc Natl Acad Sci U S A. 1996;93(19):10078–10083. doi: 10.1073/pnas.93.19.10078. [PubMed] [Cross Ref]
Santos SD, Verveer PJ, Bastiaens PI. Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate. Nat Cell Biol. 2007;9(3):324–330. doi: 10.1038/ncb1543. [PubMed] [Cross Ref]
Mackeigan JP, Murphy LO, Dimitri CA, Blenis J. Graded mitogen-activated protein kinase activity precedes switch-like c-Fos induction in mammalian cells. Mol Cell Biol. 2005;25(11):4676–4682. doi: 10.1128/MCB.25.11.4676-4682.2005. [PMC free article] [PubMed] [Cross Ref]
Harding A, Tian T, Westbury E, Frische E, Hancock JF. Subcellular localization determines MAP kinase signal output. Curr Biol. 2005;15(9):869–873. doi: 10.1016/j.cub.2005.04.020. [PubMed] [Cross Ref]
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):4776–4784. doi: 10.1091/mbc.E08-04-0407. [PMC free article] [PubMed] [Cross Ref]
Ferrell JE Jr, Machleder EM. The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science. 1998;280(5365):895–898. doi: 10.1126/science.280.5365.895. [PubMed] [Cross Ref]
McAdams HH, Arkin A. Stochastic mechanisms in gene expression. Proc Natl Acad Sci U S A. 1997;94(3):814–819. doi: 10.1073/pnas.94.3.814. [PubMed] [Cross Ref]
Bar-Even 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):636–643. doi: 10.1038/ng1807. [PubMed] [Cross Ref]
Cohen-Saidon 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):885–893. doi: 10.1016/j.molcel.2009.11.025. [PubMed] [Cross Ref]
Niepel M, Spencer SL, Sorger PK. Non-genetic cell-to-cell variability and the consEquationuences for pharmacology. Curr Opin Chem Biol. 2009;13(5–6):556–561. [PMC free article] [PubMed]
Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009;459(7245):428–432. doi: 10.1038/nature08012. [PMC free article] [PubMed] [Cross Ref]
Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S. Transcriptome-wide noise controls lineage choice in mammalian progenitor cells. Nature. 2008;453(7194):544–547. doi: 10.1038/nature06965. [PubMed] [Cross Ref]
Raser JM, O'Shea EK. Noise in gene expression: origins, consEquationuences, and control. Science. 2005;309(5743):2010–2013. doi: 10.1126/science.1105891. [PMC free article] [PubMed] [Cross Ref]
Wilkinson DJ. Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009;10(2):122–133. doi: 10.1038/nrg2509. [PubMed] [Cross Ref]
Pedraza JM, Paulsson J. Effects of molecular memory and bursting on fluctuations in gene expression. Science. 2008;319(5861):339–343. doi: 10.1126/science.1144331. [PubMed] [Cross Ref]
Raj A, Peskin CS, Tranchina D, Vargas DY, Tyagi S. Stochastic mRNA synthesis in mammalian cells. PLoS Biol. 2006;4(10):e309. doi: 10.1371/journal.pbio.0040309. [PMC free article] [PubMed] [Cross Ref]
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):1600–1603. doi: 10.1126/science.1119623. [PubMed] [Cross Ref]
Cai L, Friedman N, Xie XS. Stochastic protein expression in individual cells at the single molecule level. Nature. 2006;440(7082):358–362. doi: 10.1038/nature04599. [PubMed] [Cross Ref]
Colman-Lerner A, Gordon A, Serra E, Chin T, Resnekov O, Endy D, Pesce CG, Brent R. Regulated cell-to-cell variation in a cell-fate decision system. Nature. 2005;437(7059):699–706. doi: 10.1038/nature03998. [PubMed] [Cross Ref]
Perez OD, Nolan GP. Phospho-proteomic immune analysis by flow cytometry: from mechanism to translational medicine at the single-cell level. Immunol Rev. 2006;210:208–228. doi: 10.1111/j.0105-2896.2006.00364.x. [PubMed] [Cross Ref]
Taniguchi Y, Choi PJ, Li GW, Chen H, Babu M, Hearn J, Emili A, Xie XS. Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science. 2010;329(5991):533–538. doi: 10.1126/science.1188308. [PMC free article] [PubMed] [Cross Ref]
Shahrezaei V, Swain PS. Analytical distributions for stochastic gene expression. Proc Natl Acad Sci U S A. 2008;105(45):17256–17261. doi: 10.1073/pnas.0803850105. [PubMed] [Cross Ref]
Cohen AA, Kalisky T, Mayo A, Geva-Zatorsky 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. doi: 10.1371/journal.pone.0004901. [PMC free article] [PubMed] [Cross Ref]
Paulsson J, Berg OG, Ehrenberg M. Stochastic focusing: fluctuation-enhanced sensitivity of intracellular regulation. Proc Natl Acad Sci U S A. 2000;97(13):7148–7153. doi: 10.1073/pnas.110057697. [PubMed] [Cross Ref]
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]
Boykevisch S, Zhao C, Sondermann H, Philippidou P, Halegoua S, Kuriyan J, Bar-Sagi D. Regulation of ras signaling dynamics by Sos-mediated positive feedback. Curr Biol. 2006;16(21):2173–2179. doi: 10.1016/j.cub.2006.09.033. [PubMed] [Cross Ref]
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):337–351. doi: 10.1016/j.cell.2008.11.051. [PMC free article] [PubMed] [Cross Ref]
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):528–533. doi: 10.1073/pnas.0805981105. [PubMed] [Cross Ref]
Kholodenko BN, Birtwistle MR. Four-dimensional dynamics of MAPK information processing systems. Wiley Interdiscip Rev Syst Biol Med. 2009;1(1):28–44. doi: 10.1002/wsbm.16. [PMC free article] [PubMed] [Cross Ref]
Markevich NI, Tsyganov MA, Hoek JB, Kholodenko BN. Long-range signaling by phosphoprotein waves arising from bistability in protein kinase cascades. Mol Syst Biol. 2006;2:61. [PMC free article] [PubMed]
Sauro HM, Kholodenko BN. Quantitative analysis of signaling networks. Prog Biophys Mol Biol. 2004;86(1):5–43. doi: 10.1016/j.pbiomolbio.2004.03.002. [PubMed] [Cross Ref]
Birtwistle MR, Kolch W. Biology using engineering tools: the negative feedback amplifier. Cell Cycle. 2011;10(13):2069–2076. doi: 10.4161/cc.10.13.16245. [PMC free article] [PubMed] [Cross Ref]
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. doi: 10.1126/scisignal.2001212. [PubMed] [Cross Ref]
Fritsche-Guenther 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. [PMC free article] [PubMed]
Borisov N, Aksamitiene E, Kiyatkin A, Legewie S, Berkhout J, Maiwald T, Kaimachnikov NP, Timmer J, Hoek JB, Kholodenko BN. Systems-level interactions between insulin-EGF networks amplify mitogenic signaling. Mol Syst Biol. 2009;5:256. [PMC free article] [PubMed]
Kholodenko BN. Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur J Biochem. 2000;267(6):1583–1588. doi: 10.1046/j.1432-1327.2000.01197.x. [PubMed] [Cross Ref]
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):4334–4340. doi: 10.1074/jbc.270.9.4334. [PubMed] [Cross Ref]
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):1115–1127. doi: 10.1111/j.1600-0854.2009.00943.x. [PMC free article] [PubMed] [Cross Ref]
Goentoro L, Shoval O, Kirschner MW, Alon U. The incoherent feedforward loop can provide fold-change detection in gene regulation. Mol Cell. 2009;36(5):894–899. doi: 10.1016/j.molcel.2009.11.018. [PMC free article] [PubMed] [Cross Ref]
Nakakuki T, Birtwistle MR, Saeki Y, Yumoto N, Ide K, Nagashima T, Brusch L, Ogunnaike BA, Okada-Hatakeyama M, Kholodenko BN. Ligand-specific c-Fos expression emerges from the spatiotemporal control of ErbB network dynamics. Cell. 2010;141(5):884–896. doi: 10.1016/j.cell.2010.03.054. [PMC free article] [PubMed] [Cross Ref]
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 immediate-early gene expression via CArG elements. Mol Cell Biol. 1990;10(8):4243–4255. [PMC free article] [PubMed]
Kholodenko BN. Cell-signalling dynamics in time and space. Nat Rev Mol Cell Biol. 2006;7:165–176. doi: 10.1038/nrm1838. PMID: 16482094. [PMC free article] [PubMed] [Cross Ref]
Kiyatkin A, Aksamitiene E, Markevich NI, Borisov NM, Hoek JB, Kholodenko BN. Scaffolding protein Grb2-associated binder 1 sustains epidermal growth factor-induced mitogenic and survival signaling by multiple positive feedback loops. J Biol Chem. 2006;281(29):19925–19938. doi: 10.1074/jbc.M600482200. [PMC free article] [PubMed] [Cross Ref]
Aksamitiene E, Achanta S, Kolch W, Kholodenko BN, Hoek JB, Kiyatkin A. Prolactin-stimulated activation of ERK1/2 mitogen-activated protein kinases is controlled by PI3-kinase/Rac/PAK signaling pathway in breast cancer cells. Cell Signal. 2011;23(11):1794–1805. doi: 10.1016/j.cellsig.2011.06.014. [PMC free article] [PubMed] [Cross Ref]
Hartigan JA, Hartigan PM. The Dip Test of Unimodality. Ann Stat. 1985;13(1):70–84. doi: 10.1214/aos/1176346577. [Cross Ref]
Hartigan PM. Computation of the Dip Statistic to Test for Unimodality. Appl Stat J R Stat Soc Ser C. 1985;34(3):320–325.

Articles from BMC Systems Biology are provided here courtesy of
BioMed Central