A mathematical model for Vibrio-phage interactions

1.   Introduction

Cholera is a gastrointestinal illness that presents with diarrhea and vomiting which, if left untreated, may result in death. Death occurs as a result of dehydration, metabolic acidosis and uremia ^[1,2]. A cholera infection is contracted upon consumption of water contaminated by the Vibrio cholerae bacterium. An infected individual may then release additional pathogenic bacteria (V. cholerae) into the environment by means of defecating or vomiting into or near a water source (this is referred to as shedding). As consuming contaminated water is the main route of transmission, areas without proper water sanitation are effected the most. These regions experience endemic cholera which appears with seasonality. Regions that have been significantly impacted by cholera include Haiti, Dominican Republic, Tanzania, Democratic Republic of the Congo, Somalia, Nigeria, South Sudan, Mozambique, Iraq and Afghanistan ^[3,4]. With 1.3 to 4 million cases per year, cholera remains a global public health concern ^[2].

While V. cholerae is the pathogen that creates a potentially devastating infection, it is not the only microbe at play. In the environment, V. cholerae interact with bacteriophage (which we will refer to as just phage). Phage are viruses that insert their genetic material into the bacterial cell. These phage may be either lytic or lysogenic. A lysogenic phage integrates it genetic material into the bacterial cell's genome. When the bacterium replicates by means of binary fission, the phage genetic material remains present in the bacterial daughter cells. The lysogenic cycle is a non destructive means for phage replication. Some lysogenic phage may be beneficial to the survival of a bacterial cell. In the case of cholera, the toxin production is attributed to the CTX$\varphi$ phage. A lytic phage inserts its genetic material into the bacterial cell and replicates within the host cell. This causes the cell to lyse (burst) resulting in cell death. The lytic process releases additional phage into the environment. While some phage aid in the survival of the bacterial cell, such as the phage that are responsible for the pathogenicity of $V. cholerae$, others do quite the opposite. The interaction of the bacteria and lytic phage resembles a predator-prey type relationship.

While V. cholerae may exist in biofilms which aid in the survival of the bacterium, another mechanism for cell survival is the ability of the bacterium to enter a viable but not culturable state (VBNC). A bacterial cell enters a VBNC state when conditions are harsh. These bacteria are viable to produce infection but do not replicate on routine media and so are not culturable. When more favorable conditions present, the bacterial cell may resuscitate ^[5]. Such favorable conditions may be the gut of an individual who has consumed water contaminated with bacteria in the VBNC state. Upon resuscitation, the bacteria replicate rapidly within the gut, resulting in a cholera infection.

While cholera has plagued the developing world since the 1800s, there has been a significant amount of research in the spread of the disease. Despite this, there has been only a few publications in modelling the impact of the bacteria and phage interactions. Since the pioneering work of Codeço in 2001 ^[6], mathematical modelling has been used to better understand cholera and how it might be controlled in regions that experience endemic cholera. The model developed by Codeço is an ordinary differential equation model consisting of a human and a bacterial population. The human population is divided amongst two explicit compartments (susceptible and infectious) and the bacteria are accounted for in one bacterial compartment (bacterial reservoir). The bacterial compartment included the human contribution of shedding new pathogenic bacteria into the reservoir. Codeço also explored seasonality of the cholera by the use of a periodic function for the contact and shedding terms. While the model Codeço developed included only three rather simple differential equations, it laid the foundation for all models to come.

After Codeço's work, many researchers have sought to understand the dynamics of cholera in a more complicated setting. As a large proportion of individuals infected with cholera produce only mild symptoms, King et al. incorporated asymptomatic infections in their model ^[7]. As shown in their paper, although there might be relatively few cholera infections reported in a host population, free-living bacteria in the environmental reservoir (aka, a reservoir environment) may be critical to the disease's endemicity. Models such as the one developed by Jensen et al. ^[8] take a different approach to study the environmental reservoir and explore the disease dynamics influenced by the bacteria and phage interaction. Specifically, the authors developed a model consisting of a human population, a bacterial population and a phage population. The human population was divided into four compartments; a susceptible, a recovered and two infectious compartments. The infectious compartments included those who were phage "positive" and those who were phage "negative." This model showed that the presence of phage effectively reduce the bacterial concentration and thus infections, which is consistent with the biological understanding of the bacteria-phage interactions. Kong, Davis and Wang ^[9] further incorporated an immunological threshold (a minimum dose of bacteria is needed to yield an infection) and showed that oscillating trajectories could exist in both the microbial and population scales. Notice that different functional responses were employed for vibrio-phage interaction: Holling type Ⅰ in ^[8] while Holling type Ⅱ in ^[9]. It remains unclear whether vibrio-phage interaction of Holling type Ⅰ could drive oscillation in cholear dynamics.

In this paper we investigate the cholera disease dynamics incorporating the bacteria-phage interaction of Holling type Ⅰ. Specifically, the model consists of a human population (SIRS), a bacteria (B) and phage population (P) to exhibit these interactions. We are able to show that no oscillating trajectories exist in a non-reservoir environment, while oscillating trajectories could exist in both the microbial and population scales in a reservoir environment. Our results demonstrate that a reservoir environment could be another driven force for cholera periodicity. In addtion, strategies for controlling cholera are discussed by conventional means (water purification, sanitary waste disposal, etc.) as well as unconventional means (presence of phage).

2.   Model development

