Complex mathematical SIR model for spreading of COVID-19 virus with Mittag-Leffler kernel.

This paper investigates a new model on coronavirus-19 disease (COVID-19), that is complex fractional SIR epidemic model with a nonstandard nonlinear incidence rate and a recovery, where derivative operator with Mittag-Leffler kernel in the Caputo sense (ABC). The model has two equilibrium points when the basic reproduction number R 0 > 1 ; a disease-free equilibrium E 0 and a disease endemic equilibrium E 1 . The disease-free equilibrium stage is locally and globally asymptotically stable when the basic reproduction number R 0 < 1 , we show that the endemic equilibrium state is locally asymptotically stable if R 0 > 1 . We also prove the existence and uniqueness of the solution for the Atangana-Baleanu SIR model by using a fixed-point method. Since the Atangana-Baleanu fractional derivative gives better precise results to the derivative with exponential kernel because of having fractional order, hence, it is a generalized form of the derivative with exponential kernel. The numerical simulations are explored for various values of the fractional order. Finally, the effect of the ABC fractional-order derivative on suspected and infected individuals carefully is examined and compared with the real data.

Introduction

The coronavirus is a new viral pneumonia, deadly and rapidly spreading infection that has put great panic around the globe since the break out late 2019. Within a short span of time, it affected every country all over the world, which urges every nation take severe action to control the wild spread of the virus as the severity of the disease will harm human life badly.

It is well known that infectious diseases are a massive threat for humans and also for national economies. Proper understanding of a disease dynamics is necessary, which plays a major role for controlling and eventually exterminate the infection in a community. Implementation of the correct applicable strategy against the disease transmission is another challenge. Mathematical modeling of the infectious diseases is one of the key tools in order to handle these challenges. A number of disease models have been established in the past, which enables us to investigate and control the spread of infectious diseases in a better way [1-3].

Regarding COVID-19, from the initial spread of COVID-19, many of the known mathematical models were used, or newly established by many researchers, For example Giordano et al. in [4] proposed the so called SIDARTHE model and concluded that restrictive social-distancing measures will need to be combined with widespread testing and contact tracing to end the ongoing COVID-19 pandemic. At a same time, Alshammari [5] proposed the SEYNHR compartmental model and applied real data for Saudi Arabia and obtained good predictions on a short term. The role of quarantining and isolation to control the spread of COVID-19 is studied by Memon et al. [6], where they used an extended SIR model. They show that quarantine and isolation are effective measures to control the outbreak.

From real data of the infectious disease, we know that the outbreak of the disease within the country or state for time is generally nonlinear, which tells us to design the system where we can study those dynamic nonlinear phenomena. By this system, we can be able to define the transmission of such a contagious disease, which helps us to interpret the remedial measures to stop or contain the spread of the contagious disease.

Most of the present day studies on mathematical modeling of infectious disease are based on integer-order differential equations (IDEs). However, more recently, many authors state that fractional-order models (FDEs) can describe natural phenomena more precisely than the integer-order differential equations because fractional-order derivatives and integrals describe the memory and hereditary properties of different phenomena [7]. Further, the classical IDEs cannot provide the result for non-integer values. In order to overcome the above limitation of integer-order derivatives, different types of fractional-order operators were defined, some of the research can be found in the existing literature [8] and applications of these fractional-order operators are noted in [9], also Sene in [10, 11] consider the Liu et al. [12] chaotic system with Caputo fractional derivative and show the Lyapunov exponent to characterize the nature of chaos and prove the dissipativity of the considered chaotic system.

However, it is known that the classical form of the fractional derivatives with singular kernel may not be suitable to characterize the non-local dynamics in an appropriate manner. In order to overcome this problem, Atangana and Baleanu [13] establish a new model, the so called generalized ABC differential operator involving the Mittag-Leffler (ML) and non-local type kernel, which has an anti-derivative fractional integral operator [14]. This new definition of derivative is shown to be more efficient for the SIR model with generalized incidence rate compared to the other existing fractional models [15] and an exothermic reactions model having a constant heat source in porous media with power problem [16]. More recently, Atangana and Atangana [17] solved the system of ABC fractional partial differential equations to see the side effects of facemasks, which is necessary to protect from COVID-19 and concluded that appropriate facemasks that could help avoid re-inhaling enough CO2 should be used every time one is in open air even when alone especially in a windy environment. Later, Khan et al. [18] explored the dynamics of COVID-19 with quarantine and isolations with real statistical cases reported in the mainland China, where the model is developed by using a fractal-fractional derivative in the Atangana–Baleanu sense. They indicate that their results are useful in the early eradication of the disease in the community. Similarly, Sane in [19] also consider the SIR epidemic model with Mittag-Leffler fractional derivative and delay, he uses a Lyapunov direct method and analyzes the global asymptotic stability of the disease-free equilibrium and the endemic equilibrium. Naik et al. [20] proposed an Atangana–Baleanu–Caputo operator for the transmission of COVID-19 epidemic, where they used real data from Pakistan. They show that fractional operators have many advantages over the existing non-integer-order types.

