Keywords: équations différentielles stochastiques, chemostat, processus de saut, approximation diffusion, méthode "tau-leap", méthode de Monte Carlo, algorithme de Gillespie

Modèles stochastiques du chemostat

Résumé : Nous reprenons la modélisation de la dynamique du chemostat à sa source. Le chemostat est classiquement représenté par un système d'équations différentielles. Notre objectif est d'établir un modèle stochastique qui est valable à l'échelle qui précède immédiatement celle qui correspond au modèle déterministe. Partant d'une échelle microscopique, nous présentons un modèle stochastique de sauts purs qui conduit, à l'échelle macroscopique, au modèle d'équation différentielle. À une échelle intermédiaire, une approximation diffusion nous permet de proposer un modèle sous la forme d'un système d'équations différentielles stochastiques. Nous détaillons les techniques qui permettent de passer d'une échelle à une autre ainsi que de simuler ces différents modèles. Nous décrivons également les domaines de validité des différents modèles.

Introduction

The evolution of the state of a single species/single substrate chemostat is usually described by a set of ordinary differential equations (ODE) derived from a mass balance principle, see [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]. More precisely, if s(t) denotes the concentration of nutrient (substrate) and b(t) the concentration of the organism (biomass) at time t(expressed in g/L), then the couple x(t) = (b(t), s(t)) * is the solution of the following ODE [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]:

ḃ(t) = [µ(s(t)) -D] b(t) , (1a)
ṡ(t) = -k µ(s(t)) b(t) + D [s in -s(t)] (1b)
where D > 0 is the dilution rate, s in > 0 the substrate concentration in the influent, and k > 0 the stoichiometric coefficient. The initial condition lies in the positive orthant, that is b(0) ≥ 0 and s(0) ≥ 0. Equation (1) will also be denoted:

ẋ(t) = f (x(t)) .
The specific growth rate function µ(s) is non-negative; we suppose that µ(0) = 0, µ(s) > 0 for s > 0, µ(s) ≤ µ max < ∞ and that it is continuous at 0. Commonly used models are the Monod model (uninhibited growth) and the Haldane model (inhibited growth) that reads respectively:

µ(s) = µ max s k s + s , µ(s) = µ max s k s + s + s 2 ki . (2)
This approach relies on the fact that the stochastic effects can be neglected, thanks to the law of large numbers, or at least can be averaged out. Although this level of description is sufficient for a number of applications of interest, it could be a valuable way of accounting for the stochastic nature of the system. Indeed, at small population sizes the chemostat could present stochastic behaviors, also the accumulation of small perturbations in the context of multi-species could not be neglected. Moreover, whereas the experimental results observed in well mastered laboratory conditions match closely the ODE theoretical behavior, a noticeable difference may occur in operational conditions. In these cases, stochastic features may not be neglected. We aim to build a model that still relies on a mass balance principle and that encompasses the useful stochastic information.

Many works [START_REF] Stephanopoulos | A stochastic analysis of the growth of competing microbial populations in a continuous biochemical reactor[END_REF][START_REF] Thomas | Asymptotic and transient analysis of stochastic core ecosystem models[END_REF][START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF] propose to superpose a stochastic term on Equation [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF] in order to model the uncertainty on the phenomenon, principally due to imprecise experimental conditions. Paradoxically, this amounts to the addition of an ad hoc perturbation to a model that has been obtained by neglecting these perturbations. We propose instead to consider the stochastic aspect at the very beginning of the modeling process, and to determine the conditions under which it is insignificant. This approach is not individual-based per se, as it starts from the macroscopic model [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF]. However, the first stochastic model proposed will be described at the individual level. This method will allow for a justification of the specific structure of the stochastic perturbation that affects the mean behavior. More generally, we will outline a modeling strategy based on many available tools, either stochastic or deterministic, depending on the regularity of the phenomenon to be modeled. In this paper we focus on the modeling and simulation process rather than on the mathematical developments; moreover we make use of known mathematical results. Our goal is to establish a stochastic model that is valid at the scale immediately preceding the one corresponding to the deterministic model [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF].

The paper is organized as follows: in Section 2, we recall the origin of model [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF] and the assumptions ensuring its validity. We show that since different timescales naturally appear in the problem, these assumptions need to be checked at each scale. Section 3 is devoted to the different models: the pure jump description that will be considered as the reference model is introduced in Section 3.1; the discrete time approximation, Poisson and normal, are presented in Section 3.2; the discrete-time normal approximation appeared to be a time discretization of a diffusion process given by a stochastic differential equation presented in Section 3.3. In Section 3.4 we describe the asymptotic results that bridge these different models. Section 4 is devoted to the associated simulation algorithm, Section 5 to numerical tests.

Scale and geometry issues in ODE model

An individual-based model should keep track of the position in space of each cell, together with their current biological states, it should also account for discrete events such as the division of a cell. Such a description of the system at the finest level could be of interest but unnecessary in view of our goals, namely to set a macroscopic model that gives account for stochastic phenomena. At this scale, the system is reduced to a R 2 -vector and its dynamics.

Model [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF] is obtained according to the classical approach, by choosing a small time interval ∆t on which a mass balance principle is applied to the state. However, ∆t should be large enough as we do not describe the dynamic at the timescale of jumps of one unit of substrate or bacteria but rather at the timescale of jumps of packet units. Such an interval could be called macroscopically infinitesimal [START_REF] Daniel | The chemical Langevin equation[END_REF].

Mass balance

Let (B t , S t ) denote the true concentrations at time t, assumed to be constant throughout the medium. The mass balance on interval [t, t + ∆t) reads

B t+∆t -B t = ∆B bio t + ∆B out t , (3a)
S t+∆t -S t = ∆S bio t + ∆S in t + ∆S out t (3b) where

• ∆B bio t and ∆B out t are the increments of biomass due to natural growth and to the outflow respectively, within [t, t + ∆t),

• ∆S bio t , ∆S in t and ∆S out t are the increments of substrate due to the consumption by the biomass, the inflow and the outflow respectively, within [t, t + ∆t).

Since we want to obtain an ODE, we now assume that the stochastic fluctuations are negligible relative to the increments. Again this requires ∆t to be large enough, so that sufficiently many discrete events have occurred. Moreover, ∆t should be taken even larger in case of inhomogeneity of the dynamics.

We denote by ( b(t), s(t)) for t = 0, ∆t, 2 ∆t, . . . the deterministic sequence constructed by using the mean increments of (B t , S t ):

b(t + ∆t) -b(t) = E[∆B bio t + ∆B out t ] , s(t + ∆t) -s(t) = E[∆S bio t + ∆S in t + ∆S out t ] .
Next, using the mass action law for the biomass we have

E[∆B bio t ] ≃ µ(s(t)) b(t) ∆t , E[∆S bio t ] ≃ -k µ(s(t)) b(t) ∆t (4)
where µ(s) is the specific growth rate and k > 0 the stoichiometric coefficient. Note that we again require ∆t to be large enough, since µ(s) and k make sense only for a sufficiently large population of bacteria. Now, since we have assumed perfect homogeneity of the medium, we get:

E[∆B out t ] ≃ -D b(t) ∆t , E[∆S in t ] ≃ -D s in ∆t , E[∆S out t ] ≃ -D s(t) ∆t . ( 5
)
Note that (4) and ( 5) are approximations because we have used a constant value for b(t) and s(t) within [t, t + ∆t). For this approximations to be correct, none of the quantities involved should vary significantly within [t, t + ∆t). We finally obtain the construction of the sequence

( b(t), s(t)) by b(t + ∆t) -b(t) = [µ(s(t)) -D] b(t) ∆t, (6a)
s(t + ∆t) -s(t) = (-k µ(s(t)) b(t) + D [s in -s(t)]) ∆t, (6b)
Model ( 1) is obtained by letting ∆t → 0 in System [START_REF] Stewart | Markov Processes -Characterization and Convergence[END_REF]. However, since ∆t is bounded from below, some care should be taken when this limit is achieved. System (6) can be understood as the discretization of (1) using an explicit Euler scheme with time-step ∆t. Whenever there exists ∆t sufficiently small, the deterministic sequence ( b(t), s(t)) will be close to model [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF], sampled at time 0, ∆t, 2 ∆t, . . .

Geometry and scales

The mass balance established in (3) features five terms that can be gathered according to the three sources of variations. This gives rise to a geometric structure that can be emphasized by writing (1) under the form:

d dt b(t) s(t) = µ(s(t)) b(t) 1 -k biology +D 0 s in inflow -D b(t) s(t) outflow (7)
However, whereas the geometry is well captured, the timescale of the five original terms is not readable in (1) nor [START_REF] Thomas | Asymptotic and transient analysis of stochastic core ecosystem models[END_REF]. Indeed, the fact that the approximations in ( 4) and ( 5) may be of different quality for each individual term is not exploited at all.

Models at different scales

In the previous section, we mentioned that the lower bound for ∆t is related to the size of the population and to the regularity of the phenomenon. Often, the experimental conditions are such that this bound is low enough, so that System (6) is correctly approximated by (1) sampled with period ∆t. If a smaller period is to be considered, then the conditions under which (6) has been obtained are not fulfilled. Particularly, the stochastic fluctuations should be accounted for.

