Modeling and simulation of the novel coronavirus in Caputo derivative.

The Coronavirus disease or COVID-19 is an infectious disease caused by a newly discovered coronavirus. The COVID-19 pandemic is an inciting panic for human health and economy as there is no vaccine or effective treatment so far. Different mathematical modeling approaches have been suggested to analyze the transmission patterns of this novel infection. this paper, we investigate the dynamics of COVID-19 using the classical Caputo fractional derivative. Initially, we formulate the mathematical model and then explore some the basic and necessary analysis including the stability results of the model for the case when R 0 < 1 . Despite the basic analysis, we consider the real cases of coronavirus in China from January 11, 2020 to April 9, 2020 and estimated the basic reproduction number as R 0 ≈ 4.95 . The present findings show that the reported data is accurately fit the proposed model and consequently, we obtain more realistic and suitable parameters. Finally, the fractional model is solved numerically using a numerical approach and depicts many graphical results for the fractional order of Caputo operator. Furthermore, some key parameters and their impact on the disease dynamics are shown graphically.

Introduction

Coronavirus disease (COVID-19) is a highly contagious viral disease that affects million of humans and caused reasonable death cases around the world. The infection was named by the World Health Organization (WHO) as COVID-19. The infection for the first time was identified in Wuhan China, which is the capital city of Hubei Province in Central China. Later, on 11th March 2020, WHO announced the COVID-19 disease as a pandemic spreading around the globe. According to the recent WHO statistics reported on 11 May 2020, there are 3976043 confirmed COVID-19 cases and it has taken the lives of more than 277708 people in the whole world [1]. About 215 countries, areas or territories are targeted by this novel infection so far and it is still posing a serious threat to the other countries on the globe. There are no specific vaccines or treatment for COVID-19 up till now and most of the patients have recovered without specific treatment or hospitalization. This infection affected people indiscriminately of all ages, nonetheless, people suffering from co-morbidities like diabetes, cardiovascular disease, cancer and chronic respiratory disease and with older age are at higher risk to develop serious illness leading to death.

This infection is extremely contagious, spreading from an infected person through droplets discharge from the nose while spitting, coughing, or sneezing. Although, the coronavirus affected people in different ways but the most common symptoms of this infection include dry cough, fever, aches and pains, headache and excessive fatigue, etc., while some of the serious symptoms are considered to the loss of speech or movement, chest pain or pressure, and the breath shortness. Normally, it takes five to six days from initial infection to show symptoms however, it can take up to fourteen days [1], [2]. The people once infected with this virus should transmit the virus to other people by shaking hands, so, social distances, and restricted itself at home, and using masks and gloves during in open exposure when desired, should be followed to get safe.

Due to the severity and life-threatening situation of the COVID-19 pandemic on health and society, different researchers, medical doctors, scientists and other health departments are adopting different approaches in order to explore the transmission patterns and possible control of this infection. Despite much scientific improvements, the WHO stresses that an urgent understating of the epidemiology and the evolution of the COVID-19 the outbreak is needed for the control strategies to stop the transmission. In this regards, one of the essential tools is the infectious mathematical modeling approach, which has adopted in recent days in order to explore the insights into the transmission, severity, and particularity of the disease, which could readily inform decision-making concerns in mitigating this infection [3], [4]. In the last couple of months, the researchers around the world of different field published their research related to coronavirus. Regarding the mathematical models for the coronavirus, many researchers developed mathematical models to understand the transmission dynamics of coronavirus infection in different regions. In [5], Li et al., an epidemic COVID-19 model with the latency period has been developed and the model has been fitted to the incidence data reported in the mainland of China. In [6], the authors formulated SIR based model for the  COVID-19 infection and implemented the proposed model to explore and predict the transmission dynamics of this pandemic for the countries with high a number of cases such as France, China and Italy. In [7], Ribeiro et al. employed regression models to predict COVID-19 infected cases in Brazil. While in [8], Ndairou et al. developed a transmission model with super-spreader infected class and implemented their model to the reported infected case of Wuhan. They have shown that the suggested mathematical model provides a better agreement to the reported COVID-19 cases and deaths that have accounted due to coronavirus infection. Recently, the role of non-pharmaceutical intervention on COVID-19 dynamics by considering the data of Pakistan has been studied in [9]. The authors used the real data of coronavirus of Pakistan and implemented a mathematical model with optimal control analysis.