The studies from [13] to [20] provide us with enough motivation to study and analyze a new fractional SIR model involving the ABC operator with complete memory effects. To the best of our knowledge, this is the first time to study a non-local and non-singular derivative operator for the model of the SIR, which has nonlinear incidence and recovery rates, where we also examine the existence and uniqueness results with the Atangana–Baleanu derivative by using a fixed-point method.

The rest of this paper is organized as follows. Some preliminary results are given in Sect. 2. Section 3 is devoted to the new fractional SIR model and its main characteristics. We first derived the solution of the problem and then provided the existence and uniqueness result results in Sect. 4. In the last section, we demonstrate the numerical solutions for Eq. (10) by using the ABC-derivative in a graphical method. We also discussed the necessity of maintaining enough number of hospital beds for controlling of the infectious diseases.

Preliminary results

We note some background material for the Caputo and Atangana–Baleanu fractional derivatives and Laplace transform are presented [13, 21–23], which we use later.

Definition 2.1

Mittag-Leffler function is defined as

Definition 2.2

Let be a function, then the fractional Caputo derivative of order α is defined as [13] where and . Obviously, as , tends to the .

Definition 2.3

Let , , then the Atangana–Beleanu–Caputo (ABC) derivative is defined as [13] where the term is the normalization function with . The associated fractional integral is given by

If , the initial function is recovered in Eq. (3), we have the ordinary integral for .

Proposition 2.1

If , we have

Theorem 2.1

The Laplace transform of fractional Caputo derivative and Atangana–Beleanu–Caputo (ABC) derivatives is given by and

Property 2.1

The inverse Laplace transform of specific functions can be obtained:

Mathematical model and discussion

Let us consider the total population to consist of three sub-categories: susceptible , infected and recovered individuals , where as seen in Fig. 1. There are many factors involved in modeling the infectious diseases and these factors substantially affect the dynamical behavior of the models such as incident rate and recovery rate. In the classical SIR model, incident rate and linear recovery rate are used. However, as for COVID-19, many infectious diseases show multiple peaks or periodic oscillations during the outbreak or in progress. Therefore, several nonlinear incidence rates have been suggested since they can produce rich dynamics for the epidemic models [24, 25]. Liu et al. [26] proposed the following form of nonlinear saturated incidence rate to include the effect of behavioral changes: where . This model is due to Ruan et al. [27, 28]. More recently, a more general form was considered by Rao et al. [29].

Figure 1

Time development of susceptible (a) and infected (b) individuals for integer of order derivatives, and

On the other hand, in view of the diversity of public resources, it is not enough to study the effect of treatment capacity on a large scale at this time, if we want to further explore the influence of hospital bed number on disease transmission. Then we need to define the nonlinear recovery rate introduced by [30] where is the maximum per capita recovery rate, when health care resources suffice and there are few infectious individuals, is the minimum per capita recovery rate due to lack of clinical resources, b is the ratio between the number of hospital bed and the population, all the above parameters are non-negative. It has the following properties: (i)  and , i.e., is a decreasing function of I. (ii) and , i.e., is an increasing function of b.

Here, we assume individuals move from compartment S to I with transmission rates , and when an individual is infected, then the individual either recovers at a rate or dies at the rate of γI. Then this complex SIR model can be represented by following nonlinear system of equations: where refer to total number of people and n is the birth rate, μ is a positive constant representing the natural death of the population, γ is a positive constant representing the disease induced death rate, is the minimum per capita recovery rate due to the function of basic clinical resources and is the maximum per capita recovery rate due to the sufficient health care resource and few infectious individuals as well as the inherent property of a specific disease.

