Vincent Lemaire
email: vincent.lemaire@upmc.fr

Michèle Thieullen
email: michele.thieullen@upmc.fr

Nicolas Thomas
email: nicolas.thomas@upmc.fr

Thinning and Multilevel Monte Carlo for Piecewise Deterministic (Markov) Processes. Application to a stochastic Morris-Lecar model

Keywords: Piecewise Deterministic (Markov) Processes, Multilevel Monte Carlo, Thinning, Strong error estimate, Weak error expansion, Morris-Lecar model Mathematics Subject Classification: 65C05, 65C20, 60G55, 60J25, 68U20

In the first part of this paper we study approximations of trajectories of Piecewise Deterministic Processes (PDP) when the flow is not explicit by the thinning method. We also establish a strong error estimate for PDPs as well as a weak error expansion for Piecewise Deterministic Markov Processes (PDMP). These estimates are the building blocks of the Multilevel Monte Carlo (MLMC) method which we study in the second part. The coupling required by MLMC is based on the thinning procedure. In the third part we apply these results to a 2-dimensional Morris-Lecar model with stochastic ion channels. In the range of our simulations the MLMC estimator outperforms the classical Monte Carlo one.

Introduction

In this paper we are interested in the approximation of the trajectories of PDPs. We establish strong error estimates for a PDP and a weak error expansion for a PDMP. Then we study the application of the Multilevel Monte Carlo (MLMC) method in order to approximate expectations of functional of PDMPs. Our motivation comes from Neuroscience where the whole class of stochastic conductance-based neuron models can be interpreted as PDMPs. The response of a neuron to a stimulus, called neural coding, is considered as a relevant information to understand the functional properties of such excitable cells. Thus many quantities of interest such as mean first spike latency, mean interspike intervals and mean firing rate can be modelled as expectations of functionals of PDMPs.

PDPs have been introduced by Davis in [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF] as a general class of stochastic processes characterized by a deterministic evolution between two successive random times. In the case where the deterministic evolution part follows a family of Ordinary Differential Equations (ODEs) the corresponding PDP enjoys the Markov property and is called a PDMP. The distribution of a PDMP is thus determined by three parameters called the characteristics of the PDMP: a family of vector fields, a jump rate (intensity function) and a transition measure. We consider first a general PDP (x t ) which is not necessarily Markov on a finite time interval [0, T ] for which the flow is not explicitly solvable. Approximating its flows by the classical Euler scheme and using our previous work [START_REF] Lemaire | Exact simulation of the jump times of a class of Piecewise Deterministic Markov Processes[END_REF], we build a thinning algorithm which provides us with an exact simulation of an approximation of (x t ) that we denote (x t ). The process (x t ) is a PDP constructed by thinning of a homogeneous Poisson process which enjoys explicitly solvable flows. Actually this thinning construction provides a whole family of approximations indexed by the time step h > 0 of the Euler scheme. We prove that for any real valued smooth function F the following strong estimate holds

∃ V 1 > 0, V 2 > 0, E[|F (x T ) -F (x T )| 2 ] ≤ V 1 h + V 2 h 2 . ( 1
)
Moreover if (x t ) is a PDMP the following weak error expansion holds

∃ c 1 > 0, E[F (x T )] -E[F (x T )] = c 1 h + o(h 2 ). (2)
The estimate ( 1) is mainly based on the construction of the couple (x t , x t ) and on the fact that the Euler scheme is of order 1 this is why it is valid for a general PDP and its Euler scheme. On the contrary, the estimate (2) relies on properties which are specific to PDMPs such as the Feynman-Kac formula.

The MLMC method relies simultaneously on estimates (1) and ( 2) that is why we study its application to the PDMP framework instead of the more general PDP one. MLMC extends the classical Monte Carlo (MC) method which is a very general approach to estimate expectations using stochastic simulations. The complexity (i.e the number of operations necessary in the simulation) associated to a MC estimation can be prohibitive especially when the complexity of an individual random sample is very high. MLMC relies on repeated independent random samplings taken on different levels of accuracy which differs from the classical MC method. MLMC can then greatly reduces the complexity of the classical MC by performing most simulations with low accuracy but with low complexity and only few simulations with high accuracy at high complexity. MLMC have been introduced by S. Heinrich in [START_REF] Heinrich | Multilevel Monte Carlo methods[END_REF] and developed by M. Giles in [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF]. The MLMC estimator has been efficiently used in various fields of numerical probability such as SDEs [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF], Markov chains [START_REF] Anderson | Multilevel Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics[END_REF], [START_REF] Anderson | Complexity of Multilevel Monte Carlo tauleaping[END_REF], [START_REF] Glynn | Exact estimation for markov chain equilibrium expectations[END_REF], Lévy processes [START_REF] Ferreiro-Castilla | Multilevel Monte Carlo simulation for Lévy processes based on the Wiener-Hopf factorisation[END_REF], jump diffusions [START_REF] Xia | Multilevel Path Simulation for Jump-Diffusion SDEs[END_REF], [START_REF] Dereich | Multilevel Monte Carlo algorithms for Lévy-driven SDEs with Gaussian correction[END_REF], [START_REF] Dereich | A Multilevel Monte Carlo algorithm for Lévy-driven Stochastic Differential Equations[END_REF] or nested Monte Carlo [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF], [START_REF] Giorgi | Théorémes limites pour estimateurs Multilevel avec et sans poids[END_REF]. See [START_REF] Giles | Multilevel Monte Carlo methods[END_REF] for more references. To the best of our knowledge, application of MLMC to PDMPs has not been considered.

For the sake of clarity, we describe here the general improvement of MLMC. We are interested in the estimation of E[X] where X is a real valued square integrable random variable on a probability space (Ω, F , P). When X can be simulated exactly the classical MC estimator (1/N ) N k=1 X k with X k , k ≥ 1 independent random variables identically distributed as X, provides an unbiased estimator. The associated L 2 -error satisfies Y -E[X] 2 2 = Var(Y ) = 1 N Var(X). If we quantify the precision by the L 2 -error, then a user-prescribed precision ǫ 2 > 0 is achieved for N = O(ǫ -2 ) so that in this case the global complexity is of order O(ǫ -2 ). Assume now that X cannot be simulated exactly (or cannot be simulated at a reasonable cost) and that we can build a family of real valued random variables (X h , h > 0) on (Ω, F , P) which converges weakly and strongly to X as h → 0 in the following sense

∃ c 1 > 0, α > 0, E[X h ] -E[X] = c 1 h α + o(h 2α ), (3)
and

∃ V 1 > 0, β > 0, E[|X h -X| 2 ] ≤ V 1 h β . ( 4
)
Assume moreover that for h > 0 the random variable X h can be simulated at a reasonable complexity (the complexity increases as h → 0). The classical MC estimator now consists in a sequence of random variables

Y = 1 N N k=1 X k h , (5)
where X k h , k ≥ 1 are independent random variables identically distributed as X h . The bias and the variance of the estimator [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF] are respectively given by E

[Y ] -E[X] = E[X h ] -E[X] ≃ c 1 h α and Var(Y ) = 1
N Var(X h ). From the strong estimate (4) we have that Var(X h ) → Var(X) as h → 0 so that Var(X h ) is asymptotically a constant independent of h. If as above we quantify the precision by the L 2 -error and use that

Y -E[X] 2 2 = (E[Y ] -E[X]
) 2 + Var(Y ), we obtain that the estimator (5) achieves a user-prescribed precision ǫ 2 > 0 for h = O(ǫ 1/α ) and N = O(ǫ -2 ) so that the global complexity of the estimator is now O(ǫ -2-1 α ). The MLMC method takes advantage of the estimate (4) in order to reduce the global complexity. Let us fix L ≥ 2 and consider for l ∈ {1, . . . , L} a geometrically decreasing sequence (h l , 1 ≤ l ≤ L) where h l = h * M -(l-1) for fixed h * > 0 and M > 1. The indexes l are called the levels of the MLMC and the complexity of X h l increases as the level increases. Thanks to the weak expansion (3), the quantity E[X hL ] approximates E[X]. Using the linearity of the expectation the quantity E[X hL ] can be decomposed over the levels l ∈ {1, . . . , L} as follows

E[X hL ] = E[X h * ] + L l=2 E[X h l -X h l-1 ]. (6)
For each level l ∈ {1, . . . , L}, a classical MC estimator is used to approximate

E[X h l -X h l-1 ] and E[X h * ].
At each level, a number N l ≥ 1 of samples are required and the key point is that the random variables X h l and X h l-1 are assumed to be correlated in order to make the variance of X h l -X h l-1 small. Considering at each level l = 2, . . . , L independent couples (X h l , X h l-1 ) of correlated random variables, the MLMC estimator then reads

Y = 1 N 1 N1 k=1 X k h * + L l=2 1 N l N l k=1 (X k h l -X k h l-1 ), (7)
where (X k h * , k ≥ 1) is a sequence of independent and identically distributed random variables distributed as X h * and (X k h l , X k h l-1 ), k ≥ 1 for l = 2, . . . , L are independent sequences of independent copies of (X h l , X h l-1 ) and independent of (X k h * ). It is known, see [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] or [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF], that given a precision ǫ > 0 and provided that the family (X h , h > 0) satisfies the strong and weak error estimates (4) and (3), the multilevel estimator [START_REF] Dereich | Multilevel Monte Carlo algorithms for Lévy-driven SDEs with Gaussian correction[END_REF]

achieves a precision Y -E[X] 2 2 = ǫ 2 with a global complexity of order O(ǫ -2 ) if β > 1, O(ǫ -2 (log(ǫ)) 2 ) if β = 1 and O(ǫ -2-(1-β)/α ) if β < 1.
This complexity result shows the importance of the parameter β. Finally, let us mention that in the case β > 1 it possible to build an unbiased multilevel estimator, see [START_REF] Glynn | Unbiased Estimation with Square Root Convergence for SDE models[END_REF].

Estimates (1) and (2) suggest to investigate the use of the MLMC method in the PDMP framework with β = 1 and α = 1. Letting X = F (x T ) and X h = F (x T ) for h > 0 and F a smooth function, we define a MLMC estimator of E[F (x T )] just as in [START_REF] Dereich | Multilevel Monte Carlo algorithms for Lévy-driven SDEs with Gaussian correction[END_REF] (noted Y MLMC in the paper) where the processes involved at the level l are correlated by thinning. Since these processes are constructed using two different time steps, the probability of accepting a proposed jump time differs from one process to the other. Moreover the discrete components of the post-jump locations may also be different. This results in the presence of the term V 1 h in the estimate [START_REF] Anderson | Multilevel Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics[END_REF]. In order to improve the convergence rate (to increase the parameter β) in (1), we show that for a given PDMP (x t ) we have the following auxiliary representation

E[F (x T )] = E[F (x T ) RT ]. (8)
The PDMP (x t ) and its Euler scheme are such that their discrete components jump at the same times and in the same state. ( Rt ) is a process which depends on (x t , t ∈ [0, T ]). The representation ( 8) is inspired by the change of probability introduced in [START_REF] Xia | Multilevel Path Simulation for Jump-Diffusion SDEs[END_REF] and is actually valid for a general PDP (Proposition 2.2) so that

E[F (x T )] = E[F (x T ) RT ]
where (x t ) is the Euler scheme corresponding to (x t ) and ( Rt ) is a process which depends on (x t , t ∈ [0, T ]). Letting X = F (x T ) RT and X h = F (x T ) RT we define a second MLMC estimator (noted Ỹ MLMC ) where now the discrete components of the Euler schemes (x t ) involved at the level l always jump in the same states and at the same times. To sum up, the first MLMC estimator we consider (Y MLMC ) derives from [START_REF] Davis | Markov Models and Optimization[END_REF] where the corrective term at level

l is E[F (x h l T )-F (x h l-1 T
)] whereas the corrective

term of the second estimator ( Ỹ MLMC ) is E[F (x h l T ) Rh l T -F (x h l-1 T ) Rh l-1 T ].
For readability, we no longer write the dependence of the approximations on the time step. For the processes (F (x t ) Rt ) and (F (x t ) Rt ) we show the following strong estimate

∃ Ṽ1 > 0, E[|F (x T ) RT -F (x T ) RT | 2 ] ≤ Ṽ1 h 2 ,
so that we end up with β = 2 and the complexity goes from a

