Bursting noise in gene expression dynamics: linking microscopic and mesoscopic models

The dynamics of short-lived mRNA results in bursts of protein production in gene regulatory networks. We investigate the propagation of bursting noise between different levels of mathematical modelling and demonstrate that conventional approaches based on diffusion approximations can fail to capture bursting noise. An alternative coarse-grained model, the so-called piecewise deterministic Markov process (PDMP), is seen to outperform the diffusion approximation in biologically relevant parameter regimes. We provide a systematic embedding of the PDMP model into the landscape of existing approaches, and we present analytical methods to calculate its stationary distribution and switching frequencies.

Transcription and translation in the process of gene expression occur at the molecular level and in environments of relatively small copy numbers.The discreteness of the molecular dynamics and the inherent randomness with which reactions occur are known as 'intrinsic noise'.It is now widely accepted that intrinsic noise plays an important role in gene regulatory networks [1][2][3] .It promotes epigenetic diversity and enhances the adaptability of a single phenotype in changing environments 4,5 .To investigate the effects of intrinsic noise, mathematical models at different levels have been constructed, ranging from microscopic models 3,[6][7][8][9][10] describing the finer origins of intrinsic noise to mesoscopic models [11][12][13][14][15] .While the former capture the biological processes in more detail, the latter are computationally scalable and constructed to model more complex networks.These models all capture some signatures of intrinsic noise, but the detailed implementation of stochasticity varies from model to model.It is then important to consider how noise propagates between different levels of mathematical modelling.At present coarse-grained models are often proposed ad hoc and not derived from the more detailed lower-scale models.Is this always mathematically appropriate?What statistics of noise should modellers use at the different levels of coarse graining?What are the consequences of the choice of noise statistics, and what are the pitfalls in deriving models on the meso-level from finer models on smaller scales?These are some of the questions we aim to address in this work.
The above difficulties in transitioning between different levels of modelling can nicely be illustrated in the context of biological switches.These are systems with different metastable states and the possibility to 'switch' between those states.Biological organisms with such behaviour include the Lac switch 6 in Escherichia coli and the Enterobacteria phage λ switch 7 .Computational and mathematical models of these range from very detailed descriptions 6,7 over individual-molecule approaches [8][9][10]16 to mesoscopic models 11,12,14,15 .
The difficulties in connecting these different levels of modelling biological switches are amplified by the recent recognition that the mRNA populations are essential to describing the statistics of regulatory processes 16,17 .Biologically, mRNA molecules are a relatively short-lived source compared to the proteins into which they ultimately translate.Protein production from a given mRNA molecule proceeds while it exists, but ceases after the mRNA decays.This leads to a production of protein in bursts-that is, the production is active for a relatively short and random period of mRNA lifetime, and during that time a random number of proteins is generated.This phenomenon is termed translational bursting 1 and it can be observed in single-molecule experiments 18 .While some mesoscopic models account for such bursting 14,15 , the theoretical investigation of these processes is often limited to their stationary distribution and frequently does not include dynamic features such as switching times.
The aim of our work is to investigate the effects bursting noise in gene regulatory networks [9][10][11]13,16,[19][20][21] , and to construct connections between individual-based models and mesoscopic approaches. Specificall, we start from microscopic and individual-molecule-based models of a toggle switch and set out to construct coarse grained, mathematically tractable models without systematically biasing the outcomes.

