Coinfection in a stochastic model for bacteriophage systems

A system modeling bacteriophage treatments with coinfections in a noisy context is analyzed. We prove that in a small noise regime, the system converges in the long term to a bacteria free equilibrium. Moreover, we compare the treatment with coinfection with the treatment without coinfection, showing how the coinfection affects the dose of bacteriophages that is needed to eliminate the bacteria and the velocity of convergence to the free bacteria equilibrium.


Introduction
The emergence of pathogenic bacteria resistant to most currently available antimicrobial agents has become a critical problem in medicine. The development of alternative antiinfection modalities has become a priority. Bacteriophage therapies are one of these alternatives. Prior to the discovery and widespread use of antibiotics, it has been suggested that bacterial infections could be treated by the administration of bacteriophages, but early clinical studies with bacteriophages were not pursued in the United States and Western Europe. Nowadays, these therapies are reemerging and attracting the attention of the scientific community.
Let us explain the (lytic) bacteriophage mechanism: after attachment, the virus' genetic material penetrates into the bacteria and uses the host's replication mechanism to self-replicate. After some time τ , the bacteria encounters death releasing some new viruses, ready to attack other bacteria.
When a bacteria has absorbed a virus particle it is called an infection of a bacteria by a bacteriophage (virus). An infected bacteria, that is, one that has already been previously absorbed and infected, can have a second or additional adsorption. Some authors have used the term superinfection as a synonym for what we have described as secondary adsorption. On the other hand, the word coinfection denotes, generally, the infection of a single cell by more than one virus. Particularly, it means the infection of a single bacteria by two or more bacteriophages of the same type. Coinfection is similar but not identical to superinfection, since coinfection does not necessarily follows multiple adsorptions of a single bacteria since later phages can be blocked from entering the bacteria.
There is a long history of mathematical modelling of phage dynamics. One of the first papers was the work of Campbell [8] where he proposed a model based on a system of delay differential equations. Deterministic models can be found, for instance, in [4], [5], [6], [10], [11], [13], [15]. The literature about stochastic models is scarce. For instance, in [14] the authors give a stochastic model allowing multiple bacteriophage adsorption to host. On the other hand, [7] was one of the first papers dealing with co/infection and superfinfection models in evolutionary epidemiology, as previous models took only first infections into account. A general discussion about how superinfections and coinfections have been modeled in evolutionary epidemiology can be found in [2].
In [3] we have considered a stochastic model with a constant injection of phages into the system. This variant corresponds to a treatment for cattle against Salmonella, which was brought to our attention by the Molecular Biology Group of the Department of Genetics and Microbiology at Universitat Autònoma de Barcelona. We modeled a bacteriophage system to a kind of predator-prey equation.
Set S(t) (resp. Q(t)) for the non-infected 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 the model is as follows: where α is the reproducing rate of the bacteria, k 1 is the adsorption rate, d also stands for the quantity of bacteriophages inoculated per unit of time, m is their death rate, b is the number of bacteriophages which are released after replication within the bacteria cell, τ is the delay necessary for 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). In fact, α = β − µ where β is the bacteria's reproduction rate. A given initial condition {S 0 (t), Q 0 (t); −τ ≤ t ≤ 0} is also specified.
Given a large enough M we show that when k 1 d/m > α, there exists a unique stable steady state E 0 = (0, d/m) (bacteria have been eradicated), while when k 1 d/m < α, the point E 0 is still an equilibrium but it becomes unstable and there exists another coexistence equilibrium. This paper only studies results regarding the bacteria-free equilibrium E 0 , since it corresponds to the main practical situation, where high doses of phages are usually introduced in cattle feed. Actually we also consider a small random perturbation of the form (2) where ε is a small positive coefficient and W = (W 1 , W 2 ) is a 2-dimensional Brownian motion and with Stratonovich differentials, denoted by • dW . We get a concentration result for the perturbed system around E 0 .
Our aim is to study the problem of coinfection in the model we have presented in [3]. Due to the ambiguity in terminology we have to specify what we understand by coinfection. After the first infection by a bacteriophage and before the death of the bacteria (we assumed that there is a period of time τ ) when the bacteriophages can infect the bacteria more times. These later adsorptions will not affect the behaviour of the bacteria but they cause the destruction of the attacker bacteriophages. Thus, we will loose some bacteriophages.
In order to describe this situation, we introduce a new process I(t), that gives the infected bacteria concentration at time t. Thus, we transform the model (1) into the next one Notice that the variation on I(t) increases with the number of infected bacteria k 1 σ(Q(t))S(t) and decreases with the rate of dead infected bacteria µI(t) and the infected bacteria that disappear after the latency time. On the other hand, in dQ(t) the term −k 2 σ(Q(t))I(t) appears that explains the bacteriophages that the system looses when they try to infect infected bacteria (coinfection). Now k 1 is the adsorption rate for noninfected bacteria and k 2 is the adsorption rate by infected ones. Moreover, when k 2 = 0 we get the system without coinfection. Thus, our model (3) is an extension of the model presented in [3]. We will show that under certain conditions a unique stable steady state E 0 = (0, 0, d/m) (free-bacteria equilibrium) exists and we obtain also a concentration result around E 0 for a perturbed system. These results are analogous to those obtained in [3]. Furthermore, we will compare both models to determine the role of the coinfection in the behaviour of the system.
Our article is structured as follows: Section 2 is devoted to giving the main results with the corresponding biological discussion. Then we show the detailed mathematical analysis in Section 3.

