Analysis of Atangana–Baleanu fractional-order SEAIR epidemic model with optimal control

We consider a SEAIR epidemic model with Atangana–Baleanu fractional-order derivative. We approximate the solution of the model using the numerical scheme developed by Toufic and Atangana. The numerical simulation corresponding to several fractional orders shows that, as the fractional order reduces from 1, the spread of the endemic grows slower. Optimal control analysis and simulation show that the control strategy designed is operative in reducing the number of cases in different compartments. Moreover, simulating the optimal profile revealed that reducing the fractional-order from 1 leads to the need for quick starting of the application of the designed control strategy at the maximum possible level and maintaining it for the majority of the period of the pandemic.


Introduction
Epidemiological mathematical models provide several aspects for understanding the dynamics of the spread of an epidemic and suggestions of effective control strategies. The insight that the transmission dynamics of endemic diseases can be formulated using mathematical language dates back to 1766 when Daniel Bernoulli published a paper where he described the effects of smallpox variolation on life expectancy [1]. Other early mathematical models in epidemiology were introduced in 1927 when Kermack and McKendrick published a series of papers that described the dynamics of disease transmission in terms of a system of differential equations [2].
Since the beginning of the widespread use of mathematical models for public health making, a large number of studies and publications have been made on modeling and analysis of epidemiological diseases. However, the majority of the studies were restricted to integer-order differential equations. For instance, the global stability of the SEIR and SEIAHR epidemic models with integer derivatives and different saturating contact rates were investigated in [3][4][5][6]. A delayed SIRS epidemic model with integer-order derivative involving saturation incidence and temporary immunity is studied in [7]. An epidemic stochastic mathematical model with integer order is used to predict the spread of coronavirus in [8]. A mathematical model is used to study correlation between the weather conditions and the COVID-19 pandemic in India by Borah et al. [9].
It has recently been turned out that fractional differential equations can successfully be used to model several phenomena in different fields including epidemiology [1]. Fractional calculus is a branch of mathematical analysis that studies calculus of derivatives and integrals of arbitrary orders.
The advantage of describing mathematical models using fractional derivatives is their nonlocal property in the sense that the nth derivative of a function g(x) at x is a local property only when n is an integer, whereas a noninteger fractional derivative of g(x) at x = b depends on all values of g(x) including those far away from b. In other words, in the case of fractional-order derivative, since it involves derivatives and integrals of arbitrary real or complex orders, the future state depends on the present and previous states. Using fractional-order derivatives in modeling, a dynamic system helps to describe the hereditary properties and efficacy (effectiveness, usefulness) of the memory as essential features in many biological mechanisms [10].
Many authors contributed to the development of fractional calculus starting from 1695 when L'Hospital asked Leibniz, what if the order of a derivative is n = 1/2? Some of the other contributors include Euler in 1730 and Lagrange in 1849. In 1812, Laplace defined a fractional derivative using an integral, and in 1819 the derivative of arbitrary order appears in a text by Lacroix. Other prominent names in fractional derivatives and integrals include Reimann-Liouville, Hadamard, Caputo, Caputo-Fabrizio, and more recently Atangana-Baleanu. More historical notes and the nature of fractional calculus can be obtained in [11,12] and the references therein.
In his recent publication, Atangana [13] exposed several previous mistakes and made corrections to the use of fractional calculus related to the fundamental theorem of calculus. Baleanu et al. [14] developed a fractional model for a tumor-immune surveillance mechanism and investigated the effect of chemotherapy treatment on the model. The result showed that the optimal control strategy was efficient. Sene [15] considered a fractional diffusion equation in the context of the fractional operator with Rabotnov fractional exponential kernel and determined the form of the analytical solution of the equation.
In [16] the Black-Scholes equation with Caputo-Fabrizio fractional derivative and a Mittag-Leffler fractional derivative is used to determine the value of an option, which plays an important role in finance. Analysis of a four-dimensional hyperchaotic system with Caputo-Liouville fractional derivative addressing the chaotic, hyperchaotic, and periodic behaviors of the system is considered in [17]. Kahan et al. [18] investigated a COVID-19 mathematical model with a fractal-fractional model in the sense of Atangana-Baleanu fractional operator. The issue of image processing with an Atangana-Baleanu fractional derivative in the sense of Caputo is pondered in [19]. New understandings in the existence of solution for Atangana-Baleanu Willis Aneurysm system and singular perturbation of boundary value problems for the nonlinear fuzzy differential equation are discussed by Panda et al. [20]. An insight on the existence and uniqueness of the solution to a COVID-19 mathematical model using fractional and fractional-fractal operators and fixed point theorem is performed in [21]. In [22,23] a fractional-order mathematical model for COVID-19 is developed, and an investigation of the dynamics of the pandemic is performed. In [24] a fractal-fractional differentiation mathematical model is used for the analysis of diarrhea that occurred in Ghana during 2008-2018. Kahan et al. [10] analyzed HIV-TB coinfected mathematical model in the sense of Atangana-Baleanu fractional derivative.
Numerical computations of different ordinary and fractional derivatives were considered in [25][26][27][28][29]. Kolade and Owolabi [30] performed analysis and numerical simulation of a system using a two-step family of Adams-Bashforth method to approximate the Atangana-Baleanu fractional derivative. A SEIR fractional model with its stability analysis is considered in [31]. The Atangana-Baleanu fractional derivative operator involving the Mittag-Leffler kernel is used to analyze SEIRA mathematical model in [32].
From the above-surveyed works of the literature we can say that fractional derivatives have many applications in mathematical modeling and analysis of real phenomena. In particular, the recently developed Atangana-Baleanu fraction operator has earned popularity and respect due to its immense applications in biological, physical, medical engineering, and several other nonlinear analyses.
Motivated by the aforementioned arguments, in this paper, we study a SEAIR (Susceptible-Exposed-Asymptomatic-Symptomatic-Recovered cases) mathematical model involving a saturating contact rate. The recently developed Atangan-Baleanu fractional derivative and Toufic-Atangana numerical scheme [27] are used to develop the fractional derivative version of the model and estimate its numerical solution. To the best of the authors' knowledge, the SEAIR endemic model applying the Atangana-Baleanu fractional derivative is not yet investigated. The authors also argue that optimal control analysis of mathematical models in the sense of Atangana-Baleanu fractional operators is uncommon in the existing literature. As a result, in this study, we consider an optimal control analysis of the SEAIR model. The rest of this paper is organized as follows. In Sect. 2, we accomplish the description and formulation of the model. In Sect. 3, we establish the existence and uniqueness of the solution of the model including the positivity and boundedness in the sense of the Atangana-Baleanu fractional operator. Section 4 deals with the diseases-free and endemic equilibrium points and the corresponding global stability analysis. The numerical solution of the SEAIR model via Atangana-Baleanu numerical scheme and simulation is detailed Sect. 5. In the last section, we consider an optimal control analysis of the fractional model by incorporating a control parameter in the model. Moreover, in this section, we perform a numerical simulation verifying the effect of the designed control strategy for different values of fractional order and different compartments of the model.

