Analysis of a Fractional-Order COVID-19 Epidemic Model with Lockdown

The outbreak of the coronavirus disease (COVID-19) has caused a lot of disruptions around the world. In an attempt to control the spread of the disease among the population, several measures such as lockdown, and mask mandates, amongst others, were implemented by many governments in their countries. To understand the effectiveness of these measures in controlling the disease, several mathematical models have been proposed in the literature. In this paper, we study a mathematical model of the coronavirus disease with lockdown by employing the Caputo fractional-order derivative. We establish the existence and uniqueness of the solution to the model. We also study the local and global stability of the disease-free equilibrium and endemic equilibrium solutions. By using the residual power series method, we obtain a fractional power series approximation of the analytic solution. Finally, to show the accuracy of the theoretical results, we provide some numerical and graphical results.


Introduction and Preliminaries
The study of mathematical models of infectious diseases has attracted the attention of many researchers since they provide a better understanding of their evaluations, existence, stability, and control [1][2][3][4][5][6]. Most of the mathematical models of infectious diseases are composed of a system of integer-order differential equations. However, in the last few decades, fractional-order differential equation has been used in the modeling of biological phenomena because they provide a greater degree of accuracy than the integer-order models [7][8][9][10][11][12][13]. The theory of fractional calculus, which involves differentiation and integration of non-integer orders, is as old as the classical calculus of integer orders. However, this theory has gained considerable attention from many researchers in recent years due to the numerous applications found in the sciences, engineering, economics, control theory, and finance, amongst others. The increasing interest in using fractional calculus in the modeling of real-world phenomena is due to its various properties which are not found in classical calculus. Unlike the integer-order derivatives, which are local in nature, most of the fractional-order derivatives (for example, the Caputo and Riemman-Liouville fractional order derivatives [14,15]) are non-local and possess the memory effects which make it more superior because in many situations the future state of the model depends not only upon the current state but also on the previous history. For this realistic property, the usage of fractional-order systems is becoming popular to model the behavior of real systems in various fields of science and engineering. Many researchers have therefore expanded the integer-order models to fractional-order models via various mathematical techniques. Recently, several authors have used fractional-order models to understand the dynamics of the COVID-19 epidemic to achieve the best strategies to stop the spread of the disease, chose a better effective immunization program, allocate scarce resources to control or prevent infections and also predict the future course of the outbreak [16][17][18][19][20][21][22][23][24][25][26][27].
Motivated by the above works and the high interest in finding the best solution to control the COVID-19 pandemic, we study a fractional-order model to understand the dynamics of the COVID-19 infections when the population is under lockdown. Most of the COVID-19 models in the literature, to the best of our knowledge, do not include the impact of the preventive measures adopted by many governments around the world to control the virus. Lockdown was one of the preventive measures and that makes our model quite unique and interesting.
Several definitions of fractional derivatives and integrals have been provided in the literature. For the purpose of this work, we present the definitions and some properties of the Riemann-Liouville fractional integral, Caputo fractional derivative and other related fractional-order derivatives. For more information about the Riemann-Liouville fractional integral and the Caputo derivative, we refer the interested reader to [14,15] and the references therein. Definition 1 ([14]). The Riemann-Liouville fractional integral of order α > 0 of a function f is given by where Γ denotes the Gamma function defined by Definition 2 ([14]). The Caputo fractional derivative of order α > 0 is given by D α a f (t) = 1 Γ(n − α) t a f (n) (z) (t − z) α−n+1 dz, t > a, n − 1 < α ≤ n, n ∈ N. (2) The Caputo derivative satisfies the following properties: D α a (J α a f (t)) = f (t). 3.
D α a (t − a) k = 0, for k = 0, 1, 2, · · · , n − 1. Most fractional-order differential equations do not have a closed-form solution, and thus approximation and numerical methods are extensively used. One way to find an approximated solution of a system of fractional-order differential equations is by using the technique of fractional power series. The following is the definition of the fractional power series.

