The effects of intrinsic noise on the behaviour of bistable cell regulatory systems under quasi-steady state conditions

The effects of intrinsic noise on the behaviour of bistable cell regulatory systems under quasi-steady state conditions

Roberto de la Cruz Centre de Recerca Matemàtica. Edifici C, Campus de Bellaterra, 08193 Bellaterra (Barcelona), Spain. Departament de Matemàtiques, Universitat Atonòma de Barcelona, 08193 Bellaterra (Barcelona), Spain.    Pilar Guerrero Department of Mathematics, University College London, Gower Street, London WC1E 6BT, UK.    Fabian Spill Department of Biomedical Engineering, Boston University, 44 Cummington Street, Boston MA 02215, USA. Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA.    Tomás Alarcón Centre de Recerca Matemàtica. Edifici C, Campus de Bellaterra, 08193 Bellaterra (Barcelona), Spain. Departament de Matemàtiques, Universitat Atonòma de Barcelona, 08193 Bellaterra (Barcelona), Spain.
July 3, 2019

We analyse the effect of intrinsic fluctuations on the properties of bistable stochastic systems with time scale separation operating under1 quasi-steady state conditions. We first formulate a stochastic generalisation of the quasi-steady state approximation based on the semi-classical approximation of the partial differential equation for the generating function associated with the Chemical Master Equation. Such approximation proceeds by optimising an action functional whose associated set of Euler-Lagrange (Hamilton) equations provide the most likely fluctuation path. We show that, under appropriate conditions granting time scale separation, the Hamiltonian can be re-scaled so that the set of Hamilton equations splits up into slow and fast variables, whereby the quasi-steady state approximation can be applied. We analyse two particular examples of systems whose mean-field limit has been shown to exhibit bi-stability: an enzyme-catalysed system of two mutually-inhibitory proteins and a gene regulatory circuit with self-activation. Our theory establishes that the number of molecules of the conserved species are order parameters whose variation regulates bistable behaviour in the associated systems beyond the predictions of the mean-field theory. This prediction is fully confirmed by direct numerical simulations using the stochastic simulation algorithm. This result allows us to propose strategies whereby, by varying the number of molecules of the three conserved chemical species, cell properties associated to bistable behaviour (phenotype, cell-cycle status, etc.) can be controlled.

I Introduction

The networks of interacting genes and proteins that are responsible for regulation, signalling and response, and which, ultimately, orchestrate cell function, are under the effect of noise Kepler and Elston (2001); Kaern et al. (2005); Maheshri and O’Shea (2007); Losick and Desplan (2008); Raj and van Oudenaarden (2008). This randomness materialises in the form of fluctuations of the number of molecules of the species involved, subsequently leading to fluctuations in their activity. Besides external perturbations, biochemical reactions can be intrinsically noisy, especially when the number of molecules is very low.

Far from necessarily being a mere disturbance, fluctuations are an essential component of the dynamics of cellular regulatory systems which, in many instances, are exploited to improve cell function Cai, Dalal, and Elowitz (2008); Eldar and Elowitz (2010). For example, randomness has been shown to enhance the ability of cells to adapt and increase their fitness in random or variable environments Kussell and Leibler (2005); Acar, Mettetal, and van Oudenaarden (2008); Guerrero et al. (2015). Random noise also serves the purpose of assisting cell populations to sustain phenotypic variation by enabling cells to explore the phase space Maheshri and O’Shea (2007); Raj and van Oudenaarden (2008); Losick and Desplan (2008); MacArthur, Please, and Oreffo (2008); Eldar and Elowitz (2010); Balazsi, van Oudenaarden, and Collins (2011).

One of the mechanisms that allows noise-induced phenotypic variability relays on multi-stability Cinquin and Demongeot (2005); Jaeger and Monk (2014). The basis of this mechanism was first proposed by Kauffman Kauffman (1993), who associated phenotypes or differentiated states to the stable attractors of the dynamical systems associated to gene and protein interaction networks. In the presence of noise, the corresponding phase space generates an epigenetic landscape, where cells exposed to the same environment and signalling cues coexist in different cellular phenotypes Huang (2012).