Model description and formulation
In this section, we develop the Atangana-Baleanu fractional derivative representation of the SEAIR endemic mathematical model. Let us first recall the basic definitions of Atangana-Baleanu fractional operators.

Lemma 1 ([34])
The AB fractional derivative and AB fractional integral of the function g ∈ C 1 (a, b) satisfies the Newton-Leibniz equality . 33,35]) For two functions f , g ∈ 1 (a, b), b > a, the AB fractional derivative satisfies the following inequality: Now we proceed with the formulation of the model. The following flow-diagram ( Fig. 1) is used in constructing the mathematical model of this study.
Based on the flow diagram, the mathematical model with integer order used in this study is expressed by the equation system [36] is the saturating contact rate, α and b are positive constants, and the total population is given by  Disease induces death rate α 5 The recovery rate of asymptomatic cases The recovery rate of symptomatic cases α 7 The natural death rate is proportional to the population size N ; the death rate term is α 3 N . Thus in the absence of disease the differential equation of the total population size N is dN/dt =α 3 N , and thus lim x→∞ N(t) = /α 3 , which implies that the carrying capacity of the demographic structure under consideration in this study is /α 3 . Moreover, we can see that C(N) ≈ αN for small N and C(N) ≈ α/b for large N ; C(N) is nondecreasing, and C(N)/N is nonincreasing. We also consider the disease-induced death rate. We assume that once a patient is recovered, he/she develops a permanent immune, and there is no chance of returning to the susceptible group.
Thus the mathematical model taking into account the assumptions, the saturating contact rate, the flow diagram ( Fig. 1), and the AB derivative is described by the system of differential equations where the kernels are given by 1+bN, and α 1 is the probability per unit time of transmitting the infection between two individuals taking part in contact.
In the presence of endemic, we have dN/dt =α 3 Nα 5 (A + I), which indicates that the population size is not constant. The parameters used in the model are indicated in Table 1.

