# The Fractional SIRC Model and Influenza A

Authors: Moustafa El-Shahed; Ahmed Alsaedi
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This paper deals with the fractional-order SIRC model associated with the evolution of influenza A disease in human population. Qualitative dynamics of the model is determined by the basic reproduction number,R0. We give a detailed analysis for the asymptotic stability of disease-free and positive fixed points. Nonstandard finite difference methods have been used to solve and simulate the system of differential equations.


---

## Body

## 1. Introduction

Influenza is transmitted by a virus that can be of three different types, namely A, B, and C [1]. Among these, the virus A is epidemiologically the most important one for human beings, because it can recombine its genes with those of strains circulating in animal populations such as birds, swine, horses, and so forth [2, 3]. Over the last two decades, a number of epidemic models for predicting the spread of influenza through human population have been proposed based on either the classical susceptible-infected-removed (SIR) model developed by Kermack and McKendrick [4].Casagrandi et al. [5] have introduced SIRC model by adding a new compartment C, which can be called cross-immune compartment, to the SIR model. This cross-immune compartment (C) describes an intermediate state between the fully susceptible (S) and the fully protected (R) one. They have studied the dynamical behaviors of this model numerically [6]. Jódar et al. [7] developed two nonstandard finite difference schemes to obtain numerical solutions of a influenza A disease model presented by Casagrandi et al. [5]. Very recently Samanta [6] considered a nonautonomous SIRC epidemic model for Influenza A with varying total population size and distributed time delay.The notion of fractional calculus was anticipated by Leibniz, one of the founders of standard calculus, in a letter written in 1695. In recent decades, the fractional calculus and fractional differential equations have attracted much attention and increasing interest due to their potential applications in science and engineering (see [8, 9]).In this paper, we consider the fractional order SIRC model associated with the evolution of influenza A disease in human population. Qualitative dynamics of the model is determined by the basic reproduction number,R0. We give a detailed analysis for the asymptotic stability of disease-free and positive fixed points. Numerical simulations are presented to verify the obtained results.


## 2. Model Derivation

There are many definitions of fractional derivatives [8, 9]. Perhaps the best known is the Riemann-Liouvile definition. The Riemann-Liouville derivative of order α is defined as(2.1)RLD0+αf(t)=1Γ(n-α)(ddt)n∫0tf(s)(t-s)α-n+1ds,n=[α]+1,
where Γ is the gamma function and n is an integer. An alternative definition was introduced by Caputo as follows, which is a sort of regularization of the Riemann-Liouville derivative:(2.2)Dtαf(t)=1Γ(n-α)∫0tf(n)(s)(t-s)α-n+1ds.
The most common definition is the Caputo definition, since it is widely used in real applications. The initial conditions for the fractional order differential equations with the Caputo’s derivative are in the same form as for the integer-order differential equations. The Grunwald-Letnikov (GL) definition is given as(2.3)GLDtαf(t)=Limh→0h-α∑j=0[(t-a)/h](-1)j(αj)f(t-jh).
This formula can be reduced to(2.4)0Dtαy(tm)=h-α∑j=0mωj(α)ym-j,
where h is the time step and ωjα are the Grunwald-Letnikov coefficients defined as ωjα=(1-(1+α)/j)ωj-1α,j=0,1,2,…,and ω0α=h-α. The model presented in [5] for the spread influenza disease in the human population classifies the population in four groups or classes: S(t) is the proportion of susceptibles at time t (individuals that do not have the virus), I(t) is the proportion of infected at time t (individuals that have the virus and can infect), R(t) is the proportion of recovered at time t (individuals recovered from the virus and have a temporary immunity), and C(t) is the proportion of cross-immune individuals at time t. One of the main assumptions of this model is that the per capita birth rate is a constant μ>0 and the birth rate is the same as death rate. Using the above assumptions Casagrandi et al. [5] introduced the following system:
(2.5)dSdt=μ(1-S)-βSI+γC,dIdt=βSI+σβCI-(μ+θ)I,dRdt=(1-σ)βCI+θI-(μ+δ)R,dCdt=δR-βCI-(μ+γ)C,
where the parameter β is the contact rate for the influenza disease also called the rate of transmission for susceptible to infected, γ-1 is the cross-immune period, θ-1is the infectious period, δ-1is the total immune period and σ is the fraction of the exposed cross-immune individuals who are recruited in a unit time into the infective subpopulation [5, 7].Recently great considerations have been made to the models of FDEs in different area of researches. The most essential property of these models is their nonlocal property which does not exist in the integer order differential operators. We mean by this property that the next state of a model depends not only upon its current state but also upon all of its historical states. Now we introduce fractional order into the ODE model by Casagrandi et al. [5]. The new system is described by the following set of fractional order differential equations:
(2.6)DtαS=μ(1-S)-βSI+γC,DtαI=βSI+σβCI-(μ+θ)I,DtαR=(1-σ)βCI+θI-(μ+δ)R,DtαC=δR-βCI-(μ+γ)C,
whereDtα is the Caputo fractional derivative. Because model (2.6) monitors the dynamics of human populations, all the parameters are assumed to be nonnegative. Furthermore, it can be shown that all state variables of the model are nonnegative for all time t≥0 (see, for instance, [7, 10]).Lemma 2.1.
The closed setΩ={(S,I,R,C)∈R+4:S+I+R+C=1}is positively invariant with respect to model (2.6).Proof.
The fractional derivative of the total population, obtained by adding all the equations of model (2.6), is given by
(2.7)DtαN(t)=μ-μN(t).
The solution to (2.7) is given by N(t)=N(0)Eα,1(-μtα)+μtαEα,α+1(-μtα), where Eα,βis the Mittag-Leffler function. Considering the fact that the Mittag-Leffler function has an asymptotic behavior [9, 11],
(2.8)E∝,β(z)~-∑k=1pz-kΓ(β-αk)+O(|z|-1-p),(|z|⟶∞,απ2<|arg⁡(z)|≤π).
One can observe that N(t)converges to 1 when t→∞. Therefore, all solutions of the model with initial conditions in Ω remain in Ω for all t>0. Thus, region Ω is positively invariant with respect to model (2.6).In the following, we will study the dynamics of system (2.6).


