Bacteria–bacteriophage cycles facilitate Cholera outbreak cycles: an indirect Susceptible-Infected-Recovered-Bacteria- Phage (iSIRBP) model-based mathematical study

Cholera is an acute enteric infectious disease caused by the Gram-negative bacterium Vibrio cholerae. Despite a huge body of research, the precise nature of its transmission dynamics has yet to be fully elucidated. Mathematical models can be useful to better understand how an infectious agent can spread and be properly controlled. We develop a compartmental model describing a human population, a bacterial population as well as a phage population. We show that there might be eight equilibrium points, one of which is a disease free equilibrium point. We carry out numerical simulations and sensitivity analyses and we show that the presence of phage can reduce the number of infectious individuals. Moreover, we discuss the main implications in terms of public health management and control strategies.


Introduction
Cholera is an acute enteric infectious disease caused by the Gram-negative bacterium Vibrio cholerae (V. cholerae), which spreads by water or food contaminated with feces. It causes diarrhea, which can lead to dehydration that might lead to death within hours if it untreated [7]. Its effect is more dangerous among individuals who have insufficient access to safe water and proper sanitation; it is also more dramatic in areas where basic environmental infrastructures are disrupted or have been destroyed. The World Health Organization (WHO) estimated that there are from 1.3 to 4 million cholera cases and 21,000 to 143,000 deaths due to cholera yearly [19]. However, in countries with adequate health care, cholera imposes a much less dramatic burden.
The precise nature of the cholera transmission dynamics is not clear yet, despite a huge body of research and the many efforts of the scientific community. For example, [4] discussed the response of men to infection with V. cholerae. In [17], Snow model of cholera was proposed. [9] discussed the role of blue-green algae in maintaining endemicity and seasonality of cholera in Bangladesh. [6] proposed that a viable but nonculturable V. cholerae O1 strain can revert to a cultivable state in the human intestine [6], and in [8], the authors discussed a critical element in the ability of V. cholerae to cause outbreaks and epidemics. It is of paramount importance to understand how V. cholerae can transmit through men. This would enable researchers and public health authorities to monitor the evolution of the outbreak and help put into effect adequate intervention measures that limit the effects of cholera spreading.
Mathematical models of infection transmission dynamics are extremely useful in better understanding how an infectious agent can spread and how can be effectively managed and controlled. More specifically, mathematical modelling of cholera has a long history; through the decades, many models have been devised and proposed. However, few of them have incorporated indirect transmission of the disease. Among these, the major models that employ ordinary differential equations are those built upon improved and branched off versions of the Cappasso and Paveri-Fontana model [3] and the Codeco model [5]. For example, these models include those proposed by [8,[10][11][12]14,18].
The model employed in this project is an extension of the work by [12], which considers the role of the bacteriophage, host and pathogen on the infection's sustainability. The bacteria population is assumed to be divided into two subgroups: non-pathogenic and pathogenic strains. This assumption is due to the fact that the pili protein is produced only by pathogenic strains and is necessary to anchor the bacteria to the wall of the gastrointestinal tract, so that infection can develop. It is during this crucial process that the phage targets the pili protein. The pili and toxin genes are activated by the same bacterial switch, so that toxin is produced during bacterial colonization of the gastrointestinal tract. Noting that only virulent strains of V. cholerae can produce toxin, this biological process is insured thanks to the phage and the bacterium transfer CTX (a genetic sequence that produces the toxin) through a very specific mechanism, see [1].
In the present paper, we find and discuss eight potentially equilibria, one of which is a free equilibrium point. We discuss the major epidemiological properties of the model, perform mathematical simulations and carry out sensitivity analysis to test the robustness of our proposed model.

Model formulation
Our model describes the cholera transmission taking into account host, pathogen and bacteriophage dynamics. Here, the phage population is divided into two groups: nonpathogenic and pathogenic strains. In the compartmental model presented, the host population follows a Susceptible-Infectious-Recovered (SIR) framework with S + I + R = N, where N is the total population size, which is kept constant. The bacteria population is described by the compartments B p and B np , while P represents the bacteriophage population. Observe that the interaction between bacteria and bacteriophage follows the wellknown predator-prey relationship. Hence, we use Holling 1 predation term θ(B np , P) = θ PB np , where θ is a constant of proportionality, θ(B np , P) is the conversion rate from nonpathogenic strain to pathogenic strain, induced by phage. K 1 is the half saturation constant of predation (the bacteria level at which predation occurs at half the maximum rate).
The populations' dynamics is represented by the following system of ordinary differential equations:Ṡ α(B p )S defines the incidence term, where α(B p ) is the pathogenic bacteria densitydependent component. The indirect part of the incident term is defined as Human contamination of the water supply through infected feces contributes to both bacteria and phage levels and is called shedding and it depends on both bacteria shedding of the infected individuals (a biological quantity) and the level of sanitation in the environment (an environmental assessment). Bacteria and phage shedding rates may not be the same, so the rate for the bacteria is ξ and that for the phage is ξφ, where φ is constant. The list of model parameters values and their description are seen in Table 1. The ranges of the parameters values were obtained from [4] and [10].

