Molecular Distributions in Gene Regulatory Dynamics

Molecular Distributions in Gene Regulatory Dynamics

Michael C. Mackey Departments of Physiology, Physics & Mathematics and Centre for Nonlinear Dynamics, McGill University, 3655 Promenade Sir William Osler, Montreal, QC, CANADA, H3G 1Y6 Marta Tyran-Kamińska Institute of Mathematics, University of Silesia, Bankowa 14, 40-007 Katowice, POLAND Romain Yvinec Université de Lyon CNRS Université Lyon 1 , Bât Braconnier 43 bd du 11 nov. 1918 F-69622 Villeurbanne Cedex France

We show how one may analytically compute the stationary density of the distribution of molecular constituents in populations of cells in the presence of noise arising from either bursting transcription or translation, or noise in degradation rates arising from low numbers of molecules. We have compared our results with an analysis of the same model systems (either inducible or repressible operons) in the absence of any stochastic effects, and shown the correspondence between behaviour in the deterministic system and the stochastic analogs. We have identified key dimensionless parameters that control the appearance of one or two steady states in the deterministic case, or unimodal and bimodal densities in the stochastic systems, and detailed the analytic requirements for the occurrence of different behaviours. This approach provides, in some situations, an alternative to computationally intensive stochastic simulations. Our results indicate that, within the context of the simple models we have examined, bursting and degradation noise cannot be distinguished analytically when present alone.

Stochastic modelling, inducible/repressible operon.

1 Introduction

In neurobiology, when it became clear that some of the fluctuations seen in whole nerve recording, and later in single cell recordings, were not simply measurement noise but actual fluctuations in the system being studied, researchers very quickly started wondering to what extent these fluctuations actually played a role in the operation of the nervous system.

Much the same pattern of development has occurred in cellular and 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 are increasingly focussed on trying to access the role of various types of fluctuations on the operation and fidelity of both simple and complex gene regulatory systems. Recent reviews (Kaern et al., 2005; Raj and van Oudenaarden, 2008; Shahrezaei and Swain, 2008b) give an interesting perspective on some of the issues confronting both experimentalists and modelers.

Typically, the discussion seems to focus on whether fluctuations can be considered as extrinsic to the system under consideration (Shahrezaei et al., 2008; Ochab-Marcinek, 2008, 2010), 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 Elowitz et al. (2002) have used an elegant experimental technique to distinguish between the two, see also Raser and O’Shea (2004), while Swain et al. (2002) and Scott et al. (2006) have laid the groundwork 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” (Blake et al., 2006) or “detrimental” (Fraser et al., 2004) 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.

Here, we consider in detail 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. Our work is motivated by the well documented production of mRNA and/or protein in stochastic bursts in both prokaryotes and eukaryotes (Blake et al., 2003; Cai et al., 2006; Chubb et al., 2006; Golding et al., 2005; Raj et al., 2006; Sigal et al., 2006; Yu et al., 2006), and follows other contributions by, for example, Kepler and Elston (2001), Friedman et al. (2006), Bobrowski et al. (2007) and Shahrezaei and Swain (2008a).

In Section 2 we develop the concept of the operon and treat simple models of the classic inducible and repressible operon. Section 4 considers the effects of bursting alone in an ensemble of single cells. Section 5 then examines the situation in which there are continuous white noise fluctuations in the dominant species degradation rate in the absence of bursting.

2 Generic operons

2.1 The operon concept

The so-called ‘central dogma’ of molecular biology is simple to state in principle, but complicated in its detail. Namely through the process of transcription of DNA, messenger RNA (mRNA, ) is produced and, in turn, through the process of translation of the mRNA proteins (intermediates, ) are produced. There is often feedback in the sense that molecules (enzymes, ) whose production is controlled by these proteins can modulate the translation and/or transcription processes. In what follows we will refer to these molecules as effectors. We now consider both the transcription and translation process in more detail.