The fractional epidemic model is another strong tool to explore the dynamics of infectious disease in a better way than the ordinary integer-order models. These models are based on fractional order differential operators, which generalizes the classical integer-order derivatives. Epidemic models with fractional derivatives give a greater degree of accuracy and provide a better fit to the real data in compression with classical integer-order models [10], [11]. A diverse variety of fractional operators were introduced time to time in the literature with a different kernel. Some of the frequently  used fractional orders (FO) operators are Caputo [12], Caputo-Fabrizio (CF) [13] and Atangana-Baleanu operator (AB) [14]. Although most of the COVID-19 models developed so far are based on classical integer-order derivative but only a few can be found with fractional operators. For example, in [15] Khan and Atangana formulated a fractional model of COVID-19 using Atangana-Bleanu in Caputo sense (ABC) operator and provided a better fit to the reported case in Wuhan. Recently, in [16] Abdo et al. extended a classical COVID-19 model in the literature to fractional order using Mittag–Leffler kernel and carried out the Ulam stability analysis and simulation results. Atangana studied a new COVID-19 transmission model using different fractal-fractional operators [17]. The results are obtained both mathematically and statistically considering the real cases. Further, the lockdown effect in the model has been analyzed. The authors in [18] formulated a mathematical model for coronavirus from natural to human host. Some recent mathematical studied on corona and other infectious diseases are studied by authors see [19], [20], [21], [22]. A novel results for the corona virus model has been discussed in [19]. The application of fractal-fractional differential equations to partial differential equations has been analyzed in [20]. The authors in [21] presented a theoretical results for fractal fractional differential equations. Fractal fractional differential equations and its solutions in detailed with novel analysis is presented in [22].

Inspired by the above discussion, in the present study, we reformulate the model [15] using the fractional operator in the Caputo sense. In [15], they considered the reported cases for a short period of time from January 21, 2020 to January 28, 2020, in Wuhan city China. In this study, we consider the reported case in the mainland of China from January 21, 2020 to April 9, 2020 and explore a better theoretical and graphical analysis of the model. The data is taken from [23]. Using the reported data we will obtain the real parameter values for the model and will perform the simulation to predict the dynamics of coronavirus in the community. We will show that the model provides better fitting to the reported cases. Some important parameters that can help the disease under control will be explored in detail. The rest of the results are organized is as follows: The concepts and basic results regarding fractional and Caputo fractional derivative are recalled in Section 2. Formulation of the model, the integer case, the model, and their fitting to the coronavirus cases and parameters are discussed in Section 3. The COVID-19 model with Caputo fractional derivative, its basic mathematical properties as well as stability results are given in Section 4. The numerical iterative scheme and the simulation results of the fractional model are plotted and discussed in Section 5. Finally, a brief concluding remarks are given in Section 6.

Basics of the fractional calculus

In the following, we first recall some basic details of fractional calculus [12].

Let be a function then the fractional derivative in Caputo sense having order where is given as:

Clearly, approaches to as .

A fixed point say, satisfying , then it is called the equilibrium point of the model formulated via Caputo derivative given by:

As proved in [24], the point will be locally asymptotically stable if all eigenvalues of the Jacobian matrix calculated around the equilibrium point satisfy the following condition:

For the global stability results of the fractional system in Caputo derivative considering the Lyapunov method, we need to have the following important results [25], [26]. and and . Where and denote the continuous positive definite functions on . Then, will be uniformly asymptotically stable equilibrium point of the (1).

Suppose be an equilibrium point of the system (1) and be a domain containing , and let be a continuously differentiable function satisfying

