A Well-Posed Fractional Order Cholera Model with Saturated Incidence Rate

Abstract

A fractional-order cholera model in the Caputo sense is constructed. The model is an extension of the Susceptible–Infected–Recovered (SIR) epidemic model. The transmission dynamics of the disease are studied by incorporating the saturated incidence rate into the model. This is particularly important since assuming that the increase in incidence for a large number of infected individualsis equivalent to a small number of infected individualsdoes not make much sense. The positivity, boundedness, existence, and uniqueness of the solution of the model are also studied. Equilibrium solutions are computed, and their stability analyses are shown to depend on a threshold quantity, the basic reproduction ratio ($R_0$). It is clearly shown that if $R_0<1,$ the disease-free equilibrium is locally asymptotically stable, whereas if $R_0>1$, the endemic equilibrium exists and is locally asymptotically stable. Numerical simulations are carried out to support the analytic results and to show the significance of the fractional order from the biological point of view. Furthermore, the significance of awareness is studied in the numerical section.

1. Introduction

Cholera is a prolific diarrheal disease that leads to death in a short period of time if treatment measures are not taken. The disease is estimated to cause about 21,000 to 143,000 deaths from the 1,300,000 to 4,000,000 cholera cases annually. This represents about 1.62% to 3.58% of the reported cases [1]. People who live in slums and refugee camps due to natural disasters, social conflicts, climate change, and economic meltdowns are the most affected. The camps often possess poor drinking water, which serves as a cholera-causing factor [2]. Many of the symptoms of cholera include vomiting, profuse rice–water stool, sunken eyes, cramps, shock, and severe dehydration. Vibrio cholerae carriers are those people that are exposed to an incomplete cholera-causing dose, therefore, the disease may not manifest any symptoms in their body [3]. Acute cholera leads to death in a short period that varies from hours up to three days. Exposed individuals have only half a chance of being infected with the disease if the concentration of Vibrio cholerae is 105 cells per milliliter. A minimum of 1 liter per day is the least daily consumption of untreated water that may cause the disease [4].

Up to 75% of Vibrio cholerae carriers do not show any symptoms of the disease, potentially spreading it through their feces [5]. This causes a big risk of cholera disease and outbreaks. The cholera infection falls into one of three classes: asymptomatic, mild, or severe. In a population of infected individuals, 80% have mild or moderate symptoms and only 20% develop severe watery diarrhea [6].

Mathematical models help in studying the dynamics of a given infectious disease [7,8,9]. The investigation of cholera models using the SIR epidemic model introduced in 1927 [10] provides a cavernous understanding of the transmission mechanisms of the disease, while other models use time-dependent coefficients [11]. This is why many researchers put more effort into constructing and analyzing mathematical models of a cholera epidemic. Cholera dynamical models involve the study of the disease transmission among humans and V. cholerae concentration in contaminated water. Most of the time, the direct and indirect pathway transmission of cholera gives rise to a basic framework for investigating the dynamics of the disease [12,13,14,15,16,17,18].

Due to their hereditary properties and memory description abilities, many researchers prefer to use fractional order derivatives and integrals as tools in the study of mathematical modeling. Nowadays, a fractional differential equation is used to study biological phenomenaby developing related mathematical models [19,20]. This is due to the fact that fractional calculus can be used to explain the retention and heritage properties of many materials more accurately compared to its corresponding integer-order analog [21,22]. In this paper, we use the Caputo fractional order to model the dynamics of cholera in a homogeneous setting.

The mechanism of transmission of any transmissible disease is controlled by a certain function that depends on the subpopulations of infected individuals called the incidence rate. In epidemic models, the most frequently used incidence rates are the standard incidence rate $\frac{\beta SI}{N}$ and the bilinear incidence rate $\beta SI,$ where β is the contact rate, S stands for the susceptible population, I stands for the infected population, and N is the total population. Assuming the increase in incidence for a large number of infected individualsis equivalent to a small number of infected individuals does not make much sense, hence, there is a need for a more realistic incidence rate of the form $g\left(I\right)S$ [23] where g tends to a saturation level, as it is a nonnegative function, such that $g\left(0\right)=0.$ To incorporate the effect of the behavioral changes of the susceptible individuals, a more general incidence rate of the form $\frac{\beta S{I}^m}{1+\delta I^n}$, where $m$ and $n$ are positive constants and $\delta$ is a nonnegative term, was proposed [24]. These types of incidence rates that allow for the possibility of the introduction of psychological effects are called saturated incidence rates; $\delta$ and $1+\delta I^n$ determine the amount of psychological effect and the inhibition effect, respectively. As the number of infective individuals increases, the rate of infection spread may decrease due to public awareness, potentially leading to contact reduction [25].

Many mathematical models of cholera transmission exist in the literature, and they study the dynamics of the disease. For example, Leo [26] developed an ML reference cholera model that can be used to overcome the existing complexities of modeling the disease. His results indicate, at an average of 87%, that the developed measures can integrate a large number of datasets, including environmental factors, for the timely prediction of cholera epidemics in Tanzania. Daudel et al. [27] constructed and studied a compartmental malaria model. Their results show that the higher implementation of strategies combining awareness programs and therapeutic treatments increases the efficacy of control measures. In [28], a stochastic norovirus epidemic model with a time delay and random perturbations was explored. In [29], a mathematical model for cholera considering vaccination effects was proposed. In [30], Capasso and Paveri-Fontana suggested a mathematical model for the 1973 cholera epidemic in the European Mediterranean region. In 2017, the transmission dynamics of cholera in Yemen were investigated by Nishiura et al. [31]. Lastly, a model containing optimal intervention strategies for cholera control was formulated in [32].

Luchko and Yamamoto [33] proposed a new differential operator with a general kernel function. Due to the existence of flexibility in choosing the kernel, they provide a basis for a broad range of applications [34]. Changing the kernel in the general derivative leads to the discovery of various asymptotic behaviors. Hence, the hidden features of real-world systems are more accurately observed than in the classical sense. However, both properties and applications regarding this operator must be studied in practical situations. There is also the need to state and prove some theorems in order to study models using this operator.

For cholera, the saturated incidence rate is more reasonable than the bilinear incidence rate. This is because it includes the behavioral change and crowding effect of the cholera-infective individuals and prevents the unboundedness of the contact rate by choosing suitable parameters. Motivated by the above discussions, we construct a novel mathematical model that studies the transmission dynamics of cholera in the Caputosense. The model is novel because, to the best of our knowledge, no previous study has analyzed a mathematical model with a saturated incidence rate and fractional derivative for the cholera disease in detail. The main contributions of this research and the new achievements obtained within this manuscript are summarized as follows:

• This paper addresses a new mathematical model of cholera disease which involves the Caputo fractional derivative.

• The fundamental characteristics of the new model are discussed in detail.

• A numerical scheme is developed to carry out numerical simulations.

• The effect of awareness is studied.

• Comparative results in this research show an obvious linkage between the mathematical and biological mechanisms.

Hence, in this paper, we study a fractional-order cholera model with a saturated incidence rate. The main contribution of the paper is the introduction of a more reasonable incidence rate of the form $\frac{\beta S{I}^m}{1+\delta I^n}$, which makes more sense when considering that assuming the increase in incidence for a large number of infected individualsis equivalent to a small number of infected individualsdoes not make much sense. It is also the aim of the paper to consider the effect of awareness in the control of cholera and study all the properties of well-posedness.

The paper is arranged as follows: Section 1 gives an introduction, Section 2 gives important definitions and preliminaries, Section 3 gives model formulation, Section 4 gives the well-posedness properties of the model, and Section 5 gives the numerical simulation results to support the analytic results. In Section 6, a discussion and conclusion are given.

2. Preliminary Definitions and Theorems

Definition 1

[35]. The Caputo fractional derivative of order $\alpha \in \left(n-1, n\right]$ of $f\left(x\right)$ is defined as,