Forward invariance
We would like to identify a forwardly invariant set in which solutions of system (1) will be bounded. From the first equation of system (1), we see thatṠ > 0 if S = 0. Thus S(t) > 0 for t > 0. From the second equation of system (1), we see that if I = 0, thenİ = α(B p )S ≥ 0. Since α(B p ) ≥ 0 by definition, and S(t) > 0 for t > 0, we get thatİ ≥ 0. If R = 0, then form the third equation of system (1), we getṘ = δI, and hence R(t) ≥ 0. As S + I + R = N, we must have S, I, R ≤ N.
Half saturation bacteria predation density To find an upper bound for B np and B p , Moreover, F(B) = 0 has two real roots, namely The graph of F(B) is shown in Figure 1. Since B(t) ≥ 0 for t > 0, we conclude that 0 ≤ B(t) ≤ B 1 for t > 0. Thus B(t) is bounded above.
Since B = B np + B p , and both B np and B p are non-negative, we get B np , B p ≤ B 1 . Finally, to find an upper bound for P(t), the following lemma is needed.
Then for all values of B, the following is true: To show that P(t) is bounded above, let y(t) = βB + P. Then Now, by solving the ordinary differential equationẏ(t) = −Ly(t) + W, we get that We summarize the previous results is in the following proposition:

Equilibria with no shedding and their stability
In this section, we determine the equilibrium points of the model system when ξ = 0, and then perform stability analysis of the equilibria. Since α(B) is a non-smooth function, to ease our analysis, we study the dynamics of the system when α(B p ) = 0 α(B p ) = 0 separately.

Existence of equilibria
Substituting ξ = 0 in system (1) and noticing that R = N−S−I, the third equation is then not necessary. Thus system (1) reduces tȯ To analyse the types of solutions that this model could produce, a steady-state analysis was performed. The algebraic analysis is provided in the Supplementary Material. Case 1: If the pathogenic bacteria level is below or equal to the minimum infectious dose, then α(B p ) = 0. Hence system (2) becomeṡ Then system (3) has four disease-free equilibrium points which are listed below: , exits if the following conditions hold: Note that the existence of B * p and B p 2 is contrary.
Case 2: If the pathogenic bacteria level is above the minimum infectious dose, then α(B p ) = 0, leaving us with the following system: Solving the last three equations of system (4), we get the following equilibrium points: In this case, The interior point E * * = (S 1 , I 1 , B * * p , B * * np , P * * ), exits if the following conditions are hold: Note that the existence of B * * p and B p 4 is contrary.

Stability of equilibria
Depending on the pathogenic bacteria level, the linearization of system (2) has two forms, one for system (3)  Let Then, we have Then the Jacobian matrix for model system (3) is given as follows: Let J i be the Jacobian matrix evaluated at E i , i ∈ {0, 1, 2, * }. From the Jacobian matrix J 0 , it is found that its eigenvalues are −(μ + η), −(μ + δ), −d < 0 and r. Since r > 0, then E 0 is always unstable.
When considering the equilibrium point E 12 , we found that the eigenvalues corresponding to J 1 are: −(μ + η), −(μ + δ), −r, 0 and βγ 1 K K+K 1 − d. So, E 12 might be stable if But if E 2 exists, then K > dK 1 βγ 1 −d = B p 2 and hence, E 12 might be unstable whenever E 2 exists. From the Jacobian matrix J 2 , it is found that none of the eigenvalues is zero, and one of the eigenvalue is λ = r − But in order for P 2 to be positive, we must have K > B p 2 . Thus we conclude that λ > 0, and hence E 2 is unstable whenever it exists.

Equilibria with shedding and their stability
In this section, we determine the equilibrium points of the model system when ξ = 0, and then perform stability analysis of the equilibria.

Existence of equilibria
Noting that as R = N−S−I the third equation of system (1) is not necessary, leaving us with the following system: If α(B p ) = 0, then we will have the same equilibrium points E 0 , E 1 , E 2 and E * as in Section 4.1, with the same conditions. If B P > c, then α(B p ) = 0, and hence at steady state system (5) reduces to Then system (6) has one equilibrium point, namely E * * * = (S 1 , I 1 , B * * * p , B * * * np , P * * * ) which exists if the following condition holds: The proof of the existence of E * * * is provided in the Supplementary Material.

Stability of equilibria E * , E * * , E * * *
For the equilibrium points E * , E * * , E * * * we lack exact expressions for the equilibrium quantities, and so the local stability is difficult to be found analytically.