Definition 3 ([28]
). The fractional power series about t = t 0 is defined as where n − 1 < α ≤ n for some n ∈ N and c m for m = 0, 1, 2, · · · are the coefficients of the power series.
The following result hold for the Caputo derivative.
Theorem 1 ([28]). Suppose the fractional power series representation of f is of the form If D mα a f (t) for m = 0, 1, 2, 3 · · · are continuous on the interval (t 0 , t 0 + ρ), then , where D mα a = D α a D α a · · · D α a (m-times) and ρ is the radius of convergence.

Remark 1.
We note that the kernel function in the definition of the Caputo derivative (t − z) n−1−α has a singularity at z = t. Recently, some new fractional derivatives with non singular kernels were introduced in the literature and we present them here for the readers' reference.
In 2015, Caputo annd Frabrizio [29] proposed the following generalization of the Caputo derivative as follows: Definition 4. The Caputo-Fabrizio fractional-order derivative (CF) of order α is defined as follows: where M(α) is a normalization function such that M(0) = M(1) = 1. The Caputo-Fractional derivative of a constant function is zero, but unlike the Caputo derivative in (2), the kernel does not have a singularity at z = t.
The Caputo-Fabrizio fractional-order had also been generalized by Atangana and Baleanu in 2016 as follows: 30]). Let f ∈ H 1 (a, b), a < b and α ∈ [0, 1]. The Atangana-Baleanu fractionalorder derivative in the Caputo sense is given by where B(α) is a normalization function such that B(0) = B(1) = 1, E α (·) is the one parameter Mittag-Leffler function (see Remark 3) and H 1 (a, b) is the Sobolev space of order one defined as follows: We refer the interested reader to [29,30] and the references therein for more information about these fractional derivatives.

Remark 2.
It is worth pointing out that recently, Hattaf [31] also introduced some new generalizations of the Attangana-Baleanu fractional-order derivatives by introducing a weight function and thus giving rise to the definition of several fractional-order derivatives with non-singular kernels.
Next, we present the definition of the two-parameter Mittag-Leffler function, which will be utilized later in the paper.

Model Formulation
In order to understand the impact of lockdown in preventing the spread of COVID-19 infection, Baba et al. [33] proposed the following model governed by a system of nonlinear ordinary differential equations: subjected to the initial conditions In this model, the authors considered a population of total size N(t) at time t with a constant recruitment rate of Λ. The population is divided into four compartments denoted by S(t), S L (t), I(t) and I L (t). The S(t) class denotes the susceptible individuals that are not under lockdown. The group S L (t) contains those individuals who are susceptible and are under lockdown. The group who are infected and are not under lockdown are represented by I(t), while those individuals who are infected and are under lockdown are denoted by I L (t). Finally, the cumulative density of the lockdown program is denoted by L(t). The authors have studied the local stability of the equilibrium solutions in relation to the basic reproduction number.
In order to include the memory effects and the past history to get a better understanding of the dynamics of COVID-19 infections under lockdown, we reformulate the model (7) by using the Caputo fractional derivative as follows: where D α 0 denotes the Caputo derivative for 0 < α ≤ 1; under the initial conditions We explore some interesting results of the fractional-order COVID-19 model (8), description of the parameters shown in Table 1. In particular, after proving the existence of a unique positive global solution of the model, we study the local and global stability of the various equilibrium points of the fractional-order model using the comparison theory of fractional differential inequality and fractional La-Salle invariance principle. Moreover, we apply the method of residual power series technique to approximate the solution of the fractional-order system (8). Finally, we provide numerical simulations to illustrate some of the theoretical results and error analysis to show how good the approximated solution is.

Remark 4.
In [17], the authors also discussed a slightly different extension of the model in (7) to fractional calculus where they include the order of the fractional derivative as powers on the parameters. They established the existence and uniqueness of solutions to their model using the Schauder and Banach fixed point theorems.

