Modeling and analysis of COVID-19 epidemics with treatment in fractional derivatives using real data from Pakistan.

Coronaviruses are a large family of viruses that cause different symptoms, from mild cold to severe respiratory distress, and they can be seen in different types of animals such as camels, cattle, cats and bats. Novel coronavirus called COVID-19 is a newly emerged virus that appeared in many countries of the world, but the actual source of the virus is not yet known. The outbreak has caused pandemic with 26,622,706 confirmed infections and 874,708 reported deaths worldwide till August 31, 2020, with 17,717,911 recovered cases. Currently, there exist no vaccines officially approved for the prevention or management of the disease, but alternative drugs meant for HIV, HBV, malaria and some other flus are used to treat this virus. In the present paper, a fractional-order epidemic model with two different operators called the classical Caputo operator and the Atangana-Baleanu-Caputo operator for the transmission of COVID-19 epidemic is proposed and analyzed. The reproduction number R 0 is obtained for the prediction and persistence of the disease. The dynamic behavior of the equilibria is studied by using fractional Routh-Hurwitz stability criterion and fractional La Salle invariant principle. Special attention is given to the global dynamics of the equilibria. Moreover, the fitting of parameters through least squares curve fitting technique is performed, and the average absolute relative error between COVID-19 actual cases and the model's solution for the infectious class is tried to be reduced and the best fitted values of the relevant parameters are achieved. The numerical solution of the proposed COVID-19 fractional-order model under the Caputo operator is obtained by using generalized Adams-Bashforth-Moulton method, whereas for the Atangana-Baleanu-Caputo operator, we have used a new numerical scheme. Also, the treatment compartment is included in the population which determines the impact of alternative drugs applied for treating the infected individuals. Furthermore, numerical simulations of the model and their graphical presentations are performed to visualize the effectiveness of our theoretical results and to monitor the effect of arbitrary-order derivative. © Società Italiana di Fisica and Springer-Verlag GmbH Germany, part of Springer Nature 2020.

Introduction

In the present world, more and more attention has been given to the research of epidemic diseases such as HIV, HBV, Ebola, H1N1 and malaria, and it is a big challenge to control the spread of epidemic diseases among the population. On the one hand, while the world continues to fight existing infectious diseases, on the other hand, changing world conditions causes the emergence of different types of viruses. The most recent and newest of these viruses is the novel coronavirus, which is called COVID-19, that appeared in early 2020 and is still not fully controlled. While up to now the biological origin of the disease is unclear, the first cases were traced back to December 2019 in the city of Wuhan in China. It is a virus that can cause lung disease and, when left untreated, causes diseases such as severe acute respiratory failure syndrome. Coronaviruses are viruses that most people encounter instantly in their lives. Human coronaviruses often cause mild to moderate upper respiratory diseases. Coronaviruses have three subgroups known as alpha, beta and gamma; there is also a fourth new group called delta coronaviruses, SARS-CoV. Human coronaviruses were first detected in the mid-1960s. Until 2020, the virus, which appeared only in Saudi Arabia, Qatar and Jordan, caused the death of three people. The case, which was reported in Wuhan city with 11 million population in Hubei province in China on December 31, 2019, has been found to be infected with a novel coronavirus that has never been seen before. According to the World Health Organization (WHO) reports [1], this virus is thought to be transmitted from animals to humans, such as SARS-CoV and MERS-CoV. Nowadays, the disease has been transmitted from person to person, and on August 31, 2020, the number of confirmed infected cases has reached 26,622,706 with almost 874,708 deaths worldwide so far. When the patients who were died were examined, the majority of them were found to be elderly patients or patients diagnosed with chronic heart, lung and kidney, Parkinson’s and diabetes. Coronaviruses can cause diseases in many different creatures. However, SARS-CoV can infect humans, monkeys and animals such as Himalayan civet cat, raccoon dog, cat, dogs and rodents. It can be transmitted easily like flu through the removal of viruses that come into contact with the mouth and nose after touching the infected material.

