Guilherme C P Innocentini

Arran Hodgkinson

Fernando Antoneli

Arnaud Debussche

Ovidiu Radulescu

Push-forward method for piecewise deterministic biochemical simulations

A biochemical network can be simulated by a set of ordinary differential equations (ODE) under well stirred reactor conditions, for large numbers of molecules, and frequent reactions. This is no longer a robust representation when some molecular species are in small numbers and reactions changing them are infrequent. In this case, discrete stochastic events trigger changes of the smooth deterministic dynamics of the biochemical network. Piecewise-deterministic Markov processes (PDMP) are well adapted for describing such situations. Although PDMP models are now well established in biology, these models remain computationally challenging. Previously we have introduced the push-forward method to compute how the probability measure is spread by the deterministic ODE flow of PDMPs, through the use of analytic expressions of the corresponding semigroup. In this paper we provide a more general simulation algorithm that works also for non-integrable systems. The method can be used for biochemical simulations with applications in fundamental biology, biotechnology and biocomputing.

Introduction

Stochastic simulation is a powerful tool in biology. Moreover, stochasticity is a general property of biochemical networks, having multiple origins. Noise is generated intrinsically by these molecular systems when neither the law of large numbers nor the averaging theorem can be applied, for instance when some of the constituent molecular species are present in small numbers and trigger relatively slow reactions [START_REF] Crudu | Hybrid stochastic simplifications for multiscale gene networks[END_REF]. There are also extrinsic sources of noise, resulting from a fluctuating cellular environment. The extrinsic noise paradigm also applies to stochasticity resulting from manipulating biochemical networks in an artificial environment such as lab on a chip or simply in a titration device [START_REF] Abate | Experimental biological protocols with formal semantics[END_REF].

Stochastic simulation is used in single cell experimental studies, where the amounts of mRNA [START_REF] Thattai | Stochastic gene expression in fluctuating environments[END_REF][START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF][START_REF] Cai | Stochastic protein expression in individual cells at the single molecule level[END_REF][START_REF] Tantale | A single-molecule view of transcription reveals convoys of RNA polymerases and multi-scale bursting[END_REF] and protein [START_REF] Elowitz | Stochastic gene expression in a single cell[END_REF][START_REF] Ferguson | Reconciling molecular regulatory mechanisms with noise patterns of bacterial metabolic promoters in induced and repressed states[END_REF] products of a gene can be determined for each cell. By double-or multiple-fluorophore fluorescence techniques, products from several genes can be quantified simultaneously and one can have access to multivariate probability distributions of mRNA or proteins. The stochastic dynamics of promoters and gene networks can have important consequences for fundamental biology [START_REF] Eldar | Functional roles for noise in genetic circuits[END_REF] but also for HIV [START_REF] Razooky | A hardwired HIV latency program[END_REF] and cancer research [START_REF] Gupta | Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells[END_REF]. Predicting the probability distributions of biochemical networks' products is also instrumental for lab-on-a-chip applications, when one wants to optimize and control the functioning of these networks. For these reasons we aim to develop effective methods for computing time-dependent distributions for stochastic models. Our main objective is the reduction of computation time which is prerequisite for parameter scans and machine learning applications [START_REF] Herbach | Inferring gene regulatory networks from single-cell data: a mechanistic approach[END_REF].

The traditional stochastic simulation protocol uses the chemical master equation and the Gillespie algorithm. In such a protocol all the chemical reactions are simulated individually as discrete stochastic events. Simulations on relevant experimental times have to cope with 10 6 -10 9 such events and generate 10 2 -10 4 samples in order to get statistically significant estimates of molecular species distributions. Altogether, these simulations are extremely costly in terms of execution time.

An important reduction of the simulation time is obtained by noticing that the sources of noise can be localized in small parts of the system that behave discretely, or in the discrete environmental variables in the case of intrinsic or extrinsic noise, respectively. The remaining, large part of the system consists of molecular species in large numbers that evolve continuously. This leads to piecewise-deterministic Markov process (PDMP) approximations of the biochemical dynamics, coupling discrete state Markov chain dynamics with continuous state ordinary differential equations (ODE) dynamics. The justification of the PDMP approximations of the chemical master equation can be found in [START_REF] Crudu | Hybrid stochastic simplifications for multiscale gene networks[END_REF][START_REF] Crudu | Convergence of stochastic gene networks to hybrid piecewise deterministic processes[END_REF]. Although simpler than the chemical master equation, direct simulation of the PDMP remains time consuming because the Markov chains can still have a very large number of states [START_REF] Crudu | Hybrid stochastic simplifications for multiscale gene networks[END_REF].

In the CMSB2019 proceedings paper we have introduced new methods for simulating PDMPs for gene networks [START_REF] Innocentini | Effective computational methods for hybrid stochastic gene networks[END_REF]. A gene network PDMP model can be simulated by numerical integration of ODEs satisfied by the mRNA and the protein species, coupled to the Markov chain describing the successive transitions of the gene promoters [START_REF] Zeiser | Simulation of genetic networks modelled by piecewise deterministic Markov processes[END_REF][START_REF] Crudu | Hybrid stochastic simplifications for multiscale gene networks[END_REF][START_REF] Riedler | Almost sure convergence of numerical approximations for piecewise deterministic Markov processes[END_REF][START_REF] Ting Lin | Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic markov processes[END_REF]. The simulation becomes particularly effective when analytic solutions of the ODEs are available [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF].

Probability distributions of PDMP are solutions of the Liouville-master partial differential equations (PDEs). Numerical integration of these PDEs is an interesting alternative to direct simulation, combining precision and speed for small models. Finite difference methods, however, are of limited use in this context as they can not cope with high dimensional models (for instance, extant gene networks applications are restricted to the dimension 2, corresponding to a single promoter, with or without self-regulation see [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF][START_REF] Kurasov | Stochastic hybrid models of gene regulatory networks -a PDE approach[END_REF]).

Another interesting method for computing time dependent distributions is the push-forward method. For gene networks, this method has been first introduced in [START_REF] Innocentini | Protein synthesis driven by dynamical stochastic transcription[END_REF] and further adapted for continuous mRNA variables in [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF]. It is based on the idea to compute the probability distribution of gene products as the push-forward measure of the semigroup defined by the ODEs. This method is an approximation, as one has to consider that the discrete PDMP variables are piecewise constant on a deterministic time partition. The transition rates between promoter states were computed using a mean field approximation in [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF]. In the CMSB2019 proceeding paper, the mean field approximation was replaced by the next order approximation taking into account the second moment of the protein distribution [START_REF] Innocentini | Effective computational methods for hybrid stochastic gene networks[END_REF]. In this paper we extend the method to general PDMP models, thus covering all biochemical networks, with both intrinsic and extrinsic noise simulation protocols.

Models

ODE models of biochemical networks

A biochemical network is defined by a set of chemical reactions and a set of chemical species. The amounts of different chemical species form a vector x = (x 1 , x 2 , . . . , x N ) ∈ R N + , where x i is the concentration of the species i. Each reaction is characterized by a stoichiometric vector ν i ∈ Z N and a reaction rate function R i : R N + → R. When all chemical species are present in large numbers, the biochemical network is well described by a systems of ODEs

dx dt = r i=1 ν i R i (x), (1)
For example, a simple gene transcription model can be obtained when a large number, let's say G, of a specific gene is present in the cell and, moreover, the gene promoter has just two states: active and inactive. In this simple model the gene product is the mRNA concentration. The vector of species concentration is given by: x 1 = P (inactive promoter concentration), x 2 = P * (active promoter concentration) and x 3 = X (mRNA concentration). The concentration of inactive promoters is measured by the ratio between the number of inactive promoters (G I ) and the cell volume V cell , so P = G I /V cell ; the same holds for the concentrations of active promoters and mRNA, P = G A /V cell , where G A is the number of active promoters. The reactions are P → P * , P * → P , P * → P * + X, X → ∅. Thus, the stoichiometric vectors are ν 1 = (-1, 1, 0), ν 2 = (1, -1, 0), ν 3 = (0, 0, 1), ν 4 = (0, 0, -1) and reaction functions

R 1 = k 1 P , R 2 = k 2 P * , R 3 = k 3 P * , R 4 = k 4 X, where the reaction rate constants k i , 1 ≤ i ≤ 4, have dimension [1/time].
For large copy numbers of P , P * and X the simple transcription model reads:

dP dt = -k 1 P + k 2 P * , dP * dt = k 1 P -k 2 P * , dX dt = k 3 P * -k 4 X . (2)
Note that dP dt + dP * dt = 0 and so one of the first two equations can be discarded. We must emphasize that gene transcription is only one of the many possible examples of our formalism. The formalism defined by (1) covers all biochemical network models and have large applicability in cell physiology.

PDMP models

When some, but not all, species are present in low numbers and there are also slow reactions, a stochastic representation, such as a piecewise deterministic Markov process (PDMP), is more appropriate than the deterministic equations [START_REF] Abate | Experimental biological protocols with formal semantics[END_REF].

For instance, in the simple transcription model (2), assume that there is only one copy of the gene, meaning G = G A + G I = 1. In this case, it does not make sense to describe the time evolution of promoter variables by differential equations anymore. In this scenario, the promoter variables are discrete numbers. As the copy number of the gene is one, is more suitable to replace the concentration P by the discrete variable G A assuming values G A = 1 or G A = 0. The same holds for P * , that is replaced by G I ∈ {0, 1}. Moreover, since G A + G I = 1 one of these variables can be discarded.

