Which of the following is correct regarding the regulation of lac operon in E. coli?

[M˙]=kM[D]ppc[[Ge]] ℘R[[A]]−[γM+μ]M.

[1]

[E˙]=kE[M]−[γE+μ][E].

[2]

[L˙]=kLβL[[Le]]βG[[ Ge]][Q]−2ϕMM[[L]][B]−μ[L].

[3]

[A]=[L].

[4]

[Q]=[E].

[5]

[B]=[E]/4.

[6]

μ=εJG[[Ge]]+JL[[L]] .

[7]

ppc[[Ge]]=pp1+pc[[Ge ]]kpc−11+pppc[[Ge]] kpc−1.

[8]

pcp[[Ge]]=pc1+pp[[Ge ]]kpc−11+pppc[[Ge]] kpc−1.

[9]

pc[[Ge]]=KGnhKGnh +[Ge]nh.

[10]

℘R[[A]]=1+ξ2ρ1[[A]]1+ξ3ρ1[[A]]+ξ23ρ2[[A]]∏i=1,2,31+ξiρ1[[A]]+∑σ∈P[1,2,3]σ2 L*I, and S[t] = 0 otherwise, where L*I is the value of L at the unstable saddle node located between the uninduced and induced states. Then, the times the system spends in the uninduced and induced states are calculated from the time series S, as well as the histograms of residence times in both states. In all cases, these histograms can be fitted by exponential distributions of the form H[t] = A exp[–kt], where k is the corresponding transition rate constant. In what follows, the transition rate constants from the uninduced [induced] to the induced [uninduced] state shall be denoted as k+ [k−].

Results

Comparison with experimental results

Ozbudak et al. [

] carried out experiments in which E. coli cultures were grown in the M9 minimal medium, with succinate as the main carbon source, supplemented with varying amounts of glucose and TMG. They engineered a DNA segment in which the gfp gene was under the control of the wild-type lac promoter, and inserted this segment into the chromosome of the cultured E. coli bacteria, at the λ-insertion site. In these mutant bacteria, Ozbudak et al. estimated the lac operon expression level in each bacterium by simply measuring the intensity of green fluorescence.

Ozbudak et al. observed that the histograms of fluorescence intensities were unimodal, and that the mean value corresponded to low induction levels of the lac operon, when the bacterial growth medium had low TMG levels. After the TMG concentration surpassed a given threshold, the histograms became bimodal, which can be viewed as evidence for bistability: the original [new] mode corresponds to the uninduced [induced] steady state. With further increments of the TMG concentration, the mode corresponding to the uninduced state disappeared, and the histogram became unimodal again. When the experiment was repeated by decreasing the concentrations of TMG, the opposite behavior was observed. Ozbudak et al. measured the range of TMG concentrations for which bistability was obtained, for several concentrations of external glucose. When they repeated the same experiments with the natural inducer [lactose], they were unable to find analogous evidence for bistability, even when lactose was given at saturation levels. In these last experiments, they employed glucose concentrations in the same range as in the experiments with TMG.

Following Santillán et al. [

], we stated by simulating the Ozbudak et al. experiments. For this, we set μ = 0.02 min−1 to account for the presence of a reliable carbon source [succinate] and ϕM = 0 min−1 [to simulate induction with TMG which is not metabolized by β-galactosidase]. Then, we calculated the bifurcation points and plotted them in the Le versus Ge parameter space, using the procedure described in Methods. We took KA as a free parameter, and found that KA = 9.0 × 105 mpb gives a reasonable fit to the experimental points of Ozbudak et al. Both the model bifurcation diagram and the experimental points are presented in Fig. 3 A. Note that the bistability region predicted by the model is wider than the experimental one. Three possible explanations for this discrepancy are: 1], the lac promoter-gfp fusion employed by Ozbudak et al. as a reporter lacks Operators O2 and O3; 2], the difficulty in measuring exactly the Le values at which the bimodal histograms appear and disappear; and 3], the phase diagram of Fig. 3 A is based upon a mean-field analysis and so, biochemical noise can change the phase boundaries [

].

