MULTIPLE RECURRENT OUTBREAK CYCLES IN AN AUTONOMOUS EPIDEMIOLOGICAL MODEL DUE TO MULTIPLE LIMIT CYCLE BIFURCATION

Multiple recurrent outbreak cycles have been commonly observed in infectious diseases such as measles and chicken pox. This complex outbreak dynamics in epidemiologicals is rarely captured by deterministic models. In this paper, we investigate a simple 2-dimensional SI epidemiological model and propose that the coexistence of multiple attractors attributes to the complex outbreak patterns. We first determine the conditions on parameters for the existence of an isolated center, then properly perturb the model to generate Hopf bifurcation and obtain limit cycles around the center. We further analytically prove that the maximum number of the coexisting limit cycles is three, and provide a corresponding set of parameters for the existence of the three limit cycles. Simulation results demonstrate the case with the maximum coexisting attractors, which contains one stable disease free equilibrium and two stable endemic periodic solutions separated by one unstable periodic solution. Therefore, different disease outcomes can be predicted by a single nonlinear deterministic model based on different initial data.


Introduction
Periodic outbreaks of infectious diseases have been observed worldwide from the beginning of the twentieth century. In particular, viral infectious diseases such as measles, rubella, and mumps could exhibit both regular annual or multiannual cycles and irregular dynamical patterns before [21] and after [25] the introduction of vaccination programs. The unpredicted erratic periodicities also occur in the epidemics of bacterial infections like whooping cough and walking pneumonia [25].
Understanding the dynamics of the re-emergence is crucial for the disease intervention strategies. Many mechanisms of sustain oscillations have been proposed and studied by using mathematical models [2,6,14,18]. The basic framework of these models consists of susceptible, infected, and recovered population classes. Many proposed mechanisms such as environmental changes, school calendar, and international travels are modeled as various types of contact rates. Based on compartmental epidemiological models and dynamical system theory, the irregular oscillatory phenomena are explained as chaotic behaviors or stochastic-driven switching among multiple attractors [3,9,28].
Seasonality has been well recognized as one of the driven forces for the irregular periodicities [1,11,30]. The seasonal mechanisms include the external influences on the host contact rate. This new contact rate can be periodically varying around a baseline with a seasonal sinusoidal forcing term [1,5,8,10,12,13,23,24,26]. Another type of contact rate is modeled as a periodic piecewise function [27] to mimic the school calender. Among these works, the non-autonomous deterministic epidemiological models exhibit multiple attractors leading to chaos through the well-known period-doubling scenario [4,7,8,22,23,29]. An advanced theoretical work using Melnikov's method [12] proves the existence of chaotic motion and homoclinic bifurcations.
Autonomous epidemiological models with nonlinear incidence rates have been proved to generate complex dynamical behaviors [17] including multiple equilibria, periodic motions and homoclinic loops [19,20]. However, the multiple attractors, specially the coexisting multiple stable periodic solutions, are barely studied over the past twenty years. One of the reasons is that the multiple periodic attracts are generated from small-amplitude limit cycle bifurcations, whose techniques rely heavily on tedious symbolic computation. The main contribution of this paper provides analytical proofs and numerical simulations of the coexistance of multiple periodic attractors in a classical epidemiological model with a general nonlinear incidence rate.
We consider the standard SIR (Susceptible-Infectious-Recovered) model [20]  (1.1) Under the measle season, we assume the constant recruitment rate of new susceptibles is λ [9]. Due to the advanced public health care, the infected individuals do not suffer heavy mortality. The three population classes share the same per capita death rate d [9]. F (I, S) is the incidence rate with a single infected individual. The duration of infection is v −1 . All parameters take positive real values. Since the infection of measle viruses gives lifelong immunity [15], model (1.1) is reduced to the following 2-dimensional system Here, we adopt the assumption in [20] that the incidence rate increases with the growth of susceptible population. That is F (I, 0) = 0 and ∂F/∂S > 0, for all I. Further F (I, 0) takes the same form in [20] as F (I, S) = S q G(I) and q = 1 [17], while G(I) is modified as where m > 0 and p > 1. We assume a nonzero infection rate in the absence of infected individuals, that is, lim I→0 G(I) = β 0 . The incidence rate is monotonically increasing in terms of infected populations, since dG/dI > 0 for p > 1, and saturated at (β 0 + k/m) I S, since lim I→+∞ G(I) = β 0 + k/m. System (1.2) can then be rewritten as We further introduce and obtain the following dimensionless system, where the new positive parameters are defined as In the next section, we present general solution properties and stability of equilibria of system (1.4). In section 3, we give certain conditions under which system (1.4) does not exhibit bifurcation of limit cycles. Then, in section 4, we first derive the condition for which system (1.4) has an isolated center and then study the limit cycle bifurcation around the center. Further, the existence of three limit cycles is proved and verified by simulation. Finally, conclusion is drawn in section 5.