RESULTS
Different scales of individual-based models of a toggle switch network.We compare four individualbased models and investigate the effect of bursting noise in a toggle switch network.The first model we consider describes both the mRNA and the protein population dynamics 16 .Fig. 1a illustrates the Markovian model of the regulatory network.Genes X and Y are transcribed into mRNA X and mRNA Y, respectively, which in turn are translated to produce proteins X and Y.The transcription of each of the two genes is suppressed by pro- teins of the respective other type via a Hill function 2,3 H(N ) = K [r 0 + r/[1 + (N/K) n ]], where N stands for the number of suppressing proteins.The model parameter K represents a typical population scale of the proteins, and the parameters r and r 0 set the minimal (r 0 K) and maximal transcription rates ((r 0 +r)K).The parameter n > 0 is the so-called Hill coefficient which models the cooperative binding of the repressors 3 .More details of the reaction scheme can be found in the Supplementary Information.Proteins of either type, and the mRNA molecules degrade with constant rates γ 0 and γ respectively.Biologically, mRNA molecules degrade much faster than the proteins do (γ γ 0 ) 2,14,18 .The translation rate of the mRNA is parametrised by γB where the parameter B is the relative frequency of protein production to mRNA degradation.In this parametrisation, the number of proteins one single mRNA molecule produces during its lifetime is a geometrically distribution random variable with mean B (see Supplementary Information).Biologically the parameter B varies depending on the type of product protein 24 .We assume B 10 in this work 2,8 to investigate the effect of translational bursting.Together with the relatively short lifetime of mRNA molecules, this constitutes the origin of 'translational bursting' in the model 14,25 : a relatively large number of protein molecules is synthesized in a relatively short period of time.
For simplicity, the process in Fig. 1a is assumed to be symmetric with respect to X and Y, but the analysis is easily generalised to asymmetric circuits.In Table I we list a set of estimated values of the parameters for the model organism E. coli, along with relevant references.
In the context of this work the model just described constitutes the most detailed model we will investigate and compare against.It serves as a starting point for the derivation of more coarse grained models, and for these purposes we will refer to it as the 'full model' (FM) in the following.
The FM describes both the mRNA and the protein populations, hence it constitutes a relatively highdimensional system which complicates the mathematical analysis.Notably, the only role of mRNA in the FM is to generate proteins, and so mRNA can be left out, so long as the correct statistics of protein production is retained.The timescale separation between the mRNA and protein lifetimes leads to the following reduced model describing only the protein dynamics.In the limit of infinitely-fast mRNA degradation (γ γ 0 ), proteins are generated instantaneously in bursts of geometrically distributed sizes with a mean B, and in between bursting events protein populations decay with rate γ 0 .We will refer to the reduced model as the GB model (geometrically distributed bursts), see Fig. 1b 17,24 .In the GB model, the transcription rates are regulated via the Hill function exactly as before in the FM.
A further reduction of the GB model involves replac-    n Hill coefficient 3.0 Dimensionless 9,13,19,23 a In 13 r = 1.8 and the time unit is defined as the inverse of the protein degradation rate.In our full model we use this value, normalized by to the mean burst size B = 30 molecules (r = 1.8/30 = 0.06).b In 13 r 0 = 0.2.After normalising with respect to the burst size 30, we obtain 1/150.In 23 r 0 = .05r,which is of the same order as 13 .c In 13 K is set to be 200 molecules.In 23 only the deterministic dynamics are provided and r + r 0 = 4.0.To match the protein population scale ≈ 400 in 13,22 , we impose rK = 400, resulting in a typical population scale of the proteins K ∼ 100 molecules, which is of the same order as that of 13 .

