## Abstract

We show theoretically and experimentally a mechanism behind the emergence of wide or bimodal protein distributions in biochemical networks with nonlinear input–output characteristics (the dose–response curve) and variability in protein abundance. Large cell-to-cell variation in the nonlinear dose–response characteristics can be beneficial to facilitate two distinct groups of response levels as opposed to a graded response. Under the circumstances that we quantify mathematically, the two distinct responses can coexist within a cellular population, leading to the emergence of a bimodal protein distribution. Using flow cytometry, we demonstrate the appearance of wide distributions in the hypoxia-inducible factor-mediated response network in HCT116 cells. With help of our theoretical framework, we perform a novel calculation of the magnitude of cell-to-cell heterogeneity in the dose–response obtained experimentally.

## 1. Introduction

Thermal fluctuations inherently affect all biochemical reactions, and the variability in molecular copy numbers owing to stochastic effects decreases, as the mean number of molecules in the system rises. Gene expression is inevitably stochastic, as usually there are not more than two gene copies in a cell. It is because of this biochemical noise that cells within an isogenic population, be it a bacterial colony or mammalian cells, at any given point in time exhibit a distribution of measurable properties rather than exhibit a precise value. Cells in genetically identical populations vary in, for example, size, the concentrations of proteins and mRNA, and the stage of progression along the cell cycle.

Quantification of biochemical noise requires single-cell techniques [1–5] capable of reconstructing the whole histogram (a distribution) of molecule counts over the entire cellular population. A bimodal (a double-peaked) distribution is particularly interesting from the physiological point of view, as it indicates existence of two subpopulations with (possibly) distinct phenotypic features. For instance, a bimodal distribution of antibiotic-resistance protein, ZeoR, was suggested to increase viability of the yeast population treated with large doses of Zeocin [6]. In another study, heterogeneous response to treatment with a chemotherapeutic drug, camptothecin, resulted in a bimodal distribution of 24 protein accumulation slopes. Additionally, for the RNA helicase DDX5 and the replication factor RFC1, peaks of their bimodal distributions correlated with different cell fates: fast protein accumulation was found in cells that survived the treatment, whereas slow accumulation existed in cells undergoing changes associated with cell death [7,8].

In biological systems, bimodal distributions of protein concentrations at the population level can be generated by several mechanisms (figure 1). A genetic switch is the simplest mechanism. For a large difference between the rate of gene on/off switching and the turnover of the gene product, for example, mRNA [9,17], two distinct peaks of product abundance arise (figure 1*a*). Such infrequent gene expression bursts leading to bimodality have been observed in *Escherichia coli* [18], yeast [19] and mammalian cells [2].

In the context of signalling networks, bistability is a well-known mechanism that can generate bimodal responses. It usually requires positive feedback that may be hidden or explicit (figure 1*b*) [10–12,20] or double negative feedback, as in synthetic toggle switch [21]. While stochastic switching requires non-deterministic analysis in order to show the existence of two states, bistable systems can be analysed deterministically with classical ODEs. Importantly, bistability in the deterministic setting does not necessarily lead to bimodal distributions. Dynamic switching between two stable states has to be promoted by fluctuating biochemical reactions or by external noise [22], or by cell-to-cell variability in the thresholds for switching between two steady states [14,23].

A much less appreciated mechanism capable of generating population-level bimodality relies on independently oscillating cells [13]. Cell-to-cell variability in protein abundances affects the period and/or phase of cells' oscillations which in turn generates a bimodal protein distribution (figure 1*c*). The existence of this distribution depends on the extent to which phases of oscillations vary between the cells (*phase mixing*) as well as the functional form of oscillations, for example, pure triangle wave will fail to generate bimodality. Notably, if cells oscillate with a slightly variable period, then the emergence of a bimodal distribution is inevitable with time.

The nature of the mechanism leading to bimodality that we discuss here differs from that of stochastic switching, bistability or superposition of oscillations. Here, two distinctive peaks in steady-state distributions of protein concentrations that arise on the population level stem from the interplay of two factors: (i) nonlinear input–output characteristic (the response) of the biochemical network at the single-cell level [24,25], and (ii) cell-to-cell variability in the abundance of network components which translates to different input–output characteristics of that network across the population (figure 1*d*).

