Stability Analysis and Optimal Control of a Fractional Cholera Epidemic Model

Abstract

In this paper, a fractional model for the transmission dynamics of cholera was developed. In invariant regions of the model, solutions were generated. Disease-free and endemic equilibrium points were obtained. The basic reproduction number was evaluated, and the sensitivity analysis was performed. Under the support of Pontryagin’s maximum principle, the fractional order optimal control was obtained. Furthermore, an optimal strategy was discussed, which minimized the total number of infected individuals and the costs associated with control. Treatment, vaccination, and awareness programs were regarded as three means to reduce the number of infected. Finally, numerical simulations and cost-effectiveness analysis were presented to show the result that the best strategy was the combination of treatment and awareness programs.

1. Introduction

During the last few decades, mathematics models were used to analyze different disease dynamics. An appropriate mathematical model plays a very important role in describing complex biological problems and giving reasonable suggestions. As John Snow studied in 1854, the ingestion of contaminated water can cause an outbreak of cholera [1]. Capasso established a mathematical model with two equations to describe the dynamics of cholera disease in Italy in 1979 [2]. Codeco added an additional equation to the Capasso model and studied a new cholera model [3]. R. L. M. Neilan proposed an SIR model and considered two bacterial concentrations, high and low, respectively, and considered two infectious individuals, asymptomatic individuals and symptomatic individuals [4]. Many researchers have also developed other models of cholera disease transmission since the 1980s [3,5,6,7].

Fractional derivatives and integrals have non-local properties and are very suitable for considering biological systems [8,9,10]. Ref. [11] added the memory effect on the classic SIR model and studied the role of fractional derivatives in disease transmission. Ref. [12] provides a fractional differential model of a class of biological systems with memory. As a result of the memory property of fractional derivative, its theory and application are widely used in the modeling process of the engineering field [13]. Ref. [14] introduced fractional order into an HIV infection model and analyzed the stability of the model in detail. The dynamic model with the Atangana–Baleanu fractional derivative is very effective for studying natural phenomena [15]. The current state and all previous states of the fractional order model determine the next state together [16]. The optimal control problem for the fractional tuberculosis infection model whose fractional derivative is defined in the Atangana–Baleanu–Caputo (ABC) sense including the impact of diabetes and resistant strains is studied [17]. A fractional optimal control of corruption model in the ABC sense based on generalized Mittag–-Leffler is investigated [18]. The fractional COVID-19 epidemic model in the ABC sense is considered [19].

Ref. [20] proposed a model of infectious diseases and took the vaccine into account. Ref. [21] studied the therapeutic effect of a single dose of azithromycin on adults infected with cholera. Ref. [22] summarized the therapeutic effect of different antibiotics on cholera. Ref. [23] studied the control effect of the awareness program on cholera. Ref. [24] explored a combination of treatment and vaccine strategies to control malaria and solved the cost–benefit analysis of different strategies. In order to overcome the epidemic, Ref. [25] incorporated vaccines, treatments, and public health campaigns into their control strategies. In simple terms, we are studying a controlled dynamic system, and in many feasible control schemes, we met the goals the best solution required. FOCPs (fractional order optimal control problems) are generalizations of the classical optimal control problem, and its differential equations are fractional differential equations. In all FOCPs, the author has obtained the necessary conditions for FOCPs optimality.

Cholera is an acute diarrheal infectious disease caused by the contamination of food or water by Vibrio cholerae. Every year, it is estimated that there are 3–5 million cholera cases, resulting in 100,000–120,000 deaths. The peak of the disease is in summer, which can cause diarrhea, dehydration, and even death within a few hours. Most disease outbreaks are caused by Vibrio cholerae O1, while O139, first identified in Bangladesh in 1992, is limited to Southeast Asia. The serotypes of Vibrio cholerae O1 and O139 can cause disease outbreaks. In addition to these two, other kinds of Vibrio cholerae can cause mild diarrhea, but it will not cause an epidemic. Recently, new mutant strains have been found in some parts of Asia and Africa. These strains can cause more serious cholera diseases and higher mortality. Vibrio cholerae exists in water. The most common cause of infection is drinking water contaminated by patients’ feces. Vibrio cholerae can produce cholera toxins and cause secretory diarrhea. It will continue to diarrhea whether eating or not. “Rice water” feces is one of the characteristics of cholera. Therefore, the study of cholera transmission dynamics can predict the spread of cholera and give reasonable countermeasures in time to reduce the loss of people’s lives and property.

In this paper, a control problem is proposed where the control function represents the corresponding cost of preventing the further spread of cholera. One aims to determine the optimal control strategy to minimize costs and maximize benefits. So far, few people have considered the awareness program as a strategy to control cholera. At the same time, few people have combined measures such as vaccines, treatments, and awareness programs. In addition, few people apply fractional order to the optimal control of cholera epidemic models.

The rest of the paper is organized as follows. Section 2 describes the SEIRS fractional epidemic model and briefly reviews the definition of fractional order calculus. Invariant regions of the model are presented. Then, the existence of the solution as well as obtained the disease-free and endemic equilibrium points are proved. The sensitivity analysis of the basic reproduction number $R_0$ is studied. In Section 3, the necessary conditions of the model (2) and the fractional optimal control are derived. In Section 4, numerical simulation of a fractional optimal control problem using different control strategies is presented. In Section 5, a cost-efficiency analysis of the FOCP is provided. Finally, we draw conclusions in Section 6.

In this paper, the parameters are

A:

Susceptible population growth rate;

S:

The amount of susceptible population;

E:

The amount of exposed population;

I:

The amount of infected population;

R:

The amount of recovered population;

$V_I$:

The amount of vibrio cholerae in human intestine;

$V_E$:

The amount of vibrio cholerae in the environment;

a:

Effective inoculation rate;

m:

Mortality of infected people due to illness rate;

u:

The rate of the infected population that receives the treatment;

v:

The rate of the susceptible population who was vaccinated;

w:

