images
Computer Modeling in Engineering & Sciences
images

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

Abstract: This study employs a semi-analytical approach, called Optimal Homotopy Asymptotic Method (OHAM), to analyze a coronavirus (COVID-19) transmission model of fractional order. The proposed method employs Caputo's fractional derivatives and Reimann-Liouville fractional integral sense to solve the underlying model. To the best of our knowledge, this work presents the first application of an optimal homotopy asymptotic scheme for better estimation of the future dynamics of the COVID-19 pandemic. Our proposed fractional-order scheme for the parameterized model is based on the available number of infected cases from January 21 to January 28, 2020, in Wuhan City of China. For the considered real-time data, the basic reproduction number is R0 ≈ 2.48293 that is quite high. The proposed fractional-order scheme for solving the COVID-19 fractional-order model possesses some salient features like producing closed-form semi-analytical solutions, fast convergence and non-dependence on the discretization of the domain. Several graphical presentations have demonstrated the dynamical behaviors of subpopulations involved in the underlying fractional COVID-19 model. The successful application of the scheme presented in this work reveals new horizons of its application to several other fractional-order epidemiological models.

Keywords: Novel COVID-19; semi-analytical scheme; fractional analysis

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.

Figure 2: Exposed population vs. time t

5.3 Dynamics of Infected (Symptomatic) Population

The values of auxiliary constant C1 and C2, present in Eq. (46), have been calculated by minimizing the total residual functional given by:

JIp=∫01⁡ℜIp2dt

where, the residual of infected (symptomatic) population, denoted by ℜIp, is defined as under:

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

Solving the following equations for C1 and C2we get their optimum values that are presented in Tab. 6 for various values of α. The corresponding OHAM scheme based second order approximate solution for infected population is shown in Tab. 7. The graphical view of evolution of infected population under the influence of various values of order of the fractional derivative has been exhibited in Fig. 3.

∂JIp∂C1=0 and ∂JIp∂C2=0

Figure 3: Infected population vs. time t

5.4 Dynamics of Asymptotically Infected Population

The values of auxiliary constant C1 and C2 involved in Eq. (47) of asymptotically infected population have been calculated as follows:

JAp=∫01⁡ℜAp2dt

where, the residual of asymptotically infected population, denoted by ℜAp, is defined as under:

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

The values of C1andC2 presented in Tab. 8 have been obtained by solving following equations:

∂JAp∂C1=0 and ∂JAp∂C2=0

Tab. 9 displays the second order solution for Ap. Fig. 4 demonstrates the evolutionary behavior of asymptomatically infected population over time.

Figure 4: Asymptomatically infected population vs. time t

5.5 Dynamics of the Recovered or the Removed Population

The following total error function corresponding to recovered population is obtained from Eq. (48):

JRp=∫01⁡ℜRp2dt

The residual ℜRp, for the recovered or the removed population, is defined given as:

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

Solving the following necessary conditions we obtain the values of C1 and C2 presented in Tab. 10.

∂JRp∂C1=0 and ∂JRp∂C2=0

The closed form solution for recovered population is presented in Tab. 11 using different values of α and the relevant graphs are exhibited in Fig. 5.

Figure 5: Recovered population vs. time t

5.6 Dynamics of the Reservoir Compartment

Applying the least square approach for the auxiliary constant of Eq. (49) we get following total error function:

JM=∫01⁡ℜM2dt

The residual ℜM of the class M is defined below:

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

The solutions of following equations give the values of C1 and C2 shown in Tab. 12:

∂JM∂C1=0 and ∂JM∂C2=0

The approximate solutions for M have been presented in Tab. 13 whereas the impact of variation in order of the fractional derivative has been shown in Fig. 6.

Figure 6: Reservoir saturation vs. time t

Conclusion

In this study, closed form semi-analytical approximate solution has been presented for fractional order COVID-19 model in Caputo sense of derivative operator. The underlying model involves parameters that were extracted by fitting real data into the dynamical model. The fitted parameters are responsible for a high reproduction number R0≈2.48293 that indicates the fact that the infection will persist. The sensitivity of R0 reveals that the infection can be controlled by reducing the contact rates ηp and ψ of susceptible population to symptomatically and asymptomatically infections, the contribution rates Ωp and Qp of infected populations to environment and the reservoir's disease transmissibility rate ηω. Such reductions can be directly achieved by implementing more quarantining, hospitalization and lockdown strategies by policymakers. Moreover, the availability of a healthy environment to increase population immunity and effective remedies will increase the recovery rates τap and τa which will result in a reduced reproductive number. Taking into account of available data and the fitted parameters, the dynamics of human subpopulations have been efficiently examined by considering various values of order of fraction derivatives. Analytical and graphical results demonstrate that an increase in the order α of fractional derivative shows increase in the exposed, symptomatically infected, asymptomatically infected and recovered subpopulations sizes whereas the number of susceptible individuals and environmental saturation decrease. The results obtained are promising providing us great inspiration to use our proposed methodology further to explore more advanced models of COVID-19 disease, especially, delayed and stochastic models. It is expected that this latitude study would provide researchers a good viewpoint of more advanced research that could be put forward on closed form solutions of epidemiological models.

