You are viewing a javascript disabled version of the site. Please enable Javascript for this site to function properly.
Go to headerGo to navigationGo to searchGo to contentsGo to footer
In content section. Select this link to jump to navigation

The utility of simple mathematical models in understanding gene regulatory dynamics

Abstract

In this review, we survey work that has been carried out in the attempts of biomathematicians to understand the dynamic behaviour of simple bacterial operons starting with the initial work of the 1960’s. We concentrate on the simplest of situations, discussing both repressible and inducible systems and then turning to concrete examples related to the biology of the lactose and tryptophan operons. We conclude with a brief discussion of the role of both extrinsic noise and so-called intrinsic noise in the form of translational and/or transcriptional bursting.

1Introduction

The operon concept for the regulation of bacterial genes, first put forward in [1], has had an astonishing and revolutionary effect on the development of understanding in molecular biology. It is a testimony to the strength of the theoretical and mathematical biology community that modeling efforts aimed at clarifying the implications of the operon concept appeared so rapidly after the concept was embraced by biologists. Thus, to the best of our knowledge the first analysis of operon dynamics appeared in [2] and in [3]. These first attempts were swiftly followed by Griffith’s analysis of a simple repressible operon [4] and an inducible operon [5], and these and other results were beautifully summarized in [6].

Since these modeling efforts in the early days of development in molecular biology, both our biological knowledge and level of sophistication in modeling have proceeded apace to the point where new knowledge of the biology is actually driving the development of new mathematics. This is an extremely exciting development and one which many have expected–that biology would act as a driver for mathematics in the 21st century much as physics was the driver for mathematics in the 19th and 20th centuries. However, as this explosion of biological knowledge has proceeded hand in hand with the development of mathematical modeling efforts to understand and explain it, the difficulty in comprehending the nature of the field becomes ever deeper due to the sheer volume of work being published.

In this highly idiosyncratic review we discuss work from our group over the past few years directed at the understanding of really simple operon control dynamics. We start this review in Section 2 by discussing transcription and translation kinetics (Section 2.1) and then pass to general dynamics considerations in Section 2.2 which is largely a recap of earlier work with additional insights derived from the field of nonlinear dynamics. We then pass to the role of transcriptional and translational delays in Section 2.3 and finish with a short consideration of fast and slow variables in Section 2.4. Following this, in Sections 3.1 and 3.2 we pass from the realm of mathematical nicety to biological reality by looking at realistic models for the lactose and tryptophan operons respectively. These two examples, two of the most extensively experimentally studied systems in molecular biology, and for which we have vast amounts of data, illustrate the reality of dealing with experimental biology and the difficulties of applying realistic modeling efforts to understand that biology.

Finally, in Section 4 we turn to one of the more interesting and challenging aspects of understanding operon dynamics. In the last few years with the advent of ever improved imaging techniques combined with rapid data acquisition techniques experimentalists have acquired the ability to peer ever more closely into the details of these dynamics at virtually the single molecule level. This means, therefore, that all manner of interesting statistical behaviours are being uncovered–behaviours that reveal many interesting types of ‘random’ behaviour not well understood from a mathematical perspective. We explore aspects of this in Section 4.1 where we consider the effects of transcriptional and/or translational bursting, and in Section 4.2 where we look at the effects of fluctuations in degradation rates. The review ends with a brief discussion in Section 5.

2Generic deterministic models of prokaryoticgene regulation

The central tenet of molecular biology was put forward some half century ago, and though modified in detail still stands in its basic form. Transcription of DNA produces messenger RNA (mRNA, denoted M here). Then through the process of translation of mRNA, intermediate protein (I) is produced which is capable of controlling metabolite (E) levels that in turn can feedback and affect transcription and/or translation. A typical example would be in the lactose operon of Section 3.1 where the intermediate is β-galactosidase and the metabolite is allolactose. These metabolites are often referred to as effectors, and their effects can, in the simplest case, be either stimulatory (so called inducible) or inhibitory (or repressible) to the entire process. This scheme is often called the ‘operonconcept’.

2.1Kinetic aspects of regulation of transcription and translation

We first outline the relatively simple molecular dynamics of both inducible and repressible operons and how effector concentrations can modify transcription rates. If transcription rates are constant and unaffected by any effector, then this is called a ‘no control’ situation.

2.1.1Inducible regulation

The lac operon considered below in Section 3.1 is the paradigmatic example of inducible regulation. In an inducible operon when the effector (E) is present then the repressor (R) is inactive and unable to bind to the operator (O) region so DNA transcription can proceed unhindered. E binds to the active form R of the repressor and we assume that this binding reaction is

R+nEk1k1+REn,
in which k1+ and k1- are the forward and backward reaction rate constant, respectively. The equilibrium equation for the reaction above is
(2.1)
K1=REnR·En,
where K1=k1-/k1+ is the reaction dissociation constant and n is the number of effector molecules required to inactivate repressor R. The operator O and repressor R are also assumed to interact according to
O+Rk2k2+OR,
which has the following equilibrium equation:
K2=ORO·R,K2=k2-k2+.
The total operator O tot is given by
Otot=O+OR=O+K2O·R=O(1+K2R),
while the total repressor is R tot
Rtot=R+K1R·En+K2O·R.
Furthermore, by definition the fraction of operators free to synthesize mRNA (i.e., not bound byrepressor) is
f(E)=OOtot=11+K2R.
If the amount of repressor R bound to the operator O is small
RtotR+K1R·En=R(1+K1En)
so
R=Rtot1+K1En,
and consequently
(2.2)
f(E)=1+K1En1+K2Rtot+K1En=1+K1EnK+K1En,
where K = 1 + K 2 R tot . Maximal repression occurs when E = 0 and even at that point mRNA is produced (so-called leakage) at a basal level proportionalto K -1.

Assume that the maximal transcription rate of DNA (in units of time - 1) is φ¯m. Assume further that transcription rate φ in the entire population is proportional to the fraction of unbound operators f. Thus we expect that φ as a function of the effector level will be given by φ=φ¯mf, or

(2.3)
φ(E)=φ¯m1+K1EnK+K1En.

2.1.2Repressible regulation

The tryptophan operon considered below in Section 3.2 is the classic example of a repressible system. This is because the repressor is active (capable of binding to the operator) when the effector molecules are present which means that DNA transcription is blocked. Using the same notation as before, but realizing that the effector binds the inactive form R of the repressor so it becomes active and take this reaction to be the same as in Equation 2.1. However, we now assume that the operator O and repressor R interaction is governed by

O+REnk2k2+OREn,
with the following equilibrium equation
(2.4)
K2=OREnO·REn,K2=k2-k2+.
The total operator is
Otot=O+OREn=O+K1K2O·R·En=O(1+K1K2R·En),
so the fraction of operators not bound by repressor is
f(E)=OOtot=11+K1K2R·En.
Assuming, as before, that the amount of R bound to O is small compared to the amount or repressor gives
f(E)=1+K1En1+(K1+K1K2Rtot)En=1+K1En1+KEn,
where K = K 1 (1 + K 2 R tot ). In this case we have maximal repression when E is large, and even when repression is maximal there is still be a basal level of mRNA production (again known as leakage) which is proportional to K 1 K -1 < 1. Variation of the DNA transcription rate with effector level is given by φ=φ¯mf or
(2.5)
φ(E)=φ¯m1+K1En1+KEn.

Both (2.3) and (2.5) are special cases of

(2.6)
φ(E)=φ¯m1+K1EnA+BEn=φ¯mf(E).
The constants A, B ≥ 0 are defined in Table 2.1.

2.2General dynamic considerations

The Goodwin model for operon dynamics [2] considers a large population of cells, each of which contains one copy of a particular operon, and we use that as a basis for discussion. We let (M, I, E) respectively denote the mRNA, intermediate protein, and effector concentrations. For a generic operon with a maximal level of transcription b¯d (in concentration units), the dynamics are given by [2, 4, 5, 7, 8]

(2.7)
dMdt=b¯dφ¯mf(E)-γMM,
(2.8)
dIdt=βIM-γII,
(2.9)
dEdt=βEI-γEE.
It is assumed here that the rate of mRNA production is proportional to the fraction of time the operator region is active, and that the rates of protein and metabolite production are proportional to the amount of mRNA and intermediate protein respectively. All three of the components (M, I, E) are subject to random degradation, and the function f is as determined in Section 2.1 above.

To simplify things we formulate Equations (2.7)–(2.9) using dimensionless concentrations. To start we rewrite Equation (2.6) in the form

φ(e)=φmf(e),
where φ m (which is dimensionless) is defined by
(2.10)
φm=φ¯mβEβIγMγEγIandf(e)=1+enΛ+Δen,
Λ and Δ are defined in Table 2.1, and a (dimensionless) effector concentration (e) is defined by
E=ηewithη=1K1n.
We continue and define dimensionless intermediate protein (i) and mRNA concentrations (m):
I=iηγEβEandM=mηγEγIβEβI,
so Equations (2.7)–(2.9) take the form
dmdt=γM[κdf(e)-m],didt=γI(m-i),dedt=γE(i-e),
with the dimensionless constants
κd=bdφmandbd=b¯dη.

To continue our simplifications we rename the dimensionless concentrations through (m, i, e) = (x 1, x 2, x 3), and subscripts (M, I, E) = (1, 2, 3) to finally obtain

(2.11)
dx1dt=γ1[κdf(x3)-x1],
(2.12)
dx2dt=γ2(x1-x2),
(2.13)
dx3dt=γ3(x2-x3).
In all of these equations, γ i for i = 1, 2, 3 denotes a degradation rate (units of inverse time), and thus Equations (2.11)–(2.13) are not in dimensionless form. The dynamics of this classic operon model have been fully analyzed [9], the results of which we simply summarize here. We set X = (x 1, x 2, x 3) and let S t  (X) be the flow generated by the system (2.11)– (2.13), i.e., the function t ↦ S t  (X) is a solution of (2.11)– (2.13) such that S 0 (X) = X. For both inducible and repressible operons, for all initial conditions X0=(x10,x20,x30)3+ the flow St(X0)3+ for t > 0.

The steady state solutions of (2.11)– (2.13) are given by the solutions of

(2.14)
xκd=f(x)
and for each solution x * of Equation (2.14) there is a steady state X*=(x1*,x2*,x3*) of (2.11)– (2.13) which is given by
x1*=x2*=x3*=x*.
Whether there is a single steady state X * or there are multiple steady states will depend on whether we are considering a repressible or inducible operon.

2.2.1No control

No control simply means f (x) ≡1, and in this case there is a single steady state x * = κ d that is globally asymptotically stable.

2.2.2Inducible regulation

Single versus multiple steady states. For an inducible operon (with f given by Equation (2.2)) there may be one (X1* or X3*), two (X1*,X2*=X3* or X1*=X2*,X3*), or three (X1*,X2*,X3*) steady states, with the ordering 0<X1*X2*X3*, corresponding to the possible solutions of Equation (2.14) (cf. Figure 2.1). The smallest steady state (X1*) is typically called the un-induced state, while the largest steady state (X3*) corresponds to the induced state. The steady state values of x are easily obtained from (2.14) for given parameter values, and the dependence on κ d for n = 4 and a variety of values of K is shown in Fig. 2.1. Figure 2.2 shows a graph of the steady states x * versus κ d for various values of the leakage parameter K.

Analytic conditions for the existence of one or more steady states come from Equation (2.14) in conjunction with the observation that the delineation points are marked by the values of κ d at which x/κ d is tangent to f (x) (see Figure 2.1). Differentiation of (2.14) yields a second condition

(2.15)
1κdn(K-1)=xn-1(K+xn)2.
From Equations (2.14) and (2.15) the values of x at which tangency will occur are given by:

(2.16)
x±=K-12{[n-K+1K-1]±n2-2nK+1K-1+1}n.
The corresponding values of κ d at which a tangency occurs are given by
(2.17)
κd±=xK+xn1+xn.

A necessary condition for the existence of two or more steady states is obtained by requiring that the radical in (2.16) is non-negative:

(2.18)
K(n+1n-1)2.
Thus a second necessary condition follows:
(2.19)
κdn+1n-1n+1n-1n.
Further, from Equations (2.14) and (2.15) we can find the boundaries in (K, κ d ) space in which there are one or three locally stable steady states as shown in Fig. 2.3. There, we have given a parametric plot (x is the parameter) of κ d versus K, using
K(x)=xn[xn+(n+1)](n-1)xn-1andκd(x)=[K(x)+xn]2nxn-1[K(x)-1],
for n = 4 obtained from Equations (2.14) and (2.15). As is clear from the figure, when leakage is appreciable (small K, e.g for n = 4, K < (5/3) 2) then the possibility of bistable behaviour is lost.

We can make some general comments on the influence of n, K, and κ d on the appearance of bistability from this analysis. First, the degree of cooperativity (n) in the binding of effector to the repressor plays a significant role and n > 1 is a necessary condition for bistability. If n > 1 then a second necessary condition for bistability is that K satisfies Equation (2.18) so the fractional leakage (K -1) is sufficiently small. Furthermore, κ d must satisfy Equation (2.19) which is quite interesting. For n→ ∞ the limiting lower limit is κ d  > 1 while for n → 1 the minimal value of κ d becomes quite large. This simply tells us that the ratio of the product of the production rates to the product of the degradation rates must always be greater than 1 for bistability to occur, and the lower the degree of cooperativity (n) the larger the ratio must be.

If n, K and κ d satisfy these necessary conditions then bistability is only possible if κ d  ∈ [κ d-, κ d+] (c.f. Figure 2.3). The locations of the minimal (x -) and maximal (x +) values of x bounding the bistable region are independent of κ d . And, finally, (x + - x -) is a decreasing function of increasing n for constant κ d , K while (x + - x -) is an increasing function of increasing K for constant n, κ d .