Basic properties of solutions of system (1.4)
We first consider the basic properties of the solutions of system (1.4). For convenience, define the parameter space as where R 4 + means that all the four parameters take positive real values. We present a summary of some early results obtained in [32][33][34] without proof. Define the following trapping region, where 0 < ε ≪ 1. Then, we have the following result [32]. System (1.4) has two equilibrium solutions: the disease-free equilibrium E 0 and the Endemic equilibrium E 1 , given by given in the form of , (2.5) in which Note that the E 0 is a boundary equilibrium (on the X-axis), while the E 1 is an interior (positive) equilibrium. To show the bifurcation diagram, we may choose one of the parameters as bifurcation parameter and treat others as control parameters. Moreover, it is easy to define the reproduction number as R 0 = B D . So according to R 0 < 1, R 0 = 1 and R 0 > 1 (or respectively, B < D, B = D and B > D),we have the typical three bifurcation diagrams shown in Figure 1, where the parameter A is chosen as the bifurcation parameter. It should be noted that in the bifurcation diagrams only the part in the first quadrant and below the red line are biologically meaningful.   When no Hopf bifurcation exists in system (1.4), the endemic equilibrium E 1− is asymptotically stable. But whether or not it is globally stable depends upon additional conditions. The conditions for no-existence of limit cycles will be discussed in the next section. There may exist limit cycles enclosing the equilibrium. When A ≤ BC, the disease-free equilibrium E 0 is a saddle, and so single (unstable) limit cycle can not exist since all trajectories converge to the trapping region Ω, which has only one stable equilibrium inside. Application of Poincaré-Bendixson theory implies the possible existence of two limit cycles enclosing the E 1− , which will be shown in section 4.

Conditions for non-existence of limit cycles
It has been shown in [32] that system (1.4) can exhibit limit cycles for all three cases, R 0 < 1, R 0 = 1 and R 0 > 1. Now we want to ask under what conditions on the parameters there do not exist limit cycles. First of all, it is clear that if a limit cycle exists, it must enclose the equilibrium E 1− since the E 1+ is a saddle and the E 0 is a boundary equilibrium. Further, the existence of biologically meaningful equilibria implies Y 1 ≥ 0, yielding 0 < X 1− < 1 D . Hence, when ∆ < 0, there is no real solution for the equilibrium E 1 and so no limit cycles can exist. A necessary condition for the above inequality to hold is B < D. When ∆ ≥ 0, but X 1− ≥ 1 D , i.e., A+B−D−BC+ √ ∆ ≤ 0 which is equivalent to D ≥ max{B, A+B−BC}, the equilibrium E 1− is not in the first quadrant and thus no physically meaningful solutions exist for limit cycles. Summarizing the above discussions gives the following result. Theorem 3.1. For system (1.4), if one of the following conditions holds: then the system either has no equilibrium solutions or has no physically meaningful equilibrium solutions, and thus no limit cycles exist.

