Bimodality in E. coli gene expression: Sources and robustness to genome-wide stresses

2.1. Empirical distributions of single-cell protein levels

A previous work [37] investigated single-cell expression levels of genes during the exponential and the stationary growth phases as a function of promoter σ factor preference. That study proposed a predictive model of changes in expression level during growth-phase transition as a function of the promoter sequence. Whereas mean and variability of single-cell expression levels was investigated, the modes of the single-cell distributions of expression levels were not studied. Subsequent data analysis revealed that two genes (metK and pgi) exhibited distributions of single-cell expression levels that could potentially be bimodal.

To find additional genes with significant chances of exhibiting bimodal dynamics, we considered the known mean and noise of single-cell protein numbers reported for 1018 genes of E. coli during standard growth conditions [24] (Fig 1A). Then, we selected 56 genes along a large portion of the space state of the mean over the noise in single-cell protein levels. Of these 56 genes, 70% were selected for having relatively high noise (i.e., are above the best fitting line in Fig 1A) reported in [24]. The remaining genes were selected for the opposite reason, for comparison. Relevantly, 6 out of the 7 genes found to have bimodal distributions are located above the line. Fig A in S1 Text shows how our measurements of these genes’ expression levels compare to the levels reported in [24].

Next, using flow cytometry, we obtained single-cell distributions of protein levels for each of the 56 genes, in standard growth conditions (Section 1.1 and 1.2 in S1 Text). We found that 7 of the 56 genes exhibited bimodal distributions, according to at least one of the criteria in Section 1.4 in S1 Text (Fig 1B).

Fig 1C shows the best fits to the merged single-cell distributions of protein levels for each of the seven genes and corresponding strains (Section 1.5, histograms plotted in Fig B in S1 Text). Specifically, for each gene, the distribution results from merging the data from three biological replicates. Also shown is the data corresponding to the wild-type (WT) strain. Data for each replicate is provided in S1 Text. From here onwards, for simplicity, we refer to the fits as “distributions” as well.

Visibly, the distributions differ significantly between genes. For example, arguably, whereas for elaB, metK and pyrH, the peaks of the two modes are very distant from each other, large distance is not evident for manX, paaX, and lldD. Also, whereas in elaB, metK and pgi, the “choice of mode” by each cell is highly asymmetric, for the other genes the distribution looks symmetric (same chances for both modes). Meanwhile, the distributions from each replicate are shown in Figs C and D in S1 Text. Visibly, there is little difference between replicates, suggesting that the bimodality is stable under stable conditions.

Notably, all distributions are asymmetric (i.e., right sided), with the ‘proportion’ of the total areas related to the right mode being larger than the proportion related to the left mode. Whereas ‘higher peaks’ in the left modes could potentially be explained by the boundary at zero, this asymmetry likely cannot.

One could hypothesize that the low-fluorescence intensity values of the first mode of the distributions in Fig 1C result from unhealthy or dead cells, rather than from healthy cells with low expression levels. To test this hypothesis, we compared the normalized single-cell distribution of fluorescent levels of each strain (FITC-H, commonly used as proxies for protein numbers) with the corresponding normalized single-cell distributions of pulse width, SSC-H, and FSC-H, respectively (parameters commonly used as proxies of cell size) [38]. From Figs E–G in S1 Text, we found little to no similarity between the distributions. Visibly, the peaks are not aligned, and the mode weights are not similar. Only rarely is there similarity, but never in all biological replicates of a given strain. These results are line with [24], which reported that these strains are as healthy as the other strains of the YFP library. This is supported by the density scatter plots between the single-cell protein fluorescence levels and the pulse width, SSC-H, and FSC-H, respectively (Figs H–J in S1 Text). In these, while some correlations are observed between YFP intensity and the cell morphology parameters (which is expected, since cells often regulate their protein concentrations rather than total numbers), no concurrent bimodality is evident in these channels.

Also noteworthy, the shapes of the first modes of the distributions differ widely between strains (Fig 1C), e.g. in mean values and in fraction of cells in that mode. We use this as evidence that the fractions of the cells with low or absent fluorescence levels differ and the levels themselves also differ, suggesting that they do not result from the same cause. Moreover, from microscopy images, we did not find differences in cell morphology (i.e., length, width and area) that could explain the bimodal gene expression such as asymmetric in cell division (Fig K in S1 Text). Also, we found no biased spatial distributions of YFP tagged proteins (Fig L in S1 Text).

