Noise-induced transitions in gene circuits: a perturbative approach for slow noise

Gerardo Aquino¹ Gerardo.Aquino@apha.gov.uk Andrea Rocco^2,3 A.Rocco@surrey.ac.uk ¹Department of Computing, Goldsmiths, University of London, London, UK ²Department of Physics, University of Surrey, GU2 7XH Guildford, UK ³Department of Microbial Sciences, University of Surrey, GU2 7XH, Guildford, UK

Abstract

We consider a generic class of gene circuits affected by nonlinear extrinsic noise. To address this nonlinearity we introduce a general perturbative methodology based on assuming timescale separation between noise and genes dynamics, with fluctuations exhibiting a large but finite correlation time. We apply this methodology to the case of the toggle switch, and by considering biologically relevant log-normal fluctuations, we find that the system exhibits noise-induced transitions. The system becomes bimodal in regions of the parameter space where it would be deterministically monostable. We show that by including higher order corrections our methodology allows one to obtain correct predictions for the occurrence of transitions even for not so large correlation time of the fluctuations, overcoming thereby limitations of previous theoretical approaches. Interestingly we find that at intermediate noise intensities the noise-induced transition in the toggle switch affects one of the genes involved, but not the other one.

I Introduction

Extrinsic noise has been investigated thoroughly over the last decades in a variety of physical and chemical systems Horsthemke84 . It has emerged as an active dynamical player, which can produce so-called noise-induced transitions to bimodal dynamics in systems otherwise deterministically monostable.

In chemical and physical systems, the search and the possible emergence of such transitions relies on the multiplicative nature of the corresponding stochastic dynamics, which is mathematically well characterized when the noise appears linearly in the system’s dynamical equations, and it is gaussian and white. In this case the Stratonovich prescription to interpret the corresponding stochastic integral leads to the emergence of the so-called Stratonovich drift in the Langevin and Fokker-Planck dynamics, which is responsible for dramatic changes in the stability properties of stationary states, and possibly in their number Gardiner85 .

In biological systems, however, none of the assumptions on the noise being linear, Gaussian and white is usually fulfilled. In gene regulatory networks, for instance, the noise is often non-Gaussian Cai06 ; Bengtsson05 , it is usually characterized by correlation times that can easily exceed cell cycles times Rosenfeld05 , and it appears mostly in a nonlinear fashion in the dynamics determining gene expression.

While it can be very hard to treat exactly the corresponding stochastic dynamics under these conditions, a geometric construction based on noise filtering has been introduced Ochab-Marcinek10 ; Ochab-Marcinek17 . In essence, this approach relies on the transformation properties of probability distributions of stochastic variables related by a given input/output transformation. Conservation of probability allows one to express the output probability distribution in terms of the input probability distribution, which results in non-trivial dynamics when the relation between the stochastic variables is nonlinear.

This approach has been rederived dynamically in Aquino20 by extending to nonlinear noise the so-called switching-curve approximation introduced in Arnold78 ; Reichl82 . The dynamical derivation highlights the crucial role played by the timescales of the system. The approximation works well when the input noise is much slower (in the limit, infinitely slower) than the dynamics of the system of interest. In a gene expression process, this translates in extrinsic noise being much slower that protein degradation and synthesis Aquino20 .

In this paper we take up this issue again, and develop a perturbative framework that allows one to go beyond the approximation defined in Aquino20 by extending that approach further to large but finite correlation times of the noise. We show that the inclusion of higher order perturbative terms can modify dramatically the analysis of the emergence of noise-induced transitions, by removing them when present at lower orders, and providing therefore a qualitative and quantitative more accurate description of the stochastic dynamics.

We apply our formalism to the case of the toggle switch Gardner00 , and show that when log-normal noise is considered the system undergoes a transition in a region of the phase space where the deterministic system would be monostable. We also find that tuning the noise intensity may give rise to a situation in which the marginalized probabilities of the genes involved show different qualitative features, with overlapping unimodal and bimodal behaviours. We name this feature ’partial bimodality’, and discuss its origin in the context of noise propagation across the network.

The paper is organized as follows. In Section II we present the genetic circuits considered here. In Section III we introduce our perturbative approach, and develop it up to the second perturbative order for the generic circuit introduced in Section II (with the third perturbative order described for the toggle switch in the Appendix). In Section IV we introduce the toggle switch as an application of the methodology, and assess the perturbative framework by comparing the contributions of different orders to direct stochastic simulations. In Section V we draw conclusions on the methodology introduced, and make some final remarks.

II Genetic circuits

We consider a generic circuit of interacting genes. We make the assumption that the half-life of mRNAs is much shorter than that of the corresponding proteins, so that mRNA dynamics can be adiabatically eliminated, and gene expression can be described in a single step of combined transcription and translation Rosenfeld02 . We then describe the system by decomposing the rate equations for each synthesized protein into a production term and a degradation term as follows:

\frac{d{\bf x}}{dt}={\bf g}({\bf x},R)-{\bf K}{\bf x}\ .

(1)

Here ${\bf x}$ is an $n$ -dimensional vector whose components are the $x_{i}$ species involved, while ${\bf K}$ is a diagonal matrix whose diagonal elements are the degradation rates $k_{i}=1/\tau_{x_{i}}$ of each species $x_{i}$ , with $\tau_{x_{i}}$ the respective half-lifes. The function ${\bf g}({\bf x},R)$ is a $\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$ Hill-like function of the concentrations of the $n$ species ${\bf x}$ , and $R$ is a factor that exerts control on the dynamics of ${\bf x}$ . The factor $R$ can be thought of as an external control parameter (for instance a binding constant), or any other species at high enough concentrations so that it is unaffected by the dynamics of ${\bf x}$ (for instance an external signalling molecule, or a transcription factor). For the sake of simplicity, we here limit ourselves to a single factor $R$ , but the procedure can be further generalised to multiple factors.

III Perturbative calculation

Let us now consider fluctuations acting on $R$ . We describe them by adding a generic Langevin-type dynamics to equations (1),

	$\displaystyle\frac{d{\bf x}}{dt}$	$\displaystyle=$	$\displaystyle{\bf g}({\bf x},R)-{\bf K}{\bf x}\ ,$		(2)
	$\displaystyle\frac{dR}{dt}$	$\displaystyle=$	$\displaystyle\frac{\mu(R)}{\tau}+\sqrt{\frac{D}{\tau}}\nu(R)\xi(t)\ ,$		(3)

with $\tau$ the timescale of the process $R$ . In Eq. (3) the functions $\mu(R)$ and $\nu(R)$ are kept generic for now. Their specification allows us to reproduce different noise distributions. Also, the noise $\xi(t)$ is assumed Gaussian and white, with correlator $\langle\xi(t)\xi(t^{\prime})\rangle=2\delta(t-t^{\prime})$ .

By rescaling time as $t\to t/\tau_{M}$ , with $\tau_{M}=1/k_{M}$ the largest degradation time of the species involved (we assume that $\tau_{x_{i}}\lessapprox\tau_{M}$ for $x_{i}\neq x_{M}$ ) we obtain

	$\displaystyle\frac{d{\bf x}}{dt}$	$\displaystyle=$	$\displaystyle\varepsilon{\bm{\gamma}}({\bf x},R)-{\bf\bar{K}}{\bf x}\ ,$		(4)
	$\displaystyle\frac{dR}{dt}$	$\displaystyle=$	$\displaystyle\varepsilon\mu(R)+\sqrt{\varepsilon D}\nu(R)\xi(t)\ ,$		(5)