TABLE I. Parameter set.
ing the geometrically distributed burst sizes by a constant size B. We will call this the CB model (constant bursts) 8 .While the reduction of the full model to the model with geometrically distributed bursts is well controlled and exact in the limit γ γ 0 , the effects of introducing constant burst sizes are unclear at this stage, and require a detailed analysis (see below).
An even more reduced model is a model with no bursts 8-10 , we will refer to this as the NB model.The reaction scheme is illustrated in Fig. 1c.In this model, only one single protein is synthesized when a transcription event occurs.We assume a B-fold increased transcription rate so that the average number of proteins synthesized per unit time is consistent with the FM, GB, and CB models.
Only the GB model approximates the stationary distribution of the FM.Numerical simulations of each of the models are carried out using standard methods 26,27 .In the following we present statistical properties of the models, leaving typical sample paths to the Supplementary Information.Fig. 2 displays the numerically computed stationary distributions for the FM, GB, CB and NB models.They illustrate that the profiles of protein expressions in different model settings are quite distinct.This is due to the different representation of the underlying intrinsic noise.
While the stationary distributions of the FM and the GB model are in good agreement with each other, substantial discrepancies from the full model are found in the CB and NB models.In the CB model the stationary distribution of protein numbers is very localised compared to the FM and the GB model.In the NB model the probability distribution is even more sharply concentrated.This is because the NB model misses out two pertinent sources of noise.Bursting production in the CB model amplifies the stochasticity of transcription events and leads to a broadening of the protein distribution.Adding randomly distributed burst sizes (GB model) in-troduces further stochasticity, and diversifies of protein numbers even further.
Based on these results, we conclude that the bursting noise introduced by the mRNA populations significantly broadens the stationary distribution.In addition, the GB model approximates the FM model significantly better than the CB and NB models do.We can effectively discard the CB and NB models as faithful representation of the FM, and our subsequent discussion hence focuses mostly on the GB model.
The GB model approximates the mean first switching time of the FM.The toggle switch has two dynamic attractors, one in which protein X is highly expressed and where protein Y has a low concentration, and the other with inverted roles by symmetry.Starting from one attractor the switch can be driven to the other attractor by fluctuations.The timescale of such a transition quantifies the dynamical stability: the longer the timescale, the more stable the system is at the initial position.As we will study next, the way in which the bursting production of protein is implemented significantly affects the timescale of these switching processes.
Starting from initial condition N X (0) = n x,0 and N Y (0) = n y,0 , we define the first switching time as the time it takes a sample path to reach the symmetric boundary N X = N Y .Mathematically, the first switching time is a random variable.The mean first switching time (MFST) is then the average value of the random first switching time.The MFST depends on the initial condition (n x,0 , n y,0 ).Sweeping across the space of possible initial configuration, the MFST of the FM and of the GB model are measured in simulations and presented in Fig. 3.We show the MFST of the CB in the Supplementary Information.As with the stationary distributions, the data in Fig. 3 indicates that the GB model approximates the switching times of the full model to a good accuracy.We remark that the MFST of the CB model is almost as twice as long as that of the GB and FM models, and the switching time in the NB model is longer than 1000 cell cycles (Supplementary Information).
Diffusion approximation of the GB model.The evolution of the protein population in the GB model is described by a master equation (Supplementary Information).Solving master equations mathematically is however difficult and mostly limited to linear dynamics 24,28 .The only realistic way forward for a theoretical analysis is often the so-called diffusion approximation.
In the diffusion approximation, the discrete-molecule process is approximated by a Gaussian process for continuous concentrations-numbers of the different types of molecules normalized by a typical population scale.The Gaussian process satisfies a diffusion equation (the Fokker-Planck equation) 29,30 .Based on these methods, it is often possible to calculate or approximate the stationary behaviour and switching times of model gene networks.For existing studies in the context of toggle switches see [11][12][13] .
Deriving the diffusion approximation of the GB model requires modest modifications to the standard Kramers-Moyal expansion 29,30 .These modifications are necessary to account for the randomness induced by the geometrically distributed burst size.Details of the derivation can be found in the Supplementary Information, we here only report the final outcome.The expansion results in two coupled Itō stochastic differential equations for the concentrations x t = N X (t)/K and y t = N Y (t)/K.These are valid in the limit of large but finite populations 31 and of the form with drift v and diffusion D given by v(w, z) The quantities dW and dW (y) t represent independent Wiener processes.
The diffusion approximation can only be expected to be accurate when molecule numbers are large so that the concentations x t and y t are effectively continuous.In principle, a similar analysis can also be applied to the master equations of the full model.In the FM mRNA numbers are rather small though (typically < 5, see Supplementary Information), so the Gaussian approximation does not capture the statistics of the intrinsic noise faithfully.Similarly further analysis of the CB and NB models can be carried out based on the diffusion approximation.Given that CB and NB models fail to reproduce the behaviour of the FM, these results are relegated to the Supplementary Information.
Results from simulating the Gaussian process of equations (1) are shown in Fig. 4.While the data for the stationary distribution (Fig. 4a) looks similar to that of the full model (Fig. 2a), noticeable discrepancies are manifest in the mean first switching times (compare Fig. 4b and Fig. 3a).In Fig. 4c and d, we show the differences between simulation outcomes of the full model and those of the diffusion approximation of the GB model.Although the GB model itself approximates the full model well (Figs. 2 and 3), we conclude that the diffusion approximation fails to capture the relevant model outcomes.
Constructing a mesoscopic piecewise deterministic Markov process.We have seen that the diffusion approximation of the GB model fails to reproduce the statistics of the full model.This underlines the need to construct coarse-grained models directly from the full model and without the intermediate step of a proteinonly dynamics.We now proceed to introduce such a model.As before we describe protein concentrations by continuous variables, x and y.The mRNA dynamics are captured by introducing three 'states': The 0-state describes phases in which no mRNA is present.In the X-state there is one mRNA of type X and protein X is generated with rate γb.The quantity b = B/K is the mean burst size in the unit of protein concentration.No proteins of type Y are produced in the X-state.Similarly, in the Y-state protein Y is generated with rate γb.Both types of protein are subject to natural degradation with rate γ 0 in any of the three states.This is described by the following deterministic differ- The rates with which the system transits between the states are based on the dynamics of the FM: No transitions occur directly between the X and Y states.The kinetic scheme is illustrated in Fig. 1d.The stochasticity and discreteness of the mRNA populations is reflected in the random transitioning between the 0-, X-and Y-states.Between these Markovian events the protein concentrations evolve deterministically.We will refer to this model as the piecewise deterministic Markov process (PDMP).
Notably, at most one mRNA molecule of either type can be present in PDMP at any time.Although the model can be generalised to allow more than one mRNA molecule, the analysis below shows that the lowest-order approximation is sufficient to capture the relevant fluctuations of the mRNA dynamics.
The PDMP approximation outperforms the diffusion approximation of the GB model.As in the GB model, we work in the limit of infinitely fast degrading mRNA (γ → ∞).Simulations of the PDMP model in this limit can be carried out using a minor modification of a previously proposed algorithm 15 .We measure the stationary distribution of the PDMP model and the mean first switching times for different initial protein numbers.Results are shown in Fig. 5a and b, and we compare the outcome against that of the full model in Fig. 5c and d.
The simulation data indicate that the PDMP approximation outperforms the diffusion approximation of the GB model, and it provides a more faithful approximation to the FM.This is because the diffusion approximation introduces Gaussian noise.It retains some information about the variance of protein production and degradation, but it does not capture the geometrically distributed burst sizes in the GB model well enough.The PDMP approximation, on the other hand, models exponentially distributed bursts in protein concentration.
The exponential distribution in the PDMP model is the analogue of the geometric distribution in the discrete-molecule GB model.While the PDMP model is an approximation as well, it retains the typical characteristics of the stationary distribution and switching times of the original model.At the same time the PDMP model is suitable for further mathematical analysis (see below).
When does the PDMP outperform the diffusion approximation?We now investigate the robustness of these findings.In Fig. 6 we vary two essential parameters, the mean burst size B, and the population scale K, while keeping the other parameters fixed.We measure the Jensen-Shannon distance 32,33 between the resulting stationary distributions of the PDMP and that of the the full model.Data is shown in Fig. 6a and c.We also compare the mean first switching times starting from one of the stable modes, see Fig. 6b and d.The figure also shows results from the diffusion approximation of the GB model.
Results indicate that PDMP model outperforms the diffusion approximation of the GB model for mean burst sizes of B 5. We conclude that the bursting noise has to be considered in this biologically relevant regime 24 .The PDMP model incorporates only the bursting noise and neglects the demographic noise from random degradation of the proteins.The strength of this demographic noise is proportional to 1/ √ K.The results in Fig. 6 c and d indicate that the difference in describing intrinsic noise propagates to physical observables even when the noise is weak (K ≈ 1000 for fixed B = 30).
Analytic investigation of the PDMP process.The simplicity of the PDMP approach allows us to proceed with a mathematical analysis.We here only outline the main steps, further details are reported in the Supplementary Information.We denote the probability density that the system is in the 0-state and with protein densities x, y at time t by p 0 (x, y, t).Similarly we write p X (x, y, t) and p Y (x, y, t) when the system is in the Xor Y-states.The evolution of these distributions then follows the forward equation where L † d and L † s drive the deterministic flow and the random switching between states respectively.These operators are of the form with The differential operators ∂ x and ∂ y act on all that follows to their right, including the probability densities p 0 , p X and p Y outside the matrix notation in equations ( 5) and ( 6).
The PDMP approximation applies in the limit γ → ∞, i.e., fast return into the 0-state.The resident time in the X-and Y-states is exponentially distributed and scales as γ −1 .It formally tends to zero as γ → ∞.On the other hand the translation rate γB tends to infinity in this limit.Combining the limiting behaviours of resident time and translation rate results in an exponentially distributed increment of protein concentration in each cycle of switching from the 0-state to the X-or Y-state, and then returning to the 0-state.As a consequence the PDMP converges to previously proposed continuousstate bursting models 14,15,18 in the limit γ → ∞, and p 0 (x, y, t) satisfies as detailed in the Supplementary Information.
Analytic investigation of the mean first switching time.One of the strengths of the PDMP formulation (equations ( 5) and ( 6)) is the relative ease with which mean first switching times can be obtained.We first proceed by computing mean escape time from an arbitrary open domain Ω.The mean first switching time can be calculated by setting Ω = {(x, y) : x < y}, recognising that the process can only exit this domain by crossing the boundary x = y.Suppose, the system is initially at (x, y) ∈ Ω, and in state Z ∈ {0, X, Y}.We write T Z (x, y) for the mean first time at which the process exits the domain Ω.The quantities T Z then satisfy the following backward equation 29,30 where L d and L s are adjoint to the operators in equations (6).They are given by with In the infinitely-fast degrading mRNA limit (γ → ∞), and using appropriate boundary conditions (Supplementary Information) we arrive at This is the adjoint equation 29 of the expression in equation (8) on the open domain Ω.Equation ( 12) is solved by a finite difference method, noting that it is self-consistent and no boundary condition needs to be specified.The solution is shown in Fig. 7, and reproduces the simulation outcome of the FM well.We remark that equation ( 12) is only valid for the half-plane Ω.A detailed discussion can be found in the Supplementary Information.
Analytic investigation of the weak-noise limit.The analytical calculation of the stationary distributions of the PDMP model can be pursued further using the socalled Wentzel-Kramers-Brillouin (WKB) method.This technique is based on the ansatz where b = B/K 1.One proceeds by considering x, y) = 0 order-by-order in b.To leading order we find the Hamilton-Jacobi equation where h(z) := H(Kz)/K.This equation is then numerically solved using the algorithm of Heymann and Vanden-Eijnden 34 .Results are shown in Fig. 8.Even though this only provides a first-order approximation and despite the fact that we have used b = 0.15 (which is not very small) we obtain a reasonable agreement with the stationary distribution in Fig. 5a.
For completeness we have also carried out a WKB analysis of the diffusion approximation of the GB, CB and NB models.These are presented in the Supplementary Information.
The leading order function S 0 (x, y) is the so-called 'rate function' which quantifies the rare-event statistics of the process in the weak-noise limit b 1 35,36 .Several studies have suggested that S 0 (x, y) is a suitable candidate for a 'landscape' of the non-equilibrium random processes in models of gene regulatory networks [9][10][11]16,17 . The Hmilton-Jacobi equation ( 14) contains cubic terms such as (∂ x S 0 ) 2 (∂ y S 0 ), while diffusion equations are quadratic in derivatives of S 0 .This illustrates the fundamental difference between the statistics of intrinsic noise in the diffusion approximation and the bursting noise in PDMP.Further more rigorous mathematical investigations into these differences would be very welcome in our view.
We compare the functions S 0 of the PDMP and the diffusion approximation of the GB model in Fig. 8.One observes a much 'shallower' rate function in the PMDP model, especially at larger protein numbers (N X , N Y ≈ 700).This is due to the long tails in the exponential bursting kernel of the PDMP model, which are not present in the diffusion approximation of the GB model.Such a fat-tail bursting kernel enhances the probability for the system to evolve to high protein concentrations.We identify this as the origin of the qualitatively distinct rare-event statistics in the two models.Effects of bursting noise in a multi-switch network.Recently, multi-switch systems have gained interest 13,21,37 .A schematic diagram of the three-way switch network proposed by 13 is shown in Fig. 9a.It is obtained from the classical toggle switch network by including a self-enhancing autoregulation.Our computational and mathematical setup requires only minor modifications to include generalisation to this case.Specifically, we replace the earlier Hill functions by with parameters 13 q 0 = 4, r 1 = −4/5, r 2 = 7/3, n 1 = 3, n 2 = 1, K 1 = 160, and K 2 = 320.The rest of the parameters follows Table I.The negative value of r 1 reflects the positive autoregulation.To evaluate the effects of bursting noise on this multi-switch model, we consider again the full model, the diffusion approximation of the GB model, as well as the CB and NB models of the extended network.Fig. 9 displays the stationary distribution to illustrate the effects of the bursting noise in the multi-switch network.The model without bursts (NB, panel f) has a stationary distribution consisting of three modes, as reported earlier 13 .Inclusion of constant bursts (CB, panel e) diversifies the protein expression and reduces the stability of the mode located at N X = N Y ≈ 230.In the full model (panel b) there is no discernible concentration of probability in the symmetric mode, hence the three-way switching capability appears to be absent.We also notice that the saddle of the distribution in the FM is located at a state with much lower number of proteins compared to the NB and CB models.The most likely switching path 9 from one of the asymmetric modes to the other will differ significantly between the different variants of the model.The diffusion approximation of the GB model (panel c) does not capture the outcome of the FM either.
Overall these findings confirm again that the inclusion of bursty noise statistics has significant effects on the model outcome.Finally, we observe in Fig. 9d that the PDMP model approximates the full model of the threeway switch well.We conclude that randomly distributed burst sizes are again the predominant form of intrinsic noise in the multi-switch network.

