Fractional order mathematical modeling of COVID-19 transmission.

In this article, the mathematical model with different compartments for the transmission dynamics of coronavirus-19 disease (COVID-19) is presented under the fractional-order derivative. Some results regarding the existence of at least one solution through fixed point results are derived. Then for the concerned approximate solution, the modified Euler method for fractional-order differential equations (FODEs) is utilized. Initially, we simulate the results by using some available data for different fractional-order to show the appropriateness of the proposed method. Further, we compare our results with some reported real data against confirmed infected and death cases per day for the initial 67 days in Wuhan city.

Introduction

Coronavirus disease 2019 (COVID-19) is one of the infectious diseases caused by Coronavirus-2 SARS-CoV-2, a serious acute respiratory syndrome. The disease was first detected in Wuhan city, China, in December 2019 and has since spread globally, leading to a continuing pandemic outbreak in 2020 [9]. It is declared that the COVID-19 pandemic is the biggest global threat in 2020 which has affected 212 countries and territories around the world. According to the data reported by Worldometer [9] and WHO (World Health Organization) [10], [11], as of May 03, 2020, it has been noticed that more than 3.5 million people were infected with 0.247 million deaths. Even in some countries like Italy and Spain, the death rate is as high as almost 0.066. This verifies the severity and high infectivity of 2019-nCoV. It is confirmed that most people infected with 2019-nCoV will experience mild to moderate respiratory illness, such as breathing difficulty, fever, sickness, cough, and other symptoms. However, other symptoms such as gastroenteritis and neurological diseases of varying severity have also been reported by [12], [31], [32], [33]. The 2019-nCoV transmits mainly through droplets from the nose when an infected person coughs or sneezes. Once one person inhales the droplets from infected people in the air, he will be exposed to the danger of getting the infection. As a result, the best way to prevent the virus is to avoid mixing up with the people. The severity of this pandemic attracted the researchers and scientists throughout the world [13], [14]. It was observed that more and more countries started to ban international traveling, close schools, shopping malls, and companies. The 2019-nCoV pandemic has led to serious economic damage in the whole world. A large number of doctors and researchers also devoted themselves to the anti-pandemic war and conducted researches in their areas of expertise. They looked into 2019-nCoV from various points of view, such as virology, infectious diseases, microbiology, public environmental occupational health, veterinary sciences, and sociology, media studies, political economics, etc. China, USA, and Korea are the leading countries on the 2019-nCoV research because the early outbreak of the virus urged them to start relevant research immediately. A group of researchers studied the origin of 2019-nCoV. It is noted that by adding the class of super-spreaders a number of free equilibrium points for any compartmental model were provided for analysis of the disease [15], [16]. Ndaȯrou et al., amalgamate the models proposed by Kim et al. [17] and Alasmawi et al. [18], and suggested a new epidemiological compartment model that would take into account the super-spreading phenomenon of some individuals [19]. They further take a death-related compartment due to the virus infection. By doing this they adopted the new model as given below:

In this model they divided the total population N into eight epidemiological classes: S susceptible class, E exposed class, I symptomatic and infectious class, P super-spreaders class, A infectious but asymptomatic class, H hospitalized, R recovery class, and F fatality class. By quantifying the transmission coefficient of human-to-human, per unit time per person, where, 0 and 1, quantifies the said transmission coefficient of super-spreaders and hospitalized patients respectively. The rate at which an infectious person becomes symptomatic, super-spreader, or asymptomatic leaves the exposed class. The exposed individuals become super-spreaders at a very low rate. Individuals belonging to the symptomatic and super-spreaders classes are hospitalized at the rate of a; i and r are the recovery rate of hospitalized, and without being hospitalized patient. The death rates induced infected disease are i, p, and h of super spreaders, and individuals hospitalized, respectively. At any moment in time,die due to the disease. This model the transmissibility from asymptomatic individuals, as their behavior was not evident. This problem is still controversial for epidemiologists at present.