Multi-stability is also an essential element in the control of cell response and function via signalling pathways Tyson, Chen, and Novak (2003). In particular, bi-stability as a means to generate reliable switching behaviour is widely utilised in numerous pathways such as the apoptosis Legewie, Bluthgen, and Herzel (2006), cell survival Legewie, Bluthgen, and Herzel (2007), differentiation Kalmar et al. (2009), and cell-cycle progression Ferrel and Xiong (2001); Tyson and Novak (2001) pathways. For example, bi-stability is used to regulate such critical cell functions such as the transition from quiescence to proliferation through bistable behaviour associated with the Rb-E2F switch within the regulatory machinery of the mammalian cell-cycle Gerard and Goldbeter (2009, 2012); Yao et al. (2012); Yao (2014); Gerard and Goldbeter (2014); Bedessem and Stephanou (2014).

A common theme which appears when trying to model cell regulatory systems is separation of time scales, i.e. the presence of multiple processes evolving on widely diverse time scales. When noise is ignored and systems are treated in terms of deterministic mean-field descriptions, such separation of time scales and the associated slow-fast dynamics are often exploited for several forms of model reduction, of which one of the most common is the so-called quasi-steady state approximation (QSSA) Keener and Sneyd (1998). This approximation is ubiquitously used whenever regulatory processes involve enzyme catalysis, which is a central regulation mechanism in cell function Tyson, Chen, and Novak (2003). In this paper, we investigate the effects of intrinsic noise on the bi-stability of two particular systems, namely, an enzyme-catalysed system of mutual inhibition and a gene regulatory circuit with self-activation. The mean-field limit of both these systems has been shown to exhibit bi-stability Tyson and Novak (2001); Frigola et al. (2012). The aim of this paper is to analyse how noise alters the mean-field behaviour associated to these systems when they operate under quasi-steady state conditions.

We note that this work does not concern the subject of noise-induced bifurcations García-Ojalvo and Sancho (1999). Such phenomenon has been studied in many situations, including biological systems. An example which is closely related to the systems we analyse here is the so-called enzymatic futile cycles. Samoilov et al. Samoilov, Plyasunov, and Arkin (2005) have shown that noise associated to the number of enzymes induce bistability. In the absence of this source of noise, i.e. in the mean-field limit, the system does not exhibit bistable behaviour. The treatment of this phenomena would require to go to higher orders in the WKB expansion, which we do not explore here.

The issue of separation of time scales in stochastic models of enzyme catalysis has been addressed using a number of different approaches. Several such analysis have been carried out in which the QSSA is directly applied to the master equation by setting the fast reactions in partial equilibrium (i.e. the probability distribution corresponding to the fast variables remains unchanged), and letting the rest of the system to evolve according to a reduced stochastic dynamic Rao and Arkin (2003); Turner, Schnell, and Burrage (2004). Other approaches have been proposed such as the QSSA to the exact Fokker-Planck equation that can be derived from the Poisson representation of the chemical master equation Thomas, Straube, and Grima (2010). Approaches based on enumeration techniques have also been formulated Dóka and Lente (2012). Furthermore, Thomas et al.Thomas, Grima, and Straube (2012) have recently formulated a rigorous method to eliminate fast stochastic variables in monostable systems using projector operators within the linear noise approximation Thomas, Grima, and Straube (2012). Methods for model reduction based on perturbation analysis have been developed in Alarcón (2014); Bruna, Chapman, and Smith (2014). Additionally, driven by the need of more efficient numerical methods, there has been much activity regarding the development of numerical methods for stochastic systems with multiple time-scales Burrage, Tian, and Burrage (2004); Vlachos (2005); MacNamara, Burrage, and Sidje (2008). Several of these methods are variations of the stochastic simulation algorithm Rao and Arkin (2003); Cao, Gillespie, and Petzold (2005a, b); Samant and Vlachos (2005); E, Liu, and Vanden-Eijnden (2007); Sanft, Gillespie, and Petzold (2011) or the -leap method Rathinam et al. (2006) where the existence of fast and slow variables is exploited to enhance their performance with respect to the standard algorithms. Another family of such numerical methods is that of the so-called hybrid methods, where classical deterministic rate equations or stochastic Langevin equations for the fast variables are combined with the classical stochastic simulation algorithm for the slow variables Haseltine and Rawlings (2002); Salis and Kaznessis (2005). Other related methods were studied in Assaf and Meerson (2006); Newby and Chapman (2013); Kang and Kurtz (2013).

