Dynamics of a Fractional-Order COVID-19 Epidemic Model with Quarantine and Standard Incidence Rate

: In this paper, we propose a fractional-order COVID-19 epidemic model with a quarantine and standard incidence rate using the Caputo fractional-order derivative. The model consists of six classes: susceptible ( S ), exposed ( E ), infected ( I ), quarantined ( Q ), recovered ( R ), and deceased ( M ). In our proposed model, we simultaneously consider the recovery rate and quarantine rate of infected individuals, which has not been considered in other fractional-order COVID-19 epidemic models. Furthermore, we consider the standard incidence rate in the model. For our proposed model, we prove the existence, uniqueness, non-negativity, and boundedness of the solution. The model has two equilibrium points: disease-free equilibrium and endemic equilibrium. Implementing the spectral radius of the next-generation matrix, we obtain the basic reproduction number ( R 0 ). The disease-free equilibrium always exists and is locally and globally asymptotically stable only if R 0 < 1. On the other hand, endemic equilibrium exists and is globally asymptotically stable if R 0 > 1. Our numerical simulation conﬁrms the stability properties of the equilibrium. The smaller the order of the derivative, the slower the convergence of the solution of the model. Both the recovery rate and quarantine rate of the infected class are important parameters determining the stability of the equilibrium point. Based on parameter estimation from COVID-19 data in Indonesia, the fractional-order model has better performance than the ﬁrst-order model for both the calibration and 20-day forecasting of conﬁrmed daily active cases of COVID-19.


