Bacterial Competition in the Presence of a Virus in a Chemostat

: We derive a mathematical model that describes the competition of two populations in a chemostat in the presence of a virus. We suppose that only one population is affected by the virus. We also suppose that the substrate is continuously added to the bioreactor. We obtain a model taking the form of an “SI” epidemic model using general increasing growth rates of bacteria on the substrate and a general increasing incidence rate for the viral infection. The stability of the steady states was carried out. The system can have multiple steady states with which we can determine the necessary and sufﬁcient conditions for both existence and local stability. We exclude the possibility of periodic orbits and we prove the uniform persistence of both species. Finally, we give some numerical simulations that validate the obtained results.


Introduction
A bioreactor ( Figure 1) is a tank (a lake) in which microorganisms multiply (yeasts, bacteria, fungi, algae, etc.) which consume substrates or feed on other organisms to develop, and which use precursors and activators to produce biomass, synthesize metabolites, or even bioconvert molecules of interest (e.g., depollution). Thanks to the bioreactor, it is possible to control the culture conditions (temperature, pH, aeration, etc.) and, therefore, to collect relatively reliable experimental data for monitoring bacterial growth and/or the chemical reaction of interest. If we consider a competition between two species for an essential substrate, a classical postulate, known as the competitive exclusion principle, suggests that, at most, one species can survive and the other species disappear. This principle has been frequently demonstrated mathematically and validated experimentally (see, for example, [1][2][3][4][5]). Several works [6][7][8][9][10][11][12][13][14] have tried to explain coexistence of bacterial competitors using several approaches. Increasingly interested in aquatic environments, researchers are discovering that the organisms colonized by viruses are much more varied than the bacterial species anticipated [15]. The shape and size of some viruses are also surprising. Finally, we are beginning to measure the impact of viral diversity on the living world. Through various mechanisms, such as the destruction of a dominant species to the benefit of rarer species [16] or the transfer of viral genes to the host, viruses (bacteriophage) maintain the biodiversity of aquatic ecosystems and facilitate genetic mixing [17]. Several works [17,18] confirm that viruses have a significant role in aquatic bacterial diversity. Therefore, the role of viruses in aquatic ecosystems cannot be neglected and should be taken into account when modeling bacterial competition in an aquatic ecosystem. S X s X i X 2 DS in D(S, X s , X i , X 2 ) Figure 1. A chemostat is a well-stirred bioreactor [19] where a limiting substrate (S in ) is continuously added to a liquid culture containing two competitors (X s , X i , X 2 ) in the presence of a virus affecting only the first competitor.
Note that a viral infection can be modeled with an epidemiological model, using either a deterministic, delayed, or stochastic approach [20][21][22][23][24]. One of the basic models for the spread of a disease was proposed in [25], dividing the population into three compartments, namely the infected population compartment (I), the susceptible population compartment (S), and the recovered population compartment (R), known as the "SIR" model. An extension of the "SIR" models is given by the "SEIR" (Susceptible, Exposed, Infected, Recovered) ones [26][27][28][29]. The "SEIR" epidemic models were extended to "SVEIR" models (Susceptible, Vaccinated, Exposed, Infected, Recovered), taking into account the proportion of immigrants who have been vaccinated [30][31][32][33]. Most of these works investigated the proposed models by giving the basic reproduction number and the local and global stability of the steady states using local linearisation and Lyapunov theory.
An important question has been asked in [34]: does the presence of a virus induce the stable coexistence of bacterial competitors in an aquatic-like system?
The response was given by proposing and analyzing a mathematical model of exploitative competition in a continuous reactor containing a virus [34]. The authors assume that only the species which have the best affinity with the substrate are affected by the virus. They proved under certain conditions that the coexistence of competitor bacteria is possible. Mestivier et al. [35] and Weitz et al. [36] proposed some mathematical models where the virus dynamics are given explicitly. It is shown that the coexistence between two competitor bacteria is possible in the presence of a virulent virus. Similarly, in [37], the authors considered a mathematical model where the virus behaviour is given explicitly and they give some conditions satisfying the coexistence of all competitors.
In this paper, we propose a generalized model of the one given in [34] by considering general increasing growth rates of bacteria on the substrate and a general increasing incidence rate for the viral infection. We introduce the model in Section 2 and we give some general results. In Section 3, we discuss the case where there is no viral infection where the competitive exclusion principle is valid. In Section 4, we reduce the system to a three-dimensional one which facilitates the mathematical analysis. We discuss the local analysis in Section 4.1, we prove that there is no periodic orbits on the faces in Section 4.2, and then we conclude on the persistence in Section 4.3 and the uniform persistence in Section 4.4. Then, we return in Section 5 to the main model where we discuss the uniform persistence. Then, we give some numerical simulations in Section 6. Finally, we summarize the main results and discuss certain implications in Section 7.

