UMYU Scientifica

A periodical of the Faculty of Natural and Applied Sciences, UMYU, Katsina

ISSN: 2955 – 1145 (print); 2955 – 1153 (online)

ORIGINAL RESEARCH ARTICLE

A Deterministic Model of Cholera Transmission with Vaccination, Environmental Reservoirs, and Vector Dynamics

Abdullahi Ahmed^1,2,3, Adamu Ishaku^2,3, Ahmadu Bappah Muhammad^2,3 and Ayuba Sanda^2,3

¹Department of Mathematics and Computer Science, College of Education Billiri, Gombe, Nigeria

²Department of Mathematical Sciences, Gombe State University, Gombe, Nigeria

³GSU-Mathematics for Innovative Research (GSU-MIR) Group, Department of Mathematics, Gombe State University, Gombe, Nigeria

Corresponding Author: Abdullahi Ahmed [email protected]

Abstract

This research develops an extended deterministic model to investigate the transmission dynamics of Cholera by incorporating three critical components: vaccinated human population, an exposed human compartment, and vector recruitment, which was modelled using logistic growth that contributes to environmental contamination. The model was analyzed by studying the stability of equilibrium points. The results of the analysis show that the basic reproduction number, R₀, controls disease dynamics. It is shown that the disease-free equilibrium is locally asymptotically stable when R₀ < 1 and unstable when R₀ > 1, leading to endemic persistence. In an attempt to examine the effect of some parameters of the disease dynamic, parameter estimation was performed using available epidemiological data, and model validation was conducted using root mean square error (RMSE) and coefficient of determination (R²); sensitivity analysis was employed using Partial Rank Correlation Coefficients (PRCC) revealing that transmission rates and environmental contamination parameters significantly influence disease spread. Finally, numerical simulations are also performed to verify the analytic results. Analytical results show that cholera transmission is strongly influenced by the combined effects of vaccination coverage, environmental contamination, and vector density. Although increasing vaccination reduces transmission, the presence of reinfection and environmental reservoirs may sustain the disease even when the reproduction number is close to unity. Numerical simulations are performed to support the theoretical findings and to illustrate the impact of key parameters on disease dynamics. Generally, this research highlights the importance of integrating control strategies by incorporating vaccination, access to clean water, proper sanitation, vector control, and also early treatment in achieving effective and sustainable cholera control in endemic settings.

Keywords: Cholera; Mathematical modelling; Stability analysis; Sensitivity analysis; Numerical simulations.

Introduction

Cholera is a short-term, life-threatening disease caused by a bacterium called Vibrio cholerae. The disease affects the intestine and causes diarrhoea (WHO et al., 2011). Cholera exists in different serogroups, but those responsible for outbreaks are subgroups 01 to 0139 (CDC, 2019). Diarrhea the main symptom of Cholera and is often referred to as acute diarrheal disease. It is accompanied by severe dehydration, which can lead to death within hours if not treated. The disease remains asymptomatic for 12 hours to 5 days after ingestion of contaminated Food or water. Nonetheless, infected individuals may still shed the bacterium, contaminating the environment and posing risks to others (WHO et al., 2011). Cholera is transmitted primarily through the ingestion of Food or water contaminated with V. cholerae, making unhygienic environments highly vulnerable to outbreaks. Furthermore, direct human-to-human contact, such as handshakes, can facilitate transmission (Yang, Wang, Gao, & Wang, 2017). People with low immunity, such as malnourished children or people living with HIV, have a higher tendency to develop Cholera (Hove-Musekwa et al., 2011). To prevent the spread of Cholera, proper environmental sanitation is necessary, along with access to clean water and safe Food. Vaccination also plays a significant role in reducing the prevalence of the disease (WHO et al., 2011). Treatment primarily (CDC, 2019) involves rapid replacement of lost fluids and salts using oral rehydration solutions (ORS) (WHO et al., 2011). Individuals who recover from Cholera acquire immunity that can protect them for many years (Harris, LaRocque, Kadri, Ryan, & Calderwood, 2012).

Globally, Cholera affect between 1.3 and 4.0 million cases are estimated and 21,000 to 143,000 deaths from Cholera each year (WHO et al., 2011). The current pandemic began in 1961 in Indonesia and spread to Europe, the South Pacific and Japan by the late 1970s and to South America by the 1990s. Notable outbreaks have occurred in India (2007), Congo, Zimbabwe, and Iraq (2008), Vietnam and Zimbabwe (2009), and Nigeria and Haiti (2010). In Nigeria, Cholera is a recurring disease, often peaking during the rainy season. Major epidemics were recorded in 1970, 1990, and between 1992–1997. In 2010, the Federal Ministry of Health reported 37,289 cases and 1,434 deaths between January and October. In 2011, there were 22,797 cases and 728 deaths. The Nigeria Centre for Disease Control (NCDC) reported 42,466 suspected cases and 830 deaths in 2018 (CDC, 2019). In 2021, Cholera led to 103,589 infections and 3,566 deaths—surpassing the COVID-19 death toll of 2,977 in the same year (CDC, 2021). As of August 11, 2024, a total of 5,951 suspected cholera cases and 176 deaths had been reported across Nigeria’s 36 states (CDC, 2024). Cholera is a food-borne enteric disease caused by the bacterium Vibrio cholerae. Globally, it accounts for an estimated 1.6% to 3.5% of food-borne disease-related deaths and 1.3 to 4.0 million cases annually. In developing countries, food-borne diseases like Cholera remain prevalent due to poor hygiene, lack of education, inadequate drainage systems, and the abundance of carriers (Misra, Mishra, Pathak, & et al., 2013). Promoting awareness about safe food storage and enhancing environmental sanitation can reduce the incidence of food-borne diseases (Lata, Mishra, & Misra, 2020). These diseases are also transmitted by vectors such as houseflies, ticks, and mites. Houseflies, in particular, are known mechanical vectors for various pathogens, including bacteria (Levine & Levine, 1991). Studies suggest that houseflies contribute significantly to cholera transmission (Keiding, 1986), especially during seasons of high fly density (Levine & Levine, 1991). Although humans rarely interact directly with houseflies, the fly population in the environment still significantly influences cholera transmission dynamics (Lata et al., 2020). Studies estimate that infected vectors can shed approximately one thousand Vibrio per gram of stool for one day (Nelson, Harris, Morris, Calderwood, & Camilli, 2009). Therefore, vectors such as houseflies play a crucial role in the transmission of Cholera.

Several mathematical models have been developed to study cholera transmission (Fister, Gaff, Lenhart, & Schaefer, 2016). A common model structure is the SIRB framework, which incorporates bacterial concentration in contaminated water as an environmental class alongside human compartments. Extensions of the SIRB model include asymptomatic infectious compartments and hyper-infectious bacteria classes (Das, Mukherjee, & Sarkar, 2005), while others include partial immunity and person-to-person transmission (Pascual, Koelle, & Dobson, 2006). Another study (Edward & Nyerere, 2015) explored the impact of vaccination, treatment, sanitation, and awareness campaigns on cholera control. While many models focus on direct human-to-human transmission and indirect transmission via contaminated water, fewer studies have addressed the role of vectors in cholera spread, such as (Kaka & Anteneh, 2024; Lata et al., 2020; Al-Shanfari, ELmojtaba, & Alsalti, 2019; Omondi, 2016; Fotedar, Banerjee, Samantray, & Shriniwas, 1992; Levine & Levine, 1991). Vectors are increasingly recognized as major contributors to the global transmission of Cholera, especially in developing countries (Das et al., 2005). A model proposed in (Kaka & Anteneh, 2024) introduced a compartmental framework involving human hosts, vectors, and environmental reservoirs. However, the model did not consider vaccination and exposure compartments for humans. In this study, we extend the model in (Kaka & Anteneh, 2024) by incorporating a vaccinated class and an exposed human compartment. We also model vector recruitment using a logistic growth function instead of a constant rate, for a more realistic representation. Contaminated and uncontaminated zones remain. The inclusion of the exposed human class targets the model’s realism, capturing the latent period before individuals become infectious and helping to better understand cholera transmission dynamics.

