Bertrand Cloez
email: bertrand.cloez@supagro.inra.fr

Coralie Fritsch
email: coralie.fritsch@inria.fr

GAUSSIAN APPROXIMATIONS FOR CHEMOSTAT MODELS IN FINITE AND INFINITE DIMENSIONS

Keywords: chemostat model, central limit theorem on Hilbert-space, individual-based model, weak convergence, Crump-Young model, stationary and quasi-stationary distributions Mathematics Subject Classification (MSC2010): 60F05, 92D25, 60J25, 60G57, 60B10, 60H10

In a chemostat, bacteria live in a growth container of constant volume in which liquid is injected continuously. Recently, Campillo and Fritsch introduced a massstructured individual-based model to represent this dynamics and proved its convergence to a more classic partial differential equation.

In this work, we are interested in the convergence of the fluctuation process. We consider this process in some Sobolev spaces and use central limit theorems on Hilbert space to prove its convergence in law to an infinite-dimensional Gaussian process.

As a consequence, we obtain a two-dimensional Gaussian approximation of the Crump-Young model for which the long time behavior is relatively misunderstood. For this approximation, we derive the invariant distribution and the convergence to it. We also present numerical simulations illustrating our results.

Introduction

The chemostat is a biotechnological process of continuous culture developed by [START_REF] Monod | La technique de culture continue, théorie et applications[END_REF] and [START_REF] Novick | Description of the chemostat[END_REF] in which bacteria live in a growth container of constant volume in which liquid is continuously injected.

From a mathematical point of view, beyond classic models based on systems of ordinary differential equations (see for instance [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]) or integro-differential equations (see for instance [START_REF] Fredrickson | Statistics and dynamics of procaryotic cell populations[END_REF]; [START_REF] Ramkrishna | Statistical models of cell populations[END_REF][START_REF] Ramkrishna | Population Balances: Theory and Applications to Particulate Systems in Engineering[END_REF]), several stochastic models were introduced in the literature. The first-one seems to be the one developed by [START_REF] Crump | Some stochastic features of bacterial constant growth apparatus[END_REF] and is a birth and death process for the biomass growth coupled with a differential equation for the substrate evolution. This one is the main object of interest in Section 3 below. Recently, [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF] and Collet et al. (2013a) studied some extensions of this model. In particular, [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF] propose some stochastic differential equations to model the demographic noise from the microscopic interactions.

Other stochastic models were introduced by [START_REF] Stephanopoulos | A stochastic analysis of the growth of competing microbial populations in a continuous biochemical reactor[END_REF]; [START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF]. Let us also mention [START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF]; [START_REF] Mirrahimi | Evolution of species trait through resource competition[END_REF][START_REF] Mirrahimi | Direct competition results from strong competition for limited resource[END_REF] or, for individual-based models, Campillo et al. (2016,b); [START_REF] Champagnat | Adaptation in a stochastic multiresources chemostat model[END_REF]; [START_REF] Fritsch | A numerical approach to determine mutant invasion fitness and evolutionary singular strategies[END_REF] which model the evolutionary dynamics of the chemostat.

We focus here on the individual-based model developed by [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF] and [START_REF] Fritsch | A modeling approach of the chemostat[END_REF]. In this mass-structured model, the bacterial population is represented as a set of individuals growing in a perfectly mixed vessel of constant volume. This representation combines discrete mechanisms (birth and death events) as well as continuous mechanisms (mass growth and substrate dynamics). [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF] set the exact Monte Carlo simulation algorithm of this model and its mathematical representation as a stochastic process. They prove the convergence of this process to the solution of an integro-differential equation when the population size tends to infinity. In the present work, we investigate the study of the fluctuation process; namely the difference between the measure-valued stochastic process and its deterministic approximation. We first show that, conveniently normalized, this fluctuation process converges to some superprocess. Our proof is based on a classic tightness-uniqueness argument in infinite dimension. In contrast with [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF]; [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF]; [START_REF] Haskovec | Convergence of a stochastic particle approximation for measure solutions of the 2D Keller-Segel system[END_REF], one difficulty is that the main process is a signed measure and we have to find a suitable space in which it, as well as its limit, are to be immersed (because the space of signed measures endowed with the weak convergence is not metrizable). Inspired by [START_REF] Meleard | Convergence of the fluctuations for interacting diffusions with jumps associated with Boltzmann equations[END_REF] and [START_REF] Tran | Modèles particulaires stochastiques pour des problèmes d'évolution adaptative et pour l'approximation de solutions statistiques[END_REF], we consider the fluctuation process as an element of some Sobolev space (see Section 2.1 for a description of this space). This type of spaces takes the advantage to be Hilbertian and one can use martingale techniques on Hilbert spaces to obtain the tightness (and then the convergence of this process); see for instance [START_REF] Métivier | Convergence faible et principe d'invariance pour des martingales à valeurs dans des espaces de Sobolev[END_REF]. The limit object that we obtain is then an infinite dimensional degenerate Gaussian process.

We illustrate the interest of this result applying it in finite dimension. More precisely, for particular parameters, the mass-structured model of [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF] can be reduced to the two-dimensional Crump-Young model. As pointed out by Collet et al. (2013a), the long time behavior of this model is complex and misunderstood; only few properties are known about the behavior before extinction. The convergence developed by [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF] induces an approximation by an ordinary differential equation of the Crump-Young model, whereas our main result allows to obtain a stochastic differential approximation for which we are able to plainly describe the long-time behavior.

Our main results are described in section that follows: Theorems 1.2 and 1.3 are the central limit theorems (convergence of the fluctuation processes) in infinite and finite dimensions. Theorem 1.4 gives the long time behavior of a stochastic differential approximation of the Crump-Young model. Section 2 is devoted to the proofs of the two central limit theorems. We first, introduce all the notations and preliminaries that we need from Section 2.1 to Section 2.3, then Theorem 1.2 is proved in Section 2.4. The main steps of the proof of Theorem 1.3 are given in Section 2.5. The finite-dimensional case is studied in Section 3. We prove the convergence in time of the stochastic differential approximation of the Crump-Young model in Section 3.1. We present numerical simulations and discussion illustrating our results in Section 3.2. In particular, we discuss about the validity of the approximation and introduce another diffusion process, obtained from Theorem 1.3, whose numerical behavior seems to have a better mimic of the Crump-Young model in some particular situations. The extinction time of this new process is studied in Section 3.3.

Main results.

Let us be more precise on our main results before to introduce all the machinery (notations, Sobolev spaces, ...) that we will use.

We consider the following mass-structured chemostat model : each individual is characterized by its mass x P r0, M s, where M is the maximal mass of a bacterium. At each time t ě 0, the system is characterized by the random variable pS n t , ν n t q, where S n t is the substrate concentration and ν n t " ř N n t i"1 δ X i t is the population of the N n t individuals with mass X 1 t , ¨¨¨, X

N n t t
in the chemostat at time t. The parameter n represents a scaling parameter.

We assume that one individual with mass x P r0, M s ' divides at rate bpS, xq into two individuals with masses α x and p1 ´αq x where α is distributed according to a probability distribution Qpdαq on r0, 1s; ' is withdrawn from the chemostat at rate D, with D the dilution rate of the chemostat; ' grows at speed gpS, xq : 9

x t " gpS n t , x t q, where the substrate concentration S n evolves according to the following equation 9 S n t " D ps in ´Sn t q ´k n V

N n t ÿ i"1 gpS n t , X i t q , S n 0 " S 0 ,
where s in is the input substrate concentration in the chemostat, S 0 is a deterministic initial substrate concentration, n V is the volume of the chemostat and k is a stoichiometric coefficient. Note that the scale parameter n is only involved in front of the volume and the initial number of individuals. The approximations below then holds when the volume and the initial population become larger and larger. In this context, let us do a small remark on modelling. Parameter D corresponds to a dilution rate, which is usually defined as the ratio between the flow and the volume. As we assume that the dilution rate is constant, approximations below only hold when also the flow became larger and larger. A more complete description of the stochastic process is given in Section 2.1 in term of martingale problem. To have a better understanding of the dynamics let also see (Campillo and Fritsch, 2014, Section 2.2).

For every n ě 1, we consider the renormalized process ps ν n t q tě0 defined by

s ν n t " 1 n ν n t , t ě 0 (1)
and we make the following assumptions.

Assumptions 1.1 (Regularity of the division rate and the growth speed).

(1) (2) The function g P C 1,1 pR `ˆr0, M sq is such that gps, 0q " gps, M q " 0 .

(3) In absence of substrate the bacteria do not grow, i.e. gp0, xq " 0 for all x P r0, M s.

Note that due to the form of the differential equation satisfied by the substrate concentration pS n t q tě0 , one can see that it remains in the compact set r0, maxpS 0 , s in qs. As a consequence, the regularity of the functions g and b, induced by Assumptions 1.1, implies that the division rate and the growth speed are bounded : 0 ď bps, xq ď s b, 0 ď gps, xq ď s g, s ě 0, x P r0, M s.

With these assumptions, [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF] show that if the sequence ps ν n 0 q n converges in distribution towards a deterministic, finite and positive measure ξ 0 then, under Assumptions 1.1, the following limit holds in distribution (see Section 2.1 for details about the topology), lim nÑ8 pS n t , s ν n t q tPr0,T s " pS t , ξ t q tPr0,T s ,

for any horizon time T ą 0, where pS t , ξ t q tPr0,T s is the solution of the deterministic system of equations

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % S t " S 0 `şt 0 " D ps in ´Su q ´k V ş r0,M s gpS u , xq ξ u pdxq  du , ξ t pf q " ξ 0 pf q `şt 0 ş r0,M s " bpS u , xq ş 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ı Qpdαq ´D f pxq `gpS u , xq f 1 pxq  ξ u pdxq du , (3)
for any f P C 1 pr0, M sq, t ě 0, with, for any h P Cpr0, M sq and ν in the set M F pr0, M sq of finite (positive) measures,

ν phq :" ż r0,M s hpxq νpdxq . (4)
Let us finally introduce the main object of the present article, that is the fluctuation process s η n t " pη n t , R n t q defined by η n t :"

? n ps ν n t ´ξt q , R n t :"

? n pS n t ´St q .

(5)

Our main result is Theorem 1.2 below. For presentation convenience, we don't detail here the topology of Dpr0, T s, Hq, H or H 0 but all details are given in Section 2.1, in particular H and H 0 are defined in (16) (see also Remark 2.2). Briefly Dpr0, T s, Hq is the Skohorod space associated to an appropriately chosen Sobolev space H and H 0 Ĺ H. Theorem 1.2 (Convergence of the fluctuation process). Under Assumption 1.1 and if sup ně1 E ´}s η n 0 } 2 H 0 ¯ă 8 and ps η n 0 q ně1 converges to some s η 0 " pη 0 , 0q in H then, for any horizon time T ą 0, the sequence of process ps η n q ně1 converges in distribution in Dpr0, T s, Hq towards s η " pη, Rq solution of the system

$ ' ' ' ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' ' ' ' % η t pf q " η 0 pf q `şt 0 ş r0,M s " bpS u , xq ş 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ‰ Qpdαq ´D f pxq `gpS u , xq f 1 pxq  η u pdxq du `şt 0 R u ş r0,M s " B s bpS u , xq ş 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ‰ Qpdαq `Bs gpS u , xq f 1 pxq  ξ u pdxq du `Gt pf q R t " ´şt 0 " D R u `k V η u pgpS u , ¨qq `Ru k V ξ u pB s gpS u , ¨qq  du (6)
where Gpf q is a centred Gaussian process with quadratic variation

xGpf qy t " ż t 0 ż r0,M s bpS u , xq ż 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ‰ 2 Qpdαq ξ u pdxq du `D ż t 0 ż r0,M s f 2 pxq ξ u pdxq du, (7)
for any f P C 1 pr0, M sq and t P r0, T s.

Recall that the notations

η u pgpS u , ¨qq " ż r0,M s gpS u , xq η u pdxq , ξ u pB s gpS u , ¨qq " ż r0,M s B s gpS u , xq ξ u pdxq
have been defined in (4).

This theorem may look complicated but let us illustrate the interest of this type of result with a finite dimensional application. Let us choose M " 8, gps, xq " µpsqm, bps, xq " µpsq,

where s Þ Ñ µpsq is the specific growth rate of the population that we will assume to be Lipschitz. Even though the previous assumptions are not included in the set of assumptions of Theorem 1.2 (see however remark 2.8), we can obtain the same result for specific functions f , in particular, when f " 1 (see Theorem 1.3 below and its proof in Section 2.5).

For parameters (8), the substrate concentration of the stochastic model satisfies