Existence and uniqueness of solutions
In this section, existence and uniqueness, nonnegativity, and boundedness of the solutions of the fractional-order model (3b) is deliberated. To show the existence of the solution to model (3b), we use the famous theorem referred to as the Banach fixed point theorem. For an extensive study on fixed points and contractions, we refer the reader to [21] and the references therein.
To show the existence and uniqueness of the solution, we proceed as follows. Applying the AB fractional integral to model (3b), we obtain
As a result, we have Analogously, the remaining expressions of (7) can be reduced to the following expressions: Theorem 2 The AB fractional model given in (3b) has a solution if we can find M 0 satisfying the inequality Proof From (8) and (9) we have

The existence of the solution (the existence of a fixed point) is confirmed by Theorem 1, and we have to show that the functions S(t), E(t), A(t), I(t), R(t) are solutions of model (3b).
Let us assume that the following are satisfied: From (12) we obtain Repeating the process recursively leads to Applying the limit to both sides of (13) as n → ∞, we see that a 1n (t) → 0 for Similarly, we can show that a 2n (t) → 0, a 3n (t) → 0, a 4n (t) → 0, a 5n (t) → 0, Theorems 1 and 2 guarantee the existence of the solution of model (3b) by the Banach fixed point theorem. The uniqueness of the solution is proved in Theorem 3.

Theorem 4 The epidemiologically feasible region of AB model (3b) is given by
The existence and uniqueness of the solution of model (3b) are now proved, and it remains to show that the set defined in (15) is positively invariant. The following lemma will be used for the proof of Theorem 4. To show that the set is positively invariant, using Lemma 3, we have ABC 0 It follows from (16) that each of the solutions of (3b) is nonnegative and remains in R 5 + , and hence the set defined in (15) is positively invariant for model (3b).
Finally, to establish the boundedness of the solutions of the fractional model (3b), taking into account that all the parameters are positive, we continue by summing all the equations of the model, which gives ABC 0 Applying the Laplace transform leads to Following the work [38] and applying the inverse Laplace transform, the solution is given by where E α,β refers to the Mittag-Leffler function. Taking into account the fact that the Mittag-Leffler function has the asymptotic behavior it is not difficult to observe that N(t) → /α 3 as t → ∞. Hence (15) is the biologically feasible region of model (3b).

I. The disease-free equilibrium point (DFE)
The disease-free equilibrium of (3b) is given by N 0 = ( /α 3 , 0, 0, 0, 0) ∈ ∂ . The global stability of the disease-free equilibrium point is proved in Theorem 5 after defining the basic reproductive number.

