Bistable Behavior of the Lac Operon in E. Coli When Induced with a Mixture of Lactose and TMG

In this work we investigate multistability in the lac operon of Escherichia coli when it is induced by a mixture of lactose and the non-metabolizable thiomethyl galactoside (TMG). In accordance with previously published experimental results and computer simulations, our simulations predict that: (1) when the system is induced by TMG, the system shows a discernible bistable behavior while, (2) when the system is induced by lactose, bistability does not disappear but excessively high concentrations of lactose would be required to observe it. Finally, our simulation results predict that when a mixture of lactose and TMG is used, the bistability region in the extracellular glucose concentration vs. extracellular lactose concentration parameter space changes in such a way that the model predictions regarding bistability could be tested experimentally. These experiments could help to solve a recent controversy regarding the existence of bistability in the lac operon under natural conditions.

a mutually inhibitory double-negative-feedback loop, and it is responsible for the decision between the lytic and the lysogenic pathways made by bacteria after being infected by the phage λ virus. In the lac operon, a positive-feedback loop is responsible for bistability, and this phenomenon has been adduced to explain how E. coli feeds on lactose but only when glucose (its favorite carbon source) is absent from the growing medium.
Recently, Ozbudak et al. (2004) quantified the lac operon bistable behavior by carrying out a series of clever experiments with genetically modified E. coli bacteria. Their results confirmed the existence of bistability in this system, when it is induced by the non-metabolizable inducer thiomethyl galactoside (TMG), but failed to demonstrate bistability when the operon was induced with lactose.
The results of Ozbudak et al. (2004) have open a controversy regarding the existence of bistability in the lac operon of E. coli under natural conditions. On the one hand, van Hoek and Hogeweg (2006) carried out in silico simulation of the lac operon evolution in bacterial populations. In these simulations, the parameters that control the expression of the lac operon genes evolve in such a way that the system is not bistable when the bacterial population grows in a lactose rich medium, but it does exhibit bistability when artificial inducers are employed. Thus, van Hoek and Hogeweg (2006) argue from their computational experiments that the wildtype lac operon is not a bistable switch under natural conditions. On the other hand, Santillán and colleagues (Santillán et al., 2007;Santillán, 2008) developed a detailed mathematical model for the lac operon regulatory pathway (estimating the model parameters from reported experimental data) and examined the effect of lactose metabolism. According to their results, lactose metabolism

IntroductIon
Mathematical models and computer simulations of genetic networks have demonstrated collective behaviors (commonly called systems, or emergent properties) that were not apparent from examination of the network components. Among the various patterns of complex behavior associated with non-linear kinetics, multistability is noteworthy. A system is said to be multistable when it can rest in two or more stable steady states for a given parameter set. By smoothly changing the parameter values one can make the stable steady state count discontinuously change, via saddle-node bifurcations. This allows a graded signal to be turned into a discontinuous evolution of the system (usually accompanied by hysteresis) along different discrete states.
The two paradigmatic gene-regulatory networks concerning bistability are the phage λ switch and the lac operon of Escherichia coli (Chung and Stephanopoulos, 1996;Ferrell, 2002;Angeli et al., 2004). As a matter of fact, bistability in biochemical networks was first discovered in the lac operon of E. coli, when induced with non-metabolizable inducers (Novick and Weiner, 1957;Cohn and Horibata, 1959). In the phage λ switch, bistability arises through greatly modifies the system bistable region, in the external lactose vs. external glucose parameter space, but does not make this behavior disappear.
The current understanding of the lac operon alludes to bistability to explain how E. coli consumes glucose and lactose in the most efficient possible way (Lehninger et al., 2005). Namely, the genes necessary to uptake and metabolize lactose are only turned on when this sugar is present (at high levels) in the growing environment, while glucose (the bacterium favorite carbon source) is absent. If, as some studies suggest, bistability is lacking in the lac operon under natural conditions, this would open the question of how the regulatory pathway in this system allows an efficient consumption of alternative carbon sources.
Interestingly, the gene networks responsible for the alternate consumption of glucose and tryptophan in E. coli (Konan and Yanofsky, 1997), and of glucose and galactose in S. cerevisiae (Acar et al., 2005;Bhat and Iyer, 2009), have the same architecture as the lac operon of E. coli. That is, the regulatory pathway controlling the expression of the genes involved in the transport and metabolism of carbon sources other than glucose seems to be conserved over different metabolic networks and over different species. We see from these facts that the issue of whether the lac operon is bistable or not in its natural habitat transcends this system, and the clarification of this question may help to understand the performance of a whole class of gene networks. Moreover, experiments similar to the ones here proposed may help to study bistability in other gene networks that allow the cell to transport and metabolize carbon sources other than glucose.
As mentioned above, the modeling studies of Santillán and colleagues (Santillán et al., 2007;Santillán, 2008) predict that the lac operon shows bistability when induced with lactose. Unfortunately, as asserted by the same authors, is quite hard to experimentally test this prediction. With this in mind, the purpose of the present work is to extend the Santillán et al. (2007) model (by considering the operon induction with a mixture of lactose and TMG) in order to propose feasible experiments whose results will validate or disprove it, as well as its predictions regarding bistability. Similar experiments may help to elucidate whether bistability exists or not in other gene networks with the same architecture. This work is in the line of Acar et al. (2005), who suggested methods for experimentally studying bistability in S. cerevisiae.