d dt S n t " D ps in ´Sn t q ´k V n m µpS n t q N n t
where the process pN n t q t , depicting the number of individuals, is a birth-death process with non-homogeneous birth rate µpS n t q and death rate D. It is exactly the Crump-Young model as studied in [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF]Collet et al. (2013a); [START_REF] Crump | Some stochastic features of bacterial constant growth apparatus[END_REF]. In particular, the long time behavior of this process is investigated in Collet et al. (2013a). It is shown that this process extincts after a random time and, under suitable assumptions (µ increasing,...), admits (at least) a quasi-stationary distribution. This distribution describes the behavior of the process before the extinction (when it is unique and there is convergence to it); see for instance [START_REF] Collet | Quasi-stationary distributions[END_REF].

Let us define

s N n t " s ν n t p1q " N n t n , N t " ξ t p1q, Q n t " ? np s N n t ´Nt q.
We have then the following result.

Theorem 1.3 (Convergence of the Crump-Young fluctuation process).

If sup ně1 Ep| s N n 0 | 2 |Q n 0 | 2 q ă `
8 and the sequence of random variables p s N n 0 , Q n 0 q ně1 converges in distribution towards pN 0 , Q 0 q then the sequence of processes pp s N n t , S n t , Q n t , R n t q tě0 q ně1 converges in distribution in Dpr0, T s, R 4 q towards pN, S, Q, Rq solution of the following system of stochastic differential equations:

$ ' ' ' & ' ' ' % dN t " pµpS t q ´Dq N t dt , dS t " " D ps in ´St q ´k V m µpS t q N t ‰ dt , dQ t " rpµpS t q ´DqQ t `µ1 pS t qR t N t s dt `apµpS t q `Dq N t dB t , dR t " ´"D R t `k V µpS t q m Q t `k V R t µ 1 pS t q m N t ‰ dt . ( 9
)
for all t ě 0, where pB t q tě0 is a classic Brownian motion.

The previous theorem suggests, if n is sufficiently large, that

N n t « p N n t :" n N t `?n Q t , S n t « p S n t :" S t `1 ? n R t , with p p N n t , p S n t q tě0 solution of # d p N n t " " pµpS t q ´Dq p N n t `µ1 pS t qp p S n t ´St q n N t ‰ dt `apµpS t q `Dq n N t dB t d p S n t " " D ps in ´p S n t q ´k V n m µpS t q p N n t ´k V n µ 1 pS t qp p S n t ´St q n N t ‰ dt. ( 10
)
Note that another (SDE type) approximation is given in Section 3.2. This is a Fellerdiffusion type approximation (see [START_REF] Bansaye | ume 1 of Mathematical Biosciences Institute Lecture Series[END_REF]) and it is closer to the SDE introduced in [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF].

The two first equations of (9) are, up to a factor m{V in front of N t , the classic differential equations for representing the chemostat (see [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]). The four-component process is a non-elliptic diffusion time-homogeneous process whose long time behavior is given by Theorem 1.4 below.

Theorem 1.4 (Long time behavior of the Crump-Young SDE). Assume that µ is strictly increasing on r0, s in s, µp0q " 0 and µps in q ą D. There exists a unique pN ˚, S ˚q such that

µpS ˚q " D , N ˚" V km ps in ´S˚q ,
and for any initial condition in R ˚ˆR `ˆR ˆR, the process ppN t , S t , Q t , R t q T q tě0 (T designs the transpose of the vector) converges in distribution to a Gaussian random variable with mean m and variance Σ defined by m "

¨N S0 0 ‹ ‹ ' , Σ " ¨0 0 0 0 0 0 0 0 0 0 α ´µpS ˚q µ 1 pS ˚q 0 0 ´µpS ˚q µ 1 pS ˚q β ‹ ‹ ‹ ' , (11)
where

α " `k V m µ 1 pS ˚q N ˚`3 2 D ˘2 ´5 4 D 2 k V m µ 1 pS ˚q pD `k V m µ 1 pS ˚q N ˚q (12)
and

β " k V m D 2 µ 1 pS ˚q pD `k V m µ 1 pS ˚q N ˚q . ( 13
)
Some extensions of this Theorem are given in Section 3 such as the rate of convergence and non-monotonic growth rate. The last Theorem gives the heuristic that, until extinction and if n and t are sufficiently large, the discrete model is almost distributed as a normal distribution:

pN n t , S n t q « N ˜ˆn N S˚˙, ˜n α ´µpS ˚q µ 1 pS ˚q ´µpS ˚q µ 1 pS ˚q β n ¸¸. (14)
As one would expect, the number of bacteria is negatively correlated to the substrate rate: more individuals implies less food and vice versa. Recall that pS n t `km

V n N n t q tě0 is a martingale (i.e. the total mass is in mean conserved in the container).

This theorem can be understood as a first step to fully describe the Crump-Young model such as in the case of the logistic model described in [START_REF] Chazottes | Sharp asymptotics for the quasistationary distribution of birth-and-death processes[END_REF]. Indeed, in Section 3.2, we will see that, in large population, the quasi-stationary distribution of the Crump-Young model matches with the stationary distribution of its approximation. This is not trivial (and also not proved) because, for instance, the limits when n Ñ 8 and t Ñ 8 do not even commute! An example with a non-monotonic µ with different behaviors (several invariant measures, behavior depending on the initial conditions) is also presented.

Central limit theorems

2.1. Functional notations. For any n ě 1 and t ě 0, the population of bacteria is represented by the punctual measure ν n t "

ř N n t i"1 δ X i t .
We denote by Mpr0, M sq the set of such measures (punctual measures), it is a subset of the set M F pr0, M sq of finite (positive) measures.

For any n ě 1 and T ą 0, the process pν n t q tPr0,T s is a càd-làg process. It (almost-surely) belongs to the space Dpr0, T s, M F pr0, M sqq of càd-làg functions of r0, T s with values in M F pr0, M sq, endowed with the (usual) Skohorod topology; see for instance [START_REF] Billingsley | Convergence of Probability Measures[END_REF]; [START_REF] Ethier | Markov Processes -Characterization and Convergence[END_REF] for an introduction. In contrast, pS n t q tPr0,T s is (almost surely) a continuous function. We denote by Cpr0, T s, R `q the set of continuous functions from r0, T s to R `endowed with its usual topology.

The convergence (2) corresponds to a convergence in distribution in the product space Cpr0, T s, R `q ˆDpr0, T s, M F pr0, M sqq. Roughly this convergence is proved by a compactness/uniqueness argument. The compactness (or tightness) is proved by the (well-known) Aldous criterion which is a stochastic generalisation of the Arzelà-Ascoli Theorem. One of the key assumption of this theorem is to work in metric space. Considering the fluctuation process pη n t q tě0 , such arguments cannot be used to establish any convergence. Indeed, in contrast to the measure νn t , the measure η n t is not a positive measure but it is a signed measure. The set of signed measures being not metrisable [START_REF] Varadarajan | Weak convergence of measures on separable metric spaces[END_REF], one has to consider η n t as an operator acting on a different space than those of continuous and bounded functions. As [START_REF] Meleard | Convergence of the fluctuations for interacting diffusions with jumps associated with Boltzmann equations[END_REF] and [START_REF] Tran | Modèles particulaires stochastiques pour des problèmes d'évolution adaptative et pour l'approximation de solutions statistiques[END_REF], we will use some Sobolev spaces that are defined as follows: for every integer j, we let C j pr0, M sq be the set of functions being j times continuously differentiable endowed with the norm } ¨}C j , defined for all f P C j pr0, M sq by }f } C j " ř j i"0 }f piq } 8 (with } ¨}8 the infinity norm). Let now }¨} W j be the norm defined, for any f P C j pr0, M sq, by }f } 2 W j :"

ż r0,M s j ÿ i"0
pf piq pxqq 2 dx.

Let W j " W j pr0, M sq be the completion of C j pr0, M sq with respect to this norm (note that it is the classical Sobolev space associated to the classical norm }¨} L 2 ). Contrary to the Banach space C j pr0, M sq, the Sobolev space W j is Hilbertian. Let W j be its dual space, classically endowed with the norm

}µ} W j " sup }f } W j ď1 |µ pf q | .
Another useful property is given by the Sobolev-type inequalities: there exist universal constants

C j , C 1 j such that }f } W j ď C j }f } C j , }f } C j ď C 1 j }f } W j`1 . (15)
See for instance (Meleard, 1998, Equations (3.5) and (3.6)) or [START_REF] Adams | Sobolev spaces[END_REF], Theorem V-4). In particular, W j`1 is continuously embedded in W j . Moreover, this embedding is a Hilbert-Schmidt embedding (see (Meleard, 1998, Equation (3.7)) or (Adams, 1975, Theorem VI-53)). Therefore bounded and closed sets of W j`1 are compact for the W j 's topology.

Let us illustrate an application of inequalities (15) that will be useful in the proof of Theorem 1.2.

Lemma 2.1 (Useful bound on the basis). Let pe k q kě0 be an orthonormal basis of W 2 , we have,

K 1 :" sup xPr0,M s ÿ kě0 e 2 k pxq ă `8 .
Proof. By definition, for any x P r0, M s,

}δ x } W 2 " sup }f } W 2 ď1 |δ x pf q | " sup }f } W 2 ď1 |f pxq| ď sup }f } W 2 ď1 }f } 8 .
Moreover, we have the Sobolev-type inequalities, }f } 8 " }f } C 0 ď C }f } W 2 , for some C ą 0. Hence, by the Parseval identity,

C 2 ě sup xPr0,M s }δ x } 2 W 2 " sup xPr0,M s ÿ kě0 e 2 k pxq .
Finally, contrary to the models of [START_REF] Meleard | Convergence of the fluctuations for interacting diffusions with jumps associated with Boltzmann equations[END_REF]; [START_REF] Tran | Modèles particulaires stochastiques pour des problèmes d'évolution adaptative et pour l'approximation de solutions statistiques[END_REF], the fluctuation process (as the empirical measure) is not here a Markov process by itself. We have to consider the couple population/substrate to have a homogeneous dynamics. As a consequence, we will use a slightly larger space than those of [START_REF] Meleard | Convergence of the fluctuations for interacting diffusions with jumps associated with Boltzmann equations[END_REF]; [START_REF] Tran | Modèles particulaires stochastiques pour des problèmes d'évolution adaptative et pour l'approximation de solutions statistiques[END_REF]. Let

H 0 " W 2 ˆR, H " W 3 ˆR, (16)
be the Hilbert spaces endowed with the following norms : for pµ, Rq P H 0 or pµ, Rq P H,

}pµ, Rq} H 0 " b }µ} 2 W 2 `|R| 2 , }pµ, Rq} H " b }µ} 2 W 3 `|R| 2 .
Remark 2.2 (Weaker assumptions on pη n 0 q ně1 ). More generally, Theorem 1.2 holds for H 0 " W j ˆR and H " W j`1 ˆR with j ě 2 (the entire proof holds replacing W 2 , W 2 , W 3 and W 3 by W j , W j , W j`1 and W j`1 ). For j ą 2, the assumptions on pη n 0 q ně1 are weaker than for H 0 and H defined by (16), however the convergence result is also weaker.

Martingale properties.

The sequence of processes pps ν n t q tě0 q ně1 , defined by (1), can be rigorously defined as a solution of stochastic differential equations involving Poisson point processes; see (Campillo and Fritsch, 2014, Section 4). Instead of using this characterisation, we only need that it is solution to the following martingale problem.

Lemma 2.3 (Semi-martingale decomposition, [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF]). We assume that Eps ν n 0 p1q 2 q ă 8. Let f P C 1 pr0, M sq, then, under Assumptions 1.1, for all t ą 0:

s ν n t pf q " s ν n 0 pf q `ż t 0 ż r0,M s " bpS n u , xq ż 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ‰ Qpdαq ´D f pxq `gpS n u , xq f 1 pxq  s ν n u pdxq du `Zn t pf q (17) where d dt S n t " Dps in ´Sn t q ´k V ż r0,M s gpS n t , xq s ν n t pdxq
and pZ n t pf qq tě0 is a martingale with the following predictable quadratic variation:

xZ n pf qy t " 1 n ż t 0 ż r0,M s bpS n u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs 2 Qpdαq s ν n u pdxq du `1 n D ż t 0 ż r0,M s f pxq 2 s ν n u pdxq du .
Let us recall that s η n t " pη n t , R n t q is defined by (5). From the last lemma, we deduce that s η n t " s

