Structure Preserving Algorithm for Fractional Order Mathematical Model of COVID-19

Zafar Iqbal，Muhammad Aziz-ur Rehman，Nauman Ahmed，Ali Raza，Muhammad Rafiq Ilyas Khanand Kottakkaran Sooppy Nisar

1Department of Mathematics，University of Management and Technology，Lahore，Pakistan

2Department of Mathematics and Statistics，The University of Lahore，Lahore，Pakistan

3Stochastic Analysis&Optimization Research Group，Department of Mathematics，Air University，Islamabad，44000，Pakistan

4Department of Mathematics，National College of Business Administration and Economics，Lahore，Pakistan

5Department of Mathematics，F(xiàn)aculty of Sciences，University of Central Punjab，Lahore，54500，Pakistan

6Faculty of Mathematics and Statistics，Ton Duc Thang University，Ho Chi Minh City，72915，Vietnam

7Department of Mathematics，College of Arts and Science at Wadi Aldawaser，Prince Sattam Bin Abdulaziz University，Alkharj，11991，Kingdom of Saudi Arabia

Abstract:In this article，a brief biological structure and some basic properties of COVID-19 are described.A classical integer order model is modified and converted into a fractional order model with ξ as order of the fractional derivative.Moreover， a valued structure preserving the numerical design，coined as Grunwald–Letnikov non-standard finite difference scheme， is developed for the fractional COVID-19 model.Taking into account the importance of the positivity and boundedness of the state variables，some productive results have been proved to ensure these essential features.Stability of the model at a corona free and a corona existing equilibrium points is investigated on the basis of Eigen values.The Routh–Hurwitz criterion is applied for the local stability analysis.An appropriate example with fitted and estimated set of parametric values is presented for the simulations.Graphical solutions are displayed for the chosen values of ξ (fractional order of the derivatives).The role of quarantined policy is also determined gradually to highlight its significance and relevancy in controlling infectious diseases.In the end， outcomes of the study are presented.

Keywords: Coronavirus pandemic model; deterministic ordinary differential equations; numerical methods; convergence analysis

1 Introduction(COVID-19)

