Dynamical Behaviors of an SI R Epidemic Model with Discrete Time

: Analytically and numerically, the study examines the stability and local bifurcations of a discrete-time SIR epidemic model. For this model, a number of bifurcations are studied, including the transcritical, ﬂip bifurcations, Neimark–Sacker bifurcations, and strong resonances. These bifurcations are checked, and their non-degeneracy conditions are determined by using the normal form technique (computing of critical normal form coefﬁcients). We use the MATLAB toolbox M ATCONT M, which is based on the numerical continuation method, to conﬁrm the obtained analytical results and specify more complex behaviors of the model. Numerical simulation is employed to present a closed invariant curve emerging from a Neimark–Sacker point and its breaking down to several closed invariant curves and eventually giving rise to a chaotic strange attractor by increasing the bifurcation parameter.


Introduction
There is a great deal that can be done to minimize the impact of infectious diseases through research. With relevant knowledge about the dynamics of an infection, disease transmission can often be prevented. The transmission and dynamics of most infectious diseases are greatly influenced by seasonal factors, including climatic factors and human phenomena [1]. It appears that intense seasonality causes erratic patterns based on some empirical data. The presence of chaotic oscillations in response to seasonal forces has been demonstrated in many studies focusing on seasonal influenza and measles [2,3]. When vaccination programs are not in place, many recurrent infectious diseases exhibit strong annual, biennial, or irregular oscillations in response to seasonality [4].
A mathematical model called the SIR allows us to estimate the number of people infected with a disease in a closed population over time. Susceptibility, infection, and recovery models are included in the SIR group.
The behavior of epidemic diseases is being studied in an effort to detect and control them. One of them is the dynamical epidemic model, which is used to study epidemics [5][6][7][8][9][10][11][12][13][14]. The dynamical nature of the measles epidemic model is analyzed in [14], and the dynamical nature of the disease is also strongly influenced by migration processes. The study in [15] demonstrated that chickenpox prevalence is inversely related to the size of the population on an annual cycle.
The mathematical modeling of infectious diseases leads to detect the dynamical behavior of epidemics and provide sufficient disease control measures. The treatment of epidemics is studied using a dynamic model [5][6][7][8][9]. An SIR (Sensitive-Infected-Improved) model has been used to study the dynamics of the measles epidemic [14]. According to research, the disease is very sensitive to migration, and its onset was accompanied by an epidemic. Based on [15], chickenpox prevalence is inversely proportional to population size on an annual cycle. In this paper, we aim to analyze model (3) for different types of bifurcations. The novelty of our work is that we studied the bifurcation results of discretetime SIR epidemic model [1,4,16,17], as none of the studies in the literature has studied the complex dynamic behavior of the model.
According to classical infectious disease transmission models (Kermack and McKendrick [18], Hethcote [19]), the population is divided into three classes derived from S(ι), I(ι), and R(ι) indicating the number of susceptibles, infected individuals, and recovered or removed individuals at time ι, respectively. Here, we develop a model of SIR epidemics based on a modified saturated incidence rate as follows: where the saturated contact rate is indicated by 1+α S(ι) and more information about parameters can be found in Table 1. Model (1) focuses on the following model because the first two equations are not dependent on R; see [20,21]: Euler's method is applied to (2) to obtain the following discrete time model: where h is step size.

Stability of Fixed Points
First, we present a lemma that describes the dynamics of the model (1).

Lemma 1.
In the first quadrant, the plane S + I + R = Λ δ is an invariant manifold of the system of the model (1).
Proof. When we combine the three equations in (1) and denote N(τ) = S(τ) + I(τ) + R(τ), we obtain For any N(τ 0 ) ≥ 0, is the general solution. So, we have the result can be drawn from this.
There are two fixed points for model (3) as follows: can be selected as the Jacobian matrix of (3) at E 0 and E * Theorem 1. If R 0 < 1, the fixed point E 0 is asymptotically stable when − δ (δ+σ) Proof. See [22,23].

Bifurcation Analysis of the Boundary FIXED Point E 0
Our goal in this part is to investigate the bifurcations of model (3) at the trivial fixed point E 0 by computing corresponding critical normal form coefficients; see [24][25][26].
In this section, the parameter Λ is considered as a bifurcation parameter. Proof. For Λ = Λ LP , A 0 has the following multipliers: When h = 2 δ , the Jacobian matrix A 0 has a single multiplier +1 and no other multiplier with |λ| = 1. So, the model (3) at β = β LP,0 can be reduced to its normal form E 0 develops a transcritical bifurcation because it is always the fixed point and will never disappear, and ϕ LP,0 = 0. Proof. For Λ = Λ PD,0 , A 0 has the following eigenvalues: When λ PD,0 2 = ±1, the Jacobian matrix A 0 has a single multiplier −1 and no other multiplier with |λ| = 1. So, the model (3) at Λ = Λ PD,0 can be reduced to its normal form As long as φ PD,0 > 0 (φ PD,0 < 0), the flip bifurcation is super-critical (sub-critical, resp.); moreover, the two-period emerging cycle is stable (unstable, resp.).

Bifurcation Analysis of the Positive Fixed Point E *
Based on the equations given in [24][25][26], we will determine the critical normal form coefficient at the bifurcation points of the model (3).

One Parameter Bifurcations
Λ is referred to as a bifurcation parameter.