O(ǫ -2 (log(ǫ)) 2 ) to a O(ǫ -2 ).
As an application we consider the PDMP version of the 2-dimensional Morris-Lecar model, see [START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with application to neuron models[END_REF], which takes into account the precise description of the ionic channels and in which the flows are not explicit. Let us mention [START_REF] Benaïm | Quantitative ergodicity for some switched dynamical systems[END_REF] for the application of quantitative bounds for the long time behavior of PDMPs to a stochastic 3-dimensional Morris-Lecar model. The original deterministic Morris-Lecar model has been introduced in [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF] to account for various oscillating states in the barnacle giant muscle fiber. Because of its low dimension, this model is among the favourite conductance-based models in computational Neuroscience. Furthermore, this model is particularly interesting because it reproduces some of the main features of excitable cells response such as the shape, amplitude and threshold of the action potential, the refractory period. We compare the classical MC and the MLMC estimators on the 2-dimensional stochastic Morris-Lecar model to estimate the mean value of the membrane potential at fixed time. It turns out that in the range of our simulations the MLMC estimator outperforms the MC one. It suggests that MLMC estimators can be used successfully in the framework of PDMPs. As mentioned above, the quantities of interest such as mean first spike latency, mean interspike intervals and mean firing rate can be modelled as expectations of path-dependent functional of PDMPs. This setting can then be considered as a natural extension of this work.

The paper is organised as follows. In section 2, we construct a general PDP by thinning and we give a representation of its distribution in term of the thinning data (Proposition 1). In section 3, we establish strong error estimates (Theorems 1-2). In section 4, we establish a weak error expansion (Theorem 3). In section 5, we compare the efficiency of the classical and the multilevel Monte Carlo estimators on the 2-dimensional stochastic Morris-Lecar model.

Piecewise Deterministic Process by thinning 2.1 Construction

In this section we introduce the setting and recall some results on the thinning method from our previous paper [START_REF] Lemaire | Exact simulation of the jump times of a class of Piecewise Deterministic Markov Processes[END_REF]. Let E := Θ × R d where Θ is a finite or countable set and d ≥ 1. A piecewise deterministic process (PDP) is defined from the following characteristics

• a family of functions (Φ θ ) θ∈Θ such that Φ θ : R + × R d → R d for all θ ∈ Θ, • a measurable function λ : E →]0, +∞[, • a transition measure Q : E × B(E) → [0, 1].
We denote by x = (θ, ν) a generic element of E. We only consider PDPs with continuous ν-component so that for A ∈ B(Θ) and B ∈ B(R d ), we write

Q(x, A × B) = Q(x, A)δ ν (B). ( 9
)
If we write x = (θ x , ν x ), then it holds that

Q((θ x , Φ θx (t, ν x )), dθdν) = Q((θ x , Φ θx (t, ν x )), dθ)δ Φ θx (t,νx) (dν).
Our results do not depend on the dimension of the variable in R d so we restrict ourself to R (d = 1) for the readability. We work under the following assumption Assumption 2.1. There exists λ * < +∞ such that, for all x ∈ E, λ(x) ≤ λ * .

In [START_REF] Lemaire | Exact simulation of the jump times of a class of Piecewise Deterministic Markov Processes[END_REF] we considered a general upper bound λ * . In the present paper λ * is constant (see Assumption 2.1). Let (Ω, F , P) be a probability space on which we define 1. an homogeneous Poisson process (N * t , t ≥ 0) with intensity λ * (given in Assumption 2.1) whose successive jump times are denoted (T * k , k ≥ 1). We set T * 0 = 0. 2. two sequences of iid random variables with uniform distribution on [0, 1], (U k , k ≥ 1) and (V k , k ≥ 1) independent of each other and independent of (T * k , k ≥ 1). Given T > 0 we construct iteratively the sequence of jump times and post-jump locations (T n , (θ n , ν n ), n ≥ 0) of the E-valued PDP (x t , t ∈ [0, T ]) that we want to obtain in the end using its characteristics (Φ, λ, Q). Let (θ 0 , ν 0 ) ∈ E be fixed and let T 0 = 0. We construct T 1 by thinning of (T * k ), that is

T 1 := T * τ1 , (10)
where

τ 1 := inf {k > 0 : U k λ * ≤ λ(θ 0 , Φ θ0 (T * k , ν 0 ))} . ( 11
)
We denote by |Θ| the cardinal of Θ (which may be infinite) and we set Θ = {k 1 , . . . , k |Θ| }. For j ∈ {1, . . . , |Θ|} we introduce the functions a j defined on E by

a j (x) := j i=1 Q(x, {k i }), ∀x ∈ E. ( 12
)
By convention, we set a 0 := 0. We also introduce the function H defined by

H(x, u) := |Θ| i=1 k i 1 ai-1(x)<u≤ai(x) , ∀x ∈ E, ∀u ∈ [0, 1].
For all x ∈ E, H(x, .) is the inverse of the cumulative distribution function of Q(x, .) (see for example [START_REF] Devroye | Non-uniform random variate generation[END_REF]). Then, we construct (θ 1 , ν 1 ) from the uniform random variable V 1 and the function H as follows

(θ 1 , ν 1 ) = H (θ 0 , Φ θ0 (T * τ1 , ν 0 )), V 1 , φ θ0 (T * τ1 , ν 0 ) , = (H ((θ 0 , Φ θ0 (T 1 , ν 0 )), V 1 ) , φ θ0 (T 1 , ν 0 )) .
Thus, the distribution of (θ 1 , ν 1 ) given (τ 1 , (T

* k ) k≤τ1 ) is Q((θ 0 , Φ θ0 (T * τ1 , ν 0 )), .) or in view of (9), k∈Θ Q (θ 0 , Φ θ0 (T * τ1 , ν 0 )), {k} δ (k,φθ 0 (T * τ 1 ,ν0)) .
For n > 1, assume that

τ n-1 , (T * k ) k≤τn-1 , (θ n-1 , ν n-1 ) is constructed. Then, we construct T n by thinning of (T * k ) conditionally to τ n-1 , (T * k ) k≤τn-1 , (θ n-1 , ν n-1 )
, that is

T n := T * τn ,
where

τ n := inf k > τ n-1 : U k λ * ≤ λ(θ n-1 , Φ θn-1 (T * k -T * τn-1 , ν n-1 )) .
Then, we construct (θ n , ν n ) using the uniform random variable V n and the function H as follows

(θ n , ν n ) := H (θ n-1 , Φ θn-1 (T * τn -T * τn-1 , ν n-1 )), V n , Φ θn-1 (T * τn -T * τn-1 , ν n-1 ) = H (θ n-1 , Φ θn-1 (T n -T n-1 , ν n-1 )), V n , Φ θn-1 (T n -T n-1 , ν n-1 ) .
We define the PDP x t for all t ∈ [0, T ] from the process (T n , (θ n , ν n )) by

x t := (θ n , Φ θn (t -T n , ν n )) , t ∈ [T n , T n+1 [. (13)
Thus, x Tn = (θ n , ν n ) and x - Tn = (θ n-1 , ν n ). We also define the counting process associated to the jump times N t := n≥1 1 Tn≤t .

Approximation of a PDP

In applications we may not know explicitly the functions Φ θ . In this case, we use a numerical scheme Φ θ approximating Φ θ . In this paper, we consider schemes such that there exits positive constants C 1 and C 2 independent of h and θ such that

sup t∈[0,T ] |Φ θ (t, ν 1 ) -Φ θ (t, ν 2 )| ≤ e C1T |ν 1 -ν 2 | + C 2 h, ∀θ ∈ Θ, ∀(ν 1 , ν 2 ) ∈ R 2 . ( 14
)
To the family (Φ θ ) we can associate a PDP constructed as above that we denote (x t ). We emphasize that there is a positive probability that (x t ) and (x t ) jump at different times and/or in different states even if they are both constructed from the same data (N * t ), (U k ) and (V k ). However if the characteristics (Φ, λ, Q) of a PDP (x t ) are such that λ and Q depend only on θ, that is λ(x) = λ(θ) and Q(x, .) = Q(θ, .) for all x = (θ, ν) ∈ E, then its embedded Markov chain ( Tn , ( θn , νn ), n ≥ 0) is such that ( θn , n ≥ 0) is an autonomous Markov chain with kernel Q and ( Tn , n ≥ 0) is a counting process with intensity λt = n≥0 λ( θn )1 Tn≤t< Tn+1 . In particular, both ( θn ) and (τ n ) do not depend on Φ. The particular form of the characteristics λ and Q implies that the PDP (x t ) and its approximation (x t ) are correlated via the same process (τ n , θn ). In other words, these processes always jump exactly at the same times and their θ-component always jump in the same states. Such processes (x t ) are easier theoretically as well as numerically than the general case. They will be useful for us in the sequel.

The following lemma (which is important for several proofs below) gives a direct consequence of the estimate [START_REF] Glynn | Exact estimation for markov chain equilibrium expectations[END_REF]. Lemma 2.1. Let (Φ θ ) and (Φ θ ) satisfying [START_REF] Glynn | Exact estimation for markov chain equilibrium expectations[END_REF]. Let (t n , n ≥ 0) be an increasing sequence of non-negative real numbers with t 0 = 0 and let (α n , n ≥ 0) be a sequence of Θ-valued components. For a given ν ∈ R let us define iteratively the sequences (β n , n ≥ 0) and (β n , n ≥ 0) as follows

β n = Φ αn-1 (t n -t n-1 , β n-1 ), β 0 = ν,
and

β n = Φ αn-1 (t n -t n-1 , β n-1 ), β 0 = ν.
Then, for all n ≥ 1 we have

|β n -β n | ≤ e C1tn nC 2 h,
where C 1 and C 2 are positive constants independent of h.

Proof of Lemma 2.1. Let n ≥ 1. From the estimate (14), we have for all k ≤ n

β k -β k ≤ e C1(t k -t k-1 ) |β k-1 -β k-1 | + C 2 h,
and therefore

e -C1t k β k -β k ≤ e -C1t k-1 |β k-1 -β k-1 | + C 2 h.
By summing up these inequalities for 1 ≤ k ≤ n and since β 0 = β 0 we obtain

β n -β n ≤ e C1tn nC 2 h.

Application to the construction of a PDMP and its associated Euler scheme

In this section we define a PDMP and its associated Euler scheme from the construction of the section 2.1. For all θ ∈ Θ, we consider a family of vector fields (f θ , θ ∈ Θ) satisfying Assumption 2.2. For all θ ∈ Θ, the function f θ : R → R is bounded and Lipschitz with constant L independent of θ.

If we choose Φ θ = φ θ in the above construction where for all x = (θ, ν) ∈ E, we denote by (φ θ (t, ν), t ≥ 0) the unique solution of the ordinary differential equation (ODE)

dy(t) dt = f θ (y(t)) , y(0) = ν, ( 15
)
then the corresponding PDP is Markov since φ satisfies the semi-group property which reads φ θ (t + s, ν) = φ θ (t, φ θ (s, ν)) for all t, s ≥ 0 and for all (θ, ν) ∈ E. In this case, the process (x t ) is a piecewise deterministic Markov process (see [START_REF] Davis | Markov Models and Optimization[END_REF] or [START_REF] Jacobsen | Point Process Theory and Applications, Marked Point and Piecewise Deterministic Processes[END_REF]).

Let h > 0. We approximate the solution of (15) by the Euler scheme with time step h. First, we define the Euler subdivision of [0, +∞[ with time step h, noted (t i , i ≥ 0), by t i := ih.

Then, for all x = (θ, ν) ∈ E, we define the sequence (y i (x), i ≥ 0), the classical Euler scheme, iteratively by

y i+1 (x) = y i (x) + hf θ (y i (x)), y 0 (x) = ν,
to emphasize its dependence on the initial condition. Finally, for all x = (θ, ν) ∈ E, we set

φ θ (t, ν) := y i (x) + (t -t i )f θ (y i (x)), ∀t ∈ [t i , t i+1 ]. (16)
We construct the approximating process (x t ) as follows. Its continuous component starts from ν 0 at time 0 and follows the flow φ θ0 (t, ν 0 ) until the first jump time T 1 that we construct by ( 10) and ( 11) of section 2.1 where we replace

Φ θ0 (T * k , ν 0 ) by φ θ0 (T * k , ν 0 ). At time T 1 the continuous component of x T 1 is equal to φ θ0 (T 1 , ν 0 ) := ν 1 since
there is no jump in the continuous component. The discrete component jumps to θ 1 . We iterate this procedure with the new flow φ θ1 (t -T 1 , ν 1 ) until the next jump time T 2 given by ( 10) and ( 11) with φ θ1 (T * k -T 1 , ν 1 ) and so on. We proceed by iteration to construct (x t ) on [0, T ]. Consequently, the discretisation grid for (x t ) on the interval [0, T ] is random and is formed by the points T n + kh for n = 0, . . . , N T and k = 0, . . . , ⌊(T n+1 ∧ T -T n )/h⌋. This differs from the SDE case where the classical grid is fixed. By classical results of numerical analysis (see [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF] for example), the continuous Euler scheme ( 16) (also called Euler polygon) satisfies estimate [START_REF] Glynn | Exact estimation for markov chain equilibrium expectations[END_REF]. If we choose Φ θ = φ θ in the above construction then the corresponding PDP (x t ) is not Markov since the functions φ θ (., ν) do not satisfy the semi-group property (see [START_REF] Jacobsen | Point Process Theory and Applications, Marked Point and Piecewise Deterministic Processes[END_REF]).

Thinning representation for the marginal distribution of a PDP

The sequence (T n , (θ n , ν n ), n ≥ 0) is an R + × E-valued Markov chain with respect to its natural filtration F n and with kernel K defined by

K (t, θ, ν), dudjdz := 1 u≥t λ(θ, Φ θ (u -t, ν))e - u-t 0 λ(θ,Φ θ (s,ν))ds Q((θ, Φ θ (u -t, ν)), djdz)du . (17)
For n ≥ 0, the law of the random variable T n -T n-1 given F n-1 admits the density given for t ≥ 0 by

λ(θ n-1 , Φ θn-1 (t, ν n-1 ))e - t 0 λ(θn-1,Φ(s,νn-1))ds . ( 18
)
Classically the marginal distribution of x t is expressed using [START_REF] Giorgi | Théorémes limites pour estimateurs Multilevel avec et sans poids[END_REF], the intensity λ via (18) and the kernel K (see [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]). Indeed for fixed x 0 = x ∈ E and for any bounded measurable function g we can write,

E [g(x t )] = n≥0 E [g(θ n , Φ θn (t -T n , ν n ))1 Nt=n ] = n≥0 E g(θ n , Φ θn (t -T n , ν n ))1 Tn≤t E[1 Tn+1>t |F n ] = n≥0 E g(θ n , Φ θn (t -T n , ν n ))1 Tn≤t e - t-Tn 0 λ(θn,Φ θn (u,νn))du = n≥0 t 0 E g(θ, Φ θ (t -s, ν))e - t-s 0 λ(θ,Φ θ (u,ν))du K n ((0, x), dsdθdν)
where K 0 := δ and

K n = K • . . . • K n times, that is t 0 E K n ((0, x), dsdy) = t 0 E (R+×E) n-1 K((0, x), dt 1 dy 1 ) . . . K((t n-1 , y n-1 ), dsdy).
However since we have constructed (x t ) by thinning, we would prefer to express the distribution of x t using the upper bound λ * , the Poisson process (N * t , t ≥ 0) and the sequences (

U k , k ∈ N), (V k , k ∈ N). Proposition 2.1. Let (x t , t ∈ [0, T ]) be a PDP with characteristics (Φ, λ, Q) constructed in section 2.1 and let n ∈ N. Then E[g(x t )1 {Nt=n} ] = 1≤p1<p2....<pn≤m θ∈Θ E[Q(x - T * p n-1 , θ) g(θ, Φ θ (t -T * pn , ν n ))1 {τi=pi,1≤i≤n,N * t =m} m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , ν n )) λ * )].
The following proposition and its corollaries will be useful in section 3. In their statements (x t , t ∈ [0, T ]) and (x t , t ∈ [0, T ]) are PDPs constructed in section 2.1 using the same data (N * t ), (U k ), (V k ) and the same initial point x ∈ E but with different sets of characteristics. The following results are inspired by the change of probability introduced in [START_REF] Xia | Multilevel Path Simulation for Jump-Diffusion SDEs[END_REF] where the authors are interested in the application of the MLMC to jump-diffusion SDEs with state-dependent intensity. In our case, we need a change of probability which guarantees not only that the processes jump at the same times but also in the same states.

Proposition 2.2. Let us denote by

(Φ, λ, Q) ( resp. (Φ, λ, Q)) the characteristics of (x t ) (resp. (x t )).
Let us assume that λ and Q depend only on θ, that Q is always positive and 0 < λ(θ) < λ * for all θ ∈ Θ. For all integer n, let us define on the event { Ñt = n},

Zn = Q(x - T * τn , θn ) Q( θn-1 , θn ) 1 - λ( θn ) λ * N * t -τn -1 N * t q=τn+1 1 - λ( θn , Φ θn (T * q -T * τn , νn )) λ * , the product being equal to 1 if τn = N * t and for all 1 ≤ ℓ ≤ n -1, Zℓ = Q(x - T * τℓ , θℓ ) Q( θℓ-1 , θℓ ) λ( θℓ ) λ * 1 - λ( θℓ ) λ * τℓ+1 -τ ℓ -1 -1 λ( θℓ , Φ θℓ (T * τℓ+1 -T * τℓ , νℓ )) λ * τℓ+1 -1 q=τ ℓ +1 1 - λ( θℓ , Φ θℓ (T * q -T * τℓ , νℓ )) λ * , Z0 = λ( θ0 ) λ * 1 - λ( θ0 ) λ * τ1-1 -1 λ( θ0 , Φ θ0 (T * τ1 , ν0 )) λ * τ1-1 q=1 1 - λ( θ0 , Φ θ0 (T * q , ν0 )) λ * , Rn = Zn n-1 ℓ=0
Zℓ .

Then, for all n ≥ 0 we have

E[g(x t ) Rn 1 { Ñt=n} ] = E[g(x t ) 1 {Nt=n} ].
Corollary 2.1. Under the assumptions of Proposition 2.2, setting Rt = R Ñt , we have

E[g(x t ) Rt ] = E[g(x t )].
Remark 2.1. Proposition 2.2 looks like a Girsanov theorem (see [START_REF] Palmowski | A technique for exponential change of measure for Markov processes[END_REF]) however we do not use the martingale theory here.

Remark 2.2. We have chosen to state Proposition 2.2 with a PDP (x t ) whose intensity and transition measure only depend on θ for readability purposes. Actually the arguments of the proof are valid for non homogeneous intensity and transition measure of the form λ(x, t) and

Q((x, t), dy) for x = (θ, ν) ∈ E. A possible choice of such characteristics is λ(x, t) = λ(θ, Φθ (t, ν)) and Q((x, t), dy) = Q((θ, Φθ (t, ν))
, dy) for Φ a given function. This remark will be implemented in section 5.4.

Corollary 2.2. Let (Φ, λ, Q) (resp. ( Φ, λ, Q)) be the set of characteristics of (x t ) (resp. (x t )).
We assume that Q is always positive and that 0 < λ(x) < λ * for all x ∈ E. Let (µ n ) be the sequence defined by µ 0 = ν and

µ n = Φθn-1 (T n -T n-1 , µ n-1 ) for n ≥ 1.
For all integer n, let us define on the event {N t = n},

Zn = Q (θ n-1 , µ n ), θ n Q (θ n-1 , ν n ), θ n   N * t q=τn+1 1 - λ θ n , Φ θn (T * q -T * τn , ν n ) λ *   -1 N * t q=τn+1 1 - λ θ n , Φθn (T * q -T * τn , µ n ) λ * , the products being equal to 1 if τ n = N * t and for all 1 ≤ ℓ ≤ n -1, Zℓ = Q (θ ℓ-1 , µ ℓ ), θ ℓ Q (θ ℓ-1 , ν ℓ ), θ ℓ λ θ ℓ , Φ θ ℓ (T * τ ℓ+1 -T * τ ℓ , ν ℓ ) λ * τ ℓ+1 -1 q=τ ℓ +1 1 - λ θ ℓ , Φ θ ℓ (T * q -T * τ ℓ , ν ℓ ) λ * -1 λ θ ℓ , Φθ ℓ (T * τ ℓ+1 -T * τ ℓ , µ ℓ ) λ * τ ℓ+1 -1 q=τ ℓ +1 1 - λ θ ℓ , Φθ ℓ (T * q -T * τ ℓ , µ ℓ ) λ * , Z0 = λ θ 0 , Φ θ0 (T * τ1 , ν 0 ) λ * τ1-1 q=1 1 - λ θ 0 , Φ θ0 (T * q , ν 0 ) λ * -1 λ θ 0 , Φθ0 (T * τ1 , µ 0 ) λ * τ1-1 q=1 1 - λ θ 0 , Φθ0 (T * q , µ 0 ) λ * , Rn = Zn n-1 ℓ=0
Zℓ .

Then, for all n ≥ 0 we have

E[g θ n , Φθn (t -T n , µ n ) Rn 1 {Nt=n} ] = E[g(x t ) 1 { Ñt=n} ]. Proof of Proposition 2.1. It holds that {N t = n, τ i = p i , 1 ≤ i ≤ n} ⊂ {N * t ≥ p n }. Then E[g(x t )1 {Nt=n} ] = 1≤p1<p2<...<pn≤m E[g(x t )1 {Nt=n,τi=pi, 1≤i≤n,N * t =m} ]. The set {N t = n, τ i = p i , 1 ≤ i ≤ n, N * t = m} is equivalent to the following -N * t = m, -among the times T * ℓ , 1 ≤ ℓ ≤ m
exactly n are accepted by the thinning method they are the T * pi , 1 ≤ i ≤ n, all the others are rejected. We proceed by induction starting from the fact that all the T * q , p n + 1 ≤ q ≤ m are rejected which corresponds to the event

∀ p n + 1 ≤ q ≤ m, U q > λ(θ n , Φ θn (T * q -T * pn , ν n )) λ * .
The random variable

1 {τi=pi, 1≤i≤n} depends on (θ ℓ , ν ℓ , 1 ≤ ℓ ≤ n -1, T * i , 1 ≤ i ≤ p n , U j , 1 ≤ j ≤ p n ) where by construction ν ℓ = φ θ ℓ-1 (T * p ℓ -T * p ℓ-1 , ν ℓ-1 ), θ ℓ = H((θ ℓ-1 , ν ℓ ), V ℓ ) which implies that (θ ℓ , ν ℓ , 1 ≤ ℓ ≤ n -1) depend on (T * i , 1 ≤ i ≤ p n-1 , U j , 1 ≤ j ≤ p n-1 , V k , 1 ≤ k ≤ n -1). Thus V n is independent
* i , 1 ≤ i ≤ m + 1, U j , 1 ≤ j ≤ m, V k , 1 ≤ k ≤ n -1)
is therefore an expectation indexed by this vector as parameters. Since the law of

H(x, V n ) is Q(x, •) for all x ∈ E we obtain for p 1 < p 2 < ... < p n ≤ m, E[g(x t )1 {Nt=n,τi=pi, 1≤i≤n, N * t =m} ] = E[ θ∈Θ Q(x - T * p n-1 , θ) g(θ, Φ θ (t -T * pn , ν n )) F (θ, U j , 1 ≤ j ≤ m, T * ℓ , 1 ≤ ℓ ≤ m + 1, V k , 1 ≤ k ≤ n -1)], (19)
with

F (θ, U j , 1 ≤ j ≤ m, T * ℓ , 1 ≤ ℓ ≤ m + 1, V k , 1 ≤ k ≤ n -1) = 1 {N * t =m,τi=pi, 1≤i≤n} m q=pn+1 1 Uq> λ(θ,Φ θ (T * q -T * pn ,νn )) λ * .
In [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] the random variables (U q , p n + 1 ≤ q ≤ m) are independent of the vector (T * i ,

1 ≤ i ≤ m + 1, U j , 1 ≤ j ≤ p n , V k , 1 ≤ k ≤ n -1).
Conditioning by this vector we obtain

E[g(x t )1 {Nt=n,τi=pi, 1≤i≤n, N * t =m} ] = θ∈Θ E[Q(x - T * p n-1 , θ) g(θ, Φ θ (t -T * pn , ν n ))1 {N * t =m,τi=pi, 1≤i≤n} m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , ν n )) λ * )].
We can iterate on the latter form by first conditioning V n-1 by all the other r.v. and then conditioning (U q , p n-1 + 1 ≤ q ≤ p n ) by all the remaining ones and so on. However the terms that appear do not have the same structure since the U q correspond to a rejection for p n-1 + 1 ≤ q ≤ p n -1 whereas U pn corresponds to an acceptation. So that the next step yields