Introduction
Coronavirus Disease 2019, or COVID-19, is an acute respiratory disease that was first indicated in Wuhan Province, China, in 2019 [1]. COVID-19 is an infectious disease caused by infection with the SARS-CoV-2 virus [2]. Symptoms of this disease are fever, flu, cough, pneumonia, and other symptoms that cause complications such as asthma [1,3]. COVID-19 can spread quickly, and the total number of cases reached 624 million cases in the world [4]. The number of cases is also influenced by the variant of the COVID-19 virus. In general, there are five variants of COVID-19 virus: Alpha, Beta, Gamma, Delta, and Omicron [5,6]. Various studies on how COVID-19 is transmitted have been described in [7][8][9]. The mechanism of transmission most often occurs through droplets, directly or indirectly [2,10,11]. To prevent the spread of the virus, quarantine is an effective intervention [11][12][13]. Quarantine causes infected individuals to lose contact with other individuals, which can wipe out the transmission of the disease [14]. Thus, the mechanism of transmission with quarantine interventions is important to understand and study further.
Mathematical modeling is an approximation used to understand the mechanism of transmission and current state of COVID-19 [15]. If COVID-19 data are sufficient, then the forecasting of COVID-19 via parameter estimation is very helpful for policy makers, as was conducted in [16][17][18]. Otherwise, the dynamics of an epidemic model can provide conditions that lead to an endemic disease. Some researchers have studied the COVID-19 pandemic via the dynamics of mathematical models, such as in [15,19,20]. The mathematical model for a disease epidemic is basically represented by first-order differential equations.
Some researchers have developed COVID-19 epidemic models. Zeb et al. [21] proposed a COVID-19 epidemic model with a quarantined class. The paper of Zeb et al. [21] focuses on the dynamics of the COVID-19 epidemic model. In the model, there is an assumption that exposed individuals can transmit the COVID- 19 virus. The current literature shows that an exposed individual is an infected individual but does not always transmit the virus [22,23]. Furthermore, the mortality of COVID-19 is not considered in the model of Zeb et al. [21]. Another COVID-19 epidemic model was proposed by López and Rodó [24] with a quarantined class and protected class. They focus on the calibration of COVID-19 data in Italy and Spain. In their model, there is an assumption that the death class only comes from the quarantined class. Musafir et al. [15] stated that infected individuals can also cause death from COVID-19. In both models, the incidence rate used is bilinear, with a maximum proportion of the total of living classes. A recent model has been proposed to resolve the weak assumption for these models. The proposed model by Darti et al. [25] assumes that the exposed class is not contagious and the infected class may induce death from COVID-19. These assumptions improve the proposed models by Zeb et al. [21] and López and Rodó [24]. All three models are represented by nonlinear first-order differential equations.
In recent decades, fractional calculus has been developed to provide some integral and derivative options [26,27]. A fractional derivative is a non-local operator, which means the rate of the function depends on previous histories from the initial value, which is called the memory effect [26]. In epidemiology, the future spread of a disease depends on recent and previous conditions because of the impact of both heredity and history on individual habits [28]. Hence, the fractional derivative is more realistic than the integer derivative. Moreover, some researchers, such as [2,20,29,30], state that fractional derivatives provide better accuracy and flexibility than integer derivatives. There are two types of fractional operators obtained by generalizing the integer-order derivative, namely, the Riemann-Liouville operator and the Caputo operator [31]. Li et al. [32] show that the Caputo operator is better than the Riemann-Liouville operator in terms of applied mathematics.
Recently, many researchers have proposed a fractional-order COVID-19 epidemic model. Yadav et al. [33] introduced a fractional-order SEIR-type model by dividing the infected class into two subclasses: the infected subclass without symptoms and the infected subclass with symptoms. The proposed model was analyzed dynamically and solved using q-homotopy analysis with the Sumudu transform method. With the increasing incidences of COVID-19, various measures have been taken by many countries to control the spread of COVID-19. One of the well-known preventive measures is a lockdown. The impact of lockdowns has been modeled by [34], i.e., by considering a fractional-order SI-type epidemic model where the population is split into the subpopulation not under lockdown and the subpopulation under lockdown. In [35], the authors considered a fractional SEIR model incorporating public risk awareness and the spread of disease caused by both humanto-human and zoonotic transmissions. Then, in [36], the authors extended the model of [33] by considering the classes of quarantine, isolation, and environmental impacts.
We notice that although the fractional-order COVID-19 epidemic models above were derived based on the Caputo fractional derivative operator, they did not pay attention to the inconsistency in the time dimension. Such dimensional inconsistency is caused by replacing the first-order derivative with a Caputo fractional-order derivative. Moreover, none of those models have considered the effects of fractional-order derivatives. In this study, we propose a fractional-order COVID-19 epidemic model by including the dimensional issue. We implement the proposed model using real COVID-19 data. Since the available data only consist of the confirmed active COVID-19 cases (the infected subclass), without noting whether this is with or without symptoms, we consider a simpler model than proposed by [36]. For this aim, we modify the COVID-19 epidemic model proposed by Darti et al. [25] to a fractional-order model with a Caputo fractional operator. The proposed COVID-19 epidemic model by Darti et al. [25] is given as follows, with the parameter definitions listed in Table 1: =δI(t) +κQ(t).
(1) In this study, we first propose the fractional-order COVID-19 epidemic model in Section 2. The basic properties of the proposed model are presented in Sections 3 and 4. The equilibrium points of the model and their stability are given in Sections 5 and 6, respectively. To measure how sensitive the model parameters are to COVID-19 transmission, in Section 7, we provide a sensitivity analysis. In Section 8, we numerically simulate the properties of the dynamic analysis of our proposed model. The numerical solutions are obtained by adopting the predictor-corrector scheme proposed by [37]. In Section 9, we implement our proposed model using COVID-19 data in Indonesia. Parameter estimation is performed using the nonlinear least square method. This can be achieved using MATLAB's built-in function lsqcurvefit. Here, we minimize the squared error between numerical solutions for infected subpopulations and confirmed active daily data. Finally, we remark upon the conclusions of this study in Section 10.