Mathematical model

In this section, we briefly explain the model formulation using classical integer-order derivative. To develop the model total human population denoted by is further classified into five mutually-exclusive epidemiological classes, which are the susceptible , exposed , infected (symptomatic) , asymptotically infected showing no clinical symptoms , and the recovered people . The class denotes the COVID-19 in reservoir or the seafood place or market. The recruitment rate is shown by the while the natural death rate is represented by . The symptomatic and asymptomatic infected people could export the virus into M at the rate and , respectively. The virus in M will subsequently leave the M class at a rate of , where accounts the lifetime of the COVID-19 virus. The susceptible individuals acquired infection after effective contacts with the people in I and A compartments at the rate , where the parameter is the disease transmission coefficient and is the transmissibility multiple of A to I and . Further, we assume that susceptible people could be infected after the interaction with M, given by where, is the transmission rate from M to S. The exposed people develop the COVID-19 infection after the completion of incubation period and then move to either class I or A at the transmission rate and respectively. The parameter is the proportion of asymptomatic infection. The infectious period of symptomatic I and asymptomatic A individuals are defined as and respectively after which they enter to recovered class R. The nonlinear differential equations that describes the dynamics of the COVID-19 disease are given by [15]:

The corresponding initial conditions are

For simplicity, let us denote

Then, the above model can be written as

Next, we estimate the parameters of the model and provide the method with brief explanation.

Data fitting and estimation of the parameters

We investigating the parameters estimations and curve fitting for the reported data of coronavirus cases from January 21, 2020 to April 9, 2020, in the mainland of China to the model (5). We utilized the approach first by selecting some of available parameters from the literature and then the rest of the parameters are fitted to the coronavirus cases. We consider the time unit per day. In the following, we provide a brief explanations of the parameters estimation:

(i) Natural mortality rate : The life span an average in China for the year 2019 is years [27], so, the estimated value of natural death rate is per year.

(ii) Recruitment rate : Population of China for the year 2019 is about =1.43 billion [27], so, we can obtain the recruitment rate parameter by computing , and also under the assumption that in the disease absence it is the limiting population, so we have per day.

Using the above parameters values and the rest of the parameters are simulated and fitted to the actual cases that is depicted in Fig. 1 using the method of least-square curve fitting briefly discussed and utilized in [28], [29]. The best fitted model parameters are calculated by minimizing the error between the reported infected data and the infected cases obtained from the model (5). The objective function implemented in the estimation procedure is as follows:where denotes the reported cumulative COVID-19 confirmed case and shows the model solution for infected class at time , while m is the number of available reported data points. The cumulative reported infected cases in the mainland of China are shown in Fig. 1 and the resulting best fit by the proposed model to field data is shown in Fig. 2 . One can observe that the model provide very accurate fitting to the reported cases of coronavirus and considered to be more reasonable in order to do prediction based on the fitted parameters. Table 1 contains the fitted and the estimated parameters and the approximated value of the basic reproduction number using these parameter values is .

Fig. 1

Cumulative confirmed COVID-19 cases notified in China from January 21, 2020 to April 9, 2020.

Fig. 2

Model fitting (blue solid curve) to the cumulative reported COVID-19 confirmed cases using model (5).

Table 1

Parameters values and their descriptions.

Parameter	Description	Value	Source
Λ	The recruitment rate	μ×N(0)	Estimated
η	Contact rate	0.134	Fitted
μ	Natural mortality rate	1/(76.79×365)	[27]
ψ	Transmissibility multiple	0.0002	Fitted
η_w	Disease transmission coefficient	0.0000000080082	Fitted
θ	The proportion of asymptomatic infection	0.41003	Fitted
ω	Incubation period	0.0000213	Fitted
ρ	Incubation period	0.480322	Fitted
τ	Rate of recovery of I	0.000033	Fitted
τ_a	Rate of recovery of A	0.59	Fitted
ϱ	Transmission of the virus to M by I	0.0101	Fitted
ϖ	Transmission of the virus to M by A	0.1314	Fitted
ν	Removal rate of virus from M	0.5	Fitted

