X Bardina

D Bascompte

C Rovira

S Tindel

A STOCHASTIC MODEL FOR BACTERIOPHAGE THERAPIES

Keywords: 2010 Mathematics Subject Classification. Primary 60H35; Secondary 60H10, 65C30, 92C60 Bacteriophage, competition systems, Brownian motion, large deviations

In this article, we analyze a system modeling bacteriophage treatments for infections in a noisy context. In the small noise regime, we show that after a reasonable amount of time the system is close to a sane equilibrium (which is a relevant biologic information) with high probability. Mathematically speaking, our study hinges on concentration techniques for delayed stochastic differential equations.

Introduction

In the last years Bacteriophage therapies are attracting the attention of several scientific studies. They can be a new and powerful tool to treat bacterial infections or to prevent them applying the treatment to animals such as poultry or swine. Very roughly speaking, they consist in inoculating a (benign) virus in order to kill the bacteria known to be responsible of a certain disease. This kind of treatment is known since the beginning of the 20 th century, but has been in disuse in the Western world, erased by antibiotic therapies. However, a small activity in this domain has survived in the USSR, and it is now re-emerging (at least at an experimental level). Among the reasons of this re-emersion we can find the progressive slowdown in antibiotic efficiency (antibiotic resistance). Reported recent experiments include animal diseases like hemorrhagic septicemia in cattle or atrophic rhinitis in swine, and a need for suitable mathematical models is now expressed by the community.

Let us be a little more specific about the (lytic) bacteriophage mechanism: after attachment, the virus' genetic material penetrates into the bacteria and use the host's replication mechanism to self-replicate. Once this is done, the bacteria is completely spoiled while new viruses are released, ready to attack other bacteria. It should be noticed at this point that among the advantages expected from the therapy is the fact that it focuses on one specific bacteria, while antibiotics also attack autochthonous microbiota. Roughly speaking, it is also believed that viruses are likely to adapt themselves to mutations of their host bacteria.

At a mathematical level, whenever the mobility of the different biological actors is high enough, bacteriophage systems can be modeled by a kind of predator-prey equation. Namely, set S t (resp. Q t ) for the bacteria (resp. bacteriophages) concentration at time t. Consider a truncated identity function σ : R + → R + , such that σ ∈ C ∞ , σ(x) = x whenever 0 ≤ x ≤ M and σ(x) = M + 1 for x > M + 1. Then a model for the evolution of the couple (S, Q) is as follows:

dS t = [α -kσ(Q t )] S t dt dQ t = d -mQ t -kσ(Q t )S t + k b e -µζ σ(Q t-ζ )S t-ζ dt, (1)
where α is the reproducing rate of the bacteria and k is the adsorption rate. In equation [START_REF] Calsina | Optimal latent period in a bacteriophage population model structured by infection-age[END_REF], d also stands for the quantity of bacteriophages inoculated per unit of time, m is their death rate, we denote by b the number of bacteriophages which is released after replication within the bacteria cell, ζ is the delay necessary to the reproduction of bacteriophages (called latency time) and the coefficient e -µζ represents an attenuation in the release of bacteriophages (given by the expected number of bacteria cell's deaths during the latency time, where µ is the bacteria's death rate). A given initial condition (S 0 , Q 0 ) is also specified. When modeling biological phenomena, one usually assumes models like [START_REF] Calsina | Optimal latent period in a bacteriophage population model structured by infection-age[END_REF], where σ is replaced by the identity function. We have considered here the truncation of the identity σ in order to manipulate bounded coefficients in our equations, but our parameter M can also be interpreted as a maximal infection rate of bacteria by bacteriophages. One should also be aware of the fact that the latency time ζ (which can be seen as the reproduction time of the bacteriophages within the bacteria) cannot be neglected, and is generally of the same order (about 20mn) as the experiment length (about 60mn).

According to the values of the different parameters of the system and of the initial conditions, different types of equilibriums for equation (1) might emerge. We shall focus in the sequel on the simplest of these regimes, namely when d is large enough (the exact condition is kd/m > α). This makes the mathematical analysis easier, and it corresponds to the existence of a unique stable steady state E 0 = (0, d/m) for our system (in particular bacteria have been eradicated). Notice however that we can perfectly assume the regime kd/m > α since the treatment allows to inject high quantities of viruses. One should also mention a natural generalization of our problem: Consider the action of several varieties bacteriophages, which is an option widely considered among practitioners. We have restricted our analysis here to a simplified situation for sake of readability.

It is perfectly assumable that noise will appear when collecting data from laboratory tests. Moreover, when one wishes to go from in vitro to in vivo modeling, it is commonly accepted that noisy versions of the differential systems at stake have to be considered. This program has been carried out e.g. for HIV dynamics in [START_REF] Dalal | Xuerong: A stochastic model for internal HIV dynamics[END_REF] and for bacteriophages in marine organisms in [START_REF] Carletti | Mean-square stability of a stochastic model for bacteriophage infection with time delays[END_REF]. In those references it is always assumed that the noise enters in a bilinear way, which is quite natural in this situation and ensures positivity of the solution. We shall take up this strategy here, and consider system (1) with a small random perturbation of the form

dS ε t = [α -kσ(Q ε t )] S ε t dt + εσ(S ε t ) • dW 1 t dQ ε t = d -mQ ε t -kσ(Q ε t )S ε t + k b e -µζ σ(Q ε t-ζ )S ε t-ζ dt + εσ(Q ε t ) • dW 2 t , (2)
where ε is a small positive coefficient and W = (W 1 , W 2 ) is a 2-dimensional Brownian motion defined on a complete probability space (Ω, F , P). Our aim will then be to prove that for a time τ 0 within a reasonable range, the couple Z ε τ 0 := (S ε τ 0 , Q ε τ 0 ) is not too far away from its stable equilibrium E 0 . Note that reasonable range is meant here as a time which corresponds to the order of both the latency delay and the time when the immune system of the animal can cope with the remaining bacteria.