${}_a{}^{C}D{}_x^{\alpha }f\left(x\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{{\displaystyle \int }}_a^x{\left(x-t\right)}^{n-\alpha -1}{f}^n\left(t\right)dt, n=\left[\alpha \right]+1.$

Definition 2

[36]. The linearity of the fractional derivative.

Let $f, g$be continuous and $b, c$ be scalars, then

${}_a{}^{RL}D{}_x^{\alpha }\left[bf\left(x\right)+dg\left(x\right)\right]=b{}_a{}^{RL}D{}_x^{\alpha }f\left(x\right)+d{}_a{}^{RL}D{}_x^{\alpha }g\left(x\right),$

${}_a{}^{C}D{}_x^{\alpha }\left[bf\left(x\right)+dg\left(x\right)\right]=b{}_a{}^{C}D{}_x^{\alpha }f\left(x\right)+d{}_a{}^{C}D{}_x^{\alpha }g\left(x\right).$

Definition 3

[37]. Contraction.

An operator $f:X \to X$ that maps a metric space onto itself is said to be contractive if for $0<q<1$.

$d\left(f\left(x\right), f\left(y\right)\right)=qd\left(x, y\right),\hspace{1em}\forall x, y\in X.$

Theorem 1

[37]. Picard–Banach fixed point.

Any contractive operator that maps a metric space onto itself has a unique fixed point. Furthermore, if $f:X \to X$ is a contractive operator that maps a metric space onto itself and $a$ is its fixed point: $f\left(a\right)=a;$ then for any iterative sequence

$x_0,\hspace{1em}{x}_1=f\left(x_0\right),\hspace{1em}{x}_2=f\left(x_1\right), \dots , x_{n+1}=f\left(x_n\right), \dots$

converges to $a$.

In other words, $a$is a solution or equilibrium for a continuous dynamical system and a fixed point for a discrete dynamical system.

Theorem 2

[36]. The equilibrium solutions $x^{\ast }$ of a system of fractional order differential equation is locally asymptotically stable if all the eigenvalues ${\lambda }_i$ of the Jacobian matrix $\frac{\partial f}{\partial x_i}$ evaluated at the equilibrium points satisfy

$\left|\arg \left({\lambda }_i\right)\right|>\frac{\alpha \pi }{2},\hspace{1em}0<\alpha <1.$

Theorem 3

[36]. Let $x\left(t\right)\in {ℝ}^{+}$ be a continuous and derivable function. Then, for any time instant $t\ge t_0$ and $\alpha \in \left(0, 1\right)$

${}_0{}^{C}D{}_t^{\alpha }[x\left(t\right)-x^{\ast }-x^{\ast }\ln \left(\frac{x\left(t\right)}{x^{\ast }}\right)]\le \left(1-\frac{x\left(t\right)}{x^{\ast }}\right){}_0{}^{C}D{}_t^{\alpha }x\left(t\right),\hspace{1em}{x}^{\ast }\in {ℝ}^{+}.$

3. Model Formulation

The model consists of the following classes: susceptible S, latently infected individuals $E$, infectious individuals $I,$and those that recovered from the infection $R$. We also define $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$. Assuming homogeneous mixing of the population, the model is given as:

(1)$$\begin{gathered} {}_0{}^{C}D{}_t^{\alpha }S\left(t\right)={\Lambda }^{\alpha }-\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}+{\gamma }^{\alpha }E\left(t\right)+{\xi }^{\alpha }I\left(t\right)-{\mu }^{\alpha }S\left(t\right), \\ {}_0{}^{C}D{}_t^{\alpha }E\left(t\right)=p\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}-\left({\gamma }^{\alpha }+{\eta }^{\alpha }+{\varphi }^{\alpha }+{\mu }^{\alpha }\right)E\left(t\right), \\ {}_0{}^{C}D{}_t^{\alpha }I\left(t\right)=\left(1-p\right)\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}+{\eta }^{\alpha }E\left(t\right)-\left({\xi }^{\alpha }+{\delta }^{\alpha }+q^{\alpha }+{\mu }^{\alpha }\right)I\left(t\right), \\ {}_0{}^{C}D{}_t^{\alpha }R\left(t\right)={\delta }^{\alpha }I\left(t\right)-{\mu }^{\alpha }R\left(t\right), \end{gathered}$$

subject to the initial conditions;

$S\left(0\right)>0, E\left(0\right)\ge 0, I\left(0\right)\ge 0, R\left(0\right)\ge 0.$

where the function $f\left(I\right)$ is given as;

$f\left(I\right)=1+I^2.$

We then modify the fractional operator via an auxiliary parameter ${\rm Y}>0$ to avoid dimensional mismatching [38] to obtain

(2)$$\begin{gathered} {{\rm Y}}^{\alpha -1}{}_0{}^{C}D{}_t^{\alpha }S\left(t\right)={\Lambda }^{\alpha }-\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}+{\gamma }^{\alpha }E\left(t\right)+{\xi }^{\alpha }I\left(t\right)-{\mu }^{\alpha }S\left(t\right), \\ {{\rm Y}}^{\alpha -1}{}_0{}^{C}D{}_t^{\alpha }E\left(t\right)=p\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}-\left({\gamma }^{\alpha }+{\eta }^{\alpha }+{\varphi }^{\alpha }+{\mu }^{\alpha }\right)E\left(t\right), \\ {{\rm Y}}^{\alpha -1}{}_0{}^{C}D{}_t^{\alpha }I\left(t\right)=\left(1-p\right)\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}+{\eta }^{\alpha }E\left(t\right)-\left({\xi }^{\alpha }+{\delta }^{\alpha }+q^{\alpha }+{\mu }^{\alpha }\right)I\left(t\right), \\ {{\rm Y}}^{\alpha -1}{}_0{}^{C}D{}_t^{\alpha }R\left(t\right)={\delta }^{\alpha }I\left(t\right)-{\mu }^{\alpha }R\left(t\right), \end{gathered}$$

subject to the initial conditions;

$S\left(0\right)>0, E\left(0\right)\ge 0, I\left(0\right)\ge 0, R\left(0\right)\ge 0.$

The meaning of the parameters in the model are given in Table 1 as follows.

4. Well-Posednessof the Model

In this section, the mathematical properties of the model are explored. This consists of the positivity and boundedness, the existence and uniqueness, the computation of equilibrium solutions and basic reproduction ratio, and the stability analysis of the solutions.

4.1. Positivity and Boundedness

The positivity of solutions means that the population thrives, while boundedness means that the population growth is restricted naturally due to limited resources.

To show positivity, consider Equation (1),

$$\begin{gathered} {}_0{}^{C}D{}_t^{\alpha }S\left(t\right){|}_{S=0}={\Lambda }^{\alpha }+{\gamma }^{\alpha }E\left(t\right)+{\xi }^{\alpha }I\left(t\right)>0, \\ {}_0{}^{C}D{}_t^{\alpha }E\left(t\right){|}_{E=0}=p\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}\ge 0, \\ {}_0{}^{C}D{}_t^{\alpha }I\left(t\right){|}_{I=0}={\eta }^{\alpha }E\left(t\right)\ge 0, \\ {}_0{}^{C}D{}_t^{\alpha }R\left(t\right){|}_{R=0}={\delta }^{\alpha }I\left(t\right)\ge 0. \end{gathered}$$

Hence, we can observe that the solution of (1) is non-negative.

For the boundedness, observe that,

$N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right).$

Then, by Definition 3, we have

${}_0{}^{C}D{}_t^{\alpha }N\left(t\right)={}_0{}^{C}D{}_t^{\alpha }S\left(t\right)+{}_0{}^{C}D{}_t^{\alpha }E\left(t\right)+{}_0{}^{C}D{}_t^{\alpha }I\left(t\right)+{}_0{}^{C}D{}_t^{\alpha }R\left(t\right).$

Then,

