# On Fractional SIRC Model withSalmonella Bacterial Infection

Authors: Fathalla A. Rihan; Dumitru Baleanu; S. Lakshmanan; R. Rakkiyappan
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose a fractional order SIRC epidemic model to describe the dynamics ofSalmonella bacterial infection in animal herds. The infection-free and endemic steady sates, of such model, are asymptotically stable under some conditions. The basic reproduction number ℛ0 is calculated, using next-generation matrix method, in terms of contact rate, recovery rate, and other parameters in the model. The numerical simulations of the fractional order SIRC model are performed by Caputo’s derivative and using unconditionally stable implicit scheme. The obtained results give insight to the modelers and infectious disease specialists.


---

## Body

## 1. Introduction

During the past three decades, the subject of fractional calculus has gained popularity and importance, mainly due to its demonstrated applications in numerous diverse and widespread fields of science and engineering. For example, fractional calculus has been successfully applied to system biology, physics, chemistry and biochemistry, hydrology, medicine, and finance (see, e.g., [1–10] and the references therein). In many cases, the fractional order differential/integral equations models are more consistent with the real phenomena than the integer-order models because the fractional derivatives and integrals enable the description of the memory and hereditary properties inherent in various materials and processes. Hence, there is a growing need to study and use the fractional order differential and integral equations. However, analytical and closed solutions of these types of fractional equations cannot generally be obtained. As a consequence, approximate and numerical techniques are playing important role in identifying the solution behavior of such fractional equations and exploring their applications (see, e.g., [9, 11, 12] and the references therein).We recall that theSalmonella infection is a major zoonotic disease which is transmitted between humans and other animals. Most persons infected withSalmonella develop diarrhea, fever, and abdominal cramps 12 to 72 hours after infection. The illness usually lasts 4 to 7 days, and most persons recover without treatment. However, in some persons, the diarrhea may be so severe that the patient needs to be hospitalized.Salmonella live in the intestinal tracts of humans and other animals, including birds.Salmonella are usually transmitted to humans by eating foods contaminated with animal feces. Contaminated foods usually look and smell normal. Contaminated foods are often of animal origin, such as beef, poultry, milk, or eggs, but any food, including vegetables, may become contaminated [13]. Therefore,Salmonella is considered as a serious problem for the public health throughout the world. There are no doubts that mathematical modeling ofSalmonella bacterial infection plays an important role in gaining understanding of the transmission of the disease in specific environment and in predicting the behavior of any outbreak. Furthermore, mathematical analysis leads to determining the nature of equilibrium states and to suggesting recommended actions to be taken by decision makers to control the spreading of the disease. The objective of this work is to adopt the fractional order epidemic model to describe the dynamics ofSalmonella infections in animal herds.A large amount of work done on modelling biological systems has been restricted to integer-order ordinary (or delay) differential equations (see, e.g., [14–19]). In [20], the authors proposed the classicalSusceptible-Infected-Recovered (SIR) model. The authors in [21] introduced a new compartment into SIR model, which is called cross-immune compartment to be called SIRC model. The new compartment cross-immune C(t) describes an intermediate state between the fully susceptible S(t) and the fully protected R(t) one. Recently, the fractional order SIRC model of influenza, a disease in human population, was discussed in [22]. In the present paper, we consider the fractional order SIRC model associated with evolution ofSalmonella bacterial infection in animal herds. However, we will take into account the disease-induced mortality rate m in the model. Qualitative behavior of the fractional order SRIC model is then investigated. Numerical simulations of the fractional order SRIC model are provided to demonstrate the effectiveness of the proposed method by using implicit Euler’s method.We first give the definition of fractional order integration and fractional order differentiation [23–25]. Let L1=L1[a,b] be the class of Lebesgue integrable functions on [a,b], a<b<∞.Definition 1.
The fractional integral of orderν∈ℝ+ of the function f(t), t>0 (f:ℝ+→ℝ) is defined by
(1)Iaνf(t)=1Γ(ν)∫at(t-s)ν-1f(s)ds‍,t>0.
However, the fractional derivative of order α∈(n-1,n) of f(t) is defined by two ways.(i)
Riemann-Liouville fractional derivative: take fractional integral of order(n-α) and then take nth derivative:
(2)Daαf(t)=DanIan-αf(t),D*n=dndtn,n=1,2,….(ii)
Caputo’s fractional derivative: takenth derivative and then take a fractional integral of order (n-α):
(3)Daαf(t)=Ian-αDanf(t),n=1,2,….We notice that the definition of time-fractional derivative of a functionf(t) at t=tn involves an integration and calculating time-fractional derivative that requires all the past history, that is, all the values of f(t) from t=0 to t=tn. Caputo’s definition, which is a modification of the Riemann-Liouville definition, has the advantage of dealing properly with initial value problems. For more properties of the fractional derivatives and integrals, we refer to [8, 9, 24, 25] and references therein. The generalized mean value theorem is defined in the following Remark [26].Remark 2.
(i) Suppose thatf(t)∈C[a,b] and D*αf(t)∈C(a,b] for 0<α≤1; then we have
(4)f(t)=f(a)+1Γ(α)D*αf(ξ)(t-a)α,hhhwitha<ξ<t∀t∈(a,b].
(ii) If (i) holds andD*αf(t)≥0∀t∈[a,b], then f(t) is nondecreasing for each t∈[a,b]. If D*αf(t)≤0∀t∈[a,b], then f(t) is nonincreasing for each t∈[a,b].