It is very important to study the mathematical models of infectious diseases for a better understanding of their evaluation, existence, stability, and control [1], [2], [3], [37], [40], [41]. As the classical approaches of mathematical models do not determine the high degree of accuracy to model these diseases, fractional differential equations were introduced to handle such problems, which have many applications in applied fields like production problems, optimization problem, artificial intelligence, medical diagnoses, robotics, cosmology and many more. In the last few decades, the fractional differential has been used in mathematical modeling of biological phenomena [4], [5], [6], [7], [8], [38], [39]. This is because fractional calculus can explain and process the retention and heritage properties of various materials more accurately than integer-order models. We include [20], [21], [22], [23] for further applications about fractional calculus. Hence, the aforementioned area has been investigated from various angles such as qualitative theory, numerical analysis, etc. (see [24], [25], [26]). Researchers, therefore, expanded the classical calculus to the fractional-order via fractional-order modeling in ([27], [28], [29], [30]) using different mathematical techniques. Since mathematical models are powerful tools to investigate infectious disease. The mentioned area has been explored very well. Recently some authors have considered mathematical models of COVID-19 under fractional order derivatives and produced very good results (see some [34], [35], [36]).

Therefore, motivated from the above mentioned work, here we study the model (1) under Caputo fractional derivative of order γ under initial conditions,For the proposed model, we first derive existence results by using fixed point theory. Then, we extend the famous modified Euler method for numerical simulations. The concerned method is a powerful technique for the computation of numerical results. Hence, we first simulate the results against the available data taken from [2]. Then, we compared simulated data at different fractional order with real data.

We organized the article as follows:

Section 1 presents the introduction to the pandemic disease COVID-19, mathematical models, importance, and development of fractional calculus and fractional differential equations. In Section 2 basic definitions of fractional calculus are given. Section 3 deals with the existence and uniqueness of the proposed model and provides the proper procedure of finding the general solution of the considered model by using the modified Euler’s method. Section 4 describes the graphical representation of the proposed model. In Section 5, we present the conclusion of the manuscript.

Preliminaries

Podlubny [26] Let be a continuous function on a fractional integral in Riemann-Liouville sense corresponding to t is defined as:

Kilbas et al. [20] Let be a continuous function on [0, T]. The Caputo fractional derivative may be expressed as:where and ⌊β⌋ represents the integers part of β.

Main work

In this section, we will discuss the existence and uniqueness of the proposed model. Also, by using fractional Euler’s method we will derive the numerical solution of (2).

Equilibrium points and stability analysis

To determine the equilibrium points of the proposed model equate the right hand side of (2) to zero.The disease-free equilibrium points is given byThe endemic equilibrium points is denoted by and can be obtained by solving the system of equations in (3) simultaneously, taking into the consideration the fact that

Since, the dimension of the proposed model is much higher, therefore, it is very difficult to discuss the stability of the proposed model concerning equilibrium points. Thus, we will study the stability analysis based on the basic reproduction number. The basic reproduction number is computed by using the next-generation matrix approach. For more detail see [19] where and .

The diseases free equilibrium of system(2), i.e. (N, 0, 0, 0, 0, 0, 0) is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

The proof of the theorem follows from Ndarou et al. [19]. □

Existence and uniqueness

Here we will discuss the existence and uniqueness of the considered model. Applying fractional integral to (2) and using initial conditions we obtained,where the functions under the integral signs in (5) are defined as:Assume that S(t), E(t), I(t), P(t), A(t), H(t), R(t) and F(t) are non-negative bounded functions. i.e, there exists some positive constants such that

The functions Ksatisfy Lipschitz’s conditions and are contraction mappings, if the following condition holds,

First we consider the function K 1. For any S and S 1 consider where . Thus, K 1satisfy Lipschitz condition. Similarly, it can be shown that, we can find ς, for so that K for the Lipschitz’s conditions are satisfied. Moreover, under the conditionthe functions are contractions. □

Now, we can write (5) recursively as:

The initial components of the above equations are determined by the given initial conditions. The difference between two terms can be represented as :whereConsiderIn similar fashion, we can obtain

The functions defined in(9)are exist and smooth.

If there exist t0 > 1 such that for then, at least one solution of the system exist.

Since, the functions S(t), E(t), I(t), P(t), A(t), H(t), R(t) and F(t) are bounded and each kernels K for fulfill Lipschitz’s conditions, thus, we obtain the following relations:

The system (10) shows the existence and smoothness of the function defined in (9).

We will show that S(t), E(t), I(t), P(t), A(t), H(t), R(t) and F(t) converge to system of solutions of (2).

