Stability analysis and persistence of a phage therapy model

1.   Introduction

Phage therapy is the use of bacteriophage as a therapeutic agent for the treatment of bacterial infections and has existed since the early 20th century. A bacteriophage is a virus that infects and lyses bacteria. This type of virus does not harm other cells and can kill bacteria. The killing ability of a bacteriophage is a remarkable feature that can be employed in treating bacterial disease. Bacteriophages also have many advantages over antibiotics. Accordingly, bacteriophages are used as potential agents to substitute for antibiotics for curing bacterial diseases. Accordingly, numerous bacteriophage models have been established in recent decades, and many valuable results have been obtained.

Many papers on predator–prey models are available for the study of the stability and other dynamic behavior of various mathematical models. Some significant studies on the Leslie–Gower predator–prey models can be found in ^[1,2,3]. In ^[1,2], the authors investigated the global stability of the interior equilibrium with help from the Lyapunov function. In ^[3], the authors explored the positivity, boundedness, existence, and stability of various equilibria in addition to the Hopf bifurcation. Kar ^[4,5] investigated the stability and the different dynamic behavioral characteristics of the prey–predator model as regards the Holling type II functional response. Other prey–predator fishery models with various types of functional responses and which investigate the feasible steady states together with their existence and stability were discussed in ^[6,7,8]. The authors in ^[9,10,11] presented some studies related to the teaming approach of the prey–predator model. The persistence, permanence criteria of the system, and local and global behavior of the different equilibrium solutions were depicted in these articles. Many other variants for studying the persistence of the system, existence, and local and global dynamics at all the possible equilibria can be found in ^{[12,13,14,15,16,17,18,19]}.

Campbell ^[20] studied the predator–prey association between bacteriophage and bacteria. This association was developed by Levin et al. ^[21] to propose a general chemostat model on resource-limited growth, predation, and competition. Aviram and Rabinovitch ^[22] presented a mathematical model concerning the coexistence of bacteria and bacteriophage in a chemostat. Refer to ^[23,24] for the study of the persistence and extinction of bacteria and their resistant strain in the chemostat. The authors in ^[25,26] introduced the other significant work that examines the boundedness, permanence, existence, and local and global stability of the chemostat model of bacteria and virulent bacteriophage. Sahani and Gakkhar ^[27] provided an impulsive phage therapy model. They identified two important parametric conditions for curing bacterial disease under certain conditions. The authors in ^{[28,29,30,31,32]} analyzed mathematical models for the interactions in marine bacteriophage infections. Gakkhar and Sahani ^[33] established the coexistence of bacteria, bacteriophage, and infected bacteria. They provided the conditions for the existence and stability of susceptible bacteria-free equilibrium and also considered a simple Hopf-bifurcation for non-zero equilibrium point. Calsina et al. ^[34] introduced a structured cell-population model for the interaction of bacteria and phages, and computed the optimal lysis timing (latent period).

The research on the interactions among bacteria, phages, and the immune system is vital for the reasonable use of bacteriophage treatments. Meanwhile, bacterial elimination using bacteriophage is potentially beneficial and can be used in curing bacterial infection ^[35]. Therefore, mathematical models regarding the phage therapy that combine the nonlinear interactions of bacteria, phages, and the immune system have attracted more attention from authors ^{[36,37,38,39,40]}. For example, Wang ^[38] extended the basic mathematical model of bacteria and phages in a chemostat that was proposed by ^[23] to include host innate immunity. The author investigated the effects of the host immune response on the dynamics of the model. Shu et al. ^[39] investigated a bacteriophage model based on the adaptive immune system in the bacteria to examine the stability analysis and bifurcation of equilibria. Leung and Weitz ^[40] proposed a mathematical model for phage therapy that incorporates the interactions of bacteria, phages, and the immune system to identify a synergistic regime whereby the phage and immune system jointly contribute to the elimination of bacteria. They also show that the synergy between the phage and the immune system is crucial for effective phage therapy in eliminating bacterial pathogens.

The model in ^[40] was presented by a system of nonlinear ordinary differential equations as follows:

with the initial conditions