## 2. The SIRC Epidemic Model

Assume that theSalmonella infection spreads in animal herds which are grouped as four compartments, according to their infection status: S(t) is the proportion of susceptible individuals at time t (individuals that do not have the bacterial infection), I(t) is the proportion of infected individuals (that have the bacterial infection), R(t) is the proportion of recovered individuals (that recovered from the infection and have temporary immunity), and C(t) is the proportion of cross-immune individuals at time t. The total number of animals in the herd is given by N=S+I+R+C. We consider that initially all the animals are susceptible to the infection. Once infected, a susceptible individual leaves the susceptible compartment and enters the infectious compartment where it then becomes infectious. The infected animals pass into the recovered compartment. The individuals who have recovered from the disease have temporary immunity and grouped into C(t) compartment. Therefore, we consider that the disease transmission model consists of nonnegative initial conditions together with system of equations:
(5)S˙(t)=μN+ηC(t)-(βI(t)+μ)S(t),I˙(t)=βS(t)I(t)+σβC(t)I(t)-(θ+m+μ)I(t),R˙(t)=(1-σ)βC(t)I(t)+θI(t)-(μ+δ)R(t),C˙(t)=δR(t)-βC(t)I(t)-(η+μ)C(t).
Here, parameter μ denotes the mortality rate in every compartment and is assumed to be equal to the rate of newborns in the population. β is the contact rate and also called transmission from susceptible to infected. η-1 is the cross-immune period, θ-1 is the infectious period, δ-1 is 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 [21, 27]. We also assume that the disease induces mortality rate m; see the diagram of Figure 1.Figure 1
Schematic diagram of SIRC model.


### 2.1. Fractional Order of SIRC Epidemic Model

Although a large number of work has been done in modeling the dynamics of epidemiological diseases, it has been restricted to integer-order (delay) differential equations. In recent years, it has turned out that many phenomena in different fields can be described very successfully by models usingfractional order differential equations (FODEs) [1, 5, 28]. Now, we introduce fractional order into model (5) and assume that s(t)=S(t)/N, i(t)=S(t)/N, r(t)=R(t)/N, c(t)=C(t)/N, where N is the total number of population. Then the model takes the form
(6)Dα1s(t)=μ+ηc(t)-(βi(t)+μ)s(t),Dα2i(t)=βs(t)i(t)+σβc(t)i(t)-(θ+m+μ)i(t),Dα3r(t)=(1-σ)βc(t)i(t)+θi(t)-(μ+δ)r(t),Dα4c(t)=δr(t)-βc(t)i(t)-(η+μ)c(t)
with initial conditions s(0)=s0, i(0)=i0, and r(0)=r0.


### 2.2. Stability Criteria for the Fractional Order SIRC Model

In model (6), assume that α1=α2=α3=α4=α. Then, to find the equilibria, we put Dαs=Dαi=Dαr=Dαc=0.Consider(7)ℰ0=(1,0,0,0),ℰ+=(s*,i*,r*,c*),
where
(8)s*=θ+m+μβ-σ(δθi*(μ+δσ)βi*+(μ+δ)(μ+η)),r*=θi*(βi*+η+μ)(μ+δσ)βi*+(μ+δ)(μ+η),c*=θδi*(μ+δσ)βi*+(μ+δ)(μ+η).
The positive endemic equilibrium ℰ+=(s*,i*,r*,c*) satisfies (6) and i* is the positive root of A1i*2+A2i*+A3, where
(9)A1=-β2[m(μ+δσ)+μ(θ+μ+δσ)],A2=β[βμ(μ+δσ)+ηθδ-(θ+m+μ)hhβμ×[(μ+δ)(μ+η)+μ(μ+δσ)]+μσδθ],A3=βμ(μ+δ)(μ+η)[1-(θ+m+μβ)].
The Jacobian matrix of model (6) is(10)J=(-βi(t)-μ-βs(t)0ηβi(t)βs(t)+σβc(t)-(θ+m+μ)0σβi(t)0(1-σ)βc(t)+θ-(μ+δ)(1-σ)βi(t)0-βc(t)δ-βi(t)-(η+μ)).


