Hybrid Method for Simulation of a Fractional COVID-19 Model with Real Case Application

In this research, we provide a mathematical analysis for the novel coronavirus responsible for COVID-19, which continues to be a big source of threat for humanity. Our fractional-order analysis is carried out using a non-singular kernel type operator known as the Atangana--Baleanu--Caputo (ABC) derivative. We parametrize the model adopting available information of the disease from Pakistan in the period 9th April to 2nd June 2020. We obtain the required solution with the help of a hybrid method, which is a combination of the decomposition method and the Laplace transform. Furthermore, a sensitivity analysis is carried out to evaluate the parameters that are more sensitive to the basic reproduction number of the model. Our results are compared with the real data of Pakistan and numerical plots are presented at various fractional orders.


Introduction
The novel coronavirus SARS-CoV-2, responsible for COVID-19, which is member of the family of Severe Acute Respiratory Syndrome (SARS) viruses, has been recognized as the most dangerous virus of this decade [1]. This virus has become the new novel strain of the SARS family, which was not recognized in humans before [2]. COVID-19 has not just affected humans, but also a number of animals have been infected by the virus. The SARS-CoV-2 virus has been transmitted from human to human and similarly in animals, but its origin is still a controversy [3]. Infected humans and different species of various animals are recognized as active causes of spreading of the virus [1]. In the past, some similar viruses, like the Middle East Respiratory Syndrome Coronavirus (MERS-CoV), were spread out from camels to human population, and for SARS-CoV-1 the civet was recognized as the source of spreading into humans. For COVID-19, the main source or the major reason of spreading is human-to-human interaction, where the virus transmission is easily made by an infected person to a susceptible one. Currently, thousands of research studies have been proposed and many predictions have been given on COVID-19 dynamics, see in [4][5][6][7][8] and references therein. Our paper is, however, different from those in the literature. In [4], special focus is given to the transmissibility of the so-called superspreaders, with numerical simulations being arXiv:2111.01031v1 [math.DS] 1 Nov 2021 given for data of Galicia, Spain, and Portugal. It turns out that, for each region, the order of the Caputo derivative takes a different value, always different from one, showing the relevance of considering fractional-order models to investigate COVID- 19. The work in [5] studies the COVID-19 pandemic in Portugal until the end of the three states of emergency, describing well what has happen in Portugal with respect to the evolution of active infected and hospitalized individuals. In [6,7], a non-fractional but stochastic time-delayed model for COVID-19 is given, with the aim to study the situation of Morocco. In [8], the authors provide a S-E-I-P-A-H-R-F model while here we propose a much simpler P-I-Q model (our model has only three state variables while the model in [8] is much more complex, with eight state variables). In [8], the authors use the classical operator of Caputo; differently, here we use the more recent ABC derivatives that, in contrast, use non-singular kernels that allow us to consider a much simpler model. While the main result in [8] is the proof of the global stability of the disease free equilibrium; in contrast, here we prove Ulam-Hyers stability. We also construct a practical algorithm to compute numerically the solution of the model (see Section 5), while such algorithmic approach is not addressed in [8]. Moreover, we do a sensitivity analysis to the parameters of the model. Such sensitivity analysis is also not addressed in [8].
In contrast with the work in [8], that investigates the realities of Wuhan, Spain, and Portugal, we study the case of Pakistan.
COVID-19 generally transfers by interaction with humans in close contact for a particular time period, with most common symptoms of sneezing and coughing. The virus droplets stay on the layer of matters and when they come to the contact with any susceptible human, then the virus symptoms easily transfer to the individuals. Such infected humans can pass the infection to others by touching their mouth, eyes, or nose. This virus has the strength to be alive on different surfaces, like cardboard and copper, for many hours up to some days. As the time passes, the amount of the virus symptoms decreases over a time span and might not be alive in sufficient amount for spreading the infection. It has been recorded that the symptoms appearance and COVID-19 infection initial stage lies between 1 to 14 days [1]. Several countries have prepared and implemented a COVID-19 vaccine program and are trying to protect their populations. However, to date, there is yet no treatment available. At present, the most effect way to protect ourselves from the virus remains the quarantine or isolation, effective use of mask, following the guidelines that have been passed by governments of all countries along with the World Health Organization (WHO).
Modeling of infectious diseases has a rich literature and a number of research articles have been developed, both using classical dynamical systems as well as fractional models [4,8]. Fractional-order derivatives can be useful and helpful as compared to classical derivatives, because the dynamics of real phenomena can be comprehensively understood by fractional-order derivatives due to its special properties, i.e., hereditary and memory [9][10][11][12][13][14][15][16]. For a comparison between classical (integer-order) and fractional-order models, see in [4,8]. Roughly speaking, ordinary derivatives cannot distinguish the phenomenon at two distinct closed points. To sort out this problem of ordinary derivatives, generalized derivatives have been introduced in the framework of fractional calculus [17]. The first concept of fractional-order derivative was given by Leibniz and L'Hôpital in 1695. Aiming quantitative analysis, optimization, and numerical estimations, many number of attempts have been made by employing fractional differential equations (FDEs) [1,2,. The growing interest in the modeling of complex real-world issues with the use of FDEs is due to its numerous properties that can not be found in the ordinary sense. These characteristics allow FDEs to model effectively not only non-Markovian processes but also non-Gaussian phenomena [39]. Different non-classical fractional-order derivatives and different kinds of FDEs were proposed [40][41][42]. Among them, one has the Atangana-Baleanu-Caputo (ABC) derivative, which is a nonlocal fractional derivative with a non-singular kernel, connected with various applications. For a discussion of the ABC and related operators see in [43,44], and for their use on contemporary modeling we refer the interested reader to the works in [45][46][47][48].
The famous method of decomposition was developed from the 1970s to the 1990s by George Adomian, to analytically handle nonlinear problems. After that, the Adomian decomposition method became a powerful tool to simulate analytical or approximate solutions for various problems of an applied nature. Many mathematical models have been studied by the applications of homotopy, Laplace Adomian Decomposition Method (LADM), and variational methods [49][50][51]. To the best of our knowledge, no one has studied a variable order epidemic model with ABC derivatives by the LADM. Motivated by this fact, here we study a fractional-order COVID-19 epidemic model with ABC derivatives by the Laplace Adomian decomposition algorithm. In particular, we use Banach and Krassnoselskii fixed point theorems to define some sufficient conditions to prove existence and uniqueness of solution. As stability is important for the estimated solution, we consider Ulam type stability through nonlinear functional analysis. The aforementioned stability is investigated for ordinary fractional derivatives in many research papers, see, e.g., in [52][53][54], but research on Ulam type stability regarding ABC derivatives is a rarity. At the end of the paper, our results are illustrated with real data based on Pakistan COVID-19 cases in March 2020.
The paper is organized as follows. Section 2 is devoted to the model formulation. Section 3 is concerned with some preliminary results on fractional differential equations. Existence and uniqueness are carried out in Section 4. Section 5 deals with the solution of the COVID-19 model using the LADM. Some plots are given in Section 6, showing the simplicity and reliability of the proposed algorithm. In Section 7, a sensitivity analysis is given to find the most sensitive parameter with respect to the basic reproduction number. We end with Section 8 of conclusions, including some possible future directions of research.