As we shall see in the sequel, the treatment of equation ( 2) involves the introduction of some rather technical assumptions on our coefficients. For sake of readability, we have thus decided to handle first the following system without delay:

dS ε t = [α -kσ(Q ε t )] S ε t dt + εσ(S ε t ) • dW 1 t dQ ε t = [d -mQ ε t + k(b -1)σ(Q ε t )S ε t ] dt + εσ(Q ε t ) • dW 2 t , (3)
where we notice that the only difference between (2) and ( 3) is that we have set ζ = 0 in the latter.

The main advantage of equation ( 3) lies into the fact that we are able to work under the following rather simple set of assumptions: Hypothesis 1.1. We will suppose that the coefficients of equation (3) satisfy: (i) The initial condition (S 0 , Q 0 ) of the system lies into the region

R 0 := 0, mM -d k(b -1)M × [d/m, M].
(ii) The coefficient γ = kd/m -α is strictly positive and M > d/m.

We shall also use extensively the following notations:

Notation 1.2. The letters c, c 1 , c 2 , . . . will stand for universal constants, whose exact value is irrelevant. For a continuous function f , we set f ∞,I = sup x∈I |f (x)|.

Then the previous loose considerations about convergence to E 0 can be summarized in the following theorem, which is the main result of our paper for our bacteriophage system without delay: Theorem 1.3. Given positive initial conditions, equation (3) admits a unique solution which is almost surely an element of C(R + , R 2 + ). Assume furthermore Hypothesis 1.1, set η = m/2 ∧ γ and consider 3 constants 1 < κ 1 < κ 2 < κ 3 . Then there exists ρ 0 such that for any ρ ≤ ρ 0 and any interval of time of the form I = [κ 1 ln(c/ρ)/η, κ 2 ln(c/ρ)/η], we have

P ( Z ε -E 0 ∞,I ≥ 2ρ) ≤ exp - c 1 ρ 2+λ ε 2 , ( 4
)
where λ is a constant satisfying λ > κ 3 /η.

Remark 1.4. Relation (4) can be interpreted in the following manner: assume that we observe a noise with intensity ε. Then the kind of deviation we might expect from the noisy system (3) with respect to the equilibrium E 0 is of order ε ϑ with ϑ = 2η/κ 3 . This range of deviation happens at a time scale of order ln(ρ -1 )/η.

A second part of our analysis is then devoted to the more realistic delayed system, for which we obtain a result which is analogous to Theorem 1.3: Theorem 1.5. Equation (4) still holds for the delayed system (2), under some slightly more restrictive conditions on the initial condition which shall be specified at Hypothesis 2.7.

Theorem 1.5 can be seen as the main result of the current paper, and deserves some additional comments:

(1) We have produced a concentration type result instead of a large deviation principle for equation [START_REF] Dalal | Xuerong: A stochastic model for internal HIV dynamics[END_REF], because it seemed more adapted to our biological context. Indeed, in the current situation one wishes to know how far we might be from the desired equilibrium at a given fixed time, instead of producing asymptotic results as in the large deviation theory. At a technical level however, we rely on large deviation type tools, and in particular on an extensive use of exponential inequalities for martingales.

(2) Let us compare our result with [START_REF] Carletti | Mean-square stability of a stochastic model for bacteriophage infection with time delays[END_REF][START_REF] Dalal | Xuerong: A stochastic model for internal HIV dynamics[END_REF], which deal with closely related systems. The interesting article [START_REF] Carletti | Mean-square stability of a stochastic model for bacteriophage infection with time delays[END_REF] is concerned with a predator-prey system similar to ours, but it assumes that a linearization procedure around equilibrium in the highly nonlinear situation (3) can be performed. The analysis relies then heavily on this unjustified step. As far as [START_REF] Dalal | Xuerong: A stochastic model for internal HIV dynamics[END_REF] is concerned, it roughly shows that if the noise intensity of the system is high enough, then HIV epidemics can be kept under control (in terms of exponential stability). This is a valuable information, but far away from our point of view which assumes a low intensity for the noise. We should also mention the related thorough deterministic studies [START_REF] Calsina | Optimal latent period in a bacteriophage population model structured by infection-age[END_REF][START_REF] Gourley | A delay reaction-diffusion model of the spread of bacteriophage infection[END_REF][START_REF] Kuang | Delay differential equations with applications in population dynamics[END_REF] (3) Mathematically speaking, it would certainly be interesting to play with the rich picture produced by equation ( 1) and its perturbed version in terms of stable and instable equilibrium. We have not delved deeper into this direction because it did not seem directly relevant to the biological problem at stake. It should be pointed out however that the analysis of our random dynamical system (3) is non standard due to the coefficient d, which accounts for the bacteriophage inoculation. Many of our considerations below will be devoted to handle this problem.

Our article is structured as follows: Section 2 is devoted to some preliminary considerations (existence and uniqueness results for our stochastic systems, convergence to equilibrium for the corresponding deterministic equations). Then we show our concentration results at Section 3. Finally some simulations are lead at Section 4 in order to illustrate the theoretical results.

Preliminaries

In this section, we give some basic results concerning our competition system. We first establish existence and uniqueness for the solution to the perturbed system (2), starting from the simpler system (3). Then we deduce some properties for the equilibria of the deterministic counterpart of both systems (2) and (3).

2.1. Existence, uniqueness and positivity of solution. Recall that we are considering the perturbed problem (2), with a coefficient σ and some initial conditions of the following form: Hypothesis 2.1. The coefficients of our differential systems satisfy the following assumptions:

(i) The function σ : R + → R + is such that σ ∈ C ∞ , and satisfies σ(x) = x for 0 ≤ x ≤ M and σ(x) = M + 1 for x > M + 1. We also assume that 0 ≤ σ ′ (x) ≤ C for all x ∈ R + , with a constant C such that C > 1.
(ii) As far as the initial condition is concerned, we assume that it is given as continuous positive functions {S 0,τ , Q 0,τ ; -ζ ≤ τ ≤ 0}. In case of the non delayed system (1), it is simply given by two positive constants (S 0 , Q 0 ).

Due to the fact that we have assumed a bounded coefficient σ, the existence and uniqueness of the solution to our differential system is a matter of standard considerations.

Theorem 2.2 (Global existence of solution). For any positive initial condition there exists a unique solution of (2), which is defined for all t ≥ 0.

Proof. It is readily checked that the coefficients of the equation are locally Lipschitz with linear growth. The existence and uniqueness of the solution is then a direct consequence of classical results (see e.g. [6, Section 5.2] for the non delayed system and [START_REF] Mohammed | Stochastic functional differential equations[END_REF] for the delayed one).

Positivity of the solution is also an important feature, if we want the quantities S t , Q t to be biologically meaningful. Moreover, part of our analysis will rely on this property, that we label for further use: Proposition 2.3 (Positivity). If we take positive initial conditions S 0,t ≥ 0, Q 0,t ≥ 0 for all t ∈ [-ζ, 0] for the system (2), then the solution fulfills S ε t ≥ 0, Q ε t ≥ 0 for all t > 0. Proof. Let us first consider the system with σ(x) = x for all x, namely:

dS ε t = [α -kQ ε t ] S ε t dt + εS ε t • dW 1 t dQ ε t = d -mQ ε t -kQ ε t S ε t + k b e -µζ Q ε t-ζ S ε t-ζ dt + εQ ε t • dW 2 t , (5)
with initial condition (S 0,t , Q 0,t ). Assuming existence and uniqueness of the solution to (5), we shall prove that S ε t , Q ε t ≥ 0 for all t ≥ 0 almost surely. Indeed, after the change of variables x t = e -εW 1 t S ε t , y t = e -εW 2 t Q ε t , we can recast (5) into the following system of differential equations with random coefficients:

   x ′ t = α -ke εW 2 t y t x t y ′ t = de -εW 2 t -my t -ke εW 1 t x t y t + k b e -µζ-ε(W 2 t -W 2 t-ζ -W 1 t-ζ ) y t-ζ x t-ζ , (6)
with initial conditions x 0 (t) = S 0,t ≥ 0, y 0 (t) = Q 0,t ≥ 0 for all t ∈ [-ζ, 0]. Then, the positivity of x t is immediate from the representation

x t = x 0 (0) exp t 0 (α -ke εW 2 s y s )ds ≥ 0.
In order to see the positivity of y t let us observe that for y 0 (0) = 0 we have

y ′ 0 = d + k b e -µζ-ε(W 2 0 -W 2 -ζ -W 1 -ζ ) y -ζ x -ζ > 0.
Therefore, for all initial condition y 0 ≥ 0 there exists δ > 0 such that y t > 0 for all t ∈ (0, δ). Let us suppose now that y t < 0 for some t > 0, and let t 0 = inf{t > 0 | y t < 0}. Due to the continuity of the solution we have that y t 0 = 0. Then

y ′ t 0 = de -εW 2 t 0 + k b e -µζ-ε(W 2 t 0 -W 2 t 0 -ζ -W 1 t 0 -ζ ) y t 0 -ζ x t 0 -ζ > 0,
which is impossible since it would yield y t > 0 for t ∈ (t 0 , t 0 + δ) for δ small enough. This contradiction means exactly that y t ≥ 0 for all t ≥ 0. Now that we have the positivity for system (5), we can prove the positivity for (2) in the following way. Let us first handle the case of S ε t , and assume that the initial condition is such that S 0,0 ≥ M. Set then τ 0 M,S = inf{t ≥ 0 such that S ε t ≤ M/2}, and observe that τ 0 M,S is a stopping time for the natural filtration of the Brownian motion W , such that S ε has remained positive until τ 0 M,S . Furthermore, the strong Markov property for (S ε , Q ε ) entails that the process

S ε τ 0 M,S +t , Q ε τ 0
M,S +t ; t ≥ 0 also satisfies (2) on the set Ω M,S = {ω ∈ Ω; τ 0 M,S < ∞}, with an initial condition S 0,0 = M/2. With these considerations in mind, we can assume that the initial condition of our differential system satisfies S 0,0 < M.

With such an initial condition we can conclude the positivity of S ε t until the stopping time τ 0 M,S = inf{t ≥ 0 such that S ε t ≥ M} as we have done for the system (5), since up to time τ 0 M,S we have σ(S ε t ) = S ε t . Then, invoking again the strong Markov property, we can also guarantee positivity until time τ 1 M,S = inf{t ≥ τ 0 M,S such that S ε t ≤ M/2} as above. We are now in a position to obtain the positivity of S ε t until time τ 1 M,S = inf{t ≥ τ 1 M,S such that S ε t ≥ M}, once again with the same reasoning than for the system (5). The global positivity of S ε t on any interval of the form [τ k M,S , τ k+1 M,S ] for k ≥ 0 now follows by iteration of this reasoning.

It remains to show that lim k→∞ τ k M,S = ∞. This is easily obtained by combining the following two ingredients: (ii) Owing to the specific coefficients we have for equation [START_REF] Carletti | Mean-square stability of a stochastic model for bacteriophage infection with time delays[END_REF], it can be checked that for any η 2 > 0 one can find η 1 > 0 small enough such that P(τ 1 M,S > η 1 ) ≥ 1 -η 2 . Details of this assertion are omitted for sake of conciseness.