### 2.3. The Reproduction Numberℛ0

The basic reproduction number (the number of individuals infected by a single infected individual placed in a totally susceptible population)ℛ0 that includes the indirect transmission may be obtained using next-generation matrix method. The spectral radius of the next-generation matrix (FV-1), which is the dominant eigenvalue of the same matrix, gives the value of ℛ0 [29]. Then, the basic reproductive number ℛ0 is obtained by the form
(11)ℛ0=ρ(FV-1),
where the matrices F=[∂ℱi(x)/∂xj]x=x0 and V=[∂𝒱i(x)/∂xj]x=x0. ℱi(x), where x is the set of all disease free states in the compartment i, is the rate of appearance of new infections in compartment i, and Vi(x) is net transfer rate (other than infections) of compartment i. The net transfer rate is given by 𝒱i=𝒱i--𝒱i+, where 𝒱i- is the rate of transfer of individuals out of compartment i, and 𝒱i+ is the rate of transfer of individuals into compartment i by all other means. Therefore, the disease transmission model consists of nonnegative initial conditions, xi(0), together with the following system of equations:
(12)xj′=fj(x)≡ℱj(x)-𝒱j,j≥1.
From model (6), we have
(13)F=(∂ℱ1∂i(t)∂ℱ1∂r(t)∂ℱ2∂i(t)∂ℱ2∂r(t))=(βs000),V=(∂𝒱1∂i(t)∂𝒱1∂r(t)∂𝒱2∂i(t)∂𝒱2∂r(t))=(θ+m+μ0-θμ+δ).
Since we have only two distinct stages, namely, I(t) and R(t), it follows that both F and V are 2×2 square matrices. Further, it can be noticed that F is nonnegative and V is nonsingular. The basic reproductive number ℛ0 is the dominant eigenvalue of the matrix FV-1, which is obtained by solving the characteristic equation (FV-1)I-ΛI=0, where Λ is the eigenvalue and I(t) is the identity matrix. At the disease-free equilibrium ℰ0=(1,0,0,0), we have
(14)ℛ0=βθ+m+μ.
The following theorem states that ℛ0 is a threshold parameter for the stability of model (6).Theorem 3.
The disease-free equilibrium is locally asymptotically stable and the infection will die out ifℛ0<1 and is unstable if ℛ0>1.Proof.
The disease-free equilibrium is locally asymptotically stable if all the eigenvalues,λi, i=1,2,3,4, of Jacobian matrix J(ℰ0) satisfy condition [30]:
(15)|arg(λi)|>απ2,
where
(16)J(ℰ0)=(-μ-β0η0β-(θ+m+μ)0000-(μ+δ)000δ-(η+μ)).
Figure 2 depicts the stability region of the fractional order system, according to condition (15). The eigenvalues of Jacobian matrix J(ℰ0) are λ1=-μ, λ2=β-(θ+m+μ), λ3=-(μ+δ), λ4=-(η+μ). Hence, ℰ0 is locally asymptotically stable if ℛ0<1 and is unstable if ℛ0>1.
Now, we extend the analysis to endemic equilibriumℰ+. Jacobian matrix J(ℰ+) evaluated at the endemic equilibrium is(17)J(ℰ+)=(-βi*-μ-βs*0ηβi*βs*+σβc*-(θ+m+μ)0σβi*0(1-σ)βc*+θ-(μ+δ)(1-σ)βi*0-βc*δ-βi*-(η+μ)),with characteristic equation
(18)λ4+a1λ3+a2λ2+a3λ+a4=0,
where
(19)a1=(D1+D3+D5),a2=(D1D3-D4δ+D1D5+D3D5+β2i*s*+σβ2c*i*),a3=(D1D3D5-D1D4δ+D3β2i*s*+D5β2i*s*+β2c*ηi*h-D2σβδi*+σβ2D1c*i*+σD3β2c*i*),a4=D3D5β2i*s*-D2βδηi*+D3β2c*ηi*-D4β2δi*s*-σβδD1D2i*+σD1D3β2c*i*,D1=βi+μ,D2=(1-σ)βc*+θ,D3=(μ+δ),D4=(1-σ)βi*,D5=βi*+(η+μ).
If D(Φ) denotes the discriminant of the polynomial: Φ(λ)=λ4+a1λ3+a2λ2+a3λ+a4, then denote
(20)D(Φ)=|1a1a2a3a40001a1a2a3a40001a1a2a3a443a12a2a3000043a12a2a3000043a12a2a3000043a12a2a3|.
From [31], we have the proposition.
Proposition 4.Assume thatℰ+exists inR+4.(1)
Letc1,c2,c3be the Routh-Hurwitz determinants:c1=a1,c2=|a11a3a2|,c3=|a110a3a2a10a4a3|. 											Therefore, whenα=1, the equilibrium pointℰ+is locally asymptotically stable if(21)c1>0,c2>0,c3=0,a4>0.However, conditions  (21)  are sufficient (not necessary) conditions forℰ+to be locally asymptotically stable for allα∈[0,1).(2)
IfD(Φ)>0,a1>0,a2<0, andα>2/3, then the equilibrium pointℰ+is unstable.(3)
IfD(Φ)<0,a1>0,a2>0,a3>0,a4>0, andα<1/3, then the equilibriumℰ+is locally asymptotically stable. Also, ifD(Φ)<0,a1<0,a2>0,a3<0,a4>0, then the equilibrium pointℰ+is unstable.(4)
IfD(Φ)<0,a1>0,a2>0,a3>0,a4>0, anda2=(a1a4/a3)+(a3/a1), then the equilibrium pointℰ+is locally asymptotically stable for allα∈(0,1).(5)
a4>0is the necessary condition for the equilibrium pointℰ+to be locally asymptotically stable.Figure 2
Stability region of the fractional order system withα∈(0,1].