While mathematical models do not provide a cure for a given infectious disease, they can be used to replicate possible scenarios of the dynamic at hand. At present, all sectors in the world such as medical bodies, politicians, armies, law enforcement, business, chemists, physicists, engineers and many others are putting their efforts to help stop the spread of COVID-19; mathematicians are not left behind. New mathematical models that could be used for simulation, with the aim to predict the future behavior of the spread and flatten the curve of infection and deaths, are developed. Zu et al. [2] proposed and studied transmission patterns of COVID-19 in the mainland of China and the efficacy of different control strategies: a data- and model-driven study. They constructed a compartmental model and based on reported data from the National Health Commission of PR China during January 10–February 17, 2020; they estimated the model parameters. At the end, they predicted the epidemic trend and transmission risk of COVID-19. With the help of sensitivity analysis method, they also estimated the efficacy of several intervention strategies. In their study, they have found that the quarantine measures adopted by the Chinese government since January 23, 2020, were necessary and effective. Postponing the relaxation of isolation, early diagnosis, patient isolation, broad close-contact tracing and strict monitoring of infected persons could effectively control COVID-19 epidemic. Atangana [3], in his paper, modeled the spread of COVID-19 with new fractal–fractional operators with the aim of proposing the question: Can the lockdown save mankind before vaccination? He suggested a mathematical model taking into account the possibility of transmission of COVID-19 from dead bodies to humans and the effect of lockdown. In his paper, he considered three cases. The first case suggested that there is transmission from dead to the living (medical staffs as they perform postmortem procedures on corpses and direct contacts with during burial ceremonies). This case has no equilibrium points except for disease free equilibrium, a clear indication that care must be taken when dealing with corpses due to COVID-19. In the second case, he removed the transmission rate from dead bodies. This case showed an equilibrium point, although the number of deaths, carriers and infected grew exponentially up to a certain stability level. In the last case, he incorporated a lockdown and social distancing effect, using the next-generation matrix. He achieved a zero-reproduction number, with number of deaths, infected and carriers decaying very rapidly. This is a clear indication that if lockdown recommendations are observed, the threat of COVID-19 can be reduced to zero in few months. He used Italy’s data to guide the construction of the mathematical model. To include non-locality into mathematical formulas, differential and integral operators were suggested. Properties and numerical approximations were presented in details. Finally, the suggested differential and integral operators were applied to the model. Tang et al. [4] studied the effectiveness of quarantine and isolation that determined the trend of COVID-19 epidemic in the final phase of the current outbreak in China. In their study, they have seen that the uncertainty analyses reveal that the epidemic is still uncertain, and it is important to continue enhancing the quarantine and isolation strategy and improving the detection rate in mainland China that helped the country to tackle with COVID-19 outbreak during its peak time. On the one side, researchers are continuously working toward the development for the cure of COVID-19, while on the other side mathematicians proposed many models for the spread and control of COVID-19 that have been used for some decision making. Ahmet et al. [5] analyzed a mathematical model of coronavirus disease (COVID-19) by using numerical approaches and logistic model. In their paper, they have reviewed and introduced some models for COVID-19 that included important questions about the global health care and suggested important notes. They suggested three well-known numerical techniques for the solution of proposed equations; these are Euler’s method, Runge–Kutta method of order two (RK2) and of order four (RK4). Results based on the suggested numerical techniques and provided approximate solutions gave important key answers to this global issue. They have obtained the results for two countries, namely Turkey and Iraq. More interestingly, for both countries, Turkey and Iraq, the basic reproduction numbers have been reported recently by several groups; a research estimation by April 9, 2020, revealed that for Turkey is 7.4 and for Iraq is 3.4, which are noticeably increased from the beginning of the pandemic. Thus, they investigated the forecasting epidemic size for Turkey and Iraq using the logistic model. They concluded that the suggested model is a reasonable description of this epidemic disease. One key driver of the spread is direct contact with infected patients, object or corpses from COVID-19. While the world is still waiting for a possible vaccine, measures have been initiated in countries around the world such as lockdown, self-isolation and social distancing. Mathematical models rely on mathematical tools called differential and integral operators. Several have been suggested in the last decades as researchers recognized the complexity of nature and inadequacies of existing differential and integral operators. Chen et al. [6] developed a mathematical model which investigates the transmission dynamics of COVID-19. They considered in their model four compartments as bats, hosts, reservoir and people network to point out the infection source. Munster et al. [7] examined the key questions for effect analysis of COVID-19. Corman et al. [8] designed a workflow diagram to model the process in the absence of available COVID-19 isolates or original patient specimens. Sookaromdee et al. [9] developed a quadratic model that can help better understand how the disease can be controlled and managed. Shen et al. [10] presented a mathematical model to analyze COVID-19 processes and evaluated the basic reproduction number. They also compared the results with the severe acute respiratory syndrome (SARS) and Middle East respiratory syndrome (MERS) results and pointed out that the fatality rate (FR) of COVID-19 is lower than the FR of the SARS and MERS. Kucharski et al. [11] estimated the median reproduction number between 1.6 and 2.9 and concluded that COVID-19 has still important potential for ongoing human-to-human transmission which is the most significant issue. Zhou et al. [12] estimated and evaluated the basic reproduction number of COVID-19 in their study in the range (2.2–3.0). This has concluded that the disease is controllable with moderate–high transmissibility. Ming et al. [13] proposed an extended form of basic SIR epidemic model to overcome the burdens faced by healthcare system during COVID-19 transmission in Wuhan, China. Meanwhile, many studies have also dealt with how the virus first appeared in some regions where a huge number of people can be affected seriously such as USA [14], Vietnam [15], Thailand, Japan and South Korea [16]. Tang et al. [17] estimated the transmission risk of COVID-19 and its implication for public health interventions. Their estimations, which are based on likelihood and model analysis, show that the control reproduction number is high as 6.47 (95% CI 5.71–7.23). They have performed the sensitivity analysis of their study and suggested that the disease can be controlled by reducing the contact among the people, using quarantine procedure and isolation of the population. Peng et al. [18] provided epidemic analysis of COVID-19 in China by dynamical modeling. Wu et al. [19] provided a modeling approach for the nowcasting and forecasting the potential domestic and international spread of COVID-19 outbreak originated in Wuhan, China. They have estimated the size of the epidemic in Wuhan on the basis of the number of cases exported from Wuhan to cities outside mainland China and forecasted the extent of the domestic and global public health risks of epidemics, accounting for social and non-pharmaceutical prevention interventions. Tang et al. [20] provided an updated estimation of the risk of transmission of the novel coronavirus (COVID-19). By using time-dependent contact and diagnose rates, they refitted their previously proposed dynamics transmission model [17] to the data available and re-estimated the effective daily reproduction ratio that better quantifies the evolution of the interventions. They have estimated that the effective daily reproduction ratio has fallen below 1 and when the epidemics will peak. They have shown with their updated results that disease can be controlled by using the self-isolation strictly and persistence. Nadima et al. [21] proposed a COVID epidemic model for its transmission and control. They have calculated the and analytically. They have investigated the detailed stability analysis of their proposed model. Further, with their results they have found that the burden of disease spread can be reduced if the quarantined individuals can be managed properly than the isolated individuals.