$$\begin{gathered} {}_0{}^{C}D{}_t^{\alpha }N\left(t\right)={\Lambda }^{\alpha }-{\mu }^{\alpha }N\left(t\right)-\left({\varphi }^{\alpha }E\left(t\right)+q^{\alpha }I\left(t\right)\right), \\ {}_0{}^{C}D{}_t^{\alpha }N\left(t\right)\le {\Lambda }^{\alpha }-{\mu }^{\alpha }N\left(t\right). \end{gathered}$$

We apply the Laplace transform method to solve the Gronwall’s like inequality with initial condition $N\left(t_0\right)\ge 0$ to find

$L\left\{{}_0{}^{CF}D{}_t^{\alpha }N\left(t\right)+{\mu }^{\alpha }N\right\}\le L\left\{{\Lambda }^{\alpha }\right\}.$

By the linearity of the Laplace transform, we obtain

$L\left\{{}_0{}^{CF}D{}_t^{\alpha }N\left(t\right)\}+{\mu }^{\alpha }L\{N\left(t\right)\right\}\le L\left\{{\Lambda }^{\alpha }\right\},$

we find,

$S^{\alpha }L\left\{N\left(t\right)\right\}-{\displaystyle \sum }_{k=0}^{n-1}{S}^{\alpha -k-1}{N}^k\left(t_0\right)+{\mu }^{\alpha }L\left\{N\left(t\right)\right\}\le \frac{{\Lambda }^{\alpha }}{S}.$

Simplifying, we obtain

$L\left\{N\left(t\right)\right\}\le {\Lambda }^{\alpha }\left(\frac{1}{S}-\frac{1}{S}\frac{1}{\left(1+\frac{{\mu }^{\alpha }}{S^{\alpha }}\right)}\right)+{\displaystyle \sum }_{k=0}^{n-1}\frac{1}{S^{k+1}}\frac{1}{(1+\frac{{\mu }^{\alpha }}{S^{\alpha }})}{N}^k\left(t_0\right).$

Using the Taylor series expansion, we have

$\frac{1}{\left(1+\frac{{\mu }^{\alpha }}{S^{\alpha }}\right)}={\displaystyle \sum }_{n=0}^{\infty }{\left(\frac{-{\mu }^{\alpha }}{S^{\alpha }}\right)}^n.$

Therefore,

$L\left\{N\left(t\right)\right\}\le {\Lambda }^{\alpha }\left(\frac{1}{S}-\frac{1}{S}{\displaystyle \sum }_{n=0}^{\infty }{\left(\frac{-{\mu }^{\alpha }}{S^{\alpha }}\right)}^n\right)+{\displaystyle \sum }_{k=0}^{n-1}\frac{1}{S^{k+1}}{N}^k\left(t_0\right){\displaystyle \sum }_{n=0}^{\infty }{\left(\frac{-{\mu }^{\alpha }}{S^{\alpha }}\right)}^n.$

Taking the Laplace inverse, we find

$N\left(t\right)\le {{\rm Y}}^{\alpha }-{\Lambda }^{\alpha }{\displaystyle \sum }_{n=0}^{\infty }\frac{-{\left({\mu }^{\alpha }{t}^{\alpha }\right)}^n}{\Gamma \left(\alpha n+1\right)}+{\displaystyle \sum }_{k=0}^{n-1}{\displaystyle \sum }_{n=0}^{\infty }\frac{-{\left({\mu }^{\alpha }{t}^{\alpha }\right)}^n}{\Gamma \left(\alpha n+k+1\right)}{t}^k{N}^k\left(t_0\right).$

Substituting the Mittag–Leffler function, we have

$N\left(t\right)\le {\Lambda }^{\alpha }\left[1-E_1\left(-{\mu }^{\alpha }{t}^{\alpha }\right)\right]+{\displaystyle \sum }_{k=0}^{n-1}{E}_{k+1}\left(-{\mu }^{\alpha }{t}^{\alpha }\right)t^k{N}^k\left(t_0\right).$

where $E_1\left(-{\mu }^{\alpha }{t}^{\alpha }\right), E_{k+1}\left(-{\mu }^{\alpha }{t}^{\alpha }\right)$ are the series of the Mittag–Leffler function which converges for any argument, hence, we say that the solution to the model is bounded.

Thus, we define,

$\omega =\left\{\left(S\left(t\right), E\left(t\right),I\left(t\right),R\left(t\right)\right)\in R_{+}^4:S\left(t\right), E\left(t\right),I\left(t\right),R\left(t\right)\le {\Lambda }^{\alpha }\left[1-E_1\left(-{\mu }^{\alpha }{t}^{\alpha }\right)\right]+{\displaystyle \sum }_{k=0}^{n-1}{E}_{k+1}\left(-{\mu }^{\alpha }{t}^{\alpha }\right)t^k{N}^k\left(t_0\right)\right\}$

hence, all solutions of (1) commencing in $\omega$ stays in $\omega$ for all $t\ge 0.$

4.2. Existence and Uniqueness

In this section, we study the existence and uniqueness properties of the solution of Equation (1). First, we consider the following theorem to show Lipschitz continuity.

Theorem 4.

The kernels of Equation (1) satisfy the Lipschitz continuity for $L_i\ge 0, i=1,2,\dots , 4$.

Proof. 

Let the kernel be defined as,

(3)$$\begin{gathered} g_1\left(t,S\right)={\Lambda }^{\alpha }-\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}+{\gamma }^{\alpha }E\left(t\right)+{\xi }^{\alpha }I\left(t\right)-{\mu }^{\alpha }S\left(t\right), \\ g_2\left(t,E\right)=p\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}-\left({\gamma }^{\alpha }+{\eta }^{\alpha }+{\varphi }^{\alpha }+{\mu }^{\alpha }\right)E\left(t\right), \\ {g}_3\left(t,I\right)=\left(1-p\right)\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}+{\eta }^{\alpha }E\left(t\right)-\left({\xi }^{\alpha }+{\delta }^{\alpha }+q^{\alpha }+{\mu }^{\alpha }\right)I\left(t\right), \\ {g}_4\left(t,R\right)={\delta }^{\alpha }I\left(t\right)-{\mu }^{\alpha }R\left(t\right). \end{gathered}$$

Now,

$$\begin{gathered} \left|g_1\left(t, S\right)-g_1\left(t, S^{\ast }\right)\right|=\left|\left(\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)I\left(t\right)}{f\left(I\right)}+{\mu }^{\alpha }\right)\left(S-S^{\ast }\right)\right| \\ \le \left(\left|{\mu }^{\alpha }\right|+\left|\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)I\left(t\right)}{f\left(I\right)}\right|\right)\Vert S-S^{\ast }\Vert \\ \le \left(\left|{\mu }^{\alpha }\right|+{\beta }^{\alpha }\left(1-k^{\alpha }\right)\underset{t\in \left[0, h^{\ast }\right]}{\max }\left|\frac{\left(1-k^{\alpha }\right)I\left(t\right)}{f\left(I\right)}\right|\right)\Vert S-S^{\ast }\Vert \\ \le L_1\Vert S-S^{\ast }\Vert , L_1=\left|{\mu }^{\alpha }\right|+{\beta }^{\alpha }\left(1-k^{\alpha }\right)\underset{t\in \left[0, h^{\ast }\right]}{\max }\left|\frac{I\left(t\right)}{f\left(I\right)}\right|. \end{gathered}$$

Hence,

(4)$\left|g_1\left(t, S\right)-g_1\left(t, S^{\ast }\right)\right|\le L_1\Vert S-S^{\ast }\Vert .$

In a similar way, we obtain

(5)$$\begin{gathered} \left|g_2\left(t, E\right)-g_2\left(t, E^{\ast }\right)\right|\le L_2\Vert E-E^{\ast }\Vert , \\ \left|g_3\left(t, I\right)-g_3\left(t, I^{\ast }\right)\right|\le L_3\Vert I-I^{\ast }\Vert , \\ \left|g_4\left(t, R\right)-g_4\left(t, R^{\ast }\right)\right|\le L_4\Vert R-R^{\ast }\Vert . \end{gathered}$$