Here, we advance the formalism developed in Alarcón (2014), in which a method based on the semi-classical approximation of the Chemical Master Equation allows to evaluate the effects of intrinsic random noise under quasi-steady conditions. In our analysis of the Michaelis-Menten model of enzyme catalysis in Alarcón (2014), we showed that the semi-classical quasi-steady state approximation reveals that the velocity of the enzymatic reaction is modified with respect to the mean-field estimate by a quantity which is proportional to the total number of molecules of the (conserved) enzyme. In this paper, we extend this formalism to show that, associated to each conserved molecular species, the associated (constant) number of molecules is a bifurcation parameter which can drive the system into bi-stability beyond the predictions of the mean-field theory. We then proceed to test our theoretical results by means of direct numerical simulation of the Chemical Master Equation using the stochastic simulation algorithm Gillespie (1976). We should note the Hamiltonian formalism derived from the semi-classical approximation is formulated on a continuum of particles, which requires the number of particles to be large enough. This must hold true for all the species in our model, both fast and slow. Since this separation between fast and slow species is based on their relative abundance, one must be careful that the scaling assumptions are consistent, particularly in the case of the model of self-activating gene regulatory circuit where the number of binding sites is typically small. This assumption, however, has been used in previous studies Assaf, Roberts, and Luthey-Schulten (2011). Also we show that our simulation results of the full stochastic processes agree with our analysis and, therefore, our re-scaled equations are able to predict the behaviour of the system. We note that the mean-field limit, which is conventionally obtained by ignoring noise in the limit of large particle numbers, is obtained by setting the momenta in our phase-space formalism to .

The approximation we develop in this paper falls within the general framework of the optimal fluctuation path theory Bressloff (2014). This framework is a particular case of the large deviation theory which allows us to study rare events (i.e. events whose frequency is exponentially small with system size). Within these framework we will show that, upon carrying out the QSSA, the only source of noise in the system is associated to the random initial conditions of the species whose numbers are conserved. We therefore predict that a population of cells, each having a random number of conserved molecules, will have a bimodal distribution.

This paper is organised as follows. Section 2 is devoted to a detailed exposition of the semi-classical quasi-steady state approximation for stochastic systems. In Sections 3 and 4, we apply this formalism to analyse the behaviour of a bistable enzyme-catalysed system and a gene regulatory circuit of auto-activation, respectively. We will show that our semi-classical quasi-steady state theory allows us to study the effect of intrinsic noise on the behaviour of these systems beyond the predictions of their mean-field descriptions. We also verify our theoretical predictions by means of direct stochastic simulations. Finally in Section 5, we summarise our results and discuss their relevance.

Ii Semi-classical quasi-steady state approximation

Our aim in this paper is to formulate a stochastic generalisation of the quasi-steady state approximation for enzyme-catalysed reactions and simple circuits of gene regulation and use such approximation to determine if the presence of noise has effects on the behaviour of the system beyond the predictions of the corresponding mean-field models. Specifically, we analyse stochastic systems for which the mean-field models predicts bi-stability and investigate how such behaviour is affected by stochastic effects. Our analysis is carried out in the context of Markovian models of the corresponding reaction mechanisms formulated in terms of the so-called chemical master equation (CME) Kampen (2007). Two example of such stochastic systems, a bistable enzyme-catalysed system and a gene regulatory circuit of auto-activation, are formulated and analysed in detail in Sections III and IV, respectively. Following Alarcón (2014), we formulate the QSS approximation for the asymptotic solution of the CME obtained by means of large deviations/WKB approximations Kubo, Matsuo, and Kitahara (1973); Alarcón and Page (2007); Touchette (2009). The CME is given:


where is the transition rate corresponding to reaction channel and is a vector whose entries denote the change in the number of molecules of each molecular species when reaction channel fires up, i.e. .

An alternative way to analyse the dynamics of continuous-time Markov processes on a discrete space of states is to derive an equation for the generating function, of the corresponding probabilistic density:


where is the solution of the Master Equation (1). satisfies a partial differential equation (PDE) which can be derived from the Master Equation. This PDE is the basic element of the so-called momentum representation of the Master Equation Doi (1976); Peliti (1985); Assaf and Meerson (2006); Assaf, Meerson, and Sasorov (2010); Kang and Kurtz (2013).