We formulate a cholera model incorporating the interaction of bacteria and phage so that we could investigate their influence on the disease dynamics. Specifically the model couples the classic Susceptible-Infectious-Recovered-Susceptible (SIRS) model of the human disease with a predator-prey model of a Holling type Ⅰ functional response for lytic phage (P) and bacteria (B). Furthermore, we include only natural deaths (natural death rate $\mu$) and neglect disease related deaths for simplicity. This is justified as cholera is easily treated with either oral or intravenous fluids in conjunction with antibiotics (for severe cases), and thus cholera has a relatively low death rate. The probability of contracting an infection from consuming V. cholera contaminated water is given by the dose response $f(B) = \frac{B}{B+H}$. Here, $H$ is the bacterial concentration that yields a $50\%$ chance of infection. This dose response is the same function used in the Codeço model ^[6]. We assume that there are constant births $\Lambda$. The lytic phage interact with bacteria and result in bacterial death at a rate of $b$ while the phage have a gain from such bacterial death represented by $\chi$. The bacteria replicate in the environment and the bacterial "birth" rate is given by $\nu$, while the natural bacterial death rate is $\delta$. An environment is called as a reservoir environment if $\nu > \delta$ (bacteria grow in the absence of shedding from infectious individuals); otherwise, if $\delta > \nu$, then it is called as a non-reservoir environment (bacteria die out in the absence of the shedding from infectious individuals). Phage require a host cell (bacterium) to replicate and so there is no replication independent of bacterial cell death. The phage deactivate at a rate $m$. Note that we use the term deactivate instead of die, as a phage is a virus which does not fit the biological definition of a living organism. In contrast to the model developed by Kong et al. ^[9], we include a temporary disease induced immunity. This inclusion is important to consider as the immunity obtained from a cholera infection may range from a few weeks to greater than 3 years, depending on severity of infection ^[10,11,12]. This immunity is lost at a rate $\gamma$. With these considerations in mind, we have the following system of ordinary differential equations:

with non-negative initial conditions. A visual representation of the model is given by the model diagram in Figure 1. Note that our model can be further extended to include direct transmission ^[13] and the logistic growth of bacteria ^[14], which are ignored here, so that the impact of bacteria-phage interaction can be highlighted. The proposed model (2.1) can be regarded as a simplification of the model in ^[11]. This simplification allows for a more theoretical and deeper analysis to include global stability results which were not included in ^[11].

A list of parameters and their descriptions can be found in Table 1, all parameters are assumed to be positive, unless indicated otherwise. In addition, we include parameter values which are commonly used in the literature and also will be used for our simulations in Section 9. It is important to make a note regarding the parameter $\alpha$ before proceeding. In the literature, $\alpha$ takes on values smaller than one; for example, $\alpha = 0.1$ in ^[8]. However, $\alpha$ may be viewed as the product of the incubation rate and the number of contaminated consumptions per day. Given that the mean incubation for cholera is about 3 days, assuming an individual consumes a considerable amount of contaminated water 3 times per day, we set $\alpha = \frac{1}{3}\cdot 3 = 1$.

3.   Feasible region

Set $N = S+I+R$, and so $\dot{N} = \dot{S}+\dot{I}+\dot{R}$. It follows from the model (2.1) that

Now, using the method of the integrating factor, we have the equation $\frac{d}{dt}\{e^{\mu t}N\} = \Lambda e^{\mu t}$. Integrating this yields

where $N_0 = N(0).$ Letting $t\to \infty$, we see that $N(t)\to \frac{{\Lambda}}{{\mu}}$. If $N_0 > \frac{{\Lambda}}{{\mu}}$, we see that $\dot{N} < 0$ and so $N$ monotonically decreases from $N_0$ to $\frac{{\Lambda}}{{\mu}}$. Hence $N{\le}N_0$ in this case. On the other hand, if $N_0 < \frac{{\Lambda}}{{\mu}}$, $\dot{N} > 0$ and so $N$ monotonically increases from $N_0$ to $\frac{{\Lambda}}{{\mu}}$. So, in this case, $N{\le}\frac{{\Lambda}}{\mu}$. Lastly, if $N_0 = {\frac{{\Lambda}}{\mu}}$, $\dot{N} = 0$ and so $N = \frac{{\Lambda}}{\mu}$ for all time. From this, we conclude that for any initial population size $N_0$, $N$ is bounded. More precisely, $N{\le}\max\{N_0, \frac{{\Lambda}}{\mu}\}$ for all time.

Now, consider the case of a non-reservoir environment (i.e., $\delta > \nu$), and set $\epsilon = \min\{\delta-\nu, m\} > 0$. Let $E = \chi B + P$, and so $\dot{E} = \chi \dot{B}+\dot{P}$. From (2.1), we have

This gives $\dot{E}+\epsilon E \le \eta \chi\frac{\Lambda}{\mu}.$ Again using the method of the integrating factor, we obtain $\frac{d}{dt}\{Ee^{\epsilon t}\}\le \eta \chi \frac{\Lambda}{\mu}e^{\epsilon t}$, which upon integrating gives $E(t) \le \frac{\eta{\chi}\Lambda}{{\epsilon}\mu}+Ce^{-\epsilon{t}},$ where $C$ is a constant. Thus, $\limsup_{t\to\infty} E(t)\le \frac{\eta{\chi}\Lambda}{{\epsilon}\mu}$. This motivates the following lemma.

Lemma 3.1. Given the condition $\delta > \nu,$ we define the feasible region by

which is positively invariant with respect to (2.1).

For the case $\delta > \nu$, dynamics of (2.1) within $\Gamma$ will be carried out analytically in Sections 4 and 5.

For the case of a reservoir environment (i.e., $\nu > \delta$), bacteria can grow in the absence of cholera infection, and dynamics of (2.1) could be more complicated. For example, its solutions with initial value $P(0) = 0$ could be unbounded, and numeral simulations in Section 9 also show the existence of sustained oscillations.