Define, B(t), C(t), D(t), G(t), L(t), M(t), N(t), O(t), as remainder terms after n-iterations, such thatUsing triangle inequality along with the Lipschitz condition of K 1, we obtain:Applying the above process recursively, we getThen, at t 0 one hasTaking limit as n tends to infinityUsing hypothesis Eq. (12) becomesSimilarly, on using limit as n tends to infinity, we obtainThus, at least one solution of the system exist.

□

If the conditionforthen the system(2)has a unique solution.

Assume that {S 1(t), E 1(t), I 1(t), P 1(t), A 1(t), H 1(t), R 1(t), F 1(t)} is another set of solution of system (2) then,Rearranging the terms, we getusing hypothesis the last Eq. (13) gets the formIt means that . Applying the same procedure to each solution for we obtainThus, the theorem is proved. □

Procedure for solution

Here we derive the general procedure of fractional order Euler method for our considered model. Reformulate (2) as follows:where

To procure an iterative scheme, we go ahead with first equation of the model (2) as follows:

Let [0, b] be the set of points which we want to find the solution of the Eq. (16). Actually, we cannot evaluate the function S(t) which will be the solution of Eq. (16). Instead of this, a set points is generated from which the points are used for our iterative procedure. For this, we divide the interval [0, b] into k subintervals of equal width using the nodes for Assume that and are continous on [0, b]. By generalized Taylor formula expand S(t) about For each value t there is a value C 1 so thatwhen and are substituted into (17) the result is an expression for

If the step size h is chosen small enough, then we may neglect the second-order term (h 2) and getOn repeating the same fashion, a sequence of points that approximates the solution is formed. A general formula about is Eq. (19) represents the general formula for fractional Euler’s method. Now, we shall derive the fundamental algorithm for the numerical solution of the Eq. (16). Applying fractional integral to both sides of (16), we haveTo obtain the solution point (t 1, S(t 1)), we substitute into (20)and we getNow if the modified trapezoidal rule is used to approximate with the step size then Eq. (21)becomesNotice that the formula on the right hand side of Eq. (22)involves the term S(t 1). So, we use an estimate for S(t 1). Fractional Euler’s method will suffice for this purpose. Substituting (18) into (22), yieldsThe process is repeated to generate a sequence of points that approximate the solution S(t). The general formula for our algorithm is

Using the same procedure, we obtain the numerical scheme for the other compartments of the model (2)

Graphical presentations and discussion

In the present section, we consider the model (2) and using the parameters values shown in Table 1 to obtain the graphical results. We conduct numerical simulations to equate the proposed model results with the real data obtained from various reports from Worldometer and WHO, started from Jan 4, 2020, when the Chines authorities confirmed 6 cases in one day. Slowly and gradually the number rose to 1460 on 20th January followed by 26 deaths. On the next day, the number of confirmed cases increased to 1739 followed by 38 deaths. This number rapidly increased to 3892 with 254 deaths on 4th February 2020 according to Worldometer. To control the disease, Chines authorities put lockdown in Wuhan city. Therefore, the spread of the disease was reduced. As the total population of Wuhan city is 11 million take, . This denominator was measured in the first days of the outbreak and was later proven to be a reasonable value; it is an acceptable value for the restriction of individual movements according to the actual data reported by the WHO. As for the preconditions, fix theses values: .

Table 1

Parameter values for the numerical simulations of the proposed model [2].

Name	Description	Value	Units
β	“Transmission coefficient from infected individuals”	2.55001	perday
l	“Relative transmissibility of hospitalized patients”	1.56001	dimensionless
β′	“Transmission coefficient due to super-spreaders”	7.5001	perday
k^^′	“Rate at which exposed become infectious”	0.25001	perday
ρ1	“Rate at which exposed people become infected I”	0.58001	dimensionless
ρ2	“Rate at which exposed people become super-spreaders”	0.00100	dimensionless
γα	“Rate of being hospitalized”	0.94001	perday
γi	“Recovery rate without being hospitalized”	0.27001	perday
γr	“Recovery rate of hospitalized patients”	0.50001	perday
δi	“Disease induced death rate due to infected class”	3.50001	perday
δp	“Disease induced death rate due to super-spreaders”	1.00001	perday
δh	“Disease induced death rate due to hospitalized class”	0.30001	perday