Theorem 2.
For any t ≥ 0 and X(0) ∈ Ω there exist a unique positive solution X(t) = (S(t), S L (t), I(t), I L (t), L(t)) of system (8). Moreover, the set Ω attracts all solutions of system (8), and thus it is positively invariant.
Proof. Denote the vector field associated with system (8) by f (X(t)). Then f (X(t)) can be written as f (X(t)) = A + ( Clearly, f (X(t)) and ∂ f (X) ∂X are continuous for all X(t) ∈ Ω. Moreover, where k 1 = ||A|| and k 2 = max{1, Λ d }||B + C + D|| are two positive constants. Thus by Theorem (3.1) and remark (3.2) of [34], the system (8) has a unique solution. Now, to prove that Ω is positively invariant, let N(t) = S(t) + S L (t) + I(t) + I L (t). Then adding the first four equations of system (8), we have Taking the Laplace transform on both sides yields Simplifying this equation, we have the following inequality Taking the inverse Laplace transform and using the fact that Consequently, inequality (9) and the last equation of system (8), implies that . Similar to the above discussion, taking the Laplace transform and using the Mittag-Leffler function, it follows that L(t) ≤ µΛ φd for any t > 0. In conclusion, for any X(0) ∈ Ω we have N(t) ≤ Λ d and L(t) ≤ µΛ φd for any t ≥ 0, and thus Ω is positively invariant.

Equilibrium Points and Basic Reproduction Number
One way to see what will happen to the population eventually is to explore when the system is at equilibrium. By setting we get the following equilibrium points and The disease-free equilibrium E 0 is the case when the pathogen has suffered extinction and, in the long run, everyone in the population is susceptible. The endemic equilibrium E 1 is the state where the disease cannot be totally eradicated and remains in the population without a lockdown, while E 2 is the endemic equilibrium point in the presence of lockdown. The basic reproduction number of the model, denoted by R 0 , is a constant that is used to approximate the expected number of cases directly generated by one case in a population where all individuals are susceptible to infection. One way to calculate R 0 is using the next generation matrix approach [35]. The rate at which secondary infections are produced is given by Similarly, the transfer of infection from compartment i to j is given by where i, j ∈ {1, 2} and x j ∈ {I, I L } for j = 1, 2. Now define the next generation matrix G as Hence R 0 is the dominant eigenvalue of G and it is given by

Local and Global Asymptotic Stability of the Disease-Free and Endemic Equilibrium Points
We use the following theorem to prove the local asymptotic stability of the disease-free equilibrium point E 0 .

Theorem 3 ([36]
). Given a fractional-order system of differential equation Let x 0 be an equilibrium point of the given system, and let A = D( f (x 0 )) be the Jacobian matrix of f evaluated at x 0 . Then x 0 is locally asymptotically stable if and only if |arg(λ i )| > απ 2 , for all eigenvalues λ i of the matrix A.
Theorem 4. The disease-free equilibrium point E 0 is locally asymptotically stable if R 0 < 1.

Proof.
After calculating the associated Jacobian matrix of system (8), it can be shown that the characteristic equation of the Jacobian matrix satisfies As a result, all the eigenvalues are negative real numbers, and hence |arg(λ)| = π > απ 2 for any 0 < α ≤ 1. Thus by Theorem 3 the disease-free equilibrium point E 0 is locally asymptotically stable.

Remark 5.
Since Mittage-Leffler stability implies global asymptotic stability (see Remark 4.4 in [38]), it is sufficient to prove the disease-free equilibrium point E 0 is Mittage-Leffler stable on Ω.
Note that a similar result is provided in Definition 5 [39] for a generalized Hattaf fractional (GHF) derivative which encloses the popular forms fractional derivatives with non-singular kernels.
Proof. Define the positive definite Lyapunov function where c 1 and c 2 are two positive constants that will be determined later in the proof. By linearity of the Caputo fractional derivative, we have Now since then Equation (15) can be rewritten as If R 0 = βΛ d(γ 1 +α 1 +d) < 1 then we can choose c 1 , c 2 > 0 such that Thus the inequality in (16) becomes where Now taking the Laplace transform of inequality (17) results Thus using (2), we have Let ||x(t)|| = |S(t)| + |S L (t)| + |I(t)| + |I L (t)| + |L(t)| be the norm defined on the solution x(t) = (S(t), S L (t), I(t), I L (t), L(t)) of system (8). Then from Equation (14) it follows that where c 4 = min{1, c 1 , c 2 } and c 5 = max{1, c 1 , c 2 }. As a result, from inequality (19) and (20) it follows that ||x(t) , where M(S, S L , I, I L , L) := c 4 c 5 (S(t) + S L (t) + I(t) + I L (t) + L(t)). In conclusion, the disease-free equilibrium point is Mittage-Leffler stable, and thus it is globally asymptotically stable. Remark 6. The proof of the local stability of the endemic equilibrium point E 1 is similar to Theorem 4. Thus we will provide proof of the global stability of E 1 by constructing a Lyapunov function.

Theorem 6 ([40])
. Let x 0 ∈ Γ be an equilibrium point for the non-autonomous fractional-order and D α 0 L(t, x(t))) ≤ −W 3 (x) for all α ∈ (0, 1) and all x ∈ Γ, where W 1 , W 2 , W 3 are continuous positive definite functions on Γ. Then the equilibrium point x 0 is globally asymptotically stable. Now using Theorem 6 we will show that the endemic equilibrium point E 1 of system (8) is globally asymptotically stable under some conditions. Theorem 7. If R 0 > 1, then E 1 is globally asymptotically stable in Ω. where and k is a positive constant to be determined later in the proof. Using Lemma 3.1 and Remark 3.1 of [41] we have and by the linearity property of the fractional Caputo derivative, it follows that Now using Equations (10), (22)- (24) and the fact that Choose any k > S 1 λ 2 φ and note that if R 0 > 1 then is the singleton set containing the endemic equilibrium E 1 . Thus by Theorem 6, we conclude that the endemic equilibrium is globally asymptotically stable in Ω.
Since the proof of the local stability of the endemic equilibrium point with lockdown E 2 is similar to the proof of (4), we will provide the proof of the global stability of the endemic equilibrium point with lockdown E 2 .