A Caputo fractional model

We extend the classical integer model by replacing the derivative with the Caputo fractional-order derivative in order to provide a better understanding of the COVID-19 infection. The parameter in the model is used for the order of Caputo derivative. Hence, the model is given through the following system of nonlinear fractional differential equations:with the initial conditions given in (6).

Existence and positivity of the Caputo model solution

To proceed the model analysis, first we explore the basic necessary mathematical properties. For the non-negativity of the model solution, let us construct the following set

We follow the approach used in [30] in order to prove the desired results and provide the following theorem.

[30] Consider, and , then with .

[30] Suppose that and , where , then the following conclusions are drawn if

For the fractional model (9) there exists a solution, say, which is unique and will remain in and non-negative.

Utilizing the approach presented in [31], it is easy to show the existence of the Caputo model (9). Further, the uniqueness of the solution can be obtained by making use of the Remark 3.2 in [31] for all values of . For the positivity of the solution, it is necessary to show that on each hyperplane bounding the positive orthant, the vector field points to . From the system (9), we deduced that

Hence, based on the above corollary, we concluded that the solution will remain in , for all .

Finally, to confirm the boundedness of the fractional model, we proceed by summing the first five equations of the model (9), which give

Applying the Laplace transform of (10), we obtainedApplying Laplace inverse, we obtained

The Mittag–Leffler function describe by infinite power series i.e;and laplace transform of Mittag–leffler function is Thus keeping the fact in mind that the Mittag–Leffler function has an asymptotic behavior [12], we obtain

Hence, the biologically feasible region is constructed as:

Next, we investigate the equilibria and the associated basic reproduction number .

Model equilibria and basic reproduction number

The equilibrium points of the Caputo fractional model (9) can be evaluated by setting

The disease free equilibrium (DFE) denoted by is obtained by solving the system (12) at disease free state

Further, using next generation approach we evaluate the following expression of the basic reproduction number:

Existence of endemic equilibrium

The endemic equilibrium (EE) of the Caputo model (9) is represented by and is given as follows:andsatisfying the following equation

The coefficients in (13) are as follows:

Clearly, and if , and . Thus, no endemic equilibrium will exist if .

Stability of the DFE

In order to confirm the local stability of DFE it is enough to show that the eigenvalues of the matrix of dynamics given (16) lie outside the closed angular sector as proved in [24].

The following theorem deals with the desired result.

DFE of the model is locally asymptotically stable ( LAS ) if for all eigenvalues of the Jacobian matrix of the model (9) satisfy the following condition :

The Jacobian of linearization matrix of model is given as:

The characteristic equation in term of for is given below:where,

From (17) the argument of the first two eigenvalues , satisfy the condition given in (15) for all . Further, it is clear that if , then all , and it is easy to show that . Thus, following [24], [25], the arguments of all eigenvalues satisfy the necessary condition given in (15). This completes the proof.

Global asymptotic stability of the DFE

In this part, we prove the global asymptotical stability (GAS) of the fractional model (9). For this purpose, the following result is presented.

For the DFE of the proposed model (9) is GAS on if 1.

To present the proof, let us define the following suitable Lyapunov function:where , for j= , which are positive constants to be consider later. Taking the time Caputo fractional derivative of we obtain

Using system (9), we get

Let us choose the constants values asand then after some simplification, we have,

It is clear that when then is negative, because all the model parameters are non-negative. Thus, it follows from the results given in Theorem 0.1, that the DFE is GAS in .

Numerical solution of fractional model

This section investigates the numerical solution of the fractional model (9) and to present the simulation results for various values of model parameters and . In order to do this, we utilize the Euler’s type approach presented discussed in the recent literature and references therein [32], [33]. To obtain the iterative scheme, let us express the fractional model (9) in the following simple form:where , is used for a continuous real valued vector function, which additionally satisfies the Lipschitz condition and stands for initial state vector. Taking integral on both sides of (18) we get