We now introduce a stochastic process built on the same premise, that is a mean mass balance principle at a given ∆t. This model will have (1) as a fluid limit as ∆t goes to 0. This latter model suitably features the geometry of the chemostat but, as a limit model, cannot feature all its natural scales. The proposed stochastic models will respect both the geometry and the natural scales of the chemostat. We first establish a pure jump process representation of the chemostat at a microscopic scale, then we derive a diffusion process representation which will be valid at mesoscopic and macroscopic scales.

3.1 Pure jump model X t = (B t , S t ) * Even if do not aim at deriving an individual-based model, we try to preserve the discrete feature in the dynamics. We achieve this by considering only aggregated jumps obtained by adding up small and frequent jumps resulting from individual events. The resulting stochastic process will be a pure jump process X t = (B t , S t ) * , fully determined by its jumps and the corresponding jump rates; the state variable will be denoted x = (b, s) * .

In view of (3), we are led to consider five jumps: biology term: biomass increase of size ν 1 (x) at rate λ 1 (x); biology term: substrate decrease of size ν 2 (x) at rate λ 2 (x); inflow term: substrate inflow of size ν 3 (x) at rate λ 3 (x);

RR n°7458

outflow term: biomass outflow of size ν 4 (x) at rate λ 4 (x); outflow term: substrate outflow of size ν 5 (x) at rate λ 5 (x);

(see Figure 1). It remains to set the jump size rates so as to comply with the mass balance principle and the stochastic mass action law. For a macroscopically infinitesimal ∆t, denote by ∆X We first focus on the first two expressions. The stochastic mass action law [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF] requires

E[∆X b,bio t |X t = x] ≃ µ(s) b ∆t 0 , E[∆X s,bio t |X t = x] ≃ 0 -k µ(s) b ∆t .
Now notice that, for small ∆t, the number of jumps of type (resp. ) within [t, t + ∆t) is approximately P(λ 1 (x) ∆t) (resp. P(λ 2 (x) ∆t)), so that

E[∆X b,bio t |X t = x] ≃ λ 1 (x) ∆t ν 1 (x) , E[∆X s,bio t |X t = x] ≃ λ 2 (x) ∆t ν 2 (x) .
So we are looking for (λ i (x), ν i (x)) satisfying:

λ 1 (x) ν 1 (x) = µ(s) b 0 and λ 2 (x) ν 2 (x) = 0 -k µ(s) b . (8)
We therefore introduce the scale parameters K 1 and K 2 and we choose:

λ 1 (x) def = K 1 µ(s) b , ν 1 def = 1 K 1 0 , λ 2 (x) def = K 2 k µ(s) b ν 2 def = - 0 1 K 2 .
this choice is not unique and will be explain later in Section 3.5.

Here by "scale" we mean that jumps due to i will be of magnitude 1 Ki and the corresponding rates will be of magnitude K i . Large K i yields frequent and small jumps. Using the Poisson argument mentioned above, we see that these scale parameters K i do not act on the mean values of the increments but on their variances (large K i will correspond to small variances). The K i 's can thus be regarded as tuning parameters quantifying the uncertainty or regularity of the corresponding source of variation.

Reproducing this discussion with the three other types of jumps, and considering only admissible jumps (in the positive orthant), we obtain a pure jump Markov process with rate coefficients λ i (x) and associated jumps ν i (x) defined in Table 1

λi(x) K1 µ(s) b K2 k µ(s) b K3 D s in K4 D b K5 D s jump νi(x) 1 K 1 0 - 0 1∧K 2 s K 2 0 1 K 3 - 1∧K 4 b K 4 0 - 0 1∧K 5 s K 5
Table 1: Rates and jumps of the five basic mechanisms of the pure jump process. Note that the jumps ν i (x) essentially do not depend on x except for the negative jumps near the border {x = (b, s) * ∈ R 2 + ; b = 0 or s = 0}.

About scales parameters

Let m b and m s denote the representative masses of a single bacteria and of a single molecule of substrate. Typically m b ≫ m s (e.g. m b ≃ 10 6 m s ). Hence:

0 < K i ≤ 1 m b , for i = 1, 4 and 0 < K i ≤ 1 m s , for i = 2, 3, 5 .
In most cases K i ≪ K j for i = 1, 4 and j = 2, 3, 5, but it is possible to adjust the coefficients K's to the specific application considered. For example K 2 will be large in laboratory experimental conditions, but for a real implementation the substrate inflow concentration could have a large variance. Also for the outflow, in regular conditions K 4 and K 5 could be large, but in bad mixing conditions they could be smaller. Finally K 1 could be smaller than K 4 , as the biomass concentration increase presents more variance than the substrate decrease (which is more regular as it is related to the diffusion of substrate across cell membranes).

As proved later in Lemma A.1, the jumps ν i (x) are essentially constant and equal to:

ν 1 def = 1 K1 0 , ν 2 def = - 0 1 K2 , ν 3 def = 0 1 K3 , ν 4 def = - 1 K4 0 , ν 5 def = - 0 1 K5 . (9)
Representation of X t

The constructive description of the process X t that has been just presented would be used for simulation purposes, see Section 4.1. Nevertheless, it should be completed by a more comprehensible and synthetic representation. This will require some mathematical developments which we summarize now and that are detailed in Appendix A.

First we should notice that the jump process X t can be represented as the following (jump) SDE:

X t = X 0 + 5 i=1 (0,t]×[0,∞) ν i (X u -) 1 {v≤λi(X u -)} N i (du × dv) ( 10
)
where N i are independent random Poisson measures with intensity measure du × dv (Lebesgue measure). The process X t can be described as a Markov process with infinitesimal generator Aφ(x) = lim t→0 1 t [Eφ(X x t )-φ(x)] where X x t is the process X t starting from x. This operator completely characterizes the law of the process X t . In Appendix A we prove that this process is non explosive, i.e. it is defined for all t ≥ 0; that it admits moments as soon as X 0 does; and that it is solution of (10) see Proposition A.3.

RR n°7458

Still representation [START_REF] Daniel | Improved leap-size selection for accelerated stochastic simulation[END_REF] is opaque. We can establish that the process X t essentially admits the following representation:

dB t = µ(S t ) B t -D B t dt + d m1 t √ K1 + d m4 t √ K4 , (11a)
dS t = -k µ(S t ) B t + D (s in -S t ) dt + d m2 t √ K2 + d m3 t √ K3 + d m5 t √ K5 (11b)
where mi t are independent square integrable martingales with zero mean. The exact representation (41) differs from [START_REF] Grasman | Breakdown of a chemostat exposed to stochastic noise volume[END_REF] only through terms (1 ∧ K i s) and (1 ∧ K i b) that are equal to 1 except on a very limited neighborhood of the axes.

The martingales mi t are of mean 0 and they are explicitly known, see (40), as well as their quadratic variation, see (42). From Equation [START_REF] Grasman | Breakdown of a chemostat exposed to stochastic noise volume[END_REF], the deterministic part of the process X t , its drift coefficient, appears to be essentially the classical ODE (1); and the stochastic part of this dynamics, the martingale terms, are of order 1/ √ K i .

Discrete time approximations

Poisson approximation Xtn = ( Btn , Stn ) *

For any small ∆t > 0 given, let t n = n ∆t. We propose a discrete time Poisson approximation ( Xtn ) n≥0 of (X t ) t≥0 : on the interval [t n , t n+1 ) we froze the rate functions λ i (X t ) to λ i (X tn ) so that we get a Poisson distribution. The jumps ν i (X t ) are also frozen to ν i (X tn ). Let X0 = X 0 , the approximation is defined by:

Xtn+1 = Xtn + 5 i=1 ν i ( Xtn ) P i n (∆t λ i ( Xtn )) (12)
where

(P i n (ρ)) n∈N,i=1•••5 are independent Poisson variables with intensities ρ. We have: E[ Xtn+1 | Xtn = x] = x + 5 i=1 ν i (x) E[P i n (∆t λ i ( Xtn ))| Xtn = x] = x + ∆t 5 i=1 ν i (x) λ i (x) (13)
and let

f K (x) def = 5 i=1 ν i (x) λ i (x) , (14)
f K (x) is "essentially" the r.h.s. function f (x) of the O.D.E. (1), more precisely f k (x) = f (x)
except near the axes (see. Lemma A.1). In other words, the infinitesimal increments of the conditional mean follow the O.D.E. [START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF]. Also:

cov[ Xtn+1 | Xtn = x] = 5 i=1 cov[ν i (x) P i n (∆t λ i ( Xtn ))| Xtn = x] = Σ2 1 0 0 Σ2 2 (15a) with Σ2 1 = 1 K 2 1 cov[P 1 n (∆t λ 1 (x))] + 1 K 2 4 (1 ∧ K 4 b) 2 cov[P 4 n (∆t λ 4 (x))] = ∆t 1 K1 µ(s) b + 1 K4 (1 ∧ K 4 b) 2 D b , (15b)
Σ2 2 = 1 K 2 2 (1 ∧ K 2 s) 2 cov[P 2 n (∆t λ 2 (x))] + 1 K 2 3 cov[P 3 n (∆t λ 3 (x))] + 1 K 2 5 (1 ∧ K 5 s) 2 cov[P 5 n (∆t λ 5 (x))] = ∆t 1 K2 (1 ∧ K 2 s) 2 k µ(s) b + 1 K3 D s in + 1 K5 (1 ∧ K 5 s) 2 D s . ( 15c
)
Diffusion approximation ξtn = ( βtn , σtn ) *

