Many examples (particularly models of gene regulation) have been studied wherein a deterministic treatment of the kinetic scheme describing the relevant processes has two stable steady states in a certain parameter regime (1–3, 13, 14). In such systems, stochastic effects can lead to bimodality (e.g., populated “on” and “off” states) in the parameter range where bistability is predicted by the deterministic equations as well as outside this range where there is a single stable steady state (1–3, 13, 14). The latter phenomenon is a consequence of stochastic fluctuations enabling the system to sample parameters (e.g., rate constants) that effectively fall within the range where two deterministically stable fixed points are present. In these examples, the existence of bistability in the deterministic analysis in some parameter range underlies the observation of bimodal behavior in the stochastic treatment.

The model we study exhibits a different feature. The deterministic dynamical equations yield a single steady state in all parameter ranges; i.e., there is no bistability. Yet, stochastic fluctuations result in a bimodal long-time response with neither mode corresponding to the steady state obtained deterministically. Upon increasing the copy numbers of molecules, the stochastic description ultimately converges to the deterministic behavior. Thus, we find a purely stochastically driven instability when none exists in the deterministic treatment in any parameter range. When fluctuations are important, we find that average quantities scale with parameters “anomalously” compared with the corresponding mean-field behavior. Our analyses suggest that the necessary and sufficient conditions for this phenomenon to occur are quite common.

Our simple (“toy”) model is inspired by ideas proposed recently to describe T cell responses to diverse stimuli (7–9, 15–17). T cell receptor (TCR) molecules expressed on the surface of T cells can bind complexes of peptides (p) bound to MHC proteins on the surface of antigen-presenting cells (APCs). TCR can potentially bind strongly to pMHC molecules where the peptide is derived from a pathogen's proteins (agonists). In contrast, thymic selection ensures that TCR bind weakly to “self” or endogenous pMHC molecules that are also expressed on APCs (18). The binding of TCRs to pMHC molecules can initiate signaling cascades that result in T cell activation and an immune response. T cells are as good a sensory apparatus as any in biology, and can detect as few as three agonists in a sea of tens of thousands of endogenous pMHC molecules, and it has been suggested that this extraordinary sensitivity is mediated by cooperative interactions between self pMHC and agonists (7–10, 19).

Another interesting response of T cells to pMHC molecules is called antagonism (15, 20). Antagonists are pMHC molecules obtained by mutating agonist peptide residues. When present on APC surfaces in sufficient numbers, they can shut down intracellular signaling stimulated in response to agonists. Recent experimental results (15) have suggested that this phenomenon may be mediated by dueling positive and negative feedback loops (Fig. 1). One of the earliest steps in downstream signaling initiated by the binding of the TCR to pMHC molecules is the phosphorylation of cytoplasmic domains of the TCR complex by a kinase called Lck. It has been proposed that Lck also activates its own inhibitor, a phosphatase called Shp (negative feedback). This inhibitory interaction is prevented by a product (ERK) of signaling downstream of phosphorylation of the TCR complex that protects Lck by phosphorylating one of its sites (positive feedback). It has been proposed, and supported by detailed calculations (16, 17), that the positive feedback is dominant when T cells are stimulated by agonists (and synergistic endogenous ligands), and negative feedback shuts down signaling when sufficient numbers of antagonists are present.