(i) The increments {τ k+1 M,S -τ k M,S ; k ≥ 0} form a i.
We let the reader check that the positivity of Q ε t can be obtained along the same lines, which ends the proof.

Remark 2.4. Using the a priori positivity properties stated above, we could have also obtained existence and uniqueness of the solution for system [START_REF] Gourley | A delay reaction-diffusion model of the spread of bacteriophage infection[END_REF]. We did not include those developments for sake of conciseness.

2.2.

Analysis of the deterministic non delayed system. This section is devoted to the analysis of the non perturbed system corresponding to [START_REF] Dalal | Xuerong: A stochastic model for internal HIV dynamics[END_REF]. Namely, we shall consider the following dynamical system:

dS t = [α -kσ(Q t )] S t dt dQ t = [d -mQ t + k(b -1)σ(Q t )S t ] dt. (7)
We will give some sufficient conditions for the existence of a unique stable equilibrium E 0 and then show exponential convergence to this equilibrium.

Let us start with the basic results we shall need about equilibria of (7).

Theorem 2.5. If either M + 1 < α k or M > α k and kd m ≥ α, system (7) has a unique (positive) steady state E 0 = (0, d m ). Moreover, the bacteria-free equilibrium E 0 is asymptotically stable for kd m > α and M > d m . Proof. To obtain the equilibria, we have to find the solutions of the following equation:

0 = (α -kσ( Q)) Ŝ 0 = d -m Q + k(b -1)σ( Q) Ŝ, (8)
where Ŝ, Q are positive constants.

Owing to the first equation we have either Ŝ = 0 or α -kσ( Q) = 0. Since Ŝ = 0 and the second equation imply Q = d m , we have that bacteria-free equilibrium E 0 exists for any value of the parameters. In the case M + 1 < α k one can observe that no other equilibrium exists (since α -kσ( Q) > 0 for any Q).

Taking M > α k , α -kσ( Q) = 0 if and only if Q = α k .
Then, using the second equation in [START_REF] Mao | Stochastic differential delay equations of population dynamics[END_REF], we have

0 = d -m α k + (b -1)α Ŝ =⇒ Ŝ = mα -kd k(b -1)α ,
which is positive only for α > kd m . So we have proved the first part of the result. For the second part, the Jacobian matrix of system ( 7) at E 0 is

A 0 := α -kσ( d m ) 0 k(b -1)σ( d m ) -m .
The eigenvalues of this matrix are easily shown to be λ 0 = α -kσ( d m ) and λ 1 = -m, which are negative for kd m > α and M > d m .

We now wish to study the rate of convergence towards the E 0 equilibrium in the stable case (i.e., when kd/m > α and M > d m ). The main result we obtain to this respect is: Theorem 2.6. Under Hypothesis 1.1 and 2.1, the solution of system [START_REF] Kuang | Delay differential equations with applications in population dynamics[END_REF] with initial condition

(S 0 , Q 0 ) ∈ 0, mM -d k(b -1)M × [d/m, M]
exponentially converges to the equilibrium E 0 :

|(S t , Q t ) -E 0 | ≤ c e -ηt , with η = γ ∧ m 2 , (9)
where we recall that γ = kd m -α > 0. Proof. In order to prove our claim, we first have to show that the region

R := [0, mM -d k(b-1)M ]× [ d m , M] ⊂ [0, M] 2 is
left invariant by equation [START_REF] Kuang | Delay differential equations with applications in population dynamics[END_REF]. Towards this aim, we can invoke the same method we will use in Proposition 2.8, and we let the reader check the details. Now, since we have Q t ≤ M for all t, we can consider σ(x) = x in equation [START_REF] Kuang | Delay differential equations with applications in population dynamics[END_REF]. We will consider a version of this system centered at E 0 by means of the change of variables S = S, Q = Q -d/m. This leads to the system

   S′ t = -γ St -k Qt St Q′ t = -m Qt + kd m (b -1) St + k(b -1) Qt St . (10)
Notice that, according to our set of assumptions concerning the initial conditions, we have S0 ≥ 0 and Q0 ≥ 0. Thus the solution to (10) will remain positive for all t > 0 (it can be deduced from R being invariant, or can be proved just like in Proposition 2.3). Now, from the first equation in [START_REF] Priouret | Remarques sur les petites perturbations de systèmes dynamiques[END_REF], we have that S′ t ≤ -γ St . This implies St ≤ S0 e -γt , proving that St exponentially converges to zero.

Owing to the second equation in [START_REF] Priouret | Remarques sur les petites perturbations de systèmes dynamiques[END_REF] and using positivity properties of the solution, we also get

Q′ t ≤ -m Qt + k(b -1) S0 e -γt d m + Qt .
Finally, the variation of constants method will lead to the stated result, following the same steps we will detail later in the proof of Theorem 2.10.

2.3.

Analysis of the deterministic delayed system. We now try to generalize the results of Section 2.2 to our deterministic delayed system (1). To this aim, we shall work under the following assumptions.

Hypothesis 2.7. We will suppose that the coefficients of equation ( 3) satisfy the following conditions, valid for any t ∈ [-ζ, 0]: (i) The initial condition (S 0,t , Q 0,t ) of the system lies into the region

R 0 := [0, M] × d m , M .
(ii) We have b e -µζ Q 0,t S 0,t > d m S 0,0 , and b e -µζ > 1. (iii) The condition S 0,t < mM -d kbe -µζ M is satisfied. A first step towards exponential stability is then the invariance of a certain region under our dynamical system: Proposition 2.8. Under Hypothesis 1.1, 2.1 and 2.7, the region