Given the above, we conclude that the bimodal distributions of single-cell protein fluorescence levels cannot be explained by health or morphology features of these strains, which would cause a significant fraction of cells to behave similarly to wild-type cells. Instead, the bimodal distributions likely result from each cell population being able to express the tagged proteins at two distinct rates (which differ between strains).

In this regard, from Fig 1C, in some strains, particularly pyrH, paaX and manX, cells composing the first mode of distribution appear to have low expression levels (Table A in S1 Text). In the other strains, the cells exhibit fluorescence levels similar to the WT strain, suggesting that the genes are not active.

2.2. Stresses responsiveness and robustness

Next, we investigated the robustness of the bimodality of each of the seven genes to five different stresses. Noteworthy, these stresses have direct, as well as indirect effects on gene expression. For example, cold shock is expected to have multiple effects, from influencing DNA supercoiling to cellular energy, etc. [10].

First, we studied the effects of the transitioning from exponential to stationary growth. During this transition, cells actively increase the number of RpoS (σ³⁸), activating ~10% of the genome [39–42]. Moreover, since RNA polymerase (RNAP) numbers are lower than σ factor numbers, many genes not recognized by σ³⁸ become less expressed [43–45]. Meanwhile, lldD, manX, paaX, and pyrH are σ⁷⁰-dependent [46] and, thus, are expected to be less expressed during this growth phase.

Results are shown in Fig 2 (orange distributions). Compared to the control condition (green distributions), we find that most genes (lldD, manX, paaX, pyrH) were strongly repressed, as expected. Namely, their distributions of expression levels appear to become unimodal, in the region originally occupied by the lower peak of the bimodal distribution. This would imply that the lower expression “state”, already existent during standard growth, now became dominant.

Another possible explanation is that the system remains bistable, but its two states do not differ sufficiently in expression levels to be distinguishable. This would occur if, e.g., the differences between the mean protein numbers of the two states are smaller than the variabilities.

Meanwhile, metK and pgi are not repressed, maintaining bimodality, which is consistent with them being dual σ⁷⁰ and σ³⁸-dependent [37]. Finally, elaB was also repressed during the growth phase transition, albeit being reported to be σ³⁸-dependent [46]. Relevantly, high expression of this gene during exponential growth was also reported in [24]. This may explainable by its negative regulation by rob, whose expression might be enhanced during the stationary growth phase [47]). Another explanation would be the dual regulation by oxyR [46].

Second, we perturbed the processes regulating DNA supercoiling. We used Novobiocin, an antibiotic that binds to the GyrB subunit of DNA gyrase, which is needed to resolve positive supercoiling buildup (PSB), blocking ATPase activity [48,49]. We expect the responsiveness of the genes, and, thus, changes in their bimodality, to differ with their supercoiling sensitivity. A past study reported that lldD, manX, metK, and pgi are highly sensitive to positive supercoiling buildup during standard growing conditions [10]. Meanwhile, from Fig 2, these are the genes whose activity was, comparatively, more reduced, with pyrH being the only one maintaining bimodality.

Third, we subjected cells to cold shock, which hampers gene expression in multiple ways, including the disruption of PSB [10,50,51]. In addition, cold shock is expected to influence transcription directly, particularly the initiation process [50,52]. Other effects include hampering diffusion in the cytoplasm [53] and indirectly causing significant decreases in cellular ATP levels [10]. As such, effects on individual genes are not easily predictable. To evaluate the consistency of the results we compared them with the results under Novobiocin. Past studies showed that nearly half of cold shock repressed genes are also highly responsive to gyrase inhibition [10]. In agreement, we find that the genes show similar responses, except pgi and metK. For example, the genes remaining bimodal under cold shock and under novobiocin are elaB and pyrH. Meanwhile, only metK and pgi remains bimodal under cold shock alone, and only manX remains bimodal under novobiocin alone.

Fourth, we hampered transcription rates. For this, we used rifampicin, an antibiotic that binds to the β sub-unit of RNAP [54] and prevents RNAP from escaping beyond 2–3 nucleotides away from the TSS [55,56]. In agreement, we find that all genes have weaker expressions than during standard growth.