## 2.1. Fractional Order of SIRC Epidemic Model

Although a large number of work has been done in modeling the dynamics of epidemiological diseases, it has been restricted to integer-order (delay) differential equations. In recent years, it has turned out that many phenomena in different fields can be described very successfully by models usingfractional order differential equations (FODEs) [1, 5, 28]. Now, we introduce fractional order into model (5) and assume that s(t)=S(t)/N, i(t)=S(t)/N, r(t)=R(t)/N, c(t)=C(t)/N, where N is the total number of population. Then the model takes the form
(6)Dα1s(t)=μ+ηc(t)-(βi(t)+μ)s(t),Dα2i(t)=βs(t)i(t)+σβc(t)i(t)-(θ+m+μ)i(t),Dα3r(t)=(1-σ)βc(t)i(t)+θi(t)-(μ+δ)r(t),Dα4c(t)=δr(t)-βc(t)i(t)-(η+μ)c(t)
with initial conditions s(0)=s0, i(0)=i0, and r(0)=r0.


## 2.2. Stability Criteria for the Fractional Order SIRC Model

In model (6), assume that α1=α2=α3=α4=α. Then, to find the equilibria, we put Dαs=Dαi=Dαr=Dαc=0.Consider(7)ℰ0=(1,0,0,0),ℰ+=(s*,i*,r*,c*),
where
(8)s*=θ+m+μβ-σ(δθi*(μ+δσ)βi*+(μ+δ)(μ+η)),r*=θi*(βi*+η+μ)(μ+δσ)βi*+(μ+δ)(μ+η),c*=θδi*(μ+δσ)βi*+(μ+δ)(μ+η).
The positive endemic equilibrium ℰ+=(s*,i*,r*,c*) satisfies (6) and i* is the positive root of A1i*2+A2i*+A3, where
(9)A1=-β2[m(μ+δσ)+μ(θ+μ+δσ)],A2=β[βμ(μ+δσ)+ηθδ-(θ+m+μ)hhβμ×[(μ+δ)(μ+η)+μ(μ+δσ)]+μσδθ],A3=βμ(μ+δ)(μ+η)[1-(θ+m+μβ)].
The Jacobian matrix of model (6) is(10)J=(-βi(t)-μ-βs(t)0ηβi(t)βs(t)+σβc(t)-(θ+m+μ)0σβi(t)0(1-σ)βc(t)+θ-(μ+δ)(1-σ)βi(t)0-βc(t)δ-βi(t)-(η+μ)).


## 2.3. The Reproduction Numberℛ0