The rate of promoted awareness;

$\alpha$:

The growth rate of Vibrio Cholerae in human intestine;

$\psi$:

The shedding rate of Vibrio Cholerae in human intestine;

$\beta$:

The rate of susceptible people becoming exposed;

${\gamma }_1$:

Natural recovery rate;

$\delta$:

Recovered crowd loses immunity rate;

$\sigma$:

The rate of exposed population becoming infected;

$\chi$:

The rate of virus eliminate from human intestine;

$\lambda$:

Recovery rate;

${\mu }_d$:

Natural mortality rate of individuals;

${\mu }_{V_E}$:

Natural mortality rate of Vibrio cholerae.

2. Preliminaries and Model Description

In this section, some basics of fractional calculus and some necessary lemmas are given; meanwhile, the factional system is introduced.

2.1. Model Description

Prabir proposed a SEIRS cholera transmission model to study the transmission of cholera among people [26]. It is well known that the incubation period plays a crucial role in controlling the cholera. One divides Vibrio cholerae into Vibrio cholerae in the intestine and the environment. Suppose that the susceptible population becomes exposed because of insufficient cognition. After vaccinating, the proportion of susceptible people who become recovered people is $av$. The increase in the number of people recovering is due to the recovery of infected individuals. There are two reasons: one is natural recovery and the other is treatment. The rate of treatment is $u{\gamma }_1(\lambda -1)$. The number of Vibrio cholerae in human intestines increased because of the treatment of infected people, but the treatment rate is only $\psi u\alpha$, and the rate of no cure is $(1-u)\alpha$, relatively. The model is given by [26]

$$\left\{\begin{array}{cc}\hfill & \frac{dS}{dt}=A-{\mu }_{d}S-(1-w)\beta S{V}_E+\delta R-avS,\hfill \\ \hfill & \frac{dE}{dt}=(1-w)\beta S{V}_E-{\mu }_{d}E-\sigma E,\hfill \\ \hfill & \frac{dI}{dt}=\sigma E-{\mu }_{d}I-mI-\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]I,\hfill \\ \hfill & \frac{dR}{dt}=\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]I-{\mu }_{d}R-\delta R+avS,\hfill \\ \hfill & \frac{d{V}_I}{dt}=[\psi u+(1-u)]\alpha I-\chi V_I,\hfill \\ \hfill & \frac{d{V}_E}{dt}=\chi V_I-{\mu }_{V_E}{V}_E,\hfill \end{array}\right.$$

(1)

where $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I\left(0\right)\ge 0,R\left(0\right)\ge 0,V_I\left(0\right)\ge 0,V_E\left(0\right)\ge 0$.

The total population is $N_p\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$. The total number of Vibrio cholerae is $N_c\left(t\right)=V_I\left(t\right)+V_E\left(t\right)$. The above model does not include the memory effect, and the memory effect is very important for accurately describing the biodynamic model. Atangana and Baleanu introduced the $ABC$ operator. This operator describes the memory effect more effectively, the model with $ABC$ operator ${}_0^{ABC}{D}_t^{\alpha }$ is given as follows:

$$\left\{\begin{array}{cc}\hfill & {}_0^{ABC}{D}_t^{\alpha }S=A-{\mu }_{d}S-(1-w)\beta S{V}_E+\delta R-avS,\hfill \\ \hfill & {}_0^{ABC}{D}_t^{\alpha }E=(1-w)\beta S{V}_E-{\mu }_{d}E-\sigma E,\hfill \\ \hfill & {}_0^{ABC}{D}_t^{\alpha }I=\sigma E-{\mu }_{d}I-mI-\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]I,\hfill \\ \hfill & {}_0^{ABC}{D}_t^{\alpha }R=\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]I-{\mu }_{d}R-\delta R+avS,\hfill \\ \hfill & {}_0^{ABC}{D}_t^{\alpha }{V}_I=[\psi u+(1-u)]\alpha I-\chi V_I,\hfill \\ \hfill & {}_0^{ABC}{D}_t^{\alpha }{V}_E=\chi V_I-{\mu }_{V_E}{V}_E,\hfill \end{array}\right.$$

(2)

where $S\left(t\right)=S_0\left(t\right),E\left(t\right)=E_0\left(t\right),I\left(t\right)=I_0\left(t\right),R\left(t\right)=R_0\left(t\right),V_I\left(t\right)=V_{I0}\left(t\right),V_E\left(t\right)=V_{E0}\left(t\right)$.

2.2. Preliminaries

In this section, some preliminaries of fractional derivatives are introduced. The fractional derivative is a generalization of the integer derivative. Several definitions are provided here, such as the Riemann Liouville and Caputo fractional derivative. Start with the Liouville–Caputo fractional derivative.

Definition 1.

The fractional derivative of the Atangana–Baleanu type in the Liouville–Caputo sense ($ABC$) is defined as [10]

$${}_0^{ABC}{D}_t^{\alpha }\left\{f\left(t\right)\right\}=\frac{B\left(\alpha \right)}{1-\alpha }{\int }_0^t{E}_{\alpha }\left[-\alpha \frac{{(t-\theta )}^{\alpha }}{1-\alpha }\right]\dot{f}\left(\theta \right)d\theta ,$$

where $0<\alpha <1$. $B\left(\alpha \right)$ is the normalization constant and satisfies $B\left(0\right)=B\left(1\right)=1$. The one-parameter Mittag–Leffler functions $E_{\alpha }\left(z\right)$ is defined as

$$E_{\alpha }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }(\alpha k+1)},$$

where $z\in \mathbb{C},\Re \left(\alpha \right)>0$. $\mathrm{\Gamma }(·)$ is the gamma function, $\mathrm{\Gamma }\left(x\right)={\int }_0^{+\infty }{t}^{x-1}{e}^{-t}\mathrm{d}t$, and satisfies $\mathrm{\Gamma }(x+1)=x\mathrm{\Gamma }\left(x\right)$.

Definition 2.