Model Formulation

The model considers two host populations: human and vector populations. It incorporates both direct and indirect environmental transmission by modeling the dynamics of free-living Vibrio cholerae concentrations in safe and unsafe environments. The human population is divided into five compartments: Susceptible individuals S_h(t), Vaccinated individuals V_h(t), Exposed individuals E_h(t), Infected individuals I_h(t), Recovered individuals R_h(t). Susceptible individuals are vaccinated at a constant rate ω, while vaccine protection wanes at rate ψ. The environmental bacterial concentration is denoted by P (t) (safe environment) and C(t) (unsafe environment). We define: β₁: Environment-to-human transmission rate; β₂: Human-to-human transmission rate; k: Half-saturation constant of the bacterial concentration. The force of infection for humans is defined by:

\(\lambda_{h} = \frac{\beta_{1}P}{k + P} + \beta_{2}I_{h}\) (1)

The vector population in the model is divided into two compartments: susceptible vectors S_v(t), and exposed vectors E_v(t). In this population, exposure occurs when susceptible vectors come into contact with infected human feces in an unsafe environment C(t), or with a previously safe environment P (t) that has been contaminated. The rate at which susceptible vectors are exposed, denoted by λ_v, is modeled as a nonlinear function that incorporates saturation effects:

\(\lambda_{v} = \frac{\tau_{1}C}{k_{1} + C} + \frac{\tau_{2}P}{k_{2} + P}\) (2)

Here, τ₁ and τ₂ represent the contact rates with unsafe and safe environments, respectively, while k₁ and k₂ are the corresponding half-saturation constants.

Figure 1: Model Diagram

The model equation is given by

\[{\frac{dS_{h}}{dt} = \Lambda_{h} - \lambda_{h}S_{h} - (\mu_{h} + \omega)S_{h} + \psi V_{h} }{\frac{dE_{h}}{dt} = \lambda_{h}S_{h} + \sigma\lambda_{h}V_{h} + \varepsilon\lambda_{h}R_{h} - (\mu_{h} + \pi_{h})E_{h} }{\frac{dV_{h}}{dt} = \omega S_{h} - (\mu_{h} + \psi)V_{h} - \sigma{\widehat{\lambda}}_{h}V_{h} }{\frac{dI_{h}}{dt} = \pi_{h}E_{h} - (\mu_{h} + d + \gamma)I_{h} }{\frac{dR_{h}}{dt} = \gamma I_{h} - \mu_{h}R_{h} }{\frac{dS_{v}}{dt} = r_{v}N_{v}\left( 1 - \frac{N_{v}}{L} \right) - \lambda_{v}S_{v} - \mu_{v}S_{v} }{\frac{dE_{v}}{dt} = \lambda_{v}S_{v} - \mu_{v}E_{v} }{\frac{dC}{dt} = \eta I_{h} + \alpha_{1}E_{v} - \delta_{1}C }\]\(\frac{dP}{dt} = \alpha_{2}E_{v} - \delta_{2}P\) (3)

with initial conditions

S_h(0) > 0, E_h(0) ≥ 0, V_h(0) ≥ 0, I_h(0) ≥ 0, R_h(0), S_v(0) > 0, E_v(0) ≥ 0, C(0) ≥ 0, P (0) ≥ 0

Basic Properties of the Model

Positivity of Solutions

The model (3) monitors the population, which cannot be negative; thus, we verify the positivity of the solutions to the model system defined in (3), assuming non-negative initial conditions. This step confirms that the populations and concentrations involved do not attain negative values at any time t ≥ 0. Let (S_h = x₁, E_h = x₂, V_h = x₃, I_h = x₄, R_h = x₅, S_v = x₆, E_v = x₇, C = x₈, P = x₉) in system of equations (3); we get

\[{\frac{dx_{1}}{dt} = \Lambda_{h} - \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} - (\mu_{h} + \omega)x_{1} + \psi x_{3} }{\frac{dx_{2}}{dt} = \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} + \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} + \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)R_{h} - (\mu_{h} + \pi_{h})E_{h} }{\frac{dx_{3}}{dt} = \omega x_{1} - (\mu_{h} + \psi)x_{3} - \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} }{\frac{dx_{4}}{dt} = \pi_{h}x_{2} - (\mu_{h} + d + \gamma)x_{4} }{\frac{dx_{5}}{dt} = \gamma x_{4} - \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{5} - \mu_{h}x_{5} }{\frac{dS_{v}}{dt} = r_{v}(x_{6} + x_{7})\left( 1 - \frac{(x_{6} + x_{7})}{L} \right) - \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6} - \mu_{v}x_{6} }{\frac{dx_{7}}{dt} = \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6} - \mu_{v}x_{7} }{\frac{dx_{8}}{dt} = \eta x_{4} + \alpha_{1}x_{7} - \delta_{1}x_{8} }\]\(\frac{dx_{9}}{dt} = \alpha_{2}x_{7} - \delta_{2}x_{9}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4)\)

with initial conditions

x_i(0) ≥ 0 i = 1, 2, . . . , 9

Lemma 3.1. If the initial conditions satisfy x_i(0) ≥ 0 for all i = 1, 2, . . . , 9, then the corresponding solutions x_i(t), for i = 1, 2, . . . , 9, remain positive for all t ≥ 0.

Proof. We proof the lemma by contradiction, following the method used for establishing positivity of solutions in systems of nonlinear equations as discussed in (Ishaku & Hussaini, 2020; Hai-Feng & Na-Na, 2012).

Suppose, contrary to the claim, that there exists a time t_j such that

\(x_{j}(t_{j}) = 0,\mspace{6mu}\frac{dx_{j}(t_{j})}{dt}\) < 0, for some \(j = 1,2,3,...,9\) (5)

\(x_{i}(t)\) > 0 for \(i \neq j,\quad i = 1,2,3,...,9,\quad\) for 0 < \(t\) < \(t_{j}\) (6)

Since all the parameters of the model are positive, we evaluate x^′ (t_j) for each j under equation

(5), using the system of equations in model (3):

For \(j = 1\) then \(x_{1} = 0,\) \(\frac{dx_{1}(t_{1})}{dt} = \Lambda_{h} + \psi x_{3}\) > 0 (7)

For \(j = 2\) then \(x_{2} = 0,\)

\(\frac{dx_{2}(t_{2})}{dt} = \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} + \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} + \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{5}\) > 0 (8)

For \(j = 3\) then \(x_{3} = 0,\) \(\frac{dx_{3}(t_{3})}{dt} = \omega x_{1}\) > 0 (9)

For \(j = 6\) then \(x_{6} = 0,\) \(\frac{dx_{6}(t_{6})}{dt} = r_{v}x_{7}\left( 1 - \frac{x_{7}}{L} \right)\) > 0 (10)

For \(j = 7\) then \(x_{7} = 0,\) \(\frac{dx_{7}(t_{7})}{dt} = \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6}\) > 0 (11)

Similarly, one can show that x^′ (t_j) > 0 for j = 3, 4, 5, 6, 7, 8, 9, using the fact that x_i(t) > 0 for i ̸= j and all parameters are positive.

This contradicts the assumption that x_j(t_j) = 0 and x^′ (t_j) < 0, since we have shown that

x^′ (t_j) > 0 in all cases.

j

Therefore, the assumption must be false, and the solution satisfies

x_i(t) > 0 for all i = 1, 2, 3, . . . , 9, and t ≥ 0.

Invariant Region

Theorem 3.1. the close set \(\Omega_{h} = \left\{ (S_{h},\mspace{6mu} E_{h},\mspace{6mu} V_{h,}\mspace{6mu} I_{h},\mspace{6mu} R_{h},) \in \mathfrak{R}_{+}^{5}:N_{h} \leq \frac{\Lambda_{h}}{\mu_{h}} \right\}\)and \(\Omega_{v} = \left\{ (S_{v},\mspace{6mu} E,) \in \mathfrak{R}_{+}^{2}:N_{v} \leq L \right\}\)

+

µ_h

+

