Analytic solutions for stochastic hybrid models of gene regulatory networks
AANALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS OF GENEREGULATORY NETWORKS
PAVEL KURASOV, DELIO MUGNOLO, AND VERENA WOLF
Abstract.
Discrete-state stochastic models are a popular approach to describe the inherentstochasticity of gene expression in single cells. The analysis of such models is hindered by thefact that the underlying discrete state space is extremely large. Therefore hybrid models, inwhich protein counts are replaced by average protein concentrations, have become a popularalternative.The evolution of the corresponding probability density functions is given by a coupled systemof hyperbolic PDEs. This system has Markovian nature but its hyperbolic structure makes itdifficult to apply standard functional analytical methods. We are able to prove convergencetowards the stationary solution and determine such equilibrium explicitly by combining abstractmethods from the theory of positive operators and elementary ideas from potential analysis. Introduction
Very small copy numbers of genes, RNA, or protein molecules in single cells give rise toheterogeneity across genetically identical cells [RvO08]. During the last decades discrete-statestochastic models have become a popular approach to describe gene expression in single cellssince they adequately account for discrete random events underlying such cellular processes (see[SSG17] for a review). Exact solutions for such models are available if they obey detailed balance[Lau00] or restrict to monomolecular intracellular interactions [JH07]. However, since generegulatory networks typically contain feedback loops, second-order interactions are necessary foran adequate description. Moreover, neither detailed balance nor linear dynamics are realisticassumptions even for simple regulatory networks. Recently, analytical solutions for single-genefeedback loops have been presented [GSN12, HSI +
05, KPK14, LYWZ16, VB13, VAE08].The underlying state space of discrete-stochastic models that describe gene regulatory net-works is typically extremely large due to the combinatorial nature of molecule counts for differ-ent types of chemical species. Moreover, realistic upper bounds on protein counts are often notknown. Therefore, hybrid models have become popular in which for highly-abundant speciesonly average counts are tracked while discrete random variables are used to represent specieswith low copy numbers. These hybrid approaches allow for faster and yet accurate Monte Carlosampling that stochastically selects counts of species with low copy numbers and numericallyintegrates average counts of all other species [HGK15, CDR09, HH12, PK04].Hybrid or fluid approaches have also been investigated in the context of stochastic Petri nets asa kind of mean field approximation and gave rise to fluid stochastic Petri nets [TK93, HKNT98].
Mathematics Subject Classification.
Key words and phrases.
Petri networks, systems of PDEs, positive C -semigroups. a r X i v : . [ m a t h . A P ] D ec PAVEL KURASOV, DELIO MUGNOLO, AND VERENA WOLF
L¨uck and the present authors have shown in [KLMW18] how this formalism can be successfullyadapted to study gene regulatory networks. More precisely, a stochastic hybrid approach for generegulatory networks has been proposed therein, in which the state of the genes is representedby a discrete-stochastic variable, while the evolution of the protein numbers is modeled by anordinary differential equation. The main aim of [KLMW18] was to show that fluid stochasticPetri nets allow for more efficient numerical simulations than common, purely discrete masterequations, while providing solutions that are as accurate as those provided in reference studies,like [GSN12]. The scope of the present paper is more analytic: in particular, we discuss indetails properties of the evolution equation – in fact, a coupled system of hyperbolic PDEs –that is at the core of the theory of fluid stochastic Petri nets.The model based on a one-fluid stochastic Petri net is mathematically described by the systemof first-order partial differential equations(1.1) ∂∂t (cid:18) f ( x, t ) f ( x, t ) (cid:19) = − ∂∂x (cid:20)(cid:18) a − bx c − dx (cid:19) (cid:18) f ( x, t ) f ( x, t ) (cid:19)(cid:21) + (cid:18) − λ ( x ) µ ( x ) λ ( x ) − µ ( x ) (cid:19) (cid:18) f ( x, t ) f ( x, t ) (cid:19) . The basic setting is outlined in Section 2.This system has to be complemented by suitable boundary and initial conditions: we willshow that the requirement of the solution to be a probability distribution – and in particular apositive-valued L -function – determines a specific choice of boundary conditions. Furthermore,we are interested in calculating analytically the equilibrium distribution. These distributionscan be computed numerically in a few special cases but in order to understand their propertiesit is useful to have explicit formulas for the solutions: this is done in Section 3.The analytic features of this differential equations are subtler than one may believe at afirst glance: the innocent looking coupling turns a system of hyperbolic equations (which canbe explicitly solved and whose solutions would otherwise become extinct in finite time) into atime-continuous Markov chain, which we are going to study by the classical theory of strongly-continuous semigroups of operators on an appropriate L -space. Proving well-posedness of theassociated Cauchy problem is straightforward. Yet, determining the long-time behavior of this C -semigroup turns out to be more elusive than one could naively conjecture: one expects thesystem to converge to its steady-state solution, and this guess is indeed correct. More precisely,it is a natural goal to prove that the spectrum of the differential operator matrix acting on avector-valued L -space that appears in (1.1) has no eigenvalue on the imaginary axis other than0, thus excluding oscillatory behavior. However, the most classical techniques we have tried toapply to achieve this task fall short off the mark: in Section 4 we develop a method that seemsto be new in the literature and of independent interest. It is based on a combination of classicalPerron–Frobenius theory for positive semigroups and some recent compactness results on kerneloperators on L -spaces.In Section 5 we finally compare our abstract results with numerical simulations, finding thatthey are in good agreement. NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 3 Formulation of the model
On the domain(2.1) x ∈ (cid:104) ab , cd (cid:105) , t ≥ , consider the evolution problem determined by equation (1.1). To ensure that the interval [ ab , cd ]is non-degenerate, we assume that(2.2) ad < bc and that b (cid:54) = 0 , d (cid:54) = 0. (In equations deriving from biological models it is natural to interpret a, b, c, d as strictly positive rates. However, this is not relevant for the purpose of our analysis,and in fact, in Section 5 we are going to study some toy models involving negative rates, too.)The functions λ, µ are assumed to be continuous and positive, i.e.,(2.3) λ, µ ∈ C (cid:16)(cid:104) ab , cd (cid:105)(cid:17) ,λ ( x ) > µ ( x ) > x ∈ (cid:104) ab , cd (cid:105) . Since λ, µ are positive on a compact interval, it follows that their minimum is larger than 0, too.This assumption on λ, µ may be relaxed, see Remark 4.4, but we impose it in order to simplifyour presentation.We assume that initial conditions at t = 0 are provided(2.4) f j ( x,
0) = f j ( x ) , j = 1 , . We are going to show that the evolution equation is governed by a positivity preserving C -semigroup, so that positive initial data determine positive solutions ( f ( · , t ) , f ( · , t )) for any t > . Moreover, it is straightforward to prove that if the positive vector-valued function ( f , f )solves (1.1), then (cid:90) c/da/b ( f ( x, t ) + f ( x, t )) dx, t ≥ , is an integral of motion, i.e., it is preserved over time under the evolution of (1.1): this isequivalent to saying that the semigroup is stochastic and can be proven by integration by parts,see Proposition 4.3 below. Hence the function ( f , f ) can be interpreted as (vector-valued)density assuming the normalization condition(2.5) (cid:90) c/da/b ( f ( x, t ) + f ( x, t )) dx = 1 for all t ≥ . Therefore it is natural to look for solutions in the space L ( ab , cd ) . Then the maximal domain for the semigroup generator A given by(2.6) A f ( x ) := − ∂∂x (cid:20)(cid:18) a − bx c − dx (cid:19) (cid:18) f ( x, t ) f ( x, t ) (cid:19)(cid:21) + (cid:18) − λ ( x ) µ ( x ) λ ( x ) − µ ( x ) (cid:19) (cid:18) f ( x, t ) f ( x, t ) (cid:19) PAVEL KURASOV, DELIO MUGNOLO, AND VERENA WOLF consists of functions from L ( ab , cd ) such that ∂∂x ( a − b · ) f , ∂∂x ( c − d · ) f ∈ L (cid:16) ab , cd (cid:17) . In particular this implies that if ( f , f ) is a solution to (1.1) for initial data f , f ∈ W , (cid:16)(cid:104) ab , cd (cid:105)(cid:17) , then both f and f belong to the Sobolev space W , on any compact sub-interval, moreprecisely(2.7) f ∈ W , (cid:16) ab + (cid:15), cd (cid:17) , f ∈ W , (cid:16) ab , cd − (cid:15) (cid:17) , ∀ (cid:15) > . Observe that for all x ∈ ( ab , cd ) the coefficient matrix(2.8) (cid:18) bx − a dx − c (cid:19) is indefinite, since bx − a > dx − c < ab , cd ). Indeed, (1.1) isa hyperbolic system: Because for all x in the relevant interval bx − a > dx − c < f , and on the left endpoint on the interval for f . In the following, we are going toshow that this is, in fact, necessarily the case. We are going to obtain these boundary conditionsby investigating the possible stationary state, see Theorem 3.2.3. Search for the stationary state
Numerical simulations of the system (1.1) show that it always tends to equilibrium, there-fore it looks natural to start our analysis from investigating a possible stationary solution andcalculating it explicitly. A stationary solution should satisfy the system of ordinary differentialequations:(3.1) − ddx (cid:20)(cid:18) a − bx c − dx (cid:19) (cid:18) ψ ( x ) ψ ( x ) (cid:19)(cid:21) + (cid:18) − λ ( x ) µ ( x ) λ ( x ) − µ ( x ) (cid:19) (cid:18) ψ ( x ) ψ ( x ) (cid:19) = 0 . Under our assumptions any stationary solution has to be a continuously differentiable functioninside the interval ( ab , cd ) but may have singularities at the boundary points. Nonetheless thesingularities cannot be too strong, since the densities have to be integrable functions. Let usprove the following elementary statement to be used in what follows to deduce certain necessaryboundary conditions. Lemma 3.1.
Let ψ be a positive continuous integrable function defined on the interval (0 , .Then there exists a sequence ( x n ) n ∈ N ⊂ (0 , such that (3.2) lim n →∞ x n = 0 and lim n →∞ x n ψ ( x n ) = 0 . NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 5
Proof.
Let us present such a sequence explicitly. Consider the intervals [ n , n ] and denote by x n ∈ [ n , n ] one of the minimum points for ψ in the interval: so ψ ( x ) ≥ ψ ( x n ) whenever x ∈ [ 12 n , n ]If x n ψ ( x n ) tends to zero, then we have such a sequence. Let us now assume that there is nosubsequence of ( x n ψ ( x n )) n ∈ N tending to zero: then there is a positive number δ > x n ψ ( x n ) > δ for all sufficiently large n . But then it follows that ψ satisfies the lower estimate ψ ( x ) ≥ δ n whenever x ∈ [ 12 n , n ]and hence cannot be integrable since (cid:90) f ( x ) dx > ∞ (cid:88) n =1 δ n n +1 = ∞ . This contradiction proves the lemma. (cid:3)
Lemma 3.1 applied to ψ and ψ implies that there exist sequences x − n → ab , x + n → cd suchthat the following two limits hold(3.3) lim n →∞ ( a − bx − n ) ψ ( x − n ) = 0 , lim n →∞ ( c − dx + n ) ψ ( x + n ) = 0 . Let us now return back to the differential system (3.1) and sum up the two equations to get(3.4) ddx [( a − bx ) ψ ( x ) + ( c − dx ) ψ ( x )] = 0 , which implies that(3.5) ( a − bx ) ψ ( x ) + ( c − dx ) ψ ( x ) ≡ K for some K ∈ R . Let us determine this constant by evaluating the function near x = ab and x = cd using the sequences x ± n introduced above: K = lim n →∞ (cid:0) ( a − bx − n ) ψ ( x − n ) + ( c − dx − n ) ψ ( x − n ) (cid:1) = 0 + lim n →∞ (cid:0) ( c − dx − n ) ψ ( x − n ) (cid:1) ≥ K = lim n →∞ (cid:0) ( a − bx + n ) ψ ( x + n ) + ( c − dx + n ) ψ ( x + n ) (cid:1) = lim n →∞ (cid:0) ( a − bx + n ) ψ ( x + n ) (cid:1) + 0 ≤ . PAVEL KURASOV, DELIO MUGNOLO, AND VERENA WOLF
It follows that K = 0 and therefore the functions ψ , (continuous inside the open interval) tendto zero limits at the right and left boundary points respectively(3.6) ψ (cid:16) cd (cid:17) = ψ (cid:16) ab (cid:17) = 0 . Let summarize our studies.
Theorem 3.2.
Every positive integrable solution ( ψ , ψ ) to the system (3.1) is continuous onthe open interval ( ab , cd ) and satisfies the following conditions at the boundary points: (3.7) lim x → a/b ( a − bx ) ψ ( x ) = 0 , lim x → c/d ψ ( x ) = 0 , lim x → a/b ψ ( x ) = 0 , lim x → c/d ( c − dx ) ψ ( x ) = 0 . Proof.
To accomplish the proof we just need to show that the limits lim x → a/b ( a − bx ) ψ ( x ) = 0and lim x → c/d ( c − dx ) ψ ( x ) = 0 hold, not just along subsequences as in Lemma 3.1. This followsdirectly from (3.5) and the fact that K = 0( a − bx ) ψ ( x ) + ( c − dx ) ψ ( x ) = 0 . Taking into that ψ ( x ) tends to zero at x = c/d and ψ ( x ) - at x = a/b we arrive at (3.7). (cid:3) We observe that it is not convenient to keep working with two density functions, due to anexplicit relation between them. Let us namely consider a solution ( ψ , ψ ) to the system (3.1)and introduce a new positive continuously differentiable function(3.8) h ( x ) := ( bx − a ) ψ ( x ) = ( c − dx ) ψ ( x ) . Hence the function h is continuous in the closed interval [ ab , cd ] and satisfies Dirichlet conditionsat both endpoints:(3.9) h (cid:16) ab (cid:17) = 0 = h (cid:16) cd (cid:17) . Taking the difference between the two equations in (3.1) we get the following single differentialequation on the function h introduced above:(3.10) ddx h ( x ) = (cid:18) λ ( x ) bx − a − µ ( x ) c − dx (cid:19) h ( x ) . We can solve this equation analytically whenever we can integrate the function in brackets onthe right hand side(3.11) h ( x ) = K exp (cid:26)(cid:90) xx (cid:18) λ ( y ) by − a − µ ( y ) c − dy (cid:19) dy (cid:27) , where x ∈ ( a/b, c/d ) - arbitrary and K ∈ (0 , ∞ ) is an arbitrary positive constant. The functions µ and λ are positive definite on a compact interval and therefore µ ( x ) > C > , λ ( x ) > C > C > −∞ at both endpointsof the interval. To see this, let us split the integral as (cid:82) xx (cid:16) λ ( y ) by − a dy − (cid:82) xx µ ( y ) c − dy (cid:17) dy. The second
NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 7 integral is bounded near x = a/b , while the integrand in the first integral can be estimated as λ ( y ) by − a ≥ Cby − a . Hence the following limit holdslim x → a/b (cid:90) xx λ ( y ) by − a dy = − lim x → a/b (cid:90) x x λ ( y ) by − a dy < − lim x → a/b (cid:90) x x Cby − a dy = −∞⇒ lim x → a/b (cid:90) xx λ ( y ) by − a dy = −∞ . It follows that h given by (3.11) satisfies Dirichlet condition at x = a/b. Similarly, it satisfiesDirichlet condition at the opposite endpoint as well.Summing up, we have obtained that the densities of the stationary solution of (1.1) can beexplicitly computed as follows.
Theorem 3.3.
Assume that the positive continuous functions µ and λ on the interval [ ab , cd ] aregiven. Then up to the parameter K there is a unique strongly positive solution to the stationaryequation (3.1) : it is given by ψ ( x ) = Kbx − a exp (cid:26)(cid:90) x (cid:18) λ ( y ) by − a − µ ( y ) c − dy (cid:19) dy (cid:27) ,ψ ( x ) = Kc − dx exp (cid:26)(cid:90) x (cid:18) λ ( y ) by − a − µ ( y ) c − dy (cid:19) dy (cid:27) , x ∈ (cid:16) ab , cd (cid:17) . (3.12) The solution always satisfies boundary conditions (3.7) . Note that despite h satisfies Dirichlet conditions at both endpoints, the densities ψ , mayhave singularities there (see Section 5 for illuminating examples).4. Analysis of the time-dependent system
Let us now finally turn to the study of the system of time-dependent partial differentialequations (1.1). Each of these two equations models the time evolution of a two-state continuous-time Markov chain, the vector-valued function f thus representing a probability distribution.This equation is meaningful for all t ∈ R , but we will study in particular the evolution for t → ∞ in the dependence from the configuration of the system for t = 0.As in the rest of the article, throughout this section we are still imposing the conditions statedin Section 2.4.1. Preliminaries.
The partial differential equation (1.1) describes transport phenomena: In-deed, (probability) mass is shifted to the left or to the right depending on whether the coefficientof the first derivative is positive or negative, respectively.On the other hand, (1.1) arises as a stochastic model and in the biologically relevant case of λ > µ > ddt (cid:18) f ( x, t ) f ( x, t ) (cid:19) = A (cid:18) f f (cid:19) ( x, t ) := (cid:18) − λ ( x ) µ ( x ) λ ( x ) − µ ( x ) (cid:19) (cid:18) f ( x, t ) f ( x, t ) (cid:19) A function is called strongly positive if it is positive outside a set of measure zero [RS78].
PAVEL KURASOV, DELIO MUGNOLO, AND VERENA WOLF that competes with the transport (with space-dependent speed) given by the vector-valuedpartial differential equation ∂∂t (cid:18) f ( x, t ) f ( x, t ) (cid:19) = A (cid:18) f f (cid:19) ( x, t ) := (cid:18) bx − a dx − c (cid:19) ∂∂x (cid:18) f ( x, t ) f ( x, t ) (cid:19) + (cid:18) b d (cid:19) (cid:18) f ( x, t ) f ( x, t ) (cid:19) . (4.2) Remark 4.1.
In general the two operator matrices A , A do not commute, but by [EN00,Exer. III.5.11] the solution f to the complete time-dependent problem (1.1) can be given bymeans of the Lie–Trotter product formula (4.3) f ( t ) = e t A f = lim n →∞ (cid:16) e tn A e tn A (cid:17) n f , where (cid:0) e t A (cid:1) t ≥ and (cid:0) e t A (cid:1) t ≥ are the C -semigroups generated by the operators A , A and ofcourse f is the initial data. The operator splitting in (4.3) seems to be numerically interesting,since the linear dynamical system (4.1) is solved by (4.4) e t A : f (cid:55)→ λ e − t ( µ + λ ) + µµ + λ µ ( − e − t ( µ + λ ) ) µ + λλ ( − e − t ( µ + λ ) ) µ + λ µ e − t ( µ + λ ) + λµ + λ (cid:18) f f (cid:19) , t ≥ , where f , f : ( ab , cd ) → R are the two components of the initial data f . Also the solutions to thepure transport equation (4.2) can be determined, at least in principle: in the case b > , d > that is relevant for our model one solution of (4.2) is given by (cid:18) f f (cid:19) : ( x, t ) (cid:55)→ (cid:32) e tb f (cid:0)(cid:0) x − ab (cid:1) e tb + ab (cid:1) e td f (cid:0)(cid:0) x − cd (cid:1) e td + cd (cid:1)(cid:33) . (4.5) While this formula does not bother to take into account possible boundary conditions, and indeedit is not defined for all ( x, t ) , it gives a hint about the qualitative behavior of the solution: namely,the profile of the initial data is rescaled and concentrated as the time passes by and the spacialargument is deformed according to the laws x (cid:55)→ (cid:16) x − ab (cid:17) e tb + ab , x (cid:55)→ (cid:16) x − cd (cid:17) e td + cd , t ≥ . (Observe that the intervals ( ab , + ∞ ) and ( −∞ , cd ) are left invariant under these monotone trans-formations for all t ≥ . We also note that the speed of propagation approaches 0 as x → ab or x → cd , respectively.)Accordingly, variance diminishes as the total probability is conserved but its profile gets shiftedto the left in the first and to the right in the second equation of (4.2) , respectively; this is insharp contrast with the case of b = 0 = d for the characteristics of the hyperbolic systems arestraight lines and hence mass is steadily leaving the system. As already mentioned, the system of partial differential equations in (1.1) describes the lim-iting case of a Markov chain. For this reason, one expects convergence to an equilibrium givenby the solution to the stationary equation (3.1) which we have investigated in Section 3. But
NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 9 to begin with, let us first analyze the time-dependent problem. Our analysis will be based onproperties of the Banach spaces L ( ab , cd ; R ) and L ∞ ( ab , cd ; R ), normed as (cid:107) f (cid:107) := (cid:90) cdab ( | f ( x ) | + | f ( x ) | ) dx and (cid:107) f (cid:107) ∞ := ess sup ab ≤ x ≤ cd max {| f ( x ) | , | f ( x ) |} ;they are, in fact, Banach lattices whenever endowed with the natural pointwise ordering inducedby the ordering in R . This will be important, as our analysis on the long-time asymptoticsof (1.1) will be based on properties of positivity preserving operator semigroups.It is natural to perform the analysis of the evolution equation (1.1) in an L -space, since weare interested in a stochastic model and the unknown ( f , f ) represents a probability density.Recall that a semigroup on an L -space is called irreducible if its generator A satisfies0 (cid:54)≡ g ≥ λf − Af = g ⇒ f > λ > s ( A ). Lemma 4.2.
Let α (cid:54) = β be real numbers, α < β . Let p ∈ W , ( α, β ) such that(i) either p ( s ) > for all s ∈ ( α, β ) and p ( α ) = 0 (ii) or p ( s ) < for all s ∈ ( α, β ) with p ( β ) = 0 and such that p ∈ L ( α, β ) . Then the operator A : f (cid:55)→ ddx ( pf ) with domain • either D ( A ) := { f ∈ L ( α, β ) : ( pf ) (cid:48) ∈ L ( α, β ) and f ( β ) = 0 } if ( i ) holds, • or D ( A ) := { f ∈ L ( α, β ) : ( pf ) (cid:48) ∈ L ( α, β ) and f ( α ) = 0 } if ( ii ) holds,generates an irreducible semigroup of positivity preserving contractions on L ( α, β ) .The embedding of D ( A ) in L ( α, β ) is not compact, hence A does not have compact resolvent. Observe that A is well-defined, since if pf ∈ W , ( α, β ) and hence pf is continuous on [ α, β ],so in particular its boundary values can be considered. Proof.
First of all, it is clear that the operator A is closed and densely defined. In view of [Nag86, § C-II.1] and [EN00, Cor. II.3.17], it suffices to show that A is dispersive , i.e., (cid:104) Af, φ (cid:105) ≡ (cid:90) βα Af ( x ) φ ( x ) dµ ≤ f ∈ D ( A ) and φ = { f ≥ } , and that A is m -dissipative, i.e., • for all g ∈ L ( ab , cd ) there exists a solution f ∈ D ( A ) of(4.6) (1 + p (cid:48) ) f + p dfdx = g • and additionally (cid:104) Af, ϕ (cid:105) ≡ (cid:90) βα Af ( x ) ϕ ( x ) dµ = 0 for all f ∈ D ( A ) and ϕ = sgn( f ) . Indeed, the solution to (4.6) can be found by the variation of constants formula and is given by f ( x ) = e − (cid:82) x ∗ p (cid:48) p ds (cid:90) x ∗ e (cid:82) s ∗ p (cid:48) p dr g ( s ) p ( s ) ds, x ∈ ( α, β )where ∗ = β or ∗ = α depending on whether assumption ( i ) or ( ii ) is satisfied. This solutionsatisfies the prescribed boundary conditions and, in fact, pf ∈ W , ( α, β ), so f ∈ D ( A ). Thisexplicit formula also shows that f > g (cid:54)≡ A generates an irreduciblesemigroup.Furthermore, there holds(4.7) (cid:90) βα ddx ( p ( x ) f ( x )) sgn( f ( x )) dx = (cid:104) ( p | f | ) ( x ) (cid:105) x = βx = α = 0as well as (cid:90) βα ddx ( p ( x ) f ( x )) { f ≥ } dx = (cid:90) { f ≥ } ddx ( p ( x ) f ( x )) dx = (cid:104) ( pf ) + ( x ) (cid:105) x = βx = α = 0 , respectively, since for a W , -function g one has(4.8) ( g + ) (cid:48) = g (cid:48) { g ≥ } , cf. [GT01, Lemma 7.6].To conclude, let us show that the embedding of D ( A ) in L ( α, β ) is not compact: we provethis assertion only for the case (i), the case (ii) being completely analogous. Pick a real sequence( β n ) n ∈ N ⊂ [ α, β ] with lim n →∞ β n → β and let f n ( x ) := (1 − β n ) x − β n , n ∈ N , x ∈ ( α, β ) . Then f n ∈ D ( A ) for all n ∈ N and indeed (cid:107) f n (cid:107) = 1 and also ( (cid:107) Af n (cid:107) + (cid:107) f n (cid:107) ) n ∈ N is bounded.Therefore, if the embedding of D ( A ) into L ( α, β ) was compact, then ( f n ) n ∈ N would have aconvergent subsequence, say ( f n k ) k ∈ N ; let us denote its limit by f . Since ( f n ) n ∈ N has anothersubsequence that converges to f almost everywhere. But lim n →∞ f n ( x ) = 0 for a.e. x ∈ ( α, β ),so f = 0 as well, a contradiction to the fact that f is the limit of a sequence with unit L -norm. (cid:3) In fact, investigating the stationary state we have already seen that one needs to imposeDirichlet conditions f (cid:16) cd (cid:17) = 0 = f (cid:16) ab (cid:17) . This is in accordance with the setting of Lemma 4.2.
Proposition 4.3.
Consider the operator A := A + A defined by (4.9) A (cid:18) f f (cid:19) ( x ) := − ddx (cid:18) a − b · c − d · (cid:19) (cid:18) f f (cid:19) ( x ) + (cid:18) − λ ( x ) µ ( x ) λ ( x ) − µ ( x ) (cid:19) (cid:18) f f (cid:19) ( x ) NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 11 with domain D ( A ) := (cid:26) f ≡ (cid:18) f f (cid:19) ∈ L (cid:16) ab , cd ; R (cid:17) : (cid:18)(cid:18) a − b · c − d · (cid:19) (cid:18) f f (cid:19)(cid:19) (cid:48) ∈ L (cid:16) ab , cd ; R (cid:17) and f (cid:16) cd (cid:17) = 0 , f (cid:16) ab (cid:17) = 0 (cid:111) . (4.10) Then A generates an irreducible C -semigroup of positivity preserving, stochastic contractionson L (cid:0) ab , cd ; R (cid:1) .Proof. Because λ, µ are L ∞ -functions, the multiplication operator A defined by the family ofmatrix functions M ( x ) := (cid:18) − λ ( x ) µ ( x ) λ ( x ) − µ ( x ) (cid:19) , x ∈ (cid:16) ab , cd (cid:17) , is bounded on L (cid:0) ab , cd ; R (cid:1) . In view of Lemma 4.2 the unperturbed operator (corresponding to A ≡
0) generates an irreducible semigroup of positivity preserving operators on L (cid:0) ab , cd ; R (cid:1) ,hence also the full operator A generates a semigroup on L (cid:0) ab , cd ; R (cid:1) .Positivity and irreducibility of the semigroup generated by A follows from the analogousproperty of the semigroup generated by A , see (4.4), and by a product formula analogous tothat in (4.3) [EN00, Exer. III.5.11].Finally, take a positive initial data f ∈ D ( A ) and observe that ddt (cid:107) e t A f (cid:107) L = ddt (cid:90) cdab ( e t A f ( x ) | ) R dx = (cid:90) cdab ( A e t A f ( x ) | ) R dx = (cid:90) cdab ( A e t A f ( x ) | ) R dx + (cid:90) cdab ( A e t A f ( x ) | ) R dx = 0 , with respect to the inner product ( ·|· ) of R , where in the last step the first integral vanishesintegrating by parts as in (4.7) and the second integral vanishes because lies in the null spaceof each matrix M ( x ) T . By density, this shows that (cid:107) e t A f (cid:107) L = (cid:107) f (cid:107) L for all t ≥ L -functions, i.e., e t A is stochastic and hence also contractive. (cid:3) Remark 4.4.
Our motivating model in (1.1) is based on one gene with two modes of expres-sion. Analogous models also exist that describe ensembles of genes with three or more modesof expression, leading to R n -valued functions with n ≥ . More precisely, (1.1) can begeneralized to a system of differential equations driven by the operator (4.11) A : f (cid:55)→ − ∂∂x (cid:104) diag (( a i − b i · ) f i ) i =1 ,...,n (cid:105) + M f, on L ( ab , cd ; R n ) , where M is an L ∞ -function taking values in the spaces of symmetric n × n -matrices whose rows sum up to 0 and with off-diagonal entries that are bounded below awayfrom 0. It is easy to see that our generation result extend to this more general setting. We leavethe details to the interested reader. By [Nag86, Cor. C.III.4.3] the eigenvalues that lie on i R form an additive cyclic group, i.e., σ ( A ) = i R or σ ( A ) = ik Z for some k ∈ N . Thus, the semigroup ( e t A ) t ≥ itself may a priori still converge towards a rotation group (i.e.,a linear combination of terms ( e tλ ) t ≥ , λ ∈ σ ( A ) ∩ i R ) as t →
0, as the subset σ ( A ) ∩ i R of thespectrum of A that lies on the imaginary axis may contain non-zero elements. An additionalcompactness argument is needed to rule out this case, but this turns out to be more delicatethan expected.First of all, observe that elements in the domain of A are W , , and in particularly L ∞ , oneach compact sub-interval of ( ab , cd ), hence in particular the resolvent operator R ( λ, A ) of A maps L ( ab , cd ; R ) to L ∞ loc ( ab , cd ; R ) for all λ >
0: by [Are08, Cor. 2.4], this implies that R ( λ, A )is a kernel operator for all λ >
0; we already know that it is positive as well. Now, let us recalla notion from the theory of operators on Banach lattices, based on the concept of
AM-spaces ,cf. [Sch74, § II.7]: a positivity preserving operator on L ( ab , cd ; R ) is said to be AM-compact ifit maps order intervals[ η, θ ] := (cid:110) g ∈ L (cid:16) ab , cd ; R (cid:17) : η ≤ g ( x ) ≤ θ for a.e. x ∈ (cid:16) ab , cd (cid:17)(cid:111) into precompact sets of L (cid:0) ab , cd ; R (cid:1) for all 0 ≤ θ ∈ L ( ab , cd ; R ). It is known that all positivekernel operators are AM-compact, cf. [GG, Prop. A.1]. Lemma 4.5.
The only eigenvalue of A = A + A along the imaginary axis is .Proof. By the Theorem of Perron–Frobenius for irreducible semigroups, cf. [EN06, Prop. VI.3.4],there is a unique element of ker A which is strongly positive and has norm 1; let us denote itby ( ψ , ψ ). Assume for a contradiction that ( f , f ) is a (normalised) eigenvector of A foran eigenvalue iβ where β ∈ R \ { } . Since λ and µ are continuous, it follows from ( f , f ) ∈ ker( iβ − A ) that f and f are C loc -functions. We use a standard argument from Perron–Frobenius theory to show that | ( f , f ) | = ( | f | , | f | ) ∈ ker A : indeed, we have | ( f , f ) | = | e t A ( f , f ) | ≤ e t A | ( f , f ) | for every t ≥
0. Since each semigroup operator e t A is contractive and since the norm on L ( ab , cd ; C ) is strictly monotone, it follows that actually | ( f , f ) | = e t A | ( f , f ) | for all t ≥ | ( f , f ) | ∈ ker A . From this we conclude that ( | f | , | f | ) = ( ψ , ψ ). We can thereforefind functions γ , γ : ( ab , cd ) → R such that f k = ψ k e iγ k for k = 1 ,
2. Since e iγ k is in C loc (as ψ and ψ are > ab , cd )), we can choose γ k to be in C loc , too.Now, it follows from [Nag86, Thm. C-III-2.2] that ( ψ e inγ , ψ e inγ ) ∈ ker( inβ − A ) for allintegers n ∈ Z . Using the definition of A , this yields after a short computation that β (cid:18) ψ ψ (cid:19) = − (cid:18) γ (cid:48) ( a − bx ) ψ γ (cid:48) ( c − dx ) ψ (cid:19) + in (cid:18) ddx [( a − bx ) ψ ] ddx [( c − dx ) ψ ] (cid:19) − in (cid:18) e − inγ e − inγ (cid:19) (cid:18) − λ µλ − µ (cid:19) (cid:18) ψ e inγ ψ e inγ (cid:19) (4.12) NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 13 for all n ∈ Z . By letting n → ∞ we thus obtain that β (cid:18) ψ ψ (cid:19) = (cid:18) γ (cid:48) ( a − bx ) ψ γ (cid:48) ( c − dx ) ψ (cid:19) . (4.13)Since ψ and ψ are > ab , cd ), it follows that γ (cid:48) = βa − bx and γ (cid:48) = βc − dx .On the other hand, by plugging (4.13) into (4.12) and using that A ( ψ , ψ ) = 0, we obtain (cid:18) − λ µλ − µ (cid:19) (cid:18) ψ ψ (cid:19) = (cid:18) e − inγ e − inγ (cid:19) (cid:18) − λ µλ − µ (cid:19) (cid:18) ψ e inγ ψ e inγ (cid:19) . A short computation thus yields (cid:18) ψ ψ (cid:19) = (cid:18) ψ e in ( γ − γ ) ψ e in ( γ − γ ) (cid:19) for all n ∈ Z . Plugging in n = 1 and using again that ψ and ψ are > ab , cd ), we thusconclude that γ − γ is constant. Hence, γ (cid:48) = γ (cid:48) and by using the expressions we derived for γ (cid:48) and γ (cid:48) above, we thus obtain βa − bx = βc − dx .Now we use the assumption that β (cid:54) = 0: this implies that a − bx = c − dx (for all x ∈ ( ab , cd )),and therefore a = c and b = d . Hence ab = cd , which is a contradiction. (cid:3) Now, we are finally in the position of stating our main result in this section.
Theorem 4.6.
The semigroup ( e t A ) t ≥ converges strongly towards the projector onto ker A . Our main step in the proof is to check that the semigroup is relatively compact in the strongoperator topology. Relative compactness of the orbits merely in the weak operator topologywould by [EN00, Cor. V.4.6] already imply mean ergodicity of the semigroup: but our result inTheorem 4.6 is significantly stronger, since convergence to equilibrium is actually achieved forindividual orbits not only in a time-averaged sense.
Proof.
In view of Lemma 4.5, the claim is an immediate consequence of [EN00, Thm. V.2.14] ifwe can prove that the orbits of the semigroup are relatively compact.To this purpose, let f be a strongly positive vector in the null space of A (and hence alsoa fixed point of e t A for all t ≥ f = f − A f , i.e., R (1 , A ) f = f .Introduce the set L f := (cid:110) L (cid:16) ab , cd ; R (cid:17) : ∃ c ∈ R + : | f | ≤ c | f | (cid:111) , which is dense in L (cid:0) ab , cd ; R (cid:1) , thus implying that R (1 , A ) (cid:0) L f (cid:1) is dense in R (1 , A ) (cid:0) L (cid:0) ab , cd ; R (cid:1)(cid:1) = D ( A ) and hence in L (cid:0) ab , cd ; R (cid:1) .Hence, in order to prove relative compactness of the orbits we can take without loss ofgenerality some f = R (1 , A ) g that additionally satisfies | g | ≤ cf . Then by using positivity of e t A we see that | e t A f | ≤ e t A | f | ≤ e t A R (1 , A ) | g | ≤ R (1 , A ) ce t A f = R (1 , A ) cf for all t ≥
0, and therefore e t A f ∈ [ − cR (1 , A ) f , cR (1 , A ) f ] ⊂ R (1 , A )[ − cf , cf ] for all t ≥ . All in all, we have proved that the orbits of the semigroup are contained in the image under R (1 , A ) of some order interval: because of AM-compactness of R (1 , A ), these sets are hencerelatively compact. (cid:3) Remark 4.7.
A close inspection of the proofs of Proposition 4.3 and Lemma 4.5 shows thatTheorem 4.6 holds not only for A as in (4.9) , but more generally under the assumptions that A is a bounded linear operator on L ( ab , cd ; R ) such that • A generates a positive, irreducible semigroup with A T = 0 ; and • DA − A D = Φ( e inγ , e inγ ) J for some injective bounded linear operator J and somefunction Φ : C → C that only vanishes along the diagonal { ( x, x ) : x ∈ C } , where D := (cid:18) e inγ e inγ (cid:19) . (The first conditions is equivalent to requiring the semigroup generated by A to be sto-chastic and irreducible.) We sum up our findings as follows.
Corollary 4.8.
The null space of A is spanned by a strongly positive function ψ and the semi-group generated by A on L ( ab , cd ; R ) converges strongly towards the orthogonal projector f (cid:55)→ (cid:90) cdab ( f ( x ) | φ ( x )) R dx · ψ , where ϕ is the strongly positive function that spans ker A (cid:48) and such that (cid:82) cdab ( ψ ( x ) | φ ( x )) R dx = 1 . Remark 4.9.
One may expect that the long-time behavior of (1.1) can be investigated by apply-ing classical results in the theory of positivity preserving operators like [EN00, Exer. II.4.21.(2)] , [EN06, Thm. VI.3.5] or [Dav05, Thm. 12] . However, we have not been able to prove that thesemigroup that governs our problem is • eventually norm continuous (cf. [EN00, Def. II.4.17] ), • or quasi-compact (cf. [EN06, Def. V.4.4] ), • or it has the Feller property (i.e., the solutions to (1.1) are continuous for all t > andeach initial data in L );in fact, we doubt that these properties hold at all. We also observe that if the semigroup couldbe proved to satisfy inf { e t A , e s A } > for some t > s ≥ , then by [Nag86, Cor. C-IV.2.10] theconvergence stated in Corollary 4.8 would hold in operator norm, too. Needless to say, the function ψ in Corollary 4.8 is nothing but the function ψ obtained inTheorem 3.3: this stationary solution is uniquely determined by the initial data f and theparameters λ, µ . We will present more explicit formulae for the strongly positive function ψ inthe following section, for special choices of λ, µ . NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 15 Analytical solutions and numerical examples
Let us discuss the simplest case: λ and µ are assumed to be affine functions, i.e.,(5.1) λ ( x ) = lx + k ,µ ( x ) = mx + n , for some k, l, m, n ∈ R chosen arbitrarily subject to guarantee positivity of λ and µ : then allintegrals in (3.12) can be explicitly calculated.We then have lx + kbx − a + mx + ndx − c = lb + lb ab + k/lx − ab + md + md cd + n/mx − cd and we get(5.2) h ( x ) = Ke ( lb + md ) x ( x − ab ) alb + kb ( cd − x ) cmd + nd . The solution to the system (3.1) is given by(5.3) ψ ( x ) = Kb e ( lb + md ) x (cid:16) x − ab (cid:17) alb + kb − (cid:16) cd − x (cid:17) cmd + nd ,ψ ( x ) = Kd e ( lb + md ) x (cid:16) x − ab (cid:17) alb + kb (cid:16) cd − x (cid:17) cmd + nd − . We see that the functions satisfy the boundary conditions (3.3), since(5.4) cmd + nd = 1 d µ (cid:16) cd (cid:17) > alb + kb = 1 b λ (cid:16) ab (cid:17) > . We have thus obtained an explicit, analytic solution to the stationary differential equation 3.1.This allows us to analyze the behavior of the solutions depending on the values of the parameters: • if al + kb < b , then ψ is singular at x = ab , • if al + kb = b , then ψ attains a nonzero value at x = ab , • if al + kb > b , then ψ tends to zero at x = ab ; • if cm + nd < d , then ψ is singular at x = cd , • if cm + nd = d , then ψ attains a nonzero value at x = cd , • if cm + nd > d , then ψ tends to zero at x = cd . Let us check where the maxima of ψ , ψ are situated in the case where these functions arenot singular at one of the endpoints. We calculate first the derivative of ψ ψ (cid:48) ( x ) = Kb e ( lb + md ) x (cid:16) x − ab (cid:17) ( alb + kb − ) (cid:16) cd − x (cid:17) ( cmd + nd − ) ·· (cid:20)(cid:18) lb + md (cid:19) (cid:16) x − ab (cid:17) (cid:16) cd − x (cid:17) − (cid:16) x − ab (cid:17) (cid:16) cmd + nd (cid:17) + (cid:16) cd − x (cid:17) (cid:18) alb + kb − (cid:19)(cid:21)(cid:124) (cid:123)(cid:122) (cid:125) := P ( x )=: Kb e ( lb + md ) x (cid:16) x − ab (cid:17) ( alb + kb − ) (cid:16) cd − x (cid:17) ( cmd + nd − ) · P ( x ) , (5.5)where P ( x ) is a quadratic polynomial. In the general position case al + kb > b and cm + nd ≥ d the function ψ is a non-negative function vanishing in the boundary points ab , cd , hence it musthave an odd number of local maxima in the open interval ( ab , cd ). But its critical points comefrom the zeroes of the quadratic polynomials P ( x ) , having at most two zeroes. Hence thefunction ψ has exactly one local maximum inside the interval. This maximum is necessarilyglobal.In the border case al + kb = b and cm + nd ≥ d the polynomial takes the form P ( x ) = (cid:0) lb + md (cid:1) (cid:0) x − ab (cid:1) (cid:0) cd − x (cid:1) − (cid:0) x − ab (cid:1) (cid:0) cmd + nd (cid:1) + (cid:0) cd − x (cid:1) (cid:18) alb + kb − (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) =0 = (cid:2)(cid:0) lb + md (cid:1) (cid:0) cd − x (cid:1) − (cid:0) cmd + nd (cid:1)(cid:3) (cid:0) x − ab (cid:1) . One of the roots coincides with the endpoint x = ab , hence ψ has at most one critical pointinside the interval. This point cannot be minimum point, hence the function ψ has exactly onemaximum inside the semi-closed interval [ ab , cd ) . The maximum is global and may coincide withthe left endpoint ab . This occurs if (cid:18) lb + md (cid:19) (cid:16) cd − ab (cid:17) − (cid:16) cmd + nd (cid:17) = lcbd − lab − mabd − nd ≤ ψ ψ (cid:48) ( x ) = Kd e ( lb + md ) x (cid:16) x − ab (cid:17) alb + kb − (cid:16) cd − x (cid:17) cmd + nd − P ( x ) , where P is quadratic polynomial as well. In the general position case al + kb ≥ b and cm + nd > d there is precisely one maximum inside the open interval. In the border case al + kb ≥ b and cm + nd = d , NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 17 there is one maximum in the semi-closed interval ( ab , cd ] , Note that the maximum may coincidewith the right endpoint.Let us turn to concrete examples of one-fluid systems. We compared our analytic solutionto numerical simulations for several examples (some use parameters which are biologically notmeaningful). For each example we provide the values of the parameters, analytic formula for thesolutions and plot the results. Since the numerical and analytical solutions are indistinguishable,we show only a single plot for each example. Moreover, we provide the values of the parameters,formulas for the normalized solutions and plots of these functions. Each time the functions ψ and ψ are plotted in green and blue, respectively, while their sum – the target probabilitydistribution of the system – is plotted in red. All computations were done using Matlab R (cid:13) . Example 1
Parameter values: a = − , b = d = c = 1 , λ ≡ , µ ≡ ψ ( x ) = (1 − x ) / ,ψ ( x ) = ( x + 1) / ,ψ ( x ) + ψ ( x ) = 1 / . x -1 -0.5 0 0.5 1 ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x -1 -0.5 0 0.5 1 ψ ( x ) Figure 1.
Solution of Example 1.
Example 2
Parameter values: a = − , b = d = c = 1 , λ ≡ , µ ≡ ψ ( x ) = 3( x + 1) (1 − x ) / ψ ( x ) = 3( x + 1)(1 − x ) / ψ ( x ) + ψ ( x ) = 3(1 − x ) / Example 3
Parameter values: a = − , b = d = c = 1 , λ ≡ / , µ ≡ / x -1 -0.5 0 0.5 1 ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x -1 -0.5 0 0.5 1 ψ ( x ) Figure 2.
Solution of Example 2.Normalized solution ψ ( x ) = 12 π √ − x √ x + 1 ψ ( x ) = 12 π √ x + 1 √ − xψ ( x ) + ψ ( x ) = 1 π √ − x x -1 -0.5 0 0.5 1 ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x -1 -0.5 0 0.5 1 ψ ( x ) Figure 3.
Solution of Example 3.Let us now turn to the case of affine functions λ ( x ) = lx + k , µ ( x ) = mx + n . Example 4
Parameter values: a = − , b = d = c = 1 , l = 0 , k = 1 , m = − , n = 3 NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 19
Normalized solution ψ ( x ) = e e + 1) e − x (1 − x ) ψ ( x ) = e e + 1) e − x (1 − x ) ψ ( x ) + ψ ( x ) = ee + 1 e − x (1 − x ) x -1 -0.5 0 0.5 1 ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x -1 -0.5 0 0.5 1 ψ ( x ) Figure 4.
Solution of Example 4.
Example 5
Parameter values: a = 0 , b = 1 , c = 2 , d = 1 , l = − , k = 4 , m = 1 , n = 0Normalized solution ψ ( x ) = e − e ) e − x x (2 − x ) ψ ( x ) = e − e ) e − x x (2 − x ) ψ ( x ) + ψ ( x ) = e − e ) e − x x (2 − x ) Example 6
Parameter values: a = − , b = c = d = 1 , l = 0 , k = 1 , m = 2 , n = 1Normalized solution ψ ( x ) = ( − − e x (1 − x ) ψ ( x ) = ( − − e x ( x + 1)(1 − x ) ψ ( x ) + ψ ( x ) = ( − − ( e x (1 − x ) + e x ( x + 1)(1 − x ) ) x ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x ψ ( x ) Figure 5.
Solution of Example 5. x -1 -0.5 0 0.5 1 ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x -1 -0.5 0 0.5 1 ψ ( x ) Figure 6.
Solution of Example 6.
Example 7
Parameter values: a = − , b = c = d = 1 , l = 0 , k = 1 , m = 1 , n = 2Normalized solution ψ ( x ) = 8 e − e + 48 e + e e x (1 − x ) ψ ( x ) = 8 e − e + 48 e + e e x ( x + 1)(1 − x ) ψ ( x ) + ψ ( x ) = 8 e − e + 48 e + e ( e x (1 − x ) + e x ( x + 1)(1 − x ) )In all considered examples the functions λ and µ were assumed to be affine for simplicity.Even this choice provided us with a rather rich class of models. Note that any other simpleanalytic expression will do the job. One should only assure that the integrals in (3.12) can becalculated analytically. The last two examples were chosen to illustrate the power of analyticcalculations: in both cases ψ does not approach zero and is not singular at the left endpoint. In NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 21 x -1 -0.5 0 0.5 1 ψ ( x ) ψ (x) ψ (x) ψ (x)+ ψ (x) x -1 -0.5 0 0.5 1 ψ ( x ) Figure 7.
Solution of Example 7.one case it has a maximum inside the interval, in the other case it is monotonically decreasing.It was easy to find parameters leading to such behavior since analytic formulas were available,this would be a challenging task if only computer simulations were available.6.
Conclusion
We considered a hybrid model of a self-regulating gene, which is a common motif in generegulatory networks. Our model describes the evolution of the discrete random state (mode)of the gene (“on” or “off”) and the corresponding continuous protein concentration. The latterevolves according to an ordinary differential equation and leads to a system of PDEs for theevolution of two probability densities (one for each mode). Assuming that the rate functions λ and µ for mode changing are known explicitly, we analyzed the properties of the PDE system andstudied well-posedness in an L -setting. Exploiting the theory of positive operator semigroupswe rigorously proved convergence towards stationary solutions in strong operator topology andderived an analytic expression for such stationary densities. Our solution is valid for a large classof protein production and degradation rates and structurally much simpler and easier to evaluatethan the solution of the corresponding fully discrete master equation model given by Grima etal. [GSN12]. As future work, we plan to investigate the extension of this gene feedback loop twoor more interacting gene and their corresponding proteins such as the exclusive or toggle switch[LLBB06, LLBB07]. In these cases, the support of the stationary solution of the correspondinghybrid model has a more complex shape. For instance, in the case of two genes and three modes(as one mode is not reachable), the density is only non-zero within a triangle whose endpointsare determined by the mode-conditional equilibria of the two protein concentrations. Althoughthis is straightforward to see from Monte-Carlo simulations of the model, proving this and otherproperties for the corresponding PDE system is challenging.7. Acknowledgements
This work is the outcome of a collaboration between two mathematicians working on evolu-tionary equations (DM) and partial differential equations (PK) and a computer biologist (VW).It has partially been carried out at the ZiF (Center for Interdisciplinary Research) in Bielefeld in the framework of the Cooperation Group
Discrete and continuous models in the theory ofnetworks . The authors are indebted to the ZiF for financial support and hospitality. The workof PK was also partially supported by the Swedish Research Council (Grant D0497301). Thework of VW was partially supported by the German Research Foundation (Grant 391984329).We warmly thank Jochen Gl¨uck (Ulm) for suggesting to us the proofs of Lemma 4.5 and The-orem 4.6 and for many helpful discussions.
References [Are08] W. Arendt. Positive semigroups of kernel operators.
Positivity , 12:25–44, 2008.[CDR09] Alina Crudu, Arnaud Debussche, and Ovidiu Radulescu. Hybrid stochastic simplifications for mul-tiscale gene networks.
BMC Syst. Biol. , 3(1):89, 2009.[Dav05] E.B. Davies. Triviality of the peripheral point spectrum.
J. Evol. Equ. , 5:407–415, 2005.[EN00] K.-J. Engel and R. Nagel.
One-Parameter Semigroups for Linear Evolution Equations , volume 194of
Graduate Texts in Mathematics . Springer-Verlag, New York, 2000.[EN06] K.-J. Engel and R. Nagel.
A Short Course on Operator Semigroups . Universitext. Springer-Verlag,Berlin, 2006.[GG] M. Gerlach and J. Gl¨uck. Convergence of positive operator semigroups. arxiv:1705.01587.[GSN12] R. Grima, D.R. Schmidt, and T.J. Newman. Steady-state fluctuations of a genetic feedback loop:An exact solution.
J. Chem. Phys. , 137:035104, 2012.[GT01] D. Gilbarg and N. Trudinger.
Elliptic Partial Differential Equations of Second Order . Classics inMathematics. Springer-Verlag, Berlin, 2001.[HGK15] Benjamin Hepp, Ankit Gupta, and Mustafa Khammash. Adaptive hybrid simulations for multiscalestochastic reaction networks.
The Journal of chemical physics , 142(3):034118, 2015.[HH12] Mostafa Herajy and Monika Heiner. Hybrid representation and simulation of stiff biochemical net-works.
Nonlinear Analysis: Hybrid Systems , 6(4):942–959, 2012.[HKNT98] G. Horton, V.G. Kulkarni, D.M. Nicol, and K.S. Trivedi. Fluid stochastic petri nets: Theory, appli-cations, and solution techniques.
Eur. J. Oper. Res. , 105:184–201, 1998.[HSI +
05] JEM Hornos, D Schultz, GCP Innocentini, JAMW Wang, AM Walczak, JN Onuchic, andPG Wolynes. Self-regulating gene: An exact solution.
Phys. Rev. E , 72(5):051907, 2005.[JH07] Tobias Jahnke and Wilhelm Huisinga. Solving the chemical master equation for monomolecularreaction systems analytically.
J. Math. Biol. , 54(1):1–26, 2007.[KLMW18] P. Kurasov, A. L¨uck, D. Mugnolo, and V. Wolf. Stochastic hybrid models of gene regulatory net-works.
Math. Biosci. , 305:170–177, 2018.[KPK14] Niraj Kumar, Thierry Platini, and Rahul V Kulkarni. Exact distributions for stochastic gene ex-pression models with bursting and feedback.
Phys. Rev. Lett. , 113(26):268105, 2014.[Lau00] Ian J. Laurenzi. An analytical solution of the stochastic master equation for reversible bimolecularreaction kinetics.
J. Chem. Phys. , 113(8):3315–3322, 2000.[LLBB06] Azi Lipshtat, Adiel Loinger, Nathalie Q Balaban, and Ofer Biham. Genetic toggle switch withoutcooperative binding.
Physical review letters , 96(18):188101, 2006.[LLBB07] Adiel Loinger, Azi Lipshtat, Nathalie Q Balaban, and Ofer Biham. Stochastic simulations of geneticswitch systems.
Phys. Rev. E , 75(2):021904, 2007.[LYWZ16] Peijiang Liu, Zhanjiang Yuan, Haohua Wang, and Tianshou Zhou. Decomposition and tunability ofexpression noise in the presence of coupled feedbacks.
Chaos , 26(4):043108, 2016.[MS13] J. Mi¸ekisz and P. Szyma´nska. Gene expression in self-repressing system with multiple gene copies.
Bull. Math. Biol. , 75:317–330, 2013.[Nag86] R. Nagel, editor.
One-Parameter Semigroups of Positive Operators , volume 1184 of
Lect. NotesMath.
Springer-Verlag, Berlin, 1986.
NALYTIC SOLUTIONS FOR STOCHASTIC HYBRID MODELS 23 [PK04] Jacek Pucha(cid:32)lka and Andrzej M Kierzek. Bridging the gap between stochastic and determinis-tic regimes in the kinetic simulations of the biochemical reaction networks.
Biophysical Journal ,86(3):1357–1372, 2004.[RS78] M. Reed and B. Simon.
Methods of Modern Mathematical Physics - IV: Analysis of Operators .Academic Press, San Diego (CA), 1978.[RvO08] Arjun Raj and Alexander van Oudenaarden. Nature, nurture, or chance: stochastic gene expressionand its consequences.
Cell , 135(2):216–226, 2008.[Sch74] H.H. Schaefer.
Banach Lattices and Positive Operators , volume 215 of
Grundlehren der mathema-tischen Wissenschaften . Springer-Verlag, Berlin, 1974.[SSG17] David Schnoerr, Guido Sanguinetti, and Ramon Grima. Approximation and inference methods forstochastic biochemical kinetics—a tutorial review.
Journal of Physics A: Mathematical and Theo-retical , 50(9):093001, 2017.[TK93] K.S. Trivedi and V.G. Kulkarni. FSPNs: fluid stochastic Petri nets. In
Application and Theoryof Petri Nets (14th international conference, Chicago, IL) , Lect. Notes Comp. Sci., pages 24–31,Berlin, 1993. Springer.[VAE08] Paolo Visco, Rosalind J Allen, and Martin R Evans. Exact solution of a model DNA-inversiongenetic switch with orientational control.
Phys. Rev. Lett. , 101(11):118104, 2008.[VB13] Yves Vandecan and Ralf Blossey. Self-regulatory gene: an exact solution for the gene gate model.
Phys. Rev. E , 87(4):042705, 2013.
Pavel Kurasov, Dept. of Mathematics, Stockholm Univ., 106 91 Stockholm, SWEDEN
E-mail address : [email protected] Delio Mugnolo, Chair of Analysis, University of Hagen, 58084 Hagen, GERMANY
E-mail address : [email protected] Verena Wolf, Computer Science Department, Saarland University, 66123 Saarbr¨ucken, GER-MANY
E-mail address ::