Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.015375

ARTICLE

Fractional Analysis of Dynamical Novel COVID-19 by Semi-Analytical Technique

S. Iqbal1, D. Baleanu2,3, Javaid Ali4, H. M. Younas5 and M. B. Riaz6,7,*

1School of Systems and Technology, University of Management and Technology, Lahore, 54770, Pakistan
2Department of Mathematics, Cankaya University, Ankara, 06530, Turkey
3Institute of Space Sciences, Magurele, Bucharest, 077125, Romania
4Department of Mathematics, Govt. College Township, Punjab Higher Education Department, Lahore, 54770, Pakistan
5Department of Mathematics, Islamia University of Bahawalpur, Bahawalpur, 63100, Pakistan
6Department of Mathematics, University of Management and Technology, Lahore, 54770, Pakistan
7Institute for Groundwater Studies, University of the Free State, Bloemfontein, 9300, South Africa
*Corresponding Author: M. B. Riaz. Email: bilalsehole@gmail.com
Received: 14 December 2020; Accepted: 24 June 2021

1  Introduction

Humans have invented many scientific methods so far and have set in motion several steps to avoid, even to cure, some of the lethal diseases. Although they believed they had conquered nature, corona-virus appeared killing thousands of people in China. Coronavirus has also been spread in many countries from Africa to Europe. Coronavirus 2019 (COVID-19) is a contagious virus causing infection to the respiratory system and is widely spread from humans to humans. The first infected case of this new COVID-19 disease was identified on December 31, 2019 in the city of Wuhan, China, the capital of Hubei province [1]. Reportedly, main symptom was development of the pneumonia without any diagnosable cause and the available therapies and vaccines were not found effective [2]. It was also noticed that the transfer of the virus among human takes place due to their mutual contact [3]. Following the spread of COVID-19 virus in Wuhan City of China, the virus was also spread to other Chinese cities rapidly. In turn, the virus was proliferated to other areas including Asia, Europe, and North America. It is known that it takes 2 to 10 days for signs to emerge. The symptoms include troubled breathing, fever and coughing. A total number of 4593 infected cases and 132 deaths were reported until January 28, 2020 and more details on the latest infections and deaths due to the virus were yet to come.

Mathematical models play a vital role not only in understanding the dynamics of infection but also in investigating the recommendable conditions under which the disease will persist or wiped out. Presently, governments and researchers have shown great concerns to COVID-19 because of high transmission rate and noteworthy disease induced death rate. The COVID-19 virus is generally transferred when an infected person releases droplets generated by sneezing, coughing or exhaling. The confirmed cases of COVID-19 have reached nearly fifty four million all around the globe, and more than 1.3 million deaths have been caused by this virus. As of November 13, 2020, according to Worldometers [4] USA, Brazil and India are the top countries with maximum death tolls amounting to 2.48, 1.64 and 1.28 million, respectively.

Keenly tracking the corona virus transmission, researchers have organized to speed up the diagnostic processes, and several types of vaccines are under investigation against COVID-19. For example, Cao et al. [5,6] studied and discussed the outcomes features of victims of COVID-19 in ICUs. Nesteruk [7] proposed a susceptible-infected-recovered (SIR) epidemic model and extracted some valuable guidelines through statistical analysis of the model parameters. An amended SIR model of COVID-19 was proposed in another study [8] to predict the exact count of infections and the additional burden on ICUs and isolation units. The count of coronavirus cases appeared higher than the number of cases predicted by February 2020. It lead to exert more emphasize on focusing the research on modernizing the corona virus predictions in the view of latest data and considering more complex mathematical models. Presently, no patent therapeutic agents or vaccines for treatment or prevention from coronavirus are available, however the research based investigations into vaccine candidates and potential antivirals are ongoing in several countries. Drug development is a comparatively shorter process as compared to development, testing and distribution of vaccine and any cure for COVID-19 is not expected to be available before 2021. The dense places are ideal environments where the virus can spread easily. It was agreed that human close contact is one of the possible causes of COVID-19 outbreaks. Therefore, the quarantine of the COVID-19 victims can lessen the threat of further disease proliferation. The important measures taken to reduce the spread of virus include small contact rate and social distancing. The impact of such measures was investigated by Zeb et al. [9] by proposing and analyzing a SEIIR (susceptible-exposed-infected-isolated-recovered) model. They concluded that the individuals infected from coronavirus infectious disease must be referred to isolated compartment at various rates. A logistic growth model of COVID-19 was studied by Batista [10] and was employed to predict the ultimate volume of the epidemic. The dynamical behavior of coronavirus has been studied by several researchers through various COVID-19 transmission models [7–11].