Figure 3Bifurcation diagrams in the Le versus Ge parameter space, calculated with the present model under various conditions. In all cases, the bistability region is shaded, while the monostable induced [top] and uninduced [bottom] regions are uncolored. In the left column graphs we used the parameter values tabulated in Table 3, and set KA = 8.2 × 105 mpb. In the right column graphs all parameters ξi and ξij were reduced to 0.15 times the value reported in Table 3, and KA was set to 2.8 × 106 mpb. To calculate the graphs in the first row [A and B], we set ϕM = 0 min−1 [to simulate induction of the lac operon with the nonmetabolizable TMG] and μ = 0.02 min−1 [to simulate a bacterial culture growing in a medium that contains a reliable carbon source not affecting the lac operon induction]. The KA values referred to above were chosen to fit the experimental results of Ozbudak et al. , which are shown with solid diamonds. In the second row graphs we keep μ = 0.02 min−1, but increase ϕM to simulate the usage of lactose to induce the operon. In panel C, ϕM 3.8 × 104 min−1, while in panel D, ϕM 3.0 × 104 min−1. Finally, in the third row graphs, a variable growth rate was considered, together with increased ϕM values. The open circles in graphs C–F denote the Ge and Le values for which the stochastic simulations in Figure 4, Figure 5, Figure 6, Figure 7 were carried out.

• View Large Image
• Figure Viewer
• Download [PPT]

A fourth possible explanation for the disagreement between the model and the experimental results is that our estimated parameter values differ from those corresponding to the E. coli strain used by Ozbudak et al. To account for this possibility, we explored the parameter space looking for a better fit. We found that it can be obtained by decreasing the parameters ξi and ξij to 20% of the values reported in Table 3, and by setting KA = 2.7 × 106 mpb. The results are shown in Fig. 3 B.

Effect of using lactose as inducer on the operon bistable behavior

The effect of using a metabolizable inducer [lactose], instead of a nonmetabolizable one [TMG], can be simulated by increasing the value of parameter ϕM[

]. We did this for the two parameter sets described in the previous paragraph, and the results are shown in Fig. 3, C and D. Note that, although quantitatively different, the effect of lactose metabolism on the bistable behavior of the lac operon is qualitatively similar in both cases. Thus, bistability is not eliminated by the inclusion of lactose metabolism, but the region of bistability moves upward toward higher [Le] values. If the value of ϕM is high enough, it can render the operon induction impossible, even with extremely high [Le] values. Moreover, the bistable region collapses into a cusp catastrophe for extremely low [Ge] values. These results agree completely with those in Santillán et al. [

].

Although ϕM is the parameter that changes the most when TMG is substituted by lactose, other parameters like the affinity of the allolactose-repressor complex [KA], the parameters associated to the lactose transport kinetics [kL and κL], and κM, can change as well. We tested the sensitivity of the model bistable behavior to variations on all those parameters, and found that KA is, by far, the most important in that respect.

There is evidence that TMG is a weaker inducer than lactose [

], and so KA should be decreased when induction of the operon with lactose is simulated. After recalculating the bifurcation diagrams of Fig. 3, C and D, with decreased KA values, we found that the upper boundary moves downwards along the Le axis quite a bit, while the lower border also moves downwards but to a lesser extent. Thus, the overall effect is that the width of the bistability region decreases when KA is diminished. For instance, when KA is reduced to 25%, the value used in the TMG simulations and ϕM is set such that the bistability-region lower border lies at ∼Le = 1000 μM, the bistability-region width is reduced by approximately one order of magnitude.

Variable growth rate and bistability

The presence of a reliable carbon source, such as succinate, guarantees the growth of the bacterial culture at a steady rate. Although useful when designing an experiment, this is not a likely situation under natural conditions. When the growth medium contains a mixture of glucose and lactose, and no other carbon source, the bacterial growth rate is not constant but a function of the glucose uptake and the lactose metabolism rates. To investigate the influence of a variable growth rate on the bistable behavior of the lac operon, we recalculated the bifurcation diagrams in the Le versus Ge parameter space, using the model equations in Table 2, with the parameter sets of Fig. 3, C and D. The results are shown in Fig. 3, E and F.