Model Formulation
Mathematical modeling plays a major role in investigating and thus controlling the dynamics of a disease, particularly in the vaccination privation or at the initial phases of the epidemic. Several mathematical models can be found in [12][13][14][15]. We formulated a fractional COVID-19 epidemic model, similar to other diseases [49,50], and predict its future behavior. Inspired by FDEs using the ABC derivative, we aim to simulate the COVID-19 transmission in the form of along with initial conditions where ABC D θ t is the ABC fractional derivative of order 0 < θ ≤ 1 (see Definition 1 in Section 3). In this model, P(t) represents the amount of susceptible humans, I(t) stands for the population of infected humans, and Q(t) represents the population of quarantined humans at time t. The meaning of the parameters of model (1) are given in Table 1. We take the below assumptions to the given system: The basic reproduction number R 0 , which represents the secondary cases for the model (1), is easily demonstrated to be given by .
(3) In addition, where N represents the total population.

Preliminary Results
For completeness, here we recall necessary definitions and results from the literature.
Definition 1 (See [11,48]). If x is an absolutely continuous function and 0 < θ ≤ 1, then the ABC derivative is given by where ABC(θ) is a normalization function such that 1 = ABC(1) = ABC(0) and M θ is a special Mittag-Leffler function.
Remark 2. Let x(t) be a function having fractional ABC derivative. Then, the Laplace transform of ABC D θ 0 x(t) is given by Lemma 1 (See [52]). The solution to Theorem 1 (See [54]). Let X = C[0, T] and consider the Banach space defined by Z = X × X × X with the norm-function A = (P, . Consider B to be a convex subset of Z and F, and G be operators such that F is a contraction, and 3.
G is compact and continuous.
Then, Fu + Gu = u possesses at least one solution.

Qualitative Analysis of the Proposed Model
Here, we rewrite the right-hand sides of (1) as By using (5), we have In view of Lemma 1, (6) yields Using (7) and (8), we define the two operators F and G as follows: For existence and uniqueness, we assume some basic axioms and a Lipschitz hypothesis: Proof. The theorem is proved in two steps, with the help of Theorem 1. (i) ConsiderĀ ∈ B, where B = {A ∈ Z : A ≤ ρ, ρ > 0} is a closed and convex set. Then, for F in (9), we have Therefore, F is a contraction. (ii) We want G to be relatively compact. For that it suffices that G is equicontinuous and bounded. Obviously, G is continuous as Φ is continuous and for all A ∈ B one has Thus, (11) shows the boundedness of G. For equi-continuity, we assume t 1 > t 2 ∈ [0, τ], so that The right-hand side in (12) goes to zero at t 1 → t 2 . Remembering that G is continuous, Having the boundedness and continuity of G, we conclude that G is uniformly continuous and bounded. According to the theorem of Arzelá-Ascoli, G is relatively compact and therefore entirely continuous. It follows from Theorem 1 that the integral Equation (7) has at least one solution. Now, we show uniqueness.
Proof. Let the operator T : Z → Z be defined by (13) Let A,Ā ∈ Z. Then, one can take where Thus, T is a contraction from (14). Therefore, (7) possesses a unique solution.
Next, in order to investigate the stability of our problem, we consider a small disturbance φ ∈ C[0, T], with φ(0) = 0, that depends only on the solution.
Proof. Assume A ∈ Z and letĀ ∈ Z be the unique solution of (7). Then, From (18), we can write that The proof is complete.

Construction of an Algorithm for Deriving the Solution of the Model
Herein, we derive a general series-type solution for the proposed system with ABC derivatives. Taking the Laplace transform in model (1), we transform both sides of each equation and we use the initial conditions to obtain that Now, considering each solution in the form of series, we separate the nonlinear term P(t)I(t) in terms of Adomian polynomials as Therefore, from (21) and (22), we obtain from (20) that Now, comparing the terms on both sides of (23), one has . . .
Applying the inverse Laplace transform to (24), we obtain that

Numerical Interpretation and Discussion
To illustrate the dynamical structure of our infectious disease model, we now consider a practical case study under various numerical observations and given parameter values. The concrete parameter values we have used are shown in Table 2. Quarantined population at t = 0 0.0011 millions We assume that the initial susceptible, infected, and isolated populations are 10, 0.01, and 0.0011 million, respectively. Among the 21,000 selected population, the density of susceptible population is about 0.6 percent, the infected population is 0.2 percent, and the isolated population is 0.2 percent.
By using the parameter values in Table 2, we computed the first three terms of the general series solution (25) with ABC(θ) = 1 as We have utilized the numeric computing environment MATLAB, version 2016, and plotted the solution (25) in Figure 1 by considering the first fifteen terms of the series (21). Figure 1 shows the dynamics of each one of the state variables in the classical sense when θ = 1 (green curves). Similarly, by considering the model in ABC sense, that is, for θ ∈ (0, 1), we plot in Figure 1 each of the state variables to analyze the changes in comparison with the classical case. From Figure 1(a), we see that as we increase the order θ of the fractional ABC derivative, the susceptibility increases. Further, all fractional order derivatives shows no effect after about 60 days, i.e., the susceptible population stabilizes. Figure 1(b) shows that infected individuals tend to increase, with different rates, when we decrease the fractional-order: the smaller the fractional order θ, faster the increase rate, and vice versa. All obtained curves for infected individuals, for different values of the fractional order derivatives, approach towards a non-zero steady state, which shows that the disease will persist in the community if not properly managed. On the other hand, Figure 1(c) shows that during the first month the disease progress with more and more people getting quarantined, irrespective of the order of the derivative. However, after that, the quarantined population tends to decline and, at the end, there will be no quarantined individuals in the community.
In Section 6.1, we show that the fractional model (1) with ABC derivatives has the ability to describe effectively the dynamics of transmission of the current COVID-19 outbreak.

Case Study with Real Data: Khyber Pakhtunkhawa (Pakistan)
The Khyber Pakhtunkhawa Province, like other provinces of Pakistan and the rest of the world, is also being affected by COVID-19. We decided to calibrate our model with real data of COVID-19 from Khyber Pakhtunkhawa, Pakistan, from 9th April to the 2nd of June 2020. For that, we have used the minimization method of MATLAB taking the initial weights P(0) = 35, 525, 047, I(0) = 10, 485, Q(0) = 18000, determined from the work in [55], and θ = 1, from which we arrived to the values of the parameters shown in Table 3. Figure 2 shows the total number of individuals infected by COVID-19 as registered from 9th April to the 2nd in June 2020, which corresponds to the period of one month and 24 days used to calibrate our model. Figure 3 compares the actual/real data of COVID-19 with the curve of infected given by model (1), clearly showing the appropriateness of our model to describe the COVID-19 outbreak.   Figure 4 projects the long-term behavior of the COVID-19 outbreak during a period of eight months. We can see the data matches during the first 1.8 months and, additionally, we observe that the long-term behavior consists on a rise of infected individuals with time. This means that if the government did not apply proper strategies, the incidence could increase drastically in the coming months.

Sensitivity Analysis
Here, we conduct a sensitivity analysis to evaluate the parameters that are sensitive in minimizing the propagation of the ailment. Although its computation is tedious for complex biological models, forward sensitivity analysis is recorded as an important component of epidemic modeling: the ecologist and epidemiologist gain a lot of insight from the sensitivity study of the basic reproduction number R 0 [56]. In Definition 2, we assume that the basic reproduction number R 0 is differentiable with respect to parameter ω. Given (3), this means that Definition 2 makes sense for ω ∈ {γ, λ, d 0 , µ, σ, h, η}. Definition 2. The normalized forward sensitivity index of R 0 with respect to parameter ω is defined by As we have an analytical form for the basic reproduction number, recall (3), we apply the direct differentiation process given in (27). Not only do the sensitivity indexes show us the impact of various factors associated with the spread of the infectious disease, but they also provide us with valuable details on the comparative change between R 0 and the parameters. Moreover, they also assist in the production of control strategies [57]. Table 4 demonstrates that γ, h, and σ parameters have a positive effect on the basic reproduction number R 0 , which means that the growth or decay of these parameters by 10% would increase or decrease the reproduction number by 10%, 6.36%, and 0.31%, respectively. On the other hand, d 0 -, µ-, and η-sensitive indexes indicate that increasing their values by 10% would decrease the basic reproduction number R 0 by 14.89%, 0.09%, and 1.68%, respectively. Table 4. Sensitivity indexes of the basic reproduction number R 0 (3) (see Definition 2) for relevant parameters of model (1).

Parameters Sensitivity Value Parameters Sensitivity Value
The sensitivity of the basic reproduction number R 0 is also seen graphically in Figure 5.

Conclusions and Future Work
In this manuscript, we studied a COVID-19 disease model providing a detailed qualitative analysis and showed its usefulness with a case study of Khyber Pakhtunkhawa, Pakistan. Our sensitivity analysis shows that the transmission rate γ has a huge effect on the model as compared to other parameters: the basic reproduction number varies directly with the transmission rate γ. The sensitivity analysis also showed that the death rate parameter µ has no effect on spreading the infection, which seems biologically correct. The transmission rate will be small by keeping a social distancing and self-quarantine situation that causes a decrease in the infection. In this way, one can control COVID-19 infection from rapid spreading in the community. In the future, we plan to analyze optimal control techniques to reduce the population of infected individuals by adopting a number of control measures. A modification of the given model is also possible by introducing more parameters for analyzing the early outbreaks of COVID-19 and then transmission and treatment aspects can be recalled. The given system can be also simulated by adding exposed and hospitalized classes and taking a stochastic fractional derivative. Here, we have provided a case study with real data from Pakistan, but other case studies can also be done.