The utility of simple mathematical models in understanding gene regulatory dynamics

2.1 Kinetic 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.

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]dM dt = b ¯ d φ ¯ m f ( E ) - γ M M ,(2.7) dI dt = β I M - γ I I ,(2.8) dE dt = β E I - γ E E .(2.9) 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 ) = φ m f ( e ) , where φ _m (which is dimensionless) is defined byφ m = φ ¯ m β E β I γ M γ E γ I and f ( e ) = 1 + e n Λ + Δ e n ,(2.10) Λ and Δ are defined in Table 2.1, and a (dimensionless) effector concentration (e) is defined by E = η e with η = 1 K 1 n . We continue and define dimensionless intermediate protein (i) and mRNA concentrations (m): I = i η γ E β E and M = m η γ E γ I β E β I , so Equations (2.7)–(2.9) take the form dm dt = γ M [ κ d f ( e ) - m ] , di dt = γ I ( m - i ) , de dt = γ E ( i - e ) , with the dimensionless constants κ d = b d φ m and b d = b ¯ d η .

To continue our simplifications we rename the dimensionless concentrations through (m, i, e) = (x ₁, x ₂, x ₃), and subscripts (M, I, E) = (1, 2, 3) to finally obtain dx 1 dt = γ 1 [ κ d f ( x 3 ) - x 1 ] ,(2.11) dx 2 dt = γ 2 ( x 1 - x 2 ) ,(2.12) dx 3 dt = γ 3 ( x 2 - x 3 ) .(2.13) 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 ₁, x ₂, x ₃) 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 ₀ (X) = X. For both inducible and repressible operons, for all initial conditions X 0 = ( x 1 0 , x 2 0 , x 3 0 ) ∈ ℝ 3 + the flow S t ( X 0 ) ∈ ℝ 3 + for t > 0.

The steady state solutions of (2.11)– (2.13) are given by the solutions ofx κ d = f ( x )(2.14) and for each solution x ^* of Equation (2.14) there is a steady state X * = ( x 1 * , x 2 * , x 3 * ) of (2.11)– (2.13) which is given by x 1 * = x 2 * = x 3 * = 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.3 The 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 dM dt = b ¯ d φ ¯ m f ( e - μ τ M E τ M ) - γ ¯ M M ,(2.20) dI dt = β I e - μ τ I M τ I - γ ¯ I I ,(2.21) dE dt = β E I - γ ¯ E E .(2.22) 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 = η e with η = 1 e - μ τ M K 1 n , I = i η γ E β E and M = m η γ E γ I β E β I e - μ τ I so Equations (2.20)–(2.22) take the form dm dt = γ ¯ M [ κ d f ( e τ M ) - m ] , di dt = γ ¯ I ( m τ I - i ) , de dt = γ ¯ E ( i - e ) , withκ d = b ¯ d φ m β I β E γ ¯ M γ ¯ I γ ¯ E η e - μ τ I .(2.23)

To finish our simplifications, as before rename the dimensionless concentrations (m, i, e) = (x ₁, x ₂, x ₃), and subscripts (M, I, E) = (1, 2, 3) to obtain dx 1 dt = γ ¯ 1 [ κ d f ( x 3 , τ 1 ) - x 1 ] ,(2.24) dx 2 dt = γ ¯ 2 ( x 1 , τ 2 - x 2 ) ,(2.25) dx 3 dt = γ ¯ 3 ( x 2 - x 3 ) .(2.26) 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 τ ₁ and τ ₂) 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 (τ ₁, τ ₂).

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

P ( λ ) = ( γ ¯ 1 + λ ) ( γ ¯ 2 + λ ) ( γ ¯ 3 + λ ) and ϑ = - κ d f ′ ( X * ) γ ¯ 1 γ ¯ 2 γ ¯ 3 > 0(2.27) and τ = τ ₁ + τ ₂. Writing out g (λ) we haveg ( λ ) = λ 3 + a 1 λ 2 + a 2 λ + a 3 + ϑ e - λ τ ,(2.28) where a 1 = ∑ i = 1 3 γ i , a 2 = ∑ i ≠ j = 1 3 γ i γ j , a 3 = ∏ i = 1 3 γ i .