Over the last thirty years, fractional derivatives have captivated the numerous researchers after recognition of the fact that in comparison to the classical derivatives, fractional derivatives are more reliable operators to model the real world physical phenomenons. In dynamical problems, fractional calculus (FC) based modeling is receiving a rapid popularity nowadays. The mathematical modeling of many physical and engineering models based on the idea of FC exhibits highly precise and accurate experimental results as compared to the models based on conventional calculus. The non-integer differential operators such as Caputo, Caputo-Fabrizio and ABC are fractional operators that transform the ordinary model to generalized model. In this article, we present a novel research on a fractional order dynamical model that underpins the propagation of coronavirus infectious disease and provide some forecasting with real world data. We extend an integer-order model formulation to a fractional order model by adding the Caputo sense of fractional derivative. The reason of using the Caputo fractional derivative is that it possesses several basic characteristics of fractional calculus. Moreover, the transmission behavior described in the model can be better defined by using Caputo operator. Research based on derivatives in Caputo sense and its applications to different models emerging in various disciplines of engineering and other sciences can be observed in several past studies [12–16]. Few additional linked researches on Caputo and other fractional order derivative implementations can be found in literature [17–32]. Since the number of infected cases reported from January 21 to January 28, 2020 is comparatively higher than those of on initial days, therefore we have considered this period for the formulation of the mathematical model as a parameterized model. In Section 2, transitory details of various mathematical models showing evolution of COVID-19 from bats to humans have been provided. Section 3 presents basic properties of the fractional order COVID-19 model. The Section 4 consists of the development of a semi-analytical scheme for the solution of the considered fractional order COVID-19 model. The simulation results of the proposed scheme based upon model fitted data are presented in Section 5. In the end, some fruitful conclusions have been presented.

2  Mathematical Modeling of COVID-19 Evolution

Presuming that the transmission occurs primarily within the population of bats and afterwards the transmission occurs to wild animals usually termed as hosts. Hunting of these carriers and then their transportation to the supply markets of seafood are considered as virus reservoirs. By exposing to the market, people get the risk of infection. In the following subsections, we revisit the evolution of the COVID-19 evolutionary tracks from bats to humans in the form of three mathematical models.

2.1 Model 1: Transmission of Corona Virus-19 among Bats and Hosts Populations

From the mathematical modeling point of view, let us denote the size of entire population of bats by Nb and further classify it as four subclasses, namely, the subpopulation of susceptible bats denoted by Sb, the subpopulation of exposed bats denoted by Eb, the subpopulation of bats infected by the virus symbolized as Ib and the subpopulation of removed/recovered bats denoted by Rb. The size of unknown host population is denoted by Nh and is further categorized into four subgroups Sh; Eh; Ih and Rh, respectively representing the susceptible subpopulation, exposed subpopulation, infected subpopulation and the recovered or removed subpopulation of hosts. Therefore, Nb=Sb+Eb+Ib+Rb and Nh=Sh+Eh+Ih+Rh.

The susceptible bat population is hired via birth rate Πb where death rate is μb for all types of bats. The exposed bats get infected at rates θb after having completed their incubation period and therefore, get included in the infected subpopulation Ib. The rate of removal or recovery of infected subpopulation of bats is denoted by τb at which the subclass Rb is populated. The infection due to the contact of infected population with the susceptible population takes place at a rate ηb and is denoted by the route ηbSbIb/Nb. We denote the birth rate of unknown hosts by Πh and the natural mortality rate by Ih. The exposed members of host subpopulation get infection at a rate θh and enter into the infected subpopulation (Ih), whereas τh denotes the infected host's rate of removal/recovery. The path ηbhShIb/Nh ion subpopulation Ih, models the interactions taking place among the infected bats and susceptible hosts, with the assumption that ηbh is the transmission coefficient of disease towards healthy hosts from infected population of bats. As a result of infection caused by infected members of bats population, the virus is potentially capable of dispersing inside the host's population and is mathematically represented by the term ηhShIh/Nh, where ηh denotes the disease transmission coefficient between the subgroups Sh and Ih of the host population.