Local and global stability. Although the local stability analysis of the inducible operon is possible [9], the thing that is interesting is that the global stability is possible to determine.

Theorem 2.1. [7, 10, Proposition 2.1, Chapter 4] For an inducible operon with φ given by Equation (2.3), define II = [1/K, 1]. There is an attracting box BI3+ defined by

BI={(x1,x2,x3):xiII,i=1,2,3}

such that the flow St is directed inward everywhere on the surface of BI. Furthermore, all X* ∈ BI and

  • 1. If there is a single steady state, i.e. X1* for κd ∈ [0, κd-), or X3* for κd+ < κd, then it is globally stable.

  • 2. If there are two locally stable nodes, i.e. X1* and X3* for κd ∈ (κd-, κd+), then all flows St (X0) are attracted to one of them. (See [8] for a delineation of the basin of attraction of X1* and X3*.)

2.2.3Repressible regulation

As is clear from a simple consideration of our dynamical equations the repressible operon has a single steady state corresponding to the unique solution x * of Equation (2.14). Again, rather remarkably, we can characterize the global stability of this single steady state through the following.

Theorem 2.2. [10, Theorems 4.1 & 4.2, Chapter 3] For a repressible operon with φ given by Equation (2.5), define I R  = [K 1/K, 1]. There is a globally attracting box BR3+ defined by

BR={(x1,x2,x3):xiIR,i=1,2,3}
such that the flow S t is directed inward everywhere on the surface of B R . Furthermore there is a single steady state X * ∈ B R . If X * is locally stable it is globally stable, but if X * is unstable then a generalization of the Poincare-Bendixson theorem [10, Chapter 3] implies the existence of a globally stable limit cycle in B R .

2.3The appearance of cell growth effects and delays due to transcription and translation

The considerations of the previous sections must, however, be tempered by the realization that sometimes cell growth has to be taken into account as well as the fact that significant delays may enter into the dynamical equations [11]. The effects of growth are obvious in that if a cell increases its volume then there is an effect on concentrations. But where do these delays come from? Their origin is simple to understand and arises from the fact that the transcription and translation processes take place at a finite velocity and therefore require a non-zero time for completion. The existence of these delays has been known for some time by modelers [12] and whether the incorporation of the delays will potentially change the qualitative nature of the model dynamics will depend on the type of regulation. Generally when the regulation is that of an inducible operon there will be no change, but if the system is a repressible one then the inclusion of the transcriptional and translational delays may lead to the prediction of limit cycle behaviour.

Once we take growth and these transcriptional and translational delays into account, our basic dynamical equations are modified to the form

(2.20)
dMdt=b¯dφ¯mf(e-μτMEτM)-γ¯MM,
(2.21)
dIdt=βIe-μτIMτI-γ¯II,
(2.22)
dEdt=βEI-γ¯EE.
In Equations (2.20)–(2.22) there are several changes to be noted. The first is the appearance of the terms e -μτ M and e -μτ I which account for an effective dilution of the mRNA (M) and intermediate protein (I) because the cell is growing at a rate μ (time - 1). The second is the alteration of the decay rates γ i to γ¯iγi+μ because the cell growth leads to an effective increase in the rate of destruction. The last is the altered notation E τ M  (t) = E (t - τ M ) and M τ I  (t) = M (t - τ I ) indicating that both E and M are now to be evaluated at a time in the past due to the non-zero times required for transcription and translation. From a dynamic point of view, the presence of these delays can have a dramatic effect.

Equations (2.20)–(2.22) can be put in a simpler form, just as we did for (2.7)– (2.9), but by now setting

E=ηewithη=1e-μτMK1n,I=iηγEβEandM=mηγEγIβEβIe-μτI
so Equations (2.20)–(2.22) take the form
dmdt=γ¯M[κdf(eτM)-m],didt=γ¯I(mτI-i),dedt=γ¯E(i-e),
with
(2.23)
κd=b¯dφmβIβEγ¯Mγ¯Iγ¯Eηe-μτI.

To finish our simplifications, as before rename the dimensionless concentrations (m, i, e) = (x 1, x 2, x 3), and subscripts (M, I, E) = (1, 2, 3) to obtain

(2.24)
dx1dt=γ¯1[κdf(x3,τ1)-x1],
(2.25)
dx2dt=γ¯2(x1,τ2-x2),
(2.26)
dx3dt=γ¯3(x2-x3).
Again Equations (2.24)–(2.296) are not in dimensionless form.

It is important to realize that the appearance of the delays τ M and τ I (or τ 1 and τ 2) plays absolutely no role in the determination of the steady state(s) of inducible and repressible systems as discussed above.

For an inducible operon in which f′ (X *) >0 a simple extension of the proof in [10, Proposition 6.1, Chapter 6] shows that the global stability properties are not altered by the presence of the delays (τ 1, τ 2).

However, for a repressible operon there are, at this point in time, no extensions of the global stability results of [10, Theorem 4.1 & Theorem 4.2, Chapter 3] for inducible systems. The best that we can do is to linearize Equations (2.24)–(2.26) in the neighborhood of the unique steady state X * to obtain the eigenvalue equation g (λ) = P (λ) + ϑe -λτ  = 0 wherein

(2.27)
P(λ)=(γ¯1+λ)(γ¯2+λ)(γ¯3+λ)andϑ=-κdf(X*)γ¯1γ¯2γ¯3>0
and τ = τ 1 + τ 2. Writing out g (λ) we have
(2.28)
g(λ)=λ3+a1λ2+a2λ+a3+ϑe-λτ,
where
a1=i=13γi,a2=ij=13γiγj,a3=i=13γi.

Let λ (τ) = α (τ) +  (τ) be the root of Equation (2.28) satisfying α (τ 0) =0 and ω (τ 0) = ω 0, and set p=a12-2a2, q=a22-2a1a2, r=a32-ϑ2, and let h (z) = z 3 + pz 2 + qz + r. [13, Theorem 2.4] gives the conditions for X * to be locally stable and for the existence of a Hopf bifurcation.

Theorem 2.3. [13, Theorem 2.4]

  • 1. If r ≥ 0 and Δ = p 2 - 3q < 0, then all roots of Equation (2.28) have negative real parts for all τ ≥ 0.

  • 2. If r < 0 or r ≥ 0, z 1 > 0 and h (z 1) <0, then all roots of Equation (2.28) have negative real parts when τ ∈ [0, τ 0).

  • 3. If the conditions of (2) are satisfied, τ = τ 0 and h(ω02)0, then ± 0 is a pair of simple purely imaginary roots of Equation (2.28) and all other roots have negative real parts. Moreover,

    dReλ(τ0)dτ>0.

2.4Fast and slow variables

Identifying fast and slow variables can give considerable simplification and insight into the long term behaviour of the system. A fast variable in a given dynamical system relaxes much more rapidly to an equilibrium than a slow one [14]. Differences in degradation rates in chemical and biochemical systems lead to the distinction that the slowest variable is the one that has the smallest degradation rate.

Typically the degradation rate of mRNA is much greater than the corresponding degradation rates for both the intermediate protein and the effector (γ 1 ⪢ γ 2, γ 3) so in this case the mRNA dynamics are fast and we have the approximate relationship

x1κdf(x3).
If γ 1 ⪢ γ 2 ⪢ γ 3 so that the effector is the slowest variable, then we have
x2x3
and the three variable system (2.11)– (2.13) describing the generic operon reduces to a one dimensional system
(2.29)
dx3dt=γ3[κdf(x3)-x3]
for the relatively slow effector dynamics. If instead the effector qualifies as a fast variable (as for the lac operon) so that γ 1 ⪢ γ 3 ⪢ γ 2 and
x3x2
then the intermediate protein is the slowest variable described by the one-dimensional equation
(2.30)
dx2dt=γ2[κdf(x2)-x2].
Consequently both Equations (2.30) and (2.29) are of the form
(2.31)
dxdt=γ[κdf(x)-x],
where γ is either γ 2 for protein (x 2) dominated dynamics or γ 3 for effector (x 3) dominated dynamics.

Eliminating fast variables, also known as the adiabatic elimination technique [14], has been extended to stochastically perturbed systems when the perturbation is a Gaussian distributed white noise, c.f. [15, Section 11.1], [16, 17], and [18, Section 6.4]. For the case of perturbation being a jump Markov process we refer to [19].

3Specific examples in various systems

3.1The lactose operon

Glucose is the favourite carbon and energy source for E. coli, as well as for many other organisms. Although this bacterium can also feed on other sugars, it only does so when glucose is absent. A typical population of E. coli doubles its size approximately every hour in presence of a pure sugar, like glucose or lactose. The existence of the lactose operon was conjectured by Jacob and Monod after observing that a population of E. coli is initially unable to digest lactose, when it is fed with a mixture of the glucose and lactose sugars.

Monod [20] observed in his PhD work that in the presence of a mixture of glucose and lactose the exponential growth begins as usual, then it pauses for about one hour before resuming at a similar pace. The bacterial growth curve shows two distinctive phases, as can be seen in Fig. 3.1. The key observation was that the timing of the pauses was controlled by the ratio of the initial amounts of glucose and lactose: the larger initial amount of glucose the later the pause would begin. Monod realized that E. coli is initially unable to digest lactose, so that the bacteria initially feeds exclusively on glucose, until it is totally consumed and the bacteria then needs to change its internal metabolism to consume lactose. It is worth mentioning at this point that diauxic growth only occurs in batch cultures, and simultaneous usage of sugars is often observed in continuous cultures [21].

Jacob and co-workers [1] proposed the lactose operon model as a mechanism for explaining these features. Thus, the lac genes that encode the enzymes necessary for lactose absorption and hydrolysis are all controlled by a single mechanism, and they are all turned off in the presence of glucose or the absence of lactose. Properly speaking, the lactose operon is a DNA segment composed of a promoter/operator region, followed by the structural genes lacZ, lacY, and lacA, and finally by the corresponding terminator. The promoter/operator is the DNA region where the transcription factors (RNA polymerase, lactose repressor, cyclic-AMP receptor protein, et cetera) bind in order to initiate the transcription of a corresponding mRNA strand or to regulate the corresponding transcription process. The gene lacZ codes for the enzyme β-galactosidase (LacZ) that in E. coli cleaves the disaccharide lactose into glucose and galactose. The gene lacY codes for the enzyme β-galactoside permease (LacY), an inner membrane-bound symporter that pumps lactose into the cell using a proton gradient. Finally, lacA encodes the enzyme β-galactoside transacetylase (LacA) that transfers an acetyl group from acetyl-sides. Nevertheless, it is not completely understood what its precise function is.

The β-galactosidase enzyme. Few genes have a history of study as long and distinguished as lacZ. The lacZ gene encodes an open reading frame of 1024 amino acids and was one of the first large genes to be completely sequenced. In E. coli, the biologically active β-galactosidase protein exists as a tetramer of four identical subunits and has a molecular weight of approximately 480– 500 kDa. The primary enzymatic function of β-galactosidase relevant to its role as a biotechnological tool is to cleave the chemical bond between the anomeric carbon and glycosyl oxygen of appropriate substrates; see for example [22].

lacZ was chosen as the target of a very extensive early analysis, in part owing to specific experimental advantages accompanying work with β-galactosidase. These advantages continue to provide a rationale for using this protein in biotechnological applications today.

The β-galactoside permease protein. Active transporters (pumps) require a cellular energy source (i.e. ATP hydrolysis) to catalyze the transport of charged components against an electrochemical gradient. Depending on their energy source, active transporters are classified as primary or secondary. In particular, secondary transporters use the free energy stored in a given electrochemical ion gradient, as shown in [23]. β-galactoside permease is a secondary transporter that couples free energy released from downhill translocation of protons to drive the uphill translocation of galactosides against a concentration gradient. This protein is composed of 417 amino acid residues and has 12 helices that transverse the membrane in a zigzag fashion, connected by relatively hydrophilic loops with both N and C termini on the cytoplasm side. β-galactoside permease is encoded by the lacY gene, the second structural gene in the lactose operon. lacY was the first gene encoding a membrane transport protein to be cloned into a recombinant plasmid, over-expressed and sequenced; see for example [24] and the references therein. This success in the early days of molecular biology opened the study of secondary active transport at the molecular level. Thus, β-galactoside permease was the first protein of its class to be solubilized and purified in a completely functional state, thereby demonstrating that this single gene product is solely responsible for all the translocation reactions catalyzed by the galactoside transport system in E. coli. [24] has also shown that this protein is both structurally and functionally a monomer in the membrane.

The lactose operon regulatory pathway. The lactose operon plays two main important roles in the E. coli metabolism: It controls the production of the enzymes necessary for lactose absorption and hydrolysis, but it also closes a positive feedback loop, the so called lactose operon regulatory pathway. Once the disaccharide lactose is pumped inside the bacteria by the β-galactoside (lac) permease, the second enzyme β-galactosidase has the dual role of transforming the lactose into allolactose and hydrolyzing both (lactose and allolactose) into the monosaccharides galactose and glucose. The positive feedback loop is closed when the intermediary sugar allolactose interacts with the control mechanisms of the lactose operon. Thus the allolactose binds to the lactose repressor lacI reducing its ability to repress the transcription and expression of the structural genes lacZ, lacY, and lacA. We refer the reader to the cartoon in Fig. 3.5 for a better understanding. Consequently an increment in the concentration of lactose or allolactose inside the bacteria enhances the production of the enzymes β-galactosidase and β-galactoside permease, via the expression of the structural genes lacZ and lacY. This incremental enzyme production enhances the absorption of more external lactose and its transformation into allolactose, closing the feedback loop.