In the following sections, we investigate analytically and numerically how bimodal steady-state protein distributions can emerge in nonlinear biochemical networks. We explore the effect of heterogeneity in the input signals (extrinsic noise) and/or variability in the input–output relationship owing to stochastically expressed network components (intrinsic noise) on the population-level distributions of protein concentrations. We apply our theoretical framework to infer variability in dose–responses across the cellular population from flow-cytometric measurements in the hypoxia response signalling system. We use a mathematical model of this system to demonstrate possible contributions to variability of the dose–response.

## 2. Results

### 2.1. Variability of response thresholds

The dose–response characteristics is a relationship between the magnitude of the input stimulus (e.g. growth factor concentration) and the output target, which can be the concentration of an active kinase (such as doubly phosphorylated MAPK) or a transcription factor. It is a function that depends both on the concentration of network components and kinetic parameters of biochemical reactions that involve these components. This input–output relationship, *R*(*x*), is commonly approximated by a Hill curve parametrized by four coefficients (*β*, *x*_{50}, *R*_{max} and *H*):
2.1where *x* is the input signal (e.g. concentration of a drug or a toxin), *β* is basal response level, *R*_{max} is the maximum output level and *H* is the Hill coefficient that determines the steepness of the response. The parameter *x*_{50} is the response threshold, also known as the half maximal effective concentration (EC_{50}). It is the concentration of the stimulus (the input *x*) for which the response assumes half of the maximum output (minus basal level).

We assume that cells within a population experience the same level of stimulation *x*. This condition holds, during our experiments, when cells are grown in a monolayer in a tissue culture and are subjected to equal levels of treatment (input). What sets cells apart is the threshold *x*_{50} of the response function *R*. Such variability may stem from differences in the molecule numbers of the same network components and fluctuations in their activation levels. In mathematical terms, we describe this variability by the probability density function (pdf) of thresholds, . The resulting distribution of output *y* for a given input value *x* in the presence of variable thresholds *x*_{50} of the response function *R* is given by a general expression (electronic supplementary material, equations (S1)–(S3)),
2.2where is the inverse of the response function calculated with respect to threshold *x*_{50}. Note the ‘minus’ sign, which reflects the fact that the response function is monotonically decreasing with *x*_{50}.

### 2.2. Existence of bimodality

Intuitive conditions for the existence of bimodality can be obtained for a model system with Hill-type response function and lognormally distributed threshold *x*_{50} across the population. It can be demonstrated numerically that such a parameter distribution is indeed a reasonable approximation for signalling cascades with concentrations of network components lognormally distributed across the cellular population (electronic supplementary material, figure S4). The assumption behind the distribution of protein concentrations is justified in the light of experimental evidence. Measurements have shown that protein distributions are lognormal around the mean with additional power-law tails that may arise from feedbacks in biochemical networks [26–28].

The first condition, necessary but not sufficient, relates the steepness *H* of the dose–response with the shape parameter of the lognormal distribution of *x*_{50} (electronic supplementary material, equation (S5)),
2.3where is related to the squared coefficient of variation of *x*_{50} distribution, , through: . Thus, steep (large *H*) steady-state response function *R* prompts a bimodal distribution even for very narrowly distributed thresholds (small ). And vice versa, bimodality may result from very heterogeneous but graded (small *H*) response function *R*.

Once *H* and satisfy equation (2.3), a bimodal output distribution may arise but only when the input stimulus, *x*, is within the steep region of the mean response curve. Parameters *H* and determine the width of that range. Bimodality will therefore ensue as long as the ratio of the input to the median of the threshold distribution, , satisfies
2.4where depends only on *H* and (electronic supplementary material, equation (S15)). The range of admissible ratios widens for a steep dose–response and/or large threshold variability (electronic supplementary material, figure S2).