Approximate Solution
In this section, we present a solution of the fractional-order model by using the residual power series method, which consists of expressing the solution as fractional power series expanded about the initial point t = t 0 .
Consider the following fractional-order model: where D α 0 denotes the Caputo derivative for 0 < α ≤ 1 with h = d + θ 1 , p = γ 1 + α 1 + d and r = d + θ 2 + γ 2 + α 2 ; under the initial conditions Similar to the procedure in [42], we do the following steps in order to obtain a fractional power series solution for the nonlinear fractional-order model in (26).
Step 1: Suppose that the fractional power series of S(t), S L (t), I(t), I L (t) and L(t) around t = 0 are as follows: where 0 ≤ t < η for some η > 0. Now, we let S n (t), S L,n (t), I n (t), I L,n (t) and L n (t) denote the n-th truncated power series approximation of S(t), S L (t), I(t), I L (t) and L(t), respectively. That is, Step 2: Next, we define the residual functions for the model in (26) as follows: Hence, the n-th residual functions for S(t), S L (t), I(t), I L (t) and L(t), respectively, are as follows: Res S n (t) = D α 0 S n (t) − Λ + βS n (t)I n (t) + λ 1 S n (t)L n (t) + dS n (t) − γ 1 I n (t) − γ 2 I L,n (t) − θ 1 S L,n (t) Res S L,n (t) = D α 0 S L,n (t) − λ 1 S n (t)L n (t) + hS L,n (t) Res I n (t) = D α 0 I n (t) − βS n (t)I n (t) + pI n (t) + λ 2 I n (t)L n (t) − θ 2 I L,n (t) (30) Res I L,n (t) = D α 0 I L,n (t) − λ 2 I n (t)L n (t) + rI L,n (t) Res L n (t) = D α 0 L n (t) − µI n (t) + φL n (t).
Now, we observe that Res S n (t) = Res S L,n (t) = Res I n (t) = Res I L,n (t) = Res L n (t) = 0, and lim n→∞ Res S n (t) = Res S (t), lim n→∞ Res S L,n (t) = Res S L (t), lim Since the Caputo derivative of any constant is zero, it is straightforward to see that for Res L n (0).
Step 3: To obtain the coefficients a k , b k , c k , d k and e k , for k = 1, · · · , n, we substitute the n-th truncated series of S(t), S L (t), I(t), I L (t) and L(t) into (30) and then apply the Caputo fractional derivative operator D (n−1)α 0 on S n (t), S L,n (t), I n (t), I L,n (t) and L n (t), and evaluate the result at t = 0. We obtain the following algebraic system of equations. for n = 1, 2, 3, · · · .
Step 4: We solve the algebraic system (31) for the values of a k , b k , c k , d k , e k for k = 1, 2, 3, · · · , n to get the n-th residual power series approximate solution of the system in (26).
Step 5: We repeat the procedure to obtain a sufficient number of coefficients. Higher accuracy for the solution can be achieved by evaluating more terms in the series solution.
By following the steps outlined in above, we derive a recursive formula for the coefficients as follows: First, we note that the coefficients a 0 , b 0 , c 0 , d 0 and e 0 are given by the initial conditions. That is, The first truncated power series approximations will have the forms: Thus, the first residual functions are: Evaluating Res S 1 (t), Res S L,1 (t), Res I 1 (t), Res I L,1 (t) and Res L 1 (t) at t = 0, we have By solving the equations Res S 1 (0) = 0, Res S L,1 (0) = 0, Res I 1 (0) = 0, Res I L,1 (0) = 0 and Res L 1 (0) = 0, we have that Now, the second truncated power series approximations will have the forms: t 2α , and We apply the operator D α 0 to Res S 2 (t), Res S L,2 (t), Res I 2 (t), Res I L,2 (t) and Res L 2 (t) and then evaluate the result at t = 0 to get By solving the equations we have Similarly, the third truncated power series approximations are given by t 3α , and So, the third residual functions are: We apply the operator D 2α 0 to Res S 3 (t), Res S L,3 (t), Res I 3 (t), Res I L,3 (t) and Res L 3 (t) and then evaluate the result at t = 0 to get By solving the equations we have By continuing in this direction, we deduce that the coefficients a n , b n , c n , d n and e n for n ≥ 2 are given by the recursive formula a n = −βΓ(1 + (n − 1)α) Remark 7. It is worth noting that these recursive formulas for the coefficients will ensure that we obtain higher-order approximate solutions as compared with similar results in the literature with a smaller order of the approximate solutions. Thus, due to this recursive nature of the coefficients, we can calculate a n , b n , c n , d n and e n for any large n values if necessary. For the purpose of the numerical simulation given in Section 6, we used n = 25. That is, the coefficients a n , b n , c n , d n , e n for n = 1, 2, 3, . . . 25 are computed and the approximated fractional power series solutions of system (8) are plotted in section 6, numerical simulation and examples.