Finally, we targeted translation by adding Streptomycin, an antibiotic that binds to the ribosome’s 30S subunit [57] and affects the ability of ribosomes to decode the RNA [58], which results in the formation of non-functional proteins [59]. Given the strong correlation between RNA and protein numbers in E. coli [60], one would expect that this antibiotic should cause similar effects to rifampicin, which agrees with what we observe (Fig 2).

Overall, qualitatively, we observed a wide variety of responses to the stresses. Different genes lose their bimodality to different perturbations, with 4 of 7 of the genes losing their bimodality to all perturbations. Interestingly, no perturbation forced all genes to lose bimodality. These results suggest that the genes of E. coli with bimodal gene expression studied here exhibit responsiveness to some stresses but also robustness to other stresses targeting gene expression.

2.3. Changes in the shapes of the bimodal distributions with the stresses

To quantify how the perturbations targeting gene expression affected the shape of the bimodal distributions, we considered the following shape parameters:

• (i)

  Distance between peaks 1 and 2 (|pk₁ - pk₂|) as a measure of the difference between the two possible ‘states’ of gene expression (high and low expression).

• (ii)

  Difference between peak heights (i.e., between the PDF values at the peaks’ positions, PDF₂ - PDF₁) as a measure of the asymmetry in the proportions of subpopulations in different states.

• (iii)

  Fraction of cells whose fluorescence level is in between the two peaks (n_overlap), is used as an indirect measure (positively correlated) of the frequency with which cells may be switching between states and/or for the stochasticity of either (or both) states. However, it needs to be considered that, instead, they may result from broad distributions with complex shapes.

To compare between genes, we next defined relative shape parameters (illustrated in Fig M in S1 Text). We calculated the relative distance between peaks (d) by normalizing the above difference with the entire range of the distribution (Eq. 1). Likewise, we calculated the relative difference between peaks height (h) by normalizing the difference by the maximum height of the distribution (Eq. 2). Finally, we obtained a relative overlapping region (o) from the number of observations between the two peaks divided by the total amount of observations (n_total) (Eq. 3):

Noteworthy, these parameters should only be considered if the distributions are bimodal. Fig N in S1 Text shows d, h, and o, for each distribution in Fig 2 (when a distribution for a given perturbation was not classified as bimodal, it was discarded). Visibly, the perturbations alter each of the features characterizing the shapes and do it differently between genes and between perturbations. These results suggest that it should be possible to externally fine-tune the bimodal distributions since d, h, and o exhibit many different values within a significant range of values rather than, e.g., being binary-like, which would hamper fine-tuning.

2.4. Bimodality lost due to cells entering stationary growth phase can reemerge if cells reenter the exponential growth phase

We explored whether the bimodality, after becoming undetectable, can again become detectable. We tested this in the case of the transition to stationary growth phase, since this perturbation (along with streptomycin) was the one causing most genes to lose bimodality. For this, we tracked the same cell culture for three generations, while placing the cells in new media once they reach the late stationary growth phase at each generation. Specifically, at each generation, we measure gene expression by flow-cytometry at four growth phases: early and late exponential, followed by early and late stationary growth (Section 1.2 in S1 Text). With respect to time, the consecutive phases are separated by ~50 min, ~150 min, ~200 min, respectively.

In Fig 3A–3G, for each strain, we plotted the single-cell distributions of expression levels at each of the four growth phases (early exponential, exponential, early stationary, and stationary), for the three consecutive generations. First, for all genes and in all generations, the shapes of the distributions differ significantly between growth phases. Second, for most genes and phases, the behaviors observed in the first generation (unimodality and bimodality) were also observed in the next generations. For example, for manX, the bimodality observed in the late exponential phase (and only in that phase) during the first generation, is again observed in that phase (and only in that phase) in the subsequent generations.

Nevertheless, the shapes of the distributions do, in some cases, differ significantly. There are many possible reasons for this. One explanation may be in the methodology used for the measurements, in that the original media from which we removed cells into fresh media at each generation, was not renewed overtime. Finally, arguably, since the changes in the shape of the distributions occurred over short periods of time (implying that many cells changed state), these modifications can be used as evidence that the genes can “switch state” (between high- and low-expression rates) relatively frequently.

2.5. Evidence for hysteresis

