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.
Introduction
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.
Generic 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’.
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.
Inducible 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
Assume that the maximal transcription rate of DNA (in units of time - 1) is . 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 , or
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
Both (2.3) and (2.5) are special cases of
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 (in concentration units), the dynamics are given by [2, 4, 5, 7, 8]
To simplify things we formulate Equations (2.7)–(2.9) using dimensionless concentrations. To start we rewrite Equation (2.6) in the form
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
The steady state solutions of (2.11)– (2.13) are given by the solutions of
No control simply means f (x) ≡1, and in this case there is a single steady state x * = κ d that is globally asymptotically stable.
Inducible regulation
Single versus multiple steady states. For an inducible operon (with f given by Equation (2.2)) there may be one ( or ), two ( or ), or three () steady states, with the ordering , corresponding to the possible solutions of Equation (2.14) (cf. Figure 2.1). The smallest steady state is typically called the un-induced state, while the largest steady state 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
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:
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.
such that the flow St is directed inward everywhere on the surface of BI. Furthermore, all X* ∈ BI and
If there is a single steady state, i.e. for κ
d
∈ [0, κ
d-), or for κ
d+ < κ
d
, then it is globally stable. If there are two locally stable nodes, i.e. and for κ
d
∈ (κ
d-, κ
d+), then all flows S
t
(X
0) are attracted to one of them. (See [8] for a delineation of the basin of attraction of and .)
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.
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
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
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
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
Let λ (τ) = α (τ) + iω (τ) be the root of Equation (2.28) satisfying α (τ 0) =0 and ω (τ 0) = ω 0, and set , , , 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.
If r ≥ 0 and Δ = p
2 - 3q < 0, then all roots of Equation (2.28) have negative real parts for all τ ≥ 0. 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). If the conditions of (2) are satisfied, τ = τ
0 and , then ±iω
0 is a pair of simple purely imaginary roots of Equation (2.28) and all other roots have negative real parts. Moreover,
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
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].
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
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
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, 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 accounts for dilution of mRNA concentration due to cell growth in the time interval . Finally, the notation stands for the delayed function . The parameters k 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 is identically equal to , 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
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
The dynamical equation for the concentration of allolactose A is much simpler:
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
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
The function 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) 2, 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 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 is equal to , 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
The functions p
pc
(G
e
) and 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),
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 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.
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:
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
:
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:
The reaction rate constants for repressor binding are proportional to the concentration of active repressors R
A
(T). That is,
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
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.
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.
Dynamics with bursting
Generalities
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
The corresponding master equation for the evolution of the density u (t, x) of x (t) is given by
We consider the situation in which the function f in the burst frequency φ = γκ
b
f is given [103] by
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
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
or
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
or
There are two different bifurcation patterns that are possible. In what will be referred as
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
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
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
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 *.
A 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
The equation for the steady state of (4.14) is of the form
In particular, if condition (4.12) holds and h is geometric as in (4.13) then as in (4.15) is given by
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 . Consequently, we replace Equation (2.31) with a stochastic differential equation in the form
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
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
±
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.
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.
Footnotes
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.