Here, all the parameters are positive. $B(t)$, $P(t)$, and $I(t)$ denote the population densities of the bacteria, phages, and innate immune response at time $t$, respectively. $r$ and $\alpha$ are the intrinsic growth rates of the bacteria and the innate immune response. $K_C$ and $K_I$ are the carrying capacities of the bacteria and the innate immune response. $\beta$, $\phi$, and $\omega$ indicate the burst size, adsorption rate, and decay rate of the phage, respectively. $\epsilon$ is the killing rate of the innate immune response, $K_N$ is the bacterial concentration when the innate immune response growth rate is half its maximum, and $K_D$ is the bacterial concentration when the innate immune response is half saturated. In model (1.1), the immune system that kills bacteria is activated when the bacteria are persistent. The phage can infect and lyse the bacterial population reproduction. Moreover, phage particles can decompose at the outside of cells.

The model formulated by the authors in ^[40] is useful in clinical trials for bacterial infections. However, in this study, the conditions for the existence of all the equilibria and their stability behavior and the persistence and extinction of model (1.1), were not discussed. However, these features are biologically and ecologically crucial. When the equilibrium points are stable, their description summarizes the biologically significant aspects of the model. When the said points are unstable, it must be established whether stable cycles of population fluctuation occur or whether the instability leads to ever increasing fluctuations with the eventual extinction of at least one of the populations ^[21]. The persistence of populations must be sustained to maintain the balance of an ecosystem in the real world. Note that the persistence of a system indicates that all the species are present and none of them will face extinction.

Given the above observations, we examine the dynamic behavior of the phage therapy model, including the existence of all feasible equilibria, their stability, persistence and nonpersistence of system (1.1) and verify the theoretical results through numerical simulation. From these studies, the necessary and sufficient conditions that have biologically compelling interpretations for bacterial persistence and phage extinction are obtained. Moreover, we numerically confirmed the effect of intrinsic growth rate of bacteria and the immune killing rate on the persistence and extinction of the bacteria and phages.

The study is organized as follows. In Section 2, we describe the positivity and boundedness of the solutions of system (1.1). In Section 3, we provide all the feasible equilibria of the system, their existence, and local stability analysis. In Section 4, we establish the global stability analysis of the equilibrium solutions $E_3$, $E_4$, and $E_5$ under certain conditions. Moreover, we examine the uniform persistence of the model system (1.1) and phage extinction. In Section 5, we present and discuss the numerical simulation results. In Section 6, the conclusion of this study is presented.

2.   Positivity and boundedness

In this section, we analyze the positivity and boundedness of the solutions of model (1.1).

Theorem 2.1. (i) All solutions $(B(t), P(t), I(t))$ of system (1.1) with initial condition (1.2) exist in the interval $\left[0, \infty\right)$ and satisfy $B(t)\geq 0$, $P(t)\geq 0$, $I(t)\geq 0$, and $\forall t\geq 0$.

(ii) All the solutions of system (1.1) with initial condition (1.2) are bounded for all $t\geq 0$.

Proof. (ⅰ) Let us define the right side of model (1.1) as function $g$. Clearly $g$ is continuous. Therefore, $g$ is locally Lipschitz on $\mathbb{R}_+ ^{3} = \left\{(B, P, I):B\geq 0, P\geq 0, I\geq 0\right\}$. Thus, the solution $(B(t), P(t), I(t))$ of system (1.1) with (1.2) exists and is unique on $\left[0, \zeta\right)$, where $0 < \zeta \leq +\infty$. Under model system (1.1) with (1.2), we obtain

which completes the proof.

(ⅱ) We consider $W(t) = \beta B(t)+P(t)$. Then, differentiating $W$ w.r.t $t$ along the trajectories of model (1.1) yields

Hence,

We obtain the following expression through the theorem of differential inequality ^[41]:

This expression shows that $B(t)$ and $P(t)$ are bounded. Now, we will ascertain the boundedness of $I(t)$. Given that variables $B$ and $I$ are positive, we obtain the following expression through the third equation of (1.1):

Changing variable $u(t) = \frac{1}{I(t)}$ yields

Hence, both sides of Eq (2.1) are multiplied by the integrating factor $\frac{{\rm e}^{\alpha t}}{K_I}$. Then, we obtain the following expression after the integration:

Consequently,

Therefore, $I(t)$ is bounded in its maximal domain. This expression completes the proof.

3.   Existence and local stability analysis of equilibria

In this section, the existence and local behavior of all the possible equilibrium points of system (1.1) are considered.

3.1.   Existence of equilibrium points

