# Stability Analysis of COVID-19 Epidemic Model of TypeSEIQHR with Fractional Order

Authors: Sid Ahmed Ould Beinane; Mohamed Reda Lemnaouar; Rabie Zine; Younes Louartassi
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this article, we consider a fractionalSEIR model, denoted by the SEIQHR model, which aims to predict the outbreak of infectious diseases in general. In particular, we study the spread of COVID-19. The fractional order offers a flexible, appropriate, and reliable framework for pandemic growth characterization. Firstly, we analyze some elementary results of the model (boundedness and uniqueness of solutions). In addition, we establish certain conditions to ensure the local stability of the disease-free and endemic equilibrium points. Based on analytical and numerical results, we conclude that coronavirus infection (COVID-19) remains endemic, which requires long-term prevention and intervention strategies.


---

## Body

## 1. Introduction

COVID-19 represents the disease caused by a virus in theCoronaviridae family, and SARS-CoV-2 appeared at the end of 2019 in Wuhan, China [1]. It spread rapidly around the world, causing a worldwide epidemic [2, 3].COVID-19 is a respiratory illness that can be fatal for old patients or other chronic diseases. It is transmitted through close contact with infected people. The disease could also be spread by asymptomatic patients [4], but scientific data are lacking to attest with certainty. In addition, the international health organization has imposed security measures to control the spread including isolation, quarantine, increased home confinement, promotion of wearing face masks, travel restrictions, the closure of public space, and the cancellation of events. The number of confirmed cases increased rapidly to reach more than 517 million cases, and approximately 6.25 million deaths were recorded worldwide as of May 2022.There are severalSEIR (susceptible-exposed-infectious-recovered) models that study infectious diseases in general [5–7]; since the appearance of the virus, the SEIR models proposed aim to study the behavior of this epidemic [8].In the last decade, fractional models have been used to model diseases in general. For example, Singh [9] proposed a new fractional model of blood alcohol model using the Hilfer fractional operator. Kumar et al. [10] proposed a model to study the transmission dynamics of dengue, and they considered the generalized Caputo fractional operator. Since the apparition of COVID-19, researchers have used fractional derivatives used in the modeling of this infectious disease. As an example, the work in [11–15] has developed the classic SEIR model known by fractional models (SEIQHR and SEIQRd) for epidemic analysis of COVID-19 worldwide. This development is founded on the establishment of new quarantine conditions and the hospitalization of confirmed cases, which are considered as epidemic parameters for COVID-19. Recently, many types of fractional derivatives have been employed to model the propagation of COVID-19. For instance, Danane et al. proposed a fractional-order model of the disease (COVID-19) with government action and individual response by the Caputo operator. Bonyah et al. used the Atangana–Baleanu operator for investigating the fractional optimal control dynamics of a coronavirus model. Yadav et al. studied the dynamics of the fractional-order COVID-19 model with memory effect using the Liouville–Caputo operator, so they employed the Adams–Bashforth–Moulton approach to find an approximate solution. Consequently, this development shows us more precisely the behavior of this epidemic.Recently, fractional calculus has shown wide applicability in many fields, which can be obtained by extracting a dynamic behavior of biological systems shown by a mathematical formulation of integer derivatives. Fractional models have been proposed to study several phenomena involving the memory effect, including epidemic behavior [16], and they offer more flexibility than classical integer-order models to fit the data accurately [17]. Several researchers presented the differential fractional theory [18–21] for this reason, and many papers studied fractional biological models [22, 23]. Recently, new types of fractional derivatives have been developed. For instance, Khalil et al. [24] proposed a new fractional derivative and its properties. Subsequently, Abdeljawad [25] developed the properties of the conformable fractional derivative. After that, Benmakhlouf et al. [26] studied the finite time stability (FTS) and finite time boundedness (FTB) of the conformable fractional derivative. Currently, Hattaf [27, 28] introduced a new definition of the fractional derivative with a non-singular kernel in the sense of Caputo to generalize the various types by making these properties, and he proposed an approach for studying the stability of the latter. After that, Hattaf et al. [29] proposed a new numerical method for approximating the generalized Hattaf fractional derivative involving a non-singular kernel based on Lagrange polynomial interpolation.Inspired by the aforementioned work and previous literature, we provide the compartmental model of COVID-19 with a standard incidence rate explored in [30] considering fractional Caputo derivatives for a better insight into the disease. The aim of the epidemic model (4) is to describe the dynamic behavior of the disease and predict the tendency for the disease to spread. Therefore, we have proposed a SEIQHR model. First, we analyze the qualitative properties of this model, including the existence and uniqueness of the disease-free and endemic equilibrium points. Second, we establish the conditions to ensure local asymptotic stability of disease-free and endemic equilibrium points. Third, we study the contribution of parameters to reproductive numbers. At the end, a numerical simulation is proposed to show the behavior of the different compartments of our model (4), and we will give some clarification into the interpretation and role of fractional derivatives.