Since bimodality in gene expression reemerged when cells were again placed in fresh media (albeit with differences in the shapes of the distributions), we searched for potential hysteresis. Its existence could be evidence for the influence of the systems’ history in its response. This can occur in gene expression patterns across generations when, e.g., there were transcriptome changes during the first adaptation that still remain during the second adaptation [61] and has been previously linked to multistable networks [62].

We tracked the bimodal distributions shape parameters d and o, to search for hysteresis, of genes with bimodal behaviors in more than one growth condition (elaB, metK, pgi, pyrH). We did not track the parameter h, as it depends on the fraction of cells extracted between generations and, thus, on the experimental procedure. Nevertheless, as we inoculate the same volume from one generation to the next, this fraction should be relatively similar.

From Fig O in S1 Text, the phase transitions have distinct paths for d and o, implying that the distribution depends on its history, in agreement with the existence of hysteresis. Nevertheless, in most genes, there is no ‘loop closing’. In these cases, there is no evidence that the cell population returned to its original state, implying the potential existence of only 'partial hysteresis'. Moreover, there are other potential explanations for the observations that will need to be considered in future works.

2.6. Bimodality in single cell promoter activity levels

To investigate which mechanisms might be causing bimodality in the expression levels, we next made use of a library of transcriptional fusions of GFP of promoters of E. coli [63]. Each promoter fusion is on a low-copy plasmid, and it controls the expression of GFP alone (Section 1.2 in S1 Text). The library includes plasmids with the promoters of 6 out of the 7 genes studied, namely, elaB, manX, metK, paaX, pgi, and pyrH, respectively.

We measured the single-cell fluorescent levels of each of the six strains under standard growth conditions. Then, we compared them with the strains of the YFP library (in Fig 1C). For this, we first considered that GFP and YFP differ in fluorescence intensity. Thus, we ‘mapped’ the data from the promoter fusion library into the data from the YFP strain library. For this, for each of the 6 promoters measured, we calibrated 4 data points of the single-cell distributions of fluorescence levels using the strain of the GFP library to the corresponding data points using the strain from the YFP library (Fig 4A). The points were the fluorescence levels of the cell with lowest signal, the cell with highest signal, the cells at the mean of the weak mode, and the cells at the mean of the high mode. Finally, using these points, we estimated the calibrated distribution by linear interpolation. We plotted the calibrated distributions in Fig 4B. Using the peak detection method (Section 1.4 in S1 Text), we find that all calibrated distributions from the GFP library are bimodal, confirming the conclusions using the strains of the YFP library.

Next, we considered that: (i) The promoters controlling GFP expression are in low-copy plasmids, not in the chromosome. As such, being at the chromosome is not essential for the generation of bimodality in single-cell protein levels. Given this, we argue that supercoiling buildup is not an a essential source of bimodality of the genes studied, since the plasmids used lack topological constrains that could allow for strong positive supercoiling buildup [64] (ii) The plasmids only express GFP, not the native protein, and many plasmids of the GFP library do not exhibit bimodal distributions. Thus, events during transcription elongation or termination of the native proteins in the YFP library are not necessary to generate the bimodality observed (albeit these events could influence the shapes of the single-cell distributions of fluorescence levels). (iii) Finally, for the same reasons, we argue that the processes of RNA and protein degradation of the native proteins studied using the YFP library are not necessary to generate the bimodality, as the bimodality still emerges in the corresponding cells of the GFP library.

Overall, since the only structure common to both libraries are the promoter regions, where transcription initiation occurs, we conclude that the bimodality in the single-cell distributions of protein expression levels in cells of both libraries is generated by events that occur at the promoter region.

2.7. Reduced stochastic model of bimodal gene expression

We developed stochastic models of bimodal gene expression to mimic the responsiveness to the perturbations studied above. We use the models to explore the state space, in order to identify influential parameter values and find saturation points.

First, we define a ‘detailed’ illustrative model based on the experiments conducted, specifically, cold shock, novobiocin, stationary growth-phase transition, rifampicin, and streptomycin (Fig 5). Shortly, the model assumes σ factor regulation (pink area, Fig 5) and a two rate-limiting-step process of transcription initiation: a step modelling promoter binding and finding of the transcription start site by RNAP holoenzymes, followed by a step modeling the processes that occur until promoter escape, such as the open-complex formation (blue area, Fig 5) [65,66]. Whereas all steps are reversible [66], we do not model reversibility, as we lack empirical data on in vivo parameter values for the events. These steps would increase noise in gene expression, but not mean rates.