Comparing Fig. 3, C and E, as well as Fig. 3, D and F, we conclude that a variable growth rate has no appreciable effect on the lower boundary of the bistable region. However, the upper border is affected. Instead of being nearly constant for moderate to high values of Ge, and rapidly decreasing for low Ge levels as in Fig. 3, C and D, the upper border increases steadily as Ge increases when μ is not constant. Moreover, for low to moderate Ge levels, the bistable region corresponding to a constant growth rate is wider than that corresponding to a variable μ, while the reverse happens for high external glucose concentrations. Another important difference is that the location of the cusp catastrophe moves to the right when a variable growth rate is considered.

We also tested the effect that decreasing the parameter KA has on the system bistable behavior and found, again, that the width of the bistability region rapidly decreases as KA decreases, while its lower boundary remains with little change.

Intrinsic noise and its effect on the operon dynamics

The lacZ mRNA degradation rate is so high that, according to the model, the average number of mRNA molecules per bacterium is ∼0.75 when the operon is fully induced. Therefore, there must be quite high levels if intrinsic noise, which in turn may have important effects on the system dynamic behavior. To test this, we carried out stochastic simulations [following the process outlined in the Methods section] for several conditions of external lactose and glucose. We repeated these simulations for the four parameter sets corresponding to Fig. 3, C–F.

In all cases, the simulations were carried out for 200,000 min to guarantee that the system reached the steady state and, if the [Ge, Le] point lies inside the bistability region, that the system state switches back and forth between the uninduced and induced states several times. Following Santillán et al. [

], we assumed that ϕM is ∼3–4 × 104 min−1. This implies that even if Le is as high as 1000 μM, the [Ge, Le] point lies very close to the lower border of the bistable region for moderate to high Ge levels.

Histograms of lactose molecule counts were calculated for each stochastic simulation. We observed that whenever the [Ge, Le] point lies below the bistability region, the histograms’ maximum is very close to zero, and they quickly decay as the number of lactose molecules increase. We also carried out simulations for Ge and Le values inside the bistable region, and even a few with the [Ge, Le] point located in the induced monostable region. The histograms we obtained are presented in Figure 4, Figure 5, Figure 6, Figure 7, which respectively correspond to the parameter sets of Fig. 3, C–F.

Figure 4Histograms of lactose molecules per bacterium, calculated from stochastic simulations carried out with the set of parameter values corresponding to Fig. 3 C [constant growth rate], for different combinations of Ge and Le values. All simulations were run for 200,000 min to ensure that, if the [Ge, Le] point lies in the bistability region, the system gates several times between the uninduced and the induced states. In each histogram, the transition rates from the uninduced to the induced [k+], and from the induced to the uninduced [k–] states are shown.

• View Large Image
• Figure Viewer
• Download [PPT]

Figure 5Histograms of lactose molecules per bacterium, calculated from stochastic simulations carried out with the set of parameter values corresponding to Fig. 3 D [constant growth rate], for different combinations of Ge and Le values. In each histogram, the transition rates from the uninduced to the induced [k+] and from the induced to the uninduced [k–] states are shown.

• View Large Image
• Figure Viewer
• Download [PPT]

Figure 6Histograms of lactose molecules per bacterium, calculated from stochastic simulations carried out with the set of parameter values corresponding to Fig. 3 E [constant growth rate], for different combinations of Ge and Le values. In each histogram, if bistability exists, the transition rates from the uninduced to the induced [k+], and from the induced to the uninduced [k–] states are shown.

• View Large Image
• Figure Viewer
• Download [PPT]

Figure 7Histograms of lactose molecules per bacterium, calculated from stochastic simulations carried out with the set of parameter values corresponding to Fig. 3 F [constant growth rate], for different combinations of Ge and Le values. In each histogram, if bistability exists, the transition rates from the uninduced to the induced [k+], and from the induced to the uninduced [k–] states are shown.

• View Large Image
• Figure Viewer
• Download [PPT]