are positively invariant and attracting with respect to the model (3). The biologically feasible region is given by the set Ω = Ω_h × Ω_V .

Proof. Adding the equations in model (3) where N_h = S_h + V_h + E_h + I_h + R_h, N_v = S_v + E_v

denote the total human and vector populations, respectively.

\[\Omega_{h} = \left\{ (S_{h},\mspace{6mu} E_{h},\mspace{6mu} V_{h,}\mspace{6mu} I_{h},\mspace{6mu} R_{h},) \in \mathfrak{R}_{+}^{5}:N_{h} \leq \frac{\Lambda_{h}}{\mu_{h}} \right\}\]

\[\frac{dN}{dt} = \Lambda_{h} - \mu_{h}N_{h} - dI_{h}\]

Since all the model parameters are positive for all t > 0, then the following inequality holds,

\(\frac{dN}{dt} \leq \Lambda_{h} - \mu_{h}N_{h}\) (12)

\(\frac{dN}{dt} + \mu_{h}N_{h} \leq \Lambda_{h}\) (13)

Solving the above differential inequality using the integrating factor procedure, we get

\[{\frac{dx_{1}}{dt} = \Lambda_{h} - \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} - (\mu_{h} + \omega)x_{1} + \psi x_{3} }{\frac{dx_{2}}{dt} = \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} + \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} + \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)R_{h} - (\mu_{h} + \pi_{h})E_{h} }{\frac{dx_{3}}{dt} = \omega x_{1} - (\mu_{h} + \psi)x_{3} - \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} }{\frac{dx_{4}}{dt} = \pi_{h}x_{2} - (\mu_{h} + d + \gamma)x_{4} }{\frac{dx_{5}}{dt} = \gamma x_{4} - \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{5} - \mu_{h}x_{5} }{\frac{dS_{v}}{dt} = r_{v}(x_{6} + x_{7})\left( 1 - \frac{(x_{6} + x_{7})}{L} \right) - \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6} - \mu_{v}x_{6} }{\frac{dx_{7}}{dt} = \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6} - \mu_{v}x_{7} }{\frac{dx_{8}}{dt} = \eta x_{4} + \alpha_{1}x_{7} - \delta_{1}x_{8} }\]\(\frac{dx_{9}}{dt} = \alpha_{2}x_{7} - \delta_{2}x_{9}\) (14)

So, the solution becomes

As \(t \rightarrow \infty\) \(\lim_{t \rightarrow \infty}\)sup\(N_{h}(t) \leq \frac{\Lambda_{h}}{\mu_{h}}\) (15)

Thus, the total human population N_h(t) is bounded above by \(\frac{\Lambda_{h}}{\mu_{h}}\)as t → ∞. Further, \(N_{h}(t) \leq \frac{\Lambda_{h}}{\mu_{h}}\)

Thus, Ω is positively invariant.

For Vector population

\(\Omega_{v} = \left\{ (S_{v},\mspace{6mu} E,) \in \mathfrak{R}_{+}^{2}:N_{v} \leq L \right\}\) (16)

solving the logistic growth equation for vector population

\(\frac{dN_{v}}{dt} = r_{v}N_{v}\left( 1 - \frac{N_{v}}{L} \right) - \mu_{h}N_{v}\) (17)

Let \(y = N_{v}\) solving the logistic equation

\(\frac{dy}{dt} \leq r_{v}y\left( 1 - \frac{y}{L} \right)\) (18)

Using Variable separable

\(\frac{1}{y\left( 1 - \frac{y}{L} \right)}dy \leq r_{v}dt\) (19)

Putting \(y(t) = N(t)\)

(20)

As \(t \rightarrow \infty\)

\(\lim_{t \rightarrow \infty}\)sup\(N_{v}(t) \leq \frac{L}{1 + 0} = L\) (21)

Thus, the total vector population N_v(t) is bounded above by L (carrying capacity) as t → ∞.

Further, N_v(t) ≤ L. Thus, Ω is positively invariant.

However, for the pathogen, the boundedness is shown as follows:

\(\frac{dC}{dt} = \eta I_{h} + \alpha_{1}E_{v} - \delta_{1}C\) (22)

But

\[I_{h} \leq \frac{\Lambda_{h}}{\mu_{h}},\quad E_{v} \leq L\]

Then

\(\frac{dC}{dt} \leq \eta\frac{\Lambda_{h}}{\mu_{h}} + \alpha_{1}L - \delta_{1}C\) (23)

By integrating this inequality using the integrating factor , the solution becomes

(24)

Where a is a constant, then

\(\lim_{t \rightarrow \infty}\) \(C(t) \leq \frac{\eta\Lambda_{h}}{\delta_{1}\mu_{h}} + \frac{\alpha_{1}}{\delta_{1}}L\) (25)

By integrating this inequality using the integrating factor , the solution becomes

(26)

Where A is a constant, then

\(\lim_{t \rightarrow \infty}\) \(P(t) \leq \frac{\alpha_{2}}{\delta_{2}}L\) (27)

Hence

\[0 \leq N_{h}(t) \leq \frac{\Lambda_{h}}{\mu_{h}},\quad\mspace{6mu} 0 \leq N_{v}(t) \leq L,\quad 0 \leq C(t) \leq \frac{\eta\Lambda_{h}}{\delta_{1}\mu_{h}} + \frac{\alpha_{1}}{\delta_{1}}L,\quad 0 \leq P(t) \leq + \frac{\alpha_{2}}{\delta_{2}}L\]

which implies that N_h, N_v, and all other variables (S_h, V_h, E_h, I_h, R_h, S_v, E_v, C, P ) of model equations (3) are bounded and all the solutions starting in Ω approach, enter, or stay in Ω. Hence, the set Ω is positively invariant, that is, all the solutions in Ω remain in Ω for t ≥ 0.

Disease-Free Equilibrium

The Disease-Free Equilibrium (DFE) of the cholera model is a steady-state in which no infection exists in the population or the environment. It is defined by setting the right hand side of system (3) to zero;

Thus, the disease-free equilibrium of the model is

\[E_{0} = \left( S_{h}^{0},E_{h}^{0},V_{h}^{0},I_{h}^{0},R_{h}^{0},S_{v}^{0},E_{v}^{0},C^{0},P^{0} \right)\]

\(E_{0} = \left( \frac{\Lambda_{h}(\mu_{h} + \psi)}{\mu_{h}(\mu_{h} + \omega + \psi)},0,\frac{\omega\Lambda_{h}}{\mu_{h}(\mu_{h} + \omega + \psi)},0,0,L\left( 1 - \frac{\mu_{h}}{r_{v}} \right),0,0,0 \right)\) (28)

Basic Reproduction Number

In epidemiology, a key threshold known as the basic reproduction number is defined as the average number of secondary infectious cases generated by a single primary infectious case introduced into an entirely susceptible population. To compute R₀, we use the next-generation matrix approach described in (van den Driessche & Watmough, 2002), where R₀ is obtained as the largest eigenvalue (spectral radius) of the matrix

\(R_{0} = \rho\left( FV^{- 1} \right)\) (29)

Evaluating both F and V by taking the partial derivatives of the terms in F and V at the disease free equilibrium we have.

\(F = \begin{bmatrix} 0 & \left( S_{h}^{0} + \sigma V_{h}^{0} \right)\beta_{2} & 0 & 0 & \left( S_{h}^{0} + \sigma V_{h}^{0} \right)\frac{\beta_{1}}{k} \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & S_{v}^{0}\frac{\tau_{1}}{k_{1}} & S_{v}^{0}\frac{\tau_{2}}{k_{2}} \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix},\quad\mspace{6mu}\quad V = \begin{bmatrix} a_{1} & 0 & 0 & 0 & 0 \\ \pi_{h} & a_{2} & 0 & 0 & 0 \\ 0 & 0 & \mu_{v} & 0 & 0 \\ 0 & - \eta & - \alpha_{1} & \delta_{1} & 0 \\ 0 & 0 & - \alpha_{2} & 0 & \delta_{2} \end{bmatrix}\) (30)

