Dynamics of a fractional order mathematical model for COVID-19 epidemic

In this work, we formulate and analyze a new mathematical model for COVID-19 epidemic with isolated class in fractional order. This model is described by a system of fractional-order differential equations model and includes five classes, namely, S (susceptible class), E (exposed class), I (infected class), Q (isolated class), and R (recovered class). Dynamics and numerical approximations for the proposed fractional-order model are studied. Firstly, positivity and boundedness of the model are established. Secondly, the basic reproduction number of the model is calculated by using the next generation matrix approach. Then, asymptotic stability of the model is investigated. Lastly, we apply the adaptive predictor–corrector algorithm and fourth-order Runge–Kutta (RK4) method to simulate the proposed model. Consequently, a set of numerical simulations are performed to support the validity of the theoretical results. The numerical simulations indicate that there is a good agreement between theoretical results and numerical ones.

posed class), I (infected class), Q (isolated class), and R (recovered class). This model is a generalization of a well-known ODE model formulated in [35]. In the proposed fractionalorder model, we use the Caputo fractional derivative instead of the integer-order one because when modeling processes and phenomena arising in the real world, fractional-order models are more accurate than integer-order ones. In particular, fractional-order models provide more degrees of freedom in the model while an unlimited memory is also guaranteed in contrast to integer-models with limited memory [11,24,27].
Our main aim is to study the dynamics and numerical approximations for the proposed fractional-order model. Firstly, the positivity and boundedness of the model are investigated based on standard techniques of mathematical analysis. Secondly, the basic reproduction number of the model is calculated by using the next generation matrix approach. Then, asymptotic stability of the model is investigated based on the Lyapunov stability theorem for fractional dynamical systems. Lastly, we apply the adaptive predictor-corrector algorithm and fourth-order Runge-Kutta (RK4) method formulated in [20] to simulate the proposed model. Consequently, a set of numerical simulations is performed to support the validity of the theoretical results. The numerical simulations indicate that there is a good agreement between theoretical results and numerical ones.
This work is organized as follows. The fractional-order differential model is formulated in Sect. 2. Dynamics of the model is investigated in Sect. 3. Numerical simulations by the adaptive predictor-corrector algorithm are performed in Sect. 4. The last section present some remarks, conclusions, and discussions.

Model formulation
The total population is divided into five compartments: susceptible (S), exposed (E), infected (I), isolated (Q), and recovered (R) from the disease. We assume that all the classes are normalized. This leads to the mathematical model formulation in which the interaction of the exposed population and infected population is linked to the susceptible population. In this model, we assume that the natural death rate includes the disease death rate. When there is no symptom of the disease, the exposed class moves with a certain rate to the isolated class, but in the case when symptoms are developed, it moves to the infected class. Keeping in mind the above assumptions, we obtain the following ODE model (see [35]): where the parameters and variables used are described in Table 1. Let us introduce the notation which implies that the total population in this case we take as constant because Λ = μN .
To include into the model (1) the past history or hereditary properties, we replace the first derivative by the Caputo fractional derivative. More precisely, we propose the following system of fractional differential equations: where 0 < α < 1, and C 0 D α t denotes the fractional derivative in the Caputo sense. We refer the readers to [2-4, 10, 25] for the definition of the fractional Caputo derivative.

