Molecular Distributions in Gene Regulatory Dynamics
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.
keywords:Stochastic modelling, inducible/repressible operon.
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
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
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.
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
are dimensionless constants.
For notational simplicity, henceforth we denote dimensionless concentrations by , and subscripts . Thus we have
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 .
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 .
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
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
Some general observations on the influence of , , and on the appearance of bistability in the deterministic case are in order.
The degree of cooperativity in the binding of effector to the repressor plays a significant role. Indeed, is a necessary condition for bistability.
If then a second necessary condition for bistability is that satisfies Equation 19 so the fractional leakage is sufficiently small.
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.
If , and satisfy these necessary conditions then bistability is only possible if (c.f. Figure 3).
The locations of the minimal and maximal values of bounding the bistable region are independent of .
is a decreasing function of increasing for constant
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:
such that the flow is directed inward everywhere on the surface of . Furthermore, all and
If there is a single steady state, i.e. for , or for , then it is globally stable.
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
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
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 .
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
Case 2. Effector (enzyme) dominated dynamics. Alternately, if then the intermediate is a fast variable relative to the effector and we have
for the relatively slow effector dynamics.
4 Distributions with intrinsic bursting
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
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
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
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
for all initial densities .
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
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
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
and are positive solutions of equation
We explicitly have
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:
Unimodal type 1: and is decreasing for and
Unimodal type 2: and has a single maximum at
at for and
Bimodal type 1: and has a single maximum at for
Bimodal type 2: and has two maxima at , for and
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.