In the transcription process an amino acid sequence in the DNA is copied by the enzyme RNA polymerase (RNAP) to produce a complementary copy of the DNA segment encoded in the resulting RNA. Thus this is the first step in the transfer of the information encoded in the DNA. The process by which this occurs is as follows.

When the DNA is in a double stranded configuration, the RNAP is able to recognize and bind to the promoter region of the DNA. (The RNAP/double stranded DNA complex is known as the closed complex.) Through the action of the RNAP, the DNA is unwound in the vicinity of the RNAP/DNA promoter site, and becomes single stranded. (The RNAP/single stranded DNA is called the open complex.) Once in the single stranded configuration, the transcription of the DNA into mRNA commences.

In prokaryotes, translation of the newly formed mRNA commences with the binding of a ribosome to the mRNA. The function of the ribosome is to ‘read’ the mRNA in triplets of nucleotide sequences (codons). Then through a complex sequence of events, initiation and elongation factors bring transfer RNA (tRNA) into contact with the ribosome-mRNA complex to match the codon in the mRNA to the anti-codon in the tRNA. The elongating peptide chain consists of these linked amino acids, and it starts folding into its final conformation. This folding continues until the process is complete and the polypeptide chain that results is the mature protein.

The lactose (lac) operon in bacteria is the paradigmatic example of this concept and this much studied system consists of three structural genes named lacZ, lacY, and lacA. These three genes contain the code for the ultimate production, through the translation of mRNA, of the intermediates -galactosidase, lac permease, and thiogalactoside transacetylase respectively. The enzyme -galactosidase is active in the conversion of lactose into allolactose and then the conversion of allolactose into glucose. The lac permease is a membrane protein responsible for the transport of extracellular lactose to the interior of the cell. (Only the transacetylase plays no apparent role in the regulation of this system.) The regulatory gene lacI, which is part of a different operon, codes for the lac repressor, which is transformed to an inactive form when bound with allolactose, so in this system allolactose functions as the effector molecule.

2.2 The transcription rate function

In this section we examine the molecular dynamics of both the classical inducible and repressible operon to derive expressions for the dependence of the transcription rate on effector levels. (When the transcription rate is constant and independent of the effector levels we will refer to this as the no control situation.)

2.2.1 Inducible regulation

For a typical inducible regulatory situation (such as the lac operon), in the presence of the effector molecule the repressor is inactive (is unable to bind to the operator region preceding the structural genes), and thus DNA transcription can proceed. Let denote the repressor, the effector molecule, and the operator. The effector is known to bind with the active form of the repressor. We assume that this reaction is of the form


where is the effective number of molecules of effector required to inactivate the repressor . Furthermore, the operator and repressor are assumed to interact according to

Let the total operator be :

and the total level of repressor be :

The fraction of operators not bound by repressor (and therefore free to synthesize mRNA) is given by

If the amount of repressor bound to the operator is small


and consequently


where . There will be maximal repression when but even then there will still be a basal level of mRNA production proportional to (which we call the fractional leakage).

If the maximal DNA transcription rate is (in units of inverse time) then, under the assumption that the rate of transcription in the entire population is proportional to the fraction of unbound operators, the variation of the DNA transcription rate with the effector level is given by , or


2.2.2 Repressible regulation

In the classic example of a repressible system (such as the trp operon) in the presence of the effector molecule the repressor is active (able to bind to the operator region), and thus block DNA transcription. We use the same notation as before, but now note that the effector binds with the inactive form of the repressor so it becomes active. We assume that this reaction is of the same form as in Equation 1. The difference now is that the operator and repressor are assumed to interact according to

The total operator is now given by

so the fraction of operators not bound by repressor is given by

Again assuming that the amount of repressor bound to the operator is small we have

where . Now there will be maximal repression when is large, but even at maximal repression there will still be a basal level of mRNA production proportional to . The variation of the DNA transcription rate with effector level is given by or


