Stability Analysis and Application of a Mathematical Cholera Model

In this paper, we conduct a dynamical analysis of the deterministic cholera model proposed in [9]. We study the stability of both the disease-free and endemic equilibria so as to explore the complex epidemic and endemic dynamics of the disease. We demonstrate a real-world application of this model by investigating the recent cholera outbreak in Zimbabwe. Meanwhile, we present numerical simulation results to verify the analytical predictions.

1. Introduction.Cholera, a water-borne disease caused by the bacterium Vibrio cholerae, continues to represent a significant public health burden in developing countries.For example, 142,311 cholera cases were officially reported in 2002 by 52 nations worldwide, among which Africa accounts for the majority.In 2006, the number of cholera cases climbed to almost 240,000 [42].Mostly recently, the 2008-2009 cholera outbreak in Zimbabwe ended up with 98,592 reported cases and 4,288 deaths [40][41][42].In addition, many believe the actual infection and death numbers have been far more than those reported due to poor surveillance and incomplete records.
Despite a large body of clinical, experimental and theoretical studies (see [1,21,25,27,29,36] and the references therein), the fundamental mechanism of transmission for cholera is not well understood at present, which has hindered effective prevention and control strategies of the disease.The difficulty stems from the complex, multiple transmission pathways which include both direct human-to-human and indirect environment-to-human modes, and which distinct cholera from many other infectious diseases.
Mathematical modeling, simulation and analysis provide a promising way to look into the nature of the cholera dynamics, and many efforts have been devoted to this topic [3,6,9,21,[29][30][31].Among these, Capasso and Paveri-Fontana [3] proposed a The remainder of this paper is organized as follows.In Section 2, we analyze the local and global stability of the disease-free equilibrium.This analysis naturally yields the basic reproduction number, R 0 , which provides a sharp threshold for the epidemic dynamics.In Section 3, we study the endemic equilibrium with regard to the long-term dynamics of the disease.We prove the existence and local asymptotic stability of the positive endemic equilibrium when R 0 > 1 , and briefly discuss its linear global asymptotic stability.In Section 4, we apply the model to investigate the Zimbabwean cholera outbreak and to verify our analytical results.Finally, we close the paper by some discussion in Section 5.
2. The disease-free equilibrium.The deterministic model of Hartley, Morris and Smith [9] consists of the following ordinary differential equations (ODE) Written in a vector form, the above equations become with Here S , I and R denote the susceptible, the infected, and the recovered numbers, respectively; N = S + I + R is the total population which is assumed to be a constant.B H and B L denote the concentrations of the hyperinfectious (HI) and less-infectious (LI) vibrios.The parameters β H and β L represent the HI and LI ingestion rates, κ H and κ L the HI and LI half saturation rates (i.e., ID 50 , the infectious dose sufficient to produce disease in 50% of those exposed), b the natural human birth/death rate, χ the bacterial transition rate, ξ the shedding rate, δ L the bacterial death rate, and γ the recovery rate.To be biologically feasible, all these parameters are positive.
It is straightforward to see that equations (1-5) have a unique disease-free equilibrium (DFE) The local stability of the DFE, which is directly related to the disease epidemics [6,12], is analyzed as follows.
The Jacobian of the ODE system (1-5) is After substituting the values for the DFE: S = N , I = R = B H = B L = 0 , the above matrix becomes The characteristic polynomial of the matrix J B is found as The equilibrium ( 8) is locally asymptotically stable if and only if all roots of the above polynomial have negative real parts.Obviously λ = −b is a negative root of multiplicity 2 .To analyze the three roots of the cubic polynomial inside the square brackets, we set Based on the Routh-Hurwitz criterion [13,28], the sufficient and necessary condition for stability is Note that the first inequality is automatically satisfied since all the model parameters are positive.The second inequality, b 3 > 0 , holds if and only if In addition, we have It is thus clear to see that b 1 b 2 −b 3 > 0 as long as the inequality (11) or, equivalently, (12), holds.The condition (12) provides a threshold for the total population (which is assumed to be completely susceptible initially): When N is below S c , the DFE is stable and no epidemicity would occur.In contrast, if N is above this critical value, the DFE becomes unstable and any infection entering the population would persist and lead to an epidemic.We define the basic reproduction number, R 0 , of this model by Biologically speaking, R 0 represents the average number of secondary infections that occur when one infective is introduced into a completely susceptible host population [12,37,38].The condition (12) is then equivalent to R 0 < 1 Thus, we have established the result below Theorem 2.1.The disease-free equilibrium of the model ( 6) is locally asymptotically stable if R 0 < 1 , and unstable if R 0 > 1 .
We mention that the basic reproduction number, given in equation ( 14), can also be derived by the next generation matrix analysis [37].
To study the global asymptotic stability of the DFE, one common approach is to construct an appropriate Lyapunov function [16,18]).We have found, however, that it is simpler to apply the following result introduced by Castillo-Chavez et al. [4].Lemma 2.2.[4] Consider a model system written in the form where X 1 ∈ R m denotes (its components) the number of uninfected individuals and X 2 ∈ R n denotes (its components) the number of infected individuals including latent, infectious, etc; X 0 = (X * 1 , 0) denotes the disease-free equilibrium of the system.
Also assume the conditions (H1) and (H2) below: ) is an M -matrix (the off diagonal elements of A are nonnegative) and Ω is the region where the model makes biological sense.
Theorem 2.3.The disease-free equilibrium of the model ( 6) is globally asymptotic stable if R 0 < 1 .
Proof.We only need to show that the conditions (H1) and (H2) hold when R 0 < 1 .In our ODE system (1-5), X 1 = (S, R) , X 2 = (I, B H , B L ) , and X * 1 = (N, 0) .We note that the system is linear and its solution can be easily found as Clearly, R(t) → 0 and S(t) → N as t → ∞ , regardless of the values of R(0) and S(0) .Thus X * 1 = (N, 0) is globally asymptotically stable.Next, we have We can then obtain which is clearly an M -matrix.Meanwhile, we find 3. Endemic dynamics.

