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

## Abstract

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.

## 1Introduction

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 [1]. 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 [6]. For example, randomness has been shown to enhance the ability of cells to adapt and increase their fitness in random or variable environments [8]. Random noise also serves the purpose of assisting cell populations to sustain phenotypic variation by enabling cells to explore the phase space [3].

One of the mechanisms that allows noise-induced phenotypic variability relays on multi-stability [13]. The basis of this mechanism was first proposed by Kauffman [15], 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 [16].

Multi-stability is also an essential element in the control of cell response and function via signalling pathways [17]. In particular, bi-stability as a means to generate reliable switching behaviour is widely utilised in numerous pathways such as the apoptosis [18], cell survival [19], differentiation [20], and cell-cycle progression [21] 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 [23].

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) [29]. This approximation is ubiquitously used whenever regulatory processes involve enzyme catalysis, which is a central regulation mechanism in cell function [17]. 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 [22]. 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 [31]. 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. [32] 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 [33]. 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 [35]. Approaches based on enumeration techniques have also been formulated [36]. Furthermore, Thomas et al.[37] have recently formulated a rigorous method to eliminate fast stochastic variables in monostable systems using projector operators within the linear noise approximation [37]. Methods for model reduction based on perturbation analysis have been developed in [38]. 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 [40]. Several of these methods are variations of the stochastic simulation algorithm [33] or the -leap method [48] 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 [49]. Other related methods were studied in [51].

Here, we advance the formalism developed in [38], 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 [38], 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 [54]. 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 [55]. 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 [56]. 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.

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) [57]. 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 Section 3 and Section 4, respectively. Following [38], we formulate the QSS approximation for the asymptotic solution of the CME obtained by means of large deviations/WKB approximations [58]. 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 (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 [61].

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 [63]. 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 (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 (Equation 1) by and summing up over all the possible values of

From the mathematical point of view, Eq. (Equation 2) 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 [58]. This approach is based on the WKB-like Ansatz that . By substituting this Ansatz in Eq. (Equation 2) we obtain the following Hamilton-Jacobi equation for the function :

Instead of directly tackling the explicit solution of Eq. (Equation 3), we will use the so-called semi-classical approximation. We use the Feynman path-integral representation which yields a solution to Eq. (Equation 2) of the type [62]:

where indicates integration over the space of all possible trajectories and is given by [58]:

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. (Equation 4) 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 conditions[67] , , where is the initial number of molecules of species .

Eqs. (Equation 7)-( ?) are the starting point for the formulation of the semi-classical quasi-steady state approximation (SCQSSA) [38]. In order to proceed further, we assume, as per the Briggs-Haldane treatment of the Michealis-Menten model for enzyme kinetics [69], 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. (Equation 5):

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 2, , 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 6, , 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. (Equation 7)-( ?) 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. (Equation 8). The QSS approximation consists on assuming that and in Eqs. (Equation 12)-( ?),

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

## 3Bistable enzyme-catalysed systems

As a prototype of a bistable enzyme-catalysed system, we analyse a stochastic system proposed in [38], 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 [22]. Tyson & Novak [22] 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 [22] 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 Figure 1 are given in Table Table 2. 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 Figure 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 2): , , 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 2, we assume that and .

The corresponding stochastic Hamiltonian, , which is derived by applying the methodology of Section 2 to the Master Equation associated to the chemical kinetics described in Table 2, 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 2):

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

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

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 2. We first need to determine which of the variables are slow variables and which ones are fast variables. As shown in Table 3, 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:

where

with

The rescaled parameters are given in Table 3. Last, by rescaling time and defining the dimensionless time variable as (Table 3), the SCQSSA equations (Equation 11)-( ?) and (Equation 13)-( ?) lead to (see [38] 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 [38], the parameter values are determined by comparing the corresponding mean-field approximation, which is obtained by taking [67], 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 [22]. In Eq. (Equation 20) we have redefined in order to make explicit the dependence on the bifurcation parameter, , as used by Tyson & Novak [22]. The parameter values are shown in Table 4.

Upon rescaling of the variables (Table 3) and the Hamiltonian (Eq. (Equation 9)), 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. (Equation 21), , reduces to:

where, as per the SCQSSA, and are determined by Eqs. ( ?) and ( ?), 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 [38]. 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 [38], we have shown that, applying a Laplace-type asymptotic method [71] 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 [38]:

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

### 3.1Bifurcation analysis

Fig. ? shows results regarding the bifurcation behaviour of the SCQSS approximation of the stochastic bistable enzyme-catalysed system Eqs. (Equation 20)-( ?). 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. (Equation 25). 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 [54]. Figure 2 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.

Regarding the dependence on , we have checked the predictions of the SCQSS approximation by means of simulations with different values of . Figure ? shows the bi-stability region of system Eqs. (Equation 20)-( ?) 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 3 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 [73]. Finally, as continues to increase, the system becomes trapped into G-fixed point (see Figure 3). These results fully reproduce the behaviour predicted by our SCQSSA stability analysis.

The aforementioned behaviour regarding unbounded increase of fluctuations close to a bifurcation [73] is used to locate the critical value of the associated control parameter, i.e. and for the simulations shown in Figs. Figure 2 and Figure 3, 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. ?(a) (associated to the simulations shown in Figure 2), 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. ?(b) with ). The results shown in in Fig. ?(b) (corresponding to the simulations shown in Figure 3), the critical value of , , is approximately . The prediction of our asymptotic analysis (see Fig. ?(b) with ) is .

## 4Auto-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 [55] in the context of the quasi-steady regime. Many instances of genetic switches, i.e. bistable gene regulatory circuits, have been identified [75]. 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 [78]. 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. [30] 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.

We study the stochastic system of the simple self-activating gene regulatory circuit schematically represented in Figure 4. 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 6 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 6). 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 6 and Fig. Figure 4), 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 Table 6 (see Section 2):

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

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

The re-scaled Hamilton equations are thus given by:

From these equations, we observe that for , where , to hold we must have that for all . Imposing this condition on Eq. ( ?) implies that , which, in turn, together with Eqs. ( ?) and ( ?), imply that const. Finally, applying the QSS approximation to remaining equations, Eqs. (Equation 29)-( ?), we obtain:

As for the bistable enzyme-catalysed system, the parameter values are determined by matching the mean-field limit of our stochastic model, which is obtained by setting for all [67] and (i.e. the number of binding sites exactly equal to its average), to the mean-field system proposed by Frigola et al. [30]. The mapping of our parameters to those of reference [30] and their associated values are given in Table Table 8.