Both (3) and (4) are special cases of the function


where are given in Table 1.

parameter inducible repressible
Table 1: Definitions of the parameters , , , and . See the text and Section 2.2 for more detail.

2.3 Deterministic operon dynamics in a population of cells

The reader may wish to consult Polynikis et al. (2009) for an interesting survey of techniques applicable to this approach.

We first consider a large population of cells, each of which contains one copy of a particular operon, and let denote mRNA, intermediate protein, and effector levels respectively in the population. Then for a generic operon with a maximal level of transcription (in concentration units), we have dynamics described by the system (Griffith, 1968a, b; Othmer, 1976; Selgrade, 1979)


Here we assume that the rate of mRNA production is proportional to the fraction of time the operator region is active, and that the rates of intermediate and enzyme production are simply proportional to the amount of mRNA and intermediate respectively. All three of the components are subject to random loss. The function is calculated in the previous section.

It will greatly simplify matters to rewrite Equations 6-8 by defining dimensionless concentrations. To this end we first rewrite Equation 5 in the form


where (dimensionless) is defined by


and are defined in Table 1, and we have defined a dimensionless effector concentration through

Further defining dimensionless intermediate () and mRNA concentrations () through

Equations 6-8 can be written in the equivalent form



are dimensionless constants.

For notational simplicity, henceforth we denote dimensionless concentrations by , and subscripts . Thus we have


In each equation, for denotes a net loss rate (units of inverse time), and thus Equations 12-14 are not in dimensionless form.

The dynamics of this classic operon model can be fully analyzed. Let and denote by the flow generated by the system (12)-(14). For both inducible and repressible operons, for all initial conditions the flow for .

Steady states of the system (12)-(14) are in a one to one correspondence with solutions of the equation


and for each solution of Equation 15 there is a steady state of (12)-(14) given by

Whether there is a single steady state or there are multiple steady states will depend on whether we are considering a repressible or inducible operon.

2.3.1 No control

In this case, , and there is a single steady state that is globally asymptotically stable.

2.3.2 Inducible regulation

Single versus multiple steady states. For an inducible operon with given by Equation 2, there may be one ( or ), two ( or ), or three () steady states, with the ordering , corresponding to the possible solutions of Equation 15 (cf. Figure 1). The smaller steady state is typically referred to as an uninduced state, while the largest steady state is called the induced state. The steady state values of are easily obtained from (15) for given parameter values, and the dependence on for and a variety of values of is shown in Figure 1. Figure 2 shows a graph of the steady states versus for various values of the leakage parameter .

Figure 1: Schematic illustration of the possibility of one, two or three solutions of Equation 15 for varying values of with inducible regulation. The monotone increasing graph is the function of Equation 10, and the straight lines correspond to for (in a clockwise direction) , , , , and . This figure was constructed with and for which and as computed from (18). See the text for further details.
Figure 2: Full logarithmic plot of the steady state values of versus for an inducible system, obtained from Equation 15, for and (left to right) illustrating the dependence of the occurrence of bistability on . See the text for details.

Analytic conditions for the existence of one or more steady states can be obtained by using Equation 15 in conjunction with the observation that the delineation points are marked by the values of at which is tangent to (see Figure 1). Simple differentiation of (15) yields the second condition


From equations (15) and (16) we obtain the values of at which tangency will occur:


The two corresponding values of at which a tangency occurs are given by


(Note the deliberate use of as opposed to .)

A necessary condition for the existence of two or more steady states is obtained by requiring that the square root in (17) be non-negative, or


From this a second necessary condition follows, namely


Further, from Equations 15 and 16 we can delineate the boundaries in space in which there are one or three locally stable steady states as shown in Figure 3. There, we have given a parametric plot ( is the parameter) of versus , using

for obtained from Equations 15 and 16. As is clear from the figure, when leakage is appreciable (small , e.g for , ) then the possibility of bistable behaviour is lost.