Let λ (τ) = α (τ) + iω (τ) be the root of Equation (2.28) satisfying α (τ ₀) =0 and ω (τ ₀) = ω ₀, and set p = a 1 2 - 2 a 2, q = a 2 2 - 2 a 1 a 2, r = a 3 2 - ϑ 2, and let h (z) = z ³ + pz ² + 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 ² - 3q < 0, then all roots of Equation (2.28) have negative real parts for all τ ≥ 0.

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 (γ ₁ ⪢ γ ₂, γ ₃) so in this case the mRNA dynamics are fast and we have the approximate relationship x 1 ≃ κ d f ( x 3 ) . If γ ₁ ⪢ γ ₂ ⪢ γ ₃ so that the effector is the slowest variable, then we have x 2 ≃ x 3 and the three variable system (2.11)– (2.13) describing the generic operon reduces to a one dimensional systemdx 3 dt = γ 3 [ κ d f ( x 3 ) - x 3 ](2.29) for the relatively slow effector dynamics. If instead the effector qualifies as a fast variable (as for the lac operon) so that γ ₁ ⪢ γ ₃ ⪢ γ ₂ and x 3 ≃ x 2 then the intermediate protein is the slowest variable described by the one-dimensional equationdx 2 dt = γ 2 [ κ d f ( x 2 ) - x 2 ] .(2.30) Consequently both Equations (2.30) and (2.29) are of the formdx dt = γ [ κ d f ( x ) - x ] ,(2.31) where γ is either γ ₂ for protein (x ₂) dominated dynamics or γ ₃ for effector (x ₃) 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].