One would expect the histograms to be bimodal whenever the [Ge, Le] point is located within the bistability region; each mode corresponding to one of the two available steady states. This is clearly observed in some of the plots in Figure 4, Figure 5, Figure 6, Figure 7, but not in as many as expected. The reason for this is that both modes are remarkably wide, and thus merge such that they are hard to distinguish separately. In general, the bimodal distributions are more prominent as the value of Le increases and the [Ge, Le] point moves upwards in the bistability region, and as the value of Ge decreases and the [Ge, Le] point approaches the cusp catastrophe.

The effect of a variable growth rate can be analyzed by comparing Figure 4, Figure 6, as well as Figure 5, Figure 7. Notice that the changes in the behavior of the constant and variable growth-rate histograms are consistent with the cusp catastrophe point moving right along the Ge axis, in the variable μ case. Interestingly, the fact that the upper border of the bistable region behaves differently in the constant and variable μ cases seems to have no effect on the histograms, except for the points that lie in the induced monostable region in the later case.

We also calculated, following the procedure detailed in the Methods section, the transition rates from the uninduced to the induced state [k+], and from the induced to the uninduced one [k–]. We did this for each of the stochastic simulations described in the previous paragraph, and the results are shown next to the corresponding histograms in Figure 4, Figure 5, Figure 6, Figure 7. Notice that, in all cases, k+ increases and k− decreases as Le increases; and that k– increases and k+ decreases as Ge increases. This behavior indicates that the induced [uninduced] state strengthens [weakens] its stability as the [Ge, Le] point approaches the bistability region upper boundary, and as it gets away from its lower border.

Other interesting observations that can be derived from the reported transition rates are that the residence times predicted by the model are all of ∼103 min, and that the bimodal distributions can only be clearly noticed if k+⪆2k−.

Effect of the media conditions on the operon induction level and the bacterial growth rate

To better understand the dynamic performance of the lac operon, we numerically solved the model equations with the parameter set of Fig. 3 E, in which the growth rate is variable for several different levels of Ge and Le. In all cases, the numeric simulations were carried out for 2000 min to ensure that the system reaches a steady state, starting with the initial values M = 0 mpb, E = 0 mpb, and L = 0 mpb. Given this initial condition, the system reaches the uninduced steady state whenever it exists; otherwise, it evolves toward the induced steady state. For every one of these simulations, we recorded the β-galactosidase steady-state value, which we use as an indication of the operon induction level, and computed the corresponding steady-state growth rate. The results are shown in Fig. 8.

Figure 8Steady-state operon induction level [A], defined as proportional to the concentration of LacZ polypeptides [E], and steady-state bacterial growth rate [B], calculated for different values of Ge and Le. For given values of Ge and Le, the E and μ steady-state values were computed by numerically solving the model equations for 2000 min, to ensure that the system reaches a steady state. This starts with the initial condition M = 0 mpb, E = 0 mpb, and L = 0 mpb, so the system reaches the uninduced steady state whenever it exists. Black and white, respectively, represent the minimum and maximum values of the corresponding variables in both graphs.

• View Large Image
• Figure Viewer
• Download [PPT]

In Fig. 8 A, note that, for Le⪆500 μM, the system remains in the uninduced state whenever Ge⪆50μM [the exact value depends on Le]. However, there is a sudden transition to the induced state when the external glucose level decreases below a certain value. From the way we carried out the simulations, this transition corresponds to the saddle-node bifurcation in which the uninduced state disappears by colliding with the intermediate saddle node. Once the system is in the induced state, the operon induction level increases continuously until it reaches its maximum value for very low Ge values. However, when 200μM⪅Le⪅500μ M, there is a smooth transition between the uninduced and the induced states as Ge decreases, in agreement with the bifurcation diagram of Fig. 3 E.

The operon can also be induced by increasing Le, given a fixed Ge value. We see in Fig. 8 A that the lac operon induction is smooth for 30μM⪅Ge, that the transition to the induced state is abrupt for 30μM⪅Ge⪅50μM, and that the operon never gets induced if Le ≤ 1000 μM and 50μM⪅Ge. These results are again in agreement with the bifurcation diagram of Fig. 3 E. In particular, the border between the smooth and abrupt transitions to the steady state corresponds to the location of the cusp catastrophe in the Le versus Ge parameter space.

