SubmittedtotheAnnalsofAppliedProbability ADIABATIC REDUCTION OF A PIECEWISE DETERMINISTIC MARKOV MODEL OF STOCHASTIC GENE EXPRESSION WITH BURSTING TRANSCRIPTION 2 By Jinzhi Lei , Michael C. Mackey , Romain 1 Yvinec∗ and Changjing Zhuge 0 2 Tsinghua University, McGill University, Universit´e de Lyon b e This paper considers adiabatic reduction in a model of stochas- F tic gene expression with bursting transcription. We prove that an 4 adiabatic reduction can be performed in a stochastic slow/fast sys- 2 tem with a jump Markov process. In thegene expression model, the production of mRNA(thefast variable) is assumed to follow acom- ] pound Poisson process (the phenomena called bursting in molecular R biology)andtheproductionofprotein(theslowvariable)islinearas P afunctionofmRNA.WhenthedynamicsofmRNAisassumedtobe . faster than the protein dynamics (due to a mRNA degradation rate h t largerthanfortheprotein)weprovethat,withtheappropriatescal- a ing,theburstingphenomenacanbetransmittedtotheslowvariable. m Weshow that thereduced equation is eithera stochastic differential [ equation with a jump Markov process or a deterministic ordinary differential equation depending on the scaling that is appropriate. 1 These results are significant because adiabatic reduction techniques v seem to have not been developed for a stochastic differential system 1 1 containingajumpMarkovprocess.Wehopethattheyaregeneraliz- 4 abletogiveinsightintoadiabaticmethodsinmoregeneralstochastic 5 hybridsystems.Inparticular,asimilarproofgivesusamathematical . justificationoftheburstingofmRNAthroughanappropriatescaling 2 of an “ON-OFF” switching gene model. Last but not least, for our 0 2 particular system, the adiabatic reduction allows us to understand 1 whatarethenecessary conditionsfortheburstingproduction-likeof : mRNAand protein to occur. v i X r 1. Introduction. The adiabatic reduction technique gives results that a allow one to reduce the dimension of a system and justifies the use of an effective set of reducedequations in lieu of dealing with a full,higher dimen- sional, model if different time scales occur in the system. Adiabatic reduc- tion results for deterministic systems of ordinary differential equations have ∗ corresponding author 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 1 imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 2 R.YVINECET AL. been available since thevery precise results of Tikhonov (1952) andFenichel (1979). The simplest results, in the hyperbolic case, give an effective con- struction 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.Suchpreciseandgeometric resultshavebeengeneralized 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 knowl- edge, analogous results for stochastic differential equations with jumps have not been obtained. The present paper gives a theoretical justification of an adiabatic reduc- tion of a particular piecewise deterministic Markov process (Davis, 1984). The results we obtain do not give a bound on the error of the reduced sys- tem, but they do allow us to justify the use of a reduced system in the case of a piecewise deterministic Markov process. In that sense, the results are close to the recent ones by Crudu et al. (2012) and Kang and Kurtz (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 approxima- tiondependingonthedifferentscalingchosen(seeBalletal.(2006)forsome examples in a reaction network model). After the scaling, the limiting mod- els can be deterministic (ordinary differential equation), stochastic (jump Markov process), or hybrid (piecewise deterministic process). For illustra- tive and motivating examples given by a simulation algorithm, see Haseltine and Rawlings (2002); Rao and Arkin (2003); Goutsias (2005). Ourparticular modelis meanttodescribestochastic geneexpressionwith explicit bursting (Friedman, Cai and Xie, 2006). The variables evolve under the action of a continuous deterministic dynamical system interrupted by positive jumps of random sizes that model the burst production. In that sense, the convergence theorems we obtain in this paper can be seen as an example in which there are an infinite number of reactions (one for each “burstsize”) and give complementary results to those of Cruduet al. (2012) and Kang and Kurtz(2012). We hopethat 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 is defined as the produc- tion of several molecules within a very short time, indistinguishable within imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 ADIABATICREDUCTION FOR PDMP 3 the time scale of the measurement. In the biological context of models of stochasticgeneexpression,explicitmodelsofburstingmRNAand/orprotein 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-Kaminska and Yvinec, 2011) as more and more experimental evidence from single-molecule visualization techniques has re- vealed 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, andsometimes theDNAstate). Theuseofareducedone-dimensionalmodel (that has theadvantage that it can besolved 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 butvery different behavior (for an example in that context, see Mackey, Tyran-Kaminska and Yvinec (2011)). Hence, our results are of importance to rigorously prove the validity of using a reduced model. 1.1. Continuous-state burstingmodel. Themodelsreferredtoabovehave explicitly assumed the production of several molecules instantaneously, through a jump Markov process, in agreement with experimental obser- vations. 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. The intensity of the jumps can be a bounded function to allow for self-regulation. The assumption that the intensity is bounded is crucial, and makes the stochastic production term a compound Poisson process. LetX andY denotetheconcentrationsofmRNAandproteinrespectively. A simple model of single gene expression with bursting in transcription is given by dX (1) = −γ X +N˚(h,ϕ(Y)), 1 dt dY (2) = −γ Y +λ X. 2 2 dt Here γ and γ are the degradation rates for the mRNA and protein re- 1 2 spectively, λ is the mRNA translation rate, and N˚(h,ϕ(Y)) describes the 2 transcription that is assumed to be a compound Poisson white noise occur- ring at a rate ϕ with a non-negative jump size ∆X distributed with density imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 4 R.YVINECET AL. h. The equations (1)-(2) are a short hand notation for t t ∞ ∞ (3) Xt = X0− γ1Xs−ds+ 1{r≤ϕ(Ys−)}zN(ds,dz,dr) Z0 Z0 Z0 Z0 t t (4) Yt = Y0− γ2Ys−ds+ λ2Xsds. Z0 Z0 where Xs− = limt→s−X(t), and N(ds,dz,dr) is a Poisson random measure on (0,∞)×[0,∞)2 with intensity dsh(z)dzdr, where s denotes the times of thejumps,r isthestate-dependency inanacceptance/rejection fashion,and z the jump size. Note that (X ) is a stochastic process with almost surely t finite variation on any bounded interval (0,T), so that the last integral is well defined as a Stieltjes-integral. The following discussion is valid for general rate functions ϕ and density functions h(∆X) that satisfy (5) ϕ ∈ C1, ϕ and ϕ′ are bounded, ie ϕ ≤ ϕ,ϕ′ ≤ ϕ ∞ (6) h ∈ C0 and xnh(x)dx < ∞, ∀n≥ 1. Z0 For a general density function h, the average burst size is given by ∞ (7) b = xh(x)dx. Z0 Remark 1. Hill functions are often used to model gene self-regulation, so that ϕ is given by 1+yn ϕ(y) = κ Λ+∆yn whereκ,Λ,∆arepositiveparametersandnisapositiveinteger(seeMackey, Tyran-Kaminska and Yvinec (2011) for more details). An exponential dis- tributionoftheburstingtranscriptionisoften usedinmodelinggeneexpres- sion, in accordance with experimental findings, so that the density function h is given by 1 h(∆X) = e−∆X/b, b with b the average burst size. If ϕ is independent of the state Y, the average transcription rate is λ = 1 imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 ADIABATICREDUCTION FOR PDMP 5 bϕ, and the average mRNA and protein concentrations are bϕ (8) X = . eq γ 1 λ bϕλ 2 2 (9) Y = X = . eq eq γ γ γ 2 1 2 1.2. Statement of the results. Inthefollowingdiscussion,weconsiderthe situation whenmRNAdegradation is afast process,i.e.γ is ‘large enough’, 1 but the average protein concentration Y remains unchanged. Thus, we eq always assume one of the following three scaling relations: (S1) ϕ/γ is independentof γ (= Y γ2 if ϕ is a constant) as γ → ∞ and 1 1 eqbλ2 1 h,λ remain unchanged; or 2 (S2) h → 1 h(∆X) (so that b = Y γ2γ1 ) as γ → ∞ and ϕ,λ remain γ1 γ1 eqϕ(Y)λ2 1 2 unchanged; or (S3) λ /γ is independent of γ (= Y γ2γ1 if ϕ is a constant) as γ → ∞ 2 1 1 eq bϕ 1 and h, ϕ remain unchanged. In this paper we determine an effective reduced equation for equation (2) for each of the three scaling conditions (S1)-(S3). In particular, we show that under assumption (S1), equation (2) can be approximated by the de- terministic ordinary differential equation dY (10) = −γ Y +λ ψ(Y) 2 2 dt where (11) ψ(Y) = bϕ(Y)/γ . 1 We further show that under the scaling relations (S2) or (S3), equation (2) can be reduced to the stochastic differential equation dY (12) = −γ Y +N˚(h¯(∆Y),ϕ(Y)). 2 dt where h¯ is a suitable density function in the jump size ∆Y (to be detailed below). We first explain, using some heuristic arguments, the differences between the three scaling relations and the associated results. When γ → ∞, apply- 1 ing a standard quasi-equilibrium assumption we have dX ≈ 0 dt imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 6 R.YVINECET AL. which yields 1 (13) X(t) ≈ N˚(h(.),ϕ(y)) = N˚(γ h(γ ·),ϕ(y)), 1 1 γ 1 and therefore the second equation (2) becomes dY λ ≈ −γ Y + 2N˚(h(.),ϕ(y)) 2 dt γ 1 γ γ ≈ −γ Y +N˚ 1h( 1),ϕ(y) . 2 λ λ (cid:18) 2 2 (cid:19) Hence in (12), h¯(∆Y) = (λ /γ )−1h((λ /γ )−1∆Y) under the scaling (S2) 2 1 2 1 and (S3). Furthermore, wenote that the scaling (S2) also implies γ h(γ ·)= 1 1 h(·), and therefore we can also write h¯(∆Y) = (λ )−1h((λ )−1∆Y) under 2 2 the scaling (S2). Finally, it is intuitively understandable that scaling (S1) gives (10), as the jumps become more frequent and smaller. In this paper, we provide three different proofs of the results mentioned above. In particular, we prove the results using a Master equation approach (the Kolmogorov forward equation) as well as starting from the stochastic differential equation. Note that both techniques have been used in the past, in particular within the context of discrete models of stochastic reaction networks. For the Master equation approach, see Haseltine and Rawlings (2005); Zeron and Santill´an (2010) while for the stochastic differential equa- tion approach, we refer to Crudu et al. (2012); Kang and Kurtz (2012). 1.3. Outline of the paper. This paper is organized as follows. In Section 2, we consider the situation without auto-regulation so the rate ϕ is inde- pendent of protein concentration Y. In this case the two equations (1)-(2) form a set of linear stochastic differential equations. We use the method of characteristic functionals to give a very short proof of the adiabatic reduc- tion when the white noise process is independent of the stochastic variables (external regulation). In Section 3 we consider the situation in which the rate ϕ is dependent on the protein concentration Y. We first give a result on the evolution equation on densities (Section 3.2), then we prove rigor- ously a weak convergence of the stochastic process using the machinery of generators (Section 3.3). Finally, in Section 4, we apply the same formalism with generators to prove that a binary “ON-OFF“ model converges to a bursting model, thus illustrating that these techniques can be applicable in other contexts. imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 ADIABATICREDUCTION FOR PDMP 7 2. The case without auto-regulation. First, we consider the simple case with no auto-regulation in the mRNA transcription so the rate ϕ is independent of protein concentration Y. In this case, equation (1) is inde- pendent of (2) which makes the analysis easier. Consider the equations dX (14) = −γ X +N˚(h,ϕ), X(0) = X ≥ 0, 1 0 dt dY (15) = −γ Y +λ X, Y(0) = Y ≥ 0, 2 2 0 dt where N˚(h,ϕ) is a compound Poisson white noise, and ϕ is a constant. The solutions X(t) and Y(t) of (14)-(15) are stochastic processes uniquely determined by the equation parameters and the stochastic process N˚. We first state main results for each of the scaling relations we have con- sidered. Remark 2. In all scaling situations that we consider in the three the- orems below, the stochastic process X converges only in finite-dimensional t law. Here finite-dimensional law means that for any finite number of fixed times (t ,t ,···t ) the random variable (X(t ),X(t ),···X(t )) converges 1 2 n 1 2 n inlawtowardstheappropriatedistribution.Convergenceinfinite-dimensional law is weaker than convergence in law of the stochastic process. We will see in the proofs of the theorems that this is because a lack of compactness for the process X(t). Theorem 1. Consider the equations (14)-(15). If the scaling (S1) is γ γ 2 1 satisfied, i.e., ϕ = Y , then when γ → ∞, eq 1 bλ 2 1. The stochastic process X(t) converges in finite-dimensional law to the (deterministic) steady-state value X ; eq 2. The stochastic process Y(t) converges in law towards the deterministic solution of the ordinary differential equation dY (16) = −γ Y +λ X , Y(0) = Y ≥ 0. 2 2 eq 0 dt Theorem 2. Consider the equations (14)-(15). If the scaling (S2) is γ γ satisfied, i.e., h → 1 h(∆X) (so b = Y 2 1), then when γ → ∞, γ1 γ1 eqϕλ 1 2 1. The stochastic process X(t) converges in finite-dimensional law to the compound Poisson white noise N˚(h,ϕ); imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 8 R.YVINECET AL. 2. The stochastic process Y(t) converges in law to the stochastic process defined by the solution of the stochastic differential equation dY (17) = −γ Y +N˚(h¯,ϕ), Y(0) = Y ≥ 0 2 0 dt where h¯(∆Y) = (λ )−1h((λ )−1∆Y) is the density function for the 2 2 jump size ∆Y. Theorem 3. Consider the equations (14)-(15). If the scaling (S3) is γ γ 2 1 satisfied, i.e., λ = Y , then when γ → ∞, 2 eq 1 bϕ 1. The stochastic process X(t) converges in finite-dimensional law to the (deterministic) fixed value 0; 2. The stochastic process Y(t) converges in law to the stochastic process determined by the solution of the stochastic differential equation dY (18) = −γ Y +N˚(h¯,ϕ), Y(0) = Y ≥ 0 2 0 dt where h¯(∆Y) = (λ /γ )−1h((λ /γ )−1∆Y) is the density function for 2 1 2 1 the jump size ∆Y. Remark 3. Note that scalings (S2) and (S3) give similar results for the equation governing the protein variable Y but very different results t for the asymptotic stochastic process related to the mRNA. In particular, in Theorem 2, very large bursts of mRNA are transmitted to the protein, where in Theorem 3, very rarely is mRNA present but when present it is efficiently synthesized into a burst of protein. These theorems will be proved using the method of characteristic func- tionals. For a stochastic process ξ (t ≥ 0), the characteristic functional t C : Σ → R is defined as ξ ∞ (19) C [f]= E eR0 if(t)ξtdt ξ h i for any function f in a suitable function space Σ so that the integral ∞ if(t)ξ dt is well defined. Before continuing, we need to introduce some 0 t topological background as well as properties of the Fourier transform in R nuclear spaces (see Gel’fand and Vilenkin (1964)) imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 ADIABATICREDUCTION FOR PDMP 9 2.1. Stochastic process as a distribution. We are going to recall here the continuous correspondence between a stochastic process and a distribution. All stochastic processes we consider here are right continuous with finiteleft limits. It is naturalto definea stochastic processas arandomelement in the space D = D (0,∞),R , the space of all real-valued functions on (0,∞) that are right(cid:16)continuou(cid:17)s with finite left limits. Such a space is classically endowed with the Skorokhod topology (see for instance Jacod and Shiryaev (1987, Chap. 6)). We define D(R+), the space of smooth functions with compact support, with the inductive limit topology given by the family of semi-norms (k = 0,1,2...) p (f) = sup|f(k)| on every D([0,n]), n ∈ N (c.f. k (Schaefer, 1971, Example 2, page 57)). Let f ∈ D(R+), and define x˜ in the dual space D′(R+) such that ∞ (20) x˜(f) = x(t)f(t)dt Z0 Lemma 4. The map D (0,∞),R → D′(R+) (21) (cid:16) (cid:17) (x ) 7→ x˜, t t≥0 where x˜ is defined by equation (20), is continuous. Proof. Itisaclassical resultthatx ∈ D hasatmostacountablenumber of discontinuity points so that x is locally integrable, the integral in Eq (20) is well defined for all f ∈ D(R+), x˜ ∈ D′(R+) and T |x˜(f)|≤ | x(s)| ds k f k ∞ (cid:16)Z0 (cid:17) for any f with support in [0,T] (Rudin, 1991, Section 6.11, page 142). We conclude by noticing that | x˜(f)|≤ sup |x(s) |k f k T ∞ s≤T and x 7→ sup | x(s) | is continuous for the Skorohod topology (Jacod s≤T and Shiryaev, 1987, Proposition 2.4, page 339) for all T such that T is not a discontinuity point. imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012 10 R.YVINECET AL. 2.2. Bochner-Minlos theorem for a nuclear space. Let E be a nuclear space. Westate akey resultthatwillallow usto uniquelyidentify ameasure on the dual E′ of E. Bochner-Minlos Theorem. (Gel’fand and Vilenkin, 1964, Theorem 2, page 146) For a continuous functional C on a nuclear space E that sat- isfies C(0) = 1, and for any complex z and elements x ∈A, j,k = 1,...,n, j j n n (22) z z¯ C(x −x )≥ 0, j k j k j=1k=1 XX there is a unique probability measure µ on the dual space E′, given by (23) C(y) = eihx,yidµ(x). ZE′ Note that the space D(R+) is a nuclear space (Schaefer, 1971, Example 2, page 107). 2.3. The characteristic functional of a Poisson white noise. The use of the characteristic functional allows us to definea generalized stochastic pro- cess that does not necessarily have a trajectory in the usual sense (like in D for instance). Indeed a (compound) Poisson white noise is seen as a ran- dom measure on the distribution space D, associated with the characteristic functional (given in Hida and Si (2004), here f ∈ D(R+)) ∞ ∞ (24) C [f]= exp ϕ (eizf(t) −1)h(z)dzdt , N˚ (cid:20) Z0 Z0 (cid:21) where ϕ is the Poisson intensity and h the jump size distribution. It is not hard to see that C [f − g] and C [g − f] are conjugate to each other, N˚ N˚ C [0] = 1 and C [.] is continuous for h ∈ L1(R+), so the conditions in the N˚ N˚ Bochner-Minlos Theorem2.2aresatisfiedandthereforeC uniquelydefines N˚ a measure on D′(R+). We refer to Prokhorov (1956); Hida and Si (2004, 2008) for further ma- terial on characteristic functionals and generalized stochastic processes. 2.4. Proofs of Theorems 1 through 3. The proofs of Theorems 1 to 3 are based on the idea of Levy’s continuity theorem. However in the infinite- dimensional case, the convergence of the Fourier transform does not imply convergence inlawoftherandomvariable,andoneneedstoimposemorere- strictions, namelyacompactness condition.Wewillusethefollowinglemma imsart-aap ver. 2011/11/15 file: adiabatic_aap_finalpreprint_24feb12.tex date: February 27, 2012