DISCUSSION AND CONCLUSION
Explicitly including mRNA dynamics in gene regulatory models inevitably introduces more complexity.We have quantitatively studied the effects of bursting noise 1 in a biologically relevant regime or the model organism E. coli.To our knowledge, this is one of the first which attempts to build a rigorous connection between existing individual-based models 3,8-10 and more coarse-grained models 11,12,14,15 .Results of our simulations indicate that the bursting statistics of transcription and translation are essential ingredients of models of gene regulation.Coarse grained models need to account for bursting to retain correct statistics of noise-driven phenomena such as the switching between different dynamic attractors.
The implications of our observations are relevant to the abstract modelling of regulatory networks in different ways.We are now in a better position to address our opening question, and to say how noise propagates between different levels of modelling.Perhaps more importantly, our study may ultimately help to decide what level of modelling is most appropriate to study gene regulatory circuits computationally.The answer will of course depend on the question in the focus of the investigation.We have examined different levels of coarse graining, and we have identified the steps in these reduction procedures at which significant alternations to different model outcomes are introduced.
Systematically choosing a suitable level of coarsegraining also facilitates the mathematical analysis of regulatory networks.The high dimensionality of full regulatory network effectively makes them intractable.Model reduction is needed to make progress, and our analysis demonstrates that the PDMP formulation is a powerful way forward, and that it can be more suitable than the conventional diffusion approximation.The PDMP model explicitly retains the bursting noise originating from the mRNA dynamics.Even though it effectively disregards the demographic noise from random degradation of the proteins, it delivers accurate predictions for stationary distributions and switching times.
As another strength, the PDMP formulation can relatively easily be generalised to accommodate more complex reactions.For example, in the Enterobacteria phage λ switch it is not the monomer of the synthesized proteins which acts as the repressor to regulate transcription, but instead their dimer.Modelling these processes requires the inclusion of dimerization further downstream after transcription and translation 7,10 .Preliminary results not shown here reveal that the PDMP approximates such dynamics well.
The fact that the piecewise deterministic Markov process is successful in approximating the full model opens a relatively new type of modelling paradigm.We acknowledge that we are not the first to propose this 28,[38][39][40][41] .Our contribution consists in a first analytical treatment of PDMP models and in an systematic embedding into a wider landscape of modelling approaches.The bursting phenomenon is ubiquitous whenever there is a separation of time scales between the source and the product of a biological process.These are mRNA and protein in models of gene regulation, but we expect that these ideas can be applied to other biological problems with similar time-scale separation.