E[g(x t )1 {Nt=n,τi=pi, 1≤i≤n, N * t =m} ] = α∈Θ θ∈Θ E[Q(x - T * p n-2 , α)Q ((α, ν n ), θ) g(θ, Φ θ (t -T * pn , ν n ))1 {N * t =m,τi=pi, 1≤i≤n-1} λ(α, Φ α (T * pn -T * pn-1 , ν n-1 )) λ * pn-1 q=pn-1+1 (1 - λ(α, Φ α (T * q -T * pn-1 , ν n-1 )) λ * ) m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , ν n )) λ * )], (20)
where we write ν n for simplicity keeping in mind that

ν n = Φ θn-1 (T * pn -T * pn-1 , ν n-1 ) = Φ θn-1 (T * pn - T * pn-1 , Φ θn-2 (T * pn-1 -T * pn-2 , ν n-2 )) = Φ α (T * pn -T * pn-1 , Φ θn-2 (T * pn-1 -T * pn-2 , ν n-2 )). Moreover the previous arguments apply to E(g(x t )f (θ i , ν i , 1 ≤ i ≤ n -1, θ n , ν n , T * k , 1 ≤ k ≤ m) 1 {Nt=n,τi=pi, 1≤i≤n, N * t =m} ) and provide E[g(x t )f (θ i , ν i , 1 ≤ i ≤ n -1, θ n , ν n , T * k , 1 ≤ k ≤ m) 1 {Nt=n,τi=pi, 1≤i≤n, N * t =m} ] = θ∈Θ E[Q(x - T * p n-1 , θ)g(θ, Φ θ (t -T * pn , ν n ))f (θ i , ν i , 1 ≤ i ≤ n -1, θ, ν n , T * k , 1 ≤ k ≤ m) 1 {N * t =m,τi=pi, 1≤i≤n} m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , ν n )) λ * )]. (21)
We prove below Proposition 2.2. The other statements can be proved analogously.

Proof of Proposition 2.2. By assumption the (jump) characteristics ( λ, Q) of (x t ) depend only on θ. Let p 1 < p 2 < ... < p n ≤ m. Applying the same arguments as in [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF] to (x t ) and using the definitions of Zℓ , 0 ≤ ℓ ≤ n and Rn we obtain,

E[g(x t ) Rn 1 { Ñt=n,τi=pi,1≤i≤n,N * t =m} ] = θ∈Θ E[ Q( θn-1 , θ) g(θ, Φ θ (t -T * pn , νn )) Zn n-1 ℓ=0 Zℓ 1 {N * t =m,τi=pi, 1≤i≤n} ] (1 - λ(θ) λ * ) m-pn = θ∈Θ E[ Q( θn-1 , θ) g(θ, Φ θ (t -T * pn , νn )) n-1 ℓ=0 Zℓ 1 {N * t =m,τi=pi, 1≤i≤n} (1 - λ(θ) λ * ) m-pn (1 - λ(θ) λ * ) m-pn -1 Q(x - T * p n-1 , θ) Q( θn-1 , θ) m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , νn )) λ * )] = θ∈Θ E[Q(x - T * p n-1 , θ) g(θ, Φ θ (t -T * pn , νn )) Zn-1 n-2 ℓ=0 Zℓ 1 {N * t =m,τi=pi, 1≤i≤n} m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , νn )) λ * )].
We iterate the previous argument based on the use of ( 21) and we use the definition of Zn-1 to obtain