MaterIals and Methods the lactose operon regulatory pathway
The genes in the lac operon code for the proteins necessary to transport and metabolize lactose. These genes, which span a DNA segment around 6000 bp long, are: lacZ, lacY, and lacA. The protein coded by lacZ forms an homotetrametic enzyme called β-galactosidase. This enzyme metabolizes lactose into either allolactose (approximately half of it) or into glucose and galactose (the other half). The same enzyme is also responsible for metabolizing allolactose into glucose and galactose. On the other hand, lacY codes for a permease involved in the uptake of lactose from the environment, while lacA codes for thiogalactoside transacetylase, which plays no role in the operon regulatory pathway. Lactose permease is also capable of transporting gratuitous inducers (like TMG or IPTG), but they cannot be metabolized by β-galactosidase.
The elements involved in the lac operon regulation are located as shown in Figure 1A. There, P represent the promoter, C is the binding site for the catabolite activator protein (which enhances transcription initiation), and O1, O2, and O3 stand for active repressor protein binding sites (Oehler et al., 1990).

Figure 1 | (A)
Schematic representation of the regulatory elements located in lac operon DNA. P denotes the promoter, O1, O2, and O3 correspond to the three operators (repressor binding sites), and C is the binding site for the cAMP-CRP complex. The different ways in which a repressor molecule can interact with the operator sites are represented in (B-e). Namely: a free repressor molecule (B), one with a single subunit bound by an inducer molecule (D), or one with the two subunits in the same side bound by allolactose (e) can bind a single operator. Moreover, a free repressor molecule can bind two different operators simultaneously (C). Figure adapted from Santillán (2008). increase its affinity for mRNA polymerase. On the other hand, the presence of glucose in the medium prevents the enhancement of transcription initiation because cAMP is not synthesized, and thus CRP proteins are not activated. This regulatory mechanism is known as catabolite repression.
Another regulatory mechanism operates through transcriptional repression and involves a repressor molecule. The lac repressor is a homo-tetramer (consisting of two functional homo-dimers) of LacI polypeptides (Lewis, 2005;Wilson et al., 2007). Each functional dimer can bind Operators O1, O2, and O3. Whenever a repressor is bound to operator O1, it avoids transcription initiation by blocking the mRNA polymerase. On the other hand, DNA can fold in such a way that a single repressor binds two operators simultaneously, one per dimer. The effect of these loops is to enhance the stability of the corresponding complex.
Each monomer in the lac repressor can be bound by an allolactose or a TMG molecule, inhibiting the capability of the corresponding dimer to bind an operator. This means that free repressors can bind one operator ( Figure 1B), or two of them simultaneously ( Figure 1C); that repressors with three free monomers can bind one but not two operators ( Figure 1D); that repressors with two free monomers can bind one operator, if the bound monomers belong to the same dimer ( Figure 1E), or none at all; and that repressors with only one free monomer, as well as repressors with all four monomers bound by an inducer molecule (Narang, 2007), are unable to bind any operator. Therefore, active repressors are inhibited, and so repression is relieved, by increasing the intracellular inducer concentration.
The last regulatory mechanism in the lac operon is the so-called inducer exclusion. In it, external glucose decreases the efficiency of lac permease to transport lactose (Reznikoff, 1992), and by doing so negatively affects the induction of the operon genes.