A n t `Ď M n t with s A n t pf q :" ˆAn t pf q R n t ˙Ď M n t pf q :" ˆM n t pf q 0 ˙,
where the processes pA n t q tPr0,T s and pR n t q tPr0,T s have finite variations and are defined by

A n t pf q :" η n 0 pf q `ż t 0 "ż r0,M s bpS n u , xq ż 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ‰ Qpdαq η n u pdxq ´D η n u pf q `ηn u `gpS n u , .q f 1 ˘ du `?n ż t 0 " ξ u `pgpS n u , .q ´gpS u , .qq f 1 żr0,Ms pbpS n u , xq ´bpS u , xqq ż 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ı Qpdαq ξ u pdxq * du (18)
and

R n t " ż t 0 " ´D R n u ´k V η n u pgpS n u , .qq ´?n k V ξ u pgpS n u , .q ´gpS u , .qq  du. ( 19
)
The process pM n t pf qq r0,T s , defined by

M n t pf q :" ? n Z n t pf q , ( 20
)
is a martingale with predictable quadratic variation

xM n pf qy t " ż t 0 ż r0,M s bpS n u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs 2 Qpdαq s ν n u pdxq du `D ż t 0 ż r0,M s f pxq 2 s ν n u pdxq du . (21)
2.3. Preliminary estimates. For any n ě 1, let pxxM n yy t q tPr0,T s be the trace of the process pM n t q tPr0,T s , that is the process such that p}M n t } 2 W 2 ´xxM n yy t q tPr0,T s is a local martingale [START_REF] Joffe | Weak convergence of sequences of semimartingales with applications to multitype branching processes[END_REF][START_REF] Métivier | of de Gruyter Studies in Mathematics[END_REF][START_REF] Métivier | Convergence faible et principe d'invariance pour des martingales à valeurs dans des espaces de Sobolev[END_REF]. Let pe k q kě0 be an orthonormal basis of W 2 , we have

xxM n yy t " ÿ kě0 xM n pe k qy t . (22)
Indeed, on the one hand, the Parseval's identity entails that

}M n t } 2 W 2 " ÿ kě0 pM n t pe k qq 2 .
On the other hand, by definition of the predictable quadratic variation, pM n t pe k q 2 xM n pe k qy t q tě0 is a martingale for any k ě 0. Therefore, the uniqueness of the trace implies (22).

Lemma 2.4 (Uniform moment). If sup ně1 E ´}s η n 0 } 2 H 0 ¯ă 8 then for any T ą 0 sup ně1 E ˆsup uďT }s η n u } 2 H 0 ˙ă 8 .
Proof. By definition of }.} H 0 and the triangular inequality, we have

}s η n t } 2 H 0 " }η n t } 2 W 2 `|R n t | 2 ď 2 p}A n t } 2 W 2 `}M n t } 2 W 2 q `|R n t | 2 .
By the Fubini-Tonelli theorem,

E " sup uďt }M n u } 2 W 2  " E « sup uďt ÿ kě0 M n u pe k q 2 ff ď ÿ kě0 E " sup uďt M n u pe k q 2  .
Applying the Doob inequality (see for instance (Jacod and Shiryaev, 2003, Theorem 1.43)) to the martingale pM n u pe k qq uě0 for any k ě 0, we have

E " sup uďt }M n u } 2 W 2  ď 4 ÿ kě0 E " M n t pe k q 2 ‰ " 4 ÿ kě0 E rxM n pe k qy t s " 4 E « ÿ kě0 xM n pe k qy t ff .
By ( 21) and by Lemma 2.1,

ÿ kě0 xM n pe k qy t " ÿ kě0 ż t 0 ż r0,M s ! bpS n u , xq ż 1 0 re k pα xq `ek pp1 ´αq xq ´ek pxqs 2 Qpdαq `D e 2 k pxq ) s ν n u pdxq du ď `5 s b `D˘K 1 T sup uďt s ν n u p1q
.

By assumption sup n E `}s η n 0 } H 0 ˘ă 8, which implies that sup n E ps ν n 0 p1qq ă `8. Then, from (Campillo and Fritsch, 2014, Lemma 5.4)

sup n E " sup uďT s ν n u p1q ‰ ă `8. Hence K 2 :" sup n E " sup uďT }M n u } 2 W 2  ă `8 . ( 23
)
Using the Gronwall inequality, we easily check that

C ξ T :" sup tPr0,T s ξ t p1q ă `8 ; ( 24
)
it is however proved in [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF], Proof of Theorem 5.2). Therefore, from Assumptions 1.1 and the Sobolev-type inequalities,

}A n t } W 2 ď }η n 0 } W 2 `C ż t 0 }η n u } W 2 `s g `3s b `D˘d u `4 C ż t 0 |R n u | K gb C ξ T du ď }η n 0 } W 2 `K3 ż t 0 ´}η n u } W 2 `|R n u | ¯du , with K 3 :" C max s g `3s b `D; 4 K gb C ξ T (
and C depends on the constants C 1 and C 1 1 of the Sobolev-type inequalities. Hence

}A n t } 2 W 2 ď 2 ˆ}η n 0 } 2 W 2 `K2 3 T ż t 0 }s η n u } 2 H 0 du ˙.
Moreover,

|R n t | ď ż t 0 ˆD |R n u | `C k V s g }η n u } W 2 `k V |R n u | K gb C ξ T ˙du hence |R n t | 2 ď K 2 4 T ż t 0 }s η n u } 2 H 0 du , with K 4 " max ! C k V s g ; D `k V K gb C ξ T )
. Therefore, by the Fubini-Tonelli theorem

E " sup sďt }s η n s } 2 H 0  ď p4 K 2 3 `K2 4 q T ż t 0 E " sup sďu }s η n s } 2 H 0  du `4 E }s η n 0 } 2 H 0 `2 K 2 .
By Gronwall lemma, we finally get

sup ně1 E " sup sďt }s η n s } 2 H 0  ď e p4 K 2 3 `K2 4 q T t p2 K 2 `4 sup ně1 E }s η n 0 } 2 H 0 q ă `8 .
2.4. Proof of the Theorem 1.2. We divide the proof in two steps. The first one is devoted to the proof of the tightness of the sequence of processes ps η n q ně1 in Dpr0, T s, Hq.

In the second one, we prove that the limit of the process is unique and given by (6-7).

Step 1 : Tightness of ps η n q in Dpr0, T s, Hq. From (Meleard, 1998, Lemma C) or [START_REF] Joffe | Weak convergence of sequences of semimartingales with applications to multitype branching processes[END_REF], the sequence of processes pps η n q tPr0,T s q ně1 is tight in Dpr0, T s, Hq if the two following conditions hold:

rT s for all t ď T , sup

ně1 E " }s η n t } 2 H 0 ‰ ă `8;
rAs for any ε ą 0, α ą 0, there exist θ ą 0 and n 0 such that for any sequence pσ n , τ n q n of pairs of stopping times with σ n ď τ n ď σ n `θ,

sup něn 0 Pp › › s A n τn ´s A n σn › › H ě αq ď ε , sup něn 0
Pp|xxM n yy τn ´xxM n yy σn | ě αq ď ε .

Indeed recall that the embedding H 0 Ă H is Hilbert-Schmidt and then, using Markov inequality, rT s implies that the sequence ps η n t q n almost surely belongs to a bounded set of H 0 (which is compact in H). In short, rT s implies the tightness of ps η n t q ně0 for every t ě 0 in H.

In order to prove the tightness of ps η n q in Dpr0, T s, Hq, we have to prove the conditions rT s and rAs. Condition rT s is a direct consequence of Lemma 2.4. Let us now prove rAs. By the Markov inequality,

P `› › s A n τn ´s A n σn › › H ě α ˘ď E › › s A n τn ´s A n σn › › H α .
By (18), we have for any

f P W 3 Ă W 2 such that }f } W 3 ď 1, |pA n τn ´An σn q pf q | ď C ż σn`θ σn }η n u } W 2 ps g `3 s b `Dq du `C ż σn`θ σn |R n u | 4 K gb C ξ T du ď C ż σn`θ σn p}η n u } W 2 `|R n u |q du ď C θ sup uďT }s η n u } H 0
where the constant C can be different from a line to another. By the same way, |R n τn

´Rn σn | ď C θ sup uďT }s η n u } H 0 then Ep › › s A n τn ´s A n σn › › H q ď C θ sup ně1 E ˆsup uďT }s η n u } H 0 ˙.
By Lemma 2.4, the first condition of rAs is then satisfied.

In the same way, Step 2 : Identification of the accumulation points. From Step 1, the sequence ps η n q ně1 is tight in Dpr0, T s, Hq. Therefore, by Prokhorov's theorem, it is relatively compact and then we can extract, from ps η n q ně1 , a subsequence that converges weakly to a limit ps η t q tPr0,T s " pη t , R t q tPr0,T s P Dpr0, T s, Hq. We want to prove, in this step, that this limit is unique and defined by ( 6). Then, the theorem will follow (see for example [START_REF] Billingsley | Convergence of Probability Measures[END_REF], Corollary p.59)).For a better simplicity in the notations, we assume, without loss of generality that the entire sequence ps η n q ně1 converges towards the limit s η " pη, Rq.

|xxM
Lemma 2.5 (Convergence of the martingale part). The sequence of martingale processes pM n q n converges in distribution in Dpr0, T s, W 3 q towards a process G with values in Cpr0, T s, C 0,˚p r0, M sqq Ă Dpr0, T s, W 3 q, where C 0,˚p r0, M sq is the dual of C 0 pr0, M sq. For any f P C 0 pr0, M sq, the process G pf q is a continuous centred Gaussian martingale with values in R with quadratic variation defined by (7).

Proof. Let f P C 0 pr0, M sq and xG pf qy t be the quadratic variation defined by ( 7), then by (21),

|xM n pf qy t ´xG pf qy t | ď ż t 0 ż r0,M s " bpS n u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs 2 Qpdαq `D f pxq 2  ˆ|s ν n u pdxq ´ξu pdxq| du `ż t 0 ż r0,M s |bpS n u , xq ´bpS u , xq| ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs 2 Qpdαq ξ u pdxq du ď p9 s b `Dq }f } 2 8 ż t 0 ż r0,M s |s ν n u pdxq ´ξu pdxq| du `9 K gb }f } 2 8 C ξ T ż t 0 |S n u ´Su | du
with C ξ T defined by ( 24). Therefore, by (2) and the dominated convergence theorem, xM n pf qy converges in distribution towards xG pf qy.

Moreover, a discontinuity of t Þ Ñ ν t only happens during a birth or death event and the jump of the population number is ˘1. Then from (1) and ( 17), for any f P C 0 pr0, M sq,

|∆Z n t pf q| ď }f } 8 n . Therefore, from (20) sup tPr0,T s |∆M n t pf q | ď }f } 8 ? n , (25)
and then sup tPr0,T s |∆M n t pf q | converges in probability towards 0. Hence, according to (Jacod and Shiryaev, 2003, Theorem 3.11 page 473), for each f , the sequence of processes ppM n t pf qq tPr0,T s q ně1 converges to pG t pf qq tPr0,T s . To have an (infinite dimensional) convergence of the sequence of (operator valued) processes ppM n t q tPr0,T s q ně1 , it then suffices to prove its tightness in Dpr0, T s, W 3 q. To do it, it is enough to use (23) and arguments of the step 1.

Lemma 2.6 (Limit equation for pR t q). The limit process pR t q tPr0,T s satisfies

R t " ´ż t 0 " D R u `k V η u pgpS u , .qq `Ru k V ξ u pB s gpS u , .qq  du .
Proof. By definition, pR t q tPr0,T s is a limit point of the sequence of processes ppR n t q tPr0,T s q n , then by ( 19), we have the following limit in distribution : for any t P r0, T s,

R t " ´lim nÑ8 ż t 0 " D R n u `k V η n u pgpS n u , .qq `?n k V ξ u pgpS n u , .q ´gpS u , .qq  du . ( 26
)
By definition of s η " pη, Rq as a limit of ps η n q n and as the function which converges towards 0, in distribution, by Lemma 2.4. Hence, from Lemma 2.4 again and the dominated convergence theorem in (26), the conclusion follows.

s Þ Ñ gps, .q is contin- uous, from Lemma 2.4, D R n u `k V η n u pgpS n u , .qq converges in distribution towards D R u k V η u pgpS u ,
Lemma 2.7 (Semi-martingale decomposition). The process pM t q tPr0,T s defined for any f P W 3 by

M t pf q " η t pf q ´η0 pf q ´ż t 0 ż r0,M s " bpS u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs Qpdαq ´D f pxq `gpS u , xq f 1 pxq  η u pdxq du ´ż t 0 R u ż r0,M s " B s bpS u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs Qpdαq `Bs gpS u , xq f 1 pxq  ξ u pdxq du ( 27
)
has the same law as the process G defined in Lemma 2.5.