Critically, the model assumes that the promoter can switch between two “transcription states” that differ in the rate of RNA production (blue area, Fig 5) [52]. This (and only this) allows one gene to have two distinct RNA production rates at different time points. For simplicity, we implement this difference in the rate of promoter binding by RNAP holoenzymes. However, in some genes with bimodal expression levels, it could be caused by different rates of promoter escape instead (or by differences in both rates). Note that many mechanisms could tune these rates, such as TFs. Here, we do not model explicitly the mechanisms making this possible.

Next, the model assumes that the promoter can be locked due to PSB and, thus, be influenced by the antibiotic novobiocin, as well as by cold shock, among other phenomena (green area, Fig 5) [10]. The effects of PSB and novobiocin can be set to differ with the promoter state. For simplicity, as a past study suggests that this is the step most influenced by cold shock, we do not model its influence on other steps [10].

The model also considers rifampicin. Separating transcription into two rate-limiting steps allows modelling that rifampicin specifically hampers promoter escape alone [55,56] (red area, Fig 5).

Meanwhile, translation, and the effects of streptomycin, are modeled by two separate (competing) sets of reactions. The first (translation, orange area), produces functional proteins. The second (yellow area), models the effect of streptomycin, namely the distortion of Ribosomes, which then produce nonfunctional proteins [57–59].

Finally, the model includes RNA and protein decay events (grey area). These decays account for degradation as well as dilution due to cell division.

From the detailed model in Fig 5, we then developed a reduced stochastic model of bimodal gene expression. This is necessary, since we lack empirical data on several rate constants necessary to implement the detailed model. This model, shown in Fig 6A, also allows promoters to have two states (low ‘L’ and high ‘H’ transcription rates). Switching between states occurs at the rates given by constants k_H and k_L, respectively. Transcription occurs at the corresponding rates constants, k_bind^H, k_esc^H, k_bind^L and k_esc^L. Once transcribed, RNAs are translated at a rate k_tr independent from the promoter state. Finally, RNAs and proteins decay due to degradation and dilution due to cell division (k_d^RNA or k_d^P), also independently from the promoter state.

Using this model, the effects of rifampicin are simulated by reducing the rates of promoter escape [55,56], k_esc^L and k_esc^H, depending on the promoter state. Meanwhile, the effects of shifting to stationary growth phase are modeled by tuning the numbers of RNAP.σ⁷⁰ and of RNAP.σ³⁸ (for simplicity, we assume a promoter with σ⁷⁰ preference). Also, to study the effects of different σ factor preferences, we tuned the rate of RNAP-promoter binding k_bind^L or k_bind^H. Further, the effects of Novobiocin are modeled by changing k_L or k_H, which tunes the expected time that a promoter is in each state [64], while the effects of streptomycin are modeled by reducing the rate of translation, k_tr [57]. Finally, we modeled the effects of changing decay rates of RNA and proteins, to simulate changes in cell division rates. In addition to this, one could model the effects of cold shock by tuning according to different degrees for all rates above (not done here as empirical data are not abundant).

From genome-wide measurements, we obtained reference values for each rate constant of the reduced model (Table B in S1 Text). The bimodal distribution of single-cell protein numbers obtained from simulations using these reference values is shown in Fig 6B. The distribution is bimodal, and exhibits a weak asymmetry, with a higher peak in the first mode.

It is noteworthy that the distribution produced by the model in Fig 6B has far less variability than any of the empirical distributions in Fig 4B. This is expected since the only source of variability in the model is the stochastic nature of the chemical reactions. Meanwhile, the sources of variability of the empirical distributions include not only the stochastic nature of the chemical reactions, but also several other sources, such as differences in cell components between individual cells, variability in the moments when the DNA replicates, variability in the asymmetries in partitioning of components during cell division, measurement uncertainties, among other.

Finally, we note that, in the model above, unimodal distributions emerge when the dynamics of two, distinct transcription events are sufficiently similar to make these events dynamically indistinguishable. For comparison, we next considered models that, instead, have the two-state system reduced into a single-state system. As such, these models’ reactions are the same as the reactions in Fig 6A (except that one of the two transcription reactions no longer exists). We considered three such models, which differ in the expression rate of the remaining transcription state, equaling: (i) the rate of the state (of the two original ones) with lower-expression, or (ii) the average of the rates of the two original states, or (iii) the rate of the original state with higher expression.