The Hill function is linearly dependent only on *β* and *R*_{max}, hence variability in neither of these parameters alone can generate a bimodal output distribution (figure 2*b*). Variability in *H* can also introduce bimodality, however, for inputs around the midpoint of the dose–response the distribution reverts to unimodal (figure 2*c*). Because protein distributions for a range of input stimuli that we observe experimentally correspond to those shown in figure 2*a*, we assume that variability in response threshold *x*_{50} contributes the most to the population-level protein variability.

Variability in the input stimulus *x* is mathematically equivalent to variability in *x*_{50}; therefore, it is equally capable of generating bimodal distributions as has been argued previously in the context of transcriptional regulation [15,29–31]. In those scenarios, variable input signal corresponds to a noisy concentration of a regulatory protein, for example, a transcription factor. By fixing threshold *x*_{50} in equation (2.1) and by treating *x* as a variable subjected to fluctuations described by a lognormal distribution, we obtain conceptually similar results as previously. The first condition is the same as previously stated in equation (2.3) with the only difference that *σ* relates to the variability of the input stimulus rather than to the threshold. The interpretation of the second condition changes accordingly. Function bounds the ratio of the input distribution's median *m* and now fixed threshold level *x*_{50}.

Equations (2.3) and (2.4) that warrant bimodal system response are derived for a specific case of a lognormal input or response threshold distribution and a Hill-type dose–response. However, the interpretation of the two conditions also holds for other types of distributions and responses, for example, gamma distribution, sigmoidal and logistic response characteristics. Equation (2.3) simply states that in order to facilitate a bimodal distribution both the steepness of the response and variability of thresholds must be coupled; reduction of the former requires an increase in the latter and vice versa. Equation (2.4) prescribes how far the midpoint of the threshold distribution (its ‘centre of mass’) can be from the midpoint of the response (the steepest point) in order to maintain bimodality. These two conditions hold for other cases as well, although their exact mathematical form is not as simple as for lognormal distribution and Hill-type response. Notably, bimodality may arise only for nonlinear response characteristics. A linear dose–response would result in a trivial result where the output distribution is simply a rescaled distribution of a given parameter in that response function.

Additionally, in order to generate bimodality, the dose–response cannot equal the cumulative distribution (CDF) of the variable parameter or, equivalently, the distribution cannot be a derivative of the response function [32]. If this is the case, then the output distribution becomes uniform on the interval between the lowest and highest response. Consider a cell population with the response function and the threshold distribution across the ensemble equal to the derivative of *R* with respect to *x*_{50} (with a minus sign). Because depends also on the input stimulus *x*, the population-level response to any input value is a uniform distribution between 0 and 1. It implies that any signal evokes the exact same population-level response. However, from the biological perspective, the scenario in which the dose–response always remains equal to the CDF of the threshold distribution is intriguing, yet unlikely. It requires *x*_{50} distribution to change with *x*, which in turn would affect the dose–response curve. In a more plausible scenario, the system adapts its dose–response to the distribution of inputs that it experiences. It has been shown that a number of neural systems use the strategy of matching the dose–response to CDF of inputs in the environment [32,33] to fully use neuron's limited response range and maximize the information transfer: the most frequently perceived stimuli evoke outputs from the steepest region of the dose–response. Arguably, cellular signalling networks as discussed throughout this paper undergo this type of adaptation.

According to equation (2.4), the range of input/median ratios for which bimodality arises depends on the width of the threshold distribution, (electronic supplementary material, figure S2). The larger it is, the stronger the separation of the heterogeneous population into two groups with distinct levels of protein concentrations. Figure 3*a–c* depicts this intriguing property that runs counter to the conventional assumption that cellular variability destroys robust signalling. Here, we consider a system with a mildly ultrasensitive, *H* = 3, dose–response. Compare this with *H* ≈ 5 reported for the MAPK cascade [24], *H* = 5 … 9 observed for Rap–GTP responding to cannabinoid-1 receptor signal [35], or *H* = 1 … 4 measured for a synthetic system with multiple autoinhibitory modules [36]. For , protein distributions become significantly wider for input stimuli in the steepest part of the dose–response. For , the responses tend to concentrate around basal and saturation values, and two peaks emerge for intermediate stimuli. Such bimodality may facilitate further decision-making, which is not entirely random but is performed based on two well-defined options instead.

