A BACTERIOPHAGE MODEL BASED ON CRISPR/CAS IMMUNE SYSTEM IN A CHEMOSTAT

. Clustered regularly interspaced short palindromic repeats (CRISPRs) along with Cas proteins are a widespread immune system across bacteria and archaea. In this paper, a mathematical model in a chemostat is proposed to investigate the eﬀect of CRISPR/Cas on the bacteriophage dynam- ics. It is shown that the introduction of CRISPR/Cas can induce a backward bifurcation and transcritical bifurcation. Numerical simulations reveal the co- existence of a stable infection-free equilibrium with an infection equilibrium, or a stable infection-free equilibrium with a stable periodic solution.


1.
Introduction. Bacteriophage or virulent phage is a virus which can grow and replicate by infecting bacteria. Once residing in bacteria, phage grow quickly, which result in the infection of bacteria and drive the bacteria to die [17]. Thus, we can view them as bacteria predators and use them to cure the diseases induced by infecting bacteria [9,35]. Phage therapy has become a promising method because of the emergence of antibiotic resistant bacteria [9,12]. Indeed, as a treatment, phages have several advantages over antibiotics. Phages replicate and grow exponentially, while antibiotics are not [47]. Generally speaking, a kind of phages infect only particular classes of bacteria, and this limitation of their host is very beneficial to cure the diseases. Moreover, phages are non-toxic, and cannot infect human cells. Hence, there are fewer side effects as compared to antibiotics [12,35].
It is important to understand the interaction dynamics between bacteria and phage to design an optimal scheme of phage therapy. There have been a number of papers that study the mathematical models of bacteriophages (see, for example, [1-4, 8-11, 21, 23, 28, 42] and the references cited therein). Campbell [11] proposed a deterministic mathematical model for bacteria and phage, which is a system of differential equations containing two state variables, susceptible bacteria S and phage P , and incorporating a delay τ that represents the time period for phage replication in the bacteria. The model was developed by [28,32,36,40] to consider the effect of predation and competition with the limited resources. It was further studied by [3,4,21] to study the infections of marine bacteriophages, where bacteria populations from 5-6% up to 70% are infected by bacteriophages [38]. Smith and Thieme [42] formulate two models: a system of delay differential equations with infinite delay and a more general infection-age model that leads to a system of integro-differential equations to study the persistence of bacteria and phage.
It was observed that phage can exert pressure on bacteria to make them produce resistance through loss or modification of the receptor molecule to which a phage binds with an inferior competition ability for nutrient uptake [25]. More recently, biological evidences are found that there exists an adaptive immune system across bacteria, which is the Clustered regularly interspaced short palindromic repeats (CRISPRs) along with Cas proteins [6,13,14,18,26,27,31,41]. In this system, phage infection is memorized via a short invader sequence, called a proto-spacer, and added into the CRISPR locus of the host genome. And, the CRISPR/Cas system admits heritable immunity [13,14,18,26,31]. The replication of infecting phage in bacteria is aborted if their DNA matches the crRNAs (CRISPR RNAs) which contains these proto-spacer. On the other hand, if there is no perfect pairing between the protospacerand the foreign DNA (as in the case of a phage mutant), the CRISPR/Cas system is counteracted and replication of the phage DNA can occur [19,33,34,39]. Therefore, the CRISPR/Cas system participates in a constant evolutionary battle between phage and bacteria [7,13,27,47].
Mathematical models are powerful in understanding the population dynamics of bacteria and phages. Han and Smith [23] formulated a mathematical model that includes a phage-resistant bacteria, where the resistant bacteria is an inferior competitor for nutrient. Their analytical results provide a set of sufficient conditions for the phage-resistant bacteria to persist. Recently, mathematical models have been proposed to study the contributions of adaptive immune response from CRISPR/Cas in bacteria and phage coevolution [27,29]. In these papers, numerical simulations are used to find how the immune response affects the coexistence of sensitive strain and resistance strain of bacteria. In the present paper, we extend the model in [23] by incorporating the CRISPR/Cas immunity on phage dynamics. Following [23], we focus on five state variables: R, S, M , I and P which represent the concentrations of nutrient, sensitive bacteria and resistant bacteria, productively infected bacteria and phage in a chemostat. The following assumptions will be used in the formulation of mathematical model: (i) all sensitive bacteria have the potential to get the acquired immunity; (ii) only the resistance from the CRISPR/Cas immunity is considered, which is incomplete because phage can mutate under the pressure of CRISPR/Cas immunity; the resistant bacteria suffer the cost of competitive ability [28]; (iii) the infected bacteria don't uptake nutrient, and cannot be infected again [23]; (iv) the latent time for phage infection is omitted because it is small; (v) the loss of phage due to multiple phage attachment to bacteria is omitted.
Let k be the valid phage adsorption rate to both susceptible and resistant bacteria, which is reasonable because the CRISPR/Cas immunity does not interfere with phage adsorption. Following [23], we let the loss rates of phages due to adsorption to susceptible bacteria and resistant bacteria be described by kSP and kM P respectively, and the infection rate of sensitive bacteria be given by kSP . Note that the CRISPR/Cas immunity arises during the process of phage invasion. If ε is the probability that the CRISPR/Cas immunity is successfully established in bacteria, then the transition rate of invaded bacteria to the resistance class is εkSP and the transition rate to the productively infected class is (1 − ε)kSP . Moreover, we let k = κk denote the infection coefficient of resistant bacteria where 0 < κ ≤ 1 describes the mutant probability of phage to escape CRISPR/Cas suppression. Finally, by incorporating the terms discussed above to describe the effects of CRISPR/Cas immunity into the model in [23], we arrive at the following model: where a dot denotes the differentiation with respect to time t; D and R 0 represent the dilution rate and nutrient supply concentration in the chemostat, respectively; f (R) is the growth rate function which is taken to be Monod type: f (R) = mR/(a + R), where m is the maximal growth rate and a is the halfsaturation constant; parameter µ with 0 < µ < 1 is the cost of bacteria's resistance; parameter b represents the total number of free virus particles released by each productively infected cell over its lifespan; parameter δ is the death rate of infected bacteria.
The paper is organized as follows. In the next section we present the mathematical analysis of the model that include the stability and bifurcation of equilibria. Numerical simulations are provided in Section 3 to show the complicated dynamical behaviors of the model. A brief discussion concludes in Section 4. The proofs of a few theorems are given in Appendix.
2. Mathematical analysis. In this section, we present the mathematical analysis for the stability and bifurcations of equilibria of (1). We start with the positivity and boundedness of solutions. Proof. We examine only the last conclusion of the proposition. First, we claim that R(t) is positive for all t > 0 in the existence interval of the solution. Suppose not. Then there is a time t 1 > 0 such that R(t 1 ) = 0 and R(t) > 0 for t ∈ (0, t 1 ). Since R(t 1 ) = 0, we haveṘ(t 1 ) = DR 0 > 0. Thus, there is a sufficiently small ε > 0 such that R(t) < 0 for all t ∈ (t 1 − ε, t 1 ), which is a contradiction. Hence, R(t) is positive for t > 0. It is evident that In a similar way we can show the positivity of M (t), I(t) and P (t).