Riemann–Liouville fractional order integral (p order) is defined as

$${}_a^{RL}{D}_t^{-p}f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(p\right)}{\int }_a^t{(t-\tau )}^{p-1}f\left(\tau \right)\mathrm{d}\tau ,p>0.$$

Definition 3.

The Laplace transform of the Riemann–Liouville fractional derivative is

$$L\left\{{}_0^{RL}{D}_t^{p}f\left(t\right);s\right\}=s^{p}F\left(s\right)-\sum _{k=0}^{n-1}{s}^k{\left[{}_0^{RL}{D}_t^{p-k-1}f\left(t\right)\right]}_{t=0},$$

(3)

where $n-1⩽p<n$.

Definition 4.

The Laplace transform of the Caputo fractional derivative is

$$L\left\{{}_0^C{D}_t^{p}f\left(t\right);s\right\}=s^{p}F\left(s\right)-\sum _{k=0}^{n-1}{s}^{p-k-1}{f}^{\left(k\right)}\left(0\right),$$

(4)

where $n-1<p⩽n$.

2.3. Invariant Region

In this section, one will discuss the boundary problem of the model. The concept of invariant region is very important in the dynamic system and is the key research object. The positive invariant region is for any set D, and the solution on D is still in D.

The model (2) contains two parts: the human and the Vibrio cholerae. Now, one will consider the feasible region of model (2).

$$D=D_p\cup D_{V_f}\cup D_{V_E}\subset R_{+}^4\times R_{+}^1\times R_{+}^1,$$

with

$$D_p=\left\{(S,E,I,R)\in R_{+}^4:S+E+I+R\le \frac{{\mathrm{\Lambda }}_l}{{\mu }_l}\right\},$$

$$D_{V_I}=\left\{V_I\in R_{+}^1:V_I\le \frac{[\psi u+(1-u\left)\right]\alpha A}{\chi {\mu }_d}\right\},$$

$$D_{V_E}=\left\{V_E\in R_{+}^1:V_E\le \frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_{V_E}{\mu }_d}\right\}.$$

In order to obtain an invariant region, one needs to perform the following five steps.

Step 1: Adding the first four equations of model (2) gives:

$${}_0^{ABC}{D}_t^{\alpha }{N}_p\left(t\right)=A-{\mu }_d{N}_p\left(t\right)-mI.$$

Step 2: By definition, one can know $mI\le 0$; then

$${}_0^{ABC}{D}_t^{\alpha }{N}_p\left(t\right)\le A-{\mu }_d{N}_p\left(t\right).$$

(5)

Step 3: After Laplace transform and inverse transform (4) on the above Formula (5), one can observe

$$N_p\left(t\right)\le \left(N_p\left(0\right)-\frac{A}{{\mu }_d}\right)E_{\alpha }\left(-{\mu }_d{t}^{\alpha }\right)+\frac{A}{{\mu }_d},$$

Calculating the above with $N_p\left(0\right)\le \frac{A}{{\mu }_d}$ and $E_{\alpha }\left(-{\mu }_a{t}^{\alpha }\right)\ge 0$, one has

$$N_p\left(t\right)\le \frac{A}{{\mu }_d}.$$

(6)

Step 4: In order to eliminate ambiguity, one will temporarily record $\alpha$ in the fractional order equation as $\beta$ and bring the Formula (6) into it. According to the above analysis, one can obtain

$${}_0^{ABC}{D}_t^{\beta }{V}_I\left(t\right)=[\psi u+(1-u)]\alpha I\left(t\right)-\chi V_I\left(t\right)\le \frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_d}-\chi V_I\left(t\right).$$

After Laplace transform and inverse transform on the above formula, one can observe

$$V_I\left(t\right)\le \left(V_I\left(0\right)-\frac{[\psi u+(1-u\left)\right]\alpha A}{\chi {\mu }_d}\right)E_{\beta }\left(-\chi t^{\beta }\right)+\frac{[\psi u+(1-u\left)\right]\alpha A}{\chi {\mu }_d},$$

After doing the same operation as Step 3, it gives

$$V_I\left(t\right)\le \frac{[\psi u+(1-u\left)\right]\alpha A}{\chi {\mu }_d}.$$

(7)

Step 5: According to the above Formula (7), one can obtain

$${}_0^{ABC}{D}_t^{\beta }{V}_E=\chi V_I-{\mu }_{V_E}{V}_E\le \frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_d}-{\mu }_{V_E}{V}_E,$$

After Laplace transform and inverse transform on the above formula, one can observe

$$V_E\left(t\right)\le \left(V_E\left(0\right)-\frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_{{\mu }_g}{\mu }_d}\right)E_{\beta }\left(-{\mu }_{V_s}{t}^{\beta }\right)+\frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_{V_s}{\mu }_d},$$

It can easily be

$$V_E\left(t\right)\le \frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_{V_S}{\mu }_d}.$$

So far, one has successfully obtained an invariant region:

$$\begin{array}{c}D=\left\{\left(S,E,I,R,V_I,V_E\right)\in R_{+}^6:S+E+I+R\le \frac{{\mathrm{\Lambda }}_l}{{\mu }_l},V_I\left(t\right)\le \frac{[\psi u+(1-u\left)\right]\alpha A}{\chi {\mu }_d},V_E\left(t\right)\le \frac{[\psi u+(1-u\left)\right]\alpha A}{{\mu }_{V_E}{\mu }_d}\right\}\end{array}.$$

(8)

Lemma 1.

If the initial condition is non-negative, the region D (8) of the model (2) is an invariant region on $R_{+}^6$.

3. Existence of the Solution

In this section, one will investigate the existence of the solutions for model (2) by using the fixed points theorem.

3.1. Disease-Free Equilibrium Point

It is easy to see that the disease-free equilibrium points are:

$$E_0=\left(S^0,E^0,I^0,R^0,V_I^0,V_E^0\right)=\left(\frac{A\left({\mu }_d+\delta \right)}{{\mu }_d\left({\mu }_d+av+\delta \right)},0,0,\frac{Aav}{{\mu }_d\left({\mu }_d+av+\delta \right)},0,0\right),$$

and the endemic equilibrium points are

$$S^{*}=\frac{{\mu }_{V_E}\left({\mu }_d+\sigma \right)\left({\mu }_d+m+\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]\right)}{\sigma \beta (1-w)[\psi u+(1-u\left)\right]\alpha },$$

$$E^{*}=\frac{\left({\mu }_d+m+\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]\right)I^{*}}{\sigma },$$

$$I^{*}=\frac{\left({\mu }_d+m+\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]\right){\mu }_{{\gamma }_s}\left({\mu }_d^2+{\mu }_d\delta +av{\mu }_d\right)}{\alpha \beta \left(1-w\left)\right({\mu }_d+\sigma \right)\left({\mu }_d+\delta \right)\left({\mu }_d+m\right)(\psi u+1-u)\left(m{u}_d^2+{\mu }_d\delta +{\mu }_d\sigma \right)\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]},$$

$$V_I^{*}=\frac{\alpha [\psi u+(1-u)]I^{*}}{\chi },$$

$$V_E^{*}=\frac{\alpha [\psi u+(1-u)]I^{*}}{{\mu }_{V_E}}.$$

In the ongoing study of infectious diseases, $R_0$ is defined as the average number of infections produced by an infected person during the period of infection.

Theorem 1.

Disease-free equilibrium point $E_0$ is locally asymptotically stable at $R_0<1$. Endemic equilibrium point $E_{*}$ is locally asymptotically stable at $R_0>1$.

Proof. 

The Jacobi matrix $J\left(E_0\right)$ of disease-free equilibrium points for model (2) is

$$J\left(E_0\right)=\left(\begin{array}{cccccc}-{\mu }_d-av& 0& 0& \delta & 0& -\frac{\beta A(1-\omega )\left({\mu }_d+\delta \right)}{{\mu }_d\left({\mu }_d+av+\delta \right)}\\ 0& -\left({\mu }_d+\sigma \right)& 0& 0& 0& \frac{\beta A(1-\omega )\left({\mu }_d+\delta \right)}{{\mu }_d\left({\mu }_d+av+\delta \right)}\\ 0& \sigma & -\left({\mu }_d+m+\left((1-u){\gamma }_1+u{\gamma }_1\lambda \right)\right)& 0& 0& 0\\ av& 0& 0& -\left({\mu }_d+\delta \right)& 0& 0\\ 0& 0& [\psi u+(1-u\left)\right]\alpha & 0& -\chi & 0\\ 0& 0& 0& 0& \chi & -{\mu }_{V_E}\end{array}\right)$$

The eigenvalue of $J\left(E_0\right)$ is difficult to solve, so one will solve its eigenequation. $a_{ij}$ denotes the elements of row i column j of $J\left(E_0\right)$. The characteristic equation is abbreviated as

$$T_6{X}^6+T_5{X}^5+T_4{X}^4+T_3{X}^3+T_2{X}^2+T_1{X}^1+T_0=0,$$

It is calculated that $T_i>0(i=1,2,3,4,5,6)$. The Routh table of the above equation is

$$\begin{array}{cc}& x^6{T}_6{T}_4{T}_2{T}_0\hfill \\ & x^5{T}_5{T}_3{T}_{1}0\hfill \\ & x^4{S}_1{S}_2{S}_3\hfill \\ & x^3{R}_1{R}_{2}0\hfill \\ & x^2{Q}_1{Q}_2\hfill \\ & x^1{P}_{1}0\hfill \\ & x^0{O}_1\hfill \end{array}$$

After analysis, it can be concluded that the first column elements of the Routh table of the above polynomial and the ordinal principal subscripts are all positive. Thus, based on the Routh–Hurwitz discriminant, it is possible to conclude that the disease-free equilibrium point is locally asymptotically stable. □

3.2. $R_0$ Sensitivity Analysis

Here, one-at-a-time (OAT) is used to test parameter sensitivity, and the effect of variables on results is obtained by using partial derivatives of variables. It is easy to find that

$$R_0=\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B},$$

Let us evaluate the sensitivity of $R_0$ to each parameter.

$$\frac{\partial R_0}{\partial A}=\frac{\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}>0,$$

$$\frac{\partial R_0}{\partial \alpha }=\frac{A\beta \sigma (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}>0,$$

$$\frac{\partial R_0}{\partial \beta }=\frac{A\sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}>0,$$

$$\frac{\partial R_0}{\partial w}=-\frac{A\beta \sigma \alpha \left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}<0,$$

$$\frac{\partial R_0}{\partial \psi }=\frac{A\beta \sigma \alpha u(1-w)\left({\mu }_d+\delta \right)}{{\mu }_{V_E}B}>0,$$

$$\frac{\partial R_0}{\partial \sigma }=\frac{A\beta \alpha {\mu }_d{\mu }_{V_E}(1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}^2({\mu }_d+\sigma )B}>0,$$

$$\frac{\partial R_0}{\partial \delta }=\frac{A\beta \sigma \alpha av{\mu }_d{\mu }_{V_E}(1-w)[\psi u+(1-u)]}{{\mu }_{V_E}^{2}B({\mu }_d^2+av{\mu }_d+{\mu }_d\delta )}>0,$$

$$\frac{\partial R_0}{\partial v}=-\frac{A\beta \sigma \alpha a{\mu }_d(1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}\left({\mu }_d^2+\alpha v{\mu }_d+{\mu }_d\delta \right)B}<0,$$

$$\frac{\partial R_0}{\partial m}=-\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}{B}^2}<0,$$

$$\frac{\partial R_0}{\partial u}=\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)C}{{\mu }_{V_E}{B}^2}<0,$$

$$\frac{\partial R_0}{\partial {\mu }_d}=\frac{A\beta \sigma \alpha (1-w)[\psi u+(1-u\left)\right]D}{{\mu }_{V_E}B}>0,$$