Model developMent
In Section "Model Development" in Appendix we develop a mathematical model for the lac operon regulatory pathway when it is induced with a mixture of lactose and the artificial inducer TMG. This model is based on a previous model introduced by Santillán (2008), which is modified to take into account the presence of TMG. Like the original model, the present one accounts for all of the regulatory mechanisms described in the previous subsection.
The model here introduced consists of four ordinary differential equations that account for the dynamics of mRNA concentration  Table 1 and are briefly described in the forthcoming paragraphs. Equation (1) accounts for the mRNA concentration dynamics. There, k M stands for the transcription initiation rate of the active lacP promoter; p pc ([G e ]) is a function of the extracellular glucose concentration ([G e ]) -defined in (9) -is the probability that a mRNA polymerase is bound to promoter lacP, taking into account cooperativity with a CAP protein bound to its corresponding site; function P R ([A],[T]) -defined in (12)is the probability that the lac promoter is not repressed and so that a polymerase bound to it can initiate transcription; [A] is defined in (5) and represents the intracellular allolactose concentration; parameter γ M represents the mRNA degradation rate; and μ Glucose is the favorite carbon and energy source for E. coli, as well as for many other organisms. Although this bacterium can also feed on other sugars (like lactose), it only does so when glucose is absent. Thus, if a bacterial culture grows in a medium containing a mixture of glucose and another sugar, it will exclusively feed on the former until it is exhausted, and only then will it switch to the second one. A consequence of this behavior is that the bacterial growth curve shows two distinctive phases. This phenomenon was originally studied by Monod (1941), who described it as diauxic growth. As seen below, the regulatory mechanisms in this operon are responsible for diauxic growth in E. coli. The regulatory mechanisms responsible for this behavior are reviewed below. We refer the reader to the schemes in Figure 2 for a better understanding.
When glucose is absent from the bacterial growing medium, cyclic AMP (cAMP) is synthesized at high rates. In turn, cAMP binds CRP proteins turning them active. Finally, active CRP proteins (also known as catabolite activator proteins, CAP) bind an specific site right upstream of promoter lacP and considerably Figure 2 | Schemes representing the lac operon regulatory pathways when the system is induced with lactose (A), and when it is induced with an artificial inducer like TMG (B). Labeled rectangles represent chemical species, arrows with empty heads denote processes through which one chemical species is transformed into another, arrows with solid heads indicate interactions that enhance the process they point to, and finally, lines ending in perpendicular bars denote interactions that diminish the process they point to.

Equation
(2) stands for the LacE polypeptide concentration dynamics. The meaning of parameters in it is as follows: k E is the translation initiation rate at the ribosome binding site of gene lacE and γ E is the degradation rate of lacE polypeptides.
Equation (3) governs the intracellular lactose concentration dynamics. There, k L is the maximal lactose uptake rate per permease molecule; [Q] -defined in (6) (19) -corresponds to the activity of β-galactosidase enzymes, as a function of intracellular lactose concentration.
Analogously, Equation (4) governs the dynamics of the intracellular TMG concentration. Parameter k T in this equation corresponds to the maximal TMG uptake rate per permease, while β([L e ],[T e ]) is the permease activity for TMG transport as a function of [L e ] and [T e ]. Unlike lactose, TMG is not metabolized by β-galactosidase and so it is known as gratuitous inducer. That is, TMG allows the functioning of a part of the regulatory pathway, but disables the functioning of other one (see Figures 2A,B). Thus, a comparison of the system performance under different lactose and TMG concentration would allow us to gain insight into the performance of the whole system.
As mentioned before, Equation (9) defines function p pc ([G e ]), which stands for the probability that a mRNA polymerase is bound to promoter lacP, taking into account cooperativity with a CAP protein bound to its corresponding binding site. Similarly, Equation (10) defines function p cp ([G e ]), which denotes the probability that a CAP protein is bound to its site, taking into consideration cooperativity with a polymerase bound to promoter lacP. Both, Equations (9) and (10) invoke function p c ([G e ]) -defined in (11)that accounts the probability that a CAP proteins binds its binding site non-cooperatively as a function of [G e ].
The function accounting for the repression regulatory mechanism (12)  In this work we have paid special attention to the estimation of all the model parameters from reported experimental data. The estimation procedure is described in detail in Section "Parameter Estimation" in Appendix and the resulting parameter values are tabulated in Table 2.