In [START_REF] Hamza | Conditions for integrability of Markov chains[END_REF], the variable P i n (∆t λ i (x)) is Poisson distributed with parameter ∆t λ i (x). When this parameter is large (greater than 10 or 20) then this last distribution is very close to the normal distribution of mean ∆t λ i (x) and variance ∆t λ i (x). Hence, we get a (discrete time) normal approximation ( ξtn ) n≥0 of (X t ) t≥0 by letting ξ0 = X 0 and, conditionally on ξtn-1 = x:

ξtn+1 = x + 5 i=1 ν i (x) N i n
where N i n are 5 independent Gaussian random variables :

N i n ∼ N λ i (x) ∆t , λ i (x) ∆t
So conditionally on ξtn = x, ξtn+1 is normal with mean ( 13) and covariance matrix [START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF]. Let ξtn = ( βtn , σtn ) * , given βtn = b and σtn = s:

βtn+1 = b + µ(s) -(1 ∧ K 4 b) D b ∆t + ∆t µ(s) b K1 w 1 n + ∆t (1∧K4 b) 2 D b K4 w 4 n (16a) σtn+1 = s + -(1 ∧ K 2 s) k µ(s) b + D s in -(1 ∧ K 5 s) D s ∆t + ∆t (1∧K2 s) 2 k µ(s) b K2 w 2 n + ∆t D s in K3 w 3 n + ∆t (1∧K5 s) 2 D s K5 w 5 n (16b)
where w i n are i.i.d. N (0, 1) random variables.

Boundary conditions

In both approximations ( 12) and ( 16), no mechanism prevents the processes Xtn or ξtn from staying within the positive orthant R 2 + . An ad hoc solution is to set the concentration to 0 whenever it becomes negative, see Section 4. 16) is the Euler-Maruyama time discretization of the diffusion process ξ t = (β t , σ t ) * solution of the following SDE:

Diffusion model ξ

t = (β t , σ t ) * A stochastic differential equation model System (
dβ t = µ(σ t ) -(1 ∧ K 4 β t ) D β t dt + µ(σt) βt K1 dW 1 t + (1∧K4 βt) 2 D βt K4 dW 4 t dσ t = -(1 ∧ K 2 σ t ) k µ(σ t ) β t + D s in -(1 ∧ K 5 σ t ) D σ t dt + (1∧K2 σt) 2 k µ(σt) βt K2 dW 2 t + D s in K3 dW 3 t + (1∧K5 σt) 2 D σt K5 dW 5 t
where W i t are independent standard Wiener processes. Note that this result can be obtained directly from the process X t without the help of the discrete-time approximation. Indeed the infinitesimal generator of process X t given by ( 26) is a difference operator, and by Taylor development, it can be approximated by a second order differential operator corresponding to a diffusion process [START_REF] Stewart | Markov Processes -Characterization and Convergence[END_REF].

For small K i 's the last system is equivalent considering:

dβ t = µ(σ t ) β t -D β t dt + µ(σt) βt K1 dW 1 t + D βt K4 dW 4 t dσ t = -k µ(σ t ) β t + D s in -D σ t dt + k µ(σt) βt K2 dW 2 t + D s in K3 dW 3 t + D σt K5 dW 5 t RR n°7458
then we can group the Brownian motions in the following way:

dβ t = µ(σ t ) -D β t dt + µ(σt) βt K1 + D βt K4 dW b t (17a) dσ t = -k µ(σ t ) β t + D (s in -σ t ) dt + k µ(σt) βt K2 + D s in K3 + D σt K5 dW s t (17b)
where W b and W s are independent standard Wiener processes.

Behavior of the system of SDE's near the axes System ( 16) is the Euler-Maruyama time discretization of the SDE (17) (for large K i 's). Even if the diffusion approximation is only valid for large values of the biomass and the substrate, we can study the behavior of ( 17) near the axes. As for the discrete-time normal approximation, we should clarify the boundary conditions. As we well see, the component β t given by (17a) will remain positive, but the component σ t given by (17b) could become negative. We must first require that µ(s) = 0 for s < 0. Then, note that each equation of ( 17) is related to the following well-known CIR model for interest rates: Remark 3.1 (Cox-Ingersoll-Ross model) Consider the one-dimensional SDE:

dX t = (a + b X t ) dt + σ X t dW t , X 0 = x 0 ≥ 0 . ( 18
)
with a ≥ 0, b ∈ R, σ > 0.
According to [START_REF] Lamberton | Introduction to Stochastic Calculus Applied to Finance[END_REF]Prop. 6.2.4], for all x 0 ≥ 0, X is a continuous process taking values in R + , and let τ = inf{t ≥ 0, X t = 0}, then:

(i) If a ≥ σ 2 /2, then τ = ∞ P x -a.s.; (ii) if 0 ≤ a < σ 2 /2 and b ≤ 0 then τ < ∞ P x -a.s.;
(iii) if 0 ≤ a < σ 2 /2 and b > 0 then P x (τ < ∞) ∈ (0, 1).

In the first case, X never reaches 0. In the second case X a.s. reaches the state 0, in the third case it may reach 0. If a = 0 then the state 0 is absorbing.

It is clear that β = 0 is an absorbing state for (17a), and when β = 0, (17b) reduces to

dσ t = D (s in -σ t ) dt + D s in K3 + D σt K5 dW s t
and from Remark 3.1, the solution of this SDE will stay on the half-line [-K5 K3 s in , ∞) and:

(i) if s in ( 1 K3 + 1 K5 ) ≥ 1 2 K 2 5 then σ t never reaches -K5 K3 s in ; (ii) if s in ( 1 K3 + 1 K5 ) < 1 2 K 2 5
then σ t reaches -K5 K3 s in in finite time and is reflected.

Indeed, it is enough to apply Itô formula to σt = D s in K3 + D σt K5 and to use the Remark 3.1. Note that, as K 5 is large, condition (i) is more realistic than condition (ii).

To extend the definition of (17) for negative value of σ, let suppose that µ(σ) = 0 for σ ≤ 0. As we seen, β t will stay non-negative and β = 0 is an absorbing state. Also σ t ≥ -K5 K3 s in and for large K 5 this state will be repulsive. Note that for small values of σ t , as the K i are large, the diffusion term in (17b) will be small and the drift part will be dominated by D s in so that σ t will increase fast and its probability to be negative will be small.

The fact that the substrate concentration could be "negative" is due to the normal approximation. This approximation is valid for large values of concentration and the validity of the diffusion system ( 17) is questionable for small concentration. Nonetheless we can study its properties. 17b). This solution is rather artificial and will not be retained.

RR n°7458

So the solution of the system [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF]

remains in the domain D = [0, ∞) × [-K5 K3 s in , ∞).
The proof that this system admits a strong solution with pathwise uniqueness is presented in Appendix B.

Asymptotic analysis

The convergence of the pure jump model [START_REF] Daniel | Improved leap-size selection for accelerated stochastic simulation[END_REF] or of the diffusion approximation [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] to the deterministic model (1) as all the K i → ∞ can be rigorously established.

Let X K t be the pure jump model defined at the beginning of Section 3.1, or as the solution of the Equation (10) for a given

K def = (K 1 , K 2 , K 3 , K 4 , K 5 ). Let ξ K
t be the solution of the SDE [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF]. Let x(t) be the EDO model solution of Equation (1). Then X K t converges toward x(t) in the following way: for all T > 0 and all δ > 0,

P sup 0≤t≤T X K t -x(t) ≥ δ -→ 0 (19)
as

K i → ∞ for all i = 1 • • • 5.
This result is not surprising if we consider the representation (41) of (X t ) t≥0 ; it was obtained in a context of martingale convergence theorems in [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF][START_REF] Thomas | Limit theorems for sequences of jump Markov processes approximating ordinary differential processes[END_REF] or in a more general context of convergence of sequences of infinitesimal generators in [START_REF] Stewart | Markov Processes -Characterization and Convergence[END_REF]. We can also prove the same type of convergence for the process ξ K t . Indeed, in Equation ( 17) the scale coefficients appears as 1/ √ K i in the diffusion part of the SDE, and the convergence clearly holds as all the K i tends to infinity.

Other models

As already noticed, the choice of (λ i (x), ν i (x)) satisfying ( 8) is not unique. We choose not to make the jump sizes depend on the state value x (except for the boundary conditions), only the jump rates depend on x. Another possibility is to choose jump sizes that depend on the state value x. For example instead of the choice of Table 1, we can choose:

λ1(x) = K1 µ(s) , λ2(x) = K2 k µ(s) , λ3(x) = K3 D , λ4(x) = K4 D , λ5(x) = K5 D and ν1(x) = b K 1 0 , ν2(x) = - 0 b K 2 , ν3(x) = 0 s in K 3 , ν4(x) = - b K 4 0 , ν5(x) = - 0 s K 5 RR n°7458
(if we neglect the boundary condition). Then in place of [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] we have the following set of equations:

