Limits of noise for autoregulated gene expression

Gene expression is influenced by extrinsic noise (involving a fluctuating environment of cellular processes) and intrinsic noise (referring to fluctuations within a cell under constant environment). We study the standard model of gene expression including an (in-)active gene, mRNA and protein. Gene expression is regulated in the sense that the protein feeds back and either represses (negative feedback) or enhances (positive feedback) its production at the stage of transcription. While it is well-known that negative (positive) feedback reduces (increases) intrinsic noise, we give a precise result on the resulting fluctuations in protein numbers. The technique we use is an extension of the Langevin approximation and is an application of a central limit theorem under stochastic averaging for Markov jump processes (Kang et al. in Ann Appl Probab 24:721–759, 2014). We find that (under our scaling and in equilibrium), negative feedback leads to a reduction in the Fano factor of at most 2, while the noise under positive feedback is potentially unbounded. The fit with simulations is very good and improves on known approximations.


Introduction
It is now widely accepted that gene expression is a stochastic process. The reason is that a single cell is a system with only one or two copies of each gene and of the order tens for mRNA molecules Elowitz et al. 2002;Raj and van Oudenaarden 2008). Experimentally, this stochasticity can even be observed directly by single-cell measurements such as flow cytometry and fluorescence microscopy, which show the inherent fluctuations of protein numbers arising from cell to cell (Li and Xie 2011).
Usually, noise in gene expression is divided into an intrinsic and an extrinsic part Raser and O'Shea 2005). While the intrinsic part leads to variation of protein numbers from cell to cell in the same environment, the extrinsic part is attributed to the different environmental conditions of the cell. In practice, ensemble averages eliminate intrinsic noise, while single-cell measurements over time can be thought of having a constant environment, thus eliminating extrinsic noise (Singh and Soltani 2013;Singh 2014).
Stochasticity in gene expression is not only interesting per se. Today, its role in evolution, development and cell fate decisions is under discussion (Kaern et al. 2005;Maamar et al. 2007;Fraser and Kaern 2009;Eldar and Elowitz 2010;Balázsi et al. 2011;Silva-Rocha and Lorenzo 2010;Wang and Zhang 2011). For instance, noisy gene expression can be detrimental for the survival of cells under harsh conditions (Mitosch et al. 2017;Fraser and Kaern 2009). Still, many cells have to function constantly. Therefore, mechanisms reducing and controlling the level of noise are beneficial for most of real systems.
Under the central dogma of molecular biology, modeling stochasticity of gene expression is straight-forward (see Paulsson 2005 for a review). A gene, which is either turned on or off, is transcribed into mRNA, which is translated into protein.
Both, mRNA and protein are degraded at constant rates. Since the resulting chemical reaction network is linear, the master equation can be solved and all moments can be derived analytically. Most interestingly, the variance can be decomposed into the effects of switching the gene on and off, noise due to the finite life-time of mRNA, and random fluctuations in the production of protein (Paulsson 2005). It is often stated that gene expression tends to occur in bursts, which occur due to the short life-time of the on-state of the gene and due to the short life-time of mRNA (Kumar et al. 2015).
We are interested in the effect of self-regulation on gene expression noise. It is known that a negative feedback loop, i.e. a protein suppressing its own transcription (or translation) leads to a reduced noise, while positive feedback is attributed to increase noise (Lestas et al. 2010;Hornung and Barkai 2008). Although these findings are widespread, a complete mathematical analysis is lacking. At least, for negative feedback, Thattai and Oudenaarden (2001) and in more generality Swain (2004) quantify the effect of negative feedback using a linearization argument. The latter paper further analyzes different feedback models differing between translational and transcriptional autoregulation. Moreover, Dessalles et al. (2017) derive the equilibrium distribution using a multi-scale approach under negative feedback.
Most analyses of noise in unregulated gene expression rely on the master equation (e.g. Paulsson 2005). By the linearity of this equation, a solution can be given explicitly.
Using the approximation that gene switching is so fast that it is effectively constantly transcribed to mRNA, this linearity can as well be used under negative feedback (Thattai and Oudenaarden 2001;Swain 2004). Our approach and also the one performed in Dessalles et al. (2017) differs in two ways. First, we are using martingale methods from stochastic analysis in order to describe the chemical system (Ethier and Kurtz 1986). Second, we can relax the assumption that the gene is transcribed effectively constantly, and therefore derive a more general result. Consequently, we are able to analyze noise in a truly non-linear system under a quasi-steady-state assumption.
The explicit expression we derive for the noise in the number of proteins is also the main difference between our findings and the results obtained in Dessalles et al. (2017). Since the authors of that paper are interested in the case of not very abundant proteins they compute a stationary distribution for the protein using martingale techniques in the context of birth-death processes as opposed to our stochastic diffusion setting.
While the full model of regulated gene expression (or any other chemical reaction network) is usually hard to study, considering an ODE approach instead, which approximates the full model, leads to new insights. Formally, a law of large numbers-usually referred to as a fluid limit-can be obtained connecting the stochastic and deterministic model (Kurtz 1970b;Darling 2002). While such a law of large numbers gives a deterministic limit, fluctuations are studied using central limit results; see Kurtz (1970a). The special situation for gene expression is that the gene and mRNA only have a few copies, while the protein is often in large abundance. Such multi-scale models are often studied under a quasi-steady-state assumption (Seegel and Slemrod 1989). Here, the species in low abundance are assumed to evolve fast, such that the slow, abundant, species only sense their time-average. For such a stochastic averaging, not only a law of large numbers is given e.g. by Ball et al. (2006), but also a central limit result has recently been obtained by Kang et al. (2014).
While a multi-scale approach to stochastic gene expression is not new (see Bokes et al. 2012;Dessalles et al. 2017), the analysis of fluctuations for such systems is not finished yet. In the case of multi-scale diffusion systems, Veretennikov (2001, 2003) derive a limit result for the slow components using a Poisson-equation. The results by Kang et al. (2014) are similar but are based on Markov jump processes instead of a diffusion limit framework. We apply the techniques of Kang et al. (2014) on the chemical reaction network of (un-)regulated gene expression. As our results show, fluctuations take into account all sources of noise and we give explicit formulas for the reduction of noise under negative feedback and the increase in noise under positive feedback.