The switching between the two discrete states of G A can be described by a continuous time Markov chain with transition rates k 1 and k 2 . That is, if

G A (0) = 1, then G A (t) = 1 for 0 ≤ t < T 1 where T 1 is a random time such that P[T 1 > t] = exp(-k 2 t) and G A (t) = 0 for T 1 ≤ t < T 1 + T 2 where T 2 is a random time such that P[T 2 > t] = exp(-k 1 t).
Moreover, suppose that the switching constants, k 1 and k 2 , are small compared to k 4 , and also that k 3 /k 4 1. Thus, the variable X has a switching behavior alternating accumulation and degradation periods when G A = 1 and G A = 0, respectively, each one described by ODEs:

dX dt = k 3 /V cell -k 4 X when G A = 1 -k 4 X when G A = 0 , (3)
where V cell is the cell volume. Using methods from [START_REF] Crudu | Convergence of stochastic gene networks to hybrid piecewise deterministic processes[END_REF] it is possible to show that the "hybrid system" (G A /V cell , X), with X given by (3), converges to (2), as G → ∞.

The PDMP formalism can be extended to a rather general class of biochemical network models as follows:

A PDMP biochemical model is a stochastic process ζ t having states in a set of the form E = R N + × S, where S is a finite set encoding discrete states of the biochemical model, and ζ t = (x t , s t ), where x t is a vector in R N + whose components x i t (1 ≤ i ≤ N ) encode the dynamics of i continuous biochemical species, and s t describes the jump Markov process between the discrete states. The PDMP ζ t = (x t , s t ) is determined by three characteristics:

1) For all fixed s ∈ S, a vector field

V s : R N + → R N determining a unique global flow Φ s (t, x) in R N + , such that, for t > 0, dΦ s (t, x) dt = V s (Φ s (t, x)), Φ s (0, x) = x. (4)
The flow Φ s (t, x) represents a one parameter semigroup fulfilling the properties

(i) Φ s (0, x 0 ) = x 0 , (ii) Φ s (t + t , x 0 ) = Φ s (t , Φ s (t, x 0 )).
2) A transition rate matrix

H : R N + → M Ns×Ns (R), such that H r,s (x) is the (r, s) element of the matrix H, N s = #S is the number of states. If (s = r), H r,s (x)
is the rate of probability to jump to the state r from the state s. Furthermore, H s,s (x) =r =s H r,s (x) for all s ∈ S and for all x ∈ R N .

3) A jump rate λ : E → R + . The jump rate can be obtained from the transition rate matrix λ(x, s) =

r =s H r,s (x) = -H s,s (x). (5)
From these characteristics, right-continuous sample paths {x t : t > 0} starting at ζ 0 = (x 0 , s 0 ) ∈ E can be constructed as follows. Define

x t (ω) := Φ s0 (t, x 0 ) for 0 ≤ t ≤ T 1 (ω), (6)
where T 1 (ω) is a realisation of the first jump time of s t , with the distribution

F (t) = P[T 1 > t] = exp - t 0 λ(Φ s0 (u, x 0 ), s 0 ) du , t > 0, (7)
and ω is the element of the probability space for which the particular realisation of the process is given. The pre-jump state is ζ T - 1 (ω) (ω) = (Φ s0 (T 1 (ω), x 0 ), s 0 ) and the post-jump state is ζ T1(ω) (ω) = (Φ s0 (T 1 (ω), x 0 ), s), where s has the distribution

P[s = r] = H r,s0 (Φ s0 (T 1 (ω), x 0 ), s 0 ) λ(Φ s0 (T 1 (ω), x 0 ), s 0 )
, for all r = s 0 .

We then restart the process ζ T1(ω) and recursively apply the same procedure at jump times T 2 (ω), etc.. Note that between each two consecutive jumps x t follow deterministic ODE dynamics defined by the vector field V s . At the jumps, the values x t are continuous. More general definitions of PDMPs include jumps in the continuous variables but will not be discussed here.

We define multivariate probability density functions p s (t, x). These functions satisfy the Liouville-master equation which is a system of partial differential equations:

∂p s (t, x) ∂t = -∇ x .(V s (x)p s (t, x)) + s H s,s (x)p s (t, x). (9)
This general formalism covers all biochemical models that have discrete variables. As salient examples we can cite gene networks [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF], ion channels dynamics for neuron networks [START_REF] Anderson | Stochastic Representations of Ion Channel Kinetics and Exact Stochastic Simulations of Neuronal Dynamics[END_REF], lab on chip biochemical devices [START_REF] Abate | Experimental biological protocols with formal semantics[END_REF].

PDMP models of ON/OFF gene networks

A particular example of PDMP biochemical model is represented by the gene networks. Because these models are extensively used in systems biology, we provide their details in this subsection.

Two state ON/OFF gene networks generalize the simple two state single gene transcription model, by considering more interacting genes and by using separate variables for the two gene products: mRNA and protein.

To each gene we associate a discrete variable σ i with two possible values σ i = 0 for a non-productive state (OFF) and σ i = 1 for a productive state (ON). Furthermore, a gene i produces proteins and mRNAs in the amounts y i and r i , respectively. A gene network state is described by the N -dimensional vector x t = (r 1 (t), y 1 (t), r 2 (t), y 2 (t), . . . , r N/2 (t), y N/2 (t)). The gene products amounts satisfy ODEs:

dy i dt = b i r i -a i y i , dr i dt = k i (σ i ) -ρ i r i (10)
The coupling between genes is modeled at the level of discrete state transitions. The elements of the matrix H are functions of products from various genes. For instance, if a gene i inhibits a gene j the transition rate from ON to OFF of the gene j is an increasing function of the protein concentration y i . As a first example that we denote as model M 1 , let us consider a two genes network; the expression of the first gene being constitutive and the expression of the second gene being activated by the first. We consider that the transcription activation rate of the second gene is proportional to the concentration of the first protein f 2 y 1 . All the other rates are constant f 1 , h 1 , h 2 , representing the transcription activation rate of the first gene, and the transcription inactivation rates of gene one and gene two, respectively. For simplicity, we consider that the two genes have identical protein and mRNA parameters b 1 = b 2 = b, a 1 = a 2 = a, ρ 1 = ρ 2 = ρ. We further consider that k i = k 0 if the gene i is OFF and

k i = k 1 > k 0 if the gene i is ON.
The gene network has four discrete states, in order (0, 0), (1, 0), (0, 1), and (1, 1). Then, the transition rate matrix for the model M 1 is

H(y 1 , y 2 ) =     -(f 1 + f 2 y 1 ) h 1 h 2 0 f 1 -(h 1 + f 2 y 1 ) 0 h 2 f 2 y 1 0 -(f 1 + h 2 ) h 1 0 f 2 y 1 f 1 -(h 1 + h 2 )     . ( 11
)
The Liouville-master equation for the model M 1 reads

∂p 1 ∂t = - ∂[(br 1 -ay 1 )p 1 ] ∂y 1 - ∂[(k 0 -ρr 1 )p 1 ] ∂r 1 - ∂[(br 2 -ay 2 )p 1 ] ∂y 2 - ∂[(k 0 -ρr 2 )p ] ∂r 2 + h 2 p 3 + h 1 p 2 -(f 1 + f 2 y 1 )p 1 , ∂p 2 ∂t = - ∂[(br 1 -ay 1 )p 2 ] ∂y 1 - ∂[(k 1 -ρr 1 )p 2 ] ∂r 1 - ∂[(br 2 -ay 2 )p 2 ] ∂y 2 - ∂[(k 0 -ρr 2 )p ] ∂r 2 + f 1 p 1 + h 2 p 4 -(h 1 + f 2 y 1 )p 2 , ∂p 3 ∂t = - ∂[(br 1 -ay 1 )p 3 ] ∂y 1 - ∂[(k 0 -ρr 1 )p 3 ] ∂r 1 - ∂[(br 2 -ay 2 )p 3 ] ∂y 2 - ∂[(k 1 -ρr 2 )p ] ∂r 2 + h 1 p 4 + f 2 y 1 p 1 -(h 2 + f 1 )p 3 , ∂p 4 ∂t = - ∂[(br 1 -ay 1 )p 4 ] ∂y 1 - ∂[(k 1 -ρr 1 )p 4 ] ∂r 1 - ∂[(br 2 -ay 2 )p 4 ] ∂y 2 - ∂[(k 1 -ρr 2 )p ] ∂r 2 + f 1 p 3 + f 2 y 1 p 2 -(h 1 + h 2 )p 4 .
(12) The model M 2 differs from the model M 1 by the form of the activation function. Instead of a linear transcription rate f 2 y 1 we use a Michaelis-Menten model f 2 y 1 /(K 1 + y 1 ). This model is more realistic as it takes into account that the protein p 1 has to attach to specific promoter sites which become saturated when the concentration of this protein is high.

The transition rate matrix for the model M 2 is

H(y 1 , y 2 ) =      -(f 1 + f2y1 K1+y1 ) h 1 h 2 0 f 1 -(h 1 + f2y1 K1+y1 ) 0 h 2 f2y1 K1+y1 0 -(f 1 + h 2 ) h 1 0 f2y1 (K1+y1) f 1 -(h 1 + h 2 )      . (13)
The Liouville-master equation for the model M 2 reads