Next, define
Then, the remaining case is under the conditions: There are two cases: R 0 ≥ 1 and R 0 < 1. When R 0 ≥ 1 (i.e., B ≥ D), it has been shown in [32] that there must exist limit cycles if H 2 > 0. Thus, for R 0 ≥ 1 and H 2 ≤ 0, there may exist certain conditions on the parameters such that no limit cycles exist. Note that for this case, there are only two equilibrium solutions: a saddle point on the boundary at the E 0 : ( 1 D , 0), and a positive equilibrium at the E 1− , which may be a stable focus or node, or a center. The E 1+ is outside the first quadrant when R 0 > 1 and coincides with E 0 when R 0 = 1 for which E 0 becomes a degenerate saddle. So by Lemma 2.1, there exists a trapping region in the first quadrant, attracting all trajectories. Hence, the unstable manifold of the saddle point E 0 must converge either to a stable equilibrium or to a stable limit cycle which encloses the E 1− . In order for the system to have no limit cycles to exist, one needs to find the conditions such that the unstable manifold of the saddle E 0 converges to the stable equilibrium E 1− . It should be noted that if the unstable manifold of the saddle E 0 converges to a stable limit cycle, there must exist two limit cycles enclosing E 1− with inner one unstable and outer one stable. A single (unstable) limit cycle can not exist if R 0 ≥ 1 and H 2 ≤ 0, but it is certainly possible when H 2 > 0 since the limit cycle becomes stable.
The case for ∆ ≥ 0 and R 0 < 1 (B < D) is more complicated, since now the equilibrium E 0 is a stable node, and both the two positive equilibria E 1+ (which is a saddle point) and E 1− (which may be a focus or node, or a center) are located in the first quadrant. In this case, the unstable manifold of the saddle E 1+ may converge to the stable node E 0 , or to the stable equilibrium E 1− if H 2 < 0, or to a stable limit cycle which encloses the equilibrium E 1− . It should be pointed that unlike the case R 0 ≥ 1 which has limit cycles if H 2 > 0, while for this case R 0 < 1, the system may have no limit cycles for H 2 > 0 since the unstable manifold of the saddle E 1+ can converge to the stable node E 0 and the stable manifold of the saddle E 1+ connects the unstable equilibrium E 1− , called a heteroclinic orbit.
Next, we present a more general condition applicable for both cases R 0 ≥ 1 and

Theorem 3.2. When ∆ ≥ 0 and
there does not exist bifurcation of limit cycles in system (1.4).
Proof. By using the Dulac's Criterion, we choose the following function and let the vector filed of system (1.4) be denoted as (P, Q) T . Then, we have Moreover, since all trajectories are attracted to the trapping region where Y < 1, the condition can be further improved as

Bifurcation of multiple limit cycles around a center
It has been shown in [32] that maximal two limit cycles can bifurcate from in the equilibrium E 1− near a Hopf critical point. In this section, we will show that there exist conditions for which system (1.4) is integrable with a positive center, and three limit cycles can bifurcate from the center. To achieve this, we first apply the method of normal form computation to find the center conditions. In general, suppose the amplitude equation of normal form of a general dynamical system associated with a Hopf bifurcation is given by where the focus values v i 's are expressed in terms of the system parameters. We have the following lemma which can be used to find bifurcation of multiple limit cycles form Hopf critical point. The proof can be found, for example, in the book [16].
Then, for any given µ near µ c , the dynamical system can have at most k small limit cycles bifurcating from the origin; and for some µ near µ c , the system can have k limit cycles around the origin.
Now we derive the condition under which system (1.4) is integrable, and have the following theorem.