Fig P in S1 Text shows how the original bimodal distributions shift to unimodal in all three cases, after changing the transcription reactions of the system. All models take a similar time to exhibit unimodal single-cell distributions of protein numbers. Meanwhile, panel A of Fig Q in S1 Text shows that, as expected, these models differ in mean single-cell protein numbers, with only model (ii) matching the original model. For this, and more reasons, they also differ in the CV² of single cell protein numbers (panel B of Fig Q in S1 Text).

Interestingly, for example, model (i) is in good agreement with the observed dynamics of manX and paaX under many of the perturbations (Fig 2). Meanwhile, model (ii) is in good agreement with the observed dynamics of, e.g., metK subjected to novobiocin (Fig 2). Finally, model (iii) describes well the dynamics of, e.g., pgi under cold shock.

2.8. Exploring the robustness of the bimodality using model simulations

We explored the reduced model by tuning, independently, each rate constant in Fig 6A, assuming that they can, at most, become 10× larger or 10x smaller than the reference value in Table B in S1 Text, respectively. This assumption is based on the empirical interquartile range of mean protein numbers reported in [24] (Fig R in S1 Text). The in-silico distributions of single-cell protein numbers are shown in Fig 7. For each bimodal distribution, let the ‘L peak’ be the peak corresponding to the lower expression rate and the ‘H peak’ be the peak corresponding to the higher expression rate.

Comparing the distributions qualitatively, first, changes in k_d^RNA and k_d^P cause the strongest increase in the distance between peaks and, at the same time, the strongest broadening of the modes. This relates with the influence of these rate constants on the noise in RNA and protein numbers, due to their direct control over the mean in RNA and protein numbers (as the mean numbers decrease, noise, as defined by CV², should increase [67]).

Other parameters allowing large changes in the distance between the two modes are k_bind^H and k_esc^H. This implies that the time that promoters remain occupied, i.e., unavailable for new RNAP bindings is critical in determining the distance between modes. On the other hand, k_bind^L and k_esc^L do not cause similar changes, albeit also controlling the time that promoters remain occupied. This relates to their smaller values relative to all other rate constants of the model (from Table B in S1 Text, k_bind^L and k_esc^L are one order of magnitude smaller than k_bind^H and k_esc^H).

Finally, only a few conditions transform the distribution from bimodal into unimodal. Expectedly from above, only k_bind^H, k_esc^H and k_tr, at their lowest values, as well as k_d^RNA and k_d^P at their highest values achieve this transition. Contrarily, only when increasing both k_L as well as k_H, does one observe a significant overlapping between the two modes of the distribution without loss of bimodality. This is in line with the robustness to perturbations of the bimodality of the genes observed above.

In addition to the simulations above, where we increased or decreased one rate constant at a time, we next considered a more complex transition. Namely, we modeled the events known to occur when cells shift to stationary growth. First, during this transition, the average lifetime of mRNAs increase to 7.8 min [68] while protein degradation rates increase by 8%/h [69]. Meanwhile, dilution due to cell division decreases by approximately 91.5% (as cells stop dividing). We model these changes by reducing decay rates k_d^RNA and k_d^P significantly. Additionally, translation activity also decreases 40 times when cells enter stationary phase [70]. Finally, there are changes in σ³⁸ amounts, but also in the formation rate of RNAPs holoenzymes, due to changes in the levels of anti-sigma factors [71], 6S RNA [72–74], ppGpp/DksA [75], Crl [76], among others (this is approximately accounted for by changing the ratio of available holoenzymes, Section 1.7.3 in S1 Text).

We tuned the bimodal gene expression dynamics assuming gradual changes from the values in Table B in S1 Text to the values in Table C in S1 Text. Fig 7K shows gradual changes in the bimodal distribution of single cell protein numbers, with the mean proteins of the H and L peaks becoming smaller due to the reduction in the numbers of holoenzymes RNAP.σ⁷⁰ and the increased decay rates of RNAs and proteins. Nevertheless, gene expression levels remain bimodal. This suggests that the genes measured whose distributions became unimodal in the stationary growth phase are subjected to even stronger changes than what is assumed in this model.