Theorem 5. The presence of
causes a flip bifurcation of E * .
Proof. For Λ = Λ PD, * , A * has the following multipliers: The Jacobian matrix A SIR * has a simple eigenvalue −1 and no other eigenvalue with |λ| = 1 if λ PD, * 2 = ±1. So, the model (3) at Λ = Λ PD, * can be reduced to its normal form An indication of the type of flip bifurcation is given by the sign of φ PD, * . The bifurcation is supercritical (sub-critical) if it is positive (negative).

Continuation Method
The numerical bifurcation analysis is performed using the MATLAB package MAT-CONTM; see [27].

Numerical Continuation of E 0
Taking into account the following fixed parameters which will lead to a numerical continuation of E 0 : We consider Λ as a bifurcation parameter.
By varying Λ, the flip bifurcation occurs at E 0 for Λ = 0.040541 where φ PD,0 = −5.448258 × 10 −2 . The sign of φ PD,0 determines the sub-critical flip bifurcation. The stability region of E 0 near E 0 is presented in Figure 1.

Numerical Continuation of E *
Taking into account the following fixed parameters which will lead to a numerical continuation of E * : We consider Λ as a bifurcation parameter.
By varying Λ, we can obtain the following one-parameter bifurcations: The following bifurcations can be obtained with two parameters, based on the selected the Neimark-Sacker point and the continuation with two free parameters (Λ, σ): 1.

2.
The resonance 1:3 bifurcation occurs at E * for Λ = 4.999711 and σ = 0.392641 where 3 4 (2e i 4π 3 υ R 3 , * − |γ R 3 , * | 2 ) = 4.403109 × 10 −2 . If we compute the convergent orbits from initial point (S, I) = (2.8495, 5.9841) with respect to Λ and σ, a two-dimensional bifurcation diagram in the neighborhood of the R3 point can be displayed with the period number of the corresponding orbits; see Figure 5. In addition to the parameter region with a period-3 cycle, there only exist regions with fixed points and a period-2 cycle. Here, a stable period-3 cycle occurs when (Λ, σ) = (4.9, 0.38) and one of the period-3 cycle is (2.796052631578960,5.972329721362207).

3.
The resonance 1:2 bifurcation occurs at E * for Λ = 6.321289 and σ = 0.357760 where υ R 2 , * = −1.731233 × 10 −1 and γ R 2 , * = 1.575539 × 10 −2 . If we compute the convergent orbits from initial point (S, I) = (2.7021, 8.3779) with respect to Λ and σ, a twodimensional bifurcation diagram in the neighborhood of the R2 point can be displayed with the period number of the corresponding orbits; see Figure 6. In addition to the parameter region with a period-2 cycle, there only exist regions with fixed points and period-4, -6, and -8 cycles. Here, a stable period-2 cycle occurs when (Λ, σ) = (5.6, 0.32) and one of the period-2 cycles is (2.232027014018290, 8.330641606985276); see Figure 7.      The stability region of E * in space (Λ, σ) and the bifurcation curves of the flip and the Neimark-Sacker are shown in Figures 8 and 9. Figure 9 confirms the results of Theorems 7-9.  The bifurcation curves of the second and third iterates of (3) are presented in Figure 10a,b.

Discussion
We presented a discrete-time SIR epidemic model with detailed complex dynamics in this study. Using analytical and numerical methods, we analyzed the bifurcation of the boundary and positive fixed points E SIR 0 and E SIR * . By incorporating the Neimark-Sacker bifurcation into the model, we can infer that susceptible and infective individuals can fluctuate around some mean values of the population recruitment rate, and these fluctuations remain stable as well as constant if υ NS, * < 0. According to biological theory, an invariant curve bifurcates from a fixed point, which allows susceptible and infected individuals to coexist and produce their densities. Periodic or quasi-periodic dynamics may be observed on an invariant curve. It appears that the susceptible and infected individuals change from one period to the next in this model based upon the period-doubling bifurcation. On the other hand, the strong resonances bifurcation of the model suggests susceptible and infected individuals coexist in stable high period cycles around some mean values of the rate of population recruitment rate and infection rate of infected individuals. Some two-dimensional bifurcation diagrams in the neighborhood of two-parameter bifurcation point are computed and displayed to show possible periodic dynamics.

Conclusions
The dynamics of a system can be identified and predicted using bifurcation theory. In this sense, bifurcation theory is an important branch of dynamical systems theory. In this paper, we provide a standard research format of bifurcation analysis. The existence and stability of fixed points are provided in Section 2. In Sections 3 and 4, one-parameter bifurcations and two-parameter bifurcations are analyzed, respectively. Detailed instructions are given in Section 5 regarding the computation of fixed point curves. According to Sections 3-5, the numerical observations and the analytical predictions are in excellent agreement. Discussions are summarized in Section 6. Both analytical and numerical aspects of bifurcations are considered in dynamic models. Some methods are more efficient than others for studying bifurcations in each of these two aspects. The computation of the critical normal form coefficients is a very effective analytical method in bifurcation theory. One can see different kinds of methods employed in bifurcation analysis [30][31][32][33][34][35][36]. There are many dynamical systems that are prone to notice this method, discrete or continuous; see [37][38][39][40][41][42][43]. An analytical computation is performed in this paper, and the results are also validated numerically using MATCONTM. More details can be found in Kuznetsov and Meijer (2005) [25] and Govaerts et al. (2007) [27]. The paper also provides a robust analytical and numerical method that can be applied to different discrete-time models. Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable
Acknowledgments: The authors express their sincere thanks to three referees for careful reading and valuable suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.