Sensitivity analysis
The objective of this section is to discuss the sensitivity of the peak value and time to model parameters. The ranges of the parameters values were obtained from [4] and [10]. We assume that recovered patient immediately become susceptible. We equally assume that within the period of our studies, deaths and births are equal. Thus our system reduces to the following system: Using the normalized forward sensitivity index, we study the sensitivity of the parameters of System (7) to cholera outbreak peak and peak time. The normalized forward sensitivity index: where x * is the quantity being considered (peak value or peak time) and p is a parameter which x * depends upon. Sensitivity indices can be positive or negative which indicates the nature of the relationship, and it is the magnitude that ranks the strength of the relationship as compared to the other parameters. Because we do not have an explicit formula for the quantities we are interested in (outbreak peak, outbreak peak time and the endemic steady state), we estimate ∂x * ∂p using the central difference approximation: We choose p = 1% of p. plugging all these in (8) we get Figure 2 contains the sensitivity of the outbreak peak and the outbreak peak time to the parameters of the model. Figure 3 shows reservoir related parameters: phage attack rate and phage burst size have the strongest relationships to the magnitude of the outbreak peak and outbreak peak time. The figure shows that if the phage attack rate and burst sizes are large, the magnitude of the outbreak peak value will be small and the time at which it occurs will be large. This indicates that phage-bacteria predator prey like cycles may have a key role to play in driving cholera outbreak cycles. Looking at the two parameters that can be controlled by policy-makers, these are K and ξ . The carrying capacity K has a large influence on the dynamics of the system. It has one of the largest sensitivity indices, being many times greater than that of the shedding rate ξ . Figure a shows the sensitivity indices of the outbreak to these two important parameters. The figure shows that the outbreak peak value increase as K but almost has no response to the variation in ξ . This further emphasizes the importance of paying attention to the dynamics in the reservoir in order to control and prevent outbreaks. The figure shows that if there are good sanitation practices that prevent shedding, but without taking care of the dynamics in the reservoir the bacteria-phage dynamics may still be able to cause outbreak in human population. This is clearly illustrated in the figure. In Figure b we plot the time course dynamics of infected population (left y-axis) and pathogenic bacteria population (right y-axis). The parameter values are chosen in such a way that we have two outbreaks every year. In such a situation, the figure illustrates that the bacteria population peaks before the human population. This indicates that the bacteria-bacteriophage cycles in the reservoir may be driving the outbreak cycles in the human population.
Our finding that the bacteria-bacteriophage cycles in the reservoir may be driving the outbreak cycles in the human population has important practical implications in terms of public health policies. Our model confirms and expands empirical and theoretical proofs that lytic bacteriophages are one of the major drivers of cholera infection, shedding light on the complex interplay between the host, pathogen and phages [15]. Cholera is one of the infectious diseases for which ecological factors and indirect transmission routes play a major role, together with those infections caused by pathogens like rotavirus, hantavirus, Legionella, Schistosoma, Cryptosporidium or Giardia. Devising mathematical models which incorporate bacteriophages can shed light on the epidemiological features of cholera and provide a better understanding of it. The seasonal self-limiting pattern of cholera epidemics can be explained taking into account the phage predation behaviour of the V. cholerae during the late stages of the outbreak, which are usually characterized by significant phage shedding in stools. Bacteriophages are able to tune/modulate the bacterial infectivity and, as such, to lead to the collapse/quenching of the outbreak, with a very low probability of a phage-resistant variant to take-over [16]. Advancements in the understanding of phages/bacteria interactions could potentially lead to the development of tools for monitoring and even predicting (forecasting/now-casting) of cholera outbreaks. For instance, it has been suggested that phage monitoring in acquatic environments and reservoirs could be a feasible strategy for controlling and mitigating against the burden imposed by cholera [16]. Moreover, recent discoveries have shown that cocktails based on a number of virulent bacteriophages can effectively prevent V. cholerae infection [20]. Phage therapy appears to be a promising pharmacological tool to counteract cholera transmission, given the organizational difficulties to implement mass vaccination campaigns in low-resource contexts as well as the lack of drugs and the issue imposed by multi-drug resistance [1,13]. In conclusion, mathematical models could play a significant role in elucidating the mechanisms of such a therapy, helping physicians optimize it [10]. As such, further research is warranted.

Conclusion
Cholera is an acute enteric infectious disease which imposes a relevant societal and clinical burden in developing countries, characterized by inadequate access to water, sanitation and lacking of proper healthcare facilities. Despite a huge body of research, the precise nature of cholera transmission dynamics has yet to be fully elucidated. Mathematical models can be useful to better understand how an infectious agent can spread and be properly controlled. In the present paper, we developed a compartmental model describing a human population, a bacterial population as well as a phage population. We showed that there might be up to eight equilibrium points; one of which is a disease free equilibrium point. We carried out numerical simulations and sensitivity analyses and we showed that the presence of phages can significantly reduce the number of infectious individuals. Moreover, we discussed the main implications in terms of public health management and control strategies. Our results agree with empirical and theoretical proofs that lytic bacteriophages are one of the major drivers of cholera infection [15]. Our findings support [16]'s suggestion that phage monitoring in acquatic environments and reservoirs could be a feasible strategy for controlling and mitigating cholera outbreaks. Moreover, our results suggest that phages could be used to efficiently treat V. cholerae infections. This is in line with the findings in [1,10,13,20]. Our proposed model can inform public health policies in low-resource settings, helping health policy-and decision-makers implement evidence-informed programs aimed at targeting cholera outbreaks and epidemics. Moreover, our model can be further expanded and refined, inspiring further research in the field.