Uncertainty of a distribution is often quantified by Shannon entropy [37]. Its low value indicates that a small amount of information is required to describe the varying quantity, for example, protein concentration across the population. Therefore, a peaked unimodal distribution requires a twice shorter description than a bimodal distribution with two narrow peaks. Any wide distribution in between requires more information than either of the two, which is reflected by the high value of entropy. Figure 3*d* quantifies distributions induced by inputs equal to *x*_{50} of the dose–response (red dots in figure 3*a–c*). The entropy peaks at and decreases owing to uncertainty introduced by threshold variability. However, for biologically realistic values of between 0.4 and 4, the distributions change their shape from wide to bimodal, whereas their entropy changes only slightly. This suggests that the quantification of uncertainty using entropy may be misleading as even protein distributions with large entropy may reveal a physiologically relevant widening of a population into two groups. Even though this division is not complete and a number of responses still appear between the peaks, such distributions have been demonstrated to be physiologically advantageous to the population [6].

### 2.3. Variability in the hypoxia response network

Here, we study experimentally the cellular hypoxia response network, analyse the emergence of wide protein distributions and estimate the response variability from experimental data. When cellular oxygen demand exceeds supply, the cells enter a phase of hypoxia. As a consequence, stabilization of the hypoxia-inducible factor (HIF) ensues. Stabilized HIF mediates transcription of genes to adapt to the hypoxic insult [38]. Central to the response is the action of prolyl-hydroxylases (PHD), enzymes that hydroxylate HIF at residues Pro-402 and Pro-564 [39] and target it for ubiquitination–degradation via the von Hippel–Landau protein (VHL) [40]. Figure 4 depicts a simplified scheme of the network.

For our experimental set-up, we used a stable HCT116 cell line expressing a fragment of the HIF protein containing residues 403–603, termed the oxygen-dependent degradation (ODD) domain [41] tagged to GFP (cells courtesy of Prof. E. Gottlieb [42]). The ODD–GFP is our readout of the hypoxic response. We activate the system using the hydroxylase inhibitor dimethyloxalylglycine (DMOG), which mimics the condition of low oxygen levels in the HIF system [43]. Cells in tissue culture were grown up to 70% confluency at the end of the treatment, which minimized the effect of cell contact and maintained cells in a monolayer such that all of them were exposed to equal levels of DMOG. Hence, any variability in the response can be attributed to intrinsic variations of network components in individual cells, which facilitates our aim of measuring dose–response variability while assuming a fixed input. The condition, however, may not hold in general, especially when cells are embedded in tissue and/or subjected to different microenvironments.

### 2.4. Hypoxia-inducible factor responses to dimethyloxalylglycine averaged over the cell population

Using flow cytometry, we first identify a sigmoidal dose–response. For each DMOG condition, we calculate the median of the single-cell ODD–GFP fluorescence across a population of a minimum of 10 000 cells. In doing so, we derive the dose–response curve, i.e. the ODD–GFP response to DMOG. We normalize the median to control, unstimulated case. The median response averaged over a minimum of three biological repeats at 4, 8 and 16 h following DMOG addition is shown in figure 5*a* where we also fit the Hill model from equation (2.1) (electronic supplementary material, table S1). The increase with time of the steepness, *H*, and the maximum response, *R*_{max}, is consistent with earlier measurements using Western blots [44], which support the qualitative correspondence between population-averaged flow-cytometric and bulk measurements.

### 2.5. Hypoxia-inducible factor response distributions across the cell population