$$\frac{\partial R_0}{\partial {\mu }_{V_E}}=-\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}^{2}B}<0,$$

where

$$B=\left({\mu }_d^2+av{\mu }_d+{\mu }_d\delta \right)\left({\mu }_d+\sigma \right)\left({\mu }_d+m+\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]\right)>0,$$

$$C=(\psi -1)\left({\mu }_d^2+av{\mu }_d+{\mu }_d\delta \right){\mu }_{V_E}\left({\mu }_d+\sigma \right)+{\gamma }_1(\psi -\lambda )<0,$$

$$D=1-\frac{\left({\mu }_d+\delta \right)}{{\mu }_d+m+[(1-u){\gamma }_1+u{\gamma }_1\lambda ]}-\frac{({\mu }_d+\delta )}{{\mu }_d+\sigma }>0.$$

Hence, $R_0$ is increasing with $A,\alpha ,\beta ,\psi ,\sigma ,\delta$ and is decreasing with $w,v,m,u,{\mu }_d,{\mu }_{V_E}$.

4. Fractional Optimal Control

Here, one will discuss the fractional optimal control problem ($FOCP$) and the optimal control conditions.

4.1. Fractional Optimal Control Problem

We aim to reduce the number of infected individuals and at the same time to reduce the cost of the treatment (u), vaccine (v), and awareness program (w). Now, the feasible control function of the above model on $R^6$ is considered.

$$\mathsf{\Psi }=\left\{\left(u(·),v(·),w(·)\right)\in \left(L^{\infty }(0,T)\right),0\le u(·),v(·),w(·)\le 1,\forall t\in [0,T],\right\}$$

Here, one takes into account the overdose of drugs. Therefore, the square of the respective cost is used [27]. The objective function is defined as

$$J(u,v,w)={\int }_0^T\left(p_{0}I\left(t\right)+p_1{u}^2\left(t\right)+p_2{v}^2\left(t\right)+p_3{w}^2\left(t\right)\right)dt,$$

(9)

where $p_0,p_1,p_2,$and$p_3$ are the relative costs of stopping the spread of cholera. One attempts to find an optimal control $u^{*},v^{*},w^{*}$ in order to obtain

$$J\left(u^{*},v^{*},w^{*}\right)=\underset{\mathsf{\Psi }}{min}{\int }_0^T\left(p_{0}I\left(t\right)+p_1{u}^2\left(t\right)+p_2{v}^2\left(t\right)+p_3{w}^2\left(t\right)\right)dt,$$

Next, the Lagrangian is given

$$L(I,u,v,w)=p_{0}I\left(t\right)+p_1{u}^2\left(t\right)+p_2{v}^2\left(t\right)+p_3{w}^2\left(t\right),$$

Furthermore, the Hamiltonian of the model is $H(S,E,I,R,V_I,V_E,u,v,w,{\lambda }_S,{\lambda }_E,{\lambda }_I,{\lambda }_R,{\lambda }_{V_I},{\lambda }_{V_E})$, such that ${\lambda }_S,{\lambda }_E,{\lambda }_I,{\lambda }_R,{\lambda }_{V_I}$ and ${\lambda }_{V_E}$ are the adjoint representation. Then, substituting model (2) into H, one can obtain

$$\begin{array}{cc}\hfill H=& p_{0}I\left(t\right)+p_1{u}^2\left(t\right)+p_2{v}^2\left(t\right)+p_3{w}^2\left(t\right)+{\lambda }_S\left\{A-(1-w)\beta S{V}_E-{\mu }_{d}S+\delta R-avS\right\}\hfill \\ & +{\lambda }_E\left\{(1-w)\beta S{V}_E-{\mu }_{d}E-\sigma E\right\}\hfill \\ & +\lambda \left\{\sigma E-{\mu }_{d}I-mI-\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]I\right\}\hfill \\ & +{\lambda }_R\left\{\left[(1-u){\gamma }_1+u{\gamma }_1\lambda \right]I-{\mu }_{d}R-\delta R+avS\right\}\hfill \\ & +{\lambda }_{V_I}\left\{[\psi u+(1-u)]\alpha I-\chi V_I\right\}+{\lambda }_{V_B}\left\{\chi V_I-{\mu }_{V_B}{V}_E\right\},\hfill \end{array}$$

From the above, one can derive

$${}_0^{ABC}{D}_t^{\alpha }{\lambda }_S=-\frac{\partial H}{\partial S},{}_0^{ABC}{D}_t^{\alpha }{\lambda }_E=-\frac{\partial H}{\partial E},{}_0^{ABC}{D}_t^{\alpha }{\lambda }_I=-\frac{\partial H}{\partial I},$$

$${}_0^{ABC}{D}_t^{\alpha }{\lambda }_R=-\frac{\partial H}{\partial R},{}_0^{ABC}{D}_t^{\alpha }{\lambda }_{V_I}=-\frac{\partial H}{\partial V_I},{}_0^{ABC}{D}_t^{\alpha }{\lambda }_{V_E}=-\frac{\partial H}{\partial V_E},$$

moreover

$$\frac{\partial H}{\partial u}=0,\frac{\partial H}{\partial v}=0,\frac{\partial H}{\partial w}=0.$$

In addition, there is a necessary condition that ${\lambda }_S,{\lambda }_E,{\lambda }_I,{\lambda }_R,{\lambda }_{V_I},$and${\lambda }_{V_E}$ are zero.

So far, the above equation gives the necessary conditions for the previously defined Hamiltonian of the $FOCP$, resulting in different sets of differential equations, state variables $(S,E,I,R,V_I,V_E)$, controls $(u,v,w)$ and the Lagrangian.

4.2. Fractional Optimal Control Conditions

Theorem 2.

Model (2) has an optimal control such that the objective function on Ψ is minimal and its accompanying variables satisfy the following equation.

