Optimal Control Analysis of Cholera Dynamics in the Presence of Asymptotic Transmission

Abstract

Many mathematical models have explored the dynamics of cholera but none have been used to predict the optimal strategies of the three control interventions (the use of hygiene promotion and social mobilization; the use of treatment by drug/oral re-hydration solution; and the use of safe water, hygiene, and sanitation). The goal here is to develop (deterministic and stochastic) mathematical models of cholera transmission and control dynamics, with the aim of investigating the effect of the three control interventions against cholera transmission in order to find optimal control strategies. The reproduction number $R_p$ was obtained through the next generation matrix method and sensitivity and elasticity analysis were performed. The global stability of the equilibrium was obtained using the Lyapunov functional. Optimal control theory was applied to investigate the optimal control strategies for controlling the spread of cholera using the combination of control interventions. The Pontryagin’s maximum principle was used to characterize the optimal levels of combined control interventions. The models were validated using numerical experiments and sensitivity analysis was done. Optimal control theory showed that the combinations of the control intervention influenced disease progression. The characterisation of the optimal levels of the multiple control interventions showed the means for minimizing cholera transmission, mortality, and morbidity in finite time. The numerical experiments showed that there are fluctuations and noise due to its dependence on the corresponding population size and that the optimal control strategies to effectively control cholera transmission, mortality, and morbidity was through the combinations of all three control interventions. The developed models achieved the reduction, control, and/or elimination of cholera through incorporating multiple control interventions.

1. Introduction

Cholera is a life threatening but easily preventable and treatable infection, which is known as an acute diarrhoeal infection caused by the intake of food or water contaminated with the bacterium Vibrio cholerae, which can kill within a short range of time if left untreated [1,2,3,4]. It remains a global threat to public health [5]. Cholera is a major economic and social problem that slows down the growth of the African economy yearly [2]. There are roughly 1.3 to 4.0 million cases and 21,000 to 143,000 deaths worldwide due to cholera every year [1]. The number of cholera cases reported to WHO (World Heath Organization) has continued to increase over the last few years [1,5]. In 2015, about 172,454 cases were reported from 42 countries including 1304 deaths [1,5]. Symptoms of cholera ranges from mild or moderate symptoms, while many develop acute watery diarrhoea with severe dehydration, which may lead to death if left untreated [1,5]. Other studies revealed that it takes between 12 h and 5 days for someone to show symptoms after ingesting contaminated food or water [1,6].

Cholera is generally transmitted through water or food that has been contaminated with faeces from an infected individual [7,8,9,10]. Its transmission is closely related to inadequate access to clean water and sanitation facilities [1]. Transmission occurs in a vulnerable or susceptible person when in contact with an infected person, which is called direct transmission, or by ingesting food contaminated with vibrio cholerae, which is called indirect transmission [8,9,11,12,13].

Cholera is known as an endemic or epidemic disease but there are control interventions recommended by WHO to control its morbidity and reduce death due to cholera [1,2]. It has been recommended by the WHO that a multifaceted approach be used to control morbidity and reduce mortality due to the disease through a combination of surveillance, water sanitation and hygiene, social mobilization, treatment, and oral cholera vaccines [1,2,14].

Mathematical modeling has been a useful tool in answering questions and making effective decisions for policies in epidemiology, public health, and biological sciences [3,15]. There are a number of mathematical models for cholera that have been developed over the years, yet they are still scant [3,15,16,17,18,19,20]. Codeco in [21] proposed a basic model for the dynamics of cholera. His model was extended by Hartley et al. [22] in order to explain the role played by hyper infectious. Capasso and Paveri-Fontana [17] proposed a mathematical model for cholera epidemics occurring in the European Meditteranean region in 1973. In Capasso’s version, two equations describe the dynamics of infected people in the community and the dynamics of the aquatic population of pathogenic bacteria. After that, several cholera models were formulated and analyzed.

In 2001, Codeco extended the cholera model of [17], he presented a basic model for the dynamics of cholera with the inclusion of susceptible individuals in the host population and explored the role of the aquatic reservoir on the persistence of endemic cholera. In 2006, Hartley et al. [22] extended Codeco’s work with an addition of the hyper infectious vibrios. Liao and Wang [10] explored and presented the dynamical analysis that is the stability of equilibria of the system of the deterministic cholera model.

The work in [23] investigated a cholera model with vaccination, and analyzed the local and global asymptotic stability of the disease-free and endemic equilibria of their system. In 2010, [24] presented a water-borne disease model which includes the duel transmission pathways, with bilinear incidence rates employed for both the environmental-to-human and human-to-human infection routes. The did not consider the effect of saturation in their work.

The work of [25] proposed a mathematical model to study how lytic bacteriophage specific for Vibrio cholerae affects cholera outbreaks. Their model explored the vibrios (V), the phage (P), the infection $\left(I_{-}\right)$ solely caused by Vibrio cholerae, and the infection $\left(I_{+}\right)$ caused by both Vibrio cholerae and phage.

The work of [26] also proposed a model to estimate the reproductive number for the 2008–2009 cholera outbreak in Zimbabwe, their model includes both environmental-to-human and human-to-human transmission pathways. Ref. [23] presented a cholera epidemic model with a periodic transmission rate. They defined the reproductive number and showed that the disease-free-equilibrium is globally asymptotically stable and the cholera eventually disappears if the basic reproduction number is greater than one, there exits a positive periodic solution which is globally asymptotically stable. They finally provided numerical simulation to illustrate their analytic results. Ref. [27] presented a cholera transmission dynamic model for public health practitioners. In their work, it was discussed that a basic model (ordinary differential equation models of cholera transmission dynamics) adapted from Codeco (2001) and how it can be modified to incorporate different hypotheses, including the importance of asymptomatic or in apparent infections, and hyper infectious Vibrio cholerae and human-to-human transmission. He highlighted three important challenges of cholera models: (1) Model mis-specification and parameter uncertainty; (2) modeling the impact of water, sanitation, and hygiene intervention; and (3) model structure. He also emphasized that the choice of models should be dictated by the research questions in mind. Fiser et al. [28] developed an age-structured cholera model with optimal control of vaccination using both partial differential equations (PDEs) and ordinary differential equations (ODEs). They carried out the existence of a solution to their non-linear state system using the method of characteristic and a fixed point argument. Numerical results indicate a clear dependence on age in the number of individuals vaccinated, and suggest less vaccination for individuals in young and middle-aged adults, which is expected from the choices of our rates depending on age. Ref. [29] carried out a mathematical modeling study on the impact evaluation of vaccination programs.

Data analysis to forecast the global cholera burden has been carried out. It was found that the impact of vaccination will likely vary depending on local epidemiological conditions including age distribution of cases and the relative contribution of different transmission routes. Sun et al. [30] proposed a mathematical model to describe the transmission of cholera in China. It was found that it may be reasonable to increase the immigration coverage rate and make efforts to improve environmental management, especially for drinking water in order to stop the progress of cholera in China. Hailemariam et al. [31] carried out an analysis on cholera epidemic control using mathematical modeling. They presented a mathematical model for the transmission dynamics of cholera and its preventive measure as a cohort of individuals, along with the Susceptible-Infected-Recovered (SIR) class of individual with Vibrio cholerae concentration. Their results showed that the disease dies out in areas with adequate preventive measures but spreads more and has increased mortality in areas with inadequate preventive measures. Other literatures on the recent development within the framework of biological modeling can be found in ([32,33,34,35]).