Distributions of ODD–GFP fluorescence intensity (figure 5*b*) increase their width for DMOG = 1 and 2 mM. At these DMOG concentrations, the response curve is the steepest; it is the most sensitive region of the cellular response. Widening of the ODD–GFP distribution for DMOG concentrations that induce the steepest region of the dose–response hints at the existence of cell-to-cell variability in response thresholds (cf. figure 2). If there were no threshold variability, then cells treated with DMOG would exhibit a peak of ODD–GFP distribution shifting towards higher values without increase in its width as shown in figure 3*a*. This is not what we observe experimentally. In the presence of threshold variability and steep nonlinear response, the distribution may widen or even become ‘split’ between the low and high response—the fact quantified by equations (2.3) and (2.4), and illustrated by figure 3*b,c*. Consistently with our predictions, the observed widening of the ODD–GFP distribution becomes more pronounced the steeper the dose–response, and the larger, the maximum response level. This is the case at 8 and 16 h post-treatment.

As a next step, we estimate distributions of parameters *β*, *R*_{max} and *x*_{50} in the Hill model, equation (2.1), from a single experiment. Parameter *H* is assumed to be constant across the population, and we calculate it by fitting the Hill curve to population-averaged dose–responses (electronic supplementary material, figure S5*a* and table S2). The distribution of basal level *β* is obtained from the ODD–GFP distribution measured for the unstimulated case, DMOG = 0 mM (electronic supplementary material, figure S5*b* and table S3). Based on our fits of the Hill function to the dose–response, we recognize that even at DMOG = 4 mM the curve does not saturate. However, the ODD–GFP distribution at 8 and 16 h post-DMOG treatment becomes left-skewed, which indicates that the system is not far from saturation. Because higher DMOG levels would be toxic to cells, we decide to use the response at DMOG = 4 mM as a proxy of *R*_{max} variability. At that DMOG level, the response is also affected by an independent contribution from variability in the basal level *β*. We therefore subtract the mean and variance of *β* from the distribution owing to DMOG = 4 mM and use this new distribution as an estimate of *R*_{max} variability (electronic supplementary material, figure S5*c* and table S4).