The basic reproduction number (the number of individuals infected by a single infected individual placed in a totally susceptible population)ℛ0 that includes the indirect transmission may be obtained using next-generation matrix method. The spectral radius of the next-generation matrix (FV-1), which is the dominant eigenvalue of the same matrix, gives the value of ℛ0 [29]. Then, the basic reproductive number ℛ0 is obtained by the form
(11)ℛ0=ρ(FV-1),
where the matrices F=[∂ℱi(x)/∂xj]x=x0 and V=[∂𝒱i(x)/∂xj]x=x0. ℱi(x), where x is the set of all disease free states in the compartment i, is the rate of appearance of new infections in compartment i, and Vi(x) is net transfer rate (other than infections) of compartment i. The net transfer rate is given by 𝒱i=𝒱i--𝒱i+, where 𝒱i- is the rate of transfer of individuals out of compartment i, and 𝒱i+ is the rate of transfer of individuals into compartment i by all other means. Therefore, the disease transmission model consists of nonnegative initial conditions, xi(0), together with the following system of equations:
(12)xj′=fj(x)≡ℱj(x)-𝒱j,j≥1.
From model (6), we have
(13)F=(∂ℱ1∂i(t)∂ℱ1∂r(t)∂ℱ2∂i(t)∂ℱ2∂r(t))=(βs000),V=(∂𝒱1∂i(t)∂𝒱1∂r(t)∂𝒱2∂i(t)∂𝒱2∂r(t))=(θ+m+μ0-θμ+δ).
Since we have only two distinct stages, namely, I(t) and R(t), it follows that both F and V are 2×2 square matrices. Further, it can be noticed that F is nonnegative and V is nonsingular. The basic reproductive number ℛ0 is the dominant eigenvalue of the matrix FV-1, which is obtained by solving the characteristic equation (FV-1)I-ΛI=0, where Λ is the eigenvalue and I(t) is the identity matrix. At the disease-free equilibrium ℰ0=(1,0,0,0), we have
(14)ℛ0=βθ+m+μ.
The following theorem states that ℛ0 is a threshold parameter for the stability of model (6).Theorem 3.
The disease-free equilibrium is locally asymptotically stable and the infection will die out ifℛ0<1 and is unstable if ℛ0>1.Proof.
The disease-free equilibrium is locally asymptotically stable if all the eigenvalues,λi, i=1,2,3,4, of Jacobian matrix J(ℰ0) satisfy condition [30]:
(15)|arg(λi)|>απ2,
where
(16)J(ℰ0)=(-μ-β0η0β-(θ+m+μ)0000-(μ+δ)000δ-(η+μ)).
Figure 2 depicts the stability region of the fractional order system, according to condition (15). The eigenvalues of Jacobian matrix J(ℰ0) are λ1=-μ, λ2=β-(θ+m+μ), λ3=-(μ+δ), λ4=-(η+μ). Hence, ℰ0 is locally asymptotically stable if ℛ0<1 and is unstable if ℛ0>1.
Now, we extend the analysis to endemic equilibriumℰ+. Jacobian matrix J(ℰ+) evaluated at the endemic equilibrium is(17)J(ℰ+)=(-βi*-μ-βs*0ηβi*βs*+σβc*-(θ+m+μ)0σβi*0(1-σ)βc*+θ-(μ+δ)(1-σ)βi*0-βc*δ-βi*-(η+μ)),with characteristic equation
(18)λ4+a1λ3+a2λ2+a3λ+a4=0,
where
(19)a1=(D1+D3+D5),a2=(D1D3-D4δ+D1D5+D3D5+β2i*s*+σβ2c*i*),a3=(D1D3D5-D1D4δ+D3β2i*s*+D5β2i*s*+β2c*ηi*h-D2σβδi*+σβ2D1c*i*+σD3β2c*i*),a4=D3D5β2i*s*-D2βδηi*+D3β2c*ηi*-D4β2δi*s*-σβδD1D2i*+σD1D3β2c*i*,D1=βi+μ,D2=(1-σ)βc*+θ,D3=(μ+δ),D4=(1-σ)βi*,D5=βi*+(η+μ).
If D(Φ) denotes the discriminant of the polynomial: Φ(λ)=λ4+a1λ3+a2λ2+a3λ+a4, then denote
(20)D(Φ)=|1a1a2a3a40001a1a2a3a40001a1a2a3a443a12a2a3000043a12a2a3000043a12a2a3000043a12a2a3|.
From [31], we have the proposition.
Proposition 4.Assume thatℰ+exists inR+4.(1)
Letc1,c2,c3be the Routh-Hurwitz determinants:c1=a1,c2=|a11a3a2|,c3=|a110a3a2a10a4a3|. 											Therefore, whenα=1, the equilibrium pointℰ+is locally asymptotically stable if(21)c1>0,c2>0,c3=0,a4>0.However, conditions  (21)  are sufficient (not necessary) conditions forℰ+to be locally asymptotically stable for allα∈[0,1).(2)
IfD(Φ)>0,a1>0,a2<0, andα>2/3, then the equilibrium pointℰ+is unstable.(3)
IfD(Φ)<0,a1>0,a2>0,a3>0,a4>0, andα<1/3, then the equilibriumℰ+is locally asymptotically stable. Also, ifD(Φ)<0,a1<0,a2>0,a3<0,a4>0, then the equilibrium pointℰ+is unstable.(4)
IfD(Φ)<0,a1>0,a2>0,a3>0,a4>0, anda2=(a1a4/a3)+(a3/a1), then the equilibrium pointℰ+is locally asymptotically stable for allα∈(0,1).(5)
a4>0is the necessary condition for the equilibrium pointℰ+to be locally asymptotically stable.Figure 2
Stability region of the fractional order system withα∈(0,1].


