|Home | About | Journals | Submit | Contact Us | Français|
There is increasing experimental evidence that cells can utilize biochemical noise to switch probabilistically between distinct gene-expression states. In this paper, we demonstrate that such noise-driven switching is dominated by tails of probability distributions and is therefore exponentially sensitive to changes in physiological parameters such as transcription and translation rates. Exponential sensitivity limits the robustness of noise-driven switching, suggesting cells may use other mechanisms in order to switch reliably. We discuss our results in the context of competence in the bacterium Bacillus subtilis.
Recent experiments indicate that cells can use noise in gene expression to switch probabilistically between distinct gene-expression states . For example, in the phenomenon of competence, the soil-dwelling bacterium Bacillus subtilis utilizes noise in gene expression to switch probabilistically between a vegetatively growing state and a competent state where a bacterium can take up exogenous DNA [2–6]. Noise is also likely responsible for the phenotypic heterogeneity found in mycobacteria  and the bistable behavior of B. subtilis in nutrient-rich media where genetically identical cells are found in two distinct phenotypes: motile swimmers and immobile chains . Recent experiments also suggest that a noise-driven switch is at least in part responsible for the ability of both yeast  and the human fungal pathogen Candida albicans [10, 11] to exist in two distinct phenotypes. The role of stochastic fluctuations in many other biological phenomenon such as the lytic-lysogeny decision during bactriophage infection is still the subject of current debate [12–14].
In this paper, we demonstrate that in standard models of noise-driven switching [15–26] many biologically relevant quantities, such as the switching rate, exhibit an exponential sensitivity to physiological parameters. The root cause of this exponential sensitivity is that switching properties are dominated by rare concatenations of biochemical events [23, 24]. Consequently, great care must be taken when calculating switching properties using coarse-grained models . In particular, we show that sources of noise often ignored in theoretical treatments, such as protein bursting [28–31], are extremely important for noise-driven switches. However, provided mRNA lifetimes are short, we show that switching can still be accurately simulated using protein-only models of gene expression.
The exponential sensitivity of switching may also have important biological consequences. For example, under many physiological conditions, B. subtilis cells can become competent only during a window of opportunity at the end of exponential growth [2, 3]. Thus, in standard models of switching, small changes in physiological parameter are likely to lead to large variations in the percentage of cells that become competent. It is then natural to ask if and how cells can achieve robust switching. We address this issue by discussing several alternative switching models that are more robust to changes in kinetics parameters.
Perhaps the simplest example of a genetic switch is a genetic network in which a protein induces its own expression (see figure 1). If this feedback is nonlinear, e.g. if multiple copies of the protein are required to induce expression, such a network may be bistable with two distinct steady states, a low-expression state where the protein is present in small numbers and a high-expression state where the protein is present in large numbers. Noise from the inherent stochasticity in biochemical reactions can probabilistically switch a cell between these two gene-expression states [32, 33]. Such a simple bistable network is at the core of many noise-driven systems in biology [8, 10, 11]. For example, in competence, the master-regulator protein ComK promotes its own expression through a positive-feedback loop with noise in gene-expression switching bacteria from a low-ComK vegetative state to a high-ComK competent state in response to stress. It is worth noting that under certain physiological conditions, the high-ComK competent state is actually transient. However, since we are concerned only with entrance into competence this will not be important to our arguments below.
The average number of mRNAs, , proteins, , in such a genetic network can be described by the deterministic differential equations
where fm(p) is of the Hill form and is given by
with α0m being the basal rate of mRNA transcription and the second term representing positive feedback. Equation (1) describes the change in the mean number of mRNAs due to mRNA transcription and degradation. Equation (2) describes the change in the mean number of proteins due to protein translation and degradation. In the bistable regime, these equations have two stable fixed points divided by a separatrix containing an unstable intermediate fixed point (see figure 1).
The deterministic differential equations (1) and (2) ignore fluctuations in mRNA and protein numbers. To capture the effects of these fluctuations, we employ the exact Master Equation for the probability, P(m, p), of finding m mRNA molecules and p protein molecules . Define shift operators, , which act on an arbitrary function g(m, p) according to and . In terms of these operators, the Master Equation describing our simple bistable switch is
The terms in equation (4) have the same biological meaning as the corresponding terms in equations (1) and (2). However, the Master Equation describes the full probability distributions P(m, p) and not just the mean numbers of mRNAs and proteins. In writing this equation, we have assumed that the noise in the system is due to intrinsic noise in transcription and translation  and we have ignored other sources of noise. In particular, we have modeled both protein and mRNA production with a simple Poisson model. This approximation may not hold in many biological systems where, in particular, transcription of mRNA may occur in bursts [35–37].
An important property of biological switches is the mean first-passage time (MFPT)—in our case, the average time it takes to switch gene-expression states. The MFPT is also inversely related to the mean switching rate for a population. As discussed in the introduction, in competence the MFPT is directly related to the number of cells in a bacterial population that become competent during a window of opportunity at the end of exponential growth [2, 3]. We have calculated the MFPT for the simple positive-feedback network described by equation (4) by performing 5000 independent Gillespie simulations for each set of parameters and averaging the first-passage times (FPTs). (The Gillespie algorithm is a dynamic Monte Carlo method that numerically generates trajectories of a stochastic Master Equation .) The FPT for each simulation was calculated by initializing the system at the low-protein-number fixed point and calculating the time the trajectory took to reach the deterministic separatrix. Figure 2 shows the MFPT for two different ratios of the mRNA and protein lifetimes, holding the mean protein number fixed. When τm/τp = 1/16 1, mRNAs are short-lived and typically many proteins are produced in a rapid burst from each mRNA molecule. In contrast, when τm/τp = 1, proteins are produced nearly continuously. In the two simulations, the mean number of proteins at the fixed points is identical. This is achieved by varying αp with τm so that = αpτm, the average number of proteins made per mRNA, is fixed.
Note that the MFPT is exponentially sensitive to the transcription-rate coefficient αm. Similar sensitivity is observed with respect to other parameters such as the protein-degradation rate (data not shown). The MFPT in figure 2 also shows a strong dependence on the ratio of mRNA and protein lifetimes though the mean protein number is held fixed. In particular, the MFPT is significantly shorter when proteins are produced in bursts (τm/τp 1). Bursty protein production increases protein fluctuations, and as a result lowers the MFPT for switching. This dramatic reduction in MFPT due to bursting has also been observed in the context of a dimerizing genetic switch . The increase in protein fluctuations due to bursts of production can be quite large. For example, the variance of the number of proteins for a fixed rate of mRNA transcription is given in the bursting limit by [28–31]
where = αpτm 1 is the mean burst size, i.e. the mean number of proteins made from each mRNA. This should be contrasted with the case where proteins are produced continuously and δp2/2 = 1/ (1 + )/.
In the bursting limit (τm τp), gene regulation is often approximated using protein-only models. In these approximations, time derivatives of the mRNA species are set equal to zero, and the resulting mean mRNA concentrations are substituted into the equations describing protein dynamics. For our network, such a procedure yields the reduced equation
It is worth noting that since these approximations fail when τm ≥ τp, switching properties cannot be accurately described by protein-only models in the continuous protein synthesis regime.
Even in the bursting regime, commonly used protein-only models often fail to accurately reproduce switching kinetics since these models neglect or oversimplify noise sources such as protein bursting . Stochastic fluctuations are often included in protein-only models by assuming that proteins are produced probablistically one at a time. This assumption leads to the master equation
We have calculated the MFPT for switching between the low-expression and high-expression states for this Master Equation using the Gillespie algorithm and the result is shown in figure 3. Note that the MFPTs calculated in this way can be nearly four orders of magnitude longer than those calculated using the full mRNA-protein model, demonstrating that protein bursting significantly increases the rate of switching.
In fact, it is straightforward to include protein bursting in protein-only models in the bursting regime . Since mRNA degradation and protein translation from an mRNA molecule are independent Poisson processes, it can be easily shown that the protein-burst size b (defined as the number of proteins produced from a single mRNA molecule) is geometrically distributed with mean burst size = αpτm and distribution
We consider two different approximations for incorporating bursting in a protein-only model. In the first, once an mRNA is produced, translation proceeds deterministically and proteins are produced. In the second, b proteins are produced from a single mRNA molecule with the geometric distribution, G(b). Both these approximations assume τm τp since all the proteins from a single mRNA molecule are considered to be produced instantaneously. The MFPTs computed within these two approximations are shown in figure 3. As for the full mRNA-protein model, the protein-only MFPTs are exponentially sensitive to changes in the transcription-rate coefficient αm. The MFPT computed using the geometric distribution of protein-burst sizes agrees well with the MFPT from the full mRNA-protein model. However, the protein-only model where proteins are produced with a fixed burst size is inadequate for calculating the MFPT even though this approximation correctly yields equation (5) for the variance in protein number at a fixed rate of mRNA transcription. Note that the MFPT computed within the geometric-burst-size approximation is smaller than that computed using the fixed-burst-size approximation. This is expected since geometrically distributed burst sizes increase the fluctuations in protein production.
A qualitative understanding of our simulation results can be obtained from the Fokker–Planck approximation to the Master Equation (7). Gaussian approximations such as the Fokker–Plank and Langevin equations generally fail to accurately capture switching properties since they underestimate rare fluctuations , but nonetheless they allow us to identify the regime of exponential sensitivity. The Fokker–Plank equation that approximates our genetic network is 
where and . In deriving (9), we assumed that the state of the system is described by a single reaction coordinate, protein number, and have ignored fluctuations such as transcription-factor binding and unbinding. Additionally, we have made the simplifying approximation that proteins are produced in bursts of size . We can calculate the MFPT using equation (9) and obtain 
where we have defined an effective potential U(p) using the equation and where A(B) is the mean number of proteins at the low-expression stable fixed point (intermediate-expression, unstable fixed point). The result is derived assuming the constant-diffusion approximation where we replace a2(p) by its value at the low-expression fixed point a2(A) and performing a saddle-point approximation. As expected, this approximate Fokker–Plank equation correctly determines the regime of exponential sensitivity but does not quantitatively reproduce the MFPT. Note that the MFPT resembles the FPT for barrier crossing by a particle moving in a potential U(p) with the effective diffusion constant D = a2(A) and the barrier height UBA . The analogy between barrier crossings and biochemical switches was discussed earlier in [19, 20]. We have extended this analogy to include protein bursting and shown that the burst size changes both the barrier height and the effective diffusion constant. When UBA/D 1, i.e. when the barrier height is large, the MFPT is exponentially sensitive to changes in the kinetic parameters that determine U(p) and D. Thus, exponential sensitivity indicates that switching is dominated by the tails of a probability distribution and the central limit theorem no longer applies. Thus, aspects of gene regulation, such as protein bursting which may not be important when considering small fluctuations around mean protein levels, may still play an important role in governing switching rates.
There is increasing evidence that the dynamic properties of bistable genetic switches are extremely sensitive to the details of gene regulation. This work, as well as others [25–27], shows that small changes in parameters can lead to large changes in the dynamics of switching. We have argued that this sensitivity occurs because the dynamics of switching depends exponentially on the effective barrier height and the effective diffusion coefficient for switching. Such exponential sensitivity has important theoretical implications for coarse-grained descriptions of biochemical networks . For example, many sources of noise, such as mRNA bursting [35–37] and protein dimerization [16, 25, 26], likely play important roles in determining switching rates in some biological systems. Thus, great care must be taken theoretically when integrating out sources of noise . This work complements and extends earlier theoretical studies by showing that when mRNA lifetimes are short, protein fluctuations can be correctly incorporated using protein-only models. In addition, it provides a simple understanding of exponential sensitivity using intuition from the statistical mechanics of barrier crossings.
The exponential sensitivity of genetic switches also has important physiological consequences. In recent experiments, Maamar et al [2, 3] showed that protein bursting is a major source of the noise driving the induction of competence. Furthermore, they were able to change the noise-driven switching rate (the inverse of the MFPT) in the B. subtilis competence system, while keeping mean protein numbers fixed, by manipulating transcription and translation rates. This suggests that evolution can tune switching rates while keeping mean protein levels fixed. (Note, however, that switching rates are also likely to be influenced by additional factors such as temperature, pH, fluctuations in the number of cellular components important for transcription and translation such as polymerases and ribosomes, and extra-cellular factors such as chemicals present in the environment.)
By our analysis, switching—because it is like barrier crossing in statistical mechanics—is likely to be exponentially sensitive to all these factors. Indeed, this may be one of the causes for the large variations observed in experiments on the number of swimming cells and chains in a cultures of B. subtilis . However, it is then natural to ask if, and how, cells can achieve robust switching. For example, robust switching may be achieved using alternative switching mechanisms such as all-or-none events that switch the cell between the low- and high-expression states . Such an all-or-none event could correspond to the disassembly of an entire repressor complex bound to the promoter of a gene, allowing multiple mRNA molecules to be made rapidly. If switching is dominated by these all-or-none events and not the tails of probability distributions, then switching is no longer necessarily exponentially sensitive to changes in parameters. This scenario may be relevant to competence since a repressor complex is known to be important in regulating the master competence regulator ComK . Alternatively, cells could switch robustly by actively tuning noise to control switching rates. Recent experiments suggest that this may indeed be the case in mycobacteria . A final possibility is that large variations in the switching rate may be acceptable or even desirable for cells living in certain environments [44, 45]. Thus, cells may exploit the exponential sensitivity of noise-driven switching to increase their fitness. It will be interesting to better understand the ecological and evolutionary implications of noise-driven switching.
We would like to acknowledge D Dubnau, H Maamar, A Raj, M Elowitz, G Suel, S Goyal and A Sengupta for useful discussions. This work was supported in part by US National Institutes of Health (NIH) grants R01 GM073186 and PSO GM071508, and by the Defense Advanced Research Project Agency (DARPA) under grant HR0011-05-0057.