dβ t = µ(σ t ) β t -D β t dt + µ(σt) K1 β t dW 1 t + D K4 β t dW 4 t (20a) dσ t = -k µ(σ t ) β t + D s in -D σ t dt + k µ(σt) K2 β t dW 2 t + D K3 s in dW 3 t + D K5 σ t dW 5 t (20b)
where W i are independent standard Wiener processes.

A three components model

Instead of the five components to , we can consider a case with three independent sources of jump variation. This example strictly preserves the geometry (7) by considering three independent sources of jump variation:

′ biology term: biomass increase and substrate decrease at scale K ′ 1 ; ′ inflow term: substrate inflow at scale K ′ 2 ; ′ outflow term: biomass and substrate outflow at scale K ′ 3 see Figure 2.

Again the jump sizes and rates should be chosen so as to satisfy the mass balance principle and the stochastic mass action, with no canonical choice. An ad hoc choice is given in the following table:

′ ′ ′ biology inflow outflow rate λ ′ i (x) K ′ 1 µ(s) b K ′ 2 D K ′ 3 D jump ν ′ i (x)   1 K ′ 1 - k∧K ′ 1 s K ′ 1   0 s in K ′ 2   - 1∧K ′ 3 x K ′ 3 b - 1∧K ′ 3 x K ′ 3 s   Table 2:
Rates and jumps of an ad hoc choice for three mechanisms of the pure jump process. Note that the third jump ν ′ 3 (x) is radial.

These jumps are now essentially equal to

ν ′ 1 def = 1 K ′ 1 1 -k , ν ′ 2 def = 1 K 2 0 s in , ν ′ 3 def = 1 K ′ 3 -b -s .
This setting forces the jumps to be directed along the corresponding vector field, which is a strong constraint. In particular, the stoichiometry is strictly respected: the production of 1 unit of biomass requires exactly k units of substrate. Moreover the outflow jump is always radial, so that the increments of biomass and substrate are again strongly linked. Notice that for this particular choice of λ ′ 3 and ν ′ 3 , the jump rate is constant but the jump size is not. In other words, the jump carries information both in the direction and the intensity of the variation. This will affect the qualitative behavior of the process and of its diffusion approximation, regarding extinction for example.

As for our canonical model, we obtain a SDE for the diffusion approximation of the jump process :

dβ t = µ(σ t ) -(1 ∧ K ′ 3 β t ) D β t dt + µ(σt) βt K ′ 1 dW 1 t + D K ′ 3 (1 ∧ K ′ 3 ξ t ) β t dW 3 t dσ t = -(k ∧ K ′ 1 σ t ) µ(σ t ) β t + D s in -(1 ∧ K ′ 3 ξ t ) D σ t dt + µ(σt) βt K2 (k ∧ K ′ 1 σ t ) dW 1 t + D K ′ 2 s in dW 2 t + D K ′ 3 (1 ∧ K ′ 3 ξ t ) σ t dW 5 t .
Notice that since W 1 and W 3 affect both components of ξ t , the quadratic variation process ξ t will not be a diagonal matrix. In comparison with [START_REF] Thomas | Asymptotic and transient analysis of stochastic core ecosystem models[END_REF], we write the vector form of the SDE for large K i 's:

d βt σt = µ(σ t ) β t 1 -k biology +D 0 s in inflow -D βt σt outflow dt + µ(σt) βt K ′ 1 1 -k dW 1 t biology + D K ′ 2 0 s in dW 2 t inflow + D K ′ 3 βt σt dW 3 t outflow .
The diffusion term appears as the conjunction of three perturbations acting along the three vector fields determined by the sources of variation. Moreover, the intensity of the noise could be different for each type of perturbation. Considering this model could therefore be of interest, if the geometric interpretation of the noise is meaningful, see [START_REF] Joannides | On geometry and scale of a stochastic chemostat[END_REF].

Comparison with the Imhof-Walcher model [START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF] We finally mention that the diffusion model appearing in [START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF], is obtained from (20) by letting K 1 = K 2 = K 3 = 0 which leads to:

dβ t = µ(σ t ) -D β t dt + c b β t dW b t (21a) dσ t = -k µ(σ t ) β t + D (s in -σ t ) dt + c s σ t dW s t (21b)
The choice of these coefficients is justified in [START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF] by constructing an approximating Markov chain, and then taking the limit as the sampling rate goes to 0. This model will be compared to the diffusion approximation [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] model on a simulation test in Section 5.4.

Simulation algorithms

We presented several models for the chemostat system: the pure jump model (X t ) t≥0 could be considered as a detailed model at the microscopic scale. The Poisson approximation ( Xtn ) n∈N given by ( 12) and the normal approximation ( ξtn ) n∈N given ( 16) are constant time step approximation of the pure jump process. Finally the diffusion process (ξ t ) t≥0 solution of the SDE ( 17) is a continuous time approximation of the pure jump process. The now present the three associated simulation algorithms that will be valid at different scales.

Pure jump model

The pure jump model in continuous time described in Section 3.1 can be exactly simulated thanks to the Gillespie algorithm, also called stochastic simulation algorithm, described in Algorithm 1.

When the rate coefficients λ i (x) are large the time increment will be small and the Gillespie algorithm is impractical. As the scale coefficients K i are large, the λ i (x), i = 3, are large only when β and σ are small; λ 3 (x) will remain large as it does not depend on x.

t ← 0, x ← x 0 save (t, x) while t ≤ T max do compute λ i (x) % see Table 1 λ = 5 i=1 λ i (x) ∆t ∼ Exp(λ) % exponential distribution u ∼ U[0, 1] % uniform distribution t ← t + ∆t if u ≤ λ 1 (x)/λ then x ← x + ν1 % biomass reproduction else if u ≤ {λ 1 (x) + λ 2 (x)}/λ then x ← [x -ν2 ] + % consumption else if u ≤ {λ 1 (x) + λ 2 (x) + λ 3 (x)}/λ then x ← x + ν3 % substrate inflow else if u ≤ {λ 1 (x) + λ 2 (x) + λ 3 (x) + λ 4 (x)}/λ then x ← [x -ν4 ] + % biomass outflow else x ← [x -ν5 ] + % substrate outflow end if save (t, x) end while Here ν1 = 1/K1 0 , ν2 = 0 1/K2 , ν3 = 0 1/K3 , ν4 = 1/K4 0 , ν5 = 0 1/K5 and [x] + is the projection on the positive quadrant: [x] + = [( β σ )] + = β∨0 σ∨0
. Algorithm 1: Gillespie algorithm (or stochastic simulation algorithm). The simulation of the previous model could be cumbersome for very high rates of event. In this case it is desirable to use the fixed time step Poisson approximation method (12) also called tau-leap [START_REF] Daniel | Approximate accelerated stochastic simulation of chemically reacting systems[END_REF]. Recently many papers have addressed the numerical analysis of this approximation scheme [START_REF] Rathinam | Consistency and stability of tau-leaping schemes for chemical reaction systems[END_REF][START_REF] Li | Analysis of explicit tau-leaping schemes for simulating chemically reacting systems[END_REF][START_REF] Anderson | Error analysis of tauleap simulation methods[END_REF]. In this method the time step should be small enough so that it fulfills the following "leap condition": the state change in any leap should be small enough that no rate function λ i (x) will experience a macroscopically significant change in its value, that is:

Poisson approximation

t ← 0, x ← x 0 save (t, x) while t ≤ T max do compute λ i (x) % see Table 1 λ = 5 i=1 λ i (x) compute m i (x), v i (x) % see (23) ∆t ← min i=1•••5 ε λ/|m i (x)| , ε 2 λ 2 /v i (x) t ← t + ∆t P i ∼ Poisson(λ i (x) ∆t) for i = 1 • • • 5 x ← [x + ν1 P 1 -ν2 P 2 + ν3 P 3 -ν4 P 4 -ν5 P 5 ] + save (t, x) end while Here ν1 = 1/K1 0 , ν2 = 0 1/K2 , ν3 = 0 1/K3 , ν4 = 1/K4 0 , ν5 =
λ i x + i ′ ν i ′ (x) P i ′ n (∆t λ i ′ (x)) -λ i (x) ≤ ε λ(x) (22)
for i = 1 • • • 5, where 0 < ε ≪ 1 is an error control parameter.

For this method to be practicable [START_REF] Daniel | Improved leap-size selection for accelerated stochastic simulation[END_REF] proposed an automatic and simple way of determining the largest time step ∆t compatible with the leap condition. Define:

m i (x) def = 5 i ′ =1 λ i ′ (x) ∇λ i (x) • ν i ′ , (23a)
v i (x) def = 5 i ′ =1 λ i ′ (x) ∇λ i (x) • ν i ′ 2 (23b) for i, i ′ = 1 • • • 5, and let ∆t = min i=1•••5 ε λ(x) |m i (x)| , ε 2 λ 2 (x) |v i (x)| ( 24
)
where ε is an error control parameter (0 < ε ≪ 1), see Algorithm 2. Note that in the original context the jumps ν i (x) do not depends on x, but in our situation they do not essentially depend on x, the dependence on x was introduced to handle the jump near the axes in order to avoid negative concentration.