## 3. Implicit Euler’s Scheme for FODEs

Since most of the FODEs do not have exact analytic solutions, so approximation and numerical techniques must be used. In addition, most of resulting biological systems are stiff (one definition of the stiffness is that the global accuracy of the numerical solution is determined by stability rather than local error and implicit methods are more appropriate for it). Thestiffness often appears due to the differences in speed between the fastest and slowest components of the solutions and stability constraints. In addition, the state variables of these types of models are very sensitive to small perturbations (or changes) in the parameters which occur in the model. Therefore, efficient use of a reliable numerical method for dealing with stiff problems is necessary.Consider the following fractional order differential equation:(22)Dαy(t)=f(t,y(t)),t∈[0,T],y(k)(0)=y(k)(0),k=0,1,2,…,m-1,0<α≤1.
Here, y(t)=[y1(t),y2(t),…,yn(t)]T and f(t,y(t)) satisfy the Lipschitz condition in variable y:
(23)∥f(t,y(t))-f(t,x(t))∥≤K∥y(t)-x(t)∥,K>0,
where x(t) is the solution of the perturbed system.Theorem 5.
Problem (22) has a unique solution provided that Lipschitz condition (23) is satisfied and M-=KTα/Γ(α+1)<1.Proof.
Using the definitions of Section1, we can apply a fractional integral operator to the differential equation (22) and incorporate the initial conditions, thus converting the equation into the equivalent equation:
(24)y(t)=∑k=0m-1y0(k)tkk!+1Γ(α)∫0t(t-s)α-1f(s,y(s))ds‍‍,
which also is a Volterra equation of the second kind. Define operator ℒ, such that
(25)ℒy(t)=∑k=0m-1y0(k)tkk!+1Γ(α)∫0t‍(t-s)α-1f(s,y(s))ds‍.
Then, we have
(26)∥ℒy(t)-ℒx(t)∥≤1Γ(α)∫0t‍(t-s)α-1∥f(s,y(s))-f(s,x(s))∥ds≤KΓ(α)∫0t(t-s)α-1sup⁡s∈[0,T]|y(s)-x(s)|ds‍≤KΓ(α)∥y-x∥∫0t‍sα-1ds≤KTαΓ(α+1)∥y-x∥Tα.
So, we obtain
(27)∥ℒy(t)-ℒx(t)∥≤M-∥y-x∥.
By Banach contraction principle [32], we can deduce that ℒ has a unique fixed point which implies that our problem has a unique solution.Several numerical methods have been proposed to solve the FODEs [11, 33]. Recently, the predictor-corrector algorithm is an efficient and powerful technique for solving the FODEs, which is a generalization of the Adams-Bashforth-Moulton method. The modification of Adams-Bashforth-Moulton algorithm is proposed by Diethelm [34, 35] to approximate the fractional order derivative. However, converted Volterra integral equation (24) is with a weakly singular kernel, such that a regularization is not necessary any more. It seems that there exist only a very small number of software packages for nonlinear Volterra equations. In our case, the kernel may not be continuous, and therefore the classical numerical algorithms for the integral part of (24) are unable to handle the solution of (22). Therefore, we implement the implicit Euler’s scheme to approximate the fractional order derivative.Given model (22) and mesh points 𝒯={t0,t1,…,tN}, such that t0=0 and tN=T, then a discrete approximation to the fractional derivative can be obtained by a simple quadrature formula, using Caputo’s fractional derivative (3) of order α, 0<α≤1 and usingimplicit Euler’s approximation as follows (see [12]):
(28)D*αxi(tn)=1Γ(1-α)∫0tdxi(s)ds(tn-s)-αds‍≈1Γ(1-α)×∑j=1n‍∫(j-1)hjh[xij-xij-1h+O(h)](nh-s)-αds‍=1(1-α)Γ(1-α)×∑j=1n{[xij-xij-1h+O(h)]‍hhhhh[xij-xij-1h+O(h)]×[(n-j+1)1-α-(n-j)1-α]}h1-α=1(1-α)Γ(1-α)1hα×∑j=1n[xij-xij-1][(n-j+1)1-α-(n-j)1-α]‍+1(1-α)Γ(1-α)×∑j=1n[xij-xij-1]‍hhhh×[(n-j+1)1-α-(n-j)1-α]O(h2-α).
Setting
(29)𝒢(α,h)=1(1-α)Γ(1-α)1hα,ωjα=j1-α-(j-1)1-α,(whereω1α=1),
then the first-order approximation method for the computation of Caputo’s fractional derivative is then given by the expression
(30)D*αxi(tn)=𝒢(α,h)∑j=1nωjα(xin-j+1-xin-j)‍+O(h).
From the analysis and numerical approximation, we also arrive at the following proposition.Proposition 6.
The presence of a fractional differential order in a differential equation can lead to a notable increase in the complexity of the observed behavior, and the solution continuously depends on all the previous states.


