Mathematical model of SIR epidemic system (COVID-19) with fractional derivative: stability and numerical analysis

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0}$\end{document}R0 is less than unity and unstable when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0} > 1$\end{document}R0>1. The analysis shows that the phenomenon of backward bifurcation occurs when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0}<1$\end{document}R0<1. 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.


Introduction
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.

The dimensional model
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

Positivity of solution
In this section, we prove that under nonnegative conditions, the model solutions are positive.
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.

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

Equilibria
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.

Theorem 2 System (3)-(5):
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.

Stability
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 .

Local stability of the disease-free equilibrium
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.

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

Stability of equilibria E 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.

Numerical simulations
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.

The model with fractional derivative
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.

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

Conclusion
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.
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.