$$\begin{array}{cc}\hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_S\left(t\right)& =\left(\left({\mu }_d+(1-w)\right)\beta V_E+\alpha v\right){\lambda }_S-(1-w)\beta V_E{\lambda }_E-{\lambda }_{R}av,\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_E\left(t\right)& =\left({\mu }_d+\sigma \right){\lambda }_E-{\lambda }_I\sigma ,\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_l\left(t\right)& =-p_0+{\lambda }_l\left({\mu }_d+m+\left({\gamma }_1+u{\gamma }_1(\lambda -1)\right)\right)-{\lambda }_R\left({\gamma }_1+u{\gamma }_1(\lambda -1)\right)-\alpha {\lambda }_{V_I}(1+u(\psi -1)),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_R\left(t\right)& =\left({\mu }_d+\delta \right){\lambda }_R-\delta {\lambda }_S,\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_{V_I}\left(t\right)& =\chi \left({\lambda }_{V_I}-{\lambda }_{V_E}\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_{V_E}\left(t\right)& =(1-w)\beta S\left({\lambda }_S-{\lambda }_E\right)+{\lambda }_{V_E}{\mu }_{V_E,}\hfill \end{array}$$

with transversal conditions

$${\lambda }_S\left(T\right)=0,{\lambda }_E\left(T\right)=0,{\lambda }_I\left(T\right)=0,{\lambda }_R\left(T\right)=0,{\lambda }_{V_I}\left(T\right)=0,{\lambda }_{V_E}\left(T\right)=0,$$

and the optimal controllers

$$\begin{array}{cc}\hfill u^{*}& =min\left\{max\left\{\frac{\left(\left({\gamma }_1-{\gamma }_1\lambda \right)\left({\lambda }_R-{\lambda }_I\right)+\alpha (1-\psi ){\lambda }_{V_I}\right)I}{2p_1},0\right\},1\right\},\hfill \\ \hfill v^{*}& =min\left\{max\left\{\frac{a\left({\lambda }_S-{\lambda }_R\right)S}{2p_2},0\right\},1\right\},\hfill \\ \hfill w^{*}& =min\left\{max\left\{\frac{\beta \left({\lambda }_E-{\lambda }_S\right)S{V}_E}{2p_3},0\right\},1\right\}.\hfill \end{array}$$

Proof. 

With the help of the Pontryagin maximum principle, one could obtain the adjoint equation and transversal condition [17,18,19]. The necessary conditions for the optimality of Equation (9) are model (2) and

$$\begin{array}{cc}\hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_S\left(t\right)& =-\frac{\partial H}{\partial S}=\left(\left({\mu }_d+(1-w)\right)\beta V_E+av\right){\lambda }_S-(1-w)\beta V_E{\lambda }_E-{\lambda }_{R}av,\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_E\left(t\right)& =-\frac{\partial H}{\partial E}=\left({\mu }_d+\sigma \right){\lambda }_E-{\lambda }_I\sigma ,\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_l\left(t\right)& =-\frac{\partial H}{\partial I}=-p_0+{\lambda }_l\left({\mu }_d+m+\left({\gamma }_1+u{\gamma }_1(\lambda -1)\right)\right)-{\lambda }_R\left({\gamma }_1+u{\gamma }_1(\lambda -1)\right)-\alpha {\lambda }_{V_I}(1+u(\psi -1)),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_R\left(t\right)& =-\frac{\partial H}{\partial R}=\left({\mu }_d+\delta \right){\lambda }_R-\delta {\lambda }_S,\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_{V_I}\left(t\right)& =-\frac{\partial H}{\partial V_I}=\chi \left({\lambda }_{V_I}-{\lambda }_{V_E}\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{\lambda }_{V_E}\left(t\right)& =-\frac{\partial H}{\partial V_E}=(1-w)\beta S\left({\lambda }_S-{\lambda }_E\right)+{\lambda }_{V_E}{\mu }_{V_E},\hfill \end{array}$$

The initial conditions are

$$S\left(0\right)=S_0,E\left(0\right)=E_0,I\left(0\right)=I_0,R\left(0\right)=R_0,V_I\left(0\right)={V_I}_0,V_E\left(0\right)={V_E}_0,$$

It is obtained from the condition

$$\frac{\partial H}{\partial u}=0,\frac{\partial H}{\partial v}=0,\frac{\partial H}{\partial w}=0,$$

that

$$\begin{array}{cc}\hfill u^{*}& =\frac{\left(\left({\gamma }_1-{\gamma }_1\lambda \right)\left({\lambda }_R-{\lambda }_I\right)+\alpha (1-\psi ){\lambda }_{V_I}\right)I}{2p_1},\hfill \\ \hfill v^{*}& =\frac{a\left({\lambda }_S-{\lambda }_R\right)S}{2p_2},\hfill \\ \hfill w^{*}& =\frac{\beta \left({\lambda }_E-{\lambda }_S\right)S{V}_E}{2p_3},\hfill \end{array}$$

After applying the boundary conditions to the controllers, optimal control can be written in the following form

$$\begin{array}{cc}\hfill u^{*}& =min\left\{max\left\{\frac{\left(\left({\gamma }_1-{\gamma }_1\lambda \right)\left({\lambda }_R-{\lambda }_I\right)+\alpha (1-\psi ){\lambda }_{V_I}\right)I}{2p_1},0\right\},1\right\},\hfill \\ \hfill v^{*}& =min\left\{max\left\{\frac{a\left({\lambda }_S-{\lambda }_R\right)S}{2p_2},0\right\},1\right\},\hfill \\ \hfill w^{*}& =min\left\{max\left\{\frac{\beta \left({\lambda }_E-{\lambda }_S\right)S{V}_E}{2p_3},0\right\},1\right\}.\hfill \end{array}$$

□

5. Numerical Simulation

In this section, the main purpose is to provide the solution to the practical model above, using the Adams–Moulton numerical method of the Atangana–Baleanu fractional integration, while the values of this part are given in the paper [26].

$$\begin{array}{cc}\hfill {}_0^{ABC}{D}_t^{\alpha }S& =f_1\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right),u\left(t_j\right),v\left(t_j\right),w\left(t_j\right),t\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }E& =f_2\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right),u\left(t_j\right),v\left(t_j\right),w\left(t_j\right),t\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }I& =f_3\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right),u\left(t_j\right),v\left(t_j\right),w\left(t_j\right),t\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }R& =f_4\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right),u\left(t_j\right),v\left(t_j\right),w\left(t_j\right),t\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{V}_I& =f_5\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right),u\left(t_j\right),v\left(t_j\right),w\left(t_j\right),t\right),\hfill \\ \hfill {}_0^{ABC}{D}_t^{\alpha }{V}_E& =f_6\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right),u\left(t_j\right),v\left(t_j\right),w\left(t_j\right),t\right),\hfill \end{array}$$