## 3. Equilibrium Points and Stability

To evaluate the equilibrium points let(3.1)DtαS=0,DtαI=0,DtαR=0,DtαC=0.
Then E0=(1,0,0,0). By (2.6), a positive equilibrium E1=(S1,I1,R1,C1) satisfies (3.2)S1=μ+θβ-σ(δθI1(μ+σδ)βI1+(μ+γ)(μ+δ)),R1=θI1(βI1+μ+γ)(μ+σδ)βI1+(μ+γ)(μ+δ),C1=δθI1(μ+σδ)βI1+(μ+γ)(μ+δ),
andI1 is the positive root of g(I)=A1I2+A2I+A3(1-R0), where R0=β/(μ+θ) and
(3.3)A1=βμ(θ+μ+δσ),A2=βμ(-βμ+γ(δ+θ+μ)+(θ+μ)(δ+2μ)+δ(-β+μ)σ),A3=μ(μ+γ)(μ+δ)(μ+θ).
The Jacobian matrix J(E0) for the system given in (2.6) evaluated at the disease-free equilibrium is as follows:(3.4)J(E0)=(-μ-β0γ0β-μ-θ000θ-μ-δ000δ-μ-γ).Theorem 3.1.
The disease-free equilibriumE0 is locally asymptotically stable if R0<1 and is unstable if R0>1.Proof.
The disease-free equilibrium is locally asymptotically stable if all the eigenvalues,λi,i=1,2,3,4of the Jacobian matrix J(E0)satisfy the following condition [12, 13]:
(3.5)|arg(λi)|>απ2.
The eigenvalues of the Jacobian matrix J(E0)areλ1=-μ,λ2=-(μ+δ),λ3=-(μ+γ) and λ4=(μ+θ)(R0-1). Hence E0is locally asymptotically stable if R0<1and is unstable if R0>1. We now discuss the asymptotic stability of the endemic (positive) equilibrium of the system given by (2.6). The Jacobian matrix J(E1)evaluated at the endemic equilibrium is given as:
(3.6)J(E1)=(-μ-βI1-βS10γβI1βS1+σβC1-(μ+θ)0σβI10(1-σ)βC1+θ-(μ+δ)(1-σ)βI10-βC1δ-βI1-(μ+γ)).
The characteristic equation of J(E1)is
(3.7)(λ+μ)(λ3+a1λ2+a2λ+a3)=0,
where
(3.8)a1=γ+δ+2μ+2βI1,a2=(γ+μ)(δ+μ)+βI1(γ+δ+θ+3μ+δσ+βI1),a3=βI1(Ψ((δ+μ)(θ+μ)+γ(δ+θ+μ)-δθσ)+β(θ+μ+δσ)I1(2Ψ+β(μ+δσ)I1))Ψ+β(μ+δσ)I1,
where Ψ=(μ+γ)(μ+γ). Let D(Φ) denote the discriminant of a polynomial f. If Φ(x)=x3+a1x2+a2x+a3then. Denote
(3.9)D(Φ)=-|1a1a2a3001a1a2a332a1a20003a1a200032a1a2|=18a1a2a3+(a1a2)2-4a3a13-4a23-27a32.
Following [12, 14–16], we have the proposition.Proposition 3.2.
One assume thatE1 exists in R+3. (i)
If the discriminant ofΦ(x), D(Φ) is positive and Routh-Hurwitz are satisfied, that is, D(Φ)>0,a1>0, a3>0, a1a2>a3, then E1 is locally asymptotically stable.(ii)
IfD(Φ)<0, a1>0, a2>0, a1a2=a3, α∈[0,1) then E1 is locally asymptotically stable.(iii)
IfD(Φ)<0, a1<0, a2<0,α>2/3, then E1 is unstable.