The approximate solutions given in (23)-(30) are presented by graphs corresponding to different fractional order as:

In Figs. 1 –8, we have presented the approximate solutions of different compartments against the given data and corresponding to different fractional orders. We see that initially the people are assumed uninfected (susceptible). When the outbreak started the population of susceptible was going on decreasing as in Fig. 1. Since they were exposed to infection, therefore the population density of exposed, symptomatic and infectious, super-spreaders, infectious, hospitalized classes were increasing rapidly as presented in Figs. 2 –6 . This increase resulted in the rise of death rate, and many people got rid of the infection which led to an increase in the population of recovered and fatality classes. The rate of decay and growth is different due to fractional order. The smaller the order the faster the concerned process and vice versa. Hence, when the fractional order γ → 0, the corresponding solutions also approach to the solution at integer order. As fractional differential operator has a greater degree of freedom which provides a complete spectrum of the geometry, we have taken only a few fractional orders to observe the dynamical behaviors of the model under consideration. Further, in Fig. 9 , we compare our simulated results with the available real data published in [9] from 4th January 2020 to 8th March 2020 for 67 days as [6, 12, 19, 25, 31, 38, 44, 60, 80, 131, 131, 259, 467, 688, 776, 1776, 1460, 1739, 1984, 2101, 2590, 2827, 3233, 3892, 3697, 3151, 3387, 2653, 2984, 2473, 2022, 1820, 1998, 1506, 1278, 2051, 1772, 1891, 399, 894, 397, 650, 415, 518, 412, 439, 441, 435, 579, 206, 130, 120, 143, 146, 102, 46, 45, 20, 31, 26, 11, 18, 27, 29, 39, 39].

Fig. 1

Dynamical behavior of susceptible class at various fractional order of the considered model.

Fig. 8

Dynamical behavior of fatality class at various fractional order of the considered model.

Fig. 2

Dynamical behavior of exposed class at various fractional order of the considered model.

Fig. 6

Dynamical behavior of hospitalized class at various fractional order of the considered model.

Fig. 9

Comparison of simulated and real data at different fractional order for the confirmed reported cases per day of the proposed model.

Dynamical behavior of symptomatic and infectious class at various fractional order of the considered model.

Dynamical behavior of super-spreaders class at various fractional order of the considered model.

Dynamical behavior of infectious but asymptomatic class at various fractional order of the considered model.

Dynamical behavior of Recovered class at various fractional order of the considered model.

We see that the graphs of the curves of simulated data and real data are very close to each other at the order of 0.97. Hence, is the best suitable fractional-order value. Further, the confirmed reported death in [9] as [0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 4, 8, 15, 15, 25, 26, 26, 38, 43, 46, 45, 57, 64, 66, 73, 73, 86, 89, 97, 108, 97, 254, 121, 121, 142, 106, 106, 98, 115, 118, 109, 97, 150, 71, 52, 29, 44, 37, 35, 42, 31, 38, 31, 30, 28, 27, 23, 17, 22, 11, 7, 14, 10, 14, 13, 13] from 4th January 2020 to 8th March 2020 for 67 days are compare with the simulated data against different fractional order in Fig. 10 .

Fig. 10

Comparison of simulated and real data at different fractional order for the confirmed reported death per day of the proposed model.

Again from Fig. 10, we see that at fractional-order the simulated data and real are very close to each other; hence, the best choice of the fractional-order value is .

Conclusion

This paper investigates the fractional-order mathematical modeling of COVID-19 transmission. We show the existence and uniqueness of the model by using nonlinear analysis. We obtain the numerical solution of the proposed model through fractional Euler’s method. In the end, numerical simulation is also provided. We observe that fractional differential operators provide the global dynamics of the model we have considered. It is also observed that the smaller the fractional-order faster the decay or growth process and hence stability occurs rapidly on smaller fractional orders. Further, the results we have simulated in various Figs. 1–8 for different fractional order. As the order is increasing the solution approaches to the result at integer order 1. Also, we have compared our simulated results at different fractional-order against real data for reported cases of infection and death. We see that numerical results are close to real data solutions. The curve very well coincides with real data at . So this is the best value of fractional order.

Authors contribution

All the authors have equal contribution in this work.

Funding

No source exist of funding this work.

Declaration of Competing Interest

None.