## 2. Preliminary Results

We begin by giving some definitions of fractional calculus.Definition 1.
(see [31]). The Caputo fractional derivate of function f order is defined by(1)dβft=d−n−βdndtnft=1Γn−β∫0tfnxt−xn−β−1dx,wheren−1≤β<n∈ℕ∗.Definition 2.
(see [32]). Let the fractional system(2)dβXt=ft,X,t>0,Xt0=X0>0,whereβ∈0,1,t0>0, andf:t0,+∞×Ω⟶ℝn,Ω⊂ℝn. Ifft,Xfulfill the local Lipschitz condition with respect toX, there exists a unique solution of the above system.Lemma 1 (see [33]).
Let the fractional-order system(3)dβxt=fx,t>0,xt0=x0>0,with0≤β<1,x∈ℝ. The equilibrium points of the above system are calculated by solving the following equation: fx=0.x∗are locally asymptotically stable if all eigenvaluesλof the Jacobian matrixJ=∂f/∂xatx∗evaluated of the equilibrium points satisfy Matignon’s condition (see [34]) argλ>βπ/2.


## 3. Description of Model

In this section, we present a mathematical formulation to model the behavior of compartments; first we will split our population into six categoriesSt Susceptible, Et Exposed, It Infected, Qt Quarantined, Ht Hospitalized and RtRecovered. The description of the parameters is indicated in Table 1, so that the description of the interactions between the compartments is illustrated in Figure 1. The incidence rate plays an important role in describing the evolution of an infectious disease. Based on the spread of different diseases, there are many forms of incidence rates [30, 35–37]. In our model, we assume the bilinear incidence rates [30] (β1SI and β2SE). Thus, our proposed model (4) is in the following form:(4)dβSt=Λ−β1SI−β2SE−dS−asS,dβEt=β1SI+β2SE−αE−dE,dβIt=αE−δI−ηI−μI−dI,dβQt=δI−kQ−εQ−dQ,dβHt=εQ+ηI−rH−dH,dβRt=rH+μI−dR−arR,with initial conditions(5)S0≥0,E0≥0,I0≥0,Q0≥0,H0≥0andR0≥0,where dβ is in the sense of Caputo fractional derivative and 0<β≤1.Table 1
Description of variables and parameters.
VARS and PRMExplanationSSusceptible individualsEExposed individualsIInfected individualsQQuarantined individualsHHospitalized individualsRRecovered individualsΛInflow number of susceptible individualsβ1Infection rates of the infected individualsβ2Infection rates of the exposed individualsαIncubation rateδRate at which symptomatic infections are diagnosed and quarantineddNatural rate mortalityμRate at which symptomatic infections are diagnosed and recoveredηRate at which symptomatic infections are diagnosed and hospitalizedεTransition rate of quarantined individuals to the hospitalized infected classrTransition rate of hospitalized individuals to the recovered classkThe death rate caused by the diseaseFigure 1
Schematic diagram of (4).Our model (4) is based on the following hypothesis. H1: the number of susceptible individuals added is constant each day. H2: the susceptible individuals S move into the exposed E class before the infected class. H3: when an individual is infected I, he is either quarantined Q, hospitalized H (severe case), or recovered. H4: a confined individual can be hospitalized if his situation is critical.The detailed description of (4) is presented in the following schema.


