Modeling and analysis of the dynamics of novel coronavirus (COVID-19) with Caputo fractional derivative

The new emerged infectious disease that is known the coronavirus disease (COVID-19), which is a high contagious viral infection that started in December 2019 in China city Wuhan and spread very fast to the rest of the world. This infection caused million of infected cases globally and still pose an alarming situation for human lives. Pakistan in Asian countries is considered the third country with higher number of cases of coronavirus with more than 200,000. Recently, many mathematical models have been considered to better understand the coronavirus infection. Most of these models are based on classical integer-order derivative which can not capture the fading memory and crossover behavior found in many biological phenomena. Therefore, we study the coronavirus disease in this paper by exploring the dynamics of COVID-19 infection using the non-integer Caputo derivative. In the absence of vaccine or therapy, the role of non-pharmaceutical interventions (NPIs) is examined on the dynamics of theCOVID-19 outbreak in Pakistan. First, we construct the model in integer sense and then apply the fractional operator to have a generalized model. The generalized model is then used to present the detailed theoretical results. We investigate the stability of the model for the case of fractional model using a nonlinear fractional Lyapunov function of Goh-Voltera type. Furthermore, we estimate the values of parameters with the help of least square curve fitting tool for the COVID-19 data recorded in Pakistan since March 1 till June 30, 2020 and show that our considered model give an accurate prediction to the real COVID-19 statistical cases. Finally, numerical simulations are presented using estimated parameters for various values of the fractional order of the Caputo derivative. From the simulation results it is found that the fractional order provides more insights about the disease dynamics.