The model considered in this paper differs from previous work done and from the work done by [23], in that time dependent control interventions and constant control interventions are incorporated and asymptotic transmission is involved. The aim of this work is to gain insight into how the various control interventions can control the morbidity of cholera and reduce its mortality in the presence of asymptotic transmission in order to guide public health officials to make proper decisions on the control strategies that should be put in place to stop the progression of the spread of cholera and reduce its prevalence in the community. In this paper, a cholera disease transmission mathematical model with asymptotic transmission is formulated using ordinary differential equation and stochastic differential equations and their analysis was also performed. The stochastic differential equation was developed to determine errors in measurement, variability in the population, and other factors that introduce uncertainties. Stability properties and qualitative analysis of the optimal control theory of the cholera model were carried out. The necessary conditions for optimal control of the disease using Pontryagin’s Maximum Principle in order to determine optimal strategies for controlling the spread of the disease and existence of optimal control and state variables are performed. Numerical simulations to illustrate the analytical findings are carried out. The developed model predicted the reduction of cholera in the presence of various scenarios on the choice of the control interventions.The rest of this work is organized as follows. In the next section, the cholera model formulation and its analysis is presented. In Section 3, we present the sensitivity and elasticity analysis, while in Section 4, the stochastic model of the transmission dynamics of cholera is presented. The optimal control problem is formulated and analyzed in Section 5. In Section 6, we discuss the numerical simulations to validate our analytical findings. In Section 7, the numerical simulation of the stochastic model is discussed. The paper ends with discussions of the results and findings in the conclusion in Section 8.

2. Derivation of the Models and Their Analysis

2.1. Deterministic Model of the Transmission Dynamics of the Cholera Model and Its Biological Description

The standard model of the type SIRS (Susceptible-Infected-Recovered-Susceptible) will be used in line with [4,23]. We had earlier discussed in [4] the impact of constant control on cholera transmission. We have in this paper discussed the effect of some time-dependent control intervention in the presence of asymptotic transmission. These time-dependent controls are:

(i)

The use of hygiene promotion and social mobilization, $u_1$;

(ii)

The use of treatment by drug/oral re-hydration solution (ORS), $u_2$;

(iii)

The use of safe water, hygiene, and sanitation, $u_3$.