Fractional calculus is a common field trying to understand the real-world phenomena that is modeled with non-integer-order derivatives. Using these types of operators, more effective and up-to-date studies have been revealing over time. In this context, fractional calculus theory and its illustrative applications are attracting attention all over the world day by day. New fractional operators that have different features have been defined and have been used extensively to model real-life problems. The emergence of the new operators in the literature can be considered as a result of the reproduction of new problems that model different types of real-life events. Fractional derivative operators that address the kind of nonlinear differential equations can be stated as non-local. There exist nowadays many types of fractional derivatives with and without singular kernels. The fractional derivative begins with Leibniz’s question in 1695. The list of the existing fractional derivatives is long. With singular kernels, we have the Caputo derivative [22], Riemann–Liouville derivative [22] and the Katugampola derivative [23]. Without singular kernels, we have two types: the fractional derivative with exponential kernel known as Caputo–Fabrizio fractional derivative (CF) [24] and the fractional derivative with Mittag–Leffler kernel known as the Atangana–Baleanu fractional derivative (ABC) [25]. Because of effective properties, fractional calculus has found wide applications to model dynamics processes in many well-known fields of science, finance  [26, 27], engineering, biology, medicine and many others [28-46]. The importance of dealing with fractional-order derivatives is the involvement of memory and hereditary properties that gives a more realistic way to model COVID-19 epidemics. Due to the memory effect, the non-integer models integrate all previous information from the past that makes it able to predict and translate the epidemic models more accurately. Saeedian et al. [47] formulated SIR epidemic model with the inclusion of memory effect and studied its behavior along the memory effect on the disease spread with the help of fractional derivatives. Ucar et al. [48] provided the dynamics of a kind of smoking model and its community health by considering the ABC fractional operator. Going by the antecedents, we have seen clearly that modeling of physical and real-life scenarios with the fractional-order derivatives is much more accurate when compared with the integer-order cases. This assertion has been demonstrated in a number of research papers, monographs and books, see, for example, [49-58]. In view of these achievements, we are motivated in this research work to model and analyze COVID-19 epidemics for disease transmission using the Caputo and ABC fractional-order operators. The choice of using the Caputo derivative is due to the fact that if the given function is a constant, then the Caputo derivative of that function gives zero. Primarily, the Caputo operator computes an ordinary differential equation, followed by a fractional integral to obtain the desired order of fractional derivative. More importantly, the Caputo fractional differential equation (FDO) permits the use of local initial conditions to be included in the derivation of the model. Furthermore, models with integer-order derivatives can be used to capture dynamical systems of infectious disease, when only the initial conditions are used to forecast future behaviors of the spread. However, when the scenario is unpredictable or cannot be described adequately, maybe due to some uncertainties which are inherent to many physical and real-world processes, integer-order derivatives and integrals are both deficient. In the case of COVID-19, there are many uncertainties, many unknowns and much misinformation that make it very difficult to really provide a suitable mathematical model with classical differentiation. In general, non-local operators are more suitable for such situations, as they are able to capture non-localities and some memory effects depending on whether power law, fading memory or crossover effects are included.