E[g(x t ) Rn 1 { Ñt=n,τi=pi, 1≤i≤n, N * t =m} ] = α∈Θ θ∈Θ E[Q(x - T * p n-2 , α)Q((α, νn ), θ) g(θ, Φ θ (t -T * pn , νn )) n-2 ℓ=0 Zℓ 1 {N * t =m,τi=pi, 1≤i≤n-1} m q=pn+1 (1 - λ(θ, Φ θ (T * q -T * pn , νn )) λ * ) λ(α, Φ α (T * pn -T * pn-1 , νn-1 )) λ * pn-1 q=pn-1+1 (1 - λ(α, Φ α (T * q -T * pn-1 , νn-1 )) λ * )],
where for short νn = φ α (T * pn -T * pn-1 , νn-1 ) and νn-1 = φ θn-2 (T * pn-1 -T * pn-2 , νn-2 ). Comparing the latter expression to [START_REF] Jacobsen | Point Process Theory and Applications, Marked Point and Piecewise Deterministic Processes[END_REF] and using an induction we conclude that

E[g(x t ) Rn 1 { Ñt=n,τi=pi, 1≤i≤n, N * t =m} ] = E[g(x t ) 1 {Nt=n,τi=pi,1≤i≤n,N * t =m} ].
It remains to sum up on p i , 1 ≤ i ≤ n and m.

Strong error estimates

In this section we are interested in strong error estimates. Below, we state the main assumptions and theorems of this section, the proofs are given in sections 3.2, 3.3 respectively. Assumption 3.1. For all θ ∈ Θ and for all A ∈ B(Θ), the functions ν → λ(θ, ν) and ν → Q((θ, ν), A) are Lipschitz with constants L λ > 0, L Q > 0 respectively independent of θ. Theorem 3.1. Let Φ θ and Φ θ satisfying [START_REF] Glynn | Exact estimation for markov chain equilibrium expectations[END_REF] and let (x t , t ∈ [0, T ]) and (x t , t ∈ [0, T ]) be the corresponding PDPs constructed in section 2.1 with x 0 = x 0 = x for some x ∈ E. Assume that Θ is finite and that λ and Q satisfy Assumption 3.1. Then, for all bounded functions F : E → R such that for all θ ∈ Θ the function ν → F (θ, ν) is L F -Lipschitz where L F is positive and independent of θ, there exists constants V 1 > 0 and V 2 > 0 independent of the time step h such that

E |F (x T ) -F (x T )| 2 ≤ V 1 h + V 2 h 2 . Remark 3.1. When the numerical scheme Φ θ is of order p ≥ 1, which means sup t∈[0,T ] |Φ θ (t, ν 1 )- Φ θ (t, ν 2 )| ≤ e C1T |ν 1 -ν 2 | + C 2 h p we have E |F (x T ) -F (x T )| 2 ≤ V 1 h p + V 2 h 2p .

Assumption 3.2.

There exist positive constants ρ, λmin , λmax such that for all (i, j)

∈ Θ 2 , ρ ≤ Q(i,
F : E → R such that for all θ ∈ Θ the function ν → F (θ, ν) is L F -Lipschitz (L F > 0), there exists a positive constant Ṽ1 independent of the time step h such that E |F (x T ) RT -F (x T ) RT | 2 ≤ Ṽ1 h 2 ,
where RT has been defined in Corollary 2.1.

We now introduce the random variable τ † which will play an important role in the strong error estimate of Theorem 3.1 as well as in the identification of the coefficient c 1 in the weak error expansion in section 4 (see the proof of Theorem 4.1 in section 4.2).

Definition 3.1. Let us define τ

† := inf k > 0 : (τ k , θ k ) = (τ k , θ k ) .
The random variable τ † enables us to partition the trajectories of the couple (x t , x t ) in a sense that we precise now. Consider the event

{min(T τ † , T τ † ) > T } = N T = N T , (T 1 , θ 1 ) = (T 1 , θ 1 ), . . . , (T NT , θ NT ) = (T N T , θ N T ) , ( 22
)
where (T n ) and (T n ) denote the sequences of jump times of (x t ) and (x t ). On this event {min(T τ † , T τ † ) > T } the trajectories of the discrete time processes (T n , θ n ) and (T n , θ n ) are equal for all n such that T n ∈ [0, T ] (or equivalently

T n ∈ [0, T ]). Moreover the complement i.e {min(T τ † , T τ † ) ≤ T } contains the trajectories for which (T n , θ n ) and (T n , θ n ) differ on [0, T ] (there exits n ≤ N T ∨ N T such that T n = T n or θ n = θ n ).

Preliminary lemmas

In this section we start with two lemmas which will be useful to prove Theorems 3.2 and 3.3.

Lemma 3.1. Let K be a finite set. We denote by |K| the cardinal of K and for i = 1, . . . , |K| we denote by k i its elements. Let (p i , 1 ≤ i ≤ |K|) and (p i , 1 ≤ i ≤ |K|) be two probabilities on K. Let a j := j i=1 p i and a j := j i=1 p i for all j ∈ {1, . . . , |K|}. By convention, we set a 0 = a 0 := 0. Let X and X be two K-valued random variables defined by

X := G(U ), X := G(U ),
where

U ∼ U([0, 1]), G(u) = |K| j=1 k j 1 aj-1<u≤aj and G(u) = |K| j=1 k j 1 aj-1<u≤aj for all u ∈ [0, 1]. Then, we have P(X = X) ≤ |K|-1 j=1 |a j -a j |.
Proof of Lemma 3.1. By definition of X and X and since the intervals ]a j-1 , a j ]∩]a j-1 , a j ] are disjoints for j = 1, . . . , K, we have

P(X = X) = |K| j=1 P U ∈ ]a j-1 , a j ]∩]a j-1 , a j ] .
Moreover, for all 1 ≤ j ≤ |K|, we have

P U ∈ ]a j-1 , a j ]∩]a j-1 , a j ] = 0 if ]a j-1 , a j ]∩]a j-1 , a j ] = ∅, a j ∧ a j -a j-1 ∨ a j-1 if ]a j-1 , a j ]∩]a j-1 , a j ] = ∅.
Thus, denoting by x + := max(x, 0) the positive part of x ∈ R and using that x + ≥ x, we obtain

P(X = X) ≥ |K| j=1 (a j ∧ a j -a j-1 ∨ a j-1 ).
Adding and subtracting a j ∨ a j in the the above sum yields

P(X = X) ≥ |K| j=1 (a j ∨ a j -a j-1 ∨ a j-1 ) + |K| j=1 (a j ∧ a j -a j ∨ a j ).
The first sum above is a telescopic sum. Since a |K| = a |K| = 1 and a 0 = a 0 = 0, we have

P(X = X) ≥ 1 - |K|-1 j=1 |a j -a j |.
Lemma 3.2. Let (a n , n ≥ 1) and (b n , n ≥ 1) be two real-valued sequences. For all n ≥ 1, we have

n i=1 a i - n i=1 b i = n i=1 (a i -b i ) n j=i+1 a j i-1 j=1 b j
Proof of Lemma 3.2. By induction.

Proof of Theorem 3.1

First, we write

E |F (x T ) -F (x T )| 2 = E 1 min(T τ † ,T τ † )≤T |F (x T ) -F (x T )| 2 + E 1 min(T τ † ,T τ † )>T |F (x T ) -F (x T )| 2 =: P + D,
where τ † is defined in Definition 3.1. The order of the term P is the order of the probability that the discrete processes (T n , θ n ) and (T n , θ n ) differ on [0, T ]. The order of the term D is given by the order of the Euler scheme squared because the discrete processes (T n , θ n ) and (T n , θ n ) are equal on [0, T ]. In the following we prove that P = O(h) and that D = O(h 2 ).

Step 1: estimation of P . The function F being bounded we have P ≤ 4M 2

F P min(T τ † , T τ † ) ≤ T where M F > 0. Moreover, for k ≥ 1, τ † = k = τ † > k -1 (τ k , θ k ) = τ k , θ k . Hence P min(T τ † , T τ † ) ≤ T = k≥1 E 1 min(T k ,T k )≤T 1 τ † =k = k≥1 E 1 min(T k ,T k )≤T 1 τ † >k-1 1 (τ k ,θ k ) =(τ k ,θ k ) ≤ k≥1 J k + 2I k
where

J k := E 1 min(T k ,T k )≤T 1 τ † >k-1 1 τ k =τ k 1 θ k =θ k , I k := E 1 min(T k ,T k )≤T 1 τ † >k-1 1 τ k =τ k . ( 23
) We start with J k . First note that, for k ≥ 1, {τ k = τ k } = {T k = T k } and that on the event {T k = T k }, we have min(T k , T k ) = T k , so that J k = E 1 T k ≤T 1 τ † >k-1 1 τ k =τ k 1 θ k =θ k . We emphasize that it makes no difference in the rest of the proof if we choose min(T k , T k ) = T k . Since {τ † > k -1} = k-1 i=0 {(τ i , θ i ) = (τ i , θ i )}, we can rewrite J k as follows 1≤p 1 <...<p k α 1 ,...,α k-1 ∈Θ E[1 {τi=τ i =pi,1≤i≤k} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T 1 θ k =θ k ]. (24)
By construction we have

θ k = H((θ k-1 , ν k ), V k ) and θ k = H((θ k-1 , ν k ), V k ). The random variable 1 {τi=τ i=pi ,1≤i≤k} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k
≤T depends on the vector (U i ,

1 ≤ i ≤ p k , T * j , 1 ≤ j ≤ p k , V q , 1 ≤ q ≤ k -1) which is independent of V k .
Conditioning by this vector in [START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF] and applying Lemma 3.1 yields

E[1 {τi=τ i =pi,1≤i≤k} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T 1 θ k =θ k ] ≤ E   1 {τi=τ i=pi ,1≤i≤k} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T |Θ|-1 j=1 |a j (α k-1 , ν k ) -a j (α k-1 , ν k )|   .
From the definition of a j (see [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF]), the triangle inequality and since

Q is L Q -Lipschitz, we have |Θ|-1 j=1 |a j (α k-1 , ν k ) -a j (α k-1 , ν k )| ≤ (|Θ|-1)|Θ| 2 L Q |ν k -ν k |.
Since we are on the event

{τ i = τ i = p i , 1 ≤ i ≤ k} {θ i = θ i = α i , 1 ≤ i ≤ k -1}, the application of Lemma 2.1 yields |ν k -ν k | ≤ e LT * p k kCh. Thus J k ≤ C 1 hE[1 T k ≤T k] where C 1 is a constant independent of h. Moreover, k≥1 1 T k ≤T k = NT k=1 k ≤ N 2 T and E[N 2 T ] ≤ E[(N * T ) 2 ] < +∞ so that k≥1 J k = O(h).
From the definition of I k (see ( 23)), we can write

I k = E 1 min(T k ,T k )≤T 1 τ † >k-1 (1 τ k <τ k + 1 τ k >τ k ) = E 1 T k ≤T 1 τ † >k-1 1 τ k <τ k + E 1 T k ≤T 1 τ † >k-1 1 τ k >τ k =: I (1) k + I (2)
k .

The second equality above follows since {τ

k < τ k } = {T k < T k } and {τ k > τ k } = {T k > T k }.
We only treat the term I k can be treated similarly by interchanging the role of (τ k , T k ) and (τ k , T k ). Just as in the previous case, we can rewrite I (1)