4.   Disease free equilibrium and stability

Letting $I = 0$ and setting $\dot{S} = \dot{I} = \dot{R} = \dot{B} = \dot{P} = 0$ in (2.1), it is clear that $Q^0 = (\frac{\Lambda}{\mu}, 0, 0, 0, 0)$ is the disease free equilibrium (DFE). This equilibrium always exists and lies on the boundary of $\Gamma$. We now discuss both local stability as well as global stability of the DFE.

To explore the local stability of the DFE, we will examine the Jacobian matrix. However, before proceeding we exchange the $I$ and $R$ equations to obtain a matrix with a nice structure. Note that this arrangement does not effect the component order of the DFE as both the $I$ and $R$ components at DFE are zero. The Jacobian matrix at the DFE is then given by

This matrix is block upper triangular and so three eigenvalues are $-\mu, -(\gamma+\mu), -m < 0$ and the remaining two eigenvalues are given by the eigenvalues of the $2\times 2$ matrix

Now, we may conclude that the remaining two eigenvalues have negative real part provided $\mathrm{Tr}(M) < 0$ and $\mathrm{Det}(M) > 0$. We have $\mathrm{Det}(M) = -(r+\mu)(\nu-\delta)-\frac{{\Lambda}{\alpha}{\eta}}{H{\mu}} > 0$ implies $\frac{\Lambda{\alpha}\eta}{H{\mu}\delta(r+\mu)}+\frac{\nu}{\delta} < 1$. For this inequality to hold, it is required that $\delta > \nu$. Now, $\mathrm{Tr}(M) = -(r+\mu)+\nu-\delta < 0$. Thus, under the condition $\frac{\Lambda{\alpha}\eta}{H{\mu}\delta(r+\mu)}+\frac{\nu}{\delta} < 1$, the DFE is locally stable and unstable otherwise. Now we define the following threshold parameter

Biologically, $\mathcal{R}_B$ describes the average number of next generation of bacteria, and accounts of both reproductions due to bacteria shedding from infectious individuals ($\frac{\Lambda{\alpha}\eta}{H{\mu}\delta(r+\mu)}$) and due to its natural growth ($\frac{\nu}{\delta}$). The following result shows that $\mathcal{R}_B$ serves as a threshold value determining whether the disease dies out (when $\mathcal{R}_B < 1$). The relation between $\mathcal{R}_B$ and the basic reproduction number $\mathcal{R}_0$ (a more commonly used threshold value in the literature) will be discussed in Section 6.

Theorem 4.1. The disease free equilibrium is globally asymptotically stable if $\mathcal{R}_B < 1$, and is unstable if $\mathcal{R}_B > 1$.

Proof. It suffices to show that the disease free equilibrium is globally asymptotically stable when $\mathcal{R}_B < 1$. Consider the function

Clearly $L \ge 0$. Now, differentiating with respect to (2.1), we have

So, it is clear that when $\mathcal{R}_B < 1$, $\dot{L} \le 0$. We have that $L$ is a Lyapunov function. Now, $\dot{L} = 0$ implies $B = P = 0$. Since $B = P = 0$, we have $\dot{B} = 0$ and so from the equation for $\dot{B}$ in (2.1), we see that $I = 0$. Now, since $I = 0$, $\dot{I} = 0$ and so by the equation for $\dot{I}$ in (2.1), we may conclude $S = \frac{\Lambda}{\mu}$. Since $S = \frac{\lambda}{\mu}$, $\dot{S} = 0$ and so the equation for $\dot{S}$ in (2.1) gives that $R = 0$. So, the largest invariant set such that $\dot{L} = 0$ is given by $\{Q^0\}$. Thus, by LaSalle's invariance principle ^[15], we conclude that $Q^0$ is globally asymptotically stable in $\Gamma$.

Biologically, Theorem 4.1 implies that when the critical value $\mathcal{R}_B < 1$, the disease fails to persist for any initial conditions within the feasible region. Alternatively, if $\mathcal{R}_B > 1$, the disease free equilibrium is unstable. Simulations with different initial conditions have been performed, as depicted in Figure 2. The biological meaning and importance of the critical value $\mathcal{R}_B$ is discussed further in section 6.

5.   Endemic equilibria and their stability

In this section, we provide detailed analysis for disease dynamics when $\mathcal{R}_B > 1$. First, model (2.1) admits potentially two distinct endemic equilibria. One of the equilibria, when it exists, exists in the absence of phage. This equilibrium will be referred to as the phage free endemic equilibrium (PFEE) and will be denoted by $Q_{PFEE} = (S_*, I_*, R_*, B_*, 0)$. Another endemic equilibrium, when it exists, exists in the presence of phage and so this equilibrium is referred to as the phage present endemic equilibrium (PPEE) and will be denoted by $Q_{PPEE} = (S^*, I^*, R^*, B^*, P^*).$

5.1.   Phage free endemic equilibrium and stability

For the PFEE, set $P = 0$, and $\dot{S} = \dot{I} = \dot{R} = \dot{B} = \dot{P} = 0$. It follows from (2.1) that

So, (5.3) gives that $R_* = \frac{r}{\gamma+\mu}I$, (5.4) gives $B_* = \frac{\eta}{\delta-\nu}I_*$ and so it is necessary that $\delta > \nu$, and (5.2) gives $S_* = \frac{r+\mu}{\alpha}\big(\frac{\delta-\nu}{\eta}H+I_*\big)$. We now must determine when $I_*$ is feasible. We use $S_*, R_*, B_*$ in (5.1) and obtain

