Bistability in the Lactose Operon of Escherichia coli: A Comparison of Differential Equation and Boolean Network Models
Raina Robeva* and Necmettin Yildirim†, *Department of Mathematical Sciences, Sweet Briar College, Sweet Briar, VA 24595, USA, †Division of Natural Sciences, New College of Florida, FL 34243, USA
2.1 Introduction
Cellular systems are complex and involve many interacting components, which are ultimately responsible for cellular functions. Such complexity is capable of producing interesting dynamic behaviors. Among such behaviors, bistability is extremely important and is a recurring pattern in biological systems. A system is called bistable if it is capable of resting in two stable steady states separated by an unstable steady state. A bistable system can reside in one of the two states for a long time unless a strong external perturbation kicks in. Bistability also provides a perfect discontinuous switching among the stable steady states and can convert graded inputs into switch-like outputs.
Bistability can be either reversible or irreversible. A system is irreversibly bistable if it cannot go back to its earlier state once it settles down to a state after a perturbation. However, a reversible bistable system can go back and forth between the stable states as the system parameters vary. Figure 2.1a shows a reversible bistable switch and Figure 2.1b is for an irreversible bistable switch. As seen in Figure 2.1a, the value of the parameter 
for the system to switch from a lower steady state to a higher steady state is different from the value of the parameter 
for which a switch from a higher steady state to a lower steady state occurs. This is called hysteresis. Bistability comes with hysteresis, which can be defined as dependence of the current state of a system not only on its current state but also on its past state. To predict future state of a bistable system, both its current state and its history must be known.
Figure 2.1 Bistability and hysteresis emerging in cellular reaction networks. Bistability can be reversible (a) or irreversible (b). In reversible bistable systems, the steady states can change from low to high or high to low state as the parameter p is tuned. However, the steady state changes from a low state to a high state at a higher value of p than that for which the change from a high steady state to a low steady state occurs. In an irreversible bistable system, once a system shifts from a low state to a high state, it stays there even when the value of the parameter p is set back to a value much lower than that where an earlier shift took place.
In the case of the lactose (lac) operon in Escherichia coli (E. coli), it has been known since the 1950s [1] that it exhibits a hysteresis effect: In the absence of glucose, the operon is uninduced for low concentrations below a threshold 
of extracellular lactose and fully induced at high concentrations above a threshold of 
. In the interval 
, both induced and uninduced cells can be observed and their status depends on the cell history (the system response is hysteresis). The interval 
is defined as a region of bistability and is referred to as maintenance concentration. Cells grown in an environment poor in extracellular lactose will have low levels of internal lactose and allolactose and may remain uninduced for lactose levels in the interval 
. In contrast, cells grown in a lactose-rich environment remain induced for concentrations in the interval 
. A recent discussion of the bistability feature of the lac operon together with green-fluorescence and inverted phase-contrast images of a cell population showing a bimodal distribution of lac expression levels can be found in [2]. A survey of the quantitative approaches to the study of bistability in the lac operon is given in [3].
In this chapter we examine two types of mathematical models of the lac operon: differential equation models and Boolean network models. In both cases the focus will be on bistability and on the capability of the models to capture the bistable behavior of the system. The chapter is organized as follows: Section 2.2 presents a short description of the regulatory mechanism of the lac operon. In Section 2.3 we give a step-by-step tutorial on how to model biochemical reactions using differential equations. These approaches are then used in Section 2.4 to justify several widely used differential equation models of the lac operon. Section 2.5 extends the material from Chapter 1 on Boolean networks, describing Boolean models with varying elimination times and Boolean models involving delays. In Section 2.6 we build Boolean model approximations of the differential equation models from Section 2.4 and demonstrate that those models capture the bistability of the system. Section 2.7 contains some closing comments and conclusions.
2.2 The Lactose Operon of Escherichia Coli
The lac operon is a well-known example of an inducible genetic circuit. It has been serving as a model system for understanding many aspects of gene regulations since its discovery in the late 1950s. The lac operon encodes the genes for internalization of extracellular lactose and then its conversion to glucose. A cartoon that depicts the major regulatory components of this system is shown in Figure 2.2. The lac operon consists of a promoter/operator region 
and three structural genes lacZ, lacY, and lacA. A regulatory gene lacI (I) preceding the lac operon is responsible for producing a repressor 
protein. The molecular mechanism of this operon is as follows: In the presence of allolactose, a complex between allolactose and the repressor is formed. This complex is unable to bind to the operator region, which allows the RNA polymerase to bind to the promoter site, starting initiation of transcription of the structural genes to produce mRNA
. However, in the absence of allolactose 
the repressor protein R binds to the operator region (O) and blocks the RNA polymerase from transcribing the structural genes. The translation process takes place once the mRNA has been produced. The lacZ gene encodes the portion of the mRNA for the production of the enzyme 
-galactosidase 
and translation of the lacY gene produces the mRNA that is needed for the production of the enzyme permease 
. The final portion of mRNA produced by the lacA gene is responsible for the production of the enzyme thiogalactoside transacetylase, which does not play a role in the regulation of the lac operon [4].
Figure 2.2 Schematic representation of the lactose operon regulatory system. See the text for details. Reprinted from Methods in Enzymology, v. 484, Yildirim, N. and Kazanci, C., Deterministic and stochastic simulation and analysis of biochemical reaction networks the lactose operon example, p. 371–395, (2011), with permission from Elsevier.
A bistable system often involves a positive feedback loop or a double negative feedback. Such motifs are very common in control mechanisms in gene networks. The lactose operon and the arabinose operon of Escherichia coli both have positive feedback mechanisms [5,6]. In the lactose operon, once external lactose is converted into allolactose, the allolactose feeds back and forms a complex with the lactose repressor, which allows the transcription process to start. This eventually leads to synthesis of more enzymes to transport more external lactose into the cell and then convert it into allolactose. More information about this mechanism and about the lac operon in general can be found in Chapter 1.
2.3 Modeling Biochemical Reactions with Differential Equations
Differential equations are widely used to study reaction kinetics. In this section we outline the basics of reaction modeling with differential equations. This approach will be applied in Section 2.4 to justify two well-known differential equation models of the lac operon.
Two or more species can react if they come together and collide. Collision alone is not enough for a reaction to occur since the molecules have to collide the right way and with enough energy. Even if the molecules are orientated properly when they collide, a reaction may not happen if the molecules do not have enough energy. On the other hand, if the molecules collide with enough energy but their orientations are not appropriate, then the reaction still may not occur. Reactions with a single species, such as degradation reactions, do not require orientation of collisions. Only if they have enough energy will those reactions take place. Most of the reactions involve one or two species. Reactions with three or more species are extremely uncommon.
The rate of a reaction reflects how fast or slow it takes place. There are different approaches and methodologies to studying reaction rates. Mass-action kinetics is a microscopic approach to reaction kinetics and it is one of the most common methods. Modeling based on mass-action kinetics results in a system of differential equations as a mathematical description of the reaction rates. This approach is fully deterministic and it is only appropriate under certain conditions, such as constant temperature, well-mixed solution, large number of molecules, etc. In this section, we briefly describe how to write differential equation models from a reaction mechanism.
Reactions can be classified based on their reactants. The reaction 
is called uni-molecular, since one reactant A becomes a product P. The reaction 
is called bi-molecular, since there are two reactants A and B becoming a product P. A tri-molecular reaction looks like 
, and has three reactants A, B, and C.
According to the mass-action kinetics, a reaction rate is proportional to the probability of collision of the reactants involved. At higher concentration, the collisions occur more often. This probability can be taken to be proportional to each reactant concentration.
Now consider a uni-molecular reaction in which A becomes P,
(2.1)
Here k is called a kinetic rate constant describing how likely it is for this reaction to occur and produce the product P. According to mass action kinetics, the rate of this reaction can be written as
Here 
and 
represent concentrations of A and P, respectively; k is a first order rate constant. Units for 
and 
are concentrations and unit of 
is concentration per time. Therefore, unit of the rate constant k has to be 
.
A bi-molecular reaction looks like,
(2.2)
In this reaction the reactants A and B react and become a product P with a rate constant k. The rate constant quantifies how likely it is that the collision of A and B ends up with the product P. The rate of this reaction is
(2.3)
Here k is a second order rate constant and its unit becomes 
.
The reactions given in both Eq. (2.1) and Eq. (2.2) are uni-directional reactions. In theory, all chemical reactions are reversible. Now let’s consider a two-directional reaction like
(2.4)
In this reaction, A and B react and become P with an associated rate constant 
. On the other hand, P can break down and produce A and B. The rate constant for the backward reaction is 
. Here P is produced at a rate 
and consumed at a rate 
. Hence, how the concentration of P changes over time is given by the difference of these two rates:
(2.5)
(2.6)
In this equation 
is a second order rate constant and 
is a first order rate constant.
2.3.1 Enzymatic Reactions and the Michaelis-Menten Equation
Enzymes are specific proteins that catalyze reactions. An enzyme can increase the rate of a reaction up to 
-fold [7], compared to the spontaneous reaction without the enzyme. The enzyme first binds to its substrate (reactant), forms a complex (an enzyme-substrate complex), and performs a chemical operation on it. Then it releases from the complex, resulting in conversion of the substrate into a product. Enzymes stay unchanged after the reaction. Some enzymes bind to a single substrate while others can bind to multiple substrates and combine them to produce a final product.
Consider now the enzyme catalyzed reaction in Eq. (2.7), which involves three individual reactions.
(2.7)
The first reaction is 
, where an enzyme E binds to a substrate S and forms an enzyme-substrate complex ES with an associated rate constant 
. Since this is a reversible reaction, ES can break down into E and S 
. The associated rate constant for this backward reaction is 
. The third reaction is 
in which the enzyme E releases from ES, producing a product P with a rate constant 
.
The differential equation describing the dynamics of the concentration of the enzyme-substrate complex ES is the difference between the gain and loss terms. ES is produced with the first reaction and consumed with the second and third reactions. Hence, we have
(2.8)
The dynamics of the product P are modeled by Eq. (2.9). This equation has a single term, since P is produced by the third reaction and is not consumed by any of the reactions.
(2.9)
If the total concentration of the enzyme stays constant over the duration of this reaction, we can write,
(2.10)
where 
represents the initial enzyme concentration.
Equations (2.8)–(2.9) now describe the dynamics of the single enzyme single substrate reaction in Eq. (2.7). However, with an additional assumption, these two equations can further be simplified. In an enzymatic reaction, the enzyme-substrate complex reaches a steady state much earlier than the product does. This allows us to take 
. Therefore,
(2.11)
Solving Eq. (2.11) for 
yields
(2.12)
After plugging 
from Eq. (2.12) into the enzyme conversation equation given by Eq. (2.10) and solving the resultant equation for 
, we obtain 
in terms of only 
:
(2.13)
Substituting 
given by Eq. (2.13) into Eq. (2.9) yields
(2.14)
where
(2.15)
are two positive parameters. This equation is called the Michaelis-Menten equation. The right-hand side of this equation 
is an increasing function of the substrate concentration with two properties: (i) 
, and (ii) 
. From a biological point of view, this reaction never occurs at a rate greater than 
. 
is the substrate concentration at which the reaction rate is equal to half of its maximum value 
. Furthermore, since 
, the reaction rate is a linear increasing function of the initial enzyme concentration 
, which means that higher initial enzyme concentration makes the reaction go faster and this relationship between the rate and 
is linear.
If the total substrate concentration is conserved throughout the course of the reaction, we can write
(2.16)
Since the initial concentration of the enzyme is usually much smaller than the initial substrate concentration, the concentration of ES is negligibly small compared to either the substrate or the product concentration. Hence, (2.16) simplifies to
(2.17)
Differentiating both sides of this equation with respect to t yields 
. Hence, Eq. (2.14) becomes
(2.18)
The evolution of the substrate concentration can be simulated by solving this differential equation after fixing 
and 
and assigning a value to the initial concentration of the substrate. Then, from Eq. (2.17), one can compute the progress curve for the concentration of the product.
Exercise 2.1
Consider the reactions where two substrates 
and 
compete for binding to an enzyme E to produce two different products 
and 
. Under the assumption that each reaction follows the Michaelis-Menten kinetics, derive rate equations for 
and 
in this system and explain the effects of the competition occurring.
2.3.2 Multi-Molecule Binding and Hill Equations
Some enzymes react with more than one substrate molecule, e.g., the reaction in Eq. (2.19), in which 
-molecules of a substrate react with an enzyme E. For simplicity, let’s assume that the binding sides on the enzyme are all identical, the bindings take place simultaneously and the bindings of 
-molecules are independent. Otherwise, the influence of each binding step to other binding steps has to be considered.
(2.19)
This reaction consists of three individual reactions. In the first reaction, 
-molecules of a substrate S interact with one molecule of an enzyme E to form an enzyme-substrate complex 
with a rate constant 
. 
is the second reaction, in which the enzyme-substrate complex 
can break down into the enzyme E and 
-substrate S molecules with a rate constant 
. In the third reaction, 
, the enzyme E releases unchanged from the complex and a product P is produced with a rate constant 
. The differential equation modeling dynamics of the concentration of the enzyme-substrate complex 
is the difference between the production and consumption rates of this complex. 
is produced with the first reaction and it is used up with both the second and the third reactions. Hence,
(2.20)
Since the product P is produced by the third reaction and it is not consumed, the dynamics of the product P concentration can be modeled as in Eq. (2.21).
(2.21)
Suppose that the complex 
reaches a steady state faster than the product P. Then, Eq. (2.20) becomes
(2.22)
Solving Eq. (2.22) for 
gives
(2.23)
Assuming that the total amount of enzyme is conserved throughout the reaction, we have
(2.24)
where 
stands for the total enzyme concentration. Substituting Eq. (2.23) into Eq. (2.24) and then solving it for 
, we obtain
Then the dynamics of the concentration of P given by Eq. (2.21) becomes
(2.25)
where 
and 
are positive constants. The parameters 
and 
are as in Eq. (2.15), and n is a Hill coefficient. This equation is known as the Hill equation. Unlike the Michaelis-Menten equation, the Hill equation has three adjustable parameters. The right-hand side of the equation 
is a strictly increasing function of 
and satisfies the following two conditions: (i) 
, and (ii) 
. Biologically, this reaction can never take place at a rate larger than 
. The value 
is the substrate concentration necessary for this reaction to occur at a rate equal to half of the maximum rate, 
. Since 
, increasing the initial enzyme concentration linearly increases the reaction rate. When 
, the Hill function behaves like a step function. For example, for large n values, 
when 
is small. On the other hand, 
when 
is big (see Figure 2.3). When the initial substrate concentration is conserved throughout the course of the reaction and the initial enzyme concentration is much smaller than that of the substrate concentration, Eq. (2.25) can be written in terms of the substrate consumption rate. That is, when
(2.26)
where 
is the initial substrate concentration, Eq. (2.25) becomes
Notice that the Hill function reduces down to a Michaelis-Menten function when 
. Both the Hill function and the Michaelis-Menten function are commonly used in modeling of biochemical reaction networks.
Exercise 2.2
The Hill equation is an approximation for multi-molecule binding and it assumes simultaneous binding of 
-molecules of a substrate S to the enzyme E above. Now suppose that two molecules of the substrate S are undergoing a reaction with an enzyme in an ordered manner as follows:
(2.27)
Derive a rate equation under the steady state assumption and compare it with the Hill equation given in Eq. (2.25). When do these two equations become roughly the same?
Figure 2.3 Representative curves from the Hill equation given by Eq. (2.25) for various values of the Hill coefficient 
when both 
and 
are constant. The curves are all sigmoidal and they behave like a step function when n is big (solid curve).
2.4 The Yildirim-Mackey Differential Equation Models for the Lactose Operon
In this section, we talk about two delay differential equation models [8,9] recently developed to study the lac operon. Both models employ experimentally estimated parameters and are compared against independently collected, published experimental data from this genetic circuit. The first model is a 3 variable model and only takes into account mRNA (M), 
-galactosidase (B), and allolactose (A) dynamics. Internal lactose (L) is assumed to be readily available and is included in the model as a parameter [8]. In the second model, besides these three dynamic variables, the internal lactose (L) and permease (P) are also explicitly modeled [9].
2.4.1 Model Justification
We need a few more preliminaries before presenting the differential equation models.
Modeling dilution in protein concentration due to bacterial growth: When developing a mathematical model for bacterial systems like the lac operon, it is important to know how fast the bacteria actually grow. Depending on the environmental conditions, an E. coli population can double in size as fast as every 20 min [10]. This increase in the cellular volume results in a dilution in the protein concentrations, which can be critical in developing reliable mathematical models. Dilution in protein concentrations due to bacterial growth can be modeled as follows: Suppose that V is the average volume of a bacterial cell and x represents the number of molecules of a protein X in that cell. Assume that the cell volume increases exponentially in time. Then we have,
(2.28)
where 
denotes the growth rate. Also, assume that the degradation of X is exponential. Then we can write,
(2.29)
where 
represents this decay rate. Concentration of the protein X is equal to the number of molecules x of X divided by volume of the cell, namely
(2.30)
where 
stands for concentration of the protein X. Differentiating both sides of Eq. (2.30) with respect to time t results in
(2.31)
Substituting Eq. (2.28) and Eq. (2.29) into Eq. (2.31) and making necessary simplifications gives us
(2.32)
As seen in this equation, an increase in the volume over time can successfully be modeled by adding the growth rate to the protein degradation rate.
Modeling of the lactose repressor dynamics: The repressor protein plays a major role in the regulation of the lac operon. In both absence and presence of the external lactose, the repressor protein is produced at a constant rate. A binary complex between allolactose and the repressor protein is formed. This complex cannot bind to the operator site (O). Yagil and Yagil [11] have modeled the dynamics of the repressor protein (R) as follows:
(2.33)
Here n represents the effective number of allolactose molecules required to form this complex. If we assume this reaction is at equilibrium, we have
(2.34)
where 
is the association constant for this reaction. The repressor molecules can also bind to the operator region (O) in the absence of allolactose and block the transcription process. We assume this interaction is of the following form:
(2.35)
At equilibrium, we have
(2.36)
Here 
is also an association constant of this reaction. Let 
be the total operator concentration. It is plausible to take the total concentration of the operator as constant. Therefore,
(2.37)
From Eq. (2.36) and (2.37), we can write
(2.38)
Now, let 
be the total concentration of the repressor protein. Since the repressor protein is not regulated by the extracellular lactose, we can assume that the total concentration of this protein also stays unchanged. Hence,
(2.39)
Since there are at most a few molecules of the operator site per cell, it is reasonable to assume 
. Then Eq. (2.39) simplifies to
(2.40)
Putting Eq. (2.34) into Eq. (2.40) and solving the resultant equation for 
, we obtain
(2.41)
By plugging the repressor protein concentration from Eq. (2.41) into Eq. (2.38), we obtain the portion of free operator site concentration in terms of the allolactose concentration as
(2.42)
where 
. This equation describes the available operator concentration for RNA polymerase binding and for initiating the transcription of the structural genes. The right-hand side of the function in Eq. (2.42) is a sum of two functions: (i) a decreasing Hill function, 
, and (ii) an increasing Hill function 
of the form given in Eq. (2.25). Denote 
and notice that f takes a positive value 
when 
. Since f is an increasing function of 
and the concentration of A is always non-negative, this is the lowest value f can take. Biologically this can be interpreted as the process of ongoing transcription of the structural genes occurring at a constant rate in the absence of allolactose, and providing basal concentrations of the proteins encoded by the structural genes. On the other hand, when allolactose is abundant, 
. This can be interpreted as having all of the operator sites available to mRNA polymerase for transcription of the structural genes, with transcription occurring at the maximal possible rate. Yagil and Yagil [11] have experimentally shown that Eq. (2.42) can accurately fit their observations and estimated the parameters K, 
, and n from the data.
The Yildirim-Mackey models of the lac operon: The dynamics of mRNA is modeled with the following equation
(2.43)
where 
and 
. The rate of change in mRNA concentration is the difference between its gain and loss terms. We assume that the production rate of mRNA is proportional to the fraction of free operators as described in Eq. (2.42). Here 
represents such a proportionality constant. 
represents the allolactose concentration at time 
. The production of mRNA from DNA via transcription is not an instantaneous process. In fact, it requires a period of time 
for RNA polymerase to transcribe the first ribosomes binding site. Hence, the production rate of mRNA is not a function of the available allolactose concentration at time t, but a function of the available concentration of allolactose at time 
. The exponential prefactor 
accounts for the dilution of the allolactose through bacterial growth during the transcriptional period, which can be derived as follows: Suppose that a protein P decays exponentially. If p represents the concentration of this protein, then we have
(2.44)
where 
quantifies how fast this degradation happens. When this equation is integrated from 
to t, we obtain
(2.45)
where 
. This derivation explains the prefactor 
in Eq. (2.43).
in Eq. In the 5-variable model, the mRNA equation has an extra parameter 
that accounts for the basal transcription rate in the absence of the extracellular lactose. Therefore, we use Eq. (2.46) for the dynamics of mRNA concentration in the 5 variable model.
(2.46)
The basal transcription rate for mRNA is not explicitly modeled in Eq. (2.43) but is included in the parameters 
and K, instead. The second term in Eq. (2.43) models the loss in mRNA concentration. It is a sum of two terms and given by 
. The term 
accounts for the loss due to mRNA degradation and the term 
is the effective loss due to the bacterial growth, as explained in Eq. (2.32).
Equation (2.47) models the dynamics of the 
-galactosidase concentration. It can be assumed that the rate of production of 
-galactosidase is directly proportional to the mRNA concentration at time 
with a proportionality constant 
.
(2.47)
Here 
is the time required for the translation of mRNA, 
and 
.
Equation (2.48) governs the allolactose dynamics. The first term on the right-hand side of this equation gives the gain in allolactose due to the conversion of lactose mediated by 
-galactosidase. We assume that this conversion follows the Michaelis-Menten equation as derived in Eq. (2.14).
(2.48)
The second term denotes the loss of allolactose due to its conversion to glucose and galactose by 
-galactosidase, which is also assumed to follow a Michaelis-Menten kinetics. The last term accounts for the loss in the allolactose concentration due to its degradation and dilution. The equations for the 3-variable model are presented in Table 2.1.
In addition to the equations in Table 2.1, the 5-variable model includes specific equations for the dynamics of permease and internal lactose. The equation for the dynamics of the permease is similar to the 
-galactosidase equation and it is given in Eq. (2.50).
(2.50)
Here 
and 
. The first term models the gain in the permease concentration due to the transcriptional induction of lacY gene by the external lactose, and it is assumed to be proportional to the mRNA concentration at a time 
ago. Here 
is the proportionality constant for this process, and 
and 
are the times required for the translation of functional 
-galactosidase and permease, respectively. The term 
models the dilution in mRNA concentration due to cellular growth during the translation of such enzymes. The last term is for the loss in the permease concentration due to cellular growth and degradation of this enzyme.
The equation governing the dynamics of internal lactose is given in Eq. (2.51). The first three terms in this equation are in the Michaelis-Menten form. The first term models the gain in the concentration of lactose due to the transport of external lactose by the permease.
(2.51)
Here 
denotes the concentration of external lactose, which is a parameter for the 5-variable model. The permease-mediated transport of lactose through the cellular membrane is reversible [12]. The second term in Eq. (2.51) accounts for the loss in the internal lactose concentration due to the reversible nature of this transport process. The third term describes the loss in the lactose level due to its conversion to allolactose, which is mediated by 
-galactosidase. The last term models the loss in the lactose concentration due to the cellular growth and its degradation. The equations for the 5-variable model are presented in Table 2.2.
2.4.2 Numerical Simulation of the Yildirim-Mackey Models and Bistability
In this section, we will perform steady state analysis and numerical simulations for both models. The models consist of delay differential equations with discrete time delays due to the transcription and translation processes. Besides delay parameters, they also involve a number of unknown parameters that need to be fixed in order to perform steady state analysis and numerical simulations. Yildirim et al. [9,8] did an extensive literature search to estimate these parameters from published articles. The values and details on estimation of the parameters can be found in [9,8,13]. In our simulations in this chapter, we have used two different 
values, as was done in the original papers. For simulations with the 3 variable model we took 
. For the simulations with the 5 variable model, a smaller 
value was used, 
. These values were estimated by fitting the differential equation models to experimental data.
Experimental results showing that the lac operon in Escherichia coli is capable of showing bistable behavior for a range of extracellular lactose concentrations have been available since the late 1950s [1,14] and more recently in [2]. To see if the 3 variable model can show bistability, steady state analysis is performed, which can be studied by setting the left-hand side of each equation in the 3 variable model given in Table 2.1 to zero and solving it for a range of L concentrations after keeping all the other parameters fixed at their estimated values. The result is shown in (Figure 2.4). The 3 variable model predicts that there is a range for the internal lactose concentration, which corresponds to the S-shaped curve in the figure. When the lactose concentration is in this range, the lac operon can have three coexisting steady states.
Figure 2.4 Bistability arises in 
space in the 3 variable lac operon model (
-axis is in logarithmic scale). For a range of L concentrations there are three coexisting steady states for the allolactose concentration. We estimated this range to be (0.039,0.055) mM of L concentration.
Figure 2.5 shows how the bistable behavior arises in the time series simulation of the 3 variable model. For this simulation, all the parameters are kept constant at their estimated values when 
. As seen from the steady state curve in (Figure 2.4), there are three distinct steady states for this particular concentration of L. We calculated these steady state values numerically as in Table 2.3. The steady state values were calculated by solving the nonlinear system of equations obtained by setting time derivatives to zero in the model equations in Table 2.1 and solving it for the concentrations. Then six distinct initial values were chosen for the protein concentration around the unstable steady state concentration (steady state II in Table 2.3). Three of them were below the unstable steady state and the other three were above it. As expected, three initials converged to the lower steady state (steady state I in Table 2.3), the other three ended up settling at the higher stable steady state (steady state III in Table 2.3).
Exercise 2.4
The Lac repressor protein is a tetramer of identical subunits [11]. It has been shown experimentally that two allolactose molecules on average bind to this protein and effectively block it so that transcription of new proteins can take place in the presence of external lactose. Now suppose that three molecules of allolactose are needed for effective blockage of the repressor protein. Numerically study how this will affect the bistability range in this system. Use the 3 variable model and the parameter values given in TableI in the paper [8]. You should take the Hill coefficient as 
.
MATLAB starter code is provided for this exercise. Open the file Code_for_Ex_2_ 4_Starter.m and add the appropriate lines of code. Note that this exercise requires the use of the MATLAB’s Global Optimization Toolbox.
Figure 2.5 Time series simulation of the mRNA, 
-galactosidase and allolactose concentrations. These plots were produced by numerically solving the 3 variable model when 
. For this value of the internal lactose concentration, there exist three coexisting steady states (see Figure 2.4). The (
)’s in these plots represent the location of the low and high stable steady states. See the text for selection of the initials.
Table 2.3
The steady state values calculated by solving the nonlinear system of equations by setting time derivatives to zero in the 3 variable model when 
for which there exist three steady states.
The 5 variable model was analyzed in a similar way. Figure 2.6 shows a plot of the allolactose steady state values as responses to the external lactose concentration in the 5 variable model. To produce this plot, the system of nonlinear equations in Table 2.2 was solved for a range of 
values after setting each of the time derivatives to zero. As seen in (Figure 2.6), there is a range of 
concentrations for which there are three coexisting steady states for the allolactose concentration.
Figure 2.6 Bistability arises in 
space in the 5 variable lac operon model. Notice that the 
-axis is in the logarithmic scale and there exists a range of 
concentrations for which there are three coexisting steady states for the allolactose concentration. We estimated this range to be (0.027,0.062) mM of 
.
Table 2.4 shows the three steady states calculated from the 5 variable model when the extracellular lactose concentration is 
. To compute those steady states, the time derivatives in the model equations given in Table 2.2 are set to zero and the system of nonlinear equations is solved numerically while all the parameters are kept constant at their estimated values.
Table 2.4
The steady state values calculated from the 5 variable model when 
for which there exist three steady states (see Figure 2.4).
Figure 2.7 shows time-series simulations for the time evolution of mRNA, 
-galactosidase, allolactose, lactose, and permease concentrations in the 5 variable model when the extracellular lactose concentration is 
. There are three coexisting steady states for this value of 
(see Figure 2.6 and Table 2.4), two of which are locally stable and the third one is unstable [9,15]. For this simulation, we have chosen six different initial values for the protein concentrations in the neighborhood of the unstable steady state (steady state II in Table 2.4) when all the parameters were held at their estimated values in Table 2.4. Since we started our simulations from a point close to the unstable steady state (but not exactly from it), three of the initials converged to the higher steady state (steady state III in Table 2.4) and and the other three initials converged to the lower steady state (steady state I in Table 2.4).
Exercise 2.5
Consider the 5 variable model and the parameter values given in Table I in the paper [9] except for the bacterial growth rate 
. Bacterial growth rate can change depending on the environmental condition. In our analysis, we estimated this parameter to be about 30 min. Assume that you have a bacteria population that can double in size every 100 min then compute numerically the range for the extracellular lactose concentration and produce the bistability plot in 
space. Furthermore, compute and estimate the value for the bacterial growth rate from this model for which the lactose operon is no longer capable of showing bistable behavior.
MATLAB starter code is provided for this exercise. Open the file Code_for_Ex_2_ 5_Starter.m and add the appropriate lines of code. Note that this exercise requires the use of the MATLAB’s Global Optimization Toolbox.
Figure 2.7 Time series simulation of the mRNA, 
-galactosidase, allolactose, lactose and permease concentrations produced from the 5 variable model when the extracellular lactose concentration is 
. For this value of the external lactose concentration, there are three coexisting steady states (see Figure 2.6). Two of the three states are locally stable and the third one is unstable [15,9]. In the plots, 
’s represent the location of the low and high stable steady states. See the text for selection of the initials and details for this simulation.
2.5 Boolean Modeling of Biochemical Interactions
Boolean models were introduced and discussed in detail in Chapter 1 of this volume and we refer the reader to Chapter 1 for a primer on Boolean networks. We repeat some of the basics here for easy reference. In a Boolean model, all model variables are discretized to take values 1 or 0 (we also say that a Boolean variable can be On or Off, respectively). A certain threshold of discretization is chosen for each variable and a value of 0 typically represents the case when only “trace” (baseline) values of a substance are available. A value of 1 refers to concentrations larger than the threshold level.
As in the differential equation models, each variable in a Boolean model represents the concentration of a molecular species (e.g., enzyme, substrate, protein, or mRNA). A Boolean model of n-variables 
consists of n transition Boolean functions (also called update functions) 
describing the dynamical evolution of the model variables, where 
. Equations (2.53) present a simple example.
(2.53)
Boolean models are time-discrete and the dynamical evolution of the model is determined by iterating the transitions defined by the update functions. Starting from an initial condition 
at time 
, the state of the system at time 
will be 
where 
for 
. The transition between time t and 
is given by 
. As an example, consider the set of transition functions defined in Eqs. (2.53), with the initial condition 
. Substituting these values into Eqs. (2.53) yields
(2.54)
Next, the values 
are used to evaluate the transition functions again, producing 
. A subsequent iteration of the functions 
returns the same values 
. We say that we have computed the trajectory 
, and that (0,0,0,0) is a fixed point for the Boolean network. A visual representation of this process for all 16 different sequences 
composed of 0s and 1s as initial states is given by the directed graph in Figure 2.8. This graph depicts the transition diagram of the Boolean model. Loops of length larger than one on the transition diagram correspond to limit cycles. Figure 2.8 shows that the Boolean model from Eqs. (2.53) has two fixed points (0,0,0,0) and (1,1,1,1) and no limit cycles. The fixed points in a Boolean model are equivalent to the steady states of a differential equation model.
Figure 2.8 The state space diagram for the Boolean model from Eqs. (2.53).
In Chapter 1 of this volume we considered several models of the lac operon, specifically examining their transition functions, state-space diagrams and fixed points. A common feature of those models was the assumption that all major processes in the system’s regulation (e.g., dilution and degradation, translation and transcription, and the effects of interaction between the system’s components) take place in a single time step. Although models constructed under such assumptions can play an important role in the qualitative analysis of a system, such simple models lack the capability of describing delayed interactions and systems for which not all degradation rates are equal. Further, since those models assume only high and low levels of internal or external lactose, they effectively bypass the bistability phenomenon occurring at intermediate lactose concentrations.
We next consider Boolean models that can incorporate time delays, account for different degradation rates, and can exhibit bistability. The approach we outline below follows in most parts the general methodology proposed by Hinkelmann et al. in [16]. We begin with some background considerations.
1. Assume that a variable Y regulates the production of X.1 As in Chapter 1, we use 
and 
to denote the values of the Boolean variables X and Y at time t, where t is discrete, 
Assume that 
implies 
; that is, the presence of Y at time t implies that X will be present at time 
. Assume that the loss of X due to dilution and degradation occurs over the course of several time steps. To account for this process, n additional Boolean variables 
,…, 
, are introduced (where n is a positive integer) with the property that:
i. If 
and 
, then 
. A value of 1 for 
indicates that the amount of X present at time 
is already reduced once by dilution and degradation.
ii. If 
and 
, for some 
, then 
. A value of 1 for 
indicates that the amount of X present at time 
is already reduced i times by dilution and degradation.
iii. The number of required “old” variables is determined by the number of time steps needed to reduce the concentration of X below the discretization threshold when there is no new production of X. For instance, assume that in the absence of new production, the concentration of X needs to be diluted n times before falling below the discretization threshold. Thus, 
when either 
(that is, a new amount of X will be produced by time 
) or when 
(that is, previously available amounts of X are still available at time t and have not yet been reduced n-fold). In other words, 
. We should note that our approach here differs somewhat from that in [16]. We highlight the differences in Section 2.7.
2. To be able to model bistability we need to be able to distinguish between low, medium, and high concentrations of lactose L. This can be achieved by introducing an additional variable 
such that 
implies 
. In this framework, 
means that the concentration of L is at least medium, and when 
and 
, the concentration of L is only medium. 
and 
stands for low concentration of 
and 
represents high concentration of L.
3. In a Boolean model, delays can be incorporated by introducing additional variables. The number of additional variables depends on the magnitudes of the delays and on the choice of the time step for the Boolean model. Assume a variable R regulates the production of X but the effect that R exerts on X is delayed by time 
.2 If the delay 
is commensurable with n time steps in the model, n additional Boolean variables 
, can be introduced to represent the delayed regulation. The following motif will then be present in the set of transition functions: 
, …, 
. Expressed in terms of the transition functions, the same can be stated as 
.
We now use these techniques to create Boolean models that approximate the Yildirim-Mackey differential equation models from Section 2.4.
2.6 Boolean Approximations of the Yildirim-Mackey Models
2.6.1 Boolean Variants of the 3-Variable Model
We want to build a Boolean model based on the assumptions used in the 3-variable differential equation model from Table 2.1. Recall that 
. The delays are estimated in [8] to be 
. Various estimates are available from the literature for the loss rates 
, and 
and in some cases the range for the estimates may be rather wide. For instance in Yildirim and Mackey [9], the degradation rate of A is estimated at 
, in Yildirim et al. [8] the estimate is 
, and in Wong et al. [17] an estimate as low as 
is considered. The degradation rates 
and 
are estimated in both [8,9] to be 
. The effective loss in the concentrations due to dilution is proportional to the growth rate 
and is estimated to be between 
and 
.
As Boolean models are qualitative in nature, the actual values of these constants are not essential and the only consideration of importance is the comparative order of their magnitudes. For our first Boolean model below, we have selected the following values by considering middle-of-the-range estimates: 
and 
. Thus the loss terms in the 3-variable model in Table 2.1 are estimated to be 
, and 
. From here the times needed to reduce the concentrations of 
, and A by half are calculated to be 
, 
, and 
.
Combining this information with the methodology outlined in Section 2.5 leads to the following choice of variables for the Boolean model:
1. We will model the dynamics of the Boolean variables 
, and A, denoting mRNA, 
-galactosidase, and allolactose.
2. We assume that glucose is absent and intracellular lactose is present at all times. The intracellular lactose concentration L is a parameter for the model. From Section 2.4, we know that the model in Table 2.1 has multiple steady states when the internal lactose concentration is in a certain intermediate (maintenance) range, estimated numerically to be between (0.039, 0.055) mM (see Figure 2.4). To distinguish between low, medium, and high lactose concentrations in the Boolean case, we introduce an additional Boolean parameter 
. The value 
implies intracellular lactose concentrations are within the maintenance range or higher while 
implies intracellular lactose concentrations are below the maintenance range. When 
, internal lactose concentration is higher than the maintenance range.
3. We select a discrete time step of about 10 min for the Boolean model. The delays 
can then be ignored since they are much smaller than the time step. Similarly, since 
min, in the absence of new mRNA production, any available amounts of mRNA would be reduced below the discretization threshold in one time step. Thus there is no need for introducing 
.
4. We define additional variables 
, and 
to model the different degradation rates of allolactose and 
-galactosidase.
We emphasize that the choices regarding the size of the discrete time step and the exact number of “old” variables chosen for this model represent just one of many possibilities. Alternative models are considered later in the chapter and in the exercises.
Combining these assumptions with the assumption that translation and transcription happen in one time step, we obtain the following Boolean model:
(2.55)
A justification for the transition functions in Eqs. (2.55) follows:
Transition Equation for M: In the presence of allolactose at time t, mRNA will be produced and present at the next time step 
. Availability of M at time t does not affect its availability at time 
since M is assumed to degrade completely within one time step.
Transition Equation for B: When mRNA is available at time t, translation and transcription of 
-galactosidase will take place and 
at step 
. If amounts of B produced earlier are still available in high enough concentrations to not fall below the discretization threshold by time 
(that is, 
), B will still be available at time 
.
Transition Equation for 
: When no mRNA is available at time t, no newly produced 
-galactosidase will be available at time 
. By time 
, the amounts of B produced at time t will be reduced once due to dilution and degradation and 
.
Transition Equation for 
: When no mRNA is available at time t, no newly produced 
-galactosidase will be available at time 
. By time 
, the amounts of B, already reduced once at time t, will be further reduced due to dilution and degradation and 
.
Transition Equation for A: There are three possible ways for A to be available at time 
: (i) At time t, 
-galactosidase above the discretization threshold is available together with at least medium concentration of lactose. Under those conditions, 
-galactosidase will convert lactose into allolactose by time 
. (ii) At time t, the high concentration of intracellular lactose ensures that available trace amounts of 
-galactosidase will convert enough lactose molecules into allolactose to bring the concentration of allolactose at time 
above the discretization threshold. (iii) At time t, the concentration of A is high enough not to be reduced below the threshold level at time 
due to dilution and it will not be lost via conversion into glucose and galactose (mediated by 
-galactosidase).
Transition Equation for 
: At time t, the conditions for producing A by time 
are not met. Amounts of A available at time t will be reduced once by time 
due to dilution and degradation and 
.
We analyze the Boolean model defined by Eqs. (2.55) using the web-based DVD software [18], considering all four possible combinations for the parameter values. When 
, and 
, the cell is producing all necessary proteins for the metabolism of lactose. In this case we say that the operon is On. We say that the lac operon is Off when those proteins are not being produced (
, and 
). The results are presented in Table 2.5. Regardless of the parameter values, the system has only fixed points and no limit cycles. As expected, low concentration of lactose 
drives the system to a single steady state in which the operon is Off, while for high concentrations of lactose 
the system settles in a fixed point at which the operon is On. For intermediate lactose concentrations 
the model approximates the bistable nature of the operon, qualitatively replicating the bistability results from Section 2.4. If we start at the fixed point 1 in Table 2.5 (corresponding to low inducer concentration and an Off state for the operon) and increase the lactose concentration to medium, the operon remains Off (fixed point 3). If we start at fixed point 1 and increase the lactose concentration to high, the operon turns On (fixed point 4). In a similar way, if we start at fixed point 2 (corresponding to high inducer concentration and an On state for the operon) and reduce the lactose concentration to medium, the operon remains On (fixed point 4). Thus, in intermediate inducer concentrations the long-term behavior of the system depends on its history (hysteresis).
Table 2.5
Fixed points for the Boolean model from Eqs. (2.55). The system always settles in a fixed point. There are no limit cycles. Two steady states are possible for medium lactose concentrations (the system is bistable).
By considering the middle-of-the-range value for 
, we essentially made the assumption that allolactose degrades much slower than mRNA and that the delays 
and 
are negligible in comparison with the time for degradation and dilution of allolactose. However, considering estimates for 
among the largest reported in the literature (e.g., 
), a different Boolean model would be more appropriate since in that case the half-life for A, estimated at 
, is similar to the half-life of mRNA. This situation calls for assumptions different from those used to build the model given by Eqs. (2.55): (i) If a larger time step is considered (e.g., t = 15 min), we can eliminate the variables 
and 
. (ii) If we consider a much smaller time step (e.g., 
), more additional variables will be needed in order to account for the system delays and for the fact that multiple time steps will be needed for dilution and degradation to bring M and A below the discretization threshold levels.
The Boolean model defined by Eqs. (2.56) would be appropriate under assumption (i).
(2.56)
Exercise 2.6
(a) Justify the model presented by Eqs. (2.56) by explaining the logical expression defining each of the transition functions as was done above for the Boolean model in Eqs. (2.55); (b) Explain why the model would be consistent with choosing a time step 
min.; (c) Use DVD to analyze the model and obtain the state space diagram presented in Figure 2.9; (d) Does the result from (c) imply bistability for the system?
Figure 2.9 The state space diagram displaying the points (M,B,
, A) for the Boolean model from Eqs. (2.56). Panel A: Low level of internal lactose, corresponding to 
. Single fixed point, the operon is Off; Panel B: High level of lactose, corresponding to 
. Single fixed point, the operon is On; Panel C: Medium level of lactose, corresponding to 
. Two possible fixed points exist, the operon can be On or Off.
The Boolean model defined by Eqs. (2.57) on the other hand would be appropriate under the assumption (ii) above 
. New variables 
and 
are introduced to model the delayed effect (by 
) of mRNA on the production of 
-galactosidase and 
is needed to represent the delayed action of A on the production of mRNA by 
. As before, 
and 
track the loss of M and A due to dilution and degradation. Since the loss of B is much slower, several “old” variables are needed. We have used two such variables here but one would think that a much larger number of such variables will be needed to represent the time scales accurately. In the discussion we argue that this is unnecessary.
(2.57)
Analysis of the long-term behavior of the models from Eqs. (2.25) and (2.57) leads to results similar to those for the Boolean model from Eqs. (2.55). Since the state space for the model defined by Eqs. (2.56) is relatively small, we have presented it in Figure 2.9 (see Exercise 2.6). Regardless of the initial conditions, for low lactose 
the system settles in a state where 
, and 
and the operon is Off. For high lactose 
the system settles in a state where 
, and 
and the operon is On. For intermediate levels of lactose 
, there are two fixed points and the state of the operon depends on the initial conditions (its history): If the lactose concentration is raised to medium for cells grown in low lactose concentrations (those cells have thus settled in the Off state), the operon remains Off. If the concentration of lactose is decreased to medium for cells grown in lactose-rich environment (those cells have thus settled in the On state) the operon remains On.
For the model defined by Eqs. (2.57) the entire state space is too large to display but we present its fixed points in Table 2.6.
Table 2.6
Fixed points for the Boolean model from Eqs. (2.57). The system always settles in a fixed point. There are no limit cycles. Two steady states are possible for medium lactose concentrations (the system is bistable).
So far, we have used designated “old” variables to separate the time scales of the dilution and degradation processes from those of synthesis. As the next model shows, this could also be done implicitly through the logical assumptions built into the transition functions. If we choose to work with the estimate for 
from Wong et al. [17], the degradation time for A will be the slowest among the three variables 
, and A. The transition function for A in the next model accounts for this by allowing for 
either because of new production or because previously produced amounts of A are still present and have not been lost due to conversion to glucose and galactose.
(2.58)
2.6.2 Boolean Variants of the 5-Variable Model
We next consider Boolean variants of the model from Table 2.2 that uses five dynamic variables, 
, and P, representing mRNA, 
-galactosidase, allolactose, intracellular lactose, and lac permease. Recall that extracellular glucose is assumed to be absent and extracellular lactose 
present at all times. 
is a parameter for the model.
The degradation constants for L and P are estimated to be 
(meaning the dilution loss of L is fully due to the growth rate) and 
is an estimate for the delay 
in the equations for P[9].
As before, we need to be able to distinguish between high, medium, and low concentrations of external lactose, where medium concentration would roughly correspond to the maintenance range estimated to be (0.027, 0.062) mM (see Figure 2.6). We introduce an additional parameter 
. Assuming that 
stands for at least medium external lactose, the combination 
indicates medium external lactose, while 
and 
represent low and high external lactose, respectively.
As in the case of the 3 variable differential equation model from Table 2.1, it was shown in [9] that the delays in the equations do not affect the bistable behavior of the system. With this motivation, our first Boolean model for the 5-variable system ignores the delays. We also choose to introduce a single “old” variable for each of the model variables (see Section 2.7 for more on this). The transition functions for this Boolean model are:
(2.59)
The justification for the transition functions of the “old” variables is the same as before: no new quantities are made and previously available amounts have not already been reduced once by dilution and degradation. Also note that the justifications for the transition functions for M and B are as those for the Boolean models from the previous section.
Transition Equation for A: There are three possible ways for A to be available at time 
: (i) At time t, 
-galactosidase is available together with at least medium concentration of lactose. Under those conditions 
-galactosidase will convert lactose into allolactose by time 
. (ii) At time t, internal lactose is already present and the concentration of extracellular lactose is high; this will ensure that by time 
, the intracellular lactose concentration is high enough to find available trace amounts of 
-galactosidase and initiate enough conversions into allolactose to bring the concentration of allolactose above the discretization threshold. (iii) At time t the concentration of A is high enough not to be reduced below the threshold level at time 
due to dilution or degradation and that it will not be lost to conversion into glucose and galactose (mediated by 
-galactosidase) by time 
.
Transition Equation for L: There are three possible ways for L to be available at time 
: (i) At time t, permease above the discretization threshold is available together with at least medium levels of external lactose. Under those conditions, the permease will facilitate the transport of external lactose into the cell. (ii) At time t, the concentration of P may be below the discretization threshold but the high concentration of external lactose at time t will guarantee that the external lactose enters the cell and brings, by time 
, internal lactose above the threshold. (iii) At time t there is still high enough concentration of L that will not be lost due to the reversible nature of permease-facilitated lactose transport and also not lost to conversion into allolactose mediated by B.
Transition Equation for P: When mRNA is available at time t, translation and transcription of permease will take place to make P available at time 
. If high enough amounts of P produced earlier are still available by the time t (that is, by time t they have not already been reduced enough to be below the threshold at time 
) P will still be available at 
.
When the variables M, A, L, B, and P all have values 1, the operon is On. When they all have values 0, the operon is Off. Analyzing this model with DVD shows no limit cycles, a single fixed point for the cases of low and high concentrations of external lactose and two fixed points for medium concentrations. Regardless of the initial conditions, the operon is On for high external lactose and Off for low external lactose. When the external lactose concentration is medium, the operon can be On or Off based on the initial conditions (Table 2.7).
Exercise 2.9
Assuming a time step t = 1 min, expand the Boolean model from Eqs. (2.59) to include additional variables 
, and 
accounting, as in the model from Eqs. (2.57), for the delays in the 5-variable model from Table 2.2. Use DVD to show that introducing the delays does not affect bistability.
Exercise 2.10
Reduce the Boolean model from Eqs. (2.59) by removing all “old” variables except for 
. Explain why this would be justified. Use DVD to find the fixed points of the system and show that the qualitative bistable behavior of this model remains the same as for the Boolean model from Eqs. (2.59).
Exercise 2.11
In Santillán et al., 2007 [3], the authors consider a minimal ODE model of the lac operon involving three variables: the intracellular concentration of mRNA (M), lacZ polypeptide (E), and intracellular lactose (L). Since 
-galactosidase is a homo-tetramer made up of four identical lacZ polypeptides and the translation rate of the lacY transcription is assumed to be the same as the rate for the lacZ transcript, the following holds for the intracellular concentrations of 
-galactocidase (B) and permease (Q): 
and 
. The model also assumes that the concentrations of internal lactose (L) and allolactose (A) are the same. The three ODEs for the model are given by Eqs. (2.60). External lactose 
and external glucose 
are parameters for the model. The full model is presented and justified in [3].
(2.60)
In Chapter 1 of this volume we examined and analyzed a Boolean model built under the same assumptions and defined by the set of transition functions in Eqs. (2.61). Since this Boolean model does not have the ability to distinguish between low, medium, and high levels of lactose, it does not exhibit bistability.
(2.61)
1. Give justification for the transition functions presented by Eqs. (2.61). Use DVD to calculate the system’s fixed points and present a table with the fixed points for all possible combinations of the parameter values.
2. Now assume that in the model from Eqs. (2.61), 
stands for external concentration of lactose that is at least medium. Then introduce a new parameter 
to denote high levels of external lactose. The combination of parameter values 
and 
now stands for medium external lactose. Modify the transition functions from Eqs. (2.61) to make the modified model exhibit bistable behavior for medium lactose concentrations.
Table 2.7
Fixed points for the Boolean model from Eqs. (2.59). The system always settles in a fixed point. There are no limit cycles. Two steady states are possible for medium lactose concentrations (the system is bistable).
2.7 Conclusions and Discussion
In this chapter we developed and compared a number of mathematical models of the lactose regulatory mechanism of E. coli. The focus was on bistability. The continuous models were previously published in [8,9], while the Boolean models we have considered are new. Since the discovery of the lac operon by Jacob and Monod in 1960 ([19,20]), both differential equation models and Boolean models have been introduced to describe operon system dynamics, beginning with the work of Goodwin [21] and Kauffman [22], respectively. Bistability for the lac operon system has been experimentally observed and simulated by a number of continuous models, including those considered in this chapter. However, establishing that Boolean models are capable of capturing the system’s bistability has only been done recently [16,23]. A necessary condition for bistability when modeling with Boolean networks is to make possible distinguishing between at least three levels of inducer concentrations: low, medium, and high. Both deterministic models (e.g., [16]) and models involving stochasticity (e.g., [23]) have been proposed.
In this chapter we confirmed that Boolean models can approximate delayed differential equation models, preserving critical qualitative features such as bistability. It was noted in [8] that the model in [9] does not consider inducer exclusion or catabolite repression (both are external-glucose-dependent mechanisms), indicating that bistability is independent from the presence of glucose in the extracellular medium. The 3-variable model from Table 2.1 in [8] ignores the lactose permease in the operon regulation and considers only the role of 
-galactosidase. This shows that of all feedback loops in the system, the 
-galactosidase regulation is the one responsible for the bistable behavior of the operon. The fact that the Boolean models that approximate this continuous model preserve bistability further underscores this result. In [8,9] the authors also show that bistability of the system is preserved when the delays in the regulatory equations from Tables 2.1 and 2.2 are ignored. The Boolean models replicate this finding, indicating once again that bistability is a robust phenomenon that arises from the network structure.
The Santillán model from Exercise 2.11[3], considers external lactose as a parameter of the model and the function 
(an increasing function of L and a decreasing function of 
) in the equation for M in Eqs. (2.60) is a possible way to model the dependence of the system on these parameters. In [24] the authors provide a detailed comparison of a number of ODE models of the lac operon that differ in the numbers of variables and type (stochastic or deterministic), including the three ODE models considered in this chapter. Many of those models differ in the ways the dependence of mRNA on external glucose and internal lactose is modeled. It is noted in [24] that in some cases heuristic reasoning is used to propose Hill-type equations for the function 
with a different level of detail. Establishing that Boolean approximations of the minimal model from Eqs. (2.60) can describe and explain bistability indicates that such differences are nonessential with regard to bistability and do not impact this feature of the system.
It has been shown that bistability requires a direct positive feedback loop or an indirect positive feedback loop, such as a double negative feedback loop [25]. However, a feedback loop alone is not enough for a system to exhibit bistable behavior. It must also possess some type of nonlinearity within the feedback circuit. That is, some of the proteins in the feedback circuit must respond to their upstream regulators in an ultrasensitive manner [26–28]. The Hill coefficient is used to quantify the steepness of this response. A Hill coefficient larger than one is considered to be the ultrasensitive response [29]. In Boolean models this condition is automatically satisfied since the On-Off switches in such models can be viewed as responses with very large Hill coefficients.
The analysis of the Boolean models was done using the web-based software DVD [18]. When the number of variables is small, DVD can be used to produce the entire state space for the model. When the number of variables grows, obtaining the entire state space is not feasible, but DVD computes and returns the model’s fixed points together with the number of connected components in state space. DVD’s successor ADAM [30], which handles a broader class of discrete models, including logical models and Petri nets in addition to Boolean and polynomial models, can also be used.
We found that introducing multiple “old” variables into the Boolean models does not appear to change the long-term behavior of the system and to impact bistability. This is consistent with findings reported in [31] where the author presents a method for reducing Boolean networks and their wiring diagrams while preserving the set of fixed points. The reduction is done by deleting vertices with no self-loops (that is, vertices whose transition functions do not include them as inputs). If a vertex is to be removed from the network, its transition function is substituted for the variable representing this vertex in all of the other transition functions in the model. The idea is to reduce the size of the network by providing a way to delete a vertex and “pass on” its functionality to other variables in the network. A reduced model with the same fixed points as the original model would be indicative of characteristics that emerge not as a result of some specific system interactions but, instead, from the core network topology. For instance, in [23] Veliz-Cuba and Stigler present examples of Boolean models of the lacoperon incorporating inducer inclusion and catabolite repression. The networks for those models are then reduced by this method to smaller models that no longer include those regulatory mechanisms but preserve the bistable behavior of the system.
In the Boolean models proposed in this chapter we use “old” variables to track the dilution and degradation of concentrations that require multiple steps of time for reduction below the discretization threshold. In the context of the notation used in Section 2.5, if a variable Y regulates the production of X, which (if no new amounts are produced) would degrade below threshold levels in n steps, the Boolean model will include the following motif:
(2.62)
We can reduce the network by eliminating 
. To do this, substitute its transition function in place of 
in the transition function of X: 
, which simplifies to 
, resulting in exactly the same structural motif with one less “old” variable. Variables 
can be eliminated similarly, leading to the reduction 
. When applied to the Boolean models in this chapter, this reduction process indicates that the bistability property of the system does not depend on the number of “old” variables used in the model. We finally note that if we choose to also eliminate 
from the model, this leads to the transition equations 
, reflecting a situation in which X is considered completely stable (a situation in which the combined degradation and dilution rate for X is infinitely small and thus, the growth rate is negligible), which is biologically unrealistic. Thus, at least one “old” variable is needed.
Exercise 2.12
Confirm that 
.
The use and treatment of “old” variables here differ from the approach introduced in [16]. We stipulate that an “old” variable has value 1 at time 
only when conditions for new production are not met at time tand when previously produced amounts available at time t have not already been reduced by a certain factor due to dilution and degradation. In [16], an “old” variable has value 1 at time 
when conditions for new production are not met at time t. Our approach provides a mechanism to track the level of reduction of X: since only one 
could be equal to 1 at each time step, 
means that at time t the concentration of X has already been reduced exactly k times 
.
Thus, when delay variables are added to our Boolean model from Eqs. (2.59) (see Exercise 2.9) the resulting model provides a Boolean approximation of the 5-variable differential equation model from Table 2.57 that differs from the Boolean approximation of the same model presented in [16]. Both models capture the bistable behavior of the lactose operon but the model in [16] generates several fixed points that are not biologically meaningful. Another difference in comparison with the Boolean model in [16] is that the Boolean models in this chapter factor in the bacterial growth rate 
.
The Boolean modeling carried out here is strictly deterministic in nature and uses synchronous updates for all variables. It would be interesting to examine what impact stochasticity and asynchronous update schedules may have on the model outcomes.
Acknowledgment
Raina Robeva gratefully acknowledges the support of NSF under the DUE award 0737467.
2.8 Supplementary Materials
All supplementary files and/or computer code associated with this article can be found from the volume’s website http://booksite.elsevier.com/9780124157804
References
1. Novick A, Wiener M. Enzyme induction as an all-or-none phenomenon. Proc Natl Acad Sci USA. 1957;43:553–566.
2. Ozbudak EM, Thattai M, Lim HN, Shraiman BI, Oudenaarden AV. Multistability in the lactose utilization network of Escherichia coli. Nature. 2004;427:737–740.
3. Santillán M, Mackey MC, Zeron E. Origin of bistability in the lac operon. Biophys J. 2007;92:3830–3842.
4. Beckwith J. The lactose operon. In: Washington, DC: American Society for Microbiology; 1987;1444–1452. Neidhardt FC, Ingraham JL, Low KB, Magasanik B, Umbarger HE, eds. Escherichia coli and Salmonella: Cellular and molecular biology. vol. 2.
5. Schleif R. Regulation of the L-arabinose operon of Escherichia coli. Trends in Genetics. 2000;16(12):559–565.
6. Lewin Benjamin. Genes. 9th ed. Jones and Bartlett Publishers 2008.
7. Klipp E, Herwig R, Kowald A, Wierling C, Lehrach H. Systems biology in practice. 1st ed. Wiley-VCH 2005.
8. Yildirim N, Santillán M, Horike D, Mackey MC. Dynamics and bistability in a reduced model of the lac operon. Chaos. 2004;14(2):279–292.
9. Yildirim N, Mackey MC. Feedback Regulation in the Lactose Operon: A Mathematical Modeling Study and Comparison with Experimental Data. Biophysical J. 2003;84:2841–2851.
10. Watson JD. Molecular Biology of the Gene. 3rd ed WA, New York: Benjamin Inc; 1977.
11. Yagil G, Yagil E. On the relation between effector concentration and the rate of induced enzyme synthesis. Biophys J. 1971;11:11–27.
12. Saier MH, Ramseier TM, Reizer J. Regulation of carbon utilization. In: Washington, DC: American Society for Microbiology; 1996;1325–1343. Neidhardt FC, Curtiss R, Ingraham JL, eds. Escherichia coli and Salmonella: cellular and molecular biology. vol. 1.
13. Mackey MC, Santillán M, Yildirim N. Modeling operon dynamics: the tryptophan and lactose operons as paradigms. C R Biol. 2004;327:211–224.
14. Cohn M, Horibata K. Inhibition by glucose of the induced synthesis of the 
-galactosidase-enzyme system of Escherichia coli: Analysis of maintenance. J Bacteriol. 1959;78:613–623.
15. Yildirim N, Kazanci C. Deterministic and stochastic simulation and analysis of biochemical reaction networks the lactose operon example. Methods Enzymol. 2011;487:371–395.
16. Hinkelmann F, Laubenbacher R. Boolean Models of Bistable Biological Systems. Discrete and Continuous Dynamical Systems. 2011;4:1443–1456.
17. Wong P, Gladney S, Keasling JD. Mathematical model of the lac operon: Inducer exclusion, catabolite repression, and diauxic growth on glucose and lactose. Biotechnol Prog. 1997;13:132–143.
18. Vastani H, Jarrah A, Laubenbacher. Visualization of Dynamics for Biological Networks. http://dvd.vbi.vt.edu/dvd.pdf.
19. Jacob F, Perrin D, Sanchez C, Monod J. L’Operon: groupe de gène à expression par un operatour. C R Seances Acad Sci. 1960;250:1727–1729.
20. Jacob F, Monod J. Genetic regulatory mechanisms in the synthesis of proteins. J Mol Biol. 1961;3:318–356.
21. Goodwin B. Oscillatory behaviour in enzymatic control process. Adv Enz Regul. 1965;3:425–438.
22. Kauffman S. Metabolic stability and epigenetics in randomly constructed gene nets. J Theor Biol. 1969;22:437–467.
23. Veliz-Cuba A, Stigler B. Boolean models can explain bistability in the lac operon. J Comput Biol. 2011;18:783–794.
24. Santillán M, Mackey M. Quantitative approaches to the study of bistability in the lac operon of Escherichia coli. J R Soc Interface. 2008;5:S29–S39.
25. Laurent M, Kellershohn N. Multistability: a major means of differentiation and evolution in biological systems. Trends Biochem Sci. 1999;24:418–422.
26. Koshland DE, Goldbeter A, Stock JB. Amplification and adaptation in regulatory and sensory systems. Science (New York, N.Y.). 1982;217:220–225.
27. Ferrell JE. Self-perpetuating states in signal transduction: positive feedback, double-negative feedback and bistability. Curr Opin Cell Biol. 2002;14:140–148.
28. Ferrell JE. Tripping the switch fantastic: how a protein kinase cascade can convert graded inputs into switch-like outputs. Trends Biochem Sci. 1996;21:460–466.
29. Ferrell JE. Building a cellular switch: more lessons from a good egg. Bioessays. 1999;21:866–870.
30. Hinkelmann F, Brandon M, Guang B, et al. ADAM: Analysis of discrete models of biological systems using computer algebra. BMC Bioinformatics. 2011;12:437–467.
31. Veliz-Cuba A. Reduction of Boolean network models. J Theor Biol. 2011;289:167–172.