Where

\(a_{1} = \mu_{h} + \pi_{h},\quad a_{2} = \mu_{h} + d + \gamma\) (31)

\(S_{h}^{0} = \frac{\Lambda_{h}\left( \mu_{h} + \psi \right)}{\mu_{h}\left( \mu_{h} + \psi + \omega \right)},\quad S_{v}^{0} = L\left( 1 - \frac{\mu_{v}}{r_{v}} \right)\) (32)

Thus

\(FV^{- 1} = \begin{bmatrix} \frac{\beta_{2}(S_{h}^{0} + \sigma V_{h}^{0)}}{a_{1}a_{2}} & \frac{\beta_{2}(S_{h}^{0} + \sigma V_{h}^{0)}}{a_{2}} & \frac{\beta_{2}(S_{h}^{0} + \sigma V_{h}^{0)}}{\delta_{2}k\mu_{2}} & 0 & \frac{\beta_{2}(S_{h}^{0} + \sigma V_{h}^{0)}}{\delta_{2}k} \\ 0 & 0 & 0 & 0 & 0 \\ \frac{S_{v}^{0}\eta\pi_{h}\tau_{1}}{a_{1}a_{2}\delta_{1}k_{1}} & \frac{S_{v}^{0}\eta\tau_{1}}{a_{2}\delta_{1}k_{1}} & \frac{S_{v}^{0}\alpha_{1}\tau_{1}}{\delta_{1}k_{1}\mu_{v}} + \frac{S_{v}^{0}\alpha_{2}\tau_{2}}{\delta_{2}k_{2}\mu_{v}} & \frac{S_{v}^{0}\tau_{1}}{\delta_{1}k_{1}} & \frac{S_{v}^{0}\tau_{2}}{\delta_{2}k_{2}} \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}\) (33)

The characteristic polynomial of \(FV^{- 1}\) result to

\(\lambda^{2}(\lambda^{2} - (A + D)\lambda + (AD - BC)) = 0\) (34)

Which gives

\(\lambda_{1} = 0,\quad\lambda_{2} = 0,\quad\lambda^{2} - (A + D)\lambda + (AD - BC) = 0\)

\(R_{0} = \rho\left( FV^{- 1} \right) = \frac{1}{2}((A + D) + \sqrt{(A - D)^{2} + 4BC)}\) (35)

Where

\[A = \frac{\beta_{2}\pi_{h}\Lambda_{h}\left\lbrack (\mu_{h} + \psi) + \sigma\omega \right\rbrack}{a_{1}a_{2}\mu_{h}(\mu_{h} + \omega + \psi)},\quad B = \frac{\beta_{1}\Lambda_{h}\left\lbrack (\mu_{h} + \psi) + \sigma\omega \right\rbrack}{\delta_{2}k\mu_{v}\mu_{h}(\mu_{h} + \omega + \psi)},\quad C = \frac{L\left( 1 - \frac{\mu_{v}}{r_{v}} \right)\eta\pi_{h}\tau_{1}}{a_{1}a_{2}\delta_{1}k_{1}},\quad D = \frac{L\left( 1 - \frac{\mu_{v}}{r_{v}} \right)\alpha_{1}\tau_{1}}{a_{1}k_{1}\mu_{v}} + \frac{L\left( 1 - \frac{\mu_{v}}{r_{v}} \right)\alpha_{2}\tau_{2}}{\delta_{2}k_{2}\mu_{v}}\]

(36)

Following the theorem of van den Driessche & Watmough (2002), we state the following result.

Theorem 5.1. The disease-free equilibrium E₀ of the cholera model (3) is locally asymptotically stable if R₀ < 1 and unstable if R₀ > 1.

The DFE has a strong influence on the behavior of the disease transmission dynamics in the population. This implies that if we are looking into eliminating Cholera from the population, we have to establish the conditions necessary for the cholera-free equilibrium to be stable. Epistemologically, when R₀ < 1, on average each infected individual infects fewer than one individual and the disease dies out. On the other hand, if R₀ > 1, on average each infected individual infects more than one other individual and we would expect the disease to spread. However, this depends on the initial sizes of infected persons.

Backward Bifurcation

The significance of examining backward bifurcation is that the classical requirement of R₀ < 1 is necessary but not sufficient for effective disease elimination. Consequently, the dynamics of the transmission of Cholera with backward bifurcation is difficult to control. This phenomenon has been reported in several epidemiological contexts; see, for example, (Ishaku and Hussain 2024). Therefore, it is always important to examine the possibility of backward bifurcation before exploring the global stability. We then examine the presence of backward bifurcation in accordance with (Ishaku & Hussaini, 2020;) by applying Center Manifold theory. The outcome is shown as follows: let S_h = x₁, E_h = x₂, V_h = x₃, I_h = x₄, R_h = x₅, S_v = x₆, E_v = x₇, C = x₈, P = x₉ and the total population is given by N = x₁ + x₂ + x₃ + x₄ + x₅ + x₆ + x₇ + x₈ + x₉ Furthermore, using the vector representation X = (x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉), the cholera dynamical system (3) will be re-written as

\(\frac{dX}{dt} = Y(X)\) (37)

\[{\frac{dx_{1}}{dt} = \Lambda_{h} - \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} - (\mu_{h} + \omega)x_{1} + \psi x_{3} }{\frac{dx_{2}}{dt} = \left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{1} + \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} + \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)R_{h} - (\mu_{h} + \pi_{h})E_{h} }{\frac{dx_{3}}{dt} = \omega x_{1} - (\mu_{h} + \psi)x_{3} - \sigma\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{3} }{\frac{dx_{4}}{dt} = \pi_{h}x_{2} - (\mu_{h} + d + \gamma)x_{4} }{\frac{dx_{5}}{dt} = \gamma x_{4} - \varepsilon\left( \frac{\beta_{}x_{9}}{k + x_{9}} + \beta_{2}x_{4} \right)x_{5} - \mu_{h}x_{5} }{\frac{dS_{v}}{dt} = r_{v}(x_{6} + x_{7})\left( 1 - \frac{(x_{6} + x_{7})}{L} \right) - \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6} - \mu_{v}x_{6} }{\frac{dx_{7}}{dt} = \left( \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}} \right)x_{6} - \mu_{v}x_{7} }{\frac{dx_{8}}{dt} = \eta x_{4} + \alpha_{1}x_{7} - \delta_{1}x_{8} }\]\(\frac{dx_{9}}{dt} = \alpha_{2}x_{7} - \delta_{2}x_{9}\) (38)

Where

\[N_{v} = x_{6} + x_{7}\]

\[\lambda_{h} = \frac{\beta_{1}x_{9}}{k + x_{9}} + \beta_{2}x_{4}\]

\(\lambda_{v} = \frac{\tau_{1}x_{8}}{k_{1} + x_{8}} + \frac{\tau_{2}x_{9}}{k_{2} + x_{9}}\),

consider the case when R₀ = 1. Suppose, that β₂ is chosen as a bifurcation parameter, implies β₂ = β^c Given a system of ODEs The Jacobian matrix is defined as