Proof. We define, for any ζ P Dpr0, T s, W 3 q, f P W 3 , t P r0, T s,

Ψ f t pζq " ζ t pf q ´ζ0 pf q ´ż t 0 ż r0,M s " bpS u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs Qpdαq ´D f pxq `gpS u , xq f 1 pxq  ζ u pdxq du ´ż t 0 R ζ u ż r0,M s " B s bpS u , xq ż 1 0 rf pα xq `f pp1 ´αq xq ´f pxqs Qpdαq `Bs gpS u , xq f 1 pxq  ξ u pdxq du ( 28
)
where

R ζ t " ż t 0 " ´D R ζ u ´k V ζ u pgpS u , .qq ´Rζ u k V ξ u pB s gpS u , .qq  du . (29)
Following, for example, the approach of (Campillo and Fritsch, 2014, Lemma 5.8), we can prove that ζ Þ Ñ Ψ f t pζq is continuous from Dpr0, T s, W 3 q to R in any point ζ P Cpr0, T s, W 3 q. Indeed, using some rough bounds, we have the existence of a constant C ą 0, such that

|Ψ f t pζq ´Ψf t p r ζq| ď C sup sPr0,T s › › ›ζ s ´r ζ s › › › W 3
and on continuous points, the Skorohod topology coincides with the uniform topology. However, as for (25),

|η n t pf q ´ηn t´p f q | ď }f } 8 ? n
therefore, η is a continuous process and then lim nÑ8 Ψ f t pη n q " Ψ f t pηq in distribution. By Lemma 2.5, it is sufficient to prove the proposition that pΨ f t pη n qq t and M n t pf q converge in distribution towards the same limit.

As η n t pf q " A n t pf q `M n t pf q, by ( 18) and ( 28), Ψ f t pη n q ´M n t pf q " B n t pf q `Cn t pf q with B n t pf q "

ż t 0 ξ u `"? n pgpS n u , .q ´gpS u , .qq ´Rη n u B s gpS u , .q ‰ f 1 ˘du `ż t 0 ż r0,M s "? n pbpS n u , xq ´bpS u , xqq ´Rη n u B s bpS u , xq ‰ ż 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ı Qpdαq ξ u pdxq du and C n t pf q " ż t 0 " η n u `pgpS n u , .q ´gpS u , .qq f 1 żr0,Ms pbpS n u , xq ´bpS u , xqq ż 1 0 " f pα xq `f pp1 ´αq xq ´f pxq ‰ Qpdαq η n u pdxq  du .
By the same approach using in the proof of Lemma 2.6, we get ˇˇ? n pgpS n u , .q ´gpS u , .qq ´Rη n u B s gpS u , .q ˇˇď ˇˇB s gpS u , .qpR n u ´Rη n u q ˇˇ`K gb 2 ? n pR n u q 2 (30)

and

ˇˇ? n pbpS n u , xq ´bpS u , xqq ´Rη n u B s bpS u , xq ˇˇď ˇˇB s bpS u , xqpR n u ´Rη n u q ˇˇ`K gb 2 ? n pR n u q 2 .
Hence,

sup tďT |B n t pf q| ď p › › f 1 › › 8 }B s g} 8 `3 }f } 8 }B s b} 8 q C ξ T T sup tďT |R n t ´Rη n t | `p› › f 1 › › 8 `3 }f } 8 q K gb C ξ T T 2 ? n sup tďT pR n t q 2 ,
where C ξ T was defined in (24). From Lemma 2.4, the second term converges towards 0 in probability. By ( 19), ( 29) and ( 30)

|R n t ´Rη n t | ď ˆD `k V }B s g} 8 C ξ T ˙ż t 0 |R n u ´Rη n u | `Kgb k V ż t 0 |S n u ´Su | η n u p1q du `k V K gb 2 ? n T C ξ T sup tďT |R n t | 2 .
From Lemma 2.4 and by the Gronwall lemma, we deduce that pR n t ´Rη n t q tďT converges uniformly towards 0 in probability and then that pB n t pf qq tďT converges uniformly towards 0 in probability. Furthermore,

sup tďT |C n t pf q| ď p › › f 1 › › 8 `3 }f } 8 q K gb T C 1 1 sup tďT }η n t } W 2 sup tďT |S n t ´St |
where that C 1 1 was defined in (15). The sequence S n converges in distribution towards S then by Lemma 2.4, we deduce that sup tďT |C n t pf q| converges in probability towards 0. Finally, pΨ f t pη n qq t and M n t pf q have the same limit G pf q in distribution. To conclude the proof of Theorem 1.2, it rests to prove the uniqueness of the solution of ( 27) but it can be easily proved via the classic argument involving Gronwall lemma.

2.5. Proof of Theorem 1.3. The proof is quite similar to the proof of Theorem 1.2 so we do not provide all details.

Firstly, similarly to Lemma 2.4, we can use Lemma 2.3 (with f " 1), Doob's inequality, Gronwall lemma and some rough bounds to show that

sup ně1 E ˆsup tďT |N n t | 2 `|S n t | 2 `|Q n t | 2 `|R n t | 2 ˙ă `8 . (31)
Indeed, note that one can bound µpS n s q by s µ " max 0ďsďs in _S 0 µpsq for every s ě 0, because S n remains in r0, s in _ S 0 s; see for instance (Collet et al., 2013a, Proposition 2.1).

Using Equation (31) and Markov inequality, we obtain, as in the proof of Theorem 1.2 that pN n , S n , Q n , R n q ně1 satisfies the Aldous Robolledo criterion (Joffe and Métivier, 1986, Corollary 2.3.3.) and then that pN n , S n , Q n , R n q ně1 is tight in Dpr0, T s, R `ˆr0, s in _ S 0 s ˆR2 q.

It remains to show the uniqueness of the limit point. Using Lemma 2.3, we can see that each limit point of the sequence pN n , S n q is solution to the classic chemostat ODE (i.e. the two first equations of ( 9)) and then by uniqueness of the solution, it converges to pN, Sq. Let now study a limit of a convergent subsequence of pQ n , R n q n . Following the way of Lemma 2.5 with f " 1, we obtain that sup tPr0,T s |∆M n t p1q | ď 1 ? n . Moreover, we have the following convergence in distribution, lim nÑ8 xM n p1qy tPr0,T s " ˆż t 0 pµpS u q `Dq N u du ˙tPr0,T s .

Then, by (Jacod and Shiryaev, 2003, Theorem 3.11 page 473), we deduce the convergence, in distribution, of ppM n t p1qq tPr0,T s q ně1 towards p ş t 0 a pµpS u q `Dq N u dB u q tPr0,T s . The end of the proof is then as in the proof of Theorem 1.2.

Remark 2.8 (Infinite dimensional case when M " 8). According to the proof of Theorem 1.3, we see that to obtain the convergence from finite M to infinite M (i.e. non compact support for the mass) for the finite-dimensional process then it is enough to prove that the uniform bound of Lemma 2.4 remains valid (that is equation 31).

The situation is more tricky in infinite dimension. Indeed, firstly, we crucially need the convergence of s ν n to ξ in the space of positive measure endowed with the weak topology in the proof of Theorem 1.2. [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF] proved the convergence in the space of positive measure endowed with the vague topology. Although vague and weak topologies coincide on a compact space, it is not valid anymore on non-compact sets. Extending the convergence to the weak topology is not trivial; see for instance [START_REF] Cloez | Limit theorems for some branching measure-valued processes[END_REF]; [START_REF] Méléard | Slow and fast scales for superprocess limits of agestructured populations[END_REF]. Also, note that Inequalities (15) are also no longer valid in infinite dimension.

The Crump-Young model

In this section, we propose an application of the previously demonstrated central limit theorem (Theorem 1.3) to understand the Crump-Young model. In particular, Section 3.1 contains the proof of Theorem 1.4, which gives an approximation of the long-time behavior of the Crump-Young model (see ( 14)). The process pS n t , N n t q tě0 satisfies the Markov property and is generated by the following infinitesimal generator Lf ps, q " " Dps in ´sq ´km nV µpsq  B s f ps, q `µpsq pf ps, `1q ´f ps, qq `D pf ps, ´1q ´f ps, qq , for all ě 0, s ě 0 and smooth f . This model is a particular case of the general model of [START_REF] Campillo | Weak convergence of a mass-structured individualbased model[END_REF], where we suppose that division rate and the growth rate (per capita) do not depend on the mass of the bacteria. This is a rough assumption which enables us to considerably weaken the dimension of the problem (from an infinite dimension to two dimensions). Our main result implies that it can be approximated by (10) and this diffusion process will be the main object of interest. Note that, through the function µ, we introduce a parameter m which can be understood as the mean size of one bacterium induced by the mass-structured model. Indeed, if we consider the integrodifferential equation ( 3) with parameters given by ( 8) and we set

@t ě 0, N t " ż r0,M s ξ t pdxq, Y t " ż r0,M s x ξ t pdxq,
where N t represents the number of individuals at time t and Y t the biomass. As pointed out by (Campillo and Fritsch, 2014, Section 5.4), one can prove that these two quantities can be described as a solution to the classic chemostat equations. Moreover d dt

Y t N t " µpS t q ˆm ´Yt N t ˙,
and then Yt Nt converges to m when t tends to infinity. Before studying rigorously the behavior of the system (9), let us end this section by a remark on the modelling.

Remark 3.1 (Reinforced process for indirect interactions). Consider the system (9), with starting points N 0 " N ˚, S 0 " S ˚, R 0 " 0, Q 0 P R, where pN ˚, S ˚q ‰ p0, s in q is some equilibrium of the two first equations. As µpS ˚q " D, the system then reduces to

# dQ t " µ 1 pS ˚q R t N ˚dt `?2 D N ˚dB t , dR t " ´"D R t `k V µpS ˚q m Q t `k V R t µ 1 pS ˚q m N ˚‰ dt.
In particular the second equation became a simple linear (ordinary) differential equation and then by the variation of constants method, we have

R t " ´k V m µpS ˚q ż t 0 e ´pD`k V m µ 1 pS ˚q N ˚q pt´sq Q s ds.
Hence

dQ t " ´µ1 pS ˚q N ˚k V m µpS ˚q ż t 0 e ´pD`k V m µ 1 pS ˚q N ˚q pt´sq Q s ds dt `?2 D N ˚dB t .
The solution of this equation then represents the evolution of the population around an equilibrium under an indirect competition (presence of substrate). This process belongs to the large class of self-interacting diffusions; see [START_REF] Gadat | Long time behaviour and stationary regime of memory gradient diffusions[END_REF]; [START_REF] Gadat | A stochastic model for speculative dynamics[END_REF] and reference therein. These processes are not Markov and if, more generally,

dQ t " ´ż t 0 κpt ´sq Q s ds dt `?2 D N ˚dB t ,
for some function κ, then κ represents the memory of the substrate consumption. In a different context than the chemostat, one can imagine a different function κ to model an indirect interaction which can influence the size of the population.

3.1. Proof of theorem 1.4. In this section, we study the solution of the system of equations ( 9) under the assumptions of Theorem 1.4. Firstly, let us see that the two first equations of (9) forms a homogeneous system of ODE. It is the classic chemostat equations; see [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]. In particular, as the specific growth rate µ is supposed to be increasing, the couple pN t , S t q tě0 admits only two equilibria that are p0, s in q, which is usually called the washout and corresponds to the extinction of the population, and another pN ˚, S ˚q corresponding to the unique solution of µpS ˚q " D and N ˚" V km ps in ´S˚q .

Moreover, we have

d ˆSt `mk V N t ˙" D ˆsin ´ˆS t `mk V N t ˙˙dt,
and then

lim tÑ8 S t `mk V N t " s in .
Also a calculus of the Jacobian at these two points shows that pN ˚, S ˚q is stable while p0, s in q is unstable. As a consequence, the Poincaré-Bendixson theorem (see for instance (Smith and Waltman, 1995, Page 9)) entails that, whatever the initial condition pN 0 , S 0 q P R ˚ˆR `, the following deterministic convergence holds: lim tÑ8 pN t , S t q " pN ˚, S ˚q.

Now, let us study the dynamics of pQ t , R t q tě0 . We set Z t " pQ t , R t q T , then

dZ t " A t Z t dt `Ct dB t , (32)
where

A t " ˆµpS t q ´D µ 1 pS t q N t ´k V µpS t q m ´pD `k V µ 1 pS t q N t mq ˙, C t " ˆapµpS t q `Dq N t 0 ˙.
In particular, one can think pN t , S t , Q t , R t q tě0 as a homogeneous-time Markov process or only pQ t , R t q tě0 as an inhomogeneous one. Equation ( 32

Therefore, for all t ě 0, the law of Z t is a Gaussian distribution of mean e ş t 0 As ds EpZ 0 q and variance matrix Σ t given by Σ t :"

ż t 0 e ş t s Au du C s C T s e ş t s A T u du ds.
To prove the convergence in law of Z t to a Gaussian variable, it is then enough to study the convergence of its mean and its variance. Note that the eigenvalues of A s are

λ 1 s " µpS s q ´D ´k V m µ 1 pS s q N s , λ 2 s " ´D,
which are (at least for large s because pN s , S s q Ñ pN ˚, S ˚q) negative because µpS ˚q " D and µ 1 ą 0. Nevertheless, it does not directly imply the convergence of the mean; see for instance (Amato, 2006, Example 2.2). However, we have

A 8 :" lim tÑ8 A t " ˆ0 µ 1 pS ˚q N k V µpS ˚q m ´pD `k V µ 1 pS ˚q N ˚mq ˙,
whose eigenvalues are λ `" ´k V m µ 1 pS ˚q N ˚and λ ´" ´D, and then, by a Cesàro-type theorem, [START_REF] Amato | Robust control of linear systems subject to uncertain time-varying parameters[END_REF], Theorem 2.9), we have lim tÑ8 e ş t 0 As ds " 0.

We have then obtained the convergence of the mean, it rests to prove the convergence of the variance matrix to

Σ 8 " ż 8 0 e A8 u C 8 C T 8 e A T 8 u du,
where

C 8 :" ˆ?2DN 0 ˙" lim tÑ8 C t .
Again by [START_REF] Amato | Robust control of linear systems subject to uncertain time-varying parameters[END_REF], Theorem 2.9), there exist C, α ą 0 such that

@t ě 0, › › ›e ş t 0 As ds › › › `› › ›e ş t 0 A8 ds › › › ď
› › › ds `ż 8 τ › › ›e A8 s C 8 C T 8 e A T 8 s › › › ds ď ż τ 0 › › ›e ş t t´s Au du C t´s C T t´s e ş t t´s A T u du ´eA8 s C 8 C T 8 e A T 8 s › › › ds `Ke ´2 α τ .
The introduction of the variable τ allows us to obtain an integral whose integration interval does not depend on t. As the second term of the last member is negligible for large τ , by dominated convergence, it then remains to prove that the last integrand vanishes when t Ñ 8. This is a direct application of the convergences of pA t q tě0 , pC t q tě0 and the continuity of the various applications (exponential, product...). Also note that ş t t´s A u du " ş s 0 A t´u du. This concludes the proof of the convergence of pΣ t q tě0 and then of the convergence of pN t , S t , Q t , R t q tě0 . Let us finally express the calculus of Σ 8 . We have e A8 s " 1 pD ´Lq ˆpD e ´L s ´L e ´D s q µ 1 pS ˚q N ˚pe ´L s ´e´D s q k V m D pe ´D s ´e´L s q D e ´D s ´L e ´L s ẇhere L " k V m µ 1 pS ˚q N ˚" ps in ´S˚q µ 1 pS ˚q, and then

e A8 s C 8 C T 8 e A T 8 s " 2 D N pD ´Lq 2 pD e ´L s ´L e ´D s q 2 k V m D pD e ´L s ´L e ´D s qpe ´D s ´e´L s q k V m D pD e ´L s ´L e ´D s qpe ´D s ´e´L s q `k V m D ˘2 pe ´D s ´e´L s q 2 ¸.
As a consequence, the term

ş 8 0 e A8 s C 8 C T 8 e A T 8 s ds is equal to 2 D N pD ´Lq 2 ˆ˜D 2 2 L `L2 2 D ´2 D L D`L k V m D " 1 ´D 2 L ´L 2 D ‰ k V m D " 1 ´D 2 L ´L 2 D ‰ p k V m Dq 2 " 1 2 D `1 2 L ´2 D`L ı ¸. Finally Σ 8 " D N L ˆ˜L 2 `3 L D`D 2 D pD`Lq ´k V m ´k V m p k V mq 2 D pD`Lq " ¨p k V m µ 1 pS ˚q N ˚`3 2 Dq 2 ´5 4 D 2 k V m µ 1 pS ˚q pD`k V m µ 1 pS ˚q N ˚q ´D µ 1 pS ˚q ´D µ 1 pS ˚q k V m D 2 µ 1 pS ˚q pD`k V m µ 1 pS ˚q N ˚q '.
Remark 3.2 (Rate of convergence). Due to the simple form of (32), one can give some estimates on the rate of convergence. Let W 2 be the (second order) Wasserstein distance, defined for any probability measure µ, ν by

W 2 pµ, νq " inf Er}X ´Y } 2 2 s 1{2
, where } ¨}2 is the classic Euclidean norm in R 2 , and the infimum runs over all random vectors pX, Y q with X " µ and Y " ν. Using (Givens and Shortt, 1984, Proposition 7), we find, for any t ě 0,

