BACTERIOPHAGE-RESISTANT AND BACTERIOPHAGE-SENSITIVE BACTERIA IN A CHEMOSTAT

In this paper a mathematical model of the population dynamics of a bacteriophage-sensitive and a bacteriophage-resistant bacteria in a chemostat where the resistant bacteria is an inferior competitor for nutrient is studied. The focus of the study is on persistence and extinction of bacterial strains and bacteriophage.

1. Introduction.Mathematical models of viral predation on bacteria (bacteriophagy) have been extensively studied since the pioneering work of Campbell [7] and Levin, Stewart and Chao [16].Recent work appears in both the ecological [15,28,19,6] and the biomathematical literature [1,2,3,4,20,22,17].However, a mathematically rigorous study providing sharp conditions for persistence/extinction of both the virus and the bacteria has only very recently appeared and this study, see [22], considered only the simplest scenario consisting of a single phage-sensitive bacterial population and a population of virulent phages in the setting of a chemostat.
Experimental studies of the interaction of bacteria and phages typically result in the appearance of various mutants of both bacteria and phage.For example, experimental studies by Chao, Levin, and Stewart [8] of bacteria-phage interaction in continuous culture (chemostat), using a strain of E. coli bacteria and phage T7 resulted in two distinct outcomes.In one of these, a mutant strain of E. coli, resistant to phage infection, evolved within a few hundred hours and, later, a mutant phage evolved that was capable of infecting both the original E. coli strain and the resistant strain.In other replicates of the experiment, the same mutants appeared as before, and in addition, a mutant bacterial strain evolved that was resistant to both phage strains.Pairwise competition experiments between the phage-sensitive bacterial strain and the resistant mutant strains in a virus-free chemostat, described in [8], showed that the resistant mutants were inferior competitors relative to the susceptible strain.
Bohannan and Lenski [6] note that bacterial resistance is generally due to loss or modification of the receptor molecule to which a phage binds and that often this receptor is involved in bacterial metabolism.This explains the observed tradeoff 1. Sharp criteria for persistence/extinction of the sensitive bacteria, obtained for a model that did not include a resistant mutant in [22], are unaffected by the presence of the resistant bacteria.The same is true for persistence/extinction of the phages.Therefore, the evolution of a resistant bacteria cannot cause extinction of sensitive bacteria or phages.2. A necessary condition for the resistant mutant to persist is that the growth rate advantage of the sensitive organism be, on average, balanced by the phage adsorption rate to sensitive bacteria.3. Sufficient conditions for persistence of the resistant organism in the chemostat are that, (i) when it is absent, the phages and sensitive bacteria coexist in a globally stable equilibrium, and (ii) the resistant bacteria can grow at the equilibrium nutrient level set by this equilibrium.4. If the resistant organism suffers no growth rate disadvantage relative to the sensitive bacteria, then the phages are eliminated.
Obviously, our sufficient conditions for persistence of the resistant bacteria are not optimal and they are difficult to verify.Numerical simulations suggest that these conditions are satisfied for a large set of parameter values.We show that the conditions for persistence of the resistant bacteria are met when the latent period is sufficiently short.
It is well-known that predator-prey models, like those used to model phages and bacteria, can lead to coexistence of predator and prey in a stable periodic solution [1,2,4,20] as well as a stable equilibrium.We find this to be the case in our model as well.Numerical simulations and Hopf bifurcation calculations reveal attracting periodic orbits with and without the resistant bacteria.These may oscillate strongly such that bacteria and phages reach perilously low densities in parts of the cycle and quite large levels at other times.Such strong oscillations indicate that while phages and bacteria may persist in a technical mathematical sense, they are unlikely to survive over long periods due to demographic stochasticity.
Our numerical simulations reveal interesting phenomena.We simulate the behavior of our ecosystem when all parameters except for the cost of resistance born by the resistant bacteria and the virulence of the phages, as measured by the burst size, are held constant.Then, fixing the cost of resistance, we find that the largest interval of burst size for which the resistant bacteria can persist in a stable equilibrium with phages and sensitive bacteria occurs for an intermediate value of the cost of resistance.For burst sizes beyond the upper limit of this interval, oscillations develop and their amplitude increases.Furthermore, if the cost of resistance is fixed at a low value (20%) and virulence allowed to vary over a wide range, we find the chemostat is dominated by bacteria, a majority being resistant and a minority being sensitive, whose combined density is essentially at levels typical of the phage-free system.Bacterial densities are controlled by nutrient levels which are reduced to very low levels.If the cost of resistance is fixed at a relatively large value (60%) and virulence ranges over a wide range, we observe that bacteria are controlled by a much more numerous phage population, resistant bacteria are non-existent or rare depending on virulence, resource is high and goes unused, and strong oscillations may be present.
2. The model.Our model, an elaboration on the one in [22], includes phage P , phage-sensitive bacteria S that are uninfected, infected bacteria I, phage-resistant bacteria M , and nutrient R supporting bacterial growth in a well-stirred chemostat with dilution rate D and nutrient supply concentration R 0 .The resistant organism is assumed to enjoy complete resistance to phage infection.
Yield constants γ S and γ M may be scaled out by using auxiliary variables S = γ S S, M = γ M M , I = γ S I, P = γ S P and k = k/γ S .Therefore, we will hereafter assume that: f S (R) and f M (R) are the nutrient uptake functions for microbes S and M .They are assumed to be continuously differentiable functions, vanishing at zero nutrient, with positive derivative.We assume, for simplicity, that infected bacteria do not uptake nutrient.As we do not wish to waste effort on trivial cases, we make the following hypotheses: Therefore, phage-sensitive bacteria can survive in the chemostat with nutrient feed concentration R 0 and dilution rate D in the absence of phages and resistant organisms.The resistant bacteria suffers a cost of resistance to phage infection in the form of a reduced growth rate.It is assumed that its growth rate is not so reduced as to be unable to survive in the chemostat with nutrient feed concentration R 0 and dilution rate D in the absence of the sensitive strain because then it cannot survive in the presence of sensitive bacteria, with or without the presence of phages.Indeed, if either of the inequalities in (F1) are reversed, the corresponding bacteria cannot survive in the chemostat.
A popular choice of f S and f M are Michaelis-Menten type functions given by f i (R) = viR ui+R where v i , u i > 0 for i = S, M .Bohannan and Lenski [6] note that bacterial resistance is generally due to loss or modification of the receptor molecule to which a phage binds and that often this receptor is involved in bacterial metabolism.It is plausible that u S < u M or that v M < v S or both.For mathematical simplicity, we will take f M = (1 − ε)f S for our numerical simulations, where ε ∈ (0, 1) represents the cost of resistance.
Following [22], we describe variation in the latent period by a cumulative probability distribution η(τ ).To be precise, for τ > 0, η(τ ) is the probability that an infected bacterium lyses during the time period [0, τ ] following infection.Mathematically, η(τ ) = ν([0, τ ]) where ν is a probability measure on [0, ∞) so ∞ 0 dν(s) = 1.b(τ ) is the average burst size of phage with latent period τ .Generally an infected bacterium releases more than 50 phage particles [9].The average number of new phages eventually released by an infected bacterium, B, (taking into account that it could be washed out before doing so) is given by the Laplace transform of the measure bν evaluated at D, B = ∞ 0 e −Dτ b(τ )dν(τ ). ( We also assume there exists b 0 > 1 such that b(τ ) ≤ b 0 for all τ ≥ 0. Our theoretical results are derived for a general latent period distribution η and burst-size function b(•); however, for our numerical simulations, we choose a gammadistributed latent period and a constant burst size, independent of the length of the latent period.This allows the reduction of our model to ordinary differential equations by the "linear chain trick".
As shown in [22], the initial data for I must satisfy a constraint in order for solutions to be positive.Moreover, the differential equation for I can be integrated directly: where F is the sojourn function for the latency stage [27, Sec.12.1], i.e., is the probability that an infected bacterium has not yet lysed s time units after infection.As a consequence of (4), the equation for I can be dropped from our system.However, it will occasionally be useful to include it.For example, and consequently Since (1) involves a potentially infinite distributed delay, the phase space must be chosen very carefully.See e.g., [13] and [11].Here we let Here C 0 is the space of all constant functions on (−∞, 0].And C γ is defined by e γs ϕ(s) exists}, where γ > 0 is a fixed number.
To make the integral in the differential equation of P (t) in (1) valid, we also need γ < D 2 .A norm on C γ is defined as |ϕ| γ = sup{e γs |ϕ(s)| : −∞ < s ≤ 0}.B is given the maximum norm • .All biologically reasonable solutions should be non-negative so we take the state space to be the positive cone B + in B.More precisely, the initial data is given by x = (R(0), S 0 , M (0), P 0 ) ∈ B + .
More importantly, by using the framework developed in [13] and [11], we can prove the local existence and uniqueness of solutions.With these preliminaries and Theorem 9.1 in [12, Chap.12], we can prove (1) has a compact global attractor K: Theorem 2.1.Solutions of (1) corresponding to nonnegative initial data in B + exists for all t ≥ 0 and are nonnegative.Moreover, there exists a maximal compact invariant set K ⊂ B + such that K attracts all bounded sets in B + .
3. Equilibria and local stability.Equilibria of (1) and their stability properties are studied in this section.As the infected cell density is determined by the other densities via (4), the state vector is taken as (R, S, M, P ).
Our assumptions (F1)-(F2) ensure that three equilibria without phage always exist: E 0 = (R 0 , 0, 0, 0), E S = (R S , S, 0, 0), E M = (R M , 0, M , 0, ), where f i (R i ) = D, i = S, M , S = R 0 − R S , and M = R 0 − R M .E 0 and E M are unstable; E 0 is unstable to colonization by either S or M and E M is unstable to invasion by S which can out compete M in the absence of phage by our assumptions (F1)-(F2).
As noted in [22], the stability of E S is determined by the "Phage Reproduction Number".Define the "Phage Reproduction Number at bacteria density S" by: The Phage Reproduction Number is P RN = P RN (S), the evaluation of P RN (S) at S = S.
Lemma 3.1.E S is asymptotically stable if P RN < 1 and unstable if P RN > 1.
There are two more possible equilibria of (1) which include phage: It is intuitively evident that the total bacteria population cannot exceed the level S supported in E S : M + S * ≤ S. See Lemma 7.2.Necessary and sufficient conditions for existence and positivity of equilibria are summarized in the table below.