where $\varepsilon=\tau_{M}/\tau$ , ${\bm{\gamma}}=\tau{\bf g}$ and ${\bf\bar{K}}={\bf K}/k_{M}$ with ${\bf\bar{K}}$ still diagonal and ${\bf\bar{K}}_{ii}=\kappa_{ii}=k_{x_{i}}/k_{M}=O(1)$ . In order for the dynamics of $R$ to be slower than that of all the species $x_{i}$ , we consider the condition $\tau_{M}\ll\tau$ , implying that $\varepsilon$ is a small parameter, with $\varepsilon\ll 1$ .

This naturally leads to seeking a solution of the system (4) and (5) in terms of a perturbative expansion in orders of $\varepsilon$ . This perturbative approach is inspired by the one presented in Aquino20 , with the important difference that the latter was based on assuming the limit $\tau\to\infty$ , while here we consider instead the limit $\tau_{M}\to 0$ (i.e. $\varepsilon\to 0$ ) with $\tau$ large but fixed and finite. This different perspective keeps the function ${\bm{\gamma}}=\tau{\bf g}$ in (4) fixed in the limit $\varepsilon\to 0$ , and allows us to obtain a systematic perturbative expansion in $\varepsilon$ . As a result, the approach presented here allows us to go beyond the lowest order calculation of Aquino20 , accounting for the noise filtering approach of Ochab-Marcinek10 ; Ochab-Marcinek17 , and to obtain all perturbative corrections analytically at any order in $\varepsilon$ .

Let us define the Fokker-Planck operator $L_{\rm FP}$ for the $R$ process as

L_{\rm FP}(R)=\frac{\partial}{\partial R}\left[-(\mu+D\nu\nu^{\prime})+D\frac{\partial}{\partial R}\left(\nu^{2}\right)\right]\ ,

(6)

where we have assumed the Stratonovich interpretation for the multiplicative noise equation (5).

The total Fokker-Planck equation for the time-dependent probability density $W_{t}({\bf x},R)$ is therefore:

	$\displaystyle\frac{\partial}{\partial t}W_{t}({\bf x},R)$	$\displaystyle=$	$\displaystyle\left[-\mathbf{\nabla}\cdot\left(\varepsilon{\bm{\gamma}}({\bf x},R)-{\bf\bar{K}}{\bf x}\right)\right.$		(7)
			$\displaystyle\qquad\quad\left.+\varepsilon L_{\rm FP}(R)\right]W_{t}({\bf x},R)\ .$		(7)

We now aim to compute the stationary solution of Eq. (7), namely the time-independent probability density $W({\bf x},R)$ satisfying

\displaystyle\left[\mathbf{\nabla}\cdot\left(\varepsilon{\bm{\gamma}}({\bf x},R)-{\bf\bar{K}}{\bf x}\right)-\varepsilon L_{\rm FP}(R)\right]W({\bf x},R)=0\ ,

(8)

obtained by imposing $\frac{\partial}{\partial t}W_{t}({\bf x},R)=0$ in (7). We seek a solution of the type

W({\bf x},R)=\bar{W}(R)\psi({\bf x},R)

(9)

where $\bar{W}$ is chosen so that $L_{\rm FP}\bar{W}(R)=0$ , since the dynamics of the variable $R$ is independent of $x,y$ (but not viceversa). In general it can be shown Gardiner85 that $\bar{W}(R)$ can be written as

\bar{W}(R)=\frac{\mathcal{N}}{\nu(R)}\exp\left(\frac{1}{D}\int^{R}dR^{\prime}\frac{\mu(R^{\prime})}{\nu^{2}(R^{\prime})}\right)\ ,

(10)

and it must be as well

\int{\bf dx}W({\bf x},R)=\bar{W}(R)\Rightarrow\int{\bf dx}\psi({\bf x},R)=1\ .

(11)

III.1 Zeroth order solution

We begin by assuming the following perturbative expansion for $W({\bf x},R)$ of generic order $m$ in powers of $\varepsilon$

W({\bf x},R)=\sum^{m}_{j=0}\frac{\varepsilon^{j}}{j!}W^{(j)}({\bf x},R)+O(\varepsilon^{m+1})\ .

(12)

By replacing this expansion in the lefthand side of (8), and by keeping only the zeroth order terms in $\varepsilon$ we obtain

\left(\mathbf{\nabla}\cdot{\bf x}\right)W^{(0)}({\bf x},R)=0\ ,

(13)

which immediately leads to the zeroth order solution

W^{(0)}({\bf x},R)=\delta^{n}({\bf x})\bar{W}(R)=\bar{W}(R)\prod_{i=1}^{n}\delta(x_{i})\ .

(14)

III.2 First order correction

Having evaluated $W^{(0)}$ we can replace again the expansion (12) in the lefthand side of (8), and keep this time all terms up to first order in $\varepsilon$ . We obtain:

\displaystyle-\mathbf{\nabla}\cdot\left({\bm{\gamma}}({\bf x},R)W^{(0)}-{\bf\bar{K}}{\bf x}W^{(1)}\right)+L_{\rm FP}(R)W^{(0)}=0\ ,

which, using the property $x\frac{d^{n}\delta(x)}{dx}=-n\frac{d^{n-1}\delta(x)}{dx}$ and noticing that $L_{\rm FP}(R)W^{(0)}=0$ , leads to the solution

W^{(1)}({\bf x},R)=-\bar{W}(R)\sum_{i}\left(\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}\gamma_{i}({\bf x},R)\prod_{j\neq i}\delta(x_{j})\right)\ .

(15)

Therefore, up to the first order in $\varepsilon$ we can write

	$\displaystyle W({\bf x},R)$
	$\displaystyle\quad=W^{(0)}({\bf x},R)+\varepsilon W^{(1)}({\bf x},R)$
	$\displaystyle\quad=\bar{W}(R)\left[\vphantom{\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}}\delta^{n}({\bf x})-\varepsilon\sum_{i}\left(\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}\gamma_{i}({\bf x},R)\prod_{j\neq i}\delta(x_{j})\right)\right]$
	$\displaystyle\quad\simeq\bar{W}(R)\;\delta^{n}\left({\bf x}-{\frac{\varepsilon}{\bf\bar{K}}}{\bm{\gamma}}({\bf x},R)\right)\ ,$		(16)

where the last equality is meant in the sense of distributions, and is also valid up to the first order in $\varepsilon$ :

	$\displaystyle\int d{\bf x}\delta^{n}\left({\bf x}-\varepsilon\frac{1}{\bf\bar{K}}{\bm{\gamma}}({\bf x},R)\right)h({\bf x})=\int d{\bf x}\left[\vphantom{\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}}\delta^{n}({\bf x})\right.$
	$\displaystyle-\left.\varepsilon\sum_{i}\left(\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}\gamma_{i}({\bf x},R)\prod_{j\neq i}\delta(x_{j})\right)\right]h({\bf x})+O(\varepsilon^{2})\,$		(17)