Modeling Bacterial Competition in the Presence of a Virus
Consider a bio-reactor in which the bacterial competition of two species in the presence of a virus that affects only species 1 was studied (see Figures 1 and 2). Therefore, species 1 is present in two compartments, susceptibles (X s ) and infectives or bacteriophage (X i ); however, species 2 is present in a single form (X 2 ). We know that a virus requires a host to replicate; that is why we do not explicitly model the virus dynamics and we assume that the virus spreads when an infected species takes contact with a susceptible one, which is the case of the classical epidemic models ("SI", "SIS", "SIR", "SIRS", "SEIR", "SEIRS", "SVEIR"). The limiting substrate (S) was added instantaneously to the reactor with a flow rate D and a concentration S in . The culture liquid, containing the substrate, species 1(either infected or not), and species 2, is continuously mixed and removed at the same flow rate, D. Note that the viral infection concerns only species 1. We neglected all natural mortality rates compared to the dilution rate. Figure 2. Competition diagram of the competition of the two species in the presence of a virus inside a bioreactor. Compartments S, X s , X i , and X 2 are described by circles and transition rates between compartments are described by arrows and labels.
We proposed a mathematical model describing the competition of two species for a single non-reproducing growth-limiting substrate in a continuous reactor that it is wellstirred in the presence of a virus that affects only species 1 (species 2 is not susceptible to the virus attack, Figure 2). The mathematical model takes the form of an "SI" epidemic model where the main goal is to find under what conditions the coexistence of all species is possible. This model is a generalization of the model proposed in [34] by considering generalized growth rates for all species and also a generalized incidence rate for the viral infection. The model is given by the following fourth-dimensional system of ordinary differential equations: Here, S denotes the concentration of the resource with S(0) ≥ 0, whereas X s , X i , and X 2 stand for the concentrations of susceptible species 1, infected species 1, and species 2, respectively, with initial conditions satisfying X s (0) > 0, X i (0) > 0 and X 2 (0) > 0. Note that D and S in describe the dilution rate and the substrate input concentration, respectively, and are assumed to be constant and positive. Y 1 and Y 2 denote the yield coefficients, commonly referred to as the substrate-to-species-1 (either infected or not) and substrate-tospecies-2 yields, respectively. The significance of the variables and parameters is shown in Table 1. By making the following change of variable, we obtain a more simplified model. Let . Then, the model takes the form: Let us define some operating parameters as follows: D s = µ s (s in ), D 2 = µ 2 (s in ), and D i = µ i (s in ). All the mentioned parameters are positive. Through the paper, we will consider the most important case by using the following assumption: Assumption 1. The growth rates µ s , µ i , and µ 2 are increasing, non-negative, Let us define the valuess 1 ands 2 as the solutions of µ s (s) = D and µ 2 (s) = D, respectively. Remark 1.

1.
Assumption 1 expresses that species 1 has the best affinity with the substrate and then it wins the competition in the absence of the infection. Once the infection is present, Assumption 1 expresses that the non-infected species 1 (x s ) still has the best affinity with the substrate; however, infected species 1 (x i ) has a growth rate (µ i ) smaller than both growth rates (µ s and µ 2 ) of the non-infected species 1 and species 2.