\[J(x) = \left( \begin{matrix} - \lambda - (\mu_{h} + \psi) & 0 & \psi & - x_{1}\beta_{2} & 0 & 0 & 0 & 0 & - x_{1}\beta_{1}\frac{k}{(k + x_{9})^{2}} \\ \lambda_{h} & - (\mu_{h} + \pi_{h}) & \sigma\lambda_{h} & \beta_{2}(x_{1} + \sigma x_{3} + \varepsilon x_{5}) & \varepsilon\lambda_{h} & 0 & 0 & 0 & \beta_{1}\frac{k}{(k + x_{9})^{2}}(x_{1} + \sigma x_{3} + \varepsilon x_{5}) \\ \omega & 0 & - (\mu_{h} + \psi) - \sigma\lambda_{h} & - \sigma x_{3}\beta_{2} & 0 & 0 & 0 & 0 & - \sigma x_{3}\beta_{1}\frac{k}{(k + x_{9})^{2}} \\ 0 & \pi_{h} & 0 & - (\mu_{h} + d + \gamma) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \gamma - \varepsilon x_{5}\beta_{2} & - (\varepsilon\lambda_{h} + \mu_{h}) & 0 & 0 & 0 & - \varepsilon x_{5}\beta_{1}\frac{k}{(k + x_{9})^{2}} \\ 0 & 0 & 0 & 0 & 0 & r_{v}(1 - \frac{2N_{v}}{L}) - \lambda_{v} - \mu_{v} & r_{v}(1 - \frac{2N_{v}}{L}) & - x_{6}\tau_{1}\frac{k_{1}}{(k + x_{9})^{2}} & - x_{6}\tau_{2}\frac{k_{2}}{(k + x_{9})^{2}} \\ 0 & 0 & 0 & 0 & 0 & \lambda_{v} & - \mu_{h} & - x_{6}\tau_{1}\frac{k_{1}}{(k + x_{9})^{2}} & x_{6}\tau_{2}\frac{k_{2}}{(k + x_{9})^{2}} \\ 0 & 0 & 0 & \eta & 0 & \alpha_{1} & 0 & - \delta_{1} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \alpha_{2} & - \delta_{1} \end{matrix} - \right)\]

The Jacobian of f = (f₁, f₂, f₃, f₄, f₅, f₆, f₇, f₈, f₉)^T evaluated around disease free equilibrium (E_o) denoted by J(E_o) is given by

\[J(E_{0})\beta^{c} = \begin{pmatrix} - (\mu_{h} + \psi) & 0 & \psi & - x_{1}\beta_{2} & 0 & 0 & 0 & 0 & - x_{1}\beta_{1}\frac{k}{(k)^{2}} \\ 0 & - (\mu_{h} + \pi_{h}) & 0 & \beta_{2}(x_{1} + \sigma x_{3}) & 0 & 0 & 0 & 0 & \beta_{1}\frac{k}{(k)^{2}}(x_{1} + \sigma x_{3}) \\ \omega & 0 & - (\mu_{h} + \psi) & - \sigma x_{3}\beta_{2} & 0 & 0 & 0 & 0 & - \sigma x_{3}\beta_{1}\frac{k}{(k)^{2}} \\ 0 & \pi_{h} & 0 & - (\mu_{h} + d + \gamma) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \gamma - \varepsilon x_{5}\beta_{2} & - \mu_{h} & 0 & 0 & 0 & - \varepsilon x_{5}\beta_{1}\frac{k}{(k)^{2}} \\ 0 & 0 & 0 & 0 & 0 & r_{v}(1 - \frac{2N_{v}}{L}) - \lambda_{v} - \mu_{v} & r_{v}(1 - \frac{2N_{v}}{L}) & - x_{6}\tau_{1}\frac{k_{1}}{(k_{1})^{2}} & - x_{6}\tau_{2}\frac{k_{2}}{(k)^{2}} \\ 0 & 0 & 0 & 0 & 0 & \lambda_{v} & - \mu_{v} & - x_{6}\tau_{1}\frac{k_{1}}{(k_{1})^{2}} & x_{6}\tau_{2}\frac{k_{2}}{(k)^{2}} \\ 0 & 0 & 0 & \eta & 0 & \alpha_{1} & 0 & - \delta_{1} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \alpha_{2} & - \delta_{1} \end{pmatrix}\]

The Jacobian J(E₀)β^c of the linearized system has a simple zero eigenvalue (meaning that the remaining eigenvalues have negative real part). Hence, by center manifold theorem developed in (Castillo-Chavez, Feng, & Huang, 2002), it is possible to analyze the dynamics of the system within β₂ = β^c. Eigenvectors of J(E₀)β^c: Let w = [w₁, w₂, w₃, w₆, w₈, w₉]^T and v = [u₁, u₂, u₃, u₄, u₅, u₆, u₇, u₈, u₉]^T be the right and left eigenvectors, respectively, associated with the zero eigenvalue of the Jacobian of the linearized system J(E0)β^c, defined by:

\(w_{1} = \frac{a_{3}w_{3}}{\psi}\) > 0, \(w_{2} = \frac{\pi_{h}w_{4}}{a_{2}}\) > 0, \(w_{3} = \frac{\psi w_{1}}{a_{3}}\) > 0, \(w_{4} = \frac{a_{4}w_{1} + a_{5}w - a_{6}w_{3} + \eta w_{8}}{a_{3}}\) < 0,

\(w_{5} = 0\), \(w_{6} = \frac{\mu_{v}w_{7}}{a_{10}}\) > 0, \(w_{7} = 0,\) \(w_{8} = \frac{a_{4}w_{1} - a_{5}w_{2} + a_{6}w_{3} + a_{7}w_{4} - \gamma w_{5}}{\mu}\) > 0,

\(w_{9} = - (\frac{\delta_{1}w_{8} + a_{11}w_{6} + a_{12}w_{7}}{a_{2}})\) < 0, (39)

And

\(u_{1} = \frac{- (a_{3}a_{4} + \psi a_{6})u_{4} + (a_{3}a_{12} + \psi a_{14})u_{9}}{a_{1}a_{3} - \psi\omega}\) < 0, \(u_{2} = \frac{a_{7}u_{4}}{\pi_{h}}\) > 0

\(u_{3} = \frac{- (\omega a_{4} + a_{1}a_{6})u_{4} + (\omega a_{12} + a_{1}a_{14})u_{9}}{a_{1}a_{3} - \psi\omega}\) < 0, \(u_{4} = \frac{\pi_{h}u_{2}}{a_{7}}\) > 0, \(u_{5} = \frac{- \gamma u_{4}}{a_{8}}\) < 0,

\(u_{6} = \frac{- (a_{10}u_{7} + a_{11}u_{8} + a_{15}u_{9})}{a_{9}}\) < 0, \(u_{7} = \frac{a_{12}u_{8} + a_{16}u_{9}}{\mu_{v}}\) > 0, \(u_{8} = \frac{\delta_{2}u_{9}}{\alpha_{2}}\) > 0

\(u_{9} = \frac{\alpha_{2}u_{8}}{\delta_{2}}\) > 0, (40)

Following the Castillo-Chavez and Song theorem (Castillo-Chavez et al., 2002), we compute the associated backward bifurcation parameters A and B that are responsible for the direction of bifurcation,

\(A = \sum_{k,i,j = 1}^{n}{u_{k}w_{i}w_{j}\frac{\partial^{2}f_{k}(0,0)}{\partial x_{i}\partial x_{j}}},\quad\quad\quad B = \sum_{k,i, = 1}^{n}{u_{k}w_{i}\frac{\partial^{2}f_{k}(0,0)}{\partial x_{i}\partial\beta^{c}}},\) (41)

After computing the non-zero second-order partial derivatives of f as defined in (3.16), and evaluating around diseases free E0, we obtain

\[{A = \frac{2\beta_{1}}{k^{2}}(S_{h}^{0}(u_{1}w_{9}^{2} - u_{9}^{2}) + \sigma V_{h}^{0}(u_{3}w_{9}^{2} - u_{9}^{2})) + \frac{2\beta_{1}}{k}\left( (u_{2} - u_{1})w_{1}w_{9} + \sigma(u_{2} - u_{3})w_{3}w_{9} \right) }{\quad + 2\beta_{2}\left( (u_{2} - u_{1})w_{1}w_{4} + \sigma(u_{2} - u_{3})w_{3}w_{4} \right) - 2\left( u_{6}w_{6}^{2} + u_{7}w_{6}w_{7} \right)\left( \tau_{2} + \frac{\tau_{1}}{L_{0}} \right) }\]\(\quad + \frac{2}{k_{1}^{2}}\left( \tau_{1}S_{v}^{0}u_{6}w_{8} - \frac{\tau_{1}}{k_{2}}S_{v}^{0}u_{7}w_{8} \right) + \frac{2}{k_{2}^{2}}\left( \tau_{2}S_{v}^{0}u_{6}w_{8} - \frac{\tau_{2}}{k_{1}}S_{v}^{0}u_{7}w_{9} \right) + \frac{2u_{7}w_{8}\tau_{1}}{k_{2}} + \frac{2v_{7}w_{6}w_{9}\tau_{2}}{k_{1}k_{2}}\) (42)