Model Formulation
We modify the model proposed by Darti et al. [25] to a fractional-order model with the Caputo operator. The Caputo fractional operator is given as follows.
Definition 1 ([38]). Let t ≥ 0 and α ∈ (0, 1). The Riemann-Liouville integral of order α is defined by where x ∈ C n ([0, ∞), R) and Γ(·) is a Gamma function. Definition 2 ([38]). Let t ≥ 0 and α ∈ (0, 1). The Caputo derivative of order α is defined by Our proposed model has five population classes: the susceptible class (S), the exposed class (E), the infected class (I), the quarantined class (Q), the recovered class (R), and the death class (M). We define variable N as the total of living individuals, which is denoted as N = S + E + I + Q + R. The individuals of the susceptible class (S) who are in close contact with infected individuals have the potential to be exposed individuals. The exposed class is infected but not yet infectious. After the incubation period of the virus, exposed individuals move to the infected class (I), which is composed of the infected individuals. The infected individuals can be either quarantined, recovered, or dead from illness. The quarantined individuals (Q) can be either recovered or dead. The interaction of all classes is depicted in Figure 1. The proposed fractional-order SEIQRM model with the Caputo derivative operator is derived by replacing the first-order derivatives in the SEIQRM model (1) with the Caputo fractional derivatives. However, if we perform this step, then there is a time-dimensional inconsistency on both sides of all the equations because the fractional derivatives of order α have the dimension of (time) −α . To handle the dimension inconsistency, we rescale the parameters and obtain the following fractional-order SEIQRM model: where Λ =Λ α , β =β α , µ =μ α , γ =γ α , σ =σ α , η =η α , δ =δ α , ϑ =θ α , and κ =κ α .

Existence and Uniqueness of Solution
In this section, we explicitly prove the existence and uniqueness of the solution of model (2). (2) be non-negative. The model (2) has a unique solution in R 5 + for all t ≥ 0.

Theorem 1. Let the initial values of model
Proof. If the right-hand side of model (2) is written by In the same way, we can show that . . , 5, then all kernels of f i contract and satisfy the Lipschitz condition.
If both sides of model (2) are integrated by Definition 1, then we obtain Hence, we obtain the iterative scheme Next, we construct a recursive formulation based on the Equation (4) as follows: The norm of the first equation of (5) can be expressed as where (X 1 , X 2 , X 3 , X 4 , X 5 ) = (S, E, I, Q, R). Assume that there is t = T such that It is clear that for every j = 1, 2, . . . , 5, ψ jn (t) → 0 as n → ∞. Therefore, the solution of model (2) exists; that is, it is constructed from the Equation (3).
Next, we show that every solution with the same initial value is unique. Suppose the model (2) has two different solutions with the same initial value, where the other solutions are denoted byŜ,Ê,Î,Q, andR. We obtain The norm of the above equation is so that by using the Lipschitz condition, we obtain The inequality (6) is satisfied if and only if S(t) −Ŝ(t) = 0, which means S(t) =Ŝ(t). This leads to a contradiction, so that S(t) is unique. In a similar way, we can show that E(t), I(t), Q(t), and R(t) are unique.

Non-Negativity and Boundedness of the Solution
In this section, we show the non-negativity and boundedness of the solution of model (2). We also obtain the feasible region of model (2). The following theorem guarantees the non-negativity and boundedness of the solution. (2) subject to non-negative initial values are non-negative and ultimately bounded. The feasible region of model (2) is given by

Theorem 2. All solutions of model
Proof. We first show that S is non-negative. Assume S < 0 for some t > 0. Then, there is at least one τ > 0 such that From model (2), we obtain According to the corollary of the non-decreasing function in [39], S > 0 in (τ, τ + ε) for arbitrary small positive constant ε. This leads to a contradiction. Hence, S ≥ 0 for all t > 0. Similarly, we can prove E, I, Q, and R are non-negative. So, all solutions of model (2) are non-negative.
We next show the boundedness of the solution. From model (2), we have According to comparison theorem in [39], we obtain Hence, all solutions of model (2) are ultimately bounded. Then, we obtain the feasible region of model (2) as follows:

The Equilibrium Points and Basic Reproduction Number
To obtain the equilibrium point of model (2), we set where X = (S, E, I, Q, R) T . Based on the third equation to fifth equation of (7), we obtain If E is substituted to the second equation of (7), then we have either For I = 0, we have E = 0, Q = 0, R = 0, N = Λ µ , and S = Λ µ , from which we obtain disease-free equilibrium point X 0 = ( Λ µ , 0, 0, 0, 0). From X 0 , we establish the basic reproduction number. The basic reproduction number, denoted R 0 , is the average number of secondary individuals infected by one primary infected individual in the susceptible class [40]. The basic reproduction number can be determined by the spectral radius of the next-generation matrix.
Let Z = (E, I, Q) T . From model (2), we have By calculating the spectral radius of next generation matrix, we obtain see [25] for more details. Notice that R 0 is inversely proportional to the recovery rate of infected class η and quarantine rate σ. Therefore, if θ or σ increase, then R 0 decreases. We next determine the second equilibrium point. For S = (γ+µ)(σ+η+δ+µ) βγ N = N R 0 , based on the first equation of (7), we obtain from which we obtain the endemic equilibrium point (S * , E * , I * , Q * , R * ) with which exists if R 0 > 1. X 0 = (S 0 , 0, 0, 0, 0) as the disease-free equilibrium with S 0 = Λ µ , which always exists; 2.