## 4. Boundedness and Uniqueness of Solutions

This section is dedicated to demonstrate the boundedness of the solutions of (4).Lemma 2.
The setΩ=S,E,I,Q,H,R∈ℝ+6:Nt≤N0+Λ/dis a region of attraction for all solutions initiating in the interior of the positive octant, whereN=S+E+I+Q+H+R.Proof 1.
We poseNt=St+Et+It+Qt+Ht+Rt; then,(6)dβN+dN=Λ−kQ≤Λ.
Using the theory of fractional inequality (see [31]), we obtain(7)Nt≤N0Eβ−dtβ+Λd1−Eβ−dtβ,where Eβz=∑k=0∞zk/Γβk+1 is Mittag-Leffler function [31], Γz=∫0∞xz−1e−xdx is Euler’s gamma function, and 0<Eβ−dtβ≤1 if t⟶∞, and we have 0<Nt≤N0+Λ/d, proving this lemma.
Otherwise, model (4) is equivalent to the model:(8)dβX=FX,where(9)X=SEIQHR,FX=Λ−β1SI−β2SE−dSβ1SI+β2SE−αE−dEαE−δI−ηI−μI−dIδI−kQ−εQ−dQεQ+ηI−rH−dHrH+μI−dR=F1XF2XF3XF4XF5XF6X.Theorem 1.
Assume that the initial conditionsX0≥0. Then, there exists a unique solution of system (4) defined on0,+∞.
The sufficient condition for the existence and uniqueness of the solution of system (4) in the region Ω×t0,Twith initial conditions X0 and L=maxd+2β1+β2M,d+2α+β2M,d+2η+μ+δ+β1M,d+k+2ε,d+2r,d is ‖F(X)−F(X’)‖‖1≤L‖X−X’‖1..Proof 2.
To prove the global existence and uniqueness of (4), consider the region Ω×t0,T, where Ω=S,E,I,Q,H,R∈ℝ+6:maxS,E,I,Q,H,R≤M,M>0.
For anyX=S,E,I,Q,H,RT, X′=S′,E′,I′,Q′,H′,R′T∈Ω,(10)FX−FX′1=∑i=16FiX−FiX′=β1SI−S′I′+β2SE−S′E′+dS−S′+β1SI−S′I′+β2SE−S′E′−α+dE−E′+αE−E′−δ+η+μ+dI−I′+δI−I′−k+d+εQ−Q′+εQ−Q′+ηI−I′−r+dH−H′+rH−H′+μI−I′−dR−R′≤d+2β1+β2MS−S′+d+2α+β2ME−E′+d+2η+μ+δ+β1MI−I′+d+k+2εQ−Q′+d+2rH−H′+dR−R′≤LX−X′1,where(11)L=maxd+2β1+β2M,d+2α+β2M,d+2η+μ+δ+β1M,d+k+2ε,d+2r,d.
Thus,F satisfies Lipschitz’s condition (see [32] and Definition 2) with respect to X.


## 5. Equilibria and Local Stability