In this subsection, we present the existence of various equilibrium solutions of model (1.1). Straightforward calculations reveal that the possible equilibria of system (1.1) are as follows:

1. Trivial equilibrium: $E_0 = (0, 0, 0)$.

2. Axial equilibrium: (i) $E_1 = (K_C, 0, 0)$ and (ii) $E_2 = (0, 0, K_I)$.

3. Planar equilibrium:

(ⅰ) $E_3 = (\bar{B}, \bar{P}, 0)$, where $\bar{B} = \frac{\omega}{\beta\phi}$, $\bar{P} = \frac{r}{\phi}\left(1-\frac{\omega}{\beta\phi K_C}\right)$ with

(ⅱ) $E_4 = (B^{\prime}, 0, I^{\prime})$, where

with

4. Interior equilibrium: $E_5 = (B^*, P^*, I^*)$, where

with

3.2.   Local stability of equilibria

In this subsection, we analyze the local stability of all possible equilibria of model (1.1) by calculating the corresponding variational matrices of each equilibrium point. The findings are shown by the following theorem:

Theorem 3.1. (i) The trivial equilibrium $E_0 = (0, 0, 0)$ of system (1.1) is unstable;

(ii) The axial equilibrium $E_1 = (K_C, 0, 0)$ is unstable;

(iii) The axial equilibrium $E_2 = (0, 0, K_I)$ is stable if $r < \epsilon K_I$;

(iv) The planar equilibrium $E_3 = (\bar{B}, \bar{P}, 0)$ is unstable;

(v) The planar equilibrium $E_4 = (B^{\prime}, 0, I^{\prime})$ is locally asymptotically stable if

(vi) The positive interior equilibrium $E_5 = (B^*, P^*, I^*)$ is locally asymptotically stable if it exists, and

i.e., the intrinsic growth rate of bacteria exceeds a threshold value.

Proof. (ⅰ) The variational matrix of system (1.1) at $E_0 = (0, 0, 0)$ is presented by

The roots of the characteristic equation of $J(E_0)$ are $\lambda_1 = r > 0$ and $\lambda_2 = -\omega < 0$. This concept implies that $E_0$ is unstable.

(ⅱ) The variational matrix of system (1.1) at $E_1 = (K_C, 0, 0)$ is expressed by

Then, the roots of the characteristic equation of $J(E_1)$ are $\lambda_1 = -r < 0$, $\lambda_2 = \beta\phi K_C-\omega$, and $\lambda_3 = \frac{\alpha K_C}{K_C+K_N} > 0$. $\lambda_3$ is positive, an outcome which implies that $E_1$ is unstable in the $B \rm{-}P \rm{-}I$ space. If $\lambda_2 > 0$ (i.e., $\omega < \beta\phi K_C$), then $E_1$ is unstable in the $P \rm{-}I$ plane and stable in the $B$ direction. However, if $\lambda_2 < 0$ (i.e., $\omega > \beta\phi K_C$), then $E_1$ is stable in the $B \rm{-}P$ plane and unstable in the $I$ direction. We can describe this outcome biologically as follows: if the carrying capacity of bacteria is less than a threshold value, then equilibrium without the phage and immune system is locally asymptotically stable in the $B \rm{-}P$ plane. This threshold value depends on the phage parameters (decay rate, burst size, and adsorption rate). Given the existence condition of $E_3$ (i.e., $\lambda_2 = \beta\phi K_C-\omega > 0$), $E_1$ is unstable in $B \rm{-}P$ plane. Thus, $E_1$ is unstable when $E_3$ exists.

(ⅲ) The variational matrix of system (1.1) at $E_2 = (0, 0, K_I)$ is denoted by

Thus, the roots of the characteristic equation of $J(E_2)$ are $\lambda_1 = r-\epsilon K_I$ and $\lambda_2 = -\omega < 0$. Hence, $E_2$ is stable if $\lambda_1 < 0$, a condition which implies that $r < \epsilon K_I$. This outcome can be described biologically as follows: if the maximum bacterial growth rate $r$ is less than the maximum per capita immune killing rate $\epsilon K_I$, then equilibrium without the bacteria and phage is stable. Given the existence condition of $E_4$ (i.e., $\lambda_1 = r-\epsilon K_I > 0$), $E_2$ is unstable. Hence, $E_2$ is unstable if $E_4$ exists.