BIfurcatIon dIagraM calculatIon
The first step to calculate bifurcation diagrams is to identify the system steady states. Given the model complexity it was not possible to get an analytic expression whose roots rendered the steady states. Thus we proceeded as follows: corresponds to the bacterial growth rate. In this work we consider a variable growth rate, which depends on the glucose uptake flux (20) and on the lactose metabolism rate (21) as given by Equation (8); see Santillán (2008) for a derivation of this equation.
molecule count of all chemical species in one (the other) of the chosen initial states is smaller (larger) than that of the uninduced (induced) steady state, we can assure that these initial conditions lie in different basins. The method accuracy depends on the chosen ∆E threshold and on the total simulation time; the smaller the threshold and the longer the simulation time, the more accurate the method. We carried out a numerical convergence analysis in a few selected cases and found that decreasing the threshold and increasing the simulation time beyond the reported values does not change the obtained bifurcation diagrams.

stochastIc sIMulatIons
To carry out stochastic simulations of the lac operon dynamics we employed a variation of the Gillespie algorithm introduced by Cao et al. (2005), which we implemented in Octave. The algorithm proceeds as follows. First, we assumed that transcription, mRNA degradation, translation, and enzyme degradation are slow reactions as compared with all the other reactions in the operon regulatory pathway. The stoichiometry of these reactions is while the reaction propensities are given by Once reaction scheme is determined, the following steps are iteratively carried out: values, the steady state solutions calculated in the previous steps shall be the same, otherwise the system is bistable. In particular, we considered that the system is monostable whenever the difference between the stationary enzyme counts resulting from the simulations in steps 1 and 2 (∆E) is smaller than 10 −2 mpb. This method detects bistability whenever it exists for the following reasons. If the system is monostable, it will converge to the only available steady state regardless of the initial condition. However, if the system is bistable, each of the available stable steady states will have its own basin of attraction. Therefore, if we pick two initial conditions, one in each basin, and solve the model equations, the system will evolve to two different steady states. Finally, since the We further simulated induction of the operon with lactose, solely. To do that we set [T e ] = 0 and repeated the calculations described in the previous paragraph. The resulting bifurcation diagram is plotted in Figure 3B. Note that the bistable region disappears via a cusp bifurcation for extremely low [G e ] values. On the other hand, when compared with that obtained with TMG induction, much larger [L e ] values are needed to enter the bistability region. As a matter of fact, the lower [L e ] bifurcation value slowly grows as [G e ] increases until about 200 μM, and reaches a plateau afterward. Conversely, the upper [L e ] boundary of the bistability region increases steadily for all the analyzed [G e ] values. These results corresponds to those reported in Figure 3F of Santillán (2008).
Next, we investigated the effect of using a mixture of lactose and TMG (to induce the lac operon) on the system bistable behavior. We These diagrams are qualitatively similar to that in Figure 4A in the sense that the cusp catastrophe disappears prolonging the bistable region down to [L e ] = 0. However, the bistability region for [L e ] < 100 μM gets wider as [T e ] increases. These results are in agreement with those of Narang and Pilyugin (2008), who studied induction of the lac operon when the bacterial population grow in a mixture of lactose and glucose, and concluded that the lac operon is much more prone to bistability if the medium contains carbon sources that cannot be metabolized by the lac enzymes.
To account for the intrinsic noise caused by the stochastic timing of transcription, translation, and mRNA and enzyme degradation, we carried out stochastic simulations of the operon dynamics as described in Section "Materials and Methods." First, we set [T e ] = 130 μM, as in Figure 4B, and performed simulations for the [G e ] and [L e ] values corresponding to the 12 different points shown there. For each simulation we calculated the histogram of the number of permease molecules and the results are shown in Figure 5. In all cases, we took four million simulation steps to guaranty that the system switches several times between the uninduced and the induced states, whenever the ([G e ],[L e ]) point lies inside the bistability region. In fact, we chose such simulation length after performing a numerical convergence analysis and noticing that the shape of all histograms does not change after one million steps.
We results Following Santillán (2008), we started by simulating the experiments reported by Ozbudak et al. (2004), who quantified the lac operon bistable behavior when it is induced with TMG. To calculate the corresponding bifurcation points we set [L e ] = 0 and varied [G e ] and [T e ] as described in the Methods. The resulting bifurcation diagram is plotted in Figure 3A in the [T e ] vs. [G e ] parameter space, together with the experimental points of Ozbudak et al. (2004). This diagram is equivalent to that in Figure 3B of Santillán (2008). The agreement between our results and those of Ozbudak et al. (2004) validate both the model and the numerical procedure to calculate the system bifurcation diagrams proposed in Section "Materials and methods." On the other hand, observe in Figure 3A that bistability is present for all G e values, and that both the low and high [T e ] boundaries increase as [G e ] increases, reaching a plateau at about 200 μM of glucose. In particular, we have that the system becomes bistable for very low TMG levels when [G e ] < 5 μM, as reported by Ozbudak et al. (2004). In order to contribute to the clarification of the above referred controversy, we have extended the model developed by Santillán (2008) for the lactose operon in E. coli to account for the presence of both lactose and TMG in the bacterial growing medium. Our aim is to make quantitative predictions regarding the bistable behavior of this system, which can be experimentally tested to validate the model itself. We have paid special attention to the estimation of all the model parameters from reported experimental data.
To test the extended model we reproduced the bifurcation diagrams calculated by Santillán (2008) for when the lactose operon is induced with either TMG or lactose. Our results are completely equivalent to those of Santillán (2008), as one can verify by comparing the bifurcation diagrams in Figures 3A,B with the corresponding diagrams in Santillán (2008).
Afterward, we investigated the effect of inducing the lac operon of E. coli with a mixture of lactose and TMG. We did this by using both a deterministic model and stochastic simulations. The results are summarized in Figures 4-6. When comparing with the bistability diagram corresponding to induction by lactose - Figure 3B, we can see that the addition of enough TMG makes the cusp catastrophe disappear and prolongs the bistability region down to 0 in the [L e ] axis -see We can see in Figures 5 and 6 that, as expected, all the histograms corresponding to ([G e ], [L e ]) points lying outside of bistability regions are monomodal. On the other hand, not all the histograms corresponding to ([G e ], [L e ]) points within the bistability region are clearly bimodal. In some cases the histogram bimodality is blurred by the wide dispersion of one or the two modes. Finally, we observe that all histograms corresponding to bistable points located close to the uninduced monostable region have a distinct bimodal shape. This suggests that the uninduced steady state is more strongly stable in all cases.