where

$L_2=\left|{\gamma }^{\alpha }+{\eta }^{\alpha }+{\varphi }^{\alpha }+{\mu }^{\alpha }\right|, L_3=\left|{\xi }^{\alpha }+{\delta }^{\alpha }+q^{\alpha }+{\mu }^{\alpha }\right|+{\beta }^{\alpha }\left(1-k^{\alpha }\right)\left(1-p\right)\underset{t\in \left[0, h^{\ast }\right]}{\max }\left|\frac{I\left(t\right)}{f\left(I\right)}\right|,$

$\mathrm{and} L_4=\left|{\mu }^{\alpha }\right|.$□

The following Lemma converts the system to Volterra integral equations.

Lemma 1.

The continuous system (1) can be transformed to equivalent Volterra integral equations.

Proof. 

Consider

${}_0{}^{C}D{}_t^{\alpha }S\left(t\right)=g_1\left(t, S\left(t\right)\right).$

On integrating fractionally, we find

(6)$$\begin{gathered} {}_0{}^{C}I{}_t^{\alpha }\left[{}_0{}^{C}D{}_t^{\alpha }S\left(t\right)\right]={}_0{}^{C}I{}_t^{\alpha }[g_1\left(t, S\left(t\right)\right)] \\ S\left(t\right)-S\left(0\right)=\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_1\left(\tau , S\left(\tau \right)\right)dt \\ S\left(t\right)=S\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_1\left(\tau , S\left(\tau \right)\right)dt. \end{gathered}$$

Similarly,

(7)$$\begin{gathered} E\left(t\right)=E\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_2\left(\tau , E\left(\tau \right)\right)dt, \\ I\left(t\right)=I\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_3\left(\tau , I\left(\tau \right)\right)dt, \\ R\left(t\right)=R\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_4\left(\tau , R\left(\tau \right)\right)dt. \square \end{gathered}$$

The following theorem provides the existence of the unique solution.

Theorem 5.

Let $0<\alpha <1, I=\left[0,{\mathrm{h}}^{\ast }\right]\subseteq ℝ \mathit{and}\mathrm{J}=\left|\mathrm{S}\left(\mathrm{t}\right)-\mathrm{S}\left(0\right)\right|\le {\mathrm{k}}_1.$

Let ${\mathrm{g}}_1:\mathrm{I}\times \mathrm{J}\to ℝ$be continuous bounded function, that is $\exists !\mathrm{M}>0$ such that $\left|{\mathrm{g}}_1\left(\mathrm{t},\mathrm{S}\right)\right|\le {\mathrm{M}}_1.$

Assume that ${\mathrm{g}}_1$ satisfies the Lipschitz conditions. If $L_1{k}_1<M_1$, then there exist unique $S\in C\left[0, h^{\ast }\right], where h^{\ast }=min[h, \left(\frac{k_1\Gamma \left(\alpha +1\right)}{M_1}){ }^{\frac{1}{\alpha }}\right]$ that holds the equation.

Proof. 

$Let T=\left\{S\in C\left[0, h^{\ast }\right]: \Vert S\left(t\right)-S\left(0\right)\Vert \le k_1\right\}.$

Since $T\subseteq ℝ$ and its closed set, then $T$ is a complete metric space.

Recall that,

(8)$S\left(t\right)=S\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_1\left(\tau , S\left(\tau \right)\right)dt.$

Define operator $F$ in $T$ by,

(9)$F\left(S\left(t\right)\right)=S\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_1\left(\tau , S\left(\tau \right)\right)dt.$

To show that (6) satisfies Theorem 1, we have

$$\begin{gathered} \left|F\left(S\left(t\right)\right)-S\left(0\right)\right|=\left|\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{g}_1\left(\tau , S\left(\tau \right)\right)dt\right| \\ \le \frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}\left|g_1\left(\tau , S\left(\tau \right)\right)\right|dt \\ \le \frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{M}_{1}dt \\ =\frac{M_1}{\Gamma \left(\alpha +1\right)}{t}^{\alpha } \\ =\frac{M_1}{\Gamma \left(\alpha +1\right)}{\left(h^{\ast }\right)}^{\alpha } \\ \le \frac{M_1}{\Gamma \left(\alpha +1\right)}\frac{k_1\Gamma \left(\alpha +1\right)}{M_1} \\ =k_1. \end{gathered}$$

Hence,

(10)$\left|F\left(S\left(t\right)\right)-S\left(0\right)\right|\le k_1.$

Similarly,

(11)$$\begin{gathered} \left|F\left(E\left(t\right)\right)-E\left(0\right)\right|\le k_2, \\ \left|F\left(I\left(t\right)\right)-I\left(0\right)\right|\le k_3, \\ \left|F\left(R\left(t\right)\right)-R\left(0\right)\right|\le k_4. \end{gathered}$$

Therefore $F$ maps $T$ onto itself.

Secondly, to show that $T$ is a contraction, we have

$F\left(S\right)-F\left(S^{\ast }\right)=S\left(0\right)-S^{\ast }\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}\left[g_1\left(\tau , S\left(\tau \right)\right)-g_1\left(\tau , S^{\ast }\left(\tau \right)\right)\right]d\tau .$

Since $S\left(0\right)=S^{\ast }\left(0\right)$

$$\begin{gathered} \left|F\left(S\right)-F\left(S^{\ast }\right)\right|=\left|\frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}\left[g_1\left(\tau , S\left(\tau \right)\right)-g_1\left(\tau , S^{\ast }\left(\tau \right)\right)\right]d\tau \right| \\ \le \frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}\left|g_1\left(\tau , S\left(\tau \right)\right)-g_1\left(\tau , S^{\ast }\left(\tau \right)\right)\right|d\tau . \\ \le \frac{1}{\Gamma \left(\alpha \right)}{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{L}_1\Vert S-S\Vert d\tau \\ =\frac{L_1}{\Gamma \left(\alpha \right)}\Vert S-S\Vert {{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\alpha -1}{\tau }^{0}d\tau \\ =\frac{L_1}{\Gamma \left(\alpha \right)}\Vert S-S^{\ast }\Vert \frac{\Gamma \left(\alpha \right)}{\Gamma \left(\alpha +1\right)}{t}^{\alpha } \\ =\frac{L_1}{\Gamma \left(\alpha +1\right)}\Vert S-S^{\ast }\Vert t^{\alpha } \\ \le \frac{L_1}{\Gamma \left(\alpha +1\right)}\Vert S-S^{\ast }\Vert {\left(h^{\ast }\right)}^{\alpha } \\ \le \frac{L_1}{\Gamma \left(\alpha +1\right)}\Vert S-S^{\ast }\Vert \frac{k_1\Gamma \left(\alpha +1\right)}{M_1}. \end{gathered}$$

Hence,

(12)$\Vert FS-F{S}^{\ast }\Vert \le \frac{L_1{k}_1}{M_1}\Vert S-S^{\ast }\Vert .$

Since, by hypothesis $\frac{L_1{k}_1}{M_1}<1$, $T$ is contractive and has a unique fixed point.

In a similar way, we obtain

(13)$$\begin{gathered} \Vert F\left(E\right)-F\left(E^{\ast }\right)\Vert \le \frac{L_2{k}_2}{M_2}\Vert E-E^{\ast }\Vert , \\ \Vert F\left(I\right)-F\left(I^{\ast }\right)\Vert \le \frac{L_3{k}_3}{M_3}\Vert I-I^{\ast }\Vert , \\ \Vert F\left(R\right)-F\left(R^{\ast }\right)\Vert \le \frac{L_4{k}_4}{M_4}\Vert R-R^{\ast }\Vert . \end{gathered}$$

Thus, Equation (1) has a unique solution. □

4.3. Existence of Equilibrium Solutions

Since $R=N-\left(S+E+I\right),$we can consider the first three equations in Equation (1) for this analysis.

Setting the first three equations in Equation (1) to zero and solving simultaneously, we find the following equilibrium solutions;

• Disease-free equilibrium $\left(E_0\right)$ is given as;

  $E_0=\left(S^0,0,0\right)=\left(\frac{{\Lambda }^{\alpha }}{{\mu }^{\alpha }},0,0\right).$

• Endemic equilibrium $\left(E_1\right)$ is given as;

  $E_1=\left(S^1,E^1,I^1\right),$

  where

  $$\begin{gathered} S^1=\frac{A_1{A}_{2}f\left(I\right)}{{\beta }^{\alpha }\left(1-k^{\alpha }\right)\left[\left(1-p\right)A_1+{\eta }^{\alpha }p\right]}, \\ {E}^1=\frac{p{A}_2{I}^1}{{\beta }^{\alpha }\left(1-k^{\alpha }\right)\left[\left(1-p\right)A_1+{\eta }^{\alpha }p\right]}, \end{gathered}$$

  and $I^1$ can be obtained by solving,

  ${\Lambda }^{\alpha }-\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)SI}{f\left(I\right)}-{\mu }^{\alpha }S+{\gamma }^{\alpha }E+{\xi }^{\alpha }I=0,$

  where

  $A_1={\gamma }^{\alpha }+{\eta }^{\alpha }+{\varphi }^{\alpha }+{\mu }^{\alpha }, and A_2={\xi }^{\alpha }+{\delta }^{\alpha }+q^{\alpha }+{\mu }^{\alpha }.$