k as follows 1≤p 1 <...<p k α 1 ,...,α k-1 ∈Θ E[1 {τi=τ i =pi,1≤i≤k-1} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T 1 τ k =p k 1 p k <τ k ]. (25)
In [START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with application to neuron models[END_REF] we have

{τ k = p k } ∩ {p k < τ k } ⊆ {λ(α k-1 , Φ α k-1 (T * p k -T * p k-1 , ν k-1 )) < U p k λ * ≤ λ(α k-1 , Φ α k-1 (T * p k -T * p k-1 , ν k-1 ))}. The random variable 1 {τi=τ i=pi ,1≤i≤k-1} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T depends on (U i , 1 ≤ i ≤ p k-1 , T * j , 1 ≤ j ≤ p k , V q , 1 ≤ q ≤ k -1) which is independent of U p k .
Conditioning by this vector in [START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with application to neuron models[END_REF] yields

E[1 {τi=τ i=pi ,1≤i≤k-1} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T 1 τ k =p k 1 p k <τ k ] ≤ E[1 {τi=τ i=pi ,1≤i≤k-1} 1 {θi=θi=αi,1≤i≤k-1} 1 T * p k ≤T |λ(α k-1 , Φ α k-1 (T * p k -T * p k-1 , ν k-1 )) -λ(α k-1 , Φ α k-1 (T * p k -T * p k-1 , ν k-1 ))|]
Using the Lipschitz continuity of λ then Lemma 2.1 we get that

I (1) k ≤ C 2 hE[1 T k ≤T k] where C 2 is a constant independent of h. Concerning the term I (2)
k , we will end with the estimate

I (2) k ≤ C 2 hE[1 T k ≤T k].
We conclude in the same way as in the estimation of J k above that

k≥1 I k = O(h).
Step 2: estimation of D. Note that for n ≥ 0 we have

{N T = n} ∩{min(T τ † , T τ † ) > T } = {N T = n} ∩ {N T = n} ∩ {τ † > n},
where we can interchange the role of {N T = n} and {N T = n}. Thus, using the partition {N T = n, n ≥ 0}, we have

D = n≥0 E 1 NT =n 1 NT =n 1 τ † >n F (θ n , Φ θn (T -T n , ν n )) -F (θ n , Φ θn (T -T n , ν n )) 2
The application of the Lipschitz continuity of F and of Lemma 2.1 yields

F (θ n , Φ θn (T -T n , ν n )) -F (θ n , Φ θn (T -T n , ν n )) ≤ L F e LT (n + 1)Ch.
Then, we have

D ≤ C 3 h 2 n≥0 E 1 NT =n (n + 1) 2 where C 3 is a constant independent of h. Since n≥0 E 1 NT =n (n + 1) 2 = E[(N T + 1) 2 ] ≤ E[(N * T + 1) 2 ] < +∞, we conclude that D = O(h 2 ).

Proof of Theorem 3.2

First we reorder the terms in RT . We write RT = QT ST HT where

QT = ÑT l=1 Q(x - T * τl , θl ) Q( θl-1 , θl ) , ( 26
) ST = ÑT l=1 λ( θl-1 , Φ θl-1 (T * τl -T * τl-1 , νl-1 )) λ * τl k=τ l-1 +1 (1 - λ( θl-1 , Φ θl-1 (T * k -T * τl-1 , νl-1 )) λ *
) ( 27)

N * T l=τ ÑT +1 (1 - λ( θ ÑT , Φ θ ÑT (T * l -T * τ ÑT , ν ÑT )) λ * ), HT = ÑT l=1 λ( θl-1 ) λ * (1 - λ( θl-1 ) λ * ) τl -τ l-1 -1 -1 (1 - λ( θ ÑT ) λ * ) N * T -τ ÑT -1 . ( 28
)
Likewise we reorder the terms in RT writing RT = QT ST HT where QT and ST are defined as ( 26) and ( 27) replacing x and Φ by x and Φ. Since the processes ( θn ) and (τ n ) do not depend on Φ or Φ, the term H is the same in R and R . To prove Theorem 3.2, let us decompose the problem and write

|F (x T ) RT -F (x T ) RT | = |(F (x T ) -F (x T )) RT + ( RT -RT )F (x T )| ≤ |F (x T ) -F (x T )|| RT | + | RT -RT ||F (x T )|, so that E |F (x T ) RT -F (x T ) RT | 2 ≤ 2E |F (x T ) -F (x T )| 2 | RT | 2 + 2E | RT -RT | 2 |F (x T )| 2 =: 2D + 2C.
In the following we show that C = O(h 2 ) and that D = O(h 2 ).

Step 1: estimation of C. The function F being bounded we have

C ≤ M 2 F E | RT -RT | 2 where M F is a positive constant. Moreover, for all θ ∈ Θ, we have (1 -λ(θ)/λ * ) -1 ≤ (1 -λmax /λ * ) -1 and ( λ(θ)/λ * ) -1 ≤ ( λmin /λ * ) -1 . Thus, HT ≤ λmin λ * (1 -λmax λ * ) -N *
T and using the definition of R and R (see [START_REF] Palmowski | A technique for exponential change of measure for Markov processes[END_REF], ( 27) and ( 28)) we can write

| RT -RT | ≤ λmin λ * (1 - λmax λ * ) -N * T | QT -QT | ST + | ST -ST | QT .
We set J = | QT -QT | ST and I = | ST -ST | QT . To provide the desired estimate for C, we proceed as follows. First, we work ω by ω to determine (random) bounds for J and I from which we deduce a (random) bound for | RT -RT |. Finally, we take the expectation. We start with I. For all (θ, ν) ∈ E and for all t ≥ 0 we have, from Assumption 2.1, that 1λ(θ, Φ θ (t, ν))/λ * ≤ 1 and λ(θ, Φ θ (t, ν))/λ * ≤ 1. Then, using Lemma 3.2 (twice) we have

| ST -ST | ≤ 1 λ * ÑT +1 l=1 τl ∧N * T k=τ l-1 +1 |λ( θl-1 , Φ θl-1 (T * k -T * τl-1 , νl-1 )) -λ( θl-1 , Φ θl-1 (T * k -T * τl-1 , νl-1 ))|.
Using the Lipschitz continuity of λ and Lemma 2.1, we find that, for all l = 1, . . . , ÑT + 1 and

k = τl-1 + 1, . . . , τl ∧ N * T , |λ( θl-1 , Φ θl-1 (T * k -T * τl-1 , νl-1 )) -λ( θl-1 , Φ θl-1 (T * k -T * τl-1 , νl-1 ))| ≤ e LT Chl.
Moreover, for all l = 1, . . . , ÑT + 1 we have τl

∧ N * T -τl-1 ≤ N * T so that | ST -ST | ≤ N * T (N * T + 1) 2 C 1 h where C 1 is a positive constant independent of h. Finally, since QT ≤ ρ -N * T we have I ≤ ρ -N * T N * T (N * T + 1) 2 C 1 h. ( 29
)
Now, consider J. Note that from Assumption 2.1 we have ST ≤ 1. We use the same type of arguments as for I. That is, we successively use Lemma 3.2, the Lipschitz continuity of Q and Lemma 2.1 to obtain

J ≤ ρ -N * T (N * T ) 2 C 2 h, ( 30
)
where C 2 is a positive constant independent of h. Then, we derive from the previous estimates (29) and (30) that

| RT -RT | ≤ Ξ 1 (N * T )C 3 h, where Ξ 1 (n) = ρ λmin λ * (1 -λmax λ * ) -n n(n + 1) 2 and C 3 = max(C 1 , C 2 ). Finally, we have E[| RT - RT | 2 ] ≤ C 3 h 2 E[Ξ 1 (N * T ) 2 ]. Since E[Ξ 1 (N * T ) 2 ] < +∞ we conclude that C = O(h 2 ).
Step 2: estimation of D.

Recall that xT = ( θ ÑT , Φ θ ÑT (T -T ÑT , ν ÑT )) and xT = ( θ ÑT , Φ θ ÑT (T - T ÑT , ν ÑT )).
Then, using the Lipschitz continuity of F , Lemma 2.1 and since ÑT ≤ N * T we get

|F (x T ) -F (x T )| ≤ L F e LT ( ÑT + 1)Ch ≤ L F e LT (N * T + 1)Ch. Moreover, | RT | ≤ ρ λmin λ * (1 -λmax λ * ) -N * T so that D ≤ C 4 h 2 E[Ξ 2 (N * T ) 2 ] where C 4 is a positive constant independent of h and Ξ 2 (n) = (n + 1) ρ λmin λ * (1 -λmax λ * ) -n . Since E[Ξ 2 (N * T ) 2 ] < +∞ we conclude that D = O(h 2 ).

Weak error expansion

In this section we are interested in a weak error expansion for the PDMP (x t ) of section 2.3 and its associated Euler scheme (x t ). First of all, we recall from [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF] that the generator A of the process (t, x t ) which acts on functions g defined on R + × E is given by

Ag(t, x) = ∂ t g(t, x) + f (x)∂ ν g(t, x) + λ(x) E (g(t, y) -g(t, x))Q(x, dy), ( 31
)
where for notational convenience we have set ∂ ν g(t, x) := ∂g ∂ν (t, θ, ν), ∂ t g(t, x) := ∂g ∂t (t, x) and f (x) = f θ (ν) for all x = (θ, ν) ∈ E. Below, we state the assumptions and the main theorem of this section. Its proof which is inspired by [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] (see also [START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF] or [START_REF] Graham | Stochastic simulation and Monte Carlo[END_REF]) is delayed in section 4.2. Assumption 4.1. For all θ ∈ Θ and for all A ∈ B(Θ), the functions ν → Q ((θ, ν), A), ν → λ (θ, ν) and ν → f θ (ν) are bounded and twice continuously differentiable with bounded derivatives.

Assumption 4.2. The solution u of the integro differential equation

Au(t, x) = 0, (t, x) ∈ [0, T [×E, u(T, x) = F (x), x ∈ E, ( 32
)
with F : E → R a bounded function and A given by ( 31) is such that for all θ ∈ Θ, the function

(t, ν) → u(t,
E[F (x T )] -E[F (x T )] = hc 1 + O(h 2 ). ( 33
)
Remark 4.1. If (x t ) is a PDMP whose characteristics λ, Q satisfy the assumptions of Proposition 2.2 and (x t ) is its approximation we deduce from Theorem 4.1 that

E[F (x T ) RT ] -E[F (x T ) RT ] = hc 1 + O(h 2 ). ( 34
)

Further results on PDMPs: Itô and Feynman-Kac formulas

Definition 4.1. Let us define the following operators which act on functions g defined on R + ×E.

T g(t, x)

:= ∂ t g(t, x) + f (x)∂ ν g(t, x), Sg(t, x) := λ(x) E (g(t, y) -g(t, x))Q(x,

dy).

From Definition 4.1, the generator A defined by (31) reads Ag(t, x) = T g(t, x) + Sg(t, x). We introduce the random counting measure p associated to the PDMP (x t ) defined by p([0, t]×A) := n≥1 1 Tn≤t 1 Yn∈A for t ∈ [0, T ] and for A ∈ B(E). The compensator of p, noted p ′ , is given from [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF] by

p ′ ([0, t] × A) = t 0 λ(x s )Q(x s , A)ds.
Hence, q := pp ′ is a martingale with respect to the filtration generated by p noted (F p t ) t∈[0,T ] . Similarly, we introduce p, p ′ , q and (F p t ) t∈[0,T ] to be the same objects as above but corresponding to the approximation (x t ). The fact that p ′ is the compensator of p and that q is a martingale derives from arguments of the marked point processes theory, see [START_REF] Brémaud | Point Processes and Queues, Martingale Dynamics[END_REF].

Definition 4.2. Let us define the following operators which act on functions g defined on R + ×E.

T g(t, x, y) := ∂ t g(t, x) + f (y)∂ ν g(t, x), Ag(t, x, y) := T g(t, x, y) + Sg(t, x). Remark 4.2. For all functions g defined on R + ×E, T g(t, x, x) = T g(t, x), so that Ag(t, x, x) = Ag(t, x).

The next theorem provides Itô formulas for the PDMP (x t ) and its approximation (x t ). For all s ∈ [0, T ], we set η(s) := T n + kh if s ∈ [T n + kh, (T n + (k + 1)h) ∧ T n+1 [ for some n ≥ 0 and for some k ∈ {0, . . . , ⌊(T n+1 -T n )/h⌋}. Theorem 4.2. Let (x t , t ∈ [0, T ]) and (x t , t ∈ [0, T ]) be a PDMP and its approximation respectively constructed in section 2.3 with x 0 = x 0 = x for some x ∈ E. For all bounded functions g : R + × E → R continuously differentiable with bounded derivatives, we have

g(t, x t ) = g(0, x) + t 0 Ag(s, x s )ds + M g t , ( 35
)
where M g t :=

t 0 E (g(s, y)g(s, x s-))q(dsdy) is a true F p t -martingale, and

g(t, x t ) = g(0, x) + t 0 Ag(s, x s , x η(s) )ds + M g t , ( 36
)
where,

M g t := t 0 E (g(s, y) -g(s, x s-))q(dsdy) is a true F p t -martingale. Proof of Theorem 4.2.
The proof of ( 35) is given in [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF]. We prove (36) following the same arguments. Since q = pp ′ , we have

M g t = k≥1 1 T k ≤t g(T k , x T k ) -g(T k , x - T k ) - t 0 Sg(s, x s )ds.
Consider the above sum. As in [START_REF] Davis | Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models[END_REF], we write, on the event {N t = n}, that

k≥1 1 T k ≤t g(T k , x T k ) -g(T k , x - T k ) = g(t, x t ) -g(0, x) -g(t, x t ) -g(T n , x T n ) + n-1 k=0 g(T k+1 , x - T k+1 ) -g(T k , x T k ) .
For all k ≤ n-1, we decompose the increment g(T k+1 , x -

T k+1 )-g(T k , x T k ) as a sum of increments on the intervals [T k +ih, (T k +(i+1)h)∧T k+1 ] ⊂ [T k , T k+1 ].
Without loss of generality we are led to consider increments of the form g(t, θ, φ θ (t, ν))g(ih, θ, y i (x)) for some i ≥ 0, t ∈ [ih, (i + 1)h] and for all x = (θ, ν) ∈ E where we recall that φ is defined by [START_REF] Graham | Stochastic simulation and Monte Carlo[END_REF]. The function g is smooth enough to write

g(t, θ, φ θ (t, ν)) -g(ih, θ, y i (x)) = t ih (∂ t g + f θ (y i (x))∂ ν g) (s, θ, φ θ (s, ν))ds.
Then, the above arguments together with definition 4.2 yields

g(t, x t ) -g(T n , x T n ) + n-1 k=0 g(T k+1 , x - T k+1 ) -g(T k , x T k ) = t 0 T g(s, x s , x η(s) )ds.
The following theorem gives us a way to represent the solution of the integro-differential equation (32) as the conditional expected value of a functional of the terminal value of the PDMP (x t ). It plays a key role in the proof of Theorem 4.1.

Theorem 4.3 (PDMP's Feynman-Kac formula [START_REF] Davis | Markov Models and Optimization[END_REF]). Let F : E → R be a bounded function. Then the integro-differential equation (32) has a unique solution u : R + × E → R given by

u(t, x) = E[F (x T )|x t = x], (t, x) ∈ [0, T ] × E.

Proof of Theorem 4.1

We provide a proof in two steps. First, we give an appropriate representation of the weak error

E[F (x T )] -E[F (x T )]
. Then, we use this representation to identify the coefficient c 1 in (33).

Step

1: Representing E[F (x T )] -E[F (x T )].
Let u denote the solution of (32). From Theorem 4.3 we can write

E[F (x T )] -E[F (x T )] = E[u(T, x T )] -u(0, x).
Then, the application of the Itô formula (36) to u at time T yields

u(T, x T ) = u(0, x) + T 0 Au(s, x s , x η(s) )ds + M u T .
Since (M u t ) is a true martingale, we obtain

E[u(T, x T ) -u(0, x)] = E T 0 Au(s, x s , x η(s) )ds . For s ∈ [0, T ] we have Au(s, x s , x η(s) ) = ∂ t u(s, x s ) + f (x η(s) )∂ ν u(s, x s ) + Su(s, x s ) (see Definition 4.2).
From the regularity of λ, Q and u (see assumptions 4.1 and 4.2), the functions ∂ t u, ∂ ν u and Su are smooth enough to apply the Itô formula (36) between η(s) and s respectively. This yields

∂ t u(s, x s ) = ∂ t u(η(s), x η(s) ) + s η(s) A(∂ t u)(r, x r , x η(r) )dr + M ∂tu s -M ∂tu η(s) , ∂ ν u(s, x s ) = ∂ ν u(η(s), x η(s) ) + s η(s) A(∂ ν u)(r, x r , x η(r) )dr + M ∂ν u s -M ∂ν u η(s) , Su(s, x s ) = Su(η(s), x η(s) ) + s η(s) A(Su)(r, x r , x η(r) )ds + M Su s -M Su η(s) .
Moreover, since η(r) = η(s) for r ∈ [η(s), s], we have

f (x η(s) )∂ ν u(s, x s ) = f (x η(s) )∂ ν u(η(s), x η(s) ) + s η(s) f (x η(r) )A(∂ ν u)(r, x r , x η(r) )dr + f (x η(s) )(M ∂ν u s -M ∂ν u η(s) ), so that Au(s, x s , x η(s) ) = Au(η(s), x η(s) , x η(s) ) + s η(s) Υ(r, x r , x η(r) )dr + M ∂tu s -M ∂tu η(s) + f (x η(s) )(M ∂ν u s -M ∂ν u η(s) ) + M Su s -M Su η(s) ,
where,

Υ(t, x, y) := A(∂ t u) + f (y)A(∂ ν u) + A(Su) (t, x, y). ( 37
)
Since Au(t, x, x) = Au(t, x), the first term in the above equality is 0 by Theorem 4.3. By using Fubini's theorem and the fact that (M ∂tu t ) and (M Su t ) are true martingales, we obtain

E T 0 M ∂tu s -M ∂tu η(s) ds = E T 0 M Su s -M Su η(s) ds = 0. Moreover, since (M ∂ν u t
) is a F p t -martingale, we have

E T 0 f (x η(s) )(M ∂ν u s -M ∂ν u η(s) )ds = T 0 E f (x η(s) )E[M ∂ν u s -M ∂ν u η(s) |F p η(s) ] ds = 0.
Collecting the previous results, we obtain

E[F (x T )]-E[F (x T )] = E T 0 s η(s) Υ(r, x r , x η(r)
)drds . We can compute an explicit form of Υ in term of u, f , λ, Q and their derivatives. Indeed, Υ is given by (37), and we have

A(∂ t u)(t, x, y) = ∂ 2 tt u(t, x) + f (y)∂ 2 tν u(t, x) + S(∂ t u)(t, x), f A(∂ ν u) (t, x, y) = f (y) ∂ 2 tν u(t, x) + f (y)∂ 2 νν u(t, x) + S(∂ ν u)(t, x) , A(Su)(t, x, y) = ∂ t (Su)(t, x) + f (y)∂ ν (Su)(t, x) + S(Su)(t, x).
The application of the Taylor formula to the functions

∂ 2 tt u, ∂ 2 tν u, ∂ 2 νν u, S(∂ t u), S(∂ ν u), ∂ t (Su), ∂ ν (Su) and S(Su) at the order 0 around (η(r), x η(r) ) yields Υ(r, x r , x η(r) ) = Υ(η(r), x η(r) , x η(r) )+ O(h). Setting Ψ(t, x) = Υ(t, x, x) and recalling that for r ∈ [η(s), s], η(r) = η(s) and that |s -η(s)| ≤ h, we obtain E[F (x T )] -E[F (x T )] = E T 0 (s -η(s))Ψ(η(s), x η(s) )ds + O(h 2 ).
Consider the expectation in the right-hand side of the above equality. We decompose the integral into a (finite) sum of integrals on the intervals [T n +kh, (T n +(k+1)h)∧T n+1 ] where Ψ is constant. Without loss of generality, we are led to consider integrals of the form t kh (skh)Cds for some k ≥ 0, t ∈ [kh, (k + 1)h] and C a bounded constant. We have t kh (skh)Cds = t-kh 2 t kh Cds moreover adding and subtracting h in the numerator of (tkh)/2 yields

t kh (s -kh)Cds = h 2 t kh Cds + t -(k + 1)h 2 t kh
Cds.

Since C is bounded we deduce that

t kh (s -kh)Cds = h 2 t kh Cds + O(h 2
). Since Ψ is assumed bounded and E[N T ] < +∞, the above arguments yields the following representation

E[F (x T )] -E[F (x T )] = h 2 E T 0 Ψ(η(s), x η(s) )ds + O(h 2 ). ( 38
)
Step 2: From the representation (38) to the expansion at the order one. In this step, we show that

E T 0 Ψ(η(s), x η(s) )ds = E T 0 Ψ(s, x s )ds + O(h).
First, we introduce the random variables Γ and Γ defined by Γ := T 0 Ψ(η(s), x η(s) )ds and Γ := T 0 Ψ(η(s), x η(s) )ds and write

E[|Γ -Γ|] = E 1 min(T τ † ,T τ † )≤T |Γ -Γ| + E 1 min(T τ † ,T τ † )>T |Γ -Γ| ,
where τ † is defined in Definition 3.1. Since Ψ is bounded and P(min

(T τ † , T τ † ) ≤ T ) = O(h) (see the proof of Theorem 3.1), we have E |Γ -Γ|1 min(T τ † ,T τ † )≤T = O(h). Now, recall from (22) that, on the event {min(T τ † , T τ † ) > T }, we have T k = T k and θ k = θ k for all k ≥ 1 such that T k ∈ [0, T ]. Thus, for all n ≤ N T and for all s ∈ [T n , T n+1 [ we have x η(s) = (θ n , φ θn (η(s) -T n , ν n )) and x η(s) = (θ n , φ θn (η(s) -T n , ν n )). Consequently, on the event {min(T τ † , T τ † ) > T } we have |Γ -Γ| ≤ N T n=0 T n+1∧T T n |Ψ(η(s), θ n , φ θn (η(s) -T n , ν n )) -Ψ(η(s), θ n , φ θn (η(s) -T n , ν n ))|ds.
From the regularity assumptions 4.1 and 4.2, the function ν → Ψ(t, θ, ν) is uniformly Lipschitz in (t, θ) with constant L Ψ as sum and product of bounded Lipschitz functions. Thus, from this Lipschitz property and the application of Lemma 2.1, we get

|Ψ(η(s), θ n , φ θn (η(s) -T n , ν n )) -Ψ(η(s), θ n , φ θn (η(s) -T n , ν n ))| ≤ L Ψ Ce LT (n + 1)h.
From the above inequality, we find that

E 1 min(T τ † ,T τ † )>T |Γ -Γ| ≤ L Ψ Ce LT T hE[N T (N T +1)]. Since N T ≤ N * T and E[N * T (N * T + 1)] < +∞ we conclude that E 1 min(T τ † ,T τ † )>T |Γ -Γ| = O(h). We have shown that E T 0 Ψ(η(s), x η(s) )ds = E T 0 Ψ(η(s), x η(s) )ds + O(h).
Secondly, from the regularity assumptions 4.1 and 4.2, the function (t, ν) → Ψ(t, θ, ν) is uniformly Lipschitz in θ. Moreover, for all s ∈ [0, T ] there exits k ≥ 0 such that both s and η(s) belong to the same interval

[T k , T k+1 [ so that x s = (θ k , φ θ k (s -T k , ν k )) and x η(s) = (θ k , φ θ k (η(s) -T k , ν k )).

Numerical experiment

In this section, we use the theoretical results above to apply the MLMC method to the PDMP 2-dimensional Morris-Lecar (shortened PDMP 2d-ML).

The PDMP 2-dimensional Morris-Lecar

The deterministic Morris-Lecar model has been introduced in 1981 by Catherine Morris and Harold Lecar in [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF] to explain the dynamics of the barnacle muscle fiber. This model belongs to the family of conductance-based models (just as the Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]) and takes the following form

dv dt = 1 C I -g Leak (v -V Leak ) -g Ca M ∞ (v)(v -V Ca ) -g K n(v -V K ) , dn dt = (1 -n)α K (v) -nβ K (v), (39)
where

M ∞ (v) = (1 + tanh[(v -V 1 )/V 2 ])/2, α K (v) = λ K (v)N ∞ (v), β K (v) = λ K (v)(1 -N ∞ (v)), N ∞ (v) = (1 + tanh[(v -V 3 )/V 4 ])/2, λ K (v) = λ K cosh((v -V 3 )/2V 4 ).
In this section we consider the PDMP version of (39) that we denote by (x t , t ∈ [0, T ]), T > 0, whose characteristics (f, λ, Q) are given by

• f (θ, ν) = 1 C I -g Leak (ν -V Leak ) -g Ca M ∞ (ν)(ν -V Ca ) -g K θ N K (ν -V K ) , • λ(θ, ν) = (N K -θ)α K (ν) + θβ K (ν), • Q (θ, ν), {θ + 1} = (N K -θ)α K (ν) λ(θ,ν) , Q (θ, ν), {θ -1} = θβ K (ν) λ(θ,ν) .
The state space of the model is E = {0, . . . , N K } × R where N K ≥ 1 stands for the number of potassium gates. The values of the parameters used in the simulations are

V 1 = -1.2 , V 2 = 18, V 3 = 2, V 4 =

Classical and Multilevel Monte Carlo estimators

In this section we introduce the classical and multilevel Monte Carlo estimators in order to estimate the quantity E [F (x T )] where (x t , t ∈ [0, T ]) is the PDMP 2d-ML and F (θ, ν) = ν for (θ, ν) ∈ E so that F (x T ) gives the value of the membrane potential at time T . Note that other possible choices are F (θ, ν) = ν n or F (θ, ν) = θ n for some n ≥ 2. In those cases, the quantity E [F (x T )] gives the moments of the membrane potential or the number of open gates at time T so that we can compute statistics on these biological variables. Let X := F (x T ). In the sequel it will be convenient to emphasize the dependence of the Euler scheme (x t ) on a time step h. We introduce a family of random variables (X h , h > 0) defined by X h := F (x T ) where for a given h > 0 the corresponding PDP (x t ) is constructed as in section 2.3 with time step h. In particular, the processes (x t ) for h > 0 are correlated through the same randomness (U k ), (V k ) and (N * t ). We build a classical Monte Carlo estimator of E[X] based on the family (X h , h > 0) as follows

Y MC = 1 N N k=1 X k h , (40)
where (X k h , k ≥ 1) is an i.i.d sequence of random variables distributed like X h . The parameters h > 0 and N ∈ N have to be determined. We build a multilevel Monte Carlo estimator based on the family (X h , h > 0) as follows

Y MLMC = 1 N 1 N1 k=1 X k h * + L l=2 1 N l N l k=1 (X k h l -X k h l-1 ), (41)
where (X k h l , X k h l-1 ), k ≥ 1 for l = 2, . . . , L are independent sequences of independent copies of the couple (X h l , X h l-1 ) and independent of the i.i.d sequence (X k h * , k ≥ 1). The parameter h * is a free parameter that we fix in section 5.4. The parameters L ≥ 2, M ≥ 2, N ≥ 1 and q = (q 1 , . . . , q L ) ∈]0, 1[ L with L l=1 q l = 1 have to be determined, then we set N l := ⌈N q l ⌉, h l := h * M -(l-1) . We also set X := F (x T ) RT where RT is defined as in Proposition 2.2 with an intensity λ and a kernel Q that will be specified in section 5.4 and let ( Xh , h > 0) be such that Xh := F (x T ) RT for all h > 0 . By Proposition 2.2, we have E[X] = E[ X] and E[X h ] = E[ Xh ] for h > 0. Consequently, we build likewise a multilevel estimator Ỹ MLMC based on the family ( Xh , h > 0).

The complexity of the classical Monte Carlo estimator Y MC depends on the parameters (h, N ) and the one of the multilevel estimators Y MLMC and Ỹ MLMC depends on (L, q, N ). In order to compare those estimators we proceed as in [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF] (see also [START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF]), that is to say, for each estimator we determine the parameters which minimize the global complexity (or cost) subject to the constraint that the resulting L 2 -error must be lower than a prescribed ǫ > 0.

As in [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF], we call V 1 , c 1 , α, β and Var(X) the structural parameters associated to the family (X h , h > 0) and X. We know theoretically from Theorem 3.1 (strong estimate) and Theorem 4.1 (weak expansion) that (α, β) = (1, 1) whereas V 1 , c 1 and Var(X) are not explicit (we explain how we estimate them in section 5.3). Moreover, the structural parameters Ṽ1 , c1 , α, β and Var( X) associated to ( Xh , h > 0) and X are such that α = α, c1 = c 1 (see (34)), β = 2 (see Theorem 3.2) and Ṽ1 , Var( X) are not explicit. The classical and the multilevel estimators defined above are linear and of Monte Carlo type in the sense described in [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF]. The optimal parameters of those estimators are then expressed in term of the corresponding structural parameters as follows (see [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF] or [START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF]). For a user prescribed

L 1 + log(|c1| 1 α h * ) log(M) + log(A/ǫ) α log(M) , A = √ 1 + 2α q q 1 = µ * (1 + ρ(h * ) β 2 ) q j = µ * ρ(h * ) β 2 n -β 2 j-1 +n -β 2 j √ nj-1+nj , j = 2, . . . , L; µ * = 1/ 1≤j≤L q j N 1 + 1 2α Var(X) 1+ρ(h * ) β 2 L j=1 n -β 2 j-1 +n -β 2 j √ nj-1+nj 2 ǫ 2 L j=1 qj (nj-1+nj )
(ǫ) = (1 + 2α) -1 2α ǫ |c 1 | 1 α , N (ǫ) = 1 + 1 2α Var(X) 1 + ρh β/2 (ǫ) 2 ǫ 2 , ( 42
)
where ρ = V 1 /Var(X). The parameters of the estimator Y MLMC are given in Table 1 where n l := M l-1 for l = 1, . . . , L with the convention n 0 = n -1 0 = 0. The parameters of Ỹ MLMC are given in a similar way using Ṽ1 , β and Var( X). Finally, the parameter M (ǫ) is determined as in [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF] section 5.1.

Methodology

We compare the classical and the multilevel Monte Carlo estimators in term of precision, CPUtime and complexity. The precision of an estimator Y is defined by the

L 2 -error Y -E[X] 2 = (E[Y ] -E[X]
) 2 + Var(Y ) also known as the Root Mean Square Error (RMSE). The CPUtime represents the time needed to compute one realisation of an estimator. The complexity is defined as the number of time steps involved in the simulation of an estimator. Let Y denote the estimator (40) or (41). We estimate the bias of Y by

b R = 1 R R k=1 Y k -E[X],
where Y 1 , . . . , Y R are R independent replications of the estimator. We estimate the variance of

Y by v R = 1 R R k=1 v k ,
where v 1 , . . . , v R are R independent replications of v the empirical variance of Y . In the case where Y is the crude Monte Carlo estimator we set

v = 1 N (N -1) N k=1 (X k h -m N ) 2 , m N = 1 N N k=1 X k h . If Y is the MLMC estimator, we set v = 1 N 1 (N 1 -1) N1 k=1 (X k h -m (1) N1 ) 2 + L l=2 1 N l (N l -1) N l k=1 (X k h l -X k h l-1 -m (l) N l ) 2 ,
where m

(1)

N1 = 1 N1 N1 k=1 X k h and for l ≥ 2, m (l) N l = 1 N l N l k=1 X k h l -X k h l-1 . Then, we define the empirical RMSE ǫ R by ǫ R = b 2 R + v R . ( 43
)
The numerical computation of (43) for both estimators (40) and (41) requires the computation of the optimal parameters given by (42) and in table 1 of section 5.2 which are expressed in term of the structural parameters c 1 , V 1 and Var(X). Moreover the computation of the bias requires the value E[X]. Since there is no closed formula for the mean and variance of X we estimate them using a crude Monte Carlo estimator with h = 10 -5 and N = 10 6 . The constants c 1 and V 1 are not explicit, we use the same estimator of V 1 as in [START_REF] Lemaire | Multilevel Richardson-Romberg extrapolation[END_REF] section 5.1, that is

V 1 = (1 + M -β/2 ) -2 h -β E |X h -X h/M | 2 , ( 44
)
and we use the following estimator of c 1

c 1 = 1 -M -α -1 h -α E X h/M -X h . ( 45
)
The estimator of c 1 is obtained writing the weak error expansion for the two time steps h and h/M , summing and neglecting the O(h 2 ) term. In (44) we use (h, M ) = (0.1, 4) and in (45), we use (h, M ) = (1, 4) and the expectations are estimated using a classical Monte Carlo of size N = 10 4 on (X h/M , X h ). We emphasize that we interested in the order of c 1 and V 1 so that we do not need a precise estimation here.

Numerical results

In this section we first illustrate the results of Theorems 3.1 and 3.2 on the Morris-Lecar PDMP, then we compare the MC and MLMC estimators. The simulations were carried out on a computer with a processor Intel Core i5-4300U CPU @ 1.90GHz × 4. The code is written in C++ language. We implement the estimator Ỹ MLMC (see section 5.2) for the following choices of the parameters ( λ, Q).

Case 1: λ(θ) = 1 and Q θ, {θ + 1} = N K -θ N K , Q θ, {θ -1} = θ N K . Case 2: λ(x, t) = λ(θ, v(t)) and Q((x, t), dy) = Q((θ, v(t))
, dy) where v denotes the first component of the solution of (39). Cases 1 and 2 correspond to the application of Proposition 2.2. Based on Corollary 2.2 we also consider the following case.

Case 3: Consider the quantity E[F (x T ) -F (x T )] where (x t ) and (x t ) are PDPs with characteristics (Φ, λ, Q) and ( Φ, λ, Q) respectively. By Corollary 2.2, we have

E[F (x T )] = E[F (y T ) RT ]
where (y t ) is a PDP whose discrete component jumps in the same states and at the times as the discrete component of (x t ) do and ( Rt ) is the corresponding corrective process. Thus, we consider the quantity E

[F (x T ) -F (y T ) RT ] instead of E[F (x T ) -F (x T )].
The case 3 implies to use the following MLMC estimator which is slightly different from (41).

Ỹ MLMC = 1 N 1 N1 k=1 X k h * + L l=2 1 N l N l k=1 X k h l -Xk h l-1 ,
where (X k h l , Xk h l-1 ), k ≥ 1 for l = 2, . . . , L are independent sequences of independent copies of the couple (X h l , Xh l-1 ) = (F (x T ), F (y T ) RT ) where (y t ) is a PDP whose discrete component jumps in the same states and at the same times as the Euler scheme (x t ) with time step h l do, whose deterministic motions are given by the approximate flows with time step h l-1 and ( Rt ) is the corresponding corrective process (see Corollary 2.2). The figure 2

confirms numerically that E[|X h l -X h l-1 | 2 ] = O(h l ) and that E[| Xh l -Xh l-1 | 2 ] = O(h 2
l ) for the cases 1,2 and 3 (see Theorems 3.1 and 3.2 respectively). Indeed, for T = 10 (see figure 2a), we observe that the curve corresponding to the decay of E[|X h l -X h l-1 | 2 ] as l increases is approximately parallel to a line of slope -1 and that the curves corresponding to the decay of E[| Xh l -Xh l-1 | 2 ] in the cases 1,2 and 3 are parallel to a line of slope -2. We also see that the curves corresponding to the cases 2 and 3 are approximately similar and that for some value of l those curves go below the one corresponding to

E[|X h l -X h l-1 | 2 ]
. The curve corresponding to the case 1 is always above all the other ones, this indicates that the L 2 -error (or the variance) in the case 1 is too big (w.r.t the others) and that is why we do do not consider this case in the sequel. As T increases (see figures 2b and 2c), the theoretical order of the numerical schemes is still observed. However, for T = 20, a slight difference begin to emerge between the cases 2 and 3 (the case 3 being better) and this difference is accentuated for T = 30 so that we do not represent the case 2. For the Monte Carlo simulations we set T = 30, λ * = 10 and the time step involved in the first level of the MLMC is set to h * = 0.1. We choose this value for h * because it represents (on average) the size of an interval [T * n , T * n+1 ] of two successive jump times of the auxiliary Poisson process (N * t ). The estimation of the true value and variance leads E[X] = -31.4723 and Var(X) = 335. Note that v(30) = -35.3083 where v is the deterministic membrane potential solution of (39) so that there is an offset between the deterministic potential and the mean of the stochastic potential. We replicate 100 times the simulation of the classical and multilevel estimators to compute the empirical RMSE so that R = 100 in (43). As an example, the first line of table 3 reads as follows: for a user prescribed ǫ = 2 -1 = 0.5, the MLMC estimator Y MLMC is implemented with L = 2 levels, the time step at the first level is h * = 0.1, this time step is refined by a factor n l = M l-1 with M = 2 at each levels and the sample size is N = 2600. For such parameters, the numerical complexity of the estimator is Cost(Y MLMC ) = 28200, the empirical RMSE ǫ 100 = 0.389 and the computational time of one realisation of Y MLMC is 0.362 seconds. We also reported the empirical bias b 100 and the empirical variance v 100 in view of (43). However, the difference between the complexity of these two MLMC estimators increases more slowly than the one between a MC and a MLMC estimator. Recall that the computational benefit of the MLMC over the MC grows as the prescribed ǫ decreases. Both classical and multilevel estimators provide an empirical RMSE which is close to the prescribed precision (see tables 2, 3 and 4). We can conclude that the choice of the parameters is well adapted. For the readability, figures 3a, 3b show the ratios of the complexities and the CPU-times of the three estimators Y MC , Y MLMC and Ỹ MLMC as a function of ǫ.

2 E T 0 Ψ

 20 Thus, from the Lipschitz continuity of Ψ, from the fact that |sη(s)| ≤ h and sincef θ is uniformly bounded in θ we have |Ψ(s, x s ) -Ψ(η(s), x η(s) )| ≤ Ch where C is a constant independent of h. Then, we obtain sup s∈[0,T ] |E[Ψ(s, x s )] -E[Ψ(η(s), x η(s) )]| ≤ Ch from which we deduce that E T 0 Ψ(η(s), x η(s) )ds -E T 0 Ψ(s, x s )ds ≤ CT h.Finally, the weak error expansion reads E[F (x T )] -E[F (x T )] = h (s, x s )ds + O(h 2 ).

  30, λ K = 0.04, C = 20, g Leak = 2, V Leak = -60, g Ca = 4.4, V Ca = 120, g K = 8, V K = -84, I = 60, N K = 100. Membrane potential as a function of time. Red curves: stochastic potential. Black curve: deterministic potential. Proportion of opened gates as a function of time.Red curves: stochastic gates (θ/N K ). Black curve: deterministic gates (n).

  Probability of opening a gate (Q(xt, {θt + 1})) as a function of time. Jump rate (λ(xt)) as a function of t.

Figure 1 :

 1 Figure 1: 10 trajectories of the characteristics of the PDMP 2d-ML on [0, 100].

Table 3 :

 3 50e-01 2.59e-01 1.69e-01 3.87e-02 1.55e+00 8.47e+03 3.15e-02 2.69e+05 3 1.25e-01 1.17e-01 6.25e-02 9.78e-03 8.80e+003.34e+04 1.58e-02 2.12e+06 4 6.25e-02 5.67e-02 2.73e-02 2.47e-03 5.62e+011.32e+05 7.88e-03 1.68e+07 5 3.12e-02 2.50e-02 -1.78e-03 6.21e-04 3.93e+02 5.24e+05 3.94e-03 1.33e+08 Table 2: Results and parameters of the Monte Carlo estimator Y MC . Estimated values of the structural parameters: c 1 = 4.58, V 1 = 7.25. Results and parameters of the Multilevel Monte Carlo estimator Y MLMC . Estimated values of the structural parameters: c 1 = 4.58, V 1 = 7.25. Xh l -Xh l-1 ) 2 ] E[( Xhl -Xhl-1 ) 2 ]: Case 1 E[( Xhl -Xhl-1 ) 2 ]: Case 2 E[( Xhl -Xhl-1 ) 2 ]: Case 3 (a) T=10.

  Xh l -Xh l-1 ) 2 ] E[( Xhl -Xhl-1) 2 ]: Case 2 E[( Xhl -Xhl-1)2 ]: Case 3 (b) T=20.

  Xh l -Xh l-1 ) 2 ] E[( Xhl -Xhl-1 ) 2 ]: Case 3 (c) T=30.

Figure 2 :

 2 Figure 2: The plots (a),(b) and (c) show the decay of E[(X h l -X h l-1 ) 2 ] and E[ Xh l -Xh l-1 ) 2 ] (y-axis, log M scale) as a function of l with h l = h × M -(l-1) , h = 1, M = 4, for different values of the final time T . For visual guide, we added black solid lines with slopes -1 and -2.

Table 4 :

 4 The results indicate that the MLMC outperforms the classical MC. More precisely, for small values of ǫ (i.e k = 1, 2, 3) the complexity and the CPU-time of the classical and the multilevel MC estimators are of the same order. As ǫ decreases (i.e as k increases) the difference in complexity and CPU-time between classical and multilevel MC increases. Indeed, for k = 5 the complexity of the estimator Y MC is approximately 13 times superior to the one of Y MLMC and 19 times superior to the one of Ỹ MLMC . The same fact appears when we look at the complexity ratio of the estimators Y MLMC and Ỹ MLMC (i.e Cost(Y MLMC )/Cost( Ỹ MLMC )) as ǫ decreases. Results and parameters of the Multilevel Monte Carlo estimator Ỹ MLMC (case 3). Estimated values of the structural parameters: c1 = 3.91, Ṽ1 = 34.1.

  Y MLMC Ỹ MLMC : case 3 (a) Ratio of the complexities.

  Y MLMC Ỹ MLMC : case 3 (b) Ratio of the CPU-times.

Figure 3 :

 3 Figure 3: The plots (a) and (b) show the complexity and CPU-time ratios w.r.t the complexity and CPU-time of the estimator Ỹ MLMC as a function of the prescribed ǫ (log 2 scale for the x-axis, log scale for the y-axis).

  of all the other random variables of thinning that are present in g(x t )1 {Nt=n,τi=pi, 1≤i≤n, N * t =m} . The conditional expectation of g(x t )1 {Nt=n,τi=pi, 1≤i≤n,N *

	w.r.t. the vector (T	t =m}

Table 1 :

 1 Optimal parameters for the MLMC estimator (41). ǫ > 0, the classical Monte Carlo parameters h and N are h


Jean-Claude Bermond
email: jean-claude.bermond@inria.fr

Takako Kodate
email: kodate@lab.twcu.ac.jp

Joseph Yu
email: joseph.yu@ufv.ca

Gossiping with interference in radio chain networks

In this paper, we study the problem of gossiping with interference constraint in radio chain networks. Gossiping (or total exchange information) is a protocol where each node in the network has a message and wants to distribute its own message to every other node in the network. The gossiping problem consists in finding the minimum running time (makespan) of a gossiping protocol and efficient algorithms that attain this makespan.

Transmission model

The radio chain network is modeled as a symmetric dipath P n , where the vertices represent the nodes and the arcs represent the possible communications. A call (s, r) is defined as the transmission from the node s to the node r, in which s is the sender and r is the receiver and (s, r) is an arc of the dipath. The network is assumed to be synchronous and the time is slotted into steps. We suppose that each device is equipped with a half duplex interface; so, a node cannot both receive and transmit during a step.

Interference model Furthermore, communication is subject to interference constraints. We use a binary asymmetric model of interference based on the distance in the communication digraph like the ones used in [START_REF] Bermond | Optimal gathering protocols on paths under interference constraints[END_REF][START_REF] Bermond | Hardness and approximation of gathering in static radio networks[END_REF][START_REF] Klasing | On the complexity of bandwidth allocation in radio networks[END_REF]. Let d(s, r) denote the distance, that is the length of a shortest directed path, from s to r in P n and d I be a non negative integer. We assume that when a node s transmits, all nodes v such that d(s, v) ≤ d I are subject to the interference from s transmission. So two calls (s, r) and (s ′ , r ′ ) do not interfere if d(s, r ′ ) > d I and d(s ′ , r) > d I . During a given step only non interfering (or compatible) calls can be done and we will define a round as a set of such compatible calls. We focus here on the case where d I = 1.

Main result

This problem has been studied in general in [START_REF] Gasieniec | Gossiping with unit messages in known radio networks[END_REF] where approximation results are given (see also the survey [START_REF] Gasieniec | On efficient gossiping in radio networks[END_REF]). In [START_REF] Bermond | Gossiping with interferences in radio ring networks presented at JCDCG 3 2017[END_REF] we solved completely the gossiping problem in radio ring networks within this model. Here we determine exactly the minimum number of rounds R needed to achieve a gossiping when transmission network is a dipath P n on n nodes and the interference distance is d I = 1. We first prove the lower bound and then give gossiping algorithms which meet this lower bound.

Theorem 1

 1 The minimum number of rounds R needed to achieve a gossiping in a chain network P n (n ≥ 3), with the interference model d I = 1 is :R = { 3n -5 n ≥ 4 5 n = 3


Laurent Monasse

Virginie Daru

Christian Mariotti

Serge Piperno

A Conservative Coupling Method for Fluid-Structure Interaction in the Compressible Case

We propose a general explicit coupling method between a Finite Volume method for compressible fl w and a rigid body. The coupling strategy is based on the idea of Embedded Boundary methods (Pember et al., J. Comput. Phys. 120:278-304, 1995). The flu es are computed everywhere in the Cartesian grid, and are modifie at the solid boundaries to enforce flui mass conservation. The coupling between the flui and the solid is designed to ensure a balance in momentum and energy. We prove the exact numerical conservation of several simple uniform fl ws. An illustrative example of the liftoff of a cylinder by a shock wave is presented and compared with existing results.

Introduction

For fluid-structur interaction problems, two main types of methods were developed. The Arbitrary Lagrangian-Eulerian method (ALE) was firs developed in the late seventies [START_REF] Donea | An arbitrary Lagragian Eulerian finit element method for transient dynamic fluid-structur interactions[END_REF]. However, as a body-fitte method, the ALE technique requires remeshing of the flui domain when the solid is subjected to large deformations or breaking, which can be computationally demanding. In order to avoid bodyfittin and remeshing, Peskin [START_REF] Peskin | Numerical analysis of blood fl w in the heart[END_REF] proposed the Immersed Boundary method for the coupling of incompressible biological flui fl ws with moving elastic boundaries. Penalization type method have been applied to compressible fluid-structur interaction [START_REF] Aquelet | Euler-lagrange coupling with damping effects: Application to slamming problems[END_REF], but the stiff condition on the time-step as the penalization parameter is increased is unsatisfactory. Direct forcing methods, which modify the values of the flui cells in the vicinity of the solid boundary, have been widely applied [START_REF] Fadlun | Combined immmersed-boundary finite-di ference methods for three-dimensional complex fl w simulations[END_REF][START_REF] Mohd-Yusof | Combined immersed-boundary/b-spline methods for simulation of fl w in complex geometries[END_REF][START_REF] De Palma | An immersedboundary method for compressible viscous fl ws[END_REF], but still do not ensure the conservation of physical quantities at the interface.

L. Monasse (B)

Université Paris-Est, CERMICS, 77455 Marne-la-Vallée, France e-mail: monassel@cermics.enpc.fr Conservative methods for Finite Volumes with complex geometries have been developed, such as the Embedded Boundary method [START_REF] Pember | An adaptive Cartesian grid method for unsteady compressible fl w in irregular regions[END_REF]. Several procedures were proposed to avoid the stability condition restriction due to cut cells with small volume. We adopt here the solution consisting in merging the small cut-cells with their stable neighbors, that were successfully applied for compressible fluid-structur interaction [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF].

In the sequel the equations for flu modificatio remain similar to [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF]. The main difference lies in the coupling approach. Ultimately, our aim is to couple the compressible w with a Discrete Element method for the solid. As the Discrete Element method is computationally expensive, we use an explicit coupling method. In this paper, we restrict ourselves to rigid solids, but the method also applies for deformable bodies. The algorithm is carried out in a partitioned way: the flui pressure makes the solid move, and the updated position of the solid induces a modificatio of the flui flu es at the new boundary. The flui solver is based on the time splitting high order scheme developed in [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible fl w calculations[END_REF]. However the theoretical analysis of the coupling algorithm does not depend on the scheme that is used.

Description of the Method

The position of the solid in the flui domain is taken into account using the Embedded Boundary method [START_REF] Dragojlovic | An embedded boundary method for viscous, conducting compressible fl w[END_REF][START_REF] Falcovitz | A two-dimensional conservation laws scheme for compressible fl ws with moving boundaries[END_REF][START_REF] Hu | A conservative interface method for compressible fl ws[END_REF][START_REF] Miller | A conservative three-dimensional Eulerian method for coupled solidflui shock capturing[END_REF][START_REF] Pember | An adaptive Cartesian grid method for unsteady compressible fl w in irregular regions[END_REF]. At time t, the solid occupies a volume fraction α i of cell i, and all variables are assumed to be uniform in the cell. The conservative quantities contained in the cell are therefore equal to their value at the center of the cell times the volume of the cell and the volume fraction of flui 1 -α i . In the same way, the computed flu es are assumed to be constant on cell faces. If we denote κ i j the surface fraction of the face between cells i and j occupied by solid, we set the effective flu between i and j as the computed flu times the surface of their interface and the surface fraction of flui 1 -κ i j . Additional flu es come from the presence of the solid boundary, and are computed in order to yield exact conservation of flui mass and of total momentum and energy of the system.

We now describe the algorithm we developed, which preserves the fully conservative properties of the Finite Volume method combined with the Embedded Boundary method. At the beginning of a time step, at time n t, the known quantities are:

• The position of the center of mass of the solid particle X n and the rotation matrix Q n • The velocity of the center of mass of the solid particle V n and the angular quantity of movement matrix P n • The flui state: density n , velocity u n and pressure p n . The general algorithm is as follows:

SOLID FLUID COUPLING X n , Q n , V n , P n n , u n , p n (1) Computation of fluxes F num ˜, ũ, p (2)
Pressure is transferred to the solid boundary (3) Solid step (using boundary pressure)

X n+1 , Q n+1 , V n+1 , P n+1
(4) Update of the boundary position, computation of the α n+1 and κ n+1

(5) Fluid update :

⎧ ⎨ ⎩ n+1 = n + Δ n+1 u n+1 = n u n + Δ( u) n+1 e n+1 = n e n + Δ( ) n+1 , u n+1 , p n+1
The pressures p ref x and p ref y transferred to the solid boundary are chosen to be the temporal average of the reference pressures used during the partial steps of the flu computation (1) in the x and y directions (in our case, due to the splitting method; these steps could also be the steps of a Runge-Kutta method). The solid is advanced in step (3) using a classical second order Verlet scheme for translation and a second order RATTLE scheme for rotation. Details can be found in [START_REF] Monasse | An energy-preserving discrete element method for elastodynamics[END_REF]. The solid rigid body is supposed to be of general polygonal shape, and the computation of intersections in step (4) is straightforward.

Let the indices l, r , b and t denote the respective left, right, bottom and top faces of a cell C . Following [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF], the cut-cells computation in step ( 5) is written as follows in the 2D case:

• For each solid boundary F , compute the quantity w n F swept by boundary F , such that :

F w n F = C α n+1 C -α n C w n C ( 1
)
where w denotes alternately density , momentum u or energy e. In practice, we compute w n F as the integral of w n in the quadrangle define by the positions of F at time n t and (n + 1) t. -

If w = , F num F = 0 -If w = u, F num F = 1 S F F p ref x n x F 3 -If w = v, F num F = 1 S F F p ref y n y F -If w = e, F num F = 1 S F V n+ 1 2 F • F p ref x n x F p ref y n y F • Compute w C for each cell C : 1 -α n+1 C w C = t 1 -κ n+1 C l x F num C l - 1 -κ n+1 C r x F num C r + 1 -κ n+1 C b y F num C b - 1 -κ n+1 C t y F num C t + F ∈C t S F x y F num F + F ∈C w n F (2)
• Update the value of w C in every cell :

w n+1 C = w n C + w C (3)
To ensure stability of the method for small cut-cells, we follow the procedure described in [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF], with minor changes due to the impossibility to defin a normal vector in a cell occupied by two boundaries. We defin small cells as α C > 0.5. If we mix cell C with a target cell C t , so that the fina value w of the two cells is equal, we have to exchange the quantities

M C C t = α C t α C +α C t w C t -w C and M C t C = α C α C +α C t w C -w C t ,

and we easily check that w

C + M C C t = w C t + M C t C .
In the 2D case, we have to make a choice for the target cell C t . We f x C t to be the fully-flui cell (α C t = 0) nearest to cell C , such that the path between the two cells does not cross a solid boundary. A recursive subroutine find such a target cell after few iterations.

Let us note that the mixing procedure is entirely conservative, and ensures that the significan volume for a cell is consistent with the usual CFL condition based on standard cell size.

Theoretical Results

The following results were theoretically proven:

• Mass, momentum and energy conservation:

When there is no infl w into or outfl w from the domain, conservation of flui mass is ensured. For periodic boundary conditions, the momentum is exactly balanced between flui and solid during each time-step. For periodic or reflectin boundary conditions, the energy received by the flui from the solid is exactly balanced by the work of flui pressure forces on the solid during each time-step.

• Frame indifference:

Let an arbitrary shaped rigid body moving at constant velocity and no rotation, be immersed in a uniform flui fl wing at the same velocity. The uniform motion of both the solid and the flui is preserved by the coupling algorithm.

• Free slip along a straight boundary:

A uniform fl w parallel to a rigid semi-infinit half-plane is preserved by the coupling algorithm.

The last result shows that no numerical boundary layer or artificia boundary roughness appear at the solid boundary, even when it is not aligned with the Cartesian mesh.

Numerical Results

Due to space limitations, we only present here a moving boundary benchmark that was firs proposed in [START_REF] Falcovitz | A two-dimensional conservation laws scheme for compressible fl ws with moving boundaries[END_REF] and also treated in [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF]. A rigid cylinder of density 7.6 kg m -3 , initially resting on the lower wall of a two-dimensional channel fille with air at standard conditions, is driven and lifted upwards by a Mach 3 shock wave. The results obtained on a 1,600 × 320 grid are shown in Fig. 1. We observe good agreement with the results shown in [START_REF] Arienti | A level set approach to Eulerian-Lagrangian coupling[END_REF][START_REF] Hu | A conservative interface method for compressible fl ws[END_REF]. Small differences in the position of the shock waves can be noticed but no reference solution exists for this case. We also observe a strong vortex under the cylinder which is much weaker in [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF]. This vortex seems to be associated with a Kelvin-Helmholtz instability originating at the contact discontinuity present below the cylinder.

In addition, when considering the fina position of the center of mass of the cylinder, we observe a fast grid convergence of the method. The position we obtain on a 400 × 80-grid is comparable to that obtained on a 1,600 × 320-grid in [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF]. With increasing resolution on grids 400 × 80, 800 × 160 and 1, 600 × 320 [START_REF] Hu | A conservative interface method for compressible fl ws[END_REF] gives positions (0.659, 0.132), (0.649, 0.145) and (0.641, 0.147); on the same grids, we obtain (0.64375, 0.1463), (0.64278, 0.1471) and (0.64253, 0.1471).

•

  Compute the numerical flu es F num F at the solid boundary using p ref x and p ref y stored in the boundary F . The flu es are given by:

Fig. 1

 1 Fig. 1 60 contours of flui pressure from 0 to 28 at different times, x = y = 6.25 × 10 -4