In order to formulate an iterative scheme, we consider a uniform grid on with is the step size and . Thus, the Eq. (19) gets the structure as follows after make use of the Euler method [34]

Thus, utilizing the above scheme (20), we deduced the following iterative formulae for the corresponding classes of the model (9) where

Further, we used the estimated parameters to obtain the simulations for the Caputo model (9) by varying fractional order p. The parameter values used in graphical interpretation are given in Table 1. In plots, the time level is taken up to 150 days. The behavior of the model (9) for different values of fractional order p is described in Fig. 3 . By making decrease in the value of arbitrary order p of the Caputo operator, a decrease in the infected compartments while increase in the susceptible compartment is observed. Fig. 4, Fig. 5 respectively describe the impact of disease transmission rate on the symptomatic infected () and asymptomatic infected () people for various cases of fractional order p. It is observed that the population in both infected classes is decreasing significantly with the decrease in as shown in subplots (a-d) of Fig. 4, Fig. 5. Further, the same behavior is observed for all values of fractional order p but the decrease in infected individuals is more significant for the smaller values of p.

Fig. 3

Dynamics of fractional COVID-19 model when .

Fig. 4

Influence of (disease transmission coefficient) on symptomatic infected individuals where (a) , (b) , (c) , (d) .

Fig. 5

The effect of (disease transmission coefficient) on asymptomatic infected individuals where (a) , (b) , (c) , (d) .

Conclusion

The dynamics of novel coronavirus (COVID-19) pandemic is presented using a mathematical model. The model was first analyzed in [15] with Atangana-Baleanu fractional operator. In this study, we reformulated the model using the Caputo fractional operator and provides a better theoretical and graphical analysis of the model. Additionally, in the present paper, we considered the reported COVID-19 confirmed cases from the early stage until the end of pandemic in the mainland China. The findings of the present study show a better agreement of the model prediction to the real reported infected cases and ultimately, a more suitable model parameters are estimated. Initially, the model is briefly explored with classical integer-order derivative and then the Caputo type fractional derivative is applied to reformulate the model to best described the disease dynamics. The basic mathematical results of the fractional model such as solution positivity, feasibility, equilibria and the most important biological threshold parameter basic reproduction number is analyzed. Further, we have shown that the COVID-19 model with Caputo operator is locally and globally stable if . Besides the basic analysis of the model, we further estimated the biological parameters using the reported confirmed infected cases in the mainland of China for the period January, 21 to April 9, 2020. It is worth mentioning that the present model fits very well with the real data of daily confirmed cases, as shown in Fig. 2 which reflects the reality of this infection in China. Based on the fitted parameters, the estimated value of the reproduction number is obtained as . Furthermore, using the real estimated parameters we carried out the numerical simulation and depicted the graphical results of the fractional model for various values of p. We also analyzed the influence of model parameters on the symptomatic and asymptomatic infected individuals for four different values of fractional order p. It is observed that by decreasing the parameter (the disease transmission coefficient) the infected individuals decrease significantly. The same behavior is observed for the fractional order p, but the decrease in the infected individuals is more significant for smaller values of fractional order of the Caputo operator p. We believe the work presented in this paper for the mathematical analysis of coronavirus with real cases will be more useful for the health authorities or the other decision making agencies in combating the disease. The authors wish to continue this study in the near future by extending this work using the control measures by applying the optimal control theory in order to get the possible elimination or decrease the infection of COVID-19 pandemic.

CRediT authorship contribution statement

Muhammad Awais: Conceptualization, Methodology, Writing - original draft. Fehaid Salem Alshammari: Conceptualization, Methodology, Writing - original draft. Saif Ullah: Conceptualization, Methodology, Writing - original draft. Muhammad Altaf Khan: Conceptualization, Methodology, Writing - original draft, Supervision, Writing - review & editing. Saeed Islam: Conceptualization, Methodology, Writing - original draft.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.