Remarkably, the bacterial growth rate reaches close to maximal values for almost all combinations of Ge and Le values; see Fig. 8 B. The exceptions are the region just right to the border of the monostable induced region [where the value of μ is intermediate], the region of very low Ge and medium to high Le levels [where the μ-value is intermediate as well], and the region of very low Ge and medium to low Le concentrations [in which the growth rate is very close to zero].

Conclusions

We have improved a previously published model for the lac operon in Escherichia coli. The improvements include: the new model takes into account the interaction between the repressor and its three different binding sites [operators] in a more detailed way; it accounts for the cooperativity between the activator [CAP] and DNA folding; and it considers the fact that the bacterial growth rate depends on the glucose uptake and the lactose metabolism rates.

Except for KA, all the model parameters were estimated from reported experimental data. KA was taken as a free parameter to fit the model results to the experimental data of Ozbudak et al. [

], who reported the [Ge, Le] bifurcation points corresponding to the induction of E. coli's lac operon with TMG. We can see in Fig. 3 A that the bistability region predicted by the model is wider than the experimental one. In our opinion, there are four possible reasons for this discrepancy between the data and the model predictions:

• The lac promoter-gfp fusion employed by Ozbudak et al. as a reporter lacks Operators O2 and O3.

• Experimentally, it is quite difficult to measure with precision the Le values at which the bistability region starts and ends.

• The phase diagram of Fig. 3 A is based upon a mean-field analysis and so, biochemical noise can change the phase boundaries [

  ].

• The estimated parameter values differ from those corresponding to the E. coli strain used by Ozbudak et al.

We explored the parameter space to account for the last possibility, and found that by decreasing parameters ξi and ξij to 20% of their estimated value a much better fit can be obtained [compare to Fig. 3 B]. In all the following analysis, we used both parameter sets to test the robustness of the obtained conclusions.

Ozbudak et al. found no evidence of bistability when the natural inducer, lactose, was employed. From their simulations of the lac operon evolution on in silico bacterial populations, van Hoek and Hogeweg [

] proposed that the operon parameters evolve in such a way that bistability is not available for the bacteria when they grow on lactose. Santillán et al. [

] tackled this issue from a mathematical modeling viewpoint and concluded, in opposition to van Hoek and Hogeweg [

], that bistability does not disappear because of lactose metabolism, although it is very hard to identify with the Ozbudak et al. experimental setup.

Since the present model is an improvement and a generalization of the model in Santillán et al. [

], we started by corroborating their results. In Fig. 3, C and D, we can see that with a constant growth rate the current model results are in complete agreement with those of Santillán et al. [

] for the two considered parameter sets. That is, the bistability region moves upwards in the Le versus Ge parameter space, and it disappears in the region of very low Ge values via a cusp catastrophe. Following Santillán et al. [

], we choose the value of ϕM [which accounts for the maximum lactose-to-allolactose metabolism rate] such that even if Le = 1000 μM, the [Ge, Le] point is very close to the lower border of the bistability region.

Given that a constant growth rate is not realistic under natural conditions, we considered the influence of a variable μ on the bistable behavior of the lac operon. It can be seen in Fig. 3, E and F, that having a μ that depends on the external conditions has no appreciable effect on the lower border of the bistability region. However, the upper border is greatly modified; instead of being nearly constant from moderate to high values of Ge and rapidly decreasing for low Ge levels, as in Fig. 3, C and D, the upper border increases steadily as Ge increases. Furthermore, the cusp catastrophe location moves rightwards, along the Ge axis, when a variable μ is considered.

We also carried out stochastic simulations, considering both a constant and a variable growth rate, to test how feasible the Ozbudak et al. experimental setup is to identify bistability when lactose is the inducer. According to Ozbudak et al., whenever a bimodal distribution is found, it can be considered a proof for bistability. More precisely, having a bimodal distribution is a sufficient but not a necessary condition for this behavior: for instance, if we have large standard deviations, a bimodal distribution will hardly be noticed even is bistability is present. Our numeric experiments confirm that, for the employed ϕM values, identifying bistability via expression-level distributions is very difficult, regardless of whether a constant or a variable growth rate is considered; compare to Figure 4, Figure 5, Figure 6, Figure 7. It is known that the shot noise resulting from the production of multiple protein copies per mRNA transcript also lowers the switch stability [

] and probably widens the distribution, thus further blurring the bimodal distribution. This could perhaps explain why, in the experiments by Ozbudak et al., only a unimodal distribution was observed in the presence of lactose.

