Dynamics model analysis of bacteriophage infection of bacteria

A bacteriophage (in short, phage) is a virus that can infect and replicate within bacteria. Assuming that uninfected and infected bacteria are capable of reproducing with logistic law, we investigate a model of bacteriophage infection that resembles simple SI-models widely used in epidemiology. The dynamics of host-parasite co-extinctions may exhibit four scenarios: hosts and parasites go extinct, parasites go extinct, hosts go extinct, and hosts and parasites coexist. By using the Jacobian matrix and Bendixson–Dulac theory, local and global stability analysis of uninfected and infected steady states is provided; the basic reproduction number of the model is given; general results are supported by numerical simulations. We show that bacteriophages can reduce a host density. This provides a theoretical framework for studying the problem of whether phages can effectively prevent, control, and treat infectious diseases.


Introduction
A bacteriophage (in short, phage) is a virus that can infect and replicate within bacteria. Bacteriophages were discovered by Frederick Twort in 1915 [1] and Felix d'Herelle in 1917 [2]. Bacteriophages consist of a core of genetic material (nucleic acid) surrounded by a protein capsid. Usually a phage follows one of two life cycles: lytic (virulent) or lysogenic (temperate). When lytic phages infect bacteria, they attach their tails to the bacterial surface and then inject DNA from their heads into the bacteria. The DNA that enters the bacteria, being biosynthesised of the raw material of the host bacteria, causes proliferation of new daughter phages. Reaching a certain amount, these phages give rise to bacterial cells bursting, after which they can multiply rapidly and form hundreds of daughter phage particles. Lysogenic phages incorporate their nucleic acid into the chromosome of the host cell and replicate with it as a unit without destroying the cell. Under certain conditions, lysogenic phages can be induced for which a lytic cycle takes place. Each daughter phage can infect another bacterial cell, and the process is repeated over and over again, so that the phages can kill many cells. It has long been noted that phage therapy can be used to treat pathogenic bacterial infections. In particular, phage therapy has been widely used in livestock and poultry breeding (see [3,4]), aquaculture (see [5]), food production (see [6]), and other fields (see [7,8]). Also, phage therapy proved itself as an important tool in the treatment of human diseases such as diarrheal diseases caused by e. coli, shigella or vibrio, and wound infections caused by skin facultative pathogens (for example, staphylococcus and streptococcus). In recent years, phage therapy has also been used for systemic and even intracellular infections. It has been 100 years since the discovery of phage, and the research on phage has never stopped. With the emergence of antibiotics, phage therapy was gradually ignored. However, the widespread antibiotic resistance of bacteria in recent years has made the research on bacteriophages a hot spot. To which extent phages can replace antibiotics as a new treatment for bacterial infections is still under debate, and probably it will take some time for them to appear in clinical trials. However, finding a theoretical basis for studying the relationship between phage and bacteria is a problem of great importance and formidable complexity. This paper attempts to model the dynamic stability of phage infected bacteria, to discuss the relationship between phage and bacteria, and to provide some theoretical basis for whether phages can effectively prevent, control, and treat infectious diseases.
Phages can be thought of as organisms that prey on bacteria, or parasites that host bacteria. To the best of our knowledge, there are many dynamics models of viral infection in host cells (see [9][10][11][12][13][14][15][16][17][18][19][20] and the references therein). For example, in [11], Ebert et al. considered a model of microparasite transmission for a horizontally transmitted parasite and established the host-born density-dependent cabin model by ignoring the possibility of host recovery. The model equations are where x, y are the densities of uninfected (susceptible) and infected (infective) hosts at time t, respectively; r is the maximum per capita birth rate of uninfected hosts; f is the relative fecundity of infected hosts; c measures the per capita density-dependent reduction in birth rate; μ is the parasite-independent host background mortality; β is the infection rate constant; and ν is the parasite-induced excess death rate.
This model predicts the existence of a stable equilibrium of infected and uninfected hosts, and the population is predicted to approach this equilibrium either monotonically or by damped oscillations.
In this paper, we modify Ebert's model by assuming that infected hosts (bacteria) are capable of reproducing with logistic law, and investigate a class of parasites (phages) infection models. The model is given by the following system of differential equations: Here, r 1 , r 2 are the proliferation constants of uninfected and infected hosts, M is the environmental tolerance of a host population, d 1 stands for the phage-independent bacteria background mortality, β is the proportionality coefficient of parasite infection, ε denotes the phage-induced excess death rate. In this model, we assume that the total hosts population N is composed of two population classes: the first one consists of uninfected hosts (denoted S), while the second one consists of the virus infected hosts (denoted I). Thus, N(t) = S(t) + I(t). We make the following assumption: both uninfected and infected hosts are capable of reproducing with logistic law, and the logistic growth of the uninfected bacteria and infected bacteria is given by r 1 S(t) [