For $I_*$ to be positive the left hand side of the above equation must be positive. That is, we must have ${\Lambda}\alpha{\eta}-{\mu}(r+\mu)(\delta-\nu)H > 0$ which occurs when $\mathcal{R}_B > 1$. So, under the condition $\mathcal{R}_B > 1$ and $\delta > \nu$, the phage free endemic equilibrium exists, denoted as $Q_{PFEE} = (S_*, I_*, R_*, B_*, 0)$.

Before proceeding, we define another threshold parameter

Biologically, $\mathcal{R}_P$ can be regarded as the phage invasion reproduction number, characterizing whether the phage population can succeed to invade (see Theorems 5.1 and 5.2). It is straight forward to show $\mathcal{R}_P < \mathcal{R}_B$.

Now, to explore the local stability of the PFEE, we again turn to the Jacobian matrix. We have

One of the eigenvalues is $\lambda_5 = \chi b B_*-m$ while the remaining eigenvalues are not readily apparent. This is again a block upper triangular matrix and so the remaining 4 eigenvalues are given by the eigenvalues of

We compute the characteristic polynomial of $C$ and obtain the polynomial

It can be show that $t = -\mu$ is a root of the polynomial and so it remains to show the roots of the third degree factor have negative real parts. To do this, we utilize the Routh-Hurwitz criterion which states that a third degree polynomial of the form $p(t) = t^3+a_1t^2+a_2t+a_3$ has roots that lie in the left half of the complex plane provided $a_1, a_2, a_3 > 0$ and $a_1a_2 > a_3$. Indeed, it is clear that for $\delta > \nu$, $a_1, a_2, a_3 > 0$. Now,

So, by the Routh-Hurwitz criterion, we may conclude that when $Q_{PFEE}$ exists, $\lambda_1, \lambda_2, \lambda_3, \lambda_4$ have negative real parts. All that is left to conclude local stability of the PFEE is determining when the eigenvalue $\lambda_5 = \chi b B_*-m < 0$. Using $B_*$ in terms of $I_*$ and writing $I_*$ out explicitly, we obtain the the inequality

which holds if and only if $\mathcal{R}_P < 1$. Thus we can conclude that $Q_{PFEE}$ is locally stable provided $\mathcal{R}_P < 1 < \mathcal{R}_B$. Note that the condition $\delta > \nu$ is implied by $\mathcal{R}_P < 1$.

We now explore the global stability of the PFEE when $\gamma = 0$. That is, when recovered individuals attain permanent immunity.

Theorem 5.1. If $\mathcal{R}_P < 1 < \mathcal{R}_B$ and $\gamma = 0$, then $Q_{PFEE}$ is globally asymptotically stable.

Proof. To show global stability of the PFEE, consider the function

It is clear that $L \ge 0$. Taking the derivative with respect to $t$, we have

As discussed previously, ${\chi}bB_*-m < 0$ provided $\mathcal{R}_P < 1$. It remains to show negativity of the remaining terms in the sum. Indeed,

implies $2-\frac{S_*}{S}-\frac{S}{S_*} \le 0$. We now write

Now,

implies $4-\dfrac{S_*}{S}-\dfrac{f(B_*)B}{f(B)B_*}-\dfrac{f(B)SI_*}{f(B_*)S_*I}-\dfrac{B_*I}{BI_*} \le 0$. We also have

That is, $\dot{L} \le 0$ and so $L$ is a Lyapunov function. It is clear that $\dot{L} = 0$ implies $S = S_*, I = I_*, B = B_*$ and $P = 0$. So, the largest invariant set such that $\dot{L} = 0$ is the set $\{Q_{PFEE}\}.$ Thus, it follows from LaSalle's invariance principle ^[15] that $Q_{PFEE}$ is globally asymptotically stable in $\Gamma^0$.

Theorem 5.1 implies that when the critical values $\mathcal{R}_P < 1 < \mathcal{R}_B$, the phage free endemic equilibrium is globally asymptotically stable provided $\gamma = 0$. Biologically, the infectious population persists at the endemic level $I_*$, irregardless of the initial conditions. Simulations have shown this also holds for $\gamma > 0$, as depicted in Figure 3.

5.2.   Phage present endemic equilibrium and stability

We now examine an equilibrium in which phage persist. To do this, we examine the equilibrium equations

with $P\ne 0$. Equation (5.9) implies $B^* = \frac{m}{\chi b}$. This and Eq (5.8) give $P^* = \frac{(\nu-\delta)}{b}+\frac{{\eta}\chi}{m}I^*$ which is feasible and distinct from the PFEE when $\frac{(\nu-\delta)}{b}+\frac{{\eta}\chi}{m}I^* > 0$. From Eq (5.7), it is clear that $R^* = \frac{r}{\gamma+\mu}I^*$ while it is clear from Eq (5.6) that $S^* = \frac{r+\mu}{\alpha}\big(\frac{{\chi}bH}{m}+1\big)I^*$. Using this in Eq (5.5), we obtain

So, the only condition for existence of the PPEE is $\frac{(\nu-\delta)}{b}+\frac{{\eta}\chi}{m}I^* > 0$. It can be shown that this inequality holds if and only if $\mathcal{R}_P > 1$. It is shown in Section 8 that $I^*\le I_*$, where $I_*$ is the infectious component of the PFEE as discussed in the previous section.

The following theorem establishes global stability of the PPEE when recovered individuals have permanent immunity ($\gamma = 0$).

Theorem 5.2. Suppose $\delta > \nu$ and $\gamma = 0$. If $\mathcal{R}_P > 1$, then the phage present endemic equilibrium is globally asymptotically stable.

Proof. We again make use of a Lyapunov function. Consider the function