The model in (11) is not capable of representing the influence of memory effects with respect to time derivative. In order to account for the memory effect, in addition, we modify the fractional operator by an auxiliary parameter σ, having the dimension of second, to ensure that the right- and left-hand sides of the resultant equation possess the same dimension . Consequently, we can reformulate the system in (9) by using the ABC derivatives as where . In the first two equations in (12), does not appear, it is therefore assumed the first two equations are sufficient to handle the full system. We modify the fractional operator by an auxiliary parameter σ, having the dimension of day, to ensure that the right- and left-hand sides of the resultant equation possess the same dimension

Non-negative solution

In this section, we show that the feasibility region of the system (12) is given by the closed set which is positively invariant.

Lemma 3.1

The closed set defined in (13) is positively invariant with respect to the fractional model (12).

Proof

Adding all three equations in (12), we have

Applying the Laplace transform on both sides, we obtain

Using Theorem 2.1 and the properties of the Laplace transform, we got

Hence now, applying the inverse Laplace transform and using Property 2.1, we find where , since the Mittag-Leffler function has an asymptotic behavior, we obtain as . Hence, a solution of the fractional model (12) stays in Ω for every . Consequently, the closed set is a positively invariant set regarding the fractional model (12). □

Equilibrium points, basic reproduction number and local stability of equilibrium points

The equilibrium points of the system are the zeros of second side of Eq. (12), obviously, disease-free equilibrium point i.e., and or . Therefore, the model (12) has a threshold parameter , known as the basic reproduction number, this is defined as the number of secondary infections produced by a single infection in a completely susceptible population. It is not difficult to write the reproduction number for the disease-free equilibrium point of Eq. (12) as [5]

The existence of the endemic equilibrium point(s) can be determined by the relation where is the solution of the following equation:

Here, the coefficients are given by

The solution of quadratic equation in (21) is given by

Then we have the three cases of

If , so, and , if no endemic equilibrium, if , we have no endemic equilibrium for , one positive endemic equilibrium for and two endemic equilibrium points for .

If , then, if , then system unique endemic equilibrium point.

If , then, if , then system unique endemic equilibrium point.

Theorem 3.1

The disease-free equilibrium point is locally asymptotically stable and if then it is unstable if .

The jacobian matrix corresponding to system (12) is

Then the solution of characteristic equation associated to is given by

Hence, if , the disease-free equilibrium point is locally asymptotically stable, if , it is unstable. □

Theorem 3.2

A sufficient condition for the endemic equilibrium to be locally asymptotically stable is and , where

See the book of Chou and Friedman [31], we here note that

Therefore, local asymptotic stability of the endemic equilibrium point is a function of and the number of available hospital beds . □

Theorem 3.3

For the system (12), if , i.e. , the disease-free equilibrium is globally asymptotically stable.

Consider the Lyapunov function in , Calculating the time fractional derivative of the , we obtain

Utilizing the system (12), we have

The Lyapunov–Lasalle theorem implies that solutions in Ω approach the largest positively in-variant subset of the set , i.e., the plane . In this plane, as . Thus all solutions in the plane go to the disease-free equilibrium . Therefore is globally asymptotically stable. □

Existence and uniqueness of the solution

Applying the AB-fractional integral on both sides of (12), we obtain

By using Proposition 2.1, we have

Using Eq. (3) in the above, we find and where

The above system of Eqs. (32)–(33) is the solution of the systems of Eq. (12),

Let with subnorm be a Banach space of the real valued continuous functions. with the norm and .

Theorem 4.1

The kernels, and , satisfy the Lipschitz condition and contraction if the following inequality holds:

Proof

Consider , Let and be two functions, then we have where is bounded since is bounded, hence the Lipschitz condition is satisfied. If , then becomes a contraction mapping. Similarly, for , we find where

Since and are bounded, hence the Lipschitz condition is satisfied, if , becomes a contraction mapping. □

Theorem 4.2

Suppose that the following condition holds:

Then the fractional epidemic model (12) has a unique solution for .

Recursively, the expressions in (35) can be written as

As usual taking the difference of successive approximations, we obtain

It is not difficult to see that and , taking the sup norm on both sides, and using the Lipschitz condition on kernel functions and triangle inequality, we obtain

Starting from , we obtain by back substitution from (30) and (31), respectively, and

It is clear from our hypothesis that and as . We obtain a convergent series as follows: and

Since , are a Cauchy series in , and this series is uniformly convergent, it has a unique limit (Adams). So as , in Eq. (18), and have unique fixed points and therefore the system in (7) has a solution and this solution in unique. □