R := 0, mM -d kbe -µζ M × d m , M ⊂ [0, M] 2
is left invariant by equation (1).

Proof. We separate the analysis of S and Q in two steps.

Step 1: boundedness of S. Since S is obviously positive (along the same lines as for equation ( 6)) and owing to the fact that

S ′ t = (α -kσ(Q t )) S t we obtain that S ′ t ≤ 0 whenever Q t > α k , and
S ′ t ≥ 0 whenever Q t < α k .
Furthermore, our system starts from an initial condition Q 0,0 ≥ d m > α k . Thus S is non increasing as long as Q remains in the interval [ d m , ∞).

Let us now observe what happens in the limiting case Q 0,0 = d m : recalling that our initial conditions are denoted by S 0,t , Q 0,t for t ∈ [-ζ, 0], we have

Q ′ 0 = -k d m S 0,0 + kbe -µζ σ(Q 0,-ζ )S 0,-ζ = k be -µζ Q 0,-ζ S 0,-ζ - d m S 0,0 > 0,
where we have used the fact that be -µζ Q 0,-ζ S 0,-ζ > d m S 0,0 . According to this inequality, we obtain the existence of a strictly positive ε such that Q t > d m for all t ∈ (0, ε). We thus introduce the quantity t 0 = inf{t > 0 : Q t = d m }, and notice that we have

Q ′ t 0 = -k d m S t 0 + kbe -µζ σ(Q t 0 -ζ )S t 0 -ζ .
We can now distinguish two cases:

(1) If t 0 > ζ, since S t is non-increasing in [0, t 0 ], S t 0 -ζ ≥ S t 0 and hence

Q ′ t 0 ≥ kS t 0 be -µζ σ(Q t 0 -ζ ) - d m > 0, due to the fact that be -µζ > 1, M > d m and Q t 0 -ζ > d m . (2) If t 0 ≤ ζ, since S t 0 ≤ S 0,0 we obtain Q ′ t 0 ≥ -k d m S 0,0 + kbe -µζ σ(Q 0,t 0 -ζ )S 0,t 0 -ζ = k be -µζ Q 0,t 0 -ζ S 0,t 0 -ζ - d m S 0,0 > 0,
where we have used the fact that be -µζ Q 0,t S 0,t > d m S 0,0 for all t ∈ [-ζ, 0]. This discussion allows thus to conclude that t 0 cannot be a finite time. Indeed, we should have Q ′ t 0 > 0 and hence Q increasing in a neighborhood of t 0 , while Q should be decreasing in a neighborhood of t 0 according to its very definition. We have thus reached the following partial conclusion:

Q t ≥ d m , t → S t decreasing, S t ≥ 0.
In particular, any interval of the form [0, L] for L ≥ 0 is left invariant by t → S t .

Step 2: boundedness of Q. Our claim is now reduced to prove that for (S 0,t , Q 0,t ) ∈ R we have Q t ≤ M for all t ≥ 0.

To this aim notice that, whenever Q 0,0 = M we have

Q ′ 0 = d -mM -kMS 0,0 + kbe -µζ σ(Q 0,-ζ )S 0,-ζ ≤ d -mM + kbe -µζ MS 0,-ζ < 0,
where we recall that S 0,-ζ < mM -d kbe -µζ M according to Hypothesis 2.7. This yields the existence of ε > 0 such that Q t < M for all t ∈ (0, ε).

We now define

t 1 = inf {t > 0 : Q t = M}. It is readily checked that Q ′ t 1 = d -mM -kMS t 1 + kbe -µζ σ(Q t 1 -ζ )S t 1 -ζ = d -mM -kMS t 1 + kbe -µζ Q t 1 -ζ S t 1 -ζ ≤ d -mM + kbe -µζ MS t 1 -ζ ,
and we can distinguish again two cases:

(1) If t 1 > ζ, thanks to the fact that t → S t is non-increasing on [0, t 1 ], we have

Q ′ t 1 ≤ d -mM + kbe -µζ MS 0,0 < 0, since we have assumed that S 0,0 < mM -d kbe -µζ M . (2) If t 1 ≤ ζ then Q ′ t 1 ≤ d -mM + kbe -µζ MS 0,t 1 -ζ < 0, thanks to the fact that S 0,t < mM -d kbe -µζ M for all t ∈ [-ζ, 0].
As for the discussion of the previous step, this allows thus to conclude that t 1 cannot be a finite time, due to the contradiction Q ′ t 1 < 0 and Q t < Q t 1 for all t ∈ (0, t 1 ). We have thus shown Q t ≤ M for all t ≥ 0, which finishes the proof.

Remark 2.9. Before stating the exponential convergence to the bacteria-free equilibrium result, let us observe that Theorem 2.5 still holds true for the delayed system (1). It can be easily checked using exactly the same steps we have done for the non-delayed system.

We are now ready to state our result on exponential convergence of the delayed dynamics:

Theorem 2.10. Assume Hypothesis 1.1, 2.1, and 2.7 are satisfied, and let R be the region defined at Proposition 2.8. Then the solution of system (1) with initial condition (S 0 , Q 0 ) ∈ R exponentially converges to the equilibrium E 0 :

|(S t , Q t ) -E 0 | ≤ c e -ηt , with η = γ ∧ m 2 , (11)
where we recall that γ = kd m -α > 0. Proof. According to Proposition 2.8, we have Q(t) ≤ M for all -ζ ≤ t < ∞ under our standing assumptions. Hence one can recast equation (1) as

dS t = (α -kQ t ) S t dt dQ t = d -mQ t -kQ t S t + kbe -µζ Q t-ζ S t-ζ dt

Let us perform now the change of variables