dSb(t)dt=Πb−μbSb(t)−ηbSb(t)(Ib(t))Nb(t)(1)

dEb(t)dt=ηbSb(t)(Ib(t))Nb(t)−(μb+θb)Eb(t)(2)

dIb(t)dt=(θb)Eb(t)−(τb+μb)Ib(t)(3)

dRb(t)dt=τbIb(t)−μbRb(t)(4)

dSh(t)dt=Πh−μhSh(t)−ηbhSh(t)(Ib(t))Nh(t)−ηhSh(t)(Ih(t))Nh(t)(5)

dEh(t)dt=ηbhSh(t)(Ib(t))Nh(t)+ηhSh(t)(Ih(t))Nh(t)−(μh+θh)Eh(t)(6)

dIh(t)dt=(θh)Eh(t)−(τh+μh)Ih(t)(7)

dRh(t)dt=τhIh(t)−μhRh(t)(8)

Eqs. (1)–(8) are subject to some non-negative initial conditions.

2.2 Model 2: Transmission of COVID-19 from Seafood Market to Human

Let Np(t) denote the total population of individuals, which is again categorized into five subgroups such as Sp(t);Ep(t);Ip(t);Ap(t) and Rp(t) respectively representing the susceptible subpopulation, exposed subpopulation, infected (symptomatic) subpopulation, asymptomatically infected subpopulation, and the recovered or removed subpopulations of human. Simply, the following system of coupled system of nonlinear differential equations defines the evolutionary track of COVID-19 from hosts to human via seafood market reservoir.

dSp(t)dt=Πp−μpSp(t)−ηpSp(t)(Ip(t)+ψAp(t))Np(t)−ηωSp(t)M(t)(9)

dEp(t)dt=ηpSp(t)(Ip(t)+ψAp(t))Np(t)+ηωSp(t)M(t)−(1−θp)ωpEp(t)−θpρpEp(t)−μpEp(t)(10)

dIp(t)dt=(1−θp)ωpEp(t)−(τp+μp)Ip(t)(11)

dAp(t)dt=θpρpEp(t)−(τap+μp)Ap(t)(12)

dRp(t)dt=τpIp(t)+τapAp(t)−μpRp(t)(13)

dM(t)dt=bM(t)Ih(t)Nh(t)+QpIp(t)+ΩpAp(t)−πM(t)(14)

The governing Eqs. (9)–(14) are also subject to some initial conditions. The parameter τap represents the removal or recovery rate of Ap, ηω is the coefficient of transmission of disease from M to Sp and θp is the asymptomatic infection factor.

2.3 Model 3: Transmission of COVID-19 among Human with Reservoir

This model presents changing aspects of transmission of COVID-19 among human that are due to the close contacts of human population with the contagious environment without direct contact to virus hosts. Peoples' birth and natural death rates are denoted by the parameters ρp and Ip, respectively. The prone people Sp will be infected by ample encounters with the infected people Ip through the definition provided by ηpSpIp where the ηp is coefficient of disease transmission. Transmission amongst asymptomatically infected healthy individuals may occur in the form of ΨηpSpAp, where Ψ numerous transmissibility of Ap to the Ip such that Ψ∈[0,1]; if Ψ=0, there will be no numerous transmissibility, and if Ψ=1, the same will occur as Ip infection. The θp value is the asymptomatic infection factor. The ωp and ρp parameters respectively reflect the rate of transmission after the completion of the time of incubation and becoming infected, entering the classes Ip and Ap. The individuals in symptomatic class Ip and asymptomatic class Ap join the group Rp with the rate of removal or recovery by τp and τap, respectively. Section M is the reservoir or the marketplace location where the seafood is stored. The prone individuals are affected by ηωMSp after contact with M, where ηω is the coefficient of transmission of disease from M to Sp. The number of hosts who visit the seafood market for buying products (retail purchase) is indicated by b and the corresponding infection is modeled through the relation bMIh/Nh. In view of the fact that the COVID-19 can be introduced into the seafood market in a short time with ample source of infection and hence, without lack of generality, the effect of disease transmission due to bats’ contact to hosts can be ignored. Consequently, the following model is achieved:

dSp(t)dt=Πp−μpSp(t)−ηpSp(t)(Ip(t)+ψAp(t))Np(t)−ηωSp(t)M(t)(15)

dEp(t)dt=ηpSp(t)(Ip(t)+ψAp(t))Np(t)+ηωSp(t)M(t)−

(1−θp)ωpEp(t)−θpρpEp(t)−μpEp(t)(16)

dIp(t)dt=(1−θp)ωpEp(t)−(τp+μp)Ip(t)(17)

dAp(t)dt=θpρpEp(t)−(τap+μp)Ap(t)(18)

dRp(t)dt=τpIp(t)+τapAp(t)−μpRp(t)(19)

dM(t)dt=QpIp(t)+ΩpAp(t)−πM(t)(20)

There also exist some initial conditions of this model. Considering that the rates of symptomatic and asymptomatic infections due to reservoir M are constant, we denote these constant rates by parameters Qp and ϖp, respectively. The parameter π in Eq. (20) denotes the rate at which the virus is removal from environment reservoir as a result of some treatment policies implemented by policy makers.

3  Dynamical Fractional Order Model of COVID-19

To develop fractional order COVID-19 model, we describe some basic definitions from fractional calculus, which play vital role in fractional calculus for solving fractional order system of differential equations. These definitions consist of fractional integral operator of a function f(t) in Riemann–Liouville sense and the fractional derivative of a function f(t) in Caputo sense.

The integral operator of fractional order in Riemann–Liouville sense with order α ≥ 0 of f(t)∈Cμ is defined as beneath.

Icαf(t)=1Γ(α)∫ct⁡(t−μ)α−1f(μ)dμ,α>0,t>0

where Cμ is said to be space and μ∈R, μ≥−1. The following formula describes Caputo sense based fractional order derivative operator.

Dcαf(t)=1Γ(1−α)∫ct⁡(t−η)i−α−1fi(η)dη

Fori−1≤α≤i,i∈N,f(t)∈C−1i

Replacing the time derivatives of state variables by Caputo sense fractional order derivatives in Model 3, we obtain the following generalized COVID-19 model of fractional order in Caputo sense fractional derivative operator.

dαSp(t)dtα=Πp−μpSp(t)−ηpSp(t)(Ip(t)+ψAp(t))Np(t)−ηωSp(t)M(t)(21)

dαEp(t)dtα=ηpSp(t)(Ip(t)+ψAp(t))Np(t)+ηωSp(t)M(t)−(1−θp)ωpEp(t)−θpρpEp(t)−μpEp(t)(22)

dαIp(t)dtα=(1−θp)ωpEp(t)−(τp+μp)Ip(t)(23)

dαAp(t)dtα=θpρpEp(t)−(τap+μp)Ap(t)(24)

dαRp(t)dtα=τpIp(t)+τapAp(t)−μpRp(t)(25)

dαM(t)dtα=QpIp(t)+ΩpAp(t)−πM(t)(26)

Initial conditions obeyed by the model are:

Sp(0)≥0,Ep(0)≥0,Ip(0)≥0,Ap(0)≥0,Rp(0)≥0,M(0)≥0,

By adding Eqs. (15)–(19), we obtain human population's total dynamics as under.

dNp(t)dt=Πp−μpNp(t).

Integrating over [0,t] we get:

⇒Np(t)=Πpμp(1−e−μpt)

⇒limt→∞⁡Np(t)=Πpμp

It implies that the total population has an upper bound of Πpμp at any time step. The positivity of the state variables and the upper bound of the entire population constitute the following feasible region for the model.

Ω={(Sp(t),Ep(t),Ip(t),Ap(t),Rp(t))∈R+5:Np(t)≤Πpμp,M∈R+:ΠpμpQp+wpπ}

The points of equilibrium of the above fractional order dynamical COVID-19 model are calculated by solving the nonlinear algebraic equations obtained by equating the fractional time derivatives in Eqs. (21)--(26) to zero. Disease free and endemic equilibrium points denoted by D0 and D∗ respectively are given below:

D0=(Πpμp,0,0,0,0,0)

D∗=(Sp∗(t),Ep∗(t),Ip∗(t),Ap∗(t),Rp∗(t),Mp∗(t)),

where, Sp∗(t)=Πpλ+μp ; Ep∗(t)=λSp∗(t)θpρp−θpωp+μp+ωp; Ip∗(t)=Ep∗(t)(1−θp)ωpμp+τp; Ap∗(t)=Ep∗(t)θpρpμp+τap;

Rp∗(t)=Ap∗(t)τap+Ip∗(t)τpμp; Mp∗(t)=Ap∗(t)wp+Ip∗(t)Qpπ

The solution of the polynomial equation P(λ∗)=m1(λ∗)2+m2λ∗=0 is given as:

λ∗=ηp(ψAp∗(t)+Ip∗(t))Sp∗(t)+(t)+Ip∗(t)+Ap∗(t)+Rp∗(t)+ηwMp∗(t),

where, m1=π(μp+τp)(μp+τap)(θp(ρp−ωp)+μp+ωp),

m2=πμp(μp+τp)(μp+τap)(θp(ρp−ωp)+μp+ωp)(1−R0).

Clearly m1>0 and m2≥0 whenever R0<1, so that λ∗=−m2/m1≤0, thus, no endemic equilibrium exists whenever R0<1.

The essential reproduction quantity (R0) is usually defined as the number of supplementary infections that a certain infectious person would create on the period of the infectious time period provided that everybody else is susceptible. To compute the essential reproductive number (R0) of the Model 4, we follow the work in [33].

A=[0ηpψηpηωΠpμp000000000000]

B=[θpρp+(1−θp)ωp+μp000(θp−1)ωpμp+τp00−θpρp0μp+τap00−Qp−Ωpπ]

The spectral radius γ(AB−1) gives the basic reproduction number of the Model 4.

R0=(μpηpψπ+ΠpΩpηω)θpρp(θp(ρp−ωp)+μp+ωp)πμp(μp+τap)+(μpηpπ+ΠpQpηω)(1−θp)ωp(θp(ρp−ωp)+μp+ωp)πμp(μp+τap)(27)

The threshold value for R0 is 1. If R0<1, the size of secondary infections in human population will not be large enough for persistence of disease and hence the infection will die out over course of time. In the case when R0 > 1, there will be a geometric increase in size of infected population and therefore the infection will continue to spread. The fitted values of transmission rates are presented in Tab. 1.

4  Optimal Homotopy Asymptotic Method (OHAM) Scheme for Fractional COVID-19 Model

Now we develop OHAM scheme for solving underlying fractional order COVID-19 model. OHAM technique is known for its rapid convergence as compared to other techniques. It produces an approximate closed form of the desired solutions and, hence, is known as semi-analytical approach. Procedure of our approximate OHAM scheme has been derived by following the relevant principles as described in literature [34,35]. The underlying COVID-19 fractional order model consists of six governing equations. In the following steps we carry out complete calculations of OHAM scheme for all equations of the model.

Step-1) Homotopies of the governing equations:

We construct the homotopy equations by defining the real valued functions ΦS(t,q),ΦE(t,q),ΦI(t,q),ΦA(t,q),ΦR(t,q) and ΦM(t,q) on the domain Ω×[0,1] representing the approximate solutions of state variables Sp,Ep,Ip,Ap,Rpand M, respectively.

(1−q)(dαΦS(t,q)dtα)−H(q,c)

(dαΦS(t,q)dtα−Πp−μpΦS(t,q)+ηpΦS(t,q)(Ip(t)+ψAp(t))Np(t)+ηωΦS(t,q)M(t))=0(28)

(1−q)(dαΦE(t,q)dtα)−H(q,c)

(dαΦE(t,q)dtα−ηpSp(t)(Ip(t)+ψAp(t))Np(t)−ηωSp(t)M(t)+(1−θp)ωpΦE(t)+θpρpΦE(t)+μpΦE(t))=0(29)

(1−q)(dαΦI(t,q)dtα)−H(q,c)(dαΦI(t,q)dtα−(1−θp)ωpEp(t)+(τp+μp)ΦI(t))=0(30)

