Mathematical analysis of SIRD model of COVID-19 with Caputo fractional derivative based on real data.

We discuss a fractional-order SIRD mathematical model of the COVID-19 disease in the sense of Caputo in this article. We compute the basic reproduction number through the next-generation matrix. We derive the stability results based on the basic reproduction number. We prove the results of the solution existence and uniqueness via fixed point theory. We utilize the fractional Adams-Bashforth method for obtaining the approximate solution of the proposed model. We illustrate the obtained numerical results in plots to show the COVID-19 transmission dynamics. 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

It is very important to study the mathematical models of infectious diseases for the better understanding of their evaluation, existence, stability and control [1], [2], [3], [4]. As the classical approaches of mathematical models do not determine 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 problems, artificial intelligence, medical diagnoses, robotics, cosmology, and many more. Fractional differential equations are used in mathematical models of biological phenomena during the last few decades [5], [6], [7], [8], [9], [10], [11], [12]. One of the infectious diseases is COVID-19 caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The disease was first reported in December 2019 in Wuhan city, China, which spread throughout the world and led to a continuing pandemic outbreak in 2020 [12]. It is declared that the COVID-19 pandemic represents a great global threat 2019–2020, which has affected 212 countries and different territories around the world. According to the data reported by Worldometer and WHO [13], [14], [15], [16], 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. It is reported that most people infected with COVID-19 will experience mild to moderate lung infections, such as trouble breathing, fatigue, vomiting, cough, and others physical signs. However, other symptoms of varying severity, such as gastroenteritis and neurological disorders have also been identified [15]. If an infected person coughs or sneezes, the COVID-19 transmits mostly by droplets from the nose, and once a person inhales the droplets 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 [16], [17], [18], [19], [20], [21], [22], [23], [24]. It was observed that more and more countries started to put ban on international traveling, close schools, businesses, and even shopping malls. The COVID-19 has caused considerable financial damage around the globe. A significant number of doctors and researchers dedicated themselves to the struggle against pandemics and carried out work in their areas of expertise. They studied COVID-19 from various viewpoints, such as virology, infectious diseases, public environmental occupational safety, microbiology, sociology, veterinary sciences, political economics, and media studies, etc. China, Spain, and the USA are the leading countries on COVID-19 research because the early outbreak of the virus urged them to start relevant research immediately. A group of researchers studied the origin of COVID-19. Recently, Hasan et al. [25] proposed a new compartmental SIRD model of COVID-19 as follow

In this model they divided the total population N into four epidemiological classes: (susceptible class), (infected class), (recovered class), and (death class). The description of the parameters are as follow:

is the average number of contacts per person per time,

is the recovery rate,

is the death rate.

In the last few decades, researchers have given significant attention to the study of fractional calculus. That is because the fractional calculus can more effectively describe and process the preservation and inheritance properties of different structures than the integer-order models. We include [26], [27], [28], [29], [30] for further applications about fractional calculus. Hence, the aforementioned area has been investigated from various angles such as qualitative theory, numerical analysis, etc. (see [28], [29], [30]). Researchers, therefore, expanded the classical calculus to the fractional order via fractional order modeling in [31], [32], [33] using different mathematical techniques. For example, in 1980, Adomian devised a popular method of decomposition to handle nonlinear problems analytically. Before that, the aforementioned method was gradually introduced as an important tool for calculating analytical or approximate solutions to many specie problems. The mathematical models were extensively analyzed using decomposition, homotopy and variation techniques [34], [35], [36], [37]. The methods mentioned were extensively used to treat both linear and nonlinear FODEs [38], [39], [40], [41], [42], [43]. Recently, the “residual power series method”, “Fourier transform method”, “spectral methods” and “method of collocation” as well as several new ways of computational methods have been used to manage differential equations in both fractional and classical order and their systems (see for details [45], [44], [46], [47], [48], [49]).

Since mathematical modeling of dynamical real-world processes/phenomena has been well-predicted by using fractional-order equations. The mentioned operators have appeared in various fields such as physics, chemistry, and engineering. Because in most situations, the memory and hereditary characteristics of materials and processes have not been well-predicted via ordinary differential operators. They can be explained comprehensively with real acceptable remarks by using fractional order operators. Further since the most notable definitions among the fractional differential operators are those given by Reimann-Liouville, Caputo, and recently those given by Caputo-Fabrizo and Atangana-Baleanu type fractional operators. These operators have their characteristics in using. Since in the past two decades, Caputo derivatives have adopted mostly to handle biological models of infectious disease. Therefore by using fractional differential operator, we state that in the fields of continuous-time modeling, various researchers have pointed out that fractional derivatives are better tools in describing acoustics, rheology, polymeric chemistry, linear viscoelasticity, and many more such type of sciences [50]. Therefore, motivated from the above mentioned work, here we study the model (1) under Caputo fractional-order derivative. For ,along with initial conditions