The Stability of Equilibrium Points
In this section, we first study the local asymptotic stability of the equilibrium points of model (2). The following theorem provides the local stability of the equilibrium point.
Theorem 4. The disease-free equilibrium point X 0 is locally asymptotically stable if R 0 < 1 and is unstable if R 0 > 1.
We next study the global asymptotic stability of equilibrium points. The following theorems show the condition of the global stability of the disease-free equilibrium X 0 and the endemic equilibrium X * . Theorem 5. The disease-free equilibrium point X 0 is globally asymptotically stable if R 0 < 1.
Proof. Assume that R 0 < 1. Consider a Lyapunov function L, which is defined by If we apply fractional derivative to L with respect to t [41], we obtain If R 0 < 1, then C D α t L ≤ 0. Moreover, the condition C D α t L = 0 is satisfied if and only if E = 0 and I = 0. The Lemma of the LaSalle Invariance Principle in [42] guarantees E → 0 and I → 0 as t → ∞. Hence, Q → 0 and R → 0 as t → ∞, and thus, S → Λ µ = S 0 . Therefore, X 0 is globally asymptotically stable. Theorem 6. Let the endemic equilibrium X * exists, that is, R 0 > 1. X * is globally asymptotically stable in the domain Ω.
The proof of Theorem 6 uses same Lyapunov function W as in [25] by considering the inequality of the fractional derivative of the Lyapunov function [41]: where we simply change the equal sign in [25] to an inequality sign.

Sensitivity Analysis
In this section, we perform sensitivity analysis of basic reproduction number R 0 . This analysis aims to identify the role of each parameter of model (2) to R 0 , which is represented by the sensitivity index. The normalized forward sensitivity index of R 0 for parameterp is formulated by The greater the absolute value of the index, the more influential the parameter is to R 0 . The sign of the index shows the enlargement and reduction of R 0 , that is, a positive sign enlarges R 0 while a negative sign reduces R 0 . We analytically calculate the sensitivity index for each parameter as follows: , , , , To display the sensitivity index of R 0 for each parameter, we substitute the value of the parameter given in Table 2 into the index formula for each parameter. The sensitivity indices of R 0 are depicted in Figure 2, and the values are given in Table 3. It can be seen that parameter β, σ, and η have important role in reducing R 0 . However, the sensitivity of β is clear because the definition of R 0 is very close to disease transmission β. Thus, we discuss the sensitivity indices of both σ and η, i.e., −0.4585 and −0.4458, respectively. Since both these indices are negative, the greater the values of σ or η, the smaller R 0 is. The results show both parameters σ and η have an important role in reducing the basic reproduction number and thus have an impact on the stability of the equilibrium point.    Table 2.