Differentiating with respect to $t$ yields

From which, it is clear that $\dot{L}\le 0$. So, $L$ is indeed a Lyapunov function. Now, $\dot{L} = 0$ implies $S = S^*, I = I^*, B = B^*$. Moreover, since $B = B^*$, we have $\dot{B} = 0$ and so by (5.9), we have that $P = P^*$. So the largest invariant set such that $\dot{L} = 0$ is $\{Q_{PPEE}\}$. Now, by LaSalle's invariance principle, we have that $Q_{PPEE}$ is globally attracting $\Gamma^0$. Moreover, we have that $Q_{PPEE}$ is locally stable by Lyapunov's stability theorem. These two results give that the phage present endemic equilibrium $Q_{PPEE}$ is globally asymptotically stable. in $\Gamma^0$.

Recall from Sections 4 and 5.1, the condition $\mathcal{R}_P > 1$ implies $\mathcal{R}_B > 1$. Moreover, the condition $\nu > \delta$ ensures the existence of the PFEE. We now observe two feasible endemic equilibrium and one disease free equilibrium. Theorem 5.2 states that when these conditions are satisfied, the DFE and PFEE are unstable while the PPEE is globally asymptotically stable provided $\gamma = 0$. Biologically speaking, the infection remains endemic, with the presence of phage, at the level $I^*$. Our simulations confirm this and also show this still holds for $\gamma > 0$, as depicted in Figure 4. This endemic equilibrium (PPEE) gives insight into how the phage may aid in understanding and controlling cholera infection as discussed in Section 8.

6.   Reproduction numbers

In this section we revisit the next-generation matrix method ^[16,17], which has been widely used in the literature of mathematical epidemiology, and derive various reproduction numbers for our model (2.1). The relations among these reproduction numbers and $\mathcal{R}_B$ defined in (4.2) are discussed.

We examine the infectious sub-system of (2.1), which is given by

and its linearization at the DFE

As shown in ^[14,18,19], the decomposition of $M$ into a transmission matrix $F$ and a transfer matrix $V$ (i.e., $M = F-V$) might not be unique, yielding various threshold values (aka various reproduction numbers) that all determine the stability of the DFE. We choose the "largest" matrix $F$ as possible, since all other threshold values can be regarded as the corresponding target reproduction numbers as shown in Theorem 6 in ^[18]. Hence, set

and the next-generation matrix is

The spectral radius of $K$ provides a threshold value

Now we apply the next generation matrix $K$ and the corresponding weighted digraph $\mathcal{G}$ as depicted in Figure 5 to derive other threshold values which are of more biological interests.

First, $\mathcal{R}_B$ can be regarded as the average number of the secondary bacteria of a single bacterium. Specifically, the secondary bacteria can be due to either its self reproduction (the loop $B \to B$ in $\mathcal{G}$) or the bacteria shedding ($I \to B$) from the infectious individuals caused by the bacterium ($B \to I$). Thus, the number of the secondary bacteria becomes the sum of the weight of the loop $B\to B$ and the weight of the closed walk $B\to I \to B$, that is, $\mathcal{R}_B = \frac{\nu}{\delta} + \frac{\alpha \Lambda}{\mu \delta H} \cdot \frac{\eta}{r+\mu}$. Notice that the weight of a closed walk is the product of weights of all arcs in the walk.

Next, the basic reproduction number $\mathcal{R}_0$, the average number of the secondary infections while introducing a single infection in a completely susceptible host population, can be evaluated as the sum of weights of all closed walks starting at $I$ and ending at $I$ in $\mathcal{G}$ ($I\to B \to I$, $I \to B \to B \to I$, $I \to B \to B \to B \to I$, $\cdots$), that is,

In a non-reservoir environment (i.e., $\delta > \nu$), the geometric series in (6.4) converges, which leads to

In a reservoir environment (i.e., $\delta < \nu$), the geometric series in (6.4) diverges, implying that $\mathcal{R}_0$ approaches to infinity and expression (6.5) becomes invalid. Biologically, this means that the new infections can be produced without introducing any single infection. Notice that the threshold value $\mathcal{R}_K$ defined in (6.2) has been customarily called as the "basic reproduction number" in the literature. However, here we follow the biological interpretation of the basic reproduction number, so prefer it as $\mathcal{R}_0$ defined in (6.4) (in a non-reservoir environment, which becomes (6.5)).

We will show in the next section that $\mathcal{R}_B$ and $\mathcal{R}_0$ are target reproduction numbers of $K$ with specific disease control strategies, and serve as threshold values (staying at the same side of 1). In specific, the following result holds and follows immediately from Theorems 6 and 7 in ^[18].

Lemma 6.1. (i) Either $\mathcal{R}_B > \mathcal{R}_K > 1$ or $\mathcal{R}_B = \mathcal{R}_K = 1$ or $1 > \mathcal{R}_K > \mathcal{R}_B$.

(ii) Suppose that $\delta > \nu$. Then either $\mathcal{R}_0 > \mathcal{R}_B > \mathcal{R}_K > 1$ or $\mathcal{R}_0 = \mathcal{R}_B = \mathcal{R}_K = 1$ or $1 > \mathcal{R}_K > \mathcal{R}_B > \mathcal{R}_0$.

7.   Target reproduction number and control strategies