And

\(B = \frac{w_{9}(u_{2} - u_{3})\sigma V_{h}^{0} + (u_{2} - u_{1})S_{h}^{0}}{k}\) (43)

which indicates that B > 0. Hence, backward bifurcation exists only if A > 0.

Theorem 5.2. The disease-free equilibrium E₀ of model (3) is globally asymptotically stable if R₀ < 1 and unstable if R₀ > 1.

The global stability is proved for a special case β₂ = 0 and ϵ = 0. Thus, the modified model is given as

\[{\frac{dS_{h}}{dt} = \Lambda_{h} - {\widehat{\lambda}}_{h}S_{h} - (\mu_{h} + \omega)S_{h} + \psi V_{h} }{\frac{dE_{h}}{dt} = {\widehat{\lambda}}_{h}S_{h} + \sigma{\widehat{\lambda}}_{h}V_{h} - (\mu_{h} + \pi_{h})E_{h} }{\frac{dV_{h}}{dt} = \omega S_{h} - (\mu_{h} + \psi)V_{h} - \sigma{\widehat{\lambda}}_{h}V_{h} }{\frac{dI_{h}}{dt} = \pi_{h}E_{h} - (\mu_{h} + d + \gamma)I_{h} }{\frac{dR_{h}}{dt} = \gamma I_{h} - \mu_{h}R_{h} }{\frac{dS_{v}}{dt} = r_{v}N_{v}\left( 1 - \frac{N_{v}}{L} \right) - \lambda_{v}S_{v} - \mu_{v}S_{v} }{\frac{dE_{v}}{dt} = \lambda_{v}S_{v} - \mu_{v}E_{v} }{\frac{dC}{dt} = \eta I_{h} + \alpha_{1}E_{v} - \delta_{1}C }\]\(\frac{dP}{dt} = \alpha_{2}E_{v} - \delta_{2}P\) (44)

Where

\({\widehat{\lambda}}_{h} = \frac{\beta_{1}}{k + P}\) (45)

The associated reproduction number for model (51) is presented as:

\({\widehat{R}}_{0} = \rho(FV^{- 1}) = \frac{1}{2}\left( D + \sqrt{D^{2} + 4BC} \right)\) (46)

Proof. Following the approach of (Castillo-Chavez et al., 2002), let X = (S_h, V_h, R_h, S_v)^T be non-infected classes and Z = (E_h, I_h, E_v, C, P )^T be infected classes. Thus, our model can be written in the form

\[\frac{dX}{dt} = F(X,Z)\]

\[\frac{dZ}{dt} = G(X,Z)\]

Where

\[{F(X,Z) = (\Lambda_{h} - {\widehat{\lambda}}_{h}S_{h} - (\mu_{h} + \omega)S_{h} + \psi V_{h},\mspace{6mu}\omega S_{h} - (\mu_{h} + \psi)V_{h},\mspace{6mu}\gamma I_{h} - \mu_{h}R_{h} }{\quad\quad\quad\quad r_{v}N_{v}\left( 1 - \frac{N_{v}}{L} \right) - \lambda_{v}S_{v} - \mu_{v}S_{v}) }{G(X,Z) = ({\widehat{\lambda}}_{h}S_{h} + \sigma{\widehat{\lambda}}_{h}V_{h} - (\mu_{h} + \pi_{h})E_{h},\mspace{6mu}\pi_{h}E_{h} - (\mu_{h} + d + \gamma)I_{h}, }\]\(\quad\quad\quad\quad\lambda_{v}S_{v} - \mu_{v}E_{v},\mspace{6mu}\eta I_{h} + \alpha_{1}E_{v} - \delta_{1}C,\mspace{6mu}\alpha_{2}E_{v} - \delta_{2}P)\) (47)

Condition (i). Set Z = 0 and consider the uninfected subsystem \(\frac{dX}{dt} = F(X,0)\) we obtain

\[{\frac{dS_{h}}{dt} = \Lambda_{h} - \mu_{h}S_{h} }{\frac{dV_{h}}{dt} = \omega S_{h} - (\mu_{h} + \psi)V_{h} }\]\(\frac{dS_{v}}{dt} = r_{v}N_{v}\left( 1 - \frac{N_{v}}{L} \right) - \mu_{h}S_{v}\) (48)

Following (Castillo-Chavez et al., 2002) the disease-free equilibrium point E₀ is a globally asymptotically stable (G.A.S) equilibrium of (28) provided that R₀ < 1 if the following conditions are satisfied.

• Condition i: For ^dX = F (X, 0), the equilibrium point X is globally asymptotically stable.

dt

• Condition ii: Gˆ(X, Z) = AZ − G(X, Z), where Gˆ(X, Z) ≥ 0 for all X, Z ∈ Ω.

To show Condition i we solve system (48) converges to their respective points in the DFE, E₀, as t → ∞. The equation (48) can be obtain as

\(S_{h}(t) \rightarrow S_{h}^{0} = \frac{\Lambda_{h}}{\mu_{h}},\mspace{6mu}\quad V_{h}(t) \rightarrow V_{h}^{0}(t) = \frac{\omega\Lambda_{h}}{\mu_{h}(\mu_{h} + \psi)},\quad S_{v}(t) \rightarrow S_{v}^{0}(t) = L\) (49)

Thus, the uninfected subsystem has a globally asymptotically stable equilibrium \((S_{h}^{0},V_{h}^{0},R_{h}^{0},S_{h}^{0})\)

the biologically feasible region. Hence condition (i) of the Castillo–Chavez framework (Castillo-

Chavez et al., 2002) is satisfied.

Condition (ii). We write the infected subsystem in the form

\(\widehat{G}(X,Z) = AZ - G(X,Z),\) (50)

\(G(X,Z) = \begin{pmatrix} {\widehat{\lambda}}_{h}S_{h} + \sigma{\widehat{\lambda}}_{h}V_{h} - (\mu_{h} + \pi_{h})E_{h} \\ \pi_{h}E_{h} - (\mu_{h} + d + \gamma)I_{h} \\ \lambda_{v}S_{v} - \mu_{v}E_{v} \\ \eta I_{h} + \alpha_{1}E_{v} - \delta_{1}C \\ \alpha_{2}E_{v} - \delta_{2}P \end{pmatrix}\) (51)

\[A = \begin{pmatrix} - (\mu_{h} + \pi_{h}) & 0 & 0 & 0 & (S_{h}^{0} + \sigma V_{h}^{0})\frac{\beta_{1}}{k} \\ \pi_{h} & - (\mu_{h} + d + \gamma) & 0 & 0 & 0 \\ 0 & 0 & - \mu_{v} & S_{v}^{0}\frac{\tau_{1}}{k_{1}} & S_{v}^{0}\frac{\tau_{2}}{k_{2}} \\ 0 & \eta & \alpha_{1} & - \delta_{1} & 0 \\ 0 & 0 & \alpha_{2} & 0 & \delta_{2} \end{pmatrix}\]

(52)

\[\widehat{G}(X,Z) = AZ - G(X,Z),\]

\(\widehat{G}(X,Z) = \begin{pmatrix} (S_{h}^{0} + \sigma V_{h}^{0})\frac{\beta_{1}}{k}P - {\widehat{\lambda}}_{h}S_{h} - \sigma{\widehat{\lambda}}_{h}V_{h} \\ 0 \\ S_{v}^{0}\frac{\tau_{1}}{k_{1}}C + S_{h}^{0}\frac{\tau_{2}}{k_{2}}P - \widehat{\lambda}S_{v} \\ \begin{aligned} & 0 \\ & 0 \\ & 0 \end{aligned} \end{pmatrix}\) (53)