W 2 pLpZ t q, N p0, Σ 8 qq " › › ›e ş t 0 Asds › › › 2 2 `Tr ´Σt `Σ8 ´2 Σ 1{2 t Σ 8 Σ 1{2 t ¯,
where Tr is the classic trace operator. The decay of the right-hand side depends on the rate of convergence of the two-component ODE towards pN ˚, S ˚q. However, even if we assume that pN 0 , S 0 q " pN ˚, S ˚q, one can not simplify this expression because, even in this case, Σ t " ş t 0 e A8 s C 8 C T 8 e A T 8 s ds does not necessary commute with Σ 8 . The bound of (Givens and Shortt, 1984, Proposition 7) also induces a bound in Wasserstein distance for the four-component process pS t , N t , Q t , R t q tě0 . This is not trivial because it is not the case, for example, in total variation in contrast to pZ t q tě0 . Remark 3.3 (An example of non-increasing growth rate). Let us consider the following growth rate:

µ : s Þ Ñ µ max s K `s `s2 {C ,
where µ max , K, C are some positive constants. This rate is often called Haldane kinetics in the literature and can sometimes be more realistic in application; see for instance [START_REF] Mailleret | Nonlinear adaptive control for bioreactors with unknown kinetics[END_REF]. We assume that D ą µps in q and sup sPr0,s in s µpsq ą D.

In this case there are two solutions of µpSq " D , N " V km ps in ´Sq , (see Figure 1). Let us denote by p0, s in q, pN ˚, S ˚q, pN ue , S ue q the tree equilibria for pN t , S t q tě0 , with µ 1 pS ˚q ą 0 and µ 1 pS ue q ă 0. The study of the Jacobian matrix and the Poincaré-Bendixson theorem implies here that, if the ODE system does not start from the unstable equilibrium pN ue , S ue q then it necessary converges to one of the two stable equilibria pN ˚, S ˚q or p0, s in q, depending on the initial condition. As a consequence, the set of invariant distributions of the process pN t , S t , Q t , R t q tě0 is the convex hull of the Gaussian distributions (11) (with the stable equilibrium pN ˚, S ˚q for the Haldane growth) and δ p0,s in ,0,0q . Indeed, for the stable equilibrium p0, s in q, the matrix A 8 is zero as well as the vector C 8 (when we replace pN ˚, S ˚q by p0, s in q) while Σ t explodes for pN ue , S ue q because A 8 admits as positive eigenvalue ´k V mµ 1 pS ue qN ue (again when we replace pN ˚, S ˚q Table 1. Simulation times, in seconds, of the simulations of Figure 2 (15 runs of the Crump-young model and the SDE in small, medium and large population sizes). Simulations were carried out on a laptop computer with 2.5 GHz i5 (2 cores) processor and 4 GB memory.

Population size

Small Medium Large 15 15 runs of the SDE 4.631 4.573 4.272 by pN ue , S ue q). Also, mimicking the previous proof gives the convergence to one of them according to the starting distribution.

3.2. Numerical simulations and discussion. We use a Gillespie algorithm for the simulation of the Crump-Young model (C-Y) (see Algorithm 1 of [START_REF] Fritsch | A modeling approach of the chemostat[END_REF]) and an Euler method for the simulation of the stochastic differential equations (SDE) (10).

The system of ordinary differential equations (ODE) (two first equations of ( 9)) is solved by the odeint function of the scipy.integrate module of Python.

In general, a chemostat is described by the substrate and the biomass concentrations rather than the substrate concentration and the number of individuals. The biomass concentration is obtained by multiplying the number of individuals by m V , therefore, the graphs are the same up to the multiplicative constant m V . 3.2.1. Monod growth. We use the Monod growth parameters of the Escherichia coli bacteria in glucose with a temperature equals to 30 degrees Celsius [START_REF] Monod | La technique de culture continue, théorie et applications[END_REF], i.e.

µpSq " 1.35 S 0.004 `S , k " 0.23 , and with m " 7 .10 ´13 g, D " 0.5 h ´1, S 0 " s in " 0.003 g.l

´1

The convergence, in large population size, of the Crump-Young model towards the SDE (10) is illustrated in Figure 2. In small population size, the behavior of the Crump-Young C-Y SDE ( 10) ODE small population size medium population size large population size

V " 10 ´8 l , N 0 " 5 V " 10 ´6 l , N 0 " 500 V " 10 ´5 l , N 0 " 5000 Figure 2
. Time evolution of the number of individuals (top) and the substrate concentration (bottom) for small, medium and large population sizes for the Monod growth model. The magenta lines are simulations of 15 independent runs of the Crump-Young model, the blue lines are simulations of 15 independent runs of the system of stochastic differential equations (10) and the large red line is the simulation of the system of ordinary differential equations given by the two first equations of (9).

model is different from the one of the SDE. In particular, contrary to the Crump-Young model, the SDE can not depict the population extinction. Moreover, in small population size, we observe that the number of individuals can be negative for the SDE, therefore this model is not satisfactory in this situation. Also note that the Crump-Young model is a jump model, whereas the SDE is a continuous model. However, in large population size, the jumps of the number of individuals (˘1) in the Crump-Young model become negligible with respect to the population size, then this model can be approximated by a continuous one. According to Figure 2, the SDE seems to be a good approximation of the Crump-Young model from medium population size. Moreover it is much faster to compute than the Crump-Young model (see Table 1). In very large population, both models converge to the deterministic system of ODE, given by the two first equations of (9), then the ODE model is sufficient to describe the behavior of the chemostat in this context.

Figure 3 1 compares the estimated quasi-stationary distribution of the Crump-Young model to the invariant mesure of the SDE given in Theorem 1.4 for the three population α β " ´0.942470). On the top (resp. right) of each graph, the blue histogram represents the empirical distribution of the number of individuals (resp. the substrate concentration) of the Crump-Young model, the dashed green line is this distribution regularized with a Gaussian kernel and the red curve represents the probability density function of the Gaussian law N pN ˚, αq (resp. N pS ˚, βq), with α defined by (12) (resp. β defined by ( 13)), where pN ˚, S ˚q is the non-trivial (‰ p0, s in q) equilibrium of the ODE system (see ( 14)). sizes of Figure 2. In small population size, we observe that the two laws are different. The main reason is the large probability of extinction of the Crump-Young model. Indeed, on the 1000 non-extinct populations, many are close to the extinction p0, s in q, whereas the invariant measure predicts a convergence in a neighbourhood of the non-trivial stationary state pN ˚, S ˚q ‰ p0, s in q. However, in medium and large population sizes, the invariant mesure ( 14) is a very good approximation of the quasi-stationary distribution of the Crump-Young model. The behavior of the chemostat, for Haldane growth, depends on the initial condition. Indeed, there is, for the ODE, two basins of attraction which are associated to the two stable equilibria p0, s in q and pN ˚, S ˚q (see Figure 1), contrary to Monod growth for which there is only one stable equilibrium (the washout is an unstable equilibrium).

If the initial condition is close to the boundary of the two basins of attraction, the ODE remains in its initial basin and converges to its attractor whereas, due to the randomness, the Crump-Young model can change basin of attraction. The SDE ( 10) is very depending on the ODE solution and will converge to the invariant mesure of the basin of attraction associated to the initial condition. Therefore, the SDE ( 10) is not representative of the two possible convergences for one given initial condition. The SDE ( 10) is in fact a good approximation of the Crump-Young model when the population size is sufficiently large (which depends on the distance between the initial condition and the boundary of the two basins of attraction) to ensure that the Crump-Young model does not change (with a large probability) basin of attraction (see Figure 4). Even if the approximation only holds for large population, in Figure 4 (right), both models converge to the population extinction (even if the SDE is not absorbed, it converges to 0).

However, if the object of interest is the convergence towards p0, s in q or pN ˚, S ˚q for one given initial condition (close to the boundary of basins of attraction) then we must either use the Crump-Young model (if the simulation time is reasonable) or use a model which keeps more qualitative properties than the SDE (10).

In fact, Theorem 1.3 suggests that, for n sufficiently large, the Crump-Young model can be approximated by pN n t , S n t q « p r N n t , r S n t q with r N n t :" n N t `?n rQ t `F n t s , r S n t :" S t `1 ? n rR t `Hn t s , where pF n q n and pH n q n are two sequences of processes which converge in distribution towards the process 0 in Dpr0, T s, Rq. The SDE (10) is obtained by letting F n " H n " 0. Let now consider pF n q n and pH n q n be defined by dF n t " " ? n pµpS n t q ´µpS t qq 10) (blue lines) and the ODE (red curve) for initial conditions N 0 " 1000, S 0 " 0.08 g.l ´1 (left) and N 0 " 3000, S 0 " 0.04 g.l ´1 (right). and