Numerical Simulation
In this section, we confirm our analytical results by setting the parameter values and implementing a predictor-corrector scheme in [37] to obtain the numerical solution of model (2). The values of the parameters in the first simulation are given by Table 2. Note that the parameter units shown in Table 2 are due to parameter scaling. Therefore, to maintain consistency with the physical dimensions, the units of all model parameters (2) are adjusted to the dimensions day α .
Since the quarantined class has better prevention than the infected one [45], we set the recovery rate to ϑ = 0.1 > η and the death rate to κ = 0.005. Here, we use (38,000, 2000, 1000, 3000, 50, and 3000) as initial values in the model (2). For these parameters value, we obtain R 0 = 1.1881 > 1, and thus two equilibrium points exist, namely, the unstable diseasefree equilibrium X 0 = (220,410,536.4, 0, 0, 0, 0) and the globally asymptotically stable endemic equilibrium X * = (184,244,619.17, 3,856,537.79, 2,198,361.42, 1,834,462.79, 26,769,323.88). These stability properties are verified via our simulation depicted in Figure 3. In Figure 3, there is different behavior for different orders of fractional derivative (α). It is clear that the smaller the order α, the slower the convergence of the solutions is. The peaks are also lower when α is smaller. Furthermore, for the numerical solution of R at α = 0.8, there is no peak in the time interval [0, 3500]. However, the curves of the solution for each α have the same behavior and preserve their convergence.  Table 2. For the second simulation, we execute the model using the parameter values in Table 2, except the quarantine rate is σ = 0.2. For the value of the parameters, we have R 0 = 0.81456 < 1, and therefore, only one equilibrium exists, that is, a globally asymptotically stable disease-free equilibrium X 0 = (220,410,536.4, 0, 0, 0, 0). This stability equilibrium is confirmed by our numerical simulation, depicted in Figure 4; that is, the disease-free equilibrium X 0 is globally asymptotically stable for all orders of fractional derivative (α). Again, the smaller the order of α, the slower the convergence of the solution. Therefore, endemic stability from the first simulation can be prevented by increasing the quarantine rate, which is actually supported by education and lockdown.
For the third simulation, the parameter values listed in Table 2 are used, except the recovery rate of infected individuals is η = 0.2. We have R 0 = 0.8075 < 1, and therefore, only one equilibrium point exists, that is, a globally asymptotically stable disease-free equilibrium X 0 = (220,410,536.4, 0, 0, 0, 0). This global asymptotic stability of the diseasefree equilibrium is confirmed by our numerical simulation, depicted in Figure 5. We observe that all numerical solutions are convergent to X 0 . Again, the order α affects the convergence rate of the solution. This simulation is the difference between our paper and those of Zeb et al. [21] and López and Rodó [24], both of whom did not consider the recovery rate of infected class η. We conclude that the recovery rate of the infected class is important to consider in the COVID-19 epidemic model because it can change the stability of the equilibrium point. Epidemiologically, the recovery rate can be increased by treatment based on Molnupiravir-assisted pharmacotherapy [46].