(ⅳ) The variational matrix of system (1.1) at $E_3 = (\bar{B}, \bar{P}, 0)$ is provided by

The roots of the characteristic equation of $J(E_3)$ are $\lambda_1 = \frac{\alpha\omega}{\omega+\beta\phi K_N}$ and the solutions of the quadratic equation,

where

Thus, if the existence condition of $E_3$ (3.1) holds, then $E_3$ is locally asymptotically stable in the $B \rm{-}P$ plane. Given that $\lambda_1 = \frac{\alpha\omega}{\omega+\beta\phi K_N}$ is positive, $E_3$ is unstable in the $B \rm{-}P \rm{-}I$ space. We observed that if the carrying capacity of the bacteria is greater than a threshold value, then equilibrium without the immune system is locally asymptotically stable in the $B \rm{-}P$ plane. This threshold value depends on the phage parameters (decay rate, burst size, and adsorption rate).

(ⅴ) The variational matrix of system (1.1) at $E_4 = (B^{\prime}, 0, I^{\prime})$ is defined by

The roots of the characteristic equation of $J(E_4)$ are $\lambda_1 = r-\frac{2r}{K_C}B^{\prime}-\frac{\epsilon K_I K_D^2}{(B^{\prime}+K_D)^2}$, $\lambda_2 = \beta \phi B^{\prime}-\omega$, and $\lambda_3 = -\frac{\alpha B^{\prime}}{B^{\prime}+K_N} < 0$. $E_4$ is locally asymptotically stable in the $B \rm{-}I$ plane if $\lambda_1 < 0$ (i.e., $r < \frac{\epsilon K_I K_C K_D^2}{(K_C-2B^{\prime})(B^{\prime}+K_D)^2}$). Moreover, $E_4$ is locally asymptotically stable in the $B \rm{-}P \rm{-}I$ space if $\lambda_1 < 0$ and $\lambda_2 < 0$, provided that $r < \frac{\epsilon K_I K_C K_D^2}{(K_C-2B^{\prime})(B^{\prime}+K_D)^2}$ and $\omega > \beta\phi B^{\prime}.$ This outcome can be biologically interpreted as follows: if the intrinsic growth rate and equilibrium density of the bacteria is lower than some threshold values, then equilibrium without phages is locally asymptotically stable in the $B \rm{-}P \rm{-}I$ space.

(ⅵ) The variational matrix of system (1.1) at the positive equilibrium $E_5 = (B^*, P^*, I^*)$ is provided as follows:

Hence, $\lambda_1 = -\frac{\alpha B^*}{B^*+K_N}$ is one of the roots of the characteristic equation of $J(E_5)$. The other two roots are provided by equation

where

Substituting the value of $B^*$, $P^*$, and $I^*$ in (3.5) yields

$\lambda_1 = -\frac{\alpha\omega}{\omega+\beta\phi K_N} < 0$. The existence condition (3.4) of $E_5 = (B^*, P^*, I^*)$ provided $N_2 > 0$. Thus, $E_5$ is locally asymptotically stable in the $B \rm{-}P \rm{-}I$ space whenever it exists and $N_1 > 0$, thereby implying that $r > \frac{\epsilon\beta^2\phi^2K_I K_C K_D}{(\omega+\beta\phi K_D)^2}.$ This outcome can be biologically described as follows: if the intrinsic growth rate of bacteria exceeds a threshold value, then coexistence equilibrium $E_5$ is locally asymptotically stable. This expression completes the proof.

4.   Global stability analysis and uniform persistence

In this section, the global behavior at equilibria $E_3$, $E_4$ and $E_5$ of system (1.1) is established under certain parametric conditions. We also discuss some conditions for the persistence and nonpersistence of model (1.1).

4.1.   Global stability

In this subsection, we first present the global stability of $E_3$ and $E_4$ by applying the Bendixson–Dulac criterion. Then, we explore the global stability of $E_5$ by using an appropriate Lyapunov function.

Theorem 4.1. If equilibrium without immune system $E_3$ exists and is locally asymptotically stable in the interior of the positive quadrant of the $B \rm{-}P$ plane, then $E_3$ is globally asymptotically stable in that plane.

Proof. Let us construct

$\psi(B, P) > 0$ in the interior of the positive quadrant of the $B \rm{-}P$ plane. Thus, we obtain

