Relatively slow stochastic gene-state switching in the presence of positive feedback significantly broadens the region of bimodality through stabilizing the uninduced phenotypic state

Within an isogenic population, even in the same extracellular environment, individual cells can exhibit various phenotypic states. The exact role of stochastic gene-state switching regulating the transition among these phenotypic states in a single cell is not fully understood, especially in the presence of positive feedback. Recent high-precision single-cell measurements showed that, at least in bacteria, switching in gene states is slow relative to the typical rates of active transcription and translation. Hence using the lac operon as an archetype, in such a region of operon-state switching, we present a fluctuating-rate model for this classical gene regulation module, incorporating the more realistic operon-state switching mechanism that was recently elucidated. We found that the positive feedback mechanism induces bistability (referred to as deterministic bistability), and that the parameter range for its occurrence is significantly broadened by stochastic operon-state switching. We further show that in the absence of positive feedback, operon-state switching must be extremely slow to trigger bistability by itself. However, in the presence of positive feedback, which stabilizes the induced state, the relatively slow operon-state switching kinetics within the physiological region are sufficient to stabilize the uninduced state, together generating a broadened parameter region of bistability (referred to as stochastic bistability). We illustrate the opposite phenotype-transition rate dependence upon the operon-state switching rates in the two types of bistability, with the aid of a recently proposed rate formula for fluctuating-rate models. The rate formula also predicts a maximal transition rate in the intermediate region of operon-state switching, which is validated by numerical simulations in our model. Overall, our findings suggest a biological function of transcriptional “variations” among genetically identical cells, for the emergence of bistability and transition between phenotypic states.

1 Stochastic models of central dogma 1.1 Simplest mechanism: a single DNA state A single DNA waits for an exponential distributed time window with parameter k 1 until it gives rise to a piece of mRNA, and then the life time of each mRN A also obeys an exponential distribution with parameter d 1 . Within the life time of each mRN A, protein P is being synthesized by a Poisson process with intensity k 2 , and at the same time, each molecule of P degrades with parameter d 2 .
During each single burst, the distribution of synthesized protein is geometric with parameter q = k 2 k 2 +d 1 . More precise, Its average is q/(1 − q), which is regarded as the burst size. This fact was first proved by [12] and recently confirmed by experiments [14]. And the master equation for the population distribution of protein is dP (n, t) dt = k 1 ( n j=1 G j P (n − j, t) − qP (n, t)) + d 2 ((n + 1)P (n + 1, t) − nP (n, t)), whose stationary distribution is just the negative binomial [16] P ss (n) = b n (1 + b) a+n Γ(a + n) Γ(a)n! , where a = k 1 d 2 is the burst frequency per cell cycle and b = k 2 d 1 = q/(1 − q) is just the burst size. a and b are typically determined through the fitting of the stationary distribution that measured.
Interestingly, it is not the real burst size observed if you track the protein in a single cell, because one could not observe G 0 (although P ss (n) could be observed), hence this kind of real burst sizeb = b/(1 − G 0 ) = b + 1.

Two-state model
See Fig S1. Assume the parameters for the exponential distributed "on" time of the DNA and mRNA are β, d 1 respectively, and within that "on" time, the mRNA and protein are being synthesized, by another two Poisson process with intensities k 1 and k 2 repectively.
It is easy to prove that the number of mRNA synthesized per single "on" period of DNA and the number of newly synthesized protein per single "on" period of mRNA both obey geometric distribution with parameter θ 1 = k 1 β+k 1 and θ 2 = k 2 d 1 +k 2 respectively. Therefore, the distribution of newly synthesized protein number during a single "on" period of DNA is This is a modified geometric distribution with parameter θ 1 and θ 2 . In this case, we don't think experimentally one can tell the difference between this and a real geometric distribution since all the difference is in the distributions at n = 0 and n = 1, overwhelmed in the experimental noise.
Consequently, the averaged value of this modified geometric distribution is The master equation for the population distribution of protein under the bursty condition is

Its stationary distribution is just the negative binomial
is the burst size. Notice that b could not be less than k 2 d 1 which is the mean translated protein number per mRNA molecule.
The real burst size if you track the protein in a single cellb = n /(1 − P (0)) = α, k 1 , then the two-state model can be approximated by the simplest model with only one gene state, since the rate limiting step to synthesize a new mRNA is the switching from the inactive state to the active state. In this case, the equallyk 1 = k 1 α α+β . Then the burst frequency and size in the two-state model reduce to those in the simplest model.