\(\widehat{G}(X,Z) = \begin{pmatrix} (\frac{(S_{h}^{0} - S_{h})k + S_{h}^{0}P}{k(k + P)})\beta P + (\frac{(V_{h}^{0} - V_{h})k + V_{h}^{0}P}{k(k + P)})\sigma\beta_{1}P \\ 0 \\ (\frac{(S_{v}^{0} - S_{v})k_{1} + S_{h}^{0}P}{k_{1}(k_{1} + C)})\tau_{1}C + (\frac{(S_{v}^{0} - S_{v})k_{2} + S_{v}^{0}P}{k_{2}(k_{2} + P)})\tau_{2}P \\ \begin{aligned} & 0 \\ & 0 \\ & 0 \end{aligned} \end{pmatrix}\) (54)

Since \(0 \leq S_{h}(t) \leq \frac{\Lambda_{h}}{\mu_{h}},\quad 0 \leq V_{h}(t) \leq \frac{\omega\Lambda_{h}}{\mu_{h}(\mu_{h} + \psi)}\)and \(0 \leq S_{v}(t) \leq L,\) then it follows that \(\widehat{G}(X,Z) \geq 0.\)

Thus \(E_{0}\) is GAS whenever \({\widehat{R}}_{0}\) < 1.

Validation

The model was validated by comparing simulated results with observed cholera cases data obtained from Nigeria Centre for Disease Control (NCDC) between the period of 2022–2024. The model was compared based on the monthly reported cholera cases over three years which is 36-month period (from January 2022 - December 2024), as presented in Figure 2 and Ta

\[{RMSE = \sqrt{\frac{1}{n}\sum_{i = 1}^{n}{(y_{i} - {\widehat{y}}_{i})^{2}}} }{R^{2} = 1 - \frac{\sum_{}^{}{(y_{i} - {\widehat{y}}_{i})^{2}}}{\sum_{}^{}{(y_{i} - \bar{y})^{2}}} }\] (55)

Figure 2 illustrates the comparison between observed data and estimated number of infected individuals predicted by the model (I_h) over time (monthly).

Figure 2: Observed data obtain from (NCDC) fitted with a model solution

Parameter Estimation

Model parameters were obtained from published literature and, where necessary, estimated by fitting the model to available cholera cases data . The fitting was performed using a least squares approach to minimize the difference between observed and simulated data. Some parameters were assumed due to limited data availability, but their influence on model outcomes was examined through sensitivity analysis.

The table below presents the parameter description and values for cholera model. Table 1: Parameter Description and Values for Cholera Model.

Sensitivity Analysis of the Model

In this section, Partial Rank Correlation Coefficient (PRCC) is used to study the sensitivity of the model output due to changes in some specific parameter values. The effective reproduction number R₀ is used as the response function for the PRCC. If a parameter falls on the positive side, then R₀ will increase as the parameter increases. However, if the sensitivity index of a parameter is negative, then R₀ decreases as the parameter increases. The results of the (PRCC) illustrate the degree of effect of each parameter that most strongly influences the basic reproduction number R₀ in the extended model of transmission of Cholera. Generally, the PRCC analysis identifies the key parameters driving cholera transmission in the model. The strongest contributors to an increase in R₀ include Λ_h, β₂, τ₂, π_h, and L, while µ_h, µ_v, γ, and δ₂ are the most important parameters reducing R₀. This information is valuable for guiding targeted control strategies aimed at reducing disease spread

Figure 3: Parameter affecting reproduction number

Numerical Simulation

In this section, we perform a numerical simulation of model system Equations (3) to demonstrate our simulation results and to illustrate the asymptotic behaviour of the model. The systems of differential equations were solved over a specific period of time period using Range- Kutta method embedded in MATLAB with ode45. The parameter values used in the simulations are given in Table 1

1. Effect of Re-infection

Figure 4: Infection Dynamic with and without Re-infection

The above Figure 4, illustrates the dynamics of infected individuals over a period of 120 months under two scenarios, with and without re-infection. ε = 0 (Without re-infection) Recovered indi-viduals acquire permanent immunity. Consequently, the number of infected individuals declines over time, leading to possible disease elimination. ε ̸= 0 (With re-infection) Immunity wanes and recovered individuals can be re-infected. This leads to a higher and more persistent number of infected individuals in the population. Re-infection enhances disease persistence and increases the burden of Cholera, emphasizing the importance of considering waning immunity in transmission modeling. These results indicate the critical role of re-infection in maintaining cholera endemic. Models that did not consider re-infection may underestimate the long-term disease burden, whereas incorporating re-infection provides more realistic description of cholera dynamics, particularly in regions where immunity is temporary. Consequently, sustained control strategies such as continu-ous vaccination, improved sanitation, and access to clean water are essential for effective disease management.

2. Assessing the role of Vaccine with σ = 0

Figure 5: Vaccine with σ = 0

Figure 5 illustrates the effect of varying vaccination rates on the number of Vacinated individuals over time. In the absence of vaccination, the number of infected individuals increases rapidly and stabilizes at a high endemic level, indicating persistent disease transmission. As the vaccination rate increases to 25% and 50%, a noticeable reduction in the peak and endemic level of infections is observed, although the disease remains present in the population. At higher vaccination rates (75%), the infected population declines significantly, approaching disease elimination. When the vaccination rate reaches 100%, the number of infected individuals decreases to zero, indicating com-plete eradication of the disease. These results demonstrate that vaccination substantially reduces cholera transmission and that high vaccination coverage is essential for effective disease control and elimination.

3. Assessing the role of Vaccine with \(\mathbf{\sigma \neq 0}\)

Figure 6: Vaccine with \(\sigma \neq 0\)

Figure 6 illustrates the dynamics of vaccinated individuals for different vaccination rates over time. When the vaccination rate is zero, the vaccinated population remains at zero throughout the simulation, showing the absence of vaccination. when the vaccination rate increases to 25% and 50%, the vaccinated population increases gradually and approaches higher steady-state levels. At higher vaccination rates (75% and 100%), the number of vaccinated individuals increases rapidly and stabilizes at significantly larger values. These results show that increasing vaccination coverage leads to a higher proportion of vaccinated individuals in the population, thereby reducing susceptibility to infection. This supports the effectiveness of vaccination as a key intervention strategy for controlling and eliminating Cholera.

Vaccination Rate Density Plot

Figure 7: Vaccine Coverage and effective control

The figure 7 shows that the color bar indicate that the darker blue color indicate to lower density levels, whereas brighter yellow colors indicate higher density levels. The Figure 7 illustrates that vaccination is most effective when both vaccine coverage and vaccine effectiveness are high. As the vaccine failure rate increases, the disease burden rises sharply, even at high levels of vaccine coverage. In contrast, when the vaccine failure rate is low, increasing vaccine coverage leads to a substantial reduction in disease transmission. The highest disease levels occur when vaccine coverage is low and the failure rate is high. Overall, the figure emphasizes that effective disease control depends not only on vaccinating a large proportion of the population but also on ensuring that the vaccine provides strong and reliable protection.

Discussion

The research developed an extended deterministic compartmental model to investigate the trans-mission dynamics of Cholera by incorporating three critical components: a vaccinated human popu-lation, an exposed human compartment, and vector recruitment, which was modelled using logistic growth, which contributes to environmental contamination. The human host is divided into five compartments: susceptible, vaccinated, exposed, infected, and recovered. The vector population is divided into two Compartment: susceptible, exposed vectors and the environmental reservoir divide the environment into safe and unsafe environments according to the concentration of Vibrio cholera. We established that the model is epidemiologically feasible and well-posed and we also showed the existence of the disease-free equilibrium. Furthermore, we employed the next generation matrix technique to derive the reproduction number R₀. We proved that the model has two equilibrium points; the disease-free equilibrium which is locally asymptotically stable whenever R₀ < 1, unsta-ble otherwise. We performed the sensitivity analysis on the reproductive number, R₀. Our analyses revealed that the parameters human recruitment rate Λ_h introduces more susceptible individuals into the population, thereby enhancing transmission potential. Similarly, the human-to-human transmission rate β₂, the high-level environmental shedding rate τ₂, and the progression rate from exposed to infectious individuals π_h exhibit strong positive correlations with R₀. Contrarily, the natural death rate of humans µ_h shows a strong negative correlation with R₀. Higher values of µ_h decrease the number of individuals available for infection, thereby lowering transmission. Other parameters with notable negative influence include the vector mortality rate µ_v, the human recov-ery rate γ, and the environmental pathogen decay rate δ₂. numerical simulations where carried out to confirm the theoretical analysis and explored more patterns of dynamical behaviours of our model. The analytical results showed that the basic reproduction number depends on the combined influence of human-to-environment, human-to-human, and vector-mediated transmission pathways. The model demonstrates that increasing vaccination coverage and improving vaccine effectiveness reduce cholera transmission, but these measures alone may not guarantee disease elimination due to the persistence of environmental reservoirs and reinfection. Numerical simulations further con-firmed that integrated intervention strategies, including vaccination, improved sanitation, vector control, and early treatment, are essential for achieving effective and sustainable cholera control.