Diffusion (normal) approximation

t ← 0, (β, σ) ← (β 0 , σ 0 ) save (t, β, σ) while t ≤ T max do w b ∼ N (0, 1), w s ∼ N (0, 1)

β ′ ← β + (µ(σ) -D) β ∆t + µ(σ) β K1 + D β K4 √ ∆t w b σ ′ ← σ + (-k µ(σ) β + D (s in -σ)) ∆t + k µ(σ) β K2 + D s in K3 + D σ K5 √ ∆t w s β ← [β ′ ] + % 0 is an absorbing state for β σ ← |σ ′ -σ min | + σ min % reflection at σ min = -K5
K3 s in for σ t ← t + ∆t save (t, β, σ) end while Algorithm 3: Normal approximation.

The normal approximation ( 16) can be slightly modified in order to take into account the qualitative behavior of the SDE (17) near the axes. We propose the following scheme:

βtn+1 = βtn + µ(σ tn ) -(1 ∧ K 4 βtn ) D βtn ∆t + √ ∆t µ(σt n ) βtn K1 + (1∧K4 βtn ) 2 D βtn K4 w b n + , (25a)
σtn+1 = σtn + -(1 ∧ K 2 σtn ) k µ(σ tn ) βtn + D s in -(1 ∧ K 5 σtn ) D s ∆t + √ ∆t (1∧K2 σtn ) 2 k µ(s) βtn K2 + D s in K3 + (1∧K5 σtn ) 2 D s K5 w s n -σ min + σ min . ( 25b
)
Indeed as β = 0 is an absorbing state for the component β t of the SDE, instead of the standard Euler-Maruyama (16a), we can use (25a) where [•] + is the positive part and w b n are i.i.d. N (0, 1) random variables. Also, to take into account that the component σ is reflected in σ min = -K5 K3 s in we use the scheme (25b) where w s n are i.i.d. N (0, 1) random variables. This discretization scheme was proposed in [START_REF] Diop | Sur la discrétisation et le comportement à petit bruit d'EDS unidimensionnelles dont les coefficients sont à dérivées singulières[END_REF] in the context of the CIR diffusion process. In order to get a positive substrate concentration we can consider σ+ tn = σtn ∨ 0 or let σ min = 0 in (25b). The simulation procedure is presented in Algorithm 3.

Remark 4.1 (Scales and hybrid simulation) The three algorithms proposed here are valid at different scales. In the Gillespie algorithm all the detailed microscopic jumps of the dynamics are simulated.

The idea of the Poisson approximation is to consider a time step ∆ that should be small enough so that the different event rates barely evolve in the time interval [t, t + ∆t], but large enough for the approximation to be worthwhile. Starting in x at t, the time step ∆t is given by (24) but if it is less than a few multiples of 1/λ(x) then the Gillespie algorithm should be preferred. Now the Poisson variables P i of parameter λ i (x) ∆t could be approximated by normal variables N (λ i (x) ∆t, λ i (x) ∆t) as soon as λ i (x) ∆t ≥ 20.

The simulation method can automatically switch from one algorithm to another one according to the scale. We can also imagine that different components of the state vector are simulated with different algorithms.

Simulation study