Small burst
During each single small burst, the model could be reduced to in which OR and O * R are the gene states with the repressor bound to both operators or only the weaker one. Here Therefore, the population distribution of protein is the negative binomial with param- where O * R denotes the partial dissociation state of DNA. Here , and the frequency could be computed by the first passage time starting from the state OR to O, i.e. T = Hence both the large burst frequency and size increase with intracellular inducer concentration, but the frequency will not change that much while the size would increase significantly.
When the inducer concentration is high, then Fig 1C in the main text could be approximated by an alternative pathway Hence the large burst size is nearly constant, and the frequency will slightly increase with intracellular inducer concentration with the However, as we have mentioned in the main text, the bursts would be much more indistinct when the intracellular inducer concentration is quite high.
For intermediate intracellular inducer concentration, the general expression of large burst size can be well approximated by

Calculation of parameters
First of all, based on Fig 1 of the main text, we assume that the system satisfies the thermodynamic constrain: Then we could get that where K 3 = r 3 r −3 . Therefore, Lac operon copy number and production/degradation rates Kennell measured the translation rate of LacZ as 18.8min −1 [30] and (Kennell and Riezman, 1977), which is also consistent with [14]. And since LacZ is nearly expressed 10 times higher than LacY, hence we set transcription rate of LacY as k Y = 1.8min −1 . The mRNA half-life time for LacA is 55 seconds according to (Kennell and Riezman, 1977), hence γ M = ln 2/(55/60) = 0.756min −1 . Also Tsr and permease should have similar rates. In fact, Choi, et al. [3] did a direct comparison of them (LacY-Venus replaced by Tsr-Venus) and find similar expression levels.
The cell-doubling time is about 60 minutes, hence the dilution rate for the stable protein TMG γ I = 0.012min −1 . And the Lac permease degradation rate is about 0.01min −1 , then γ Y = 0.01 + 0.012 = 0.022min −1 .
The in-vivo transcription rate of LacY k M is more complicated. Kennell's estimate of k M in vitro is closer to 18min −1 . However, for the fully induced cells, the probability of open operon is nearly 1 and the mRNA should approximately k M /γ M ≈ 23.8, which is 2 − 3 folds of that observed recently (almost every gene has an mRNA copy number of less than 10 per cell in [20]). This value also varies: for instance, k M ≈ 9min −1 in (Chung and Stephanopoulos, 1996) determined from theoretical fitting, and even simple physical or chemical data published by different groups will differ by factors of 2 − 10. So in order to become more quantitatively consistent with recent experiments, we here set k M = 8min −1 .
Parameters for operon stochastic dynamics According to Elf, et al.'s kinetic measurement in vivo (Elf et al., 2007), there are about 1−3 repressers per cell, and each repressor will spend about 5−6 minutes to rebind, hence R T = 2molec. and the association rate r 1 ≈ 0.2molec. −1 min −1 . Referred to the in-vitro measurements of Table I and III of [31], the binding constants for repressor and inducer √ K is about 10µM , and that for repressor/operon and inducer 1 √ K 3 is several mM (Matzke et al., 1992). However, in the in-vivo experiments, it seems that K = 100molec. −1 µM 2 is too low since the uninduced cells still should exist when the intracellular inducer concentration increases to nearly 50 − 100µM . The inconsistence of the in vitro and in vivo measurement might result from the much more complicated environment inside a living cell. For example, the DNA-protein interactions are highly salt-dependent, and the flexibility and persistence length of DNA in vivo does not seem to match in vitro values, probably because there are other proteins binding to the DNA that can cause bends, etc. Therefore, we modify K to be 2500molec. −1 µM 2 and set K 3 = 8 × 10 −6 µM −2 . Then according to the thermodynamic constrain, we can have K 1 K 4 = 1 KK 3 = 50molec.. Furthermore, the mRNA level/ protein level is nearly k Y /γ Y ≈ 100, which is consistent with experimental observations for 1000 induction ratio (for uninduced cells, mRNA=0.01 and protein=1; for induced cells, mRNA=10 and protein=1000). It also implies the proba-bility for unrepressed operon p O as the intracellular inducer concentration is low (≤ 50µM ) the burst frequency Inducer uptake rate Without feedback, the intracellular level will likely be equal to the extracellular, within 10 minutes at most , hence we could set c ≈ ln 10/10 = 0.23min −1 .
And finally k l = 0.25min −1 is determined to make the bistability range of our simulation roughly consistent with the experimental observation.
All parameters are summarized in S1 Table. Fig S2 Copy-number distributions for the permease protein in wild-type cells. We compare the copy-number distribution of permease with different extracellular concentration of inducers I e and show that the I e range of the bimodal distribution is much more broader than that predicted in the deterministic bifurcation diagram( Fig. 2A in the main text).    Table S1 Values of kinetic parameters in the fluctuating-rate model.