Conclusion

The results of the model highlight the significant role of environmental contamination and vacci-nation in cholera transmission dynamics. The inclusion of an exposed compartment and vector logistic growth provides a more realistic disease progression compared to existing models.

The findings are consistent with previous studies based on SIRB and SEIR frameworks, which emphasize the importance of environmental reservoirs in sustaining cholera outbreaks.

However, the model has some limitations such as:

1. seasonal variations such as rainfall, temperature, and environmental conditions that influence cholera outbreaks. Seasonal effects can significantly affect both vector population dynamics and environmental contamination levels.

2. stochastic modeling approaches could be explored to account for random variations and un-certainties in disease transmission, particularly in small populations or during the early stages of an outbreak.

3. optimal control theory to determine the cost effectiveness intervention strategies, including vaccination, treatment, sanitation improvement, and vector control.

4. fractional-order differential equations to capture memory effects in disease transmission pro-cesses, which may provide a more accurate description of cholera dynamics.

Despite these limitations, the model provides useful insights into effective intervention strategies, particularly the combined use of vaccination and environmental sanitation.

References

Al-Shanfari, S., ELmojtaba, I. M., & Alsalti, N. (2019). The role of houseflies in cholera transmission. Communications in Mathematical Biology and Neuroscience, Article 31. [Crossref]

Bertuzzo, E., Mari, L., Righetto, L., Gatto, M., Casagrandi, R., Blokesch, M., ... Rinaldo, A. (2011). Prediction of the spatial evolution and effects of control measures for the unfolding Haiti cholera outbreak. Geophysical Research Letters, 38(6). [Crossref]

Castillo-Chavez, C., Feng, Z., & Huang, W. (2002). On the computation of R₀ and its role on global stability. In C. Castillo-Chavez, S. Blower, P. van den Driessche, & et al. (Eds.), Mathematical approaches for emerging and reemerging infectious diseases: Models, methods and theory (pp. 229-250). Springer. [Crossref]

CDC. (2019). Cholera - Vibrio cholerae infection. [Link]

CDC. (2021). Premiumtimes. [Link]

CDC. (2024). NCDC. [Link]

Das, P., Mukherjee, D., & Sarkar, A. K. (2005). Study of a carrier dependent infectious disease - cholera. Journal of Biological Systems, 13(2), 233-244. [Crossref]

Edward, S., & Nyerere, N. (2015). A mathematical model for the dynamics of cholera with control measures. Applied and Computational Mathematics, 4(2), 53-63. [Crossref]

Fister, K. R., Gaff, H., Lenhart, S., & Schaefer, E. (2016). Optimal control of vaccination in an age-structured cholera model. In G. Chowell & J. M. Hyman (Eds.), Mathematical and statistical modeling for emerging and re-emerging infectious diseases (pp. 221-248). Springer. [Crossref]

Fotedar, R., Banerjee, U., Samantray, J. C., & Shriniwas. (1992). Vector potential of hospital houseflies with special reference to Klebsiella species. Epidemiology and Infection, 109(1), 143-147.

Hai-Feng, H., & Na-Na, S. (2012). Global stability for a binge drinking model with two stages. Discrete Dynamics in Nature and Society, 1-15. [Crossref]

Harris, J. B., LaRocque, R. C., Kadri, F., Ryan, E. T., & Calderwood, S. B. (2012). Cholera. The Lancet, 379(9835), 2466-2476. [Crossref]

Hartley, D. M., Morris Jr., J. G., & Smith, D. L. (2006). Hyperinfectivity: A critical element in the ability of V. cholerae to cause epidemics? PLoS Medicine, 3(1), e7. [Crossref]

Hove-Musekwa, S. D., Nyabadza, F., Chiyaka, C., Das, P., Tripathi, A., & Mukandavire, Z. (2011). Modelling and analysis of the effects of malnutrition in the spread of cholera. Mathematical and Computer Modelling, 53(9-10), 1583-1595. [Crossref]

Isere, A., Osemwenkhae, J., & Okuonghae, D. (2014). Optimal control model for the outbreak of cholera in Nigeria. African Journal of Mathematics and Computer Science Research, 7(2), 24-30. [Crossref]

Ishaku, A., & Hussaini, N. (2020). Modelling effects of HCV plus-strand RNA influx into a cell during HCV replication. Journal of the Indonesian Mathematical Society, 26(1), 37-54. [Crossref]

Kaka, R., & Anteneh, L. M. (2024). Mathematical modelling and analysis of cholera dynamics via vector transmission. Communications in Mathematical Biology and Neuroscience, Article 74. [Crossref]

Keiding, J. (1986). The house-fly: Biology and control (Technical Report). World Health Organization.

Lata, K., Mishra, S. N., & Misra, A. K. (2020). An optimal control problem for carrier dependent diseases. Biosystems, 187, 104039. [Crossref]

Levine, O. S., & Levine, M. M. (1991). Houseflies (Musca domestica) as mechanical vectors of shigellosis. Clinical Infectious Diseases, 13(4), 688-696. [Crossref]

Misra, A. K., Mishra, S. N., Pathak, A. L., & et al. (2013). A mathematical model for the control of carrier-dependent infectious diseases with direct transmission and time delay. Chaos, Solitons and Fractals, 57, 41-53. [Crossref]

Modnak, C. (2017). A model of cholera transmission with hyperinfectivity and its optimal vaccination control. International Journal of Biomathematics, 10(6), 1750084. [Crossref]

Mukandavire, Z., Liao, S., Wang, J., Gaff, H., Smith, D. L., & Morris Jr., J. G. (2011). Estimating the reproductive numbers for the 2008-2009 cholera outbreaks in Zimbabwe. Proceedings of the National Academy of Sciences, 108(21), 8767-8772. [Crossref]

Namawejje, H., Obuya, E., & Luboobi, L. S. (2018). Modeling optimal control of cholera disease under the interventions of vaccination, treatment and education awareness. Journal of Mathematics Research, 10(5), 137-152. [Crossref]

Nelson, E. J., Harris, J. B., Morris, J. G., Calderwood, S. B., & Camilli, A. (2009). Cholera transmission: The host, pathogen and bacteriophage dynamic. Nature Reviews Microbiology, 7(10), 693-702. [Crossref]

Omondi, E. O. (2016). A mathematical model for onchocerciasis and its treatment with ivermectin [Master's thesis, University of Stellenbosch].

Onuorah, M. O., Atiku, F., & Juuko, H. (2022). Mathematical model for prevention and control of cholera transmission in a variable population. Research in Mathematics, 9(1), 2018779. [Crossref]

Pascual, M., Koelle, K., & Dobson, A. P. (2006). Hyperinfectivity in cholera: A new mechanism for an old epidemiological model? PLoS Medicine, 3(6), e280. [Crossref]

van den Driessche, P., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1-2), 29-48. [Crossref]

WHO. (2011). Weekly epidemiological record, 2011, vol. 86, 36 [full issue]. Weekly Epidemiological Record, 86(36), 389-400.

Yang, C., Wang, X., Gao, D., & Wang, J. (2017). Impact of awareness programs on cholera dynamics: Two modeling approaches. Bulletin of Mathematical Biology, 79, 2109-2131. [Crossref]