It is worth noting that, regarding the switching rates between the H and L states, at the moment, there is no empirical information on these rates. One could argue that, by setting these rates slower that all other rates (Table B in S1 Text), we made them unrealistically slow. We based ourselves on the simulations data. First, (i) in Fig 7C, if the transitions were faster than assumed, the two modes of the distribution would overlap (hampering bimodality). Similarly, (ii) the model shows (Fig 7I and 7J) that if the degradation rates become fast relative to the switching frequencies, again the distribution shifts to unimodal. As such, we expect the switching frequencies to be slower than protein degradation rates in E. coli (which are in the order of tens of minutes).

Finally, observing Fig 7, note that in all cases where the bimodal distributions shifted to unimodal by tuning the parameter values, the unimodal distributions are centered at the original ‘low’ state. This occurred because of the boundaries that were set for the parameter values. To show this, we expanded this range for the parameters k_bind^L and k_esc^L (Fig 7D and 7F, respectively). The results, in Fig S in S1 Text, show that, for these wider ranges, the bimodal distribution shifts to unimodal distributions centered around the original ‘high’ state, as long as the other rate-limiting step is increased to match the value in the ‘high’ state.

2.9. Potential sources of bimodality

We first investigated potential common sources of bimodality of the 7 genes. However, since the bimodality is observed also from the genes’ promoters when on plasmids, we excluded the accumulation of positive supercoiling or RNAP collisions due to the activity of neighboring genes as the sources of bimodality. Similarly, we excluded RNA silencing or sequences enhancing transcriptional pausing, since the plasmids only code for the promoter regions of the genes of interest followed by the coding region of the fluorescent protein (i.e., they do not code for the original RNA of the gene of interest or for any membrane-associated protein). Further, we did not find any known TF regulator common to any two genes with bimodal expression levels here reported.

While we did not find a common explanation for all genes, we did find potential, diverse explanations for each but all genes. First, the bimodality could emerge from the regulation by each gene’s specific input TF(s), including a self-regulator TF. This is the case for the paaX gene, who is capable of autorepression [46]. This can lead to bimodality if the repression is nonlinear due to, e.g., feedback delays [77] or if the repression has a switch-like behavior.

While we did not find any other evidence of self-regulation, positive feedback loops are expected to be able of generating bimodality as well. To test, we modeled a positive feedback loop (Fig T in S1 Text) whose parameter values are largely based on the reduced model. From simulations (Fig T in S1 Text), using these values one observes bimodality. Moreover, akin to the reduced bimodal model, altering the transcription rate-limiting steps (RNAP binding and promoter escape rates) changes the distance between the peaks. Meanwhile, changing the ratio of binding and unbinding of the proteins to the promoter region (i.e., the rates of gene activation), changes the height of the peaks (relative to each other), similarly to when changing the rates of shifting between the low and high states in the reduced model.

Another source of bimodality could be TFs with ‘dual effects’. This may occur in the cases of elaB and lldD genes (OxyR and LldR TFs, respectively) [46]. Their TF regulators can switch between being activating or repressing based on environmental or intracellular conditions [78]. Slight differences in these signals among individual cells can drive some cells will have the TF acting as an activator, while others will have the TF acting as a repressor. This can lead to bimodal gene expression within a cell population by enabling the coexistence of two distinct subpopulations: one with high expression and the other with low expression.

Another potential mechanism generating bimodality is multiple σ factor preference. This is the case of pgi and metK genes, as their promoters can be recognize by σ⁷⁰ and by σ³⁸ [46]. This could lead to bimodal distributions if the transcription rates differ with the σ factor [65]. In these genes, cell-to-cell variability in σ factors numbers could potentially cause significant cell-to-cell variability in transcription rates, thus creating bimodal distributions of gene expression.

Finally, manX is regulated by nine TFs, including two activators and two repressors. Such a complex regulatory mechanism could be subject to significant cell-to-cell variability, particularly if small differences in the numbers of the nine transcription factors can have cumulative effects.

Finally, we were unable to identify any potential regulatory mechanism that could explain the bimodal expression observed in the pyrH gene. Given these regulatory differences, the cause of bimodal expression in these genes remains unclear, suggesting that other factors may be contributing to the observed expression patterns.