Bifurcations and chaos in a discrete SI epidemic model with fractional order

In this study, a new discrete SI epidemic model is proposed and established from SI fractional-order epidemic model. The existence conditions, the stability of the equilibrium points and the occurrence of bifurcation are analyzed. By using the center manifold theorem and bifurcation theory, it is shown that the model undergoes flip and Neimark–Sacker bifurcation. The effects of step size and fractional-order parameters on the dynamics of the model are studied. The bifurcation analysis is also conducted and our numerical results are in agreement with theoretical results.


Introduction
Mathematical modeling plays an important role in understanding the dynamics of many infectious diseases. Thus, the use of modeling is crucial in order to analyze the spread, control strategies and the mechanisms of transmission of diseases. Over the years, numerous epidemiological models have been formulated mathematically (see, e.g., [1][2][3][4][5][6]). Although most of these models have been restricted to integer-order differential equations (IDEs), in the last three decades, it has turned out that many problems in different fields such as sciences, engineering, finance, economics and in particular epidemiology can be described successfully by the fractional-order differential equations (FDEs) (see, e.g., [7][8][9][10][11][12]). A property of these fractional-order models is their nonlocal property which does not exist in IDEs. Nonlocal property means that the next state of a model depends not only upon its current state but also upon all of its historical states [13]. The transformation of an integer-order model into a fractional-order model needs to be precise as to the order of differentiation α: a small change in α may cause a big change in the final results [14]. FDEs can be used to model certain phenomena which cannot be adequately modeled by IDEs [15]. FDEs are commonly used on biological systems since they are related in a natural way to systems with memory. Similar to a nonlinear differential system, a nonlinear fractional differential system may also have complex dynamics, such as chaos and bifurcation. Studying chaos in fractional-order dynamical systems is an attractive and interesting topic [16].
Discrete models constructed by the discretization of continuous models have been used to describe some epidemic models. Many significant and meaningful types of such research can be found in [17][18][19][20][21][22][23][24][25][26][27] and the references therein. The reason for the use of dis-crete models is that statistical data on epidemics are collected in discrete times and hence comparing data with output of a discrete model may be easier [28]. In general, the dynamical behavior of discrete fractional-order epidemic models may exhibit phenomena such as the period-doubling and chaotic behavior. Iwami et al. [29] have proposed a mathematical model in a continuous-time version to interpret a model of the spread of avian-human influenza epidemic. In [29], the dynamics of the spread into bird population and between bird and human populations have been thoroughly studied. The main objective of the current study is to extend the research done in [29] by employing fractional models with discrete-time systems. In this paper, a discrete-time SI is established by SI fractional-order epidemic model. This model contains two parameters in addition to those already existing in the original SI model proposed in [29]; time step parameter and fractional-order parameter. The occurrence of bifurcations as these parameters are varied is shown in Section 7. The stability of fixed points, the emergence of flip bifurcation and Neimark-Sacker bifurcation are also studied. Using the discretized-time SI model with fractional order is a new topic, thus this paper provides a new contribution to the literature. This paper is organized as follows. In Section 2, the fractional-time SI epidemic model from its continuous-time counterpart is constructed. In Section 3, the discretized-time SI model with fractional-order parameter is established. Section 4 discusses the existence and stability of fixed points of the discretized-time SI model with fractional order. In Sections 5 and 6, flip and Neimark-Sacker bifurcations are demonstrated, respectively. In Section 7, some numerical examples are shown. A brief discussion of our results is given in Section 8.

The FDE epidemic model
Let us consider the following IDE epidemic model: This model is constructed by Iwami et al. [29] to explain the spreads of avian influenza through the bird world and describes the interactions between them. The population in this model is divided into susceptible birds with size S and infected birds with size I. The new birds birth rate is expressed by the parameter . Susceptible birds die at the rate μ and infected birds die at the rate μ + r, where r is the additional death rate mediated by avian influenza. The parameter β is the bilinear incidence rate.
There are several definitions of fractional derivatives [30,31]. One of the most common definitions is the Caputo definition [32]. This definition is often used in real applications and shown in Definition 1.

Definition 1
The fractional integral of order β ∈ R + of the function f (t), t > 0, is defined by and the fractional derivative of order α ∈ (n -1, n) of f (t), t > 0, is defined by where f (n) represents the nth order derivative of f (t), n = [α] is the value of α rounded up to the nearest integer, I β is the βth order Riemann-Liouville integer operator and (·) is Euler's Gamma function. The operator D α is called the 'αth order Caputo differential operator' . Now, the fractional-order form of the SI epidemic model (2.1) can be formulated as follows: where D α t represents the Caputo fractional derivative, t > 0, and α is the fractional order satisfying α ∈ (0, 1].

Discretization process
There are many discretization methods that have been used to construct the discretetime model using continuous-time methods such as explicit and implicit Euler's method, Runge-Kutta method, predictor-corrector method and nonstandard finite difference methods [33][34][35][36][37]. Some of them are approximation for the derivative and some for the integral. In [38,39] a discretization process was introduced to discretize FDEs. This discretization method is an approximation for the right hand side of the differential equation has the formula D α x(t) = f (x(t)), t > 0, α ∈ (0, 1). This method is now applied to the fractional-order SI model (2.2). Assume that S(0) = S 0 and I(0) = I 0 are the initial conditions of system (2.2). So, the discretization of the system (2.2) is given by the following formulas: 2) and the solution of (3.2) reduces to [ -μS 0 -βS 0 I 0 ], Second, let t ∈ [h, 2h), which makes 1 ≤ t/h < 2. Thus, we obtain which have the following solutions: Repeating the discretization process n times yields [ -μS n -βS n I n ], βS n I n -(μ + r)I n . (3.7) Remark 3.1 It should be noticed that if α → 1 in (3.7), the Euler discretization of SI model is obtained.