(1−q)(dαΦA(t,q)dtα)−H(q,c)(dαΦA(t,q)dtα−θpρpEp(t)+(τap+μp)ΦA(t))=0(31)

(1−q)(dαΦR(t,q)dtα)−H(q,c)(dαΦR(t,q)dtα−τpIp(t)+τapAp(t)+μpΦR(t))=0(32)

(1−q)(dαΦM(t,q)dtα)−H(q,c)(dαΦM(t,q)dtα−QpIp(t)+ΩpAp(t)+πΦM(t))=0(33)

In the above relations Ω=[0,∞), t∈Ω, q∈[0,1] and H(q,c) is an auxiliary function that is always nonzero for all q∈(0,1], whereas H(0,c) is necessarily zero. Auxiliary function in the above relation is defined as under:

H(q,c)=qC1+q2C2+…(34)

Step-2) The steps of OHAM scheme for finding the approximate solution have been carried out in details for the first governing Eq. (28) only. The procedure for the remaining Eqs. (29)–(33) is quite similar. We start by expanding Φs(t,q,Ci) as a Taylor's series about q.

Φs(t,q,Ci)=Sp0(t)+∑k=1∞⁡Spk(t,Ci)qk,i=1,2,…(35)

The convergence of above series depends upon auxiliary constants Ci, i=1,2,…

For q=1 we have

Sp(t,Ci)=Sp0(t)+∑k=1∞⁡Spk(t,Ci),i=1,2,…(36)

Step-3) Comparing the coefficients of like powers of ‘q’ after using the Eq. (35) in Eq. (28) we can obtain approximations of order from zero to onwards, if needed, as below:

q0:Sp0(α)(t)=0(37)

q1:C1Np(t)((−(Np(t)μp+ηpψAp0(t)+ηpIp0(t)+Np(t)ηωM0(t))Sp0(t)+Np(Πp−Spo(α)(t)))++Np(−Spo(α)(t)+Sp1(α)(t)))=0(38)

q2:−1Np(ηpψAp1(t)Sp0(t)C1+ηpIp1(t)Sp0(t)C1+NpηωM1(t)Sp0(t)C1

+NpμpSpt(t)C1+ηpψApo(t)Sp1(t)C1+ηpIpo(t)Sp1(t)C1

+NpηωMoSp1(t)C1−Np∏pC2+NpMpSpo(t)C2

+ηpψAp0(t)Sp0(t)C2+ηpIp0(t)Sp0(t)C2+

NpηωM0(t)Sp0(t)C2+NpC2Sp0(α)(t))−

(1+C1)Sp1(α)(t)+Sp2(α)(t)=0(39)

And so on.

Step-4) Using Reimann-Liouville sense of integral on governing equation of susceptible population and using the initial condition Sp(0), one can obtain following series solution.

Sp0(t)=Sp(0)(40)

Sp1(t)=tα(8065518Np(t)(50000ηω+μp)−Np(t)Πp+16131036ηp(141+100ψ))C1Np(t)αΓ(α)(41)

Sp2(t)=1Np2(t)tα⌈tα(8065518Np(t)(50000ηω+μp)−Np(t)Πp+16131036ηp(141+100ψ))(Np(t)(50000ηω+μp)+ηp(282+200ψ))C12Γ(1+2α)−16131036Np(t)tα[Np(t)ηω(141Qp+100(−250π+Ω))−ηp(141τp+μp(141+100ψ)+100((−100θpρp+τap)ψ+1000(−1+θp)ωp))]C12Γ(1+2α)Np(t)[8065518Np(t)(50000ηω+μp)−Np(t)Πp+16131036ηp(141+100ψ)](C1+C12+C2)Γ(1+α)⌉(42)

where Γ(α) denotes the gamma function. We have the solution in the form:

Sp(t)=Sp0(t)+Sp1(t)+Sp2(t)+…(43)

By substituting the values from Eqs. (40)–(42) in Eq. (43), we obtain second order approximate result of Eq. (21) as:

Sp=8065518+496439.4332413166tαC1αΓ(α)+tα(2.427851710500268×1018tαC12Γ(1+2α)+3.392009602385773×1019(C1+C12+C2)Γ(1+α))68326756000000 (44)