Novel coronavirus is a spherical or pleomorphic shaped， particle having single stranded(positive sense) RNA (Ribonucleic Acid) linked to a nucleoprotein surrounded by a special type of protein.The outer surface of the coronavirus contains the projections of the club-shaped structure.The Classification of the coronaviruses depends upon the appearance of the outer surface (whether it is crown like or halo like)， the replication mechanism and the distinct features related to the chemistry of the virus.In general， these viruses belong to OC43-like or 229Elike serotypes.Avian and mammalian species serve as hosts for them.Both types are similar with respect to morphology and chemical structure.Corona viruses present in human beings and animals are antigenically similar.These are capable of attacking on different types of tissues in animals.But， in human beings this family of viruses generally cause only the upper respiratory tract infection.This virus belons to the subclass Orthocoronavirinae， class Coronaviridae， order Nidovirales， realm Riboviria， kingdom Orthcornavirae and phylum Pisuviricota.The dimension of this virus varies from 26 to 32 kilobases， which is largest in the class of RNA viruses.They have distinct protruded club or clove shaped studs or spikes [1].Like other corona viruses， COVID-19 also contains protein in the form of spikes ejecting outside from the surface.These spikes cling with the host (human) cells then its genome bears a structural change and the viral membrane fuse with the host cell cytoplasm.After this step， the viral genes of the COVID-19 enter into the host cell for replication and multiplication of the viruses.Depending upon the protease of the host cell，cleavage reaction permits it to reach into the host cell by endocytosis or fusion.After entering into the host cell， the virus becomes uncovered and their genome attacks on the cell cytoplasm.The genome of the coronavirus works as a messenger and it is translated by the ribosomes of the host cells.These viruses are divided into four categories as alpha coronavirus， beta coronavirus， gamma coronavirus and delta coronavirus.The first two viruses infect the mammals while the last two viruses initially attack the birds.The genera and species of these viruses are described as follows:the species Alphacoronavirus 1， Human coronavirus 229E， Human coronavirus NL63， Miniopterus bat coronavirus 1， Miniopterus bat cor onavirus HKU8， Porcine epidemic diarrhea virus， Rhinolophus bat coronavirus HKU2 and Scotophilus bat coronavirus 512 belong to the Alpha coronavirus.While the species， Betacoronavirus 1(Bovine Coronavirus， Human coronavirus OC43)，Hedgehog coronavirus 1， Human coronavirus HKU1， Middle East respiratory syndrome-related coronavirus， Murine coronavirus， Pipistrellus bat coronavirus HKU5， Rousetlus bat coronavirus HKU9， Severe acute respiratory syndrome- related coronavirus (SARS-Cov， SARS-Cov-2) and Tylonycteris bat coronavirus HKU4 belong to Beta coronavirus.Furthermore， the species Avian coronavirus and Beluga whale coronavirus SW1 are the members of the Gamma coronavirus.Lastly， the Bulbul coronavirus HKU11 and Porcine coronavirus HKU15 species are the family members of the Delta coronavirus.Coronaviruses are deleterious to health with high risk factor.Some of them have more than 30% mortality rate， for instance MERS-COV.But other are not so harmful like as common cold.All types of the coronaviruses can be the causative agent of cold with prime symptoms including fever， sore throat and swollen adenoids.Moreover，they can cause primary viral pneumonia or secondary bacterial pneumonia or bronchitis in the same way as that of pneumonia [2].The SARS-COV appeared in 2003， resulted in severe acute respiratory syndrome (SARS).It effected both the upper and lower respiratory tract due to an unmatched pathogenesis.There are six classes of human coronaviruses that are known so far，each specie is categorized into two types.There are seven types of human coronaviruses.Four coronaviruses which show mild symptoms are:Human coronavirus OC43 (HCOV-OC43)，β-Cov，Human coronavirus HKU1 (HCOV-HKU1) andβ-Cov， Human coronavirus 229E (HCOV-229E)，α-Cov， Human coronavirus NL63 (HCOV-NL63) andα-Cov.Three coronaviruses which show severe symptoms are middle east respiratory syndrome-related coronavirus (MERS-COV)，β-Cov，Severe acute respiratory syndrome coronavirus (SARS-COV)，β-Cov and Severe acute respiratory syndrome coronavirus-2 (SARS-COV-2) andβ-Cov.The HCOV-OC43， HCOV-HKU1， HCOV-229E and HCOV-NL63.They periodically produce the mild symptoms of the common cold in the population all around the year.The outburst of pneumonia in Wuhan， China declared a pandemic by WHO (World Health Organization).This is considered due to a novel type of coronavirus， provisionally named as 2019-nCov by WHO， later it was renamed as SARS-COV-2 by the international committee on Taxonomy of viruses.As of 24 June 2020， 476，911 deaths and more than 9，237691 confirmed cases of COVID-19 were recorded.The Wuhan breed is recognized as a new class of Beta coronavirus， which is genetically similar to SARS-COV.Since COVID-19 has a great resemblance with the bat coronavirus so it is suspected to be initiated from bats also [3].As there is no vaccination or treatment for this disease so it has become a challenge for the scientists， health workers and policy makers to control the spreading of the infection.However， the research community and scientists are making efforts to find the treatment， vaccine or factors that are helpful in slowing down the dynamics of the disease.As a matter of fact，the virus is new however the virus family is not new for the human beings.The humanity has already faced such types of viruses on different scales in the near past.Currently， this infection has frozen all types of academic， business， sports and many other routine activities which created many problems and difficulties for the human beings as well as for the society [4].The eradication of the COVID-19 is an uphill task for the relevant authorities.The modern world is fighting against the infectious diseases on the one hand and changing environmental conditions that are favorable for the emergence of the viral diseases on the other hand.The examination of dead bodies revealed that most of the patients were diagnosed with severe heart， lungs， diabetes and some other diseases.The Disease can be communicated easily through the nasal viral secretion that is transmitted directly or indirectly to the susceptible person.Generally， the symptoms of the infection are mild but in some cases painful death is also observed.The effective forecast for the disease dynamics is a prolific study matter regarding epidemiology， mathematical modeling and simulations.There exist many classes of mathematical models that depend upon the assumptions imposed on the process of dynamics.For instance， SIS， SIR， SEIR， SEIQR and many other compartmental models are used to formulate nonlinear incidence rates and double epidemic hypothesis [5].In these types of modelsS，E，I，QandRdescribe the susceptible， exposed，infected， quarantined and recovered individuals.These cellular or compartmental systems are used for adjusted incidence rate and imperfect vaccinations [6].In these models， it is assumed that the susceptible individuals are aware of the infection’s presence [7].It is worth mentioning that most of the existing deterministic models rest on the ordinary differential equations which imply the assumption of constant diffusion in the domain population.On the other hand， the use of partial differential equations model highlights the non-constant diffusion of the infection [8–10].Although， many studies have been carried out for the providing the deeper insight into the disease dynamics.Chen et al.[11] considered a model with four compartments to investigate the dynamics of the novel infection.Shim and co-authors [12] addressed the questions relating to the effect analysis of the disease.Naveed et al.[13] developed a mathematical model to analyze the virus communication among the population and calculated the basic reproductive number.Many other approaches on the infection propagation may be seen in the sequel.Fractional calculus (FC) is the extension of the integer order calculus.On the basis of FC， the researchers are trying to understand the real world phenomenon of the infectious diseases in a more comprehensive way.They are developing mathematical models with derivatives of non-integer order.By the usage of these types of fractional operators many fruitful studies have been made in the recent scenario.New features and properties of the FC have fascinated the researchers， engineers and scientists to model the problems in the frame work of fractional calculus.The development of new non integer order operators have brought a number of essential features of many physical problems in to the lime light.The history of non-integer order calculus starts with a question posed by Leibniz in 1965.There is a long list of existing fractional order differential operators depending upon the nature of Kernels.Caputo， Riemann–Liouvelle and Katugampola fractional differential operators are developed by using singular kernels.While the fractional operators without singular kernels are of two types.In first type， exponential kernel is used for instance， Caputo–Fabrizio fractional differential operator.Whereas Mittag-Leffler kernel is used in the second type of operators e.g.，Atangana–Baleanu fractional operator in Caputo sense.Due to the salient features of the FC， it is used to model a wide range of physical and dynamical problems in various fields of physical sciences， mathematics， life sciences and engineering [14].Positivity， boundedness and the stability of the equilibrium points for the fractional order physical problems is a challenging task for the scientists and mathematicians.Some researchers are working in this line， but a lot of work is still to be done.Memory effect and hereditary properties are among the celebrated features of the non-integer order derivatives that are helpful in describing the disease dynamics.Taking in to account the memory effect， the fractional order models provide all the important informations from the past that are helpful in forecasting the dynamics of the infection more accurately and comprehensively.Saeedian and co-authors [15] designed a fractional SIR infectious disease model with memory effect and investigated the infection spread in the population.Ucar et al.[16] studied the dynamics of a fractional order smoking model.In the current study， an integer order model is considered initially then we switched on the non-integer order model by using Caputo differential operator.The fractional order model of COVID-19 can describe the complex dynamics of the physical phenomena with a more realistic approach.