It is remarkable to develop some effective numerical methods to compute the approximate solutions to both linear and nonlinear problems based on fractional-order differential equations. In this regard, various methods and algorithms have been proposed for the analytical and numerical solutions of the mentioned problems. In several research articles, Euler’s method and RK4 method have been utilized to investigate numerical solutions to various classes of fractional order differential equations. Therefore, finding the numerical solution of (2) with an appropriate technique is very important. In this paper, we use the Adams–Bashforth method [51] of fractional order to find the numerical solution of (2). The advantage of this method is that it converts the considered system into a system of Volterra type integral equations. In this way, one can apply the numerical schemes for Volterra type integral equations to find the numerical solution of (2). Keeping these points in view, we first establish some existence results, boundedness, and computation of the basic reproductive numbers for checking the stability of the considered model. Then, by utilizing fixed point theory, we derive some necessary conditions for the existence of at least a solution and its uniqueness.

We organized the paper as follows: Section 1 presents the introduction of the COVID-19 pandemic disease and the development of the fractional calculus. Section 2 includes basic notions of fractional calculus. Section 3 provides equilibrium points and stability analysis. It also deals with the existence and uniqueness results of the solution, and presents the numerical scheme used to find the approximate system solution of the model (2). Section 4 gives the graphical representation of the model, while Section 5 presents the conclusion of the article.

Preliminaries and derivation of basic reproductive number

.

[25] The Caputo fractional derivative of a continuous function on is defined as:where , and represents the integers part of .

[28] The fractional integral of a continuous function on of order corresponding to t is defined as:

The solution the problem with fractional order and where is any real number given by

Here, we discuss the feasibility region of solution and stability analysis of equilibrium points on the basis of the basic reproduction number.

The solution of the model under consideration is restricted to the feasible region given by and the pandemic will occur if .

Letthis implies that

Now adding all equation of (1), we get

Solving (3), we havebut as ,where is initial total population at , where . From this we can write as . This shows the boundedness of the considered population dynamical model.

Next, from first equation of (1) ortherefore, is always decreasing and hence no pandemic will occur. Second equation of (1) iswhere is called threshold phenomenon or critical community size for pandemic.

If which gives . So infection class is decreasing and hence no pandemic will occur. If which yields . Hence, we have , so infection class is increasing and the pandemic will occur.

The disease-free equilibrium point of (2) is

For this, we write (2) as

From second equation of (4) and using (2), we take , thenand putting values of in and equations of (4), we get

Thus the required disease-free equilibrium point is

Hence the theorem is proved.

The basic reproduction number for (2) is determined as

To find the basic reproduction number, we take the second equation of (2) as ,where

The non-linear term is F and the linear term is V. Now we will find the next generation matrix , whereandthen

Since is the leading eigen-value of the next generation matrix at disease-free equilibrium point , which can be written as

Therefore, the basic reproduction number is given by

Thus the requisite result is demonstrated.

If , the pandemic free equilibrium point of (2) is locally asymptotically stable and unstable if .

One can measure the Jacobian matrix for (2) aswhere represents right hand side of (2). Apply partial derivatives, we get

After using the values of , we obtain

Now, the characteristics equation of (7) can be determined as

After evaluation of the determinant our interested eigenvalue is , and is given byor , implies that must be negative if , hence proved our required result.

Main work

We demonstrate the existence of unique solution of the model (2) through fixed point theory. Finally, we present the numerical scheme of the proposed model via fractional Adams–Bashforth method. Let be the Banach space of all continuous and bounded real-valued function on a closed interval . To study the existence of unique solution, we consider (2) as follows:where

Applying fractional integral to both sides of Eq. (8), we have where and .

Now, we define a compact Banach space aswhereand

Also, we assume that

Define the infinite norm as:

Now, we define an operator by

First, we show that the operator is well defined. For this, we need to show that

Considerwhere