where $0<\alpha <1$ and $t>1$. Suppose one finds that the interval of the solution to the problem is $[0,n]$. The interval $[0,n]$ is divided into N subintervals, each of which has length $h=\frac{n}{N}$. The nodes are $t_j=jh$ where $j=0,1,2,\cdots N,$ and $t_{j+1}=t_j+h$. Thus, the generalized formula is

$$\begin{array}{cc}\hfill S\left(t_{j+1}\right)=S\left(t_j\right)& +\frac{1-\alpha }{B\left(\alpha \right)}{f}_1+\frac{\alpha }{B\left(\alpha \right)}\sum _{N=0}^{\infty }{(N+1)}^{1-\alpha }-{\left(N\right)}^{1-\alpha }{f}_1,\hfill \\ \hfill E\left(t_{j+1}\right)=E\left(t_j\right)& +\frac{1-\alpha }{B\left(\alpha \right)}{f}_2+\frac{\alpha }{B\left(\alpha \right)}\sum _{N=0}^{\infty }{(N+1)}^{1-\alpha }-{\left(N\right)}^{1-\alpha }{f}_2,\hfill \\ \hfill I\left(t_{j+1}\right)=I\left(t_j\right)& +\frac{1-\alpha }{B\left(\alpha \right)}{f}_3+\frac{\alpha }{B\left(\alpha \right)}\sum _{N=0}^{\infty }{(N+1)}^{1-\alpha }-{\left(N\right)}^{1-\alpha }{f}_3,\hfill \\ \hfill R\left(t_{j+1}\right)=R\left(t_j\right)& +\frac{1-\alpha }{B\left(\alpha \right)}{f}_4+\frac{\alpha }{B\left(\alpha \right)}\sum _{N=0}^{\infty }{(N+1)}^{1-\alpha }-{\left(N\right)}^{1-\alpha }{f}_4,\hfill \\ \hfill V_I\left(t_{j+1}\right)=V_I\left(t_j\right)& +\frac{1-\alpha }{B\left(\alpha \right)}{f}_5+\frac{\alpha }{B\left(\alpha \right)}\sum _{N=0}^{\infty }{(N+1)}^{1-\alpha }-{\left(N\right)}^{1-\alpha }{f}_5,\hfill \\ \hfill V_E\left(t_{j+1}\right)=V_E\left(t_j\right)& +\frac{1-\alpha }{B\left(\alpha \right)}{f}_6+\frac{\alpha }{B\left(\alpha \right)}\sum _{N=0}^{\infty }{(N+1)}^{1-\alpha }-{\left(N\right)}^{1-\alpha }{f}_6,\hfill \end{array}$$

where $f_i=f_i\left(S\left(t_j\right),E\left(t_j\right),I\left(t_j\right),R\left(t_j\right),V_I\left(t_j\right),V_E\left(t_j\right)\right),i=1,2,\cdots ,6$.

Figure 1 is based on Parameter

Table 1. Parameter valuation.

Parameter	Estimated Value	Parameter	Estimated Value
$\alpha$	0.73 [28]	A	50 day${}^{-1}$ [29]
$\beta$	0.01 [26]	a	0.6 [30]
$\sigma$	0.6 [26]	$\chi$	6 cell per infected human [26]
$\delta$	0.001 [4]	${\mu }_d$	0.2 [29]
${\gamma }_1$	0.52 [30]	${\mu }_{V_E}$	1/30 day${}^{-1}$ [28]
$\lambda$	2.3 [30]	$p_1$	2$ per infected person [4]
$\psi$	0.52 [30]	$p_2$	6$ per person [4]
m	0.005 [30]	$p_3$	1$ per person [26]

, which shows visually the local asymptotic stability of disease-free equilibrium. In Figure 1, it is clear that $S,E,I,R,V_I,$and$V_E$ are from top to bottom in which $t=50$ days, while the population label represents the amount of corresponding persons.

Next, the optimal control specifically will be applied to our model. Our goal is to reduce the number of infected populations while using different combinations of treatment, vaccine, and awareness campaigns. Four different control strategies are considered here.

Let $u,v,w=0$. The initial condition is $(S_0,E_0,I_0,R_0,V_{I0},V_{E0})=(50,100,50,90,40,30)$. One has

$$R_0={\left.\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}\right|}_{(u,v,w)=(0,0,0)}=56.6379>1.$$

(a)

Using the combination of vaccination and awareness programs only (${\mathit{A}}_{\mathbf{1}}$);

Here, v (vaccines) and w (awareness programs) are optimized to obtain the minimum value of the objective function J. Here, let $v=1,w=0.5$. That is, one supposes that $u=0$ (treatment) here. In Figure 2, one can see that after applying these measures, all populations have their own variations. The susceptible population is decreasing, the exposed population is increasing, the infected population is increasing, and the recovered population is decreasing. Vibrio cholerae is increasing in the human gut as well as in the environment. The basic reproduction number is

$$R_0={\left.\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}\right|}_{(u,v,w)=(0,1,0.5)}=7.1063>1,$$

(b)

Using the combination of vaccination and treatment only (${\mathit{A}}_{\mathbf{2}}$);

