{'fraction_non_alphanumeric': 0.040826544080833825, 'fraction_numerical': 0.01533840061899777, 'mean_word_length': 4.445353159851301, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gene expression dynamics with stochastic bursts: exact results for a coarse-grained model

We present a theoretical framework to analyze the dynamics of gene expression with stochastic bursts. Beginning with an individual-based model which fully accounts for the messenger RNA (mRNA) and protein populations, we propose a novel expansion of the master equation for the joint process. The resulting coarse-grained model reduces the dimensionality of the system, describing only the protein population while fully accounting for the effects of discrete and fluctuating mRNA population. Closed form expressions for the stationary distribution of the protein population and mean first-passage times of the coarse-grained model are derived and large-scale Monte Carlo simulations show that the analysis accurately describes the individual-based process accounting for mRNA population, in contrast to the failure of commonly proposed diffusion-type models.

Intrinsic noise originating from the discreteness of interacting particles plays an important role in genetic expression: it diversifies the distribution of protein population, promotes transition between different cellular phenotypes on a population level, and in turn enhances organisms' ability to adapt to changing environments without the need of genetic mutation [1]. There are two primary sources of intrinsic noise in the context of gene expression: transcriptional noise from the stochastic transition between active and repressed states of DNA transcription, and translational noise from the relatively fast action of mRNA to produce the proteins [1,2]. Both steps result in bursts of protein production which are experimentally observed [3,4].
This Letter presents a new mathematical framework to analyze bursting noise in gene expression. Starting from an individual-based model including both mRNA and protein populations we construct a novel coarse-grained process describing only the protein population dynamics that fully accounts for the discreteness effects and fluctuations in the mRNA population. When the mRNA degrades at a much shorter time scale, the approximating process converges to currently proposed bursting models [15,16]. In our process-based framework, mean first-passage times in a autoregulated gene circuit with stochastic bursts can be formulated and calculated. We present analytic solutions along with computational verification from large-scale Monte-Carlo simulations. A key conclusion is that the conventional diffusion approximation of the master equation fails to accurately estimate switching times of the individual-based model.
A simple individual-based model of autoregulated gene expression including both the mRNA and protein populations contains four reaction steps [2,15] as summarized in Fig. 1 − →). In the first step N P is the random number of available proteins and in this autoregulated genetic circuit, the population of proteins regulates the transcription rate. The Hill function H 0 (X) := r 0 + r 1 X n /(K n + X n ) with the Hill coefficient n approximates the autoregulated transcription rate when the gene switches between on and off states on a much shorter time scale [2].
We refer to the process in Fig. 1(a) as the individualbased (IB) model. Although the IB model provides a detail description of both the mRNA and protein populations, it is generally difficult to analyze theoretically except for linear cases [17,18]. Single-species models describing only the protein populations are often adopted, especially for more complicated genetic circuits [9,10,13,14]. However, fluctuations in the mRNA population are an important dynamical factor [11,12] and our objective is to construct a coarse-grained model describing only the protein population accounting for contributions from fluctuation in the mRNA population.
Generally mRNA's degrade much faster than proteins. In the model organism Escherichia coli for example, the mean lifetime of the mRNA is about 2 min while protein lifetimes are 45 ∼ 60 min [15]. As a consequence, a large number of proteins is produced in a relatively short period of time-a phenomenon termed translational bursting. In addition, due to small system size (the volume of E. coli are ∼ 10 −18 m 3 ), the onset of the transcription and the lifetime of the synthesized mRNA are observed in a stochastic manner [19].
Motivated by the observation of translational bursting, we propose a novel expansion to approximate the master equation of the IB process in Fig. 1(a). First, we notice that in the IB model, for any given mRNA number m, the protein population N P (t) is a birth-death processes with constant birth rate mγB 0 and constant per capita death rate γ 0 . Therefore, it is convenient to expand the process describing the protein dynamics conditioning on the mRNA population: each "state" of the system is labeled by the mRNA number m. The transition rate from state m to m + 1 mRNA molecules is the autoregulated transcription rate H 0 (N P ), and the transition rate from state m + 1 to m mRNA molecules is the mRNA degradation rate γ. Within each state of the system we perform a Kramers-Moyal expansion of the birth-death process [20,21] with respect to the system size K ≫ 1. In the lowest order approximation only the advection terms describing the mean-field dynamics are retained [22]. Formally letting the protein concentration be x := N P /K ≥ 0 and the number of mRNA molecules be m ∈ {0, 1, 2, . . .}, in each state the protein density evolves according to the deterministic equatioṅ with transition rates between different states where H(x) := H 0 (Kx) is the scaled Hill-function. Next we note that the mean lifetime of the mRNA is O(1/γ) and in the fast-degrading mRNA limit γ ≫ 1, most of the time the system has either m = 0 or m = 1. We therefore neglect states m ≥ 2 and formulate a closed forward equation for p m (x, t), the joint probability density that the system presents m ∈ {0, 1} mRNA molecules and protein density x at time t: where the forward operator [23] is defined to be where we have defined b := B 0 /K to be a dimensionless parameter characterizing the strength of the bursts. We shall refer to (3) as the piecewise deterministic Markov process (PDMP; schematic diagram Fig. 1(b)) and remark that the process in x alone is non-Markovian.
To proceed with our analytic investigation, an infinitely fast-degrading mRNA limit γ → ∞ is taken. Although in such a limit the mean duration when the system stays in m = 1 state is 1/γ → 0, the protein concentration in m = 1 state increases with a rate bγ → ∞, preserving exponentially distributed random burst with an average burst strength b. In this limit the process stays in m = 0 state almost surely (i.e. p 1 (x, t) → 0 ∀t), and the probability distribution p 0 satisfies a closed and second-order differential equation The stationary probability distribution is obtained by direct integration: where N is the normalization factor. Substituting the explicit form of the Hill function we find the analytic expression for the stationary distribution We remark that in the limit γ → ∞ the PDMP model reduces to the bursting model described by a continuous master equation, and that this result confirms [16]. In some parameter regimes the stationary distribution (7) exhibits bi-stability [24] and can be adopted to model a biological switch [11,16]. Our formulation (3) can be used to derive the mean switching time (MST) between two modes of gene expression in a straightforward way [25]. We begin by deriving the mean first-exit time to leave a domain ( If the initial protein concentration is x ∈ (x 1 , x 2 ) and the initial number of mRNA is m, then the mean time to exit the domain T m (x) satisfies the inhomogeneous equation [20] − where the generator L is the adjoint of the forward op-erator in (4), with boundary conditions The physical meaning of the boundary conditions is clear: when the system starts with the state m = 1-a state with fast production of proteins-at upper boundary x 2 , and when the state m = 0-a state with only degrading proteins-at lower boundary x 1 , immediately the flow leaves the domain (x 1 , x 2 ). Taking the limit γ → ∞ we deduce a closed secondorder differential equation for T 0 , where prime denote derivative with respect to x. The boundary conditions for (11) follow from (10): We remark that while formally deriving the backward equations of the bursting models [15,16] considering only the protein concentration is possible, imposing the correct boundary conditions (12) is not trivial.
The solution (derived in the Supplemental Material) is where the auxiliary functions M (x) and V (x) and the constant C are This solution is a generalization of results in [26,27].
When the system exhibits bi-modality, the mean switching times between two modes of protein expression can be obtained by taking appropriate limits of (13). First, we define a critical density x c separating the low and high protein abundance modes, then take x 1 → 0 and x 2 → x c for the low mode, and x 1 → x c and x 2 → ∞ for the high mode. Careful analysis is needed because (11) is singular at x = 0 (and is presented in the Supplemental Material). The analytic expressions for the mean switching times (MSTs) are where the constant C 2 is We now turn to the diffusion approximation (DA) of the IB process. To our knowledge there is no standard way to derive DA models for general bursting kernels.
In the Supplemental Material we present the straightforward Kramers-Moyal expansion [20,21] of the master equation of the IB process in the limit γ → ∞ yielding the Itô stochastic differential equation where X t is the random population density of the proteins, W t is the standard Wiener process and the scaling factor Γ = 2. An alternative and phenomenological construction the diffusion approximation is to insert the mean and variance of the bursting kernel in the individual-based process (see Supplemental Material) which yields (20) with the scaling factor Γ = 1. To avoid leaking probabilities to negative densities, we put a reflective boundary at the origin x = 0. Analytic expressions for the stationary distribution and the mean switching times of the diffusion equation are derived by standard analysis [20] and presented in the Supplemental Material. We performed numerical simulations to measure the stationary distributions and the mean switching times (MST) in all three models to verify the theoretical analysis. For the IB model, exact sample paths are generated by standard continuous time Markov chain simulations [28]. For the PDMP model, kinetic Monte Carlo simulations can be constructed by generating exact random waiting times to the next transition events [29]. In the limit γ → ∞, we adopted a previously proposed algorithm [30]. For the diffusion approximations we construct a standard Euler-Maruyama integrator of (20).
The parameters were chosen to be in a biologically relevant regime [15,16,31]: K = 200, n = 4, B 0 = 40, r 0 = 2, r 1 = 10, γ = 30, while γ 0 := 1 is chosen to normalize the unit of the time by a natural cell cycle. Fig. 3 presents the stationary probability distributions of the IB, PDMP, and DA models. Note that the lowmode is noise induced and does not exist in the mean-field dynamics (top panel of Fig 3). While the PDMP model captures the stationary distribution of the IB model extremely well, directly expanding the IB stochastic bursting model by Kramers-Moyal expansion (DA with Γ = 2) qualitatively captures the stationary distribution, and the phenomenological DA model with Γ = 1 failed to capture the stability of the low mode. Fig. 4 presents the MST between low and high proteinabundance modes in all three models. Again, the PDMP model well estimates the mean switching times of the IB model, and both the DA models fail by a large amount. When the state is initially below x c , both the DA models under-estimate the transition time because the bursting kernels of the DA model have a thinner (Gaussian) tail compared to to the geometric (for the derivation see Supplemental Material) bursting kernel of the IB model. When the initial state is above x c , the DA model with Γ = 1 over-estimate the MST because the the approximation does not capture the high probabilities of low-density bursts, and the DA model with Γ = 2 underestimate the MST because the approximation fail to capture that the bursting kernel is always positive.
The PDMP approximation works well even for models with a strong noise strength. In our example, the low-mode is of an order of 100 protein molecules, and the noise strength (per each burst) is of order 40 protein molecules. In addition, the PDMP approximation performs well even though an infinitely-fast degrading mRNA limit γ → ∞ is taken and consequently almost surely there is no mRNA presented in the system. Meanwhile we observe an average 0.3188 mRNA in the stationary state of the IB model.
The PDMP model can be easily generalized. For example, finite population and lifetime of mRNA can be considered by generalize (3) to include p m with m ∈ {0, 1, 2, . . .}. The transcriptional bursting can be included by considering multiple stages of the gene. Higher dimensional genetic circuit can be investigated by including more states of the system [32]. These generalizations merit future investigations.
We conclude that bursting originating from the discreteness of the fast-living mRNA molecules and the stochastic transcription events is the dominating noise in individual-based autoregulated gene expression model. Diffusion approximations are no longer adequate to analyze the dynamical properties of bursting systems while the novel expansion described here faithfully captures the dynamical properties of the individual-based model in a biologically realistic parameter regime and serves as a new analytic tool to investigate more complex models with bursting noise.

I. INDIVIDUAL-BASED MODEL INCLUDING THE MRNA POPULATION
In this section, we define the individual-based model including the mRNA populations. The kinetic scheme of the autoregulated process is [1,2] φ−→ 1 × mRNA with a rate H 0 (N Protein ) , mRNA−→ mRNA + 1 × Protein with a rate γB 0 , mRNA−→ φ with a rate γ, where the Hill function is defined to be We will denote the (random) population of mRNA and protein by N m and N P respectively.
Continuous time Markov chain simulations [3] are constructed to generate exact sample paths.

II. DERIVING THE PDMP IN THE LARGE POPULATION LIMIT
In the fast degrading mRNA limit (γ/γ 0 ≫ 1), the system only presents only 1 or 0 mRNA in a majority portion of the time. Our proposed approximation is to consider the process (1) conditioning on whether or not the system presents an mRNA, and truncate the probabilities associated with mRNA number greater than 1: P 0,n = − [H 0 (n) + γ 0 n] P 0,n + γP 1,n + γ 0 (n + 1) P 1,n+1 .
Next, for each of the master equations of the protein number n conditioning on the mRNA number m, we perform the conventional Kramers-Moyal expansion [4]. Denote a typical population scale of the protein by N Ω ≫ 1. Note that in the autocorrelated circuit, it is convenient to choose N Ω = K. In the continuum limit, the population density is defined by x := n/K, and the mean "burst" size is defined is b := B 0 /N Ω . The evolution of the probability distributions p 0 (x, t) := P 0,n (t)/K and p 1 (x, t) := P 1,n (t)/K is well-approximated by two coupled Fokker-Planck equations [5] The coupled Fokker-Planck equation can be expressed in a compact matrix form: with We again remind the reader that the differential operators ∂ x and ∂ 2 x act on p 0,1 too. It should be clear that the the discrete population of the proteins causes the demographic stochasticity, which is described by those terms with a prefactor 1/K. We further propose to take the limit K → ∞ [5] and leave only the advection terms in (7) to consider exclusively the bursting noise, which is a result of randomly production and degradation events of the mRNA. In such a limit, the process becomes a piecewise deterministic Markov process: in each state of m = 0 or m = 1, the process is deterministic but the switching between the states is Markovian. We emphasize that, in such a limit, the demographic noise which comes from the discreteness of the protein population does not exist-condition on a m state, the concentration of the protein on its own is always evolving in a deterministic fashion.
We notice that the duration of m = 1 state is of order O (1/γ), but the resident time of the m = 0 state does not depend on γ-it is of of order O γ 0 . As a consequence, p 1 scales O (1/γ) and as γ → ∞. In fact, it can be rigorously proved that p 1 → 0 as γ → ∞, and γp 1 in (7b) can be eliminated by (7a):

III. INDIVIDUAL-BASED BURSTING MODEL
Back to process (1), when γ ≫ γ 0 and H(N Protein ), there is a time-scale separation and the mRNA's degrade at a very rapid rate. As a consequence, when one mRNA is formed, almost surely the next happening events before its final degradation are the even more rapid productions of proteins.
Due to the time-scale separation, the production of other mRNA's is and the protein degradation are negligible in one mRNA's lifetime. In such a limit, the distribution of the total number of proteins an mRNA could ever synthesized before its final degradation can be computed. Define the total number of proteins an mRNA could ever synthesize to be N Σ , a non-negative random variable. Because the mRNA has only have two choices-either to degrade or to produce a protein-at any time before the final degradation, the probability that the mRNA produce a protein is B 0 /(1 + B 0 ) from reading off the ratio of the rates in the process (1). Therefore, the distribution of N Σ is a geometric distribution As a consequence, process (1) in the limit γ → ∞ can be re-formulated to neglect the mRNA population We remind the reader that the parameter B 0 is the mean number of the proteins an mRNA can produce. We shall refer to model (12) as the "individual-based bursting model".
Let P n to be the probability when the system has exactly n proteins. The master equation of process (12) can be deriveḋ

IV. EQUIVALENCE BETWEEN THE PDMP AND CONTINUOUS STATE BURST-ING MODELS
We now apply Kramers-Moyal system-size expansion is performed only to the degradation dynamics in (13). The expansion of (13) yields where p(x, t) := P n (t)/N Ω is the continuum probability distribution, x := n/N Ω is the population density of the protein, and W (x − y) is a kernel of the bursting process, defined by the approximating the discrete by the trapezoid rule: In the infinity population limit K → ∞, (14) reduces to which is exactly the continuous master equation in Friedman et al. [6] It is straightforward to establish the equivalence of the piecewise deterministic Markov process [i.e., (7) in the limit K → ∞] and the continuous-state bursting model (12): acting an operator 1 + b∂ x to (16) yields (10) as K → ∞.

V. DERIVING THE DIFFUSION APPROXIMATION
Often in higher dimensional systems, diffusion processes are adopted to analyze complex genetic circuits [7][8][9][10]. This section present the derivation to the diffusion approximation of the process (12).
In the fast-degrading mRNA limit [i.e. γ → ∞ in (1)], diffusion approximation can be obtained by by performing Kramers-Moyal expansion to (12). The corresponding Fokker-Planck approximation reads N Σ is a geometric distribution, and the exact expression of the first two moments are When the bursting number is large B 0 ≫ 1 (typical biological value 10 1 ∼ 10 2 in E. Coli [2,11]), we arrive at the final diffusion approximation of (12): We finally remark that in the large population limit, b scales O K −1 . A sensible population scaling suggests that near the mean-field fixed points, . It is clear that the diffusion term can then be simplified if we neglect demographic noise due to the protein degradation, when B 0 ≫ 1.
or equivalently the Itô stochastic differential equation where dW t is the Wiener process.
In the main text, we refer the diffusion process (20) to be the diffusion approximation of process (1) in the fast-degrading mRNA limit (γ → ∞), with a reflective boundary at x = 0.
A more phenomenological way is to assert the drift and the diffusion terms of the stochastic differential equation to be the mean (b) and variance (b 2 ) of the exponentially distributed burst size. It can be shown [12] that this approach corresponds to a constant-burst model. In this case, the stochastic differential equation is clearly Finally, we remark that (20) is derived from expanding the master equation (12) where effect that the bursting noise only enhance the population of the proteins, and (21) clearly over-estimate the noise. On the other hand, the phenomenological approach (22) clearly under-estimate the low-density bursts of the exponentially distributed kernel.

VI. MEAN SWITCHING TIME OF THE PIECEWISE DETERMINISTIC MARKOV PROCESS
On the domain Ω := {x : 0 < x 1 ≤ x ≤ x < ∞}, the backward equation reads and in the limit with fast-degrading mRNA γ → ∞, it converges to where we denote γ 0 x by f (x). It is elementary to eliminate the variable T 1 and obtain Define an auxiliary function M (x) and V (x) Now (25) can be expressed as and the formal solution is With two constants of integration C 0 and C 1 . Since T 0 (x 1 ) = 0, clearly C 0 = 0. The second constant C 1 can be determined by the second boundary condition After some algebra, we arrive at For the mean switching times between the high-and the low-concentration mode, we have to impose either x 1 → 0 (for initial state in the low mode) or x 2 → ∞ (for the initial state in the high mode).
A. x1 → 0 limit Note that exp [−M (x)] has a singularity at x = 0. However, V (x) is a well-behave function near x = 0, and we claim The second term is bounded. As for the first term, since y < ǫ ≪ 1, we have