Numerical results and discussion

In this section, we present the result, for which we use COVID-19 data. The approach developed in [31] is used to approximate the ABC integral. We first use the fundamental theorem on the system (6) and obtain and

At time  , we obtain following iterative scheme:

Now, the function is approximated by Lagrange interpolation on as

Substituting (55) into (53)–(54) and by integration, we have

Simulation

Case I: The world

We need to use real data for COVID-19 to determine the parameters, the current birth rate for the world in 2020 is 18.077 births per 1000 people, and the death rate is 7.612 per 1000 people [16]. The world’s population on February was , so , according the WHO report [32] the COVID-19 mortality rate . The other parameters can be found by a curve fitting technique with real reported data for COVID-19. Figure 1(a) shows the development of susceptible individuals decrease with the time as expected for integer-order derivatives, and . Figure 1(b) shows the variations of infected individuals with the time, it increases with time and reaches the plateau region as seen in each country, where we used same parameters as Fig. 1(a). We compare the prediction fractional-order derivatives in Fig. 2. It is well known that an increase in the number of active cases and reaching the plateau region are different for each country and the number of suspected population decreases and reaches the steady case. This is also different for each country. On the other hand, from Fig. 2, it is clear that the active or suspected population not only is a function of time but also a function of the fractional order. This fact proves that the fractional-order derivative can be used the adjust the modeling to the different countries. The other point for fractional-order derivative can also be used as a free parameter to fix the modeling with the data. The effect of the hospital beds for controlling of the diseases are given in Fig. 3(a)–(b). Figure 3(a) shows the difference of the results between and . In Fig. 3(a) we show the infected individuals’ difference between and could reach 400,000 daily cases, if we assume that only 20% of infected individuals need medical treatment at hospital, then, if we multiply this one by COVID-19 mortality rate 0.034, we obtain 2720 people life can be saved if we provide sufficient hospital beds. Figure 3(b) shows the difference of the suspected individual results of and , this is expected because active infected cases for always are bigger than active infected cases for . We also note that our numerical results converge to the disease equilibrium point.

Figure 2

Comparison of the results of integer order derivatives (blue line, ), fractional order of derivatives ( (red line), (black line) and (green line) for and

Figure 3

The differences of active infected cases and suspected individuals for and

Case 2: Saudi Arabia

In the second case, we provide a numerical simulation using real data for the transmission model of COVID-19 in Saudi Arabia. According to the WHO report, the total population of Saudi Arabia so , we assumed COVID-19 mortality rate . On April 8, and . In Fig. 4(b), we explored the prediction of active cases for different values of fractional order, (blue line), 0.9 (red line), 0.7 (black line) and 0.5 (green line), when we compare modeling prediction with the real data (WHO report active cases in Saudi Arabia [32]), we see that model prediction overestimate sometimes lower estimate, but again it seems that produces the best result. At this point, we note that even our suggested modeling seems simple, but we believe that prediction of this model can be used as a basis for more complicated model predictions. Furthermore, we have one more adjustable parameters in fractional-order derivatives, this gives us enormous flexibility. Further, in Fig. 5, we concentrate on the effect of half saturation constant k over the suspected and infected individual for several values of fractional order of the derivations. We can see as before the effect of the fractional order provides us with enough flexibility and suspected individuals decrease with time and eventually reach the equilibrium point as we expected again the time to reach the equilibrium point to depend on the fractional order of the derivative.

Figure 4

Comparison of between the results of integer order derivatives (blue line, ), fractional order of derivatives ( (red line), (black line) and (green line) for , and

Figure 5

Comparison of between the results of integer order derivatives (blue line, ), fractional order of derivatives ( (red line), (black line) and (green line) for , and

Conclusion

In this study, consideration is given to a complex mathematical fractional SIR model for spreading of the COVID-19 virus, we first show the feasibility region where the solution exists. The equilibrium points (disease-free and endemic equilibrium) of the system are found. The local and global asymptotic stability of the system are discussed. Using the fixed-point theorem, existence and uniqueness of solution is proven. A numerical solution of the system is performed for several values of the parameters involved in the modeling. We showed the importance of maintaining a sufficient number of hospital beds. We also discussed and compared the prediction of our modeling with real word data. Finally, we concluded that using the fractional derivative provides us with one more parameter, which can be used to simulate the real word data.