## 4. Numerical Methods and Simulations

Since most of the fractional-order differential equations do not have exact analytic solutions, so approximation and numerical techniques must be used. Several analytical and numerical methods have been proposed to solve the fractional-order differential equations. For numerical solutions of the system (2.6) one can use the nonstandard finite difference method (NFDM). The nonstandard finite difference schemes were introduced by Mickens in the 1980s as a powerful numerical method that preserves significant properties of exact solutions of the involved differential equation [17]. The concept of the nonstandard finite difference method is discussed in [18]. Applying this method, the system (2.6) can be discretized as follows [7]: (4.1)∑j=0k+1ωjαSk+1-j=μ(1-Sk+1)-βSk+1Ik+γCk,∑j=0k+1ωjαIk+1-j=βSk+1Ik+σβCk+1Ik-(μ+θ)Ik+1,∑j=0k+1ωjαRk+1-j=(1-σ)βCk+1Ik+θIk+1-μRk+1-δRk,∑j=0k+1ωjαCk+1-j=δRk-βCk+1Ik-μCk+1-γCk.
Doing some algebraic manipulation to (4.1) yields the following relations:(4.2)Sk+1=μ+γCk-∑j=1k+1ωjαSk+1-jω0α+μ+βIk,Ik+1=βSk+1Ik+σβCk+1Ik-∑j=1k+1ωjαIk+1-jω0α+μ+θ,Rk+1=(1-σ)βCk+1Ik+θIk+1-δRk-∑j=1k+1ωjαRk+1-jω0α+μ,Ck+1=δRk-γCk-∑j=1k+1ωjαCk+1-jω0α+μ+βIk.
The approximate solutions S(t), I(t), R(t), and C(t) are displayed in Figures 1, 2, 3, 4, and 5, respectively. In each figure three different values of α are considered. When α=1, system (2.6) is the classical integer-order system (2.5). Figure 1 indicates behavior of the approximate solutions for system (2.6) obtained for the values of α=0.9. In Figure 2, the variation of S(t) versus time t is shown for different values of α=1,0.9,0.8 by fixing other parameters. It is revealed that increase in α increases with the proportion of susceptible while behavior is reverse after certain value of time. Figure 3 depicts I(t) versus time t, for Figures 3, 4, and 5 showing the similar variations of I(t), R(t), and C(t) with various values of α.Figure 1
S(t), I(t), R(t), C(t) for μ=0.02, β=100, δ=1, γ=0.5, σ=0.05, θ=73, α=0.9.Figure 2
S(t) for μ=0.02, β=100, δ=1, γ=0.5, σ=0.05, θ=73.Figure 3
I(t) for μ=0.02, β=100, δ=1, γ=0.5, σ=0.05,θ=73.Figure 4
R(t) for μ=0.02, β=100, δ=1, γ=0.5, σ=0.05, θ=73.Figure 5
C(t) for μ=0.02, β=100, δ=1, γ=0.5, σ=0.05, θ=73.



---

*Source: 480378-2011-11-15.xml*