N n t n ´µ1 pS t q R t N t  dt `«c pµpS n t q `Dq
dH n t " k V m
" µpS t q pQ t ´Qn t q `Rt µ 1 pS t q N t ´?n pµpS n t q ´µpS t qq

N n t n  dt ,
with initial condition F n 0 " H n 0 " 0. Following, for example, the approach used in the proof of Lemma 2.6, we can prove that pF n q n and pH n q n converge towards 0 in distribution. We then (heuristically) obtain the following model of approximation :

$ & % d r N n t " pµp r
S n t q ´Dq r N n t dt `bpµp r S n t q `Dq r N n t dB t , d r S n t " " D ps in ´r S n t q ´k V n m µp r S n t q r N n t ı dt .

(

) 34
This new approximation model can be seen as a particular case of the model of [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF]; see (Campillo et al., 2011, Equations (17a) and (17b)) with K 1 " K 4 " 1 and K 2 " K 3 " K 5 " `8 (note that K 2 " K 3 " K 5 " `8 corresponds to a continuous approximation of the substrate equation for a large number of substrate particles which is an approximation that we made for all models in this article).

In contrast with (10), the variance of the population size r N n t depends on itself. Moreover this type of dependence is classic in population dynamics, see for instance [START_REF] Bansaye | ume 1 of Mathematical Biosciences Institute Lecture Series[END_REF].

Figure 5 represents the time evolution of the number of individuals for the four models (ODE, Crump-young models, SDE (10) and SDE (34)) in three cases. Each column represents the same simulation with all or some represented curves (the first line allows to compare the four models together, however we have split each graph in two graphs for the sake of clarity). The first one (on the left) is for initial condition close to the unstable equilibrium pN ue , S ue q for which the solution of the ODE converges towards the washout p0, s in q. Therefore the solutions of the SDE (10) also converge towards the washout. However the Crump-Young model changes basin of attraction with a large probability and converges either to the washout or to a neighbourhood of the stable equilibrium N 0 " 190, S 0 " 0.07 g.l ´1 N 0 " 210, S 0 " 0.065 g.l ´1 N 0 " N ue " 200, S 0 " S ue " 0.0656 g.l

´1

V " 10 ´9 l V " 10 ´9 l V " 10 ´9 l

Figure 5. Time evolution of the number of individuals for the Haldane growth model for 20 independent runs of the Crump-Young model (magenta lines), 20 independent runs of the SDE (34) (green lines), 20 independent runs of the SDE (10) (blue lines) and the ODE (red curve) for initial conditions N 0 " 190, S 0 " 0.07 g.l ´1 (left), N 0 " 210, S 0 " 0.065 g.l ´1 (middle) and N 0 " N ue , S 0 " S ue g.l ´1 (right).

pN ˚, S ˚q. As we can observe, the solutions of the SDE (34) mimic the behavior of the Crump-Young model and then SDE (34) seems better than the SDE (10) in this context. The second case (center) is for initial condition close to the unstable equilibrium pN ue , S ue q for which the solution of the ODE converges towards the stable equilibrium pN ˚, S ˚q. Once again, the Crump-Young model and the diffusion process (34) depict two possible convergences (towards the washout or the quasi-stationary distribution in the area of pN ˚, S ˚q) while the diffusion process (10) follows the solution of the ODE. We see an explosion of the noise for the blue curve. This comes from that for small time t, the matrix A t has large positive eigenvalue (due to the initial condition, recall that it tends to infinity when the initial condition is the unstable equilibrium) but, as in the Monod case, Σ t converges to a finite matrix.

For the last case (right), the initial condition equals the unstable equilibrium pN ue , S ue q. Therefore, the deterministic approximation stays at this equilibrium whereas the Crump-Young model and the diffusion process (34) depict over again the two possible convergences. We observe that the solutions of (10) diverge. In fact, we can write the SDE as in (33), but, as µ 1 pS ue q ă 0, the eigenvalue λ 1 s " ´k V m µ 1 pS ue q N ue of A s is positive which implies the divergence of the SDE (10).

Even if the solutions of (34) seems, for some parameters, to be a more suitable approximation for the Crump-Young model, it is nevertheless more difficult to study it mathematically. Indeed, as for the Crump-Young model, there is always extinction; see Theorem 3.4. Also, even if this process is continuous and solution to a stochastic differential equation, it is not possible to deduce a result of uniqueness (or convergence) for a quasi-stationary distribution because it is not reversible in contrast with the classic logistic diffusion process; see [START_REF] Collet | Quasi-stationary distributions[END_REF]. Also, in contrast with the solutions of (10), no explicit formula is known for solutions of (34).

3.3.

Extinction time of the diffusion process (34). In this section, we will consider a solution p r N t , r S t q tě0 of (34) for one fixed n. The notation P px,sq refers to the probability given the initial condition is p r N 0 , r S 0 q " px, sq and E px,sq is the expectation associated to this probability.

Theorem 3.4 (Extinction). Let p r N t , r S t q tě0 be a solution of (34) and T 0 " inftt ě 0 | r N t " 0u.

Then P px,sq pT 0 ă `8q " 1 for any starting distribution px, sq P R 2 `. Moreover, there exist C, α ą 0 such that for all px, sq P R 2 `and t ě 0, P px,sq pT 0 ą tq ď Ce ´αt px `s `1q.

Proof. First, we assume that for every compact set K Ă R 2 `, there exist t 0 , δ ą 0 such that δ " inf px,sqPK P px,sq pT 0 ă t 0 q ą 0.

Secondly, considering V 0 : pN, Sq Þ Ñ km V n N `S ´sin and using (34), pe Dt V 0 p r N t , r S t qq tě0 is a martingale; namely V 0 is a Lyapunov-type function. From (36) and the Lyapunov property, it is then classic to prove the statement of the Lemma. Indeed, shortly, the Lyapunov property entails that, whatever the initial position is, the process converges rapidly in a compact set and then, by (36), it will be absorbed in finite time. This standard argument to prove geometric ergodicity of general Markov processes is given, for instance, by [START_REF] Hairer | Yet another look at Harris' ergodic theorem for Markov chains[END_REF]. Nevertheless we can not directly apply this theorem because even if St Ñ s in , it does not reach it, therefore we can not obtain the convergence of St to s in in total variation. So let us prove that the Lyapunov property and (36) are sufficient to ensure (35). Let us fix a compact set K Ă R 2 `such that, for all x R K, V 0 pxq ě C 0 , for some C 0 ą 0. Moreover let us fix the associated t 0 and δ as in (36).

We divide the proof in two steps.

Step 1 : Bound on hitting time. Let τ be the hitting time of K. Using the stopping-time theorem, for any n P N and px, sq R K, we have E px,sq re Dτ ^nV 0 p r N τ ^n, r S τ ^nqs " V 0 px, sq.

Then

C 0 E px,sq re Dpτ ^nq s ď V 0 px, sq.

Indeed, if p r N 0 , r S 0 q R K then p r N t , r S t q R K for all t ď τ (by definition of τ ). Using the monotone convergence theorem, we have for every px, sq R K E px,sq re Dτ s ď C ´1 0 V 0 px, sq . Moreover, if px, sq P K, then τ " 0 hence E px,sq re Dτ s ď 1 therefore for any px, sq P R 2 Èpx,sq re Dτ s ď C ´1 0 V 0 px, sq `B (37)

with B " C ´1 0 s in `1 (because V 0 px, sq `sin ě 0). Then the Markov inequality gives P px,sq pτ ě tq " P px,sq pe Dτ ě e Dt q ď e ´Dt E px,sq re Dτ s ď pC ´1 0 V 0 px, sq `Bq e ´Dt .

(38)

Step 2 : Bound on the extinction time. Let s 0 " 0 and for every ě 0, τ `1 " infts ě s |p r N s , r S s q P Ku ´s s `1 " s `τ `1 `t0 Finally, this gives the existence of a constant C 1 ą 0 (which depends on t 0 but not on θ) such that P px,sq pT 0 ą tq ď C 1 e ´Dp1´θqt pC ´1 0 V 0 px, sq `Bq ÿ ě0 ´p1 ´δq θ pBe Dt 0 q 1´θ ¯ .

Choosing θ sufficiently close to 1 to guarantee that p1 ´δq θ pBe Dt 0 q 1´θ ă 1 ends the proof of ( 35) and then of the statement of the lemma.

It remains to prove (36) to end the proof. Let us introduce Ψ : px, sq Þ Ñ 2 ? x a µpsq `D .

We now set U t " Ψp r N t , r S t q for all t ě 0. This new process hits 0 at the same time T 0 as p r N t q ě0 and, using Itô formula, it verifies dU t " " pµp r S t q ´Dq r N t B x Ψp r N t , r S t q `1 2 pµp r S t q `Dq r N t B xx Ψp r N t , r S t q `ˆD ps in ´r S t q ´k V n m µp r S t q r N t ˙Bs Ψp r N t , r S t q  dt `bpµp r S t q `Dq r N t B x Ψp r N t , r S t qdB t " " 1 2 pµp r S t q ´DqU t ´1 4U t ´Ut 2pµp r S t q `Dq µ 1 p r S t qD ps in ´r S t q `U 3 t 8 µ 1 p r S t q k V n m µp r S t q ff dt `dB t .

One can then bound the drift term with quantities not depending on the substrate rate and use (Ikeda and Watanabe, 1981, Theorem 1.1 chapter VI) to see that U t ď r Z t for every t ě 0, where p r Z t q tě0 is the one-dimensional diffusion solution to d Zt " Cp r Z 3 t ´1 r Z t qdt `dB t for some constant C ą 0. By the Feller's test for explosions (see (Karatzas and Shreve, 1991, Chapter 5)), and a monotonicity argument, we deduce that, for all z 0 ą 0, there exists t 0 ą 0 such that inf