Similarly, we can derive the above inequality for the remaining terms. Thus, ifthen the operator is well-defined.

Second, we have to show that the operator satisfies the Lipschitz condition, i.e

For the first component, we havewhere

Similarly for the other terms, we haveso, the operator is a contraction, if

Hence is a contraction on a compact Banach space. It implies that has a unique solution. Consequently, the model (2) has a unique solution.

Numerical scheme

Here the general method to finds the solution of the considered system by fractional Adams–Bashforth method has been discussed. Consider the model (8) with the initial conditions

Applying fractional integral to both sides of (13) and using initial conditions, we get

Here, we will derive the numerical scheme for the first equation of the system (14).

Let us denote , where is the step size and is an integer and . At , the Eq. (15) gets the formand

By subtracting (17) from (16), we obtain

We can write the above equation aswhere

Using Lagrangian polynomial interpolation, we approximate the function asthus,

Solving the integrals on the right hand side of the above equation, we havesimilarly

Substituting the values of and in Eq. (18), we get the approximate solution of the first equation of the model (2).

Similarly for the remaning equations, we have

Numerical simulations

We conduct numerical simulations to represent the effects of our model with data from various world meter and WHO reports; Starting on 4 January 2020, when the Chinese authorities have confirmed 6 Cases a day. Slowly and gradually the number to 1460 on 20 January 2020 followed by 26 deaths. In next day the number confirmed cases increased to 1739 followed by 38 deaths. This number is repidly increased to 3892 with 254 deaths on 4 February 2020 data followd by worldometer. To control the disease Chinese authority lockdown in Wuhan city; the spread of the disease was reduced. As the total population of the Wuhan city is 11 million, therefore take, . In the first days of the outbreak, this denominator was determined and proved later to be a realistic value. It is an appropriate value for the restriction of individual movements according to the actual data recorded by the WHO. As for the preconditions, fixed theses values: .

In Fig. 1, Fig. 2, Fig. 3, Fig. 4 , we have provided the approximate solutions of different classes against the data given and corresponding to various fractional orders. We see that the population are initially assumed to be uninfected (susceptible). Once the epidemic began the susceptible population continued to drop as seen in Fig. 1. Hence, as they were exposed to infection, the population density of the infected class was rapidly increasing as shown in Fig. 2. This rise resulted in an increase in the death class, and many people got rid of the infection that contributed to an increase in the recovered population. Due to the fractional-order the rate of decay and growth is different. The smaller the order the quicker the process concerned, and vice versa. Furthermore, if the fractional-order the subsequent solutions always meet the integer-order solution. Since the fractional differential operator has a greater degree of freedom and offers a full geometry spectrum, we took just a few fractional orders to analyze the model’s dynamic behaviors under consideration. Further, in Fig. 5 , we compare our simulated results with the available real data published in [5] from 4th January 2020 to 8th March 2020 for 67 days as (see Table 1 ).

Fig. 1

Graphical representation of susceptible class at different fractional order.

Fig. 2

Graphical representation of infected class at different fractional order.

Fig. 3

Graphical representation of recovered class at different fractional order.

Fig. 4

Graphical representation of death class at different fractional order.

Fig. 5

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

Table 1

Parameters values for the simulation.

Name	Description	Value	Units
β	Average number of contact per person per time	5	day−−1
γ	Recovery rate	0.5	Dimensionless
κ	Death rate	3.5	day−−1

We see that the graphs of the curves of simulated data and real data are very close to each other at the order of . Hence, is the best suitable fractional-order value. Further, the confirmed reported death in [5] as from 4th January 2020 to 8th March 2020 for 67 days are compare with the simulated data against different fractional order in Fig. 6 .

Fig. 6

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

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

Conclusion and future work

In this paper, we investigated the fractional-order SIRD model of COVID-19 in the case of Caputo. Using the next-generation matrix, we calculated the basic reproduction number and proved the stability conditions based on the basic reproduction number. We established the results of the proposed model’s existence and uniqueness of the solution. We derived a proper numerical scheme using the Adams–Bashforth fractional method to find an approximate solution to the proposed model. Lastly, we presented the numerical simulation to illustrate the better dynamics of the SIRD fractional-order model COVID-19. Also, we have compared our simulated results of infected and death class with 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. We will study the SIRD model in the future under new operators called fractal-fractional operators.

Authors contribution

All the authors have equal contribution in this work.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.