van Hoek and Hogeweg [

] simulated, in silico, the evolution of the lac operon in bacterial populations. According to their results, the parameters that determine the expression of the lac operon genes evolve in such a way that, when the bacterial population is grown in a medium rich in lactose, bistability disappears. Contrarily, when the bacteria grow in the presence of a gratuitous inducer, van Hoek and Hogeweg predict that bistability would be present. They propose that this behavior explains why Ozbudak et al. [

] did not observe bistability when E. coli was grown on lactose. From the results discussed above, we propose a different explanation; namely, that bistability does not disappear when lactose is the employed inducer, but the bistability region in the Le versus Ge parameter space is modified in such a way that this behavior is masked due to the intrinsic noise.

E. coli and other bacteria can feed on both lactose and glucose. However, when they grow in a medium rich in both sugars, glucose is utilized before lactose starts being consumed. This phenomenon, known as diauxic growth [

,

], represents an optimal thermodynamic solution given that glucose is a cheaper energy source than lactose: to take advantage of lactose, the bacteria needs to expend energy in producing the enzymes needed to transport and metabolize this sugar [

].

To test whether the results discussed above are compatible with diauxic growth and, more generally, to better understand the lac operon performance under variable conditions of external lactose and glucose, we calculated the plots in Fig. 8. To do that, the model equations were numerically solved for several values of Ge and Le, with the parameter set of Fig. 3 E, and with the initial condition M = 0 mpb, E = 0 mpb, and L = 0 mpb. This initial condition was chosen so the system reaches the uninduced steady state, unless it does not exist. In this last case, the system evolves toward the induced state.

Consider a bacterial culture growing in a medium containing a mixture of glucose and lactose. We can see from the results in Fig. 8 A [where the operon induction level is plotted versus Ge and Le] that diauxic growth is predicted by the model. That is, induction of the lac operon only occurs when glucose is exhausted. However, it can take place in two different ways: if Le < 500 μM, the cells undergo a smooth transition from the uninduced to the induced state; otherwise, the lac operon shifts abruptly from the off- to the on-state when glucose is almost completely depleted. Why this difference? To answer this question, we carried out numeric experiments in which the lactose operon is induced by changing the bacterial culture from a medium rich in glucose [in which it has been growing for a long time] to a medium with no glucose, while Le is kept constant. Our results [not shown] indicate that it takes >1000 min for the lac operon to become fully induced when Le = 300 μM; while, when Le = 1000 μM, the operon achieves a 97% induction level 300 min after the medium shift. We conclude from these last results that, with a very high external lactose concentration, the lac operon can remain fully uninduced until glucose utilization is almost complete because the response time is short in these conditions.

In principle, it should be possible to induce the lac operon by increasing the Le value. However, the results in Fig. 8 A predict that, if Ge > 60 μM, this is practically impossible for any physiological concentration of external lactose. This can be interpreted as an optimal behavior; if we consider that glucose is a cheaper carbon source, it makes sense to keep the genes in the lactose operon uninduced, as long the energy necessary to grow can be obtained from that sugar. We can also observe from Fig. 8 A that, when 30 μM < Ge < 60 μM, the operon undergoes an abrupt transition to the induced state with high enough levels of Le. However, the induction level only reaches an intermediate value. Finally, if Ge < 30 μM, the operon induction level increases smoothly as Le increases. In this last case, the induction level can get very close to its maximum possible value; a higher value is obtained sooner as Ge decreases. This behavior can also be interpreted as optimal; if the amount of glucose in the growing medium is not enough to support the bacterial population, it is reasonable to allow a continuous increase of the lac-operon induction level as the amount of lactose in the environment increases.

In Fig. 8 B, the bacterial growth rate [calculated under the same conditions as Fig. 8 A] was plotted versus Ge and Le. Remarkably, the value of μ is very close to its maximum almost everywhere; the only exception is the region corresponding to very low levels of Ge and Le. We conclude from this that bistability not only helps to guarantee an efficient switching of the lac operon in E. coli, when feeding on glucose, lactose, or both sugars, but also plays an important role in maintaining a high bacterial growth rate for an extremely wide range of environmental conditions.