In this strategy, only two controls v (vaccines) and u (treatments) are used to obtain the minimum value of the objective function J. Here, let $v=1,u=1$. That is, one supposes that $w=0$ (awareness programs) here. In Figure 3, one can see that after applying these measures, all populations have their own variations. The susceptible population is increasing, the exposed population is decreasing, the infected population is decreasing, and the recovered population is increasing. Vibrio cholerae is decreasing in the human gut as well as in the environment. The basic reproduction number is

$$R_0={\left.\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}\right|}_{(u,v,w)=(1,1,0)}=3.8245>1,$$

(c)

Using the combination of treatment and awareness programs only (${\mathit{A}}_{\mathbf{3}}$);

Here, u (treatments) and w (awareness programs) are optimized to obtain the minimum value of the objective function J. Here, let $u=1,w=0.5$. That is, one supposes that $v=0$ (vaccines) here. In Figure 4, one can see that after applying these measures, all populations have their own variations. The susceptible population is increasing, the exposed population is decreasing, the infected population is decreasing, and the recovered population is increasing. Vibrio cholerae is decreasing in the human gut as well as in the environment. The basic reproduction number is

$$R_0={\left.\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}\right|}_{(u,v,w)=(1,0,0.5)}=7.6204>1,$$

(d)

Using the combination of treatment, vaccination, and awareness programs (${\mathit{A}}_{\mathbf{4}}$);

In this strategy, all controllers u (treatments), v (vaccines), and w (awareness programs) are optimized to obtain the minimum value of the objective function J. Here, let $u=0.5,v=0.5,w=0.5$. In Figure 5, one can see that after applying these measures, all populations have their own variations. The susceptible population is increasing, the exposed population is decreasing, the infected population is decreasing, and the recovered population is increasing. Vibrio cholerae is decreasing in the human gut as well as in the environment. The basic reproduction number is

$$R_0={\left.\frac{A\beta \sigma \alpha (1-w)\left({\mu }_d+\delta \right)[\psi u+(1-u)]}{{\mu }_{V_E}B}\right|}_{(u,v,w)=(0.5,0.5,0.5)}=5.8892>1.$$

Cost-Effectiveness Analysis

CEA (Cost-Effectiveness Analysis) is an economic analysis of cost, which can compare the relative cost of two or more schemes. It is often used in the medical field. In this section, the most effective strategy is identified for using one of the cholera control strategies. It is very important to compare the results of different measures, with the help of calculating the ICER (Incremental Cost Effectiveness Ratio) to seek an answer [31]. ICER is often described as the additional cost of increasing medical outcomes. Based on the numerical simulation results of the model (2), its incremental cost-effectiveness ratio is obtained.

The ICERs of the four strategies are calculated as follows

$$ICER\left(A_1\right)=\frac{\mathrm{17,705.2659}}{198.69}=89.11,$$

$$ICER\left(A_2\right)=\frac{\mathrm{18,561.1209}-\mathrm{17,705.2659}}{243.15-198.69}=19.25,$$

$$ICER\left(A_3\right)=\frac{8769.2521-\mathrm{18,561.1209}}{353.28-243.15}=-88.9119,$$

$$ICER\left(A_4\right)=\frac{\mathrm{24,893.2485}-8769.2521}{475.189-353.28}=132.2626,$$

Analyzing the results above, one can see in

Table 2. Comparison of strategic cost table.

Strategies	Number of Infections Avoided	Total Costs	ICER
Uncontrolled	0	0	0
$A_1$	198.69	17,705.2659	89.11
$A_2$	243.15	18,561.1209	19.25
$A_3$	353.28	8769.2521	−88.9119
$A_4$	475.189	24,893.2485	132.2626

that strategy $A_2$ is 19.25 cheaper than strategy $A_1$, which means that strategy $A_1$ is more expensive and less efficient than strategy $A_2$, so strategy $A_1$ is removed.

In

Table 3. Revised comparison of strategic cost table.

Strategies	Number of Infections Avoided	Total Costs	ICER
$A_2$	243.15	18,561.1209	19.25
$A_3$	353.28	8769.2521	−88.9119
$A_4$	475.189	24,893.2485	132.2626

, one can clearly see the ICER values of $A_2$ and $A_4$. The comparison shows that $A_4$ is less efficient, so strategy $A_4$ will be removed.

In

Table 4. Second revised comparison of strategic cost table.

Strategies	Number of Infections Avoided	Total Costs	ICER
$A_3$	353.28	8769.2521	−88.9119
$A_4$	475.189	24,893.2485	132.2626

, one can see that strategy $A_4$ is 132.2626 less efficient than strategy $A_3$. So, one concludes on the basis of CEA that strategy $A_3$, which is the combination of treatments and awareness programs, is the most efficient cholera control strategy of all plans.

6. Conclusions

In this paper, a fractional-order optimal control model of cholera is studied, the dynamic behavior of the outbreak is given, and a discriminant of the system stability is presented. We successfully applied the control of integer order systems to fractional dynamical systems. The total population was classified into four groups: susceptible, exposed, infected, and recovered, and there was also the classification of Vibrio cholerae into Vibrio cholerae in the intestine and the environment. The disease-free equilibrium and endemic equilibrium were determined, and sensitivity analyses were performed on the basic reproduction number. The results showed that the disease-free equilibrium point of the model is locally asymptotic stable when $R_0<1$, at which point the disease is indicated to be completely gone. The number of infections per unit time is less than 1; the endemic equilibrium point is locally asymptotic stable when $R_0>1$, when the disease is stable. The number of infections no longer continues to grow. A comparison of disease outcomes using different combinations of three control parameters, treatment, vaccination, and awareness programs, was also studied. Ultimately, it was concluded that the most effective way to control cholera was to use a combination of treatment and awareness programs ($A_3$).

It is noteworthy that there are some future directions for applying the method [32,33,34] to more complex models such as the dynamic behavior of switched memristive neural networks.