### 3.1. Stability and Convergence

We here prove that fractional order implicit difference approximation (30) is unconditionally stable. It follows then that the numerical solution converges to the exact solution as h→0. In order to study the stability of the numerical method, let us consider a test problem of linear scaler fractional differential equation
(31)D*αu(t)=ρ0u(t)+ρ1,U(0)=U0
such that 0<α≤1, and ρ0<0, ρ1>0 are constants.Theorem 7.
The fully implicit numerical approximation (30), to test problem (31) for all t≥0, is consistent and unconditionally stable.Proof.
We assume that the approximate solution of (31) is of the form u(tn)≈Un≡ζn; then (31) can be reduced to
(32)(1-ρ0𝒢α,h)ζn=ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1)+ρ1𝒢α,h,n≥2,
or
(33)ζn=ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)‍+ρ1/𝒢α,h(1-(ρ0/𝒢α,h)),n≥2.
Since (1-(ρ0/𝒢α,h))≥1, for all 𝒢α,h, then
(34)ζ1≤ζ0,(35)ζn≤ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)‍,n≥2.
Thus, for n=2, the above inequality implies
(36)ζ2≤ζ1+ω2(α)(ζ0-ζ1).
Using relation (34) and the positivity of the coefficients ω2, we get
(37)ζ2≤ζ1.
Repeating the process, we have, from (35),
(38)ζn≤ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)‍≤ζn-1,
since each term in the summation is negative. Thus, ζn≤ζn-1≤ζn-2≤⋯≤ζ0. With the assumption that ζn=|Un|≤ζ0=|U0|, which entails ∥Un∥≤∥U0∥, we have stability.Of course this numerical technique can be used both for linear and for nonlinear problems, and it may be extended to multiterm FODEs.


### 3.2. Numerical Simulations

The approximate solutions of epidemic model (6) are displayed in Figures 3, 4, and 5, and sensitivity of ℛ0 to transmission coefficients is displayed in Figure 6. The numerical simulations are performed by Euler’s implicit scheme. We choose commensurate fractional order that α1=α2=α3=α4=α, with different fractional order values and the parameter values given in Table 1.Table 1
List of parameters.
Parameter
Description
Value
Reference
μ
Replacement and exit rate (day−1)
0.011
[36]
β
Contact (transmission) rate of suspectable to be infected (animal−1 day−1)
0.15
[36]
θ
Recovery rate of infected animals day-1
0.16
Assumed
m
Disease-induced mortality rate (day−1)
0.041
Assumed
η
Cross-immune period
0.5
[36]
σ
The average reinfection probability ofC(t)
0.06
Assumed
δ
The average time of appearance of new dominant clusters
1
Assumed
N
The total number of population
345
AssumedFigure 3
Numerical simulation of fractional order epidemic model (6), when α=0.8, and ℛ0>1 (each infected individual infects more than one other member of the population and a self-sustaining group of infectious individuals will propagate), with parameter values of Table 1.Phase plane portrait for fractional order endemic model (6), in absence of C(t) and R(t) components, when α=0.7 (a) and α=0.9 (b) with ℛ0=0.5<1. We note that solution paths approach the disease-free equilibrium ℰ0=(1,0,0).
(a)
(b)Phase plane portrait for classic fractional order endemic model (6) when α=1 (a) and α=0.9 (b) with ℛ0=1.2>1. We note that solution paths approach the endemic equilibrium ℰ+ given by (7).
(a)
(b)Sensitivity ofℛ0 with respect to the transmission coefficients β and θ.
(a)
(b)