∂p 1 ∂t = - ∂[(br 1 -ay 1 )p 1 ] ∂y 1 - ∂[(k 0 -ρr 1 )p 1 ] ∂r 1 - ∂[(br 2 -ay 2 )p 1 ] ∂y 2 - ∂[(k 0 -ρr 2 )p 1 ] ∂r 2 + h 2 p 3 + h 1 p 2 -(f 1 + f 2 y 1 /(K 1 + y 1 ))p 1 , ∂p 2 ∂t = - ∂[(br 1 -ay 1 )p 2 ] ∂y 1 - ∂[(k 1 -ρr 1 )p 2 ] ∂r 1 - ∂[(br 2 -ay 2 )p 2 ] ∂y 2 - ∂[(k 0 -ρr 2 )p 2 ] ∂r 2 + f 1 p 1 + h 2 p 4 -(h 1 + f 2 y 1 /(K 1 + y 1 ))p 2 , ∂p 3 ∂t = - ∂[(br 1 -ay 1 )p 3 ] ∂y 1 - ∂[(k 0 -ρr 1 )p 3 ] ∂r 1 - ∂[(br 2 -ay 2 )p 3 ] ∂y 2 - ∂[(k 1 -ρr 2 )p 3 ] ∂r 2 + h 1 p 4 + f 2 y 1 /(K 1 + y 1 )p 1 -(h 2 + f 1 )p 3 , ∂p 4 ∂t = - ∂[(br 1 -ay 1 )p 4 ] ∂y 1 - ∂[(k 1 -ρr 1 )p 4 ] ∂r 1 - ∂[(br 2 -ay 2 )p 4 ] ∂y 2 - ∂[(k 1 -ρr 2 )p 4 ] ∂r 2 + f 1 p 3 + f 2 y 1 /(K 1 + y 1 )p 2 -(h 1 + h 2 )p 4 .

Simulation methods

Monte-Carlo method

The Monte-Carlo method utilizes the direct simulation of the PDMP based on the iteration of the following equations:

dx dt = V s0 (x(t)), dF dt = -λ(x, s 0 )F, (14)
for t ∈ [0, T 1 ) with initial conditions x(0) = x 0 , F (0) = 0, and stopping condition F (T 1 ) = U , where U is random variable, uniform in the range [0, 1), followed by the choice of the next discrete state s 1 by using the matrix H(x(T 1 )).

A large number N mc of sample traces is generated and the values of x t are stored at selected times. Multivariate time dependent probability distributions are then estimated from this data.

The complete algorithm is presented in the Algorithms 1 and 2.

The direct simulation of PDMPs needs the solutions of Eqs. ( 14) which can be obtained by numerical integration. This is not always computationally easy. Problems may arise for fast switching promoters when the ODEs have to be integrated many times on small intervals between successive jumps. Alternatively, the numerical integration of the ODEs can be replaced by analytic solutions or quadrature. Analytic expressions are always available for the gene network flow Algorithm 1 NextState Input: 1. x initial value continuous variable, 2. s initial discrete state, 3. CUM := CUM + H I,s0 (x)/ J =s0 H J,s0 (x) 9: end while 10: return ((x 0 , x 1 , . . . , x n ), (t 0 , t 1 , . . . , t n ), I) [START_REF] Ferguson | Reconciling molecular regulatory mechanisms with noise patterns of bacterial metabolic promoters in induced and repressed states[END_REF] and read