2 Preliminaries

In this section， we will present some fundamental definitions of non-integer order derivatives，their key properties and notations used in this article.

2.1 Non Integer Order Derivatives

Fractional order derivatives have been defined by many researchers in a number of ways according to the nature of the kernel used therein.Some basic fractional order operators are defined in this section.Firstly， the Riemann–Liouville non integer order derivative of order 0＜ξ≤1 is defined as，

wherek=[ξ]+1，k-1＜ξ≤k，is thekthorder derivative and Γ (.) is the extensively used gamma function presented by Euler Fractional order derivative of a function is not defined in a unique way， as in the classical calculus an integer order derivative of a function is defined.Fractional derivative is a generic name given to a class of differential operators used to find the non-integer order derivative of a function.These operators are defined by a number of researchers in different ways.Some of them are Liouville， Riemann–Liouville， Caputo， Fabrizio，Grunwald–Letnikov and many more can be studied in the sequel.Among many other questions，it is important that which differential operator best suits for the underlying model.A more generic fractional operator should be selected for the purpose.The Caputo fractional operator is defined as，

The importance of this operator， when applied to solve a system of fractional differential equations is that it can be associated with initial conditions of classical order， which results in an initial value problem in the desired form as，

A very useful definition relating to this article by using the classical finite differences on a uniform mesh partionised in [0，t] is described as follows， consider that) observes the particular smoothness constraints in every interval (0，t)，t≤Twith mesh points as

0=τ0＜τ1＜...＜τn+1=t=(n+1)h

wherehis defined as the difference ofτn+1andτn.By using the classical finite difference symbols，we have

This expression is derived from the famous Euler method.Consider the fractional differential equation

Now， by applying the G-L scheme on a uniform mesh， we obtain the following expression

whereeandγnsatisfy the following relationsand

whereμ，k∈N ∪{0}.

Furthermore，eiandγiobserve the relations as stated in the Lemma 1.