for any test function $h({\bf x})$ . We notice here that this first order result coincides with the result of the nonlinear noise filtering approach proposed in Ochab-Marcinek10 and derived dynamically in Aquino20 . In fact the delta function in Eq. (III.2) defines the (approximate) steady state solution, determined by the condition

{\bf z}({\bf x},R)={\bf x}-\varepsilon{\frac{1}{\bf\bar{K}}}{\bm{\gamma}}({\bf x},R)=0.

(18)

This solution coincides with the steady state solution of Eqs. (2) and (3) in the limit $\tau\to\infty$ of infinitely slow noise (i.e. constant $R$ on the time scale of the dynamics of ${\bf x}$ ) and corresponds trivially to a change of variable $R\to{\bf x}$ .

For instance, in the one-dimensional case, when $\gamma(x,R)=\gamma(R)$ (i.e. $g(x,R)=g(R))$ depends monotonically only on $R$ , one can use the one-dimensional version of Eq. (III.2) and after integrating over $R$ via

	$\displaystyle p(x)$	$\displaystyle=$	$\displaystyle\int dRW(x,R)=\int dR\bar{W}(R)\delta(x-\frac{\varepsilon}{\kappa}\gamma(x,R))$		(19)
			$\displaystyle=\int dR\bar{W}(R)\delta(x-\frac{\varepsilon}{\kappa}\gamma(R)),$		(19)

obtain the transformation of the probability density function

p({x})=\bar{W}(r({x}))|dr({x})/d{x}|.

(20)

In Eq. (20) $r({x})$ descends from the one-dimensional version of (18), i.e. $x=\frac{\varepsilon}{\kappa}\gamma(R)$ , and therefore $r(x)=R=\gamma^{-1}(\kappa x/\varepsilon)={g}^{-1}(kx)$ , with $\kappa=k/k_{M}$ being the only value of the one-dimensional matrix ${\bf\bar{K}}$ .

This solution corresponds to the nonlinear filtering approach introduced in Ochab-Marcinek10 ; Ochab-Marcinek17 , where the Hill-like $g$ function acts as transfer function of the regulator $R$ , assumed as an infinitely slow fluctuating variable, and produces a corresponding distribution of stationary protein concentration ${x}$ .

However, this perturbative approach can be extended to describe noise with finite $\tau$ beyond the first order, providing non-trivial corrections to the resulting dynamics as it is shown in the next subsection, as well as to generic $n$ -dimensional systems by using the implicit functions theorem.

In the multi-dimensional case, in fact, when ${\bf z}({\bf x},R)$ is a continuously differentiable function $\mathbb{R}^{n+1}\to\mathbb{R}^{n}$ , the implicit functions theorem tells us that a function $r({\bf x})=R$ exists which is locally invertible so that ${\bf z}({r^{-1}}(R),R)$ = $0$ . It follows that the probability density for ${\bf x}$ can be determined from $\bar{W}(R)$ via a change of variable as

p({\bf x})=\bar{W}(r({\bf x}))|J|.

(21)

Here $|J|$ is the Jacobian determinant of the transformation $(R=r({\bf x}),x_{2}\cdots x_{n})\to(x_{1}\cdots x_{n})$ which can be obtained, without knowing the explicit form of $r$ , from the function ${\bf z}$ in Eq. Eq. (18) using the implicit function theorem in $n$ dimensions. Eq. (21) generalises Eq. (20) to the case of $n$ species, with $r({\bf x})$ now depending on all species $x_{1}\cdots x_{n}$ and defined only implicitly via Eq. (18).

III.3 Second order correction

In order to evaluate the second order correction we need first to evaluate $L_{\rm FP}(R)W^{(1)}({\bf x},R)$ . From Eq. (15) and the property $L_{\rm FP}(R)\bar{W}(R)=0$ it follows that

	$\displaystyle-L_{\rm FP}(R)W^{(1)}({\bf x},R)$
	$\displaystyle\qquad=\sum_{i}\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}\prod_{k\neq i}\delta(x_{k})L_{\rm FP}(R)\left[\gamma_{i}({\bf x},R)\bar{W}(R)\right]$
	$\displaystyle\qquad=\sum_{i}\frac{\delta^{\prime}(x_{i})}{\kappa_{ii}}G_{i}({\bf x},R)\prod_{k\neq i}\delta(x_{k})$		(22)