In summary, the lactose operon is an excellent example of the inducible operon reviewed in Section 2.2. However, it took a while to interpret the lactose operon subtle behaviour in terms of what we now call bistability. This interpretation was first introduced by Novick and Wiener [25] and Cohn and Horibata [26], who suggested that a single cell may have two alternative states: induced, in which it can metabolize lactose, or uninduced, in which the corresponding genes are switched off and lactose metabolism does not occur. From their results, Novick and Wiener, as well as Cohn and Horibata, interpreted the so called maintenance effect as the consequence of a high permease concentration in induced cells, which would enable these cells to maintain the induced state and to transmit it to their progeny, even if placed in a medium with a low concentration of inducer. Although this interpretation accounts for the existence of two distinct phenotypes and provides an explanation of why induced cells placed in media with low inducer concentrations remain indefinitely induced, whereas cells that have never been induced stayed uninduced, it does not explain what makes the cells switch between alternative states. This switching remained a mystery that had to wait for the introduction of the concept of multistability to be fully explained.

We have seen in Section 2.2 that Griffith [5] introduced a mathematical model for a single gene controlled by a positive feedback loop, and found that, under certain conditions, two stable states may be accessible for the system simultaneously. However, Griffith did not use his model to explain the maintenance effect of the lac operon. The first models explicitly aimed at unraveling this phenomenon were due to Babloyantz and Sanglier [27] and to Nicolis and Prigogine [28], who were able to interpret the maintenance effect as the biological facet of the physical process of multistability. These models were quite complex, and took into account all the information regarding the lactose operon regulatory pathway available at the time. However, the level of detailed knowledge about the underlying molecular mechanisms has expanded greatly in the intervening decades. Thus, more detailed and sophisticated models are possible. Below, we review some of the most recent modeling studies of the lactose operon, many of which are by our group.

Transcription of the structural genes. Let P (O P ) be the probability that a polymerase is bound to the promoter/operator region of the lactose operon and it is ready to initiate transcription. The dynamical equations for the lacZ and lacY ribosome binding sites (RBSs) in the mRNA molecule are given [30–32] by

(3.1)
dMZdt=DkMPτZ(OP)-(γM+μ)MZ,
(3.2)
dMYdt=DkMPτY(OP)-(γM+μ)MY.

Variable M Z and M Y respectively denote the concentrations of lacZ and lacY RBSs. D stands for the concentration (number of molecules per average bacteria) of lactose promoters, k M is the maximum transcription initiation rate of the promoter, γ M denotes the mRNA degradation rate, and μ is the average bacterias grown rate. μ is included along with the degradation rate γ M to account for the effective loss due to dilution. Both (3.1) and (3.2) share the same parameters because the structural genes lacZ and lacY are located in tandem after the promoter, and thus they are transcribed by the same polymerase one after the other. Finally, the notation P τ Z  (*) (t) stands for P (*) (t - τ Z ), and we use it to take into account the time delay τ Z existing between transcription initiation and translation initiation. Hence, τ Z is the time interval between transcription initiation and the moment when the corresponding RBS is transcribed, so that a ribosome can bind to it and initiate the translation. Obviously, the time delay τ Y is larger than τ Z , because the structural genes lacZ are located close to the promoter and so are transcribed first. Note that the symmetry between Equations (3.1) and (3.2) implies that M Y  (t) is equal to M Z  (t - τ) for the difference τ = τ Y  - τ Z , so that we need to use only one of these equations.

Translation of mRNA. The dynamical equations for the concentration of the proteins encoded by the genes lacZ and lacY are given [30–32] by

(3.3)
dEZdt=kZe-μτZ*MZ,τZ*-(γZ+μ)EZ,
(3.4)
dEYdt=kYe-μτY*MY,τY*-(γY+μ)EY.

The variable E Z (E Y ) denotes the concentration of LacZ (LacY) polypeptides. The parameter k Z stands for the maximum translation initiation rate at the lacZ RBS, τZ* is the time necessary to fully translate a LacZ polypeptide, γ Z denotes the protein E Z degradation rate, and μ is as before. The exponential factor e-μτZ* accounts for dilution of mRNA concentration due to cell growth in the time interval [t-τZ*,t]. Finally, the notation MZ,τZ*(t) stands for the delayed function MZ(t-τZ*). The parameters k Y , τY*, and γ Y in Equation (3.4) have the same meaning as above for the dynamics of protein E Y . Since the lacY mRNA segment has its own ribosome binding site, it is translated independently from lacZ mRNA segment.

Observe that if the set of parameters (kZ,τZ*,γZ) is identically equal to (kY,τY*,γY), then the symmetry between Equations (3.3) and (3.4) implies that E Y  (t) is equal to E Z  (t - τ) for τ = τ Y  - τ Z , because we already know that M Y is equal to M Z,τ .

Lactose absorption and hydrolysis into lactose and allolactose. Once the lacZ and lacY polypeptides are produced, they pass through several biochemical processes like folding and tetramerization in order to produce the corresponding enzymes β-galactosidase and β-galactoside permease. The internal dynamics of these biochemical processes are not modeled in general (the corresponding reversible reactions are assumed to always be in equilibrium), and so one may take

(3.5)
B=EZ/4andQ=EY,
where B is the internal concentration of β-galactosidase and Q denotes the concentration of β-galactoside permease. The factor 1/4 comes from the fact that β-galactosidase is a homo-tetramer made up of four identical lacZ polypeptides. We thus assume that all the β-galactosidase monomers are incorporated into tetramers.

Dynamical equations for the concentration of intracellular lactose L in bacteria were developed in [30] and [31], and then later improved [32] to include explicitly the effects of the external glucose G e in the absorption of lactose. This latter formulation took the form

(3.6)
dLdt=kLβL(Le)βG(Ge)Q-kβ(L)Q-φMM(L)B-(γL+μ)L.
L, as before, is the concentration of intracellular lactose, while G E (L E ) denotes the concentration of extracellular glucose (lactose). The first term k L β L β G Q in (3.6) stands for the gain of intracellular lactose L obtained due to the action of the β-galactoside permease Q in the transport of extracellular lactose L; the second term k β  (LQ expresses the loss of intracellular lactose to the extracellular fluid due to the reversible nature of the permease mediated transport; the third term φMM(L)B accounts for the β-galactosidase mediated conversion of lactose into allolactose as well as the hydrolysis of lactose to glucose and galactose. The last term in (3.6) stands for the decrease in internal lactose due to degradation and dilution. β L  (L e ) is an increasing function of the extracellular lactose L e , and β G  (G e ) is decreasing with respect to the external glucose G e to take into account the negative influence of the glucose on the absorption of lactose:
(3.7)
βL(Le)=LeκL+LeandβG(Ge)=1-φGGeκG+Ge.
Furthermore, the terms β  (L) and (L) are both functions of the internal lactose
(3.8)
β(L)=Lκ+LandM(L)=LκM+L.

The dynamical equation for the concentration of allolactose A is much simpler:

(3.9)
dAdt=αφMM(L)B-φAM(A)B-(γA+μ)A,
where α is the fraction of internal lactose L transformed by β-galactosidase B into allolactose instead of being hydrolyzed into glucose and galactose. The term φAM(A)B represents the hydrolysis of allolactose into glucose and galactose mediated by β-galactosidase, while the last term in (3.9) stands for the decrease in internal allolactose due to degradation and dilution. We implicitly assume that the dynamics of lactose and allolactose hydrolysis are so similar that the same functions (L) and (A) can be used to represent both.

In particular, if αφ M  ≃ φ A holds, γ A  + μ is close to zero, and the allolactose dynamics are fast (so that Equation (3.9) is always close to equilibrium), then we conclude that A ≈ L and is independent of B.

The lactose operon control system. The system of Equations (3.1) to (3.9) gives a mathematical model of the biochemical reactions involved in the transcription and translation of the lac structural genes, the absorption of the extracellular lactose, its later transformation into allolactose, and the hydrolysis of lactose and allolactose into glucose and galactose. The one thing left to specify is an exact expression for the probability P (O P ) that a polymerase is bound to the promoter/operator region of the lactose operon and it is ready to initiate transcription. We need an explicit formula for P (O P ) in order to substitute it into Equations (3.1)(3.2) and to model how allolactose and glucose control the production of the enzymes necessary for the lactose absorption, transformation, and hydrolysis, closing in this way the positive feedback loop described previously.

The system (3.1) to (3.9) was presented in [30–32] and has not been significantly modified since the time it was originally developed. However, the probability P (O P ) has changed significantly from the original form

P(OP)=a+Anb+An
proposed by [30].

Other investigators [29, 33–35] have proposed different formulas for P (O P ) adding more and more new details on the lactose operon control system, which is quite complex as the most recent discoveries show. Thus [36] and [37] have established that the lactose operon regulatory elements (pictured in Fig. 3.3a) are distributed along the DNA chain as follows: the lactose promoter is located between bp -36 (bp stands for base pair, and positions are referred relative to the starting point of gene lacZ, bp +1) and bp -7. Operator O1 is 21 bp long and centred around bp +11. There are two additional operators, denoted O2 and O3, which are, respectively, located at 401 bp downstream and 92 bp upstream from O1. Finally, the activator (CAP)-binding site spans from bp -72 to bp -50.

The lactose repressor is a homo-tetramer (consisting of two functional homo-dimers) of lacI polypeptides, according to [38] and [39]. Each functional dimer can bind operators O1, O2 and O3. Furthermore, DNA can also fold in such a way that a single repressor binds two operators simultaneously, one per dimer. Each monomer in the lactose repressor can be bound by an allolactose molecule, inhibiting the capability of the corresponding dimer to bind an operator. This means that free repressors can bind one operator (Fig. 3.3b) or two of them simultaneously (Fig. 3.3c), repressors with three free monomers can bind one but not two operators (Fig. 3.3d), repressors with two free monomers can bind one operator, if the bound monomers belong to the same dimer (Fig. 3.3e), or none at all, and that repressors with only one free monomer are unable to bind any operator, as are repressors with all four monomers bound by allolactose; see for example [40].

Deletion experiments [41] have shown that a repressor bound to O1 inhibits transcription initiation, while a repressor bound to either O2 or O3 has almost no effect on the expression of the lactose operon structural genes. Nevertheless, O2 and O3 do have an indirect effect because the complex formed by a single repressor simultaneously bound to O1 and either O2 or O3 is far more stable than that of a repressor bound only to O1. The consequence of this is that interacting with the lactose repressor operator O1 is only capable of decreasing the expression of the operon genes 18 times; when it cooperates with O2, the repression level can be as high as 700-fold; when O1 and O3 act together, they can reduce the operon activity up to 440 times; when all three operators are present, the repression intensity can be as high as 1300-fold.

Also, in [36] it has been established that the intracellular production of cyclic AMP (cAMP) decreases as the concentration of extracellular glucose increases. cAMP further binds a specific receptor molecule (CRP) to form the so-called CAP complex. Finally, CAP binds a specific DNA site (denoted here as C) upstream from the lac promoter, and by doing so it increases the affinity of the mRNA polymerase for this promoter. This regulatory mechanism is known as cataboliterepression.

A novel source of cooperativity has been recently discovered [42] in the lactose operon: when a CAP complex is bound to site C, it bends DNA locally and increases the probability of the complex in which a repressor simultaneously binds operators O1 and O3.

The last regulatory mechanism in the lac operon is a so-called inducer exclusion. In it, external glucose decreases the efficiency of lac permease to transport lactose, and by doing so negatively affects the induction of the operon genes; see for example [36].

These regulatory mechanisms which we have briefly reviewed above are summarized in Fig. 3.3. As we have seen, the activity of the lactose operon is regulated by extracellular glucose and lactose. While extracellular glucose decreases the operon activity via catabolite repression and inducer exclusion, extracellular lactose increases the operon expression level by deactivating the repressor. Another important point is the existence of a positive feedback loop: as more molecules of lactose permease and β-galactosidase are produced, there is an elevated lactose uptake flux and an increased lactose metabolism rate; this further increases the production of allolactose and, as a consequence, diminishes the amount of active repressor. This, in turn, increases the operon activity, and thus more lactose permease and β-galactosidase molecules are produced.

The reader interested in the details of the lac operon regulatory mechanisms is referred to the excellent review [43] and the references therein. A good description of the operon regulatory elements and their location on the DNA chain can be found in [36]. The most recent discoveries regarding the cooperativity between CAP-binding site and operator O3 are [42].

Probability that a polymerase is bound to the promoter and a transcription initiates. Santillan and co-workers ([34] and [29]) have taken into account all the details of the lactose operon control system described above and deduced an explicit formula for the probability P (O P ) as a function of the allolactose A and external glucose G e concentrations. This rather complicated expression is given by

(3.10)
P(OP)=ppc(Ge)PR(A),
(3.11)
ppc(Ge)=pp1+(kpc-1)pc(Ge)1+(kpc-1)pppc(Ge),
(3.12)
pcp(Ge)=pc(Ge)1+(kpc-1)pp1+(kpc-1)pppc(Ge),
(3.13)
pc(Ge)=KGmKGm+Gem,
(3.14)
PR(A)=(1+ξ2ρ(A))(1+ξ3ρ(A))+ξ1*ρ(A)2Z(A)+j=1,2,3(1+ξjitscρ(A)),
(3.15)
Z(A)=j=1,2,3pcp(Ge)δ2j(1+ξjρ(A))ξj*ρ(A)2,
(3.16)
ρ(A)=(KAKA+A)2.
In the following few paragraphs we explain, step by step, the elements of this expression.

The function 𝒫R(A) in (3.14) accounts for the regulation of transcription initiation by active repressors, giving the probability that the lactose promoter is not repressed by an active repressor bound to Operator O1. It accounts for the interactions of the repressor and allolactose molecules, of the repressor molecules and the three different lactose operators (including DNA looping), of the CAP activator and the mRNA polymerase, and of CAP and the DNA loop involving operators O1 and O3.

Repressor molecules are tetramers formed by the union of two active dimers. Every one of the four repressor subunits can be bound by an allolactose molecule. According to [40], free repressors, repressors bound by one allolactose, and repressors bound by two allolactoses in the same dimer can bind a single operator. The fraction of repressors able to do so is denoted by ρ (A) in (3.16). Conversely, only free repressors, whose fraction is given by ρ (A2, can bind two different operators simultaneously. The function p pc  (G e ) in (3.11) denotes the modulation of transcription initiation by the cooperative interaction between a CAP activator and a polymerase, each bound to its respective site. Production of cyclic AMP (cAMP) is inhibited by extracellular glucose G e . cAMP further binds the so-called cAMP receptor protein to form the CAP complex. Finally, CAP binds a specific site near the lactose promoter and enhances transcription initiation. The probability of finding a CAP molecule bound to its corresponding site is given by p c  (G e ) in (3.13).

The probability that a CAP activator is bound to its corresponding site is given by the function p cp  (Ge) in (3.12). Its presence in the definition of 𝒫R(A) in (3.14) accounts for the fact that it affects the formation of the DNA loop in which a single repressor binds operators O1 and O3 at the same time. Note that (p cp δ 2j is equal to p cp only when j = 2 and it is equal to one in any other case.

Reduced model. The system of equations developed above can be reduced after assuming that the set of parameters (τZ,kZ,τZ*,γZ) is equal to (τY,kY,τY*,γY), because in this case Equations (3.1) and (3.2) are identical, and in the same way (3.3) is identical to (3.4). Thus, recalling Equations (3.6) and (3.10), we obtain the reduced system

(3.17)
dMZdt=DkMppc(Ge)PR(A)-(γM+μ)MZ,
(3.18)
dEZdt=kZe-μτZ*MZ,τZ*-(γZ+μ)EZ,
(3.19)
dLdt=kLβL(Le)βG(Ge)Q-kβ(L)Q-φMM(L)B-(γL+μ)L.

The functions p pc  (G e ) and 𝒫R(A) are given in Equations (3.10) to (3.16). Finally, if we assume in (3.9) that the equality αφ M  = φ A holds, the sum γ A  + μ is very small (close to zero), and the allolactose dynamics is very fast, then we can assume that A = L. Thus, we complete the model for the lactose operon by adding the Equations (3.5) to (3.8),

(3.20)
B=EZ/4,
(3.21)
Q=EZ,
(3.22)
A=L,
(3.23)
βL(Le)=LeκL+Le,
(3.24)
βG(Ge)=1-φGGeκG+Ge,
(3.25)
β(L)=Lκ+L,
(3.26)
M(L)=LκM+L.

The parameters of the model (3.10) to (3.26) are given in Table 3.1 as estimated from the experimental literature, see [32, 34] and [29].

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

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

Noting that TMG inactivates the lactose repressor, but it is not metabolizable, we simulate the Ozbudak et al. experiments. For this, we set φ M  = 0/min to account for the presence of a reliable carbon source (succinate) and induction with TMG, which is not metabolized by β-galactosidase. Then, we calculated the bifurcation points and plotted them in the L e versus G e parameter space. We took K A as a free parameter, and found that K A  = 8.2 × 105 mpb (here and thereafter mpb means molecules per average-size bacterium) gives a reasonable fit to the experimental points of Ozbudak et al. Both the model bifurcation diagram and the experimental points are presented in Fig. 3.4A. Note that the bistability region predicted by the model is wider than the experimental one. There are three possible explanations for this discrepancy: 1) the lactose promoter-gfp fusion employed by Ozbudak et al. as a reporter lacks operators O2 or O3; 2) the difficulty in measuring exactly the L e values at which the bimodal histograms appear and disappear; and 3) the phase diagram of Fig. 3.4A is based upon a mean-field analysis, and so biochemical noise can change the phase boundaries; see Section 4 below. A fourth possible explanation for the disagreement between the model and the experimental results is that our estimated parameter values differ from those corresponding to the E. coli strain used by Ozbudak et al. To account for this possibility, we explored the parameter space looking for a better fit. We found that it can be obtained by decreasing the parameters ξ j and ξj* to 15% of the values reported in Table 3.1, and by setting K A  = 2.8 × 106 mpb. The results are shown in Fig. 3.4B.