-(S(t) + I(t))/M] and r 2 I(t)[-(S(t) + I(t))/M], respectively.
To establish our results, the existence and number of steady states, as well as local stability and global stability of uninfected and infected steady states, are analyzed by the Jacobian matrix and Bendixson-Dulac theory. Under certain assumptions (reasonable from the biological viewpoint), we derive the basic reproduction number R 0 : E 1 is locally asymptotically stable if R 0 < 1, and E 1 is unstable if R 0 > 1. Our theoretical discoveries are supported by numerical simulations.
At present, there is not much work on mathematical modeling of phage infection bacteria, so our study has certain significance. In particular, the parasite (phage) population is not explicitly modeled in this model, which is a striking feature of this model.
The organization of this paper is as follows. In Sect. 2, we discuss the positively invariant set and equilibria. In Sect. 3, we give local and global stability analysis. In Sect. 4, we derive the basic reproduction number. Then, in Sect. 5, we give numerical simulations to support our main result. Section 6 ends the paper with a discussion about phage therapy that could be a new savior for patients infected with superbugs.

Positively invariant set and equilibria
To begin with, let us find a positively invariant set with respect to (2). Adding the equations in (2), one obtainṡ where r = max{r 1 , r 2 }. Clearly, is positively invariant with respect to (2). Next, we will study equilibria of system (2) (i.e. zeros of the right-hand side of (2)). It is easy to verify that system (2) admits one vanishing equilibrium and two boundary equilibria. More precisely, the following statement (its proof is straightforward) is true.

Proposition 2.2 Under the assumptions of Proposition 2.1, suppose, in addition, that
, provided that one of the following conditions holds: Proof Consider the algebraic system If r 1r 2 + Mβ = 0, then (4) Since (6) implies ω > 0 and S > I, it is easy to see that if β 1 < β < ω r 2 I or ω r 2 I < β < β 2 , then (2) admits the unique positive equilibrium, and the result follows. We assume that the hypotheses of Propositions 2.1-2.2 are satisfied. The Jacobian matrix of (2) is given by hence, By assumption, r 1d 1 > 0 and r 2d 2 > 0, therefore, E 0 = (0, 0) is always an unstable node.

Local stability of the positive equilibrium
It follows from (7) that the Jacobian matrix of (2) at the positive equilibrium point E * = (S * , I * ) is given by By direct computation, In particular, if β > r 2 -r 1 M (resp. β < r 2 -r 1 M ), then det J(E * ) > 0 (resp. det J(E * ) < 0). This way, we arrive at the following result. (i) If β > r 2 -r 1 M and ω r 1 S < β < ω r 2 I , then E * is locally asymptotically stable; (ii) If β < r 2 -r 1 M , then E * is a saddle and, as such, is unstable.

Global stability analysis
We complete this section with the following global stability result. Proof Take μ(S, I) = 1 SI to serve as a Dulac multiplier. Denote

Theorem 3.3 Under the notations and assumptions of Propositions
x -βSI and Q := βSI + r 2 I 1 - Then for all S > 0, I > 0. Thus, (2) has no periodic orbits in (see (3)). A simple application of the classical Poincaré-Bendixson theory shows that all solutions in converge to a single equilibrium, from which the result follows.