Parameter Estimation
In this section, we estimate the parameter values of model (2) using COVID-19 data from Indonesia. The COVID-19 data are adopted from the daily distribution map of confirmed active COVID-19 cases from the sixth wave, i.e., from 20 October 2022 to 10 January 2023, which is available on the official website of the Indonesian government at Distribution Map of COVID-19. Available online: https://covid19.go.id/ peta-sebaran (accessed on 15 January 2023) [47]. The number of confirmed active daily cases corresponds to the infected subpopulations of model (2). The model's parameter values are then estimated using the least square method, namely by minimizing the squared error between the numerical solution of the infected subpopulations and the data of confirmed active daily cases. To perform this process, we employ a builtin MATLAB function, that is, lsqcurvefit, which is based on the non-linear least square method. The initial value of model (2) is S(0) = 270 × 10 6 , E(0) = 50,000, I(0) = 18,919, Q(0) = 17,000, R(0) = 5 × 10 6 , M(0) = 150,000, where I(0) corresponds to first-day case of the data used. We also set the initial estimates for each parameter, that is, Λ = 270×10 6 70×365 , β = 10, µ = 0.00001, γ = 0.2, σ = 0.4, η = 0.9, δ = 0.001, ϑ = 0.9, and κ = 0.001. Each parameter value is estimated in the range from 0 to 1, except for recruitment rate Λ and disease transmission β. Ideally, recruitment rate Λ is around Λ = 270×10 6 70×365 = 10,567, where the total population of Indonesia is around 270 × 10 6 and the average age of the Indonesian population is around 70 years. Hence, we set Λ in the range of 0 to 20,000. We assume the estimated value of disease transmission β is in the range 0 to 50.
We perform parameter estimation for fractional-order α = 1, 0.97, 0.94, 0.9, 0.87, 0.84, 0.8. Notice that estimated values of the parameter for α = 1 correspond to the parameter of the first-order model proposed by Darti et al. [25]. For each value of α, model calibration is performed for 63 days of data, i.e., from 20 October 2022, to 21 December 2022, whereas the additional data from 22 December 2022, to 10 January 2023, is used for 20-day forecasting. We compare the calibration and forecasting of both the fractional-order model and the first-order model. To perform this comparison, we determine the root mean square error (RMSE), expressed by where N is the number of data points, I(t k ) is the numerical solution of the infected compartment of model (2) at t k , and I data (t k ) is the data of confirmed daily active cases at t k .
The estimated values of the parameters for different α are shown in Table 4. Notice that the estimated values of the parameter for each α are different. These cannot be compared between α because different alpha values have different units. The fitted curve for α = 1, 0.97, 0.93, 0.9, 0.87, 0.83, 0.8 is depicted in Figure 6. The confirmed daily active data used for calibration are marked by red dots, and the data for the 20-day forecasting are marked by blue dots. We discuss the results of calibration and forecasting from the viewpoint of the fitted curve and the obtained RMSE.
For the fitted curve of a first-order model, i.e., when α = 1, it can be seen that there are large deviations almost everywhere except for around the inflection point and peak point. Fitted curves for the fractional-order model typically have some deviation at the beginning and end of the period. Except at the end of the period, the fitted curves for α = 0.97 and α = 0.94 depicted in Figure 6 are good. It should be noted that not all orders α are good for data fitting. For example, the fitted curve for α = 0.9 appears to have a large deviation from the data used. Thus, the value of α also influences the data fitting performance. Nonetheless, based on the fitted curves, a fractional-order curve is typically better than a first-order curve for calibration and 20-day forecasting. These results show that a fractional-order model has greater flexibility for data fitting, and thus the best-fit curve for the data can be provided.  To confirm the observation of the fitted curves for different values of α, we also determine the RMSE for the calibration and 20-day forecasting of confirmed daily active cases given by Tables 5 and 6, respectively. It can be seen that the RMSE of calibration for the fractional-order model is smaller than the first-order model, except for α = 0.9. The smallest RMSE of calibration occurs when fractional order α = 0.97, where, in this case, the fractional-order model has a RMSE that is 2.7 times smaller than the first-order model. These results verify that a fractional-order model generally performs better than a first-order model in calibrating the COVID-19 data. Similarly, the RMSE of the 20-day forecasting for the fractional-order model is also smaller than the first-order model, except for α = 0.9. Again, a fractional-order model generally provides better 20-day forecasts than a firstorder model, especially the fractional-order model with α = 0.97, which has an RMSE that is 2.3 times smaller than the first-order model. Finally, we can say a fractional-order model has better accuracy and flexibility than a first-order model when implemented using confirmed daily active cases of COVID-19 in Indonesia. This also confirms the statement of Alrabaiah et al. [29] that a fractional-order model is more flexible than a first-order model.

Conclusions
We construct a fractional-order COVID-19 epidemic model with quarantine and standard incidence rate using the Caputo fractional-order derivative. We modify earlier models by simultaneously considering the recovery rate and quarantine rate of infected individuals. We also consider the standard incidence rate in the model. We have proven the existence, uniqueness, non-negativity, and boundedness of the solution. The model has two equilibrium points, i.e., a disease-free equilibrium and an endemic equilibrium. We obtain the basic reproduction number (R 0 ) by implementing the spectral radius of the nextgeneration matrix. The disease-free equilibrium always exists and is locally and globally asymptotically stable only if R 0 < 1. In contrast, the endemic equilibrium exists and is globally asymptotically stable if R 0 > 1. The stability properties of the equilibrium points are confirmed by numerical simulation. From the numerical simulation, we can see that the smaller the order of the derivative, the slower the convergence of the solution of the model and vice versa. In addition, we conclude that both the recovery rate and quarantine rate of the infected class have an important role in reducing R 0 , as well as for the stability of the equilibrium point. Based on parameter estimation results, a fractional-order model typically performs better than a first-order model for both calibrating and forecasting confirmed cases of COVID-19 in Indonesia. In our case, the fractional-order model with α = 0.97 has an RMSE that is twice as small as the first-order model for both the calibration and forecasting of confirmed COVID-19 cases.

Data Availability Statement:
The data used in this study are available on the official website of the Indonesian government at https://covid19.go.id/peta-sebaran, accessed on 15 January 2023.