$\Delta(B, P)$ has no sign change and is not zero in the positive quadrant of the $B \rm{-}P$ plane. When the Bendixson–Dulac criterion is applied, system (1.1) has no limit cycle and completely lies in the positive quadrant of the $B \rm{-}P$ plane. Therefore, $E_3$ is globally asymptotically stable.

Theorem 4.2. If equilibrium without phages $E_4$ exists and is locally asymptotically stable in the $B \rm{-}I$ plane, then $E_4$ is globally asymptotically stable in region $\mathbb{R}_+ ^{2}$ of the $B \rm{-}I$ plane, where

Proof. Let us construct

$\psi(B, I) > 0$ in the interior of the positive quadrant of the $B \rm{-}I$ plane. Hence, we obtain

Therefore, the theorem holds.

Theorem 4.3. If $\beta < 1$, then the positive interior equilibrium $E_5 = (B^*, P^*, I^*)$ is globally asymptotically stable in the interior of the positive octant ($i.e., \ Int.\mathbb{R}_+ ^{3}$).

Proof. Define the positive definite Lyapunov function $V(B, P, I):\mathbb{R}_+ ^{3} \rightarrow \mathbb{R}$ to validate the global stability at $E_5 = (B^*, P^*, I^*)$, such that

where $V_1(B, P, I) = (B-B^*-B^*\ \rm{ln}(B/{B^*}))$, $V_2(B, P, I) = (P-P^*-P^*\ \rm{ln}(P/{P^*}))$, and $V_3(B, P, I) = (I-I^*-I^*\ \rm{ln}(I/{I^*}))$.

$V(B, P, I)$ is continuous on $Int.\mathbb{R}_+ ^{3}$ and zero at $E_5 = (B^*, P^*, I^*)$. Then, differentiating function $V$ with respect to time $t$ along the trajectories of (1.1) yields

Moreover, the time derivatives of $V_1$, $V_2$, and $V_3$ along the solutions of system (1.1) are as follows:

$E_5 = (B^*, P^*, I^*)$ satisfies (1.1). Thus, we obtain the following expression after a straightforward computation:

Substituting the three values of (4.6) into (4.3)-(4.5) yields

Substituting (4.7)-(4.9) into (4.2) and using algebraic calculation yield

If the condition in Theorem 4.3 holds, then $\frac{dV}{dt} < 0$ along all the trajectories in $\mathbb{R}_+ ^{3}$, except $E_5 = (B^*, P^*, I^*)$. Therefore, $E_5 = (B^*, P^*, I^*)$ is globally asymptotically stable.

Note 1. Theorem 4.3 shows that burst size $\beta$ (the amount of new virions released per lysis) plays an important role in making system (1.1) globally asymptotically stable.

4.2.   Uniform persistence

In this subsection, we establish the uniform persistence and non-persistence behaviors of model (1.1). We show the uniform persistence of the model using the average Lyapunov function method.

Theorem 4.4. Assume that the hypotheses in Theorems 4.1 and 4.2 hold. Then, model (1.1) is permanent or uniformly persistent if $\omega < \beta\phi K_C$, $r > \epsilon K_I$, and $\omega < \beta\phi \left(\frac{K_C-K_D}{2}+ \sqrt{\frac{(K_C+K_D)^2}{4}-\frac{\epsilon K_I K_C K_D}{r}}\right)$.

Proof. We define the average Lyapunov function for (1.1) as follows:

where $\delta$, $\delta_1$, and $\delta_2$ are positive constants. Function $\Upsilon(X)$ is the nonnegative continuous function defined in $\mathbb{R}_+ ^{3}$. Therefore, we obtain

We assume that conditions (3.1), (3.2), and (3.3) hold. The hypotheses in Theorems 4.1 and 4.2 hold. Then, planar equilibria $E_3$ and $E_4$ exist. Furthermore, no periodic orbits are observed in the interior of the $B \rm{-}P$ plane and region $\mathbb{R}_+ ^{2}$ of the $B \rm{-}I$ plane. Hence, system (1.1) is uniformly persistent enough to prove that $\Omega(X)$ is positive for all equilibria $X$ in domain $G$ of $\mathbb{R}_+ ^{3}$ where

for the appropriate selection of $\delta$, $\delta_1$, and $\delta_2 > 0$. Specifically, the conditions described below must be satisfied for system (1.1) to persist.