Basic reproduction number
In this section, using the method of the next generation matrix (see [21]), we derive the basic reproduction number R 0 (β) of model (2). We have and the next generation matrix is Hence, the basic reproduction number is given by Clearly (see (13)), R 0 (β) is increasing in β, and R 0 (β c ) = 1, where β c = ω r 1 S . It is easy to see that (i) If R 0 (β) < 1, then E 1 is locally asymptotically stable; (ii) If R 0 (β) > 1, then E 1 is a saddle and, as such, unstable.
According to the stability results presented in Sect. 3, the above choice of parameter values guarantees that the equilibrium E 2 is asymptotically stable when β > ω r 2 I . Figure 1(2) shows that, for β = 0.052, the solution approaches the boundary equilibrium (0, I). parasite-independent hosts background mortality 0.01 [11] β humoral response coefficient of a susceptible host population 0.001-0.1 [11] ε parasite-induced excess death rates 0.2 [11]   Next, taking β = 0.039 mm 3 /cells/day and leaving other model parameters unchanged, we provide numerical simulations related to the positive equilibrium E * (see Fig. 2). These simulations are consistent with the theoretical results related to the local asymptotic stability of E * obtained in Sect. 3, see Table 2.

Discussion
A phage is a virus that lives in bacteria and can infect lysed host bacteria in a specific environment. It has been successfully used to treat infections caused by e. coli, pseudomonas aeruginosa, and staphylococcus (see [22][23][24]). In recent years, many scholars (both domestic and foreign) conducted in-depth research on phage biology and genome sequencing. It was confirmed that phages of many strains can effectively kill mycobacterium tuberculosis, see [25,26] and the references therein. Also, it was predicted that this "bacterial killer" may become an alternative to antibiotics in the future. There is a complex dynamic relationship between phage and bacteria, which has been described as a variety of dynamic behaviors at the population level, including ARD (arm race dynamics) model [27], FSD (fluctuating selection dynamics) model [27], KtW (kill-the-winner) model [28], and PtW (piggyback-the-winner) model [29]. The study of the dynamic changes of phage and bacteria is helpful to elucidate their changing rules and lay a foundation for the better application of phage in medical treatment, animal husbandry, and other industries. However, there is a lack of research on mathematical modeling of phage infection. Beretta and Kuang [30,31] proposed a mathematical model of phage infection in the ocean and analyzed its mathematical characteristics based on the observations of marine organisms made by biologist A. Okubo. In the present paper, we revised the mathematical model established by Ebert et al. [11] by assuming that (a) infected hosts (bacteria) are capable of reproducing with logistic law, and (b) the relative fecundity of an infected host is equal to zero.
In particular, we studied the nonlinear dynamic model of vertical transmission of infection bacteria. The parasite (phage) population is not explicitly modeled, which is a striking feature of this model. As such, our model is similar to SI-models widely used in epidemiology.
We established four possibilities for the host-parasitic relationships: both uninfected and infected hosts become extinct simultaneously; extinction of uninfected hosts; extinction of infected hosts; uninfected and infected hosts coexist. The local and global stability of the four equilibrium points was analyzed by using the Jacobian matrix and Bendixson-Dulac theory. In particular, the instability of the extinction equilibrium was established. We also figured out the basic reproduction number R 0 of the system: E 1 is locally asymptotically stable (resp. unstable) if R 0 < 1 (resp. R 0 > 1).

Figure 3
Here, r 1 = 0.4, r 2 = 0, d 1 = 0.01, ε = 0.2, M = 10. In Fig. 3(4), β = 0.0238, solutions tend to E * monotonically. In Fig. 3(5), β = 0.158, solutions tend to E * by damped oscillations The model simulations presented in Fig. 2 show the existence of the positive equilibrium to which population approaches monotonically. The model simulations presented in Fig. 3 show that when r 2 = 0, the population approaches the equilibrium either monotonically or by damped oscillations. Actually, it is a special case of system (1), and the results are consistent with Ebert's prediction in [11]. We established that phages and bacteria can coexist, and the concentration of bacteria can reduce. Also, phages cause the increase in the death rate of the bacterial host, keeping the bacteria at a low concentration that prevents them from becoming ill and spreading to others. Finally, we predict that phage therapy could be a new savior for patients infected with superbugs.