Although closed, analytic solutions are rarely available, the PDE for the generating function admits a perturbative solution, which is commonly obtained by means of the WKB method Assaf, Meerson, and Sasorov (2010). More specifically, the (linear) PDE that governs the evolution of the generating function can be written as:


where the operator is determined by the reaction rates of the Master Equation (1). Furthermore, the solution to this equation must satisfy the normalisation condition for all . This PDE, or, equivalently, the operator , are obtained by multiplying both sides of the Master Equation (1) by and summing up over all the possible values of

From the mathematical point of view, Eq. (3) is a Schrödinger-like equation and, therefore, there is a plethora of methods at our disposal in order to analyse it. In particular, when the fluctuations are (assumed to be) small, it is common to resort to WKB methods Kubo, Matsuo, and Kitahara (1973); Alarcón and Page (2007); Gonze, Halloy, and Gaspard (2002). This approach is based on the WKB-like Ansatz that . By substituting this Ansatz in Eq. (3) we obtain the following Hamilton-Jacobi equation for the function :


Instead of directly tackling the explicit solution of Eq. (4), we will use the so-called semi-classical approximation. We use the Feynman path-integral representation which yields a solution to Eq. (3) of the type Peliti (1985); Feynman and Hibbs (2010); Kubo, Matsuo, and Kitahara (1973); Dickman and Vidigal (2003); Elgart and Kamenev (2004); Täuber, Howard, and Vollmayr-Lee (2005):


where indicates integration over the space of all possible trajectories and is given by Kubo, Matsuo, and Kitahara (1973):


where the position operators in the momentum representation have been defined as with the commutation relation . corresponds to the action associated with the generating function of the probability distribution function of the initial value of each variable, , which are assumed to be independent random variables.

The so-called semi-classical approximation consists of approximating the path integral in Eq. (5) by


where are now the solutions of the Hamilton equations, i.e. the orbits which maximise the action :


where the pair (,) are the generalised coordinates corresponding to chemical species . These equations are (formally) solved with boundary conditionsElgart and Kamenev (2004) , , where is the initial number of molecules of species .

Eqs. (8)-(9) are the starting point for the formulation of the semi-classical quasi-steady state approximation (SCQSSA) Alarcón (2014). In order to proceed further, we assume, as per the Briggs-Haldane treatment of the Michealis-Menten model for enzyme kinetics Briggs and Haldane (1925); Keener and Sneyd (1998), that the species involved in the system under scrutiny are divided into two groups according to their characteristic scales. More specifically, we have a subset of chemical species whose numbers, , scale as:


where , whilst the remaining species are such that their numbers, , scale as:


where . Key to our approach is the fact that and must be such that:


We further assume that the generalised coordinates, , scale in the same fashion as the corresponding variable , i.e.


where . We refer to the variables belonging to this subset as slow variables. Similarly,


where , which are referred to as fast variables. Moreover, we assume that the moment coordinates, , are all independent of and , and therefore remain invariant under rescaling.

Under this scaling for the generalised coordinates, we define the following scale transformation for the Hamiltonian in Eq. (6):


where identifies the reaction with the largest order among all the reactions that compose the dynamics and is the corresponding rate constant. For example, in the case of the bistable enzyme-catalysed system whose reactions or elementary events and the corresponding transition rates are given in Table 1, , as this reaction is order 3 whereas all the others are order 0, 1, or 2. In the case of the self-activating gene regulatory circuit, Table 4, , since this reaction is order 3 whereas the remaining ones are order 1 at most. The exponents and correspond to the number of slow and fast variables involved in the transition rate , respectively.

The last step is to rescale the time variable so that a dimensionless variable, , is defined such that:


It is now a trivial exercise to check that, upon rescaling, Eqs. (8)-(9) read


for the slow variables. By contrast, rescaling of the Hamilton equations corresponding to the subset of fast variables leads to:


where is defined in Eq. (12). The QSS approximation consists on assuming that and in Eqs. (19)-(20),


resulting in a differential-algebraic system of equations which provides us with the semi-classical quasi-steady state approximation (SCQSSA).

Iii Bistable enzyme-catalysed systems

