key: cord-1051521-s4xly5fi
authors: Alqahtani, Rubayyi T.
title: Mathematical model of SIR epidemic system (COVID-19) with fractional derivative: stability and numerical analysis
date: 2021-01-04
journal: Adv Differ Equ
DOI: 10.1186/s13662-020-03192-w
sha: af583f6ac97b53314ad8aff8f85f939821891854
doc_id: 1051521
cord_uid: s4xly5fi

In this paper, we study and analyze the susceptible-infectious-removed (SIR) dynamics considering the effect of health system. We consider a general incidence rate function and the recovery rate as functions of the number of hospital beds. We prove the existence, uniqueness, and boundedness of the model. We investigate all possible steady-state solutions of the model and their stability. The analysis shows that the free steady state is locally stable when the basic reproduction number [Formula: see text] is less than unity and unstable when [Formula: see text] . The analysis shows that the phenomenon of backward bifurcation occurs when [Formula: see text] . Then we investigate the model using the concept of fractional differential operator. Finally, we perform numerical simulations to illustrate the theoretical analysis and study the effect of the parameters on the model for various fractional orders.

The spread of Covid-19 diseases is a very complex phenomenon carried out by many researchers. Many mathematical models were proposed including complex and simple mathematical models to understand the disease behavior. Faal et al. [1] proposed a model for the spread of the COVID-19 disease taking into account the superspreader, hospitalized, and fatality class. The authors analyzed the local stability of the steady-state solution and the model sensitivity. Mandal et al. [2] introduced a mathematical model taking into account a quarantine class and governmental intervention measures. In this study, the authors consider the basic reproduction number as an important parameter in analyzing the dynamics of the model. Recently, significant works were carried out to study the behavior of COVID-19 by means of mathematical models. Lin et al. [3] proposed SEIR models for the COVID-19 using data from China considering the impact of social isolation policies including governmental actions. The model successfully captures the course of the COVID-19 outbreak, whereas Wells et al. [4] and Gostic et al. [5] consider the impact of travel restrictions and border control on the global spread of the COVID-19.

The SIR model is commonly used for disease modeling, in particular, for the COVID-19 analysis [6] [7] [8] . The dynamic behavior of SIR model, including the stability, bifurcation, and chaos, has been studied over many decades [9] [10] [11] [12] . In most studies the authors assume that the recovery rate is a constant. However, in reality the recovery rate depends on time of recovering process such as the health system, including the number of hospital beds and medicines.

In recent years, many researchers have studied the systems of differential equations with fractional operators [13] [14] [15] . The epidemic models involving a fractional operator were also investigated by many authors because they deeply show biological and physical perspectives of the diseases [16, 17] .

Rao et al. [18] studied an SIRS epidemic model assuming different death rates for each subclass, and the fraction of newborn children is represented by the parameter p. In this paper, we propose and analyze the extended SIRS epidemic model presented in [18] with the concept of fractional differential operator. In fact, we propose and study a model including three nonlinear differential equations with general incidence rate function and nonlinear recovery rate depending on the health system. The main focus of this study is analyzing the basic properties of model and demonstrating the stability properties of the model.

The rest of the paper is arranged as follows. We propose a dynamical model in Sect. 2. Then we formulate and establish the existence, uniqueness, positivity, and boundedness of solutions in Sect. 3. The steady-state solutions of the model and their stability are studied in Sects. 4 and 5, whereas numerical simulations of the steady-state solution brunches has is presented in Sect. 6. Section (7) contains a detailed dynamic behavior of the model with fractional derivative. We finish this study with conclusion in Sect. 8.

In this section, we extend the model suggested in [18] to include a nonlinear incidence rate and recovery rate. The recovery rate is a function of both the hospital bed-population ratio b 1 > 0 and the infected I. Thus the recovery rate α is given by [19] 

where the parameter α 1 and α 0 are the maximum and minimum per capita recovery rates, respectively. The nonlinear incidence rate is generalized by the function

