{'fraction_non_alphanumeric': 0.04269748590461285, 'fraction_numerical': 0.009891851762675817, 'mean_word_length': 4.385197190707726, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Small protein number effects in stochastic models of autoregulated bursty gene expression

A stochastic model of autoregulated bursty gene expression by Kumar et al. [Phys. Rev. Lett. 113, 268105 (2014)] has been exactly solved in steady-state conditions under the implicit assumption that protein numbers are sufficiently large such that fluctuations in protein numbers due to reversible protein-promoter binding can be ignored. Here we derive an alternative model that takes into account these fluctuations and hence can be used to study low protein number effects. The exact steady-state protein number distributions is derived as a sum of Gaussian hypergeometric functions. We use the theory to study how promoter switching rates and the type of feedback influence the size of protein noise and noise-induced bistability. Furthermore we show that our model predictions for the protein number distribution are significantly different from those of Kumar et al. when the protein mean is small, gene switching is fast, and protein binding is faster than unbinding.


Introduction
One of the most common gene network motifs is an autoregulatory feedback loop whereby protein expressed from a gene activates or represses its own expression [1]. It has been estimated that 40% of all transcription factors in Escherichia coli self-regulate [2] with most of them participating in negative autoregulation [1]. Feedback leads to a regulation of the magnitude of intrinsic noise [3][4][5][6], and also to changes in the response time and relaxation time of transcription networks [2,7].
The predictions by stochastic models of autoregulation have also been shown to lead to considerably different dynamics than those by deterministic models. For example, while deterministic models of non-cooperative autoregulation predict monostable gene expression for all parameter values, stochastic models of the same set of reactions show switching between two distinct gene expression levels and thus lead to bistability [8,9]. It has also been shown that in the presence of noise, feedback loops can lead to sustained oscillations in protein numbers in regions of the parameter space where deterministic models predict damped oscillations [10][11][12]. These results suggest that stochastic models are necessary to understand the intracellular dynamics of biological systems utilizing a combination of positive and negative feedback loops and in which at least one molecular component is present in low copy numbers, e.g. circadian clocks [13,14].
The stochastic properties of an autoregulatory gene circuit have been explored mostly by stochastic simulations and to a lesser extent by analytical solutions of various discrete, continuous, and hybrid gene expression models (see [15] for a recent review). Discrete models are those in which gene, mRNA, and protein numbers change by discrete integer amounts when reactions occur; in continuous models, fluctuations correspond to hops on the real axis rather than on the integer axis [16][17][18][19]; in hybrid models, some fluctuations are modeled discretely (such as those of genes) while other types of fluctuations (such as those of mRNAs and proteins) are modeled in a continuous sense [20][21][22][23]. Here we focus on discrete models since these are the most realistic among the three types (continuous and hybrid models have the advantage of possessing simpler distribution solutions that can be easier to interpret).
In the literature, there are two exact solutions for the steady-state protein number distribution of stochastic autoregulatory models [9,24]. Both models assume that a single gene can exist in one of two states (with two different protein production rates) and that reversible binding of a protein molecule to a gene leads to switching from one gene state to the other. The differences between the two models are as follows. The model solved by Grima et al. [9] (henceforth referred to as the Grima model; see Fig.  1(a) for an illustration) assumes that (i) protein molecules are produced one at a time and (ii) when a protein molecule binds to a gene, the protein copy number decreases by one whereas it increases by one when unbinding occurs. In contrast, the model solved by Kumar et al. [24] (henceforth referred to as the Kumar model; see Fig. 1(b) for an illustration) assumes that (i) proteins are produced in random bursts whose size is a random number sampled from a geometric distribution and (ii) when a protein molecule binds to a gene or unbinds from it, there is no change in the protein copy number. The advantage of the Kumar model is its modeling of bursty protein expression which is in agreement with experimental data [25]; its disadvantage (unlike the Grima model) is the implicit neglect of fluctuations due to reversible protein-promoter binding which necessarily implies that it cannot precisely capture low protein number effects. Taking into account the latter is important because the copy number of DNA-binding proteins can be in the single digits [26].
In this paper, we remedy the shortcomings of the Grima and Kumar models, by deriving and solving a discrete stochastic model of autoregulatory feedback loops that takes into account both translational bursting and protein number fluctuations during the binding-unbinding process. The paper is organized as follows. In Section 2, starting from a stochastic model of autoregulation describing both mRNA and protein dynamics (henceforth referred to as the full model; see Fig. 1(c) for an illustration), we use multiscale decimation theory to derive a reduced stochastic model of autoregulation for the protein dynamics in the limit of fast mRNA degradation. The chemical master equation (CME) of this reduced model is similar to that of the Kumar model except that the propensities include protein fluctuations during the binding-unbinding process -we henceforth refer to this as the modified Kumar model which is illustrated in Fig. 1(d). In Section 3, we provide an exact analytical solution for the steadystate protein distribution of the modified Kumar model in terms of Gaussian hypergeometric functions, show its simplification in the cases of fast and slow gene switching, discuss the influence of feedback on protein noise and bistability, and verify the theory using stochastic simulations. In Section 4, we show that the modified Kumar model reduces to the Grima model under certain conditions. In Section 5, we compare the modified Kumar model with the Kumar model, showing agreement between the two under slow gene switching and disagreement under fast gene switching and strong protein-promoter interactions. We also investigate how the relative sensitivity of protein noise to model parameters differs between the two models. We conclude in Section 6.

Derivation of the modified Kumar model from a fine-grained model
Based on the central dogma of molecular biology, the gene expression kinetics for an autoregulatory gene network in an individual cell has a standard three-stage representation involving gene switching, transcription, and translation ( Fig. 1(c)) [27]. Due to autoregulation, the promoter could either be free or bound to a protein molecule to form a promoter-protein complex. Let G and G * denote the unbound and bound states of the gene, respectively, let M be the corresponding mRNA, and let P denote the corresponding protein. Then the effective reactions describing the autoregulatory gene circuit are given by: Here σ b is the binding rate of protein to the promoter which characterizes the strength of feedback; σ u is the unbinding rate of protein from the promoter; ρ u and ρ b are the transcription rates when the gene is unbound and bound to protein, respectively; u is the translation rate; v and d are the degradation rates of mRNA and protein, respectively. The reaction scheme describes a positive feedback loop if ρ b > ρ u and describes a negative feedback loop if ρ b < ρ u . We shall refer to this fine-grained model as the full model. (a) Grima model [9]. This model describes protein dynamics only. It neglects translational bursting but takes into account fluctuations in protein number due to binding to and unbinding from the gene. (b) Kumar model [24]. This model describes protein dynamics only. It takes into account translational bursting but neglects protein fluctuations during the binding-unbinding process, i.e. when a protein molecule binds to the promoter, there is no change in the protein number. There are two dotted boxes in this figure; the left one displays the reaction G + P σ b − −→ G * + P and the right one displays the reaction Fine-grained full model. This model has both mRNA and protein descriptions. It (implicitly) takes into account translational bursting (via the mRNA description) and also fluctuations in protein number during the binding-unbinding process. (d) Modified Kumar model. This model, which is derived from the full model in the limit of fast mRNA degradation, takes into account both translational bursting and protein fluctuations during the binding-unbinding process.
The microstate of the gene of interest can be represented by an ordered triple (i, m, n): the gene state i with i = 0, 1 corresponding to the unbound and bound states, respectively, the mRNA copy number m, and the protein copy number n. Let p i,m,n denote the probability of having m copies of mRNA and n copies of protein when the gene is in state i. Then the stochastic gene expression kinetics can be described by the Markov jump process illustrated in Fig. 2(a). The evolution of the Markovian model is governed by the CME 0,m,n = ρ u p 0,m−1,n + (m + 1)vp 0,m+1,n + mup 0,m,n−1 + (n + 1)dp 0,m,n+1 + σ u p 1,m,n−1 − (ρ u + mv + mu + nd + nσ b )p 0,m,n , p 1,m,n = ρ b p 1,m−1,n + (m + 1)vp 1,m+1,n + mup 1,m,n−1 + (n + 1)dp 1,m,n+1 + (n + 1)σ b p 0,m,n+1 Experimentally, it is commonly observed that mRNA decays much faster than protein. For example, mRNA lifetimes in prokaryotes are usually of the order of a few minutes, whereas protein lifetimes are generally on the order of tens of minutes to many hours [28]. Due to timescale separation of the underlying biochemical reaction kinetics, the full model can be reduced to a simpler one. The model reduction technique described below is similar to but slightly different from the one described in [29], where the author considered a different full model which neglects fluctuations in protein number during the binding-unbinding process.
Specifically, let λ = v/d denote the ratio of the degradation rates of mRNA and protein.
Here we make the classical assumption that λ 1 and u/v is strictly positive and bounded [27]. In addition, let q (i,m.n),(i ,m .n ) denote the transition rate of the Markovian model from microstate (i, m, n) to microstate (i , m , n ) and let denote the rate at which the system leaves microstate (i, m, n), which is defined as the sum of transition rates from (i, m, n) to other microstates. Since λ 1, we say that (i, m, n) is a fast state if lim λ→∞ q (i,m,n) = ∞, and we say that (i, m, n) is a slow state if lim λ→∞ q (i,m,n) < ∞.
If (i, m, n) is a fast state, then the time that the system stays in this state will be very short. Note that the limit of λ → ∞ is taken at constant u/v and d. It is easy to check that the leaving rates of all microstates are given by Therefore, all microstates (i, m, n) for m ≥ 1 are fast states and all microstates (i, 0, n) are slow states. By using a classical simplification method of multiscale Markov jump processes called decimation [29][30][31][32][33], the full Markovian model can be simplified to a reduced one by removal of all fast states. For simplicity, microstate (i, 0, n) will be denoted by (i, n) in the reduced model.  The remaining question is to determine the transition diagram and effective transition rates of the reduced Markovian model. This process is described as follows. Suppose that the full model jumps from microstate (i, m, n) to another microstate at a particular time. When λ 1, the transition probability from (i, m, n) to another microstate (i , m , n ) is given by Let (i 1 , m 1 , n 1 ), · · · , (i k , m k , n k ) be a sequence of microstates. We say that c : (i, m, n) → (i 1 , m 1 , n 1 ) → · · · → (i k , m k , n k ) → (i , m , n ) is a fast transition path from (i, m, n) to (i , m , n ) if the intermediate microstates are all fast states. Moreover, the probability weight w c of the fast transition path c is defined as According to decimation theory [29][30][31][32][33], the effective transition rateq (i,n),(i ,n ) from microstate (i, n) to microstate (i , n ) is given byq where c ranges over all fast transition paths from (i, 0, n) to (i , 0, n ). This formula indicates that the effective transition rate from (i, n) to (i , n ) is the superposition of two parts: the direct transition rate and the contribution of indirect transitions via all fast transition paths, as depicted in Fig. 2(c).
Since the intermediate microstates of a fast transition path c are all fast states, in order for the path to have a positive probability weight, all the intermediate transitions along this path must satisfy lim λ→∞ q (i1,m1,n1),(i2,m2,n2) = · · · = lim λ→∞ q (ik,mk,nk),(i ,m ,n ) = ∞. (1) By using this criterion, it is easy to see that the full model has only one type of fast transition paths with positive probability weight, which is given by (see the red arrows in Fig. 2 This is because any fast transition path from (i, 0, n) to (i, 0, n ) with positive probability weight cannot pass through microstate (i, m, k) for some m ≥ 2 and k ≥ 0. Otherwise, there must be an intermediate transition along the path from (i, m − 1, k) to (i, m, k) with transition rate being ρ u or ρ b , which does not diverge as λ → ∞. This contradicts the criterion (1). Moreover, since the intermediate transition rates along the path (2) are all given by u or v, and u, v → ∞ as λ → ∞, it follows that the path (2) satisfies the criterion (1).
To proceed, we define two constants p and q as Since λ 1, the transition probabilities along the path (2) are given by Therefore, the probability weight of this path is ρ u p n −n q if i = 0 and is ρ b p n −n q if i = 1. Since there is no direct transition, the effective transition rate from (i, n) to (i, n ) with n > n is exactly the indirect transition rate via the fast transition path (2) This formula indicates that the reduced model may produce large jumps of protein number within a very short period, which correspond to random translational bursts. Each random burst corresponds to a fast transition path of the full Markovian model. The above computations can be understood intuitively as follows.
Since v d and u/v is finite, the process of protein synthesis followed by mRNA degradation is essentially instantaneous. Once a transcript is synthesized, it can either produce a protein molecule with probability p = u/(u + v) or be degraded with probability q = 1 − p = v/(u + v). Thus, the probability that a transcript produces k copies of protein before it is finally degraded will be p k q, which follows a geometric distribution. Note that since the protein burst size is geometrically distributed, its expected value is given by So far, we have obtained the transition diagram and all effective transition rates of the reduced model, as depicted in Fig. 2(b). Let p i,n denote the probability of having n copies of protein when the gene is in state i. Then the evolution of the reduced model is governed by the coupled set of master equations ρ u p n−k qp 0,k + (n + 1)dp 0,n+1 + σ u p 1,n−1 − (ρ u p + nd + nσ b )p 0,n , p 1,n = n−1 k=0 ρ b p n−k qp 1,k + (n + 1)dp 1,n+1 + (n + 1)σ b p 0,n+1 − (ρ b p + nd + σ u )p 1,n .
This is exactly the CME of the following chemical reaction system (see Fig. 1(d) for an illustration): This reaction scheme takes into account both translational bursting resulting from short-lived mRNA and protein fluctuations during the binding-unbinding process, i.e. when a protein molecule binds to the promoter, the protein number in a single cell is decreased by one and conversely it is increased by one when unbinding occurs. In fact, the reduced model modifies a model of autoregulatory gene networks proposed by Kumar et al. [24]. In the Kumar model, the authors also take into account translational bursting but neglect protein fluctuations during the binding-unbinding process, i.e. when a protein molecule binds to or unbinds from the promoter, the protein number remains the same. In the following, we shall refer to the reduced model as the modified Kumar model. This model was first introduced (using intuitive arguments) and studied using approximation methods in [34]. A method to efficiently infer the parameters of this model from noisy experimental data by means of MCMC has been described in [35].
3 Exact analytic solution of the modified Kumar model

General case
Here we derive the exact analytic solution for the steady-state protein number distribution of the modified Kumar model. For convenience, we normalize all model parameters by the protein degradation rate as At steady-state, all probabilities are time-independent and thus Eq. (3) can be rewritten as To proceed, we define a pair of generating functions Let p n = p 1,n + p 0,n denote the probability of having n copies of protein and let F (z) = f (z) + g(z) denote its generating function. Then Eq. (5) can be converted into the following system of ordinary differential equations (ODEs): Solving these equations leads to (see Appendix A for details) is the Gaussian hypergeometric function, is a normalization constant that can be determined by solving F (1) = 1. Recall that the steady-state distribution of protein number can be obtained by taking the derivatives of the generating function F (z) at zero: Taking the derivatives of the generating function given by Eq. (8) yields We next focus on two trivial special cases. In the case of σ u σ b , ρ u ρ b , d, the gene is mostly in the unbound state and the parameters in Eq. (9) reduce to Since and thus the generating function reduces to Therefore, the protein number distribution is negative binomial and given by [27,37] p n = (r) n n! p n q r , where the gene is mostly in the bound state and the protein number distribution is negative binomial and given by Our analytic results can also be used to derive explicit expressions for several other quantities of interest. For example, the steady-state probability that the gene is in the bound state can be recovered from the generating function f (z) at z = 1, In addition, solving Eqs. (6) and (7) simultaneously for f (z) and g (z), we obtain Substituting z = 1 in this equation yields where B = p/q is the mean protein burst size. Since n = F (1), the steady-state mean of protein number is given by where p G is the steady-state probability of the gene being in the unbound state. The first term on the right-hand side is the mean protein number sB when the gene is in the bound state multiplied by the corresponding probability p G * and similarly the second term is the mean protein number rB when the gene is in the unbound state multiplied by the corresponding probability p G . We stress here that the protein mean of the stochastic model given by Eqs. (13) and (14) is not generally the same as that obtained by solving the corresponding deterministic rate equations in steady-state conditions. Similarly, the steady-state second moment of protein number is given by (see Appendix B for details) This is the sum of three terms: the first term is the second moment of the negative binomial distribution associated with the bound gene state, i.e. the second moment of Eq. (12), multiplied by the corresponding probability; the second term is the second moment of the negative binomial distribution associated with the unbound gene state, i.e. the second moment of Eq. (11), multiplied by the corresponding probability; the last term can hence be interpreted as that arising from the difference between the exact probability distribution and a mixture of two negative binomials (this term becomes negligible in the limit of slow gene switching as we show in Section 3.3).

Regime of fast gene switching
We next focus on two nontrivial special cases. Consider the limiting case when the gene switches rapidly between the unbound and bound states, i.e. σ u , σ b ρ u , ρ b , d. In this case, the modified Kumar model depicted in Fig. 1(d) can be further simplified using another classical simplification method of multiscale Markov jump processes called averaging [33,38]. Since σ u and σ b are large, for any n ≥ 1, the two microstates (0, n) and (1, n − 1) are in rapid equilibrium and thus can be aggregated into a group that is labeled by group n, as shown in Fig. 3(a). In addition, group 0 is composed of the single microstate (0, 0). In this way, the modified Kumar model can be further simplified to the Markovian model illustrated in Fig. 3(b), whose state space is given by Here we emphasize that the group index n cannot be interpreted as the protein number. This is because when n ≥ 1, group n is composed of two microstates with different protein numbers. The remaining question to address is how to calculate the effective transition rates between two groups. In the fast switching limit, the two microstates (0, n) and (1, n − 1) will reach a quasi-steadystate with quasi-steady-state distribution For convenience, we define the effective transcription rate as and the effective protein degradation rate as It is important to here emphasize that fast switching leads to effective transcription and degradation rates whereas it is customary to write reduced master equations in this parameter regime which only have effective transcription rates; this explains the discrepancies between conventional and exact master equation reduction reported in [39]. According to averaging theory [33,38], the effective transition rate from group n to group n + k is given by p qss (0,n)q (0,n),(0,n+k) + p qss (1,n−1)q (1,n−1),(1,n−1+k) = c n p k q and the effective transition rate from group n to group n − 1 is given by So far, we have obtained all effective transition rates for the group dynamics ( Fig. 3(b)). Let p group n denote the probability of being in group n. Then the evolution of the group dynamics is governed by the master equatioṅ To solve this equation, we notice that it is recursive with respect to the group index n. It involves two variables when n = 1, three variables when n = 2, and so on. At steady-state, solving this master equation by induction leads to where K is a normalization constant. Given that there are n copies of protein in an individual cell, the gene can exist in either microstate (0, n) or microstate (1, n). Since microstate (0, n) is contained in group n and microstate (1, n) is contained in group n+1, the steady-state distribution of protein number is given by In fact, this steady-state protein distribution is exactly the same as that obtained from the generating function given by Eq. (8) in the fast switching limit. To see this, we note that when the generating function reduces to To proceed, we recall the following Euler-Pfaff transformation for hypergeometric functions: Using this transformation, the generating function can be simplified as Therefore, the steady-state protein distribution can be recovered taking the derivatives of the generating function at zero: .
Straightforward computations show that This equality, together with the fact that r = c 0 /d, shows that Inserting Eqs. (20) and (21)

Regime of slow gene switching
We next consider the limiting case when the gene switches slowly between the unbound and bound states, i.e. σ u , σ b ρ u , ρ b , d. In this case, the modified Kumar model can also be simplified using the averaging method [33,38]. Since σ u and σ b are small, all unbound microstates (0, n), as well as all bound microstates (1, n), are in rapid equilibrium and thus can be aggregated into a group, as shown in Fig. 3(c). In this way, the modified Kumar model can be further simplified to the Markovian model illustrated in Fig. 3(d), which has only two states.
In the slow switching limit, it follows from Eqs. (11) and (12) that all unbound microstates (0, n) will reach a quasi-steady-state with quasi-steady-state distribution p qss (0,n) = (r) n n! p n q r , and all bound microstates (1, n) will reach a quasi-steady-state with quasi-steady-state distribution p qss (1,n) = (s) n n! p n q s .
According to averaging theory [33,38], the effective transition rate from the unbound group to the bound group is given by So far, we have obtained all effective transition rates for the two-state group dynamics (Fig. 3(d)). Let p G and p G * denote the steady-state probabilities of being in the unbound and bound groups, respectively. Clearly, we have Therefore, the steady-state distribution of protein number is given by which is a mixed negative binomial distribution. Since a negative binomial distribution is unimodal, the mixture of two negative binomials can yield a bimodal steady-state protein distribution. It is not hard to show that this steady-state protein distribution is exactly the same as that obtained from the generating function given by Eq. (8) in the slow switching limit. The details are omitted here. Note that in the slow switching limit, we have σ u p G * = rBσ b p G , implying that the third term in Eq. (15) vanishes leaving only the terms associated with the conditional negative binomials of each gene state.

Gene expression noise
In single-cell experiments, the size of fluctuations around the protein mean is often measured by the squared coefficient of variation η = σ 2 / n 2 , where n is the mean and σ 2 = n 2 − n 2 is the variance. From Eqs. (14) and (15), the steady-state variance of protein number is given by where L = σ u /σ b . Therefore, the size of protein number fluctuations can be decomposed into three terms as Here (1 + B)/ n is the noise due to unregulated bursty gene expression which has been previously derived in the literature [40]. In the presence of autoregulation, two additional terms η + and η − emerge.
In Fig. 4 we show how η + , η − , their sum η + + η − , and the total noise η vary as a function of the gene switching rates σ u and σ b for positive (upper panels) and negative (lower panels) feedback loops. Clearly, η + is always positive. Since 0 ≤ p G * , p G ≤ 1 and p G * + p G = 1, the term η + has the following lower and upper bounds: When L 1 or L 1, we have p G ≈ 0 or p G * ≈ 0. In this case, the lower bound is achieved and η + vanishes. Moreover, the upper bound is achieved when In this case, the ratio L of σ u and σ b is neither too small nor too large. In other words, in order to obtain a large η + , log σ b − log σ u must be controlled within a "narrow" belt ( Fig. 4(a)). We next focus on the term η − . In the slow switching limit, we have Lp G * = rBp G and thus η − vanishes ( Fig. 4(b)). In the fast switching limit, however, the pair of reversible reactions G + P G * are in rapid equilibrium and thus the following approximation is appropriate: This approximation, together with Eq. (14), shows that In addition, we note that since |s − r|B, i.e. the absolute difference between the protein mean in the two genetic states, is usually larger than 1 in living cells, then the signs of 1 − sB + rB and (r − s)B are typically the same. Therefore, in most biologically relevant cases, η − is the product of two terms with different signs and is thus negative in the regime of fast gene switching ( Fig. 4(b)).
In previous works, it has been shown that for an unregulated gene, the protein noise only contains the Poisson noise and mRNA noise [40]. Therefore, the sum of η + and η − characterizes the contribution of an autoregulatory feedback loop to protein fluctuations. In the slow switching limit, η − vanishes and thus the overall feedback contribution is η + . In the fast switching limit, η − is strictly negative and thus the overall feedback contribution η + + η − is less than that in the slow switching limit. From Eq. (23), in the fast switching limit, η − counteracts most of η + and the remaining term is given by which is proportional to s − r. Consequently, the overall feedback contribution is positive (negative) for positive (negative) feedback loops. This clearly shows that in the regime of fast gene switching, positive feedback enhances protein noise and negative feedback diminishes it (Fig. 4(c)). This is consistent with recent theoretical results obtained in [6]. In the fast switching limit, comparing Eqs. (22) and (24), we obtain η + ≈ (s − r)B(η + + η − ).
Since |s − r|B is the absolute difference between the protein mean in the two genetic states, the overall feedback contribution in the fast switching limit is much smaller than that in the slow switching limit when the protein mean is relatively large (Fig. 4(c)). The total noise η is then the superposition of the unregulated contribution (1 + B)/ n and the regulated contribution η + + η − . When L 1, the gene is mostly in the unbound state and thus the protein mean has a negative binomial distribution with mean rB and variance rB(B + 1). In this case, the total noise is given by Similarly, when L 1, the gene is mostly in the bound state and thus the protein mean has a negative binomial distribution with mean sB and variance sB(B + 1). In this case, the total noise is given by This explains why the upper-left and lower-right corners in Fig. 4(d) have different colors. Due to large feedback contribution, the total noise is also large in the regime of slow gene switching, as can be seen from the lower-left corner in Fig. 4(d). In the regime of fast gene switching, feedback regulation gives rise to a weaker enhancement or suppression of the total noise, depending on the sign of the feedback loop.

Numerical simulations
To validate our analytic solution given by Eq. (10), we compare it with the numerical solution obtained using the stochastic simulation algorithm (SSA) for both positive and negative feedback loops in the regimes of slow (Fig. 5(a),(c)) and fast ( Fig. 6(a),(c)) gene switching. Clearly, the analytic solution coincides perfectly with the SSA. In the regime of slow gene switching, both positive and negative feedback loops can lead to bistable gene expression. In the positive (negative) feedback case, the low expression peak becomes higher as the feedback strength σ b decreases (increases) and the high expression peak becomes higher as σ b increases (decreases) (Fig. 5(a),(c)). To gain a deeper insight into bimodal gene expression, we define the strength of bistability as where h low and h high are the heights of the low and high expression peaks, respectively, and h valley is the height of the valley between them. Obviously, κ is a quantity between 0 and 1 for bimodal distributions and is set to be 0 for unimodal distributions. In general, to display strong bistability, the following two conditions are necessary: (i) the two peaks should have similar heights and (ii) there should be a deep valley between the two peaks. The former ensures that the time periods spent in the low and high expression states are comparable while the latter guarantees that the two expression levels are distinguishable. Clearly, κ is large if the two conditions are both satisfied and κ is small if any one of the two conditions is violated. Therefore, κ captures both features of bistability and serves as an effective indicator that characterizes its strength.
Using this definition, we investigate how the strength of bistability κ varies with the gene switching rates σ u and σ b for positive ( Fig. 5(b)) and negative (Fig. 5(d)) feedback loops when the low expression mode is away from zero. It can be seen that both positive and negative feedback loops can only produce bistability in the regime of slow gene switching. When all model parameters are fixed (except a possible interchange between ρ u and ρ b ), a positive feedback loop requires a larger feedback strength to achieve bistability than a negative feedback loop.
The situation is totally different when the low expression mode lies at zero -see Fig. 6. Interestingly, we find that in this case, a positive feedback loop can produce strong bistability under fast gene switching. This is nontrivial because when gene switching is fast, the effective transcription rate c n is a Hill-like function with Hill coefficient being equal to 1. In this case, there is no cooperativity in the protein-promoter binding process and the conventional deterministic theory predicts that bistability can never occur (see Appendix C for a proof of this result). However, our stochastic model predicts that a positive feedback loop is capable of bistable behaviour when gene switching is fast even in the absence of cooperative binding. Note that this is not an artefact of the assumptions used to derive the modified Kumar model since simulations verify that it is also a property of the full model. On the other hand, according to our simulations, a negative feedback loop fails to achieve bistable behaviour in the regime of fast gene switching whether the low expression mode is away from zero ( Fig. 5(d)) or lies at zero (Fig. 6(d)). This is consistent with the deterministic prediction.

Reduction of the modified Kumar model to the Grima model
Consider the limiting case when ρ u , ρ b → ∞ and p → 0, while keeping ρ u p =ρ u and ρ b p =ρ b as constants. Since the mean protein burst size is B = p/(1 − p) and the effective protein production rates in the two gene states are ρ u B and ρ b B, the limits considered above can be interpreted as the limit of negligible mean protein burst size, i.e. B → 0, taken at constant effective protein production rates. In this case since q = 1 − p, we have ρ u pq →ρ u , ρ b pq →ρ b , ρ u p n q → 0 ρ b p n q → 0, n ≥ 2, ṗ 0,n =ρ u p 0,n−1 + (n + 1)dp 0,n+1 + σ u p 1,n−1 − (ρ u + nd + nσ b )p 0,n , p 1,n =ρ b p 1,n−1 + (n + 1)dp 1,n+1 + (n + 1)σ b p 0,n+1 − (ρ b + nd + σ u )p 1,n . This is exactly the CME of the following chemical reaction system (see Fig. 1(a) for an illustration): This reaction scheme coincides with the classical model of autoregulatory non-bursty gene networks proposed by Grima et al. [9]. The Grima model takes into account changes in protein number during the binding-unbinding process but neglects translational bursting, i.e. it assumes that protein molecules are produced one at a time. Here we have derived the Grima model rigorously as the non-bursty limit of the modified Kumar model which itself is the fast mRNA decaying limit of the full model illustrated in Fig. 1(c).
To derive the analytic distribution for the Grima model, we normalize all model parameters by the protein degradation rate asr Recall that when α 1 → ∞ and z → 0, while keeping α 1 z as a constant, the Gaussian hypergeometric function has the following limit where 1 F 1 (α; β; z) is the confluent hypergeometric function. Taking ρ u , ρ b → ∞ and p → 0 in Eq. (8) and applying the above formula, we obtain f (z) =r µK (1 + µ)β es z 1 F 1 (α + 1; β + 1; w(z − z 0 )), is a normalization constant that can be determined by solving F (1) = 1. This is fully consistent with the results obtained in [9]. The steady-state protein distribution for the Grima model can be obtained by taking the derivatives of the generating function F (z) at zero: 5 Comparison of the modified Kumar model with the Kumar model

Regime of slow gene switching
The modified Kumar model solved in the present paper can also be compared with the classical model of an autoregulatory bursty gene circuit proposed by Kumar et al. [24]: The Kumar model takes into account translational bursting but neglects protein fluctuations during the binding-unbinding process, i.e. when a protein molecule binds to the promoter, the protein number remains the same. The dynamics of this reaction scheme can be described by the Markovian model illustrated in Fig. 7(a). It has been shown in [24] that the generating function of the Kumar model is given by where K is a normalization constant and We can see that the generating function of the Kumar model has only one hypergeometric term but, as we proved earlier, the generating function of the modified Kumar model is the superposition of two hypergeometric terms. Comparing Eqs. (9) and (26), we find that the parameters β, w, and z 0 for the two models are exactly the same, while the parameters α 1 and α 2 for the two models are different. We next focus on the limiting case when the gene switches slowly between the unbound and bound states, i.e. σ u , σ b ρ u , ρ b , d. Using the averaging method, the Kumar model can be simplified to the same two-state Markovian dynamics (compare Fig. 7(c),(d) with Fig. 3(c),(d)). Thus, the Kumar model leads to the same steady-state protein distribution as the modified Kumar model in the slow switching limit: In other words, the Kumar and modified Kumar models share the same steady-state behaviour when gene switching is slow. Fig. 8(a)-(c) illustrate the steady-state protein distributions for the two models in different regimes of gene switching. It can be seen that the two models agree with each other very well when gene switching is slow, while they disagree, as expected, when gene switching rates are moderate or large.

Regime of fast gene switching
There is still another case when the Kumar model agrees with the modified Kumar model. In the case of σ u σ b , d, we have b 1 + µ and thus all the five parameters α 1 , α 2 , β, w, z 0 for the two models are the same. Moreover, since β 1, we have A ≈ 0 in Eqs. functions for the two models also coincide with each other. Hence in this case, the two models lead to the same steady-state protein distribution, which can be obtained by taking the derivative of the generating function at z = 0: We next focus on the limiting case when the gene switches rapidly between the unbound and bound states. Since σ u , σ b ρ u , ρ b , d, the two microstates (0, n) and (1, n) of the Kumar model are in rapid equilibrium and thus can be aggregated into a group, as shown in Fig. 7(a). Similarly, using the averaging method, the Kumar model can be simplified to the Markovian model illustrated in Fig. 7(b), where c n is the effective transcription rate defined in Eq. (16). There are two differences between the Kumar and modified Kumar models in the fast switching limit. For the modified Kumar model, both an effective transcription rate c n and an effective protein degradation rate d n < d should be introduced (Fig. 3(b)). For the Kumar model, however, only an effective transcription rate c n should be introduced. In addition, unlike the modified Kumar model, the group index n in the Kumar model can be interpreted as the protein number since each group is composed of two microstates with the same protein number. The steady-state protein distribution for the Kumar model in the fast switching limit is given by [6,18] p n = K (p/d) n n! c 0 (c 1 + d) · · · (c n−1 + (n − 1)d), where K is a normalization constant.
In the fast switching limit, we have seen that the two models lead to the same steady-state protein distribution given by Eq. (27) when σ b σ u . However, the two models may deviate significantly from each other when σ b σ u or when σ b and σ u are comparable. In particular, when σ b σ u , it follows from Eq. (12) that the steady-state protein distribution for the modified Kumar model is the negative binomial distribution p n = (s) n n! p n q s .
Moreover, note that the effective transcription rate c n has the following approximation when σ b σ u : Inserting these equations into Eq. (28), we find that the steady-state protein distribution for the Kumar model is given by a zero-inflated or zero-deflated negative binomial distribution where δ 0 (n) is Kronecker's delta which takes the value of 1 when n = 0 and takes the value of 0 otherwise, and is a constant. We make a crucial observation that γ > 0 for a positive feedback loop and γ < 0 for a negative feedback loop. Therefore, when σ b σ u , the Kumar model leads to a zero-inflated negative binomial protein distribution in the positive feedback case and leads to a zero-deflated negative binomial protein distribution in the negative feedback case.
To validate our theoretical analysis, we compare the steady-state protein distributions obtained using the SSA for the two models in the regime of fast gene switching (Fig. 9). Clearly, the two models agree with each other perfectly when σ b σ u , but they fail as predicted when σ b σ u . In the latter case, there is an apparent zero-inflation or zero-deflation in protein number, depending on the type of feedback loop. In the positive (negative) feedback case, the Kumar model has a much higher (lower) probability of having zero protein copy than the modified Kumar model, and the probability of having nonzero protein copies for the Kumar model is lower (higher). From Eq. (14), we can see that there are two ways of increasing the steady-state protein mean: increasing the transcription rates ρ u and ρ b or increasing the mean burst size B = p/q. In either case, the term q −s in Eq. (29) becomes larger and thus the constant γ becomes closer to zero. When the protein mean is very large, γ is almost zero and thus the zero-inflation/deflation effect becomes almost invisible. In other words, the Kumar and modified Kumar models share similar steady-state behaviour when protein mean is large.

Sensitivity analysis for protein noise
Here we investigate the relative sensitivities of protein noise to various model parameters for the Kumar and modified Kumar models. Recall that the relative sensitivity of protein noise with respect to a parameter s is defined as Λ(s) = ∂ log η/∂ log s, which means that 1% change in s leads to Λ(s)% change in protein noise [41].
Since the two models coincide with each other in the regime of slow gene switching, the relative sensitivities of protein noise for the two models should be the same. This prediction is validated by the numerical simulations in Fig. 8(a) the two models lead to the same relative sensitivities when gene switching is slow, while they fail to agree when gene switching rates are moderate or large.
In the regime of fast gene switching, we also compare the relative sensitivities of protein noise for the two models in positive and negative feedback cases. The relative sensitivities for the two models agree with each other when protein unbinding is much faster than protein binding ( Fig. 9(a),(b)), while they disagree when protein binding is much faster than protein unbinding (Fig. 9(c),(d)). In the positive feedback case, the Kumar model has larger relativity sensitivities than the modified Kumar model. In the negative feedback case, however, the modified Kumar model has larger relativity sensitivities to most of the model parameters compared to the Kumar model. An interesting observation is that the relative sensitivity with respect to a parameter can have a different sign in one model compared to the other model. One such example is the relative sensitivity with respect to B in the negative feedback case with σ b σ u ; in this case, an increased burst size enlarges protein noise for the Kumar model but diminishes protein noise for the modified Kumar model (Fig. 9(d)).

Conclusions
In this paper, starting from a stochastic model of an autoregulatory genetic circuit with both mRNA and protein descriptions (the full model), we use the multiscale decimation method to obtain a reduced model with only the protein description, which we refer to as the modified Kumar model. This model takes into account both translational bursting resulting from short-lived mRNA and protein fluctuations stemming from the binding/unbinding reactions with the promoter. Hence it generalizes two classical models of an autoregulatory feedback loop proposed in previous papers: the Kumar model [24], which takes into account translational bursting but neglects protein fluctuations during the binding-unbinding process, as well as the Grima model [9], which takes into account protein fluctuations during the binding-unbinding process but neglects translational bursting.
We have solved the CME of the modified Kumar model to obtain an exact expression for the steadystate protein number distribution in terms of Gaussian hypergeometric functions. Using the multiscale averaging method, we also obtain the analytic protein distributions in the limits of fast and slow gene switching. In the regime of slow gene switching, the steady-state protein distribution is a mixture of two negative binomials and thus can lead to bistable gene expression for both positive and negative feedback loops. We show that a positive feedback loop requires a larger feedback strength to achieve bistability than a negative feedback loop. Interestingly, we discover that a positive feedback loop can produce noise-induced bistability in the regime of fast gene switching even in the absence of cooperative binding, while a negative feedback loop cannot.
Moreover, we examine steady-state fluctuations in protein abundance in detail. Based on the exact solutions of steady-state protein mean and protein variance, we decompose protein noise, measured by the squared coefficient of variation, into three terms. The first term characterizes fluctuations due to unregulated bursty gene expression which has been reported in previous work [27,37,40]. The other two terms emerge due to the presence of autoregulation, one is positive and one is negative in most biologically relevant situations, which collectively characterize the overall feedback contribution to protein noise. When gene switching is slow, the positive term is large and the negative terms vanishes. In this situation, the overall feedback contribution is large, which results in large protein noise in the regime of slow gene switching. When gene switching is fast, however, the negative term counteracts most of the positive term and thus the overall feedback contribution is much weaker, with its sign depending on the type of feedback loop.
We further study the relationship between the modified Kumar model and the other two classical stochastic models of autoregulation. In the limit of small mean burst size, we have proved that the modified Kumar model reduces to the Grima model. In addition, we show that the modified Kumar model shares the same steady-state behaviour as the Kumar model under slow gene switching. In the regime of fast gene switching, however, the two models agree with each other when the binding rate of the protein to the promoter is much smaller than the unbinding rate, but deviate significantly from each other when the binding rate is much larger than the unbinding rate. In the latter case, the Kumar model yields an apparent zero-inflation (zero-deflation) in protein number for positive (negative) feedback loops. We have also shown that the relative sensitivities of protein noise with respect to all model parameters, as predicted by the modified Kumar and Kumar models, are typically considerably different in magnitude and in some cases also in sign. These differences as well as the artificial zeroinflation/deflation effect for the Kumar model become increasingly significant as the protein mean decreases, and hence our modified Kumar model provides a more accurate description of fluctuations in the low protein number regime.
In this paper we have obtained the steady-state solution of the modified Kumar model. Potential extensions that are currently under investigation include the derivation of approximate results for the time-dependent protein number distribution for both constant and time-varying transcription rates.

B Derivation of the expression for the second moment of protein noise
Here we shall derive the analytic expression for the steady-state second moment of protein number fluctuations given by Eq. (15). Substituting z = 1 in Eq. (6), we obtain µg (1) = bf (1).

C Proof that the deterministic theory only predicts monostability
In the main text, we have shown that a noisy autoregulatory gene network can perform bistability. Here we show that the traditional deterministic theory fails to predict bistability in any parameter region.
According to the chemical reaction scheme of the modified Kumar model (4) and the law of mass action, the deterministic theory of the modified Kumar model is given by the following set of coupled ODEs: ġ = σ b x(1 − g) − σ u g, where g is the mean number of genes in the bound state and x is the mean protein number. In steady-state conditions, the gene number equilibrates to Substituting this equation into the time-evolution equation for the protein numbers, we obtaiṅ is the effective transcription rate defined in Eq. (16). Fig. 10 shows graphs of the functions y = c(x)B and y = dx, whose intersections give the fixed points of the deterministic model. Clearly, there is only one intersection in the negative feedback case, which means that the deterministic model has only one fixed point, which is an attractor (Fig. 10(a)). In the positive feedback case, however, there are one or two intersections, depending on whether ρ u vanishes or not. If ρ u = 0, the deterministic model has only one fixed point that is an attractor (Fig. 10(b)). If ρ u = 0, the deterministic model has two fixed points: one is an attractor which is away from zero and the other is an repeller which lies exactly at zero (Fig. 10(c)). In all cases, there is only one attractor, which implies that the deterministic theory of the modified Kumar model does not allow bistability for any choice of model parameters.