Article Information

• PDF (441 kb)
• Full Text with Large Figures
• Cited by in Scopus (7)

PubMed

• Articles by M. Santillán
• Articles by E.S. Zeron

Related Articles

• Bistable Behavior in a Model of the lac Operon in Escherichi...

  Bistable Behavior in a Model of the lac Operon in Escherichia coli with Variable Growth Rate Biophysical Journal, Volume 94, Issue 6, 15 March 2008, Pages 2065-2081 M. Santillán Abstract This work is a continuation from another study previously published in this journal. Both the former and the present works are dedicated to investigating the bistable behavior of the operon in from a mathematical modeling point of view. In the previous article, we developed a detailed mathematical model that accounts for all of the known regulatory mechanisms in this system, and studied the effect of inducing the operon with lactose instead of an artificial inducer. In this article, the model is improved to account, in a more detailed way, for the interaction of the repressor molecules with the three operators. A recently discovered cooperative interaction between the CAP molecule (an activator of the lactose operon) and Operator 3 (which influences DNA folding) is also included in this new version of the model. The growth rate dependence on the rate of energy entering the bacteria (in the form of transported glucose molecules and of metabolized lactose molecules) is also considered. A large number of numerical experiments is carried out with this improved model. The results are discussed in regard to the bistable behavior of the lactose operon. Special attention is paid to the effect that a variable growth rate has on the system dynamics. Abstract | Full Text | PDF (550 kb)

• Multistability: a major means of differentiation and evoluti...

  Multistability: a major means of differentiation and evolution in biological systems Trends in Biochemical Sciences, Volume 24, Issue 11, 1 November 1999, Pages 418-422 Michel Laurent and Nicolas Kellershohn Abstract Very simple biochemical systems regulated at the level of gene expression or protein function are capable of complex dynamic behaviour. Among the various patterns of regulation associated with non-linear kinetics, multistability, which corresponds to a true switch between alternate steady states, allows a graded signal to be turned into a discontinuous evolution of the system along several possible distinct pathways, which can be either reversible or irreversible. Multistability plays a significant role in some of the basic processes of life. It might account for maintenance of phenotypic differences in the absence of genetic or environmental differences, as has been demonstrated experimentally for the regulation of the lactose operon in . Cell differentiation might also be explained as multistability. Abstract | Full Text | PDF (137 kb)