We can make $\Omega(E_0) > 0$ by increasing $\delta$. Hence, inequality (4.10) holds. Inequality $\omega < \beta\phi K_C$ implies that (4.11) also holds. If $r > \epsilon K_I$, then increasing $\delta$ is enough. We can make $\Omega(E_2) > 0$, implying that inequality (4.12) holds. $\frac{\alpha\omega}{\omega+\beta\phi K_N}$ is positive. Thus, inequality (4.13) holds. The following condition must be satisfied for inequality (4.14):

Finally, we know that model (1.1) is permanent or uniformly persistent.

Note 2. Theorem 4.4 indicates that system (1.1) is uniformly persistent if the bacteria growth rate is greater than the innate immune killing rate, and if the carrying capacity of bacteria and the equilibrium density of the bacteria are greater than the threshold values that depends upon the phage parameters.

Remark 1. For the persistence of system (1.1), when conditions (a) $\omega < \beta\phi K_C$, (b) $r > \epsilon K_I$, and (c) $\omega < \beta\phi \left(\frac{K_C-K_D}{2}+ \sqrt{\frac{(K_C+K_D)^2}{4}-\frac{\epsilon K_I K_C K_D}{r}}\right)$ hold, the equilibrium without phage and immune system $E_1 = (K_C, 0, 0)$ becomes unstable in the $B \rm{-}P$ plane, while the equilibrium without bacteria and phage $E_2 = (0, 0, K_I)$ and the equilibrium without phages $E_4 = (B^{\prime}, 0, I^{\prime})$ both become unstable in the $B \rm{-}P \rm{-}I$ space.

Now, we provide a sufficient condition under which system (1.1) is non-persistent. The following lemma must be recalled.

Lemma 4.5. (see ^[8,42]) If $k_1$, $k_2 > 0$, and $\frac{dX}{dt}\leq(\geq) X(t) (k_1-k_2X(t))$ with $X(0) > 0$, then

Theorem 4.6. If $\omega > \beta \phi K_C$, then $\lim \limits_{t\rightarrow \infty}P(t) = 0$, that is, the phage becomes extinct.

Proof. Applying the positivity of variables $B$, $P$, and $I$ and the first equation of (1.1) yields

Using Lemma 4.5 yields

Hence, $T^* \in \mathbb{R}_+$ for arbitrary $\eta > 0$. Accordingly,

Using the second equation in (1.1) and (4.16) yields

Then,

With the given hypothesis, $P(t)\to 0$ as $t\to\infty$. Specifically, the phage becomes extinct.

Note 3. Biologically, Theorem 4.6 indicates that when the carrying capacity of bacteria is less than the threshold value ($K_C < \omega/\beta \phi$), (i.e., a high decay rate, low adsorption rate, and small burst size), the phage disappears.

5.   Numerical simulation

In this section, we describe the numerical simulations to explain our analytical findings and stability results in the prior sections. For this, we consider the three sets of parameter values of system (1.1) as provided in Table 1.

For the set of parameter values in Data 1, the conditions of persistence in Theorem 4.4 are satisfied. For these parameter values, all the species, namely, $B(t)$, $P(t)$, and $I(t)$, persist, and a stable population is obtained for all the species (Figure 1). Figure 1 indicates that the densities of phage species $P$ and immune system $I$ initially decrease, while the bacteria species slightly decreases initially and then increase slowly. Finally, all three species achieve steady states and become asymptotically stable. We also plot the dynamics of the model for different initial conditions using these parameter values (Figure 2). In Figure 2(a)–(c), the populations of all the species, namely, $B(t)$, $P(t)$, and $I(t)$, tend to be in steady state. In Figure 2(d), all the solutions of system (1.1) approach $E_5 = (2.5, 2.091, 0.5)$, starting from various initial points. Thus, coexistence equilibrium point $E_5 = (2.5, 2.091, 0.5)$ becomes an attractor.