METHODS
Sample paths of the individual-based processes (FM, CB, NB, and GB) are generated using the standard Gillespie algorithm 26,27 .The PDMP process is simulated using the algorithm discussed by Bokes et al 15 .Simulations of the stochastic differential equations resulting from the diffusion approximation are performed with a standard Euler-Maruyama scheme.The geometric minimum action method 34 is implemented using MATLAB R2010a, as is the finite-difference scheme to solve the backward equation (12).Further details can be found in the Supplementary Information.

I. NOTATION
We briefly summarise the notation used in the main manuscript and in this supplement: • Discrete numbers of the two types of protein are denoted by N X and N Y .We write M X and M Y for the number of mRNA molecules of the two types.
• Variables such as x(t) denote continuous particle densities (or concentrations).Specifically x(t) and y(t) are protein densities, i.e., x = N X /K and y = N Y /K in the limit K 1.
• We denote the probabilities in the master equations by capital P (discrete particle numbers).
• The lower-case notation p is used for probability density functions in the diffusion approximation (continuous particle densities/concentrations).

II. MASTER EQUATIONS OF THE INDIVIDUAL-BASED MODELS
The different individual-based model in the main manuscript are uniquely defined by their master equations.Here we briefly summarise the master equations of the FM and of the GB, CB and NB models.