The target reproduction number provides a measure of effort required to control the disease by targeting only certain types of infectious interactions. That is, the target reproduction targets specific entries of the next generation matrix where the $(i, j)$ entry of $K$ represents the effect on infections that infectious type $j$ has on infectious type $i$. When targetting multiply entries of $K$, we denote the target set by $S$ where $S$ consists of all entries of $K$ that are targeted. Here, the notation and methods follows Shuai, Heesterbeek and van den Driessche ^[21], while the readers are referred to ^[18] for the recent update. For a target set $S$, the target reproduction number is given by $\mathcal{T}_S = \rho(K_S(I-K+K_S)^{-1})$ where $[K_S]_{i, j} = [K]_{i, j}$ if $(i, j){\in}S$ and $0$ otherwise. If targeting just one entry of $K$, we write the targets set as the entry being targeted. That is, if we are targeting only the $(i, j)$ entry, we write the target reproduction number as $\mathcal{T}_{i, j}$.

We wish to only consider feasible strategies for disease control. These strategies are regularly implemented in areas with endemic cholera. The control strategies include reducing the consumption of contaminated water and reducing the shedding of bacteria into the environment.

Reducing the consumption of contaminated water considers the effect on the disease that the bacteria have on the infectious population. Examples of strategies that reduce consumption of contaminated water include water filtration devices, chemical water treatment, and improving water sanitation infrastructure. Control strategies of this type consider the $(1, 2)$ entry of $K$ defined as in (6.1). The target reproduction number is then given by $\mathcal{T}_{1, 2} = \frac{{\Lambda}{\alpha}{\eta}}{H{\mu}(r+\mu)(\delta-\nu)} = \mathcal{R}_0$.

Reducing the amount of bacteria shed into the environment considers the effect on the disease that infectious human population has on the bacterial population. Examples of control strategies that reduce bacterial shedding include building proper latrines or other sanitary methods of human waste disposal. Control strategies of this type target the $(2, 1)$ entry of $K$. The target reproduction number is then given by $\mathcal{T}_{1, 2} = \mathcal{T}_{2, 1}.$

Both strategies for control may be used simultaneously and with equal effort. In this case, we consider the target set $S = \{(1, 2), (2, 1)\}$. In this case, $\mathcal{T}_S = \sqrt{\mathcal{T}_{1, 2}}$. So, all of the target reproduction numbers provide an equivalent threshold.

We did not consider the target reproduction numbers $\mathcal{T}_{1, 1}$ or $\mathcal{T}_{2, 2}$. The target reproduction number $\mathcal{T}_{2, 2}$ targets the $(2, 2)$ entry of $K$. That is, it considers the effect that bacteria have on bacteria (bacterial growth rate) which cannot feasibly be controlled. The target reproduction number $\mathcal{T}_{1, 1}$ targets the $(1, 1)$ entry of $K$. That is, it considers the effect that the infectious human population has on the infectious human population. Since our model does not consider human to human transmissions, we do not consider strategies for control of this type.

8.   Impact of phage on endemic equilibria

We now wish to explore how the presence of phage impacts the endemic equilibria of the model: the phage free endemic equilibrium $(S_*, I_*, R_*, B_*, 0)$ and the phage present endemic equilibrium $(S^*, I^*, R^*, B^*, P^*)$. We will show that the infectious levels may be reduced to any desirable degree with the appropriate phage growth. We do this by examining the impact environmental environmental factors have on disease dynamics.

Recall that the infectious component of the PPEE is given by

Now, considering $I^*$ as a function of $m$, we see that $I^*$ is strictly increasing in $m$. Also recalling the definition of $\mathcal{R}_P$, we see that $\mathcal{R}_P$ is strictly decreasing in $m$.

For a moment, suppose that $\mathcal{R}_B > 1$ and $\delta > \nu$. That is, the PFEE exists. Recall the infectious component of the PFEE is given by

By choosing

we see that $\mathcal{R}_P(m_0) = 1$ and $I^*(m_0) = I_*$. Here, $m_0 > 0$ is guaranteed by the conditions $\mathcal{R}_B > 1$ and $\delta > \nu$.

Now, since $\mathcal{R}_P(m)$ is strictly decreasing in $m$, we have that $\mathcal{R}_P(m) > 1$ for $0 < m < m_0$. This gives that $I^*$ exists uniquely for $m < m_0$ and since $I^*(m)$ is strictly increasing in $m$, we have that $I^*(m) < I_*$. This motivates the following theorem.

Theorem 8.1. If $\mathcal{R}_B > 1$ and $\delta > \nu$, then for any $\epsilon > 0$, there is an $0 < m_{\epsilon} < m_0$ such that $I^*(m) < \epsilon$ whenever $0 < m < m_{\epsilon}$.

This theorem show that the number of infectious individuals can be reduced by increasing the phage growth rate. Biologically, this means decreasing the phage deactivation rate $m$. While, on the surface, this may not seem as a feasible strategy for controlling the disease, this may be accomplished by an addition to the phage compartment proportional to the current phage population. Consider an addition of $kP$ to the phage growth rate where $k > 0$ is a constant. Then the phage growth rate is $\dot{P} = \chi b BP-(m-k)P$ and so $m-k < m$ effectively "reducing" the phage deactivation rate.

Recalling that $\mathcal{R}_B > 1$ implies $\mathcal{R}_0 > 1$, we see that the disease remains endemic. Despite the disease remaining endemic, the number of infectious individuals can effectively be reduced in the presence of phage. The ultimate goal with any illness should be complete eradication. However, in regions where cholera remains endemic, phage may aid in reducing the number of infectious individual to a more manageable amount.

9.   Simulations