Q = Q -d m .
This transforms the previous system into

   dS t = α -k( Qt + d m ) S t dt d Qt = d -m( Qt + d m ) -k( Qt + d m )S t + kbe -µζ ( Qt-ζ + d m )S t-ζ dt. Equivalently, our new system is:    dS t = -γS t + k Qt S t dt d Qt = -m Qt -k d m S t -k Qt S t + k d m be -µζ S t-ζ + kbe -µζ Qt-ζ S t-ζ dt.
Observe now that Proposition 2.8 asserts that Q t ≥ d m for all t ≥ 0, which means that Qt ≥ 0. With our change of variables, we have also shifted our equilibrium to the point (0, 0). We now wish to prove that S t and Qt exponentially converge to 0.

The bound on S t is easily obtained: just note that dS t ≤ -γS t dt, which yields S t ≤ S 0,0 e -γt . As far as Qt is concerned, one gets the bound

d Qt dt ≤ -m Qt + k d m be -µζ S 0,0 e -γ(t-ζ) + kbe -µζ Qt-ζ S 0,0 e -γ(t-ζ) ≤ -m Qt + kbe -µζ S 0,0 e -γ(t-ζ) d m + M - d m = -m Qt + c e -γt ,
with c = kbMS 0,0 e (γ-µ)ζ , and where we have used the fact that Q t ≤ M uniformly in t.

Invoking now the variation of constant method, it is readily checked that equation ẋt = -mx t + c e -γt with initial condition x 0 = Q0,0 can be explicitly solved as

x(t) = e -mt Q0,0 + c m -γ e (m-γ)t -1 = Q0,0 - c m -γ e -mt + c m -γ e -γt .
By comparison, this entails the inequality Qt ≤ c 1 e -ηt , where c 1 = max( Q0,0 -c m-γ , c m-γ ) and η = m ∧ γ. Our proof is now finished.

Fluctuations of the random system

Let us summarize the information we have obtained up to now in the non delayed case: we are considering the system

dS ε t = [α -kσ(Q ε t )] S ε t dt + εσ(S ε t ) • dW 1 t dQ ε t = [d -mQ ε t + k(b -1)σ(Q ε t )S ε t ] dt + εσ(Q ε t ) • dW 2 t . (12)
Under Hypothesis 1.1 and 2.1, we have shown the existence of a unique equilibrium E 0 = (0, d/m) for the deterministic system [START_REF] Kuang | Delay differential equations with applications in population dynamics[END_REF], corresponding to (12) with ε = 0. Furthermore, we have constructed a region R ∈ R 2 + such that for any initial condition (S 0 , Q 0 ) ∈ R, the solution converges exponentially to E 0 , with a rate η = γ ∧ m 2 . We now wish to obtain a concentration result for the perturbed system (12), that is give a proof of Theorem 1.3. To this aim, we shall divide our proof in several subsections. Notation 3.1. We will set Z ε t for the couple (S ε t , Q ε t ), and Z 0 t for the solution to the deterministic equation [START_REF] Kuang | Delay differential equations with applications in population dynamics[END_REF].

3.1.

Reduction of the problem. Recall that Theorem 1.3 states an exponential bound (valid for ρ small enough) of the form

P ( Z ε -E 0 ∞,I ≥ 2ρ) ≤ exp - c 1 ρ 2+λ ε 2 , (13)
on any interval of the form I = [κ 1 ln(c/ρ)/η; κ 2 ln(c/ρ)/η] and 1 < κ 1 < κ 2 < κ 3 such that λ > κ 3 /η.

A first step in this direction is to consider a generic interval of the form Î = [a, b], and write

P Z ε -E 0 ∞, Î ≥ 2ρ = P ( Z ε -E 0 ∞, Î ≥ 2ρ) ∩ ( Z 0 -E 0 ∞, Î ≥ ρ) + P ( Z ε -E 0 ∞, Î ≥ 2ρ) ∩ ( Z 0 -E 0 ∞, Î ≤ ρ) , which yields P Z ε -E 0 ∞, Î ≥ 2ρ ≤ A 1 + A 2 , with A 1 = P Z 0 -E 0 ∞, Î ≥ ρ , and A 2 = P Z ε -Z 0 ∞, Î ≥ ρ . (14)
Moreover, the term A 1 is easily handled: owing to (11), we have A 1 = 0 as soon as a = κ 1 ln(c/ρ)/η with κ 1 > 1. In order to prove (13), it is thus sufficient to check the following identity:

P Z ε -Z 0 ∞,I ≥ ρ ≤ exp - c 1 ρ 2+λ ε 2 , (15)
on any interval of the form I = [κ 1 ln(c/ρ)/η; κ 2 ln(c/ρ)/η] and 1 < κ 1 < κ 2 < κ 3 . We shall focus on this inequality in the next subsection.

3.2. Exponential concentration of the stochastic equation. We will now give a general concentration result for Z ε -Z 0 on suitable time scales as follows:

Proposition 3.2. Let Z ε be the solution to (12). Then there exists ε 0 = ε 0 (M, τ ) such that, for any ρ ≤ 1 and ε ≤ ε 0 we have

P Z ε -Z 0 ∞,[0,τ ] > ρ ≤ exp - c 2 ρ 2 e κ 2 τ ε 2 , ( 16
)
where c 2 , κ 2 are strictly positive constants which do not depend on ρ, ε, but both depend on our set of parameters α, k, σ, d, m, b, M.

Proof. For notational sake, let us abbreviate f ∞,[0,τ ] into f ∞ throughout the proof.

In order to bound Z ε -Z 0 , we first seek a bound for S ε -S 0 . To this aim we notice that for the deterministic function S 0 and thanks to relation (11), one can find a constant