Lemma 1 [17]:Let 0＜ξ ＜1， then the coefficients expressed byare positive and obey the relation

Also the following two relations are satisfied

In this section， we present the GL-NSFD hybrid scheme is formulated by combining the GL scheme for numerical approximation of the fractional order derivatives and NSFD scheme constructed by using the standard rules designed by Mickens.More details can be seen in [18].The system of equations for COVID-19 is described as follows:

2.2 GL-NSFD Scheme

In this portion， we will construct the proposed scheme.The discretization of fractional derivativeC0) is given as，

The above formula is used on the left hand side of Eq.(6) to get the following expression

After some simplifications， we have the final form as，

Similar procedure is adopted for the remaining compartments and we have the final forms as，

2.3 Positivity of the Solution

In this portion， positivity of the solution will be investigated.Positivity is an important feature of the compartmental models.Since， the state variables in these type of models describe the size of the population that cannot be negative.So， positivity is the basic requirement of the solutions at every moment of time.Following result is helpful in this regard.

Theorem:Assume that all the unknowns and parameters arose in the model are non-negative i.e.，So，Eo，Io，QoandRoare positive.Alsoλξ，μξ，kξ，rξ，αξ，d1ξ，d2ξ，β1ξ，β2ξ，q1ξ，q2ξ，qξandφ(h)ξall are ≥0.ThenSn，En，In，QnandRnall are ≥0 ?n∈Z+.

Proof:Taking in to account the Eqs.(11) to (15) forn=0， we have

From the restrictions imposed on the state variables and parameters， it is evident thatS1≥0.Similarly，E1≥0，I1≥0，Q1≥0and R1≥0.Continuing in the same way and by straight forward calculations， it is easy to conclude thatSn+1≥0，En+1≥0，In+1≥0，Qn+1≥0and Rn+1≥0.i.e.，

2.4 Boundedness

Since the state variables in the model represent the subpopulation of a certain compartment.So the sum of values of all the state variables must be less than or equal to the total population or equivalently the sum of solutions at any time must be bounded.The following result is helpful in this regard.

Theorem:LetS0，E0，I0，Q0andR0are all finite quantities andS0+E0+I0+Q0+R0=N0.Moreover， all the parameters involved in the model are positive.i.e.，λξ，μξ，qξ，qξ1，qξ2，αξ，γξ，kξandφ(h)ξare positive， then there is a constantBn+1such thatSn+1≤Bn+1，En+1≤Bn+1，In+1≤Bn+1，Qn+1≤Bn+1，Rn+1≤Bn+1.

Proof:Considering the Eqs.(11) to (15)， we have

By applying the principle of mathematical induction，

forn=0， we have

In the same wayE1≤B1，I1≤B1，Q1≤B1andR1≤B1.

Now， we calculate the expression (16) forn=1， and obtain the following relations，

The above inequalities help us to reach at

Now， let

Sm≤Bm，Em≤Bm，Im≤Bm，Qm≤Bm and Rm≤Bm

For somem∈Z+.

where，

Now forn∈Z+， we arrive at，

In the same fashion a adopted before， we conclude that

So， the given expression is true for all positive values ofn.

Hence， the solutions are bounded ?n∈Z+.

2.5 Stability of the Model

In this portion， we will investigate the stability of the model at both the points of equilibria i.e.， at a corona free equilibrium point and a corona existing equilibrium point.

The corona free equilibrium state of the model is given as

Theorem:The corona free equilibriumC1=(So，Eo，Io，Qo，Ro) =of the model is locally asymptotically stable if Ro＜1， otherwise unstable for Ro＞1.

Proof:The corona-free equilibriumC1=(So，Eo，Io，Qo，Ro) =is locally asymptotically stable (LAS) if all the Eigenvaluesλi ＜0，i=1，2，3，4，5 with condition |arg(λi)|＞For the Eigen values， the Jacobean matrix atis given as follows:

Notice that the two Eigen values are repeated asλ1=-μξ ＜0，λ2=-μξ ＜0 and third Eigen value is

By using the Routh–Hurwitz Criterion of 2ndorder polynomial as，

Hence， all Eigenvalues are negative and by Routh-Hurwitz criteria the given equilibrium pointC1is locally asymptotically stable.

Theorem:The corona existence equilibriumof the model is locally asymptotically stable if Ro＞1， otherwise unstable for Ro＜1.

Proof:The corona existence equilibriumis locally asymptotically stable (LAS) if all the Eigenvaluesλi ＜0，i=1，2，3，4，5 with condition