Before we begin the discussion of numerical results, we must first begin by discussing reasonable initial conditions for the model. Biologically, bacteria and phage can exist in the absence of human shedding. However, by examining the model in (2.1), we see that $\dot{P} = 0$ for $t\ge 0$ provided $P(0) = 0$. That is, an introduction of a single infectious individual to a setting absent of bacteria and phage will result in a the equation $\dot{P}$ "dropping out" of the model. While realistically, an infectious individual may be able to introduce bacteria and phage by means of shedding. This reality is not considered in this model and so to remedy this, we require the initial condition $P(0) > 0$. This requirement may be justified by the fact that in areas with previous outbreaks, both bacteria and phage may exist in the environment despite the absence of endemic cholera. We proceed with the two cases: a non-reservoir environment ($\delta > \nu$), and a reservoir environment ($\nu > \delta$).

9.1.   Non-reservoir environment ($\delta > \nu$)

To explore the disease dynamics and how they are effected by the vibrio-phage interactions when $\gamma\ne 0$, we turn to numerical simulations. After an individual recovers from a cholera infection, they retain immunity for a period of time. Though this serves as protection against infection from V. cholerae, it does not last forever. As discussed in Section 2, the immunity period after a cholera infection ranges from just a few weeks to more than 3 years. For the sake of numerical simulations, we let $\gamma = \frac{1}{365}$. That is, we choose an immunity period of 1 year.

We vary the phage deactivation rate $m$ and see how this impacts the infectious population. As observed in the previous analytic results, we anticipate that decreasing $m$ will reduce the infectious levels. Recall that the threshold number $\mathcal{R}_B$ and the infectious equilibrium $I_*$ are independent of the phage deactivation rate $m$ and so they are constant given parameters chosen from the ranges given in Table 1. We set $\nu = 0, \delta = 0.33, \chi = 100, \eta = 10, b = 0.00002, r = \frac{1}{5}, \alpha = 1$. So, we have $\mathcal{R}_B = 1.5149$ and $I_* = 46.4483$. Using these chosen parameter values, we first simulate the infectious population for $2000$ days with $m = 7.9$. Now, using the formula for $\mathcal{R}_P$, we have $\mathcal{R}_P = 0.6196 < 1 < \mathcal{R}_B$. We see that the trajectory for the infectious compartment settles to $I_*$. We see similar behavior for the infectious compartments when $m = 4$. Here, $\mathcal{R}_P = 0.8748 < 1 < \mathcal{R}_B$. We again see that the infectious population eventually settles to $I_*$. Again decreasing $m$, this time $m = 1$, we now observe that $1 < \mathcal{R}_P = 1.2086 < \mathcal{R}_B$ and $I^* = 21.1305$. In contrast to the previous values for $m$, we see that the infectious population settles to $I^* < I_*$ and phage population becomes persistent at the phage-present endemic equilibrium, as depicted in Figure 6. It is also seen in Figure 6(a) that the decreasing the parameter $m$, reduces the infectious peak through the simulations. That is, an environment that is favorable to phage survival may help reduce the infectious peak of an outbreak, thus preventing the healthcare system from becoming overburdened. Additionally, we observe from Figure 6(b) that the transient behavior (peak) of the phage population appears not to be monotone, while the phage component of the stable equilibrium (PPEE or PFEE) appears to be monotone decreasing with respect to $m$.

Now, we wish to numerically examine the impact a decrease in the parameter $m$ has on the disease dynamics. Recall from Section 8 that there exists a value $m_0$ such that $I^*(m_0) = I_*$. Moreover, recall that the infectious level may be reduced to any desired amount by selecting an appropriate $m < m_0$. Using the parameter choices given in Table 1, we consider $I^*$ as a function of the phage deactivation rate $m$, as depicted in Figure 7. We see that while the infectious component of the PFEE remains constant, the infectious component of the PPEE approaches zero as $m$ approaches zero. Specifically, when $m > m_0\approx 2.815$, $\mathcal{R}_P$ becomes to be less than 1 and thus the phage present endemic equilibrium vanishes.

9.2.   Reservoir environment ($\nu > \delta$)

In this section we explore the case $\nu > \delta$, that is, the net rate of bacterial growth is positive in the absence of shedding. As a consequence, the bacterial population could exist in the environment in the absence of the human contribution of shedding. As shown previously, if $\nu > \delta$, then $\mathcal{R}_P\ge 1$ (and thus $\mathcal{R}_B > 1$ as well), the DFE is always unstable (Section 4), the PFEE does not exist (Section 5.1), and the PPEE always exists (Section 5.2). Numerical simulations have shown that the model (2.1) could admit complex dynamics. Specifically, with parameter values given in Table 1 and setting $\nu = 2$, $\delta = 1$, $m = 7.9$ and $b = 1.4e-9$, simulations, as depicted in Figure 8, show the coexistence of two stable periodic solutions. Specifically, two solutions with different initial conditions are displayed in the $B$-$P$-$I$ phase space, while only trajectories for large time $t$ of these two solutions are shown in the $B$-$P$ phage plane and in the $I$-$t$ plot.

9.3.   Sensitivity analysis