The above discussion regarding the influence of lactose metabolism on the bistable behavior of the lac operon relies on the assumption that only the value of parameter ϕM is affected when the operon is induced with lactose instead of TMG. We know, being that TMG is a weaker inducer, KA should be decreased as well. We investigated the sensitivity of the system bistable behavior to modifications on this parameter and found that the width of the bistability region—in the Le versus Ge parameter-space bifurcation diagram—rapidly decreases as KA diminishes, while its lower boundary remains with little change.

These facts oblige us to reconsider our conclusions. Thus, if the change on KA is small and the bistability region is still quite wide, then all the above stated conclusions hold like that. Otherwise, our results regarding the system behavior when the [Ge, Le] point lies within the bistability region are probably no longer valid. However, not everything is lost because, even then, bistability does not disappear, and the incapability of Ozbudak et al. to observe it can still be explained if the bistability region is so high along the Le axis [due to a very high ϕM value], that the lac operon never gets induced, except for extremely low Ge values [below the cusp catastrophe].

Appendix A:. Model development

Interaction of the lac repressor with the operators

The lac repressor is a tetramer, formed by the union of two functional homo-dimers. Each monomer in a repressor can be bound by an allolactose molecule. When all four monomers are free, each dimer can bind any of the lac-operon operators. When one or the two monomers in a dimer are bound by allolactose molecules, it is unable to bind an operator; but the dimer in the opposite side still can do it, provided both their monomers are free. Finally, if either of the two monomers is on an opposite side, or three or all four monomers are bound to allolactose, none of the repressor dimers is capable of binding an operator.

Given that there are six possible configurations in which two repressor monomers are bound by allolactose molecules, and that in two of such configurations both bound monomers belong to the same dimer, we can assume that one-third of the repressors with two occupied monomers are capable of binding an operator. This assumption is supported by the fact that the binding of allolactose to a given monomer is independent of its binding to the other three monomers.

Let R, RA, and R2A, respectively, represent a free repressor, and a repressor bound by one and two allolactose molecules. They can bind an operator, O, according to the following reactions:

R+O ⇌KR/2OR,RA+O⇌KRORA,R2A+O⇌KROR 2A.

The dissociation constant of the first reaction [KR/2] is half as large as those of the second and third reactions [KR] because, in a free repressor, both dimers can bind the operator, and therefore the forward reaction rate is twice as large.

Recalling that only one-third of the R2A compounds can bind an operator, the equilibrium equations corresponding to the reactions in the former paragraph are

[R][O]=KR2[ OR],[RA][O]=KR[ORA],1 3[R2A][O]=KR[OR2A].

Assume that the total concentration of operators [[OT]] remains constant:

[O]+[OR]+[ORA]+[O R2A]=[OT].

Then, it follows after a little algebra that the concentration of free operators is given by

[O]= 11+2[R]KR+[RA]KR+[R2A ]3KR[OT].

[20]

The binding of allolactose [A] to the lac repressor [R] takes place in four consecutive independent steps as

R+A⇌KA/4RA,RA+A⇌2KA/3R2A,R2A+A⇌ 3KA/2R3A,R3A+A⇌4KAR4A.

The different factors multiplying the dissociation constants account for the fact that the forward and backward reaction rates are proportional to the number of empty and occupied allolactose binding sites, respectively.

The equilibrium equations for the repressor-allolactose binding reactions are

[R][A]=14KARA,[RA][ A]=23KA[R2A],[R2A][A]=32 KA[R3A],[R3A][A]=4KA[R4A ].

From these and the equation accounting for the conservation of repressor molecules,

[R]+[RA]+[R2A ]+[R3A]+[R4A]=[RT],

it follows that [R], [RA], and [R2A] are, respectively, given by

[R]=1[1+[A]/KA]4[ RT],

[21]

[RA]=4[A]/KA[1+[A]/KA]4[RT],

[22]

[R2A] =6[[A]/KA]2[1+[A]/KA]4[ RT].