Figure 1: Reactions for the bistable enzyme-catalysed system proposed by Tyson & Novak Tyson and Novak (2001). represents active Cdh/Apc, inactive Cdh/Apc, inactivating enzymes, activating enzymes, active Cdh/Apc-inactivating-enzyme complexes, inactive Cdh/Apc-activating-enzyme complexes, and the number of CycB-CDK complexes. The first two reactions correspond to enzyme-catalysed inactivation and activation of Cdh/APC. The third reaction corresponds to the dynamics of CycB activity: synthesis at a constant rate, , and degradation by natural decay and active Cdh/Apc-induced inactivation.
Variable Description
, Number of active and inactive (respectively) Cdh1 molecules
, Number of Cdh1-inactivating and Cdh1-activating (respectively) enzyme molecules
, Number of enzyme-active Cdh1 and enzyme-inactive Cdh1 (respectively) complexes
, Number of active cyclin molecules
Transition rate r Event
Enzyme and active Cdh1 form complex
Enzyme-active Cdh1 complex splits
Inactivation of Cdh1 and enzyme release
Enzyme and inactive Cdh1 form complex
Enzyme-inactive Cdh1 complex splits
Activation of Cdh1 and enzyme release
CycB synthesis
CycB degradation
Table 1: Random variables and transition rates of the stochastic model associated to the enzymatic reaction shown in Fig. 1.

As a prototype of a bistable enzyme-catalysed system, we analyse a stochastic system proposed in Alarcón (2014); Guerrero and Alarcón (2015), whose mean-field limit has been shown to correspond to a bistable system which is a part of a model for the G/S transition of the eukaryote cell cycle proposed in Tyson and Novak (2001). Tyson & Novak Tyson and Novak (2001) have formulated a (deterministic) model of the cell cycle such that the core of the system regulating the G/S transition is a system of two mutually-repressing proteins (Cdh1 and CycB). This system of mutual repression gives rise to a bistable system where one of the stable steady states is identified with the G phase whereas the other corresponds to a state where the cell is ready to go through the other three phases of the cell-cycle, known as S, G, and M. This central module, which is the one we focus on, is acted upon by a complex regulatory network which monitors if conditions are met for the cell to undergo this transition and accounts for its accurate timing. Presently, we ignore this network and focus on the central bistable system. It is shown in Tyson and Novak (2001) that the mean field version of the model exhibits bistable behaviour as a function of a bifurcation parameter , i.e. the mass of the cell. For very small values of , the system is locked into a high (low) Cdh1(CycB)-level stable fixed point (i.e. into the G phase). For very large values , the system has only one stable steady state corresponding to a low (high) Cdh1(CycB)-level fixed point. For intermediate values of the system exhibits bistability, i.e. both of these stable fixed points coexist with an unstable saddle point. In this section, we focus on how noise alters the behaviour of the mean-field dynamics.

The transition rates corresponding to the different reactions involved in the stochastic model associated to the enzyme-regulated kinetics shown in Fig. 1 are given in Table 1. This kinetics corresponds to the enzyme regulated activation and inhibition of Cdh1 (an inhibitor of cell-cycle progression). Cdh1 inactivation is further (up)regulated by the presence of CycB, an activator of cell-cycle progression. CycB is synthesised and degraded at basal rates and is further degraded in the presence of active Cdh1 (see Fig. 1). Therefore, the resulting dynamics leads to a system with mutual inhibition which produces bistable behaviour. It is important to note that the associated reaction kinetics exhibits three conservation laws (see Table 1): , , and . The first two of these conservation laws are associated to the conservation of the number of Cdh1-inhibiting and Cdh1-activating enzymes, respectively, whilst the latter expresses the conservation of the total number of Cdh1 molecules. The quantities and are the (conserved) number of enzymes and Cdh1, respectively. Note that, as per the methodology developed in Section II, we assume that and .

Rescaled variables Dimensionless parameters
Table 2: Dimensionless variables used in Eqs. (III). and are the average concentration of Cdh1 (active plus inactive) and the average concentration of both Cdh1-activating and Cdh1-inactivating enzymes, respectively. We further assume that the stationary concentration of active CycB also scales with .

The corresponding stochastic Hamiltonian, , which is derived by applying the methodology of Section II to the Master Equation associated to the chemical kinetics described in Table 1, can be split into three parts,