Thus the system of differential equations is given by

where the total population is split into three parts: S(t) is the susceptible population, I(t) is the infected population, and R(t) is the recovered population, so that N = S + I + R. The details and interpretation of the model can be found in [18] . We assume that all parameters are positive.

3 Basic properties of model

In this section, we prove that under nonnegative conditions, the model solutions are positive.

Proof Let x(t) = (S(t), I(t), R(t)) be the solution of system under initial conditions x 0 = (S(0), I(0), R(0)) = (S 0 , I 0 , R 0 ) ≥ 0. By the continuity of solution, for all of S(t), I(t), R(t) that have positive initial values at t = 0, we have the existence of an interval (0, t 0 ) such that S(t), I(t), R(t) ≥ 0 for 0 < t < t 0 . We will prove that t 0 = ∞.

If S(t 1 ) = 0 for t 1 ≥ 0 and other solutions stay positive at t = t 1 , then

This ensures that at any time the solution reaches the axis, its derivative increases, and the function S(t) does not cross to negative part. We can show by similar analysis that

dR dt (t = t 1 ) = pb + αI ≥ 0. 

is positively invariant and attracts all solutions in R 3 + .

Proof Let W (t) = S(t) + I(t) + R(t). Then from the system (3)-(5) we have

This implies that

Solving (10), we obtain

where W (0) is the initial condition. Thus 0 < W (t) < b μ as t reaches infinity, and hence is a positively invariant and attractive set.

We use the next-generation matrix method [24] to calculate the reproduction number R 0 of model (3)- (5):

In this section, we consider the number of equilibrium solutions of model (3)-(5). It is clear that the model has a disease-free equilibrium given by

The non-free steady state of model (3)-(5) can be obtained by setting the right sides to zero. From equations (3)-(5) we have

Substituting equations ( (14) and (15)) into equation (3), we obtain

where c 0 , c 1 , c 2 , and c 3 are defined by Table 1 Number of possible positive real roots of equation (16) . c 4 = basic reproduction number R 0 , c 5 = sign change number, c 6 = possible number of positive real roots

If R 0 = 1, then c 0 = 0, so equation (16) reduces to the equation

where I = 0 is the disease-free equilibrium. By equation (16) Table 1 .

1. has a one equilibrium if the basic reproduction number is greater than 1 and Cases 1, 5, and 7 are satisfied; 2. can have more than one equilibrium if the basic reproduction number is greater than 1 and Case 3 is satisfied; 3. can have two or more equilibria if the basic reproduction number is less than 1 and Cases 2, 4, and 6 are satisfied.

The existence of multiple steady state suggests the possibility of backward bifurcation where the phenomenon of three branches of steady-state equilibrium occurs at the same point.

In this section, we focus on analysis of the stability of the equilibrium of equations (3)-(5). We study the stabilities of two types of the disease equilibrium, that is, E 0 and E 1 .

In this section, we study the stability of the free equilibrium E 0 . The Jacobian matrix of system (3)-(5) at E 0 is

where

The eigenvalues of matrix (19) are given by

A simple calculation shows that J 22 = R 0 -1. So, we have the following result.

The free steady-state solution E 0 is locally asymptotically stable if R 0 < 1 and is unstable if R 0 > 1.

In this section, we show that the nonfree steady-state solution E 1 of system (3)- (5) is stable under specific condition. The Jacobian of the system can be written as

where J 11 = β 1 I(Ia 3 + a 1 ) (Ia 3 + a 2 S + a 1 ) 2 + μ 1 ,

J 12 = β 1 S(a 2 S + a 1 ) (Ia 3 + a 2 S + a 1 ) 2 ,

From equation (4) we get the following relations:

By simple analysis we get that the characteristics equation of J(E 1 ) is

where B 1 = J 22 + J 11 + μ 1 + μ 3 + γ 1 ,