Simplifying, we find

$G\left(I\right)={\Lambda }^{\alpha }-I\left[\frac{A_1{A}_2}{\left(1-p\right)A_1+{\eta }^{\alpha }p}-\left(\frac{{\gamma }^{\alpha }p{A}_2}{\left(1-p\right)A_1+{\eta }^{\alpha }p}+{\xi }^{\alpha }\right)\right]-\frac{A_1{A}_{2}f\left(I\right)}{{\beta }^{\alpha }\left(1-k^{\alpha }\right)\left[\left(1-p\right)A_1+{\eta }^{\alpha }p\right]}=0.$

Since,

$\frac{A_1{A}_2}{\left(1-p\right)A_1+{\eta }^{\alpha }p}-\left(\frac{{\gamma }^{\alpha }p{A}_2}{\left(1-p\right)A_1+{\eta }^{\alpha }p}+{\xi }^{\alpha }\right)>0, and f^{\prime }\ge 0,$

then G is an increasing function.

In addition,

$G\left(I\right)<{\Lambda }^{\alpha }-I\left[\frac{A_1{A}_2}{\left(1-p\right)A_1+{\eta }^{\alpha }p}-\left(\frac{{\gamma }^{\alpha }p{A}_2}{\left(1-p\right)A_1+{\eta }^{\alpha }p}+{\xi }^{\alpha }\right)\right],$

then,

$\underset{t\to \infty }{\lim }G\left(I\right)=-\infty .$

If $f\left(0\right)=1,$ clearly

$$\begin{gathered} G\left(0\right)={\Lambda }^{\alpha }-\frac{A_1{A}_2}{{\beta }^{\alpha }\left(1-k^{\alpha }\right)\left[\left(1-p\right)A_1+{\eta }^{\alpha }p\right]} \\ ={\Lambda }^{\alpha }\left[1-\left(\frac{{\mu }^{\alpha }{A}_1{A}_2}{{\eta }^{\alpha }p{\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}+\frac{{\mu }^{\alpha }{A}_2}{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}\right)\right]. \end{gathered}$$

Thus, a unique positive solution for $G$ exists $if f G\left(0\right)>0, i.e.,$ if

$\frac{{\mu }^{\alpha }{A}_1{A}_2}{{\eta }^{\alpha }p{\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}+\frac{{\mu }^{\alpha }{A}_2}{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}>1.$

We shall establish,

$\frac{{\mu }^{\alpha }{A}_1{A}_2}{{\eta }^{\alpha }p{\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}+\frac{{\mu }^{\alpha }{A}_2}{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}=R_0,$

in the subsequent section where $R_0$ is the basic reproduction ratio.

4.4. Basic Reproduction Ratio

This is defined as the number of secondary cases caused by a single infected individual in a population of Susceptible [37]. Here we use the idea of the next-generation matrix as in [39].

Let,

$Ƴ$ represents new infection and $Ƶ$ represents the remaining terms in Equation (1). Then,

$Ƴ=\left[\begin{array}{c}p\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}\\ \left(1-p\right)\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S\left(t\right)I\left(t\right)}{f\left(I\right)}\end{array}\right], and Ƶ=\left[\begin{array}{c}{A}_{1}E\\ -{\eta }^{\alpha }E+A_{2}I\end{array}\right].$

Then,

$Y\left(E_0\right)=\left[\begin{array}{cc}0& p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\\ 0& \left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\end{array}\right], and Z\left(E_0\right)=\left[\begin{array}{cc}{A}_1& 0\\ -{\eta }^{\alpha }& A_2\end{array}\right],$

where $Y\left(E_0\right)$ is the Jacobian of $Ƴ$ at $E_0$ and $Z\left(E_0\right)$ is the Jacobian at $Ƶ$ at $E_0.$

Therefore,

$X=Y\left(E_0\right)(Z{\left(E_0\right)}^{-1})=\left[\begin{array}{cc}{\eta }^{\alpha }p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0& A_{1}p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\\ {\eta }^{\alpha }\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0& A_1\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\end{array}\right].$

Note: $X$ is the next generation matrix and the dominant eigenvalue of $X$ is the basic reproduction ratio ($R_0$).

Here,

$R_0=\frac{{\mu }^{\alpha }{A}_1{A}_2}{{\eta }^{\alpha }p{\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}+\frac{{\mu }^{\alpha }{A}_2}{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right){\Lambda }^{\alpha }}.$

4.5. Stability Analysis of the Equilibria

In this section, we carry out a local stability analysis of both disease-free and endemic equilibrium points.

First, consider the following Jacobian matrix from Equation (1);

(14)$J=\left[\begin{array}{ccc}-\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)I}{f\left(I\right)}-{\mu }^{\alpha }& {\gamma }^{\alpha }& -\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S}{f\left(I\right)}\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)+{\xi }^{\alpha }\\ \frac{p{\beta }^{\alpha }\left(1-k^{\alpha }\right)I}{f\left(I\right)}& -A_1& \frac{p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S}{f\left(I\right)}\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)\\ \frac{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)I}{f\left(I\right)}& {\eta }^{\alpha }& \frac{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S}{f\left(I\right)}\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)-A_2\end{array}\right].$

Theorem 6.

The disease-free equilibrium ($E_0$) is locally asymptotically stable if $R_0<1$.

Proof. 

Consider Equation (14) at $E_0,$ then we have