κ 1 = κ 1 (α, k, σ, d, m, b) such that S 0 ∞ ≤ κ 1 . Set also J 1 t := t 0 σ(S ε s ) • dW 1 s . Then |S ε t -S 0 t | ≤ t 0 (α -kσ(Q ε s )) S ε s -α -kσ(Q 0 s ) S 0 s ds + ε J 1 t ≤ t 0 (α -kσ(Q ε s )) (S ε s -S 0 s ) ds + t 0 k σ(Q ε s ) -σ(Q 0 s ) |S 0 s |ds + ε|J 1 t | ≤ t 0 (α + kM)|S ε s -S 0 s |ds + κ 1 k t 0 |Q ε s -Q 0 s |ds + ε|J 1 t |. (17)
Analogously, setting

J 2 t := t 0 σ(Q ε s ) • dW 2 s , we obtain |Q ε t -Q 0 t | ≤ t 0 (m + k(b -1)κ 1 )|Q ε s -Q 0 s |ds + t 0 k(b -1)M|S ε s -S 0 s |ds + ε|J 2 t |. (18)
Hence, putting together (17) and (18), we get the existence of two positive constants κ 2 , κ 3 such that

|Z ε t -Z 0 t | 2 ≤ κ 2 ε 2 |J 1 t | 2 + |J 2 t | 2 + κ 3 t 0 |Z ε s -Z 0 s | 2 ds,
and by a standard application of Gronwall's lemma, we get for all t ∈ [0, τ ]:

|Z ε t -Z 0 t | 2 ≤ κ 2 ε 2 |J 1 t | 2 + |J 2 t | 2 exp(κ 3 t) ≤ κ 2 ε 2 |J 1 t | 2 + |J 2 t | 2 exp(κ 3 τ ). (19)
Let us now go back to our claim (16): thanks to inequality (19), we have

P Z ε -Z 0 ∞ > ρ = P Z ε -Z 0 2 ∞ > ρ 2 ≤ P J 1 2 ∞ + J 2 2 ∞ > ρ 2 κ 2 ε 2 exp(κ 3 τ ) ≤ T 1 + T 2 ,
with

T 1 = P J 1 ∞ > κ 4 ρ ε exp(κ 5 τ )
, and

T 2 = P J 2 ∞ > κ 4 ρ ε exp(κ 5 τ )
.

We now proceed to bound the quantity T 1 , and to this aim we first write J 1 t in terms of Itô's integrals: according to [6, Definition 3.13 p. 156],

J 1 t = t 0 σ(S ε s )dW 1 s + 1 2 σ(S ε ), W 1 t ,
where •, • stands for the bracket of two semi-martingales. Invoking equation ( 12) and ordinary rules of Stratonovich differential calculus, it is also readily checked that

σ(S ε t ) = σ(S ε 0 ) + ε t 0 σσ ′ (S ε s )dW 1 s + V t ,
where V is a process with bounded variation. We thus end up with the expression

J 1 t = M1 t + V 1 t ,
where

M1 t = t 0 σ(S ε s )dW 1 s , and V 1 t = ε 2 t 0 σσ ′ (S ε s )ds,
and decompose T 1 accordingly into T 1 ≤ T 1,1 + T 1,2 , with

T 1,1 = P M1 ∞ > κ 4 ρ ε exp(κ 3 τ ) , and
T 1,2 = P V 1 ∞ > κ 4 ρ ε exp(κ 3 τ )
.

We now bound the terms T 1,1 and T 1,2 separately. The term T 1,2 is easily bounded thanks to some deterministic arguments. Indeed, since σσ ′ (x) ≤ C(M + 1) for any x ∈ R + , we have V 1 ∞ ≤ C(M + 1)ετ , so that for any ρ ≤ 1 and ε ≤ ε 1 := (κ 4 /(C(M + 1)τ exp(κ 3 τ ))) 1/2 , we have T 1,2 = 0. As far as T 1,1 is concerned, one can apply the exponential martingale inequality (see, for instance, [START_REF] Dozzi | Stochastic processes with a multidimensional parameter[END_REF]) for stochastic integrals in order to get

T 1,1 ≤ exp - κ 4 ρ 2 M 2 exp(κ 3 τ )ε 2 .
Putting together the estimates for T 1,1 and T 1,2 , we have thus obtained

T 1 ≤ exp - κ 4 ρ 2 M 2 exp(κ 3 τ )ε 2 ,
for any ρ ≤ 1 and ε ≤ ε 1 := (κ 4 /(C(M + 1)τ exp(κ 3 τ ))) 1/2 . We let the reader check that the term T 2 can be handled along the same lines, which finishes our proof.

3.3. Deviation from equilibrium. Let us now prove inequality (13): recall that we have decomposed

P( Z ε -E 0 ∞,I ≥ 2ρ) into A 1 + A 2 defined by (14). Furthermore, A 1 = 0 when Î is of the form [a, b] with a = κ 1 ln(c/ρ)/η.
In order to complete our result, let us analyze the term A 2 in the light of inequality (16). Indeed, in order to go from ( 16) to (15), it is sufficient to choose ρ, τ, λ such that

ρ 2 exp(-κ 2 τ ) > ρ 2+λ ,
which is achieved for τ < b := λ ln(1/ρ)/κ 2 . Hence our claim is satisfied on the interval Î = [a, b]. We now have to verify that this interval is nonempty, namely that a < b. This gives a linear equation in ln(1/ρ), of the form

κ 1 η [ln(1/ρ) + ln(c)] ≤ λ κ 2 ln(1/ρ).
and the reader might easily check that the following conditions are sufficient:

(i) The linear terms satisfy κ 1 η < λ κ 2 , that is λ > κ 1 κ 2 η . (ii)
We take ρ small enough, namely ρ ≤ ρ 0 in order to compensate the term ln(c).

The proof of ( 13) is now finished.