Positivity and boundedness
Let us denote The following theorem is proved by using a generalized mean value theorem [ Proof First, it is easy to prove the existence and unique of solutions of the model thanks to results proved in [17]. For the model (3), we have By (5) and the generalized mean value theorem [21,22], we deduce that S(t), E(t), I(t), From the first equation of the system (3), we obtain By using the fractional comparison principle, we have the first estimate of (4). From the second equation of the system (3), we have Consequently, we have the second estimate of (4). From the third equation of the system (3), we get for t large enough. This implies the third estimate of (4).
From the fourth equation of the system (3), we have for t large enough. From this we get the fourth estimate of (4).

Finally, the fifth equation of (3) implies that
which implies the fifth estimate of (4). The proof is complete.

Equilibria and the reproduction number
Firstly, to find equilibria of the model (3), we consider the following algebraic system: By some algebraic manipulations, we obtain two solutions of the system (6) that are and Note that S * , E * , I * , Q * and R * are positive if and only if We now compute the reproduction number of the model (3) by using the next generation matrix approach developed by van den Driessche and Watmough [31]. Let x = (E, I, Q, R, S). We rewrite the model (3) in the matrix form where Hence, the reproduction number of the model (3) can be determined by . (11) It is easy to verify that S * , E * , I * , Q * , and R * > 0 if and only if R 0 > 1.
given by (8) for all values of the parameters. Moreover, the model has a unique disease endemic equilibrium point F * = (S * , E * , I * , Q * , R * ) given by (9) if and only if R 0 > 1. (3) is locally asymptotically stable if R 0 < 1.

Stability analysis Theorem 3 The DFE point of the model
Proof The Jacobian matrix of the model (3) at the DFE point is The characteristic polynomial of J is where τ 1 = -βS 0π α + μ α + γ ασ α + μ α , It is easy to verify that if R 0 < 1, then τ 1 > 0 and τ 2 > 0. This implies that two roots of the polynomial λ 2 + τ 1 λ + τ 2 have negative real parts. Consequently, the real parts of five eigenvalues of the matrix J(F 0 ) are all negative, or equivalently, F 0 is locally stable. The proof is complete.
We now prove the uniform asymptotical stability of the DFE point of the model (3) (3) is not only locally asymptotically stable but also uniformly asymptotically stable.
Proof Since the first three equations of the model (3) do not depend on the states Q and R, we need only consider the following subsystem: From (4), it is sufficient to consider the model (3) in its feasible set defined by Consider a Lyapunov function V : Ω → R + given by By using the linearity property of the Caputo derivative and [32, Lemma 3.1] and some algebraic manipulations, we obtain It is easy to check that if R 0 < 1, then τ 1 , τ 2 < 0. Consequently, by the Lyapunov stability theorem for fractional dynamical systems, we deduce that F 0 is uniformly asymptotically stable. The proof is complete.
Remark 1 The analysis of stability of F * is an interesting mathematical problem, but in this work, we mainly focus on the case R 0 < 1 to find an effective strategy to prevent the disease.

R 0 sensitivity analysis
Theorem 3 suggests that we should control the parameters such that R 0 < 1. This provides a good strategy to prevent and restrain the disease. More precisely, when R 0 < 1, then which means that the disease will be completely controlled and prevented. Motivated by this, we now perform an R 0 sensitivity analysis to find ways to choose suitable parameters. It is easy to verify that Equation (16) suggests some ways to choose the parameters such that R 0 < 1. Hence, based on this, we can propose suitable strategies to control and prevent the disease.

Numerical simulations by the adaptive predictor-corrector algorithm 4.1 The adaptive predictor-corrector algorithm
In this section we review the method that is proposed by [20]. The proposed algorithm is given as follows. Consider the initial value problem (IVP) governed by: where D α,ρ a+ is the proposed generalized Caputo-type fractional derivative operator given [20,Definition 4]. Initially, for m -1 < α ≤ m, a ≥ 0, ρ > 0 and y ∈ C m ([a, T]), the IVP (17) is equivalent, using [20,Theorem 3], to the Volterra integral equation: where The where h = (T ρ -a ρ ) N . Now, we are going to generate the approximations y k , k = 0, 1, . . . , N , to solve numerically the IVP (17). The basic step, assuming that we have already evaluated the approximations y i ≈ y(t j ) (j = 1, 2, . . . , k), is that we want to get the approximation y k ≈ y(t k+1 ) by means of the integral equation Making the substitution we get that is, Next, if we use the trapezoidal quadrature rule with respect to the weight function (t ρ k+1z) α-1 to approximate the integrals appear in the right-hand side of Eq. (24), replacing the function f (z 1 ρ , y(z 1 ρ )) by its piecewise linear interpolant with nodes chosen at the t ρ j (j = 0, 1, . . . , k + 1), then we obtain Thus, substituting the above approximations into Eq. (24), we obtain the corrector formula for y(t k+1 ), k = 0, 1, . . . , N -1, where The last step of our algorithm is to replace the quantity y(t k+1 ) shown on the right-hand side of formula (26) with the predictor value y p (t k+1 ) that can be obtained by applying the one-step Adams-Bashforth method to the integral equation (23). In this case, by replacing the function f (z 1 ρ , y(z 1 ρ )) with the quantity f (t j , y(t j )) at each integral in Eq. (24), we get Therefore, the present adaptive predictor-corrector algorithm, for evaluating the approximation y k+1 ≈ y(t k+1 ), is completely described by the formula where y j ≈ y(t j ), j = 0, 1, . . . , k, and the predicted value y p k+1 ≈ y p (t k+1 ) can be determined as described in Eq. (28) with the weights a j,k+1 being defined according to (27). It is easily observed that when ρ = 1, the present adaptive predictor-corrector algorithm will be reduced to the predictor-corrector approach presented in [13].  Figs. 1-5, it can be seen that the results obtained using the proposed algorithm match the results of the RK4 method very well, which implies that the presented method can predict the behavior of these variables accurately in the region under consideration. Figures 6-10 show the approximate solutions for S(t), E(t), I(t), Q(t), and R(t) obtained for different values of α using the proposed algorithm. From the graphical results given in Figs. 6-10, it is clear that the approximate solutions depend continuously on the timefractional derivative α. It is evident that the efficiency of this approach can be dramatically enhanced by decreasing the step size.

Remarks and conclusions
In this paper, we have formulated and analyzed a new mathematical model for COVID-19 epidemic. The model is described by a system of fractional-order differential equations and includes five classes, namely, S (susceptible class), E (exposed class), I (infected class), Q (isolated class), and R (recovered class). It should be emphasized that the model is a generalization of our recent work proposed in [35].
Firstly, the positivity, boundedness, and stability of the model have been established. Furthermore, the basic reproduction number of the model has been calculated by using the next generation matrix approach. Lastly, we have applied the adaptive predictor-corrector algorithm and fourth-order Runge-Kutta (RK4) method to simulate the proposed model.
A set of numerical simulations has been performed to support the validity of the theoretical results. The numerical simulations indicate that there is a good agreement between theoretical results and numerical ones.
In the near future, we will extend the results in this work to propose new mathematical models for COVID-19 epidemic. Especially, effective strategies to control and prevent the disease will be investigated.