The Fractional SIRC Model and Influenza A

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, 𝑅₀. 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.

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 𝐶, which can be called cross-immune compartment, to the SIR model. This cross-immune compartment (𝐶) describes an intermediate state between the fully susceptible (𝑆) and the fully protected (𝑅) 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, 𝑅₀. 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_RL𝐷^𝛼₀₊𝑓(𝑡)=1Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑡₀𝑓(𝑠)(𝑡−𝑠)^{𝛼−𝑛+1}𝑑𝑠,𝑛=[𝛼]+1, (2.1) where Γ is the gamma function and 𝑛 is an integer. An alternative definition was introduced by Caputo as follows, which is a sort of regularization of the Riemann-Liouville derivative:𝐷^𝛼_𝑡𝑓(𝑡)=1Γ(𝑛−𝛼)^𝑡₀𝑓^(𝑛)(𝑠)(𝑡−𝑠)^{𝛼−𝑛+1}𝑑𝑠. (2.2) 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_GL𝐷^𝛼_𝑡𝑓(𝑡)=Lim_ℎ→0ℎ^{[()−𝛼𝑡−𝑎/ℎ]}_𝑗=0(−1)^𝑗𝛼𝑗𝑓𝑡−𝑗ℎ. (2.3) This formula can be reduced to₀𝐷^𝛼_𝑡𝑦𝑡_𝑚=ℎ^𝑚−𝛼_𝑗=0𝜔_𝑗^(𝛼)𝑦_𝑚−𝑗, (2.4) where ℎ is the time step and 𝜔^𝛼_𝑗 are the Grunwald-Letnikov coefficients defined as 𝜔^𝛼_𝑗=(1−(1+𝛼)/𝑗)𝜔^𝛼_𝑗−1,𝑗=0,1,2,…,and 𝜔^𝛼₀=ℎ^−𝛼. The model presented in [5] for the spread influenza disease in the human population classifies the population in four groups or classes: 𝑆(𝑡) is the proportion of susceptibles at time 𝑡 (individuals that do not have the virus), 𝐼(𝑡) is the proportion of infected at time 𝑡 (individuals that have the virus and can infect), 𝑅(𝑡) is the proportion of recovered at time 𝑡 (individuals recovered from the virus and have a temporary immunity), and 𝐶(𝑡) is the proportion of cross-immune individuals at time 𝑡. 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: 𝑑𝑆()𝑑𝑡=𝜇1−𝑆−𝛽𝑆𝐼+𝛾𝐶,𝑑𝐼𝑑𝑡=𝛽𝑆𝐼+𝜎𝛽𝐶𝐼−𝜇+𝜃𝐼,𝑑𝑅=()𝑑𝑡1−𝜎𝛽𝐶𝐼+𝜃𝐼−𝜇+𝛿𝑅,𝑑𝐶𝑑𝑡=𝛿𝑅−𝛽𝐶𝐼−𝜇+𝛾𝐶, (2.5) where the parameter 𝛽 is the contact rate for the influenza disease also called the rate of transmission for susceptible to infected, 𝛾⁻¹ is the cross-immune period, 𝜃⁻¹is the infectious period, 𝛿⁻¹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 [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: 𝐷^𝛼_𝑡()𝐷𝑆=𝜇1−𝑆−𝛽𝑆𝐼+𝛾𝐶,_𝛼𝑡𝐷𝐼=𝛽𝑆𝐼+𝜎𝛽𝐶𝐼−𝜇+𝜃𝐼,_𝛼𝑡()𝐷𝑅=1−𝜎𝛽𝐶𝐼+𝜃𝐼−𝜇+𝛿𝑅,_𝛼𝑡𝐶=𝛿𝑅−𝛽𝐶𝐼−𝜇+𝛾𝐶, (2.6) where𝐷^𝛼_𝑡 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 𝑡≥0 (see, for instance, [7, 10]).

Lemma 2.1. The closed set Ω={(𝑆,𝐼,𝑅,𝐶)∈𝑅⁴₊∶𝑆+𝐼+𝑅+𝐶=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) The solution to (2.7) is given by 𝑁(𝑡)=𝑁(0)𝐸_𝛼,1(−𝜇𝑡^𝛼)+𝜇𝑡^𝛼𝐸_𝛼,𝛼+1(−𝜇𝑡^𝛼), where 𝐸_𝛼,𝛽is the Mittag-Leffler function. Considering the fact that the Mittag-Leffler function has an asymptotic behavior [9, 11], 𝐸_∝,𝛽(𝑧)∼−^𝑝_𝑘=1𝑧^−𝑘Γ|𝑧|𝛽−𝛼𝑘+𝑂^−1−𝑝,|𝑧|⟶∞,𝛼𝜋2<arg(𝑧)≤𝜋. (2.8) One can observe that 𝑁(𝑡)converges to 1 when 𝑡→∞. Therefore, all solutions of the model with initial conditions in Ω remain in Ω for all 𝑡>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𝐷^𝛼_𝑡𝑆=0,𝐷^𝛼_𝑡𝐼=0,𝐷^𝛼_𝑡𝑅=0,𝐷^𝛼_𝑡𝐶=0. (3.1) Then 𝐸₀=(1,0,0,0). By (2.6), a positive equilibrium 𝐸₁=(𝑆₁,𝐼₁,𝑅₁,𝐶₁) satisfies 𝑆₁=𝜇+𝜃𝛽−𝜎𝛿𝜃𝐼₁𝜇+𝜎𝛿𝛽𝐼₁+,𝑅𝜇+𝛾𝜇+𝛿₁=𝜃𝐼₁𝛽𝐼₁+𝜇+𝛾𝜇+𝜎𝛿𝛽𝐼₁+,𝐶𝜇+𝛾𝜇+𝛿₁=𝛿𝜃𝐼₁𝜇+𝜎𝛿𝛽𝐼₁+,𝜇+𝛾𝜇+𝛿 (3.2) and𝐼₁ is the positive root of 𝑔(𝐼)=𝐴₁𝐼²+𝐴₂𝐼+𝐴₃(1−𝑅₀), where 𝑅₀=𝛽/(𝜇+𝜃) and 𝐴₁,𝐴=𝛽𝜇𝜃+𝜇+𝛿𝜎₂+𝜎,𝐴=𝛽𝜇−𝛽𝜇+𝛾𝛿+𝜃+𝜇𝜃+𝜇𝛿+2𝜇+𝛿−𝛽+𝜇₃=𝜇𝜇+𝛾𝜇+𝛿𝜇+𝜃. (3.3) The Jacobian matrix 𝐽(𝐸₀) for the system given in (2.6) evaluated at the disease-free equilibrium is as follows:𝐽𝐸₀=−𝜇−𝛽0𝛾0𝛽−𝜇−𝜃000𝜃−𝜇−𝛿000𝛿−𝜇−𝛾. (3.4)

Theorem 3.1. The disease-free equilibrium𝐸₀ is locally asymptotically stable if 𝑅₀<1 and is unstable if 𝑅₀>1.

Proof. The disease-free equilibrium is locally asymptotically stable if all the eigenvalues, 𝜆_𝑖,𝑖=1,2,3,4of the Jacobian matrix 𝐽(𝐸₀)satisfy the following condition [12, 13]: arg𝜆_𝑖>𝛼𝜋2. (3.5) The eigenvalues of the Jacobian matrix 𝐽(𝐸₀)are𝜆₁=−𝜇,𝜆₂=−(𝜇+𝛿),𝜆₃=−(𝜇+𝛾) and 𝜆₄=(𝜇+𝜃)(𝑅₀−1). Hence 𝐸₀is locally asymptotically stable if 𝑅₀<1and is unstable if 𝑅₀>1. We now discuss the asymptotic stability of the endemic (positive) equilibrium of the system given by (2.6). The Jacobian matrix 𝐽(𝐸₁)evaluated at the endemic equilibrium is given as: 𝐽𝐸₁=−𝜇−𝛽𝐼₁−𝛽𝑆₁0𝛾𝛽𝐼₁𝛽𝑆₁+𝜎𝛽𝐶₁−𝜇+𝜃0𝜎𝛽𝐼₁0(1−𝜎)𝛽𝐶₁+𝜃−𝜇+𝛿(1−𝜎)𝛽𝐼₁0−𝛽𝐶₁𝛿−𝛽𝐼₁−𝜇+𝛾. (3.6) The characteristic equation of 𝐽(𝐸₁)is 𝜆+𝜇𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃=0, (3.7) where 𝑎₁=𝛾+𝛿+2𝜇+2𝛽𝐼₁,𝑎₂=𝛾+𝜇𝛿+𝜇+𝛽𝐼₁𝛾+𝛿+𝜃+3𝜇+𝛿𝜎+𝛽𝐼₁,𝑎₃=𝛽𝐼₁Ψ𝛿+𝜇𝜃+𝜇+𝛾𝛿+𝜃+𝜇−𝛿𝜃𝜎+𝛽𝜃+𝜇+𝛿𝜎𝐼₁2Ψ+𝛽𝜇+𝛿𝜎𝐼₁Ψ+𝛽𝜇+𝛿𝜎𝐼₁, (3.8) where Ψ=(𝜇+𝛾)(𝜇+𝛾). Let 𝐷(Φ) denote the discriminant of a polynomial 𝑓. If Φ(𝑥)=𝑥³+𝑎₁𝑥²+𝑎₂𝑥+𝑎₃then. Denote 𝐷(Φ)=−1𝑎₁𝑎₂𝑎₃001𝑎₁𝑎₂𝑎₃32𝑎₁𝑎₂0003𝑎₁𝑎₂00032𝑎₁𝑎₂=18𝑎₁𝑎₂𝑎₃+𝑎₁𝑎₂²−4𝑎₃𝑎³₁−4𝑎³₂−27𝑎²₃. (3.9) Following [12, 14–16], we have the proposition.

Proposition 3.2. One assume that 𝐸₁ exists in 𝑅³₊. (i)If the discriminant of Φ(𝑥), 𝐷(Φ) is positive and Routh-Hurwitz are satisfied, that is, 𝐷(Φ)>0,𝑎₁>0, 𝑎₃>0, 𝑎₁𝑎₂>𝑎₃, then 𝐸₁ is locally asymptotically stable.(ii)If 𝐷(Φ)<0, 𝑎₁>0, 𝑎₂>0, 𝑎₁𝑎₂=𝑎₃, 𝛼∈[0,1) then 𝐸₁ is locally asymptotically stable.(iii)If 𝐷(Φ)<0, 𝑎₁<0, 𝑎₂<0,𝛼>2/3, then 𝐸₁ 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]: ^𝑘+1_𝑗=0𝜔^𝛼_𝑗𝑆_{𝑘+1−𝑗}=𝜇1−𝑆_𝑘+1−𝛽𝑆_𝑘+1𝐼_𝑘+𝛾𝐶_𝑘,_𝑘+1𝑗=0𝜔^𝛼_𝑗𝐼_{𝑘+1−𝑗}=𝛽𝑆_𝑘+1𝐼_𝑘+𝜎𝛽𝐶_𝑘+1𝐼_𝑘−𝜇+𝜃𝐼_𝑘+1,_𝑘+1𝑗=0𝜔^𝛼_𝑗𝑅_{𝑘+1−𝑗}=(1−𝜎)𝛽𝐶_𝑘+1𝐼_𝑘+𝜃𝐼_𝑘+1−𝜇𝑅_𝑘+1−𝛿𝑅_𝑘,_𝑘+1𝑗=0𝜔^𝛼_𝑗𝐶_{𝑘+1−𝑗}=𝛿𝑅_𝑘−𝛽𝐶_𝑘+1𝐼_𝑘−𝜇𝐶_𝑘+1−𝛾𝐶_𝑘. (4.1) Doing some algebraic manipulation to (4.1) yields the following relations:𝑆_𝑘+1=𝜇+𝛾𝐶_𝑘−∑^𝑘+1_𝑗=1𝜔^𝛼_𝑗𝑆_{𝑘+1−𝑗}𝜔^𝛼₀+𝜇+𝛽𝐼_𝑘,𝐼_𝑘+1=𝛽𝑆_𝑘+1𝐼_𝑘+𝜎𝛽𝐶_𝑘+1𝐼_𝑘−∑^𝑘+1_𝑗=1𝜔^𝛼_𝑗𝐼_{𝑘+1−𝑗}𝜔^𝛼₀,𝑅+𝜇+𝜃_𝑘+1=(1−𝜎)𝛽𝐶_𝑘+1𝐼_𝑘+𝜃𝐼_𝑘+1−𝛿𝑅_𝑘−∑^𝑘+1_𝑗=1𝜔^𝛼_𝑗𝑅_{𝑘+1−𝑗}𝜔^𝛼₀,𝐶+𝜇_𝑘+1=𝛿𝑅_𝑘−𝛾𝐶_𝑘−∑^𝑘+1_𝑗=1𝜔^𝛼_𝑗𝐶_{𝑘+1−𝑗}𝜔^𝛼₀+𝜇+𝛽𝐼_𝑘. (4.2) The approximate solutions 𝑆(𝑡), 𝐼(𝑡), 𝑅(𝑡), and 𝐶(𝑡) 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 𝑆(𝑡) versus time 𝑡 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 𝐼(𝑡) versus time 𝑡, for Figures 3, 4, and 5 showing the similar variations of 𝐼(𝑡), 𝑅(𝑡), and 𝐶(𝑡) with various values of α.

Figure 1 

𝑆(𝑡), 𝐼(𝑡), 𝑅(𝑡), 𝐶(𝑡) for 𝜇=0.02, 𝛽=100, 𝛿=1, 𝛾=0.5, 𝜎=0.05, 𝜃=73, 𝛼=0.9.

Figure 2 

𝑆(𝑡) for 𝜇=0.02, 𝛽=100, 𝛿=1, 𝛾=0.5, 𝜎=0.05, 𝜃=73.

Figure 3 

𝐼(𝑡) for 𝜇=0.02, 𝛽=100, 𝛿=1, 𝛾=0.5, 𝜎=0.05,𝜃=73.

Figure 4 

𝑅(𝑡) for 𝜇=0.02, 𝛽=100, 𝛿=1, 𝛾=0.5, 𝜎=0.05, 𝜃=73.

Figure 5 

𝐶(𝑡) for 𝜇=0.02, 𝛽=100, 𝛿=1, 𝛾=0.5, 𝜎=0.05, 𝜃=73.

Acknowledgments

The authors gratefully acknowledge the support from the Deanship of Scientific Research (DSR) at King Abdulaziz University (KAU) for the group entitled “Nonlinear Analysis and Applied Mathematics”. They also acknowledge the referees for their useful comments.