The potential aim of present study is to design, with the help of epidemiological modeling, a mathematical model for understanding transmission dynamics of COVID-19 using actual cases of the pandemic in Pakistan. Primarily, we aim to obtain the basic reproductive number and equilibria in preventing the epidemic spread in the country. Our motivation emerges from a number of recently conducted research studies [44, 59–61] in the literature focused on deterministic modeling of different diseases in various countries. Each of these studies consists of compartmental modeling; however, none investigates the inclusion of treatment class and its effects on the control of the epidemic. The COVID-19 model, in the present study, is based on the assumption of continuous treatment of affected individuals. The research findings of the present study may help governments and public health authorities to formulate strategic plans to reduce the immunization gaps and thus prevent the outbreaks in the future. Moreover, a number of new studies related to COVID-19 mathematical modeling are appearing nowadays wherein the proposed model can be considered a good addition with consideration of treatment class within those studies, thereby increasing interest of researchers belonging to the field of fractional calculus modeling and mathematical epidemiology.

In the current paper, we consider the dynamical transmission process of novel coronavirus by using the Caputo and Atangana Baleanu fractional derivatives. In Sect. 1, we introduce the novel model, some related papers stated in the literature and fractional-order derivatives used in our paper. Section 2 provides some preliminary results required for the formulation of proposed mathematical model. In Sect. 3, we define the mathematical model of COVID-19 and present the fundamental structure of the fractional-order model. In Sect. 4, we present the existence and the uniqueness of the solutions of the fractional-order coronavirus model. Also, in this section, we discuss the mathematical analysis of the proposed fractional-order COVID-19 epidemic model along with equilibrium points and the stability of equilibrium points. In Sect. 5, we fit the parameters by considering the real data which have been reported in Pakistan and determine the exact parameters which we use in the simulations. In Sects. 6 and 7, we describe the numerical simulation technique namely Adams–Bashforth scheme and its application to the stated model according to the Caputo and Atangana–Baleanu operators, respectively. In Sect. 8, we illustrate our main results by the graphical representations and discuss the memory trace. We give the conclusions and perspectives in Sect. 9.

Preliminaries

In this section, we give the fundamental definitions that can be used throughout the paper. These definitions generally explain the fractal–fractional derivative in the power kernel sense and Mittag–Leffler kernel sense.

Definition 1

A real function is said to be in the space if there exists a real number such that where and it is said to be in the space if and only if

Definition 2

The Riemann–Liouville form of fractional integral operator of order of a function is given byorwhere and is a well-known gamma function.

Definition 3

[62] The Riemann–Liouville form of fractional derivative of order of a function is given by

Definition 4

The Caputo fractional derivative of order of a function is given bywhere the operator satisfies:

and .

Definition 5

[63] The left and right Atangana–Baleanu (ABC) fractional derivatives in the frame of Caputo are given as, respectively:andwhere is Mittag–Leffler function, , is an arrangement function and .

Definition 6

The corresponding integral of the Atangana–Baleanu fractional derivative is defined aswhere is defined in Eq. (5).

Definition 7

The Laplace transform of the Caputo fractional derivative of a function of order is defined as

Definition 8

The Laplace transform of the function is defined aswhere is the two-parameter Mittag–Leffler function with . Further, the Mittag–Leffler function satisfies the following equation [64]

Mathematical model formulation

Modeling the dynamics of infectious diseases has become a topic of much interest in recent years. Such efforts are useful in disease control and in the prevention of outbreaks. In the modeling transmission dynamics of a communicable disease, it is common to divide the population into disjoint classes, namely compartments whose sizes change with time. We formulated a fractional-order compartmental model for the recently emerged virus, namely COVID-19, and then analyzed it. For the understanding of COVID-19 transmission dynamics, the total population N(t) is divided into eight sub-population compartments, namely susceptible, exposed, quarantined, asymptomatic, symptomatic, isolated, treated and recovered such that for all t. When a person is healthy and susceptible to the disease (denoted by S), exposed, when the person is in a latent period but not yet infectious (denoted by E), quarantined, refers to the separation of COVID-19 infected individuals from the general population when the population are infected but not infectious (denoted by Q), asymptomatic, those individuals in the population who does not show the symptoms but are in incubation period (denoted by A), symptomatic, when the individual got the infection and is infectious to others (denoted by I), isolated, describes those COVID-19-infected individuals who are separated from the population become symptomatic infectious (denoted by P) and recovered population (denoted by R). Therefore, the proposed model in ordinary differential equations takes the following form [17, 21]:For the infected individuals from the groups quarantined, asymptomatic, symptomatic or isolated, the transmission coefficients are , and , respectively. We consider the as a transmission rate along with the modification factors for quarantined , asymptomatic and isolated individuals. The exposed population decreases with quarantine at a rate of , and becomes asymptomatic and symptomatic at a rate and dies with natural death at a rate . The quarantined population is reduced by growth of clinical symptom at a rate of and transferred to the isolated class. is the recovery rate of quarantine individuals, and is the natural death rate of quarantined population. The exposed individuals become asymptomatic at a rate by a proportion k. The recovery rate of asymptomatic individuals is , and the natural death rate is . The symptomatic individuals are produced by a proportion of of exposed class after the exposer of clinical symptoms of COVID-19 by exposed individuals. is the isolation rate of the symptomatic individuals; is the recovery rate and natural death at a rate . The isolated individuals are come from quarantined community at a rate and symptomatic group at a rate . The recovery rate of isolated individuals is , disease-induced death rate is , and natural death rate is . The symptomatic individuals join the treatment class at a rate and recover at a rate with natural death rate . Furthermore, quarantined, asymptomatic, symptomatic, isolated and treated individuals recover from the disease at rates , , , and , respectively, and this population is reduced by a natural death rate .