where $G_{i}({\bf x},R)=L_{\rm FP}(R)\left[\gamma_{i}({\bf x},R)\bar{W}(R)\right]$ which, using Eq.(6), leads to:

	$\displaystyle G_{i}({\bf x},R)$
	$\displaystyle=-(\mu(R)-3D\nu(R)\nu^{\prime}(R))\bar{W}(R)\partial_{R}\gamma_{i}({\bf x},R)$
	$\displaystyle\quad+D\nu(R)^{2}(2\bar{W}^{\prime}(R)\partial_{R}\gamma_{i}({\bf x},R)+\bar{W}(R)\partial_{R}^{2}\gamma_{i}({\bf x},R)\ .$		(23)

Having evaluated $W^{(1)}({\bf x},R)$ we can now repeat the same procedure as done for the first order, i.e. replace expansion (12) in the lefthand side of (8) and keep this time all terms up to second order in $\varepsilon$ . We obtain:

-\mathbf{\nabla}\!\cdot\!\left({\bm{\gamma}}({\bf x},R)W^{(1)}-{\bf\bar{K}}{\bf x}\frac{W^{(2)}}{2}\right)+L_{\rm FP}(R)W^{(1)}({\bf x},R)=0\ ,

(24)

which, using (III.3), can be rewritten as

	$\displaystyle\mathbf{\nabla}\cdot\left({\bm{\gamma}}({\bf x},R)W^{(1)}-{\bf\bar{K}}{\bf x}\frac{W^{(2)}}{2}\right)$
	$\displaystyle\qquad\qquad+\sum_{i}\frac{\delta^{\prime}(x_{i})}{k_{ii}}G_{i}({\bf x},R)\prod_{k\neq i}\delta(x_{k})=0\ .$		(25)

From the structure of this equation and from (III.2) we deduce that the solution for $W^{(2)}({\bf x},R)$ must have the following form:

	$\displaystyle W^{(2)}({\bf x},R)$	$\displaystyle=$	$\displaystyle\bar{W}(R)\sum_{i,q}\left[\left(A_{iq}\delta^{\prime}(x_{i})+B_{iq}\delta^{\prime\prime}(x_{i})\right)\delta(x_{q})\right.$		(26)
			$\displaystyle\left.\qquad+C_{iq}\delta^{\prime}(x_{i})\delta^{\prime}(x_{q})\right]\prod_{k\neq i,q}\delta(x_{k})\ .$		(26)

The coefficients $A_{iq}$ , $B_{iq}$ , $C_{iq}$ can be obtained by inserting (26) in (III.3) and result in


$\displaystyle A_{iq}$	$\displaystyle=\left(-\frac{2}{k_{ii}^{2}}\bar{G}_{i}({\bf x},R)+\tilde{A}_{i}({\bf x},R)\right)\delta_{iq}$	(27a)
	$\displaystyle\text{with}\;\;\frac{\tilde{A}_{i}({\bf x},R)}{x_{i}}-\partial_{x_{i}}\tilde{A}_{i}({\bf x},R)=-\frac{2}{\kappa_{ii}^{2}}\partial_{x_{i}}\bar{G}_{i}({\bf x},R)$
$\displaystyle B_{iq}$	$\displaystyle=\left(\frac{\gamma_{i}({\bf x},R)}{\kappa_{ii}}\right)^{2}\delta_{iq}$	(27b)
$\displaystyle C_{iq}$	$\displaystyle=\frac{\gamma_{i}({\bf x},R)\gamma_{q}({\bf x},R)}{\kappa_{qq}\kappa_{ii}}(1-\delta_{iq})$	(27c)

with $\bar{G}_{i}=G_{i}/\bar{W}(R)$ . Therefore we make the following ansatz for the stationary probability density $W({\bf x},R)$ :

W({\bf x},R)=\bar{W}(R)\delta^{n}\left({\bf x}-\frac{\varepsilon}{\bf\bar{K}}{\bm{\gamma}}({\bf x},R)-\frac{\varepsilon^{2}}{2}{\bf Q}({\bf x},R)\right){+O(\varepsilon^{3})}

(28)

with

{\bf Q}({\bf x},R)=\left(-\frac{2}{\bf\bar{K}^{2}}{\bf\bar{G}}({\bf x},R)+{\bf\tilde{A}}({\bf x},R)\right).

(29)

The solution given by the ansatz (28) in fact coincides up to second order in $\varepsilon$ with the solution derived with the perturbative expansion (12) to the same order, i.e. including Eq. (26) with coefficients (27).

The delta function in Eq. (28) is similar to the one in (III.2) and therefore, as it was seen for the first order in (18), it defines the (approximate) steady state solution of Eqs. (2) and (3). This is determined by the condition

{\bf z}({\bf x},R)={\bf x}-\frac{\varepsilon}{\bf\bar{K}}{\bm{\gamma}}({\bf x},R)-\frac{\varepsilon^{2}}{2}{\bf Q}({\bf x},R)=0,

(30)

which in the case that $\mathbf{\nabla}\cdot\frac{1}{\bf\bar{K}}{\bm{\gamma}}({\bf x},R)=const$ , implying in turn $\mathbf{\nabla}\cdot\frac{1}{\bf\bar{K}}{\bf G}({\bf x},R)=0$ and therefore (neglecting null contributions in distributional sense) ${\bf\tilde{A}}({\bf x},R)=0$ , simplifies to

{\bf z}({\bf x},R)={\bf x}-\frac{\varepsilon}{\bf\bar{K}}{\bm{\gamma}}({\bf x},R)-\varepsilon^{2}\frac{1}{\bf\bar{K}^{2}}{\bf\bar{G}}({\bf x},R)=0.

(31)

Note that both ${\bf Q}({\bf x},R)$ and ${\bf\bar{G}}({\bf x},R)$ depend linearly on ${\bm{\gamma}}=\tau{\bf g}$ . Therefore, when multiplied by their expansion factor $\varepsilon^{2}=\frac{\tau_{M}^{2}}{\tau^{2}}$ , they disappear in the limit $\tau\to\infty$ , leaving only the first order correction in Eqs. (28), (30) or (31). This can be verified to apply to all the higher order corrections in the same limit. Therefore Eq. (18), obtained as first order correction within our approach, is actually an exact solution in such limit, corresponding thereby to the noise filtering approach of Ochab-Marcinek10 ; Ochab-Marcinek17 .

Under the assumption that ${\bf z}({\bf x},R)$ is a continuously differentiable function $\mathbb{R}^{n+1}\to\mathbb{R}^{n}$ , the implicit functions theorem can be applied again. As a result, the probability density for ${\bf x}$ can be determined from $\bar{W}(R)$ via a change of variable as shown in Eq. (21), where the Jacobian can now be obtained from the function ${\bf z}$ in Eq. (30) or (31).

In the following section, we show how this approach is applied with a concrete example.

IV The toggle switch

Refer to caption — Figure 1: The toggle switch. The expression level of gene $x$ is downregulated by gene $y$ upon binding with the repressor R. We make the assumption that when active, both genes synthesize the proteins $x$ and $y$ in one single step of combined transcription and translation.

As an application of our methodology we consider the toggle switch, a positive feedback loop of two genes mutually repressing each other (Fig. 1) Gardner00 . The dynamic equations describing each of the synthesized protein have the form (1), namely:

	$\displaystyle\frac{dx}{dt}$	$\displaystyle=$	$\displaystyle g(y,R)-k_{1}x=\frac{g_{1}}{1+\rho_{1}{(}Ry{)}^{\beta_{1}}}-k_{1}x$		(32)
	$\displaystyle\frac{dy}{dt}$	$\displaystyle=$	$\displaystyle f(x)-k_{2}y=\frac{g_{2}}{1+\rho_{2}x^{\beta_{2}}}-k_{2}y$		(33)

Here, we have assumed that protein synthesis occurs in one single stage of merged transcription and translation, which, as discussed, is valid when the half-life of mRNA species is much shorter than the half-life of protein species Rosenfeld02 . Further to standard degradation dynamics for $x$ and $y$ proteins, characterized by degradation rates $k_{1}$ and $k_{2}$ respectively, the structure of the equations is based on two Hill-like functions, which account for each of the two synthesised proteins to act as a repressor on the other gene. Here parameters $\rho_{1,2}$ are the association constants for the binding of $x$ ( $y$ , respectively) to the promoter of the target gene, $g_{1,2}$ are the maximal expression rates of each gene, and $\beta_{1,2}$ are cooperativity parameters representing homodimerization of the relevant protein.

The bifurcation diagram of this system has been analyzed in detail in Gardner00 . The system exhibits deterministic multistable behavior for $\beta_{1},\beta_{2}>1$ , and it is otherwise monostable. To probe the effect of noise in a situation of deterministic monostability, we therefore set in the following $\beta_{1}=\beta_{2}=1$ . In this way any transition to bimodal dynamics will only be due to the effect of noise, accounting thereby for the possible emergence of a noise-induced transition in the system.

To explore this further, we have made the assumption that regulation of $x$ happens upon formation of a heterodimer of $y$ with the external (repressive) transcription factor $R$ . Hence fluctuations in $R$ will correspond to insertion of nonlinear extrinsic noise in the system, and will be amenable to be studied with the methodology developed here.

IV.1 Perturbative calculation

As for the more general case, let us add Eq. (3) to the equations (32) and (33) for the toggle switch and let us assume $\tau_{x}=1/k_{1}>\tau_{y}=1/k_{2}$ . We then redefine time as $t\to t/\tau_{x}$ , we obtain:

$\displaystyle\frac{dx}{dt}$	$\displaystyle=$	$\displaystyle\varepsilon\gamma(y,R)-x=\varepsilon\frac{\gamma_{1}}{1+\rho_{1}Ry}-x\ ,$	(34)
$\displaystyle\frac{dy}{dt}$	$\displaystyle=$	$\displaystyle\varepsilon f(x)-\kappa y=\varepsilon\frac{\gamma_{2}}{1+\rho_{2}x}-\kappa y\ ,$	(35)
$\displaystyle\frac{dR}{dt}$	$\displaystyle=$	$\displaystyle\varepsilon\mu(R)+\sqrt{\varepsilon D}\nu(R)\xi(t)\ ,$	(36)

where $\varepsilon=\tau_{x}/\tau$ , $\gamma_{1,2}=g_{1,2}\tau$ and $\kappa=k_{2}/k_{1}\simeq 1$ . For simplicity we here assume $\kappa=1$ .

We also assume that the timescale $\tau$ of the dynamics of $R$ satisfies $\tau\gg\tau_{x},\tau_{y}$ and therefore $\varepsilon$ is a small parameter. We are in this way assuming that the dynamics of $R$ is slower than that of $x$ and $y$ .

With $L_{\rm FP}$ defined as in Eq. (6), the Fokker-Planck (7) becomes

	$\displaystyle\frac{\partial}{\partial t}W_{t}(x,y,R)$	$\displaystyle=$	$\displaystyle\left[-\frac{\partial}{\partial x}\left(\varepsilon\gamma(y,R)-x\right)-\frac{\partial}{\partial y}\left(\varepsilon f(x)-y\right)\right.$		(37)
			$\displaystyle\left.\vphantom{-\frac{\partial}{\partial x}}+\varepsilon L_{\rm FP}(R)\right]W_{t}(x,y,R)\ ,$		(37)

whose stationary solution satisfies

	$\displaystyle\left[\frac{\partial}{\partial x}\left(x-\varepsilon\gamma(y,R)\right)+\frac{\partial}{\partial y}\left(y-\varepsilon f(x)\right)+\right.$
	$\displaystyle\left.\vphantom{\frac{\partial}{\partial x}}+\varepsilon L_{\rm FP}(R)\right]W(x,y,R)=0$		(38)

IV.2 Zeroth order solution

Assuming the perturbative expansion (12), Eq. (13) becomes

\left(\frac{\partial}{\partial x}x+\frac{\partial}{\partial y}y\right)W^{(0)}(x,y,R)=0\ ,

(39)

which immediately leads to the zeroth order solution:

W^{(0)}(x,y,R)=\delta(x)\delta(y)\bar{W}(R)\ .

(40)

IV.3 First order correction

We now proceed with the calculation of first order corrections $W^{(1)}(x,y,R)$ . Upon substitution of (12) in (IV.1) and by collecting terms up to first order in $\varepsilon$ , we obtain

	$\displaystyle-\frac{\partial}{\partial x}\left(\gamma(y,R)W^{(0)}-xW^{(1)}\right)$	$\displaystyle-$	$\displaystyle\frac{\partial}{\partial y}\left(f(x)W^{(0)}-yW^{(1)}\right)$		(41)
		$\displaystyle+$	$\displaystyle L_{\rm FP}(R)W^{(0)}=0$		(41)

which leads to the solution

	$\displaystyle W^{(1)}(x,y,R)$
	$\displaystyle\hskip 2.84544pt=-\bar{W}(R)\left(\vphantom{\sum}\delta^{\prime}(x)\delta(y)\gamma(y,R)+\delta^{\prime}(y)\delta(x)f(x)\right),$		(42)

as could be derived directly from Eq.(15). Following our general procedure, we can write to the first order in $\varepsilon$ :

$\displaystyle W(x,y,R)$	$\displaystyle=$	$\displaystyle W^{(0)}+\varepsilon W^{(1)}$	(43)
	$\displaystyle=$	$\displaystyle\bar{W}(R)\left[\left(\delta(x)-\varepsilon\gamma(y,R)\delta^{\prime}(x)\right)\delta(y)\right.$
		$\displaystyle\qquad\qquad\quad\;\;-\left.\delta^{\prime}(y)\delta(x)f(x)\right]$
	$\displaystyle=$	$\displaystyle\bar{W}(R)\delta(x-\varepsilon\gamma(y,R))\delta(y-\varepsilon f(x))\ ,$

analogous to Eq. (III.2). In fact, if we let

\tilde{\delta}(x,y)=\left(\vphantom{\sum}\delta(x)-\varepsilon\gamma(y,R)\delta^{\prime}(x)\right)\delta(y)-\varepsilon\delta^{\prime}(y)\delta(x)f(x)\ ,

(44)

then from the distributional point of view

	$\displaystyle\iint dxdy\tilde{\delta}(x,y)h(x,y)$
	$\displaystyle\hskip 5.69046pt=h(0,0)+\varepsilon\left(\gamma(y,R)\frac{\partial}{\partial x}h(0,0)+f(x)\frac{\partial}{\partial y}h(0,0)\right)$
	$\displaystyle\hskip 5.69046pt=h(\varepsilon\gamma(y,R),\varepsilon f(x))+O(\varepsilon^{2})$		(45)

and since

	$\displaystyle h(\varepsilon\gamma(y,R),\varepsilon f(x))$
	$\displaystyle\hskip 8.5359pt=\iint dxdy\delta(x-\varepsilon\gamma(y,R))\delta(y-\varepsilon f(x))h(x,y)$		(46)

for any $h(x,y)$ , the last equality in (43) follows. This result coincides at this order with the noise filtering approach proposed in Ochab-Marcinek10 with $p(x)$ stemming from $\bar{W}(R)$ via the change of variables $R\to x$ induced by the (analytically solvable in this case) condition $x-\varepsilon\gamma(\varepsilon f(x),R)=0$

IV.4 Second order correction

We can now carry out the evaluation of the second order correction by computing $L_{\rm FP}(R)W^{(1)}$ . Eq. (III.3) becomes

	$\displaystyle L_{\rm FP}(R)W^{(1)}(x,y,R)$	$\displaystyle=$	$\displaystyle-\delta^{\prime}(x)\delta(y)L_{\rm FP}(R)\left(\gamma(y,R)\bar{W}(R)\right)$		(47)
		$\displaystyle=$	$\displaystyle-\delta^{\prime}(x)\delta(y)G(y,R)\ ,$		(47)

where

$\displaystyle G(y,R)$	$\displaystyle=$	$\displaystyle L_{\rm FP}(R)\left(\gamma(y,R)\bar{W}(R)\right)$
	$\displaystyle=$	$\displaystyle-(\mu(R)-3D\nu(R)\nu^{\prime}(R))\bar{W}(R)\partial_{R}\gamma(y,R)+$
	$\displaystyle+$	$\displaystyle D\nu(R)^{2}(2\bar{W}^{\prime}(R)\partial_{R}\gamma(y,R)+\bar{W}(R)\partial_{R}^{2}\gamma(y,R)).$

Eq. (24) now becomes

	$\displaystyle-\frac{\partial}{\partial x}\left(\gamma(y,R)W^{(1)}-\frac{x}{2}W^{(2)}\right)$	$\displaystyle-$	$\displaystyle\frac{\partial}{\partial y}\left(f(x)W^{(1)}-\frac{y}{2}W^{(2)}\right)$		(49)
		$\displaystyle+$	$\displaystyle L_{\rm FP}(R)W^{(1)}=0\ ,$		(49)

which, using (47), can be rewritten as

	$\displaystyle-\frac{\partial}{\partial x}\left(\gamma(y,R)W^{(1)}-\frac{x}{2}W^{(2)}\right)$	$\displaystyle-$	$\displaystyle\frac{\partial}{\partial y}\left(f(x)W^{(1)}-\frac{y}{2}W^{(2)}\right)$		(50)
		$\displaystyle+$	$\displaystyle\delta^{\prime}(x)\delta(y)G(y,R)=0\ .$		(50)

The form of the solution $W^{(2)}(x,y,R)$ is analogous to (26), namely given by

	$\displaystyle W^{(2)}(x,y,R)$	$\displaystyle=$	$\displaystyle\bar{W}(R)\left(A_{10}\delta^{\prime}(x)\delta(y)+A_{11}\delta^{\prime}(x)\delta^{\prime}(y)\right.$		(51)
		$\displaystyle+$	$\displaystyle\left.A_{20}\delta^{\prime\prime}(x)\delta(y)+A_{02}\delta(x)\delta^{\prime\prime}(y)\right)$		(51)

Inserting this expression in (50) one obtains the following conditions on the coefficients $A_{ij}$

	$\displaystyle A_{10}$	$\displaystyle=$	$\displaystyle-2\bar{G}(y,R)\ ,\;\;\;A_{20}=\left(\gamma(y,R)\right)^{2}\ ,$
	$\displaystyle A_{11}$	$\displaystyle=$	$\displaystyle 2\gamma(y,R)f(x)\ ,\;\;\;A_{02}=\left(f(x)\right)^{2}\ .$

with $\bar{G}=G/\bar{W}$ , as could be derived as well directly from Eqs. (27) with the condition ${\bf\tilde{A}}=0$ . The ansatz for the steady state solution for the probability density $W(x,y,R)$ that coincides up to the second order in $\varepsilon$ with the perturbative evaluation shown above is therefore

	$\displaystyle W(x,y,R)$		$\displaystyle=\bar{W}(R)\delta\left(x-\varepsilon\gamma(y,R)-\varepsilon^{2}\bar{G}(y,R)\right)$		(52)
			$\displaystyle\times\delta\left(y-\varepsilon f(x)\right)+O(\varepsilon^{3}).$		(52)

IV.5 Log-normal Noise

The power of the approach that we propose here relies also on the fact that we can assume different types of noise, as long as these can be represented in terms of the Langevin dynamics (3). By following Aquino20 , we focus on log-normal noise, which in contrast to Gaussian noise preserves the positivity of parameter values on which it acts. It also describes well extrinsic fluctuations Rosenfeld05 ; Shahrezaei08 and universal behaviours in bacteria and yeast Salman12 ; Brenner15 .

The stationary log-normal distribution

\bar{w}(R)=\frac{1}{R\sqrt{2\pi}}\exp{\left[-\frac{1}{2D}\left(\log\frac{R}{\bar{R}}+\frac{D}{2}\right)^{2}\right]}

(53)

can be derived Aquino20 by defining the stochastic process

R(t)=\bar{R}e^{\eta(t)}e^{-D/2},

(54)

where $\eta(t)$ is the standard Ornstein-Uhlenbeck noise,

\frac{d\eta}{dt}=-\frac{\eta}{\tau}+\sqrt{\frac{D}{\tau}}\xi(t),

(55)

with $\xi(t)$ zero average Gaussian white noise, and $\bar{R}=\langle R(t)\rangle$ , since $\langle e^{\eta(t)}\rangle=e^{D/2}$ . Hence, the stochastic process (54) obeys

\frac{dR}{dt}=-\frac{R}{\tau}\log{\left(\frac{R}{\bar{R}}+\frac{D}{2}\right)}+\sqrt{\frac{D}{\tau}}R\xi(t),

(56)

which corresponds to (3) with the choice $\mu(R)=-R\log(R/\bar{R}+D/2)$ and $\nu(R)=R$ . Eq. (10) allows then the derivation of (53).

IV.6 Comparison with numerical simulations

The joint probability density $p(x,y)$ for the protein distributions can be obtained by marginalizing (52) with respect to $R$ :

p(x,y)=\int dRW(x,y,R)\ .

(57)

The probability density $p_{x}(x)$ for the distribution of protein $x$ can then be obtained by marginalizing (57) with respect to $y$ . In this case it is simpler to first integrate in $y$ after replacing (52) in (57), we obtain

	$\displaystyle p_{x}(x)$	$\displaystyle=$	$\displaystyle\int dR\bar{W}(R)\delta\left(x-\varepsilon\gamma(\varepsilon f(x),R)-\varepsilon^{2}\bar{G}(\varepsilon f(x),R)\right)$		(58)
		$\displaystyle=$	$\displaystyle\int dR\bar{W}(R)\delta(z(x,R))$		(58)

and further integration with respect to $R$ amounts to a change of variable $R=r(x)$ where the function $r(x)$ is defined implicitly by the equation

z(x,R)=x-\varepsilon\gamma(\varepsilon f(x),R)-\varepsilon^{2}\bar{G}(\varepsilon f(x),R)=0\ .

(59)

We notice that in this particular case we can also develop easily the third order correction. This is obtained following the same procedure used for the second order (see Appendix for details) leading to the following modified expression for the steady state solution for the probability density $W(x,y,R)$ :

	$\displaystyle W(x,y,R)=\bar{W}(R)\delta\left(x-\varepsilon\gamma(y,R)-\varepsilon^{2}\bar{G}(y,R)+\right.$
	$\displaystyle\left.-\varepsilon^{3}\bar{G_{3}}(x,y,R)\right)\delta\left(y-\varepsilon f(x)\right)+O(\varepsilon^{4})\ .$		(60)

Here, we have used (see Appendix):

\bar{G_{3}}(x,y,R)=\bar{G}_{M}(y,R)+\frac{\bar{G}_{2}(y,R)}{2x}\ ,

(61)

and therefore by integrating over $y$ we obtain

p_{x}(x)=\int dR\bar{W}(R)\delta(z(x,R))

(62)

with

	$\displaystyle z(x,R)$	$\displaystyle=$	$\displaystyle x-\varepsilon\gamma(\varepsilon f(x),R)-\varepsilon^{2}\bar{G}(\varepsilon f(x),R)+$
		$\displaystyle-$	$\displaystyle\varepsilon^{3}\bar{G_{3}}(x,\varepsilon f(x),R),$

which replaces (59).

According to the implicit functions theorem, from the condition $z(x,R)=0$ , it follows that a function $r(x)=R$ exists such that

p_{x}(x)=\bar{W}(r(x))\left|\frac{dr}{dx}\right|\ ,

(64)

where $dr/dx$ can be computed by applying the implicit functions theorem either to (59) for the calculation up to the second order, or to (IV.6) for the third order. Hence:

\frac{dr}{dx}=-\frac{\partial_{x}z(x,R)}{\partial_{R}z(x,R)}\ ,

(65)

which allows us to evaluate explicitly $p_{x}(x)$ for the toggle switch.

Finally, in order to derive

p_{y}(y)=\int dx\int dRW(x,y,R)\ ,

(66)

it suffices to notice that at steady state we have exactly $y=\varepsilon f(x)$ , and therefore $p_{y}(y)$ in this case can be derived directly from $p_{x}(x)$ via a change of variable as

	$\displaystyle p_{y}(y)$	$\displaystyle=$	$\displaystyle\frac{1}{\varepsilon}p_{x}\left(f^{-1}\left(\frac{y}{\varepsilon}\right)\right)\left\|\frac{df^{-1}}{dy}\right\|$		(67)
		$\displaystyle=$	$\displaystyle\frac{p_{x}(f^{-1}(y/\varepsilon))}{\varepsilon\left\|f^{\prime}(x=f^{-1}(y/\varepsilon)\right\|}\ .$		(67)

Equations (64) and (67) so constructed can then be compared at any given perturbative order with direct stochastic simulations, in order to highlight the predictive power of the theoretical approach here developed.

We perform direct stochastic simulations of fluctuatons in $R$ by using the standard approach to generate an Ornstein-Uhlenbeck process, converting it via Eq. (54) into log-normal noise and then by performing ordinary integration of the equations (32) and (33). The probability densities are obtained by averaging a single trajectory over long times. Parameters are set so as to reproduce biologically relevant conditions. In all simulations and in the analytical solution, we fix the maximal gene expression rates at the values $g_{1}=g_{2}=g_{0}=10^{-2}s^{-1}$ , which is representative of typical gene expression in bacteria Thattai01 , and the degradation rates $k_{1}=k_{2}=10^{-4}s^{-1}$ , which are representative of proteins’ half-life in the range of 2 hours. We also fix the dissociation constants of transcription factors to DNA $\rho_{1}=10nM^{-1}$ and $\rho_{2}=1nM^{-1}$ , and $\bar{R}=0.87nM$ Aquino20 .

We then focus on two different values of $\varepsilon$ in order to explore the contribution of different perturbative orders. Since we assume $k=10^{-4}s^{-1}$ , choosing $\tau=10^{7}s$ corresponds to assume $\varepsilon=0.001$ , while the choice $\tau=10^{5}s$ will reflect a larger value of $\varepsilon$ , namely $\varepsilon=0.1$ . It is this second regime, when the timescale separation between dynamics of $x$ and of the $R$ fluctuations is less pronounced, that we predict that higher order terms in the perturbative expansion will become relevant.

In Fig. 2 we show the comparison between perturbative expansion and stochastic simulations for $\varepsilon=0.001$ . The theoretical analysis predicts bimodality for the concentration $x$ , and this is confirmed excellently by the numerical result. In this case first order solution Eq. (43), second order solution via Eq. (59), and third order solution via Eq. (IV.6) practically coincide.

In Fig. 3 (top panel) we show a particular case where using the first order solution Eq. (43) leads erroneously to conclude that bimodality emerges, while including second and third order correction gives a quantitatively and qualitatively correct description of the properties of the density function, correctly predicting a unimodal distribution for both $x$ and $y$ protein concentrations.

IV.7 Noise-induced transitions: Partial bimodality

The analysis and simulations shown in Fig. 2 highlight the emergence of a particular phenomenon, which we name ’partial bimodality’, whereby bimodal behaviour emerges for the protein concentration $x$ , but not for $y$ . This is clearly unexpected in standard toggle switches affected by intrinsic noise only, as in this case the bimodality would instead affect simultaneosly both genes involved.

In order to obtain further insight into these dynamics, and clarify the nature of the noise-induced transition here identified, we perform simulations and solve the corresponding first order dynamics, for three different values of noise intensities.

In Fig. 4 we show simulations and analytical predictions for the toggle switch in a regime where at low noise intensity ( $D=0.01$ ) the dynamics are unimodal. In this case we observe that increasing the noise intensity up to $D=10$ bimodal behaviour emerges for both genes. However, interestingly we also observe that for intermediate values of noise, $D\simeq 0.1$ , a regime exists where bimodality emerges only for one gene and not the other.

For this case we show two exemplary trajectories for variables $x$ and $y$ in the upper panel of Fig. 5. While the variable $x$ makes complete transitions between a high value around the second mode ( $x=32.5$ ) and low values close to the first mode ( $x=0.24)$ , the variable $y$ remains more consistently around $y=0.2$ , as testified by the higher peak of its distribution at the origin as compared to that of $p_{x}$ (see inset in middle Figure 4). In correspondence of transitions of the $x$ variable to the low state, excursions to slightly larger values of $y$ are permitted, but their amplitude is not large enough to exit the basin of attraction of the low $y$ state. In the bottom panel of Fig. 5 instead, the case where bimodality occurs for both species is shown. These traejctories correspond to the bimodal distributions in the bottom panel of Fig. 4 with modes in $x,y=0$ and in $x=100,y=99$ . Since $p_{x}$ has a higher peak in $x=100$ and $p_{y}$ has a higher peak in $y=0$ these species spend most of their time around these values with sporadic transitions to values around $x=0$ and $y=99$ respectively. The amplitude of the fluctuations of the $y$ variable allows in this case the variable $y$ to exit the basin of attraction of its high state.

We attribute the different noise amplitudes experienced by the variables $x$ and $y$ to a noise propagation effect, whereby the fluctuations of the external parameter $R$ are perceived differently by the two genes. In fact, the variance of the distribution of $x$ is consistently larger of the variance of the distribution of $y$ , for all input noises on $R$ (parameterized by $D$ ), as shown in Table 1:

$D$	$0.01$	$0.1$	$10$
$\sigma_{x}$	56.62	325.31	590.67
$\sigma_{y}$	0.88	58.57	191.17

Table 1: Variances of

x

and

y

distributions corresponding to the probability distributions shown in Fig. 4. Parameter values are set as in Fig. 4.

While $x$ is under the direct influence of fluctuations of $R$ , gene $y$ feels those fluctuations, propagated through the repressive action from $x$ , at a lower level. This means that while gene $x$ can be above threshold for transitions to bimodal behavior, gene $y$ may still remain below the threshold, and therefore may still exhibit unimodal behaviour. This mechanism is in agreement with the result of some extensive studies performed in Hornung08 , which demonstrated that positive feedback loops limit the range of fluctuations of the target gene.

V Conclusions

In this paper we have introduced an efficient methodology to deal with nonlinear extrinsic noise when the correlation time of the fluctuations is large. The derivation here proposed is based on a timescale expansion of the master equation and reproduces to the first order the phenomenological approach based on the nonlinear noise filtering introduced in Ochab-Marcinek10 . Importantly, our perturbative approach allows one to go beyond the dynamical derivation presented in Aquino20 and obtain correct predictions even in a regime when the timescales characterizing fluctuations and gene dynamics are comparable, with the lowest perturbative orders failing to capture the qualitative dynamical features of the system.

We also find that the toggle switch is capable of a partial bimodality, in that while one of the two participating genes shows full bimodal behaviour, the other gene still preserves unimodality. This feature is intermediate in terms of noise intensities to the two extreme cases of no bimodality for any of the two genes, to full bimodality for both of them. As a follow-up to this observation, it will be interesting to analyze the effect of simultaneous intrinsic and extrinsic noise on the system. We expect that the interplay between intrinsic and extrinsic noise would provoke even richer dynamics, in which the balance between noise intensities will be crucial in determining the structure of the phase space.

The novel mechanism for noise-induced transitions here described and the methodology introduced are promising and worthy of further analysis and experimental validation. They will further our understanding of the fundamental and active role played by noise in biological systems.

ACKNOWLEDGMENTS

This study has been funded by the UK Biotechnology and Biological Sciences Research Council (BBSRC), under grant no. BB/L007789/1.

For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

DATA AVAILABILITY

All details of the mathematical models developed and the data deriving from their simulation are available in the article.

APPENDIX: EVALUATION OF THIRD ORDER PERTURBATIVE CORRECTION FOR THE TOGGLE SWITCH

We can evaluate the third order correction in the perturbative expansion just by carrying on the procedure introduced in the main text to include the third order. We start with $W(x,R)=W^{(0)}+\varepsilon W^{(1)}+\varepsilon^{2}W^{(2)}/2+\varepsilon^{3}W^{(3)}/6$ and this leads to

	$\displaystyle-\frac{\partial}{\partial x}\left(\gamma(y,R)W^{(2)}-\frac{xW^{(3)}}{3}\right)$
	$\displaystyle-\frac{\partial}{\partial y}\left(f(x)W^{(2)}-\frac{yW^{(3)}}{3}\right)+L_{\rm FP}(R)W^{(2)}=0.$		(68)

where we know $W^{(2)}(x,y,R)$ from previous calculations. Therefore we can evaluate $L_{\rm FP}(R)W^{(2)}$ , and we obtain

	$\displaystyle L_{\rm FP}(R)W^{(2)}(x,y,R)=G_{2}(y,R)\delta^{\prime\prime}(x)\delta(y)$
	$\displaystyle+2f(x)G(y,R)\delta^{\prime}(x)\delta^{\prime}(y)-2G_{M}(y,R)\delta^{\prime}(x)\delta(y)$		(69)

with

	$\displaystyle G_{2}(y,R)=L_{\rm FP}\left(\gamma^{2}(y,R)\bar{W}(R)\right)=$
	$\displaystyle\qquad 2(3D\nu(R)\nu^{\prime}(R)-\mu(R))\bar{W}(R)\gamma(y,R)\partial_{R}\gamma(y,R)$
	$\displaystyle\qquad+D\nu(R)^{2}\left[4\bar{W}^{\prime}(R)\gamma(y,R)\partial_{R}\gamma(y,R)\right.$
	$\displaystyle\qquad\left.+2\bar{W}(R)\left(\gamma(y,R)\partial_{R}^{2}\gamma(y,R)+(\partial_{R}\gamma(y,R))^{2}\right)\right]\ ,$		(70)

where it can be shown that

	$\displaystyle G_{M}(y,R)=L_{\rm FP}(R)\left(G(y,R)\right)=$
	$\displaystyle\alpha(R)G(y,R)+\sigma(R)\partial_{R}G(y,R)+\gamma(R)\partial_{R}^{2}G(y,R)\ ,$		(71)

with

$\displaystyle\alpha(R)$	$\displaystyle=$	$\displaystyle-\left(\mu^{\prime}(R)-D\nu^{\prime}(R)^{2}-D\nu(R)\nu^{\prime\prime}(R)\right)\ ,$	(72)
$\displaystyle\sigma(R)$	$\displaystyle=$	$\displaystyle-\left.(\mu(R)-3D\nu(R)\nu^{\prime}(R)\right)\ ,$	(73)
$\displaystyle\gamma(R)$	$\displaystyle=$	$\displaystyle D\nu^{2}(R)\ .$	(74)

From Eq. (APPENDIX: EVALUATION OF THIRD ORDER PERTURBATIVE CORRECTION FOR THE TOGGLE SWITCH) we can therefore expect $W_{3}(x,R)$ of the form:

	$\displaystyle W_{3}(x,y,R)=\bar{W}(R)\left(A_{10}\delta^{\prime}(x)\delta(y)+A_{11}\delta^{\prime}(x)\delta^{\prime}(y)\right.$
	$\displaystyle\qquad\left.+A_{12}\delta^{\prime}(x)\delta^{\prime\prime}(y)+A_{20}\delta^{\prime\prime}(x)\delta(y)+A_{21}\delta^{\prime\prime}(x)\delta^{\prime}(y)+\right.$
	$\displaystyle\qquad\left.+A_{30}\delta^{\prime\prime\prime}(x)\delta(y)+A_{03}\delta(x)\delta^{\prime\prime\prime}(y)\right)\ .$		(75)

Replacing Eq. (75) in (APPENDIX: EVALUATION OF THIRD ORDER PERTURBATIVE CORRECTION FOR THE TOGGLE SWITCH) and solving for $A_{ij}$ leads to

$\displaystyle A_{10}$	$\displaystyle=$	$\displaystyle-6\bar{G}_{M}(y,R)-(3/x)\bar{G}_{2}(y,R)$	(76)
$\displaystyle A_{11}$	$\displaystyle=$	$\displaystyle 6\bar{G}(y,R)f(x)$	(77)
$\displaystyle A_{12}$	$\displaystyle=$	$\displaystyle 3\gamma(y,R)f^{2}(x)$	(78)
$\displaystyle A_{20}$	$\displaystyle=$	$\displaystyle 6\bar{G}(y,R)\gamma(y,R)$	(79)
$\displaystyle A_{21}$	$\displaystyle=-$	$\displaystyle 3\gamma^{2}(y,R)f(x)$	(80)
$\displaystyle A_{30}$	$\displaystyle=$	$\displaystyle-\gamma^{3}(y,R)$	(81)
$\displaystyle A_{03}$	$\displaystyle=$	$\displaystyle-f^{3}(x)$	(82)

References

(1) W. Horsthemke and R. Lefever, Noise-Induced Transi- tions, Springer (1984).
(2) C.W. Gardiner, Handbook of stochastic methods, Springer (1985).
(3) L. Cai, N. Friedman and X.S. Xie, Nature 440, 358 (2006).
(4) M. Bengtsson, A. Stahlberg, P. Rorsman and Mikael Kubista, Genome Res. 5, 1388 (2005).
(5) N. Rosenfeld, J.W. Young, U. Alon, P.S. Swain, M.B. Elowitz, Science 307, 1962 (2005).
(6) A. Ochab-Marcinek and M. Tabaka, Proc. Natl. Acad. Sci., 107, 22096 (2010).
(7) A. Ochab Marcinek, J. Jedrak, M. Tabaka, Phys. Chem. Chem. Phys. 19, 22580 (2017).
(8) G. Aquino and A. Rocco, Math. Biosci. Eng. 17, 6993 (2020).
(9) L. Arnold, W. Horsthemke, R. Lefever, Z. Physics B 29, 367 (1978).
(10) L.E. Reichl and W.C. Schieve, Instabilities, bifurcation, and Fluctuations in Chemical Systems, University of Texas Press, Austin, 1982.
(11) T.S. Gardner, C.R. Cantor, J.J. Collins, Nature 403, 339 (2000).
(12) N. Rosenfeld, M.B. Elowitz, U. Alon, J. Mol. Biol. 323 785 (2002).
(13) V. Shahrezaei, J.F. Ollivier, P.S. Swain, Mol. Sys. Bio. 4, 196 (2008).
(14) H. Salman et al., Phys. Rev. Lett. 108, 238105 (2012)
(15) N. Brenner et al., Phys. Rev. E 92, 042713 (2015).
(16) M. Thattai and A. van Oudenaarden, PNAS 98, 8614 (2001).
(17) G. Hornung and N. Barkai, PLOS Comp. Biol. 4, 055 (2008).
(18) J. Paulsson, Nature 427, 415 (2004).