# Adiabatic reduction of models of stochastic gene expression with bursting

Romain Yvinec

<figure><img src="image_1.png" /></figure>

## To cite this version:

Romain Yvinec. Adiabatic reduction of models of stochastic gene expression with bursting. 2013. [hal-02806255](https://hal.inrae.fr/hal-02806255)

HAL Id: hal-02806255

https://hal.inrae.fr/hal-02806255

Submitted on 6 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Adiabatic reduction of models of
stochastic gene expression with bursting

Romain Yvinec*

Université de Lyon CNRS UMR 5208 Université Lyon 1
Institut Camille Jordan 43 blvd. du 11 novembre 1918
F-69622 Villeurbanne Cedex France
e-mail: yvinec@math.univ-lyon1.fr

**Abstract:** This paper considers adiabatic reduction in both discrete and continuous models of stochastic gene expression. In gene expression model, the concept of bursting is a production of several molecules simultaneously and is generally represented as a jump terms of random size. In a general two-dimensional birth and death discrete model, we prove that under specific assumptions and scaling (that are characteristics of the mRNA-protein system) an adiabatic reduction leads to a one-dimensional discrete-state space model with bursting production. The burst term appears then through the reduction of the first variable. In a two-dimensional continuous model, we also prove that an adiabatic reduction can be performed in a stochastic slow/fast system. In this gene expression model, the production of mRNA (the fast variable) is assumed to be bursty and the production of protein (the slow variable) is linear as a function of mRNA. When the dynamics of mRNA is assumed to be faster than the protein dynamics (due to a mRNA degradation rate larger than for the protein) we prove that, with the appropriate scaling, the bursting phenomena can be transmitted to the slow variable. We show that the reduced equation is either a stochastic differential equation with a jump Markov process or a deterministic ordinary differential equation depending on the scaling that is appropriate.

These results are significant because adiabatic reduction techniques seem
to have not been applied to a stochastic differential system containing a
jump Markov process. Last but not least, for our particular system, the
adiabatic reduction allows us to understand what are the necessary condi-
tions for the bursting production-like of protein to occur.

**AMS 2000 subject classifications:** Primary 92C45, 60Fxx; secondary 92C40,60J25,60J75.

**Keywords and phrases:** adiabatic reduction, piecewise deterministic Markov process, stochastic bursting gene expression, quasi-steady state assumption, scaling limit.

Introduction

The adiabatic reduction techniques give results that allow to reduce the dimen-
sion of a system and justify the use of an effective set of reduced equations
in lieu of dealing with a full, higher dimensional model, if different time scales
occur in the system. Adiabatic reduction results for deterministic systems of
ordinary differential equations have been available since the very precise results

<footnote>*corresponding author</footnote>
of Tikhonov (1952) and Fenichel (1979). The simplest results, in the hyperbolic case, give an effective construction of an uniformly asymptotically stable slow manifold (and hence a reduced equation) and prove the existence of an invariant manifold near the slow manifold, with (theoretically) any order of approximation of this invariant manifold. Such precise and geometric results have been generalized to random systems of stochastic differential equation with Gaussian white noise (Berglund and Gentz (2006), see also Gardiner (1985) for previous work on the Fokker-Planck equation). However, to the best of our knowledge, analogous results for stochastic differential equations with jumps have not been obtained.

The present paper gives a theoretical justification of an adiabatic reduction of a particular piecewise deterministic Markov process (Davis, 1984). The results we obtain do not give a bound on the error of the reduced system, but they do allow us to justify the use of a reduced system in the case of a piecewise deterministic Markov process. In fact, we prove limit theorems using martingale strategy, in a similar manner than in recent papers such as Crudu et al. (2012), Kang and Kurtz and Riedler, Thieullen and Wainrib (2012), where general convergence results for discrete models of stochastic reaction networks are given. In particular, these papers give alternative scaling of the traditional ordinary differential equation and the diffusion approximation depending on the different scaling chosen (see Ball et al. (2006) for some examples in a reaction network model). After the scaling, the limiting models can be deterministic (ordinary differential equation), stochastic (jump Markov process), or hybrid (piecewise deterministic process). For illustrative and motivating examples given by a simulation algorithm, see Haseltine and Rawlings (2002); Rao and Arkin (2003); Goutsias (2005). However, we emphasize that we do not consider here a continuous approximation of a discrete model. Rather, we perform adiabatic reduction on both discrete state-space and continuous state-space models. Time-scale reduction have been considered in Kang and Kurtz, but not on the kind we perform here.

Our particular model is meant to describe stochastic gene expression with explicit bursting (Friedman, Cai and Xie, 2006). In discrete state-space bursting models, the variables evolve under the action of a discrete birth and death process, interrupted by discrete positive jumps of random sizes. In continuous state-space bursting models, the variables evolve under the action of a continuous deterministic dynamical system, interrupted by positive jumps of random sizes. In both cases, the positive jumps model the burst production of several molecules instantaneously. In that sense, the convergence theorems we obtain in this paper can be seen as an example in which there is a reaction with size between 0 and $\infty$. We hope that the results here are generalizable to give insight into adiabatic reduction methods in more general stochastic hybrid systems (Hespanha, 2006; Bujorianu and Lygeros, 2004). We note also that more geometrical approaches have been proposed to reduce the dimension of such systems in Bujorianu and Katoen (2008).

Biologically, the bursting of mRNA or protein molecules is defined as the production of several molecules within a very short time, indistinguishable within
the time scale of the measurement. In the biological context of models of stochastic gene expression, explicit models of bursting mRNA and/or protein production have been analyzed recently, either using a discrete (Shahrezaei and Swain, 2008; Lei, 2009) or a continuous formalism (Friedman, Cai and Xie, 2006; Mackey, Tyran-Kamińska and Yvinec, 2011) as more and more experimental evidence from single-molecule visualization techniques has revealed the ubiquitous nature of this phenomenon (Ozbudak et al., 2002; Golding et al., 2005; Raj et al., 2006; Elf, Li and Xie, 2007; Xie et al., 2008; Raj and van Oudenaarden, 2009; Suter et al., 2011). Traditional models of gene expression are composed of *at least* two variables (mRNA and protein, and sometimes the DNA state). The use of a reduced one-dimensional model (that has the advantage that it can be solved analytically) has been justified so far by an argument concerning the stationary distribution in Shahrezaei and Swain (2008). However, it is clear that two different models may have the same stationary distribution but very different behavior (continuous or discontinuous trajectories, monostable or bistable, etc; for an example in that context, see Mackey, Tyran-Kamińska and Yvinec (2011)). Hence, our results are of importance to rigorously prove the validity of using a reduced model. Our results are based on the standard assumption that the mRNA molecules have a shorter lifetime than the protein molecules, that is widely observed in both prokaryotes and eukaryotes (Schwanhäusser et al. (2011)). Depending on the assumed scaling of other kinetic parameters within the mRNA degradation rates, different limiting models are obtained.

The paper is organized as follows. In the first section, we prove a reduction results for a discrete state-space model, that is a two-dimensional birth and death process. Assumptions on the birth and death rates are in agreement with a standard model of gene expression for the mRNA-protein system. That is both variables remain positive and birth of the second variable can occur only if the first variable is positive. Using an appropriate scaling of birth and death rates, we prove that this model converges to a general one-dimensional discrete bursting model.

In the second section, we prove a reduction for a continuous state-space model, that is a two-dimensional piecewise deterministic model of gene expression with a jump production term for the first variable. Using appropriate scaling on parameters, we prove that this model converge either to a deterministic ordinary differential equation or to a one-dimensional continuous bursting model.

# 1. A bursting model from a two-dimensional discrete model

The fact that bursting models arise as a reduction procedure of a higher dimensional model was already observed in Shahrezaei and Swain (2008) and Crudu et al. (2012). In Shahrezaei and Swain (2008), the authors show that, within an appropriate scaling, the stationary distribution of a 2-dimensional discrete model converge to the stationary distribution of a 1-dimensional bursting model. The authors used analytic methods through the transport equation on the generating function. Their result seems to be restricted to first-order kinetics. The
first variable is a fast variable that induces infrequent kicks to the second one. In Crudu et al. (2012), the authors show that, within an appropriate scaling, a fairly general discrete state space model with a binary variable converge to a bursting model with continuous state space. The authors obtained a conver- gence in law of the solution through Martingale techniques. The binary variable is a fast variable that, when switching in an "ON" state, induces kicks to the other variable.

We present below analogous result of Crudu et al. (2012) when the fast vari- able is similar to the one of Shahrezaei and Swain (2008). Our limiting model is still a discrete state space model. These results are more precise than the one of Shahrezaei and Swain (2008), and more general (some kinetics rates can be non- linear). We use martingales techniques, with a proof that is similar to Crudu et al. (2012) and also inspired by results from Kang and Kurtz. We present be- low the model, then state our result in the subsection 1.1, and divide the proof in the three next subsections 1.2-1.4.

We consider the following two-dimensional stochastic kinetic chemical reac- tion model

$$
\begin{equation}
\begin{aligned}
\emptyset &\xrightarrow{\lambda_1(X_1, X_2)} X_1, & \text{Production of } X_1 \text{ at rate } \lambda_1(X_1, X_2) \\
X_1 &\xrightarrow{\gamma_1(X_1, X_2)} \emptyset, & \text{Destruction of } X_1 \text{ at rate } \gamma_1(X_1, X_2) \\
\emptyset &\xrightarrow{\lambda_2(X_1, X_2)} X_2, & \text{Production of } X_2 \text{ at rate } \lambda_2(X_1, X_2) \\
X_2 &\xrightarrow{\gamma_2(X_1, X_2)} \emptyset, & \text{Destruction of } X_2 \text{ at rate } \gamma_2(X_1, X_2)
\end{aligned}
\tag{1}
\end{equation}
$$

with $\gamma_1(0, X_2) = \gamma_2(X_1, 0) = 0$ to ensure positivity. This model can be rep- resented by a continuous time Markov chain in $\mathbb{N}^2$, and is then a general birth and death process in $\mathbb{N}^2$. It can be described by the following set of stochastic differential equations

$$
\begin{align*}
X_1(t) &= X_1(0) + Y_1 \left( \int_0^t \lambda_1(X_1(s), X_2(s))ds \right) - Y_2 \left( \int_0^t \gamma_1(X_1(s), X_2(s))ds \right), \\
X_2(t) &= X_2(0) + Y_3 \left( \int_0^t \lambda_2(X_1(s), X_2(s))ds \right) - Y_4 \left( \int_0^t \gamma_2(X_1(s), X_2(s))ds \right),
\end{align*}
$$

where $Y_i$, for $i = 1...4$ are independent standard poisson processes. The genera- tor of this process is given by

$$
\begin{align}
\mathbb{B}f(X_1, X_2) ={}& \lambda_1(X_1, X_2)[f(X_1+1, X_2) - f(X_1, X_2)] \nonumber \\
& + \gamma_1(X_1, X_2)[f(X_1-1, X_2) - f(X_1, X_2)] \tag{2} \\
& + \lambda_2(X_1, X_2)[f(X_1, X_2+1) - f(X_1, X_2)] \nonumber \\
& + \gamma_2(X_1, X_2)[f(X_1, X_2-1) - f(X_1, X_2)], \nonumber
\end{align}
$$

for every bounded function $f$ on $\mathbb{N}^2$.
Exemple 1. We have in mind the standard mRNA-Protein system given by the following choice: $\gamma_i(X_1, X_2) = g_i X_i$ with $g_i > 0$ for $i = 1, 2$, $\lambda_1(X_1, X_2) = \lambda_1(X_2)$ and $\lambda_2(X_1, X_2) = k_2 X_1$ with $k_2 > 0$. Note however that even in the context of models of gene expression, different models have been proposed, that includes nonlinear feedback of mRNA and/or nonlinear degradation terms Bose and Ghosh (2012).

1.1. Statement of the result

We suppose the following scaling holds

$$ \gamma_1^N(X_1, X_2) = N\gamma_1(X_1, X_2) $$

$$ \lambda_2^N(X_1, X_2) = N\lambda_2(X_1, X_2) $$

where $N \to \infty$ that is degradation of $X_1$ and production of $X_2$ occurs at a faster time scale than the two other reactions. Then $X_1$ is degraded very fast, and induces also as a very fast production of $X_2$. The rescaled model is given by

$$ X_1^N(t) = X_1^N(0) + Y_1 \left( \int_0^t \lambda_1(X_1^N(s), X_2^N(s)) ds \right) - Y_2 \left( \int_0^t N\gamma_1(X_1^N(s), X_2^N(s)) ds \right), \\ X_2^N(t) = X_2^N(0) + Y_3 \left( \int_0^t N\lambda_2(X_1^N(s), X_2^N(s)) ds \right) - Y_4 \left( \int_0^t \gamma_2(X_1^N(s), X_2^N(s)) ds \right), \tag{3} $$

and the generator of this process is given by

$$ \mathbb{B}_N f(X_1, X_2) = \lambda_1(X_1, X_2) [f(X_1 + 1, X_2) - f(X_1, X_2)] \\ + N\gamma_1(X_1, X_2) [f(X_1 - 1, X_2) - f(X_1, X_2)] \\ + N\lambda_2(X_1, X_2) [f(X_1, X_2 + 1) - f(X_1, X_2)] \\ + \gamma_2(X_1, X_2) [f(X_1, X_2 - 1) - f(X_1, X_2)]. \tag{4} $$

We can prove the following reduction holds:

**Theorem 1.** We assume that

1. The degradation function on $X_2$ satisfies $\gamma_2(X_1, 0) \equiv 0$.

2. The degradation function on $X_1$ satisfies $\gamma_1(0, X_2) \equiv 0$, and

$$ \inf_{X_1 \ge 1, X_2 \ge 0} \gamma_1(X_1, X_2) = \underline{\gamma} > 0. $$

3. The production rate of $X_2$ satisfies $\lambda_2(0, X_2) = 0$.

4. The production rate function $\lambda_1$ and $\lambda_2$ are linearly bounded by $X_1 + X_2$.

5. Either $\lambda_1$ or $\lambda_2$ is bounded.
Let $(X_1^N, X_2^N)$ the stochastic process whose generator is $\mathbb{B}_N$ (defined in eq. (4)). Assume that the initial vector $(X_1^N(0), X_2^N(0))$ converges in distribution to $(0, X(0))$, as $N \to \infty$. Then, for all $T > 0$, $(X_1^N(t), X_2^N(t))_{t \ge 0}$ converges in $L^1(0, T)$ (and in $L^p$, $1 \le p < \infty$) to $(0, X(t))$ where $X(t)$ is the stochastic process whose generator is given by

$$ \mathbb{B}_{\infty}\varphi(X) = \lambda_1(0, X) \left( \int_0^{\infty} P_t(\gamma_1(1,.)\varphi(.))(X)dt - \varphi(X) \right) + \gamma_2(0, X) \left[ \varphi(X-1) - \varphi(X) \right], \quad (5) $$

where

$$ P_t g(X) = \mathbb{E}[g(Y(t, X)e^{-\int_0^t \gamma_1(1, Y(s, X))ds})], $$

and $Y(t, X)$ is the stochastic process starting at $X$ at $t = 0$ whose generator is given by

$$ Ag(Y) = \lambda_2(1, Y)(g(Y+1) - g(Y)). $$

**Remark 2.** The first three hypotheses of theorem 1 are the main characteristics of the mRNA-protein system (see example 1). Basically, they impose that quantities remains non-negative, that the first variable has always the possibility to decrease to 0 (no matter the value of the second variable), and that the second variable cannot increase when the first variable is 0. Hence these three hypotheses will guarantee that (with our particular scaling) the first variable converges to 0, and will lead to an intermittent production of the second variable. The last two hypotheses are more technical, and guarantee that the Markov chain is not explosive, and hence well defined for all $t \ge 0$, and that the limiting model is well defined too.

**Remark 3.** The above expression eq. (5) is a generator of a bursting model for a “general bursting size distribution”. For instance, for linear function $\gamma_1(X_1, X_2) = g_1X_1$, and $\lambda_2(X_1, X_2) = k_2X_1$, we have

$$ \begin{aligned} P_t(\gamma_1(\cdot)\varphi(\cdot))(p) &= g_1 P_t(\varphi)(p), \\ &= g_1 \mathbb{E}[\varphi(Y_t^y)e^{-g_1 t}], \\ &= g_1 e^{-g_1 t} \sum_{z \ge y} \varphi(z) \mathbb{P}\{Y_t^y = z\}, \\ &= g_1 e^{-g_1 t} \sum_{z \ge y} \varphi(z) \frac{(k_2 t)^{z-y} e^{-k_2 t}}{(z-y)!}. \end{aligned} $$

It follows by integration integration by parts that

$$ \int_0^\infty P_t(\gamma_1(\cdot)\varphi(\cdot))(y)dt = \frac{g_1}{g_1+k_2} \sum_{z \ge 0} \varphi(z+y) \left(\frac{k_2}{k_2+g_1}\right)^z, $$

which gives then an additive geometric burst size distribution of parameter $p = \frac{k_2}{k_2+g_1}$, as expected Shahrezaei and Swain (2008).
We divide the proof in three steps: moment estimates, tightness and identification of the limit.

## 1.2. Moment estimates

Because production rates are linearly bounded, it is straightforward that with $f(X_1, X_2) = X_1 + X_2$ in eq. (4), there is a constant $C_N$ (that depends on $N$ and other parameters) such that

$$ \mathbb{B}_N f(X_1, X_2) \leq C_N(X_1 + X_2). $$

Then $\mathbb{E}[X_1^N(t) + X_2^N(t)]$ is bounded on any time interval $[0, T]$ and

$$ f(X_1^N(t), X_2^N(t)) - f(X_1^N(0), X_2^N(0)) - \int_0^t \mathbb{B}_N f(X_1^N(s), X_2^N(s))ds $$

is a $L^1$-martingale.

## 1.3. Tightness

Clearly, from the stochastic differential equation on $X_1^N$, we must have $X_1^N(t) \to 0$. We can show in fact that the Lebesgue measure of the set $\{t \le T : X_1^N(t) = 0\}$ converges to 0. Indeed, taking $f(X_1, X_2) = X_1$ in eq. (4), we have

$$ X_1^N(t) - X_1^N(0) - \int_0^t (\lambda_1(X_1^N(s), X_2^N(s)) - N\gamma_1(X_1^N(s), X_2^N(s)))ds \quad (6) $$

is a martingale. Thanks to the lower bound assumption on $\gamma_1$, we have

$$ \underline{\gamma} \mathbb{E} \left[ \int_0^t \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds \right] \le \mathbb{E} \int_0^t \gamma_1(X_1^N(s), X_2^N(s)) ds. $$

Then, by the martingale property, we deduce from (6)

$$ \underline{\gamma} N \mathbb{E} \left[ \int_0^t \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds \right] \le \mathbb{E}[X_1^N(0)] + \int_0^t \mathbb{E}[\lambda_1(X_1^N(s), X_2^N(s))] ds. \quad (7) $$

Now for $X_2^N$ we obtain from eq. (3),

$$ X_2^N(t) \le X_2^N(0) + Y_3 \left( \int_0^t N \mathbf{1}_{\{X_1^N(s) \ge 1\}} \lambda_2(X_1^N(s), X_2^N(s)) ds \right). $$

Let us now distinguish between the two cases.

* Suppose first that $\lambda_2$ is bounded (say by $K$). Then

$$ \mathbb{E}[X_2^N(t)] \le \mathbb{E}[X_2^N(0)] + K N \mathbb{E} \left[ \int_0^t \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds \right]. $$
As $\lambda_1$ is linearly bounded (say by $K$) by $X_1^N + X_2^N$, the upper bound eq. (7) becomes

$$ \underline{\gamma} N \mathbb{E} \left[ \int_0^t \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds \right] \le \mathbb{E} [X_1^N(0)] + K \int_0^t \left( \mathbb{E}[X_1^N(s)] + \mathbb{E}[X_2^N(s)] \right) ds. $$

Finally, with eq. (6), it is clear that

$$ \mathbb{E}[X_1^N(t)] \le \mathbb{E}[X_1^N(0)] + K \int_0^t \left( \mathbb{E}[X_1^N(s)] + \mathbb{E}[X_2^N(s)] \right) ds. $$

Hence, with the three last inequalities, we can conclude by the Grönwall lemma that $\mathbb{E}[X_2^N(t)]$ is bounded on $[0, T]$, uniformly in $N$. Then

$$ N \mathbb{E} \left[ \int_0^T \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds \right] $$

is bounded and $X_1^N \to 0$ in $L^1([0, T], \mathbb{N})$. By the law of large number, $\frac{1}{N}Y_3(N)$ is almost surely convergent, and hence almost surely bounded. We deduce then there exists a random variable $C$ such that

$$ X_2^N(t) \le X_2^N(0) + N C \int_0^t \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds, $$

almost everywhere. By Grönwall lemma and Markov inequality

$$ \mathbb{P} \left\{ \sup_{t \in [0, T]} X_2^N(t) \ge M \right\} \to 0 $$

as $M \to \infty$, uniformly in $N$.

* Now suppose $\lambda_1$ is bounded (say $K$). By the martingale eq. (6) (and the same lower bound hypothesis on $\gamma_1$, it is clear that

$$ N \mathbb{E} \left[ \int_0^T \mathbf{1}_{\{X_1^N(s) \ge 1\}} ds \right] $$

is bounded and $X_1^N \to 0$ in $L^1([0, T], \mathbb{N})$. Now, let us denote $U^N(t) = \frac{1}{N}X_1^N(t)$, $V^N = \frac{1}{N}X_2^N(t)$ and $W^N = N\mathbf{1}_{\{X_1^N(t) \ge 1\}}$ (which is then bounded in $L^1([0, T])$). From eq. (3), and from the linear bound on $\lambda_2$ (say by $K$)

$$ V^N(t) \le V^N(0) + \frac{1}{N} Y_3 \left( \int_0^t N K W^N (U^N(s) + V^N(s)) ds \right). $$

Then, still by the law of the large number there exists a random variable $C$ such that

$$ V^N(t) \le V^N(0) + C \int_0^t W^N(U^N(s) + V^N(s)) ds, $$
and hence

$$X_2^N(t) \leq X_2^N(0) + C \int_0^t W^N(X_1^N(s) + X_2^N(s))ds.$$

By Grönwall lemma,

$$\sup_{[0,T]} X_2^N(t) \leq (X_1^N(0) + X_2^N(0)) \exp \left( C \int_0^t W^N(s)ds \right),$$

which is then bounded, uniformly in N.

For any subdivision of $[0, T]$, $0 = t_0 < t_1 < \dots < t_n = T$,

$$\sum_{i=0}^{n-1} | X_2^N(t_{i+1}) - X_2^N(t_i) | \leq \sum_{i=0}^{n-1} Y_3 \left( \int_{t_i}^{t_{i+1}} N \mathbf{1}_{\{X_1^N(s) \geq 1\}} \lambda_2(X_1^N(s), X_2^N(s)) ds \right) \\ \leq Y_3 \left( \int_0^T N \mathbf{1}_{\{X_1^N(s) \geq 1\}} \lambda_2(X_1^N(s), X_2^N(s)) ds \right)$$

so by a similar argument as above, we also get the tightness of the BV norm

$$\mathbb{P}\{\|X_2^N\|_{[0,T]} \geq K\} \to 0$$

as $K \to 0$, independently in N. Then $X_2^N$ is tight in $L^p([0, T])$, for any $1 \leq p < \infty$ (Giusti (1984)).

## 1.4. Identification of the limit

We choose an adherence value $(0, X_2(t))$ of the sequence $(X_1^N(t), X_2^N(t))$ in $L^1([0, T]) \times L^p([0, T])$. Then a subsequence (again denoted by) $(X_1^N(t), X_2^N(t))$ converge to $(0, X_2(t))$, almost surely and for almost $t \in [0, T]$. We are looking for test-functions such that

$$f(X_1^N(t), X_2^N(t)) - f(X_1^N(0), X_2^N(0)) - \int_0^t \mathbb{B}_N f(0, X_2^N(s)) 1_{X_1^N(s)=0} ds \\ - \int_0^t \mathbb{B}_N f(X_1^N(s), X_2^N(s)) 1_{X_1^N(s) \geq 1} ds$$

is a martingale and $\mathbb{B}_N f(X_1^N(s), X_2^N(s))$ is bounded independently of N when $X_1 \geq 1$. The following choice is inspired by Crudu et al. (2012). We introduce the stochastic process $Y_t^{x,y}$, starting at y and whose generator is

$$A^x g(y) = \lambda_2(x, y) [g(y+1) - g(y)],$$

for any $x \geq 1$. and we introduce the semigroup $P_t^x$ defined on bounded function, for any $x \geq 1$, by

$$P_t^x g(y) = \mathbb{E} \left[ g(Y_t^{x,y}) e^{-\int_0^t \gamma_1(x, Y_s^{x,y}) ds} \right]. \quad (8)$$
Then the semigroup $P_t^x$ satisfies the equation

$$ \frac{dP_t^x g(y)}{dt} = A^x P_t^x g(y) - \gamma_1(x, y) P_t^x g(y). $$

Now for any bounded function $g$, define recursively

$$ \begin{aligned} f(0, y) &= g(y), \\ f(x, y) &= \int_{0}^{\infty} P_{t}^{x}(\gamma_{1}(x, \cdot)f(x-1, \cdot))(y)dt. \end{aligned} $$

Such a test function is well defined by the assumption on $\gamma_1$. We then verify that

$$ \begin{aligned} \mathbb{B}_N f(0, y) &= \lambda_1(0, y) \left( \int_0^\infty P_t^1(\gamma_1(1, \cdot)g(\cdot))(y)dt - g(y) \right) + \gamma_2(0, y) [g(y-1) - g(y)], \\ \mathbb{B}_N f(x, y) &= \lambda_1(x, y) [f(x+1, y) - f(x, y)] + \gamma_2(x, y) [f(x, y-1) - f(x, y)]. \end{aligned} $$

Indeed, for any $x \ge 1$,

$$ \begin{aligned} &A^x f(x, y) - \gamma_1(x, y) f(x, y) \\ &= \int_0^\infty A^x P_t^x(\gamma_1(x, \cdot)f(x-1, \cdot))(y) - \gamma_1(x, y) P_t^x(\gamma_1(x, \cdot)f(x-1, \cdot))(y)dt, \\ &= \int_0^\infty \frac{d}{dt} P_t^x(\gamma_1(x, \cdot)f(x-1, \cdot))(y)dt, \\ &= \lim_{t \to \infty} P_t^x(\gamma_1(x, \cdot)f(x, \cdot))(y) - \gamma_1(x, y)f(x-1, y), \\ &= -\gamma_1(x, y)f(x-1, y). \end{aligned} $$

Then

$$ \lambda_2(x, y) [f(x, y + 1) - f(x, y)] + \gamma_1(x, y) [f(x - 1, y) - f(x, y)] = 0. $$

Hence $\mathbb{B}_N f(x, y)$ is independent of N, and, taking the limit $N \to \infty$ in

$$ f(X_1^N(t), X_2^N(t)) - f(X_1^N(0), X_2^N(0)) - \int_0^t \mathbb{B}_N f(X_1^N(s), X_2^N(s))ds, $$

we deduce

$$ g(X_2(t)) - g(X_2(0)) - \int_0^t \mathbb{B}_\infty g(X_2) $$

is a martingale where

$$ \mathbb{B}_{\infty} g(y) = \lambda_1(0, y) \left( \int_0^\infty P_t^1(\gamma_1(1, \cdot)g(\cdot))(y)dt - g(y) \right) + \gamma_2(0, y) [g(y-1) - g(y)]. $$

**Uniqueness** Due to assumption on $k_1$ and $k_2$, the limiting generator defines a pure-jump Markov process in $\mathbb{N}$ which is not explosive. Uniqueness of the martingale then follows classically.
## 2. Continuous-state bursting model

The model we consider now is a continuous state-space model that explicitly assume the production of several molecules instantaneously, through a jump Markov process, in agreement with experimental observations (Golding et al. (2005); Raj et al. (2006)). In line with experimental observations, it is standard to assume a Markovian hypothesis (an exponential waiting time between production jumps) and that the jump sizes are exponentially distributed (geometrically in the discrete case) as well (Suter et al. (2011)). The intensity of the jumps can be a linearly bounded function, to allow for self-regulation.

For simplicity, we will only consider the standard model of gene expression, that is with linear degradation rates and the production rate of the second variable is linear with respect to the first variable (as in example 1). Note that more general rates as in the previous section could be considered as well. Here, we ask the question of what is the correct scaling so that the bursting production term is transmitted from the first variable to the second one, when the first variable is eliminated through an adiabatic limit. The propagation of bursting property in a gene network is an important question in molecular biology Kaern et al. (2005).

This section is structured as follows. We first present the model in the rest of this paragraph, then state the results in subsection 2.1, and divide the proofs in the remaining three subsection 2.2,2.3,2.4.

Let $x_1$ and $x_2$ denote the concentrations of mRNA and protein respectively. A simple model of single gene expression with bursting transcription is given by

$$ \frac{dx_1}{dt} = -g_1x_1 + \hat{N}(h, \lambda_1(x_2)), \quad (9) $$

$$ \frac{dx_2}{dt} = -g_2x_2 + k_2x_1. \quad (10) $$

Here $g_1$ and $g_2$ are the degradation rates for the mRNA and protein respectively, $k_2$ is the mRNA translation rate, and $\hat{N}(h, \lambda_1(x_2))$ describes the transcription that is assumed to be a compound Poisson white noise occurring at a rate $\lambda_1(x_2)$ with a non-negative jump size $\Delta x_1$ distributed with density $h$.

The equations (9)-(10) are a short hand notation for

$$ x_1(t) = x_1(0) - \int_0^t g_1 x_1(s^-) ds + \int_0^t \int_0^\infty \int_0^\infty 1_{\{r \le \lambda_1(x_2(s^-))\}} zN(ds, dz, dr), \quad (11) $$

$$ x_2(t) = x_2(0) - \int_0^t g_2 x_2(s^-) ds + \int_0^t k_2 x_1(s^-) ds. \quad (12) $$

where $X_{s-} = \lim_{t \to s^-} X(t)$, and $N(ds, dz, dr)$ is a Poisson random measure on $(0, \infty) \times [0, \infty)^2$ with intensity $dsh(z)dzdr$, where $s$ denotes the times of the jumps, $r$ is the state-dependency in an acceptance/rejection fashion, and $z$ the
jump size. Note that $(x_1(t))$ is a stochastic process with almost surely finite variation on any bounded interval $(0, T)$, so that the last integral is well defined as a Stieltjes-integral.

**Hypothesis 2.** The following discussion is valid for general rate functions $\lambda_1$ and density functions $h(\cdot)$ that satisfy

* $\lambda_1 \in C^1$, $\lambda_1$ is globally lipschitz and linearly bounded with

$$0 \le \lambda_1(x) \le c + Kx.$$

* $h \in C^0$ and $\int_0^\infty xh(x)dx < \infty$.

For such a general density function $h$, we denote the average burst size by

$$b = \int_0^\infty xh(x)dx. \qquad (13)$$

*Remark 4.* Hill functions are often used to model gene self-regulation, so that $\lambda_1$ is given by

$$\lambda_1(x_2) = \frac{1 + x_2^\alpha}{L + D x_2^\alpha}$$

where $L, D$ are positive parameters and $\alpha$ is a positive integer (see Mackey, Tyran-Kamińska and Yvinec (2011) for more details). An exponential distribution of the bursting transcription is often used in modeling gene expression, in accordance with experimental findings (Xie et al. (2008)), so that the density function $h$ is given by

$$h(x) = \frac{1}{b}e^{-x/b},$$

with $b$ the average burst size.

If $\lambda_1(x_2) \equiv k_1$ is independent of the state $x_2$, the average transcription rate is $bk_1$, and the asymptotic average mRNA and protein concentrations are

$$x_1^{\text{eq}} := \mathbb{E}[x_1(t \to \infty)] = \frac{bk_1}{g_1}.$$

$$x_2^{\text{eq}} := \mathbb{E}[x_2(t \to \infty)] = \frac{k_2}{g_2} x_1^{\text{eq}} = \frac{bk_1 k_2}{g_1 g_2}.$$

## 2.1. Statement of the results

Although Equations (9)-(10) are simple, they are not analytically solvable. Hence, for practical use to interpret experimental data, and to avoid numerical simulations, one looks for a reduced, analytically solvable, one-dimensional equation. In the following discussion, we consider the situation when mRNA degradation is a fast process, i.e. $g_1$ is “large enough”, but the average equilibrium protein concentration $x_2^{\text{eq}}$ remains unchanged. In most organisms and for most genes, mRNA has a smaller lifetime than protein (Schwanhäusser et al.
(2011)). In what follows, we denote by $g_1^n$, $k_2^n$ sequences of parameters, $\lambda_1^n$ sequence of functions and $h^n$ sequence of density function that will replace $g_1$, $k_2$, $\lambda_1$, $h$ in (11)-(12) and satisfy hypothesis 2. We then denote $(x_1^n, x_2^n)$ its associated solution. We will always assume one of the following three scaling relations:

(S1) Frequent production rate of mRNA, namely $g_1^n = ng_1$, $\lambda_1^n = n\lambda_1$, and $k_2^n = k_2 h^n = h$ are independent of $n$;

(S2) Large burst of mRNA, namely $g_1^n = ng_1$, $h^n(z) = \frac{1}{n}h(\frac{z}{n})$ and $\lambda_1^n = \lambda_1, k_2^n = k_2$ remain unchanged;

(S3) Large production rate of protein, namely $g_1^n = ng_1$, $k_2^n = nk_2$, and $\lambda_1^n = \lambda_1 h^n = h$ are independent of $n$;

These three different scaling are then associated with different behaviors of the biological systems given by $(x_1, x_2)$. As different genes may have different kinetics, each one of the possible scaling are of importance (Suter et al. (2011); Schwanhäusser et al. (2011)).

In this section we determine an effective reduced equation for equation (10) for each of the three scaling conditions (S1)-(S3). In particular, we show that under assumption (S1), equation (10) can be approximated by the deterministic ordinary differential equation

$$ \frac{dx_2}{dt} = -g_2x_2 + \lambda_2(x_2) \quad (14) $$

where

$$ \lambda_2(x_2) = bk_2\lambda_1(x_2)/g_1. $$

We further show that under the scaling relations (S2) or (S3), equation (10) can be reduced to the stochastic differential equation

$$ \frac{dx_2}{dt} = -g_2x_2 + \dot{N}(\bar{h}, \lambda_1(x_2)). \quad (15) $$

where $\bar{h}$ is a suitable density function in the jump size $\Delta x_2$ (to be detailed below).

We first explain, using some heuristic arguments, the differences between the three scaling relations and the associated results. When $n \to \infty$, $g_1^n \to \infty$ and applying a standard quasi-equilibrium assumption we have

$$ \frac{dx_1^n}{dt} \approx 0, $$

which yields

$$ x_1^n(t) \approx \frac{1}{g_1^n} \dot{N}(h^n(\cdot), \lambda_1^n(x_2^n)) = \dot{N}(g_1^n h^n(g_1^n \cdot), \lambda_1^n(x_2^n)), $$

and therefore the second equation (10) becomes

$$ \begin{aligned} \frac{dx_2^n}{dt} &\approx -g_2x_2^n + \frac{k_2^n}{g_1^n} \dot{N}(h^n(\cdot), \lambda_1^n(x_2^n)), \\ &\approx -\gamma_2 x_2^n + \dot{N}\left(\frac{g_1^n}{k_2^n} h^n\left(\frac{g_1^n}{k_2^n}\right), \lambda_1^n(x_2^n)\right). \end{aligned} $$
Hence in (15), $\bar{h}(x_2) = (k_2/g_1)^{-1}h((k_2/g_1)^{-1}x_2)$ under the scaling (S2) and (S3). Furthermore, we note that the scaling (S2) also implies $nh^n(n \cdot) = h(\cdot)$, while in (S1), $nh^n(n \cdot) = nh(n \cdot)$ so that the jumps become more frequent and smaller.

We denote $(D[0, \infty), S)$ the cadlag function space of function defined on $[0, \infty)$ at values in $\mathbb{R}^+$ with the usual Skorohod topology (Jacod and Shiryaev (1987)). Similarly $(D[0, T], J)$ is the cadlag function space on $[0, T]$, with the Jakubowski topology (Jakubowski (1997)). Also, $L^p[0, T)$ the space of $L^p$ integrable function on $[0, T]$, with $T > 0$, which we endowed with total variation norm (Giusti (1984)), and $M(0, \infty)$ is the space of real measurable function on $[0, \infty)$ with the metric (Kurtz (1991))

$$d(x, y) = \int_{0}^{\infty} e^{-t} \max\{1, |x(t) - y(t)|\} dt.$$

Our main results can be stated as follows

**Theorem 3.** Consider the equations (9)-(10) and assume Hypothesis 2. If the scaling (S1) is satisfied, i.e., $k_1^n = nk_1$, and if $x_2^n(0) \to x_2^0$, then when $n \to \infty$,

1. The stochastic process $x_1^n(t)$ does not converge in any functional sense;

2. The stochastic process $x_2^n(t)$ converges in law in $(D[0, \infty), S)$ towards the deterministic solution of the ordinary differential equation

$$\frac{dx_2}{dt} = -g_2x_2 + \lambda_2(x_2), \quad x_2(0) = x_2^0, \qquad (16)$$

where

$$\lambda_2(x_2) = bk_2\lambda_1(x_2)/g_1.$$

**Theorem 4.** Consider the equations (9)-(10) and assume Hypothesis 2. If the scaling (S2) is satisfied, i.e., $h^n(z) = \frac{1}{n}h(\frac{z}{n})$, and if $x_2^n(0) \to x_2^0$, then when $n \to \infty$,

1. The stochastic process $\frac{x_1^n(t)}{n}$ converges in law in $L^p$, $1 \le p < \infty$ and in $(D[0, T], J)$ to the (deterministic) fixed value 0;

2. The stochastic process $x_2^n(t)$ converges in law in $L^p$, $1 \le p < \infty$ and in $(D[0, T], J)$ to the stochastic process defined by the solution of the stochastic differential equation

$$\frac{dx_2}{dt} = -g_2x_2 + \hat{N}(\bar{h}, \lambda_1), \quad x_2(0) = x_2^0 \ge 0, \qquad (17)$$

where $\bar{h}(x_2) = (k_2/g_1)^{-1}h((k_2/g_1)^{-1}x_2)$.

**Theorem 5.** Consider the equations (9)-(10) and assume Hypothesis 2. If the scaling (S3) is satisfied, i.e., $k_2^n = nk_2$, and if $x_2^n(0) \to x_2^0$, then when $n \to \infty$,

1. The stochastic process $x_1^n(t)$ converges in law in $L^p$, $1 \le p < \infty$ and in $(D[0, T], J)$ to the (deterministic) fixed value 0;
2. The stochastic process $x_2^n(t)$ converges in law in $L^p$, $1 \le p < \infty$ and in $(D[0, T], J)$ to the stochastic process determined by the solution of the stochastic differential equation

$$ \frac{dx_2}{dt} = -\gamma_2 x_2 + \dot{N}(\bar{h}, \varphi), \quad x_2(0) = x_2^0 \ge 0, $$

$$ \text{where } \bar{h}(x_2) = (k_2/g_1)^{-1}h((k_2/g_1)^{-1}x_2). $$

Remark 5. Note that scalings (S2) and (S3) give similar results for the equation governing the protein variable $x_2(t)$ but very different results for the asymptotic stochastic process related to the mRNA. In particular, in Theorem 4, very large bursts of mRNA are transmitted to the protein, where in Theorem 5, very rarely is mRNA present but when present it is efficiently synthesized into a burst of protein.

In the rest of this paper, we provide proofs of the results mentioned above, using martingale techniques. In a companion paper Yvinec et al. (2012), we use partial differential techniques to prove similar results (see also Haseltine and Rawlings (2005); Zeron and Santillán (2010); Santillán and Qian (2011)).

The proofs of the three theorems above are divided in three steps. In section 2.2 we first recall generator properties and derive moment estimates associated to (9)-(10). In section 2.3 we show the tightness result for all three theorems. We then identify the limit using a martingale approach in section 2.4.

## 2.2. General properties and moment estimates

We first summarize the important background results on the stochastic processes used in the next.

**One dimensional equation** For the one-dimensional stochastic differential equation (15) perturbed by a compound Poisson white noise, of intensity $\lambda(\cdot)$ and jump size distribution $\bar{h}(\cdot)$, the extended generator of the stochastic process $(x_2(t))_{t \ge 0}$ is (Davis, 1984, Theorem 5.5), for any $f \in \mathcal{D}(\mathcal{A}_\infty)$,

$$ \mathcal{A}_1 f(x) = -g_2 x \frac{df}{dx} + \lambda(x) \left( \int_x^\infty \bar{h}(z-x) f(z) dz - f(x) \right) $$

$$ \mathcal{D}(\mathcal{A}_1) = \{ f \in \mathcal{M}(0, \infty) : t \mapsto f(xe^{-\gamma_2 t}) \text{ is absolutely continuous for } t \in \mathbb{R}^+ \text{ and } \mathbb{E} \sum_{T_i \le t} |f(x_2(T_i)) - f(x_2(T_i^-))| < \infty \text{ for all } t \ge 0 \} $$

where $\mathcal{M}(0, \infty)$ denotes a Borel-measurable function of $(0, \infty)$ and the times $T_i$ are the instants of the jump of $x_2$. It is an extended domain containing all functions that are sufficiently smooth along the deterministic trajectories between the jumps, and with a bounded total variation induced by the jumps.
For any $f \in \mathcal{D}(\mathcal{A}_1)$, we have

$$ \frac{d}{dt} \mathbb{E} f(x_2(t)) = \mathbb{E} \mathcal{A}_1(f(x_2(t))). $$

**Two dimensional equation** Consideration the two-dimensional stochastic differential equation (9)-(10) perturbed by a compound Poisson white noise, of intensity $\lambda_1(x_2)$ and jump size distribution $h$ follows along similar lines. Its infinitesimal generator and extended domain are

$$ \mathcal{A}_2 g(x_1, x_2) = -g_1 x_1 \frac{\partial g}{\partial x_1} + (k_2 x_1 - g_2 x_2) \frac{\partial g}{\partial x_2} \\ + \lambda_1(x_2) \left( \int_{x_1}^{\infty} h(z-x_1) g(z, x_2) dz - g(x_1, x_2) \right), \quad (18) $$

$$ \mathcal{D}(\mathcal{A}_2) = \{g \in \mathcal{M}((0, \infty)^2) : t \mapsto g(\phi_t(x_1, x_2)) \text{ is absolutely continuous for } t \in \mathcal{R}^+ \text{ and } \mathbb{E} \sum_{T_i \le t} |g(x_1(T_i), x_2(T_i)) - g(x_1(T_i^-), x_2(T_i^-))| < \infty \text{ for all } t \ge 0\} \quad (19) $$

where $\phi_t$ is the deterministic flow given by the deterministic part of equations (9)-(10), namely

$$ \begin{aligned} \frac{dx_1}{dt} &= -g_1 x_1, \\ \frac{dx_2}{dt} &= -g_2 x_2 + k_2 x_1. \end{aligned} $$

For any $f \in \mathcal{D}(\mathcal{A}_2)$, we have

$$ \frac{d}{dt} \mathbb{E} f(x_1(t), x_2(t)) = \mathbb{E} \mathcal{A}_2(f(x_1(t), x_2(t))). \quad (20) $$

Using the stochastic differential equations (11)-(12), we can deduce moment estimates, needed to be able to use unbounded test function (namely $f(x_1, x_2) = x_1$ and $f(x_1, x_2) = x_2$) in the martingale formulation. By taking the mean into (11)-(12), neglecting negatives values and using hypothesis 2,

$$ \begin{aligned} 0 &\le \mathbb{E}[x_1(t)] \le \mathbb{E}[\int_0^t b\lambda_1(x_2(s))ds] \le \int_0^t b(c + K\mathbb{E}[x_2(s)])ds \\ 0 &\le \mathbb{E}[x_2(t)] \le \int_0^t k_2\mathbb{E}[x_1(s)]ds \end{aligned} $$

where we note $b = \mathbb{E}[h] = \int_0^\infty zh(z)dz$. By Gronwall inequalities, there exist a constant $C$ such that

$$ \begin{aligned} \mathbb{E}[\sup_{t \in [0,T]} x_1(t)] &\le C(\mathbb{E}[x_1(0)] + e^{CT}) \\ \mathbb{E}[\sup_{t \in [0,T]} x_2(t)] &\le C(\mathbb{E}[x_2(0)] + e^{CT}) \end{aligned} \quad (21) $$
Then we claim that $f(x_1, x_2) = x_1$ is in the domain of the generator $\mathcal{A}_2$. We only have to verify (see Eq (19))

$$ \mathbb{E} \sum_{T_i \le t} |x_1(T_i) - x_1(T_i^-)| < \infty \text{ for all } t \ge 0. $$

By equation (11)

$$ \begin{aligned} \mathbb{E} \sum_{T_i \le t} |x_1(T_i) - x_1(T_i^-)| &= \mathbb{E} \int_0^t \int_0^\infty \int_0^\infty 1_{\{r \le \lambda_1(x_2(s^{-}))\}} zN(ds, dz, dr), \\ &\le b \mathbb{E} \left[ \int_0^t c + Kx_2(s)ds \right]. \end{aligned} $$

which is finite according to the previous estimates.

### 2.3. Tightness

**S1** We first show the tightness property for the scaling (S1) corresponding to theorem 3. In such case $x_1^n$ does no converge in any functional sense because it fluctuates very fast, as more and more jumps appears of size that stay of order 1 (given by $h$). However, $\mathbb{E}[x_1^n(t)]$ remains bounded, $\frac{x_1^n}{n}$ goes to 0, and by eq. (12),

$$ |x_2^n(t)| \le |x_2^n(0)| + \int_0^t k_2 |x_1^n(s)| ds. $$

For any $n$, let $N_n$ be the compound Poisson process associated to (11), with $\{T_{n,i}\}_{i=1}^\infty$ the jump times which occur at a rate $n\lambda_1(x_2(s)^n)$, and $\{Z_{n,i}\}_{i=1}^\infty$ the jump sizes that are iid random variables with density $h$ (with the convention $T_{n,0} = 0$ and $Z_{n,0} = X_0$),

$$ N_n(t) = \sum_{T_{n,i} \le t} Z_{n,i}. $$

Then

$$ x_1^n(t) = \sum_{T_{n,i} \le t} Z_{n,i} e^{-ng_1(t-T_{n,i})} \mathbf{1}_{\{t \ge T_{n,i}\}}. $$

By integration,

$$ \int_0^t x_1^n(s) ds = \sum_{T_{n,i} \le t} Z_{n,i} \frac{1}{ng_1} (1 - e^{-g_1(t-T_{n,i})}) \mathbf{1}_{\{t \ge T_{n,i}\}}. $$

Then,

$$ x_2^n(t) \le x_2^n(0) + \int_0^t k_2 x_1^n(s) ds \le x_2^n(0) + \frac{k_2}{ng_1} \sum_{T_{n,i} \le t} Z_{n,i}. $$

Finally we deduce, by definition of the compound Poisson process,

$$ x_2^n(t) \le x_2^n(0) + \frac{k_2}{ng_1} N_n(t). $$
Now, by a time change, there exists a process Y such that

$$N_n(t) = Y\left(\int_0^t n\lambda_1(x_2^n(s))ds\right),$$

where Y is unit rate compound Poisson process of jump size iid (with density h). By the law of large number, $\frac{1}{n}Y(nt)$ is almost surely convergent, and hence almost surely bounded. We deduce then there exists a random variable C such that

$$x_2^n(t) \le x_2^n(0) + \frac{k_2}{g_1} C \int_0^t \lambda_1(x_2^n(s))ds.$$

By Gronwall lemma and Markov inequality

$$\mathbb{P}\left\{\sup_{t \in [0,T]} x_2^n(t) \ge M\right\} \to 0,$$

as $M \to \infty$ and uniformly in $n$. Similarly, for any $t_1, t_2 \in [0, T]$,

$$|x_2^n(t_2) - x_2^n(t_1)| \le \frac{k_2}{ng_1} |N_n(t_2) - N_n(t_1)|.$$

Again, $N_n(t_2) - N_n(t_1) = Y\left(\int_{t_1}^{t_2} n\lambda_1(x_2^n(s))ds\right)$ and, still by the law of large number

$$|x_2^n(t_2) - x_2^n(t_1)| \le \frac{k_2}{g_1} C \int_{t_1}^{t_2} \lambda_1(x_2^n(s))ds,$$

so that, for any $\varepsilon > 0$

$$\lim_{\theta \to 0} \limsup_n \sup_{S_1 \le S_2 \le S_1 + \theta} \mathbb{P}\left\{ |x_2^n(S_2) - x_2^n(S_1)| \ge \varepsilon \right\} = 0,$$

where the supremum is over stopping times bounded by T. Then by Aldous' tightness criteria ((Jacod and Shiryaev, 1987, thm 4.5 p 356)), $x_2^n$ is tight in $(D[0, \infty), S)$.

**S3** Now we show the tightness property for the scaling (S3) corresponding to theorem 5, with $k_2^n = nk_2$. In such case $x_1^n$ converges to 0 in $L^1$, and we get a control over $n \int_0^t x_1^n(s)ds$. Indeed using $g(x_1, x_2) = x_1$ in (18), we get

$$x_1^n(t) - x_1^n(0) - \int_0^t (-ng_1x_1^n(s) + b\lambda_1(x_2(s)^n))ds,$$

is a martingale so that due to Hypothesis 2,

$$g_1 \mathbb{E}[n \int_0^t x_1^n(s) ds] \le \mathbb{E}[x_1^n(0)] + b(ct + K \int_0^t \mathbb{E}[x_2^n(s)] ds)$$

By eq. (12),

$$x_2^n(t) \le \mathbb{E}[x_2^n(0)] + k_2 n \int_0^t x_1^n(s) ds.$$
then

$$ \sup_{t \in [0,T]} x_2^n(t) \leq \mathbb{E}[x_2^n(0)] + k_2 n \int_0^T x_1^n(s) ds. $$

Reporting into the estimates for $x_1^n$ yields

$$ g_1 \mathbb{E}[n \int_0^t x_1^n(s) ds] \leq \mathbb{E}[x_1^n(0)] + b(ct + K(\mathbb{E}[x_2^n(0)] + tk_2 n \int_0^t \mathbb{E}[x_1^n(s)] ds)), \\ \leq C_T^1 + C_T^2 \mathbb{E}[n \int_0^t x_1^n(s) ds], $$

for two constants $C_T^1$, $C_T^2$ that depends solely on $T$. Then $\mathbb{E}[n \int_0^t x_1^n(s) ds]$ is bounded uniformly in $n$ so that $x_1^n$ converges to 0 in $L^1$ and

$$ \mathbb{P}\{\sup_{t \in [0,T]} x_2^n(t) \geq M\} \to 0 $$

as $M \to \infty$ and uniformly in $n$. Now for any subdivision of $[0, T]$, $0 = t_0 < t_1 < \dots < t_n = T$,

$$ \sum_{i=0}^{n-1} |x_2^n(t_{i+1}) - x_2^n(t_i)| \leq \mathbb{E}[x_2^n(0)] + k_2 n \int_0^t x_1^n(s) ds, $$

so that we also get the tightness of the BV norm,

$$ \mathbb{P}\{\|x_2^n\|_{[0,T]} \geq M\} \to 0, $$

as $M \to 0$, independently in $n$. Then $x_2^n$ is tight in $L^p([0, T])$, for any $1 \leq p < \infty$ (Giusti (1984)) and also, by a similar criteria, in $(D[0, T], J)$ (Jakubowski (1997)).

**S2** Now we show the tightness property for the scaling (S2) corresponding to theorem 4, with $h^n = \frac{1}{n}h(\frac{1}{n})$. Remark that on such case, denoting $z^n = \frac{x_1^n}{n}$, the variables $(z^n, x_2^n)$ satisfies (11)-(12) with the (S3) scaling, so we already now that $x_2^n$ is tight in $L^p([0, T])$, for any $1 \leq p < \infty$.

For $x_1^n$, formally, note that each jumps yields a contribution for $\int x_1^n$ of $\frac{b}{g_1}$ so there's no hope for a convergence to 0 in $L^1$. However, we still have

$$ x_1^n(t) = \sum_{T_{n,i} \leq t} Z_{n,i} e^{-ng_1(t-T_{n,i})} 1_{t \geq T_{n,i}}. $$

where $T_{n,i}$ appears with rate $\lambda_1(x_2^n(s))$, and $\{Z_{n,i}\}_{i=1}^\infty$ are iid random variables with density $h^n$. Then

$$ x_1^n(t) \leq \sum_{T_{n,i} \leq t} Z_{n,i} \left( \mathbf{1}_{\{[T_{n,i}, T_{n,i}+\frac{1}{\sqrt{n}}]\}} + e^{-ng_1 \frac{1}{\sqrt{n}}} \mathbf{1}_{t \geq T_{n,i}} \right). $$
But for $M > 0$, by Markov inequality,

$$ \mathbb{P}\{Z_{n,i}e^{-\sqrt{n}g_1} > M\} \le \frac{nb}{Me^{\sqrt{n}g_1}} \le \varepsilon, $$

for any $\varepsilon$ and $n$ sufficiently large. Then, conditioning by the jump times,

$$ \int_0^t \mathbb{P}\{x_1^n(s) > M \mid T_{n,i}\} \le \sum_{T_{n,i} \le t} \frac{1}{\sqrt{n}} \mathbf{1}_{\{t \ge T_{n,i}\}} + \sum_{T_{n,i} \le t} \varepsilon(t - T_{n,i}) \mathbf{1}_{\{t \ge T_{n,i}\}} \le \varepsilon. $$

for $n$ large. Because $\int_0^t x_2^n(s)ds$ has been shown to be bounded independently of $n$, we can drop the conditioning, and $\int_0^t \mathbb{P}\{x_1^n(s) > M\}$ is arbitrary small. We show also similarly that

$$ \lim_{h \to 0} \sup_n \int_0^T \max(1, |x_1^n(t+h) - x_1^n(t)|)dt = 0, $$

so that $x_1^n$ is tight in $M(0, \infty)$ ((Kurtz, 1991, thm 4.1)).

## 2.4. Identification with the martingale problem

The three theorems below can be proved using martingale techniques, with similar spirit. For each scaling, the generator $A_2^n$ can be decomposed into a fast component, or order $n$, and a slow component, of order 1. In each case, one need to find particular condition to ensure that the fast component vanishes. For the scaling $(S1)$, the fast component acts only in the first variable, so ergodicity of this component will ensure that it vanishes, as in averaging theorems Kurtz (1992). For the other two scaling, the fast component acts on both variables, and we will have to find the particular relation between both variable that ensures this component vanishes.

**Proof of Theorem 3** For any $B \in \mathbb{B}(\mathbb{R}_+)$, $t > 0$, we define the occupation measure

$$ V_1^n(B \times [0, t]) = \int_0^t \mathbf{1}_{\{B\}}(x_1^n(s))ds, $$

and we identify $V_1^n$ as a stochastic process with value in the space of finite measure on $\mathbb{R}^+$. Because $E[x_1^n(t)]$ remains bounded uniformly in $n$ on any $[0, T]$, it is stochastically bounded and $V_1$ then satisfies Aldous criteria of tightness. Now take a test function $f$ that depends only on $x_1$, so that

$$ A_2^n f(x_1) = nC_{x_2}f(x_1), $$

with

$$ C_{x_2}f(x_1) = -g_1x_1f'(x_1) + \lambda_1(x_2) \left( \int_{x_1}^{\infty} h(z-x_1)f(z)dz - f(x_1) \right). $$
Then

$$M_t^n = f(x_1^n(t)) - f(x_1^n(0)) - n \int_{\mathbb{R}_+} \int_0^t C_{x_2(s)} f(x_1) V_1^n(dx_1 \times ds)$$

is a martingale. Dividing by $n$, for any limiting point $(V_1, x_2)$, we must have, for any $f \in C_b(\mathbb{R}_+)$,

$$\mathbb{E}[\int_{\mathbb{R}_+} \int_0^t C_{x_2(s)} f(x_1) V_1(dx_1 \times ds)] = 0.$$

Because for any $x_2$, the generator $C_{x_2}$ is (exponentially) ergodic, $V_1$ is uniquely determined by the invariant measure associated to $C_{x_2}$. In particular, for any $t > 0$

$$\int_{\mathbb{R}_+} \int_0^t x_1 V_1^n(dx_1 \times ds) \rightarrow \int_0^t \frac{b}{g_1} \lambda_1(x_2(s)) ds.$$

Then for $f$ that depends only on $x_2$,

$$f(x_2^n(t)) - f(x_2^n(0)) - \int_{\mathbb{R}_+} \int_0^t (k_2 x_1 - g_2 x_2^n(s)) f'(x_2^n(s)) V_1^n(dx_1 \times ds)$$

converges to

$$f(x_2(t)) - f(x_2(0)) - \int_0^t \left( \frac{bk_2}{g_1} \lambda_1(x_2(s)) - g_2 x_2(s) \right) f'(x_2(s)) ds$$

Due to the assumption on $\lambda_1$, there exists a unique solution associated to the (deterministic) equation 14 so $x_2$ is uniquely determined.

**Proof of Theorem 5** We've already seen that $x_1^n$ converges to 0 in $L^1([0, T])$ and $x_2^n$ is tight in $L^p([0, T])$. We then take a subsequence $(x_1^n(t), x_2^n(t))$ that converges to $(0, x_2(t))$, almost surely and for almost $t \in [0, T]$. Then we consider the fast component of the generator $A_2^n$, given in this case by

$$-g_1 x_1 \frac{\partial f}{\partial x_1} + k_2 x_1 \frac{\partial f}{\partial x_2}.$$

This defines a transport equation. Starting at $(x_1, x_2)$ at time 0, the asymptotic value of the flow associated to the transport equation is $(0, y)$ where

$$y = x_2 + \int_0^\infty k_2 x_1(s) ds = x_2 + \int_0^{x_1} \frac{k_2 z}{g_1 z} dz = x_2 + \frac{k_2}{g_1} x_1$$

We then consider

$$f(x_1, x_2) = g(x_2 + \frac{k_2}{g_1} x_1),$$

that satisfies, for any $x_1, x_2$,

$$-g_1 x_1 \frac{\partial f}{\partial x_1} + k_2 x_1 \frac{\partial f}{\partial x_2} = 0.$$
Now taking the limit $n \to \infty$ into

$$f(x_1^n(t), x_2^n(t)) - f(x_1^n(0), x_2^n(0)) - \int_0^t \mathcal{A}_2^n f(x_1^n(s), x_2^n(s))ds,$$

yields

$$g(x_2(t)) - g(x_2(0)) - \int_0^t -g_2x_2g'(x_2(s)) + \lambda_1(x_2(s)) \left( \int_0^\infty \bar{h}(z)g(x_2(s) + z)dz - g(x_2(s)) \right) ds,$$

where $\bar{h}(x_2) = (k_2/g_1)^{-1}h((k_2/g_1)^{-1}x_2)$. Hence the limiting process $x_2$ must satisfy the martingale problem associated with the generator

$$\mathcal{A}_{\infty}g(x) = -g_2x \frac{dg}{dx} + \lambda_1(x) \left( \int_x^\infty \bar{h}(z-x)f(z)dz - f(x) \right),$$

for which uniqueness holds for bounded $k_1$ (see (Crudu et al., 2012, thm 2.5)). A truncature argument allows then to conclude for linearly bounded $k_1$.

**Proof of Theorem 4** As noticed before, $(z^n, x_2^n)$ with $z^n(t) = \frac{x_1^n(t)}{n}$ satisfies the scaling (S3) so similar conclusion holds for $x_2^n$.

## References

BALL, K., KURTZ, T. G., POPOVIC, L. and REMPALA, G. (2006). Asymptotic Analysis of Multiscale Approximations to Reaction Networks. *The Annals of Applied Probability* **16** 1925–1961.

BERGLUND, N. and GENTZ, B. (2006). *Noise-Induced Phenomena in Slow-Fast Dynamical Systems, A Sample-Paths Approach*. Springer.

BOSE, I. and GHOSH, S. (2012). Binary Response in MicroRNA-mediated Regulation of Gene Expression. 1–7. arXiv:1205.0381v2.

BUJORIANU, M. L. and KATOEN, J. P. (2008). Symmetry reduction for stochastic hybrid systems. In *Decision and Control, 2008. CDC 2008. 47th IEEE Conference on* DOI - 10.1109/CDC.2008.4739086 233–238.

BUJORIANU, M. L. and LYGEROS, J. (2004). General stochastic hybrid systems: modelling and optimal control. In *Decision and Control, 2004. CDC. 43rd IEEE Conference on* DOI - 10.1109/CDC.2004.1430320 2 1872–1877 Vol.2.

CRUDU, A., DEBUSSCHE, A., MULLER, A. and RADULESCU, O. (2012). Convergence of stochastic gene networks to hybrid piecewise deterministic processes. *Ann. Appl. Probab.* **22** 1822–1859.

DAVIS, M. H. A. (1984). Piecewise-Deterministic Markov Processes: A General Class of Non-Diffusion Stochastic Models. *Journal of the Royal Statistical Society. Series B (Methodological)* **46** 353–388.
ELF, J., LI, G.-W. and XIE, X. S. (2007). Probing Transcription Factor Dynamics at the Single-Molecule Level in a Living Cell. *Science* **316** 1191-1194.

FENICHEL, N. (1979). Geometric singular perturbation theory for ordinary differential equations. *Journal of Differential Equations* **31** 53-98.

FRIEDMAN, N., CAI, L. and XIE, X. S. (2006). Linking Stochastic Dynamics to Population Distribution: An Analytical Framework of Gene Expression. *Phys. Rev. Lett.* **97** 168302–.

GARDINER, C. W. (1985). *Handbook of stochastic methods, for physics, chemistry and the natural sciences*, second ed. Springer Series in Synergetics **13**. Springer-Verlag, Berlin.

GIUSTI, E. (1984). *Minimal surfaces and functions of bounded variations*. Birkhäuser.

GOLDING, I., PAULSSON, J., ZAWILSKI, S. M. and COX, E. C. (2005). Real-Time Kinetics of Gene Activity in Individual Bacteria. *Cell* **123** 1025-1036.

GOUTSIAS, J. (2005). Quasiequilibrium approximation of fast reaction kinetics in stochastic biochemical systems. *J. Chem. Phys.* **122** 184102–15.

HASELTINE, E. L. and RAWLINGS, J. B. (2002). Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. *J. Chem. Phys.* **117** 6959–6969.

HASELTINE, E. L. and RAWLINGS, J. B. (2005). On the origins of approximations for stochastic chemical kinetics. *J. Chem. Phys.* **123** 164115–16.

HESPANHA, J. P. (2006). Modelling and analysis of stochastic hybrid systems. *Control Theory and Applications*, *IEE Proceedings - DOI* - **153** 520–535.

JACOD, J. and SHIRYAEV, A. N. (1987). *Limit theorems for stochastic processes*. Springer.

JAKUBOWSKI, A. (1997). A non-Skorohod topology on the Skorohod space. *Electronic Journal of Probability* **2** 1–21.

KAERN, M., ELSTON, T. C., BLAKE, W. J. and COLLINS, J. J. (2005). Stochasticity in gene expression: From theories to phenotypes. *Nat. Rev. Genet.* **6** 451-464.

KANG, H.-W. and KURTZ, T. G. Separation of time-scales and model reduction for stochastic reaction networks. To appear in *Annals of Applied Prob.*

KURTZ, T. G. (1991). Random time changes and convergence in distribution under the Meyer-Zheng conditions. *The Annals of Probability* **19** 1010–1034.

KURTZ, T. G. (1992). Averaging for martingale problems and stochastic approximation. *Applied Stochastic Analysis Lecture Notes in Control and Information Sciences* **177** 186-209.

LEI, J. (2009). Stochastic in single gene expression with both intrinsic noise and fluctuation in kinetic parameters. *Journal of Theoretical Biology* **256** 485-492.

MACKEY, M. C., TYRAN-KAMIŃSKA, M. and YVINEC, R. (2011). Molecular distributions in gene regulatory dynamics. *Journal of Theoretical Biology* **274** 84-96.

OZBUDAK, E. M., THATTAI, M., KURTSEV, I., GROSSMAN, A. D. and VAN OUDENAARDEN, A. (2002). Regulation of noise in the expression of a single gene. *Nat Genet* **31** 69–73.

RAJ, A. and VAN OUDENAARDEN, A. (2009). Single-Molecule Approaches to
Stochastic Gene Expression. *Annu. Rev. Biophys.* **38** 255–270.

RAJ, A., PESKIN, C. S., TRANCHINA, D., VARGAS, D. Y. and TYAGI, S. (2006). Stochastic mRNA Synthesis in Mammalian Cells. *PLoS Biol* **4** e309–

RAO, C. V. and ARKIN, A. P. (2003). Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. *J. Chem. Phys.* **118** 4999–5010.

RIEDLER, M. G., THIEULLEN, M. and WAINRIB, G. (2012). Limit theorems for infinite-dimensional piecewise deterministic Markov processes. Applications to stochastic excitable membrane models. *Electronic Journal of Probability* **17** 1–48.

SANTILLN, M. and QIAN, H. (2011). Irreversible thermodynamics in multiscale stochastic dynamical systems. *Phys Rev E Stat Nonlin Soft Matter Phys.* **83** 1–8.

SCHWANHÄUSSER, B., BUSSE, D., LI, N., DITTMAR, G., SCHUCHHARDT, J., WOLF, J., CHEN, W. and SELBACH, M. (2011). Global quantification of mammalian gene expression control. *Nature* **473** 337–42.

SHAHREZAEI, V. and SWAIN, P. S. (2008). Analytical distributions for stochastic gene expression. *Proc Natl Acad Sci U S A* **105** 17256–17261.

SUTER, D. M., MOLINA, N., GATFIELD, D., SCHNEIDER, K., SCHIBLER, U. and NAEF, F. (2011). Mammalian Genes Are Transcribed with Widely Different Bursting Kinetics. *Science* **332** 472–474.

TIKHONOV, A. N. (1952). Systems of differential equations containing small parameters in the derivatives. *Mat. Sb. (N.S.)* **31** (73) 575–586.

XIE, X. S., CHOI, P. J., LI, G.-W., LEE, N. K. and LIA, G. (2008). Single-Molecule Approach to Molecular Biology in Living Bacterial Cells. *Annu. Rev. Biophys.* **37** 417–444.

YVINEC, R., MACKEY, M. C., LEI, J. and ZHUGE, C. (2012). Adiabatic reduction of jump processes in a model of stochastic gene expression with bursting transcription. *submitted*. arXiv:1202.5411.

ZERON, E. S. and SANTILLÁN, M. (2010). Distributions for negative-feedback-regulated stochastic gene expression: Dimension reduction and numerical solution of the chemical master equation. *Journal of Theoretical Biology* **264** 377–385.