$J\left(E_0\right)=\left[\begin{array}{ccc}-{\mu }^{\alpha }& {\gamma }^{\alpha }& {\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0+{\xi }^{\alpha }\\ 0& -A_1& p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\\ 0& {\eta }^{\alpha }& \left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0-A_2\end{array}\right].$

In echelon form, we find

(15)$J\left(E_0\right)=\left[\begin{array}{ccc}-{\mu }^{\alpha }& {\gamma }^{\alpha }& {\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0+{\xi }^{\alpha }\\ 0& -A_1& p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\\ 0& {\eta }^{\alpha }& A_1\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0-A_1{A}_2+{\eta }^{\alpha }p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^0\end{array}\right].$

The eigenvalues of Equation (15) are:

${\lambda }_1=-{\mu }^{\alpha }<0, {\lambda }_2=-A_1, and {\lambda }_3=-A_1{A}_2\left(R_0-1\right)<0 if R_0<1.$

Hence, $E_0$ is locally asymptotically stable if $R_0<1.$ □

Theorem 7.

The endemic equilibrium ($E_1$) is locally asymptotically stable if $R_0>1$.

Proof. 

Consider Equation (14) at $E_1,$ then we find

$\begin{array}{l}J\left(E_1\right)\\ =\left[\begin{array}{ccc}-\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)I^1}{f\left(I\right)}-{\mu }^{\alpha }& {\gamma }^{\alpha }& -\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^1}{f\left(I\right)}\left(1-\frac{I^1{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)+{\xi }^{\alpha }\\ \frac{p{\beta }^{\alpha }\left(1-k^{\alpha }\right)I^1}{f\left(I\right)}& -A_1& \frac{p{\beta }^{\alpha }\left(1-k^{\alpha }\right)S^1}{f\left(I\right)}\left(1-\frac{I^1{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)\\ \frac{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)I^1}{f\left(I\right)}& {\eta }^{\alpha }& \frac{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)S^1}{f\left(I\right)}\left(1-\frac{I^1{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)-A_2\end{array}\right].\end{array}$

After computing the echelon form,

$\begin{array}{cc}\hfill trace\left(J\left(E_1\right)\right)=& -\left({\mu }^{\alpha }{A}_1+p\left({\mu }^{\alpha }+{\varphi }^{\alpha }\right)\right)\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)I}{f\left(I\right)}I\left(\left(1-p\right)A_1+p{\eta }^{\alpha }\right)\hfill \\ & -\frac{\left(1-p\right){\beta }^{\alpha }\left(1-k^{\alpha }\right)I}{f\left(I\right)}\left({\mu }^{\alpha }+q^{\alpha }+{\delta }^{\alpha }\right)I\hfill \\ & -A_2{\mu }^{\alpha }\left[\left(\left(1-p\right)A_1+p{\eta }^{\alpha }\right)-\left(1-p\right)A_1{A}_2\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)\right]\hfill \end{array}$

and

$\begin{array}{cc}\hfill det\left(J\left(E_1\right)\right)=& \left({\mu }^{\alpha }{A}_1+\frac{p{\beta }^{\alpha }\left(1-k^{\alpha }\right)\left({\mu }^{\alpha }+{\varphi }^{\alpha }\right)I}{f\left(I\right)}\right)[\left({\mu }^{\alpha }+q^{\alpha }+{\delta }^{\alpha }\right)\left(1-p\right)\frac{{\beta }^{\alpha }\left(1-k^{\alpha }\right)I}{f\left(I\right)}\hfill \\ & +A_2{\mu }^{\alpha }\left(\left(1-p\right)A_1+p{\eta }^{\alpha }\right)-{\mu }^{\alpha }\left(1-p\right)A_1{A}_2\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)]\hfill \\ & +\left({\mu }^{\alpha }\left(1-p\right)+{\mu }^{\alpha }\right)\left[\frac{{\beta }^{\alpha }{\mu }^{\alpha }{A}_1{A}_2\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)\left(1-k^{\alpha }\right)I^1}{\left(\left(1-p\right)A_1+p{\eta }^{\alpha }\right)+f\left(I\right)p\left({\mu }^{\alpha }+q^{\alpha }+{\delta }^{\alpha }\right)}I\right].\hfill \end{array}$

Clearly, $trace\left(J\left(E_1\right)\right)<0$ and $det\left(J\left(E_1\right)\right)>0$ if $\left(\left(1-p\right)A_1+p{\eta }^{\alpha }\right)-{\mu }^{\alpha }\left(1-p\right)A_1{A}_2\left(1-\frac{I{f}^{\prime }\left(I\right)}{f\left(I\right)}\right)>0.$ This is equivalent to $R_0>1.$ □

5. Numerical Simulations

This section is devoted to testing the performance of the proposed fractional-order model (1) under the Caputo differential operator while using a numerical explicit technique called the Adams–Bashforth–Moulton technique, also known as the fractional predictor-corrector method introduced and analyzed for its convergence and error bounds in [40,41]; $s=min\left(1+\alpha , 2\right)$ is the order of accuracy for the numerical technique. It is also worth mentioning that, unlike newly established numerical techniques for classical initial value problems, the present literature is not rich enough with numerical techniques for fractional order differential equations.

In this research, the parameter values used are from [28]; $\Lambda =6,\beta =10,\mu =0.02,\gamma =0.01,\xi =0.02,\varphi =0.01,q=0.02,\eta =0.2,\delta =0.1,p\in \left[0,1\right],k=0.2,\mathrm{and}\alpha \in \left(0,1\right].$ Consider an equi-spaced mesh $t_i=i\Delta t,i=0,1,\dots ,M,$where M is a positive integer and $\Delta t=\frac{T}{M},$ where $T$ is the upper limit of the closed interval of integration $\left[0,T\right].$ This setting leads to the following structure for the predictor part required by the numerical technique under consideration:

(16)$$\begin{gathered} S_{i+1}^s=S\left(0\right)+{\displaystyle \sum }_{i=1}^n{b}_{\alpha ,i,n+1}{g}_1\left(t_i,S_i,E_i,I_i,R_i\right), \\ {E}_{i+1}^s=E\left(0\right)+{\displaystyle \sum }_{i=1}^n{b}_{\alpha ,i,n+1}{g}_2\left(t_i,S_i,E_i,I_i,R_i\right), \\ I_{i+1}^s=I\left(0\right)+{\displaystyle \sum }_{i=1}^n{b}_{\alpha ,i,n+1}{g}_3\left(t_i,S_i,E_i,I_i,R_i\right), \\ {R}_{i+1}^s=R\left(0\right)+{\displaystyle \sum }_{i=1}^n{b}_{\alpha ,i,n+1}{g}_4\left(t_i,S_i,E_i,I_i,R_i\right). \end{gathered}$$

Similarly, for the corrector, we have

(17)$$\begin{gathered} S_{i+1}^u=S\left(0\right)+a_{\alpha ,n+1,n+1}{g}_1\left(t_i,S_i{}^s,E_i{}^s,I_i{}^s,R_i{}^s\right)+{\displaystyle \sum }_{i=1}^n{a}_{\alpha ,i,n+1}{g}_1\left(t_i,S_i,E_i,I_i,R_i\right), \\ {E}_{i+1}^u=E\left(0\right)+a_{\alpha ,n+1,n+1}{g}_2\left(t_i,S_i{}^s,E_i{}^s,I_i{}^s,R_i{}^s\right)+{\displaystyle \sum }_{i=1}^n{a}_{\alpha ,i,n+1}{g}_2\left(t_i,S_i,E_i,I_i,R_i\right), \\ I_{i+1}^u=I\left(0\right)+a_{\alpha ,n+1,n+1}{g}_3\left(t_i,S_i{}^s,E_i{}^s,I_i{}^s,R_i{}^s\right)+{\displaystyle \sum }_{i=1}^n{a}_{\alpha ,i,n+1}{g}_3\left(t_i,S_i,E_i,I_i,R_i\right), \\ {R}_{i+1}^u=R\left(0\right)+a_{\alpha ,n+1,n+1}{g}_4\left(t_i,S_i{}^s,E_i{}^s,I_i{}^s,R_i{}^s\right)+{\displaystyle \sum }_{i=1}^n{a}_{\alpha ,i,n+1}{g}_4\left(t_i,S_i,E_i,I_i,R_i\right). \end{gathered}$$

where

$a_{\alpha ,i,n+1}=\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +2\right)}, \mathrm{and} b_{\alpha ,i,n+1}=\frac{1}{\Gamma \left(\alpha +2\right)}\left[{\left(n-i+1\right)}^{\alpha }-{\left(n-i\right)}^{\alpha }\right].$

and s denote the predictor and u represent the corrector.

Additionally,