In Figure 3, we show the variation of $B$ and $P$ with time for five different values of parameter $r$, and the rest of the parameters have the same values as those in Data 1. More precisely, if $r = 0.6$, $r = 1$, and $r = 1.5$, such that $r > \epsilon K_I$, then all bacterial and phage populations converge to their equilibrium values, implying that system (1.1) becomes persistent. Meanwhile, if $r = 0.05$ and $r = 0.15$, such that $r\leq\epsilon K_I$, then system (1.1) is nonpersistent, that is, leads to the extinction of some species because the populations of $B$ and $P$ tend toward zero. Figure 3 shows that whenever the intrinsic growth rate of bacteria exceeds the immune killing rate, the bacteria and phage populations persist; otherwise, the bacteria and phage populations become extinct. This phenomenon indicates that the intrinsic growth rate of bacteria and the immune killing rate may suppress the persistence and extinction of bacteria and phage.

For the set of parameters stated in Data 1, except $\omega = 0.3$, the extinction condition of phage $P$ given in Theorem 4.6 is satisfied. We observe that if the carrying capacity of bacteria is less than a certain threshold value, phage population $P$ becomes extinct, while bacterial population $B$ and the population of innate immune response $I$ persist (see Figure 4). In Figure 4(a), the phage population significantly decreases and ultimately tends toward zero, while bacteria population $B$ increases toward the equilibrium level, immune system $I$ decreases toward the steady state level. In Figure 4(b), we also describe the solution curves starting from different initial points. This figure shows that all the phage populations ($P$) tend toward zeros, implying phage extinction for system (1.1).

In the set of parameters exhibited in Data 2, the stability condition of $E_2$ in Theorem 3.1(iii) is satisfied. We see that the equilibrium without bacteria and phage $E_2 = (0, 0, 2.48)$ is stable if the maximum bacterial birth rate is less than the maximum per capita immune killing rate, as described in Figure 5. Figure 5 indicates that species $B$ and $P$ are extinct, while species $I$ persists. The stability conditions of $E_4$ in Theorem 3.1(v) are satisfied by further changing the parameter values $\epsilon = 0.4$ and setting the other parameters similarly as in Data 2. In this case, system (1.1) has an equilibrium without phages $E_4 = (2.058, 0, 2.48)$, and we observe that $E_4$ is locally asymptotically stable if the intrinsic growth rate of bacteria and the equilibrium density of the bacteria are less than the threshold values (see Figure 6). Figure 6(a) illustrates that in the parametric values stated in Data 2, except $\epsilon = 0.4$, phage $P$ eventually becomes extinct, while bacteria $B$ and innate immune response $I$ persist and finally obtain their steady states $(2.058, 0, 2.48)$. Figure 6(b) shows that all the solutions, beginning from various initial conditions, approach $E_4 = (2.058, 0, 2.48)$.

In the set of parameters stated in Data 3 of Table 1, system (1.1) has a coexistence equilibrium $E_5 = (1.705, 1.038, 0.58)$. According to Theorem 3.1(vi), we conclude that $E_5$ is locally asymptotically stable if the intrinsic growth rate of bacteria species is greater than the threshold value (see Figure 7). Figure 7(a) shows that the bacteria and phage populations initially exhibit oscillations, and the amplitude of oscillations gradually decreases and eventually becomes stable. For species $I$, a steady state achieved. Figure 7(b) represents the stable phase portrait of system (1.1).

6.   Conclusions

This paper focuses on the dynamical analysis of the phage therapy model (1.1) proposed by Leung and Weitz in ^[40]. On the basis of this model, we develop the mathematical analysis of the model theoretically and numerically. First, we consider the basic dynamic behaviors, such as the positivity and boundedness of system (1.1). Theorem 2.1 shows that all the solutions of system (1.1) are positive and bounded, implying that the system is biologically well–behaved. Then, we discuss the sufficient conditions for the existence and local stability of all the equilibrium solutions of the system.

From Theorems 3.1 (ⅱ) and (3.1), we conclude that the instability of the equilibrium without phage and immune system $E_1$ offers a sufficient condition for the existence of equilibrium without immune system $E_3$. Similarly, we can deduce that the instability of the equilibrium without bacteria and phage $E_2$ provides an existence condition for the equilibrium without phages $E_4$ (refer to Theorems 3.1 (ⅲ) and (3.3)).