We further use the Rough-Hurtwiz criterion to show the stability of the steady state E 1 .

We have

By the Routh-Hurwitz theorem E 1 is locally asymptotically stable when B 1 > 0, B 3 > 0, and B 1 B 2 -B 3 > 0. Theses conditions are satisfied when the following condition holds:

Thus we have following results. Hence

Now let β 1 = β * 1 be the bifurcation parameter. When R 0 = 1, we have the following relation:

and the model equation has one zero eigenvalue, and the other eigenvalues are negative. The behavior of the system near β 1 = β * 1 can be studied by applied the center manifold theory. The Jacobian matrix at free steady state E 0 is

The right eigenvectors can be obtained as W = (w 1 , 

The coefficient a is obtained as

The bifurcation parameter b at E 0 is given by

and can be obtained as

Clearly, b is always positive. According to [25, Theorem 4.1] , the backward bifurcation phenomenon exists when the coefficient a is positive. Thus the condition for backward bifurcation is given by

The existence of the backward bifurcation at R 0 = 1 requires condition (44) to be satisfied. When the number of hospital beds b 1 is below the critical point b 1,cr , the number of hospital beds open to the public is below demand, and as a result, some patients fail to access to healthcare. In this situation, there remains a high infection leading to a backward bifurcation.

In this section, we carry out some numerical calculations to support our theoretical results. The values of parameters used for numerical simulations are indicated in Table 2 . We study the branch of steady state with respect to the model parameters. Figure 1 shows the curves of the infected population I for different values of b 1 , donated by the number of hospital beds and a specific value of general incidence rate (a 1 = a 2 = a 3 = 1). It shows that there is a forward bifurcation at R 0 = 1.

If we decrease the value of b 1 from 2 to 1.6, then the backward bifurcation does not occur. These values are higher than the critical value of b 1,cr = 1.64. If we decrease the value of b 1 to 0.1, less than the critical value b 1,cr = 1.64, then we can observe from Fig. 1(a) that the backward bifurcation occurs. Note that in Fig. 1 (a) the above line of the curve is a stable state and the below line of the curve is an unstable state. This result indicates that Assumed. a 3 1 Assumed. Table 2 in managing an infectious disease the number of hospital beds plays a significant role. Figure 2 shows the effect of the value of b 1 on the curve when the backward bifurcation occurs. We observe that as the value of b 1 decreases, the area of the curve increases. Figure 2 shows the infected population size I as a function of reproduction number R 0 when the parameter b 1 is varied for the case R 0 < 1. It illustrates that as the value of b 1 increases, the infected population size I decreases. It also shows the existence of a backward bifurcation, and the area of backward bifurcation curve decreases as the value b 1 increases.

We consider the model with the Caputo-Fabrizio fractional derivatives

Here we have 0 < α 3 < 1 and

We present the existence of positive solution of the system,

Then

We can similarly show that

Thus for all t ∈ [0, t], we have that S(t), I(t), and R(t) are positive.

Here we present the condition under which the system of equations has a unique solution.

To achieve this, we have

We will show that, for all i = 1, 2, 3,

:

where Table 2 f 2 (S j-1 , I j-1 , R j-1 , τ j-1 ) (n + 1j) α 3 +1 -(nj) α 3 (nj + 1 + α 3 )), R(t n+1 ) = R(0) + ( t) α 3 (α 3 + 2) n j=0 f 3 (S j , I j , R j , τ j ) (n + 1j) α 3 (nj + 2 + α 3 )

-(nj) α 3 (nj + 1 + α 3 )) f 3 (S j-1 , I j-1 , R j-1 , τ j-1 ) (n + 1j) α 3 +1 -(nj) α 3 (nj + 1 + α 3 )). 

In this paper, we considered the SIR model with general incidence rate function and nonlinear recovery rate to model the spread of disease. The nonlinear recovery rate depends on the influence of health system.