Acknowledgement: The authors are thankful to their respective departments and Universities.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Macrotrends (2020). Wuhan, China Population 1950–2020. http://www.macrotrends.net/cities/20712/wuhan/population.

 2.  The New York Times (2020). Is the World Ready for the Coronavirus? https://www.cnbc.com/2020/01/24/chinas-hubei-province-confirms-15-more-deaths-due-to-coronavirus.html.

 3.  CNBC (2020). China Virus Death Toll Rises to 41, More than 1,300 Infected Worldwide. https://www.cnbc.com/2020/01/24/chinas-hubei-province-confirms-15-more-deaths-due-to-coronavirus.html.

 4.  COVID-19 Coronavirus updates (2020). https://www.worldometers.info/coronavirus/.

 5.  Cao, J. L., Tu, W. J., Hu, X. R., Liu, Q. (2020). Clinical features and short-term outcomes of 102 patients with coronavirus disease 2019 in Wuhan, China. Clinical Infectious Diseases, 71(15), 748–755. DOI 10.1093/cid/ciaa243.

 6.  Cao, J. L., Hu, X. R., Tu, W. J., Liu, Q. (2020). Clinical features and short-term outcomes of 18 patients with corona virus disease 2019 in intensive care unit. Intensive Care Medicine, 46(5), 851–853. DOI 10.1007/s00134-020-05987-7.

 7.  Nesteruk, I. (2020). Statistics-based predictions of coronavirus epidemic spreading in mainland China. Innovative Biosystems and Bioengineering, 4(1), 13–18. DOI 10.20535/ibb.2020.4.1.195074.

 8.  Ming, W., Huang, J. V., Zhang, C. J. P. (2020). Breaking down of the healthcare system: Mathematical modelling for controlling the novel coronavirus (2019-nCoV) outbreak in Wuhan, China. bioRxiv. DOI 10.20535/ibb.2020.4.1.195074.

 9.  Zeb, A., Alzahrani, E., Erturk, V. S., Zaman, G. (2020). Mathematical model for coronavirus disease 2019 (COVID-19) containing isolation class. BioMed Research International, 2020, 7. DOI 10.1155/2020/3452402.

10. Batista, M. (2020). Estimation of the final size of the coronavirus epidemic by SIR model. https://www.researchgate.net/publication/339311383.

11. Victor, A. (2020). Mathematical predictions for COVID-19 as a global pandemic. DOI 10.2139/ssrn.3555879.

12. Sarwar, S., Zahid, M. A., Iqbal, S. (2016). Mathematical study of fractional-order biological population model using optimal homotopy asymptotic method. International Journal of Biomathematics, 9(6), 1650081. DOI 10.1142/S1793524516500819.

13. Sarwar, S., Iqbal, S. (2017). Stability analysis, dynamical behavior and analytical solutions of nonlinear fractional differential system arising in chemical reaction. Chinese Journal of Physics, 56(1), 374–384. DOI 10.1016/j.cjph.2017.11.009.

14. Iqbal, S., Mufti, M. R., Afzal, H., Sarwar, S. (2019). Semi analytical solutions for fractional order singular partial differential equations with variable coefficients. AIP Conference Proceedings, 2116(1), 3000071–3000074. DOI 10.1063/1.5114307.

15. Sarwar, S., Iqbal, S. (2017). Exact solutions of the non-linear fractional klein-gordon equation using the optimal homotopy asymptotic method. Nonlinear Science Letters A, 8(4), 340–348.

16. Sarwar, S., Alkhalaf, S., Iqbal, S., Zahid, M. A. (2015). A note on optimal homotopy asymptotic method for the solutions of fractional order heat and wave-like partial differential equations. Computers & Mathematics with Applications, 7(5), 942–953. DOI 10.1016/j.camwa.2015.06.017.

17. Friehet, A., Hasan, S., Al-Smadi, M., Gaith, M., Momani, S. (2019). Construction of fractional power series solutions to fractional stiff system using residual functions algorithm. Advances in Difference Equations, 2019, 95. DOI 10.1186/s13662-019-2042-3.

18. Oldham, K. B., Spanier, J. (1974). The fractional calculus. New York, USA: Academic Press.