Fig. 1. A schematic representation of dueling positive and negative feedback loops stimulated upon receptor binding to stimulatory or inhibitory ligands. The negative regulator can shut off signaling by inactivating the receptor-associated signaling complex (negative (more ...)

Although the specific molecular identity of positive and negative regulators involved in T cell signaling is still debated (21), the idea that dueling positive and negative feedback loops play a role in determining whether signaling is shut off (antagonism) or sustained/amplified (agonism) is of general significance to cellular decisions that lead to distinct outcomes. Furthermore, such processes are often mediated by small numbers of molecules. Therefore, we set out to study the effects of stochastic fluctuations on the following simple and general model with dueling positive and negative feedback regulation

The mean-field deterministic equations corresponding to the model described by Eqs. 1–6 can be written down following mass action kinetics [supporting information (SI) Text, SI Table 1, and SI Figs. 6–9], and yield the following solution for the steady state:

Although we have studied different parameter ranges for a stochastic description of this model (SI Text), let us first consider situations that are inspired by T cell signaling. Reactions 1, 2, and 4 represent multistep processes (22). Reactions 3 and 5, the dueling feedback loops, are thought to represent one-step phosphorylation or deactivation steps (15). So, we study situations where k₃ and k₅ are much larger than k₁, k₂, and k₄; i.e., both positive and negative feedback loops are strong. Recent studies (9, 17) with detailed models of membrane-proximal signaling in T cells suggests that k₄ could be larger than k₁, but we have taken them to be equal (k₄ > k₁ is considered in the SI Text). Changing the relative values of k₁ and k₂ would simply modify the specific value of the ratio of initial numbers of A₁ and A₂ molecules that would result in a transition from “agonism” to “antagonism.”

Fig. 2 shows results of spatially homogeneous stochastic simulations with discrete number of molecules [using the Gillespie algorithm (23)] of the model represented by Eqs. 1–6. When there are only a few molecules of A₁ and A₂, essentially all of the stochastic trajectories commit to one of two final states: all of the A₁ molecules are converted to the protected species, A₁^PROT, or are annihilated and signaling stops. This bimodality is in striking contrast to the mean-field solution that does not exhibit bistability for any parameter values. The qualitative phenomenon of finding a bimodal stochastic solution when the deterministic solution is unique for all parameter values is preserved as long as the positive and negative feedback loops are sufficiently strong (SI Text).

Fig. 2. Bimodal stochastic solutions distinct from the unique deterministic solution. Histogram showing the bimodal distribution of the protected agonists at steady state (red) for a situation when there are 10 agonist (A1) and 10 antagonists (A2). The corresponding (more ...)

The mechanism underlying this result is as follows. The species A₁ is converted to either A₁^PROT or A₁^INACT. The effective rates of production of these species can be obtained from the deterministic equations. Both rates equal zero initially and at long times, and exhibit a maximum (SI Fig. 9). The initial rise and amplitudes of the maxima depend on the values of the initial number of A₁ and A₂ molecules, and are very different if one of these quantities is much larger than the other. In these circumstances, either agonism or antagonism dominates in the deterministic and stochastic solutions. The more interesting cases are ones where the generation of positive and negative regulations is roughly balanced (Fig. 2) because it could result in a transition from agonism to antagonism. Now, the rates at short times and amplitudes of the maxima for the production of A₁^PROT and A₁^INACT are comparable in the mean-field sense, and the deterministic equations yield a single steady state solution with an intermediate value of A₁^PROT. However, stochastically, one of two reactions 1 and 2 occurs first. There is a stochastic delay, τ, before the other reaction occurs, and for this duration, the reaction propensities are effectively as in cases where A₁ A₂, or vice versa. For small numbers of A₁ and A₂ molecules, τ can be long. If τ is longer than the intrinsic time scale associated with the feedback reaction corresponding to the reaction that occurred first (e.g., reaction 3 if 1 occurred first), then the small number of A₁ molecules will all be converted to either A₁^PROT or be annihilated, depending upon whether reaction 1 or 2 occurred first. So, the stochastic trajectories partition into two classes (those that end with all A₁ molecules annihilated or protected), and the stochastic solution is bimodal.

The time delay (τ) becomes smaller as the number of molecules of A₁ and A₂ increases. This observations suggests that, for a sufficiently large number of particles, it will not be longer than the intrinsic time scale associated with the feedback loops and the stochastic solution will not be bimodal. Rather, it will be distributed around the mean-field solution. Fig. 3 shows results of simulations that demonstrate this unequivocally. Thus, for the same parameter values, as the number of molecules decreases past a threshold, the stochastic solution exhibits an instability from one solution to bimodality. This transition from unimodal to bimodal solutions is driven by stochastic effects, and occurs in the absence of any underlying deterministic bistability.

Fig. 3. A purely stochastically driven transition. Histograms showing the distribution of the protected agonists at steady state (red) and corresponding steady state solution of deterministic ODEs (blue) for different amounts of agonist and antagonist. All other (more ...)

The qualitative differences between the stochastic and deterministic descriptions due to the dominance of fluctuation effects suggests that the manner in which the response scales with different control parameters may be different. For example, we expect the steady state amount of A₁^PROT to scale with k₁A₁/k₂A₂ for the stochastic simulations. This is because the probability of conversion to A₁^PROT is essentially equal to the probability that reaction 1 occurs first, which is given by k₁A₁/(k₁A₁ + k₂A₂). Conversely, probability of annihilating all A₁ molecules is equal to k₂A₂/(k₁A₁ + k₂A₂) (equal to probability that reaction 2 occurs first). Both expressions depend only on the combination k₁A₁/k₂A₂. This implies, for example, that the amount of A₁^PROT scales linearly with k₁ (a measure of how effective the agonist is in stimulating signaling). The deterministic solution, on the other hand, is not expected to obey this linear scaling. Indeed, numerical solutions support these expectations (SI Fig. 8).

The complexity of the model described by Eqs. 1–6, however, makes it difficult to explore these differences in scaling behavior precisely. The complexity also prevents us from analyzing the necessary and sufficient conditions for purely stochastic instabilities (results in Figs. 2 and 3) in cell signaling dynamics. Therefore, we formulated a simpler model that enabled exploration of these issues.

This minimal model, which can be solved exactly, includes the following features: irreversibility, branching, and feedback. The model is described in terms of the three coupled reactions shown below:

The following Master Equation describes the stochastic time evolution of the reactions shown in (8)

Using the method of generating functions (24), Eq. 12 can be solved exactly (SI Text) to obtain:

Eq. 15 results in a steady state probability distribution that is bimodal for small numbers of molecules (SI Fig. 11a) when the deterministic solution does not exhibit bistability in any parameter range. In the more complex model that we studied (Eqs. 1–6), mean-field behavior was obtained as the numbers of A₁ and A₂ molecules increased past a threshold value even though their relative numbers were kept constant. The corresponding limit for the minimal model is k₃ → ∞, N → ∞, with the ratio N/k₃ (or the dimensionless, Nk₁/k₃) remaining constant. This is because a large value of N corresponds to a large amount of the source of a positive regulator (A₁ in Eqs. 1–6) and a large value of k₃ corresponds to greater annihilation or a big source (A₂ in Eqs. 1–6) of negative regulation. SI Fig. 11b shows that, like the more complex model, there is a purely stochastic transition as the stochastic solution is unimodal and distributed around the deterministic solution above a threshold value of N and k₃. So, these results establish that the sufficient conditions for the phenomena we report are: irreversibility, branching, and feedback loops. But, are these also necessary conditions?

The possibility of two different outcomes is obviously necessary, and branching is ubiquitous in cell signaling processes that lead to functional decisions. We have also found that removing irreversibility abolishes the phenomenon (data not shown). Ultimately, all reactions are, in principle, reversible. However, in the time scales of interest to signal propagation in cells, many steps are effectively irreversible.

Feedback regulation is also necessary as the bimodal stochastic solution does not exist if k₂ in the minimal model tends to zero (SI Fig. 15). Insight into the kind of feedback regulation that is necessary can be obtained by contrasting our studies of dueling feedback loops in cell signaling to a model for binary drift in population genetics (25, 26). Consider a population of heterozygote individuals with two forms, B₁ and B₂, for a particular allele. In the absence of mutations, the number of each type of allele can change from generation to generation, even in a population of fixed size, due to mating. The effects of binary selection on the numbers of B₁ and B₂ forms can be roughly represented as follows (25, 26):

The model described by Eq. 16 also contains branching and irreversibility. There is also an effective feedback, but unlike Eqs. 1–6 or 8, there is no separate intrinsic time scale associated with the feedback loops. A special case of this model (with no selection), k = k′, shares some features with the systems we are considering. The deterministic changes in this limit are trivially zero, and any initial condition (along the fixed line of B₁ + B₂ = population size) remains fixed. These deterministic steady states are unique, but the stochastic solutions yield a bimodal distribution. This is because the stochastic trajectories are divided into two classes: ones that terminate when the number of B₁ particles vanishes and those that terminate when the number of B₂ reaches zero.

There is an important difference, however, between the model for binary drift with k = k′ and the class of cell signaling models we have been considering. The stochastic solution of the model represented by Eq. 16 does not converge to the deterministic solution when the number of particles becomes large. The stochastic solution at t → ∞ is always bimodal. The stochastic trajectories cease to evolve when either B₁ or B₂ become zero because only then is the effective rate of conversion between these species equal to zero. The deterministic rates of formation of B₁ and B₂ equal the same constant for all times. Increasing the numbers of molecules does not eliminate this difference between the deterministic and stochastic cases. As the number of particles increases, the stochastically determined time (τ′) required for B₁ or B₂ to equal zero increases, but ultimately it always happens. There is no separate intrinsic time scale that can compete with increasing values of τ′ as the number of particles increases and prevent this from happening (i.e., a bimodal solution). Recall that, for the signaling models that we focused on, the relative values of the stochastic delay, τ, and the separate time scale associated with feedback loops determined the stochastically driven transition when the number of molecules was lowered (Fig. 3 and SI Fig. 11b). The absence of such an interplay prevents a purely stochastic instability in the binary drift model as the number of particles decreases. Correspondingly, if the rate coefficients in the model represented by Eq. 16 were time dependent with an intrinsic time scale, the phenomenon of a purely stochastic instability would be recovered.

The analyses presented above suggest that the necessary and sufficient conditions for a purely stochastic bimodality in the absence of any deterministic bistabilities are: (i) irreversibility, (ii) branching, and (iii) feedback regulation with an associated distinct and fast time scale.

The analytical solution for the probability distribution (Eq. 15) obtained for the minimal model of cell signaling that satisfies these conditions enables us to calculate average properties, such as the average number of molecules of the product, x. This allows us to examine whether x scales with parameter values in the same way as N_x determined from the mean-field equations (see above). The average value, x, is

To answer this question, as shown above, we need to consider the value of x in the limit of large values of N, M, and k₃. Consider first the limit of large values of N and k₃ for a fixed value of Nk₁/k₃. Simple algebra yields the value of x in this limit to be

The general solution (Eq. 18) for x does not allow us to explicate the non-mean-field scaling when fluctuations are important. This can be obtained analytically only in special limits. For example, consider the limit of infinitely strong feedback (k₂ → ∞). In this limit, x takes the following form (SI Text)

Fig. 4. Results from the minimal model. Non-mean-field scaling in the limit of large positive feedback (k2 → ∞): The average values of X species, x, at steady state obtained from Gillespie simulations scale with (1/log(1 + Mk (more ...)

Dueling positive and negative feedback loops are ubiquitous in biology. In many instances, these processes involve small numbers of the pertinent molecules, and hence stochastic fluctuations can be important. We report a striking result for such systems. The models we have studied correspond to unique deterministic steady states for all parameter values, and do not exhibit bistability. Yet, when there are a small number of molecules, stochastic effects result in a bimodal solution with neither solution corresponding to the mean-field result. Our analyses suggest that the necessary and sufficient conditions for this phenomenon are irreversibility, branching, and the existence of an intrinsic and relatively fast time scale associated with feedback regulation. Our studies show that for specific examples of such systems, near the transition from one phenotype to another (e.g., agonism to antagonism), mean-field scaling does not apply to the stochastic solutions. Whether the specific differences in scaling between mean-field and stochastic solutions that we report are universal for the class of models that exhibit the phenomenon revealed by our studies remains an open question.

There is a key difference between models of gene regulation and cell signaling where bimodality has been observed in stochastic limits under conditions where the deterministic equations yield monostable solutions and our results. In the former examples (double negative feedback, dimer mediated gene regulation, etc., e.g. refs. 1 and 27–32), bistable deterministic solutions exist in some other parameter regime. Stochastic bimodality displayed by binary drift models in population genetics are also different from the phenomena we report in that the stochastic solutions are always bimodal, regardless of the number of particles; i.e., there is no stochastically driven transition from a single solution to bistable solutions.

The necessary and sufficient conditions for the phenomenon that we report (branching, irreversibility, and feedback loops with distinct time scales) are quite common in cell biology. Our results suggest that these features, when combined with stochastic fluctuations, can enable cells to make binary decisions, whereas this would not be possible in a deterministic world. For instance, if gene transcription and effector function required greater than a threshold value of a downstream signaling product, in a mean-field world, cells would be unable to make decisions with a distinct functional outcome (Fig. 5). Under the same conditions, stochastic effects would result in cells being either “on” or “off” (Fig. 5), as observed in experimental studies in diverse contexts.

Fig. 5. Stochastic fluctuations can enable cellular decisions. Schematic representation showing that irreversibility, branching, and dueling feedback loops associated with intrinsic time scales, when combined with stochastic effects, can result in distinct functional (more ...)

For example, a recent study of HIV latency by Weinberger et al. (4) showed a “temporary” bimodal cell population in a time window when there is no instability in the set of rate equations used to describe the signaling events. The main difference between this study and the results we have discussed is that, in ref. 4, the observed bimodality disappears at long times. Another example is provided by T cell signaling. It has been proposed that dueling feedback regulation could underlie how antagonists shut off signaling in T cells. Experiments show a bimodal response for a downstream signaling product (Erk), with the proportion of “off” cells increasing as the number of antagonists becomes larger (15, 16). Stochastic simulations of a model of the T cell signaling network are in accord with these experimental observations (SI Text); i.e., bimodal distributions are the norm because of fluctuations, whereas the deterministic equations do not exhibit bistability in any parameter regime. However, we emphasize that a bimodal or “digital” ERK response in T cells could also result from important contributions from other molecular mechanisms (33).

We hope that the possibility of purely stochastic instabilities that lead to distinct cellular decisions will be broadly explored in the context of cell signaling processes by carrying out single cell assays for systems where the necessary and sufficient conditions we have described are naturally present or are engineered.

This research was supported by National Institutes of Health Grant 1 PO1 AI071195-01 and a National Institutes of Health Director's Pioneer Award (to A.K.C.).