zďz 0 Pp r T 0 ă t 0 | r Z 0 " zq ą 0, (39)
where r T 0 " inftt ě 0 | r Z t " 0u. More precisely, let S " inftt ě 0 | r Z t R p0, `8qu be the exit time from p0, `8q. The scale function p defined in (Karatzas and Shreve, 1991, Equation (5.42)) is given by Using standard results of asymptotic analysis, we have ż y 1 e az 4 z ´bdz " `8 e ay 4 y ´b ˆ1 4ay 3 .

p : x Þ Ñ c
Then as ş x 1 y ´3dy ă `8, we have that lim xÑ`8 vpxq is finite. Moreover, for x P p0, 1q, we have 0 ă vpxq " e a |b ´1| then lim xÑ0 `vpxq is also finite (note that even if the case b " 1 is not treated in the previous line, it works as well). As a consequence by (Karatzas and Shreve, 1991, Proposition 5.32 (i)), the stopping time S is finite (and even integrable) and by (Karatzas and Shreve, 1991, Proposition 5.22 (d)) PpS " r T 0 | r Z 0 " z 0 q ą 0, for every z 0 ą 0. Consequently, for every z 0 ą 0, there exits t 0 ą 0 such that Pp r T 0 ă t 0 | r Z 0 " z 0 q ą 0 and then using that for all z ď z 0 , Pp r T 0 ă t 0 | r Z 0 " zq ě Pp r T 0 ă t 0 | r Z 0 " z 0 q, we have proved (39). Finally (36) is a direct consequence of (39).

Remark 3.5 (Quasi-stationary distribution). Equation (35) is a necessary (but not sufficient) condition to ensure existence of a quasi-stationary distribution; see for instance [START_REF] Collet | Quasi-stationary distributions[END_REF].

Remark 3.6 (Extinction of the Crump-Young model). It is not difficult to see that (36) and the Lyapunov property also hold for the Crump-Young model and then (35) also holds for this process. In particular this gives a new proof of (Collet et al., 2013a, Theorem 3.1). Moreover, in contrast to (Collet et al., 2013a, Theorem 3.1), we obtain the speed of extinction (35); furthermore we do not assume any monotonicity on µ.

Figure 1 .

 1 Figure 1. Equilibria of the substrate concentration with respect to the Haldane specific growth rate.

1Figure 3 .

 3 Figure3. Distribution of the number of individuals and the substrate concentration, at time T " 1000 h, for small (top left), medium (top right) and large (bottom) population sizes for the Monod growth model. For each population size, the blue crosses represent the states of 1000 independent runs of the Crump-Young model. The green dashed ellipse is the 95% confidence ellipse of a bidimensional normal variable where the mean and the covariance matrix are estimated on the 1000 Crump-Young model simulations (the sample correlation between the substrate concentration and the number of individuals is indicated under each graph). The red ellipse is the 95% confidence ellipse of a normal variable with law (14) (the theoretical correlation equals ´µpS ˚q µ 1 pS ˚q ?

Figure 4 .

 4 Figure 4. Time evolution of the number of individuals for the Haldane growth model for 20 independent runs of the Crump-Young model (magenta lines), 20 independent runs of the SDE (10) (blue lines) and the ODE (red curve) for initial conditions N 0 " 1000, S 0 " 0.08 g.l ´1 (left) and N 0 " 3000, S 0 " 0.04 g.l ´1 (right).

ż x 1 y

 1 b e ´ay 4 dy, for some a, b, c ą 0, and then lim xÑ0 `ppxq and lim xÑ`8 ppxq are clearly finite. Moreover, the function v (defined in(Karatzas and Shreve, 1991, Equation (5.65))) verifies v :

  The functions ps, xq Þ Ñ gps, xq and ps, xq Þ Ñ bps, xq are Lipschitz continuous w.r.t. s uniformly in x and differentiable in s with derivative Lipschitz continuous w.r.t. s uniformly in x: for all s 1 , s 2 ě 0, x P r0, M s, |gps 1 , xq ´gps 2 , xq| ď K gb |s 1 ´s2 |; |bps 1 , xq ´bps 2 , xq| ď K gb |s 1 ´s2 |; |B s gps 1 , xq ´Bs gps 2 , xq| ď K gb |s 1 ´s2 |; |B s bps 1 , xq ´Bs bps 2 , xq| ď K gb |s 1 ´s2 |.

  By (21) and Lemma 2.1, we then get |xxM n yy τn ´xxM n yy σn | ď C K 1 yy τn ´xxM n yy σn | ď C K 1 θ sup

	ż σn`θ
		s ν n u p1q du
	σn
	therefore
	E |xxM n ně1	uďT E ˆsup	}s η n u } H 0	ȧnd
	by Lemma 2.4, the condition rAs holds.

n yy τn ´xxM n yy σn | " ˇˇˇˇÿ kě0 `xM n pe k qy τn ´xM n pe k qy σn ˘ˇˇˇď ÿ kě0 |xM n pe k qy τn ´xM n pe k qy σn |

  Ce ´αt , where }.} is the standard matrix norm. Hence, there exists a constant K such that for any τ and t ě τ ,

	}Σ t ´Σ8 } "	› › › › ż t 0	e ş t t´s Au du C t´s C T t´s e ş t t´s A T u du ds	´ż 8 0	e A8 s C 8 C T 8 e A T 8 s ds	› › › ›
	ď	ż τ 0	› › ›e ş t t´s Au du C t´s C T t´s e	ş t t´s A T u du ´eA8 s C 8 C T 8 e ´A8 s	› › › ds
		`ż 8	› › ›e	ş t t´s Au du C t´s C T t´s e	ş t t´s A T u du
			τ

  3.2.2. Haldane growth. We now use the following Haldane growth function:and k " 0.23, m " 7 .10 ´13 g, D " 0.5 h ´1, s in " 0.0978 g.l ´1.

	µpSq " 1.35	S 0.004 `S `S2 {0.04	,

  . Let θ P r0, 1s, by Hölder inequality, we have P px,sq pT 0 ą tq " P px,sq p r N On the first hand and if ě 1, by the strong Markov property, Equation (36) and an induction argument, we have P px,sq p r N ´δqP px,sq p r N s ´1`τ ‰ 0q ď p1 ´δqP px,sq p r N s ´1 ‰ 0q ď p1 ´δq . On the other hand, by the Markov property, Equation (38), the martingale properties (stopping time theorem on a truncated version of s and Fatou Lemma) and noting that τ " τ 1 , P px,sq pt P rs , s `1qq ď E px,sq " 1 těs P p r Dt 0 e ´Dt E px,sq " V 0 p r N s , r S s qe Ds ı `Be ´Dt e Dt 0 E px,sq " e Ds ‰ ď C ´1 0 e Dt 0 e ´Dt V 0 px, sq `Be ´Dt e Dt 0 E px,sq " e Ds ‰ . Moreover by (37), the Markov property, the martingale properties and an induction argument E px,sq " e Ds ‰ " e Dt 0 E px,sq " e Ds ´1 E p r Dt 0 V 0 px, sq `Be Dt 0 E px,sq Dt 0 V 0 px, sq Be Dt 0 ´1 ˙.

	Ns , r Ss q pτ ě t ´s ´t0 | s q 1 těs pC ´1 " 0 V 0 p r N s , r S s q `Bqe ´Dpt´s ´t0 q ı Ss ´1 q " e Dτ ‰ ı ď e Dt 0 E px,sq ď E px,sq ď C ´1 " e Ds ´1 pC ´1 0 V 0 p r N s ´1 , r S s ´1 q `Bq ı ı ď C ´1 0 e " e Ds ´1 ‰ 0 e Ns ´1 , r ď `Be Dt 0 ˘ ˆ1 `C´1 0 e

t ‰ 0q " ÿ ě0 P px,sq p r N t ‰ 0, t P rs , s `1qq ď ÿ ě0 P px,sq p r N s ‰ 0, t P rs , s `1qq ď ÿ ě0 P px,sq p r N s ‰ 0q θ P px,sq pt P rs , s `1qq 1´θ . s ‰ 0q " E px,sq

" 1 r N s ´1`τ ‰0 P p r N s ´1`τ , r S s ´1`τ q p r N t 0 ‰ 0q ı ď p1

Aknowledgments

The authors thank Sylvie Méléard about some discussions on tightness on Hilbert spaces.

This work was partially supported by the Chaire "Modélisation Mathématique et Biodiversité" of VEOLIA Environment, École Polytechnique, Muséum National d'Histoire Naturelle and Fondation X and by the project PIECE (Piecewise Deterministic Markov Processes) of ANR (French national research agency).


M El Gheche

J.-F Aujol

Y Berthoumieu
email: yannick.berthoumieu@u-bordeaux.fr

C.-A Deledalle

R Fablet
email: ronan.fablet@telecom-bretagne.eu

TEXTURE SYNTHESIS GUIDED BY A LOW-RESOLUTION IMAGE

Keywords: Texture synthesis, Super resolution, Wasserstein distance, Random phase texture, Proximal algorithms, Nonconvex optimization

In this paper, we aim at synthesizing a texture from a highresolution patch and a low-resolution image. To do so, we solve a nonconvex optimization problem that involves a statistical prior and a Fourier spectrum constraint. The numerical analysis shows that the proposed approach achieves better results (in terms of visual quality) than state-of-the-art methods tailored to super-resolution or texture synthesis.

INTRODUCTION

Texture synthesis consists of creating a large, coherent, and non-periodic texture image from a given sample. In this paper, we deal with the problem of texture synthesis under the assumption that a low-resolution version of the sought texture is available in addition of the given high-resolution sample. Roughly speaking, based on the down-sampling factor r, the following two cases may occur:

• If r is low, the resulting problem is similar to image reconstruction, with the addition of a high-resolution patch.

• If r is high, the resulting problem is similar to texture synthesis, with the addition of a low-resolution image.

The present paper is placed in the second context.

Related work

The problem addressed in this paper is closely related to texture synthesis. In fact, when no low-resolution is available, the two problems are identical. Existing texture synthesis algorithms can be broadly categorized into region-growing local methods and optimization-based global methods. Local approaches grow the texture one pixel (or patch) at a time, while maintaining the spatial coherence with nearby pixels by

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux -CPU (ANR-10-IDEX-03-02).

modeling the neighborhoods with Markov fields and fractal models [START_REF] Cross | Markov random field texture models[END_REF][START_REF] Kwatra | Texture optimization for example-based synthesis[END_REF]. A weakness of these methods is that the spatial coherence between pixels is enforced at a local scale. A possible approach to circumvent this limitation consists of resorting to a small patch from which to grow the texture [START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF], or by using a multiscale setting [START_REF] Han | Multiscale texture synthesis[END_REF][START_REF] Eisemann | ZIPMAPS: zoom-into-parts texture maps[END_REF]. Global methods process the entire texture as a whole, using some criteria for measuring its similarity with a small texture patch. For example, the latter can be modeled with a statistical descriptor based on histograms [START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF] and Fourier coefficients [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]. Recently, an approach has been proposed in [START_REF] Tartavel | Constrained sparse texture synthesis[END_REF][START_REF] Tartavel | Variational texture synthesis with sparsity and spectrum constraints[END_REF], which introduces a preliminary step of dictionary learning for exploiting the given patch, and (not least) the Wasserstein distance for comparing the histograms of the entire texture with an extended version of the small patch. While the Wasserstein distance is well-known in image processing and computer vision under the name "earth mover distance" [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF], it was only recently expanded to the context of texture synthesis [START_REF] Rabin | Wasserstein regularization of imaging problem[END_REF][START_REF] Tartavel | Variational texture synthesis with sparsity and spectrum constraints[END_REF][START_REF] Tartavel | Wasserstein loss for image synthesis and restoration[END_REF]. Under the assumption that the main texture characteristics are contained in their Fourier magnitude [START_REF] Julesz | Visual pattern discrimination[END_REF], many works have shown that an efficient synthesis method is achieved when the texture phase is randomized [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF][START_REF] Tartavel | Variational texture synthesis with sparsity and spectrum constraints[END_REF]. However, The resulting algorithm is limited to non-structured textures [START_REF] Desolneux | A compact representation of random phase and gaussian textures[END_REF].

Regarding the problem considered in this paper, the authors in [START_REF] Hidane | Image zoom completion[END_REF] propose to recover the high-resolution image from a pair of images: a complete low-resolution image and a highresolution but incomplete one (sample). The resulting problem is solved using an extension of the nonlocal total variation model, where a set of connections is built between the missing high-resolution pixels and a set of pixels that lies in the sample. However, to interpolate the missing data, the authors in [START_REF] Gheche | Texture reconstruction guided by the histogram of a highresolution patch[END_REF] design a new nonlocal graph that provides better connections between the missing pixels and the high-resolution pixels. Additionally, they introduce a histogram-based statistical prior modeled by a sum of Wasserstein distances between the histogram of some linear transformations of the texture. The work in [START_REF] Fablet | Joint interpolation of multi-sensor sea surface temperature fields using non-local and statistical priors[END_REF] exploits a patch-based nonlocal regularization and a Generalized Gaussian model of the texture gradients, whose parameters are learned on the high-resolution patch. Differently from [START_REF] Hidane | Image zoom completion[END_REF] and [START_REF] Gheche | Texture reconstruction guided by the histogram of a highresolution patch[END_REF], the authors in [START_REF] Fablet | Joint interpolation of multi-sensor sea surface temperature fields using non-local and statistical priors[END_REF] use a prior on the spatial covariance of the synthesized image.

In this paper, we propose a variational method that combines a histogram-based statistical prior [START_REF] Gheche | Texture reconstruction guided by the histogram of a highresolution patch[END_REF] and a constraint on the power spectral density [START_REF] Fablet | Joint interpolation of multi-sensor sea surface temperature fields using non-local and statistical priors[END_REF]. The main novelty of the proposed approach w.r.t. [START_REF] Gheche | Texture reconstruction guided by the histogram of a highresolution patch[END_REF] lies in introducing a spectral constraint instead of the nonlocal regularization, and in considering the texture synthesis problem from free noise acquisitions with a high down-sampling factor. On the other hand, the originality of our technique w.r.t. [START_REF] Fablet | Joint interpolation of multi-sensor sea surface temperature fields using non-local and statistical priors[END_REF] is twofolds: (i) the ability to consider multi-histogram priors without being constrained to a parametric model, (ii) the global optimization procedure that we carry out using a recent primal dual proximal algorithm.