In the first part of this section, we debate the existence of equilibria. It is evident that (4) has an infection-free equilibrium P0Λ/d,0,0,0,0,0. Also, other endemic equilibrium point P∗S∗,E∗,I∗,Q∗,H∗,R∗ is defined after in (13). The basic reproduction number of (4) is(12)R0=Λβ1α+β2δ+μ+d+ηd+asδ+μ+d+ηα+d.Theorem 2.
(1)
IfR0≤1, then system (4) has one infection-free equilibriumP0Λ/d,0,0,0,0,0.(2)
IfR0>1, then (4) has endemic equilibrium pointP∗S∗,E∗,I∗,Q∗,H∗,R∗, whereS∗,E∗,I∗,Q∗,H∗,R∗are defined after in (9).Proof 3.
To get the endemic equilibrium point of the (4), in the interior of the equilibrium P∗S∗,E∗,I∗,Q∗,H∗,R∗,Proof.
i.e.,FiS,E,I,Q,H,R=0 for i=1,…,6, we get(13)S∗=Λδ+μ+d+ηβ1αE∗+β2E∗+dδ+μ+d+η,I∗=αE∗δ+μ+d+η,Q∗=δαE∗k+d+εδ+μ+d+η,H∗=ηε+k+d+εδαE∗k+d+εr+dδ+μ+d+η,R∗=rηε+k+d+εδ+μr+dk+d+εαE∗k+d+εr+dδ+μ+d+ηd,E∗=δ+μ+d+ηdβ1α+β2δ+μ+d+ηR0−1.
Then, ifR0>1, then (4) has endemic equilibrium point P∗.
Now, we analyzed the local asymptotic stability of disease-free equilibrium pointP0 and endemic equilibrium point P∗ for (4).Theorem 3.
The disease-free equilibrium pointP0is locally asymptotically stable, ifδ+μ+η+2d+α>β2Λ/dandR0<1.Proof 5.
To prove the local stability of equilibria, the eigenvalues of the Jacobian matrix of (4) are evaluated by(14)J0=J11J12J21J22,where(15)J11=−d00−β2Λdβ2Λd−α+dα−β1Λdβ1Λd−δ+μ+d+η,J12=000000δημ,J21=000000000,J22=−d+ε+kε00−r+dr00−d.
The characteristic equation ofP0Λ/d,0,0,0,0,0 is(16)λ+d2λ+r+dλ+k+d+ελ2+sλ+p=0,where s=δ+μ+η+2d+α−β2Λ/d,p=δ+μ+η+dd+α−β2Λ/d−αβ1Λ/d=d+αδ+μ+η+d1−R0.
The discriminant of equationλ2+sλ+p=0 is δ+μ+η−α+β2Λ/d2+4Λαβ1/d>0. Then, the eigenvalues of matrix (10) to the equilibrium point P0 are reel roots, so λ1=−k+d+ε<0, λ2=−r+d<0, λ3=λ4=−d<0, λ5+λ6=−s, and λ5λ6=p.
Ifs>0 and R0<1, then λ5+λ6<0 and λ5λ6>0. So, λ5<0 and λ6<0. The proof is completed.Theorem 4.
The endemic equilibrium pointP∗is locally asymptotically stable ifR0>1and condition (12) are realized.Proof 6.
In the same way as the previous proof, let(17)J∗=J11∗J12∗J21∗J22∗,where(18)J11∗=−β1I∗+β2E∗+dβ1I∗+β2E∗0−β2S∗β2S∗−α+dα−β1S∗β1S∗−δ+μ+d+η,J12∗=000000δημ,J21∗=000000000,J22∗=−d+ε+kε00−r+dr00−d.
From the Jacobian matrixJ∗, the characteristic equation at P∗ is(19)λ+dλ+r+dλ+k+d+ελ3+a2λ2+a1λ+a0=0,where(20)a0=δ+μ+η+dβ1I∗+β2E∗+dα+d−dβ2S∗−αdβ1S∗,a1=δ+μ+η+dα+2d+β1I∗+β2E∗−β2S∗+β1I∗+β2E∗+dα+d−β1α+dβ2S∗,a2=β1I∗+β2E∗+3d+α−β2S∗+δ+μ+η.
Using (13), we get(21)a0=dS∗αβ1+β2η+μ+d+δR0−1>0,becauseR0>1,a1=δ+μ+η+α+2ddR0−β2ΛR0,a2=δ+μ+η+α+2d+dR0−β2ΛdR0.
Firstly,λ1=−d<0, λ2=−r+d<0, λ3=−k+d+ε<0. Then, argλ1,2,3=π>βπ/2.
On the other hand, the discriminant of algebraic equationλ3+a2λ2+a1λ+a0=0 is this form (see [38]):(22)dP=18a1a2a0+a2a12−4a13−4a23a0−27a02.
Then, the equilibrium pointP∗ is locally asymptotically stable in one of the following cases:(1)
IfdP>0,P∗is asymptotically stable ifa1,a2>0anda2a1−a0>0for allβ∈0,1.(2)
IfdP<0anda1,a2>0,P∗is asymptotically stable ifβ<2/3anda2a1−a0>0.(3)
IfdP<0,a1,a2>0anda2a1=a0,P∗is asymptotically stable for allβ∈0,1.
Then,P∗ is locally asymptotically stable if R0>1 and (12) are realized.