Adopting the same procedure presented in Steps 1–4, we find the following approximate solutions for the state variables Ep(t),Ip(t),Ap(t),Rp(t) and M(t).

Ep=200000+tα34163378000000(−1.231569373491795×1018tαC12Γ(1+2α)−1.695293749139955×1019(C1(2+C1)+C2)Γ(1+α))(45)

Ip=282+2tα(101.2574552024585tαC12Γ(1+2α)−28.001872579217828(C1(2+C1)+C2)Γ(1+α))(46)

Ap=4(50+87.0480383689366t2αC12Γ(1+2α)+11.641883908078782tα(C1(2+C1)+C2)Γ(1+α))(47)

Rp=−34.261698271281816t2αC12Γ(1+2α)−198.69662tα(C1(2+C1)+C2)Γ(1+α)(48)

M=50000+4.972599594940743t2αC12Γ(1+2α)+499.687764tα(C1(2+C1)+C2)Γ(1+α)(49)

Step-5) In this step, the values of auxiliary constants C1 and C2, present in Eqs. (44)--(49), are found by using the method of least squares [33,36,37] by minimizing the total error function of all the governing equations. Let χ={Sp,Ep,Ip,Ap,Rp,M} be the set consisting of all state variables and k be any member of this set. The corresponding residual of each governing equation is denoted by ℜk:k∈χ and total error function is denoted by Jk and is defined as under:

Jk=∫01ℜk2dt:k∈χ

The minimization process of each Jk is explained in following subsections of results and discussions.

5  Results and Discussion

This section is dedicated for presentation of the closed form semi-analytical solutions for all of the state variables and demonstration of their dynamics through graphical exhibition of related simulation results. Model 4 has total six equations. That means COVID-19 fractional model is to analyzed by observing the behaviors of six subpopulations Sp, Ep, Ip, Ap, Rp and M for the various values of α.

The considered values for model parameters are presented by Tab. 1. All figures in the following cases describe the individual behavior of model subpopulations for various values of α. The evolution time t is taken in days by setting a scale of 10 to 1 along horizontal axes.

5.1 Dynamics of Susceptible Human Population

For the fractional analysis of human susceptible population we calculate the auxiliary constants (C1,C2) which describe the order of the solution. We will find these auxiliary constants so that the series solution defined by Eq. (44) satisfies the related governing equation COVID-19 Model 4. The corresponding error function is:

JSp=∫01ℜsp2dt

where, the residual, of susceptible population, denoted by ℜSp, is defined as under:

ℜSp=dαSp(t)dtα−Πp+μpSp(t)+ηpSp(t)(Ip(t)+ψAp(t))Np(t)+ηωSp(t)M(t)

Imposing the necessary conditions for minimum residual, we can find C1 and C2, by solving the equations:

∂JSp∂C1=0 and ∂JSp∂C2=0

Tab. 2 presents the optimized values of auxiliary constants of the second order approximate solution for susceptible population considering various orders of fractional derivative. The corresponding second order approximate solution for Sp has been exhibited in Tab. 3. Fig. 1 demonstrates a clear impact of varying values of the order of the fractional derivative on the dynamics of the susceptible population.

5.2 Dynamics of Exposed Population

The auxiliary constants (C1,C2) describing the order of the solution for exposed population have been calculated first by enforcing the series solution given by Eq. (45) to satisfy the second governing equation of the fractional order COVID-19 model. The total error function is given as:

JEp=∫01ℜEp2dt

where, the residual of exposed population, denoted by ℜEp, is defined as under:

ℜEp=dαEp(t)dtα−ηpSp(t)(Ip(t)+ψAp(t))Np(t)−ηωSp(t)M(t)+

(1−θp)ωpEp(t)+θpρpEp(t)+μpEp(t)

Imposing the necessary conditions for minimum residual, we can find C1 and C2 by solving the equations:

∂JEp∂C1=0 and ∂JEp∂C2=0

Figure 1: Susceptible population vs. time t

Considering various orders of the fractional derivative, above necessary conditions provide the relevant optimum values of auxiliary constants for the exposed population and are presented in Tab. 4.

The OHAM scheme based second order solution for the exposed population Ep has been shown in Tab. 5 for various values of α. Fig. 2 describes the evolution of exposed population with respect to time.