The remainder of this paper is organized as follows. In Sec. 2, we present the texture synthesis as a nonconvex optimization problem. In Sec. 3, we propose to solve this problem by an iterative scheme based on proximal methods. Experimental results, showing the accuracy of the resulting images comparing with the state-of-the-art methods, are given Sec. 4. Finally, some conclusions are drawn in Sec. 5.

Notations

Let • be the standard Euclidean norm, and let Id be the identity matrix. The domain of a function f :

R N →] - ∞, +∞] is dom f = {x ∈ R N |f (x) < +∞}. Γ 0 (R N ) is the class of lower semi-continuous convex functions from R N to ] -∞, +∞] such that dom f = ∅. Let f ∈ Γ 0 (R N ). The conjugate of f is the function f * ∈ Γ 0 (R N ) defined by f * : R N →] -∞, +∞] : u → sup x∈R N x u -f (x). When f is Gâteaux-differentiable at y ∈ R N , ∂f (y) = {∇f (y)} where ∇f (y) is the gradient of f at y. A differentiable convex function has β-Lipschitz continuous gradient ∇f if (∀(x, y) ∈ R N × R N ) ∇f (x) -∇f (y) ≤ β x -y , where β ∈ ]0, +∞[. Let C be a nonempty subset of R N ,
then ι C is the indicator function of C, equal to 0 on C and +∞ otherwise.

PROPOSED APPROACH

Let x ∈ R N be the unknown signal to be recovered (which generally corresponds to an image of size N = N 1 × N 2 .), z (1) ∈ R Q the complete low-resolution image, and z (2) ∈ R M the given sample 1 such that such that

z (1) = Dx, z (2) = M x, (1)
where D ∈ R Q×N stands for spatial down-sampling by a dyadic factor in each direction yielding Q = 2 -r N , and M ∈ R M ×N is a selection operator that extracts the patch from x.

To recover x from the observations z (1) and z (2) , We propose a variational approach that leads to solving the following optimization problem:

minimize x∈R N Dx -z (1) 2 + ι C1 (x) + H(x, z (2) ) + ι C2 (x), (2)
1 A sample denotes a small patch of the synthesized image.

where (2) }. Beside the data fidelity terms w.r.t. the observations z (1) and z (2) , we use two additional pieces of information: a term H(•, z (2) ) conveying some histogram-based statistics, and a Fourier spectrum constraint C 2 . The histogram-based term is modeled through the Wasserstein distance discussed in Section 2.1, while the frequency constraint is presented in Section 2.2.

C 1 = {x ∈ R N | M x = z

Statistical prior

The Wasserstein distance between the histograms2 of two images u ∈ R N and v ∈ R N is defined as follows [START_REF] Rabin | Wasserstein regularization of imaging problem[END_REF][START_REF] Tartavel | Variational texture synthesis with sparsity and spectrum constraints[END_REF]:

W 2 2 (ν u , ν v ) = min σ∈Σ N u -v • σ 2 , (3)
where ν u and ν v are the normalized histograms of u and v, the symbol v • σ denotes a permutation of the vector v, and Σ N is the set of all the permutations of N -length vectors. For grayscale images, the optimal permutation σ * is computed as

σ * = σ v • σ -1 u , (4)
where σ v (resp. σ u ) denotes the permutation operator that arranges the pixels of v (resp. u) in ascending order. In our case, however, the two images have a different number of pixels (u ∈ R N and v ∈ R M , with M < N ). Hence, we replicate the patch v so as to obtain a larger image ṽ ∈ R N such that the normalized histogram ν ṽ is equal to ν v . Although the Wasserstein distance is nonconvex (due to the histogram transformation), its gradient is Lipschitz-continuous and takes the following form [START_REF] Tartavel | Constrained sparse texture synthesis[END_REF][START_REF] Tartavel | Wasserstein loss for image synthesis and restoration[END_REF]]

∇ u W 2 2 (ν u , ν ṽ ) = 2(u -ṽ • σ ṽ • σ -1 u ). (5)
This property allows us to employ the Wasserstein distance into the optimization algorithm presented in Section 3. In this work, inspired from [START_REF] Gheche | Texture reconstruction guided by the histogram of a highresolution patch[END_REF], we define the term H in (2) as

H(x, z (2) ) = 4 s=1 α s W 2 2 (ν Lsx , ν z(2) s ) (6)
where, for every ∀s ∈ {1, . . . , 4}, α s > 0, z

s = L s z (2) , z (2)
s is the extension of z

s (after the linear transformation), and for every s ∈ {1, . . . , 4}, L s ∈ R Ns×N is defined as follows • L 1 : the identity matrix (N 1 = N ).

• L 2 : the concatenation of the horizontal and vertical difference operators (N 2 = 2N ).

• L 3 : the concatenation of the diagonal difference operators (N 3 = 2N ).

• L 4 : the isotropic Laplacian operator (N 4 = N ).

Fig. 1. From top to bottom: "Street", "Wall". From the left to the right: low-resolution image (which was up-sampled by a factor 4 for visualization purposes), high-resolution patch (25% of the total image), true image, the proposed approach without the spectral constraint, the proposed approach with the spectral constraint, the approach of [START_REF] Dong | Image deblurring and super-resolution by adaptive sparse domain selection and adaptive regularization[END_REF] and the approach of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF].

Fourier Spectrum constraint

Textures having the same second-order statistics share a common auto-covariance and, therefore, a common Fourier magnitude. So it is meaningful to introduce a constraint enforcing texture characteristics while preserving the Fourier magnitude, through a random phase textures. Hence, inspired from [START_REF] Fablet | Joint interpolation of multi-sensor sea surface temperature fields using non-local and statistical priors[END_REF], we search for an image x ∈ R N such that

∀m, |x(m)| = |ẑ (2) (m)|, (7)
where x (resp. ẑ( 2) ) is the orthogonal discrete Fourier transform of x (resp. z (2) ), whose coefficients read

x(m) = 1 N 1 N 2 p∈R N x(p)e -2iπ( p 1 n 1 N 1 + p 2 n 2 N 2
) .

Consequently, Eq. ( 7) can be rewritten as a Fourier spectrum constraint defined as

C 2 = {x ∈ R N | ∀m, ∃ ϕ(m) : x(m) = e iϕ(m) ẑ(2) (m)} (9
) Since x and z (2) are real images, ϕ(m) of Eq. ( 9) must be antisymmetric modulo 2π. The projection onto the convex set C 2 consists in putting together the phase and the modulus, whose expression is given in Sec. 3.

Discussion

We observed that Gaussian textures could be accurately synthesized using the spectral constraint introduced in Eq. ( 9). However, for structural textures with important oscillations, the latter hypothesis is not valid. Nonetheless, taking into account the statistical information through histogram priors and the low frequencies provided by the low-resolution acquisition, the proposed approach (see Eq. ( 2)) infers the correct structures and details of the sample to the estimated texture, leading to improved results w.r.t. the classical spectral prior considered alone.

OPTIMIZATION

The solution of Eq. ( 2) requires an efficient algorithm for dealing with problems involving nonsmooth functions and linear operators. Recently, it has been shown experimentally that primal-dual proximal methods [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Briceño-Arias | A monotone + skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, lipschitzian, and parallel-sum type monotone operators[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms[END_REF][START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving largescale optimization problems[END_REF], which were originally designed for convex optimization, can be also applied to nonconvex problems in some circumstances [START_REF] Storath | Jump-sparse and sparse recovery using Potts functionals[END_REF][START_REF] Möllenhoff | Energy Minimization Methods in Computer Vision and Pattern Recognition: 10th International Conference[END_REF]. In the convex setting, the key tool of these methods is the proximity operator [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] of a lower semicontinuous convex function f : R N →] -∞; +∞], defined as

(∀y ∈ R N ) prox f (y) = argmin z∈R N f (z) + 1 2 z -y 2 .
Proximity operators enjoy many properties [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. In particular, they generalize the notion of projection onto a closed convex set C, in the sense that prox ι C = P C . Among the wide array of existing proximal algorithms, we employ the Forward-Backward Primal Dual method (FBPD) [START_REF] Condat | A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms[END_REF] reported in Algorithm 1.

Algorithm 1 FBPD [24] INITIALIZATION     Choose x [0] , y [0] ∈ R n × R Kn set τ > 0 and ω > 0 such that τ (β/2 + ω) < 1 FOR l = 0, 1, . . .          x [l] = ∇f (x [l] ) + y [l]
x

[l+1] = P {M • =z (2) } x [l] -τ x [l] y [l] = 2x [l+1] -x [l] y [l+1] = y [l] + ω y [l] -PC 2 y [l] + ω y [l]
The operators required by this algorithm are detailed below.

Fig. 2. From top to bottom: "Woolen", "Wood". From the left to the right: low-resolution image (which was up-sampled by a factor 8 for visualization purposes), high-resolution patch (25% of the total image), true image, the proposed approach without the spectral constraint, the proposed approach with the spectral constraint, the approach of [START_REF] Dong | Image deblurring and super-resolution by adaptive sparse domain selection and adaptive regularization[END_REF] and the approach of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF].

• The projection onto the set C 2 is given by the following expression

P C2 (x(m)) = x(m) • ẑ(2) (m) |x(m) • ẑ(2) (m)| ẑ(2) (m), (10)
where x • y denotes the hermitian product.

• The projection onto the convex set associated to the constraint Mx = z (2) is expressed as

P {M•=z (2) } (x) = x + M (z (2) -Mx). (11)
• The gradient of the sum of the remaining terms, that is

f (x) = Dx -z (1) 2 + 4 s=1 α s W 2 2 (ν Lsx , ν z(2) s ), reads ∇f = 2 D (Dx -z (1) ) + 2 4 s=1 α s L s (L s x -z(2) s • σ z(2) s • σ -1 Lsx ), ( 12
)
where ∇f is β-Lipschitz with β = 2(1 +

4 1 α s L s 2 ).
Although there is no theoretical guarantee about the estimate produced by Algorithm 1, in our experiments we observed that it always converges to a stable solution.

EXPERIMENTAL RESULTS

Algorithm and Initialization

We address the problem by considering a hierarchical approach, which consists of initializing our algorithm with the low-resolution image and iterating three consecutive steps: (i) up-sample the initialization by a factor 2 using a bicubic interpolation, (ii) apply a patch-based approach [START_REF] Efros | Image quilting for texture synthesis and transfer[END_REF], and (iii) run the FBPD algorithm. These three steps are repeated in loop until the size of the output signal is equal to the sought one (512x512 in our case).

Results

This section provides numerical results of our synthesis algorithm. We illustrate the effect of the Fourier constraint on the outcome of Eq. ( 2). We compare our results to the most classical synthesis method [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] and to that of super-resolution method [START_REF] Dong | Image deblurring and super-resolution by adaptive sparse domain selection and adaptive regularization[END_REF]. Fig. 1 gives two texture examples where the down-sampling factor is equal to 4, and the sample amounts to 25% of the total image. As illustrated by these experiments, we assess the performance achieved with and without the spectral constraint. Using the proposed formulation (2), we obtain the best visual results while the state-of-the-art super-resolution method [START_REF] Dong | Image deblurring and super-resolution by adaptive sparse domain selection and adaptive regularization[END_REF] tends to smooth the textures and the state-of-theart texture synthesis approaches [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] produce an image (by replicating the patch) without exploiting information provided by the low-resolution image. Hence, the results may be far from the ground truth. Fig. 2 gives a similar example where the down-sampling factor is equal to 8 and the sample amounts to 25% of the total image. The obtained results are better than the ones achieved without the spectral prior. Formulation (2) leads to a better texture reconstruction in the synthesized images compared to the state-of-the-art methods. This is confirmed by looking at the reconstructed textures and the low-resolution image. Even if the synthesized image of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] is close to the sample but it is far from the low-resolution guide.

CONCLUSIONS AND FUTURE WORK

In this paper, we have investigated the application of texture synthesis guided by a low resolution image. The proposed approach is adequate for various texture images. It is also able to exploit the potentials offered by multicore/GPU parallel architectures. However, one of the current limitations of the proposed method is that it is nonconvex, thus requiring a good initialization. In our future work, we therefore plan to consider a convex relaxation of the histogram prior in order to have a convex minimization problem [START_REF] Swoboda | Convex variational image restoration with histogram priors[END_REF].

978-1-5090-0746-2/16/$31.00 c 2016 IEEE

The histogram of a signal refers to a histogram of the pixel intensity values.