We proved the existence, uniqueness, and boundedness of the model solution. We studied all possible steady-state solutions of the model and details of stability and also derived the reproductive number. The analysis shows that the free steady state is locally stable when the reproductive number is less than unity and unstable otherwise. The model shows the phenomenon of backward bifurcation when R 0 < 0 and the parameter b 1 is less than the critical value given by b 1 < b 1,cr = μ 1 S 0 [a 2 (μ 3 + γ 1 )(α 1α 0 )S 0 + a 1 (μ 3 + γ 1 )(α 1α 0 )] [α 1 + μ 2 ][S 0 a 3 μ 1 (μ 3 + γ 1 ) + a 1 (α 1 μ 3 + μ 2 [γ 1 + μ 3 ])] .

When the parameter b 1 is sufficiently greater that the critical value b 1,cr , the disease infection decreases because the number of hospital beds increases. Therefore, to treat the disease in a community, the hospital resources must be improved. Finally, we applied the theory of fractional derivatives to the model for different values of fractional orders. We used the numerical technique of Atangan and Toufiq, which is very accurate for solving fractional differential equations.

Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan

A model based study on the dynamics of COVID-19: prediction and control

A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

Impact of international travel and border control measures on the global spread of the novel 2019 coronavirus outbreak

Estimated effectiveness of symptom and risk screening to prevent the spread of COVID-19

A time-dependent SIR model for Covid-19 with undetectable infected persons

Predicting the spread of COVID-19 using SIR model augmented to incorporate quarantine and testing

Simulating the progression of the COVID-19 disease in Cameroon using SIR models

Hopf bifurcation in two SIRS density dependent epidemic models

Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate

Stability and Hopf bifurcation for a delayed SIR epidemic model with logistic growth

Stability of Hopf bifurcation of a delayed SIRS epidemic model with stage structure

Numerical analysis for the fractional diffusion and fractional Buckmaster's equation by two step Adam-Bashforth method

A new numerical approximation of the fractal ordinary differential equation

Analysis of series RL and RC circuits with time-invariant source using truncated M, atangana beta and conformable derivatives

Mathematical model of Ebola and Covid 19 with fractional differential operators: non-Markovian process and class for virus pathogen in the environment

Facemasks simple but powerful weapons to protect against COVID-19 spread: can they have sides effects

Complicated endemics of an SIRS model with a generalized incidence under preventive vaccination and treatment controls

Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds

A stochastic SIRS epidemic model with infectious force under intervention strategies

Complete global stability for an SIRS epidemic model with generalized non-linear incidence and vaccination

Global analysis of an epidemic model with nonmonotone incidence rate

Modeling the spread and control of Dengue with limited public health resources

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Dynamical models of tuberculosis and their applications

New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models

The author would like to thank the anonymous referees for their valuable suggestions, which have greatly helped in improving the presentation of this paper.

≤ 3(pb) 2 + 3(μ 3 + γ ) 2 |R| 2 + 3|α| 2 |I| 2 (66) ≤ 3(pb) 2 + 3(μ 3 + γ ) 2 |R| 2 + 3 sup |α| 2 |I| 2 (67) ≤ 3(pb) 2 + 3(μ 3 + γ ) 2 |R| 2 + M 1 (68)Therefore, under the conditionthe system admits a unique solution.

In this section, we present the numerical solution of the equations. We use the numerical scheme of Atangan and Toufiq [26] . To use their scheme, we haveThe next step is converting the above toFollowing their scheme, we havef 1 (S j-1 , I j-1 , R j-1 , τ j-1 ) (n + 1j) α 3 +1 -(nj) α 3 (nj + 1 + α 3 )),n j=0 f 2 (S j , I j , R j , τ j ) (n + 1j) α 3 (nj + 2 + α 3 )

No funding available.

Not applicable.

The author declares that they have no competing interests.

The author worked in the derivation of the mathematical results and read and approved the final manuscript.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 18 November 2020 Accepted: 17 December 2020