dIscussIon and conclusIons
The current understanding of the lac operon invokes bistability to explain how E. coli consumes glucose and lactose in the most efficient possible way (Lehninger et al., 2005). Namely, the genes necessary to uptake and metabolize lactose are only turned on when this sugar is present (at high levels) in the growing environment, while glucose (the bacterium favorite carbon source) is absent. Nonetheless, recent experimental and modeling works have open a controversy regarding the existence of bistability in this system under natural conditions. While van Hoek and Hogeweg (2006) propose that the wild-type lac operon is not a bistable switch when the bacterial population grows in a lactose rich medium for several generations, Santillán (2008) argue that bistability is present under such conditions and plays an important evolutionary role.  Figures 5 and 6 were carried out. The stairlike appearance of the borderlines is due to the resolution of our numerical calculations.
consumption of glucose and tryptophan in E. coli (Konan and Yanofsky, 1997), and of glucose and galactose in S. cerevisiae (Acar et al., 2005;Bhat and Iyer, 2009), have the same architecture as the lac operon of E. coli. That is, the regulatory pathway controlling the expression of the genes involved in the transport and metabolism of carbon sources other than glucose seems to be conserved over different metabolic networks and over different species. Thus, a complete understanding of the lac operon dynamics may help to understand the performance of a whole class of gene networks. Finally, experiments similar to the ones here proposed may help to elucidate the existence of bistability in other gene networks with the same architecture.

Interaction of the lac repressor with a single operator
The lac repressor is a homo-tetramer made up of two identical functional dimers. Each monomer in the dimer can be bound by either an allolactose (A) or a TMG (T) molecule. When all The results described in the previous paragraph predict that bistability can be found at much lower [L e ] values when a mixture of lactose and TMG is employed, as compared with the situation where lactose is the only inducer. Moreover, the novel shape of the bistability region also implies that bistability could be experimentally tested with the setup of Ozbudak et al. (2004) by fixing the extracellular glucose and TMG concentrations in the growing medium and varying the level of lactose, but also by fixing the extracellular lactose and TMG levels and varying the glucose concentration. If validated experimentally, these predictions would speak in favor of the Santillán (2008) model, as well as of his conclusions regarding bistability in the lac operon of E. coli. The fact that some bistable points, according to the deterministic model, are hard to identify through stochastic simulations does not invalidate the previously proposed experiments: smaller bistability regions would be found than those predicted by the deterministic model, but they could certainly be recognized.
The importance of unveiling whether the lac operon of E. coli is bistable or not under natural conditions transcends this system. The reason is that the gene networks responsible for the alternate