t 0 initial time, 4. t max maximal time Output: 1. (x 0 , x 1 , . . . , x n ), 2. (t 0 , t 1 , . . . , t n ) 3. s new discrete state 1: x 0 := x; s 0 := s 2: U := randunif ([0, 1)) 3: solve dx dt = V s0 (x), dF dt = -λ(x, s 0 )F with initial data x(t 0 ) = x 0 ,F (t 0 ) = 1, until F (t) = U or t = t max , return x 0 = x(t 0 ), . . . , x n = x(t n ), where t 0 < t 1 < . . . < t n , t n = min(F -1 (U ), t max ). 4: V := randunif ([0, 1)
Φ y i (t, y 0 , r 0 ) = y 0 exp(-a i t) + b i r 0 - k i (σ i ) ρ i exp(-ρ i t) -1 a i -ρ i + k i (σ i ) ρ i 1 -exp(-a i t) a i , Φ r i (t, y 0 , r 0 ) = (r 0 -k i /ρ i ) exp(-ρ i t) + k i /ρ i .
Let us consider the following general expression of the jump intensity function

λ(x, s) = c 0 (s) + N i c i (s)y i + N i d i (s)f i (y i ),
where f i are non-linear functions, for instance Michaelis-Menten function

f i (p i ) = p i K i + p i or Hill functions f i (y i ) = (y i ) ni K ni i + (y i ) ni . If d i = 0 for all 1 ≤ i ≤ N ,
the cumulative distribution function of the waiting time T 1 can be solved analytically [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF], otherwise it can be obtained by quadrature. For example, for the model M 2 one has

λ(x, s) = f 1 + f 2 y 1 K 1 + y 1 δ s,1 + h 1 + f 2 y 1 K 1 + y 1 δ s,2 + (h 2 + f 1 )δ s,3 + (h 2 + h 1 )δ s,4 ,
Algorithm 2 PDMPmontecarlo Input: 1. p 0 (x, s), initial distribution of (x 0 , s 0 ). 2. τ time step for computing the time-dependent distribution. 3. δx bin size in one direction. Output: n t vectors H 0 , . . . , H nt , representing the distribution at times 0, τ, . . . , n t τ 1: for i := 1 to n t do 2:

H i := n M x dimensional null vector 3: end for 4: for m := 1 to N mc do 5:

Draw (x 0 , s 0 ) from the distribution p 0 (x, s) 6:

x := x 0 ,s := s 0 7:

t := t 0 ; i := 0 8:
while t < n t τ do 9:

((x 0 , . . . , x n ), (t 0 , . . . , t n ), s next ) := NextState(x, s, t, t M ) 10:

for j := 1 to n do 11: if t j ≤ iτ < t j+1 then 12: X i := (x j + x j+1 )/2 13:
increment by 1/((δx) N N mc ), the bin corresponding to X i in the histogram H i 14:

i := i + 1 15: end if 16:
end for where δ i,j is Kronecker's delta. In this case, the waiting time T 1 is obtained as the unique solution of the equation

-log(U ) = (f 1 + f 2 )T 1 + f 2 T1 0 1 K 1 + Φ y 1 (t , y 0 , r 0 ) dt δ s0,1 + (h 1 + f 2 )T 1 + f 2 T1 0 1 K 1 + Φ y 1 (t , y 0 , r 0 ) dt δ s0,2 + (h 2 + f 1 )T 1 δ s0,3 + (h 2 + h 1 )T 1 δ s0,4 , (15)
where U is a random variable, uniformly distributed in the range [0, 1). In our implementation of the algorithm we solve (15) numerically, using the bisection method.

Finite difference Liouville master equation method

The finite difference Liouville master equation method for a given number of genes, N g , uses a discrete approximation of the domain to compute the numerical solution across time for a given higher-dimensional system [START_REF] Herbach | Inferring gene regulatory networks from single-cell data: a mechanistic approach[END_REF], with initial conditions given by p 1 (0, y 1 , r 1 , . . .) := δ(y 1 , r 1 , . . .) and p 2 (0, r 1 , y 1 , . . . ) := . . . := p 2Ng (0, y 1 , r 1 , . . .) := 0 where δ(•) is the Dirac delta function. We compute the solution for 2N g distributions, since each gene has both an ON and OFF state, respectively.

In order to achieve the simulation of the equations given by system [START_REF] Herbach | Inferring gene regulatory networks from single-cell data: a mechanistic approach[END_REF], we begin by discretising each of the domains into n η intervals, whose centres are given by {η 1 , . . . , η nη } ∈ y meeting the condition that η 1 < η 2 < . . . < η nη and where, likewise, an arbitrary discrete value in y i would be denoted η i,j . Each ηi, j then represents a unique, discrete, j th position in the protein abundance domain spanned by y i . Likewise, mRNA would have an analogous discretisation given by all ρ i,j , ∀j ∈ {1, . . . , n ρ }. Time is then similarly discretised by a timestep τ into n t + 1 temporal locations given by {0, τ, . . . , τ n t }. The task then becomes the computation of the solution at each of these discrete locations in the domain, such that the sink and source terms may be trivially calculated but where the derivative terms warrant further explanation.

The solution to the advection equation under a uniform coefficient, χ, is given by a translation in the relevant domain. To achieve this, whilst also maintaining the stability of the system, as a whole, we implement a simultaneous forwardsand backwards-difference discrete operator scheme. This means that for an arbitrary probability density function, p k (t, y i , . . .), we write the discrete partial derivative operator as

χ ∂ ∂y i p k (t, η i,j , . . .) =        |χ| p k (t, η i,j-1 , . . .) -p k (t, η i,j , . . .) η i,j -η i,j-1 if χ > 0 |χ| p k (t, η i,j+1 , . . .) -p k (t, η i,j , . . .) η i,j+1 -η i,j otherwise.
(16) This may be evaluated as such for each term within the system of equations, given by [START_REF] Herbach | Inferring gene regulatory networks from single-cell data: a mechanistic approach[END_REF], and guarantees the probability balance of the system, as a whole.

Beyond the computation of the equation's solutions at a single time-step the solutions must be computed robustly across time. We therefore employ a McCormack predictor-corrector scheme [START_REF] Maccormack | The Effect of Viscosity in Hypervelocity Impact Cratering Frontiers of Computational Fluid Dynamics[END_REF], given explicitly for any i th timestep by

p k,i+1 = p k,i + τ F k (p 1,i , . . . , p Ng,i ) p k,i+1 = 1 2 (p k,i + p k,i+1 ) + 1 2 τ F k (p 1,i+1 , . . . , p Ng,i+1 ) ∀k ∈ {1, . . . , 2N g }, (17)
where p k,i = p k (t i , y 1 , r 1 , . . .), F k (p 1 , . . . , p 2Ng ) = ∂p k /∂t and p k,i is the discrete nomenclature for a prediction of the solution for p k at time t i .

The method for solving the problem, in totality, is then given by Algorithm 3 where, again, the partial derivative terms on the right-hand side of each equation are evaluated using [START_REF] Innocentini | Effective computational methods for hybrid stochastic gene networks[END_REF]. Algorithm 3 is a direct implementation of ( 17) as an algorithm with the concurrent calculation of the distributions P t , ∀t ∈ {0, τ, . . . , n t τ } from the individual distributions p k,i , ∀k ∈ {1, . . . , 2N g }, i ∈ {0, . . . , n t }. For each time-step, we solve the right-hand side of the equation using the predictor corrector scheme, update the value of each distribution, and calculate the total distribution, P t .

Algorithm 3 PDMPLiouvillemaster

Input: 1. p 1,0 , p 2,0 , p 3,0 , p 4,0 initial distributions in (y for j := 1 to 2N g do 3:

Use p k,i-1 and discretised right hand side of PDE in [START_REF] Herbach | Inferring gene regulatory networks from single-cell data: a mechanistic approach[END_REF] to compute the temporal gradient, ∂ t p j,i-1

4:

Compute predicted solution at next time step as p j,i = p j,i-1 +τ ∂ t p j,i-1

5:

end for

6:

Set P i = 0 7:

for j := 1 to 2N g do 8:

Use prediction p k,i and discretised right hand side of PDE in 12 to compute the temporal gradient, ∂ t p j,i

9:

Compute corrected solution at next time step as p j,i = 1 2 (p j,i-1 + p j,i ) + 1 2 τ ∂ t p j,i 10:

Set P i = P i + p j,i 11:
end for 12: end for 13: return (P 0 , . . . , P nt ) 3.3. Push-forward method

General algorithm

This method allows one to compute the multivariate probability distribution of the continuous variable x at a time τ given the probability distribution of (x, s) at time 0.

In order to achieve this we consider a deterministic partition

τ 0 = 0 < τ 1 < . . . < τ M = τ of the interval [0, τ ] such that ∆ M = max j∈[1,M ] (τ j -τ j-1 )
is small. The main approximation of this method is to assume that s t , for t ∈ [0, τ ], is piecewise constant on this partition, more precisely, that s t = s j := s τj , for t ∈ [τ j , τ j+1 ), 0 ≤ j ≤ M -1. This is rigorously true for intervals Algorithm 4 PDMPpushforward Input: 1. p 0 (x, s), initial distribution of (x 0 , s 0 ). 2. τ time step for computing the time-dependent distribution Output: n t vectors H 0 , . . . , H nt , representing the distribution at times 0, τ, . . . , n t τ 1: compute H 0 , P 0 from initial distribution 2: for i := 1 to n t do 3:

for j := 1 to n N x do 4:

for s0 := 1 to N s do 5:

solve dx0 dt = V s0 (x0), dΠ0 dt = H(x0)Π0 with initial conditions x0(0) = C j , Π0(0) = I from t = 0 to t = τ 1 6:
for s1 := 1 to N s do 7:

solve dx1 dt = V s1 (x1), dΠ1 dt = H(x1)Π1 with initial conditions x1(τ 1 ) = x0(τ 1 ), Π1(0) = I from t = τ 1 to t = τ 2 8:
for s2 := 1 to N s do 9:

solve dx2 dt = V s2 (x2), dΠ2 dt = H(x2)Π2 with initial conditions x2(τ 2 ) = x1(τ 2 ), Π2(0) = I from t = τ 2 to t = τ 3 10:
for s3 := 1 to N s do 11:

solve dx3 dt = V s3 (x3) with initial conditions x3(τ 3 ) = x2(τ 3 ) from t = τ 3 to t = τ 4 12: P M = Π2 s3,s2 Π1 s2,s1 Π0 s1,s0 P 0 s0 13: x = x3(τ 4 )
14: end for 20: end for [τ j , τ j+1 ) completely contained between two successive random jump times of s t . This situation becomes very frequent for a very fine partition (large M ). Thus, the error generated by the approximation vanishes in the limit M → ∞ (the rigorous result is Theorem 1 given in the Results section).

H i (BIN (x)) = H i (BIN (x)) + P M H i-1 (j) BIN (
For each path realization S M := (s 0 , s 1 , . . . , s M -1 ) ∈ Ω := S M of the discrete states, we can compute x(t), t ∈ [0, τ ] as the continuous solution of the following ODE with piecewise-defined r.h.s:

dx dt = V sj (x), for t ∈ [τ j , τ j+1 ), 0 ≤ j ≤ M -1 (18)
and with the initial condition x(0) = x 0 .

In order to compute the probability P[S M ] of a path realization we can use the fact that, given x t , s t is a finite state Markov process. Therefore,

s 0 =1 s 1 =1 s 2 =1 s 3 =0 s 4 =0 t x 0 τ 1 τ 2 τ 3 τ 4 τ
P[S M ] = Π s M -1 ,s M -2 (τ M -2 , τ M -1 ) . . . Π s1,s0 (τ 0 , τ 1 )P S 0 (s 0 ), (19)
where P S 0 : S → [0, 1] is the initial distribution of the discrete state, and Π(τ j , τ j+1 ) is the solution, at t = τ j+1 , of

dΠ(τ j , t) dt = H(x t )Π(τ j , t) (20)
with Π(τ j , τ j ) = I and x t is given by [START_REF] Ting Lin | Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic markov processes[END_REF]. In order to compute the probability distribution of x at time τ one has to sum the contributions of all solutions of (18), obtained for the N M s realisations of promoter state paths with weights given by the probabilities of the paths.

Suppose that we want to estimate the distribution of x(τ ), using a multivariate histogram with bin centers c (l1,...,l N ) = (c l1 1 , . . . , c l N N ), 1 ≤ l i ≤ n x , 1 ≤ i ≤ N where n x is the number of bins in each direction x i , 1 ≤ i ≤ N . The initial distribution of x at time t = 0 is given by the bin probabilities p (l1,...,l N ) 0 , 1 ≤ i ≤ N, 1 ≤ l i ≤ n x and the distribution at time τ are given by the probabilities p

(l1,...,l N ) τ , 1 ≤ i ≤ N, 1 ≤ l i ≤ n x .
Let x (l1,...,l N ) (t) be the solution of (18) with x (l1,...,l N ) (0) = c (l1,...,l N ) and let (l 1 , . . . , l N ) be the histogram bin containing x (l1,...,l N ) (τ ). The application (l 1 , . . . , l N ) → (l 1 , . . . , l N ) := ψ(l 1 , . . . , l N ) is in general many to one. Given the probability P[S M ] of a path S M ∈ Ω, the push forward distribution of x(τ ) is computed as

p (l 1 ,...,l N ) τ = S M ∈Ω ψ(l1,...,l N )=(l 1 ,...,l N ) p (l1,...,l N ) 0 P[S M ] , (21)
The push-forward method can be applied recursively to compute the distribution at times τ, 2τ, . . . , n t τ . The complete algorithm is schematized in Algorithm 4 for the choice M = 4.

Remark 1. The Algorithm 4 can be used for computing the full multivariate probability distribution of the vector x(t), but also for computing marginal distributions. For gene network the full multivariate distribution implies products from all genes, whereas a marginal distribution can select only one, or a small number of genes. For marginals, the dimension N is replaced by N m < N where N m is the number of components of interest of the vector x(t).

Remark 2. For gene networks, the numerical integration of the ODEs (steps 7,9,11 in the Algorithm 4) can be replaced by symbolic solutions. For each path realization S M := (s 0 , s 1 , . . . , s M -1 ) ∈ Ω := {0, 1, . . . , 2 N -1} M of promoter states, we can compute the protein and mRNA levels, y t and r t , respectively, of all genes i ∈ {1, N }, at t = τ :

r i τ = r i 0 e -ρτ + k 0 ρ (1 -e -ρτ ) + k 1 -k 0 ρ M -1 j=1 e -ρτ (e -ρτj+1 -e -ρτj )s i j (22)
y i τ = y i 0 e -aτ + br i 0 a -ρ (e -ρτ -e -aτ ) + bk 0 ρ 1 -e -aτ a + e -aτ -e -ρτ a -ρ + b(k 1 -k 0 ) ρ e -aτ M j=1 s i j-1 w j , (23)
for i ∈ {1, . . . , N }. Here, w j = e (a-ρ)τ -e (a-ρ)τj a -ρ (e ρτj -e ρτj-1 ) -e (a-ρ)τj -e (a-ρ)τj-1 a -ρ e ρτj-1 + e aτj -e aτj-1 a with s i j := 0 if promoter i is OFF for t ∈ [τ j , τ j+1 ) and s i j := 1 if promoter i is ON for t ∈ [τ j , τ j+1 ).

Complexity issues

The complexity of the push forward algorithm scales as n t n N x N M s (N + N 2 s ), because there are N M s possible paths S M , n N

x histogram bins, and n t time complexity for solving [START_REF] Ting Lin | Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic markov processes[END_REF], [START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF] is O(N + N 2 s ). The complexity of computing marginal distributions of N m < N variables is lower and scales as n t n Nm x N M s (N + N 2 s ).

Mean field push-forward method

A way to mitigate the computational burden of the Algorithm 4 is to use the mean field approximation. In the mean field approximation, the probabilities P[S M ] are computed from averaged equations that are identical for each histogram bin.

More precisely, equation ( 20) is replaced by

dΠ(t , t) dt = E[H(x t )] Π(t , t) (24)
where t ≤ t, Π(t , t ) = I.

Using a Taylor expansion of H(x t ) around the expectation E[x t ], one gets

dΠ(t , t) dt = H(E[x t ]) + 1 2 H (E[x t ]) : Var(x t ) Π(t , t), (25)
where Var(x t ) is the variance/covariance matrix of x t , H is the element-wise second derivative matrix of H and : stands for the double dot product. In [START_REF] Tantale | A single-molecule view of transcription reveals convoys of RNA polymerases and multi-scale bursting[END_REF] the moments E[x t ] and Var(x t ) are either available analytically or are solutions of ODEs obtained using moment closures such as in [START_REF] Singh | Stochastic hybrid systems for studying biochemical processes[END_REF].

The main advantage of the mean-field approximation is that the transition matrix elements can be computed outside the discrete variables loop. The complete algorithm is presented in the Algorithm 5. The time complexity in this case is n t n N x N M s N , because in the innermost loop one has to solve only [START_REF] Ting Lin | Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic markov processes[END_REF], whose time complexity is O(N ).

For gene networks, the mean field push-forward procedure can be applied also to ODE dynamics of the individual genes. This is possible because during ON or OFF periods, the ODE dynamics of one gene is uncoupled from that of another gene. Furthermore, the N s × N s transition matrix H can be replaced by N/2 2 × 2 transition matrices of one gene with elements averaged over the values of the other genes. This approximation reduces the complexity of the calculations to 2n t (n x ) 2 N 2 M , which is linear in the number of genes. The mean field push-forward algorithm for gene networks is presented in the Algorithm 6.

This method has already been used for particular models. In [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF] we have replaced the regulation term f 2 y 1 (t) occurring in the transition matrix of the gene network model M 1 by its mean f 2 E[y 1 (t)]. This means that the gene 2 switches between its ON and OFF states with rates given by the mean of the regulatory protein y 1 . In this case both H and Π can be computed analytically, which leads to a drastic reduction in the execution time. This simple mean field approach is suitable for the model M 1 , which contains only linear regulation terms. For non-linear regulation terms, Π can not generally be computed analytically. Moreover, the naive mean field approach introduces biases. For instance, in the case of the model M 2 , the approximation

f 2 y 1 (t)/(K 1 +y 1 (t)) ≈ f 2 E[y 1 (t)] /(K 1 + E[y 1 (t)]
) is poor. In the CMSB2019 paper we proposed a better approximation [START_REF] Innocentini | Effective computational methods for hybrid stochastic gene networks[END_REF], in which we replace f 2 y 1 (t)/(K 1 + y 1 (t)) by its mean value and use

E f 2 y 1 (t) K 1 + y 1 (t) ≈ f 2 E[y 1 (t)] (K 1 + E[y 1 (t)]) - f 2 (K 1 + E[y 1 (t)]) 3 Var(y 1 (t)), (26)
in order to correct the bias. This approach is generalized by [START_REF] Tantale | A single-molecule view of transcription reveals convoys of RNA polymerases and multi-scale bursting[END_REF].

Algorithm 5 PDMPpushforwardmeanfield

Input: 1. p g 0 (x, s), initial distribution of (x 0 , s 0 ). 2. τ time step for computing the time-dependent distribution Output: n t vectors H 0 , . . . , H nt representing the distribution of x at times 0, τ, . . . , n t τ 1: compute H 0 , P 0 from initial distribution 2: for i := 1 to n t do 3:

for s0 := 1 to N s do 4:

for s1 := 1 to N s do 5:

for s2 := 1 to N s do 6:

for s3 := 1 to N s do 7:

Compute P M (s 0 , s 1 , s 2 , s 3 ) = Π s3,s2 (τ 2 , τ 3 )Π s2,s1 (τ 1 , τ 2 )Π s1,s0 (τ 0 , τ 1 )P 0 s0 where Π(t , t) is the solution of (25). for s0 := 0 to 1 do 13:

for j := 1 to n N x do 14:

solve dx0 dt = V s0 (x0), with initial conditions x0(0) = C j , from t = 0 to t = τ 1 15:
for s1 := 0 to 1 do for s3 := 0 to 1 do 20: end for 28: end for Algorithm 6 PDMPpushforwardmeanfieldgenenetwork Input: 1. for each gene g ∈ {1, 2, . . . , N/2}, p 0 (x, s), initial distribution of (x 0 , s 0 ), where x = (r, p). 2. τ time step for computing the time-dependent distribution Output: for each gene g ∈ {1, 2, . . . , N/2}, n t vectors H 0 , . . . , H nt representing the distribution at times 0, τ, . . . , n t τ 1: for g := 1 to N/2 do 2:

solve dx3 dt = V s3 (x3) with initial conditions x3(τ 3 ) = x2(τ 3 ) from t = τ 3 to t = τ 4 21: x = x3(τ 4 ) 22: H i (BIN (x)) = H i (BIN (x))+P M (s 0 , s 1 , s 2 , s 3 )H i-1 (j) where BIN (x)
compute H 0 , P 0 from initial distribution for gene g 3:

for i := 1 to n t do 4:

for s0 := 0 to 1 do 5:

for s1 := 0 to 1 do 6:

for s2 := 0 to 1 do 7:

for s3 := 0 to 1 do 8:

Compute P M (s 0 , s 1 , s 2 , s 3 ) = Π s3,s2 (τ 2 , τ 3 )Π s2,s1 (τ 1 , τ 2 )Π s1,s0 (τ 0 , τ 1 )P 0 s0
where Π(t , t) is the solution of ( 25) for gene g.

9:

end for for s0 := 1 to 2 do 14:

for j := 1 to n 2 x do 15:

solve dx0 dt = V s0 (x0), with initial conditions x0(0) = C j , from t = 0 to t = τ 1 , and V s is defined by [START_REF] Ferguson | Reconciling molecular regulatory mechanisms with noise patterns of bacterial metabolic promoters in induced and repressed states[END_REF] for gene g for s2 := 1 to N s do 19:

solve dx2 dt = V s2 (x2) with initial conditions x2(τ 2 ) = x1(τ 2 ) from t = τ 2 to t = τ 3 20:
for s3 := 1 to N s do 21:

solve dx3 dt = V s3 (x3) with initial conditions x3(τ 3 ) = x2(τ 3 ) from t = τ 3 to t = τ 4 22: x = x3(τ 4 ) 23: H i (BIN (x)) = H i (BIN (x)) + P M (s 0 , s 1 , s 2 , s 3 )H i-1 (j)
where BIN (x) is the bin containing x end for 30: end for As in [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF] we can use analytic expressions for E[y 1 (t)], but also for Var(y 1 (t)). These expressions can be found in the Appendix A. Although the elements of matrix H have analytic expressions, the elements of the matrix Π contain integrals that must be computed numerically. For the model M 2 , we have

Π 1 (τ, τ ) = (1 -p 1,on ) + p 1,on e -1(τ -τ ) (1 -p 1,on )(1 -e -1 (τ -τ ) ) p 1,on (1 -e -1 (τ -τ ) ) p 1,on + (1 -p 1,on )e -1 (τ -τ ) , ( 27
)
for the transition rates of the first gene, where

p 1,on = f 1 /(f 1 + h 1 ), 1 = (f 1 + h 1 )/ρ, and
Π 2 (τ, τ ) = K(τ, τ ) + h 2 τ τ K(t, τ ) dt h 2 τ τ K(t, τ ) dt 1 -K(τ, τ ) -h 2 τ τ K(t, τ ) dt 1 -h 2 τ τ K(t, τ ) dt , ( 28
)
for the transitions of the second gene, where K(τ, τ ) = e -τ τ (h2+F2(t)) dt and

F 2 (t) = f 2 E y1(t)
K1+y1(t) .

Results

Convergence of the push-forward method

The probability distribution obtained with the push-forward method converges to the exact PDMP distribution in the limit M → ∞. This is a consequence of the following theorem: Theorem 1. Let Φ S M (t, x) be the flow defined by the formulas [START_REF] Ting Lin | Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic markov processes[END_REF], such that

x t = Φ S M (t, x(0)) for t ∈ [0, τ ]. Let µ M t : B(R N ) → R + be the probability meas- ure defined as µ M t (A) = S M ∈Ω P[S M ] µ 0 (Φ -1 S M (t, A))
, where µ 0 : B(R N ) → R + is the probability distribution of x at t = 0, P[S M ] is given by [START_REF] Maccormack | The Effect of Viscosity in Hypervelocity Impact Cratering Frontiers of Computational Fluid Dynamics[END_REF], and B(R N ) is the σ-algebra of Borel sets on R N . Let µ t , the exact distribution of x t for the PDMP defined in Section 2.2, with initial values (x 0 , s 0 ) distributed according to µ 0 × P S 0 . Assume that the vector fields V s , s ∈ S, and the transition matrix H are C 1 functions and that there exists a bounded set of R N such that all flows Φ s , s ∈ S, leave B invariant. Assume that

|τ i -τ i-1 | < C/M for all i ∈ [1, M ],
where C is a positive constant. Then, for all t ∈ [0, τ ], µ M t converges in distribution to µ t , when M → ∞. More precisely, for all Lipschitz functions ϕ on R N , there exists a constant κ depending on the data of the PDMP, τ and the Lipschitz constant of ϕ such that:

B ϕ(x)µ M t (dx) - B ϕ(x)µ t (dx) = |E(ϕ(Φ S M (t, x))) -E(ϕ(x t ))| ≤ κ/M.
Also for all Borel set A in B(R N ) such that µ t (∂A) = 0 we have µ M t (A) → µ t (A) when M → ∞. Here, ∂A denotes the boundary of A.

The proof of this theorem is given in Appendix B. It is inspired by the classical proof of weak order of convergence for the Euler scheme for a stochastic differential equation (see [START_REF] Talay | Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution[END_REF]). Remark 3. If the flows Φ s are not known explicitly, a numerical scheme can be used. This introduces another source of error. Our proof easily extends and provide a similar result of convergence. If one wants to investigate the convergence of p (l1,...,l N ) τ

given by ( 21), this follows from the above theorem only under the assumption that the probability that the PDMP is on the boundary of the bins at time t is 0. Note this is not restrictive and happens only in pathological situations. Also, if the initial distribution has a smooth density, we can prove that the error estimate above holds for borelian bounded functions ϕ, thus we can choose ϕ to be an indicator function and obtain error bounds for these probabilities without this restriction. More precisely, we have for any Borel set A:

µ M t (A) -µ t (A) ≤ κ/M
, where now κ depends on the Lipschitz constant of the initial density, see Remark 4 in Appendix B.

4.2.

Testing the performance of various methods 4.2.1. Testing accuracy and speed of push-forward method compared to the Monte-Carlo method In order to test the push-forward method, we compared the resulting probability distributions with the ones obtained by the Monte-Carlo method using the direct simulation of the PDMP. We considered the models M 1 and M 2 with the following parameters:

ρ = 1, p 1 = f f +h = 1/2, a = 1/5, b = 4, k 0 = 4, k 1 =
40 for the two genes. For the parameter = f +h ρ we used two values: = 0.5 for slow genes and = 5.5 for fast genes. We tested the slow-slow and the fast-fast combinations of parameters.

The initial distribution of the promoters state was P S 0 ((0, 0)) = 1 where the state (0, 0) means that both promoters are OFF. The initial probability measure µ 0 was a delta Dirac distribution centered at r 1 = r 2 = 0 and y 1 = y 2 = 0. This initial condition is obtained by always setting the direct simulation of the PDMP to start from r 1 (0) = r 2 (0) = 0, y 1 (0) = y 2 (0) = 0, and s 1 0 = s 2 0 = 0. The simulations were performed between t 0 = 0 and t max = 20 for fast genes and between t 0 = 0 and t max = 90 for slow genes. In order to estimate the distributions we have used M C = 50000 samples for the highest sampling. The push-forward method was implemented with M = 10 equal length sub-intervals of [0, τ ]. The time step τ was chosen to be τ = 2 for fast genes and τ = 15 for slow genes. The procedure was iterated 10 times for fast genes (up to t max = 20) and 6 times for slow genes (up to t max = 90).

The execution times are shown in the Table 1. The comparison of the probability distributions are illustrated in the Figures 2 and3. In order to quantify the relative difference between methods we use the L 1 distance between distributions. More precisely, if p(x) and p(x) are probability density functions to be compared, the distance between them is Table 1: Execution times for different methods. All the methods were implemented in Matlab R2013b running on a single core (multi-threading inactivated) of a Intel i5-700u 2.5 GHz processor. The Monte-Carlo method computed the next jump waiting time (Algoritm 1) using the analytical solution of ( 15) for M 1 and the numerical solution of ( 15) for M 2 . The pushforward method used Algorithm 6 and analytic solutions for mRNA and protein trajectories from ( 23), ( 18) and numerical computation of the integrals in ( 28), for both models. This distance was computed for distributions resulting from the push-forward method and the Monte-Carlo method with the highest sampling. We have also used a reduced sampling Monte-Carlo scheme whose execution time is similar to the one of the push-forward method. The distributions resulting from low sampling and high sampling Monte-Carlo were compared using the same distance. Figures 2 and3 clearly show that, for the same execution time, the push-forward method outperforms the Monte-Carlo method.

d = |p(x) -p(x)| dx. (29)

Comparing the Monte Carlo, push-forward and Liouville-master methods

Using Python, we have implemented Algorithm 4 for the one gene model, to compute mRNA (r(t)) and protein (y(t)) probability densities P (r, t) and P (y, t), respectively. For the one gene model, we have considered two switching regimes: slow ( = 0.5) and fast ( = 5.5). Also in Python, we have implemented Algorithm 5 for the two genes model, to compute the probability densities for mRNA and protein associated to gene one (P (r 1 , t) and P (y 1 , t)) and for mRNA and protein associated to gene two (P (r 2 , t) and P (y 2 , t)). For the two genes model, we have considered four different switch configurations: slow-slow, fastfast, slow-fast and fast-slow. For all implementations each time interval [0, τ ] has been partitioned into four sub-intervals of equal sizes (M = 4) resulting in the sequence {s 0 , s 1 , s 2 , s 3 }, representing the state of the switch inside each sub-interval (s = s j := s tj for t ∈ [t j , t j+1 ), j = 0, ..., M -1), leading to 2 4 path realizations. For all slow genes we have set t 0 = 0, t 4 = 9 and evaluated the solution up to t max = 90, using the composition rule between successive time intervals. For all fast genes we have set t 0 = 0, t 4 = 1 and evaluated the solution up to t max = 20, using the composition rule between successive time intervals, as well.

Model

Monte We have compared the results obtained by the implementation of the push forward method (PF), as described above, with Monte-Carlo simulation (MC) (Algorithms 1 and 2) and with numerical solution of the Liouville-master equation (LME) (Algorithm 3). The comparison between execution times for each model and each method can be found in Table 2 (all the execution times are expressed in minutes).

In order to illustrate our results we have produced Figures 4,5, 6 and 7, for selected models and for methods: Push forward (PF), Monte-Carlo simulation (MC) and Liouville-master equation (LME).

The one gene model (Figures 4 and5) has as initial conditions p 0 (0) = 1, where p 0 (0) is the probability to find the switch in state 0 at time t = 0, and r(0) = y(0) = 0, for both switch regimes: slow and fast. The two genes model (Figures 6 and7) has as initial conditions p 1 0 (0) = p 2 0 (0) = 1 (meaning that at time t = 0 both genes are in state 0 with probability one), r 1 (0) = y 1 (0) = r 2 (0) = y 2 (0) = 0, for both switch configurations: slow-slow and slow-fast. We used the same parameters for all the models and for all the genes: 1 = 2 = 0.5 (slow switch), 1 = 2 = 5.5 (fast switch), p 1 0 = p 2 0 = 0.5,

ρ 1 = ρ 2 = 1, k 1 (σ 1 = 0) = k 2 (σ 2 = 0) = 4, k 1 (σ 1 = 1) = k 2 (σ 2 = 1) = 40, a 1 = a 2 = 1/5 and b 1 = b 2 = 4.
We have quantified the relative difference between the distributions obtained by two distinct methods using the L 1 distance (29). The distance between the distributions obtained by the Push Forward method and the Monte-Carlo method (high sampling) is indicated in the plots by d P F , and the distance between the distributions obtained by the Liouville-master equation and Monte-Carlo method (high sampling) is indicated by d LM E .

Discussion and conclusion

Combining direct simulation of PDMP gene network models and analytic formulae for the ODE flow provides an effective, easy to implement method for computing time dependent, multivariate probability distributions of these models. However, the precision of the Monte-Carlo estimates of the distributions increases with √ M C, where M C is the number of Monte-Carlo samples. For this reason, the execution time of the Monte-Carlo method, although smaller when compared to PDMP simulation methods, which implement numerical resolution of the ODEs, such as reported in [START_REF] Ting Lin | Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic markov processes[END_REF] (data not shown), is large when compared to the push-forward method. The push-forward method represents an effective alternative to Monte-Carlo and Liouville-master equation methods, ensuring reduced execution time.

With respect to an earlier implementation of this method for gene networks in [START_REF] Innocentini | Protein synthesis driven by dynamical stochastic transcription[END_REF] we used promoter states instead of mRNA copy numbers as discrete variables of the PDMP. As a consequence, the number of discrete states is lower and we can afford to increase the number M of temporal subdivisions. Compared to the similar work done in [START_REF] Innocentini | Time dependent stochastic mRNA and protein synthesis in piecewise-deterministic models of gene networks[END_REF] we used second moments of the protein distribution, which took into account the correlation of the promoter states and lead to increased accuracy in the case of nonlinear regulation. Although the protein moments and the exponential transition rate matrix Π can be computed numerically, the effectiveness of the push-forward method is increased when analytic expressions are available for these quantities. In this paper, these expressions were computed for particular cases.

The push-forward method is an approximate method, and its accuracy relies on the careful choice of the temporal and spatial step densities, namely of the integers M , n t , n x . We have rigorously proved that the convergence of the push-forward method is of order one, with errors that scale with 1/M . We situate our findings in the broader effort of the community to produce new effective simulation algorithms for computational biology. Although in this manuscript we restricted ourselves to examples of modeling gene networks, the push forward method and the algorithms of this paper can be applied to a broader class of phenomena where some continuous quantity is perturbed by a discrete stochastic process. A very interesting and important example is the stochastic representation of ion channel kinetics in neuron models [START_REF] Anderson | Stochastic Representations of Ion Channel Kinetics and Exact Stochastic Simulations of Neuronal Dynamics[END_REF].

The system of equations describing the kinetics of an ion channel is mathematically equivalent to the system of equations describing our one gene model. More complex neuron models with multiple types of ion channels (stochastic versions of Hodgkin-Huxley or Morris-Lecar deterministic models, for instance) are also covered by our PDMP formalism and algorithms. Thus, the push forward method as presented in Algorithm 4 can be used to obtain the histograms for the membrane voltage for rather general neuron models. A detailed comparison between the push forward method and the one presented in [START_REF] Anderson | Stochastic Representations of Ion Channel Kinetics and Exact Stochastic Simulations of Neuronal Dynamics[END_REF] will be the subject of a future work.

Appendices Appendix A. Expectation and variance of the protein

For the sake of completeness, we give explicit formulas for the expectation and the variance of the protein synthesized by a constitutive promoter (gene 1 of models M 1 and M 2 ). The details of the calculation can be found in [START_REF] Innocentini | Effective computational methods for hybrid stochastic gene networks[END_REF].

The expectation is given by

E[y t ] = M 0 + M 1 e -at + M 2 e -ρt + M 3 e -t , (A.1)
where

M 0 = b(k 0 + (k 1 -k 0 )p 1 ) a , M 1 = E[x 0 ] - bE[y 0 ] a -ρ + bk 0 a(a -ρ) + b(k 1 -k 0 )(p 10 -p 1 ) (a -ρ )(a -ρ) + b(k 1 -k 0 )p 1 a(a -ρ) , M 2 = bE[y 0 ] a -ρ - bk 0 a -ρ - b(k 1 -k 0 )(p 10 -p 1 ) ρ(1 -)(a -ρ) - b(k 1 -k 0 )p 1 a -ρ , M 3 = b(k 1 -k 0 )(p 10 -p 1 ) ρ(1 -)(a -ρ) .
The variance is given by

Var(y t ) = Var(y 0 )e -2at + b 2 Var(r 0 ) e -ρt -e -at a -ρ 2 - (p 10 -p 1 )(k 1 -k 0 )b ρ(1 -) e -ρ t -e -at a -ρ - e -t -e -at a -ρ 2 + p 1 (1 -p 1 )(k 1 -k 0 ) 2 b 2 ρ 2 (V 0 + V 1 e -(a+ρ )t + V 2 e -ρ(1+ )t + V 3 e -2at + V 4 e -(a+ρ)t + V 5 e -2ρt ) + (1 -2p 1 )(p10 -p 1 )(k 1 -k 0 ) 2 b 2 ρ 2 (V 6 e -ρ t + V 7 e -(a-ρ )t + V 8 e -ρ( +1)t + V 9 e -(a+ρ )t + V 10 e -ρt + V 11 e -2ρt ), (A.2)
and the forward Kolmogorov equation:

du dt (t, x, s) = Lu(t, x, s), x ∈ R, s ∈ S, t > 0.
When u(0, x, s) = ϕ(x, s) for x ∈ R, s ∈ S, it is well known that u(t, x 0 , s 0 ) = E(ϕ(x(t, ζ 0 ), s(t, ζ 0 ))) (see [START_REF] Davis | Markov Models and Optimization[END_REF]). Therefore, Following [START_REF] Crudu | Convergence of stochastic gene networks to hybrid piecewise deterministic processes[END_REF] we apply Gronwall's lemma to show that, when ϕ is Lipschitz with respect to x, so is u, with a Lipschitz constant less than κ ϕ 0,1 . Here, κ depends on the characteristics of the PDMP and the time horizon τ ; ϕ 0,1 is the sum of the supremum of ϕ with its Lipschitz constant. Moreover, it follows from the forward Kolmogorov equation that its time derivative is also bounded. where I is the identity matrix. Since H is bounded, it follows from Gronwall's lemma that Π is bounded and Π(τ j , t) -I ≤ κ|t -τ j | for some constant κ, where • is a norm on the space of matrices. For instance we may take the supremum of the modulus of the coefficients of a matrix. Hence, Π(τ j , τ j+1 ) -I - When ϕ is only Lipschitz, we obtain the same result by choosing a sequence of C 1 functions which converges uniformly to ϕ and whose Lipschitz constants do not exceed L. This finishes the proof of the first statement. It is classical that this implies that the left-handed side converges to zero when max j=0,...,M -1 |τ j+1 -τ j | → 0 for any uniformly continuous and bounded function ϕ. We conclude thanks to the Portemanteau theorem [START_REF] Billingsley | Convergence of Probability Measures[END_REF]. Remark 4. When ϕ is only borelian bounded but µ 0 has a C 1 and Lipschitz density f 0 , the proof can be adapted. We observe that the only point where the Lipschitz property of ϕ is used is to obtain (B.1) which is used to derive (B.2). To replace (B.2), we can instead write the explicit formula for E(ϕ(x(t), s(t)) and E(ϕ(x(t), s(t)). Since the flows Φ s , s ∈ S are C 1 diffeomorphisms that depend smoothly on time, we see that (B.2) still holds, but on the left-handed side, κ( ϕ 0,1 + 1) must be replaced by a constant which depends on f 0 0,1 and the supremum of ϕ.

E

) 5 :

 5 CUM := 0; I := 0 6: while CUM < V do 7:I := I + 1 8:

  end for 23: return (H 0 , . . . , H nt )

Figure 1 :

 1 Figure 1: Deterministic partition of time in the push forward method. Trajectories of the PDMP model (red line) having a jump of the discrete state inside [τ 2 , τ 3 ) are replaced by a trajectory (black line) having a jump of the discrete state at τ 3 .

17 :

 17 = V s1 (x1) with initial conditions x1(τ 1 ) = x0(τ 1 ) from t = τ 1 to t = τ 2 for s2 := 0 to 1 do 18: solve dx2 dt = V s2 (x2) with initial conditions x2(τ 2 ) = x1(τ 2 ) from t = τ 2 to t = τ 3 19:

is the bin containing x 23

 23

16 :

 16 for s1 := 1 to N s do 17: solve dx1 dt = V s1 (x1) with initial conditions x1(τ 1 ) = x0(τ 1 ) from t = τ 1 to t = τ 2 18:

ModelMonte-

  Carlo high sampling [min] Push-forward [s]

Figure 2 :

 2 Figure 2: Histograms of protein for the second gene, produced by the Monte-Carlo method (green lines) and by the Push-forward method (black lines) for the model M 1 . The green dotted line results from low sampling Monte-Carlo with similar execution time as the pushforward method, whereas the solid green line results from high sampling Monte-Carlo. The distances, defined by (29), are between low sampling and high sampling Monte-Carlo (d * ) and between push-forward and high sampling Monte-Carlo (d).

Figure 3 :

 3 Figure 3: Histograms of protein for the second gene, produced by the Monte-Carlo method (green lines) and by the Push-forward method (black lines) for the model M 2 . The green dotted line results from low sampling Monte-Carlo with similar execution time as the pushforward method, whereas the solid green line results from high sampling Monte-Carlo. The distances, defined by (29), are between low sampling and high sampling Monte-Carlo (d * ) and between push-forward and high sampling Monte-Carlo (d).

Figure 4 :

 4 Figure 4: This set of plots shows the comparison between the computed probability distribution for mRNA (P (r, t)) and protein (P (y, t)), using the push forward method (Algorithm 4), direct Monte-Carlo simulation (Algorithms 1 and 2) and numerical solution of the Liouvillemaster Equation (Algorithm 3). The model is for one gene in the slow switch regime ( = 0.5). The execution times of PF is 2.55 minutes in Python and 97.39 minutes for LME in MATLAB. The values d LM E and d P F are shown in each plot and measure the distances between the LME equation and MC simulation, PF and MC simulation, respectively, as given in (29). The value of parameters are: p 0 = 0.5, k 0 = 4, k 1 = 40, a = 1/5 and b = 4. For initial conditions we set: r(0) = y(0) = 0 and p 0 (0) = 1.

Figure 5 :

 5 Figure 5: This set of plots shows the comparison between the computed probability distribution for mRNA (P (r, t)) and protein (P (y, t)), using PF (Algorithm 4), MC simulation (Algorithms 1 and 2) and numerical solution LME (Algorithm 3). The model is for one gene producing in the fast switch regime ( = 5.5). The execution time of the PF is 5.12 minutes in Pyhton, and 21.64 minutes for LME in MATLAB. The values d LM E and d P F are shown in each plot and measure the distances between the LME equation and MC simulation, and PF and MC simulation, respectively, as given in (29). The value of parameters are: p 0 = 0.5, k 0 = 4, k 1 = 40, a = 1/5 and b = 4. For initial conditions we set: r(0) = y(0) = 0 and p 0 (0) = 1.

Figure 6 :

 6 Figure 6: This set of plots shows the comparison between the computed probability density for mRNA (P (r 2 , t)) and protein (P (y 2 , t)) associated to gene two, using PF (Algorithm 5), MC simulation (Algorithms 1 and 2) and numerical solution LME (Algorithm 3). The model is M 1 in the slow-slow regime ( = 0.5). The execution time of PF is 6.86 minutes in Python, and 97.39 minutes for LME in MATLAB. The values d LM E and d P F are shown in each plot and measure the distances between the LME equation and MC simulation, and PF and MC simulation, respectively, as given in (29). The value of parameters, for both genes are: p 0 = 0.5, k 0 = 4, k 1 = 40, a = 1/5 and b = 4. For initial conditions we set: r 1 (0) = y 1 (0) = r 2 (0) = y 2 (0) = 0 and p 0 (0) = 1.

Figure 7 :

 7 Figure 7: This set of plots shows the comparison between the computed probability density for mRNA (P (r 2 , t)) and protein (P (y 2 , t)) associated to gene two, using PF (described in Algorithm 5), MC simulation (Algorithms 1 and 2) and numerical solution LME (Algorithm 3). The model is M 1 in the fast-slow regime ( 1 = 5.5 and 2 = 0.5). The execution time of PF is 8.83 minutes in Python and 97.39 minutes for LME in MATLAB. The values d LM E and d P F are shown in each plot and measure the distances between the LME equation and MC simulation, and PF and MC simulation, respectively, as given in (29). The value of parameters, for both genes are: p 0 = 0.5, k 0 = 4, k 1 = 40, a = 1/5 and b = 4. For initial conditions we set: r 1 (0) = y 1 (0) = r 2 (0) = y 2 (0) = 0 and p 0 (0) = 1.

- 1 j=0E

 1 Hence, |u(t 1 , x 1 , s) -u(t 2 , x 2 , s)| ≤ κ ϕ 0,1 (|x 1 -x 2 | + |t 1 -t 2 |). (B.1)Now we rewrite the error as follows:E(ϕ(x(τ, ζ 0 ), s(τ, ζ 0 )) -E(ϕ(x(τ, ζ 0 ), s(τ, ζ 0 )) = E(u(τ, x 0 , s 0 ) -u(0, x(τ, ζ 0 ), s(τ, ζ 0 ))) = M (u(τ -τ j , x(τ j , ζ 0 ), s(τ j , ζ 0 )) -u(τ -τ j+1 , ζ 0 , x(τ j+1 , ζ 0 ), s(τ j+1 , ζ 0 ))On each interval [τ j , τ j+1 ] we have, by chain rule and the forward Kolmogorov equation, thatd dt [u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 ))] = -du dt (τ -t, x(t, ζ 0 ), s(t, ζ 0 )) + dx dt (t)∂ x u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 ) = -du dt (τ -t, x(t, ζ 0 ), s(t, ζ 0 )) + V s(τj ,ζ0) (x(t, ζ 0 ))∂ x u(τ -t, x(t, ζ 0 ), ζ 0 , s(τ j , ζ 0 ) = -r∈S (u(τ -t, x(t, ζ 0 ), r) -u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 )))H r,s(τj ,ζ0) (x(t, ζ 0 )).Hence,E u(τ -τ j , x(τ j , ζ 0 ), s(τ j ), ζ 0 ) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), s(τ j+1 , ζ 0 )) τ -t, x(t, ζ 0 ), r) -u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 )))H r,s(τj ,ζ0) (x(t, ζ 0 ))dt + E(u(τ -τ j+1 , x(τ j+1 , ζ 0 ), s(τ j , ζ 0 )) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), s(τ j+1 , ζ 0 ))) τ -t, x(t, ζ 0 ), r) -u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 )))H r,s(τj ,ζ0) (x(t, ζ 0 ))dt + E r∈S (u(τ -τ j+1 , x(τ j+1 , ζ 0 )), s(τ j , ζ 0 ) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), r) Π r,s(τ,ζ0j) (τ j , τ j+1 ).By eq. (20), we have that Π(τ j , τ j+1 ) = I + τj+1 τj H(x(t), ζ 0 )Π(τ j , t)dt,