where is the Hamiltonian corresponding to the CycB-regulated enzymatic inactivation of Cdh1 (reactions 1 to 3 in Table 1):


corresponds to enzymatic activation of Cdh1 (reactions 4 to 6 in Table 1):


and, finally, , which corresponds to synthesis and degradation of CycB, is given by (reactions 7 and 8 in Table 1):


We now proceed to apply the procedure explained in Section 2 in order to obtain the SCQSSA for the system determined by the transition rates given in Table 1. We first need to determine which of the variables are slow variables and which ones are fast variables. As shown in Table 2, the pairs , , and , corresponding to the active and inactive forms of Cdh1 and to CycB, respectively, are the slow generalised coordinates, as the generalised positions scale with . The remaining generalised coordinates scale as and, therefore, are fast variables. Furthermore, the rescaled Hamiltonian is given by:






The rescaled parameters are given in Table 2. Last, by rescaling time and defining the dimensionless time variable as (Table 2), the SCQSSA equations (17)-(18) and (21)-(22) lead to (see Alarcón (2014) for a detailed derivation):


where , , and are constants to be determined and and , and and . Note that for , with , to hold must be satisfied. In this case, we have


As shown in Alarcón (2014), the parameter values are determined by comparing the corresponding mean-field approximation, which is obtained by taking Elgart and Kamenev (2004), and , i.e. the total number of molecules of Cdh1 and its activating and inhibiting enzymes be exactly equal to its average, i.e. and , to the system originally proposed by Tyson & Novak Tyson and Novak (2001). In Eq. (36) we have redefined in order to make explicit the dependence on the bifurcation parameter, , as used by Tyson & Novak Tyson and Novak (2001). The parameter values are shown in Table 3.

Rescaled parameter Parameter Units Reference
min Tyson and Novak (2001)
min Tyson and Novak (2001)
min Tyson and Novak (2001)
min Tyson and Novak (2001)
min Tyson and Novak (2001)
Dimensionless Alarcón (2014)
Dimensionless Alarcón (2014)
min Rao and Arkin (2003)
Dimensionless Alarcón (2014)
Dimensionless Tyson and Novak (2001)
Table 3: Parameter values used in simulations of the stochastic bistable enzyme-catalysed system

Figure 2: (a) Bifurcation analysis for the SCQSS approximation of the stochastic bistable enzyme-catalysed system Eqs. (36)-(39). The panels on the top plot (a) shows the bifurcation diagrams for different values of the parameters , and . If and are random Poisson variables with parameter and , respectively, then (see Eq. III). In these panels solid lines correspond to , dot-dashed lines to , and dashed lines to . The bottom plot (b) shows the bi-stability boundaries in parameter space. The region between the boundaries corresponds to the bistable region of the stochastic Tyson & Novak system according to the SCQSS approximation.

Upon rescaling of the variables (Table 2) and the Hamiltonian (Eq. (15)), the action functional reads:


It is straightforward to check that in SCQSSA conditions . Furthermore, since const. and , and for the fast generalised coordinates, the SCQSS approximation of the action Eq. (40), , reduces to:


where, as per the SCQSSA, and are determined by Eqs. (38) and (39), respectively, , which implies , and , and are constants that remain to be determined. In order to do so, we resort to the method developed in reference Alarcón (2014). The quasi-steady state characteristic function, is given by:


where is the generating function of the probability distribution for the initial condition of species . In Alarcón (2014), we have shown that, applying a Laplace-type asymptotic method Murray (1984); Ablowitz and Fokas (2003) to the integrals


where, , and can be given as functions of and , i.e. the initial numbers of Cdh1 molecules and Cdh1-inactivating and Cdh1-activating enzymes, respectively:


, and are the probabilities that initially takes the value and that and have initial values and . These probabilities can be interpreted to correspond to variability in the abundance of these enzymes within a population of cells. A particularly simple case results from assuming that , and are Poisson distributions with parameter and , respectively. In this case Alarcón (2014):


Note that, in the particular case in which the total numbers of Cdh1 and enzyme molecules are random Poisson variables, we have that , , and .

iii.1 Bifurcation analysis