## 3.1. Stability and Convergence

We here prove that fractional order implicit difference approximation (30) is unconditionally stable. It follows then that the numerical solution converges to the exact solution as h→0. In order to study the stability of the numerical method, let us consider a test problem of linear scaler fractional differential equation
(31)D*αu(t)=ρ0u(t)+ρ1,U(0)=U0
such that 0<α≤1, and ρ0<0, ρ1>0 are constants.Theorem 7.
The fully implicit numerical approximation (30), to test problem (31) for all t≥0, is consistent and unconditionally stable.Proof.
We assume that the approximate solution of (31) is of the form u(tn)≈Un≡ζn; then (31) can be reduced to
(32)(1-ρ0𝒢α,h)ζn=ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1)+ρ1𝒢α,h,n≥2,
or
(33)ζn=ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)‍+ρ1/𝒢α,h(1-(ρ0/𝒢α,h)),n≥2.
Since (1-(ρ0/𝒢α,h))≥1, for all 𝒢α,h, then
(34)ζ1≤ζ0,(35)ζn≤ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)‍,n≥2.
Thus, for n=2, the above inequality implies
(36)ζ2≤ζ1+ω2(α)(ζ0-ζ1).
Using relation (34) and the positivity of the coefficients ω2, we get
(37)ζ2≤ζ1.
Repeating the process, we have, from (35),
(38)ζn≤ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)‍≤ζn-1,
since each term in the summation is negative. Thus, ζn≤ζn-1≤ζn-2≤⋯≤ζ0. With the assumption that ζn=|Un|≤ζ0=|U0|, which entails ∥Un∥≤∥U0∥, we have stability.Of course this numerical technique can be used both for linear and for nonlinear problems, and it may be extended to multiterm FODEs.


## 3.2. Numerical Simulations

The approximate solutions of epidemic model (6) are displayed in Figures 3, 4, and 5, and sensitivity of ℛ0 to transmission coefficients is displayed in Figure 6. The numerical simulations are performed by Euler’s implicit scheme. We choose commensurate fractional order that α1=α2=α3=α4=α, with different fractional order values and the parameter values given in Table 1.Table 1
List of parameters.
Parameter
Description
Value
Reference
μ
Replacement and exit rate (day−1)
0.011
[36]
β
Contact (transmission) rate of suspectable to be infected (animal−1 day−1)
0.15
[36]
θ
Recovery rate of infected animals day-1
0.16
Assumed
m
Disease-induced mortality rate (day−1)
0.041
Assumed
η
Cross-immune period
0.5
[36]
σ
The average reinfection probability ofC(t)
0.06
Assumed
δ
The average time of appearance of new dominant clusters
1
Assumed
N
The total number of population
345
AssumedFigure 3
Numerical simulation of fractional order epidemic model (6), when α=0.8, and ℛ0>1 (each infected individual infects more than one other member of the population and a self-sustaining group of infectious individuals will propagate), with parameter values of Table 1.Phase plane portrait for fractional order endemic model (6), in absence of C(t) and R(t) components, when α=0.7 (a) and α=0.9 (b) with ℛ0=0.5<1. We note that solution paths approach the disease-free equilibrium ℰ0=(1,0,0).
(a)
(b)Phase plane portrait for classic fractional order endemic model (6) when α=1 (a) and α=0.9 (b) with ℛ0=1.2>1. We note that solution paths approach the endemic equilibrium ℰ+ given by (7).
(a)
(b)Sensitivity ofℛ0 with respect to the transmission coefficients β and θ.
(a)
(b)


## 4. Conclusions

In this paper, we provided a fractional order SIRC epidemic model withSalmonella bacteria infection. We derived the sufficient conditions to preserve the asymptotic stability of infection-free and endemic steady states. The threshold parameter (reproduction number) ℛ0 has been evaluated in terms of contact rate, recovery rate, and other parameters in the model. We provided unconditionally stable method, using Euler’s implicit method for the fractional order differential system. The solution of a fractional order model at any time t* continuously depends on all the previous states at t≤t*. Fractional order dynamical models are more suitable to model biological systems with memory than their integer-orders. The presence of a fractional differential order into a corresponding differential equation leads to a notable increase in the complexity of the observed behavior and enlarges the stability region of the solutions. However, fractional order differential models have the same integer-order counterpart steady states, when α>0.5.



---

*Source: 136263-2014-04-22.xml*