Figure 3: In this figure we present a parametric plot (for ) of the bifurcation diagram in parameter space delineating one from three steady states in a deterministic inducible operon as obtained from Equations 15 and 16. The upper (lower) branch corresponds to (), and for all values of in the interior of the cone there are two locally stable steady states , while outside there is only one. The tip of the cone occurs at as given by Equations 19 and 20. For there is but a single steady state.
Remark 0.

Some general observations on the influence of , , and on the appearance of bistability in the deterministic case are in order.

  1. The degree of cooperativity in the binding of effector to the repressor plays a significant role. Indeed, is a necessary condition for bistability.

  2. If then a second necessary condition for bistability is that satisfies Equation 19 so the fractional leakage is sufficiently small.

  3. Furthermore, must satisfy Equation 20 which is quite instructive. Namely for the limiting lower limit is while for the minimal value of 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 the larger the ratio must be.

  4. If , and satisfy these necessary conditions then bistability is only possible if (c.f. Figure 3).

  5. The locations of the minimal and maximal values of bounding the bistable region are independent of .

  6. Finally

    1. is a decreasing function of increasing for constant

    2. is an increasing function of increasing for constant .

Local and global stability. The local stability of a steady state is determined by the solutions of the eigenvalue equation (Yildirim et al., 2004)



so (21) can be written as


By Descartes’s rule of signs, (22) will have either no positive roots for or one positive root otherwise. With this information and using the notation SN to denote a locally stable node, HS a half or neutrally stable steady state, and US an unstable steady state (saddle point), then there will be:

  • A single steady state (SN), for

  • Two coexisting steady states (SN) and (HS, born through a saddle node bifurcation) for

  • Three coexisting steady states (SN) for

  • Two coexisting steady states (HS at a saddle node bifurcation), and (SN) for

  • One steady state (SN) for .

For the inducible operon, other work extends these local stability considerations and we have the following result characterizing the global behaviour:

Theorem 2.

(Othmer, 1976; Smith, 1995, Proposition 2.1, Chapter 4) For an inducible operon with given by Equation 3, define . There is an attracting box defined by

such that the flow is directed inward everywhere on the surface of . Furthermore, all and

  1. If there is a single steady state, i.e. for , or for , then it is globally stable.

  2. If there are two locally stable nodes, i.e. and for , then all flows are attracted to one of them. (See Selgrade (1979) for a delineation of the basin of attraction of and .)

2.3.3 Repressible regulation

Figure 4: Schematic illustration that there is only a single solution of Equation 15 for all values of with repressible regulation. The monotone decreasing graph is for a repressible operon, while the straight lines are . This figure was constructed with and . See the text for further details.

As illustrated in Figure 4, the repressible operon has a single steady state corresponding to the unique solution of Equation 15. To determine its local stability we apply the Routh-Hurwitz criterion to the eigenvalue equation (22). The steady state corresponding to will be locally stable (i.e. have eigenvalues with negative real parts) if and only if (always the case) and


The well known relation between the arithmetic and geometric means

when applied to both and gives, in conjunction with Equation 23,

Thus as long as , the steady state corresponding to will be locally stable. Once condition (23) is violated, stability of is lost via a supercritical Hopf bifurcation and a limit cycle is born. One may even compute the Hopf period of this limit cycle by assuming that () in Equation 22 where is the Hopf angular frequency. Equating real and imaginary parts of the resultant yields or

These local stability results tell us nothing about the global behaviour when stability is lost, but it is possible to characterize the global behaviour of a repressible operon with the following

Theorem 3.

(Smith, 1995, Theorem 4.1 & Theorem 4.2, Chapter 3) For a repressible operon with given by Equation 4, define . There is a globally attracting box defined by

such that the flow is directed inward everywhere on the surface of . Furthermore there is a single steady state . If is locally stable it is globally stable, but if is unstable then a generalization of the Poincare-Bendixson theorem (Smith, 1995, Chapter 3) implies the existence of a globally stable limit cycle in .