[23]

Finally, by substituting Eqs. [21], [22], [23] into Eq. 20, we obtain the following expression for the concentration of free operators, in terms of the allolactose concentration,

with

[R′]=1[1+[A]/KA]2[RT].

[25]

DNA folding

Assume that there are two different operators [or lac repressor binding sites] that only allolactose-free repressor molecules can bind, and that DNA can fold in such a way that a single repressor can bind both operators simultaneously. Let O, OR1,OR2,O2R, and OR [12], respectively, denote the binding states in which both operators are free; the first operator is bound by a repressor and the second is free; only the second operator is bound by a repressor; each operator is bound by a different repressor molecule; and both operators are bound by the same repressor through DNA folding. The chemical reactions standing for all of the transitions between these states are

O+R⇌KR1/2OR1,O+R⇌KR2/2OR2,OR1+R ⇌KR2/2O2R,OR2+R⇌K R1/2O2R,OR1⇌K2 OR[12],OR2⇌K1OR[12].

The equilibrium equations for these equations are

[O][R]=K R12[OR1],[O][R]=KR22[O R2],[OR1][R]=KR22[O2R ],[OR2][R]=KR12[O2R], [OR1]=K2[OR[12]],[ OR2]=K1[OR[12]].

The condition that the total concentration of DNA regulatory regions remains constant is

[O]+[OR1]+[OR2]+[O2R]+[ OR[12]]=[OT].

It follows from the last two equations that the probability for both operators to be free is

[O][OT]=11+2[R]KR1+2[R]KR2+2[R]KR12[R]KR2+ [R]KR12,

[26]

where KR12=KR1K2/ 2=KR2K1/2. Similarly, the probability that the first operator is free, regardless of the second operator state, is given by

[O]+[OR 2][OT]=1+2[R]KR21+2[R]KR1+2[R]KR2+2[R]KR1 2[R]KR2+[R]KR12.

[27]

To take into account the dependence of these probabilities on the allolactose concentration, we only need to substitute R from Eq. 21.

A more realistic model should take into account the fact that all of the repressors bound by one allolactose, and one-third of those bound by two allolactose molecules, can also bind an operator, although they are unable to simultaneously bind the two operators through DNA folding. To carry out the corresponding analysis, more chemical reactions and their corresponding equilibrium equations must be included, as in the previous subsection. Although not too complicated, the whole process is quite lengthy. Therefore, we omit the details for the sake of brevity and simply give the results. For instance, the probability that both operators are free is given by

[O][OT]=11+2[R′]KR1+2[R′ ]KR2+2[R′]KR12[R′]KR 2+[R]KR12,

[28]

where [R] and [R′] are, respectively, given by Eqs. [21], [25]. Similarly, the probability the first operator is free, regardless of the second operator state, is

[O]+[OR2][OT]=1+2[R′]KR21+2[R′] KR1+2[R′]KR2+2[R′]KR1 2[R′]KR2+[R]KR12.

[29]

It is not hard to see how these results are a direct generalization of the ones presented in Eqs. [24], [26], [27].

The three operators of the lac operon

The lactose operon has three different operators, termed O1, O2, and O3. Free lac repressor molecules, as well as those bound by one allolactose, and one-third of those bound by two allolactose molecules, can bind each of these operators with different dissociation constants. Moreover, free repressor molecules can bind any two of such operators simultaneously. Denote the binding state in which all three operators are free by O; the state in which only the ith operator is bound by a repressor by ORi; that in which operators i and j are bound by one repressor each as O2Rij ; the state in which all three operators are bound by different repressor molecules as O3R; the state in which operators i and j are bound by a single repressor as OR[ij]; and that in which operators i and j are bound by a single repressor and operator k is bound by a second repressor as O2R[ij]j.

By following a similar but lengthier process to that outlined in the previous subsections, one can calculate the probabilities associated to all the binding states of this system. Given the very large number of equations involved, and of their concomitant equilibrium equations, we omit the details and give the results directly.

The probability that all three operators are free is

[O][OT ]=1∏i=1,2,31+2[R]′KRi+ ∑σ∈P[1,2,3]σ2