$\{\begin{array}{cc}{n}^{\alpha +1}-\left(n-\alpha \right){\left(n+1\right)}^{\alpha }, \hfill & i=0\hfill \\ {\left(n-i+2\right)}^{\alpha }-2{\left(n-i+1\right)}^{\alpha +1}+{\left(n-1\right)}^{\alpha +1},\hfill & 1\le i\le n\hfill \\ 1, \hfill & i=n+1.\hfill \end{array}$

The behavior of the Susceptible S(t), Exposed E(t), Infected I(t), and Recovered R(t) population can be seen in Figure 1a–d where the required simulations have been carried out for T = 100 days while varying the values of the fractional order parameter $\alpha$. For the Susceptible population, it is observed that the population decays at a faster rate for higher values of the fractional-order parameter. On the contrary, the higher values of fractional-order parameter $\alpha$ lead to an increase in the Exposed, Infected, and Recovered populations.

The behavior of the Exposed E(t) and Infected I(t) classes can be seen in Figure 2a,b where the required simulations have been carried out for T = 100 days while varying the values of the awareness parameter k. For both the Exposed and Infected classes, it can be observed that a higher value of awareness parameter $k$ leads to a decrease in the Exposed and Infected populations.

6. Summary and Conclusions

In this paper, a fractional order cholera model in the Caputo sense is constructed. The transmission dynamics of the disease are studied by incorporating the saturated incidence rate into the model. Various well-posedness properties such as positivity, boundedness, existence, and uniqueness of the solution are also studied. Equilibrium solutions were computed, and their stability analysis was shown to depend on a threshold quantity, the basic reproduction ratio ($R_0$). It was clearly shown that if $R_0<1,$ the disease-free equilibrium is locally asymptotically stable, whereas if $R_0>1$, endemic equilibrium exists and is locally asymptotically stable. The model is found to be well-posed and more realistic than many related models in the literature due to the consideration of the saturated incidence rate.

Numerical simulations were carried out to support the analytic result and to show the significance of the fractional order from the biological point of view as well as show the significance of awareness programs in curtailing the spread of cholera in a population. From the result obtained, it is observed that the Susceptible population decay at a faster rate for higher values of the fractional-order parameter. On the contrary, the higher values of fractional-order parameter $\alpha$ lead to an increase in the Exposed, Infected, and Recovered populations. In addition, it is observed that for both the Exposed and Infected classes, the increase in the value of awareness parameter $k$ leads to a decrease in the Exposed and Infected populations, respectively.

In conclusion, this paper studied a fractional order model in the Caputo sense with a saturated incidence rate. The entire well-posedness properties of the model were studied in detail. The awareness contribution in controlling cholera was studied numerically. From the findings of this research, one can see there is a need for the consideration of a saturated incidence rate in studying epidemic diseases such as cholera. Additionally, there is a need for studying well-posedness properties before making deductions on any mathematical model. From the findings of this model in particular, there is a need for relevant agencies to mount awareness programs in order to curtail the spread of cholera in a given population. Using a fractional-order to construct a mathematical model is also of paramount importance due to its hereditary properties and memory description ability. The limitation of this work is that there is a need to consider the agent-based approach of the model proposed in this paper which will help in providing more decision-making tools for the effective prevention and control of cholera and public health interventions.

Author Contributions

Conceptualization, I.A.B. and U.W.H.; methodology, I.A.B.; software, I.A.B. and F.A.R.; validation, I.A.B., U.W.H. and F.A.R.; formal analysis, I.A.B. and U.W.H.; investigation, F.A.R.; resources, I.A.B.; data curation, F.A.R.; writing—original draft preparation, I.A.B., U.W.H. and F.A.R.; writing—review and editing, I.A.B. and F.A.R.; visualization, F.A.R.; supervision, U.W.H.; project administration, U.W.H.; funding acquisition, U.W.H. and F.A.R. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

This research was supported by King Mongkut’s University of Science and Technology Thonburi’s Postdoctoral Fellowship.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. (a) The dynamics of the Susceptible population for varying fractional order [Forumla omitted. See PDF.]. (b) The dynamics of the Exposed population for varying fractional order [Forumla omitted. See PDF.]. (c) The dynamics of the Infected population for varying fractional order [Forumla omitted. See PDF.] (d) The dynamics of the Recovered population for varying fractional order [Forumla omitted. See PDF.].

Figure 2. (a) The dynamics of the Exposed classes by increasing awareness parameter k. (b) The dynamics of the Infected classes by increasing awareness parameter k.

Figure 2. (a) The dynamics of the Exposed classes by increasing awareness parameter k. (b) The dynamics of the Infected classes by increasing awareness parameter k.

Table 1

Meaning of parameters.

Parameter	Meaning
Λ	Birthrate
β	Disease contact rate
$\mu$	Natural death rate
$\varphi$	Disease-induced death rate in the Exposed class
γ	Rate at which the Exposed become Susceptible
𝜉	Rate at which the Infectious become Susceptible
𝜂	Rate at which the Exposed become Infectious
d	Disease-induced death rate in the Infectious class
$\delta$	Rate of recovery
k	Awareness parameter
$p$	Fraction of individuals joining the Exposed class

References

1. Ganesan, D.; Gupta, S.S.; Legros, D. Cholera surveillance and estimation of burden of cholera. Vaccine; 2020; 38, pp. A13-A17. [DOI: https://dx.doi.org/10.1016/j.vaccine.2019.07.036] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31326254]

2. Idoga, P.E.; Toycan, M.; Zayyad, M.A. Analysis of factors contributing to the spread of cholera in developing countries. Eurasian J. Med.; 2019; 51, pp. 121-127. [DOI: https://dx.doi.org/10.5152/eurasianjmed.2019.18334] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31258350]

3. Llanes, R.; Somarriba, L.; Hern’andez, G.; Bardaj, Y.; Aguila, A.; Mazumder, R.N. Low detection of vibrio cholera carriage in healthcare workers returning to 12 Latin American countries from Haiti. Epidemiol. Infect.; 2015; 143, pp. 1016-1019. [DOI: https://dx.doi.org/10.1017/S0950268814001782] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25016919]

4. Momba, M.; El-Liethy, M.A. Vibrio Cholerae and Cholera Biotypes; Global Water Pathogen Project: Pretoria, South Africa, 2018.

5. Javidi, M.; Ahmad, B. A study of a fractional-order cholera model. Appl. Math. Inf. Sci.; 2014; 8, pp. 2195-2206. [DOI: https://dx.doi.org/10.12785/amis/080513]

6. Lamond, E.; Kinyanjui, J. Cholera Outbreak Guidelines: Preparedness, Prevention and Control; Oxfam GB: Oxford, UK, 2012.

7. Eskandari, Z.; Avazzadeh, Z.; Khoshsiar, G.; Li, B. Dynamics and bifurcations of a discrete—Time Lotka–Volterra model using nonstandard finite difference discretization method. Math. Methods Appl. Sci.; 2022; [DOI: https://dx.doi.org/10.1002/mma.8859]

8. Li, B.; Liang, H.; Shi, L.; He, Q. Complex dynamics of Kopel model with nonsymmetric response between oligopolists. Chaos Solitons Fractals; 2022; 156, 111860. [DOI: https://dx.doi.org/10.1016/j.chaos.2022.111860]

9. Liu, X.; Arfan, M.; Ur Rahman, M.; Fatima, B. Analysis of SIQR type mathematical model under Atangana-Baleanu fractional differential operator. Comput. Methods Biomech. Biomed. Eng.; 2022; 3, pp. 1-5. [DOI: https://dx.doi.org/10.1080/10255842.2022.2047954]

10. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. A Contain. Pap. A Math. Phys. Character; 1927; 115, pp. 700-721.

11. Marinov, T.; Marinova, R. Adaptive SIR model with vaccination: Simultaneous identification of rates and functions illustrated with COVID-19. Sci. Rep.; 2022; 12, pp. 1-13.

12. Alam, A.; LaRocque, R.C.; Harris, J.B.; Vanderspurt, C.; Ryan, E.T.; Qadri, F.; Calderwood, S.B. Hyperinfectivity of human-passaged Vibrio cholerae can be modelled by growth in the infant mouse. Infect. Immun.; 2005; 73, pp. 6674-6679. [DOI: https://dx.doi.org/10.1128/IAI.73.10.6674-6679.2005]

13. Codeço, C. Endemic and epidemic dynamics of cholera: The role of the aquatic reservoir. BMC Infect. Dis.; 2001; 1, pp. 1-14. [DOI: https://dx.doi.org/10.1186/1471-2334-1-1]

14. Hartley, D.M., Jr.; Morris, J.; Smith, D. Hyperinfectivity: A critical element in the ability of V. cholerae to cause epidemics?. PLoS Med.; 2006; 3, pp. 63-69. [DOI: https://dx.doi.org/10.1371/journal.pmed.0030007] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16318414]

15. Mukandavire, Z.; Liao, S.; Wang, J.; Gaff, H.; Smith, D.; Morris, J. Estimating the reproductive numbers for the 2008–2009 cholera outbreaks in Zimbabwe. Proc. Natl. Acad. Sci. USA; 2011; 108, pp. 8767-8772. [DOI: https://dx.doi.org/10.1073/pnas.1019712108] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21518855]

16. Nelson, E.; Harris, J.; Morris, J.; Calderwood, S.; Camilli, A. Cholera transmission: The host, pathogen and bacteriophage dynamics. Nat. Rev. Microbiol.; 2009; 7, pp. 693-702. [DOI: https://dx.doi.org/10.1038/nrmicro2204] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19756008]

17. Shuai, Z.; Tien, J.H.; Van den Driessche, P. Cholera models with hyper infectivity and temporary immunity. Bull. Math. Biol.; 2012; 74, pp. 2423-2445. [DOI: https://dx.doi.org/10.1007/s11538-012-9759-4]

18. Shuai, Z.; van den Driessche, P. Modeling and control of cholera on networks with a common water source. J. Biol. Dyn.; 2015; 9, pp. 90-103. [DOI: https://dx.doi.org/10.1080/17513758.2014.944226]

19. Escalante-Martínez, J.E.; Gómez-Aguilar, J.F.; Calderón-Ramón, C.; Aguilar-Meléndez, A.; Padilla-Longoria, P. Synchronized bioluminescence behavior of a set of fireflies involving fractional operators of Liouville Caputo type. Int. J. Biomath.; 2018; 11, pp. 1-24. [DOI: https://dx.doi.org/10.1142/S1793524518500419]

20. Escalante-Martínez, J.E.; Gómez-Aguilar, J.F.; Calderón-Ramón, C.; Aguilar-Meléndez, A.; Padilla-Longoria, P. A mathematical model of circadian rhythms synchronization using fractional differential equations system of coupled van der Pol oscillators. Int. J. Biomath.; 2018; 11, 1850041. [DOI: https://dx.doi.org/10.1142/S1793524518500419]

21. Ullah, S.; Khan, M.A.; Farooq, M. A fractional model for the dynamics of TB virus. Chaos Solitons Fractals; 2018; 116, pp. 63-71. [DOI: https://dx.doi.org/10.1016/j.chaos.2018.09.001]

22. Aguilar, J.F.G. Fundamental solutions to electrical circuits of non-integer order via fractional derivatives with and without singular kernels. Eur. Phys. J. Plus; 2018; 133, pp. 1-20.

23. Capasso, V.; Serio, G. A generalization of the Kermack-Mckendrick deterministic epidemic model. Math. Biosci.; 1978; 42, pp. 43-61. [DOI: https://dx.doi.org/10.1016/0025-5564(78)90006-8]

24. Liu, W.M.; Levin, S.A.; Iwasa, Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models. J. Math. Biol.; 1986; 23, pp. 187-204. [DOI: https://dx.doi.org/10.1007/BF00276956] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3958634]

25. Xiao, D.; Ruan, S. Global analysis of an epidemic model with non monotone incidence rate. Math. Biosci.; 2007; 208, pp. 419-429. [DOI: https://dx.doi.org/10.1016/j.mbs.2006.09.025] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17303186]

26. Leo, J. Complexity of Epidemics Models: A Case-Study of Cholera in Tanzania. Digital Transformation for Sustainability; Marx Gómez, J.; Lorini, M.R. Springer: Cham, Switzerland, 2022; [DOI: https://dx.doi.org/10.1007/978-3-031-15420-1_18]

27. Tchatat, D.; Kolaye, G.; Bowong, S.; Temgoua, A. Theoretical assessment of the impact of awareness programs on cholera transmission dynamic. Int. J. Nonlinear Sci. Numer. Simul.; 2022; [DOI: https://dx.doi.org/10.1515/ijnsns-2021-0341]

28. Wang, Y.; Abdeljawad, T.; Din, A. Modeling the dynamics of stochastic norovirus epidemic model with time-delay. Fractals; 2022; 30, pp. 1-13. [DOI: https://dx.doi.org/10.1142/S0218348X22401508]

29. Lemos-Paião, A.; Silva, C.; Torres, D. A cholera mathematical model with vaccination and the biggest outbreak of world’s history. AIMS Math.; 2018; 3, pp. 448-463. [DOI: https://dx.doi.org/10.3934/Math.2018.4.448]

30. Capasso, V.; Paveri-Fontana, S. A mathematical model for the 1973 cholera epidemic in the European Mediterranean region. Rev. D’epidémiologie Et De St.; 1979; 27, pp. 121-132.

31. Nishiura, H.; Tsuzuki, S.; Yuan, B.; Yamaguchi, T.; Asai, Y. Transmission dynamics of cholera in Yemen, 2017: A real time forecasting Theoret. Biol. Med. Model.; 2017; 14, 14. [DOI: https://dx.doi.org/10.1186/s12976-017-0061-x]

32. Neilan, R.; Schaefer, E.; Gaff, H.; Fister, K.; Lenhart, S. Modeling optimal intervention strategies for cholera. Bull. Math. Biol.; 2010; 72, pp. 2004-2018.

33. Luchko, Y.; Yamamoto, M. General time-fractional diffusion equation: Some uniqueness and existence results for the initial-boundary-value problems. Fract. Calculus Appl. Anal.; 2016; 19, pp. 676-695. [DOI: https://dx.doi.org/10.1515/fca-2016-0036]

34. Jajarmi, A.; Baleanu, D.; Vahid, K.Z.; Pirouz, H.M.; Asad, J.H. A new and general fractional Lagrangian approach: A capacitor microphone case study. Res. Phys.; 2021; 31, 104950. [DOI: https://dx.doi.org/10.1016/j.rinp.2021.104950]

35. Podlubny, I. Fractional Differential Equations; Academic Press: San Diego, CA, USA, 1999.

36. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Application of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006; Volume 204.

37. van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci.; 2002; 180, pp. 29-48. [DOI: https://dx.doi.org/10.1016/S0025-5564(02)00108-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12387915]

38. Ullah, M.; Baleanu, D. A new fractional SICA model and numerical method for the transmission of HIV/AIDS. Math. Methods Appl. Sci.; 2021; 44, pp. 8648-8659. [DOI: https://dx.doi.org/10.1002/mma.7292]

39. Diethelm, K.; Ford, N.; Freed, A. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn.; 2002; 29, pp. 3-22. [DOI: https://dx.doi.org/10.1023/A:1016592219341]

40. Diethelm, K.; Ford, N.; Freed, A. Detailed error analysis for a fractional Adams method. Numer. Algorithms; 2004; 36, pp. 31-52. [DOI: https://dx.doi.org/10.1023/B:NUMA.0000027736.85078.be]

41. Khan, M.; Parvez, M.; Islam, S.; Khan, I.; Shafie, S.; Gul, T. Mathematical Analysis of Typhoid Model with Saturated Incidence Rate. Adv. Stud. Biol.; 2015; 7, pp. 65-78. [DOI: https://dx.doi.org/10.12988/asb.2015.41059]