A. Master equation of full model (FM)
We write P a,b,c,d for the probability that the system is in state The master equation of the FM is then The probability of a state is zero if any of the variables a, b, c or d are negative.
B. Infinitely fast-degrading mRNA limit: the GB model In the kinetic scheme of the full model the mRNA decays with a rate γ, but synthesizes a protein with a rate γB.Both of the rates are constants.One an mRNA is created the next event involving this mRNA particle is either the production of protein or the decay of the mRNA molecule.The probability that the next event is the synthesis of a protein is B/(B + 1), and the probability that a decay occurs next (before production of a protein) is 1/(B + 1).The random number, , of protein molecules generated by one particular mRNA molecule during its lifetime then follows a geometric distribution The lifetime of an mRNA molecule is of order O (1/γ).As a consequence, we can think of the protein-generating process as follows in the infinitely-fast decaying mRNA limit (γ → ∞): As soon as an mRNA is transcribed, it immediately releases a random number of proteins drawn from the distribution (2) and then decays.
Protein X Protein Y where is drawn from the above geometric distribution, every time one of the first two reactions fires.The master equation of this process is We have written P c,d (t) for the probability that the system is in state Again, the probability of a state is zero if c or d are negative.

D. Master equation of CB model
The master equation for the CB model is obtained by replacing g( ) → δ ,B , i.e., takes value = B with probability one.We find

E. Master equation of the model without bursts (NB)
In this case we have

III. DERIVING THE DIFFUSION APPROXIMATIONS A. Diffusion approximation of the GB model
Simulations of the full model (FM) show that the number of mRNA molecules present at any one time is typically very small (M X , M Y < 10) when biologically relevant parameters are used.The conventional diffusion approximation relies on large particle numbers, and so it is not adequate for the full model.
Instead, we perform the diffusion approximation to the master equation of the GB model, equation . This is a standard method, and we proceed along the lines of [1].The only complication is the presence of the geometrically distributed random numbers (denoted by ) in the protein-generation reactions.As we will discuss below this requires only modest modifications to the standard Kramers-Moyal expansion.
We assume that the scale of the population size, K, is large but finite, i.e., K 1, and we write x = N X /K and y = N Y /K, and replace P c,d (t) in favour of p(x, y, t).The master equation ( 4) then becomes In the last two terms we have extended the summation over to infinity, terms in which x − /K or y − /K become negative are automatically suppressed as the corresponding probabilities p(x − /K, y, t) and p(x, y − /K, t) vanish.
The above expression can then be written as where • • • l denotes an average with respect to a geometrically distributed random number , i.e., f .
We next expand the above equation in powers of 1/K, keeping only the leading and sub-leading order terms [1].We also use the explicit expressions = B and 2 = B(2B + 1) for the first two moments of the geometric distribution.
We arrive at the Fokker-Planck equation where we have written ∂ x = ∂ ∂x and similarly for ∂ y .Realisations of the random process described by the Fokker-Planck equation ( 9) can be obtained as the solutions of the coupled Itō stochastic differential equations with the defined drift v and diffusion D given by v(w, z) The quantities W (x) t and W (y) t are independent Wiener processes.