Finally, we estimate variability in threshold *x*_{50}. At a given time point, for every DMOG concentration, we record the value of ODD–GFP distribution at a half-distance between the mode (the peak) of the basal (DMOG = 0 mM and the maximum (DMOG = 4 mM) distribution. Such a curve, which is a function of DMOG level, is largely independent of *β*, *R*_{max} and *H*, and is determined only by parameters of *x*_{50} distribution. The independence is exact when *x*_{50} is distributed lognormally (electronic supplementary material, equation (S17)). After normalization to 1, the function has the interpretation of the probability density, which can be used to determine parameters of the *x*_{50} distribution (electronic supplementary material, figures S3 and S5*d* and table S5).

For spread of parameters *β* and *R*_{max} across the cell population, we fit gamma, lognormal and Weibull distributions. Based on Akaike information criterion [45], we conclude that Weibull distribution is the best approximation of experimentally measured variability in both parameters (electronic supplementary material, figure S5*b*,*c* and tables S3 and S4). By applying these distributions to equation (2.1), we numerically calculate ODD–GFP distributions shown in figure 5*c*. Figure 5*d* quantifies the agreement between experimental and predicted ODD–GFP distributions with a dimensionless robust coefficient of variation, rCV, expressed by equation (4.2) in Material and methods.

According to our model, wide distributions that we observe experimentally can arise only in the presence of variability in the response threshold, *x*_{50} (figure 2). We find that the dimensionless squared coefficient of variation (CV^{2} = variance over squared mean) of fitted *x*_{50} distribution decreases from 2.72 at 4 h, through 1.38 at 8 h, down to 0.26 at 16 h post-treatment. The decrease in over time is compensated by increasing *H*, which promotes widening of the ODD–GFP distribution at DMOG concentrations that induce the steep part of the response. Substituting *H* and into equation (2.3) reveals that the system remains in the same regime where only mild bimodality may arise (electronic supplementary material, table S6 and figure S6). The distribution is further affected by variability in *β* and *R*_{max}, which results in widening instead of the emergence of two clearly separated peaks. The widening is reflected by the increase in rCV (figure 5*d*). In the case of small *x*_{50} variability, or lack thereof, an increasing stimulus would induce a proportionally shifting peak of ODD–GFP distribution without widening as illustrated in figure 3*a* and quantified by a theoretical dotted line in figure 5*d*.

### 2.6. Response variability in a mathematical model of the hypoxia-inducible factor system

We used a published mathematical model of the HIF system [44] to search for possible molecular sources causing large dose–response variability. First, we calculate numerically the response of the system at 16 h to various oxygen tension levels, which affects the HIF level in the opposite manner to DMOG treatment. The untreated case corresponds to 21% of oxygen, i.e. normoxia, whereas high DMOG treatment corresponds to the hypoxic condition of a low oxygen level. The Hill function is then fitted to the calculated dose–response. In order to account for cell-to-cell heterogeneity, we repeat this procedure 1000 times with total levels of all protein species drawn from a lognormal distribution with CV^{2} = 0.5 and means fixed to the original ODE model (supplementary material of [44]). The chosen CV^{2} is at the upper limit of measured protein variability in mammalian cells which typically ranges between 0.1 and 0.5 [7,16,28,46,47].

The mean and variability of the dose–response obtained from a population of deterministic models are shown in figure 6*a*. By fitting the Hill function to each of the responses from the entire population, we obtain distributions of parameters *x*_{50}, *β*, *R*_{max} and *H* at 16 h post-treatment (figure 6*c–f*). The distribution of *H* has the smallest CV^{2} of all parameters, which justifies our choice of neglecting variability in *H* when estimating variability of Hill parameters from experimental data in §2.5. The distribution of *x*_{50} is well approximated by a lognormal distribution, and is comparable to 0.26 that we estimated from experimental data in figure 5*c*. Additionally, we perform local sensitivity analysis by randomly sampling total concentrations of individual protein species from a lognormal distribution with the mean as in the ODE model and CV^{2} = 0.5 while keeping the remaining concentrations fixed. Interestingly, the variability of the response threshold in the HIF model is largely insensitive to variability of the majority of network components (figure 6*b*). The variability in only two protein concentrations, nuclear FIH and HIF-1β, translates to large CV^{2} of response thresholds.

## 3. Discussion

Signalling networks can transform analogue, continuous input signals into distinct low and high response states [48]. The sigmoidal shape of that response is typically a result of layers of phospho- and dephosphorylation cycles (e.g. MAPK cascade [24]) or simply results from nonlinearity of enzymatic reactions as in the case of the HIF network. The steep sigmoidal response minimizes the region where stimuli could result in ambiguous outputs, making such a response characteristic a prevailing biological mechanism on which cellular decision-making can rely. Because of the large number of molecules that eukaryotic signalling networks typically consist of, at the level of a single cell, these systems seem to be protected from undesirable stochastic effects, which can occur at low molecule counts. Still, cells within a population differ in their molecular make-up because of the biochemical noise that affects total protein levels. This randomness affects the shape of the dose–response and causes a similarly random response of a cellular population.

Our theoretical results demonstrate that networks with sigmoidal response characteristics possess yet another fascinating feature, this time at the level of cellular ensemble. In the presence of variability in nonlinear dose–response across the population, two subpopulations can emerge with distinct response levels. The condition being, dose–response needs to be of sufficient steepness. Our analysis provides intuitive mathematical insights into conditions under which such a response in form of a bimodal protein distribution may arise.

Cell-to-cell variability of the dose–response threshold, *x*_{50}, more commonly known as EC_{50}, is essential for generating wide or bimodal distributions in response to a common input signal. In the absence of such response variability, protein distribution would not widen but shift proportionally to the stimulus (figure 3*a*), which is not what we observed in our experiments. Flow-cytometric analysis of the HIF system demonstrated that the response, measured as ODD–GFP distribution, indeed widens for inputs (DMOG treatment) in the steepest region of the Hill-type dose–response. Our theoretical framework provides a means of estimating EC_{50} variability from flow cytometry data and an explanation of protein distribution widening.

Without EC_{50} variability, the ODD–GFP distribution would not widen as is illustrated by the dotted line in figure 5*d*. Even though the squared coefficient of variation of the response threshold decreased 10 times between 4 and 16 h after DMOG treatment, the steepness *H* of the response increased more than twice over the same period. Together, these two parameters facilitate a necessary but not sufficient condition for bimodality expressed by equation (2.3). After substituting EC_{50} variability and *H*, we find that the system remains approximately in the same regime where only very mild bimodality may arise (electronic supplementary material, table S6 and figure S6). Using our framework, we hypothesize that the response steepness at 16 h post-treatment is sufficient to bring about two distinct peaks in the response, provided that CV^{2} of EC_{50} remained at the level of 2.72 (as estimated at 4 h) throughout the treatment (electronic supplementary material, figure S6).

We have shown that owing to certain features of the nonlinear input–output characteristics commonly observed in cellular signalling networks, random differences between individual cells enable the separation of these cells into subpopulations with distinct responses to the common signal. The significance of our results extends beyond cascaded signalling networks. Any network with nonlinear dose–response characteristics may exhibit bimodal behaviour as shown in our study. In particular, a simple enzymatic reaction or a gene expression network also belong to that class. For all of such nonlinear networks, variability of EC_{50} caused by cell-to-cell differences in protein concentrations is an important factor to consider when designing a perturbation, for example a drug treatment. According to our theory, one of the strategies to avoid the emergence of two subpopulations, low and high responders, is to smooth the dose–response, i.e. decrease the coefficient *H*. This in turn can be achieved by establishing a negative feedback in the network [14]. Alternatively, the condition expressed by equation (2.3) demonstrates that the output protein distribution can be rid of bimodality by decreasing the variability of dose–response threshold resulting from cell-to-cell variability of network components. With current advancements in synthetic biology and ongoing progress in mathematical modelling of biochemical networks, such interventions in *in vivo* systems will likely become more routine. An intriguing open question is to what extent bimodal population-wide protein distributions brought about by the mechanism discussed here are used to perform physiologically relevant decision-making.

## 4. Material and methods

### 4.1. Cell culture

Human colon carcinoma (HCT116) cells stably expressing GFP–ODD were a generous gift from Prof. Gottlieb and were previously used in another study [42]. The cells were subsequently subcloned as single-cell clones by growing cells in DMEM 10% FCS supplemented with 800 µg ml^{−1} G418. The clones were treated with 2 mM DMOG for 16 h. Two clones where eGFP–ODD was induced were selected. Cells were then maintained in a humidified 5% CO_{2} incubator at 37°C and cultured in Dulbecco's modified Eagle's medium supplemented with 10% FCS, 1% l-glutamine and 50 µg ml^{−1} G418.

### 4.2. Flow cytometry

HCT116 cells were treated with cell permeable pan-hydroxylase inhibitor dimethyloxalylglycine (DMOG; Cayman Chemicals, MI, USA) dissolved in DMSO (Sigma, Wicklow, Ireland) and diluted with growth medium to appropriate concentration. The cells were then lifted by trypsinization (0.05% trypsin–EDTA, Gibco) and resuspended in 0.5 ml growth medium with DMOG concentration corresponding to initial stimulation. Experiments were replicated four times for each time point. The samples were analysed with Accuri C6. Post-gating by forward and side scatter was performed to remove events corresponding to dead cells, debris and cell clusters (i.e. doublets). For each sample, 10 000 events (cells) were collected in final gating.

Dose–response from flow cytometry experiments (figure 5) was calculated using normalized median fluorescence intensity (nMFI), 4.1

In order to quantify the width of protein distributions obtained from flow cytometry and calculated numerically, we use robust CV (from FlowJo flow cytometry analysis software) defined as follows 4.2

The rCV is not as sensitive to outliers and gives a better than CV description of the distribution spread on logarithmic scales.

## Funding statement

This study is supported by Science Foundation Ireland (under grant no. 06/CE/B1129) and European Union Grant PRIMES (no. FP7-HEALTH-2011-278568).

## Acknowledgements

We thank Lisa Heiserich for kindly providing us with HCT116 cells and the anonymous reviewers for their insightful comments.

- Received April 11, 2014.
- Accepted June 5, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.