## 6. Sensitivity Analysis

Sensitivity analysis shows us the impact of each parameter on the transmission of the disease. It is used to understand which parameters have a high impact on theR0 threshold. More specifically, sensitivity indices allow us to measure the relative change in a variable when a parameter changes. If this variable is differentiable with respect to the parameter, the sensitivity index is defined as follows [39]:(23)SaR0=∂R0∂aaR0,where a represent the contribution to the basic reproductive number R0.ForΛ, SΛR0=∂R0/∂ΛΛ/R0=1.Forβ1, Sβ1R0=∂R0/∂β1β1/R0=αβ1/αβ1+β2δ+μ+η+d.Forβ2, Sβ2R0=∂R0/∂β2β2/R0=β2δ+μ+η+d/αβ1+β2δ+μ+η+d.Forμ, SμR0=∂R0/∂μμ/R0=−β1αμ/αβ1+β2δ+μ+η+dδ+μ+η+d.Forη, SηR0=∂R0/∂ηη/R0=−β1αη/αβ1+β2δ+μ+η+dδ+μ+η+d.Forδ, SδR0=∂R0/∂δδ/R0=−β1αδ/αβ1+β2δ+μ+η+dδ+μ+η+d.Forα, SαR0=∂R0/∂α⋅α/R0=αβ1d−β2δ+μ+η+d/α+dαβ1+β2δ+μ+η+d.Ford, SdR0=∂R0/∂d⋅d/R0=−β1d+αd+η+μ+δ+dd+α+dd+η+μ+δ/d+αd+η+δ+μβ1α+β1μ+δ+η+d+β22d+αd+μ+η+δ/d+αβ1α+β1μ+δ+η+d. Using the values of parameters (14), we get the following table.The sensitivity index can depend on the system parameters, but it can also be constant. For example,SΛR0 means that an increase (decrease) in Λ by a given percentage will result in an increase (decrease) of R0 by the same percentage. Concretely, an increase of the values Λ, β1, β2 will increase the basic reproduction number by 100%, 15.23%, and 84.77%, respectively, and an increase of the values μ, η, δ, α, d will decrease R0 by 6.3%, 2.53%, 6.4%, 22.52%, and 94.01%, respectively (see Table 2).Table 2
Sensitivity index of the basic reproduction number.
ParametersSensitivity indexΛ+1β1+0.1523β2+0.8477μ−0.063η−0.0253δ−0.064α−0.2252d−0.9401


## 7. Simulations