system (1.4) is integrable, with the first integral given by
.
Proof. For convenience, we use parameters A and C to solve the two equations: F 1 (X) = 0 and Trace(J) = 0, where J is the Jacobian of system (1.4), to obtain where X = X 1− is the solution given in (2.5). It follows from A > 0 and C > 0 that −1+X −BX 2 > 0 and 1−DX = Y > 0. Thus, it requires that 0 < X < 1 D and Under these conditions, the frequency at the Hopf bifurcation point is given by into system (1.4) yields a system with its linear part in Jordan canonical form, and then applying the Maple program [31] into (4.7), we obtain the first three focus values, given by where and the lengthy expressions of G 2 and G 3 are omitted here for brevity. Here, X takes the values from the equilibrium E 1− . Next, eliminating X from the three equations D), and two resultants, R 12 = R 0 R 12a and R 13 = R 0 R 13a , where R 12a and R 13a are two different polynomials which do not have common roots, and R 0 is given by Therefore, the center conditions come from R 0 = 0, which gives five solutions: By using the solution X = X 1 (B, D), we can verify that the first solution B = D(1−D) gives X = X 1 = 1 D , yielding −1+X−BX 2 = 0; and that the second, third and fourth solutions of the above result in A < 0 and C < 0. When B = D 1+D , the focus values v 1 , v 2 and v 3 have a common factor, 2−X +D(1−X) 2 , under which A and C are reduced to the conditions given in (4.2): A = 1 1+D and C = 1+D. Now under the conditions given in (4.2), system (1.4) can be rewritten as  In the following, we will perturb the integrable system (4.9) to study bifurcation of limit cycles. System (4.9) has three equilibrium solutions given by (2.3), for which X 1± are now equal to ) . (4.10) For the equilibrium E 0 , the eigenvalues of the Jacobian of system (4.9) are −D and −D 1+D , indicating that the E 0 is a stable node. Evaluating the Jacobian at the E 1± yields the characteristic polynomial ξ 2 + det ± , where ) , meaning that the equilibrium E 1+ is a saddle; and ) , implying that E 1− is a center. When D = 1 4 , E 1± = (3, 1 4 ) is a nilpotent center. In order to make the perturbation to the integrable system (4.9) have the same structure of the system (making it biologically meaningful), we consider the following perturbed system: where p(X, Y ) = a 12 XY 2 + a 11 XY + a 10 X + a 01 Y + a 00 , in which a ij and b ij are constant parameters.
Remark 4.1. The coefficients a 00 and b 00 can be used to make the perturbed system have a same equilibrium at the E 1− . It should be pointed out that if we use general cubic-degree polynomial perturbations, we will obtain at least seven smallamplitude limit cycles around the center. However, this would destroy the structure of the original system, and the result may only have mathematical interests. For the perturbations given in (4.12), there are two cases.
This makes sense because: (1) in reality this situation can happen as a patient who may not be completely recovered (implying thatȲ 0 ̸ = 0); and (2) in this study we are interested in the dynamical behaviour of this system around the disease equilibrium E 1− and so whetherȲ 0 is exactly equal to zero is not important.
In the following, we will consider these two cases separately. First, without loss of generality, introducing another time scaling τ = (Y +D+1)τ 1 into system (4.12) yields the following system: (4.13)