DNA folding
There are three different operators in the lactose operon. A repressor molecule with at least one free dimer can bind any one of these operators. Furthermore, DNA can fold in such a way that a free repressor molecule can bind two different operators simultaneously.
the monomers in the repressor are free, each dimer can bind an operator and so the repressor can bind two operators simultaneously. When one or the two monomers in one dimer of the repressor are bound by either A or T, and the opposite dimer is free, the repressor can still bind an operator. However, when two or more monomers in opposite sides are bound by A or T, the repressor cannot bind any operator.
In what follows R represents a repressor completely free from allolactose and TMG, R A (R T ) denote a repressor bound by a single allolactose (TMG) molecule, R 2A (R 2T ) indicates a repressor bound by two allolactose (TMG) molecules, and finally R AT corresponds to a repressor bound by one allolactose and one TMG molecule. We followed Santillán (2008) to found the fractional concentrations of all these complexes. That is, we wrote down all the chemical equations through which these complexes are formed, as well as the equilibrium equation for each one of them. Afterward, we wrote conservation equation for the total repressor concentration. And

( ))
Then, the transcription initiation rate at the lac promoter is: ( ) corresponds to the state where operators i and j are bound by the same repressor, while operator k is bound by another repressor. Santillán (2008) deduced expressions for the probability of all possible operator states when the lac operon is induced solely with lactose. Following an analogous procedure we found the probabilities for all the operator states, taking into consideration that each binding site in a repressor molecule can be bound by either an allolactose or a TMG molecule. In the forthcoming paragraph we present the results, but omit the algebraic procedure for the sake of brevity.  ( , , ) ; (A7) the fraction of operators in state O R i is: Since TMG is not metabolized by β-galactosidase, the corresponding term is absent.

paraMeter estIMatIon
Given that the present model is a generalization of that introduced by Santillán (2008), we took most of the parameter values estimated in such reference. It must be emphasized that such parameters values were all estimated from data reported for E. coli.
Nonetheless, we found more accurate experimental dates we regarding the maximum initiation transcription rate, k M , and thus we re-estimated this parameter as follows. Russell and Bennett (1982) report that the trp promoter is approximately five-fold stronger than the lac promoter. On the other hand, according to Santillán and Zerón (2004) the maximum initiation transcription rate of the trp operón is 5.1 min −1 . From this data we estimate k M = 1.0 min −1 .
Having a new value for k M obligates us to re-estimate parameter K T by fitting the bifurcation diagram predicted by our model, with only TMG like inductor, to the experimental results of Ozbudak et al. (2004). The value we obtained through this procedure is K T = 1.7 × 10 7 mpb.
On the other hand, since we are considering both lactose and TMG, there are parameters in this model which do no appear in the model by Santillán (2008). These parameters are the affinity of lactose to the repressor (K A ), the affinity of TMG to the lac permease (κ T ), and the maximum rates of TMG and lactose uptake per permease (k T ). Given that TMG is a chemical analogue of lactose, we simply assumed that K A = K T , κ A = κ T , and k A = k T .

acknowledgMents
The authors are thankful to Eduardo S. Zeron for his fruitful advice. This research was partially supported by Consejo Nacional de Ciencia y Tecnología (CONACyT, Mexico) under grant 55228.

Let [G e ], [L e ], [T e ], [L], [A]
, and [T] respectively represent the concentrations of external glucose, external lactose, external TMG, intracellular lactose, intracellular allolactose, and intracellular TMG. Following Santillán (2008), the efficiencies of β-permease to transport lactose and TMG respectively are written as: Lactose and TMG transport by β-permease is inhibited by external glucose (a phenomenon know as inducer exclusion). The decrease on the uptake rate as a function of [G e ] can be modeled by (Ozbudak et al., 2004): ] .
e e e = − + 1 Let k L and k T be the maximum rates of lactose and TMG uptake, respectively. Then, the corresponding uptake rates are: Once lactose is internalized, it is metabolized by β-galactosidase. Approximately half of the lactose molecules are directly converted into glucose and galactose, while the rest are turned into allolactose. Allolactose is later converted into glucose and galactose by the same enzyme (β-galactosidase). Following Santillán (2008), and taking into account the considerations above, the dynamics of intracellular lactose can be modeled by the following equation: is the rate or lactose metabolism into glucose and galactose, and the factor 2 accounts for the fact that lactose is metabolized also into allolactose at approximately the same rate. We can also assume according to Santillán (2008) that