B. Diffusion approximation of the CB and NB models
The same procedure can be applied to the master equation of the CB and NB models, and for completeness we report the resulting Fokker-Planck equations.
For the CB model one finds and for the model without bursts (NB) one has

IV. THE PIECEWISE DETERMINISTIC MARKOV PROCESS (PDMP)
A. Construction of the PDMP model In this section we outline the construction of the PDMP approximation, starting from the full model.For this purpose it is useful to introduce the notation Thus the upper indices (a, b) denote the number of mRNA molecules of either type in the system, and the arguments a, b stand for protein numbers.
The master equation (1) can then be written in matrix form We have introduced where the operator L †(m,n) describes the forward evolution of protein numbers when there are M X = m and M Y = n mRNA molecules in the system.Specifically, where E i,j are the shift operators [1] acting on functions of protein numbers.They are defined through Next we consider the limit of fast mRNA decay, that is, large values of γ.More specifically mRNA molecules of either type are generated with rates H(N Y ) and H(N X ) respectively, and we assume that γ is much larger than either of these two rates (γ H(N Y ), γ H(N X ), for any values of N X and N Y ).In this limit, the system is almost always in the state without mRNA molecules (i.e., M X = 0, M Y = 0), except for short spells during which there is either one molecule of mRNA of type X, or one of type Y .The duration of the episodes spent in these (1, 0) and (0, 1) states is of order γ −1 , then a switch back to the (0, 0) state occurs.The probability to find the system in states with M X > 1 or M Y > 1 is even smaller, specifically of order (H/γ) 2 , and we neglect contributions from these states.Equation ( 15) can then be simplified into a forward equation of a three-state model: with In the next step we consider the limit of large values of K. Formally we take the limit K → ∞.
The system can then be described by the protein concentration x = N X /K and y = N Y /K.The corresponding probability distributions in the continuum limit are p 0 (x, y) := P (0,0) (Kx, Ky) K 2 , (21a) p Y (x, y) := P (0,1) (Kx, Ky) K 2 .(21c) On the left-hand side we have introduced the notation 0, X and Y to describe the states in which there are no mRNA molecules (M X = M Y = 0), one mRNA molecule of type and one mRNA molecule of type Y respectively (M X = 0, M Y = 1).This is in-line with the notation in the main manuscript.
The time evolution of the protein concentrations between the random switching events between these three states is then taken to be deterministic.Mathematically this corresponds to expanding the discrete operators L †(a,b) in powers of K −1 , and keeping only the lowest-order advection terms.
This generates so-called Liouville operators, and leads to where L † d and L † s are the forward operators driving the deterministic flow and the random switching between states, respectively.They are given by and Applying the operator ∂ y We note that this equation is not closed in p 0 .
Next, we take the γ → ∞ limit, keeping in mind that H and γ 0 are finite.The system then almost-surely stays in the 0-state, and consequently p X , p Y → 0. Equation ( 26) then reduces to The inverse operator of 1 + b∂ z is and so equation ( 27) turns into the 'forward equation' presented in the main text: C. Equations for the mean first switching time Here we illustrate the detail derivation to the adjoint equation in the main text.We focus on initial conditions y > x and our goal is to calculate the mean time it takes the dynamics to reach states with x = y.We write T Z (x, y) for the time it takes the dynamics to reach a state in which x = y if started from initial condition x, y, and in mRNA state Z ∈ {0, X, Y}.The T Z (x, y) then satisfy the following adjoint equation where L d and L s are the adjoint operators of L † d and L † s .They are given by In the infinitely fast degrading mRNA limit, γ → ∞, equations ( 30) can be seen to converge to The V. WKB ANALYSIS

A. WKB ansatz
In order to find the quasi-stationary distribution of the PDMP model and of the diffusion approximation of the GB and CB models, one uses the ansatz where ∝ K −1 is the magnitude of the intrinsic noise in the protein dynamics.For the purposes of the WKB analysis the noise is assumed to be weak, i.e., 1.

B. DA of the GB model
In the context of the diffusion approximation of the GB model we use = B/K.To leading order (O K 0 ) one finds a Hamilton-Jacobi equation of the form The vector v denotes the deterministic flow and the (scaled) diffusion matrix D is given by with entries C. DA of the CB model As before we use = B/K.A similar leading-order calculation delivers the Hamilton-Jacobi equation, which is again of the form described in equation (36).The only differences are minor modifications in the diffusion matrix, which now has entries

D. DA of the NB model
It is now convenient to use = 1/K.Again one finds a Hamilton-Jacobi equation of the form as above.The diffusion matrix now has entries

E. PDMP model
For the PDMP model, similar calculations deliver the Hamilton-Jacobi equation 0

