Abstract
The dynamics of a discrete-time SIS epidemic model are reported in this paper. Three types of codimension one bifurcation, namely, transcritical, flip, Neimark–Sacker (N-S) bifurcations, and their intersection codimension two bifurcations including 1 : 2, 1 : 3, and 1 : 4 resonances are discussed. The necessary and sufficient conditions for detecting these types of bifurcation are derived using algebraic criterion methods. Numerical simulations are conducted not only to illustrate analytical results but also to exhibit complex behaviors which include period-doubling bifurcation in period orbits, invariant closed cycles, and attracting chaotic sets. Especially, here we investigate the parameter space of the discrete model. We also investigate the organization of typical periodic structures embedded in a quasiperiodic region. We identify period-adding, Farey sequence of periodic structures embedded in this quasiperiodic region.
1. Introduction
Throughout the history of the world, epidemic diseases have posed a lasting threat to humans and society. For instance, tuberculosis, Ebola virus, Zika virus, influenza, etc. may bring many disasters to public security and individual health. Therefore, it has become an urgent and realistic problem to study the pathogenesis, transmission mechanism, development trend, and control strategy of infectious diseases. It is of great practical significance to analyze the spread of diseases and propose appropriate control strategies by using mathematical models and their dynamical properties [1–4]. The insights gained from mathematical modeling in epidemiology are important in the fight against epidemics.
In general, the dynamics of epidemics are a highly transited topic of investigation from the epidemiological point of view. There are various families of models depending on the assumption of the mechanisms of propagation. For example, there are susceptible-infected-susceptible (SIS) systems that have been used to model nonlinear incidence rates and double epidemic hypotheses [5], susceptible-infected-recovered (SIR) models which are employed to consider the awareness of the presence of a disease [6], susceptible-exposed-infected-recovered (SEIR) systems that predict the propagation of epidemics with nonlocal reaction functions [7], susceptible-exposed-infected-quarantined-recovered (SEIQR) systems used to account for adjusted incidence and imperfect vaccinations [8], and so on. Generally, in the theory of epidemic dynamical models, there are two kinds of mathematical models: the continuous-time models described by differential equations, and the discrete-time models described by difference equations. In recent years, more and more attention has been paid to the discrete-time epidemic models [9–12].
As is well known, difference equations have many applications in epidemic models [13–15]. In real life, data collection of disease is often based on a certain interval, and discrete models are more convenient for data simulation and comparison. On the other hand, a difference equation is the discretization of the continuous model, which makes it easier to handle for the numerical approach. Especially, discrete models exhibit more complex dynamical behaviors than the corresponding continuous models. With these advantages in mind, the problem of discrete epidemic models is of great research interest [16–18]. Based on the abovementioned reasons, we decide to consider a discrete-time epidemic model, which will be meaningful work. In this paper, we will study a discrete-time SIS epidemic model given by Hu et al. [19].where and denote the numbers of susceptible and infective individuals at time , respectively; is the recruitment rate of the population; is the natural death rate of the population; is the death rate of infective individuals which includes the natural death rate and the disease-related death rate; stands for the recovery rate of the infective individuals; and is the standard incidence rate. It is assumed that initial values and all the parameters are positive constants. And, is the step size.
For the sake of simplicity, we rewrite system (1) in the following form:
The objective is to study systematically the dynamics of model (2). Some of the codimension one bifurcations of the above model are investigated in [19]. Here, all of the bifurcations including codimension one and codimension two along with the normal form coefficients are provided; for more details, see [20–22].
The remainders of the paper are organized as follows: In Section 2, we present the existence and local stability of fixed points for map (2). In Section 3, we show that there exist some values of parameters such that map (2) possesses codimension one bifurcation. In Section 4, we extract the codimension two bifurcations of the model. In Section 5, we provide some numerical simulations, which not only validate the effectiveness of the theoretical analysis but also illustrate the mode-locking structures in the quasiperiodic rotation motions. Finally, we conclude this paper with comments and discussions.
2. Existence and Stability of Fixed Points
In this section, we discuss the existence of the positive fixed points of model (2). The basic reproductive rate of model (2) is , which is the average number of secondary infections generated by an initial population of infected individuals over their lifetimes. If , model (2) has a unique fixed point . If , model (2) has two fixed points: and , where
Let us consider model (2), given by the mapwhere . Its Jacobian matrix , where is given by
And, the multilinear functions of (2) used in its Taylor expansion can be expressed as follows:where
Thus,where
3. Codimension One Bifurcation Analysis
The following theorems provide the codimension one bifurcation analysis for fixed points and .
3.1. Codimension One Bifurcation at
Theorem 1. (1)There is a transcritical bifurcation of at when (2)There is a transcritical bifurcation of at when (3)There is a transcritical bifurcation of at when
Proof. We only provide the proof for the first part, and the proof for the other parts is similar. For map (2), there is a fixed point with a multiplier , ifwhereThe center manifold at can be considered asThe restriction of map (2) to the center manifold has the following form:where , andBy invariance property of the center manifold, we conclude thatNote that is always the fixed point and it will not be destroyed and . It is concluded that the fixed point undergoes a transcritical bifurcation.
Theorem 2. (1)There is a nondegenerate flip bifurcation of at , when (2)There is a nondegenerate flip bifurcation of at , when (3)There is a nondegenerate flip bifurcation of at , when (4)There is a nondegenerate flip bifurcation of at , when
Proof. We only provide the proof for the first part, and the proof for the other parts is similar. A fixed point of model (2) with the multiplier satisfies the following algebraic system:The center manifold at can be considered asThe restriction of map (2) to the center manifold has the formwhere , andAs regards the invariance property of the center manifold, the critical normal form coefficients of the flip bifurcation can be presented as follows:The flip bifurcation is nondegenerate provided . If , the flip bifurcation is supercritical and a stable double-period cycle bifurcates from . If , the flip bifurcation is subcritical and the double-period cycle is unstable.
3.2. Codimension One Bifurcation at
Theorem 3. (1)There is a nondegenerate flip bifurcation of at , where (2)There is a nondegenerate flip bifurcation of at , where (3)There is a nondegenerate flip bifurcation of at , where (4)There is a nondegenerate flip bifurcation of at
Proof. We only provide the proof for the first part, and the proof for the other parts is similar.
It is clear that at the point , map (2) has a fixed point with a multiplier , or on the other hand, this point convinces the following equation:The corresponding characteristic equation at is given bywhereThe next proof is similar to the proof of Theorem 2.
To avoid the complexity of the computation, we setThen, we haveBy using the invariance property of the center manifold, we can obtain the critical normal form coefficient . Since , the bifurcation is supercritical and a stable double-period cycle bifurcates from .
Theorem 4. (1)There is a nondegenerate Neimark–Sacker bifurcation of at .(2)There is a nondegenerate Neimark–Sacker bifurcation of at , where .(3)There is a nondegenerate Neimark–Sacker bifurcation of at , where .(4)There is a nondegenerate Neimark–Sacker bifurcation of at .
Proof. Only the first part is proved, and the proof of the other part is similar.
Map (2) has a fixed point with a pair of complex conjugate multipliers on the unit circle ifLet us consideras the center manifold at . Map (2) can be restricted to the center manifold which at the critical value is as follows:where is a complex number. To avoid the complexity of the computation, we letThus, has a pair of a simple multiplier.This applies to the nonresonance conditions. Note thathave been used, whereAnd, is a complex number given byAnd, .
After long but basic computations, we can achieve ; therefore, . Since , the Neimark–Sacker bifurcation is supercritical and the closed curve is stable.
4. Codimension Two Bifurcation Analysis
In this section, we study the codimension two bifurcations of model (2). By using the Jacobian matrix and some strict mathematical calculations, we find that does not come into being any codimension two bifurcations.
For , the following theorems give the codimension two bifurcation analysis.
Theorem 5. (1)If and , there is a 1 : 2 strong resonance(2)If and , there is a 1 : 2 strong resonance(3)If and , there is a 1 : 2 strong resonance
Proof. Only the first part is proved, and the proof of the other parts is similar.
Map (2) experiences a 1 : 2 strong resonance if it has a fixed point with two multipliers . Then, we haveThe exact solution is as follows:Let us consideras the center manifold at . Map (2) can be restricted to the center manifold at the critical value parameters of as follows:in whichWe havewhere satisfyTo avoid the complexity of the computation, we let . After long but basic computations, we obtain the explicit solutions and , which are the nondegenerate conditions of 1 : 2 resonance. The sign of determines the criticality of this bifurcation. The bifurcation scenario is indicated by the coefficient .
Theorem 6. (1)If and , there is a 1 : 3 strong resonance(2)If and , there is a 1 : 3 strong resonance(3)If and , there is a 1 : 3 strong resonance
Proof. Only the first part is proved, and the proof of the other parts is similar.
Map (2) undergoes a 1 : 3 strong resonance if it has a fixed point with a pair of conjugate multipliers . Then,The exact solution is as follows:The restriction of model (2) to two-dimensional center manifolds is expressed as follows:LetIt is easy to check thatModel (2) can be transformed into the following normal form:where can be determined byWe haveTo avoid the complexity of the computation, we let . Then, we have , and the stability of the bifurcating closed invariant curve is determined by .
Theorem 7. (1)If and , there is a 1 : 4 strong resonance(2)If and , there is a 1 : 4 strong resonance(3)If and , there is a 1 : 4 strong resonance
Proof. Only the first part is proved, and the proof of the other parts is similar.
Map (2) experiences a 1 : 4 strong resonance if it has a fixed point with a pair of conjugate multipliers . Then,The exact solution is as follows:The center manifold at can be considered as follows:LetIt is easy to check thatModel (2) can be transformed into the following normal form:where can be achieved by the following equation:The bifurcation scenario near the 1 : 4 resonance point is determined byTo avoid the complexity of the computation, we let . Then, we have .
5. Numerical Simulations
In this section, we will give bifurcation diagrams, phase portraits, and Lyapunov exponents to illustrate the above theoretical analyses and find new interesting complex dynamical behaviors by using numerical simulations. The bifurcation parameters are considered in the following two cases.
Case 1. Vary in range , and fix . By calculation, we find that at a unique endemic equilibrium , a flip bifurcation occurs at .
The flip bifurcation mentioned in Theorem 3 is depicted in Figure 1(a), from which we see that stability of fixed point happens for and loses its stability at and period-doubling phenomena lead to chaos for . The maximum Lyapunov exponent related to Figure 1(a) is computed and shown in Figure 1(b). We observe that the period orbits occur for and chaotic sets for . The period 6 orbit occurs at and period 5 orbit occurs at within the window of the chaotic region. Here, the powerful Lyapunov exponent test is shown in Figure 1(b) to confirm the existence of periodic orbits and chaotic motions.
(a)
(b)
Case 2. Fix and choose . By calculation, we observe that a Neimark–Sacker (NS) bifurcation appears at a unique epidemic equilibrium for . Figure 2 shows the correctness of Theorem 4.
The bifurcation diagrams shown in Figures 2(a) and 2(b) demonstrate that stability of happens for and there is a loss in its stability at and an attracting invariant curve appears if . We dispose the maximum Lyapunov exponent in Figure 2(c) relating bifurcation diagrams, which confirms the existences of chaos and period windows as parameter varies. When , the sign of the maximum Lyapunov exponent is greater than zero, confirming the presence of chaos. The status of stable, periodic, or chaotic dynamics is compatible with the sign in Figure 2(a). As we show, all our analytical predications have excellent arguments with numerical results.
The dynamical complexity of model (2) can be observed when more parameters vary in Figure 3. The bifurcation diagrams of system (2) for controlling parameters and and fixing remaining parameters as in Case 2 are shown in Figures 3(a) and 3(b). It is easy to find values of control parameters for which the dynamics of system (2) are in the status of quasiperiodic, periodic, or chaotic. This shows that the parameter may play an important role in chaotic dynamics in this system.
By the mathematical results in Theorem 6, we know that there exists 1 : 3 resonance when with , whose parameters are kept fixed at . To prove it, we present a global view of the parameter space in Figure 4(a), for and . Figure 4(a) is a diagram obtained by discretizing the parameter interval in a grid of points equally spaced [23, 24]. For each point , an orbit of initial condition converges to a chaotic, or a quasiperiodic, or a periodic attractor, or an attractor at infinity (diverges). In Figure 4(a), the period solution is shown by orange, the period is plotted in light gray, periods are codified in other different colors, and period motions are indicated by numbers. Moreover, quasiperiodic regions are indicated by the letter Q. Divergence regions are indicated by the letter D. Chaotic regions are indicated by the letter C. From Figure 4(a), we can also find 1 : 3 resonance when choosing a dark cyan region near the neighbourhood of . This is entirely consistent with the results we obtain analytically.
Figure 4(b) shows the magnification of a region with quasiperiodic behaviors in Figure 4(a) with and , which are periodic regions similar to Arnold tongues [25, 26]. For points on the white-orange borderline in Figure 4(b), on which Neimark–Sacker bifurcation (NSb) occurs, period orbits disappear, to give rise to invariant closed curves, which are the characteristics of a quasiperiodic motion. Therefore, quasiperiodicity is related to parameters near the white-orange border, on the white side, while chaotic solutions are corresponding to parameters far away from this same border, also on the white side. As is known [25, 26], the periodic structures embedded in the quasiperiodic regions possess self-similarity fractal features. In addition, these periodic structures are numbered as 7, 10, 13, 16, 19, 22, and 25. As there exists 1 : 3 resonance in this region, the difference between the adjacent period is 3. Due to the period-doubling bifurcation, we can see Arnold tongue sequence 22, 14, 20, 26. This is the main message of this paper: period-adding cascades exist in profusion in the SIS model.
To our surprise, by checking the rotation number of these Arnold tongues in Figure 4(b), it is observed that the mode-locking order is organized according to the broken Farey tree, which is organized according to a well-defined rule in [26]. Between each generic tongue of rotation numbers and , there is a tongue with a period equal to and rotation number . For instance, the tongue with the rotation number 5/18 is located between the tongues with the rotation numbers 3/11 and 2/7. Also are shown in Figure 4(b) the tongues 25 between 7 and 18, 24 between 7 and 17, 27 between 17 and 10, and 23 between 10 and 13. The sequence formed by the corresponding rotation numbers, namely, 5/18, 7/25, 2/7, 7/24, 5/17, 8/27, 3/10, 7/23, and 4/13, is only a small and incomplete part of a Farey tree, the so-called broken Farey tree.
In the sequel, we will concentrate on Figure 4(b), to investigate details on the dynamics of model (2). Figure 2 shows the bifurcation diagram and the maximum Lyapunov exponent diagram when . By increasing from 7 in Figure 2(a), the system stays in a period state until , being characterized by the largest Lyapunov exponent less than zero. At , the largest Lyapunov exponent reaches the zero value and a quasiperiodic region is attained, remaining until . The quasiperiodic motion is, therefore, originated from a period orbit and is a consequence of an NSb that occurs at . As a result of this NSb, one closed invariant is created, whose typical example can be seen in Figure 5(a), which shows a phase portrait for . There are mode-locking states (Arnold tongues) embedded in the quasiperiodic regions, for which the largest Lyapunov exponent is less than zero. By increasing from 7.6148, we cross the quasiperiodic region and penetrate a period tongue at . Period-doubling bifurcation happens at , resulting in period orbits. Examples of these periodic orbits can be seen in Figures 5(b) and 5(c). By increasing more, finally, a chaotic state for which the largest Lyapunov exponent is greater than zero is created at . A typical chaotic attractor appears in Figure 5(d), for .
From Theorem 6, we know that there also exists 1 : 3 resonance when and , whose parameters are fixed as , and , respectively. To prove it, we present the global view of the and parameter planes in Figure 6, where the color representations are as the same as in Figure 4. Similar to Figure 4, there also exist periodic structures which are numbered as 7, 10, 13, 16, 19, and 22. As there exists 1 : 3 resonance in the dark cyan region, the difference between the adjacent period is 3. In addition, we can also find periodic motion, quasiperiodic motion, and chaotic motion appear alternately. And, the complex mode-locking structures can also be seen in Figure 6.
(a)
(b)
(c)
(d)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
6. Conclusions
In summary, we give a detailed analysis for the codimension one bifurcation and codimension two bifurcation of a discrete-time SIS epidemic model. Although it is a simple discrete epidemic model, we find that there are many complex dynamical behaviors. The theoretical analysis demonstrates that model (2) undergoes codimension one bifurcation (transcritical, filp, Neimark–Sacker bifurcation) and codimension two bifurcation (1 : 2, 1 : 3, 1 : 4 strong resonances). For each bifurcation, the normal forms for maps and their approximations to the flows of the corresponding differential equations are obtained to determine the bifurcation diagrams. A good correlation between the experimental and theoretical results is found. Furthermore, we numerically report a visual view of the mode-locking structures to the discrete model in two-dimensional parameter-space diagrams. And, we also show Arnold tongue structures embedded in quasiperiodic regions regardless of any combination of parameters, respectively, organized themselves in the period-adding bifurcation cascades. The results obtained in this paper are interesting and useful and may be applied to a variety of epidemic models.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest in this paper.
Acknowledgments
This work was supported by Tarim University President Fund (TDZKQN201823 and TDZKSS201904) partially.