Remark 4.2.
If one uses Melnikov function method to consider bifurcation of limit cycles in system (4.13), one has to multiply the integrating factor γ to system (4.12) and thus the resulting vector field becomes rational functions, which will cause complexity in computation. When we use normal form (or focus value) computation method, the vector field is still in polynomial form, which greatly simplifies the computation. The difference is that Melnikov function method can be used for global limit cycle bifurcation around closed orbits while the focus value computation method is only applicable for local analysis; but they are equivalent in dealing with limit cycles bifurcating from a Hopf critical point. Proof. First, we have noticed from computing the focus values and solving the polynomial equations that there are only two independent perturbation coefficients in system (4.13). Although the choice of these coefficients is not unique, we may, for example, set a 12 = a 11 = a 10 = a 01 = a 00 = 0, i.e., p(X, Y ) = 0. (4.14) In order to make the perturbed system (4.13) still have an elementary center at the equilibrium E 1− , we let and then the frequency ω c at the Hopf critical point is given by into system (4.13) we obtain where the coefficients A ij0 , A ij1 , B ij0 and B ij1 are functions of D, b 12 and b 11 . Now, applying the Maple program [31] into (4.18) for computing the normal forms of Hopf and generalized Hopf bifurcations we obtain the focus values: where v i0 = 0, i = 1, 2, . . . since the unperturbed system is integrable, and v 11 is given by (4.20) Note that the zero-order focus value v 01 can be obtained from the trace of the Jacobian of (4.13) evaluated at the equilibrium E 1− , given as , (4.21) which equals zero when the solution b 02 given in (4.15) is applied. Next, solving v 11 = 0 for b 11 we obtain (4.23) Note that all higher-order focus values are given in the form of v k1 = (· · · ) b 12 , k = 2, 3, · · · . Thus, b 12 = 0 yields b 11 = b 02 = b 01 = 0, leading to the integrable system (4.9). Also, note that the term in the square bracket of v 21 does not equal zero for D ∈ (0, 1 4 ). So v 21 ̸ = 0, implying that at most two limit cycles can be obtained by perturbing v 01 and v 11 . Further, by direct perturbations on b 11 for v 11 and on b 02 for v 01 , we can obtain two limit cycles. There are infinitely many solutions since we have a free parameter D involved in the focus values.
Since the conditions given in Lemma 4.1 are sufficient but not necessary, we would further ask if we can choose a suitable value of D such that v 21 and v 31 have opposite sign and may yield three limit cycles. Without loss of generality, let b 12 = 1. A direct computation indeed shows that v 21 < 0 and v 31 > 0, and | v31 v21 | reaches its maximal value at D = 0.2014403192 · · · . Thus, we take D = 0. Then, the truncated normal form equation, v 01 + v 11 r 2 + v 21 r 4 = 0, gives two positive solutions: r 1 ≈ 0.0994 and r 2 ≈ 0.3582 to approximate the amplitudes of the two limit cycles, with the outer one stable and inner one unstable. If we add v 31 = 0.00294239 to the above truncated normal form equation we would obtain three solutions: r 1 ≈ 0.0994, r 2 ≈ 0.4174, and r 3 ≈ 0.7188, giving one more unstable large limit cycle. In order to check if the bifurcation of the three limit cycles is robust, we further add v 41 = −0.00592834 to the normal form equation but obtain only two positive solutions: r 1 ≈ 0.0994 and r 2 ≈ 0.3911. Therefore, we can not make a definite conclusion on the existence of three limit cycles.
This finishes the proof for Theorem4.2. A numerical example for simulation is taken from the above discussion. We take where the focus E − is stable for ϵ = 0.1 and unstable for ϵ = − 0.1. ϵ = 0.1. For this case, two limit cycles are obtained, which enclose the E − with outer stable and inner unstable. Simulations are shown in Figure 3(a), where two gray trajectories converge to a same stable limit cycle (in black loop).
Since E + is stable focus, there must exist an unstable limit cycle between the E − and the stable limit cycle, restricted to an invariant center manifold. ϵ = −0.1. For this case, the system has no limit cycle. As shown in Figure 3(b), the trajectories inside the separatrix are repelling away from the unstable focus E − , and they eventually converge to the stable node E 0 ; while while the trajectories outside the separatrix converge directly to the stable node E 0 .
(B) b 00 ̸ = 0. Now we turn to the case b 00 ̸ = 0, for which we have the following theorem. Proof. Similar to the case b 00 = 0, we may assume p(X, Y ) = 0, and make the perturbed system (4.13) still have an elementary center at the equilibrium E 1− by choosing for which the frequency ω c at the critical point is the same as that given in (4.16). Then, similarly introducing the affine transformation (4.17) into system (4.13) yields a system in the form of (4.18). Then, applying the Maple program [31] for computing the normal forms of Hopf and generalized Hopf bifurcations to obtain the focus values: where v 11 and v 21 are given by where C i 's, M 1 and M 2 are given in Appendix.
It is easy to verify that C 4 ̸ = 0 and C 5 ̸ = 0 for D ∈ (0, 1 4 ). Solving M 1 = 0 (i.e. v 11 = 0) for b 12 we obtain ] . (4.26) Next, we solve M 2 = 0 (i.e. v 21 = 0) to get (4.27) Here, C 6 and C 7 are two polynomials in D, given in Appendix. It is seen that b 02 ̸ = 0. Otherwise, b 11 = b 12 = b 01 = b 00 = 0, violating the assumption b 00 ̸ = 0. It can be shown that for D ∈ (0, 1 4 ), C 2 has a root at D = 0.2416575247 · · · , C 6 has a root at D = 0.2399255707 · · · , and C 7 has a root at D = 0.1546928262 · · · ≡ D * . Now, suppose C 2 C 6 C 7 ̸ = 0. Then using the above solutions of b 11 and b 12 , we obtain v 31 as follows: It can be shown that for D ∈ (0, 1 4 ), v 31 = 0 has a unique solution at D = D * , which cancels the root D * of C 7 . So v 31 ̸ = 0 for D ∈ (0, 1 4 ), implying that when C 2 C 6 C 7 ̸ = 0, bifurcation of four limit cycles are not possible from the Hopf critical point. Thus, at most three limit cycles can bifurcate from the E 1− due to Hopf bifurcation. Moreover, it can be shown that at the critical values of b 12 Therefore, by Lemma 3, maximal three limit cycles can bifurcate from the center E 1− for D ∈ (0, 1 4 ). The remaining task is to check the cases when C 2 C 6 C 7 = 0. When C 7 = 0, D = D * , which yields the following solutions for b 11 and b 12 : = 0.0026392443 · · · ̸ = 0, indicating that when C 7 = 0, the system can also have three small-amplitude limit cycles bifurcating from the E 1− . Similarly, we can show that there exist three small-amplitude limit cycles when C 2 = 0 or C 6 = 0.
This finishes the proof of Theorem 4.3.
To end this section, we present a numerical example to show bifurcation of three limit cycles. We again choose b 02 = 1, but D = 0.1875. Next we need to take perturbations such that 0 < v 01 ≪ −v 11 ≪ v 21 ≪ −v 31  to approximate the amplitudes of three limit cycles. Since v 31 < 0, the larger and smaller limit cycles are stable, and the middle one is unstable, all of them enclose an unstable focus E − . The perturbation parameter takes ε = 0.00001, for which ω ≈ 0.7016. For these parameter values, the following three positive equilibrium points are obtained: The simulation is shown in Figure 4. The overall picture looks similar to Figure 3, but now there exist three limit cycles enclosed by the separatries of the saddle E + . The large and small stable limit cycles are shown in Figure 4(a) and (b), respectively. In Figure 4(a), two trajectories starting from two different initial points (plotted in light and dark gray colors) converge towards the large stable limit cycle. In Figure 4(b), two trajectories starting from two different initial points (plotted in black and gray colors) converge to the small stable limit cycle. An unstable limit cycle exists between the large and small stable limit cycles.