Main results and discussion
In this section we will study our coinfection model (3) and the fluctuations of its corresponding stochastic model. Before going on with the study of the deterministic model, let us present the set of hypothesis on the coefficient σ and on the initial condition. The hypothesis on σ will be the same than those in [3].

Hypothesis 2.1
The coefficients of our differential systems satisfy the following assumptions: (ii) As far as the initial condition is concerned, we assume that it is given as continuous positive functions {S 0 (s), Q 0 (s); −τ ≤ s ≤ 0} and a constant I 0 .

Deterministic model
We need to introduce a new assumption on the initial condition in order to get the global existence and positivity of solution.
Hypothesis 2.2 Assume that the initial condition satisfies: Clearly, if Q(t) = 0 for any t ∈ [−τ, 0] then the condition is only I 0 ≥ 0. It corresponds to the case when we begin to introduce bacteriophages at time t = 0.
Under Hypothesis 2.1 and 2.2 we have the positivity and the global existence of solution (see Proposition 3.1). Nevertheless, to obtain the existence and stability of equilibrium points we need more hypothesis: a set of hypothesis on the initial condition and a set of hypothesis on the coefficients. Hypothesis 2. 3 We will suppose that the coefficients satisfy the following conditions, valid for any t ∈ [−τ, 0]: (i) The initial condition (S 0 (t), I 0 (t), Q 0 (t)) of the system lies into the region Hypothesis 2. 4 We will suppose that the coefficients satisfy that d m < M and Then under Hypothesis 2.4 Notice again that when k 2 = 0 we get the same hypothesis as the model without coinfection. Moreover, when k 2 is increasing, we find that the constant dose d must increase, i.e., if we loose more bacteriophages by coinfection we need to introduce a bigger dose of them. On the other hand, the region where the initial condition Q 0 lives can have smaller lower boundary. It means that since the dose will be bigger, the concentration of viruses in the initial condition can be smaller.
Going back to the study of our model, under Hypothesis 2.1, 2.2, 2.4 and 2.3, we can prove the boundedness of the solution of the system (see Proposition 3.2) and we are able to study the existence of equilibrium points.
Let us study the equilibrium points. We have to find the solution to the following equations: ) .
As we have explained in the introduction we are interested in the behaviour of the free disease equilibrium E 0 = (0, 0, d m ). We get that if M > α k1 and α − k 1 d m ≤ 0, the bacteria-free equilibrium E 0 is the unique steady state and if M > d m it is asymptotically stable (see Proposition 3.3). We can also get the exponential convergence to the bacteria-free equilibrium point.
Theorem 2.5 Assume Hypothesis 2.4, 2.1, and 2.3 are satisfied, and let R be the region defined at Proposition 3.2. Then the solution of system 3) with initial condition (S 0 , I 0 , Q 0 ) ∈ R exponentially converges to the equilibrium E 0 : where γ = k1d m − α > 0. Summarizing, the free equilibrium point is, in some sense, the same point that in the model without coinfection. That is, the concentration of bacteria is 0 and the concentration of bacteriophages is d m . We also have exponential convergence but in our model with coinfection it will be slower or equal that in the model without coinfection. More precisely, if in [3] it was of order e −(γ∧m)t in the model with coinfection it will be of order e −(γ∧m∧µ)t .

Stochastic fluctuations
Let us analyze the stochastic case. We will introduce the random effects following the ideas we have used in [3]. We will assume that the noise enters in a bilineal way that ensures positivity of the solution. Thus, we consider system (3) with a small random perturbation of the form 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 ) equipped with the natural filtration (F t ) t≥0 associated to the Wiener process W . Recall that •dW (t) denotes a Stratonovich integral.
The existence follows from the fact that the coefficients of the equation are locally Lipschitz with linear growth (see Theorem 2.7 in [3]). The positivity holds using the same arguments that in Proposition 2.8 in [3]. Let us introduce some notation. For a continuous function f , we set f ∞,L = sup x∈L |f (x)|. Set Z ε = (S ε , I ε , Q ε ). Then we can state the result about convergence to E 0 as follows: Theorem 2.6 Given positive initial conditions and Hypothesis 2.4, 2.1, and 2.3 , equation (6) admits a unique solution which is almost surely an element of C(R + , R 3 + ). Set η = m ∧ γ ∧ µ and consider three constants 1 < κ 1 < κ 2 < κ 3 . Then there exists ρ 0 such that for any ρ ≤ ρ 0 and any interval of time of the form L = [κ 1 ln(c/ρ)/η, κ 2 ln(c/ρ)/η], we have where λ is a constant satisfying λ > κ 3 /η. Relation (7) means that the kind of deviation we might expect from the noisy system (6) with respect to the equilibrium E 0 is of order ε 2ϑ with ϑ = 2η/κ 3 . This range of deviation happens at a time scale of order ln(ρ −1 )/η. As in the exponential convergence for deterministic model, the convergence of the stochastic model with coinfection will be slower or equal that the convergence of the stochastic model without coinfection.