• Mathematical Modeling of Gene Networks

  Mathematical Modeling of Gene Networks Neuron, Volume 26, Issue 3, 1 June 2000, Pages 567-580 Paul Smolen, Douglas A Baxter and John H Byrne Full Text | PDF (380 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 11, 3830-3842, 1 June 2007

doi:10.1529/biophysj.106.101717

Biophysical Theory and Modeling

Origin of Bistability in the lac Operon

M. Santillán^*, †, ,  , M.C. Mackey^†, ‡ and E.S. Zeron^§

^* Unidad Monterrey, Centro de Investigación y Estudios Avanzados del Instituto Politécnico Nacional, Monterrey, México; Escuela Superior de Física y Matemáticas, Instituto Politécnico Nacional, México DF, México
^† Centre for Nonlinear Dynamics in Physiology and Medicine, McGill University, Montreal, Canada
^‡ Departments of Physiology, Physics, and Mathematics, McGill University, Montreal, Canada
^§ Departamento de Matemáticas, Centro de Investigación y Estudios Avanzados del Instituto Politécnico Nacional, México DF, México

Address reprint requests to M. Santillán.

At the molecular level, biological systems function using two types of information: genes, which encode the molecular machines that execute the functions of life, and networks of regulatory interactions, specifying how genes are expressed. Substantial progress over the past few decades in biochemistry, molecular biology, and cell physiology has ushered in a new era of regulatory interaction research. Recent analysis has revealed that cell signals do not necessarily propagate linearly. Instead, cellular signaling networks can be used to regulate multiple functions in a context-dependent fashion. Because of the magnitude and complexity of the interactions in the cell, it is often not possible to understand intuitively the systems behavior of these networks. Rather, it has become necessary to develop mathematical models and analyze the behavior of these models, both to develop a systems-level understanding and to obtain experimentally testable predictions. This, together with the fact that DNA micro-arrays, sequencers, and other technologies have begun to generate vast amounts of quantitative biological data, has accelerated the shift away from a purely descriptive biology and toward a predictive one.

Recent computer simulations of partial or whole genetic networks have demonstrated collective behaviors (commonly called systems, or emergent, properties) that were not apparent from examination of only a few isolated interactions alone. Among the various patterns of complex behavior associated with nonlinear kinetics, multistability is noteworthy. Multistability corresponds to a true switch between alternate and coexisting steady states, and so allows a graded signal to be turned into a discontinuous evolution of the system along several different possible pathways. Multistability has certain unique properties not shared by other mechanisms of integrative control. These properties may play an essential role in the dynamics of living cells and organisms. Moreover, multistability has been invoked to explain catastrophic events in ecology 1, mitogen-activated protein kinase cascades in animal cells 2,3,4, cell cycle regulatory circuits in Xenopus and Saccharomyces cerevisiae5,6, the generation of switchlike biochemical responses 2,3,7, and the establishment of cell cycle oscillations and mutually exclusive cell cycle phases 6,8, among other biological phenomena. On the other hand, there are also serious doubts that multistability is the dynamic origin of some biological switches 9. Nevertheless, it is generally accepted that two paradigmatic gene-regulatory networks in bacteria, the phage λ switch and the lac operon (at least when induced by nonmetabolizable inducers), do show bistability 10,11,12. In the former, bistability arises through a mutually inhibitory double-negative-feedback loop, while in the latter, a positive-feedback loop is responsible for the bistability.

The lac operon, the phage-λ switch, and the trp operon, are three of the best known and most widely studied systems in molecular biology. The inducible lac operon in E. coli is the classic example of bistability. It was first noted by Monod and co-workers more than 50 years ago, although it was not fully recognized at the time. The bistable behavior of the lac operon has been the subject of a number of studies. It was first examined in detail by Novick and Weiner 13 and Cohn and Horibata 14. Later experimental studies include those of Maloney and Rotman 15 and Chung and Stephanopoulos 10. Over the last few years, the lac operon dynamics have been analyzed mathematically by Wong et al. 16, Vilar et al. 17, Yildirim and Mackey 18, Santillán and Mackey 19, Yildirim et al. 20, Tian and Burrage 21, and Hoek and Hogeweg 22. This was possible because this system has been experimentally studied for nearly 50 years and there is a wealth of biochemical and molecular information and data on which to draw. Recently, Ozbudak et al. 23 performed a set of ingenious experiments that not only confirm bistability in the lac operon when induced with the nonmetabolizable inducer thiomethylgalactoside (TMG), but also provide new and novel quantitative data that raise questions that may be answered via a modeling approach.

Previous studies have analyzed bistability and its dynamic properties as a systems phenomenon. However, to our knowledge, none of them have dealt with questions like: “How did multistability arise within the context of gene regulatory networks?” or “What evolutionary advantages do bistable regulatory networks have when compared with monostable ones?” In this article, we address these issues from a mathematical modeling approach, basing our examination on the dynamics of the lac operon.

A mathematical model for the lac operon is developed in Appendix A Model development,Appendix B Binding-state probability distribution for a molecule with multiple, independent, binding sites,Appendix C Cooperativity between binding sites. All of the model equations are tabulated in Table 1 and they are briefly explained below. The reader may find it useful to refer to Fig. 1, where the lac operon regulatory pathway is schematically represented.

Table 1 Full set of equations for the model of lactose operon regulatory pathway depicted in Fig. 1

		Eq. No.
	(1)
	(2)
	(3)
	(4)
	(5)
	(6)
	(7)
	(8)
	(9)
	(10)
	(11)
	(12)
	(13)

Differential Eqs. (1) govern the time evolution of the intracellular concentration of mRNA (M), polypeptide (E), and lactose (L) molecules. Ge and Le, respectively, stand for the extracellular glucose and lactose concentrations. A, Q, and B represent the intracellular allolactose, permease, and β-galactosidase molecule concentrations. The functions βL, and βG, respectively, account for the negative effect of external glucose on the initiation rate of transcription (via catabolite repression), the probability that the lactose promoter is not repressed, the positive effect of external lactose on its uptake rate, and the negative effect of external glucose on lactose uptake (inducer exclusion). The expression is the rate of lactose metabolism per β-galactosidase. Finally, pc and ρ represent internal auxiliary variables.

Display large version of this figure

Figure 1

(A) Schematic representation of the lac operon regulatory mechanisms. This operon comprises genes lacZ, lacY, and lacA. Protein LacZ is a permease that transports external lactose into the cell. Protein LacY polymerases into a homo-tetramer named β-galactosidase. This enzyme transforms internal lactose (Lac) into allolactose (Allo) and galactose (Gal), and further transforms allolactose into galactose. Allolactose can bind to the repressor (R) inhibiting it. When not bound by allolactose, R can bind to a specific site upstream of the operon structural genes and thus avoid transcription initiation. External glucose inhibits production of cAMP which, when bound to protein CRP to form complex CAP, acts as an activator of the lac operon. External glucose also inhibits lactose uptake by permease proteins. (B) Graphical representation of the interactions accounted for by the lac operon mathematical model. The meaning of the variables appearing in this figure is as follows. Ge and Le stand for external glucose and lactose concentrations; M, Q, and B denote mRNA, permease, and β-galactosidase, molecule concentrations, respectively; and L, A, and G correspond to the intracellular lactose, allolactose, and galactose molecule concentrations. All of the processes underlying the lac operon regulatory pathway are represented by rectangles. Inputs (outputs) are denoted with empty arrowheads (circles). Finally, plus and minus signs stand for the effect each input variable has on every output variable.

The model consists of three differential equations (Eqs. (1)), that respectively account for the temporal evolution of mRNA (M), lacZ polypeptide (E), and internal lactose (L) concentrations.

Messenger RNA (mRNA) is produced via transcription of the lac operon genes, and its concentration decreases because of active degradation and dilution due to cell growth. The value γ_M represents the degradation plus dilution rate, k_M is the maximum transcription rate per promoter, and D is the average number of promoter copies per bacterium. The function (defined in Eq. (9)) accounts for regulation of transcription initiation by active repressors. The fraction of active repressors is proportional to ρ(A) (compare Eq. (10)); since repressors are inactivated by allolactose (A), ρ is a decreasing function of A. Furthermore, the rate of transcription initiation decreases as the concentration of active repressors increases. Concomitantly, is a decreasing function of ρ, and thus it is an increasing function of A. The function (Eq. (9)) denotes the modulation of transcription initiation by external glucose, i.e., through catabolite repression. Production of cyclic AMP (cAMP) is inhibited by extracellular glucose. cAMP further binds the so-called cAMP receptor protein (CRP) to form the CAP complex. Finally, CAP binds a specific site near the lac promoter and enhances transcription initiation. The probability of finding a CAP molecule bound to its corresponding site is given by p_c (Eq. (8)) and is a decreasing function of Ge. Moreover, is an increasing function p_c and, therefore, a decreasing function of Ge.

The translation initiation rate of lacZ transcripts is k_E, while γ_E is the dilution and degradation rate of lacZ polypeptides. Let B and Q, respectively, denote the β-galactosidase and permease concentrations. Given that the corresponding parameters for lacY transcripts and polypeptides attain similar values, that β-galactosidase is a homo-tetramer made up of four identical lacZ polypeptides, and that permease is a lacY monomer, it follows that Q=E and B=E/4 (Eqs. (5)).

Lactose is transported into the bacterium by a catalytic process in which permease protein plays a central role. Thus, the lactose influx rate is assumed to be k_Lβ_L(Le)Q, with the function β_L(Le) given by Eq. (11). Extracellular glucose negatively affects lactose uptake (so-called inducer exclusion), and this is accounted for by β_G(Ge) (Eq. (12), which is a decreasing function of Ge. Once inside the cell, lactose is metabolized by β-galactosidase. Approximately half of the lactose molecules are transformed into allolactose, while the rest enter the catalytic pathway that produces galactose. In our model, with defined in Eq. (13), denotes the lactose-to-allolactose metabolism rate, which equals the lactose-to-galactose metabolism rate. Allolactose is further metabolized into galactose by β-galactosidase. From the assumptions that the corresponding metabolism parameters are similar to those of lactose, and that the allolactose production rate is much higher than its degradation plus dilution rate, it follows that A ≈ L (Eq. (4)).

All of the parameters in the model are estimated in Appendix D . Their values are tabulated in Table 2.

We carried out stochastic simulations with the above-described model. This was done by means of Gillespie's Tau-Leap algorithm 24,25, which we implemented in Python (http://www.python.org/download/).

All of the analytical results we used to study the model dynamic behavior are explained in detail in the Appendices.

A mathematical model for the lac operon regulatory pathway was developed as explained in Theory, above. The model equations are tabulated in Table 1. These equations determine the time evolution of variables M, E, and L, which respectively stand for mRNA, lacZ polypeptide, and internal lactose concentrations. Special attention was paid to the estimation of the model parameters from reported experimental data. The model parameters are tabulated in Table 2. Below we describe the results obtained from this model.

The model steady states and their stability are analyzed in Appendix E . As seen there, depending on the values of Le and Ge, there can be up to three steady states. For very low Le values, there is a single stable steady state corresponding to the uninduced state. As Le increases, two more steady states appear via a saddle-node bifurcation. One of these new fixed points is stable, corresponding to the induced state, while the other is a saddle node. With further increases in Le the saddle node and the original stable steady-state approach until they eventually collide and are annihilated via another saddle-node bifurcation. Afterwards, only the induced stable steady state survives.

Ozbudak et al. 23 carried out a series of clever experiments regarding the bistable behavior of the lac operon in E. coli. In these experiments, the Le values at which the bifurcations take place were experimentally determined for various extracellular glucose concentrations, using a nonmetabolizable inducer (TMG). In the present model, the usage of a nonmetabolizable inducer can be modeled by setting ϕ_M=0; recall that this parameter stands for the maximum rate of lactose-to-allolactose and of lactose-to-galactose metabolism.

Using the procedure described in Appendix E , we numerically calculated the model Le bifurcation values for several Ge concentrations when a nonmetabolizable inducer is employed. The results are shown in the bifurcation diagram of Figure 2A, and compared with the experimental results of Ozbudak et al. 23. The agreement between the experimental data and the model predictions is sufficiently close to give us confidence in the other results we report here.

Display large version of this figure

Figure 2

(A) Bifurcation diagram in the Le versus Ge parameter space, calculated with ϕ_M=0. As discussed in the text, this choice of the parameter ϕ_M simulates induction of the lac operon with a nonmetabolizable inducer, like TMG. The dots correspond to the experimental results of Ozbudak et al. (2004). Panels B–D show bifurcation diagrams calculated with the following lactose metabolism rates: ϕ_M=3.6×10³min⁻¹, ϕ_M=1.5×10⁴min⁻¹, and ϕ_M=4.0×10⁴min⁻¹, respectively.

Ozbudak et al. further assert that, when their experiments were repeated with the natural inducer (lactose), they were unable to identify any bistable behavior. This result, and other studies reporting similar conclusions 22, have sparked a lively debate about whether bistability is a property of the lac operon under naturally occurring conditions, or it is only an artifact introduced by the use of nonmetabolizable inducers.

To investigate the effects of lactose metabolism, we recalculated the bifurcation diagram of Figure 2A with different values of the parameter ϕ_M. Bifurcation diagrams in the Le versus Ge parameter space, calculated with ϕ_M=0, 3.6×10³, 1.5×10⁴, and 4.0×10⁴min⁻¹, are shown in Figure 2BD. From these results we conclude that the main effects of lactose metabolism on the Le versus Ge bifurcation diagrams are

In their experiments, Ozbudak et al. 23 let a bacterial culture grow for a long time (more than four hours) in a medium with constant glucose and TMG concentrations. Then they measured the expression level of gene lacZ in every bacterium and plotted the corresponding histogram. When the Ge and TMG concentrations were set such that the induced and uninduced stable steady states were both available, a bimodal distribution for the lacZ expression level was observed. Otherwise, they obtained unimodal distributions.

Ozbudak et al. later repeated the same experiments using lactose as an inducer, and observed only unimodal distributions, even for saturating values of Le. Thus, they did not find any evidence of bistability with lactose as an inducer. According to the bifurcation diagrams in Fig. 2, these findings can be explained if In that case, even when Le is as high as 1000μM, the system should be in the monostable uninduced region, or close to it, and thus any sign of bistability will be hard to detect. To test this hypothesis we carried out stochastic simulations by means of Gillespie's Tau-Leap algorithm 24,25, which we implemented in Python. In these simulations we set the values of Ge and Le, and let the system evolve for 20,000min, recording the value of all variables every minute. We performed these numerical experiments with Le=1000μM, and Ge=4, 10, 20, 40, 100, 200, 400, and 1000μM.

In Fig. 3 we show the histograms of the internal lactose molecule count calculated from these numerical experiments. Notice that:

Display large version of this figure

Figure 3

Normalized histograms for the number of lactose molecules per average-size bacterium. This histograms were calculated from the results of stochastic simulations in which we let the system evolve, for 10,000min, with Le=1000μM and Ge=1000, 400, 200, 100, 40, 20, 10, and 4μM; see the main text for details. Notice that the lac operon remains in the uninduced state for Ge≥100μM, and then it jumps to the induced state when Ge≤20μM. Of all the histograms, only that corresponding to Ge=40μM shows a vaguely bimodal distribution characteristic of bistability. These results are in complete accord with the bifurcation diagrams presented in Fig. 2.

Following Santillán and Mackey 19, we analyze the effect of inducer exclusion on the bistable behavior of the lac operon. For this, we construct—by setting β_G(Ge)=1 in Eq. (3)an in silico mutant strain of E. coli in which this regulatory mechanism is absent, and recalculate the bifurcation diagram in the Le versus Ge parameter space. The result is shown in Fig. 4. As seen there, the absence of this mechanism moves the bifurcation region downwards, so the monostable induced region starts at Le values similar to those obtained when TMG is used as inducer (Figure 2A). If this mutant strain can be engineered, it may be possible to identify bistability with experiments like those of Ozbudak et al. 23.

Display large version of this figure

Figure 4

Bifurcation diagram, in the Le versus Ge parameter space, for an in silico mutant strain of E. coli in which the inducer-exclusion regulatory mechanism is absent. This bifurcation diagram was calculated with ϕ_M=4.0×10⁴min⁻¹.

We have developed a mathematical model of the lac operon that, though not accurate in all detail, captures the essential system behavior with respect to its bistable characteristics. The model reproduced the experimental results of Ozbudak et al. 23, who report the bifurcation diagram for the lac operon in E. coli, when induced with a nonmetabolizable inducer.

As mentioned above, there is currently a debate about whether bistability is a property of the lac operon under naturally occurring conditions. Ozbudak et al. reported that bistability was not observed when the natural inducer, lactose, was employed. Later, van Hoek and Hogeweg 22 simulated the in silico evolution of a large set of lac operons (under conditions of fluctuating external lactose and glucose), and concluded that bistability is present for artificial inducers, but not for lactose. Moreover, Narang and Pilyugin 26 argued, from a modeling study, that E. coli can exhibit diauxic growth when growing in a mixture of sugars, even if bistability is absent. To gain further insight into this matter, we analyzed the influence of lactose metabolism on the system bistable behavior, by varying the parameter ϕ_M (the maximum lactose-to-allolactose metabolism rate). Below, we discuss how our results provide a possible explanation to such apparent absence of bistability, and suggest that bistability does not disappear because of lactose metabolism, although it is highly modified.

In terms of external lactose, the bistability region moves upwards in the Le versus Ge parameter space, as ϕ_M increases. If ϕ_M>2.5×10⁴min⁻¹, bistability is not present in the limit of very low Ge, and it only appears after Ge surpasses a given threshold, via a cusp catastrophe. We suspect that this last situation is the most likely in wild-type bacteria, because it becomes practically impossible for the operon to switch into the induced state when there is an abundance of glucose. Moreover, there is no resistance, due to bistability, to activate the lac genes when glucose is absent, and so to employ the energy provided by lactose.

The fact that Ozbudak et al. 23 did not observe signs of bistability when lactose is used as an inducer suggests that or higher. In that case, even when Le is as high as 1000μM, the (Ge, Le) point lies most of the time in the monostable uninduced region of the bifurcation diagram (see Figure 2D), or it is very close to it, and thus bistability should be very hard to detect. Our stochastic simulations support this assertion given that, of all the experiments we carried out, only that corresponding to Ge=40μM and Le=1000μM rendered a bimodal like distribution.

According to Ozbudak et al., whenever a bimodal distribution is found, it is considered a proof of 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. Then, our numeric experiments confirm that—if —identifying bistability via expression-level distributions is very difficult. We have also shown that if a mutant strain of E. coli in which inducer exclusion is absent or deficient can be found, it may be possible to identify bistability with the methods of Ozbudak et al., even if the natural, induced lactose is employed.

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 27,28, represents an optimal thermodynamic solution given that glucose is a cheaper energy source than lactose since, to take advantage of lactose, the bacteria needs to expend energy in producing the enzymes needed to transport and metabolize this sugar 29.

Consider a bacterial culture growing in a medium containing a mixture of glucose and lactose. According to the bifurcation diagram in Figure 2D, induction of the lac operon, as glucose is exhausted, can take place in two different ways: if Le<600μ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.

To better understand the behaviors described in the previous paragraph, 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 the lactose concentration is kept constant. Such simulations were performed by numerically solving the model differential equations with the Runge-Kutta method, implemented in Octave's algorithm, lsode. Our results (not shown) indicate that it takes >1000min for the lac operon to become fully induced when Le=300μM, whereas when Le=1000μM, the operon achieves a 97% induction level 300min after the medium shift.

From these results we conclude that, with a very high external lactose concentration, the lac operon can remain fully uninduced until glucose utilization is almost completely given up, because the response time is short in these conditions. That is, the most efficient performance consists of an abrupt change from the off- to the on-state, driven by a small variation on the external glucose concentration, if Le is very high. An analysis of Figure 2D reveals that, if Le=1000μM, the system goes from the uninduced to the induced monostable regions when Ge decreases from 210 to 25μM.

We have argued that bistability helps to guarantee an efficient performance of the lac operon in E. coli, when feeding on glucose, lactose, or both sugars. Could a regulatory pathway involving only monostability be equally efficient? The lac operon in E. coli is an inducible operon. This means that it is subject to positive feedback regulation through lactose (allolactose). In our model, this is accounted for by the fact that in Eq. (1) is an increasing function of A. Therefore, the higher the intracellular level of lactose, the higher the transcription initiation rate of the lactose operon genes. Depending on the functional form of we can have either bistability for some parameter values, or a unique single stable steady state for all the parameters. In general, needs to be highly sigmoidal to have bistability.

Expression of the lac operon is modulated by external glucose through the function On the other hand, by analyzing the model steady state we can see that as Ge decreases from 210 to 25μM, the system jumps from the uninduced to the induced states, and the L* steady-state value increases by a factor of 56. Substituting all these values into Eqs. (7), we conclude that an incremental increase of from 0.14 to 0.26 drives the 56-fold increase of L*. That is, the following amplification relation is observed:

It is well known, since the invention of the regenerative circuit in the early decades of the 20th Century, that very large amplifications can destabilize a monostable system subject to positive feedback 30,31,32, and that the only way to get large amplifications is to approach as much as possible to the stability limit. In Appendix F we analyze the possibility of having large amplifications, with monostable regulation, in the present model of the lac operon. There we show that the stability of the steady state is highly compromised whenever the amplification exponent is >4. Since we predict an amplification exponent of 6.5, larger than four, the stability of our system would be at risk if it were controlled by a monostable pathway. In Appendix F we also demonstrate that bistability (multistability) allows both large amplification and strongly stable steady states. Hence, since a large amplification of the lac operon expression level (driven by a small change in ) is advantageous for E. coli (when Le is very high), we conclude that bistability ensures an efficient consumption of lactose and glucose without jeopardizing the system stability.

Appendix A

Model development

In this section, a mathematical model of the lac operon in E. coli is developed. The reader may find it convenient to refer to Fig. 1 of the main text, where the lac operon regulatory mechanisms are schematically represented. The model presented here accounts for the time evolution of the following variables: intracellular mRNA (M), lacZ polypeptide (E), and intracellular lactose (L) concentrations. The dynamics of these variables are modeled by the (balance) differential equations:

(14)

(15)

(16)

The meaning of the functions and parameters in these equations is as follows:

A is the intracellular allolactose concentration.
D is the concentration of lac promoter within the bacterium.
Q and B are the permease and β-galactosidase concentrations, respectively.
k_M, k_E, and k_L are, respectively, the lac promoter transcription initiation rate, the lacZ mRNA translation initiation rate, and the maximum lactose uptake rate per permease.
γ_M, γ_E, and γ_L, respectively, stand for the dilution plus degradation rates of M, E, and L.
Ge and Le denote the external glucose and lactose concentrations.
is the probability that promoter P1 is not repressed, while takes into account the effect of external glucose on the probability of having an mRNA polymerase bound to this promoter (catabolite repression).
β_L(Le) and β_G(G_e) are functions that account for the modulation of lactose uptake as functions of the external lactose and glucose (inducer exclusion) concentrations, respectively.
represents the rate of lactose metabolization per β-galactosidase molecule.

Probability of having a polymerase bound to the lac promoter, including catabolite repression

Let P and C, respectively, denote the mRNA polymerase and CAP concentrations. A polymerase can bind by itself to the lac promoter, P₁. However, the affinity of this reaction is increased when a CAP molecule is bound to its corresponding binding site in the DNA chain. By taking this cooperative behavior into account, as well as the results in Appendix B Binding-state probability distribution for a molecule with multiple, independent, binding sites,Appendix C Cooperativity between binding sites, the probability P_D can be calculated as

(17)

where K_P and K_C are the respective dissociation constants of the polymerase-promoter and CAP-DNA complexes, and k_pc>1 is a constant that accounts for the cooperative interaction between the promoter and the CAP binding site. Equation (17) can be rewritten as

(18)

In this equation, p_p and p_c(G_e), respectively, denote the probabilities that a polymerase is bound to the promoter in the absence of CAP, and that a CAP molecule is bound to its binding site in the absence of mRNA polymerase. The probability p_p is constant and given by

However,

is a function of the external glucose concentration. Experimentally p_c(G_e) must be a decreasing function of G_e since external glucose inhibits the synthesis of cAMP, which in turns implies a decrease in C and thus in p_c(G_e). Here we assume that the functional form for p_c(G_e) is given by

(19)

Probability that the promoter P₁ is not repressed

The lac operon has three different operator regions denoted by O₁, O₂, and O₃ all of which can bind active repressor and are involved in transcriptional regulation. A repressor bound to O₁ avoids transcription initiation. Conversely, a repressor bound to either O₂ or O₃ does not seem to affect transcription. Nevertheless, DNA can fold in such a way that a single repressor simultaneously binds two operators in all possible combinations. These complexes are more stable than that of a repressor bound to a single operator, and all of them inhibit transcription initiation 33. From the results in Appendix B , the probability that the lac operon promoter is not repressed can then be calculated as

(20)

where R is the concentration of active repressors, K_i (i=1, 2, 3) is the dissociation constant of the R-O_i complex, and K_ij (i, j=1, 2, 3, i<j) is the dissociation constant of the O_i-R-O_j complex.

Define ρ(A)=R/K₁, ξ_i=K₁/K_i (i=2, 3), and With these definitions, Eq. (20) can be rewritten as

(21)

Repressor molecules are tetramers and each of their subunits can be bound by an allolactose molecule, inactivating the repressor. The concentration of active repressors as a function of the allolactose concentration is given by Eq. (19),

where R_T is the total repressor concentration, and K_A is the allolactose-repressor subunit complex dissociation rate. Since ρ(A)=R/K₁, it follows that

(22)

where ρ_max=R_T/K₁.

Inducer uptake rate

From the results of Ozbudak et al. 23, in the absence of external glucose the dependence of the normalized inducer uptake rate, per β-permease molecule, on the external inducer concentration is given by

(23)

The results of Ozbudak et al. 23 also allow us to model the decrease in the inducer uptake rate caused by external glucose by

(24)

Lactose metabolism

After being transported into the bacterium, lactose is metabolized by β-galactosidase. Approximately half of the lactose molecules are transformed into allolactose, while the rest are directly used to produce galactose. The lactose-to-allolactose metabolism rate, which as explained above equals the lactose-to-galactose metabolism rate, can be modeled as a Michaelis-Menten function, typical of catalytic reactions

Therefore, the total rate of lactose metabolized per single β-galactosidase is

(25)

Allolactose metabolism

Allolactose and lactose are both metabolized into galactose by the same enzyme β-galactosidase. Under the assumption that the allolactose production and metabolism rates are much faster than the cell growth rate, we may assume that these two metabolic processes balance each other almost instantaneously so

On the other hand, allolactose is an isomer of lactose, and therefore we can expect that the parameters related to the metabolism kinetics of both sugars attain similar values: ϕ_M ≈ ϕ_A and κ_M ≈ κ_A. Then

(26)

Permease and β-galactosidase concentrations

The translation initiation and degradation rates of lacZ and lacY are slightly different. Here, we assume they are the same for the sake of simplicity. From this and the fact that β-galactosidase is a homotetramer, made up of four lacZ polypeptides, while permease consists of a single lacY polypeptide, it follows that

(27)

Appendix B

Binding-state probability distribution for a molecule with multiple, independent, binding sites

Consider a molecule D with specific, independent, binding sites for N different molecules B_i, i=1…N. Denote the binding state of molecule D by (n₁, n₂, …n_N)=({n_i}), where n_i=1 if a molecule B_i its bound to its corresponding site in molecule D, and n_i=0 otherwise. In what follows, we demonstrate by induction that the probability of any given binding state is

where [B_i] is the concentration of chemical species B_i, and K_i is the dissociation constant of the D:B_i complex.

Case N=1

The reaction leading to the formation of complex D:B₁ is

At equilibrium, the concentrations of the chemical species involved in this reaction satisfy the relation

(28)

Assume a constant total concentration for species D, i.e.,

(29)

Then, from Eqs. (28), the fractions of free and bound B₁ sites are respectively given by

From these results, the binding-state probability distribution for the current system can be written as

(30)

where n₁=0 if the B₁ binding site is empty, and n₁=1 otherwise.

Case N=2

From the N=1 case, the probability distribution for the B₁ binding site in the absence of species B₂ is

Also, the probability distribution for the B₂ binding site in the absence of species B₁ is

If both sites are independent, their joint binding-state probability distribution is

(31)

General case

Take a D molecule with N independent binding sites and assume that the binding-state probability distribution for the first N–1 sites is given, in the absence of species B_N, by

In the absence of species B₁, …B_N–1, the probability of the B_N binding site is

Since all sites are independent from each other, the binding-state probability distribution can be calculated as

(32)

which proves our original assertion.

Appendix C

Cooperativity between binding sites

We can see from Eq. (32) that the probability of finding n_i=0, 1 molecules B_i bound to its corresponding binding site is proportional to . The denominator in the fraction of Eq. (32) plays the role of a normalizing constant.

Assume that two given sites, say l and m, have a cooperative interaction in the sense that the probability of finding the two of them bound by their respective molecules is larger than the product of the individual probabilities. From this and the considerations of Appendix B , the binding-state probability distribution accounting for cooperativity between sites l and m is

(33)

where k_c>1 measures the strength of the cooperativity.

Appendix D

Parameter estimation

Growth rate, μ

The growth rate of a bacterial culture depends strongly on the growth medium. Typically, the mass doubling time varies from 20 to more than 40min 34. For the purpose of this study, we take a doubling time of 30min, which corresponds to the growth rate μ ≈ 0.02min⁻¹.

Lac promoter concentration

According to Bremmer and Dennis 34, there are ∼2.5 genome equivalents per average E. coli cell at the growth rate determined by μ. For the purpose of this work, we take

Production rates

Malan et al. 35 measured the transcription initiation rate at the lac promoter and report k_M ≈ 0.18min⁻¹.

From Kennell and Riezman 36, translation of the lacZ mRNA starts every 3.2s. According to Beckwith 37, the production rate of lac permease is smaller than that of β-galactosidase monomers even though, as Kennell and Riezman 36 report, there are similar levels of both mRNA species. This suggests that lacY mRNA values are translated at a lower rate than is lacZ mRNA. Nevertheless, to our knowledge, there are no reported measurements of the lacY mRNA translation initiation rate. Thus, we assume it is equal to that of lacZ:

According to Chung and Stephanopoulos 10, the inducer uptake rate per lac permease is

Degradation rates

Kennell and Riezman 36 measured a lacZ mRNA half-life of 1.5min. That is, its degradation rate is 0.46min⁻¹, and the corresponding dilution plus degradation rate is

According to Kennell and Riezman 36, the lac permease degradation rate is 0.01min⁻¹. Thus, its degradation plus dilution rate is

Here we assume that the lactose degradation rate is negligible and so its degradation plus dilution rate is simply

Catabolite repression parameters

Malan et al. 35 measured the polymerase-lac promoter affinity in the presence and absence of cAMP. From their data, we estimate

The parameters p_p, K_G, and n_h were estimated by fitting Eq. (5) of the main text (with p_c as given by Eq. (6)) to the experimental data reported in Ozbudak et al. 23. See Fig. 5. The values we obtained are

Display large version of this figure

Figure 5

Effect of external glucose, due to catabolite repression, on the activity of the lac promoter. The dots stand for the experimental data of Ozbudak et al. 23, while the solid line was drawn from Eq. (5), with p_c as given by Eq. (6), with the estimated parameter values.

Repression parameters

Oehler et al. 33 studied how the three operators in the lac operon cooperate in repression. According to their results, when only O₁ and either O₂ and O₃ are present, the repression level is reduced to 53.85% and 33.85% that of the wild-type operon, respectively. Moreover, when O₂ and O₃ are destroyed, repression is reduced to 1.38% that of the wild-type operon. This allowed us to estimate the parameters ξ₂, ξ₃, and ξ₁₂₃ as and

For this, we took into consideration that eliminating O₂ (O₃) is equivalent to making 1/K₂=1/K₁₂=1/K₂₃=0 (1/K₂=K₁₃=1/K₂₃=0) in Eq. (7) of the main text. Parameters ξ₂, ξ₃, and ξ₁₂₃ are modified accordingly in Eq. (8) of the main text.

The parameters ρ_max and K_A were estimated by fitting the curves in the model bifurcation diagram to the experimental results of Ozbudak et al. 23 (see Figure 2A of the main text) obtaining the values

Inducer uptake parameters

From the experimental data of Ozbudak et al. 23, the inducer uptake rate per active permease, as a function of external inducer concentration, can be fitted by a Michaelis-Menten function with a half-saturation concentration of 680μM. That is,

Ozbudak et al. 23 measured the inducer-uptake-rate decrease per active permease as a function of external glucose concentration. We found (plot not shown) that these data are well fit by Eq. (11) of the main text with

Lactose metabolism parameters

Lactose metabolism rate and saturation constant, ϕ_M and κ_M. We estimate these parameters from the data reported in Martínew-Bilbao et al. 38 as ϕ_M ≈ 3.60×10³min⁻¹, and κ_M≈7.0×10⁵mpb, where mpb stands for molecules per average bacterium.

Appendix E

Steady-state stability and bifurcation diagrams

Stability analysis

The model equations (Eqs. (1)) in the main text can be rewritten as

(34)

(35)

(36)

The steady-state values M*, E*, and L* are, respectively, given by

(37)

(38)

(39)

It is easy to show numerically that given the parameter values tabulated in Table 2 and Ge>0, Eq. (39) may have up to three roots, and therefore that the model has up to three steady states.

Notice that, except for and all terms in the right-hand side of Eqs. (34) are linear with respect to M−M*, E−E*, and L−L*. Since we are only interested in analyzing the system dynamic behavior in a small neighborhood around the steady states, the following linear approximation is employed:

With this approximation, the full model can be reduced (around the steady state) to the following system of linear differential equations,

(40)

where

and

The general solution of the linear system (Eq. (40)) is where C_k are undetermined constants, while α_k and u_k are the respective eigenvalues and eigenvectors of matrix A. Therefore, if all eigenvalues α_k have negative real parts, the corresponding steady state is stable. Conversely, the steady state is unstable if at least one eigenvalue has a positive real part.

The eigenvalues of A are the roots of its characteristic polynomial, det(sI – A). It follows after some algebra that det(sI – A)=P(s)−Q₃, with

Notice that γ₃>0 because γ_L>0 and M is an increasing function; moreover γ_M>γ_E>0. Hence, the polynomial P(s) is positive and strictly increasing for every real s≥0. Assume now that Q₃>P(0), then, P(s) – Q₃ has one positive real root since P(0) – Q₃<0 and P(s) is strictly increasing for s>0.

On the other hand, we assert that all roots of P(s) – Q₃ have negative real part whenever Q₃<P(0). It follows from the work of Strelitz 39 that all roots have negative real part if and only if the following pair of inequalities hold simultaneously:

and we obviously have that

Therefore, the linear system is stable (unstable) whenever Q₃<P(0) (Q₃>P(0)). Furthermore, since the steady-state values of L* are the roots of the equation with

(41)

we can easily calculate that a given steady state is locally stable (unstable) whenever the derivative

(42)

is strictly negative (positive).

Calculation of the bifurcation points and the bifurcation diagram

Given Ge, the values of L* and Le at which a saddle node bifurcation occurs can be calculated by simultaneously solving and with as defined in Eq. (41). The bifurcation diagram in the Le versus Ge parameter space can be calculated by repeating the above procedure for several values of Ge.

Appendix F

Positive feedback, stability and hysteresis

Consider the following simplified model of a biological switch with nonlinear feedback,

(43)

where γ>0 is a degradation parameter, β>0 is an external control parameter, and the smooth function P(x)>0 is nonlinear with respect to x. The steady states x* of Eq. (43) are the solutions to

(44)

It is easy to prove that the dynamic system given by Eq. (43) is locally stable at x* if and only if

(45)

We choose system Eq. (43) because its dynamic behavior is representative of (and similar to) the behavior of higher dimensional systems. For example, consider the following two-dimensional system, with positive parameters β_k>0 and γ_k>0:

The steady states x* and y* are the solutions of

These equations are equivalent to Eq. (44) with β=β₁β₂ and γ=γ₁γ₂. Furthermore, the two dimensional system is locally stable at x* and y* if and only if both eigenvalues of the matrix

have strictly negative real parts. Simple algebraic calculations allows us to prove that both eigenvalues of A₂ have strictly negative real part if and only if

Notice that this last result is equivalent to Eq. (45), after setting the degradation parameter γ=γ₁γ₂ and the external control parameter β=β₁β₂. Similar results holds for higher dimensions.

Consider the model equations (Eqs. (1)) with the external glucose Ge values in the interval [25, 210] μM, the external lactose Le=10³μM and ϕ_M=4×10⁴min⁻¹. It is easy to show, under the previous conditions, that the steady state L* achieves the minimum and maximum values of 4.5×10⁵ and 2.53×10⁷molecules per average bacterium, respectively. Therefore,

The above results indicate that the term β_L(L_e)β_G(G_e) inside the square brackets of Eq.