Figure 3: Simulation results for the stochastic bistable enzyme-catalysed system Table 1. We have plotted the probability where and for different values of . The initial number of Cdh1-inactivating and Cdh1-activating enzymes are fixed according to and , respectively. . We aim to check our predictions regarding the effect of the ratio on the stability properties of the system. According to our results shown in Fig. 2, decreasing the ratio between the number of Cdh1-inactivating () and Cdh1-activating () enzymes, the system should be driven away from bistability and into the stable G-phase regime (see Fig. 2(b)). The remaining parameter values are inferred from those given by Tyson & Novak Tyson and Novak (2001) as shown in Tables 2 and 3. We see that when varying , the system switches from a state of high () ro a state of low (), whereas at the intermediate levels of (e.g. and ) the system is in a bistable state. We take in all the simulations shown in this figure. Average is performed over 1000 realisations.

Fig. 2 shows results regarding the bifurcation behaviour of the SCQSS approximation of the stochastic bistable enzyme-catalysed system Eqs. (36)-(39). In particular we are interested in a comparison between the bistable behaviour of the mean-field model, corresponding to taking for all , and that of the SCQSS approximation with , and given by Eq. (III). i.e. they are determined as functions of and .

We have shown that both the ratio of and , , and alter the bistable behaviour of the system beyond the predictions of the mean-field model. In particular, we observe that decreasing the value of extends the region of stability of the G-fixed point, i.e. the fixed point corresponding to the steady-state value of , such that . By contrast, when is increased the stability region of the G-fixed point shrinks. Intuitively, given the relation between and and the number of Cdh1-inactivating and Cdh1-activating enzyme, this result is straightforward to interpret: decreasing the number of Cdh1-inactivating enzyme demands a larger value of in order to de-stabilise the G-fixed point. This is fully confirmed by direct simulation using Gillespie stochastic simulation algorithm Gillespie (1976). Fig. 3 shows simulation results in which we compute the probability for different values of . has been chosen so that the system has reached steady state conditions. We observe, that for and , the system evolves towards the -fixed point (i.e. the S-G-M fixed point). As decreases, i.e. there is more Cdh1-inactivating enzyme than Cdh1-activating enzyme, the system enters the bistable regime. If reaches low-enough values (depending upon the initial condition), we may even observe an exchange of stability, i.e. the system evolves towards the -fixed point.

Figure 4: Simulation results for the stochastic bistable enzyme-catalysed system Table 1. We have plotted the probability where and with different initial conditions and different values of . Average is performed over 1000 realisations. and and . The remaining parameter values are inferred from those given by Tyson & Novak Tyson and Novak (2001) as shown in Tables 2 and 3. We see that when varying , the system switches from a state of high () ro a state of low (), whereas at the intermediate levels of the system is in a bistable state.

Regarding the dependence on , we have checked the predictions of the SCQSS approximation by means of simulations with different values of . Figure 2 shows the bi-stability region of system Eqs. (36)-(39) in -space, where . For a fixed value of , there is a threshold value for below which the system stops being bistable to become entrapped into the the S-G-M fixed point (i.e. ). In order to validate this prediction, we have conducted stochastic simulations for different values of . Figure 4 shows simulation results for . We observe that for small values of , the system is locked into the the S-G-M fixed point, as predicted by the SCQSS approximation. As increases, the system enters a fluctuation-dominated bistable regime where, as the system goes through the bifurcation point, the system undergoes bistable behaviour. This behaviour is typical in a system undergoing a phase transition, where fluctuations unboundedly increase Goldenfeld (1992). Finally, as continues to increase, the system becomes trapped into G-fixed point (see Figure 4). These results fully reproduce the behaviour predicted by our SCQSSA stability analysis.

Figure 5: Plots showing the variance where associated to the simulation results shown in Fig. 3 (panel (a)) and in Fig. 4 (panel (b)). These plots show how changes as the control parameter (, for the simulations associated to plot (a), and for the simulations shown in plot (b)). The maximum of as a function of the control parameter helps us to quantitatively determine the corresponding critical value Goldenfeld (1992).