3.4. Extension to the delayed system. Let us deal now with the delayed case: as mentioned in the introduction, we consider the system

dS ε t = [α -kσ(Q ε t )] S ε t dt + εσ(S ε t ) • dW 1 t dQ ε t = d -mQ ε t -kσ(Q ε t )S ε t + k b e -µζ σ(Q ε t-ζ )S ε t-ζ dt + εσ(Q ε t ) • dW 2 t , (20)
where for any t ∈ [-ζ, 0] and for any ε > 0, (S ε t , Q ε t ) = (S 0 t , Q 0 t ). Under Hypothesis 1.1, 2.1 and 2.7 we have shown the existence of a unique equilibrium E 0 for the deterministic system (1), corresponding to (20) with ε = 0. Following the non-delayed case, we wish to obtain a concentration result for the perturbed system (20), as is given in Theorem 1.5.

The proof of this result can be carried out almost exactly as for Theorem 1.3. Let us only point out the main difference: how to get an equivalent of inequalities (17) and (18). To this aim, we set again

J 1 t := t 0 σ(S ε s ) • dW 1 s and J 2 t := t 0 εσ(Q ε s ) • dW 2 s .
Then in the delayed case, relations (17) and (18) become

|S ε t -S 0 t | ≤ t 0 (α + kM)|S ε s -S 0 s |ds + κ 1 k t 0 |Q ε s -Q 0 s |ds + ε|J 1 t |, (21)
and

|Q ε t -Q 0 t | ≤ t 0 (m + kκ 1 )|Q ε s -Q 0 s |ds + t 0 kM|S ε s -S 0 s |ds + ε|J 2 t | (22) + t 0 kbMe -µζ |S ε s-ζ -S 0 s-ζ |ds + t 0 kbk 1 e -µζ |Q ε s-ζ -Q 0 s-ζ |ds.
Using that for any t ∈ [-ζ, 0] and for any ε > 0, (S ε t , Q ε t ) = (S 0 t , Q 0 t ) we can write the bounds

|Z ε t -Z 0 t | 2 ≤ κ 2 |J 1 t | 2 + |J 2 t | 2 + κ 3 t 0 |Z ε s -Z 0 s | 2 ds.
From this point, the proof follows exactly as for Theorem 1.3.

Numerical simulations

This final section is devoted to a presentation of some numerical simulations for the system described by equation [START_REF] Carletti | Mean-square stability of a stochastic model for bacteriophage infection with time delays[END_REF]. We have chosen the parameters (α, k, d, m, b, ζ) according to some real data observed by the Molecular Biology Group of the Department of Genetics and Microbiology at Universitat Autònoma de Barcelona. However, since the stochastic effects we are dealing with cannot be taken into account in laboratory experiments, we did not try to match existing curves like those of [START_REF] Calsina | Optimal latent period in a bacteriophage population model structured by infection-age[END_REF]. We have chosen instead to compare theoretical and noisy dynamics in order to see that the quantities S and Q are close to their equilibrium after a reasonable amount of time, in spite of randomness.

It is worth noticing at this point that the parameters we have chosen for our simulations do not meet the conditions stated at Hypothesis 2.7. Indeed, those conditions were imposed in order to obtain our theoretical large deviations type results with a reasonable amount of efforts, but might be too restrictive to fit to real data experiments. Nevertheless, our simulations turn out to be satisfactory, since we observe that the solution (S t , Q t ) converges to E 0 for small values of ε in a reasonable amount of time, regardless of the violation of Hypothesis 2.7.

Specifically, we have simulated trajectories with parameters estimated on an experiment involving Salmonella ATCC14028 bacteria and UAB_Phi78 virus. From the experiments conducted by the mentioned group we have chosen the parameters as: We have also put M = 10, µ = 0.5, and we have taken the initial conditions S 0,t = 4.8e α(t+ζ) , Q 0,t = 0 for t ∈ [ζ, 0]. The time is expressed in days and the amount of virus and bacteria are expressed in tens of millions of units. Our simulations are summarized at Figure 1, in which different paths of the processes S and Q are computed. We have used an Euler type discretization scheme for our equations, implemented with the R software. We have then plotted the deterministic case (ε = 0) plus the curves corresponding to several values of ε (namely ε = -3, 1). As mentioned before, the fluctuations of S and Q (which are obviously due to the randomness we have introduced) do not prevent them to converge to equilibrium.

  i.d sequence by a simple application of the strong Markov property.

kbk 1 e-ζ 0 kbk 1 e 0 kbk 1 e

 10101 -µζ |Q ε s-ζ -Q 0 s-ζ |ds = t-µζ |Q ε s -Q 0 s |ds ≤ t -µζ |Q ε s -Q 0 s |dsThen, putting these last bounds in (21) and (22) we get the existence of two positive constants κ 2 , κ 3 such that

  (α, k, d, m, b, ζ) = (12.1622, 27.36, 0.1, 0.1947, 61, 0.01875).

Figure 1 .

 1 Figure1. Simulation of the trajectories of S and Q with real parameters for the Salmonella ATCC14028 bacteria and UAB_Phi78 virus for the deterministic case (ε = 0), for ε = -3 (red curve) and ε = 1 (blue curve).

Acknowledgments

We would like to thank Prof. Montserrat Llagostera and the rest of the members of the Molecular Biology Group of the Department of Genetics and Microbiology at Universitat Autònoma de Barcelona for providing us the data we used in Section 4, and for her useful comments and suggestions.

S. Tindel is member of the BIGS (X. Bardina and D. Bascompte are supported by the grant MTM2009-08869 from the Ministerio de Ciencia e Innovación. C. Rovira is supported by the grant MTM2009-07203 from the Ministerio de Ciencia e Innovación.

1