Stability of fixed points
In this section, an approach as in [40] is employed. The stability of the system (3.7) is studied around its fixed points. It is clear that the model (3.7) has always a disease free equilibrium point E 0 = ( μ , 0) and an endemic equilibrium point E 1 = ( μ+r β , μ+r -μ β ). Furthermore, the system (3.7) has basic reproduction number 0 = β μ(μ+r) > 1, then E 1 can be reformulated as ). It can be seen that the free equilibrium point E 0 always exists while the positive endemic E 1 exists only when 0 > 1. So as to analyze the dynamical properties of (3.7), we compute the Jacobian matrix J of (3.7), and we evaluate at the fixed point E = (S * , I * ) In order to study the stability of the fixed points of the system (3.7), the two following lemmas are employed. (i) If |λ 1 | < 1 and |λ 2 | < 1, then the equilibrium point of M(x * , y * ) is locally asymptotically stable (sink). (ii) If |λ 1 | > 1 and |λ 2 | > 1, then the equilibrium point of M(x * , y * ) is unstable (source). (iii) If |λ 1 | < 1 and |λ 2 | > 1 (or |λ 1 | > 1 and |λ 2 | < 1), then the equilibrium point of Based on Lemmas 4.1 and 4.2, the following results can be achieved.

Theorem 4.3 If 0 < 1, then the free-equilibrium point E 0 has at least four different topological types for all its values of parameters
where 0 < α ≤ 1 and h, h α (1+α) > 0. Hence applying the stability conditions using Lemma 4.1 the results (i)-(iv) can be achieved.
if one of the following conditions holds: (iv) E 1 is non-hyperbolic if one of the following conditions holds: , , Proof The Jacobian matrix of E 1 can be written as The characteristic equation of J(E 1 ) has the form Then the characteristic equation J(E 1 ) has two eigenvalues, which are By applying Lemmas 4.1, 4.2 and the Jury conditions [45], the stability conditions (i)-(iv) can be achieved.
When the statement (iv.2) of Theorem 4.4 holds, then the two eigenvalues of J(E 1 ) are a pair of conjugate complex numbers and the modules of each of them equals 1. The statement (iv.2) can be reformulated as follows: If the parameter h varies in the neighborhood of h 2 and (α, h, β, , μ, r) ∈ 3 , the system (3.7) may undergo a Neimark-Sacker bifurcation of equilibrium E 1 .

Numerical examples
This section shows the bifurcation diagrams, phase portraits and maximum Lyapunov exponents for the model (3.7) to confirm the above theoretical analysis and to illustrate the complex dynamics of our model using numerical continuation. Bifurcation occurs when the stability of an equilibrium point changes [48]. In general, the dynamics of a discrete SI model with integer-order has been examined by Hu et al. [42]. As discussed earlier in Section 1, this paper focuses on varying the time step size parameter h and the fractionalorder parameter α in the model (3.7), which can be seen as extension of the corresponding results given in [29,42]. Based on the previous analysis, the parameters of the model (3.7) can be examined in the following two cases.    Figure 2   In Case 2, the basic reproduction number 0 = 1.5 > 1, thus, model (3.7) has one positive endemic equilibrium point E 1 (5, 1). Since = -0.11 < 0, according to Theorem 4.4 E 1 is asymptotically stable (sink) when α < 0.826 and unstable when α > 0.826. The Neimark-Sacker (N-S) bifurcation emerges from E 1 at α = 0.826 with = -0.778 < 0 and (α, h, β, , μ, r) = (0.826, 8.1, 0.1, 1.5, 0.2, 0.3) ∈ 3 . The occurrence of an N-S bifurcation is illustrated in Figure 4(a)-(c). Figure 4 shows that E 1 is stable for α < 0.826 and loses its stability through an N-S bifurcation at α = 0.826. Attracting invariant circle appears as fractional-order parameter increases in the range of α ∈ (0.826, 0.87). The phase portraits for various α-values corresponding to Figure 4 are plotted in Figure 5(a)-(i) to illustrate these observations. Furthermore, the quasi-periodic orbits (α = 0.913) and periodic-13 orbits (α = 0.921) are observed within the chaotic regions in Figure 5(f ) and (g), respectively. Attracting chaotic sets are also seen when α increases further and these observations are plotted in Figure 5(h)-(i). Maximal Lyapunov exponents corresponding to observations in Figure 4(a)-(c) are computed in Figure 6. It is observed that some LE values are positive and some are negative. So there exist stable fixed points or stable period windows in the chaotic regions (e.g. when α > 0.921).

Discussion and conclusion
A new discrete-time SI epidemic model has been discussed in this paper. Such a discretetime model is obtained from the discretization of the fractional-time SI model. The discretization process provides crucial terms such as h (time step parameter) and α (fractional-order parameter), which are then varied in order to describe the dynamical behaviors of the model. As h and α are varied, the model exhibits several complicated dynamical behaviors including the emergence of flip and N-S bifurcations, period-2, 4, 8, 12, 13, 16 orbits, quasi-periodic orbits, attracting invariant circle and chaotic sets. Analytically, necessary and sufficient conditions on the parameters for the occurrence of flip and N-S bifurcations are derived. Moreover, numerical continuation is carried out to illustrate the validity of the analytical results and it is observed that both numerical and analytical findings are in good agreement. In conclusion, the proposed fractional-time SI model can engender more complex dynamical behaviors than its continuous model counterpart.
Studying the dynamical behaviors of the full avian-human influenza SIR epidemic model described by a discrete-time with fractional-order version will be our future work.