We present simulation results of the discretized diffusion model [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with Monod and Haldane specific growth rates (2). The ODE ( 1) is integrated with a Runge-Kutta1 scheme but the Euler scheme, corresponding to [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with K i = ∞, gives very close results.

In addition to the deterministic case (case 0 with K i = ∞ for all i), we consider 3 basic cases (see Table 3):

"Standard" scales: K 2,3,5 = 100 × K 1,4 corresponds to the "standard" case where the substrate concentration dynamics is closer to the deterministic case than the biomass concentration dynamics.

"Unstirred inflow/outflows" scales: K 1,2 = 100 × K 3,4,5 corresponds to the case where inflow and outflows are unstirred.

"Fluid substrate" scales: K 2,3,4 = ∞, in this case the substrate equation (17b) is deterministic, i.e. the substrate dynamics is in fluid limit.

"Biological only" scales: K 3,4,5 = ∞, in this case we consider that the randomness is only due to biological aspects of the system.

A first comparison of trajectories

We consider the set of parameters of Table 4 (set 1) for the Monod case in the "standard" scales K 1,4 = 10 6 and K 2,3,5 = 10 8 . In Figure 3 we present a simulation of the pure jump process with the Gillespie method. As expected, most of the events corresponds to small jumps of the substrate concentration. Before addressing the question of reliability of these algorithms, see next subsection, we first focus on the qualitative nature of the trajectories proposed by the various methods.

In Figure 4, a simulation in a short time horizon of 0.1 (h) is proposed with Gillespie method (exact simulation) and the Poisson approximation method (tau-leaping) with a very small error control parameter ε = 10 -6 . The corresponding trajectories are very similar though 10 6 events are needed for the Gillespie method and only 3500 time steps are needed for the tau-leap method. and K 2,3,5 = 10 8 and 1237928 events (blue: a realization of the jump process (X t ) 0≤t≤0.1 , green: the ODE (x(t)) 0≤t≤0.1 ) -the substrate dynamics (RIGHT) presents many small size jumps, the biomass dynamics (LEFT) less jumps but with higher amplitude. The corresponding phase-portrait is plotted in Figure 4 (LEFT).

Figure 4: LEFT : Simulation of the pure jump process with the Gillespie Algorithm 1 for K 1,4 = 10 6 and K 2,3,5 = 10 8 and 1237928 events (blue: a realization of the jump process (X t ) 0≤t≤0.1 , green: the ODE (x(t)) 0≤t≤0.1 ) -the substrate dynamics presents many small size jumps, the biomass dynamics less jumps but with higher amplitude -the final time of simulation is 0.1 (h). RIGHT: same simulation with the Poisson approximation ( Xtn ) 0≤tn≤0.1 (tau-leap) with 3456 events for ε = 10 -6 (the time step is ≃ 1.89 × 10 -5 ).

Figure 5 present a simulation on a realistic length of time 100 (h). Here only the Poisson approximation and the normal approximation are reliable. For the Poisson approximation we use ε = 10 -3 (so that ∆t is between 0.016612 and 0.034232, for 3004 time steps) and for the normal approximation we use 2000 (corresponding to ∆t = 0.025). Again, the associated trajectories are very similar.

Law of the concentrations at a given time t

In Figure 6, we propose a Monte Carlo simulation to approximate the marginal densities of the biomass concentration B t and of the substrate concentration S t at a given time t. We consider the set of parameters of Table 4 (set 1) for the Monod case in the "standard" scales K 1,4 = 10 5 and K 2,3,5 = 10 7 . We compute (S

(j) t , B (j)
t ) for t = 3 (h) for j = 1 • • • 20000 independent Monte Carlo trials of the pure jump process (with the Gillespie method), with the Poisson approximation (tau-leap method) and with the normal approximation. For the tau-leap method we choose a constant time step. For "Poisson 1" and "Normal 1" we use a step of 0.05, for "Poisson 2" and "Normal 2" we use a step of 0.5. For each test, we compute the approximate density of S t and B t from the sample with a kernel method. Hence we compare 5 probability density functions for each component S t and B t . We also compute the empirical mean and standard deviation associated with the sample obtained from the pure jump process and we plot the associated normal density.

Initial conditions are B 0 = 0.026 and S 0 = 0.26, corresponding to the case of Figure 5, which is quite far from the equilibrium state.

The conclusions are:

• The two approximations (Poisson and normal) are very close to the exact simulation of the pure jump process; the approximation with a larger step 0.5 is slightly different.

• The computation times2 are:

for the exact simulation of the pure jump process: 5 h 45 min 32.6 s;

for the Poisson approximation: 33.2 s (with the time step 0.05) and 4.6 s (with the time step 0.5);

for the normal approximation: 0.7 s (with the time step 0.05) and 0.1 s (with the time step 0.5).

In the present situation, where the parameters K i are rather high, and for non-small concentration of the biomass and the substrate, the exact simulation of the pure jump process (Gillespie method) should be avoided.

• The resulting empirical densities are very close to normal densities and the solution of the ODE coincide with the mean of these normal densities. 4 (Monod growth function) with "standard" scales K 1,4 = 10 5 and K 2,3,5 = 10 7 . We compute the substrate and the biomass concentrations at t = 3 (h) and for j = 1 • • • 20000 independent Monte Carlo trials. The empirical densities are obtained with a kernel approximation procedure. We also compute the empirical mean and standard deviation from the sample of the pure jump process and plot with a thick grey line the corresponding normal densities: the match is very good. The vertical doted line the value of the substrate and biomass concentration at time t = 3 given by the ODE; again it matches the mean of all the empirical densities (except the ones corresponding to the time step 0.5).

A second test is proposed in the case of the Haldane growth function (Set 2 of Table 5): see Figure 7. The conclusions are the same as for the Monod case.

Figure 7: Empirical densities for the substrate and the biomass concentrations S t and B t at time t = 0.2 obtained with the exact simulation of the pure jump process (S t , B t ) with the Gillespie method (blue line), with the Poisson approximation ( St , Bt ) with constant time step (red solid line for a step 0.005, red dash line for a step 0.01), with the normal approximation (σ t , βt ) with constant time step (green solid line for a step 0.005, green dash line for a step 0.01). The corresponding CPU time are respectively: 3 h 54 min 16.9 s, 13.5 s, 6.7 s, 0.1 s, 0.1 s. The simulation parameters are the Set 2 of Table 5 (Haldane growth function) with "standard" scales K 1,4 = 10 5 and K 2,3,5 = 10 7 . We compute the substrate and the biomass concentrations at t = 0.2 (h) and for j = 1 • • • 20000 independent Monte Carlo trials. The empirical densities are obtained with a kernel approximation procedure. We also compute the empirical mean and standard deviation from the sample of the pure jump process and plot with a thick grey line the corresponding normal densities: the match is very good. The vertical doted line the value of the substrate and biomass concentration at time t = 0.2 given by the ODE; again it matches with the mean of all the empirical densities. The initial condition (S 0 , B 0 ) is the stable equilibrium solution of the ODE (1) that does not correspond to the washout: hence these empirical densities could be considered as good approximations of the limit distribution of the stochastic chemostat. Table 5: Sets of simulation parameters and the corresponding Haldane growth rate functions s → µ(s). The horizontal doted line is the maximum capacity µ max and the vertical doted line is the asymptotic substrate concentration of the ODE corresponding to the non-washout case.

About the scales parameters

As we have seen, for large populations, the diffusion approximation ξtn = ( βtn , σtn ) given by ( 25) is very close to the reference pure jump model X t = (B t , S t ). So we now propose simulations of the diffusion approximation in the case of a Monod specific growth rate according to the scales scenarios of Table 3 • Case 1 ("standard"): see Figures 8 and12.

• Case 2 ("unstirred inflowoutflows"): see Figures 9 and13.

• Case 3 ("fluid substrate"): see Figures 10 and14.

• Case 4 ("biological only"): see Figures 11 and15. We can conclude that, at this level of population and scale:

• The stochasticity is negligible only in the Case 4 ("biological only") and at the highest scale level (cases m, 3).

• The ODE solution x(t) matches the (empirical) mean of the stochastic process at these scales (as the stochastic process is solution of a nonlinear equation, there is no reason for the mean of the stochastic process to coincide with the solution of the deterministic equations).

Equivalent results have been obtained for the Haldane case.

Comparison with the Imhof-Walcher model [15]

We compare the processes ξ t = (β t , σ t ) * given by the diffusion approximation model [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] with the one given by the ad hoc model [START_REF] Norris | Markov chains[END_REF]. The parameter are: K = 1, µ max = 1, D = 0.5,

RR n°7458

Figure 8: Diffusion approximation, Case 1, Table 3 ("standard case") / Simulation of (25) with Monod specific growth rate (2) with the parameters of Table 4 -Time evolution of the biomass concentration (top left), time evolution of the substrate concentration (top right), phase portrait biomass/substrate concentrations (bottom) according to 4 cases: case 0, case 1.1, case 1.2, case 1.3 (see Table 3). Cases 0 (deterministic) and 1.3 are identical. 3 ("unstirred inflow and outflows") / Simulation of [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with Monod specific growth rate (2) with the parameters of Table 4 -Time evolution of the biomass concentration (top left), time evolution of the substrate concentration (top right), phase portrait biomass/substrate concentrations (bottom) according to 4 cases: case 0, case 2.1, case 2.2, case 2.3 (see Table 3). Cases 0 (deterministic) and 2.3 are identical.

RR n°7458

Figure 10: Case 3, Table 3 ("substrate fluid limit case") / Simulation of [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with Monod specific growth rate (2) with the parameters of Table 4 -Time evolution of the biomass concentration (top left), time evolution of the substrate concentration (top right), phase portrait biomass/substrate concentrations (bottom) according to 4 cases: case 0, case 3.1, case 3.2, case 3.3 (see Table 3). Cases 0 (deterministic) and 3.3 are identical. 3 ("biological case") / Simulation of [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with Monod specific growth rate (2) with the parameters of Table 4 -Time evolution of the biomass concentration (top left), time evolution of the substrate concentration (top right), phase portrait biomass/substrate concentrations (bottom) according to 4 cases: case 0, case 4.1, case 4.2, case 4.3 (see Table 3). Cases 0 (deterministic) and 4.3 are identical. [START_REF] Norris | Markov chains[END_REF]. The parameters are chosen so that the biomass concentration evolves from 0.5 to about 7.5, and the substrate concentration from 5 to about 0.5. Also the limit distribution is lesser than 1 in the substrate and greater than 1 in the biomass. Indeed one of the main difference between ( 17) and ( 21) is than for state values less than 1 (resp. more greater than 1) the noise variance for the first model is greater (resp. lesser) than the noise variance for the second model. This example illustrates clearly that the two models differ substantially.

Discussion

We started from a reference pure jump model X t , described by rates/jumps structure of Table 1 or as a solution of the stochastic differential equation [START_REF] Daniel | Improved leap-size selection for accelerated stochastic simulation[END_REF]. The martingale decomposition (41) clearly describes that the dynamics of X t is the combination of the classical deterministic dynamics of the chemostat (1) plus martingale terms with coefficients 1/ √ K i and with explicitly known quadratic variations, see (42). These quadratic variation terms allow us to assess the difference between the stochastic model and the deterministic one.

We presented the explicit Monte Carlo simulation procedure, called Gillespie method, for the process X t . In standard cases, that is for high population levels (i.e. K i large), this procedure is not feasible as it requires us to simulate too many events. In this case, we presented the Poisson approximation (25b) and the normal approximation [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF], both in discrete-time. These approximations are valid only for large populations, i.e. about the axes, it is necessary to return to the pure jump process representation. In the application discussed here, the Poisson approximation is of little interest: it is more time-consuming than the diffusion approximation and valid only on a very limited scale range between the pure jump model and the normal approximation model.

In contrast with previous stochastic chemostat models [START_REF] Stephanopoulos | A stochastic analysis of the growth of competing microbial populations in a continuous biochemical reactor[END_REF][START_REF] Thomas | Asymptotic and transient analysis of stochastic core ecosystem models[END_REF][START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF] where the stochasticity was introduced according to an ad hoc approach, in the present work we propose a family of models where the structure of the noise emerges from the very dynamics and where the scale parameters can be tuned according to the problem under interest. In particular it allows us to propose hybrid models where the cell population dynamics features stochasticity as the substrate is in fluid dynamics (ODE), corresponding to the Case 3 of Table 3. This kind of model has already been proposed in [START_REF] Grasman | Breakdown of a chemostat exposed to stochastic noise volume[END_REF] in a three trophic levels case where the stochasticity appears only in the top level trophic as a stochastic logistic model and with fluid limit dynamics for the two other levels; it also has been proposed in [START_REF] Crump | Some stochastic features of bacterial constant growth apparatus[END_REF] with a pure jump process for the biomass dynamics and a fluid limit for the substrate. This approach can also be related to coupled slow/fast reactions in stochastic chemical kinetics [START_REF] Haseltine | Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics[END_REF][START_REF] Ball | Asymptotic analysis of multiscale approximations to reaction networks[END_REF].

The approach proposed here can be applied to any model of population dynamics especially in cases of difference of scale between the different dynamics (e.g. cell/substrate). The dynamics of interacting populations cannot be modeled by a single model but rather by a family of models whose domain of validity depends on the scale at which the dynamics are considered. For example the normal approximation model represented as stochastic differential equations [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] or the ODE model (1) are valid in high population levels, hence using such models to infer extinction characteristics like extinction time and extinction probabilities is not valid. This was already noticed by [START_REF] Pollett | Diffusion approximations for ecological models[END_REF] and [START_REF] Wilcox | Do life history traits affect the accuracy of diffusion approximations for mean time to extinction?[END_REF].

In most standard population scales of the chemostat the ODE model is justified. Also, the ODE framework proposes analysis, control and optimization tools that are more accessible than the one of the SDE context. Though, as seen, the stochasticity cannot be neglected in many situations. This stochasticity could be of small intensity in the present single species/single substrate situation but could deeply perturb multiple species/multiple substrates situations. The SDE model could be simulated at a small extra computational cost and offers a more realistic prediction tool. Indeed, as it can account for the variability of the experiments, the simulation of the SDE offers the possibility to explore in depth the potentialities of the dynamical systems.

Figure 16: Comparison of the processes ξ t = (β t , σ t ) * given by the diffusion approximation model [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] and by the ad hoc model [START_REF] Norris | Markov chains[END_REF]. Evolution of the biomass concentration (top) and substrate concentration (center) during time; final joint density approximation of concentration (bottom). The two models differ substantially: compared to the diffusion approximation, the ad hoc model overestimates (resp underestimates) the noise variance for concentration greater (resp. lesser) than 1.

The SDE model is also more adapted for the confrontation to the data as it allows us to build a statistical model and the associated likelihood function. One of the next important steps, that we will investigate in coming work, will be to propose an adapted statistical procedure to estimate the scale parameters K i , and in a second step to estimate the parameters (D, s in ...). In the future we will also investigate the long-term behavior of these models as well as their optimal command.

Appendices

A Representation of the process X t

Infinitesimal generator

We consider the Markov process X t with infinitesimal generator:

Aφ(x) = 5 i=1 λ i (x) [φ(x + ν i (x)) -φ(x)] (26)
for all φ : R 2 + → R continuous with compact support [6, Th. 8-3.1]. The infinitesimal generator can also be understood in the following way:

P(X t+∆t = x + y|X t = x) ≃    λ i (x) ∆t + o(∆t) , if y = ν i (x) for i = 1 • • • 5, 1 - 5 i=1 λ i (x) ∆t + o(∆t) , if y = 0, o(∆t) , otherwise.
or as

Aφ(x) = lim t→0
Eφ(X x t )φ(x) t where X x t is the process X t starting from x. It can be rewritten as:

Aφ(x) = λ(x) R 2 + [φ(y) -φ(x)] ρ(x, dy) with λ(x) def = 5 i=1 λ i (x) , (27)
ρ(x, dy)

def = 5 i=1 λi (x) δ x+νi(x) (dy) with λi (x) def = λ i (x) 5 i ′ =1 λ i ′ (x) (28)
We define the jump times: τ 0 = 0 and

τ n def = inf t > τ n-1 ; X t = X τn-1
and the embedded jump chain

Y n def = X τn .
It is well know that (i) Y n is a Markov chain on R 2 + with transition probability ρ(x, dy); (ii) for all n ≥ 1, conditionally on Y 0 , . . . , Y n-1 , the holding times τ 1τ 0 , . . . , τ nτ n-1 are independent and exponentially distributed of intensity parameters λ(Y 0 ), . . . , λ(Y n-1 ), see [START_REF] Norris | Markov chains[END_REF]. These properties are at are the basis of the Gillespie simulation algorithm (see Algorithm 1).

Non-explosion and existence of moments

To study the non-explosion and existence of moments, we define the mean jump size function:

m K (x) def = R 2 + (y -x) ρ(x, dy) = 5 i=1 ν i (x) λi (x) . (29)
and, for p ≥ 1

|m K | p (x) def = R 2 + |y -x| p ρ(x, dy) = 5 i=1 |ν i (x)| p λi (x) . (30)
Note that from (14):

f K (x) = λ(x) m K (x) (31)
and if we replace ν i (x) by ν i in (31) we get:

5 i=1 ν i λ i (x) = µ(s) b -D b -k µ(s) b + D s in -D s = f (x)
where f (x) is the right-hand-side function of the ODE (1). f K (x) is the instantaneous mean of the process X t , more precisely we will show in Proposition A.3, that E(X t+∆t |X t = x) ≃ f K (x) ∆t. As proved in the next lemma f K (x) is essentially f (x), so locally in time the mean of the process X t behaves like x(t).

Lemma A.1 Consider ν i (x) defined in Table 1, ν i defined by (9) and f K (x) defined by [START_REF] Ikeda | Stochastic Differential Equations and Diffusion Processes[END_REF].

First |ν i (x)| ≤ |ν i |. Then let R K def = x = ( b s ) ∈ R 2 + ; b ≤ 1 K4 or s ≤ 1 K2 or s ≤ 1 K5 . For x ∈ R K , ν i (x) = ν i , m K (x) = m(x), f K (x) = f (x). For all x ∈ R 2 + : |ν i (x) -ν i | ≤ 1 min j=1•••5 K j , i = 1, . . . , 5 , (32)
|m K (x) -m(x)| ≤ 1 min j=1•••5 K j , ( 33
) |f K (x) -f (x)| ≤ (1 -K 2 s) + k µ(s) b + (1 -K 4 b) + D b + (1 -K 5 s) + D s (34) so that f K (x) → f (x) when K i → ∞ for all i = 1, . . . , 5
and this convergence is uniform on any compact set of (0, ∞) 2 .

Proof For

x ∈ R K it is clear that ν i (x) = ν i so f K (x) = f (x)
. For all x, we have ν 1 (x) ≡ ν 1 and ν 3 (x) ≡ ν 3 and:

ν 2 (x) -ν 2 = 1 K2 0 (1-K2 s) + , ν 4 (x) -ν 4 = 1

K4

(1-K4 b) + 0 , ν 5 (x)ν 5 = 1 K5 0 (1-K5 s) + so we get (32). The following assertions of the lemma are straightforward.

2

Non-explosion and existence of moments are given by the following result:

Theorem A.2 (Hamza and Klebaner [START_REF] Hamza | Conditions for integrability of Markov chains[END_REF]) Suppose that λ(x) > 0 for all x ∈ R 2 + and that there exist C > 0 and p ≥ 1 such that E(|X 0 | p ) < ∞ and

λ(x) |f K | p (x) ≤ C (1 + |x| p ) , ∀x ∈ R 2 + (35)
then the Markov process (X t ) t≥0 is non-explosive, that is τ n → ∞ a.s., and for all T > 0 there exists C > 0 s.t. E(|X t | p ) ≤ C for all t ≤ T . Indeed λ(x) = i λ i (x) ≥ K 3 D s in > 0 and

|f K | p (x) = 5 i=1 |ν i (x)| p λi (x) ≤ µ(s) b K p-1 1 + k µ(s) b K p-1 2 + D s in K p-1 3 + D b K p-1 4 + D s K p-1 5 ≤ 1 K p-1 1 + 1 K p-1 2 + 1 K p-1 3 + 1 K p-1 4 + 1 K p-1

5

(1 + |x|) so that (35) is fulfilled for all p ≥ 1. Hence the Markov process X t with infinitesimal generator A defined by ( 26) is non-explosive and X t admits moments of all order for all t ≥ 0.

RR n°7458

Representation for the process X t

We now give a representation for the process X t as a solution of a stochastic differential equation driven by random Poisson measures:

Proposition A.3
The process X t is defined for all t ≥ 0 and it is solution of the jump SDE [START_REF] Daniel | Improved leap-size selection for accelerated stochastic simulation[END_REF] where N i are independent random Poisson measures with intensity measure du × dv (Lebesgue measure).

Proof We first verify that the stochastic integral in [START_REF] Daniel | Improved leap-size selection for accelerated stochastic simulation[END_REF] is defined. According to [14, § II-3], this integral is defined if:

E t 0 ∞ 0 |ν i (X u )| 1 {v≤λi(Xu)} dv du < ∞ .
We have

E t 0 ∞ 0 |ν i (X u )| 1 {v≤λi(Xu)} dv du ≤ 1 K i E t 0 λ i (X u ) du ≤ C t 0 (1 + E|X u |) du
which is finite according to Theorem A.2. 2

Consider the centered random Poisson measure:

Ñ i (du × dv) def = N i (du × dv) -du × dv .
According to [14, § II-3]

M i t def = t 0 ∞ 0 ν i (X u -) 1 {v≤λi(X u -)} Ñ i (du × dv) , i = 1, . . . , 5 (36)
are five independent square-integrable martingales with finite moments of all orders. Let:

M t def = 5 i=1 M i t . (37)
From (10):

X t = X 0 + t 0 f K (X u ) du + M t . (38)
We want to study the behavior of the martingales M i t as the K i → ∞. First note that M i t = m i t 0 for i = 1, 4 and M i t = 0 m i t for i = 2, 3, 5 with

m 1 t def = 1 K 1 t 0 ∞ 0 1 {v≤λ1(X u -)} Ñ 1 (du × dv) , (39a)
m 2 t def = - 1 K 2 t 0 ∞ 0 (1 ∧ K 2 S u -) 1 {v≤λ2(X u -)} Ñ 2 (du × dv) , (39b)
m 3 t def = 1 K 3 t 0 ∞ 0 1 {v≤λ3(X u -)} Ñ 3 (du × dv) , (39c)
m 4 t def = - 1 K 4 t 0 ∞ 0 (1 ∧ K 4 B u -) 1 {v≤λ4(X u -)} Ñ 4 (du × dv) , (39d)
m 5 t def = - 1 K 5 t 0 ∞ 0 (1 ∧ K 5 S u -) 1 {v≤λ5(X u -)} Ñ 5 (du × dv) . ( 39e
) RR n°7458
As m i t is of the form m i t = t 0 ∞ 0 γ i (u -, v) Ñ i (du×dv) then the associated predictable quadratic variation is m i t = t 0 ∞ 0 [γ i (u, v)] 2 dv du so we can easily check that

m 1 t = 1 K 1 t 0 µ(S u ) B u du , m 2 t = 1 K 2 t 0 (1 ∧ K 2 S u ) 2 k µ(S u ) B u du , m 3 t = 1 K 3 D s in t , m 4 t = 1 K 4 t 0 (1 ∧ K 4 B u ) 2 D B u du , m 5 t = 1 K 5 t 0 (1 ∧ K 5 S u ) 2 D S u du . Define mi t def = K i m i t = K i [0,t]×[0,∞) γ i (u -, v) Ñ i (du × dv) . ( 40
)
So we obtained the following representation of the process X t that emphases the dependence on the K i :

dB t = µ(S t ) B t -(1 ∧ K 4 B t ) D B t dt + 1 √ K1 d m1 t + 1 √ K4 d m4 t (41a) dS t = -(1 ∧ K 2 S t ) k µ(S t ) B t + D s in -(1 ∧ K 5 S t ) D S t dt + 1 √ K2 d m2 t + 1 √ K3 d m3 t + 1 √ K5 d m5 t (41b)
where mi t are independent square integrable martingales with the following quadratic variations:

m1 t = t 0 µ(S u ) B u du , (42a)
m2 t = t 0 (1 ∧ K 2 S u ) 2 k µ(S u ) B u du , (42b)
m3 t = D s in t , (42c)
m4 t = t 0 (1 ∧ K 4 B u ) 2 D B u du , (42d)
m5 t = t 0 (1 ∧ K 5 S u ) 2 D S u du . (42e)
B Existence and uniqueness for a solution of the SDE [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] To prove that the system (17) admits a strong solution and pathwise uniqueness holds we use the results of [24, p. 134] or [START_REF] Durrett | Stochastic Calculus: A Practical Introduction[END_REF]. Let D N = [ 1 N , N ] × [-K5 K3 s in + 1 N , N ] and rewrite [START_REF] Thomas | Solutions of ordinary differential equations as limits of pure jump Markov processes[END_REF] as:

dξ t = f (ξ t ) dt + g(ξ t ) dW t ( 43
)
where ξ t = (β t , σ t ) * and W t = (W b t , W s t ) * . The coefficient g is not Lipschitz continuous, it is globally Lipschitz on D N for all N , but we can find Lipschitz continuous coefficients f N and g N such that:

f N (ξ) = f (ξ) for ξ ∈ D N ,
f N (ξ) = 0 for ξ ∈ D 2N , g N (ξ) = g(ξ) for ξ ∈ D N , g N (ξ) = 0 for ξ ∈ D 2N .

The the system dξ N t = f N (ξ N t ) dt + g N (ξ N t ) dW t admits a unique strong solution ξ N . Let:

T N = inf{t ≥ 0 ; ξ N t ∈ D N } .
Hence if N ≥ M then ξ N = ξ M for all t ≤ T M , so we can define a process ξ ∞ such that:

ξ ∞ t = ξ N t ∀t ≤ T N
and will be solution of (43) up to the explosion time T ∞ = lim N →∞ T N . Then by stability property and by the fact that the solution cannot cross the boundary of D we get T ∞ = ∞ a.s. which proves the strong existence and pathwise uniqueness defined for all t ≥ 0.

RR n°7458

Figure 1 :

 1 Figure 1: In this model, from a position x = (b, s) * the process could jump according to 5 mechanisms (2 due to the biology, 1 inflow, and 2 outflows), the basic jump i has a length 1 Ki

Figure 2 :

 2 Figure 2: In this simplified model, from a position x = (b, s) * the process could jump according to 3 mechanisms (biology, inflow, and outflow), the basic jump i ′ has a length 1 K ′ i for i = 1, . . . , 3

0 1 /

 1 K5 and [x] + is the projection on the positive quadrant: [x] + = [( β σ )] + = β∨0 σ∨0 . Algorithm 2: Poisson approximation or tau-leap method.

Figure 3 :

 3 Figure3: Simulation of the pure jump process with the Gillespie Algorithm 1 for K 1,4 = 10 6 and K 2,3,5 = 10 8 and 1237928 events (blue: a realization of the jump process (X t ) 0≤t≤0.1 , green: the ODE (x(t)) 0≤t≤0.1 ) -the substrate dynamics (RIGHT) presents many small size jumps, the biomass dynamics (LEFT) less jumps but with higher amplitude. The corresponding phase-portrait is plotted in Figure4(LEFT).

Figure 5 :

 5 Figure 5: LEFT : Simulation of the Poisson approximation with the τ -leap Algorithm 1 for K 1,4 = 10 6 and K 2,3,5 = 10 8 and 3004 events (blue: a realization of the jump process ( Xtn ) 0≤tn≤100 , green: the ODE (x(t)) 0≤t≤100 ), ε = 10 -3 , ∆t between 0.028 and 0.0344. RIGHT: same simulation with the normal approximation ( ξtn ) 0≤tn≤100 (diffusion approximation) with 2000 time steps with ∆t = 0.05.

Figure 6 :

 6 Figure6: Empirical densities for the substrate and the biomass concentrations at time t = 3 obtained with the exact simulation of the pure jump process (S t , B t ) with the Gillespie method (blue line), with the Poisson approximation ( St , Bt ) with constant time step (red solid line for a step 0.05, red dash line for a step 0.5), with the normal approximation (σ t , βt ) with constant time step (green solid line for a step 0.05, green dash line for a step 0.5). The corresponding CPU time are respectively: 5 h 45 min 32.6 s, 33.2 s, 4.6 s, 0.7 s, 0.1 s. The simulation parameters are the Set 1 of Table4(Monod growth function) with "standard" scales K 1,4 = 10 5 and K 2,3,5 = 10 7 . We compute the substrate and the biomass concentrations at t = 3 (h) and for j = 1 • • • 20000 independent Monte Carlo trials. The empirical densities are obtained with a kernel approximation procedure. We also compute the empirical mean and standard deviation from the sample of the pure jump process and plot with a thick grey line the corresponding normal densities: the match is very good. The vertical doted line the value of the substrate and biomass concentration at time t = 3 given by the ODE; again it matches the mean of all the empirical densities (except the ones corresponding to the time step 0.5).

Figures 8

 8 Figures 8 to 11 represent a simulation of a single trajectory in the 3 levels of scale: cases m, 1 to m, 3 (for m = 1 • • • 4). Figures 12 to 15 represent the result of 10000 Monte Carlo trials in the 3 levels of scale: cases m, 1 and m, 2 (for m = 1 • • • 4). We represent the mean trajectory and the empirical law of ξT at final time T .We can conclude that, at this level of population and scale:

Figure 9 :

 9 Figure 9: Diffusion approximation, Case 2, Table3("unstirred inflow and outflows") / Simulation of[START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with Monod specific growth rate (2) with the parameters of Table4-Time evolution of the biomass concentration (top left), time evolution of the substrate concentration (top right), phase portrait biomass/substrate concentrations (bottom) according to 4 cases: case 0, case 2.1, case 2.2, case 2.3 (see Table3). Cases 0 (deterministic) and 2.3 are identical.

Figure 11 :

 11 Figure 11: Case 4, Table3("biological case") / Simulation of[START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] with Monod specific growth rate (2) with the parameters of Table4-Time evolution of the biomass concentration (top left), time evolution of the substrate concentration (top right), phase portrait biomass/substrate concentrations (bottom) according to 4 cases: case 0, case 4.1, case 4.2, case 4.3 (see Table3). Cases 0 (deterministic) and 4.3 are identical.

Figure 12 :

 12 Figure 12: Cases 1.1 and 1.2 Table3/ Sampling 10000 Monte Carlo trials of the law of ( βtn , σtn ) for t n = 100 -The deterministic solution and the mean of the sampled trajectories coincide -The final law is represented by the sample and by the contour plot of the corresponding kernel approximation of the p.d.f.

Figure 13 :

 13 Figure 13: Cases 2.1 and 2.2 Table3/ Sampling 10000 Monte Carlo trials of the law of ( βtn , σtn ) for t n = 100 -The deterministic solution and the mean of the sampled trajectories coincide -The final law is represented by the sample and by the contour plot of the corresponding kernel approximation of the p.d.f.

Figure 14 :

 14 Figure 14: Cases 3.1 and 3.2 Table3/ Sampling 10000 Monte Carlo trials of the law of ( βtn , σtn ) for t n = 100 -The deterministic solution and the mean of the sampled trajectories coincide -The final law is represented by the sample and by the contour plot of the corresponding kernel approximation of the p.d.f.

Figure 15 :

 15 Figure 15: Cases 4.1 and 4.2 Table3/ Sampling 10000 Monte Carlo trials of the law of ( βtn , σtn ) for t n = 100 -The deterministic solution and the mean of the sampled trajectories coincide.

  .

	biomass	substrate	substrate	biomass	substrate
	increase	decrease	inflow	outflow	outflow
	biology	inflow		outflow
	rate

Table 3 :

 3 Simulation cases.

				Monod model
		set 1 set 2
	k	10	10	stoichiometric constant
	µmax 3	0.5	maximal growth rate (h -1 )
	D	0.12	0.4	dilution rate (h -1 )
	s in	0.5	10	input concentration (g/l)
	ks	6	1	half saturation constant (g/l)
			set 1	set 2

Table 4 :

 4 Sets of simulation parameters and the corresponding Monod growth rate functions s → µ(s). The horizontal doted line is the maximum capacity µ max and the vertical doted line the asymptotic substrate concentration of the ODE corresponding to the non-washout case.

RR n°7458

The routine ode45 of Matlab, an explicit Runge-Kutta (4,5) formula.RR n°7458

CPU time on a 2.13GHz Intel Core 2 Duo with a RAM of 2 GB.RR n°7458

Acknowledgements

The work was partially supported by the French National Research Agency (ANR) within the SYSCOMM project DISCO ANR-09-SYSC-003.