2 ≤

 2 t, ζ 0 ))(I -Π(τ j , t))dt≤ κ|τ j+1 -τ j | 2 ,for another constant κ. In particular, for r = s(τ j ):|Π r,s(τj ) (τ j , τ j+1 ) -τj+1 τj H r,s(τj ,ζ0) (x(t, ζ 0 ))dt ≤ κ|τ j+1 -τ j | 2 , Moreover, |x(τ j+1,ζ0) -x(t, ζ 0 )| ≤ κ|τ j+1 -t|,since its derivative is V sj (x(t)), which is bounded. We have shown in (B.1) that u is Lipschitz with respect to x and t, thus|u(τ -t, x(t, ζ 0 ), r) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), r)| ≤ κ( ϕ 0,1 + 1)|τ j+1 -t|. (B.2)Now we have the following estimate|E(u(τ -τ j , x(τ j , ζ 0 ), s(τ j , ζ 0 )) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), s(τ j+1 , ζ 0 ))τ -t, x(t, ζ 0 ), r) -u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 )))H r,s(τj ,ζ0) (x(t, ζ 0 ))dt + E r∈S (u(τ -τ j+1 , x(τ j+1 , ζ 0 ), s(τ j , ζ 0 )) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), r))Π r,s(τj ,ζ0) (τ j , τ j+1 ) τ -t, x(t, ζ 0 ), r) -u(τ -t, x(t, ζ 0 ), s(τ j , ζ 0 ))) + (u(τ -τ j+1 , x(τ j+1 , ζ 0 )), s(τ j , ζ 0 ) -u(τ -τ j+1 , x(τ j+1 , ζ 0 ), r))H r,s(τj ),ζ0 (x(t, ζ 0 ))dt + κ|τ j+1 -τ j | 2 ≤ κ( ϕ 0,1 + 1)|τ j+1 -τ j | 2and so|E(ϕ(x(τ, ζ 0 ), s(τ, ζ 0 )) -E(ϕ(x(τ, ζ 0 ), s(τ, ζ 0 ))| ≤ M -1 j=0 κ|( ϕ 0,1 + 1)τ j+1 -τ j | Cκτ ( ϕ 0,1 + 1)/M.The final error estimate is obtained by taking ϕ depending only on x and integrating with respect to µ 0 : [0,1] ϕ(x)µ M t (dx) -B ϕ(x)µ t (dx) = [0,1] E(ϕ(x(τ, ζ 0 ))) -E(ϕ(x(τ, ζ 0 )))dµ 0 (ζ 0 ) ≤ Cκτ ( ϕ 0,1 + 1)/M

  1 , r 1 , . . . , y Ng , r Ng ), where p i,t := p i (t, y 1 , r 1 , . . . , y Ng , r Ng ), ∀i ∈ {1, . . . , 2N g }. 2. τ time step for computing the time-dependent distribution. 3. δy, δr bin size for the protein and mRNA distributions, respectively. Output: n t + 1 distributions P 0 , . . . , P nt representing the sums of the individual distributions in (y 1 , r 1 , y 2 , r 2 ) at times 0, τ, . . . , n t τ , respectively

1: for i := 1 to n t do 2:

Table 2 :

 2 Execution times (in minutes for all cases) for different methods to compute the probability distributions for one gene model and two genes model. Methods: PF = Pushforward, LME = Liouville-master equation

		-Carlo	PF	LME
		(high sampling)
	One gene slow	-	2.55 97.39
	One gene fast	-	5.12 21.64
	Two genes slow-slow	45	6.86 97.39
	Two genes fast-fast	74	10.85 21.64
	Two genes slow-fast	243	9.22 21.64
	Two genes fast-slow	249	8.83 97.39

  (ϕ(x(t), s(t)) = R×S E(ϕ(x(t, ζ 0 ), s(t, ζ 0 )))dµ 0 (ζ 0 ) = R×S u(t, x 0 , s 0 )dµ 0 (x 0 , s 0 ).It follows from Proposition A.1 in[START_REF] Crudu | Convergence of stochastic gene networks to hybrid piecewise deterministic processes[END_REF] that u satisfies the formula:

u(t, x, s) = ϕ(Φ s (t, x), s) exp -t 0 λ(Φ s (σ, x))dσ + s∈S t 0 u(t -σ, Φ s (σ, x), r, )H r,s (Φ s (σ, x)) exp -σ 0 λ(Φ s (τ, x))dτ dσ

where

Appendix B. Proof of the Theorem 1

For simplicity, we consider only the case N = 1, the proof generalizes immediately to N > 1 at the price of cumbersome notations. We also consider the case when the process x always belongs to [0, 1]. Thus, we may assume that V s , s ∈ S and H are bounded as well as their derivatives. We first consider the case of a Lipschitz function ϕ.

We denote by (x, s) the PDMP starting with the random initial data distributed along µ 0 and by (x(t, ζ 0 ), s(t, ζ 0 )) the PDMP starting from the deterministic initial data ζ 0 = (x 0 , s 0 ).

In this proof, we use the letter κ for a constant which depends on V s , s ∈ S, H and τ .

The push-forward algorithm can be rewritten as follows: starting from s(0) = s 0 , x(0) = x 0 , we set x(τ j+1 ) = Φ s(τj ) (τ j+1 -τ j , x(τ j )) and then choose s(τ j+1 ) randomly, it is equal to r ∈ S with probability Π s(τj ),r (τ j , τ j+1 ), the latter matrix being defined in [START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF].

In fact, we consider more generally the approximation of E(ϕ(x(τ ), s(τ )), for a function ϕ : R × S → R, by E(ϕ(x(τ ), s(τ )). The result is obtained when we choose ϕ independent on s ∈ S.

We first assume that ϕ is C 1 with respect to x. Let us denote by L its Lipschitz constant with respect to x:

The main tool for the proof is the generator of the process:

Lϕ(x, s) = V s (x)∂ x ϕ + r∈S (ϕ(x, r) -ϕ(x, s))H r,s (x),