Remark 0.

There is no necessary connection between the Hopf period computed from the local stability analysis and the period of the globally stable limit cycle.

3 Fast and slow variables

In dynamical systems, considerable simplification and insight into the behaviour can be obtained by identifying fast and slow variables. This technique is especially useful when one is initially interested in the approach to a steady state. In this context a fast variable is one that relaxes much more rapidly to an equilibrium than a slow variable (Haken, 1983). In many systems, including chemical and biochemical ones, this is often a consequence of differences in degradation rates, with the fastest variable the one that has the largest degradation rate. We employ the same strategy here to obtain approximations to the population level dynamics that will be used in the next section.

It is often the case that 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

Consequently the three variable system describing the generic operon reduces to a two variable one involving the slower intermediate and effector:


In our considerations of specific single operon dynamics below we will also have occasion to examine two further subcases, namely

Case 1. Intermediate (protein) dominated dynamics. If it should happen that (as for the lac operon, then the effector also qualifies as a fast variable so

and thus from (24)-(25) we recover the one dimensional equation for the slowest variable, the intermediate:


Case 2. Effector (enzyme) dominated dynamics. Alternately, if then the intermediate is a fast variable relative to the effector and we have

so our two variable system ( 24)-(25) reduces to a one dimensional system


for the relatively slow effector dynamics.

Both Equations 26 and 27 are of the form


where is either for protein () dominated dynamics or for effector () dominated dynamics.

4 Distributions with intrinsic bursting

4.1 Generalities

It is well documented experimentally (Cai et al., 2006; Chubb et al., 2006; Golding et al., 2005; Raj et al., 2006; Sigal et al., 2006; Yu et al., 2006) that in some organisms the amplitude of protein production through bursting translation of mRNA is exponentially distributed at the single cell level with density


where is the average burst size, and that the frequency of bursting is dependent on the level of the effector. Writing Equation 29 in terms of our dimensionless variables we have

Remark 0.

The technique of eliminating fast variables described in Section 2.3 above (also known as the adiabatic elimination technique (Haken, 1983)) has been extended to stochastically perturbed systems when the perturbation is a Gaussian distributed white noise, c.f. Stratonovich (1963, Chapter 4, Section 11.1), Wilemski (1976), Titular (1978), and Gardiner (1983, Section 6.4). However, to the best of our knowledge, this type of approximation has never been extended to the situation dealt with here in which the perturbation is a jump Markov process.

The single cell analog of the population level intermediate protein dominated Case 1 above (when ) is


where denotes a jump Markov process, occurring at a rate , whose amplitude is distributed with density as given in (30). Analogously, in the Case 2 effector dominated situation the single cell equation becomes


Equations 31 and 32 can both be written as

Remark 0.

In the case of bursting we will always take in contrast to the deterministic case where .

From Mackey and Tyran-Kamińska (2008) the corresponding operator equation for the evolution of the density when there is a single dominant slow variable is given by

Remark 0.

This is a straightforward generalization of what Gardiner (1983, Section 3.4) refers to as the differential Chapman-Kolmogorov equation.

Stationary solutions of (33) are solutions of


If there is a unique stationary density, then the solution of Equation 33 is said to be asymptotically stable (Lasota and Mackey, 1994) in the sense that

for all initial densities .

Theorem 8.

(Mackey and Tyran-Kamińska, 2008, Theorem 7). The unique stationary density of Equation 34, with given by Equation 9 and given by (29), is


where is a normalizing constant such that . Further, is asymptotically stable.

Remark 0.

The stationary density (35) is found by rewriting Equation 34 in the form

using Laplace transforms and solving by quadratures. Note also that we can represent as

where is a normalizing constant.

4.2 Distributions in the presence of bursting

4.2.1 Protein distribution in the absence of control

If the burst frequency is independent of the level of all of the participating molecular species, then the solution given in Equation 35 is the density of the gamma distribution:

where denotes the gamma function. For , and is decreasing while for , and there is a maximum at .

4.2.2 Controlled bursting

We next consider the situation in which the burst frequency is dependent on the level of , c.f. Equation 5. This requires that we evaluate

where are enumerated in Table 1 for both the inducible and repressible operons treated in Section 2.2 and

Consequently, the steady state density (35) explicitly becomes


The first two terms of Equation 36 are simply proportional to the density of the gamma distribution. For we have while for , and there is at least one maximum at a value of . We have for all and from Remark 9 it follows that


Observe that if then is a monotone decreasing function of , since for all . Thus we assume in what follows that .

Since the analysis of the qualitative nature of the stationary density leads to different conclusions for the inducible and repressible operon cases, we consider each in turn.

4.2.3 Bursting in the inducible operon

For , as in the case of an inducible operon, the third term of Equation 36 is a monotone increasing function of and, consequently, there is the possibility that may have more than one maximum, indicative of the existence of bistable behaviour. In this case, the stationary density becomes

From (37) it follows that we have for if and only if


Again, graphical arguments (see Figure 5) show that there may be up to three roots of (38).

Figure 5: Schematic illustration of the possibility of one, two or three solutions of Equation 38 for varying values of with bursting inducible regulation. The straight lines correspond (in a clockwise direction) to , , (and respectively , , ), , and . This figure was constructed with , and for which and as computed from (42). See the text for further details.

For illustrative values of , , and , Figure 6 shows the graph of the values of at which as a function of . When there are three roots of (38), we label them as .

Generally we cannot determine when there are three roots. However, we can determine when there are only two roots from the argument of Section 2.3.2. At and we will not only have Equation 38 satisfied but the graph of the right hand side of (38) will be tangent to the graph of the left hand side at one of them so the slopes will be equal. Differentiation of (38) yields the second condition


We first show that there is an open set of parameters for which the stationary density is bimodal. From Equations 38 and 39 it follows that the value of at which tangency will occur is given by

and are positive solutions of equation

We explicitly have

provided that


Equation 40 is always satisfied when or when and is as in the deterministic case (19). Observe also that we have for and for . The two corresponding values of at which a tangency occurs are given by

If then and is decreasing for , while for there is a local maximum at . If then and has one or two local maximum. As a consequence, for we have a bimodal steady state density if and only if the parameters and satisfy (40), , and .

We now want to find the analogy between the bistable behavior in the deterministic system and the existence of bimodal stationary density . To this end we fix the parameters and and vary as in Figure 5. Equations 38 and 39 can also be combined to give an implicit equation for the value of at which tangency will occur


and the corresponding values of are given by


There are two cases to distinguish.

Case 1. . In this case, . Further, the same graphical considerations as in the deterministic case show that there can be none, one, or two positive solutions to Equation 38. If , there are no positive solutions, is a monotone decreasing function of . If , there are two positive solutions ( and in our previous notation, has become negative and not of importance) and there will be a maximum in at with a minimum in at .

Case 2. . Now, and there may be one, two, or three positive roots of Equation 38. We are interested in knowing when there are three which we label as as will correspond to the location of maxima in while will be the location of the minimum between them and the condition for the existence of three roots is .

We see then that the different possibilities depend on the respective values of , , , and . To summarize, we may characterize the stationary density for an inducible operon in the following way:

  1. Unimodal type 1: and is decreasing for and

  2. Unimodal type 2: and has a single maximum at

    1. for or

    2. at for and

  3. Bimodal type 1: and has a single maximum at for

  4. Bimodal type 2: and has two maxima at , for and

Remark 0.

Remember that the case cannot display bistability in the deterministic case. However, in the case of bursting in the inducible system when , if and , then and also has a maximum at . Thus in this case one can have a Bimodal type 1 stationary density.

Figure 6: Full logarithmic plot of the values of at which versus the parameter , obtained from Equation