II. The basic reproductive number R 0
The basic reproductive number is obtained using the next-generation matrix [39] and is the spectral radius (-TV -1 ), where ,
Proof Consider a Lyapunov function candidate The derivative of V in the direction of the solution of (3b) is given as Since in the set we have N = /α 3 at the DFE and C (N) = αN/(1 + bN) is nondecreasing, dV dt ≤ 0 for R 0A ≤ 1, R 0I ≤ 1 and R 0 = R 0A + R 0I ≤ 1. Moreover, the Lyapunov-Lasalle theorem [40] implies that all trajectories in approach the largest positively invariant subset of the set M where dV /dt = 0. In this work, M is the set where I = A = 0. On the boundary of , where I = A = 0, we have E = 0, R(t) = c 0 e -t → 0 as t → +∞, and Hence all the trajectories in the domain approach the disease-free equilibrium point N 0 = ( /α 3 , 0, 0, 0, 0), for R 0 ≤ 1.

III. Existence and uniqueness of endemic equilibrium point (EEP)
The proof of the global stability of the DFE ensures that when R 0 ≤ 1, there is no other equilibrium point other than N 0 = ( /α 3 , 0, 0, 0, 0). As a result, the study of the endemic equilibrium point N * = (S * , E * , A * , I * , R * ) is restricted to R 0 > 1.

Theorem 6
For R 0 > 1, model (3b) has a unique equilibrium point N * = (E * , A * , I * , R * , S * ) given by (17). The global stability of the EEP is proved in Theorem 7 using the Lyapunov function method.

Theorem 7
If R 0 > 1, then the endemic equilibrium point N * of model (3b) is globally asymptotically stable in the region .
Proof Define a Lyapunov function candidate by Then F(S, E, I, A, R) ≥ 0 and F(S * , E * , we have Note that at the EEP we have N ≤ α 3 . Hence, it follows that dF dt ≤ 0 and dF dt = 0 if and only if S = S * , E = E * , I = I * , A = A * , R = R * . Therefore the largest closed and bounded invariant set in {S, E, I, A, R ∈ :Ḟ = 0} is the set {N * : N * = (S * , E * , I * , A * , R * )}. By LaSalle's invariance principle the unique equilibrium point N * is globally asymptotically stable when R 0 > 1 in the region .

Numerical simulation I
The purpose of this section is to investigate the impact of different values of fractional order η in model (3b). To this end, we give several numerical simulations of the model using a numerical technique developed by Toufic and Atangana as shown in equations (26)- (30). We took some approximations to the real values of the parameters given as   Figure 2 shows that decreasing the fractional order of derivative η from 1 leads to flattening the curves of the suspected cases (S) with a few decreases in the number of cases in the compartment. As the value η gets smaller and smaller, the number of cases approaches the constant value S 0 .
Moreover, from Figs. 3-6 we can say that decreasing the fractional order of derivative η leads to a decrease in the number of exposed, asymptotically infected, symptomatically infected, and recovered cases significantly. In other words, the curves for each of the compartments E, A, I, and R get flattened as η reduces from 1 to 0.6 as shown in the figures. We can conclude that as η gets close to zero from the right, the number of cases in E, A, I, and R gets close to zero. It must also be noted that the EEP for fractional order and the integer order are the same and the solution of the fractional order tends to the EEP over a longer time. From Figs. 2-6 it seems that as the fractional derivative gets smaller and smaller, the time it requires to approach to the EEP gets longer and longer. That is, when the derivative order is reduced from 1, the memory effect of the dynamic system increases, and, consequently, the infection in each of the compartments increases slowly for a long time. Hence the fractional derivative order η affects the dynamics of infection of model (3b) at least in this work. The peak period of the compartments E, A, R, and I is prolonged, and the number of cases at the peak is reduced as η reduces from 1 to zero (see Figs. 3-6).  the result obtained in this study for the SEAIR model by using the Atangana-Baleanu fractional operator with results that could be obtained by applying the other fractional operators for the same model can be a gap for future work.