19. Miller, K. S., Ross, B. (1993). An introduction to the fractional calculus and fractional differential equations. New York, USA: Wiley-Interscience.

20. Podlubny, I. (1999). Fractional differential equations. New York, USA: Academic Press.

21. Younas, H. M., Mustahsan, M., Manzoor, T., Salamat, N., Iqbal, S. (2019). Dynamical study of fokker-Planck equations by using optimal homotopy asymptotic method. Mathematics, 7(3), 264. DOI 10.3390/math7030264.

22. Mustahsan, M., Younas, H. M., Iqbal, S., Nisar, K. S., Singh, J. (2020). An efficient analytical technique for time-fractional parabolic partial differential equations. Frontiers in Physics, 8, 131. DOI 10.3389/fphy.2020.00131.

23. Blake, T. D., Ruschak, K. J., Kistler, S. F., Schweizer, P. M. (1997). Liquid film coating. London, UK: Chapman & Hall.

24. Kumar, S., Atangana, A. (2020). A numerical study of the nonlinear fractional mathematical model of tumor cells in presence of chemotherapeutic treatment. International Journal of Biomathematics, 13(3), 2050021. DOI 10.1142/S1793524520500217.

25. Pandey, P., Kumar, S., Das, S. (2019). Approximate analytical solution of coupled fractional order reaction-advection-diffusion equations. The European Physical Journal Plus, 134(7), 364. DOI 10.1140/epjp/i2019-12727-6.

26. Kumar, S. (2020). Numerical solution of fuzzy fractional diffusion equation by chebyshev spectral method. Numerical Methods for Partial Differential Equations, 2020, 1–19. DOI 10.1002/num.22650.

27. Kumar, S., Gómez Aguilar, J. F., Pandey, P. (2020). Numerical solutions for the reaction-diffusion, diffusion-wave and cattaneo equations using a new operational matrix for the caputo-fabrizio derivative. Mathematical Methods in the Applied Sciences, 43(15), 8595–8607. DOI 10.1002/mma.6517.

28. Kumar, S., Pandey, P., Das, S. (2020). Operational matrix method for solving nonlinear space-time fractional order reaction-diffusion equation based on genocchi polynomial. Special Topics & Reviews in Porous Media: an International Journal, 11(1), 33–47. DOI 10.1615/SpecialTopicsRevPorousMedia.2020030750.

29. Kumar, S., Pandey, P. (2020). Quasi wavelet numerical approach of non-linear reaction diffusion and integro reaction-diffusion equation with atangana–Baleanu time fractional derivative. Chaos, Solitons & Fractals, 130, 109456. DOI 10.1016/j.chaos.2019.109456.

30. Kumar, S., Cao, J., Abdel-Aty, M. (2020). A novel mathematical approach of COVID-19 with non-singular fractional derivative. Chaos, Solitons & Fractals, 139, 110048. DOI 10.1016/j.chaos.2020.110048.

31. Pandey, P., Kumar, S., Jafari, H., Das, S. (2019). An operational matrix for solving time-fractional order cahn-hilliard equation. Thermal Science, 23, 369–369. DOI 10.2298/TSCI190725369P.

32. Ali, Z., Rabiei, F., Shah, K., Khodadadi, T. (2020). Modeling and analysis of novel COVID-19 under fractal-fractional derivative with case study of Malaysia. Fractals, 29(1). DOI 10.1142/S0218348X21500201.

33. Ali, Z., Rabiei, F., Shah, K., Khodadadi, T. (2020). Qualitative analysis of fractal-fractional order COVID-19 mathematical model with case study of wuhan. Alexandria Engineering Journal, 60(1), 477–489. DOI 10.1016/j.aej.2020.09.020.

34. Iqbal, S., Idrees, M., Siddiqui, A. M., Ansari, A. R. (2010). Some solutions of linear and nonlinear klein-gordon equations using the optimal homotopy asymptotic method. Applied Mathematics and Computation, 216(10), 2898–2909. DOI 10.1016/j.amc.2010.04.001.

35. Iqbal, S., Javed, A. (2011). Application of optimal homotopy asymptotic method for the analytic solution of singular lane-emden type equation. Applied Mathematics and Computation, 217(9), 7753–7761. DOI 10.1016/j.amc.2011.02.083.

36. Ruschak, K. J. (1985). Coating flows. Annual Review of Fluid Mechanics, 17, 65–89. DOI 10.1146/annurev.fluid.36.050802.122049.

37. Gaskell, P. H., Savage, M. D., Summers, J. L. (1995). The mechanics of thin film coatings. Proceedings of the First European Coating Symposium, pp. 19–22. Leeds University, Leeds, UK.

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.