3.1.
Existence and local stability of endemic equilibrium.The stability at the DEF determines the short-term epidemics of the disease, whereas its dynamics over a longer period of time is characterized by the stability at the endemic equilibrium.In this section we shall conduct the endemic analysis.Let us denote the endemic equilibrium of system (6) by Its components must satisfy We first show the following theorem Theorem 3.1.The positive endemic equilibrium exists and is unique if and only if R 0 > 1.
Proof.By manipulating equations (17)(18)(19)(20)(21), we obtain a cubic equation for I * : where and where S c is defined in equation (13).The zero root of equation ( 22) corresponds to the DFE.The other two (non-zero) roots, I 1 and I 2 , of the equation must satisfy: It is obvious that A < 0 since all parameters are positive.When R 0 > 1 , we have N > S c and C > 0, so that the right-hand side of equation ( 23) is negative.Hence, there is one and only one positive real root for equation (22).
On the other hand, if R 0 < 1 , then N < S c and C < 0 , so that the right-hand side of equation ( 23) is positive.Next we show B < 0 .We have when Combining ( 25) and ( 26), we obtain B < 0 .Hence, the right-hand side of equation ( 24) is negative.In this case we either have two negative real roots, or two complex conjugate roots with negative real parts, for A(I * ) 2 + BI * + C = 0 .There is no positive endemic equilibrium.
Finally, if R 0 = 1 , then C = 0 and equation ( 22) has only one nonzero root, − B A , which is negative.
We have the following result regarding the local stability of the endemic equilibrium.
Proof.Consider the Jacobian (9) at the endemic equilibrium.To make the algebraic manipulation simpler, we set ) Note that P, Q and T are all positive.The Jacobian matrix (9) then becomes Obviously, this equation has a negative root λ = −b .we expand the expression in the square brackets to obtain where To ensure that all roots of equation ( 28) have negative real parts, the Routh-Hurwitz stability criterion [13] requires Among these a 3 > 0 is obvious.The other three inequalities in (29) hold when R 0 > 1 ; the details are provided in the Appendix.
3.2.Linear global stability.The classical Poincaré-Bendixson theory [10] is a powerful tool to study global stability of nonlinear two-dimensional autonomous systems.However, this framework cannot be directly extended to higher-dimensional systems.Consequently, the global stability analysis of endemic equilibria for nonlinear higher-dimensional problems, such as the one in (1-5), is generally difficult, despite the efforts made by several authors [14][15][16][17]26,39].For some special differential equation systems, one might be able to find a suitable Lyapunov function [18,20] to prove the global stability.Unfortunately, there is no systematic way to construct or find Lyapunov functions, which hinders the application of this approach to more general model systems.
We speculate that the unique positive endemic equilibrium of system (1-5) is globally asymptotically stable when R 0 > 1 .In the present paper, we will use a special linear case (see Case one below) to illustrate this point.We will also present numerical evidence (see Section 4) for nonlinear justification.It is our plan to rigorously investigate the global endemic stability of this model in another study.
Case one.We assume the pathogen concentrations in the environment are far beyond the HI and LI half saturation rates (ID 50 ), i.e., B L >> κ L and B H >> κ H .Under this assumption, the incidence rates in the model become That is, the possibility of infection is about 100% to those exposed to pathogens.Since R = N −I −S , the model (1-5) is reduced to a two-dimensional linear system: with β = β L + β H .It is straightforward to determine the unique positive endemic equilibrium of this system: There is no DFE in this case for obvious reason.Meanwhile, the exact solution of the system (30, 31) can be quickly found as follows It is thus clear to see that S(t) → S * and I → I * as t → ∞, regardless of the initial values of S and I .Hence, the endemic equilibrium S * , I * is globally asymptotically stable.This conclusion can also be easily obtained based on the famous Bendixson Theorem [10,14], by noting that holds everywhere, so that no periodic orbit can exist.Case two.We now assume the pathogen concentrations are much lower than the half saturation rates, i.e., B L << κ L and B H << κ H . Then we have which implies the chance of getting new infection is about 0 .The model system (1-5) is then reduced to There is no endemic equilibrium in this case and the only equilibrium point of system (32,33) is a DFE: The exact solution of equations ( 32) and ( 33) is Clearly, S(t) → S 0 and I(t) → I 0 as t → ∞ , confirming the global stability of the DFE (indeed, this is a special case of Theorem 2.3).The same conclusion can be obtained by using the Bendixson Theorem.