For the equilibrium without phage and immune system $E_1$ to be locally asymptotically stable (LAS) in the $B \rm{-}P$ plane, the carrying capacity of bacteria must be less than the threshold value, which depends on the decay rate, burst size, and adsorption rate of phage. The LAS criteria of equilibrium without bacteria and phage $E_2$, the bacterial growth rate should be less than the innate immune killing rate. For the equilibrium without immune system $E_3$ to be LAS in the $B \rm{-}P$ plane, the carrying capacity of bacteria must greater than the threshold value, which depends on the decay rate, burst size, and adsorption rate of phage. For the equilibrium without phages $E_4$ to be LAS, the intrinsic growth rate of bacteria and the equilibrium density of the bacteria should be less than the threshold values. For the LAS of coexistence equilibrium $E_5$, the intrinsic growth rate must be greater than the threshold value. In Section 5, some numerical simulations are performed to verify the above theoretical results (see Theorem 3.1 and Figures 1, 2, and 4–7).

Transcritical bifurcations are directly related to the deletion or creation of a new equilibrium and its local stability nature. The system has undergone two possible transcritical bifurcations, which depend entirely on the threshold values of the carrying capacity of bacteria and the innate immune response. When the carrying capacity of bacteria ($K_C$) is less than the threshold $\omega /(\beta \phi)$, the equilibrium without phage and immune system $E_1$ is locally asymptotically stable. However, as soon as $K_C$ exceeds $\omega/(\beta \phi)$, then it leads not only $E_1$ as unstable but also this is the necessary criteria for existence of equilibrium without immune system $E_3$ and at $K_C = \omega/(\beta \phi)$, $E_3$ reduces to $E_1$. Therefore, transcritical bifurcation occurs at the threshold condition $K_C = \omega/(\beta \phi)$, around the equilibrium without phage and immune system $E_1$. Using the same argument as above, we can easily state that another transcritical bifurcation occurs at the equilibrium without bacteria and phage $E_2$ for $K_I = r/\epsilon$ and $K_I < r/\epsilon$, leading to the existence of equilibrium without phages $E_4$.

We analyze the global stability behavior for the equilibrium without immune system $E_3$ and that without phages $E_4$ by applying the Bendixson–Dulac criterion (see Theorems 4.1 and 4.2). The equilibrium without immune system $E_3$ and that without phages $E_4$ are globally asymptotically stable in the $B \rm{-}P$ and $B \rm{-}I$ planes, respectively, whenever they are locally asymptotically stable. In Theorem 4.3, we also establish the global asymptotic stability of coexistence equilibrium $E_5$ by applying the Lyapunov functional method. The role of the burst size of phage in the phage therapy model is crucial in determining the global stability behavior of the coexistence equilibrium.

We derive the sufficient conditions for the uniform persistence of system (1.1) (refer the Theorem 4.4). Biologically, if the birth rate of bacteria is greater than the immune killing rate and the carrying capacity of bacteria and the equilibrium density of the bacteria are greater than the threshold values, then the system is uniformly persistent. This is supported by some numerical examples in Figures 1, 2, and 7. In Theorem 4.6, we provide a certain condition for the extinction of the phage population. Biologically, this extinction criteria explains that if the carrying capacity of bacteria remains below the threshold value, which depends on the phage values, then the phage becomes extinct, and the system is non-persistent. This result reveals that a phage with a high value of $\omega/\beta \phi$ (i.e., a high decay rate, a low adsorption rate, and a small burst size) will become extinct. This finding is also supported by a numerical example (refer to Theorem 4.6 and Figure 4).

Numerically, the bacteria growth and innate immune killing rates affect the population of the species. Species $B$ and $P$ exist if the bacteria growth rate exceeds the innate immune killing rate; otherwise, they become extinct (see Figure 3).

Biologically, all species that co-exist exhibit an oscillatory balance behaviour. Meanwhile, a periodic solution arises in a system when the analyzed equilibrium point changes in stability as a function of its parameters. To capture the oscillating coexistence of populations, we establish the existence of Hopf bifurcation around coexistence equilibrium $E_5$ by considering the parameters in system (1.1) as a bifurcation parameter for future work. In addition, we will study model (1.1) with time delay to obtain a more realistic model. We will consider the influence of time delay on the stability of the system and the existence of a Hopf bifurcation solution. We will investigate the direction and stability of Hopf-bifurcating periodic solutions with the help of normal form theory and the center manifold theorem.

Acknowledgments

The authors would like to thank the reviewers and the editor for their careful reading, helpful comments and suggestions that greatly improved the paper.

Conflict of interest

The authors declare that they have no competing interests.