FIG. 1 .
FIG. 1. Schematic diagrams illustrating the model dynamics.(a) Full model (FM) describing both the mRNA and the protein populations; (b) Protein-only model with geometrically distributed (GB) or constant (CB) bursts.The quantity B is a geometrically distributed random number with mean B in the GB model, and B = B is a constant in the CB model; (c) Protein-only model without bursts (NB); (d) The piecewise deterministic Markov process (PDMP).

FIG. 2 .
FIG. 2. Stationary distribution of protein numbers, shown in the range 0 ≤ NX, NY ≤ 700 on a linear scale on both axes.(a) FM: Full model describing the mRNA and protein populations; (b) GB: protein-only model with geometrically distributed bursts; (c) CB: protein-only model with constant bursts; and (d) NB: protein-only model without bursts.

FIG. 3 .
FIG. 3. Mean first switching time as a function of the initial protein numbers (0 ≤ NX, NY ≤ 700, shown on a linear scale).(a) FM: Full model; (b) GB model.

4 .
FIG. 4. Diffusion approximation of the protein-only model with geometrically distributed random bursts (GB); (a) Stationary distribution as a function of the protein numbers; (b) Mean first switching time (MFST) as a function of the initial protein numbers, in the unit of cell cycles; (c) Net deviation of the stationary distribution from the full model; (d) Net deviation of the MFST from the FM.All axes show the range 0 ≤ NX, NY ≤ 700 on linear scales.
FIG. 5. PDMP approximation.(a) Stationary distribution; (b) Mean first switching time in the unit of cell cycles as a function of initial protein numbers.(c) Net deviation of the stationary distribution from the full model; (d) Net deviation of the MFST of the PDMP model from the FM.All axes are on linear scales and show the range 0 ≤ NX, NY ≤ 700.

FIG. 6 .
FIG. 6. Performance of the PDMP model and the diffusion approximation of the GB model (DA-GB).(a) Jensen-Shanon distance between the stationary distribution PDMP (and DA-GB) and the stationary distribution of the FM; (b) Mean first switching time for varying value of B at fixed K = 200; (c-d) Similar to (a-b) but now varying K at fixed B = 30.

FIG. 7 .
FIG. 7. Theoretical prediction of the PDMP model.(a) Mean first passage time as a function of initial protein numbers, calculated from the backward equation (12); (b) Stationary distribution of protein numbers calculated from the WKB method.Axes of both panels show the range 0 ≤ NX(0), NY(0) ≤ 700 on linear scales.

FIG. 8 .
FIG. 8. Rate functions S0 as functions of the protein numbers 0 ≤ NX, NY ≤ 700 on a linear scale.(a) PDMP model; (b) Diffusion approximation of the GB model.

FIG. 9 .
FIG. 9. (a) Schematic diagram illustrating the network of the three-way switch, remaining panels show stationary distribution of protein numbers in the range 0 ≤ NX , NY ≤ 700 on linear scale.(b) Full model; (c) Diffusion approximation of the GB model; (d) PDMP approximation; (e) CB model; and (f)NB model.

C.
Master equation of the GB model In the limit γ → ∞, the dynamics of the full model can effectively be coarse-grained into a single-species model Gene X H(N Y ) − −−− → × Protein X (transcription and translation of protein X), (3a) Gene Y H(N X ) − −−− → × Protein Y (transcription and translation of protein Y),
boundary conditions for the mean first exist times are determined by T Z (x b , y b ) = 0, for all locations (x b , y b ) ∈ ∂Ω at which the deterministic flow driven by L † d flows out of the domain Ω in state Z. Next, we specify a bounded domain Ω C := {(x, y) : 0 < x < y, y < C}.The boundary conditions of equations (33) are then T X (z, z) = 0 and T Y (z, C) = 0 ∀z < C. We now use these boundary conditions, and integrate the second and the third components of the expression in equation (33).Subsequently we send C → ∞ and arrive at −1 = [−γ 0 x∂ x − γ 0 y∂ y − H(Kx) − H(Ky)] T 0 (x, y)

22 IXFIG. 10 .
FIG. 10.Mean first switching times.All graphs show 0 ≤ N X , N Y ≤ 700 on linear scales.(a) Full model; (b) PDMP; (c) GB model; (d) Diffusion approximation (DA) of GB; (e) CB model; (f) Numerical solution of the adjoint equation of the PDMP.Data are plotted on the same colour scale in all panels to allow comparison.

23 XFIG. 11 .
FIG. 11. Results from the WKB analysis.Panels a, c, e, and g show the WKB rate functions S 0 (N X , N Y ), and panels b, d, f, and h the corresponding approximation for the stationary probability distribution N exp [−S 0 (N X , N Y ) / ] where N is the normalisation factor.All panels show 0 ≤ N X , N Y ≤ 700 on a linear scale.The insets show the stationary distributions viewed from (N X , N Y ) = (700, 700).(a, b) PDMP; (c, d) DA of GB ; (e, f) DA of CB; (g, h) DA of NB.