Equilibrium
existence conditions stability E R = (R 0 , 0, 0, 0) none unstable by (F1) E S = (R S , S, 0, 0) unstable by (F2) E SP = (R * , S * , 0, P * ) P RN = BkS/(D + kS) > 1 see Theorem 3.2 The phage-susceptible bacterial density S * at E SP = (R * , S * , 0, P * ) is characterized as the unique solution of P RN The resistant bacteria M can survive only in the presence of the phage so it is natural to ask whether M can invade E SP .This is easily seen to be the case if and only if where the nutrient level R * is determined by E SP .Indeed, E SM P exists and is positive if and only if M RN > 1. Similarly as for P RN , M RN is the resistant bacteria's reproductive number in the E SP environment.Notice that the nutrient density at E SM P is the same as at E M and E SM P has the same density of sensitive bacteria as E SP .Phage density is determined by It is very difficult to give criteria for stability and instability of E SP and even more so for E SM P .However, stability properties of E SP and the existence of E SM P are linked.We say an equilibrium is linearly asymptotically stable if all characteristic roots of the characteristic equation associated to the linearized system have negative real part; it is said to be linearly unstable if a characteristic root has positive real part.The next result is a consequence of the fact that the characteristic equation associated to E SP factors with one factor representing the linearization of the system without resistant bacteria and a simple linear factor measuring potential invasibility by the resistent bacteria.
Theorem 3.2.E SP is linearly asymptotically stable for (1) if it is linearly asymptotically stable for the system without M and if M RN < 1.It is linearly unstable for (1) if it is linearly unstable for the system without M or if M RN > 1.
It is well-known that E SP can lose stability through a supercritical Hopf bifurcation for the system without M (see [4,20]) for the special case of fixed latent period duration.
Figure 1 and Figure 2 illustrate the regions of stability of E SP , represented by a pentagram, and E SM P , represented by a triangle, and the existence of periodic orbits arising from Hopf bifurcation from these equilibria as two key parameters, ε and b, are varied.The burst-size distribution is assumed constant (b(τ ) ≡ b) so that b represents a measure of phage virulence and resistent bacteria are assumed to grow at rate f M = (1 − ε)f S so ε represents the cost of phage-resistance.All other parameters are fixed.The latent period obeys a gamma distribution, allowing us to use the linear chain trick to rewrite the system as a system of ODEs.See section 5 for more details of the simulations.A solid symbol (pentagram or triangle) indicates a locally asymptotically stable equilibrium and a hollow one indicates it is unstable.Circles are periodic orbits, a solid circle represents a stable orbit and  the dashed line means it is unstable.In order to more clearly display the dynamics for larger values of ε, we use two figures; Figure 2 is a continuation of Figure 1 for larger ε and the ε-scale is expanded.
We proceed to describe the curves appearing in Figure 1 and Figure 2. The horizontal line T 1 in both figures is determined by P RN = 1; for values of b below it, phage cannot survive so neither E SP nor E SM P exists.E SP exists above T 1 and it is stable in the region immediately above T 1 .The other horizontal line, i.e., H 1 , is the line along which a periodic orbit appears with M = 0 as a result of a Hopf bifurcation from E SP as b is varied with fixed ε.It is computed numerically using Mathematica to find purely imaginary roots of the characteristic equation associated with E SP for the system with M = 0, so it is independent of ε.The increasing curve T 2 above T 1 in Figure 1 is determined analytically by M RN = 1, that is, by (1 − ε)f S (R * )/D = 1.E SM P bifurcates from E SP as b is increased with fixed ε such that (ε, b) crosses above T 2 .In Figure 1, E SP is stable below T 2 and unstable above it while E SM P becomes stable, at least near the curve.Notice that in Figure 2, H 1 meets T 2 .At this point, which we designate Q = (ε c , b c ), the steady state and Hopf bifurcation coincide.At Q, the characteristic equation associated with E SP has a pair of purely imaginary roots and a zero root.E SM P may undergo a Hopf bifurcation as well and curve H 2 in Figure 1, terminating at Q, is determined by this bifurcation.It is also computed numerically using Mathematica to find, for each fixed ε, the critical value of b such that the characteristic equation associated with E SM P has purely imaginary roots.Remarkably, H 1 and H 2 , the curves along which E SP and E SM P undergo Hopf bifurcation from their respective equilibria, and the curve T 2 , along which Unfortunately, the fold-Hopf bifurcation that might be expected at this point is degenerate.
If ε < ε c is fixed and b is increased, according to the figures, bifurcations and exchanges of stability occur as follows E S → E SP → E SM P → P OM where P OM denotes a stable periodic orbit with M > 0 bifurcating from E SM P .Here, we have ignored the unstable Hopf bifurcation from the unstable E SP .If ε c < ε is fixed and b is increased, according to Figure 2, bifurcations and exchanges of stability occur as E S → E SP → P O where now P O denotes a stable periodic orbit with M = 0 bifurcating from E SP .However, as b is increased still further, P O becomes unstable to "invasion by M " and a periodic orbit with small M -component bifurcates from P O.This bifurcation occurs as the Floquet exponent of P O, changes sign.This "lift-off" bifurcation is described in Chapter 3, sec.6 of [24] for a food chain model.The curve T 3 , along which this Floquet multiplier vanishes, also meets at Q. Finally, the vertical line V in Figure 2 represents the value of For values of ε exceeding this threshold, M could not survive alone in the chemostat.
The curves described above partition parameter space into open regions in each of which numerical simulations suggest that there is a unique attractor of positive initial data indicated by solid pentagram or solid triangle or a periodic orbit with solid circle surrounding pentagram or triangle.Notable among these regions is the very large one in Figure 1 where E SM P is stable.Above it, for large values of b, a stable periodic orbit with M > 0 exists.The boundary of this oscillatory region consists of the E SM P Hopf bifurcation curve below and the curve along which the periodic orbit merges with the periodic orbit in the M = 0 subspace on the right.As one expects, when the cost of resistance is large, the region where E SP is stable is larger and it gives way to a region above it where there is a stable periodic orbit with no resistant bacteria present.As noted earlier, the oscillatory solutions may oscillate so strongly that bacteria and phage reach perilously low densities in parts of the cycle and quite large levels at other times.Such strong oscillations indicate that while phages and bacteria may persist in a technical mathematical sense, they are unlikely to survive over long periods due to demographic stochasticity.Therefore, we may regard the large region in which E SM P is stable as the region in which coexistence of sensitive-bacteria, resistant-bacteria and bacteriophages is most stable.It is notable that this region reaches maximum height at intermediate cost of resistance.
More precisely, Theorem 4.1 (a) means that there exists ε S > 0 such that We emphasize that ε S is independent of initial data satisfying S(0) > 0. Sensitive bacteria satisfy S(t) > ε S for all sufficiently large t.Theorem 4.1 (c) means that there exists P > 0 S(0) > 0 and P (0) > 0 ⇒ lim inf t→∞ P (t) > P .
Theorem 4.1 (d) gives conditions for persistence of the resistant bacteria.Obviously, resistant bacteria persist if they do not have to compete with phage-sensitive microbes since they are immune to phage infection.More interestingly, they also persist if both phage and sensitive microbes are present provided that E SP is asymptotically stable and attracts all solutions with S(0) > 0, P (0) > 0 of system (1) with M = 0. Sufficient conditions for E SP to be asymptotically stable for the system without M are known.See [4,20].In Theorem 7.4, we prove that in the special case that the latent period is of fixed duration, then E SP attracts all solutions with S(0) > 0, P (0) > 0 for system (1) with M = 0 provided that the latent period is sufficiently small.Persistence of the resistant bacteria clearly requires the presence of the phage to counter its fitness disadvantage relative to the sensitive bacteria.The following result merely phrases this in mathematical language.Theorem 4.2.Let (R(t), S(t), M (t), P (t)) be a solution of (1) with initial condition S(0) > 0, M (0) > 0 and suppose that there exists > 0 and T > 0 such that M (t) > , t > T .Then Equation ( 9) says the the growth rate advantage of the sensitive organism relative to the resistant organism (f S − f M ) is, on average, balanced by the phage infection rate (kP ) of sensitive bacteria.
It is of interest to consider the special case that the resistant organism suffers no loss of fitness relative to the phage sensitive organism.The phage cannot survive in this case.
Theorem 4.3.Suppose that (F1) holds but (F2) is replaced by has a line segment of equilibria: Every solution with S(0) + M (0) > 0 converges to an equilibrium point on L. In particular, lim t→∞ Moreover, if P RN > 1 and P (0 5. Gamma distributed latent period and numerical results.An important special case of ( 1) is the case that the latent period obeys a Gamma distribution, namely, where g m (s, a) is the probability density function with m ∈ Z + and a > 0. The mean of this distribution is m a , which also represents the average latent period.For simplicity, we assume b(τ ) = b is a constant in this section.By following the "linear chain trick" as discussed in [20] and [14], we can transform (1) into an ODE system.To perform this procedure, we introduce some new variables I j for 1 ≤ j ≤ m: Note that these equations also give the initial condition for each I j .Any solution of (1) gives rise, via (10), to a solution of and ) Moreover, any solution of (11) which exists and is bounded on R is a solution of (1).See Prop.7.3 of [20].
We use the following parameter values from [5] for the simulations: We choose m = 5 and a = 10 as gamma distribution parameters, and consequently the mean is 0.5 (the average latent period is 0.5 hour) and the variance is 0.05.Also, We also assume f M (R) = (1 − ε)f S (R), where ε ∈ (0, 1) is not so large that (F1) fails to hold.Two bifurcation diagrams, computed using XPP-AUTO, are presented using burst size b as the bifurcation parameter.In Figure 3, we fix the cost of resistance at the relatively low value ε = 0.2 and in Figure 4 it is ε = 0.61.Four plots are shown in each figure for variables R, S, M, P .Burst size b is plotted on the horizontal axis and max/min values of equilibria or periodic solutions are plotted on the vertical axis.In these diagrams, thick lines are locally asymptotically stable equilibria, thin lines are unstable ones.Hollow circles are unstable periodic orbits and solid dots are stable ones.Starting from the left (b = 0) in Figure 3, when b is small P RN < 1 and E S is stable; it has a low nutrient level R but S is large while both M and P vanish in E S .As b increases, E SP becomes a stable equilibrium, from the R-plot and Figure 1, we observe that E SP is stable when b is approximately between 10 and 14.For larger b, E SM P becomes a positive equilibrium and it is stable.Since E SP and E SM P share the same S-component, they are not distinguishable from the S-diagram.At b ≈ 80, E SP undergoes a Hopf bifurcation, an unstable periodic orbit appears.Since this orbit lies entirely in the subspace {M ≡ 0}, we see hollow circles on the b axis in the M -plot.At b ≈ 190, there is another Hopf bifurcation at E SM P resulting in a stable periodic orbit, and E SM P becomes unstable.It should be noted that although in S and P plots the periodic orbits appear to have zero minima, they are actually bounded away by uniform persistence (the lower bound is uniform for all initial data but can be extremely small).
We infer from Figure 3 that at low cost of resistance, the combined density of bacteria is essentially at levels typical of the phage-free system, a majority being resistant and a minority being sensitive.Bacterial densities are controlled by nutrient levels which are reduced to very low levels rather than by phage predation.3 but with a 61% cost of resistance.
In Figure 4, when b is small, E S is stable.Then, as b increases, E SP becomes stable, until at b ≈ 80, it undergoes a Hopf bifurcation and a stable period orbit appears.Since both E S and E SP , as well as the periodic orbit, are in {M ≡ 0}, it is easier to study these equilibria and periodic orbit in the R-plot.At b ≈ 110, E SM P becomes a positive equilibrium (the rising curve in M -diagram).Interestingly, at b ≈ 170, the periodic solution in the M = 0 hyperplane becomes unstable and a "liftoff bifurcation" occurs as a stable periodic orbit with M > 0, but small, bifurcates from the periodic solution in the M = 0 hyperplane.This lift-off bifurcation is similar to that described for the food-chain model in Chapter 3, sec.6 of [24]; it is not a small amplitude orbit.We infer from Figure 4 that at high cost of resistance, bacteria are controlled by a much more numerous phage population, resistant bacteria are non-existent or rare depending on virulence (b), available resource goes unused, and strong predator-prey oscillations may be present.
6. Summary.The mathematical modeling of bacteriophage predation on bacteria in the setting of the chemostat goes back to the classical work of Levin, Stewart and Chao [16].Unfortunately, our ability to analyze these highly nonlinear models with time-delays for virus latency in a mathematically rigorous way is lacking.Recently, one of the authors, together with H. Thieme, has provided sharp sufficient conditions for the persistence of both bacteria and phage in [22].Essentially, bacteria persist provided they can survive in the absence of phage and the phage persist provided a basic phage reproductive number exceeds unity, implying that phage can successfully invade the phage-free equilibrium.However, numerous experimental works, e.g.[8,5,6], show that the bacteria may evolve resistance to phage infection via mutation to a phage-resistant phenotype, usually at some cost to its fitness relative to the original phage-susceptible bacteria, such that the latter can out-compete the resistant organism for limiting nutrient in the absence of the phage.Our aim in this paper is to shed light on this trade-off between resistance to infection and reduced ability to compete for nutrient by seeking sufficient conditions for the persistence of the phage-resistant bacteria.Our analytical results provide a set of sufficient conditions, namely: (1) the phage reproductive number must exceed unity, (2) phage-susceptible bacteria and phage coexist in a globally stable (for positive initial data) equilibrium in the absence of the phage-resistant bacteria, and (3) that there is sufficient nutrient for the phage-resistant bacteria to grow at this equilibrium.Conditions (1) is clearly a necessary conditions for persistence of the phage-resistant bacteria since they cannot survive without the presence of the phage.Condition (2) is not a necessary condition for persistence of the phageresistant bacteria; we assume it because it allows us to verify the hypotheses of an abstract persistence theorem.Given that condition (2) holds, then (3) is necessary for persistence of the phage-resistant bacteria.Unfortunately, we are only able to verify that (2) holds in the special case of small fixed delay, a result we defer to Theorem 7.4 in Section 7. Our numerical simulations show that persistence of the phage-resistant bacteria holds for a broad set of parameter values.
While this paper was under review, we learned of recent work of Northcott, Imran, and Wolkowicz [17] which focuses on the same issue of persistence of phagesusceptible and phage-resistent bacteria in a chemostat.These authors take a different modeling approach that makes their model more mathematically tractable.They do not model phage directly; instead, they assume that phage infection of susceptible bacteria is mediated by phage-infected bacteria leading to an SIS infection model for phage-susceptible bacteria.Furthermore, as there are no phage in their model, there is no need to model virus latency, so they can avoid modeling with delay differential equations.Northcott et.al. obtain sufficient conditions for the persistence of the phage-resistant organism (and all other constituents) which are essentially the same conditions that we find.Rather surprisingly, they find a richer bifurcation scenario for their system than we do.