2.
Monod functions (or Holling's functions type II) are candidate functions that can express growth rates ( Figure 3): where k s , k i , and k 2 are Monod constants.μ s ,μ i , andμ 2 are positive constants. All constants can be chosen such that the functions µ s , µ i , and µ 2 satisfy Assumption 1. For example, we can takeμ s >μ 2 >μ i and k s = k 2 = k i .
The model (2) of the chemostat is a dynamical system defined on the non-negative cone, for which we recall some fundamental properties (see, for instance, [38]).
All solutions of system (2) are defined, non-negative, and bounded. Proof.

1.
The invariance of R 4 + is confirmed by the following points: By adding all equations of model (2), we deduce that:Ṁ and, therefore, we obtain: Since all compartments of the sum are non-negative, we can conclude on the boundedness of the solution.

2.
It can be deduced from the relation (3). Lemma 1.s 1 ands 2 exist and are unique and satisfy 0 <s 1 <s 2 < s in .
Assumption 2. The incidence rate µ is an increasing, non-negative, C 1 (R + ) concave function, such that µ(0) = 0. Furthermore, µ satisfies: and Remark 2. Monod (or Holding type II) function is a candidate function that can express the , where k is the Monod constant.μ is the maximum incidence rate.
The incidence rate µ satisfies the following lemma.
Proof. Let x, x 1 ∈ R + , and the function ϕ Let us define the basic reproduction number R 0 for system (2) using the next-generation operator approach proposed in [39] and deduced from the third equation (infected compartment) of system (2) and, therefore, given by: Here, µ i (s 1 ) + µ (0)(s in −s 1 ) describes the mean number of infective produced in a chemostat by introducing a single infective into a totally susceptible population inside the reactor. 1 D describes the average time that an infective individual passes inside the chemostat as an infective. For the rest of the paper, we consider the most important case where R 0 > 1.
Let us recall the classical 'chemostat' model in the absence of the virus.

Virus-Free Subsystem
Consider the following three-dimensional system which is the virus-free subsystem: This model is the same as (2) in the absence of the viral infection (x i = 0). This model predicts the competitive exclusion; that is, under Assumption 1, at most, the first species (which has the best affinity with the substrate) avoids extinction; however, the second species goes to extinction (see, for example, [19,38,40,41]). Let us define the steady-states of system (6) on the non-negative quadrant by SS 0 , SS 1 , and SS 2 with : wheres 1 <s 2 (according to Lemma 1). Therefore, we have: The equilibrium point SS 1 is globally asymptotically stable [38].
Note that by introducing a virus that affects only species 1, which has the best affinity with the nutriment, we aim to give a possibility of the coexistence of both competing species.

Reduction to Three-Dimensional System
Note that all solutions of the 4D-dynamics (2) converge toward Σ. Now, because we are interested by the asymptotic behavior of the dynamics (2), we will restrict the study to Σ. Thanks to Thieme's results [42], the asymptotic behavior of the reduced dynamics will be informative for the dynamics (2); see [19,43] for other applications. The reduced dynamics of (2) on Σ is given by: where the functions h 1 , h 2 , and h 3 are given by: Thus, for (7) the state-vector (x s , x i , x 2 ) belongs to the following subset of R 3 + : Formally, let F 000 , F 100 , F 001 , and F 111 be the four equilibrium points of dynamics (7) on Λ. F 000 reflects the extinction of all species and predators, and F 100 reflects the extinction of the infected first species and the second species while the non-infected first species is present. F 001 reflects the extinction of the first species (either infected or not) while the second species is present. Finally, F 111 reflects the coexistence of both species including the first species in its two forms, infected or not. F 000 , F 100 , F 001 , and F 111 are given by:

4.
F 111 = (x s ,x i ,x 2 ), where (x s ,x i ,x 2 ) is the solution of the three-dimensional system given by: From the third equation of system (9), and by Assumption 1, there exists a unique valuē Thus, x 2 = s in − x s − x i −s 2 and the system (9) is reduced to: From the second equation of (10), we have . Therefore, from the first equation of (10), we have: The derivative of ϕ i is given by: Furthermore, we have: Therefore, the equation ϕ i (x i ) = 0 admits a unique solutionx i ∈ (0, s in −s 2 ) and, thus, the existence and uniqueness of the equilibrium point F 111 corresponding to the coexistence of all species: The following equilibrium points are either not generic or not possible; that is why they are neglected.

1.
F 101 = (x s , 0,x 2 ), where (x s ,x 2 ) is the solution of the two-dimensional system given by: This case will be ignored since it is non-generic because we obtains 1 =s 2 (classical model of bacterial competition in a chemostat).

2.
F 011 = (0,x i ,x 2 ), where (x i ,x 2 ) is the solution of the two-dimensional system given by: This case will be ignored since we obtain µ 2 (s in − x i − x 2 ) = µ 2 (s 2 ) = µ i (s 2 ) = D, which is impossible because µ i (s in ) = D i < D.

3.
F 010 = (0, s in −s 3 , 0), wheres 3 is the solution of the equation µ i (s) = D. Again, this equilibrium is not possible since µ i (s in ) = D i < D.

4.
F 110 = (x s ,x i , 0), where (x s ,x i ) is the solution of the two-dimensional system given by: Let: Therefore, the isoclines are the graphs of two functions x s = ϕ 1 (x i ) and x s = ϕ 2 (x i ) and, then, ϕ 1 (0) = s in −s 1 and The derivatives of ϕ 1 and ϕ 2 are given by Note that ψ(0) = ϕ 2 (0) − s in +s 1 > 0, since R 0 > 1 by Assumption 3. Therefore, there is no equilibrium points of the form F 110 if D i < D < D 2 . Therefore, we will consider only the equilibrium points F 000 , F 100 , F 001 , and F 111 to be the four equilibrium points of dynamics (7) on Λ and we resume them in Proposition 3. (7) admit four equilibrium points F 000 , F 100 , F 001 , and F 111 .

Local stability
The Jacobian matrix at a point (x s , x i , x 2 ) solution of system (7) is given by: 1. The Jacobian matrix calculated at the steady-state F 000 is given by: J 000 admits three eigenvalues: Then, the steady-state F 000 is a saddle point.

2.
The Jacobian matrix calculated at the steady-state F 100 = (s in −s 1 , 0, 0) is given by: where µ s and µ 2 are expressed ats 1 . J 100 admits three eigenvalues: Then, the steady-state F 100 is a saddle point.

4.
The Jacobian matrix calculated at the steady-state F 111 = (x s ,x i ,x 2 ) is given by: . J 111 admits three eigenvalues: λ 1 , λ 2 , and λ 3 , roots of the characteristic polynomial given by: We can verify, by using Maple, that a 2 > 0, a 1 > 0, a 0 > 0, and a 2 a 1 > a 0 . Then, the steady-state F 111 is locally asymptotically stable once it exists.
According to Assumptions 1-3, we resume the local stability of equilibrium points in the following proposition. Proposition 4. F 000 , F 100 , and F 001 are saddle points; however, F 111 is stable node.

No Periodic Orbits on the Faces
We start by excluding the possibility of periodic trajectory in one of the faces of the invariant set Λ.

•
Consider a trajectory of dynamics (7) on the part of Λ where x 2 = 0: defined on Λ x s x i , given by: Note that the axes x s = 0 and x i = 0 are invariant. Let us apply the transformation η s = ln(x s ) and η i = ln(x i ) for x s , x i > 0. Then, one gets the following new system: Note that using Lemma 2, we have: By the criterion of Dulac [38], system (17) (and then, system (16)) has no periodic solution. Therefore, system (7) has no periodic solution in x s x i -face (x 2 = 0). • Consider a trajectory of dynamics (7) on the part of Λ where x s = 0: defined on Λ x i x 2 , given by: Note that the axes x i = 0 and x 2 = 0 are invariant. Let ua apply the transformation η i = ln(x i ) and η 2 = ln(x 2 ) for x i , x 2 > 0. Then, one gets the following new system: Note that: From the Dulac criterion [38], system (20) (and then, system (19)) has no periodic solution. Therefore, system (7) has no periodic solution in x i x 2 -face (x s = 0). • Consider a trajectory of dynamics (7) on the part of Λ where x i = 0: defined on Λ x s x 2 , given by: Note that the axes x s = 0 and x 2 = 0 are invariant. Let ua apply the transformation η s 1 = ln(x s ) and η 2 = ln(x 2 ) for x s , x 2 > 0. Then, one gets the following new system: Note that: By the criterion of Dulac [38], system (23) (and then, system (22)) has no periodic solution. Therefore, system (7) has no periodic solution in x s x 2 -face (x i = 0).

Persistence
In this subsection, we aim to prove the coexistence of both species 1 (either infected or not) and species 2 by proving the uniform persistence of dynamics (7). The saddle points F 000 , F 100 , and F 001 are the only boundary steady states for the dynamics (7). Then, we apply the proof used in [37,43,44] using the Butler-McGehee Lemma [38] frequently to prove the persistence of system (7). (7) is persistent.

Theorem 1. Dynamics
Proof. All the faces x s x i , x s x 2 , and x i x 2 are invariant. Furthermore, stable and unstable manifolds of the boundary equilibrium points are represented in Figure 4. Consider a solution z = (x s (t), x i (t), x 2 (t)) with an initial condition z(0) = (x s (0), x i (0), x 2 (0)), where x s (0) > 0, x i (0) > 0, and x 2 (0) > 0 are given data. Let us denote ω = ω(γ + ( z(0))) to be the omega limit set of γ + ( z(0)), where γ + ( z(0)) is the positive semi-orbit passing through z(0). We aim to prove that the omega limit set has no points on each of the three faces.

•
Assume that F 000 ∈ ω. Then, ∃ z * = F 000 inside ω ∩ W s (F 000 )\{F 000 }. The stable manifold W s (F 000 ) is one-dimensional and is restricted to the x i -axis. Therefore, the entire orbit passing through z * , which is inside ω, becomes unbounded, which contradicts the existence of z * . • Assume that F 100 ∈ ω. F 100 is a saddle point with a stable manifold, W s (F 100 ), of dimension two, restricted to the x s x 2 -plane. Therefore, {F 100 } is not the entire omega limit set ω. Using the Butler-McGehee Lemma [38], there exists a point z * = F 100 inside ω ∩ W s (F 100 )\{F 100 }. Since W s (F 100 ) lies entirely in the x s x 2 -plane, and since the entire orbit through z * is in ω, this orbit is unbounded, which contradicts the fact that F 100 is inside ω. • Assume that F 001 ∈ ω. Since F 001 is a saddle point where its stable manifold W s (F 001 ) is of dimension two and is restricted to the x i x 2 -plane, then {F 001 } is not the entire omega limit set ω. Therefore, using the Butler-McGehee Lemma [38], there exists a point z * = F 001 inside ω ∩ W s (F 001 )\{F 001 }. Since W s (F 001 ) lies entirely in the x i x 2plane, and since the entire orbit through z * is in ω, this orbit is unbounded, which contradicts the fact that F 001 is inside ω.
with at least one of the components x s (t), x i (t), and x 2 (t) is zero, and suppose that z ∈ ω. Thus, the entire orbit passing through z should be inside ω. However, since the orbit should lie entirely inside either x s x i , x i x 2 , or x s x 2 faces, it should converge to one of the boundary equilibrium points, since there is no periodic trajectory. Therefore, this boundary equilibrium point is inside ω, which contradicts the fact that all boundary equilibrium points are saddle points. Therefore, each of the components of the trajectory is greater than zero: lim inf t→∞ x s (t) > 0, lim inf t→∞ x i (t) > 0 and lim inf t→∞ x 2 (t) > 0, and then system (7) is persistent (see Section 4.3 in [44] for another example).

Uniform Persistence of System (7)
Persistence and uniform persistence [45] are equivalent in many examples of mathematical models. Recall a theory in [45] stating that if D is a dynamical system, such that R 3 + and ∂R 3 + are both invariant, then D is uniformly persistent if it satisfies the following statements.
∂D is isolated, where ∂D be the restriction of D to ∂R 3 + ; 4.
∂D is acyclic. Consider the dynamics D on the invariant attractor bounded set Σ. We can apply the theorem given in [45] if ∂Σ = Σ 1 ∪ Σ 2 and D is invariant on Σ 1 , but repelling into the interior of Σ on Σ 2 if Conditions 3 and 4 are satisfied when restricting D to Σ 1 . It is clear that condition 1 is satisfied. Condition 2 is also satisfied according to Theorem 1. Condition 3 is satisfied because all boundary equilibrium points are hyperbolic and then their union forms a covering of the omega limit sets of Σ 1 . Condition 4 is also satisfied because the boundary equilibrium points are not linked cyclically. Thus, we conclude on the uniform persistence of system (7). (7) is uniformly persistent, i.e., ∃ β > 0, such that:

Uniform Persistence of System (2)
Return to the main mathematical model (2) describing the competition of two bacteria in a chemostat in the presence of a virus that affects only the first bacteria. System (2) admits E 000 = (s in , 0, 0, 0), E 100 = (s 1 , s in −s 1 , 0, 0), E 001 = (s 2 , 0, 0, s in −s 2 ), and E 111 = (s in −x s −x i −x 2 ,x s ,x i ,x 2 ) as equilibrium points. E 000 , E 100 , and E 001 are saddle points; however, E 111 is locally asymptotically stable. We need to prove the uniform persistence of the main system (2). Let z 0 = (s(0), x s (0), x i (0), x 2 (0)) with s(0) ≥ 0, x s (0) ≥ 0, x i (0) ≥ 0 and x 2 (0) ≥ 0; then, ω( z 0 ) ∈ Σ. Furthermore, assume that ∃ ø ∈ R 4 + \Σ such that the the solution converges to ø. This is not possible since Σ is a global attractor according to proposition 1. Now, suppose that ω( z 0 ) contains a point on one of the faces where one of the variables x s , x i , or x 2 is zero; therefore, the entire trajectory passing through this point should be inside ω( z 0 ). Thus, the omega limit set ω( z 0 ) should be entirely inside Σ. (2) is uniformly persistent, i.e., ∃ η > 0, such that:

Numerical Simulations
We cofirm the theoretical findings by some numerical results using Monod functions (or Holling's functions type II) to express all growth rates and the incidence rate: µ s (s) =μ s s k s + s , µ i (s) =μ i s k i + s , µ 2 (s) =μ 2 s k 2 + s and µ(s) =μ s k + s where k s , k i , k 2 , and k are Monod constants.μ s ,μ i ,μ 2 , andμ are positive constants. We used Holling type-II functions as typical examples [46,47] since they are nonlinear and satisfied all our assumptions on growth rates and incidence rates. All constants are chosen such that the functions µ s , µ i , µ 2 , and µ satisfy Assumptions 1-3.
Consider the parameters values given in the following Table 2. We give two examples that satisfy Assumptions 1-3, which ensures the persistence of both species 1 (either infected or not) and species 2, as seen in Figures 5 and 6.   (2) converges to the equilibrium F 111 where the two species coexist (either infected or not).

Conclusions
Since we aim to prove that the competitive exclusion principle is not usually valid when two competitors grow on a single essential resource, we add, in this paper, an additional mechanism of competition by adding a virus in the chemostat that affects only the first species, and then the coexistence becomes possible. We propose a mathematical model that describes the competition of two species in a chemostat in the presence of a virus. We suppose that only one population is affected by the virus. We suppose also that the substrate is continuously added to the bioreactor. We obtain a model taking the form of an 'SI' epidemic model. The stability of the steady states was carried out. The system can have multiple steady states with which we can determine the necessary and sufficient conditions for both existence and local stability. We exclude the possibility of periodic orbits and we prove the uniform persistence of both species. Finally, we give some numerical simulations that validate the obtained results.
The main result of this work is that the presence of the virus allows the coexistence of the two bacterial species when the species cannot coexist unless the virus is present. A biological explanation of this result is that the virus affects the species which should win the competition and then it gives the opportunity to the second species to persist.