Conclusion
In this paper, we have rigorously proved and numerically simulated the multiple attractors in an autonomous epidemiological model. Previous work [32] has shown  starting from two different initial points (plotted in light and dark gray colors) converging towards the large stable limit cycle; and (b) two trajectories starting from two different initial points (plotted in black and gray colors) converging to the small stable limit cycle; an unstable limit cycle between the larger and smaller stable limit cycles is not shown.
that under proper perturbations on the model parameters, the model can exhibit at maximal two small-amplitude stable limit cycles around an unstable endemic equilibrium due to Hopf bifurcation. In this paper, we first find the condition on the parameters such that the system becomes integrable with the Hopf critical point becoming an isolated center. Then, we perturb the integrable system to prove the existence of at least three limit cycles, one more than that obtained in the previous work. With properly chosen parameter values and perturbations, we simulate the cases with two and three co-existing limit cycles. This indicates that an autonomous epidemiological model can indeed yield complex dynamical behaviors due to Hopf bifurcation, yielding multiple attraction regions in the phase plane.
The significant contribution of this work proves the existence of multiple periodic attractors in an autonomous epidemiological model with a nonlinear incidence rate. The incidence rate shows a cooperative effect on the contact rate in terms of infected populations. A large number of scientific and mathematical papers use non-autonomous models to generate multiple periodic attractors. In their works, incorporated with the seasonality, the contact rates are mainly considered as a period function in terms of time. Therefore, in these works, the multiple periodic attractors are generated by the periodically forcing functions, rather than an intrinsic effect from the model itself.
In a fixed set of parameter values, multiple attractors can be obtained in the autonomous epidemiological model for the case b 00 = 0 and ϵ = 0.1, as depicted in Figure 3(a), showing three attractors: the stable disease-free equilibrium E 0 , the stable endemic equilibrium E − , and the interior stable limit cycle. It is seen from Figure 3(a) that the separatrix from the saddle E + and the unstable limit cycle delimit the first quadrant of the phase plane into three regions, and the three attractors serve as the attraction basins for the three regions. The fate of the given trajectory is determined by the located region of its initial condition in the phase plane. The multiple attractors shown in Figure 4 with b 00 ̸ = 0 include the stable disease-free equilibrium E 0 , the large and small stable limit cycles. The separatrix from the saddle E + and the sandwiched unstable limit cycle delimit the phase plane into three attraction regions. The long-term behavior of a trajectory depends upon its initial condition. This indicates that for a given situation, where the corresponding set of parameter values are all fixed, the emerging diseases and regular and irregular re-emerging diseases can be explained by the bifurcation of multiple attractors and stochastic influences on the initial conditions. This work reveals that multiple attractors may be an intrinsic dynamical property causing disease outbreak, and the uncertainty may be due to stochastic influences.