3.2The tryptophan operon

Tryptophan is one of the 20 amino acids out of which all proteins are made. Arguably, tryptophan is the most expensive amino acid to synthesize, biochemically speaking. Perhaps, for this reason, humans and many other mammals do not have the enzymes necessary to catalyze tryptophan synthesis and instead they find this amino acid in their diet.

However, microorganisms like E. coli generally posses the machinery to produce tryptophan, but the production process is tightly regulated in all cases. In the particular case of E. coli, the tryptophan operon is a DNA segment containing a promoter (trpR) where transcription starts and regulation by repression takes place, a leader region (trpL) where regulation by transcriptional attenuation occurs, and five structural genes (trpE to trpA) that code for the polypeptides comprising the enzymes responsible for the catalysis of tryptophan biosynthesis. There are three different regulatory mechanisms involved in the control of the tryptophan operon dynamics: repression, transcriptional attenuation, and enzyme inhibition. The tryptophan regulatory pathway is illustrated in Fig. 3.5.

Repression occurs when an active repressor binds to one of the three available binding sites within the promoter, inhibiting the binding of a RNA polymerase, and so of transcription initiation. The repressor molecule is a homo-dimer made up of two TrpR polypeptides. Each subunit has a binding site for tryptophan, and the repressor molecules activate when both tryptophan binding sites are occupied. Of the three repressor binding sites within the promoter, the two closest to the transcription initiation site interact cooperatively. That is, when two are bound to such sites, the resulting complex is much more stable than it would be expected from the addition of the individual binding energies.

Transcriptional attenuation is regulated by the DNA leading region. The RNA strand resulting from transcription of trpL can fold into three alternative hairpin-like structures, as a result of Watson-Crick base pairing. Soon after transcription initiation, the first hairpin structure is formed, and this causes the polymerase to pause transcription. When a ribosome binds to the nascent RNA strand to start translation, it eventually disrupts the hairpin and both transcription and translation proceed together. Not long after that, the ribosome encounters two tryptophan codons in tandem. Under conditions of abundant tryptophan, there is a large amount of charged trp transfer RNA, and so the two consecutive tryptophan codons are rapidly translated. When this occurs, a second hairpin structure, that serves as a transcription termination signal, forms and transcription is prematurely aborted. Conversely, if tryptophan is scarce, the ribosome stops at the trp codons while the RNA polymerase continues transcribing the rest of the leading region. This prevents the formation of the transcription-terminating hairpin and instead promotes the formation of a third structure that allows the polymerase to go into the structural genes to transcribe them.

Tryptophan biosynthesis takes place through a series of reactions, each one catalyzed by enzymes formed from the polypeptides coded by genes trpE-A. The first of those reactions, and the slowest one, is catalyzed by the enzyme anthranilate synthase. In this reaction, anthranilate is synthesized out of chorismic acid. Being the slowest reaction of the tryptophan synthesis path, anthranilate synthesis determines the velocity of the whole process. Furthermore, anthranitale synthase is a heterotetramer made up of two TrpE and two TrpD subunits. Each TrpE subunit has a binding site for tryptophan, and when they are bound by this amino acid, the whole enzyme suffers an allosteric transformation that makes it unable catalyze the corresponding reaction. This regulatory mechanism is known as enzyme inhibition.

A deterministic model for this regulatory pathway can be constructed as follows. Consider first the dynamics of promoter switching. Denote the state of repression of the promoter as (i, j, k)—with i, j, k = 0, 1; a value of 1 means that the corresponding repressor binding site is occupied, while a value of 0 means that it is empty. If P ijk represents the average number of promoters whose state is (i, j, k), the chemical reactions through which the promoter switches between its different available states are:

P000β1α1P100,P000β2α2P010,P000β3α3P001,P100β2/kCα2P110,P100β3α3P101,P010β1/kCα1P110,P010β3α3P011,P001β1α1P101,P001β2α2P011,P110β3α3P111,P101β2/kCα2P111,P011β1/kCα1P111,
In these reactions α i represents the effective reaction rate constant for the binding of an active repressor to the i-th binding site in the promoter, β i is the corresponding unbinding reaction rate constant, and k C accounts for the cooperativity between the first two repressor binding sites.

By making use of the theory of chemical kinetics we can write the following set of differential equations governing the dynamics of variables P ijk :

(3.27)
dP000dt=-(α1+α2+α3)P000+β1P100+β2P010+β3P001,
(3.28)
dP100dt=-(β1+α2+α3)P100+α1P000+β2kCP110+β3P101,
(3.29)
dP010dt=-(α1+β2+α3)P010+β1kCP110+α2P000+β3P011,
(3.30)
dP001dt=-(α1+α2+β3)P001+β1P101+β2P011+α3P000,
(3.31)
dP110dt=-(β1kC+β2kC+α3)P110+α1P010+α2P100+β3P111,
(3.32)
dP101dt=-(β1+α2+β3)P101+α1P001+β2kCP111+α3P100,
(3.33)
dP011dt=-(α1+β2+β3)P011+β1kCP111+α2P001+α3P010,
(3.34)
dP111dt=-(β1kC+β2kC+β3)P111+α1P011+α2P101+α3P110.
These equations do not constitute a complete set because the effective binding reaction rate constants α i are directly proportional to the amount of active repressors, R A , which is, in turn, a function of the intracellular tryptophan concentration.

To complete the differential equation system let M represent the concentration of mRNA molecules resulting from transcription of the tryptophan operon, E be the concentration of anthranilate synthase enzymes, and T denote the intracellular tryptophan level. Following the development in previous sections, the differential equations accounting for the dynamics of these variables are:

(3.35)
dMdt=kMP000A(T)-γMM,
(3.36)
dEdt=kEM-γEE,
(3.37)
dTdt=kTEI(T)-γTTT+KT,
in which k M is the transcription initiation rate, 𝒜(T) represents the probability that a newly initiated transcriptional event is not prematurely aborted due to attenuation, γ M accounts for the mRNA degradation rate, k E is the enzyme synthesis rate per mRNA molecule, γ E is the enzyme degradation rate, k T represents the tryptophan synthesis rate per active enzyme, (T) is the probability that an enzyme is not inhibited by tryptophan, γ T is the maximal tryptophan consumption rate due to the cellular metabolism, and K T is the corresponding half saturationconstant.

The reaction rate constants for repressor binding are proportional to the concentration of active repressors R A  (T). That is,

(3.38)
αi=aiRA(T).
Thus, expressions for R A  (T), 𝒜(T), and (T) are required to complete the model. These functions correspond to the three known regulatory mechanisms in this system: repression, transcriptional attenuation, and enzyme inhibition, respectively. Functions R A  (T), 𝒜(T), and (T) were derived in [45] from chemical kinetics considerations by taking into account all the chemical reactions behind the corresponding regulatory mechanisms. The resulting expressionsare:
(3.39)
RA(T)=RTot(TT+KA)2,
(3.40)
A(T)=1+2αTKG+T(1+αTKG+T)2,
(3.41)
I(T)=KInKIn+Tn.
Here, R Tot is the total number of repressor molecules, K T the dissociation constant between tryptophan and one binding site of a repressor, K G is the dissociation constant between tryptophan and the corresponding transfer RNA, α the probability per unit time that a charged tRNATrp arrives at a tryptophan codon so that it is translated, K I is the dissociation constant between tryptophan and one of the TrpE subunits in anthranilate synthase, and n is a Hill coefficient.

Equations (3.27)(3.41) constitute a complete system of differential equations that model the dynamics of the tryptophan operon. However, due to its high dimensionality, this system is quite difficult to analyze. For that reason, some simplifying assumptions are useful. One which has been widely employed consists in supposing that the dynamics of repressor binding and unbinding are much faster than those of mRNA and protein synthesis and degradation, as well as those of tryptophan production and consumption. If this is the case, the subsystem given by Equations (3.27)(3.34) is much faster than that given by Equations (3.35)(3.37), and so one can make a quasi steady state approximation (also known as adiabatic elimination) for Equations (3.27)(3.34), with which the model transforms into

(3.42)
dMdt=kMP000(T)A(T)-γMM,
(3.43)
dEdt=kEM-γEE,
(3.44)
dTdt=kTEI(T)-γTTT+KT.
The concentration of non-repressed promoters is given in this case by

(3.45)
P000(T)=(1+a1β1RA(T)+a2β2RA(T)+a3β3RA(T)+kCa1β1a2β2RA2(T)+a1β1a3β3RA2(T)+a2β2a3β3RA2(T)+kCa1β1a2β2a3β3RA3(T))-1.
P 000 (T), 𝒜(T), and (T) are monotonic sigmoidally decreasing functions of T, and so is the product P000(T)A(T). This product is sometimes replaced by a decreasing Hill function [46].

As explored extensively in Section 2, an elementary classification of systems subject to feedback regulation includes those with negative feedback or, alternately, those with positive feedback. This is important because the type of feedback determines the kind of expected dynamic behaviour. Thus, positive feedback is necessary for bistability, while negative feedback is the mechanism underlying cyclic behaviour. Given that the tryptophan operon has been experimentally studied for several decades (and, thus, is one of the best known molecular systems), and that it is regulated by three different negative feedback loops, this system has become a paradigm for studying the effects of negative feedback regulation on gene expression. Below we review some of the most prominent past studies of the tryptophan operon.