The above ordinary differential model (11) is further extended to a fractional-order system of order , , including a recruitment rate of susceptible individuals into the region as per unit time and being the natural death rate. Thus, the proposed model in Caputo or Atangana–Baleanu-type fractional derivatives of order , is given bywith non-negative initial conditions

Stability analysis

In this section, the stability analysis of the proposed model (12) is discussed. Stability analysis of a system in epidemiology and immunology determines the behavior of the system in disease transmission. By stability analysis, one knows when and where the disease spreads by calculating the threshold quantity known as basic reproduction number denoted by .

Positivity and boundedness

In this subsection, the positivity and boundedness of the solution for the proposed model (12) is given, after that the basic reproduction number is obtained. Finally, the existence conditions and the stability results for both the equilibria are provided. Let and For the proof of the main theorem about the non-negativity of the solutions for model (12), we recall the following lemma.

Lemma 1

(Generalized Mean Value Theorem [65, 66]). Supposing that and Caputo fractional derivative for thenwith .

Remark 1

Suppose that and Caputo fractional derivative for It is clear from Lemma 1 that if then the function is non-decreasing and if then the function is non-increasing for all

Theorem 1

The solution of the proposed fractional-order model (12) along initial conditions (13) is unique and bounded in

Proof

The existence and uniqueness of the solution of systems (12)–(13) on the time interval can be obtained by the process discussed in the work of Lin [66]. Subsequently, we have to explain the non-negative region is positively invariant region. From model (12), we findIf then according to Eqs. (14) and Remark 1, the solution cannot escape from the hyperplanes and Also on each hyperplane bounding the non-negative orthant, the vector field points into , i.e., the domain is a positively invariant set.

In the next theorem, we will show the boundedness of the solution to the proposed model (12).

Theorem 2

The region

is a positive invariant set for system (12).

For the proof of the theorem, we have from system (12)This implies,Applying the Laplace transform to Eq. (15), we getwhich further givesFrom Eqs. (8) and (9), we infer that if thenThis shows that the total population i.e., the subpopulations and , are bounded. This proves the boundedness of the solution of system (12).

Determining the equilibria and their stabilities

The equilibrium points are obtained by equating to zero the right-hand side of system (12), asAfter simplification, the DFE, namely , where and EE , wherewithand

Basic reproduction number

For the local stability of the disease-free equilibrium, we first compute the basic reproduction number by using next-generation matrix method [67-70]. It is defined as the number of cases occurring in a population which is completely susceptible by any infectious individual. Biologically, if , then the infection will disappear, but if the infection exists and the disease persists. To determine which is considered as the spectral radius of the next-generation matrix , we assemble the compartments which are infected from system (12) and decomposing the right-hand side as , where is the transmission part, expressing the production of new infection, and is the transition part which describes the change in state. Therefore,andBy the next-generation matrix method [67-70], the matrices and at the disease-free equilibrium point are obtained by and , . This implies,Thus, we have the following expression for

Local stability of equilibria

In this subsection, we will provide the local stability results of equilibrium points in the form of theorems with proofs.

Theorem 3

The disease-free equilibrium of the proposed fractional-order COVID-19 epidemic model is locally asymptotically stable if and is unstable if .

To study the stability criterion of disease-free equilibrium , the general Jacobian matrix has been calculated at and is obtained as follows:Thus, the disease-free equilibrium is locally asymptotically stable if all the eigenvalues , of the matrix satisfy the conditionWe can evaluate these eigenvalues by solving the following characteristic equationwhere is an identity matrix and is the eigenvalue. Therefore, we get an equation of the formwhich can be arranged asAssigning thatwhere , and taking so that then we can writeTherefore, we have , which gives . Hence, when all the eigenvalues of have negative real parts. So, the DFE point is locally asymptotically stable when Moreover, when , which means ,This implies that there exists so that , which means there exists a positive eigenvalue of Eq. (19). Thus, it gives us if the DFE () is unstable.

Lemma 2

[75] Define the following characteristic equation

The following conditions make all the roots of the characteristic equation (24) satisfy Eq. (18):

For the condition for Eq. (24) is given as

For the conditions for Eq. (24) are either Routh–Hurwitz conditions or

For if the discriminant of polynomial namely is positive, the following conditions are the necessary and sufficient conditions satisfy Eq. (18): if

If then condition Eq. (18) is satisfied for all

For general is the necessary for condition Eq. (18) to be satisfied.

Theorem 4