Mathematical analysis and proofs
In this section we will state the results described in Section 2 and we will prove the Propositions and Theorems presented in the same section. Since some of the proofs are similar to those given in [3], we only will give some details of the proofs with new arguments and we will refer to those in [3] in the other cases.
The fist step to analyse the model is to get the existence of a global nonnegative solution. Proof: Let us study the positivity of the solution. Clearly On the other hand, if for some t 0 it holds that Q(t 0 ) = 0 then Q ′ (t 0 ) ≥ d > 0. So, Q(t) ≥ 0 for all t. Finally, we can write Thus, if for some t 0 it holds that I(t 0 ) = 0 then under Hypothesisi 2.2 we have that I ′ (t 0 ) = 0. It yields that I(t) ≥ 0 for all t.
In order to get the existence of global solution it is enough to check that the local solutions are bounded (see for instance [9]). Since S ′ (t) ≤ αS(t), we get that for all t > 0, S(t) ≤ S(0)e αt . On the other hand, Applying a Gronwall's lemma (see [12] Lemma A.1) we obtain that Finally, notice that I ′ (t) ≤ k 1 σ(Q(t))S(t) ≤ k 1 Q(t)S(t). So, fixed T , the local solutions are bounded in [0, T ].
For simplicity, set .
Now, we can also prove the following proposition that gives us the boundedness of the solution: is left invariant by equation (3).
Proof: We organize the proof in five steps.
Step 1: While Q ≥ ν then S ∈ R 1 and is nonincreasing. Since S is positive it is clear that On the other hand, our system starts from an initial condition Thus S is non increasing and remains in R 1 as long as Q ≥ ν.
Step 2: There exists a strictly positive ε such that Q(t) > ν for all t ∈ (0, ε). Notice that here a ε 0 exists such that I(t) < mM−d be −µτ µ for all t ∈ (0, ε 0 ). So, we have where we have used Hypothesis (ii) of 2.3.
Step 3: If S(t) is nonincreasing and I(t) remains in R 2 for any t ≤ T and Q(T ) = ν then Q ′ (T ) > 0. Let us consider what happens when Q(t 0 ) = ν. We now introduce the quantity t 0 = inf{t > 0 : Q(t) = ν}, and notice that we have We can now distinguish two cases: due to the fact that be −µτ > 1, M > ν and Q(t 0 − ζ) > ν.
where we have used again Hypothesis (ii) of 2.3.
This discussion allows thus to conclude that t 0 cannot be a finite time.
Step 4: If S(t) is nonincreasing for any t ≤ T and Q(T ) = M then Q ′ (T ) < 0.
To this aim notice that, whenever Q 0 (0) = M we have where we recall that S 0 (−τ ) < mM−d k1be −µτ M according to Hypothesis 2.3. This yields the existence of ε > 0 such that Q(t) < M for all t ∈ (0, ε). We now define 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 thanks to the fact that S 0 (t) < mM−d k1be −µτ M for all t ∈ [−τ, 0]. We have thus shown Q(t) ≤ M for all t ≥ 0, which finishes the proof.
Step 5: While S remains in R 1 and Q remains in R 3 then I lives in R 2 . Notice first that We have seen before that I is always nonnegative. Assume now that there exist t 1 such that I(t 1 ) = mM−d be −µτ µ . Then and so, I(t) ≤ mM−d be −µτ µ for all t ≥ 0. Conclusion: From the previous steps we get that there exists ε > 0 such that (S(t), I(t), Q(t)) ∈ R for any t ∈ [−τ, ε). Combining all the steps it is clear that they can not leave the region.
We can state the stability result as follows: Proof: When τ = 0, using that σ(Q(t)) = Q(t), around E 0 the differential matrix is: In order to study the system with delay, we linearize it around E 0 , i.e, S(t) = 0 + s(t), I(t) = 0 + i(t) and Q(t) = q(t) + d m and we assume that the solutions are exponential, i.e. s(t) = e λt s, i(t) = e λt i and q(t) = e λt q. We get Thus the characteristic equation is and the eigenvalues will be λ 1 = α − k 1 d m < 0, λ 2 = −µ < 0 and λ 3 = −m < 0 and so E 0 is stable under the same condition that when τ = 0.
We can now prove the exponential convergence to the bacteria-free equilibrium point.
Proof of Theorem 2.5: According to Proposition 3.2, we have Q(t) ≤ M for all t. Doing now the change of variablesQ = Q − d m we get: With our change of variables, we have also shifted our equilibrium to the point (0, 0, 0). We now wish to prove that S(t), I(t) andQ(t) exponentially converge to 0.
Following the same method, we get that I(t) ≤ c 2 e −(µ∧γ)t .
Finally, we can prove the stochastic convergence.
The proof finishes using a Gronwall's lemma type and exponential inequalities for martingales (see the proof of Proposition 3.2 in [3] for the detailed methods).