Numerical Simulations and Examples
The following table provides the description of the different parameters used in the COVID-19 model. The parameter values under column "data set 1" yield a basic reproduction number R 0 = 0.2373. Similarly, if the parameters under the column "data set 2" are used, then we have R 0 = 4.5707. In Section 6.1, we simulate both the exact and approximated solutions of system (8) using the values under "data set 1" of Table 2 for various α values (α = 0.5, 0.75, 0.95 and 1). In Section 6.2, we simulate the exact and approximated solutions of system (8) using the parameter values under "data set 2" of table Table 2. Moreover, in both cases (R 0 < 1 and R 0 > 1) the relative errors, due to approximating the exact solution of system (8) by the residual power series method, are computed for different α values as indicated in Tables 3 and 4. Using data set 1, the basic reproduction number R 0 is calculated to be 0.2373. Thus by Theorems 4 and 5, the disease-free equilibrium point is both locally and globally asymptotically stable. The disease-free equilibrium point is calculated to be E 0 = (4167, 0, 0, 0, 0). The simulations in Figures 1 and 2 support the fact that the trajectories of system (8) converges to this equilibrium point E 0 = (4167, 0, 0, 0, 0).  In Figures 1 and 2 the exact and approximated solutions of system (8) are plotted for 0 ≤ t ≤ 100, respectively. In order to show the accuracy of the residual power series method for approximating the solution of system (8), we have calculated the relative error of each state variable for various α values. A time step size ∆t = 2 −6 is considered and for each state variables x(t) ∈ {S(t), S L (t), I(t), I L (t), L(t)} the values x(i∆t) are calculated where 0 ≤ i ≤ 6400. The relative error for the state variable x(t) is computed using where Exact(x(i∆t)) and Appr(x(i∆t)) represent the exact and approximated values of the variable x(t) at t = i∆t for 0 ≤ i ≤ 6400. As can be seen from Table 3, the relative error of each state variable decreases as the value of α increases from 0.5 to 1. The basic reproduction number is calculated using the parameter values under data set 2. The value of R 0 is found to be 4.5707, and also the condition λ 2 < θ 2 is satisfied. Thus by Remark 6 and Theorem 7 the endemic equilibrium point E 1 = (1458, 0, 3370, 0, 0), which is calculated using Equation (10), is both locally and globally asymptotically stable. Note that both Figures 3 and 4 show that the trajectories of the solution of system (8) converges to the equilibrium point E 1 .  Similar to the previous discussion, the relative errors of the state variables are calculated for various α values using Equation (33).