7.
Proofs.This section contains all proofs of main lemmas and theorems in this paper.
The last part of this section is devoted to a special case of (1).In this case, we assume that the latent period is a fixed number τ , for each sufficiently small τ , we prove that M persist uniformly by using Theorem 4.1 (d) (ii).The main result of this subsection is summarized by Theorem 7.4.Proof of Theorem 2.1.Proofs of existence and positivity of solutions are omitted since they are similar to ones in [22].Therefore, the solutions of (1) generate a semiflow Φ : R + × B + → B + defined by Φ(t, x) = (R(t), S t , M (t), P t ) for x = (R(0), S 0 , M (0), P 0 ) ∈ B + .To show the existence of the compact attractor K for Φ, by Theorem 9.1 [12,Chap 12], it suffices to show positive orbits of bounded sets are bounded and Φ is point dissipative.
To show Φ is point dissipative, we prove that there exists a bounded set V that attracts all points in B + .Let x ∈ B + be an arbitrary point.Note that (5) also implies that there exists some T 1 > 0 such that X(t) ≤ R 0 + 1 for all t > T 1 .Therefore, Y (t) ≤ X(t) ≤ R 0 + 1 when t > T 1 .Also there exists T 2 > 0 such that e −γt x < 1 for all t > T 2 .Let T = max{T 1 , T 2 }, we have and the proof is complete.
Proof of Lemma 3.1.The linearization of ( 1) about E S is given by Setting (R, S, M, P ) = xe λt we find that λ and x must satisfy A(λ)x = 0 where A(λ) is given by and B = bν(λ + D) is the Laplace transform of bν.Because (R S , S) is asymptotically stable in the linear approximation for the subsystem with M, P = 0 and because it is easily seen that the stability analysis is reducible to the following scalar "phage invasion equation": The characteristic equation associated with ( 14) is obtained by inserting the ansatz P = e λt .The equation for λ is It has a positive real root if P RN > 1.To see this simply plot both sides of ( 15) and note that they intersect for positive λ precisely when P RN > 1 holds.On the other hand, if there is a root λ of ( 15) with λ ≥ 0 then it is easy to see that P RN ≥ 1.Indeed, if λ ≥ 0 then Therefore, λ < 0 for all roots of (15) if P RN < 1.The stability assertions now follow from standard results for delay equations: for finite delays, see [12], and for infinite delay see [18].
Proof of Theorem 3.2.We again calculate the linearization of ( 1) about E SP and set (R, S, M, P ) = xe λt , then λ and x satisfy A(λ)x = 0, where A(λ) is If E SP is asymptotically stable for the system without M , then the linearized stability of E SP with M is determined by the characteristic root f M (R * ) − D. In this case, M RN < 1 is equivalent to f M (R * ) − D < 0 and all roots have negative real parts.
If M RN > 1 or if E SP is unstable for the system without M , then there must be at least one characteristic root λ such that λ > 0 so E SP is unstable for the system with M .
On the other hand, the formal solution of S(t) is given by , the solution is attracted by E M , then S(t) → 0. Since P (t) → 0, we can find a δ > 0 and T > 0 such that f S (R(t)) − D − kP (t) > δ when t > T .Also, S(T ) > 0 because S(0) > 0, hence for all solutions with S(0) > 0. Thus S persists uniformly weakly, and by Theorem 4.13 in [21], S persists uniformly.Suppose S(0) > 0 but S ∞ < D Bk .Fix > 0. By suitably translating the solution, we may assume that S(t) < S ∞ + and P (t) < P ∞ + for all t ≥ 0. By the fluctuation argument, we can choose {t j } ∞ j=1 in R + such that t j → ∞, P (t j ) → 0, and P (t j ) → P ∞ .Then Bk , thus P ∞ = 0.By Lemma 7.1, S(t) → S. Therefore, either S ∞ ≥ D Bk or S(t) → S, and S ∞ ≥ min{S, D Bk } follows.To prove statement (b) and (c) in Theorem 4.1, we adapt the method used in [22] involving the Laplace transform.The Laplace transform of function f (t) is defined as for all λ ≥ 0. If f (t) is a non-negative function, then its Laplace transform f (λ) is also a non-negative function.
We take the Laplace transform of both sides of the P (t) equation to obtain Since P (λ) and SP (λ) both exist, the integral exists.And by the Fubini-Tonelli Theorem (Theorem 2.37 in [10]), we can interchange the order of the iterated integral.Thus where, as 2γ < D, We need the auxiliary estimate S ∞ ≤ S.
Proof of Theorem 4.1 (c).We will use Theorem 8.17 in [21] and follow the notation therein.Define the state space as X = {x ∈ B + : π S (x)(0) > 0}, where π S : B + → C γ is the projection map from X to C γ defined by π S (x) = S(•).Note that X is positively invariant for Φ.Because S persists uniformly and our semiflow has a compact attractor of bounded sets in B + , the restriction of Φ to X has a compact attractor of points in X.
We also need to show {E S } is compact, invariant, weakly ρ-repelling, isolated in X and acyclic.The proof of the first two properties is trivial.
Suppose {E S } is not a weak ρ-repeller, that is, there exists some x 0 ∈ X such that ρ(x 0 ) > 0 and Φ(t, x 0 ) → E S .Then, for arbitrary ∈ (0, S), we can assume, after a time shift of the solution Φ(t, x 0 ), that S − ≤ S(t) ≤ S + for all t ≥ 0. Now we apply these inequalities in the final inequality of (20) to obtain Since P (t) > 0 for t > 0, P (λ) is positive and finite for λ > 0 so we conclude that Letting λ → 0 we find that D + kS − kBS ≥ −k (1 + B).Since > 0 is arbitrary, this contradicts that P RN > 1 so {E S } is weakly ρ-repelling.
Now we show {E S } is isolated in X and acyclic in X 0 .First, {E S } is isolated in X 0 since in X 0 our system reduces to the ODEs ( 17) and E S is asymptotically stable for (17).This also shows that E S is acyclic in X 0 .
To show {E S } is isolated in X, we choose a neighborhood U of {E S } as the -ball centered at E S in X.By picking small enough, we can make U ∩ X 0 be an isolating neighborhood of {E S } in X 0 .
Suppose K is a compact invariant set in U other than {E S }.If K ∩ X 0 is not empty, by the forward invariance of X 0 , K ⊂ X 0 .However, since U ∩ X 0 is an isolating neighborhood of {E S } in X 0 , we have K ∩ X 0 is empty.So K ⊂ X\X 0 .But then we can use the same argument that established that E S is ρ-repelling to obtain a contradiction to P RN > 1.Thus, no such compact invariant set K ⊂ U , distinct from E S , exists and {E S } is isolated in X.
Let X = {x ∈ B + : π S (x)(0) > 0, π P (x)(0) > 0}.Because S and P persist uniformly, by parts (a) and (c), and Φ has a compact attractor of bounded sets in B + , it follows that the restriction of Φ to X has a compact attractor of points, as required for Theorem 8.17 in [21].
E SP is also weakly ρ-repelling.Suppose there exists x 0 ∈ X such that ρ(x) > 0 and Φ(t, so there exists some T > 0 and > 0 such that for all t > T , f M (R(t)) − D > .This implies that M (t) → ∞ as t → ∞, a contradiction.Hence there is no x 0 ∈ X with ρ(x) > 0 satisfying Φ(t, x 0 ) → E SP .
As E SP is asymptotically stable in X 0 , it is isolated and acyclic in X 0 .It is also isolated in X because if there is an entire trajectory in X, in a sufficiently small neighborhood of E SP , with M (0) > 0, then as R(t)−R * is small, so f M (R(t))−D > for some positive .This is a contradiction to the boundedness of M .By Theorem 8.17 and Theorem 4.13 in [21], M persists uniformly.To prove Theorem 4.3, we need the following lemma, also used in [22].
By the boundedness of M (t), we have S(t) → 0. Note By an application of fluctuation argument, we have 0 ≤ −DY ∞ , and by nonnegativity of solutions, Y ∞ = 0.This contradiction shows that P ∞ = 0 and consequently, I ∞ = 0.
7.4.Small delay case.In this section we consider a special case of (1), i.e, when the latent period is a single fixed delay τ for some τ ≥ 0. In this case (1) becomes where b > 1 is a constant.
For system (25), we will show that when τ > 0 is small, M persists uniformly.The main conclusion of this section is summarized below: Theorem 7.4.If bkS > D + kS, then there exists τ 0 > 0 such that for each τ ∈ [0, τ 0 ], we can find τ > 0 satisfying lim inf t→∞ M (t) > τ provided S(0) > 0, P (0) > 0 and M (0) > 0. This is a direct application of part (d) of Theorem 4.1, if we can show E SP attracts all orbits in {M ≡ 0} with positive initial data.To see that this sufficient condition is true, we first prove that E SP is "globally" stable in {M ≡ 0} when τ 0 = 0, and then analyze the stability of E SP for small τ > 0 by a perturbation argument.Without loss of generality, we assume τ 0 ≤ 1. System (25) can be obtained by letting the probability distribution of the latent period be a Delta distribution.All results regarding equilibria, boundedness, positivity and persistence obtained for (1) still hold for (25).Since components of E SP depend on τ , we write E τ SP to emphasize that it is an equilibrium of (25) with the given delay parameter τ .Since E 0 and E S remain the same for all τ ≥ 0, we skip the superscript to simplify notation.Nevertheless, it is worth mentioning that the phage reproduction number now becomes P RN = e −Dτ bkS D + kS .
Again, since I(t) can be solved by its formal solution, we omit the differential equation of I(t) in ( 25) hereafter.Our main focus is the stability of E SP given that M ≡ 0, so it suffices to consider the following system: Let the phase space of (26) be The first lemma shows that E SP attracts all orbits in {M ≡ 0} with positive initial data if τ = 0. Lemma 7.5.Suppose bkS > D + KS and τ = 0, then E 0 SP attracts all orbits of (26) with S(0) > 0 and P (0) > 0.
Part (a) and (c) of Theorem 4.1 claim that for each fixed τ > 0, S and P persist uniformly for all positive initial condition.The following lemma shows that the lower bound on the limit inferior can be chosen to be independent of the parameter τ , provided it is small.Proof.We use Theorem 5 in [25] to prove this lemma.
Let τ max = min{1, − 1 D ln D+kS bkS }.Ψ τ (t, x) is continuous in (τ, t, x) because of the continuous dependence of solutions of (26), both in parameters and initial data.Every positive orbit of Ψ τ has a compact closure in C. For each given τ and x ∈ C, let ω τ (x) be the omega-limit set of x.We claim that has a compact closure in C. Since ω τ (x) is invariant, ω τ (x) = Ψ τ (2, ω τ (x)).Thus we can rewrite W 1 as We have already noted that ω τ (x) ⊂ C τ .By Lemma 7.6, W 1 has a compact closure.We must verify hypothesis (B1) of Theorem 5 in [25].By Theorem 2.1, Ψ 0 has a global attractor.Moreover, for every x ∈ ∂U , Ψ 0 (t, x) either converges to E 0 or E S .Note that {{E 0 }, {E S }} is acyclic in ∂U and both {E 0 } and {E S } are isolated in C.Moreover, both stable manifolds W s (E 0 ) and W s (E S ) are disjoint from ρ −1 (0, ∞).Therefore, hypothesis (B1) is verified.Now we turn to hypothesis (B2) of Theorem 5 in [25].In the proof of part (a) of for all τ ∈ [0, τ max ] and x ∈ U .For E S , we would like to show there exists 0 ∈ (0, τ max ] such that for any τ ∈ (0, 0 ) and any x ∈ U , lim sup t→∞ Ψ τ (t, x) − E S ≥ 0 > 0. (29) If this is not true, then for any n ∈ N, there exists τ n ∈ (0, 1 n ) and x n ∈ U such that lim sup t→∞ After a suitable time shift, we may assume that S n (t) = πS (Ψ τn (t, x n ))(0) satisfies S n (t) > S − 1 n for all t ≥ 0. Let P n (t) = πP (Ψ τn (t, x n ))(0) for t ≥ 0. In the proof of part (c) of Theorem 4.1, we used the Laplace transform of P to show that E S is a repeller for phage P and we would like to use a similar argument to show that (30) is impossible.
For each n ∈ N, by (30), after a possible time-shift, we have

Figure 1 .
Figure 1.Equilibria, periodic orbits, and their stability in ε-b parameter space.Small ε.Recall that ε is the fitness cost of resistance to phage while b is the number of phage released when an infected cell lyses.Black fill indicates stable equilibrium, unfilled indicates unstable equilibrium; solid loop indicates stable periodic orbit, dashed loop indicates an unstable one.

Figure 2 .
Figure 2. Equilibria, periodic orbits and their stability in ε-b parameter space.Blow-up of large ε.

Figure 3 .
Figure 3. Bifurcation diagram with burst-size parameter b where the cost of resistance is 20%, i.e., ε = 0.2.Plots depict maximum and minimum values of R top left, S top right, M bottom left, P bottom right.

Figure 4 .
Figure 4. Bifurcation diagram with burst-size parameter b as in Figure 3 but with a 61% cost of resistance.
It occurs as the Floquet multiplier associated with the M -equation becomes unity as b is varied.For super-threshold values of b > 170, the bifurcating periodic orbit has small but positive M values while for sub-threshold b < 170 values the periodic orbit has small but negative M values.The particular region in the M -plot shows the periodic orbit with negative M values, see dot-dash lines, which is traced back to a Hopf point at E SM P , see dashed lines, which also has a negative M -component.In order to clarify this lift-off bifurcation phenomena, we provide two simulations, one at the sub-threshold b = 150 and one at super-threshold b = 200.Initial data are R(0) = 0.1, S(0) = 0.1, M (0) = 0.005, I 1 (0) = • • • = I 5 (0) = 0.001, P (0) = 1.At sub-threshold b = 150, the only positive periodic orbit lies in the M = 0 hyperplane and Figure 5 depicts M (t) converging to 0 as the solution approaches the periodic orbit in the M = 0 hyperplane.

Figure 5 .
Figure 5. Numerical simulation with 61% cost of resistance and burst size of b = 150.

Figure 6 .
Figure 6.Numerical simulation with 61% cost of resistance and burst size of b = 200.

7. 1 .
Preliminary results.Theorem 2.1 and local stability results are proved here.

7. 2 .Lemma 7 . 1 .
Proof of Theorem 4.1.Before proving the statements of Theorem 4.1, we need the following lemma: Let π P : B + → C γ be the projection map defined by π P which contradicts that S(t) → 0. Thus the solution cannot be attracted by E M .Similarly, it cannot converge to E 0 either.Therefore, E S attracts the solution.Now we proceed to the proof of Theorem 4.1.For simplicity, we prove each statement of Theorem 4.1 separately.Proof of Theorem 4.1 (a).In this proof we show S ∞ := lim sup t→∞ S(t) ≥ min{S, D Bk } Now we are ready to prove statement (b) and (c) of Theorem 4.1.

Lemma 7 . 3 . 0 P
Suppose h(t) is a non-negative and continuously differentiable function on [0, ∞) and|h (t)| is bounded, then h ∞ > 0 implies ∞ 0 h(s)ds = ∞.And now we are ready to prove Theorem 4.3.Proof of Theorem 4.3.Since f S (R) and f M (R) are identical, we write f (R) as the uptake function.Now let Y = I + 1 b0 P .We divide the proof into a few steps below.Step 1: Both I and P converge to 0. Suppose that P ∞ > 0, so ∞ 0 P (s)ds = ∞ by Lemma 7.3.Formal solutions of S(t) and M (t) imply that S(t) = S(0) M (0) M (t) exp −k t (s)ds .
s)ds , the integral is an increasing function of t, thus it either diverges to ∞ or converges to a finite limit.In both cases, M )(t) = R 0 − R, both lim t→∞ S(t) and lim t→∞ M (t) exist.