In our previous work [4], we developed the model of cholera dynamics with constant control intervention but we present below a time-dependent based on some assumption (and we assumed here that (i) the population is a variable one, (ii) people recovered can return to the susceptible group, (iii) people can die due to the disease, (iv) the use of hygiene promotion and social mobilization can be applied to the vulnerable group to hamper the transmission, (v) that both $u_1$ and $u_2$ can be used to prevent the disease, and (vi) that treatment $u_2$ can be used to reduce the number of infected: Thus the time-dependent model is given by:

$$\begin{array}{ccc}\frac{dS}{dt}\hfill & =& A-\frac{{\beta }_1(1-u_3\left(t\right))SB}{K+B}-{\beta }_{2}SI(1-u_3\left(t\right))-\mu S+\alpha R+\sigma I-u_1\left(t\right)S,\hfill \\ \frac{dI}{dt}\hfill & =& \frac{{\beta }_1(1-u_3\left(t\right))SB}{K+B}+{\beta }_{2}SI(1-u_3\left(t\right))-(r+\mu +d+\sigma +u_2\left(t\right))I,\hfill \\ \frac{dR}{dt}\hfill & =& rI-\alpha R-\mu R+u_1\left(t\right)S+u_2\left(t\right)I,\hfill \\ \frac{dB}{dt}\hfill & =& \zeta I-\delta B.\hfill \end{array}$$

(1)

subject to the initial conditions:

$$S\left(0\right)=S_0,I\left(0\right)=I_0,R\left(0\right)=R_0,B\left(0\right)=B_0.$$

We described the associated model variables and parameters in

Table 1. Table showing the variables in the model.

Variables	Description
S	Susceptible Human
I	Infected Human
R	Recovered Human
B	Vibrio Cholerae

and

Table 2. Descriptions of cholera control model parameters.

Parameter	Symbol
Recruitment rate of susceptible population	A
Contribution of infected individuals to
the population of Vibrio cholerae	$\zeta$
Natural death rate of human	$\mu$
Net death rate of Vibrio cholerae	$\delta$
Acquired temporary immunity	r
Disease-induced death rate	d
Concentration of Vibrio cholerae in water	K
Human spontaneous recovery rate	$\sigma$
Use of hygiene promotion and social mobilization	$u_1$
Use of treatment by drug/oral re-hydration solution (ORS)	$u_2$
Use of safe water, hygiene, and sanitation	$u_3$
Rates of ingesting vibrios from the contaminated environment	${\beta }_1$
Rates of ingesting vibrios through human-to-human interaction	${\beta }_2$
loss of immunity	$\alpha$

.

2.2. WHO Recommendation of the Current Cholera Control Interventions Used in the Model

In line with the WHO recommendation for the control and reduction of cholera mortality and morbidity in an endemic and epidemic area, it was advised that a multifaceted approach should be put in place to mitigate cholera transmission [1]. In this work, we present three different control interventions to control and reduce cholera death and morbidity. The three different control interventions introduced here are the use of hygiene promotion and social mobilization, $u_1$, the use of treatment by drug/oral re-hydration solution (ORS), $u_2$, and the use of safe water, hygiene, and sanitation, $u_3$, which are all time-dependent control interventions. Hence, we applied the principle of optimal control theory to obtain its optimal control strategies.

The cholera control model (1) is analyzed in a biologically feasible region for both human and pathogen populations. Hence, for it to be mathematically well posed, it is important to prove that all its state variables are non-negative for all time $t>0$.

2.3. Positivity of the Solution for the Deterministic Model

Theorem 1.

Suppose that the initial data $S\left(0\right)\ge 0$, $I\left(0\right)\ge 0$, $R\left(0\right)\ge 0$, and $B\left(0\right)\ge 0$ then the solutions $\left(S\right(t),I(t),R(t),B(t\left)\right)$ of the cholera control model (1) are non-negative for all $t>0$.

For the proof (see [4]).

2.4. Boundedness of the Solution for the Deterministic Model

Theorem 2.

All solutions $\left(S\right(t),I(t),R(t),B(t\left)\right)$ of the cholera control model (1) are bounded that is ∃ a $Q>0$, and $T>0$ such that $\left(S\right(t),I(t),R(t),B(t\left)\right)\le Q$ for $t\le T$.

For the proof (see [4]).

2.5. The Invariant Region for the Deterministic Model

The model system Equation (1) will be analyzed in the feasible region and all state variables and parameters are assumed to be positive since the model has to do with population. We obtain the invariant region by Theorem 3.

Theorem 3.

The solutions of the model Equation (1) are contained in the region $\Delta \in {\mathbb{R}}_{+}^4$ and $\Delta ={\Delta }_h\cup {\Delta }_b\subset {\mathbb{R}}_{+}^3\times {\mathbb{R}}_{+}.$

For the proof (see [4]).

2.6. Existence of Disease Free Equilibrium Point (DFE)

The disease-free equilibrium point of the cholera control model (1) is obtained in our paper [4].

The disease-free equilibrium is obtained in 4 by setting the LHS (left hand side) of the model (1) to zero and then we obtained the solution of the variables in the model.

Remark 1.

Obtaining the equilibrium solution implies that the long-term behavior of the cholera control model (1) is determined by the existence of the equilibrium solution of the cholera control model (1) since $\frac{A}{\mu }>0$ and $\frac{\zeta A}{\mu \delta }>0$.

2.7. The Effective Reproduction Number

The effective reproduction number $R_p$ is obtained using the next generation matrix method (see [4]).

Theorem 4.

The effective reproduction number $R_p$ is given as $R_p=\frac{A(1-u_3)(\alpha +\mu )({\beta }_{2}K\delta +{\beta }_1\zeta )}{K\delta \mu (\mu +\alpha +u_1)(r+\mu +d+\omega +u_2)}$.

For the proof (see [4]).

2.8. Global Stability of Disease Free Equilibrium Point (DFE)

Theorem 5.

If $R_p\le 1$ then the disease free equilibrium (DFE) ${\mathsf{\Omega }}_0$ is globally asymptotically stable on $\Delta .$

Proof. 

Suppose $R_p\le 1$, then ∃ only DFE:

${\mathsf{\Omega }}_0=(S^{\ast },I^{\ast },B^{\ast },R^{\ast })=(\frac{A(\alpha +\mu )}{\mu (\mu +\alpha +u_1)},0,0,\frac{u_{1}A}{\mu (\mu +\alpha +u_1)}).$

We consider the Lyapunov function $V(S,I,B,R):{\Re }^{+}⟶{\Re }^{+}$, $V(S,I,B,R)=PI$ where $P>0$ and: $\frac{dV}{dt}=P\frac{dI}{dt}.$ Let $P=\frac{{\beta }_{2}K\delta +{\beta }_1\zeta }{K\delta \mu }$ and

$$\frac{dV}{dt}=\frac{{\beta }_1(1-u_3)S^{\ast }{B}^{\ast }}{K+B}+{\beta }_2{S}^{\ast }{I}^{\ast }(1-u_3)-(r+\mu +d+\omega +u_2)I^{\ast }$$

$$\Rightarrow \frac{dV}{dt}=P\left[\frac{\left({\beta }_1(1-u_3)A(\alpha +\mu )\right)}{\mu (\mu +\alpha +u_1)}+\left\{\frac{{\beta }_{2}A(\alpha +\mu )(1-u_3)}{\mu (\mu +\alpha +u_1)}-(r+\mu +d+\omega +u_2)\right\}{I}^{\ast }\right]$$

$$\Rightarrow \frac{dV}{dt}=(r+\mu +d+\omega +u_2)\left[P\frac{{\beta }_1(1-u_3)A(\alpha +\mu )B^{\ast }}{(K+B)\mu (\mu +\alpha +u_1)}+(R_p-1)I^{\ast }\right],$$

where $\frac{dV}{dt}\le 0$. Hence $R_0\le 1$ which implies that $\frac{dV}{dt}=0$, whenever $I^{\ast }=0$.

Therefore the maximum invariant set in $\left\{(S,I,B,R)\in \Delta |\frac{dV}{dt}\le 0\right\}$ is the singleton set $\left\{{\mathsf{\Omega }}_0\right\}$. Therefore, the global stability of ${\mathsf{\Omega }}_0$ follows from the LaSalle’s invariant principle when $R_p\le 1$ [31]. □

Remark 2.

This epidemiological implication of this is that cholera can be eliminated from the community if $R_p<1$. If $R_p<1$ then the average of an infected individual produces less than one new infected individual/mosquito over the cause of its infectious period and the infection dies out. This also implies that the stability of the equilibrium of the cholera control model (1) is not just for those initial conditions that are close to it but for almost all initial conditions.

2.9. Global Stability of the Endemic Equilibrium Point

Theorem 6.

Suppose $R_p>1$, then ∃ an epidemic equilibrium ${\mathsf{\Omega }}_1$ that is globally stable.

Proof. 

Suppose $R_p>1$, then ∃ the endemic equilibrium point. We consider the Lyapunov function:

$$V(S,I,B,R)=\frac{1}{2}{(S-S^{\ast })}^2+\frac{1}{2}{(I-I^{\ast })}^2+\frac{1}{2}{(B-B^{\ast })}^2\frac{1}{2}{(R-R^{\ast })}^2$$

(2)

$$\frac{dV}{dt}=Z-Y,$$

where:

$$\begin{array}{c}Z=A(S-S^{\ast })+\frac{{\beta }_1(S-S^{\ast })(1-u_3)S{B}^{\ast }}{K+B-B^{\ast }}+\frac{{\beta }_1(S-S^{\ast })(1-u_3)S^{\ast }B}{K+B-B^{\ast }}+\hfill \\ {\beta }_2(1-u_3)(S-S^{\ast })SI+{\beta }_2(1-u_3)(S-S^{\ast })S^{\ast }I+\hfill \\ \mu (S-S^{\ast })S^{\ast }+\alpha (S-S^{\ast })R+\delta (S-S^{\ast })I+u_1(S-S^{\ast })S+\frac{{\beta }_1(I-I^{\ast })(1-u_3)SB}{K+B-B^{\ast }}+\hfill \\ \frac{{\beta }_1(I-I^{\ast })(1-u_3)S^{\ast }{B}^{\ast }}{K+B-B^{\ast }}+{\beta }_2(I-I^{\ast })(1-u_3)SI+\hfill \\ {\beta }_2(I-I^{\ast })(1-u_3)S^{\ast }{I}^{\ast }+\mu (I-I^{\ast })+d(I-I^{\ast })I+\omega (I-I^{\ast })+u_2(I-I^{\ast })I+\hfill \\ u_2(I-I^{\ast })I^{\ast }+I(B-B^{\ast })\zeta +\delta (B-B^{\ast })B^{\ast }+r(R-R^{\ast })I+\hfill \\ \alpha (R-R^{\ast })R^{\ast }+\mu (R-R^{\ast })R^{\ast }+u_1(R-R^{\ast })S+u_2(R-R^{\ast })I\hfill \end{array}$$

(3)

and

$$\begin{array}{c}Y=\frac{{\beta }_1(S-S^{\ast })(1-u_3)S{B}^{\ast }}{K+B-B^{\ast }}+\frac{{\beta }_1(S-S^{\ast })(1-u_3)S^{\ast }{B}^{\ast }}{K+B-B^{\ast }}+\hfill \\ {\beta }_2(1-u_3)(S-S^{\ast })SI+{\beta }_2(1-u_3)(S-S^{\ast })S^{\ast }{I}^{\ast }+\mu (S-S^{\ast })S+\hfill \\ \alpha (S-S^{\ast })R^{\ast }+\delta (S-S^{\ast })I^{\ast }+u_1(S-S^{\ast })S+\hfill \\ \frac{{\beta }_1(I-I^{\ast })(1-u_3)S{B}^{\ast }}{K+B-B^{\ast }}+\frac{{\beta }_1(I-I^{\ast })(1-u_3)S^{\ast }B}{K+B-B^{\ast }}+\hfill \\ {\beta }_2(I-I^{\ast })(1-u_3)S{I}^{\ast }+{\beta }_2(I-I^{\ast })(1-u_3)S^{\ast }I+I(I-I^{\ast })r+\hfill \\ r(I-I^{\ast })I^{\ast }+\mu (I-I^{\ast })I^{\ast }+d(I-I^{\ast })I^{\ast }+\omega (I-I^{\ast })I^{\ast }+\hfill \\ I^{\ast }(B-B^{\ast })\zeta +\delta (B-B^{\ast })B+r(R-R^{\ast })I^{\ast }+\alpha (R-R^{\ast })R+\hfill \\ \mu (R-R^{\ast })R+u_1(R-R^{\ast })S^{\ast }+u_2(R-R^{\ast })I^{\ast }.\hfill \end{array}$$

(4)

Therefore $\frac{dV}{dt}=Z-Y$. It follows from the above that $\frac{dV}{dt}\le 0$ if and only if $Z\le Y$. Thus, the largest compact invariant set in $\left\{(S,I,B,R)\in \Delta |\frac{dV}{dt}=0\right\}$ is a singleton $\left\{{\mathsf{\Omega }}_0\right\}$ when P is the endemic equilibrium system and since $\frac{dV}{dt}$ is a non-positive definite, the Lyapunov theorem implies that the endemic equilibrium is globally asymptotically stable in the region. □

3. Sensitivity and Elasticity Analysis

The sensitivity and elasticity analysis of model (1) was performed in order to establish the relation of the model to the model parameters [32]. To understand how to optimally reduce the transmission of cholera, human morbidity, and mortality due to cholera disease, it is essential to know the relative importance of the different factors responsible for its transmission. The numerical sensitivity analysis is carried out to detect parameters with high impact on $R_p$ and which should be targeted by the control strategies. In order to qualitatively evaluate the impact of these intervention strategies on the spread of cholera disease, an analytical sensitivity and elasticity analysis on all parameters which drives the disease dynamics was carried out with respect to reproduction number $R_p$. The sensitivity and elasticity indexes of the reproduction number that measures initial disease transmission using the approach by [32,33]. Sensitivity and elasticity indices provide the opportunity to measure the relative change in a state variable when a parameter changes. The normalized forward sensitivity index of a variable to a parameter is a ratio of the relative change in the parameter, when a variable is a differentiable function of the parameter the sensitivity index may be alternatively defined using partial derivatives.

The normalized forward sensitivity index of a variable, u, with respect to differentiating an index on a parameter, y is defined as: $V_y^u=\frac{\partial u}{\partial y}\times \frac{y}{u}$. The implicit expression for $R_p$ is been derived analytical for the sensitivity of $R_p$ as $V_y^{R_p}=\frac{\partial R_p}{\partial y}\times \frac{y}{R_p}$ to each parameter involved in $R_p$. The following result establishes the sensitivity and elasticity of the reproduction number.

Theorem 7.

The sensitivity and elasticity of the quantity $R_p$ with respect to the parameter ${\beta }_2$, is such that $V_{{\beta }_2}^{R_p}=\frac{\partial R_p}{\partial {\beta }_2}\times \frac{{\beta }_2}{R_p}$ to each parameter involved in $R_p$.

Proof. 

The normalized forward sensitivity index of a variable, b, that depends differentiably on index on a parameter, z is defined as: $V_z^b=\frac{\partial b}{\partial z}\times \frac{z}{b}$. Given the implicit expression for $R_p$ then the analysis for the sensitivity of $R_p$ is derived as $V_z^{R_p}=\frac{\partial R_p}{\partial z}\times \frac{z}{R_p}$ to each parameter involved in $R_p$.

Important indices are given thus: $V_{{\beta }_1}^{R_p}$, $V_{{\beta }_2}^{R_p}$, $V_A^{R_p}$,$V_r^{R_p}$, $V_d^{R_p}$, $V_{\delta }^{R_p}$,$V_{\sigma }^{R_p}$, $V_{\mu }^{R_p}$,$V_K^{R_p}$,$V_{\zeta }^{R_p}$.

They are given in

Table 3. Table showing numerical values of sensitivity indices for the cholera model.

Parameter	Parameter Values	Sensitivity to ${\mathit{R}}_{\mathit{p}}$
$\zeta$	100	+0.24511
${\beta }_1$	0.2143	+0.24511
${\beta }_2$	0.000002	+0.755
A	15	+0.9999
d	0.015	−0.2174
$\delta$	0.033	−0.24511
$\mu$	0.0000548	−1.0000
$\sigma$	0.05	−0.7245
r	0.004	−0.0579
K	${10}^9$	−0.24511
$\alpha$	0.025	−0.24511

:

This completes the proof. □

Interpretation of Sensitivity Indices

The most sensitive parameter is the rate of recruitment of a susceptible population, A. This is followed by the rate of ingesting vibrios through human-to-human interactions, ${\beta }_2$ and then followed by ${\beta }_1$ and $\zeta$ which are the rates of ingesting vibrios from the contaminated environment and the contribution of infected individuals to the population of vibrios cholerae.

4. Stochastic Model of the Transmission Dynamics of the Cholera Model and Its Biological Description

Here, we considered that the white noise depends on the size of the corresponding population. In this case, we introduce randomness into model (1) by considering that the white noise depends on the corresponding population size. More precisely, this amounts to consider a corresponding population of the form $e_{1}S\left(t\right)dW\left(t\right)$, where $e_1$ denotes the intensity of the random perturbation and ${\left(W\left(t\right)\right)}_{t\ge 0}$ is a one-dimensional Brownian motion defined on a complete probability space $(\mathsf{\Omega },\mathfrak{F},{\left({\mathfrak{F}}_t\right)}_{t\ge 0}\mathbb{P})$. Therefore, we obtain a stochastic model of cholera described by stochastic differential equations (SDE):

$$\begin{array}{ccc}dS\hfill & =& \left(A-\frac{{\beta }_1(1-u_3\left(t\right))SB}{K+B}-{\beta }_{2}SI(1-u_3\left(t\right))-\mu S+\alpha R+\sigma I-u_1\left(t\right)S\right)dt+e_{1}S\left(t\right)dW\left(t\right),\hfill \\ dI\hfill & =& \left(\frac{{\beta }_1(1-u_3\left(t\right))SB}{K+B}+{\beta }_{2}SI(1-u_3\left(t\right))-(r+\mu +d+\sigma +u_2\left(t\right))I\right)dt+e_{2}I\left(t\right)dW\left(t\right),\hfill \\ dR\hfill & =& \left(rI-\alpha R-\mu R+u_1\left(t\right)S+u_2\left(t\right)I\right)dt+e_{3}R\left(t\right)dW\left(t\right),\hfill \\ dB\hfill & =& \left(\zeta I-\delta B\right)dt+e_{4}B\left(t\right)dW\left(t\right).\hfill \end{array}$$

(5)

Let $X\left(t\right)=\left(\begin{array}{c}S\left(t\right)\\ I\left(t\right)\\ R\left(t\right)\\ B\left(t\right)\end{array}\right)$, system (5) can be re-written in the form of a single stochastic differential equations of the form $dX\left(t\right)=F\left(X\right(t),t)dt+G\left(X\right(t),t)dW\left(t\right)$, where the function $F:{\mathbb{R}}_{+}^2\times {\mathbb{R}}_{+}^2\to {\mathbb{R}}^2$ is defined by:

$$F=\left(\begin{array}{c}A-\frac{{\beta }_1(1-u_3\left(t\right))SB}{K+B}-{\beta }_{2}SI(1-u_3\left(t\right))-\mu S+\alpha R+\sigma I-u_1\left(t\right)S\\ \frac{{\beta }_1(1-u_3\left(t\right))SB}{K+B}+{\beta }_{2}SI(1-u_3\left(t\right))-(r+\mu +d+\sigma +u_2\left(t\right))I\\ rI-\alpha R-\mu R+u_1\left(t\right)S+u_2\left(t\right)I\\ \zeta I-\delta B\end{array}\right),$$

and the function $G:{\mathbb{R}}_{+}^2\times {\mathbb{R}}_{+}^2\to {\mathbb{R}}^2$ is given by:

$$G=\left(\begin{array}{c}{e}_{1}S\left(t\right)\\ e_{2}I\left(t\right)\\ e_{3}R\left(t\right)\\ e_{4}B\left(t\right)\end{array}\right).$$

Positivity of the Solution for the SDE Model

Lemma 1.

Suppose $(S_0,I_0,R_0,B_0)\in \mathsf{\Omega }$, then model (5) permits a unique solution $(S^{\ast },I^{\ast },R^{\ast },B^{\ast })$ on $t\ge 0$ and this solution remains in Ω with probability 1.

Proof. 

If $(S_0,I_0,R_0,B_0)\in \mathsf{\Omega }$ then the total population shows the equation:

$$\frac{d{N}_h^{\ast }}{dt}=A-\mu N_h^{\ast }-d{I}^{\ast },$$

and

$$\frac{d{B}^{\ast }}{dt}=\frac{\zeta A}{\mu }-\delta B^{\ast }.$$

Suppose $(S^{\ast },I^{\ast },R^{\ast },B^{\ast })\in {\mathbb{R}}_{+}^4$ for all $0\le t\le T$, then we almost surely obtain:

$$d{N}_h^{\ast }\le (A-\mu N_h^{\ast })dq,$$

$$d{B}^{\ast }\le (\frac{\zeta A}{\mu }-\delta B^{\ast })dq.$$

Therefore, by the comparison theorem, we obtain:

$$N_h^{\ast }\le \frac{A}{\mu }+(N_h^{\ast }\left(0\right)-\frac{A}{\mu })e^{-\mu t}<\frac{A}{\mu }$$

almost surely for all $t\in [0,T].$

$$B^{\ast }\le \frac{\zeta A}{\mu \delta }+(B^{\ast }\left(0\right)-\frac{\zeta A}{\mu \delta })e^{-\mu t}<\frac{\zeta A}{\mu \delta }$$

$a.s.$∀$t\in [0,T].$

Then, $N_h^{\ast }\le \frac{A}{\mu }$&$B^{\ast }\le \frac{\zeta A}{\mu \delta }$, so we have $(S^{\ast },I^{\ast },R^{\ast },B^{\ast })\in \mathsf{\Omega }$ a.s. ∀$t\in [0,T].$ □

5. Optimal Control Analysis

In this section we examine the optimal level of effort that will be required to control the cholera disease. We define an objective functional J, which is meant to minimize the number of infective humans I(t), the number of pathogens population at time t, B(t), and the cost of applying the controls $u_1,u_2,u_3$:

$$J(u_1,u_2,u_3)={\int }_0^{t_f}{v}_{1}I\left(t\right)+v_{2}B\left(t\right)+\frac{v_3}{2}{u}_1^2+\frac{v_4}{2}{u}_2^2+\frac{v_5}{2}{u}_3^{2}dt,$$

(6)

where $v_1,v_2,v_3,v_4$, and $v_5$ are positive weights.

Subjected to the model Equation (9):

$$\begin{array}{ccc}\frac{dS}{dt}\hfill & =& A-\frac{{\beta }_1(1-u_3)SB}{K+B}-{\beta }_{2}SI(1-u_3)-\mu S+\alpha R+\sigma I-u_{1}S\hfill \\ \frac{dI}{dt}\hfill & =& \frac{{\beta }_1(1-u_3)SB}{K+B}+{\beta }_{2}SI(1-u_3)-(r+\mu +d+\omega +u_2)I\hfill \\ \frac{dR}{dt}\hfill & =& rI-\alpha R-\mu R+u_{1}S+u_{2}I\hfill \\ \frac{dB}{dt}\hfill & =& \zeta I-\delta B\hfill \end{array}$$

(7)

with initial conditions $S\left(0\right)>0,I\left(0\right)>0,R\left(0\right)>0,B\left(0\right)>0.$

We choose a quadratic cost on the controls because this is similar to what is in the literature on epidemic controls [34,35,36,37,38]. With the given objective functional J($u_1,u_2,u_3$), our goal is to minimize the total number of infected humans, the total number of pathogen population at time (t), while minimizing the cost of control $u_1$(t), $u_2$(t), and $u_3$(t). Hence, we seek an optimal control $u_1^{\ast },u_2^{\ast }$, and $u_3^{\ast }$ such that:

$$J(u_1^{\ast },u_2^{\ast },u_3^{\ast })=min\left\{J(u_1,u_2,u_3)|u_1,u_2,u_3\in U,\right\}$$

(8)

where $U=(u_1,u_2,u_3)$ such that $u_1,u_2,u_3$ are measurable with $0\le u_1\le 1,0\le u_2\le 1$ and $0\le u_3\le 1$ for $t\in [0,t_f]$ in the control set. In the next section, we discuss the existence of optimal control in the control system.

Existence of Optimal Control

Existence of the control system (7) such that $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$, and $B\left(t\right)$ denote the state variables with the control variables with $u_1$, $u_2$, and $u_3$ as control variables. We consider a control system (7) with initial conditions for existence. Then we rewrite (7) in the following form in line with [39]:

$$F(x\left(t\right),u\left(t\right))=Cx\left(t\right)+u_{1,2}\left(t\right)x\left(t\right)+u_3\left(t\right)Z\left(x\left(t\right)\right)+H\left(x\left(t\right)\right)$$

where,

$x\left(t\right)=\left[\begin{array}{c}S\left(t\right)\\ I\left(t\right)\\ R\left(t\right)\\ B\left(t\right)\end{array}\right]$, $C=\left[\begin{array}{cccc}-\mu & 0& u_1& 0\\ 0& -(r+\mu +d+\omega )& 0& 0\\ 0& r& -(\alpha +\mu )& 0\\ 0& \zeta & 0& -\delta \end{array}\right]$,

$u_{1,2}\left(t\right)=\left[\begin{array}{cccc}-u_1\left(t\right)& 0& 0& 0\\ 0& -u_2\left(t\right)& 0& 0\\ u_1\left(t\right)& u_2\left(t\right)& 0& 0\\ 0& 0& 0& 0\end{array}\right]$, $u_3\left(t\right)=\left[\begin{array}{cccc}{u}_3\left(t\right)& u_1\left(t\right)& 0& 0\\ -u_3\left(t\right)& -u_3\left(t\right)& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$,

$H\left(x\left(t\right)\right)=\left[\begin{array}{c}A-\frac{{\beta }_{1}SB}{K+B}-{\beta }_{2}SI\\ \frac{{\beta }_{1}SB}{K+B}+{\beta }_{2}SI\\ 0\\ 0\end{array}\right]$, $Z\left(x\left(t\right)\right)=\left[\begin{array}{c}\frac{{\beta }_{1}SB}{K+B}+{\beta }_{2}SI\\ -\frac{{\beta }_{1}SB}{K+B}-{\beta }_{2}SI\\ 0\\ 0\end{array}\right]$.

Theorem 8.

Suppose the control system is uniform Lipschitz continuous, then $|F(x,u)-F(x^{\ast },u)|$$\le K|x-x^{\ast }|$ such that $K=max\{Q,\parallel C\parallel \}<\infty$.

Proof. 

Let:

$$|F(x,u\left(t\right))-F(x^{\ast },u\left(t\right))|$$

$$=|C+u_{1,2}\left(t\right)(x\left(t\right)-x^{\ast }\left(t\right))+u_3\left(t\right)(Z\left(x\left(t\right)\right)-Z\left(x^{\ast }\left(t\right)\right))+H\left(x\left(t\right)\right)-H\left(x^{\ast }\right)|$$

$$\le |C+u_{1,2}\left(t\right)\left|\right|(x\left(t\right)-x^{\ast }\left(t\right))|+u_3\left(t\right)|(Z\left(x\left(t\right)\right)-Z\left(x^{\ast }\left(t\right)\right))|+|H\left(x\left(t\right)\right)-H\left(x^{\ast }\right)|$$

$$\le (\parallel C\parallel +\underset{0,t_f}{max}|u_{1,2}\left(t\right)\left|\right)|(x\left(t\right)-x^{\ast }\left(t\right))|+\underset{0,t_f}{max}|u_3\left(t\right)\left|\right|x\left(t\right)-x^{\ast }\left(t\right)|+|x\left(t\right)-x^{\ast }|$$

$$\le (\parallel C\parallel +\underset{0,t_f}{max}|u_{1,2}\left(t\right)|+\underset{0,t_f}{max}|u_3\left(t\right)|+1)\parallel x\left(t\right)-x^{\ast }\left(t\right)\parallel \parallel F(x,u\left(t\right))-F(x^{\ast },u\left(t\right))\parallel$$

$$\le Q\parallel x\left(t\right)-x^{\ast }\left(t\right)\parallel ,$$

where,

$$Q=\parallel C\parallel +\underset{0,t_f}{max}|u_{1,2}\left(t\right)|+\underset{0,t_f}{max}\left|u_3\left(t\right)\right|+1<\infty .$$

Therefore, the function F is uniformly Lipschitz continuous. □

The necessary condition that an optimal control must satisfy comes from Pontryagin’s Maximum Principle [40]. The principle converts Equations (7) and (8) into a problem of minimizing pointwise a Hamiltonian H, with respect to $u_1$, $u_2$, and $u_3$:

$$\begin{array}{c}H=v_{1}I\left(t\right)+v_{2}B\left(t\right)+\frac{v_3}{2}{u}_1^2+\frac{v_4}{2}{u}_2^2+\frac{v_5}{2}{u}_3^2+\hfill \\ {\lambda }_S\left\{A-\frac{{\beta }_1(1-u_3)SB}{K+B}-{\beta }_{2}SI(1-u_3)-\mu S+\alpha R+\sigma I-u_{1}S\right\}+\hfill \\ {\lambda }_I\left\{\frac{{\beta }_1(1-u_3)SB}{K+B}+{\beta }_{2}SI(1-u_3)-(r+\mu +d+\omega +u_2)I\right\}+\hfill \\ {\lambda }_R\left\{rI-\alpha R-\mu R+u_{1}S+u_{2}I\right\}+\hfill \\ {\lambda }_B\left\{\zeta I-\delta B\right\},\hfill \end{array}$$

(9)

where ${\lambda }_S$, ${\lambda }_I$, ${\lambda }_R$, and ${\lambda }_B$ are the adjoint variables. When we apply Pontryagin’s Maximum Principle [40] and the existence result for the optimal control from [11].

Theorem 9.

Suppose $S^{\ast }\left(t\right)$, $I^{\ast }\left(t\right)$, $R^{\ast }\left(t\right)$, and $B^{\ast }\left(t\right)$ are optimal state solutions with associated optimal control variables $u_1^{\ast }\left(t\right)$, $u_2^{\ast }\left(t\right)$, and $u_3^{\ast }\left(t\right)$ for the optimal control problems (7) and (8). Then there exists adjoint variables ${\lambda }_S\left(t\right)$, ${\lambda }_I\left(t\right)$, ${\lambda }_R\left(t\right)$, and ${\lambda }_B\left(t\right)$ which satisfy:

$$\begin{array}{c}-\acute{{\lambda }_S}=\left(\frac{{\beta }_1(1-u_3)B}{K+B}+{\beta }_{2}I(1-u_3)\right)({\lambda }_S-{\lambda }_R)+(\mu +u_1){\lambda }_S-u_1{\lambda }_R\hfill \\ -\acute{{\lambda }_I}=-v_1-\sigma {\lambda }_S+{\beta }_2(1-u_3)S({\lambda }_S-{\lambda }_I)+(r+\mu +d+\sigma +u_2){\lambda }_I-(r+u_2){\lambda }_R\hfill \\ -\acute{{\lambda }_R}=(\alpha +\mu ){\lambda }_R-\alpha {\lambda }_S\hfill \\ -\acute{{\lambda }_B}=-v_2+\frac{{K\beta }_1(1-u_3)S}{{(K+B)}^2}({\lambda }_S-{\lambda }_I)+\delta {\lambda }_B\hfill \end{array}$$

(10)

and with transversality conditions:

$${\lambda }_S\left(t_f\right)={\lambda }_I\left(t_f\right)={\lambda }_R\left(t_f\right)={\lambda }_B\left(t_f\right)=0.$$

(11)

Furthermore, we obtain the optimal controls $u_1^{\ast },u_2^{\ast }$, and $u_3^{\ast }$:

$$\begin{array}{c}{u}_1^{\ast }=max\left\{0,min\left(1,\left(\frac{S({\lambda }_S-{\lambda }_R)}{v_3}\right)\right)\right\}\hfill \\ u_2^{\ast }=max\left\{0,min\left(1,\left(\frac{I({\lambda }_I+{\lambda }_R)}{v_4}\right)\right)\right\}\hfill \\ u_3^{\ast }=max\left\{0,min\left(1,\left(\frac{{\beta }_{1}SB({\lambda }_S-{\lambda }_I)+(K+B){\beta }_{2}SI({\lambda }_I-{\lambda }_S)}{v_5(K+B)}\right)\right)\right\}.\hfill \end{array}$$

(12)

Proof. 

We applied the Hamiltonian Equation (9) to obtain the adjoint equations and the transversality conditions. We set $S\left(t\right)=S^{\ast }\left(t\right)$, $I\left(t\right)=I^{\ast }\left(t\right)$, $R\left(t\right)=R^{\ast }\left(t\right)$, and $B\left(t\right)=B^{\ast }\left(t\right)$, differentiating the Hamiltonian Equation (9) with respect to $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$, and $B\left(t\right)$, to yield:

$$\begin{array}{c}-\acute{{\lambda }_S}=-\frac{\partial H}{\partial S}=\left(\frac{{\beta }_1(1-u_3)B}{K+B}+{\beta }_{2}I(1-u_3)\right)({\lambda }_S-{\lambda }_R)+(\mu +u_1){\lambda }_S-u_1{\lambda }_R\hfill \\ -\acute{{\lambda }_I}=-\frac{\partial H}{\partial I}=-v_1-\sigma {\lambda }_S+{\beta }_2(1-u_3)S({\lambda }_S-{\lambda }_I)+(r+\mu +d+\sigma +u_2){\lambda }_I-(r+u_2){\lambda }_R\hfill \\ -\acute{{\lambda }_R}=-\frac{\partial H}{\partial R}=(\alpha +\mu ){\lambda }_R-\alpha {\lambda }_S\hfill \\ -\acute{{\lambda }_B}=-\frac{\partial H}{\partial B}=-v_2+\frac{{K\beta }_1(1-u_3)S}{{(K+B)}^2}({\lambda }_S-{\lambda }_I)+\delta {\lambda }_B.\hfill \end{array}$$

(13)

By the optimality conditions, we obtain

$\frac{\partial H}{\partial u_1}=v_3{u}_1-{\lambda }_{S}S+{\lambda }_{R}S=0$ at $u_1^{\ast }$, ⇒ $u_1^{\ast }=\frac{S({\lambda }_S-{\lambda }_R)}{v_3}$

$\frac{\partial H}{\partial u_2}=v_4{u}_2-{\lambda }_I-{\lambda }_R=0$ at $u_2^{\ast }$, ⇒ $u_2^{\ast }=\frac{({\lambda }_I-{\lambda }_R)}{v_4}$

$\frac{\partial H}{\partial u_3}=v_5{u}_3+{\lambda }_S\frac{{\beta }_{1}SB}{K+B}+{\lambda }_S{\beta }_{2}SI-{\lambda }_I\frac{{\beta }_{1}SB}{K+B}-{\lambda }_I{\beta }_{2}SI=0$ at $u_3^{\ast }$,

⇒$u_3^{\ast }=\frac{{\beta }_{1}SB({\lambda }_S-{\lambda }_I)+(K+B){\beta }_{2}SI({\lambda }_I-{\lambda }_S)}{v_5(K+B)}.$

$$u_1^{\ast }=\left\{\begin{array}{c}0if{Z}_1^{\ast }\le 0,\hfill \\ \\ Z_1^{\ast }if0<Z_1^{\ast }<1,\hfill \\ \\ 1if{Z}_1^{\ast }\ge 1\hfill \\ \end{array}\right.$$

(14)

where

$$Z_1^{\ast }=\frac{S({\lambda }_S-{\lambda }_R)}{v_3}$$

$$u_2^{\ast }=\left\{\begin{array}{c}0if{Z}_2^{\ast }\le 0,\hfill \\ \\ Z_2^{\ast }if0<Z_2^{\ast }<1,\hfill \\ \\ 1if{Z}_2^{\ast }\ge 1\hfill \\ \end{array}\right.$$

(15)

where

$$Z_2^{\ast }=\frac{({\lambda }_I+{\lambda }_R)}{v_4}$$

and

$$u_3^{\ast }=\left\{\begin{array}{c}0if{Z}_3^{\ast }\le 0,\hfill \\ \\ Z_3^{\ast }if0<Z_3^{\ast }<1,\hfill \\ \\ 1if{Z}_3^{\ast }\ge 1\hfill \\ \end{array}\right.$$

(16)

where

$$Z_3^{\ast }=\frac{{\beta }_{1}SB({\lambda }_S-{\lambda }_I)+(K+B){\beta }_{2}SI({\lambda }_I-{\lambda }_S)}{v_5(K+B)}.$$

These can be rewritten in compact form:

$$u_1^{\ast }=max\{min\{\frac{S({\lambda }_S-{\lambda }_R)}{v_3},1\},0\},$$

(17)

$$u_2^{\ast }=max\{min\{\frac{({\lambda }_I+{\lambda }_R)}{v_4},1\},0\},$$

(18)

$$u_3^{\ast }=max\{min\{\frac{{\beta }_{1}SB({\lambda }_S-{\lambda }_I)+(K+B){\beta }_{2}SI({\lambda }_I-{\lambda }_S)}{v_5(K+B)},1\},0\}.$$

(19)

□

Therefore, according to Equations (17)–(19) for $u_1^{\ast }$, $u_2^{\ast }$, $u_3^{\ast }$ the characterization of the optimal controls. We obtained the optimal control and states by solving the optimality system which includes the state system (7) with boundary conditions, the adjoint system (10), and the characterization of the optimal control (17)–(19). We used the initial conditions, transversality conditions alongside the characterization of the optimal controls $u_1^{\ast }$, $u_2^{\ast }$, $u_3^{\ast }$ provided in (17)–(19) to solve the optimality system. Furthermore, the second derivative of the Hamiltonian with respect to $u_1$, $u_2$, and $u_3$ are positive meaning a minimum at $u_1^{\ast }\left(t\right)$, $u_2^{\ast }\left(t\right)$, and $u_3^{\ast }\left(t\right)$. By substituting the values of $u_1^{\ast }\left(t\right)$, $u_2^{\ast }\left(t\right)$, and $u_3^{\ast }\left(t\right)$ in the control model Equation (8) and we arrived at the following new system:

$$\begin{array}{ccc}\frac{dS}{dt}\hfill & =& A-\frac{{\beta }_1(1-u_3^{\ast }\left(t\right))S^{\ast }\left(t\right)B^{\ast }\left(t\right)}{K+B^{\ast }\left(t\right)}-{\beta }_2{S}^{\ast }\left(t\right)I^{\ast }\left(t\right)(1-u_3^{\ast }\left(t\right))-\mu S^{\ast }\left(t\right)+\alpha R^{\ast }\left(t\right)+\sigma I^{\ast }\left(t\right)-u_1^{\ast }\left(t\right)S^{\ast }\left(t\right)\hfill \\ \frac{dI}{dt}\hfill & =& \frac{{\beta }_1(1-u_3^{\ast }\left(t\right))S^{\ast }\left(t\right)B^{\ast }\left(t\right)}{K+B^{\ast }\left(t\right)}+{\beta }_2{S}^{\ast }\left(t\right)I^{\ast }\left(t\right)(1-u_3^{\ast }\left(t\right))-(r+\mu +d+\omega +u_2^{\ast }\left(t\right))I^{\ast }\left(t\right)\hfill \\ \frac{dR}{dt}\hfill & =& r{I}^{\ast }\left(t\right)-\alpha R^{\ast }\left(t\right)-\mu R^{\ast }\left(t\right)+u_1^{\ast }\left(t\right)S^{\ast }\left(t\right)+u_2^{\ast }\left(t\right)I^{\ast }\left(t\right)\hfill \\ \frac{dB}{dt}\hfill & =& \zeta I^{\ast }\left(t\right)-\delta B^{\ast }\left(t\right)\hfill \end{array}$$

(20)

and then together with the Hamiltonian $H^{\ast }$ at $(t,S^{\ast },I^{\ast },R^{\ast },B^{\ast },u_1^{\ast },u_2^{\ast },u_3^{\ast },{\lambda }_S,{\lambda }_I,{\lambda }_R,{\lambda }_B)$:

$$\begin{array}{c}{H}^{\ast }=v_1{I}^{\ast }\left(t\right)+v_2{B}^{\ast }\left(t\right)+\frac{v_3}{2}{u}_1^{\ast ,2}+\frac{v_4}{2}{u}_2^{\ast ,2}+\frac{v_5}{2}{u}_3^{\ast ,2}+\hfill \\ {\lambda }_S\left\{A-\frac{{\beta }_1(1-u_3^{\ast })S^{\ast }{B}^{\ast }}{K+B^{\ast }}-{\beta }_2{S}^{\ast }{I}^{\ast }(1-u_3^{\ast })-\mu S^{\ast }+\alpha R^{\ast }+\sigma I^{\ast }-u_1^{\ast }{S}^{\ast }\right\}+\hfill \\ {\lambda }_I\left\{\frac{{\beta }_1(1-u_3^{\ast })S^{\ast }{B}^{\ast }}{K+B^{\ast }}+{\beta }_2{S}^{\ast }{I}^{\ast }(1-u_3^{\ast })-(r+\mu +d+\omega +u_2^{\ast })I^{\ast }\right\}+\hfill \\ {\lambda }_R\left\{r{I}^{\ast }-(\alpha +\mu )R^{\ast }+u_1^{\ast }{S}^{\ast }+u_2^{\ast }{I}^{\ast }\right\}+\hfill \\ {\lambda }_B\left\{\zeta I^{\ast }-\delta B^{\ast }\right\}.\hfill \end{array}$$

(21)

Therefore, we obtain the optimal control and state system by numerically solving the system (20) and (21). We applied the parameter values represented in

Table 4. Table showing numerical values of parameters used in the simulations.

Parameter	Symbol	Value	Source
Recruitment rate of susceptible population	A	15/day	Assumed
Contribution of infected individuals to
the population of Vibrio cholerae	$\zeta$	100 cells/L-per day	[5]
Natural death rate of human	$\mu$	0.0000548/day	[18]
Net death rate of Vibrio cholerae	$\delta$	0.033/day	[5]
Acquired temporary immunity	r	0.004/day	[8,26]
Disease-induced death rate	d	0.015/day	[26]
Concentration of Vibrio cholerae in water	K	${10}^9$ cells/L	[26]
Human spontaneous recovery rate	$\sigma$	0.05	[18]
Use of oral cholera vaccine	$u_1$	0.9	Assumed
Use of treatment by ORS	$u_2$	0.91	Assumed
Personal hygiene/sanitation	$u_3$	0.94	Assumed
Rates of ingesting vibrios from the CO *	${\beta }_1$	0.2143	Assumed
Eates of ingesting vibrios through H-H int. *	${\beta }_2$	0.000002	Assumed
Loss of immunity	$\alpha$	0.025	Assumed

* contaminated environment = CO; human-to-human interaction = H-H int.

to get the numerical solutions for the optimality system by using numerics.

6. Numerical Simulation and Graphical Illustration of the Model

Here, we run some numerical simulations to illustrate the behavior of reproduction number and the behavior of the model. A numerical simulation of the cholera model (1) using the original system variables before normalization was conducted using MATLAB’s ode45. We observe that when an endemic equilibrium is reached the system converges to a steady state that is asymptotically stable. In order to illustrate some of the numerical results of the study, numerical simulations of the model (1) are carried out using a set of reasonable parameter values given in Table 4. Parameters were obtained from different literature like [31]. Interesting results can be found also in [41,42,43,44,45,46]. We simulate the model system by using the ODE solver coded in Matlab programming language. However, it should be known that the parameters used are theoretical so they may or may not be biologically realistic. Numerical simulation and graphical illustrations are carried out in order to verify some of the analytical results. Different initial starts have been used to perform the computer simulations for different cases and displayed graphically in the figures below:

• Strategy A: Employing hygiene promotion and social mobilisation $\left(u_1\right)$, only.
• Strategy B: Treatment of the symptomatic individuals with drug/oral re-hydration solution (ORS) $\left(u_2\right)$, only.
• Strategy C: Employing sanitation, hygiene, and safe water $\left(u_3\right)$, only.
• Strategy D: Employing the control interventions $(u_1,u_2)$.
• Strategy E: Employing the control interventions $(u_1,u_3)$.
• Strategy F: Employing the control interventions $(u_2,u_3)$.
• Strategy G: Employing all three control interventions $(u_1,u_2,u_3)$.

6.1. Strategy A: Employing Hygiene Promotion and Social Mobilization $\left(u_1\right)$, Only

Here, we use only the control measure $\left(u_1\right)$ to optimize the objective functional J, while the control interventions $\left(u_2\right)$ and $\left(u_3\right)$ are kept at zero (i.e., were utilized). It is observed that when viewing the various figures, it is of note to remember that each of the individuals with control intervention(s) are marked by undashed lines and those without control interventions are marked by dotted lines. In Figure 1, we present the plots of population of susceptible humans, infected humans, and vibrio cholerae alongside their control profile. It was observed that the solid lines denote the population of susceptible human, infected human, and vibrio cholerae with control intervention while the dotted lines denote the population of susceptible human, infected human, and vibrio cholerae without control intervention(s). By day 120, cholera is completely wiped out in the infected population in Figure 1a. The infected human population decreases sharply in the first 80 days while the maximum number of recovered humans would increase within those days as well. In Figure 1b, the population of vibrio cholerae increased within the first 20 days but dropped very sharply and wiped out completely until the end of 120 days. It also appeared that the population of vibrio cholerae decreases faster with the control intervention than without the control intervention, in place. In Figure 1c, it was observed that the susceptible human population decreases with the control intervention $u_1$ in place but increases without the control. This shows that the susceptible human population will continue to increase and more people will become more vulnerable to cholera if this control intervention is not put in place. The control profile in Figure 1d shows that the control $u_1$, remained stable from day 1 till the end of the 120 days while the others ($u_2$ and $u_3$) maintained lower bounds.

6.2. Strategy B: Treatment of the Symptomatic Individuals with Drug/Oral Re-Hydration Solution (ORS) $\left(u_2\right)$, Only

Here, we utilized only the control measure $\left(u_2\right)$ to optimize the objective functional J, while the control interventions $\left(u_1\right)$ and $\left(u_3\right)$ are kept at zero (i.e., were not utilized). In Figure 2, we present the plots of population of susceptible humans, infected humans, and vibrio cholerae and their control profile. In Figure 2a, the population of the infected human decreases sharply and by day 20, the disease is wiped out of the system in the presence of the control intervention. Figure 2b showed that the population of vibrio cholerae increases slightly and then dropped from day 20 until the end of 120 days in the presence of the intervention. It was observed in Figure 2c that the population of susceptible human with the control intervention were more than the susceptible human without the control. The control profile in Figure 2d showed that the control $u_2$, remained constant from day 1 till the end of the 120 days while the others ($u_1$ and $u_3$) maintained lower bounds.

6.3. Strategy C: Employing Sanitation, Hygiene, and Safe Water $\left(u_3\right)$, Only

Here, we use only the control measure $\left(u_3\right)$ to optimize the objective functional J, while the control interventions $\left(u_2\right)$ and $\left(u_3\right)$ are kept at zero (i.e., were not utilized). In Figure 3a, it is observed that the population of vibrio cholerae increases slightly and then drops down completely after 80 days and the disease vanishes within 80–120 days. In Figure 3b the number of susceptible humans increases faster (with the control intervention) than without the control intervention. This implies that the control $u_3\left(t\right)$ does not decrease the number of susceptible. In Figure 3c, the number of infected human population decreases completely in day 70 and vanishes by day 120. In Figure 3d, the control profile shows that the control $u_3$, remained constant from day 1 till the end of the 120 days while the others ($u_2$ and $u_3$) maintained lower bounds.

6.4. Strategy D: Employing the Control Interventions ($u_1$, $u_2$)

Here, we used control interventions ($u_1$) and ($u_2$) to optimize the objective functional J, while the control intervention ($u_3$) is kept at zero (i.e., were not utilized). In Figure 4a, the population of the infected human decreases drastically within the first 5 days and cholera was totally wiped out. In Figure 4b, the number of vibrio cholerae increases and then decreases until it vanishes. In Figure 4c, the number of susceptible humans completely decreases and vanishes from the system in the presence of the control intervention but increases without the control interventions. In Figure 4d, the control profile showed that the control ($u_1$) and ($u_2$) remained constant from day 1 till the end of the 120 days while ($u_3$) maintained lower bounds.

6.5. Strategy E: Employing the Control Interventions ($u_1$, $u_3$)

Here, we used control interventions ($u_1$) and ($u_3$) to optimize the objective functional J, while the control intervention ($u_2$) is kept at zero (i.e., were not utilized). In Figure 5a, the population of infected humans decreases essentially after 80 days until the disease was wiped out completely. Figure 5b reveals that the number of vibrio cholerae increases slightly and then decreases completely after 80 days and then vanishes completely. In Figure 5c, the number of susceptible humans completely decreases within 2 days and then vanishes from the system in the presence of the control interventions within 2–120 days. In Figure 5d, the control profile showed that the control ($u_1$), and ($u_3$) remained constant from day 1 till the end of the 120 days while ($u_2$) maintained lower bounds.

6.6. Strategy F: Employing the Control Interventions ($u_2$, $u_3$)

Here, we used control interventions ($u_2$) and ($u_3$) to optimize the objective functional J, while the control intervention ($u_1$) is kept at zero (i.e., were not utilized). In Figure 6a, the population of infected humans decreases within the first 2 days and then decreases until it vanishes completely from the system within 20–120 days. Figure 6b shows that the population of vibrio cholerae increases within the first 2 days and then decreases until it vanishes completely from the system within 20–120 days. In Figure 6c, the susceptible human population inctreases over the time period in the presence of the control intervention but was much slower within the control intervention. In Figure 6d, the control profile shows that the control ($u_2$), and ($u_3$) remained constant from day 1 till the end of the 120 days while ($u_1$) maintained lower bounds.

6.7. Strategy G: Employing All the Three Control Interventions ($u_1$, $u_2$, $u_3$)

Here, we used al the three control interventions ($u_1$), ($u_2$), and ($u_3$) to optimize the objective functional J. In Figure 7a, it was observed that the population of the infected human decreases within the first 2 days and then decreases until it vanishes completely from the system within 2–120 days in the presence of all the three interventions. In Figure 7b, the number of vibrio cholerae increases shortly within the first 20 days and then declines completely until it was wiped out from the system within 20–120 days. In Figure 7c, the number of susceptible humans begins decreasing within the first 2 days until it vanishes completely from the system. In Figure 7d, the control profile shows that the control ($u_1$), ($u_2$), and ($u_3$) remained stable from day 1 till the end of the 120 days.

7. Numerical Simulation of the SDE Model

In this section, we simulate our SDE model (5) by using the method of Milstein given the set of parameters in Table 4. We obtained results of stochastic model (5) for 200 runs, and presented and compared the mean of 200 runs of the stochastic model simulation with the results of the corresponding deterministic model as shown in Figure 8 and Figure 9, where the time series of all the variables of the model were plotted. It is observed from the figures that the mean of 200 runs of the stochastic simulation are very close to the simulation results of the deterministic model.

8. Conclusions

In this article, a nonlinear mathematical model to study the dynamics of the spread of cholera was developed and analyzed. The existence and global stability of the endemic equilibrium were discussed and it was shown that the endemic equilibrium is globally asymptotically stable whenever $R_p>1$. Sensitivity and elasticity analysis of the model were also performed and it was revealed that some parameters had major and significant influence on the spread of cholera than others. The deterministic model was converted to stochastic model to capture some variabilities in the model and the results of the stochastic model were compared with their corresponding deterministic model. It was observed that the level of infected human population was the same with the simulation result of the corresponding deterministic model. It was observed from our work that optimal control applications revealed that the combinations of the interventions in the model had a great influence in the control of cholera. The numerical experiments revealed that optimal strategies to combat effectively the cholera was by the combinations of interventions in the model. The cholera model predicted the control of the disease by the combination of the three control interventions. The models developed predicted the decrease, control, and elimination of cholera morbidity and mortality through incorporating multiple control interventions into the model. Therefore, it is advised that the use of multiple control interventions be adopted for cholera in areas where there are sufficient resources. However, in areas where there are limited or lack of resources, it is advised that treatment of the symptomatic individuals with drug and/or administration of oral re-hydration solution (ORS) to the infected should be used. We therefore advise that in a situation where there are limited resources in the area, strategy B can be applied. The result of our study recommends that: (i) Hygiene promotion and social mobilization should be organized from time to time to reduce the transmission of cholera in the community and (ii) rhe fight against cholera can as well be successful by practising sanitation, hygiene, and drinking safe water. We recognize that this work will be a useful reference for any further study on monitoring the effect of multiple control strategies for cholera disease. We also believe that this study will help society at large to have an understanding of how the disease can be controlled through multiple control strategies in order to reduce the spread of the disease. Finally, the study will also help public health officers and environmental health workers organize seminars, training, and workshops to educate the community on the impact of multiple control strategies in reducing the spread of cholera disease. In future, we will like to study the influence of migration and environmental sanitation on the dynamics of transmission of cholera. The novelty in the work is that we introduced three time-dependent control interventions based on WHO recommendations into the model in order to predict the reduction/control of the cholera in the population as the obtained result may present a good framework for planning and decision making for any national control programs on cholera especially in communities with little or no resources. We plan to extend this work in the future by considering the impact of seasonality and climate variabilities. We also plan to study the influence of host heterogeneities on the transmission of cholera in the presence of multiple control interventions and it will also be interesting to look into the aspect of using real time data to parameterize the model developed.