The model
We are dealing with the standard model of gene expression without and with feedback; see e.g. Dessalles et al. (2017). Using the terminology from Paulsson (2005), we write for the model without feedback (or the neutral model) Here, off and on refer to an inactive and an active gene, respectively. The mRNA is given by R, and the protein by P. While the first line of chemical reactions models gene switching from off to on and back, the second line encodes transcription and degradation of mRNA, and the third line gives translation and degradation of proteins. Exchanging the first line by then models a negative feedback and models a positive feedback. In all cases, we number the equations from left to right and from top to bottom by 1-6, so K = {1, . . . , 6} is the set of chemical reactions. The species counts are given by X i for i ∈ S := {off, on, R, P} for inactive and active gene, mRNA and protein, respectively. In the following we will scale the rates such that the gene switching and the mRNA production happens on a fast time-scale whereas the protein which is also present in higher abundances than the mRNA is evolving on a slower time-scale. This time-scale separation is frequently used in quantitative analyses of gene expression (see for instance Thattai and Oudenaarden 2001;Swain 2004;Ball et al. 2006;Bokes et al. 2012;Dessalles et al. 2017) since it allows to employ a quasi-steady-state assumption for the species evolving on the fast timescale, cf. Kuehn (2015). It basically means that one first solves for the stationary points in the fast sub-system which are then used to describe the dynamics in the slow sub-system. Thus, using some large constant N , we will make use of the following scaling for the abundances of chemical species Reactions are scaled for all models such that we indeed get a time-scale separation, i.e. reaction constants are such that genes and mRNAs evolve much faster than protein numbers. Note, however, that due to X P = O(N ), we need that the protein production rate needs to scale with N . We use the scaled rates κ 2 , ν 2 , κ 3 , ν 3 = O(1), which are given through For the neutral model, we also set (with κ + whereas for negative feedback (with κ + 1 , κ 1 = O(1)) and for positive feedback as the scaled number of genes, mRNA molecules and proteins, respectively, V N i (0) the corresponding initial value and for M denoting the total copy number of genes, we have in the neutral case for independent, rate 1 Poisson processes Y 1 , . . . , Y 6 . (See (Anderson and Kurtz 2015) for a general introduction on theoretical chemical reaction networks.) The first equation changes in the case of negative feedback to  and in the case of positive feedback to In the sequel, we will refer to the model without, with negative and positive feedback simply as For a summary of the parameters and their scalings, see Table 1. We will refer to λ's and μ's as the unscaled parameters and to κ's and ν's as the scaled parameters.

Results
The following results are all stated in terms of the scaled parameters (κ's and ν's and V N P ). For the corresponding formulas using unscaled parameter notation, see Appendix E.

A limiting process for the amount of protein
The following result can be obtained using a quasi-steady-state assumption. It relies on the method of stochastic averaging; see e.g. Ball et al. (2006). Basically, it says that, using the law of large numbers for Poisson processes, i.e. Y (t) ≈ t for large t, and the scaled parameter set we find where E π [·] denotes expectation with respect to the equilibrium dynamics of the fast species V N off , V N on and V N R for a fixed amount of protein, V N P . Note that we use ⇒ for weak convergence of stochastic processes (Ethier and Kurtz 1986 In particular, the equilibrium, i.e. F(v * P ) = 0, is given by Proof We apply results from Ball et al. (2006) and only sketch the proof. In order to derive Eq. (4) we need to replace V on and V R (the fast variables) in the equation for V N P by their equilibria assuming that V N P is constant. Computing these equilibria is done using the corresponding lines in (• neu ), (• neg ) and (• pos ). The resulting distribution π then is the equilibrium on the fast time-scale. For v P fixed they read Plugging this equilibrium into (3) which is the limit for large N of the corresponding equation in (• neu ), we obtain thaṫ with F as in (4). Computation of the equilibria is standard by solving F(v P ) = 0. In particular, we have to solve for the equilibrium of [neg ].

Approximate variance and Fano factor for the amount of protein
Comparing V N P and v P , where v P is the exact solution ofv P = F(v P ) = κ 3 E π [V R ]− ν 3 v P with F from Theorem 1, we assume that V N P ≈ v P + 1 √ N U for some stochastic process U . The random process U will then account for the fluctuations which are not captured by the deterministic approximation above. Hence, and therefore This approach builds on applying a quasi-steady-state assumption whenever possible, i.e. when averaging over the on/off-state of genes in order to derive the deterministic dynamics of v P using F, and the number of mRNA, which is approximated by its mean in order to derive b. Consequently, fluctuations arising from these two mechanisms cannot be accounted for in the resulting variance. As a result, fluctuations read off from (7) will be too small. In contrast, as an application of Kang et al. (2014) (see also Appendix A), we derive the following central limit result, which takes into account all fluctuations in leading order. Precisely, our next goal is to show that √ N (V N P − v P ) converges and to determine the limiting process. This limit will then provide the error due to noise between the deterministic approximation v P and the stochastic process V N P of order √ N . In the proof, we will make use of the method developed by Kang et al. (2014).
Theorem 2 (Central Limit Theorem) Let V N P , v P and F be as in Theorem 1 and assume further weak convergence of the initial conditions: with W the one-dimensional standard Brownian motion and Hence, we see that, in contrast to the Langevin approach (7) above, fluctuations arising from gene switching and RNA dynamics are also accounted for in Theorem 2; see also Sect. 3.3 for an interpretation of the individual terms. For the difference between the Langevin approximation and our result and its implications see also Sect. 4.3.
The proof of the Theorem is given in Appendix B. Briefly, we apply the stochastic averaging principle on multiple time scales developed in Kang et al. (2014). The whole approach is revisited in Appendix A. There we also state the conditions which need to be satisfied for the theory to apply. Amongst others these include solving a certain Poisson equation which enables a clean time-scale separation.
Remark 1 (Deriving the Fano factor in equilibrium) While (8) provides a dynamical result along paths of X P , we can also use this approximation and study the process in equilibrium by setting X P (0) = v * P N , where v * P is the unique solution of F(v P ) = 0 given by (5).
In order to compute the approximate variance of V N P , when started in the equilibrium v * P , we make use of the fact that the stochastic differential equation (SDE) in (8) is solved by an Ornstein-Uhlenbeck process. In particular, we obtain at late times, i.e. when the process reached its equilibrium (recall that X N P = N V N P is the total number of proteins) with F(·) from (4), v * P from (5) and c(·) from (9). Since no factor N appears on the right hand side, some authors call the Fano factor dimensionless. Empirically, it was found e.g. by Bar-Even et al. (2006), that for all classes of genes and under all conditions, the variance in protein numbers was approximately proportional to the mean, which is again reminiscent of the lacking N in the Fano factor above.
We note here that this approach of computing the Fano factor of X N P in equilibrium was achieved by an unjustified exchange of limits. Namely, for the approximate Fano factor of X N P in equilibrium, we would have to perform t → ∞ first, and only then compute N → ∞, but our approach exchanged this limit.
In order to compute the right hand side of (10), note that Plugging in the equilibrium v * P from (5) and therefore In addition, Hence, plugging these quantities into (10) gives (with X * P = N v * P ; see also (42) for the Fano factor in terms of unscaled parameters)

Interpretation of the Fano factor
The expressions above in Eq. (14) are not only a result from strict calculations but can also be interpreted in biological terms. For instance, for negative feedback, we find-as in the neutral case-contributions from randomness in gene switching, translation and transcription. Moreover, the negative feedback pushes the amount of protein faster back to its equilibrium value for a burst of gene expression. This results in the denominator in which has the biggest noise-reducing effect of negative feedback which we will study hereafter. Adjusted explanations hold in the cases of no or positive feedback.

Comparing the noise in [neu], [neg ] and [pos]
It is frequently reported that a negative feedback in gene expression results in a reduced variance (noise) of protein levels, whereas a positive feedback enhances noise (Lestas et al. 2010;Hornung and Barkai 2008). These observations can be made precise by our results from above. Here, we report some consequences on the equilibrium variance and the Fano factor, . The gene association rate κ + 1 is then chosen such that the protein mean equals 1250 or 60 in each case, respectively. Furthermore, these rates are adjusted in the cases of negative and positive feedback according to parameters as above. Then, from (4), we see that all models have v * P as their unique deterministic limit with the same c(v * P ) from (13). Setting V neu , V neg and V pos as the variance of the model without, with negative and with positive feedback, respectively, and plugging in all quantities in (14) then gives (note that the mean cancels out; see (43) for a version with unscaled parameters) In particular, we see that the variance is reduced in [neg ] and increased in [pos], as expected; see also Fig. 1. Moreover, the graphs show that the performed simulations fit our predictions well for both higher (a) and lower (b) values of v * P . Additionally, we see from Eq. (16) that the change in noise is maximal if the gene is off most of the time, while still having the same amount of protein as in the unregulated (neutral) case. This finding is reminiscent of the fact that gene expression comes in bursts. The burstiness is most extreme if the gene is on only for a short time, producing a large amount of mRNA, and afterwards off for a long period. Especially, we see that for κ 1 → ∞ the maximal reduction due to negative feedback is twofold while the increase in noise is unbounded for κ − 1 → ∞ in case of positive auto-regulation of the protein.

Negative feedback in a simpler model
In the literature simpler models of gene expression are studied as well. Here, only two molecular species are involved. Either, the gene is constitutively expressed, therefore ignoring the state of the gene, or translation is neglected and the gene is assumed to be transcribed leading to protein in one step (Hornos et al. 2005;Ramos et al. 2015Ramos et al. , 2011Shahrezaei and Swain 2008). In either case, we have the model similar to and [pos], we take ( * neg ) and ( * pos ) instead of the first line, respectively. We note that this model arises from the full model described in ( * neu ), ( * neg ) and ( * pos ), when letting ν 2 = κ 2 → ∞. Hence, we obtain the following approximation for the simpler model When comparing these equations with (15) we see that all the terms containing explicit mRNA noise disappeared. Thus, we still have contributions from gene activation and protein production and degradation as well as the scaling factor due to the autoregulation.

Refining the Fano factor
Here, we again consider the original model given by equations ( * neu ), ( * neg ) and ( * pos ), but use a different scaling for the mRNA. We assume that not only the protein evolves on the slow time-scale but also the mRNA production and degradation. This is a slightly more complex model since we cannot average the number of mRNA molecules when analyzing the protein fluctuations. This model allows us to compare our results in a more straightforward way with results obtained previously in Thattai and Oudenaarden (2001), Swain (2004); see Sects. 4.2 and 4.3.
In order to account for the new scaling we set μ 2 = O(1) and λ 3 = O(1). As scaled variables, we introduce ν 2 = μ 2 and κ 3 = λ 3 . In this model, we have X R = O(N ), but still X P = O(N ). With these assumptions it is possible to derive a refined formula for the Fano factor for X P in equilibrium. Precisely, we compute in Appendix D-see Eqs. (37), (39) and (40) with Note that this equation approximately gives (14) for ν 2 , κ 3 1. However, adding interpretations as in Sect. 3.3 is not straight-forward since the additional terms in (18) stem from interactions between RNA and protein dynamics. In practice (and in our simulations below), the life-time of proteins is much larger than the life-time of mRNA, such that (18) does not produce a better fit than (14); see also Fig. 2 and compare the dash-dotted and the solid line. Therefore, when not stated otherwise, we will use (14) in the sequel.

Comparison to previous results
Here, we compare our results in the neutral case with those obtained in Paulsson (2005), and in the case of negative feedback with the formulas for the Fano factor derived in Dessalles et al. (2017), Ramos et al. (2015), Swain (2004) and Thattai and Oudenaarden (2001).

The neutral case, Paulsson (2005)
In Paulsson (2005) the neutral model [neu] was studied without assuming any scalings of the parameter λ i , μ i or of the number of mRNA molecules or proteins. Setting as the expected life-times of a change in gene activity, mRNA and protein, respectively, we see that the Fano factor in equilibrium obeys (see equation (4) in Paulsson (2005)) The approximation in the last line corresponds to our scaling, which is exactly such that τ 3 τ 1 , τ 2 (since μ 3 μ 1 , μ 2 ), i.e. the protein is much more stable than mRNA and the state of the gene. Thus, our approximation of the Fano factor in Eq. (14) (or rather its version in unscaled parameters in Eq. (42)) is in line with (Paulsson 2005).

Negative feedback, Thattai and Oudenaarden (2001)
In Thattai and Oudenaarden (2001), a linearization of [neg ] was studied in the case of fast on-and off-switching of the gene. In particular, this will mean that both, κ − 1 , κ + 1 1. In the following, we derive their result within our framework. Therefore, consider as in the proof of Theorem 1 that Then, averaging out the gene state, we obtain the following system (compare with (• neg )) Assuming, like in Thattai and Oudenaarden (2001), that the equilibrium effective mRNA-production is a linear function, i.e. letting such that we basically model a constitutively expressing gene, we can further approximate V N R by Now, the system (19) and (21) is exactly as on p. 3 in Thattai and Oudenaarden (2001) with Plugging these variables into equation (3) of Thattai and Oudenaarden (2001) we obtain in equilibrium with that (note that η is negligible since N is large and φ is small by (8)) Considering our scaling assumption, and further assuming (20), Eq. (14) can be simplified to which equals the expression from (22). The results of this approximation are compared to our result in Fig. 2. Recalling our exact result for the Fano factor from Eq. (14), we note that due to the linearization of the mRNA expression and thus a basically constant mRNA production, the noise emerging from the random gene switches is not adequately represented in the formula obtained in Thattai and Oudenaarden (2001). To be more precise, in contrast to our formula in (14), the effect of mRNA noise due to gene switching (first term in first bracket) is not taken into account at all. Additionally, the negative feedback (last bracket in (22)) does not affect the noise in the same way as it does in the exact formula (14). As can be seen in Fig. 2, when comparing the solid and the dashed lines, these effects lead to an underestimation of the actual noise which is produced by an exact simulation of [neg ].

Negative feedback, Swain (2004)
As explained around (7), the usual Langevin approximation cannot account for all fluctuations when a quasi-steady-state assumption is made. (Precisely, it cannot account for fluctuations in the averaged variables.) In Swain (2004), a Langevin approach is carried out in order to analyze fluctuations in autoregulatory gene expression in the cases of transcriptional and translational feedback. The author considers the mRNA and the protein to evolve on the same time-scale whereas the gene (or DNA) is considered to be on a faster time-scale, see also Sect. 3.6. For transcriptional feedback (which we study here), Swain obtains in his equation (5)-see below for the transformation of his results into our parameters (unscaled, and scaled) with α neg = κ 1 v * P κ 1 v * P +κ +

. This corresponds exactly to (18)-see also Eq. (44) in
Appendix E for the unscaled version-with a missing term in the numerator, namely α 2 neg κ 2 κ 3 κ 1 v * P ( ν 2 +ν 3 ) . This term arises from fluctuations in gene activation, which was averaged out in calculations done in Swain (2004). At least, (23) arises from (18) if we assume that κ 1 v * P , κ + 1 1, i.e. fast gene switching. For a comparison of the two results we refer to the solid and dotted lines in Fig. 2 where we see that due to the missing term in the approximation obtained by Swain, his result slightly underestimates the Fano factor resulting from simulations of [neg ]. However, the difference becomes smaller for a lower number of expected proteins, see Fig. 2b. This can be explained by the difference between Eqs. (23) and (18). The missing term in Swain's derivation can be related to noise emerging from the gene switching. These processes however are given by the overall parameter configuration. Hence the larger difference between the solid and dotted lines in (a) when compared to (b) simply emerges from the corresponding term in (a) being larger than in (b). Thus, it therefore has a stronger effect on the overall fluctuations in protein numbers.

Negative feedback, Dessalles et al. (2017)
Our result for [neg ] has Theorem 2 of Dessalles et al. (2017) as a limit result. They study a similar model (with M = 1), but with a scaling such that λ 3 , μ 3 = O(1) and λ 1 = O(N ), leading to low (i.e. O(1)) abundance of protein in the system. Assuming the same time-scale separation as in our model, i.e. gene switching and mRNA processes evolve on a fast time-scale, the resulting birth-death process for P has a stationary distribution which they compute explicitly. Moreover, setting they obtain in their Corollary 4.1 that, in equilibrium, For large ρ, our model simplifies to v * P ≈ √ ρ, and then (14) gives the same limit.

Negative feedback in a simpler model, Ramos et al. (2015)
In Ramos et al. (2015), the authors derive in their equation (11) the Fano factor for the simplified model, which we introduced in Sect. 3.5. (We note that their model slightly differs since the protein binds to the gene and therefore cannot degrade in this state.) Their results connect the Fano factor to the covariance of the number of active genes and the number of proteins in equilibrium. Since they also give the limiting distribution (in terms of a confluent hypergeometric function), they can evaluate this covariance and also the Fano factor numerically. From their Figure 2, one can see that-if proteins are somewhat abundant-the Fano factor stabilizes around 1/2 for various parameter combinations; a result reminiscent of Dessalles et al. (2017) as discussed above.

Negative feedback for small amounts of protein
Theorems 1 and 2-and all subsequent calculations-only hold under the scaling described in * or in Appendix D. In these scalings, we have that there the protein is abundant. In this case, we see from (14) that the Fano factor is at least 1/2. In Sect. 4.4, we discussed the results by Dessalles et al. (2017), where X P = O(1) is used, but the limiting result for ρ → ∞ implies that proteins are abundant and leads to a Fano factor of at least 1/2. However, the scaling of Dessalles et al. (2017) also allows for smaller values for the Fano factor, which is called the infra-Fano regime in Ramos et al. (2015). Precisely, consider the model from Dessalles et al. (2017), where the scaling is used. It is shown that X P converges towards a birth-death process with death rates δ n := μ 3 n and birth rates β n := λ 3 μ 2 if there are n proteins. For this process, they compute the equilibrium distribution and Z as a normalizing constant. In this case, one can see that for λ 1 λ + 1 , π is concentrated around 1, i.e. there is a single molecule of the protein and hence, the Fano factor becomes arbitrarily small. Since this in particular means a Fano factor below 1/2, Ramos et al. (2015) call this the infra-Fano regime.

Conclusion
Quantifying noise in gene expression is essential for understanding regulatory networks in cells (Thattai and Oudenaarden 2001). Our results capture most of the previously derived results. While negative feedback is known to reduce noise under auto-regulated gene expression, we improve on the quantification of this effect, i.e. our results account for all possible sources of noise due to gene activation, mRNA fluctuations and the protein processes itself. We note that our results require that proteins are abundant. Since the infra-Fano regime described in Sect. 4.6 relies on small amounts of protein, our results do not recover this regime.
In addition, we provide the same quantification of noise also for positive feedback, where noise is increased. In particular,(14) shows that the average time the gene is off determines the reduction of noise in all cases relative to unregulated genes; see also Grönlund et al. (2013). As we saw earlier, for both, negative and positive feedback, noise difference between the non-regulated and the model with feedback is largest if the gene is off most of the time. This can be interpreted by the burstiness of gene expression. It is largest for genes which are off for long times and then turned on for a short time in which mRNA is produced. Interestingly, previous approaches mostly gave approximations for noise for negative feedback if switching the gene on and off is very fast (Thattai and Oudenaarden 2001;Swain 2004) and if the gene is on most of the time (Thattai and Oudenaarden 2001) or in a simplified model (Ramos et al. 2015). Hence, all previous papers could not have seen the effects of gene activation switching on protein noise. As in previous results also obtained in Dessalles et al. (2017), we find that in the limit where the gene is off most of the time, the negative feedback reduces noise at most by a factor of two. Additionally we find that noise can increase unboundedly for positive feedback.
Today, quasi-steady-state assumptions are frequently used when analyzing chemical reaction networks. While the intuition suggests the correct approach when approximating the system by a deterministic path, studying fluctuations is apparently much less obvious. In Kim et al. (2015), some special cases are studied when a straight-forward approximation of the fluctuations works. In our analysis, we use a new approach by Kang et al. (2014) and can also interpret all terms arising in (14), see Sect. 3.3.
Due to taking into account all potential sources of fluctuations the fit of simulations and theory (see e.g. Fig. 2) is excellent and improves on previous studies. There, noise arising from the gene switching its state has been averaged out, and only the recent approach of Kang et al. (2014) reveals the impact of these stochastic processes on the noise in protein numbers.
In their paper, Kang et al. (2014) gave as an example an approximation of noise for Michaelis-Menten kinetics and a model for virus infection. Their method relies mostly on solving a Poisson equation L 2 h = F N − F, where L 2 is the generator of the fast subsystem (gene and RNA in our example), F N and F describe the evolution of the slow system (protein) including all fluctuations and in the limit using the quasisteady-state assumption, respectively. We stress that this approach is not only useful for equilibrium situations, but also for understanding noise if the slow system has not reached equilibrium yet, e.g. after a cell division.
It was argued that complexity of gene regulatory networks leads to a reduction in the level of noise, while certain network motifs always lead to increased levels of noise (Becskei and Serrano 2000;Cardelli et al. 2016). Experimentally, gene expression noise can be used to understand the dynamics of gene regulation (Munsky et al. 2012). Our analysis should provide an approach for distinguishing between different models of gene regulation based on measurements of noise levels.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A Recalling the approach of Kang et al. (2014)
We will consider a general system of chemical reactions with S as the set of chemical species and K the set of reactions; see also Anderson and Kurtz (2015) for further reference. The chemical reactions have the form where ν k = (ν ks ) s∈S and ν k = (ν ks ) s∈S are vectors of chemical species, i.e. elements of N S . For the dynamics, we assume mass action kinetics, i.e. we set x ν ks s for the reaction rate of reaction k ∈ K. With we can then define the dynamics of the Markov process X = (X s ) s∈S through the process R k , which describes the number of occurrences of reaction k up to time t. We have for independent unit rate Poisson processes Y k , k ∈ K and therefore We will tailor the results of Kang et al. (2014) to the special case we need in our gene expression example. This means that we can make use of several simplifications, e.g. on the form of the infinitesimal generator of the full process. For some scaling parameter N , assume that and V N s = X N s /N , the system V N s is a slow (rescaled) sub-system and V N f is a fast sub-system. We assume that the generator L N of X N has the form where L N 2 describes the dynamics of V N f , i.e. L N 2 f = 0 if f only depends on v s . Our goal is to show convergence for some V s and U . Therefore, we proceed as follows.
1. We have that (with the projection π s on the slow species and F N : is a (local) martingale. For the convergence (28), we assume that for some unique process V s , which holds for where ρ N is the equilibrium of V N f for fixed slow species v s . Thus, the convergence and we have shown (28).

Note that
Assume that we (The '≈' is controlled by ε N 2 below. Note that this is a Poisson equation.) With we obtain that We assume for ε N 1 , ε N 2 from (30) that 3. In order to show convergence of M N 1 − M N 2 , we use the Martingale Central Limit Theorem; see Appendix A.1 in Kang et al. (2014). Note that since the quadratic variation of all integrals dt vanishes, we find that (recall the notation for Chemical Reaction Networks from (25)-(26), which we now equip with a superscript N to account for the scaling constant), with z ⊗2 = zz where ζ ks and ζ kf are the stochiometric changes of the kth reaction in the slow and fast subsystem, respectively. Note that in all applications, we will have that R N k either changes slowly, or changes fast and can thus be approximated by a deterministic curve, such that Now, for the equilibrium ρ N of the fast species for given concentration of slow species, v s , if where the right hand side is a deterministic, absolutely continuous, R s×s -valued function with non-negative time-derivative. Hence we know from the martingale central limit theorem (see Appendix A.1 in Kang et al. 2014) that This gives (29).

B Proof of Theorem 2
Proof We will make use of notation from Appendix A and have to show (➊)-(➍) from Appendix A in all cases. Note that the function F from Theorem 1 already satisfies (➊). In all cases, the system (V on , V R , V P ) is a Markov process with a generator of the form (27) with (1) (and different operators L N 2 ). This already implies that for all cases For [neu], From (➋) and (4), we see that we need to solve

Choosing the Ansatz
we obtain the following equation which we need to solve Solving for g 2 and then for g 1 , we obtain Then, √ N ε N 1 N →∞ ⇒ 0 since h N is bounded in N and ε N 2 = 0 by construction. Hence, we have shown (➌). For (➍), if π is the equilibrium of the fast species U, R for given value v P of the slow species as in Theorem 1, we have that For [neg ], all calculations above are the same, but with κ − 1 replaced by κ 1 v P , and for [pos], all calculations are the same with κ + 1 replaced by κ ⊕ 1 v P .

C The two-dimensional Ornstein-Uhlenbeck process
We recall results for the two-dimensional Ornstein-Uhlenbeck process; see e.g. Gardiner (2009). They are needed when refining the Fano factor in Appendix D.
Theorem 3 (Stationary variance of two-dimensional Ornstein-Uhlenbeck process) Proof Using partial integration, it is easy to see that this SDE is solved by In order to compute the right hand side, we note that, for any 2 × 2-matrix A = a b c d , we have that (for the unit matrix I 2 )

If all eigenvalues of
Hence, we can write We then write which is only possible if Combining this with (31) then gives the result.

D Refining the Fano factor
Here, we study the case , such that we have the scaling α off = α on = 0 and α R = α P = 1; compare with (1). In particular, we use that-see (2)-V R = X R /N , V P = X P /N . Here, V on = X on is fast and (V R , V P ) are slow. Note that in this case, we have that We will denote the unique solution of F(V R , V P ) = 0 by (v * R , v * P ).
Theorem 4 (Central Limit Theorem) Let V N R , V R , V N P , V P and F be as above and with W, W independent Brownian motions, and , Remark 2 In equilibrium, we have Again, we have to show (➊)-➍) from Appendix A in all cases for F as above, which already satisfies (➊). We focus on [neu] first. The system (V on , V R , V P ) is a Markov process with a generator of the form (27) with

This implies that
Hence, for (➋), we have to solve For the quadratic variation in (➍), we find that with z ⊗2 = zz Then, we will plug in c R from (35), and obtain Ornstein-Uhlenbeck processes in all cases. Since they are two-dimensional, their equilibrium (normal) distribution can be computed (see Appendix C).

E Results in unscaled parameter notation
In this section we state the formulas from Sect. 3 in terms of the unscaled parameters, i.e. λ's and μ's. The deterministic approximation or law of large numbers for the number of protein is given in the following remark (compare with Theorem 1). Throughout, we have X P ≈ N V P .
Using Theorem 2 we can now write down the Fano factor obtained in equilibrium: Remark 5 (Fano factor in unscaled parameters) From (41), we find that, with X P (0) = X * P , in equilibrium (see (14)) The comparison of the fluctuations in neutral version negative and positive feedback reads for unscaled parameters (see (16)) For completeness, we also give the unscaled version of the refined Fano factor derived in Sect. 3.6 and Appendix D. From (18), with