As discussed in Section 2.2, the first mathematical model for a repressible operon was due to Goodwin [3], who developed a model with a structure equivalent to that in Equations (3.42)(3.44), except that the regulatory functions accounting for transcriptional attenuation and enzyme inhibition were not taken into account. The repression regulatory function was modeled in the Goodwin model by a monotone decreasing Hill function. In a later paper, Goodwin [2] presented analog computer simulations of limit cycles (sustained oscillations) obtained from this model with a Hill exponent of one. However, Griffith [4] later demonstrated that the steady state is locally stable up to a Hill exponent equal to 8, making limit cycle oscillations highly unlikely for low exponent values. In a large number of simulations Griffith found limit cycles only if the steady state was unstable. Apparently, there was an error in Goodwin’s analog simulation. The controversy was finally resolved by Tyson [47], who analytically proved the existence of at least one periodic solution whenever the steady state is unstable.

As we saw in the previous paragraph, the first modeling studies on a repressible operon focused on the possibility of sustained oscillations, and ended with a negative conclusion. This question was revisited in [48], who modified the Goodwin model to include the transcriptional and translational time delays as well as the regulatory function accounting for enzyme inhibition. Bliss et al. demonstrated that time delays can induce sustained oscillations, but only when enzyme inhibition is weakened. They also presented experimental results with a mutant strain of E. coli in which the enzyme anthranilate synthase cannot be inhibited by tryptophan. This strain was first grown in a tryptophan-rich medium and then suddenly changed to a tryptophan-less medium to induce expression of the tryptophan operon genes. Both the simulations and the experiments showed periodic oscillation in both the enzyme and the tryptophan intracellular concentrations.

In later work [49] the Goodwin model was further refined by deriving a repression regulatory function from first principles, taking into consideration the underlying chemical reactions. Nonetheless, they dismissed the regulatory functions corresponding to transcription attenuation and enzyme inhibition. In [49] the possible complex behaviours the tryptophan operon can show, given the architecture of the regulatory network, was investigated. They found that the steady state, although normally stable, becomes unstable for super-repressing strains, even at low values of the cooperativity of repression. However, in order for this to happen it is necessary that the demand for end-product saturates at large end-product concentrations. Finally, in [49] it was proved that the system can also show bistability, in which a stable steady state and a stable limit cycle coexist.

In 1990, other investigators [50] introduced one more model for the tryptophan operon regulatory pathway, and used it to investigate the possibility of engineering an E. coli strain to overproduce tryptophan. The model [50] has a similar structure to that in Equations (3.42)(3.44) but, as some of the models reviewed in the former paragraphs, it ignores the transcriptional attenuation and the enzyme inhibition regulatory mechanisms. Through analytical studies and numerical simulations the authors were able to demonstrate that stable overproduction is feasible. Nevertheless, under some specific circumstances the operon may become unstable and lead to periodic synthesis. In [51] the models of [49] and [50] were further refined and employed to study the influence of periodic fluctuations in the intracellular demand for tryptophan.

In our group we have studied the dynamic behaviour of the tryptophan operon regulatory pathway for some time. In [52] we developed a mathematical model that accounts for all known regulatory mechanisms, as well as for the time delays due to transcription and translation. Moreover, we put special attention to estimating all of the model parameters from reported experimental data. Although involving one extra differential equation, a more careful analysis reveals that this model is equivalent to that in Equations (3.42)(3.44). To test the model feasibility, we compared its predictions with available dynamic experiments from wild type and two mutant strains. Later [45], we simplified our original model, but still considered all three existing regulatory mechanisms, and analyzed their influence on the system dynamic behaviour. We numerically showed that enzyme inhibition is the fastest responding mechanism. However, although it could suffice to efficiently control tryptophan biosynthesis, it would be very expensive because it would imply continuous production of enzymes. Although repression and transcription attenuation respond considerably more slowly, they allow bacteria to diminish the energy expended in enzyme synthesis when tryptophan demand is low for longer periods of time. In other words, the redundancy of feedback regulatory mechanisms allows E. coli to efficiently respond to both slow (via repression and transcription attenuation) and fast (via enzyme inhibition) fluctuations of tryptophan demand. These numerical results were analytically corroborated in [53], where we studied the global stability of the tryptophan operon model using the second Lyapunov method.

As we have seen, the first modeling studies on the tryptophan operon focused on the possibility that this system shows sustained oscillations under given circumstances. Interestingly, there is only one experimental report of such oscillatory behaviour in the tryptophan operon [48]. Taking into consideration the lack of experimental evidence, as well as recent discoveries regarding the existence of multiple repressor binding sites within the trp promoter, and of cooperativity between two of them, we have further investigated the possibility of observing sustained oscillations in this system [54]. To that end, we improved the model in [45] by incorporating the discoveries discussed above and analyzed it numerically. We found that indeed, a mutant bacterial strain lacking enzyme inhibition can behave cyclically, and that the time delays due to transcription and translation are essential for this behaviour. In Fig. 3.6 we show the model results, which show a very good agreement with the experimental results in [48].

On the other hand, regular periodic oscillations are observed in the model of [54], but only when the system intrinsic stochasticity is ignored. When the so-called intrinsic biochemical noise is taken into account, the system shows oscillations with variable periods, and this causes the global system behaviour in a cell population to be non cyclic overall. These results stress the necessity of further studying the appearance of oscillations in the tryptophan operon, both analytically and experimentally; not only to satisfy some people’s scientific curiosity, but also because answering this question may shed some light into the dynamics of gene regulation. In Fig. 3.7 we show the stochastic quasi-periodic dynamic behaviour predicted by the mathematical model, as well as the average of 100 independent cells.

All the models reviewed so far have the structure of the model represented by Equations (3.42)(3.44). This means that, either explicitly or implicitly, they assume that promoter gating between the various repressed and the non-repressed states is much faster than the transcription and translation processes. Nevertheless, recent detailed measurements of the repressor-promoter kinetics revealed that this assumption is not valid—see [55] and references therein. This further implies that the assumed separation of time scales employed to obtain the simplified model in Equations (3.42)(3.44) does not exist, and one is obliged to work with the full model: Equations (3.27)(3.37). In a recent paper [55] we studied the stochastic behaviour of such a model, but analytical and numerical studies of the deterministic counterpart are still missing.

We wish to emphasize that in our modeling studies we have followed the strategy of producing models as detailed as possible, given the available experimental evidence. This meant that not only we included all known mechanisms into the model equations, but also that we estimated all of the model parameters from reported experimental data. Understandably, this is not always possible when developing models for biological systems. However, in this particular case, the tryptophan operon of E. coli is so well studied that developing this kind of model is completely feasible. A natural consequence of having quite detailed models is the possibility of accurately reproducing dynamic experiments. In particular, we have employed the experimental results of [56] to compare with our models’ predictions. In Fig. 3.8 we show comparisons of model predictions and the Yanofsky and Horn experimental measurements for a wild type and for a enzyme-inhibition-less mutant E. coli strain. The theoretical simulations in Fig. 3.8 were carried out with our most detailed model [55], but qualitatively similar results are obtained with all the model versions previously reviewed. In our opinion, it is essential for a model to be able to reproduce existing dynamical experimental data, before it can be employed to answer dynamical questions not easily addressed experimentally.

E. coli is not the only bacterium with a tryptophan operon. Other bacteria also have an equivalent system, in particular B. subtilis. Interestingly, the structure of the regulatory pathway in both systems is very similar, although the specific mechanisms are very different. For instance, instead of repression, the tryptophan operon in B. subtilis involves a so-called TRAP molecule that promotes premature transcription termination when it is bound by 11 tryptophan molecules. Instead of transcriptional attenuation, B. subtilis has a secondary at operon that is regulated by tryptophan and produces a protein that modulates the effect of TRAP proteins. The only mechanism that E. coli and B. subtilis share in common is enzyme inhibition. A model for the tryptophan operon of B. subtilis was developed in [57] and shown that not only its regulatory pathway has a similar structure to that of E. coli, but the analogous mechanisms in both systems play similar roles from a dynamic perspective. Given that the lineages of both organisms evolved separately several millions of years ago, these similarities may be the result of evolutionary convergence.

4Noise effects in gene regulation: Intrinsicversus extrinsic

In all areas of science, when making experimental measurements it is noted that the quantity being measured does not have a smooth temporal trajectory but, rather, displays apparently erratic fluctuations about some mean value when the experimental precision is sufficiently high. These fluctuations are commonly referred to as ‘noise’ and usually assumed to have an origin outside the dynamics of the systems on which measurements are being made–although there have been many authors who have investigated the possibility that the ‘noise’ is actually a manifestation of the dynamics of the system under study. Indeed, a desire to find ways to quantitatively characterize this ‘noise’ is what led, in large part, to the development of the entire mathematical field loosely known as stochastic processes, and the interaction of stochastic processes with deterministic dynamics is of great interest since it is important to understand to what extent fluctuations or noise can actually affect the operation of the system being studied.

Precisely the same issues have arisen in molecular biology as experimental techniques have allowed investigators to probe temporal behaviour at ever finer levels, even to the level of individual molecules. Experimentalists and theoreticians alike who are interested in the regulation of gene networks increasingly focus on trying to assess the role of various types of fluctuations on the operation and fidelity of both simple and complex gene regulatory systems. Recent reviews [58–60] give an interesting perspective on some of the issues confronting both experimentalists and modelers.

As in other areas of science, in gene regulation the debate often swirls around whether the fluctuations are extrinsic to the system under consideration [61–64], or whether they are an intrinsic part of the fundamental processes they are affecting (e.g. bursting, see below). The dichotomy is rarely so sharp however, but in [65] an elegant experimental technique has been presented to operationally distinguish between the two, see also [66], while [67] and [68] have partially set the stage for a theoretical consideration of this question. One issue that is raised persistently in considerations of the role of fluctuations or noise in the operation of gene regulatory networks is whether or not they are ‘beneficial’ [69] or ‘detrimental’ [70] to the operation of the system under consideration. This is, of course, a question of definition and not one that we will be further concerned with here since it is a question without scientific meaning.

In this section we study the density of the molecular distributions in generic bacterial operons in the presence of ‘bursting’ (commonly known as intrinsic noise in the biological literature) as well as inherent (extrinsic) noise using an analytical approach. In a very real sense, the whole field of intrinsic noise behaviour owes its basis to the pioneering work of Berg [71] who first studied the statistical fluctuations of protein numbers in bacterial population (with division) through the master equation approach, and introduced the concept of what is now called bursting. Our work is further motivated by the well documented production of mRNA and/or protein in stochastic bursts in both prokaryotes and eukaryotes [72–79], and follows other mathematical contributions [80–101]. We stress, however, that we have not referenced studies in which stochasticity was studied solely using Gillespie simulations since these have become de rigeur for almost all supposed modeling efforts in spite of the fact that in and of themselves they yield little if any real insight.

Because of its relevance to the analysis of experimental data, we emphasize the behaviour of densities of gene regulatory constituents. To our knowledge, the analytical solution of the steady state density of the molecular distributions in the presence of bursting was first derived in [81]. Our approach emphasized here extends these results to show the global stability of the limiting densities and examines their bifurcation structure to give a complete understanding of the effect of bursting on molecular distributions.

4.1Dynamics with bursting

4.1.1Generalities

In this section we model the amount of the dominant protein as a Markov process {x (t)}  t≥0 with values in (0, ∞). Let x (t) denote the amount of the protein in a cell at time t, t ≥ 0. Following [73, 81] we assume that the amplitude of protein production through bursting translation of mRNA is exponentially distributed, that the frequency of bursting φ is dependent on the level of the protein, and that protein molecules undergo degradation with rate γ. We take here φ (x) = γκ b f (x) and κ b  ≡ φ m in contrast to the deterministic case where κ d  = b d φ m . If only degradation were present, then x (t) would satisfy the equation