Conclusions and Discussion
In this paper, a fractional-order COVID-19 epidemic model with lockdown is proposed and analyzed. In particular, some of the most important epidemiological constants, such as the equilibrium points (E 0 , E 1 and E 2 ) and the basic reproduction number R 0 , are calculated. The existence and uniqueness of a positive global solution of the fractional-order system (8) is established. Moreover, by implementing various techniques such as the comparison theory of fractional-order differential equations, Mittage-Leffler stability (Theorem 5), and fractional La-Salle invariance principle (Theorem 6), the local and global stabilities of the equilibrium points are discussed. Additionally, the method of residual power series is used to approximate the solution of system (8). As can be seen in Equation (32), the coefficients a n , b n , c n , d n and e n have recursive nature. That is, for any n ≥ 1, each coefficients a n , b n , c n , d n and e n can be determined from the previous ones a n−1 , b n−1 , c n−1 , d n−1 and e n−1 . Thus we can calculate as many coefficients as we need in order to get a very good approximation of the solution of system (8). Moreover, it is very interesting to note that as the value of α increases from 0.5 to 1, the trajectories of the solution of the fractional-order system are converging to the case where α = 1. That is, the solutions of the fractional-order differential equation converge to the solutions of the system of the ordinary differential equation as the order α increases towards 1. In Section 5, numerical simulations and examples are presented to show the stability of the equilibrium points when R 0 < 1 and R 0 > 1. The exact solutions of system (8) for various α values are shown in Figures 1 and 3. Similarly, the approximated solutions of system (8) using the residual power series method are included in Figures 2 and 4. Finally, in order to see the accuracy of the approximated solutions, we provided the relative errors for different α values as shown in Tables 3 and 4. The model that is studied in this paper can be used to predict and understand how COVID-19 infection spreads in the presence of a lockdown. In addition to understanding the dynamics of the infection, it is very important to determine the relative importance of the various parameters used in the model. One such technique is the study of the sensitivity of the parameters in relation to the basic reproduction number R 0 , and the two endemic equilibrium points E 1 and E 2 , which are considered to be the most important values in the study of deterministic epidemic models. One such approach is to calculate the normalized sensitivity index (NSI), denoted by SI κ , where κ is the given parameter. For example, if κ is one of the parameters in R 0 , then the NSI of κ is defined as SI κ = κ R 0 ∂R 0 ∂κ . Using this result, we have calculated the NSI of the parameters β and γ 1 , which are directly affected by the infection. As a result, we have SI β = 1, and SI γ 1 = −γ 1 γ 1 +α 1 +d < 0. Thus in order to decrease R 0 , the expected number of infections generated by one case, we have to either decrease the infection contact rate β or increase the recovery rate of the infected group γ 1 . This further implies imposing some prevention actions, such as avoiding contact with people who are infected, wearing masks, and keeping a distance from an infectious person, will help to decrease the value of R 0 , which directly leads to mitigating the infection. By modifying this model, one can study other aspects of the infection, such as the impact of vaccination, the effect of population diffusion, and other important factors that determine the transmission and persistence of the infection.