The objective of this section is to discuss the sensitivity of $\mathcal{R}_0$, $\mathcal{R}_B$, $\mathcal{R}_P$, outbreak peak value and time, and the total number of infected individuals to model parameters in (2.1). To this end, we employ the Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) method ^[20,27]. The selected parameters are those that can be used to inform policies such as: rate of bacterial consumption ($\alpha$), recovery rate ($r$), bacterial shedding rate ($\eta$), natural bacteria growth rate ($\nu$), bacteria death rate ($\delta$) and phage attack rate ($b$). We generate 3000 samples of the model parameters, using LHS and varying them between 75% and 125% of their estimated values. We then verify the monotonic relationships between the parameters and the model outcomes using the PRCC. The sign of PRCC provides a measure of the nature of the linear association and its magnitude provide a measure of the strength of the linear association. It varies between $-1$ and $1$. A relation between model output and a given parameter value is consider significant if the magnitude of the PRCC is greater than 0.5. Panels A, B and C of Figure 9 show that the bacterial death rate $\delta$ and natural growth rate $\nu$ have the strongest relationship to $\mathcal{R}_B$, $\mathcal{R}_P$ and $\mathcal{R}_0$. In contrast to the shedding rate $\eta$ which has among the lowest of PRCC, $\delta$ and $\nu$ would thus be an important parameter to control in order to reduce the harm of an outbreak. While $\mathcal{R}_0$ and $\mathcal{R}_P$ are only sensitive to $\delta$ and $\nu$, $\mathcal{R}_B$ on the other hand is significantly sensitive to $\alpha$, the rate at which humans are consuming water from the reservoirs. The total number of infected persons at the end of a year is sensitive to $\gamma$, $\alpha$, $r$, $\nu$, $\eta$, $b$, and $m$ with $r$ having one of the strongest relation (Panel D, Figure 9). The magnitude of the outbreak peak is more sensitive to $\alpha$ and $r$ (Panel E, Figure 9). The peak time like the peak magnitude is more sensitive to $\alpha$, but unlike the peak magnitude, it is as well sensitive to $H$ and $\nu$ (Panel F, Figure 9). These results suggest that control measures influencing the bacterial consumption and bacterial growth rate will be more effective in minimizing the epidemic than those concentrating on influencing the shedding rate. While improving the sanitation infrastructure of an area is the obvious step to take to control outbreaks, monitoring and controlling the bacterial levels in the reservoir itself is equally important. Improving the infrastructure would surely help control the bacterial levels in the reservoir by decreasing the amount of human contamination, but V. cholerae exist independently of humans and so other factors that influence the natural levels of bacteria in the water need to be considered as well in intervention strategies. Controlling both parameters is important and likely to be the most effective, but the bacterial consumption rate and bacterial growth rate are the more influential of the three.

10.   Human and bacterial contribution to threshold values

We now wish to examine how human behavior and natural bacterial growth influence the threshold values $\mathcal{R}_B$ and $\mathcal{R}_P$. We explore this influence numerically by means of examining heatmaps generated for $\mathcal{R}_B$ and $\mathcal{R}_P$ considering parameters that describe human contribution or involvement and natural bacterial growth. Such parameters include the consumption rate of contaminated water, $\alpha$; the rate in which humans shed bacteria, $\eta$; and the natural bacterial growth rate, $\nu$. We examine for what values these parameters give threshold values greater than one. In doing this, we make biological conclusions on how important it is to address all aspects of human behavior pertaining to shedding and consumption. This emphasizes a holistic approach to combating cholera epidemics.

We examine the heatmap in Figure 10. Here, we see that even for a small consumption rate, if the bacterial shedding is large enough, both $\mathcal{R}_B$ and $\mathcal{R}_P$ exceed the threshold value 1. Similarly, for a small rate in which humans shed bacteria, if the consumption rate is large, we again have $\mathcal{R}_B > 1$ and $\mathcal{R}_P > 1$. With this in mind, we observe that both shedding and consumption must be addressed if a disease free state is to be attained.

Now, examining the heatmaps given in Figure 11, we see that for any given bacterial replication rate $\nu$, increasing the shedding $\eta$, the threshold values $\mathcal{R}_B$ and $\mathcal{R}_P$ will cross one. We also note that we used $\delta = 0.33$. It is seen that when $\eta = 0$, as $\nu$ crosses $0.33$, the threshold values cross one. This is due to the term $\frac{\nu}{\delta}$ seen in the threshold value definitions.

11.   Discussion

A new cholera model has been developed in this paper to incorporate vibrio-phage interaction of Holling type Ⅰ, which is similar to the one in Jensen, et al. ^[8] but different as the one in Kong et al. ^[9]. Our theoretical and numerical results highlight the importance of a reservoir vs non-reservoir environment, as the former might be a new mechanism to drive cholera periodicity. As a consequence, it is of practical importance to understand the environmental growth of V. cholerae.

Our result supports the idea that phage can effectively reduce the number of infectious individuals. As our model assumes bacteria grow exponentially in the absence of phage and bacterial shedding, it is of large interests to incorporate logistic growth ^[14,22]. Further studies are also necessary to incorporate other important epidemiological features, such as the direct human-to-human transmission (i.e., a fast pathway or a short cycle of transmission) ^{[13,23,24,26]}, hyperinfectivity ^[25], and the immunological threshold (a minimum dose of bacteria is required to yield an infection) ^[9]. Infectious individuals shed not only bacteria to the environment, also a certain amount of phage. The resulting model will lose the phage-free endemic equilibrium, and more quantitative studies, incorporating real biological and epidemiological data, are needed to investigate the impact of vibrio-phage interaction on disease dynamics. Last but not least, as discussed in the Introduction, an important piece of the biological puzzle is the ability for the bacteria to enter a viable but not culturable state (VBNC). A future work on mathematical models incorporating VBNC could provide us a better understanding of cholera transmission and control.

Acknowledgments

This research was partially supported by the Natural Science and Engineering Research Council of Canada (NSERC) through Discovery Grant RGPIN-2020-03911 (HW), Accelerator Grant RGPAS-2020-00090 (HW), and by the National Science Foundation (NSF) through the grant DMS 1716445 (ZS). C. Botelho acknowledges that this work is a continuation of his M.S. thesis completed at the University of Central Florida and also acknowledges the support of the ORC Doctoral Fellowship. The authors thank the editor and reviewers for their helpful feedback.

Conflict of interest

The authors declare there is no conflict of interest.