Bifurcation diagram.
To summarize our stability analysis results, we sketch a bifurcation diagram [35] of I vs. R 0 for system (1-5) in Figure 1, which highlights the stability exchange at R 0 = 1 for the two biologically feasible equilibria: the DFE and the positive endemic equilibrium.To illustrate the positive endemic branch of the bifurcation, we write the endemic quadratic equation A(I * ) 2 + BI * + C = 0 (see equation 22) as where the constant coefficients D, E and F are given by We then obtain When I * is small or moderate, equation (35) represents approximately a straight line passing the bifurcation point R 0 = 1 , I * = 0 .
In contrast, the values of β L and β H are less well known.The authors of [9] took β L = 1.5/wk and β H as a variable; they also stated that these two parameters are  particularly "difficult to estimate".This suggests that β L and β H are sensitive in the ODE model (1)(2)(3)(4)(5).To verify it, we have carried out a sensitivity analysis for these two parameters using the Latin Hypercube Sampling (LHS) method [2,22].The LHS method, which is one of the most efficient ways to analyze sensitivity and uncertainty, calculates the Partial Rank Correlation Coefficients (PRCC) in the range of [−1, 1] for the model parameters with respect to outcome measurements.A small PRCC value (close to 0) indicates that changes in the input parameter have little impact on the outcome, whereas a bigger PRCC value suggests that the outcome is significantly influenced by, thus sensitive to, the changes of the input parameter.In this study we pick the total infection, I , as the outcome and present the results in Table 1 with a sample size n = 400 ; similar patterns are observed for differing sample sizes.The results confirm the sensitivity of the two parameters β L and β H ; in particular, we observe that the number of the infected people is very sensitive to β H .In the work of [9], the ODE model (1-5) was introduced and discussed for a hypothetical community with a total population of N = 10, 000 .In order to assess the applicability of this model to the Zimbabwean cholera outbreak, we will first need to ensure the available Zimbabwean data can be reasonably fitted by the model results, which necessitates the adjustment of parameter values.The discussion in [9] and our sensitivity analysis indicate that the parameters β L and β H are sensitive and possibly vary from country to country.Hence, in what follows we adjust these two parameters to match the reported infections in Zimbabwe.
According to the published data by WHO [41,42], the cholera outbreak started in the end of August 2008, with 11,735 cases reported by December 1, 2008, 79,613  = 0.9157 after about 5,000 weeks.
cases by February 18, 2009 and 91,164 cases by March 17, 2009; afterwards the situation had improved significantly due to various control measures and international aids.Since the total population in Zimbabwe is about 12 million [42], we scale down this number by a factor of 1,200 to match the hypothetical population of N = 10, 000 in [9].Accordingly, we now have 10, 66 and 76 (normalized) cholera infections after 12 weeks, 22.5 weeks and 26.5 weeks, respectively, from the beginning of the outbreak.By varying the two parameters, we find that when β L = 0.126 and β H = 0.0995 , the model prediction for the infected population, I , can well fit these data (see Figure 2).By substituting the values of β L , β H and other parameters into equations ( 13) and ( 14), we find the basic reproduction number and the population threshold The fact that R 0 > 1 justifies the development of the cholera epidemics.The relatively low value of R 0 for the Zimbabwean cholera outbreak is attributed to the fact that although nearly 0.1 million cases have been reported, the overall percentage of infection with respect to the total population is still pretty low (less than 0.8% ).The R 0 estimated here is regarded as a nationally averaged reproduction rate.Meanwhile, we substitute these parameter values into equations ( 22) and  Note again that we have scaled down the total population in Zimbabwe by a factor of 1, 200 to match the hypothetical population N = 10, 000 .Thus, the model predicts that the realistic endemic infection number in Zimbabwe is about 1, 099 .
To verify the model prediction, we run the numerical simulation for a much longer period of time (up to 7,000 weeks) and present the results for I , S and R in Figures 3 and 4. The first peak of the infection curve in Figure 3 represents the 2008-2009 cholera outbreak.The infection number (I) starts to decline once the susceptible population (S) falls below the threshold, S c .= 7657 .After this outbreak I drops to almost zero, meaning that the majority of the infected people have recovered and entered the R class, so that we see a significant increase of R in Figure 4. Then I stays at the zero level for the next 1,500 weeks or so (about 30 years).During this period R gradually decreases due to natural death of recovered individuals, whereas S gradually increases due to continuous birth of new susceptibles.Once the susceptible population exceeds the threshold S c , another cholera outbreak is triggered but with much lower magnitude.This pattern continues for a few more outbreaks with decaying magnitudes.After about 5,000 weeks, the infection curve rests at the endemic value, I * .= 0.9157 ; the S and R curves also converge to their endemic values, S * .
= 7666 and R * .= 2333 , respectively.Figures 5 and 6 show the results of another numerical run with different initial conditions: I(0) = 500, S(0) = 8500, R(0) = 1000, B H (0) = B L (0) = 0 .We observe very similar pattern.In particular, the I , S and R curves all approach their endemic equilibrium values after about 5,000 weeks.This pattern is also observed with various other initial conditions, which demonstrates the global asymptotic stability of the endemic equilibrium when R 0 > 1 .Meanwhile, these results also justify the instability of the disease-free equilibrium, as cholera outbreaks occur whenever the susceptible population S exceeds the critical value S c .
Finally, to validate the global asymptotical stability of the disease-free equilibrium when R 0 < 1 , we have conducted many numerical runs with differing population N and initial conditions, and the results agree with the model prediction.A typical numerical result is presented in Figure 7, where the total hypothetical population is set as N = 5, 000 (i.e., halving the current Zimbabwean population)   so that R 0 .= 0.653 .Initially I is 500 ; then the infection number quickly drops to zero and stays at zero for all time afterwards.5. Discussion.We have conducted a stability analysis of the mathematical cholera model presented in equations (1-5), originally proposed by Hartley, Morris and Smith [9].We have studied the stability property at both the disease-free equilibrium (which determines the short-term epidemic behavior) and the endemic equilibrium (which determines the long-term disease dynamics).Our analysis shows a regular stability exchange at R 0 = 1 for this combined human-environment epidemiological model.To verify the analytical theory, we have performed numerical simulations and employed the most recent Zimbabwean cholera outbreak as a realistic case study.Our simulation results demonstrate that the model (1-5) is applicable to the cholera dynamics in Zimbabwe; the results for both epidemic and endemic dynamics are consistent with the analytical predictions.
The findings in this paper have several implications to the initiation and persistence of cholera infection.As mentioned before, one distinct feature of the model (1-5) is the incorporation of an hyperinfectious (H) state of the pathogen.Our sensitivity analysis shows that the number of infections is very sensitive to the parameter β H . Together with the fact that the model can well fit the Zimbabwean data, our results that the hyperinfectious state plays an important role in the transmission of the disease, especially for the onset of an epidemic.This point of view was carefully discussed in [9].As pointed out in [30], the hyperinfectious state introduced in the model (1)(2)(3)(4)(5) can also be interpreted as a another way of modeling the direct human-to-human transmission.The model, when applied to Zimbabwean cholera, predicts that there will be a (smaller) cholera outbreak in every 20-30 years, and that the disease will remain endemic in the long term, unless there are significant changes to the current population and/or environment in Zimbabwe.We emphasize that this observation is purely based on the model system (1-5) and under the idealized condition that the population remains a constant for all the time.Thus it should be regarded as a model prediction under simplified settings, which, nonetheless, provides insight into the future evolution of cholera dynamics in this country.It is expected that similar studies can be performed to other cholera-endemic countries such as Mozambique, Zambia, and Bangladesh, with possible modifications of the model (e.g., adjustment of parameters).
The analysis and results presented in this work also enable us to study the control measures [7,12] against cholera outbreak.For a simple illustration, we consider the use of vaccination and medical therapy (which could include hydration therapy, antibiotics, etc.) to control the disease.We assume that vaccination is introduced to the susceptible population at a rate of v , so that vS individuals per time are removed from the susceptible class and added to the recovered class.Similarly, we assume that medical therapy is applied to infected people at a rate of u , so that uI individuals per time are removed from the infected class and added to the recovered class.With the incorporation of these controls into the cholera model, equations (1-3) become whereas equations ( 4) and ( 5) remain the same.Conducting similar analysis as in Sec. 2 to this control model, we can easily obtain the (modified) basic reproduction number: where R 0 is the basic reproduction number defined in equation ( 14) for the original model.Clearly R c 0 < R 0 .The result in equation (39) shows that, mathematically, each of the two types of individual controls can reduce the value of R c 0 below 1 so that the disease will be eradicated.For example, if v = 0 (no vaccination), we would just need the rate of the medical therapy to be such that u > (γ + b)(R 0 − 1) to ensure R c 0 < 1 .Practically, however, the strength and the success of each control measure would be limited by social and economic factors and available resources, and the combination of different prevention and intervention approaches would possibly achieve the best result.Other types of control measures, such as sanitation, sewage treatment, and water cleaning, can be also incorporated into the model and the resulting values of R c 0 can be similarly calculated.Such information would provide useful guidelines for the public health administrations to effectively design disease control strategies and to properly scale their efforts.
There are a few limitations in this study.As mentioned before, the model (1-5) is based on the assumption that the natural birth and death rates are always equal.This may not be true in the real world, especially when a long period of time is concerned.Next, the disease-induced mortality is not included in the model.In fact, to our knowledge, none of the published mathematical cholera models (see, e.g., [3,6,9,21,[29][30][31]) have considered this effect.This is perhaps due to the fact that most cholera infections are minor or moderate and can be effectively treated by simple hydration therapy; thus, in average, cholera-induced death rate is just about 1% [23,42].In contrast, the Zimbabwean cholera outbreak has a relatively high disease-caused death rate of 4.3% in average [23,42].It is thus expected that a model incorporating the disease mortality would yield more accurate results in modeling the Zimbabwean cholera.When the disease mortality is included and different natural birth and death rates are used, the total population N will vary in time.Specifically, it will be determined by the equation where b 1 and b 2 denote the natural birth and death rates, respectively, and α denotes the disease mortality rate.The varying population size adds some complexity to the dynamical system, though, in principle, similar epidemic and endemic analysis can be conducted.Such models have been studied for several infectious diseases with regular SIR or SEIR formulations (see, e.g., [12,16]).It is to be seen how the incorporation of the varying population size will change the dynamics of a combined human-environment epidemiological model system, particularly, the cholera model.
Acknowledgements.We would like to thank the two anonymous referees for their valuable comments to improve this paper.
Appendix.Based on Theorem 3.1, there is a unique positive endemic equilibrium when R 0 > 1 .In what follows we show the last three inequalities in (29) at the endemic equilibrium.First we rewrite a 0 as Substituting equations ( 20) and ( 21) into ( 17) and ( 18), we obtain Using equations ( 27) and (42), we have Next, we rewrite a 1 into the sum of three parts: Note that the last part in equation ( 43) is positive.After substitution of equations ( 27) and (42), the first two parts of a 1 become b Thus a 1 > 0 holds.

Figure 1 .
Figure 1.The bifurcation diagram of I vs. R0 for the DFE and the positive endemic equilibrium.The solid lines represent the stable equilibria, and the dashed line represents the unstable equilibrium.There is a stability exchange at R0 = 1 .Biologically non-feasible equilibria are not shown on the diagram.

Figure 2 .Figure 3 .
Figure 2. The data fitting for the infected population I , where the curve represents the model prediction and the squares mark the reported Zimbabwean data (normalized by a factor of 1,200).

Table 1 .
LHS sensitivity analysis for βL and βH with respect to the infection number I .The sample size is n = 400 .Results show that the infection is especially sensitive to βH .