x(t)=-γx(t),t0.
However, we interrupt the degradation at random times
t1<t2<
occurring with intensity φ, i.e.,
Pr(tk-tk-1>t|x(tk-1)=x)=e-0tφ(e-γsx)ds,t,x>0.
At each t k a random amount e k of protein molecules is produced according to an exponential distribution with density
(4.1)
h(x)=1be-x/b.
Consequently the process is given by
x(t)={e-γ(t-tk-1)x(tk-1),tk-1t<tk,e-γ(tk-tk-1)x(tk-1)+ek,t=tk,k=1,2,

The corresponding master equation for the evolution of the density u (t, x) of x (t) is given by

(4.2)
u(t,x)t-γ(xu(t,x))x=-γκbf(x)u(t,x)+γκb0xf(y)u(t,y)h(x-y)dy.
A stationary solution of Equation (4.2), which now becomes
-d(xu*(x))dx=-κbf(x)u*(x)+κb0xf(y)u*(y)h(x-y)dy,
with h given by (4.1) and nonnegative f, is of the form
(4.3)
u*(x)=Cxe-x/bexp[κbxf(y)ydy],
where 𝒞 is a normalizing constant, if u * is integrable. The next result follows from [102].

Theorem 4.1. Suppose that h is exponential as in (4.1) with b > 0 and that

C:=01xe-x/bexp[κbxf(y)ydy]dx<.
Then u * defined in (4.3) is the unique stationary density of (4.2) and the solution u (t, x) of (4.2) is asymptotically stable in the sense that
limt0|u(t,x)-u*(x)|dx=0
for all initial densities u (0, x).

4.1.2Distributions in the presence of bursting

We consider the situation in which the function f in the burst frequency φ = γκ b f is given [103] by

f(x)=1+ΘxnΛ+Δxn,
where Λ, Δ, n are positive constants and Θ ≥ 0. We take Θ = 1 to get f as defined in (2.10) for both the generic inducible and repressible operons treated in Section 2.1 with the constants Λ, Δ enumerated in Table 2.1. We have
κbxf(y)ydy=xκby[1+ΘynΛ+Δyn]dy=ln{xκbΛ-1(Λ+Δxn)θ},
where
θ=κbnΔ(Θ-ΔΛ).
Thus, the stationary density (4.3) explicitly becomes
(4.4)
u*(x)=Ce-x/bxκbΛ-1-1(Λ+Δxn)θ.
Observe that in the absence of control, i.e., if f ≡ 1 or, equivalently, Θ = Λ = Δ = 1, we obtain, as in [81], the density of the gamma distribution:
u*(x)=1bκbΓ(κb)e-x/bxκb-1,
where Γ (·) denotes the gamma function. In particular, the first two terms of Equation 4.4 are proportional to the density of a gamma distribution.

The analysis of the qualitative nature of the stationary density (4.4) leads to different conclusions for the inducible and repressible operon cases, since the parameter θ is either positive or negative. In the rest of this section we assume that Θ = 1. First note that we have u * (0) =∞ if 0 < κ b Λ -1 < 1 while u * (0) =0 for κ b Λ -1 > 1 in which case there is at least one maximum at a value of x > 0. To calculate the number of maxima we use the fact that u * (x) >0 for all x > 0 and that

u*(x)=u*(x)(κbf(x)x-1b-1x),x>0.
Consequently, we have u*(x)=0 for x > 0 if and only if
(4.5)
κb(xb+1)=1+xnΛ+Δxn.

For θ ≤ 0, as in the case of no control or a repressible operon, we have Λ = 1, Δ ≥ 1, and graphical arguments (see Fig. 4.1) easily show that Equation (4.5) may have none or one solution. Therefore, we have a stationary density which we can classify as

Unimodal type 1 if u * (0) =∞ and u * is decreasing,

or

Unimodal type 2 if u * (0) =0 and u * has a single maximum at a value of x > 0.

Observe that the stationary density u * in the case of the repressible operon is Unimodal of type 1 if 0 < κ b  < 1 and Unimodal of type 2 if 1 < κ b .

For θ > 0, as in the case of an inducible operon, the stationary density becomes

u*(x)=Ce-x/bxκbK-1-1(K+xn)θ,θ=κbn(1-K-1)
and there is the possibility that u * may have more than one maximum, indicative of the existence of bistable behaviour. Graphical arguments (see Fig. 4.2) show that there may be up to three roots of
(4.6)
1κb(xb+1)=1+xnK+xn.
There are two cases to distinguish. If 0 < κ b  < K then u * (0) =∞ and there can be none, one, or two positive solutions to equation (4.6). If 0 < K < κ b then u * (0) =0 and there may be one, two, or three positive roots of equation (4.6). If there are three we label them as x˜1<x˜2<x˜3. The values x˜1,x˜3 will correspond to the location of maxima in u * while x˜2 will be the location of the minimum between them. Consequently, the stationary density u * can be classified as Unimodal type 1, type 2, as well as

Bimodal type 1 if u * (0) =∞ and u * has a single maximum at x > 0,

or

Bimodal type 2 if u * (0) =0 and u * has two maxima at x˜1,x˜3, 0<x˜1<x˜3.

There are two different bifurcation patterns that are possible. In what will be referred as Bifurcation type 1, the maximum at x = 0 disappears when there is a second peak at x=x˜3. The sequence of densities encountered for increasing values of κ b is then: Unimodal type 1 to a Bimodal type 1 to a Bimodal type 2 and finally to a Unimodal type 2 density. Figure 4.3 illustrates Bifurcation type 1, when n = 4, K = 4, b = 1, and κ b increases from low to high values. In the Bifurcation type 2 situation, the sequence of density types for increasing values of κ b is: Unimodal type 1 to a Unimodal type 2 and then a Bimodal type 2 ending in a Unimodal type 2 density. Figure 4.4 shows Bifurcation type 2, when n = 4, K = 4, b=110, and the parameter κ b increases.

To find the analogy between the bistable behaviour in the deterministic system and the existence of bimodal stationary density u * we fix the parameters b > 0 and K > 1 and vary κ b as in Fig. 4.2. In general we can cannot determine when there are three roots of (4.6). Instead, using the argument of Section 2.2.2 one can determine when there are only two roots. Differentiation of (4.6) yields the condition

(4.7)
nxn-1(K+xn)2=1κbb(K-1).
Equations (4.6) and (4.7) can be combined to give an implicit equation for the value of x ± at which tangency will occur

(4.8)
x2n-(K-1)[n-K+1K-1]xn-nb(K-1)xn-1+K=0
and the corresponding values of κ b± are given by
(4.9)
κb±=(x+bb)(K+xn1+xn).
We see then that the different possibilities depend on the respective values of K, κ b-, κ b+, and κ b . Note that it is necessary for 0 < K < κ b in order to obtain Bimodal type 2 behaviour.

We now choose to see how the average burst size b affects bistability in the density u * by looking at the parametric plot of κ b  (x) versus K (x). Define

(4.10)
F(x,b)=xn+1nxn-1(x+b).
Then

(4.11)
K(x,b)=1+xnF(x,b)1-F(x,b)andκb(x,b)=[K(x,b)+xn]x+bb(xn+1).
Figure 4.5 presents the regions of bimodality in the presence of bursting in the (K, b · κ b ) parameter space, which should be compared to the region of bistability in the deterministic case in the (K, κ d ) parameter space ( b is the mean number of proteins produced per unit of time, as is κ d ).

4.1.3Recovering the deterministic case

The deterministic behaviour can be recovered from the bursting dynamics with a suitable scaling limit of parameters. The frequency κ b and the amplitude b are two important parameters in the bursting production, while in the deterministic production there is only κ d . Thus, if we take the limit

b0,κbwithbκbκd,
in the implicit Equations (4.5) which define the maximum points of the steady state density, then we obtain Equations (2.14) and (2.15) which define the stable steady states in the deterministic case.

Recovering Equation (2.17) in the limit implies that the bifurcations will also take place at the same points. Since we have κ b  > K when κ b → ∞, Bimodality type 1 as well as the Unimodal type 1 behaviours will no longer be present. Moreover, the steady-state density u * will became more sharply peaked as b → 0 and the mass will be more concentrated around the larger maximum of u *.

4.1.4A discrete space bursting model

The number of protein molecules in a single cell can also be described as a Markov process with values in the discrete state space {0, 1, 2, …}. Here we follow the approach of [103]. Let X (t) be the number of gene product molecules at time t. If we have X (t) = m then in a small time interval the change in the number of molecules is

mλmm + k,   mγmm1
where γ m , λ m , m ≥ 0, are constants satisfying

(4.12)
λ0>0,γ0=0,γm>0,λm0,m=1,2,,
while k is randomly chosen, independently of the actual number of molecules, according to a probability density function h = (h k k≥1, so that k=1+hk=1, h k  ≥ 0, k ≥ 1. Of particular interest is the case when h is geometric
(4.13)
hk=(1-b)bk-1,k=1,2,,
with b ∈ (0, 1), which is the discrete space analog of the exponential distribution given by (4.1). Let P m  (t) be the probability that the cell at time t has m protein molecules of the gene product. Our general master equation is an infinite set of differential equations

(4.14)
dPmdt=γm+1Pm+1-γmPm+k=1mhkλm-kPm-k-λmPm,m=0,1,,
where we use the convention that k=10=0. We supplement (4.14) with the initial condition P m  (0) = v m , m = 0, 1, …, where v = (v m m≥0 is a probability density function of the initial amount of the gene product.

The equation for the steady state p*=(pm*)m0 of (4.14) is of the form

γm+1pm+1*-γmpm*+k=1mhkλm-kpm-k*-λmpm*=0,m=0,1,,
which is uniquely solvable (up to a multiplicative constant) by
(4.15)
pm+1*=1γm+1k=0mh¯m-kλkpk*,m=0,1,,
where
h¯l=j=l+1hj,l0.
We have the following general result.

Theorem 4.2. [103, Theorem 3.1] Assume condition (4.12) and suppose that a strictly positive p*=(pm*)m0 given by (4.15) satisfies

m=0pm*=1andm=0(λm+γm)pm*<.
Then for each initial probability density functionEquation (4.14) has a unique solution and
limtm=0|Pm(t)-pm*|=0.

In particular, if condition (4.12) holds and h is geometric as in (4.13) then p*=(pm*)m0 as in (4.15) is given by

(4.16)
pm*=p0*λ0γmk=1m-1λk+bγkγk,m=1,2,.
If additionally γ m  = γm, m ≥ 1, with γ > 0 and λ m is a Hill function of the form
(4.17)
λm=λ1+ΘmnΛ+Δmn,
where Λ, Δ, n > 0 and Θ ≥ 0 are constants, then all assumptions of Theorem 4.2 are satisfied, implying that the steady-state density p*=(pm*)m0 given by (4.16) is the discrete state space analog of (4.4).

4.2Distributions with fluctuations in the degradation rate

We now examine the situation in which fluctuations appear in the degradation rate γ of the generic Equation (2.31). If the fluctuations are Gaussian distributed then it follows from standard chemical kinetic arguments [104] that the mean numbers of molecules decaying in a time dt is simply γxdt and the standard deviation of these numbers is proportional to x. Consequently, we replace Equation (2.31) with a stochastic differential equation in the form

dx=γ[κdf(x)-x]dt+σxdw,
where w is a standard Brownian motion and we use the Ito interpretation of the stochastic integral. The corresponding Fokker Planck equation for the evolution of the ensemble density u (t, x) is given by [105]

(4.18)
u(t,x)t=-[(γκdf(x)-γx)u(t,x)]x+σ222(xu(t,x))x2.

Since concentrations of molecules cannot become negative the boundary at x = 0 is reflecting and the stationary solution of Equation (4.18) is given by

u*(x)=Cxe-2γx/σ2exp[2γκdσ2xf(y)ydy].
Set κ e  = 2γκ d /σ 2. Then the stationary density is given explicitly by
(4.19)
u*(x)=Ce-2γx/σ2xκe11[+Δxn]θ,
where Λ, Δ ≥ 0 and θ are given in Table 2.1. It follows from [106, Theorem 2] that the unique stationary density of Equation (4.18) is given by Equation (4.19) and that u (t, x) is asymptotically stable.

The form of the stationary density for the situation with bursting (intrinsic noise) and extrinsic noise are identical, provided that one replaces the average burst amplitude b with b → σ 2/2γ ≡ b w and κ b  → κ e  = 2γκ d /σ 2 ≡ κ d /b w . Consequently, all of the results of the previous section can be carried over here. In particular, the regions of bimodality in the (K, κ d )-plane can be identified for a fixed value of b w . We have the implicit equation for x ±

x2n-(K-1)[n-K+1K-1]xn-nbw(K-1)xn-1+K=0
and the corresponding values of κ d are given by
κd±=(x+bw)(K+xn1+xn).
Then the bimodality region in the (K, κ d )-plane with noise in the degradation rate is the same as the bimodality region for bursting in the (K, b )-plane.

Finally, we can recover the deterministic behaviour from a limit in the extrinsic fluctuations dynamics. In this case, however, the frequency and the amplitude of the perturbation are already scaled. Then the limit σ → 0 gives the same result as in the deterministic case.

5Discussion and conclusions

Here we have attempted to give an overview of the mathematical techniques that have been used to gain understanding about the operation of bacterial operons. We have looked at generic deterministic models in a very general sense followed by more realistic considerations of both the lactose and tryptophan operons. These two examples are ones for which we have, arguably, the most extensive knowledge of the underlying biology as well as good data and if we cannot successfully understand their operation from a modeling perspective then there is little hope for more complicated situations. Finally we have discussed very recent results related to the role that noise (either extrinsic or intrinsic) may play in the steady state characteristics of a bacterial population. We have not dealt with the use of simulation techniques per se in the study of these systems as they fall far from our purpose and are a subject of study in their own right.

Acknowledgments

This work was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the State Committee for Scientific Research (Poland) Grant N N201 608240, and the Consejo Nacional de Ciencia y Tecnología (Conacyt) in México.

References

1 

Jacob F, Perrin D, Sanchez C, Monod J1960Operon: A group of genes with the expression coordinated by anoperatorC R Hebd Seances Acad Sci25017271729

2 

Goodwin BC1965Oscillatory behaviour in enzymatic control processAdv Enzyme Regul3425438

3 

Goodwin BCTemporal organization in cells1963London-New YorkAcademic Press

4 

Griffith J1968Mathematics of cellular control processes. I. Negative feedback to one geneJ Theor Biol20202208

5 

Griffith J1968Mathematics of cellular control processes. II. Positive feedback to one geneJ Theor Biol20209216

6 

Tyson J, Othmer H1978The dynamics of feedback control circuits in biochemical pathwaysRosen RProgress inBiophysicsNew YorkAcademic Press162vol. 5

7 

Othmer H1976The qualitative dynamics of a class of biochemical control circuitsJ Math Biol35378

8 

Selgrade J1979Mathematical analysis of a cellular control process with positive feedbackSIAM J Appl Math36219229

9 

Mackey MC, Tyran-Kamińska M, Yvinec R2011Molecular distributions in gene regulatory dynamicsJTheor Biol2748496

10 

Smith H1995vol. 41, American MathematicalSocietyMonotone Dynamical Systems, Mathematical Surveys and MonographsProvidence, RI

11 

Mier-y-Teran-Romero L, Silber M, Hatzimanikatis V1000The origins of time-delay in template biopolymerizationprocessesPLOS Comp Biol6e726–1e726–15

12 

Heinrich R, Rapoport TA1980Mathematical modeling of translation of mRNA in eucaryotes: Steady states,time-dependent processes and application to reticulocytesJ Theor Biol86279313

13 

Ruan S, Wei J2001On the zeros of a third degree exponential polynomial with applications to a delayed modelfor the control of testosterone secretionIMA J Math Appl Med Biol184152

14 

Haken H1983Synergetics: An introduction, Springer Series in SynergeticsBerlinSpringer-Verlagvol. 1

15 

Stratonovich RL1963Revised English edition. Translated from the Russian by Richard A. Silverman, Gordon and Breach Science PublishersTopics in the theory of random noise. Vol. I: General theory of random processes. Nonlinear transformations of signals and noiseNew York

16 

Wilemski G1976On the derivation of Smoluchowski equations with corrections in the classical theory of BrownianmotionJ Stat Phys14153169

17 

Titular U1978A systematic solution procedure for the Fokker-Planck equation of a Brownian particle in thehigh-friction casePhysica A91A321344

18 

Gardiner C1983Handbook of Stochastic MethodsSpringer VerlagBerlin, Heidelberg

19 

Yvinec R, Zhuge C, Lei J, Mackey MC2014Adiabatic reduction of a model of stochastic gene expression withjump Markov processJ Math Biol6810511070

20 

Monod J1941Ph.D. thesis, Université de ParisRecherches sur la croissance des cultures bactériennesParis

21 

Lendenmann U, Snozzi M, Egli T1996Kinetics of the simultaneous utilization of sugar mixtures by Escherichia coli in continuous cultureAppl Environ Microbiol6214931499

22 

Serebriiskii IG, Golemis EA2000Uses of lacz to study gene function: Evaluation of β-galactosidaseassays employed in the yeast twohybrid systemAnal Biochem285115

23 

Abramson J, Iwata S, Kaback HR2004Lactose permease as a paradigm for membrane transport proteins (review)Mol Membr Biol21227236

24 

Kaback HR2005Structure and mechanism of the lactose permeaseC R Biol328557567

25 

Novick A, Weiner MEnzyme induction as an all-or-none phenomenonProc Natl Acad Sci U S A431957553566

26 

Cohn M, Horibata K1959Analysis of the differentiation and of the heterogeneity within a population of Escherichia coli undergoing induced beta-galactosidase synthesisJ Bacteriol78613623

27 

Babloyantz A, Sanglier M1972Chemical instabilities of “all-or-none” type in beta - galactosidase inductionand active transportFEBS Lett23364366

28 

Nicolis G, Prigogine I1977Self-organization in nonequilibrium systems. From dissipative structures toorder through fluctuationsJohn Wiley and SonsNew York

29 

Díaz-Hernández O, Santillán MBistable behavi or of the lacoperon in E. coli when inducedwith a mixture of lactose and TMGFront Physiol1

30 

Yildirim N, Mackey MC2003Feedback regulation in the lactose operon: A mathematical modeling study and comparison with experimental dataBiophys J8428412851

31 

Yildirim N, Santillán M, Horike D, Mackey MC2004Dynamics and bistability in a reduced model of the lacoperonChaos14279292

32 

Santillán M, Mackey MC2004Inuence of catabolite repression and inducer exclusion on the bistable behorof the lac operonBiophys J8612821292avi

33 

Santillán M, Mackey MC, Zeron ESOrigin of bistability in the lac operonBiophys J92200738303842

34 

Santillán M2008Bistable behavior in a model of the lac operon in Escherichia coli with variable growthrateBiophys J920652081

35 

Santillán M, Mackey MC2008Quantitative approaches to the study of bistability in the lac operon of Escherichia coliJ R Soc Interface5S29S39

36 

Reznikoff WS1992The lactose operon-controlling elements: A complex paradigmMol Microbiol624192422

37 

Müller-Hill B1998The function of auxiliaty operatorsMol Microbiol291318

38 

Lewis M2005The lac repressorC R Biol328521548

39 

Wilson CJ, Zhan H, Swint-Kruse L, Matthews KS2007The lactose repressor system, paradigms for regulation,allosteric behavior and protein foldingCell Mol Life Sci64316

40 

Narang A2007Effect of DNA looping on the induction kinetics of the lacoperonJ Theor Biol247695712

41 

Oehler S, Eismann ER, Krämer H, Müller-Hill B1990The three operators of lac operon cooperate inrepressionEMBO J9973979

42 

Kuhlman T, Zhang Z, Saier MH, Hwa T2007Combinatorial transcriptional control of the lactose operon of Escherichia coliProc Natl Acad Sci U S A10460436048

43 

Beckwith J 1987 Escherichia coli and Salmonella thyphymurium: Cellular and molecularbiology Neidhart F, Ingraham J, Low K, Magasanik B, Umbarger H The lactose operon American Society for Microbiology Washington, DC 1439 1443 vol. 2

44 

Ozbudak E, Thattai M, Lim H, Shraiman B, van A2004Oudenaarden, Multistability in the lactose utilizationnetwork of Escherichia coliNature427737740

45 

Santillán M, Zeron ES2004Dynamic inuence of feedback enzyme inhibition and transcription attenuation onthe tryptophan operon response to nutritional shiftsJ Theor Biol231287298

46 

Santillán M2008On the use of the Hill functions in mathematical models of gene regulatory networksMathModel Nat Phenom38597

47 

Tyson JJ1975On the existence of oscillatory solutions in negative feedback cellular control processesJTheor Biol1311315

48 

Bliss RD, Painter PR, Marr AG1982Role of feedback inhibition in stabilizing the classical operonJTheor Biol97177193

49 

Sinha S1988Complex behaviour of the repressible operonJ Theor Biol132307318

50 

Sen AK, Liu W-M1990Dynamic analysis of genetic control and regulation of amino acid synthesis: Thetryptophan operon in Escherichia coliBiotechnol Bioeng35185194

51 

Giona M, Adrover A2002Modified model for the regulation of the tryptophan operon in Escherichia coliBiotechnol Bioeng80297304

52 

Santillán M, Mackey MC2001Dynamic regulation of the tryptophan operon: Modeling study and comparison withexperimental dataProc Natl Acad Sci U S A9813641369

53 

Santillán M, Zeron ES2006Analytical study of the multiplicity of regulatory mechanisms in the tryptophanoperonBull Math Biol68343359

54 

Hernández-Valdez A, Santillán M2010Cyclic expression and cooperative operator interaction in the trpoperon of Escherichia coliJ Theor Biol263340352

55 

Salazar-Cavazos E, Santillán M2013Optimal performance of the tryptophan operon of E. coli: Astochastic, dynamical, mathematical modeling approachBull Math Biol76314334

56 

Yanofsky C, Horn V1994Role of regulatory features of the trp operon of Escherichia coli in mediating aresponse to a nutritional shiftJ Bacteriol17662456254

57 

Zamora-Chimal C, Santillán M, Rodríguez-González J2012Inuence of the feedback loops in the trpoperon of B. subtilis on the system dynamic response and noise amplitudeJ Theor Biol310119131

58 

Kaern M, Elston T, Blake W, Collins J2005Stochasticity in gene expression: From theories to phenotypesNat Rev Genet6451464

59 

Raj A, van A2008Oudenaarden, Nature, nurture, or chance: Stochastic gene expression and its consequencesCell135216226

60 

Shahrezaei V, Swain PThe stochastic nature of biochemical networksCurr Opin Biotechnol192008369374

61 

Shahrezaei V, Ollivier J, Swain P2008Colored extrinsic uctuations and stochastic gene expressionMolSyst Biol4196205

62 

Ochab-Marcinek APredicting the asymmetric response of a genetic switch to noiseJ Theor Biol25420083744

63 

Ochab-Marcinek A2010Extrinsic noise passing through a Michaelis-Menten reaction: A universal response of a geneticswitchJ Theor Biol263510520

64 

Ochab-Marcinek A, Tabaka M2010Bimodal gene expression in noncooperative regulatory systemsProc NatlAcad Sci U S A1072209622101

65 

Elowitz M, Levine A, Siggia E, Swain P2002Stochastic gene expression in a single cellScience29711831186

66 

Raser J, O’Shea E2004Control of stochasticity in eukaryotic gene expressionScience30418111814

67 

Swain P, Elowitz M, Siggia E2002Intrinsic and extrinsic contributions to stochasticity in gene expressionProc Natl Acad Sci U S A991279512800

68 

Scott M, Ingallls B, Kærn M2006Estimations of intrinsic and extrinsic noise in models of nonlineargenetic networksChaos16026107–1026107–15

69 

Blake W, Balázsi G, Kohanski M, Issacs F, Murphy K, Kuang Y, Cantor C, Walt D, CollinsPhenotypic J2006consequences of promoter-mediated transcriptional noiseMol Cell24853865

70 

Fraser H, Hirsh A, Glaever G, Kumm J, Eisen M2004Noise minimization in eukaryotic gene expressionPLoS Biol28343838

71 

Berg OG1978A model for the statistical uctuations of protein numbers in a microbial populationJ TheorBiol71587603

72 

Blake W, Kaern M, Cantor C, Collins J2003Noise in eukaryotic gene expressionNature422633637

73 

Cai L, Friedman N, Xie X2006Stochastic protein expression in individual cells at the single molecule levelNature440358362

74 

Chubb J, Trcek T, Shenoy S, Singer R2006Transcriptional pulsing of a developmental geneCurr Biol1610181025

75 

Golding I, Paulsson J, Zawilski S, Cox E2005Real-time kinetics of gene activity in individual bacteriaCell12310251036

76 

Raj A, Peskin C, Tranchina D, Vargas D, Tyagi S2006Stochastic mRNA synthesis in mammalian cellsPLoSBiol417071719

77 

Sigal A, Milo R, Cohen A, Geva-Zatorsky N, Klein Y, Liron Y, Rosenfeld N, Danon T, Perzov N, AlonVariability U2006and memory of protein levels in human cellsNature444643646

78 

Suter D, Molina N, Gatfield D, Schneider K, Schibler U, Naef F2011Mammalian genes are transcribed withwidely different bursting kineticsScience332472474

79 

Yu J, Xiao J, Ren X, Lao K, Xie X2006Probing gene expression in live cells, one protein molecule at a timeScience31116001603

80 

Kepler T, Elston T2001Stochasticity in transcriptional regulation: Origins, consequences, and mathematicalrepresentationsBiophy J8131163136

81 

Friedman N, Cai L, Xie X2006Linking stochastic dynamics to population distribution: An analytical frameworkof gene expressionPhys Rev Lett97168302–1/4

82 

Morelli L, Julicher F2007Precision of genetic oscillators and clocksPhys Rev Lett9228101

83 

Bobrowski A, Lipniacki T, Pichór K, Rudnicki R2007Asymptotic behavior of distributions of mRNA andprotein levels in a model of stochastic gene expressionJ Math Anal Appl333753769

84 

Shahrezaei V, Swain P2008Analytic distributions for stochastic gene expressionProc Natl Acad Sci U S A1051725617261

85 

Iyer-Biswas S, Hayot F, Jayaprakash C2009Stochasticity of gene products from transcriptional pulsingPhys Rev E7031911

86 

Mugler A, Walczak A, Wiggins C2009Spectral solutions to stochastic models of gene expression with bursts andregulationPhys Rev E8041921

87 

Ribeiro A, Smolander O, Rajala T, Hakkinen A, Yli-Harja O2009Delayed stochastic model of transcription atthe single nucleotide levelJ Comput Biol1653953

88 

Elgart V, Jia T, Kulkarni R2010Applications of Little’s law to stochastic models of gene expressionPhys Rev E8021901

89 

Lei J2010Stochasticity in single gene expression with both intrinsic noise and uctuation in kinetic parametersErratum appears in J Theor Biol256485492Erratum appears in J Theor Biol 262(2) (2010), 381

90 

Rajala T, Hakkinen A, Healy S, Yli-Harja O, Ribeiro A2010Effects of transcriptional pausing on geneexpression dynamicsPLoS Comput Biol6e1000704

91 

Tang M2010The mean frequency of transcriptional bursting and its variation in single cellsJ Math Biol602758

92 

Bett G, Zhou Q, Rasmusson R2011Models of HERG gatingBiophys J101631642

93 

Jia T, Kulkarni R2011Intrinsic noise in stochastic models of gene expression with molecular memory andburstingPhys Rev Lett10058102

94 

Cottrell D, Swain P, Tupper P2012Stochastic branching-diffusion models for gene expressionProc NatlAcad Sci U S A10996999704

95 

Feng H, Hensel Z, Xiao J, Wang J2012Analytical calculation of protein production distributions in models ofclustered protein expressionPhys Rev E85155160

96 

Ferguson M, Le Coq D, Jules M, Aymerich S, Radulescu O, Declerck N, Royer C2012Reconciling molecularregulatory mechanisms with noise patterns of bacterial metabolic promoters in induced and repressed statesProc Natl Acad Sci U S A109155160

97 

Kuwahara H, Schwartz R2012Stochastic steady state gain in a gene expression process with mRNA degradationcontrolJ R Soc Interface915891598

98 

Singh A, Bokes P2012Consequences of mRNA transport on stochastic variability in protein levelsBiophysJ10310871096

99 

Earnest T, Roberts E, Assaf M, Dahmen K, Luthey-Schulten Z2013DNA looping increases the range ofbistability in a stochastic model of the lac genetic switchPhys Biol1026002

100 

Robinson R2013Bursting with randomness: A simple model for stochastic control of gene expressionPLoS Biol11e1001622

101 

Tian T2013Chemical memory reactions induced bursting dynamics in gene expressionPLoS One8e52029

102 

Mackey MC, Tyran-Kamińska M2008Dynamics and density evolution in piecewise deterministic growthprocessesAnn Polon Math94111129

103 

Mackey MC, Tyran-Kamińska M, Yvinec R2013Dynamic behavior of stochastic gene expression models in thepresence of burstingSIAM J Appl Math7318301852

104 

Oppenheim I, Schuler K, Weiss G1969Stochastic and deterministic formulation of chemical rate equationsJ Chem Phys50460466

105 

Lasota A, Mackey M.C1994Chaos, fractals, and noise, Applied Mathematical SciencesNew YorkSpringer-Verlagvol. 97

106 

Pichór K, Rudnicki R2000Continuous Markov semigroups and stability of transport equationsJ MathAnal Appl249668685

Figures and Tables

Fig.2.1

A schematic illustration of the possibility of one, two or three solutions of Equation (2.14) for varying values of κ d in the presence of inducible regulation. The monotone increasing graph is f of Equation (2.10), and the straight lines correspond to x/κ d for (in a clockwise direction) κ d  ∈ [0, κ d-), κ d  = κ d-, κ d  ∈ (κ d-, κ d+), κ d  = κ d+, and κ d+ < κ d . This figure was constructed with n = 4 and K = 10 for which κ d- = 3.01 and κ d+ = 5.91 as computed from (2.17). See the text for details. Taken from [9] with permission.

A schematic illustration of the
possibility of one, two or three solutions of Equation (2.14) for varying values of κ

d
 in the presence of
inducible regulation. The monotone increasing graph is f of Equation (2.10), and the straight lines correspond
to x/κ

d
 for (in a clockwise direction) κ

d
 ∈ [0, κ

d-), κ

d
 = κ

d-,
κ

d
 ∈ (κ

d-, κ

d+), κ

d
 = κ

d+, and κ

d+ < κ

d
. This figure was
constructed with n = 4 and K = 10 for which κ

d- = 3.01 and κ

d+ = 5.91 as computed from
(2.17). See the text for details. Taken from [9] with permission.
Fig.2.2

Logarithmic plot of the steady state values of x * versus κ d for an inducible operon obtained from Equation (2.14), for n = 4 and K = 2, 5, 10, and 15 (left to right) illustrating the dependence of the occurrence of bistability on K. See the text for details. Taken from [9] with permission.

Logarithmic plot of the steady state
values of x
* versus κ

d
 for an inducible operon obtained from Equation (2.14), for n = 4 and K = 2, 5, 10, and 15 (left to right) illustrating the dependence of the occurrence of bistability on K. See the text for details. Taken from [9] with permission.
Fig.2.3

This figure presents a parametric plot (for n = 4) of the bifurcation diagram in (K, κ d ) parameter space separating one from three steady states in an inducible operon as determined from Equations (2.14) and (2.15). The upper (lower) branch corresponds to κ d- (κ d+), and for all values of (K, κ d ) in the interior of the cone there are two locally stable steady states X1*,X3*, while outside there is only one. The tip of the cone occurs at (K,κd)=((5/3)2,(5/3)5/34) as given by Equations (2.18) and (2.19). For K ∈ [0, (5/3) 2) there is a single steady state. Taken from [9] with permission.

This figure presents a parametric plot
(for n = 4) of the bifurcation diagram in (K, κ

d
) parameter space separating one from three steady states
in an inducible operon as determined from Equations (2.14) and (2.15). The upper (lower) branch corresponds to
κ

d- (κ

d+), and for all values of (K, κ

d
) in the interior of the cone there are two
locally stable steady states X1*,X3*, while outside there is only one. The tip of the cone occurs at
(K,κd)=((5/3)2,(5/3)5/34) as given by Equations (2.18) and
(2.19). For K ∈ [0, (5/3) 2) there is a single steady state. Taken from [9] with permission.
Fig.3.1

Typical diauxic growth curve. Note the existence of two different exponential growth phases, separated by a short interval in which the culture does not grow. The first (second) phase corresponds to the bacterial culture feeding on glucose (lactose), while the interval with no growth corresponds to the time required for the bacteria to turn on the genes needed to metabolize lactose after glucose exhaustion.

Typical diauxic growth curve. Note the
existence of two different exponential growth phases, separated by a short interval in which the culture does not
grow. The first (second) phase corresponds to the bacterial culture feeding on glucose (lactose), while the
interval with no growth corresponds to the time required for
the bacteria to turn on the genes needed to metabolize lactose after glucose exhaustion.
Fig.3.2

Cartoon representation of the lactose operon regulatory pathway. 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 (inhibit) the process they point to. Taken from [29] with permission.

Cartoon representation of the lactose
operon regulatory pathway. 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 (inhibit) the process they point to. Taken from [29] with permission.
Fig.3.3

(A) Schematic of the regulatory elements located in lactose 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 panels B to E. Namely, a free repressor molecule (B), one with a single subunit bound by allolactose (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). Taken from [29] with permission.

(A) Schematic of the regulatory elements
located in lactose 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 panels B to E. Namely, a free
repressor molecule (B), one with a single subunit bound by allolactose (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). Taken from [29] with permission.
Fig.3.4

Bifurcation diagrams in the L e versus G e parameter space, calculated with the present model under various conditions. In both cases, the bistability region is shaded, while the monostable induced (top) and uninduced (bottom) regions are uncoloured. In the left graph we used the parameter values tabulated in Table 3.1, and set K A  = 8.2 × 105 mpb. In the right graph all parameters ξ j and ξj* were reduced to 0.15 times the value reported in Table 3.1, and K A was set to 2.8 × 106 mpb. To calculate the graphs, we set φ M  = 0/min to simulate induction of the lac operon with the non-metabolizable TMG. The K A values referred to above were chosen to fit the experimental results of [44], which are shown with solid diamonds. Taken from [34] with permision.

Bifurcation diagrams in the L

e
 versus
G

e
 parameter space, calculated with the present model under various conditions. In both cases, the bistability
region is shaded, while the monostable induced (top) and uninduced (bottom) regions are uncoloured. In the left
graph we used the parameter values tabulated in Table 3.1, and set K

A
 = 8.2 × 105
mpb. In
the right graph all parameters ξ

j
 and ξj* were reduced to 0.15 times the value reported in
Table 3.1, and K

A
 was set to 2.8 × 106
mpb. To calculate the graphs, we set
φ

M
 = 0/min to simulate induction of the lac operon with the non-metabolizable TMG. The K

A
 values referred
to above were chosen to fit the experimental results of
[44], which are shown with solid diamonds. Taken from [34] with permision.
Fig.3.5

Schematic representation of the tryptophan operon regulatory pathway. Solid lines represent the processes involved in gene expression and tryptophan synthesis, while dashed lines correspond to the operon regulatory mechanisms.

Schematic representation of the tryptophan
operon regulatory pathway. Solid lines represent the processes involved in gene expression and tryptophan
synthesis, while dashed lines correspond to the operon
regulatory mechanisms.
Fig.3.6

Oscillatory behaviour during a de-repression experiment with an enzyme-inhibition-less E. coli mutant strain, as predicted by the deterministic model in [54].

Oscillatory behaviour during a
de-repression experiment with an enzyme-inhibition-less E. coli
mutant strain, as predicted by the deterministic model in [54].
Fig.3.7

(A) Stochastic quasi-oscillatory behaviour observed during a de-repression experiment with an enzyme-inhibition-less E. coli mutant strain, as predicted by the stochastic model in [54]. (B) Average behaviour of 100 independent cells.

(A) Stochastic quasi-oscillatory
behaviour observed during a de-repression experiment with an enzyme-inhibition-less E. coli mutant strain,
as predicted by the stochastic model in
[54]. (B) Average behaviour of 100 independent cells.
Fig.3.8

Comparison of experimental results (circles and triangles) and model predictions for de-repression experiments carried out with a wild-type and an enzyme-inhibition-less mutant strain of E. coli.

Comparison of experimental results
(circles and triangles) and model predictions for de-repression experiments carried out with a wild-type and an enzyme-inhibition-less mutant strain of E. coli.
Fig.4.1

Schematic illustration that there can be one or no solution of Equation (4.5), depending on the value of κ b , with repressible regulation. The straight lines correspond (in a clockwise direction) to κ b  = 2 and κ b  = 0.8. This figure was constructed with n = 4, Δ = 10 and b = 1. See the text for further details. Taken from [9] with permission.

Schematic illustration that there can be
one or no solution of Equation (4.5), depending on the value of κ

b
, with repressible regulation. The
straight lines correspond (in a clockwise direction) to κ

b
 = 2 and κ

b
 = 0.8. This figure was
constructed with n = 4, Δ = 10 and b = 1. See the text for further details. Taken from [9] with permission.
Fig.4.2

Schematic illustration of the possibility of one, two or three solutions of equation (4.6) for varying values of κ b with bursting inducible regulation. The straight lines correspond (in a clockwise direction) to κ b  ∈ (0, κ b-), κ b  = κ b-, κ b  ∈ (κ b-, κ b+) (and respectively κ b  < K, κ b  = K, K < κ b ), κ b  = κ b+, and κ b+ < κ b . This figure was constructed with n = 4, K = 10 and b = 1 for which κ b- = 4.29 and κ b+ = 14.35 as computed from (4.9). See the text for further details. Taken from [9] with permission.

Schematic illustration of the possibility
of one, two or three solutions of equation (4.6) for varying values of κ

b
 with bursting
inducible regulation. The straight lines correspond (in a clockwise direction) to κ

b
 ∈ (0, κ

b-),
κ

b
 = κ

b-, κ

b
 ∈ (κ

b-, κ

b+) (and respectively κ

b
 < K, κ

b
 = K,
K < κ

b
), κ

b
 = κ

b+, and κ

b+ < κ

b
. This figure was constructed with n = 4,
K = 10 and b = 1 for which κ

b- = 4.29 and κ

b+ = 14.35 as computed from (4.9). See the text for further details. Taken from
[9] with permission.
Fig.4.3

In this figure we illustrate Bifurcation type 1 when intrinsic bursting is present. The stationary density u * is plotted versus x between 0 and 6. The values of the parameters used in this figure are b = 1, K = 4, and n = 4. The parameter κ b was taken to be 3, 3.7, 4, 4.5, 5, 6.5, where for κ b  = 3 we have unimodal type 1 density, and increasing κ b we obtain bimodal type 2 density for κ b  = 4.5, 5.

In this figure we illustrate Bifurcation type 1 when intrinsic bursting is present. The stationary density u
* is plotted versus x
between 0 and 6. The values of the parameters used in this figure are b = 1, K = 4, and n = 4. The parameter
κ

b
 was taken to be 3, 3.7, 4, 4.5, 5, 6.5, where for
 κ

b
 = 3 we have unimodal type 1 density, and increasing κ

b
 we obtain bimodal type 2
density for κ

b
 = 4.5, 5.
Fig.4.4

An illustration of Bifurcation type 2 for intrinsic bursting. The stationary density u * is plotted versus x between 0 and 5. The parameters used are b = 1, K = 4, and n = 4. The parameter κ b was taken to be 2.5, 22, 25, 28, 31, where we have Unimodal type 1 density for κ b  = 2.5, then Unimodal type 2 for κ b  = 22, Bimodal type 2 for κ b  = 25, 28, and back Unimodal type 2.

An illustration of Bifurcation type
2 for intrinsic bursting. The stationary density u
* is plotted versus x between 0 and 5. The parameters
used are b = 1, K = 4, and n = 4. The parameter κ

b
 was taken to be 2.5, 22, 25, 28, 31, where we have
Unimodal type 1 density for κ

b
 = 2.5, then Unimodal type 2 for κ

b
 = 22, Bimodal type 2 for
 κ

b
 = 25, 28, and back Unimodal type 2.
Fig.4.5

The presence of bursting can drastically alter regions of bimodal behaviour as shown in this parametric plot (for n = 4) of the boundary in (K, b · κ b ) parameter space delineating unimodal from bimodal stationary densities u * in an inducible operon with bursting and in (K, κ d ) parameter space delineating one from three steady states in the deterministic inducible operon. From top to bottom, the regions are for b = 10, b = 1, b = 0.1 and b = 0.01. The lowest (heavy dashed line) is for the deterministic case. Note that for b = 0.1, the two regions of bistability and bimodality coincide and are indistinguishable from one another. Taken from [9] with permission.

The presence of bursting can drastically
alter regions of bimodal behaviour as shown in this parametric plot (for n = 4) of the boundary in (K, b · κ

b
) parameter space delineating unimodal from bimodal stationary densities u
* in an inducible operon
with bursting and in (K, κ

d
) parameter space delineating one from three steady states in the deterministic
inducible operon. From top to bottom, the regions are for b = 10, b = 1, b = 0.1 and b = 0.01. The lowest (heavy
dashed line) is for the deterministic case. Note that for b = 0.1, the two regions of bistability and
bimodality coincide and are indistinguishable from one another. Taken from [9] with permission.
Table 2.1

The parameters A, B, Λ, Δ and θ for the inducible and repressible cases. See the text and Section 2.1 for more detail

ParameterInducibleRepressible
A K = 1 + K 2 R tot 1
B K 1 K = K 1 (1 + K 2 R tot )
K1K K
Λ = A K 1
Δ=BK1-1 1 KK1-1
θ=κdnΔ(1-ΔΛ) κdn·K-1K>0 κdn·K1-KK<0
Table 3.1

Values of the parameters in the lactose operon model equations. The parameter K A is the only one that we were unable to estimate. The unit mpb stands for molecules per average-size bacterium

D ≈ 2mpb k M  ≈ 0.18/min
μ ≈ 0.02/min γ M  ≈ 0.46/min
p p  ≈ 0.127 k pc  ≈ 30
K G  ≈ 2.6μM m ≈ 1.3
ξ 1 ≈ 17 ξ 2 ≈ 0.85
ξ 3 ≈ 0.17 ξ1*0
ξ2*430.6 ξ3*1261.7
K A  ≥ 0 mpb k Z  ≈ 18.8/min
γ Z  ≈ 0.01/min τ Z  ≈ 0.1min
τZ*0.42min k L  ≈ 6.0 × 104/min
k  ≈ 0/min κ  ≈ 680μM
φ M  ∈ [0, 3.8 × 104]/min κ L  ≈ 680μM
φ G  ≈ 0.35 κ G  ≈ 1.0μM
κ M  ≈ 7.0 × 105mpb