In this section, we include a numerical analysis of model (4). The parameters are estimated. We will study the impact of certain parameters (α, μ, η and ε) on the solutions and the impact β of the fractional derivation order.For this reason, we take some hypothetical data in order to illustrate the results that we have already established in the previous sections.(24)Λ=105,β1=3.8×10−6,β2=7×10−6,α=0.2657,δ=0.3352,μ=0.33029,η=0.13266,ε=0.1259,k=1.1975×10−5,r=0.0149,d=3.051×10−5,S0,E0,I0,Q0,H0,R0=11×104,0,102,102,0,0.To resolve our fractional system (4), we have used the function “fde12″ in MATLAB. It is an implementation of Adams–Bashforth–Moulton prediction correction described in [40]. The convergence and accuracy of this method are studied in [41]. The implementation of several iterations of the corrector has been proposed in [42]. In this implementation, the discrete convolution is evaluated using the FFT algorithm described in [43], which makes the computational cost proportional to N∗logN2 instead of N2 in the classical implementation; N is the number of points in time for evaluating the solution.Figure2 shows the impact of order fractional derivative β as a function of time t. We remark in Figure 2 that the increase of β leads to the growth of the speed of convergence for E and I. Also, we observe that the increase of β leads to a decrease in the speed of convergence for Q and H.Figure 2
Impact ofβ according to time of endemic equilibrium P∗.Figure3 illustrates the impact of incubation rate α according to time t. We observe that the number of susceptible, exposed, and infected individuals will initially increase, but after a period of time it will decrease. We notice that the number of susceptible cases S and recovered cases R increases when α increases. Contrarily, the cases infected I, quarantined Q, hospitalized H, and exposed E are decreased. Indeed, by increasing the rate of incubation, the number of exposed cases decreases, which reduces the number of susceptible cases who will pass into the category exposed and then will be infected. This leads to a decrease in the number of infected and hospitalized cases. As a consequence, the number of recovered cases increases.Figure 3
Impact ofα according to time of endemic equilibrium P∗.In addition, Figure4 shows that when the value of μ increases, the number of suspected S and recovered R cases increases. Contrarily, the infected cases I, quarantined cases Q, and hospitalized cases H are decreased. This explains why if the rate μ increases, many of the infected cases will recover, causing the number of infected individuals to decrease and the recovered cases to increase. Therefore, the number of hospitalized cases will decrease.Figure 4
Impact of μ according to time of endemic equilibriumP∗.We can observe in Figure5 that when the value of η increases, the number of quarantined Q and hospitalized H cases increases. On the other hand, the infected cases I and recovered cases R decrease. This implies that if the percentage of infected individuals becomes severe, they must be hospitalized. As the number of hospitalized cases grows, the number of infected cases decreases. As a result, the number of recovered cases is reduced. Moreover, in Figure 6, when the value of ε increases, the number of quarantined cases Q increases. On the other hand, the hospitalized cases H and recovered cases R decrease.Figure 5
Impact of Η according to time of endemic equilibriumP∗.Figure 6
Impact of ε according to time of endemic equilibriumP∗.Figure7 illustrates the impact of δ according to time t. We notice that the number of susceptible cases S, quarantined cases Q, and hospitalized cases H increases when δ increases. Contrarily, the cases infected I and recovered R are decreased, but it has little impact on the cases exposed E. This is explained by the fact that a significant number of infected (noncritical cases) must be quarantined. Among these confined cases, some individuals may be hospitalized at a constant ϵ rate. Consequently, the number of hospitalized cases is increasing.Figure 7
Impact of δ according to time of endemic equilibriumP∗.


## 8. Conclusion

In this paper, we have modeled the COVID-19 epidemic according to susceptible, exposed, infected, quarantined, hospitalized, and recovered compartments. An operator Caputo of fractional derivatives was used to evaluate the memory effect on the epidemic behavior. We developed a mathematical approach to prove the boundedness, uniqueness, and the existence of solutions and to demonstrate a local stability of equilibrium points.Sensibility analysis revealed that the epidemic evolution is affected by the different model parameters. Furthermore, a numerical simulation illustrates the effect of memory in a fractional derivative, and the increase of the fractional derivation order speeds up the decrease ofE, I and increase of Q, H, and R. Consequently, the calibration of this parameter provides a correct adjustment of the real data. On the other hand, they show that the susceptible cases increase with an increase in the incubation rate α, and when the rates μ and η increase, the number of infections and the exposed cases considerably decreases. Consequently, these results can help to reduce the spread of the virus and to control the epidemic. In future work, the role of the vaccination on the expansion of this disease can be incorporated into the model. Thus, we will study a fractional model using the new definition of the fractional derivative with a non-singular kernel in the sense of Caputo [27].



---

*Source: 5163609-2022-09-26.xml*