For the eigen values， the Jacobean matrix atis given as follows:

Notice that， the Eigen values areλ1=-μξ ＜0 andλ2=-(qξ+μξ+0

By using the Routh–Hurwitz Criterion of 3rdorder polynomial， we get the following expression:

and

Thus， we have concluded that all Eigenvalues are negative and by Routh Hurwitz criteria，the given equilibrium pointC2is locally asymptotically stable.Here， we will present a suitable numerical example and graphical solutions of the state variables involved in the model.This whole stuff is presented with the aid of computer simulations.

3 Numerical Example and Simulations

In this portion， an example of the fractional order COVID-19 model is provided.The parametric values are mentioned in Tab.1.Also， non-negative initial conditions are considered.

Computer aided graphs are submitted to support our assertions.These sketches support the fact that proposed numerical device is a structure preserving tool for solving the nonlinear fractional systems.The device encounters the positivity， stability and boundedness of the solutions.All the graphs in Fig.1 reveals that all the subpopulations converge at the virus free equilibrium point (with different values ofξ).E(t)→0，I(t)→0，Q(t)→0 andR(t)→0， when the population is infection free and the values of the parameters are chosen suitably as listed in Tab.1.The graphs in Fig.1 part (a) illustrate that all the curved trajectories representing the susceptible populace growing with timetapproach towards the disease-free value of the susceptible individualswhich is one in this case.Each trajectory is drawn against a certain value ofξ(the order of the fractional derivative) as mentioned in the figure.Moreover，the rate of the convergence towards the VFE of each trajectory is different， depending upon the value ofξ.Similarly， the other sketches in Fig.1 part (b)–part (e) provide the strong evidence for our declaration about the proposed numerical design.The Fig.2 exhibits the simulation results of our proposed scheme endemic equilibrium.The graphs (a)， (b)， (c)， (d) and (e) in Fig.2 provide the graphical solutions toS，E，I，QandRfor some selected values ofξ， whereξis the fractional order of the derivative.The trajectories in Fig.2 part (a) converge at the virus restricted state i.e.， endemic state for different values ofξ， while the other parametric values are kept same as mentioned in Tab.1.Each curved line in the graph Fig.2 part (a) attains the equilibrium state ofS1(t) which is represented as，

with a certain rate of convergence according to the value ofξ.Similarly， the other curved trajectories in Fig.2 part (b)–part (e) depict that they attain the virus endemic equilibrium state for various values ofξ.The endemic equilibrium is expressed as(S1(t)，E1(t)，I1(t)，Q1(t)，R1(t))and the value of each state variable is stated earlier in the section of the stability of the model.The Fig.3 describes the quarantine approach to control the infection in the populace.The values of the quarantine factorqconsidered in this figure are asq1=0.1，q2=0.3，q3=0.5 andq4=0.7 while theξ=0.9 is fixed.All the four curved representations in Fig.3 unveil the key fact that by increasing the quarantine or isolation approach， the flock of infected individuals can be minimized to a certain level.In this area， we will submit the fruitful conclusion about the current study and some future directions will be pointed out.

Table 1:Fitted and estimated parameters values for coronavirus model

Figure 1:The simulations results for all the subpopulations using proposed method at disease free equilibrium with the variation of fractional order ξ

Figure 2:The simulations results for all the subpopulations using proposed method at endemic free equilibrium with the variation of fractional order ξ

Figure 3:The effect of quarantined strategy on infected population by increasing the values of q1 for ξ=0.9

4 Conclusion

In this study， some biological and physical features of the novel corona virus-2019 are described.A classicalSIEQRmodel is converted to fractional order compartmental model withξas order of the fractional derivatives.The GL-NSFD scheme is proposed to study the propagation of the COVID-19 along with some leading properties of the system.Moreover， the numerical study is made to ensure the pre-assumed results about the numerical design.The equilibrium points of the system are also described to detect the local stability of the model.The decisive role ofR0(reproductive number) in describing the stability of the system is also discussed.Positivity and boundedness of the numerical design is also investigated to exhibit the productiveness of the scheme.The computer-aided graphs are presented via computer simulations.These solutions coincide with the exact equilibrium points for different values ofξ.As， the proposed scheme preserves the structure of the system.So， it can be used successfully to solve many other nonlinear physical systems.Moreover， this tool may be used to solve delay models， advection and diffusion reaction models in future.

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.