The endemic equilibrium of the proposed fractional-order COVID-19 epidemic model is locally asymptotically stable if and unstable otherwise.

The Jacobian matrix evaluated at the endemic equilibrium is given bywhere , , , and By considering the characteristic equation we have the following formSince the polynomial as given by Eq. (25) has all coefficients and positive (see “Appendix”). Therefore, from condition (5) of Lemma 2 the positive equilibrium point is locally asymptotically stable. This completes the proof of Theorem 4.

Global stability of equilibria

The global existence of the solution for the fractional differential equation is a most important concern and is carried out in the following section.

Theorem 5

[66] Assume that the function satisfies the following conditions in the global space: Here, A and B are two positive constants and Then, the initial value problemhas a unique solution.

The function is Lebesgue measurable with respect to t on .

The function is continuous with respect to on .

The function is continuous with respect to on .

for all most every and all .

Lemma 3

([76, 77]) Let be a continuous and derivable function. Then, for any time instant ,andwhere .

Note that for , the inequalities in (27) and (28) become equalities. Now, we provide the global stability results of the equilibria in the following theorems by considering the Lyapunov direct method.

Theorem 6

The disease-free equilibrium of the proposed model (12) is globally asymptotically stable if and unstable when .

To prove this, we can write the fractional-order system (12) aswhere denotes the number of uninfected individual compartments and denotes the number of infected individual compartments. The global stability of the disease-free equilibrium is guaranteed if the following two conditions are satisfied: where is a Metzler matrix and is the positively invariant set with respect to model (12). According to Castillo-Chavez et al. [78], we check for aforementioned conditions. Now, we have from model (12)andClearly, whenever the state variables are inside . Also, it is clear that is a globally asymptotically stable equilibrium of the system . Hence, this proves that the disease-free equilibrium of the proposed model (12) is globally asymptotically stable.

For , is globally asymptotically stable,

for ,

Theorem 7

The endemic equilibrium of the proposed model (12) is globally asymptotically stable if .

To establish the global stability of the endemic equilibrium , we define a Lyapunov function L(t) given bywhere , when , , , , , , . This function is defined, continuous and positive definite for all . It can be verified that the equality holds if and only if , , , , and . Now, we have from Lemma 3This further implies from system (12)Using the endemic conditions given asin Eq. (29), we getThis implies,It follows that if , then we have from Eq. (30), . Therefore, is bounded and non-increasing. Further, we know that , if and only if , , , , and . Also, the limit of exits as Therefore, the maximum invariant set for is the singleton set . According to the LaSalle’s invariance principle [67-70], we know that all solutions in converge to . Therefore, the endemic equilibrium of the proposed model (12) is globally asymptotically stable when This completes the proof of Theorem 7.

Parameter estimation