MENGSHI SHU, RUI FU AND WENDI WANG
Calculating the derivative of L along the solution of (1), we obtaiṅ It follows that the nonnegative solutions of (1) exist on [0, ∞) and Therefore, the nonnegative solutions of model (1) are ultimately bounded. 2.2. Infection-free equilibria. E 0 = (R 0 , 0, 0, 0, 0) is a trivial equilibrium. To find others, we assume Then It is easy to see that (3) and (4) imply Thus, the competitive exclusion in the absence of phage infection holds [24], and the boundary equilibrium E 2 is unstable. In what follows, we always assume that (3) and (4) are available. The basic reproduction number R 0 of phage in the population of sensitive bacteria is computed through the matrix F of new infection and the infective transition matrix V : and is defined as the spectral radius of F V −1 [16]. Hence, Analogously, the basic reproduction number R M 0 of phage in the population of resistant bacteria is .
The proof of Theorem 2.1 is postponed to Appendix.

Proof. Define a Lyapunov function by
where R 1 = λ 1 and S 1 = R 0 −λ 1 . Calculating the derivative of V along the solution of (1), we obtaiṅ It is easy to examine that the largest invariant set in D 0 is It follow from the LaSalle's invariance principle [22] that E 1 is globally stable.

Infection equilibria.
In this subsection, we consider the infection equilibria of system (1) which satisfy It follows that and R 3 satisfies where .
By direct calculations, we obtain where R 3 lies between λ 2 and R 0 , and E 3 does not exist when R M 0 < 1. As a result, we can state the following theorem for the existence of E 3 .
To study the local stability of infection equilibria E 3 , we set Moreover, for κ < µ, we define and is unstable when κ > µ, or The proof of Theorem 2.4 is given in Appendix.
Let us now consider the existence of coexistence equilibrium of (1). Denote such an equilibrium by E 4 = (R 4 , S 4 , M 4 , I 4 , P 4 ). By (5) we have , Since Note that where κ = k /k is the phage mutant rate and the value of κ can represents the relative resistance of bacteria for the CRISPR/Cas immune system. Though there is no definite relation between κ and µ in biology, which is presented as fluctuating selection dynamics [7], in this paper we simply consider two ways for the CRISPR/Cas immune efficiency. First, we assume κ ≥ µ, which means that the mutant rate is great than the cost of bacteria's resistance. Secondly, we consider κ < µ, which means that the mutant rate is inferior to the cost of bacteria's resistance. Note that R 4 is the positive solution of the following equation: Set Using (11) and f (R 4 ) = mR 4 /(a + R 4 ), we see that (13) holds if and only if It is easy to see that there are at most two intersection points for these two curves.
The following Theorem states the existence of infection equilibria of (1) according to the above discussions.
(iv) If 1 < R 0 , there is only one coexistence equilibrium E * * 4 . For ε < ε * we have: Obviously, F (R 4 ) > 0 always holds in case (C1) or case (C2). The conclusions of this theorem follow immediately from the preceding discussions.
This theorem indicates that the system (1) exhibits a backward bifurcation as R 0 crosses unity when ε > ε * for case (C1) or case (C2). Moreover, it admits a forward bifurcation when ε < ε * .
Notice that the equilibrium E * * 4 does not exist for case (C3). We can present another theorem for the existence of infection equilibria.
Theorem 2.6. Let (C3) hold. The existence of infection equilibria is given as follows: The proof of this Theorem is omitted because it is similar to it for Theorem 2.5. Theorem 2.6 presents the conditions for a forward bifurcation of the infection-free equilibrium and a transcritical bifurcation of the coexist equilibrium. Note that for λ 3 ≥ R 0 , there is a backward bifurcation as R 0 crosses unity, but for λ 3 < R 0 , the existence of a backward bifurcation depends on the relation between κ and µ. More specifically, if κ ≥ µ, the backward bifurcation always exists and is independent of the size of b. But for κ < µ, a critical value b * is needed to ensure the existence of backward bifurcation.