3.1 The 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 dM Z dt = Dk M P τ Z ( O P ) - ( γ M + μ ) M Z ,(3.1) dM Y dt = Dk M P τ Y ( O P ) - ( γ M + μ ) M Y .(3.2)

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 dE Z dt = k Z e - μ τ Z * M Z , τ Z * - ( γ Z + μ ) E Z ,(3.3) dE Y dt = k Y e - μ τ Y * M Y , τ Y * - ( γ Y + μ ) E Y .(3.4)

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 M Z , τ Z * ( t ) stands for the delayed function M Z ( 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 ( k Z , τ Z * , γ Z ) is identically equal to ( k Y , τ 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 takeB = E Z / 4 and Q = E Y ,(3.5) 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

dL dt = k L β L ( L e ) β G ( G e ) Q - k ℓ β ℓ ( L ) Q - φ M M ( L ) B - ( γ L + μ ) L .(3.6) 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 _ℓ β _ℓ (L) Q expresses the loss of intracellular lactose to the extracellular fluid due to the reversible nature of the permease mediated transport; the third term φ M M ( 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:β L ( L e ) = L e κ L + L e and β G ( G e ) = 1 - φ G G e κ G + G e .(3.7) Furthermore, the terms β _ℓ (L) and ℳ ( L ) are both functions of the internal lactoseβ ℓ ( L ) = L κ ℓ + L and M ( L ) = L κ M + L .(3.8)

The dynamical equation for the concentration of allolactose A is much simpler:dA dt = α φ M M ( L ) B - φ A M ( A ) B - ( γ A + μ ) A ,(3.9) where α is the fraction of internal lactose L transformed by β-galactosidase B into allolactose instead of being hydrolyzed into glucose and galactose. The term φ A M ( 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 ( O P ) = a + A n b + A n 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 P ( O P ) = p pc ( G e ) P R ( A ) ,(3.10) p pc ( G e ) = p p 1 + ( k pc - 1 ) p c ( G e ) 1 + ( k pc - 1 ) p p p c ( G e ) ,(3.11) p cp ( G e ) = p c ( G e ) 1 + ( k pc - 1 ) p p 1 + ( k pc - 1 ) p p p c ( G e ) ,(3.12) p c ( G e ) = K G m K G m + G e m ,(3.13) P R ( A ) = ( 1 + ξ 2 ρ ( A ) ) ( 1 + ξ 3 ρ ( A ) ) + ξ 1 * ρ ( A ) 2 Z ( A ) + ∏ j = 1 , 2 , 3 ( 1 + ξ jitsc ρ ( A ) ) ,(3.14) Z ( A ) = ∑ j = 1 , 2 , 3 p cp ( G e ) δ 2 j ( 1 + ξ j ρ ( A ) ) ξ j * ρ ( A ) 2 ,(3.15) ρ ( A ) = ( K A K A + A ) 2 .(3.16) 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 ρ (A) ², 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 , k Z , τ Z * , γ Z ) is equal to ( τ Y , k Y , τ 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 dM Z dt = Dk M p pc ( G e ) P R ( A ) - ( γ M + μ ) M Z ,(3.17) dE Z dt = k Z e - μ τ Z * M Z , τ Z * - ( γ Z + μ ) E Z ,(3.18) dL dt = k L β L ( L e ) β G ( G e ) Q - k ℓ β ℓ ( L ) Q - φ M M ( L ) B - ( γ L + μ ) L .(3.19)

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), B = E Z / 4 ,(3.20) Q = E Z ,(3.21) A = L ,(3.22) β L ( L e ) = L e κ L + L e ,(3.23) β G ( G e ) = 1 - φ G G e κ G + G e ,(3.24) β ℓ ( L ) = L κ ℓ + L ,(3.25) M ( L ) = L κ M + L .(3.26)

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 × 10⁵ 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 × 10⁶ mpb. The results are shown in Fig. 3.4B.

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: P 000 ⇌ β 1 α 1 P 100 , P 000 ⇌ β 2 α 2 P 010 , P 000 ⇌ β 3 α 3 P 001 , P 100 ⇌ β 2 / k C α 2 P 110 , P 100 ⇌ β 3 α 3 P 101 , P 010 ⇌ β 1 / k C α 1 P 110 , P 010 ⇌ β 3 α 3 P 011 , P 001 ⇌ β 1 α 1 P 101 , P 001 ⇌ β 2 α 2 P 011 , P 110 ⇌ β 3 α 3 P 111 , P 101 ⇌ β 2 / k C α 2 P 111 , P 011 ⇌ β 1 / k C α 1 P 111 , 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 : d P 000 dt = - ( α 1 + α 2 + α 3 ) P 000 + β 1 P 100 + β 2 P 010 + β 3 P 001 ,(3.27) d P 100 dt = - ( β 1 + α 2 + α 3 ) P 100 + α 1 P 000 + β 2 k C P 110 + β 3 P 101 ,(3.28) d P 010 dt = - ( α 1 + β 2 + α 3 ) P 010 + β 1 k C P 110 + α 2 P 000 + β 3 P 011 ,(3.29) d P 001 dt = - ( α 1 + α 2 + β 3 ) P 001 + β 1 P 101 + β 2 P 011 + α 3 P 000 ,(3.30) d P 110 dt = - ( β 1 k C + β 2 k C + α 3 ) P 110 + α 1 P 010 + α 2 P 100 + β 3 P 111 ,(3.31) d P 101 dt = - ( β 1 + α 2 + β 3 ) P 101 + α 1 P 001 + β 2 k C P 111 + α 3 P 100 ,(3.32) d P 011 dt = - ( α 1 + β 2 + β 3 ) P 011 + β 1 k C P 111 + α 2 P 001 + α 3 P 010 ,(3.33) d P 111 dt = - ( β 1 k C + β 2 k C + β 3 ) P 111 + α 1 P 011 + α 2 P 101 + α 3 P 110 .(3.34) 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: dM dt = k M P 000 A ( T ) - γ M M ,(3.35) dE dt = k E M - γ E E ,(3.36) dT dt = k T E I ( T ) - γ T T T + K T ,(3.37) 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,α i = a i R A ( T ) .(3.38) 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: R A ( T ) = R Tot ( T T + K A ) 2 ,(3.39) A ( T ) = 1 + 2 α T K G + T ( 1 + α T K G + T ) 2 ,(3.40) I ( T ) = K I n K I n + T n .(3.41) 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 tRNA^Trp 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 dM dt = k M P 000 ( T ) A ( T ) - γ M M ,(3.42) dE dt = k E M - γ E E ,(3.43) dT dt = k T E I ( T ) - γ T T T + K T .(3.44) The concentration of non-repressed promoters is given in this case by

P 000 ( T ) = ( 1 + a 1 β 1 R A ( T ) + a 2 β 2 R A ( T ) + a 3 β 3 R A ( T ) + k C a 1 β 1 a 2 β 2 R A 2 ( T ) + a 1 β 1 a 3 β 3 R A 2 ( T ) + a 2 β 2 a 3 β 3 R A 2 ( T ) + k C a 1 β 1 a 2 β 2 a 3 β 3 R A 3 ( T ) ) - 1 .(3.45) P ₀₀₀ (T), 𝒜 ( T ), and ℐ ( T ) are monotonic sigmoidally decreasing functions of T, and so is the product P 000 ( 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.

4.1 Dynamics with bursting