Fitting of parameters is one of the important features during validation of an epidemiological model. This creates confidence in the acceptance of the model in order to use it for future prediction and to better understand the transmission dynamics of the underlying epidemic. Thus, we aim to explain the fitting of parameters through least-squares curve fitting technique in this section. Since there are 19 different parameters in the proposed model for COVID-19 pandemic, we have estimated two of the parameters, whereas the rest are best fitted based upon the real cases of COVID-19 pandemic throughout Pakistan (source http://covid.gov.pk/stats/pakistan). Two demographic parameters are (recruitment rate) and (natural death rate) which have been estimated. The average natural mortality rate of a Pakistani is 66.5 years (source https://www.worldlifeexpectancy.com/pakistan-life-expectancy) so this yields per day. Further, the total population of Pakistan is 212.2M; it may be assumed that , which is the limiting population when there is no existence of the pandemic. In this way, per day. Daily cases of the pandemic are considered from March 24 to April 20, 2020 during preparation of the research paper. As far as the initial conditions are concerned, the total population of Pakistan is taken to be , the initial exposed and quarantined population is taken as and this further assists us to determine remaining initial conditions for the other state variables with the help of the identity . In this connection, we have obtained , , , , and There are 19 biological parameters which have been estimated with the aid of least-square fitting method and this leads to yield a best fit of the proposed COVID-19 model’s solution to the actual cases of the pandemic as shown in Fig. 1. The average absolute relative error between COVID-19 actual cases and the model’s solution for the infectious class is tried to be reduced and the best fitted values of the relevant parameters have been achieved. Such a value for the error is approximately . Figure 1 shows the real COVID-19 cases by blue solid circles, whereas the best fitted curve of the model is shown by the black solid line. The biological parameters included in the model are listed in Table 1 along with their best estimated values obtained via least-squares technique. These parameters have finally produced the value of the basic reproduction number equivalent to for the real COVID-19 cases in Pakistan from March 24 to April 20, 2020.

Fig. 1

The daily COVID-19 cases time series in Pakistan from March 24 to April 20, 2020, and the best fitted curve from the proposed model

Table 1

Estimated and best fitted values of the parameters used in the proposed COVID-19 model

Parameter	Meaning	Value	Sources
Λ	Recruitment rate	8.7424e+03	Estimated
λ	Natural death rate	1/(66.5∗365)	Estimated
μ_Q	Modification factor for quarantined	0.124	Fitted
μ_A	Modification factor for asymptomatic	0.956	Fitted
μ_P	Modification factor for isolated	0.076	Fitted
ρ	Rate of joining treatment class	0.2	Fitted
σ	Diseases induced mortality rate	1e-03	Fitted
k	Proportion of asymptomatic individuals	1e-04	Fitted
β_1	Recovery rate from quarantined individuals	0.0000005539	Fitted
β_2	Recovery rate from asymptomatic individuals	0.0000000196	Fitted
β_3	Recovery rate from symptomatic individuals	0.0000001257	Fitted
β_4	Recovery rate from isolated individuals	0.0000001086	Fitted
β_5	Recovery rate from treated individuals	0.6461299316	Fitted
δ	Transmission rate	0.7925264407	Fitted
α_1	Rate at which the exposed individuals are diminished by quarantine	0.0000000032	Fitted
α_2	Rate at which the symptomatic individuals are diminished by isolation	0.0000001257	Fitted
r_1	Rate at which exposed become infected	0.0000230757	Fitted
r_2	Rate at which quarantined individuals are isolated	0.0000076749	Fitted
γ	Rate at which infected individuals are treated	0.0010169510	Fitted

Numerical scheme for the Caputo fractional COVID-19 model

This section of the paper focuses on the numerical simulation for the Caputo-type fractional-order coronavirus model in Eq. (12). The Adams-type predictor–corrector method [71–73, 80] that is used to achieve numerical simulation of the nonlinear system is proposed to obtain approximate solution of the mentioned model. The following form of Cauchy ordinary differential equation is considered with respect to the Caputo operator of order :where . Equation (31) can be converted to the following Volterra equation:By using this mentioned predictor–corrector scheme associated with the Adam–Bashforth–Moulton algorithm [72] to get the numerical solutions of the fractional coronavirus model, we can take and by letting it can be discretized as follows, i.e., the corresponding corrector formula [74] whereThen, the required step is to determine the corresponding predictor formula We can calculate the mentioned predictor formula aswhere

Numerical scheme for the ABC fractional COVID-19 model

In this section, the proposed COVID-19 epidemic model under the ABC fractional derivative of order is numerically simulated [79]. The method is used to obtain approximate solutions of the proposed model. To provide the estimated solution by means of this algorithm, the following nonlinear fractional differential equation with respect to the ABC fractional derivative of order :with the following initial conditionsApplying the fundamental theorem of fractional calculus, Eq. (34) can be converted to a fractional integral equation asAt a given point , we have from Eq. (36)Within the interval , the function , using the two-step Lagrange polynomial interpolation, can be approximated as follows:Under the approximation Eq.(38), we get from Eq.(37) asWithout loss of generality, we considerandTherefore,andThis implies from Eq. (39) after substituting Eqs. (40)–(41),Equation (42) gives the numerical scheme for Atangana–Baleanu fractional derivative in the sense of Caputo. Using this scheme for the numerical solutions of the proposed fractional coronavirus model, we get

Numerical simulations and discussion

In this section, we have investigated behavior of each state variable of the proposed COVID-19 model under the two fractional-order operators called the Caputo (having singular kernel) and Atangana–Baleanu–Caputo (having non-singular kernel). Profiles for the state variables have been graphically obtained based upon variation in the fractional-order parameter and the variation in some of the more important biological parameters of the model randomly taken from Table 1. In addition to this, 3D meshes and contour plots are shown to capture behavior of the basic reproduction number in order to better understand the transmission dynamics of the pandemic.

Figure 2 shows long-term behavior of the pandemic using the Caputo and the ABC operators. It has been observed that at the end of 60 days, the total number of symptomatic individuals would reach about under the Caputo operator, whereas this number gets for the ABC operator. Thus, ABC predicts about 240 more infectious cases in comparison with the prediction obtained under the Caputo operator. This seems to be in good agreement with the ongoing situation of this pandemic. We have shown dynamical behavior of each state variable from the proposed COVID-19 model in Figs. 3, 4, 5, 6, 7, 8, 9 and 10 for varying values of the fractional-order parameter . Figure 3 shows increasing behavior of the susceptible population for decreasing values of with ABC operator predicting more susceptibility. Similar sort of behavior is observed in quarantined and asymptomatic population as depicted in Figs. 5 and 6, respectively. However, a different kind of behavior exists for rest of the state variables of the model. Figure 7 depicts significance of the fractional-order parameter for the symptomatic population of the model. This plot shows substantial decrease in the number of infectious individuals for decreasing values of under both operators with Caputo predicting 490 more cases in comparison with ABC operator assuming . For behavior of isolated population, Fig. 8 depicts that for decreasing values of in both operators, the number of isolated individuals decreases with 102 more people isolated in case of Caputo when . Under the treatment profile, Caputo operator predicts one extra case as shown in Fig. 9 for . Further, the Caputo operator suggests that 34 more symptomatic cases get recovered with in comparison with ABC operator. However, we observed previously that the Caputo generates 490 additional infectious cases. Thus, the ABC gives better agreement with real situation in this regard.

Fig. 2

Long-term prediction for infectious population in Pakistan using the proposed COVID-19 model under the a Caputo and b ABC operators with

Fig. 3

Profile for the susceptible population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 4

Profile for the exposed population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 5

Profile for the quarantined population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 6

Profile for the asymptomatic population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 7

Profile for the symptomatic population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 8

Profile for the isolated population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 9

Profile for the treated population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

Fig. 10

Profile for the recovered population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of

It can be observed from Fig. 11 that for increasing values of (transmission rate), the symptomatic individuals increase under both operators with Caputo predicting 8923 cases, whereas ABC predicts 7840 cases for at the end of the chosen time interval, thereby recommending the use of ABC operator over the Caputo. Increasing values of (rate of joining treatment class) result in decreasing number of symptomatic individuals as depicted in Fig. 12. For , Caputo in Fig. 12a predicts 3789 symptomatic cases, whereas ABC in Fig. 12b suggests 3600 cases demonstrating effectiveness of ABC operator over the Caputo.

Fig. 11

Profile for the symptomatic I(t) population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of (transmission rate) while taking

Fig. 12

Profile for the symptomatic I(t) population using the proposed COVID-19 model under the a Caputo and b ABC operators with different values of (rate of joining treatment class) while taking

Using some randomly chosen parameters from Table 1, we have shown transmission behavior of the pandemic via basic reproduction number . As shown in Fig. 13, increasing and decreasing will produce greater than 2 which is an alarming sign. Thus, reducing transmission rate of the coronavirus is the most essential strategy in order to prevent the virus further spread. Similar sort of behavior and description of the associated in Figs. 14 and 15 can be interpreted based upon the biological parameters used for the analysis.

Fig. 13

Dynamical behavior of the basic reproduction number for varying values of (transmission rate) and (recovery rate from symptomatic class)

Fig. 14

Dynamical behavior of the basic reproduction number for varying values of (rate at which exposed become infected) and (rate at which quarantined are isolated)

Fig. 15

Dynamical behavior of the basic reproduction number for varying values of (recovery rate from quarantined class) and (recovery rate from asymptomatic class)

Conclusions

In this paper, modeling and analysis of the newly emerged coronavirus (COVID-19) transmission dynamics have been provided in fractional-order derivatives with treatment class for the infected population. The basic reproduction number has been computed by the next-generation matrix method which performs as a threshold parameter in the disease transmission and determines whether the disease persists or vanishes from the population. The positivity and boundedness of the solutions have been determined. Also, the stability conditions of the equilibrium points for the proposed fractional-order system have been discussed. Meanwhile, the global dynamics of the equilibria has been obtained by the Lyapunov functional approach method. Based on the mathematical analysis, the disease-free equilibrium is locally asymptotically stable when that means the infection will die out in the population. On the other hand, the infection spreads in the population when Numerically, the endemic equilibrium tends to be locally asymptotically stable when , which means that the infection will persist in the population. Moreover, the fitting of parameters through least-squares curve fitting technique has been done and the average absolute relative error between COVID-19 actual cases and the model’s solution for the infectious class is tried to be reduced, and the best fitted values of the relevant parameters have been achieved. Finally, Adams–Bashforth–Moulton method has been applied to carry out the numerical simulations for different values of the fractional-order of the proposed model. It has been demonstrated that physical processes are better well described using the derivatives of fractional order which are more accurate and reliable in comparison with the classical-order case. Hence, we replace the integer-order time derivative with the Caputo-type fractional-order derivative and Atangana–Baleanu fractional-order derivative. It has been shown that these operators have many advantages over the existing non-integer-order types. Our choice of using the Caputo and Atangana–Baleanu derivatives is because they allow traditional initial and boundary conditions to be included in the formulation of the problem. We have formulated the model based on the available information we gathered on media and in prints about the causative agent and model of transmission of the virus disease. A number of numerical results showing the behavior of the dynamics obtained for different instances of fractional-order have been reported.

Although the fractional-order COVID-19 epidemic model based upon Caputo and ABC operators provides sufficient information to understand the epidemic transmission process and help to determine the crucial factors for its spread, for more detailed analysis one needs to have new tools to disclose unnoticed behavior of such nonlinear epidemiological systems and thus the operators known as Caputo–Fabrizio, Atangana–Gomez, Atangana beta derivative, truncated M-derivative, fractal–fractional and others can be used in future research work. In addition to this, our future research would also devise new strategies for carrying out stability analysis of the epidemiological models based upon non-autonomous nonlinear ordinary differential equations with fractional-order derivatives.