2.4.
Global analysis for full CRISPR/Cas immunity. We now explore the persistence and extinction of phages in the case where κ = 0, which leads to k = 0 and means that CRISPR/Cas immunity can fully suppress the invaded phages.
First, we claim that a positive solution of (1) admits S ∞ ≤ R 0 − λ 1 . Indeed, by (2) we getṠ where η 0 > 0 is sufficiently small. Then by similar arguments to those in [42], we obtain S ∞ ≤ R 0 − λ 1 + η 0 . Taking η 0 → 0 leads to the conclusion as claimed. As a result of the claim, we see from the last two equations of (1) that for all large t: where η 1 > 0 is sufficiently small. Let us consider a comparison system: The Jacobian matrix of (18) is Since R 0 < 1, it is easy to examine that the eigenvalues of J 1 have negative real part for small η 1 > 0. Thus, the positive solutions of (18) approach to (0, 0) as t → ∞. It follows from (17) and the comparison principle that a positive solution of (1) exhibits that (I(t), P (t)) → (0, 0) as t → ∞. As a consequence, by the theory of an asymptotically autonomous system [43], it suffices to consider the asymptotic behaviors of the limiting system: Since λ 1 < λ 2 , it follows from [24,25] that the positive solution of (1) exhibits that We wish to show that (1) is uniformly persistent with respect to (X 0 , ∂X 0 ). By Proposition 1, we see that both X and X 0 are positively invariant for model (1). Evidently, ∂X 0 is relatively closed in X. Moreover, Proposition 2 indicates that the nonnegative solutions of model (1) are point dissipative. Then we denote by J ∂ the largest positively invariant set of (1) in ∂X 0 . By similar discussions to those in [46], we see It is clear that there are three equilibria E 0 , E 1 and E 2 in J ∂ . Since λ 1 < λ 2 , by [24,25] we see that E 1 is asymptotically stable in J ∂ and the nonnegative solutions of (1) in J ∂ except for E 0 and E 2 tend to E 1 as t → ∞. Therefore, E 0 , E 1 and E 2 are isolated invariant sets in J ∂ and no subset of {E 0 , E 1 , E 2 } forms a cycle in J ∂ .
Note that (3) and (4) imply that a positive solution of (1) cannot approach E 0 as t → ∞. We claim that this is also the case for E 1 . Suppose not. Then there is a positive solution of (1) that satisfies (R(t), S(t), M (t), I(t), P (t)) → (λ 1 , R 0 − λ 1 , 0, 0, 0) as t → ∞. It follows that for all large t, we havė where η 2 > 0 is sufficiently small. Let Since R 0 > 1, it is easy to examine that J 2 has a positive eigenvalue with a positive eigenvector for small η 2 > 0. It follows that the positive solutions of the following comparison systemİ tend to infinity as t → ∞. As a result, by (20) and the comparison principle we see that the positive solution of (1) satisfies (I(t), P (t)) → (∞, ∞) as → ∞. We are led to a contradiction. Consequently, there is no positive solution of (1) that approaches E 1 as t → ∞. In a similar way, we conclude that there is no positive solution of (1) that approaches E 2 as t → ∞. Hence, W s (E i ) ∩ X 0 = ∅, i = 0, 1, 2. Finally, we claim that E 0 , E 1 and E 2 are also isolated in X 0 . Indeed, if E 1 is not isolated, then there exists a positive solution of (1) which is so close to E 1 for all time that (20) holds. Then repeating the above arguments leads to a contradiction. Similar discussions apply also to E 0 and E 2 . Consequently, we conclude from [44,48] that phage infection is uniformly persistent. By adopting the same techniques as above, we can show that the population S of sensitive bacteria is uniformly persistent. From the third equation and (2), we obtainṀ ≥ −DM + εkSP.
This, together with the uniform persistence of population S and population P , implies that the population M is uniformly persistent.
3. Numerical simulation. In this section, we implement numerical simulations to illustrate the theoretical results and explore more interesting solution patterns of model (1). Take the same parameter values as those in [23] where  Fig.3. For ε = 0.8 > ε * , there are a backward bifurcation from E 1 and a transcritical bifurcation at E 3 , which are shown in Fig.3.
With the help of MatCont package [15], we obtain more information for the dynamical behaviors of system (1), where the behaviors for κ < µ is more complicated than those for κ ≥ µ when ε and R 0 vary. As shown in Fig.2 where case (C2) holds, the coexistence equilibrium E 4 undergoes a Hopf bifurcation and   As discussed above, a bistable coexistence between the infection-free equilibrium and an infection equilibrium may occur when ε > ε * . Panel (a) of Fig.4 shows such a phenomenon. Interestingly, we find the bistable coexistence of the infection-free equilibria with a stable periodic solution, which is shown in panel (b) of Fig.4. 4. Discussions. In this paper, we have developed a bacteriophage mathematical model based on CRISPR/Cas immune system. By combining theoretical analysis and numerical simulations, we have found that the model exhibits some new dynamical behaviors than the model without the immune responses in [23]. More specifically, the introduction of the CRISPR/Cas immune system induces a backward bifurcation from the infection-free equilibrium or a transcritical bifurcation from the coexist equilibrium, which means that although the basic infection reproduction number is below unity, the phage could coexist with bacteria. The coexistence of a stable infection-free equilibrium with a stable infection equilibrium(or stable coexist equilibrium), the bistable phenomenon of a stable infection-free equilibrium and a stable periodic solution are found, which are shown in panel (a) of Fig. 4 and panel (b) of Fig. 4. They provide reasonable explanations for the complexity of phage therapy [9,29] or bacteria-phages coevolution [13], and the coexistence of bacteria with phage in the biological experiments [20,29]. In contrast, there is no the backward bifurcation or bistable phenomena in the models of previous studies [23,42] where the immune response is ignored.
For case (C3), the strain of sensitive bacteria may almost translated to resistent bacteria when the phage mutant rate (means as relative resistance) is inferior to the cost of bacteria's resistance for every values of ε, that is said, the CRISPR/Cas system have uninfluence on the bacteria and phages coexistence or the bacteria's diversity, which only influence by the released virus particles b in this case. This analysis is consistent with [27]. The value of b * increases when µ increases or κ decreases, while the effective of µ is greatly less than the change of κ. Thus, decreasing the values of phage mutant rate is contributed to sensitive bacteria survival when resistant bacteria have great bacteria's resistance cost value.
The mathematical analysis for the stability and bifurcation of equilibria of (1) in this paper present some insights into the underlying phage infection mechanisms by considering the CRISPR/Cas system in bacteria. It will be interesting to consider the analytical conditions for the Hopf bifurcation and the homoclinic bifurcation of the model and reveal how the immune response affect these bifurcations. It will be also interesting to consider the effect of latent period of infection like it in [23] or the nonlinear death rates like those in [29]. We leave these as future researches.
Note that a 1 > 0 and a 4 > 0. We see that the last two inequalities in (9)  It is easy to see F 0 (R 3 ) is positive for any value if κ ≥ µ, which leads to the instability of E 3 . Furthermore, if κ < µ, F 0 (R 3 ) is negative when which is equivalent to R 3 > λ 3 . Consequently, E 3 is asymptotically stable if (9) holds, and is unstable if (10) is available. This completes the proof.