Introduction
A highly contagious viral disease caused by a coronavirus (SARS-CoV-2), first identified in China at the end of 2019 and affecting 213 countries and territories around the world [22]. The novel Wuhan virus figures under different names such as the novel coronavirus (2019-nCoV) on January 7 and COVID-19 by the World Health Organization (WHO) on February 11 [38]. The transmission of the virus was considered initially animal-to-human and fist human-to-human transmission confirmed on January 20 in Guangdong, a coastal province of southeast China [38,12]. Over 10.27 million confirmed cases and more than 0.5 million cumulative deaths all over the world are recorded till June 30, 2020 [1]. The mortality rate is about 7% and the recovery rate is 93% globally [1,22]. The pharmaceutical companies, regulatory authorities are doing every thing possible for the availability of a safe and effective anti-viral vaccine against this novel infection. The pandemic continuously inflects severe public health, hit the laborer's community and economy throughout the world, including in Pakistan. The suggested incubation period for this viral infection is in the range 2 to 10 days by WHO [25,12]. The symptoms of a COVID-19 infected person include fever, myalgia or fatigue, dry cough, shortness of breath or dyspnoea, pneumonia, chills, sore throat, multiple organ failure, anorexia or lungs failure, acute respiratory distress syndrome (ARDS), lymphopenia (a reduction of lymphocytes in the circulating blood), loss of smell or taste and a sign of cytokine storm [25,12]. The majority of deaths have occurred in other chronic illnesses like diabetes, hypertension and cardiovascular disease. In the absence of therapies, approved vaccine and antiviral the whole world is focusing on non-pharmaceutical interventions to eradicate this COVID-19. The most common NPIs practices these days are social physical distance, exposed to be a quarantine, isolation of infected, wearing mask, more testing facility, closure of the school, non-essential business, hand washing or sanitizing, lockdown (stay at home and work from home) and protective kit for medical personal and avoiding unnecessary gathering.
Pakistan is one of the epicenters of COVID-19 in Asia with an estimated population of over 220 million [1]. Currently, Pakistan is the 3rd country in Asia and 12th in the world with high confirmed infected cases [17,1]. As on 30 June 2020, 209,337 total confirmed cases are recorded and about 4304 lost their lives. Although in Pakistan the mortality rate is very low as compared to the recovery rate, still this disease has a negative impact on the economy. The first case reported in Karachi on 26 Feb 2020, the virus spread quickly within three weeks in the entire country. Also, the accounted cases are less than the total cases due to limiting testing. Currently, the situation is worse, due to open lockdown on 9 May, the government is unable to maintain strict lockdown because of severe economic hardship. But the government imposed smart lockdown and easing restrictions, by allowing offices, business, markets with limited hours and staff capacity, five days in a week to reopen the economy of the country.
Mathematical models coupled with statistical data and its analysis are better to obtained reasonable information about the coronavirus disease and its complications and future spread and control. Recently, it is a huge challenge that without vaccines, how to curtail the pandemic with non-pharmaceutical intervention (NPIs). In this regard, different mathematical modeling approaches have been adopted to analyze the transmission pattern of COVID-19. Ferguson et al. [16] studied a COVID-19-induced mortality model with the impact of NPIs. A mathematical model for accessing the community-wide impact of mask in controlling the spread of coronavirus is developed by Eikenberry et al. [14]. The dynamics and mitigation strategies for COVID-19 in Canada have been studied in [39]. A compartmental system describing the coronavirus dynamics in three highly infected countries has been suggested in [15]. Recently, Ullah and Khan studied the coronavirus model for Pakistani data using optimal control theory [32]. Fractional calculus is the generalization of classical calculus. Mathematical models with noninteger order operators provide a better understanding of a phenomena. Further, models with fractional order derivatives are capable to capture the fading memory and crossover behavior and provides a greater degree of accuracy. Mathematical models with fractional derivative gives more insights about a disease under consideration [11,26,29,30]. Different fractional operators with singular and nonsingular kernel were suggested in literature [5,27]. The applications of these fractional operators can be found in recent literature and references therein [3,8,9,20,[33][34][35]. Recently, very few COVID-19 models based on fractional order operators are proposed. The dynamics of (2019-nCoV) fractional derivative to explore the situation of the pandemic in Wuhan is studied by Khan and Atangana [19]. Some other work about the mathematical models and its applications to fractional operators can be seen in [28,36,21,6,7,23,18].The authors in [28] considered a mathematical model for hepatitis C in fractional derivative and examined its results. The CoVID-19 model and its dynamical analysis using fuzzy Caputo differential equation is considered in [36]. The TB dynamics with relapse has been considered in [21] where the authors explored the results using conformable derivative. Some interesting results regarding fractional calculus and its connection to attract many physical structure is studied in [6]. The authors in [18] obtained the stability and numerical results for a HIV/AIDS model in fractional derivative. The oxygen diffusion problem in various fractional operators is discussed in [23]. A novel fractional-fractal operator with the singular and nonsingular kernel is applied to obtain a dynamical model for coronavirus is studied in [4]. Baleanu et al. [10] examined the dynamical investigations of coronaries model in fractional derivative using exponential kernel. A fractional order COVID-19 model in Caputo sense is recently studied in [31]. A more recent work on coronavirus considering the infected cases observed in Kingdom of Saudi Arabia is investigated in [2].
Motivated by the above discussion, currently in this work, we study the dynamical analysis of the COVID-19 model presented in [32], considering the Caputo operator in order to gain more insights about the pandemic. We consider here the real infected cases of novel coronavirus in Pakistan and obtain the results. A detailed theoretical analysis of the fractional order COVID-19 model is presented. Further, we estimate model parameters and the basic reproduction number using confirmed COVID-19 cases reported from the beginning of the outbreak till June 30, 2020 in Pakistan. The impact of important model parameters is shown for various values of the arbitrary fractional order of the Caputo operator. The rest of paper organization is as follows: The basics concepts related to fractional order derivatives are presented in Section 2. The model formulation in integer with parameter estimation and curve fitting to reported cases and its generalization to the fractional order derivative are presented in Section 3. Section 4 contains the analysis of the model and its stability results. In Section 5, we determine the important parameters and there impact on the basic reproduction number through sensitivity analysis is presented. The graphical results through biological discussion is depicted in Section 6. In Section 7, the work has been summarized with important suggestions about the disease control in the society.

Preliminaries on fractional derivative
In this section, we briefly discuss background material from fractional calculus.

Definition 1.
Consider u ∈ C n be a function then the Caputo definition with α in (r − 1, r) where r ∈ N [27] is defined as: obviously, C D α t (u(t)) approaches to u ′ (t) whenever α→1.

Definition 4.
The associated fractional integral of the ABC derivative is defined as; where α ∈ [0, 1).

Definition 5.
[37] Consider a constant point for the Caputo system say y * which is known to be its equilibrium point then, if and only if u(t, y * ) = 0.

The classical integer order model
In this section, we briefly present description of the COVID-19 model in integer case. The model was recently proposed in [32]. To formulate the mathematical model, the whole population at time t, denoted by N(t) is further split into eight mutually-exclusive sub-classes including susceptible S(t), exposed E(t), those infected people who infected after some specific incubation period with symptoms (or symptomatic infected) I(t), those have no specific disease symptoms (or asymptotically infected individuals) I a (t), quarantined class Q(t), individuals which are in hospitals or in self-isolation at home I h (t), criticallyinfected or COVID-19 infected people who were in ICU I c (t) and finally, the individuals successfully recovered from the infection is denoted by R(t). In the class I a (t) infected individual showing mild or no symptoms and the class I h (t) contains the patient admitted into the hospital and those who are self-isolating at home with clinical symptoms. The state variables of the model are described in Table 1, while the complete description of the model parameters are describe in Table 2.
The flow of among model compartments is depicted in Fig. 1. The COVID-19 transmission model rely on these subgroups and at the rate they interact, governed the non-linear system of ordinary differential equations as follows: subjected to initial conditions

Parameter estimations
The purpose of this section is to estimate the parameters with help of the least square curve fitting from the confirmed cases of coronavirus that notified in Pakistan since March 1 till 30 June 2020. The recruitment as well as the natural death rate respectively given by Λ and μ are estimated from literature and the rest of the parameter values are fitted from the reported data. The detailed estimation procedure can be found in [32]. The best fitted model predicted curve is shown graphically in Fig. 2 and the model parameters with resulting estimated or fitted numerical values are shown in the Table 2. The numerical value of the R 0 obtained via the fitted parameters is R 0 ≈ 1.8870, which is in agreement to the value estimated in recently published article [32].

Model in Caputo sense
In this subsection, we formulate the fractional order COVID-19 mathematical model. In order to observe the memory effects, the system (1), in term of integral as: where, k(t − t ′ ) act as time-dependent kernel and by power-law correlation function as: Substituting the value of kernel in Eq. (2), and then applying the Caputo fractional derivative of order α − 1, then we have; Both are the inverse operators then we will get the system (4), and Then the novel coronavirus COVID-19 model in Caputo derivative form is; subjected to the non-negative initial conditions.
In the above Caputo system (6), the model developed here has eight components, S, V, E, I, I a , Q, I h , I c and R known as COVID-19 model of disease transmission also called fractional differential equations (FDEs) of eight-dimensional dynamical system.

Properties of the model
We provide in details in the following the model important properties.

Invariant region and attractivity
The dynamics of the fractional model (6) will be analyzed in the following biologically-feasible region, Lemma 3.1. The region given by Ξ⊂R 8 + , is positively invariant for the system (6) with the initial conditions in R 8 + .
Proof. We have the following for the model (6), Hence, clearly we have, By applying the Laplace transform, Applying Laplace inverse, we obtain, The Mittag-Leffler function describe by infinite power series i.e; and laplace transform of Mittag-leffler function is So N(t) converges for t→∞ and for all t > 0 the solutions of the model with initial conditions in Ξ remain in Ξ. Thus, the region Ξ is positively invariant and attracts all solutions in R 8 + .
To show that the model has positive solution, we consider Corollary 1. [24] Suppose that

Positivity and boundedness
. Proof. In order to explore the solution non-negativity, it is require to prove that on each hyperplane bounding the positive orthant, the vector field point R 8 + . Form the system (6) we obtained: Thus, by using the above corollary, the desire goal has been obtained i.e, the solution will stay in R 8 + and so, we can have the following region that biologically feasible:

Disease-free equilibrium DFE
For the equilibrium points of the fractional model (6), we have

The basic reproduction number
To derive the important threshold parameter also known is the basic reproduction number R 0 , by using the next-generation technique, for the dynamics of the disease, whether the infection in population is decaying or growing. Let x = (E, I, I a , Q, I h , I c ) T , then we have The Jacobian of above matrices for the linearized system at disease free-state is: Hence the next-generation matrix is obtained as The basic reproduction number R 0 at DFE is obtained by taking the spectral radius of the Next-generation matrix ρ(FV − 1 ) for fractional model as follow:

Local stability of DFE
.
Proof. For local stability, the linearized system is the jacobian of system (6) at E ** 0 implies that: After the evaluation of determinant Eq. (10) implies, here we can say that the three eigen values are − μ, − μ, − K 6 have negative or negative real part and the co-efficient mentioned above is of the form Clearly, a i are all positive if R 0 < 1. The argument of the root of equation (λ m + μ) 2 = 0 and (λ m +K 6 ) = 0 are similar, that is: , where k = 0, 1, 2, …., (m − 1).
In similar fashion, we find out the arguments of the roots of equation (λ 5m +a 1 λ 4m +a 2 λ 3m +a 3 λ 2m +a 4 λ m +a 5 ) = 0 are all greater than π 2M if R 0 < 1, having an argument less than π 2M for R 0 > 1. Thus the DFE is locally asymptotically stable for R 0 < 1.

Global stability of DFE
To show global asymptotic stability (GAS) of DFE E ** 0 of the model (6) by Layapunov function approach, we proceed as in the following result.

Theorem 4.2. If R 0 ⩽1 then DFE of the Caputo COVID-19 model given by (8) is GAS.
Proof. Let the Lyapunov function be in the form below: The Caputo fractional derivative of F(E, I, I a , Q, I h ), along model (6) satisfies: Hence, it follows that C D α t F ≤ 0, for R 0 < 1 all the parameters and variables are non-negative with C D α t F = 0 iff E = I = I a = Q = I h = 0. So the infected compartments (E, I, I a , Q, I h ) approached zero whenever t→∞. Using E = I = I a = Q = I h = 0 in equations of system (6), we have S→ Λ μ ,I c →0 and R→0 as t→∞. Thus using lypaunov stability theorems for fractional case developed in [37], every solution of the Caputo model (6) with non-negative initial conditions, approaches to E ** 0 as t→∞ in Ξ. Thus, it follows that the DFE of the model (6) is GAS.
In the next subsections, we have to explore the positive equilibria of the model where at least one of the infected component is non-zero.

Existence of endemic equilibrium point
The arbitrary endemic equilibrium point (EEP) for the Caputo model (6), represented by with N ** = (S ** + E ** + I ** + I ** a + Q ** + I ** h + I ** c + R ** ). Now at steady state, after solving the equations of Caputo model gives Further, we define the force of infection at endemic state as below: substituting (12) in (13), we have after simplification, we get Hence, where,

Lemma 4.2.
There exists a unique endemic state for the model (6), whenever R 0 > 1.

Local stability of (EEP)
To formulate the Jacobian of system (6) at endemic state. The following result is obtained as: Theorem 4.3. The endemic equilibrium E ** 1 of Caputo-fractional model (6) Proof. The Jacobian of system (6) at associated endemic equilibrium (EE) follows that; where, One of the eigenvalue is − μ, that is negative and others can be obtained from the equation: Now the Hurwitz matrices are: , Here we can say that the coefficient signs b j are strictly positive for the roots of the polynomial to be negative. All roots of polynomial have negative real part iff, DetH j > 0, j = 1, 2, …., 7. Thus, the conditions ensures the endemic stability.

Global stability of EEP
For global stability at endemic state of the Caputo COVID-19 model (6) is given for the special case when infected hospitalized do not transmit infection (τ = 0), and disease-induced mortality is negligible (ξ 1 = ξ 2 = 0), then the reduced model is obtained as: where, now The reduced model shows that N(t)→ Λ μ as t→∞. Finally we have The relative expression for the reproduction number R 01 of the reduced model with (16) is The following results for model (15) at steady state obtained as: K 1 E *** = λ * 1 S *** , k 2 I *** = θωE *** , Proof. Consider the following non-linear Lyapunov function for the model (15), so that the corresponding unique EEP E * * * 1 exist by letting R 01 > 1.
The time fractional derivative of (19) is: Using the solution from (18) gives (   ( . ,  Substitute (20) to (23) in (20) Since the arithmetic mean exceeds the geometric mean, we have the following interpretation form (24): Thus every solution of the reduced model, tends to its unique endemic equilibrium for associated reproduction number as t→∞. For the GAS of EEP (E ** 1 ) whenever R 0 > 1, through the conjecture.

R 0 versus model parameters
Figs. 3-7 depict the impact of model various important parameters on the value of R 0 with the corresponding contour plots. The variation in the value of R 0 versus the effective contact rate β and quarantine rate δ is shown in Fig. 3. It is observed that R 0 decreases for values smaller than 1 with a decrease in β and an increase in δ. Similarly, R 0 can be reduced to a smaller value less than 1 with the increase in hospitalization/self-isolation rate γ and decrease in the effective contact rate β. This interpretation can be found in Fig. 4. The combined effect of the infectious rate due to hospitalized infected patients ψ and δ is analyzed in Fig. 5 while the behavior of R 0 for the variation in τ and δ is depicted in Fig. 6. This interpretation demonstrates that R 0 can be effectively reduced with a reduction in the effective infectious rates (τ and δ) and enhancing the efficacy of contact-tracing and quarantine policies of the exposed population. Finally, we analyzed the role of variation of both quarantine rate γ and hospitalization rate γ on the value of R 0 . It is worth mentioning that R 0 remains very smaller with an increase in both these parameters.

Numerical results and discussion
This section presents the numerical results of the Caputo COVID-19 model (6). The generalized predictor-corrector of Adams-Bashforth Moulton method developed in [13,11] is considered for the numerical solution of the model (6). To simulate the model, we use the real parameters fitted given in Table ( 2). To study the model qualitative behavior and its parameters effect on the system and possible control, we varied the model important parameters and the fractional order α and obtained the simulation results. We studied graphically the coro- It is observed that the susceptible population is increased while the population in all remaining classes is decreased significantly, when the fractional order α takes smaller values. The impact of effective contact rate β on the dynamics of symptomatic and asymptomatic infected individuals is shown in Fig. 10 for four different values of fractional order α of the Caputo derivative. A significant decrease in newly reported symptomatic and asymptomatic infected individuals is seen with a decrease in β. Further, with a 40% decrease in effective contacts (i.e., β = 0.4683) to its estimated baseline value, a dramatic decrease is observed in the peak of infection curves. This trend is observed for all values of α, but for smaller values of fractional order, the decrease in the infection curves is comparatively more significant. Fig. 11 illustrates the influence of effective contacts β on the behavior of hospitalized and critically infected (or ICU patients). It is seen that the total hospitalized and the critically infected population is decreased very well when β is decreased at 10%, 20%, 30% and 40% respectively to its baseline value as shown in Fig. 11. Moreover, the decrease in the cumulative hospitalized and critically infected individuals becomes more significant for smaller values of fractional order α. As a consequence of this graphical interpretation, the severity of COVID-19 infection can be reduced by following a strict Standard Operating Procedure (SOPs) in order to reduce the effective contacts among infected and susceptible individuals. The impact of quarantine rate or contact-tracing parameter δ on the dynamics of symptomatic and asymptomatic infected individuals is analyzed in Fig. 12. From this analysis one can observe that with an increase in δ (up to 50 % i.e. δ = 0.5574) to the baseline value a sharp decrease in the peak of infected curves is seen as shown in Fig. 12. We depicted the simulations for four distinct values of α and obtained identical behavior for fractional order parameter values, that the decrease in the infected individuals is comparatively faster for smaller values of α. Furthermore, Fig. 13 depicts the dynamics of the cumulative hospitalized and critically infected population for variation in the δ. With the increase in the contact-tracing policy, a decrease in hospitalized and the critically infected population is observed. The same A small variation in the proportion of symptomatic individuals occurs even if we increase γ at a significant rate. It is found that when α decreases, the variation in the infected population is slightly feasible.

Conclusion
Fractional order derivatives are more prominent and provide comparatively deeper insights about real-world problems. Moreover, the derivative with fractional order is capable to capture the fading memory and crossover behavior found in most biological phenomena including infectious diseases. In this paper, we studied the dynamics of the COVID-19 outbreak in Pakistan using a fractional order mathemat-ical model in Caputo sense which is the most commonly used operator in the modeling approach. Initially, the Caputo fractional COVID-19 model is formulated and then provided its basic and necessary mathematical properties that include the existence and positivity of the system and its solution. The detailed stability results both local and global for DFE and EE are explored using fractional stability approaches. The impact of variation in various model important parameters on the dynamics of the basic reproduction number is depicted graphically. Furthermore, the model parameters and the basic reproduction number R 0 are estimated from the COVID-19 real infective reported cases from health ministry of Pakistan from the start of the outbreak till June 30, 2020. The non-linear least square procedure is used for the estimation purpose. The present finding showed that the results of predictions obtained through our model is inline with those published in literature and thus more appropriate values of the parameters are estimated. The numerical value of the basic reproduction number evaluated using the parameters obtained from fitting is R 0 ≈ 1.8870 which is inline to the published result in [32]. Finally, we simulated the Caputo model in order to examine the impact of various NPIs on the dynamics and control of COVID-19 in Pakistan. We varied the relevant parameters at different levels to its baseline values and obtained the graphical results for distinct values of fractional order α of the Caputo derivative. It is found that the reduction in the effective contacts β (up to 40 %) and an enhancement in contact-tracing policy to quarantine the exposed individuals δ (up to 50%) to its baseline value that dramatically reduced the peak of infected curves. Moreover, the reduction in the infected population becomes more significant for smaller values of fractional order α. Thus we conclude that the results obtained for the fractional case are reliable, realistic and more biologically feasible. In future, the present study can be extended to fractional order with the non-singular kernel.