The aforementioned behaviour regarding unbounded increase of fluctuations close to a bifurcation Goldenfeld (1992) is used to locate the critical value of the associated control parameter, i.e. and for the simulations shown in Figs. 3 and 4, respectively. This property allows us to do a quantitative comparison between the simulations and asymptotic analysis. To this end, we plot how the variance, where , changes as the corresponding control parameter varies. Regarding the results shown in Fig. 5(a) (associated to the simulations shown in Fig. 3), we observe that the critical value of the control parameter , , is approximately , which, taking into account that , implies that the critical value of the renormalized mass, , . Our asymptotic analysis predicts that (see Fig. 2(b) with ). The results shown in in Fig. 5(b) (corresponding to the simulations shown in Fig. 4), the critical value of , , is approximately . The prediction of our asymptotic analysis (see Fig. 2(b) with ) is .

Iv Auto-activation gene regulatory circuit

We now proceed to analyse the effects of intrinsic noise in a model of a bistable self-activation gene regulatory circuit Assaf, Roberts, and Luthey-Schulten (2011); Frigola et al. (2012); Weber and Buceta (2013) in the context of the quasi-steady regime. Many instances of genetic switches, i.e. bistable gene regulatory circuits, have been identified Gardner, Cantor, and Collins (1999); Ozbudak et al. (2004); Yao et al. (2012); Lee et al. (2010); Yao (2014). Most of them are characterised by the presence of a positive feed-back in which one of the molecular species involved in the system up-regulates its own production. All of these systems exhibit bi-stability and hysteresis, i.e. a form of memory associated to bistable systems, and some of them are thought to exist in regimes where stochastic switching is frequent Acar, Becksei, and van Oudenaarden (2005); Lee et al. (2010). Noise effects on this kind of system has been extensively analysed and found to have both constructive and deleterious effects. For example, Frigola et al. Frigola et al. (2012) have found that noise stabilises the inactive (OFF) steady-state of a model of a bistable self-activation gene regulatory circuit by extending its stability region. In this Section, we analyse the effects of noise specifically associated to the quasi-steady state regime in the large-deviations (large number of molecules) limit.

Figure 6: Schematic representation of the self-activating gene regulatory circuit. The gene product is its own transcription factor which, upon dimerisation, binds the promoter region of the gene thus triggering gene transcription. The transition rates corresponding to this gene regulatory circuit are given in Table 4. For simplicity, we use an effective model in which the formation of the dimer and binding to the promoter region is taken into account in a single reaction, and the resulting number of promoter sites bound by two transcription factors is denoted .
Variable Description
Number of transcription factor molecules
Number of bound promoter sites in the gene promoter region
Number of unoccupied (unbound) binding sites in the gene promoter region
Transition rate r Event
Synthesis of the transcription factor
Degradation of the transcription factor
Dimer binding to the gene promoter region
Unbinding from the gene promoter region
Table 4: Random variables and transition rates associated to the stochastic dynamics of an auto-activation gene regulatory circuit Frigola et al. (2012); Weber and Buceta (2013). corresponds to the number of transcription-factor dimer/promoter binding site trimers. See Fig. 6 for an schematic representation.

We study the stochastic system of the simple self-activating gene regulatory circuit schematically represented in Fig. 6. In this circuit the gene product binds to form dimers which then act as its own transcription factor by binding to the promoter region of the gene. The rate-limiting factor is therefore the number of available binding sites within the promoter of the gene. For simplicity, our stochastic model associated to the rates shown in Table 4 does not explicitly account for dimer formation. We will assume that this process is very fast so it can be subsumed under the formation of transcription-factor dimer/promoter binding site trimers (reaction 3, Table 4). Furthermore, it is important to note that our stochastic dynamics exhibits a conservation law: at all time. This conservation law expresses the fact that the total number of binding sites, , is constant.

In order to proceed with our analysis of the stochastic model of self-activated gene regulation (see Table 4 and Fig. 6), we apply the general methodology associated to our SCQSS approximation. Following the general procedure explained in the previous sections, we start by deriving the stochastic Hamiltonian associated to the process defined by the transition rates shown in Table 4 (see Section II):


which, according to our theory (see Section II), gives raise to the re-scaled Hamiltonian, , defined by,


where and the re-scaled variables, , and re-scaled rate constants, , are defined in Table 5.

Rescaled variables Dimensionless parameters
Table 5: Dimensionless variables used in Eqs. (III). is a characteristic scale associated to the average number of molecules of transcription factor, , and is the average number of binding sites in the promoter of the self-activating gene. We further assume that .

The re-scaled Hamilton equations are thus given by: