HOW SEASONAL FORCING INFLUENCES THE COMPLEXITY OF A PREDATOR-PREY SYSTEM

. Almost all population communities are strongly inﬂuenced by their seasonally varying living environments. We investigate the inﬂuence of seasons on populations via a periodically forced predator-prey system with a nonmonotonic functional response. We study four seasonality mechanisms via a continu- ation technique. When the natural death rate is periodically varied, we get six diﬀerent bifurcation diagrams corresponding to diﬀerent bifurcation cases of the unforced system. If the carrying capacity is periodic, two diﬀerent bifurcation diagrams are obtained. Here we cannot get a “universal diagram” like the one in the periodically forced system with monotonic Holling type II functional response; that is, the two elementary seasonality mechanisms have diﬀerent ef- fects on the population. When both the natural death rate and the carrying capacity are forced with two diﬀerent seasonality mechanisms, the phenomena that arise are to some extent diﬀerent. The bifurcation results also show that each seasonality mechanism can display complex dynamics such as multiple attractors including stable cycles of diﬀerent periods, quasi-periodic solutions, chaos, switching between these attractors and catastrophic transitions. In ad- dition, we give some orbits in phase space and corresponding Poincar´e sections to illustrate diﬀerent attractors.


1.
Introduction. With almost no exception, population communities are strongly influenced by their seasonally varying living environments. Seasonal temperature variation strongly influences the reproduction rate of microorganisms, the birth rate of animals, the migration of birds, and also the breakout of some diseases among communities. Seasonal variation of light density regulates photosynthesis, further affecting the survival of herbivores. As a result, it is important to study the effect of seasons on the behaviour of population communities. One basic way to represent this effect is in ordinary differential systems with periodically varying parameters. It is well known that even a simple autonomous system can become very complex when the system is periodically forced [7]. Over the last few decades there have been extensive studies about the effect of the seasonal perturbation on population communities [10,16,6,18,19,2,12,15,20].
To study the dynamical behaviour of many predator-prey communities, the following classical system has been used for many years, where x and y are the numbers of prey and predator populations or suitable (but equivalent) measures of density or biomass; r, K, d and c are positive parameters. In the system (1), the prey population is assumed to satisfy a logistic growth, where r and K are the intrinsic growth rate of the prey population and the carrying capacity of the resources respectively. The natural death rate of the predator is denoted by d, and p(x) is the predator response function. It is assumed that the rate of the conversion of prey captured to predator, is proportional to the predator response function where c is the constant of proportionality or yield constant. Various types of response functions have been used to study the dynamics of the predator-prey system, especially Holling type response functions. According to Holling [8,9], the functional responses were originally classified into the following three types.
(i) Holling type I: linear rise to a plateau.
where a is the half-saturation constant; that is, the number or density of prey at which the per capita predation rate is half of its maximum m. (ii) Holling type II: negatively accelerated rise to a plateau.
where m and a are positive constants and a is called the half-saturation constant. (iii) Holling type III: S-Shaped rise to a plateau.
where m and a are positive constants. When p(x) = mx 2 ax 2 + bx + 1 DYNAMICS OF A SEASONALLY FORCED SYSTEM 787 with b > −2 √ a, it is called the generalized Holling type III or sigmoidal functional response. For b < 0, the predation increases to a maximum and then decreases, approaching m/a as x approaches infinity. In microbial dynamics, there are experiments indicating that an inhibitory effect on the specific growth rate may occur when the nutrient concentration reaches a high level. In population dynamics, field observations demonstrate that predation is decreased, or even prevented altogether, due to the increased ability of the prey to better defend or disguise themselves when their numbers are large enough. To describe the phenomenon of "inhibition" in microbial dynamics and "group defence" in population dynamics, the functional response is proposed [1], where a and m are positive constants; b > −2 √ a. It is called the Monod-Haldane [17] or the generalized Holling type IV [3] functional response. When b = 0, the function is called the Holling type IV functional response in the literature [17,13]. For m, a, b > 0, and x sufficiently small, p(x) resembles the Holling type II functional response while for large x the effect of inhibition can be found. If −2 √ a < b < 0 and a > 0, p(x) remains nonnegative, the inhibition effect is still observed for large x, however, for small x, p(x) resembles the Holling type III (sigmoidal) function. See [5,21] for more information about "group defence", for examples of populations that use this strategy, and for the analysis of system (1) in the case of a general unimodal response function.
In studying seasonal effect on predator-prey populations, it is natural to consider a sinusoidal perturbation of the form q = q 0 (1 + ε sin 2πt) (3) for any parameter q in system (1). q 0 is the average value of q and ε is the "degree" of seasonality (note that εq 0 is the magnitude of the perturbation). In [10] and [16], the authors studied system (1) with Holling type II response function by considering the "elementary" seasonality mechanisms, namely phenomena that involve periodic variations of a single parameter or periodic but synchronous variations of two parameters. They identified six elementary seasonality mechanisms. Through comparison of the bifurcation diagrams for each of the elementary mechanisms, the authors obtained a "universal diagram", i.e., the six seasonality mechanisms gave rise to the same phenomena. In [6], an analogous predator-prey system in a chemostat was investigated with Holling type II functional response, and bifurcation diagrams equivalent to the universal one were obtained. With the inhibition response function of generalized Holling type IV (2), system The authors in [22] carried out a complete bifurcation analysis of the unforced system (4) which revealed the complexity of the dynamics due to group defence. Besides transcritical and Hopf bifurcations (both supercritical and subcritical), the unforced system (4) undergoes homoclinic bifurcation(s), and saddle-node bifurcation of limit cycles, as well as a Bogdanov-Takens bifurcation of codimension 3. Compared to the unforced system (1) with Holling type II response function, the unforced system (4) has more complicated dynamics. Therefore, in [15] the authors investigated the complex dynamics of system (4) with periodic d. Naturally new questions arise: whether different periodic parameters have the same effect on the predator-prey communities, or whether a "universal diagram" can exist when the system (4) is periodically forced with different elementary seasonality mechanisms (periodic variations of a single parameter), and what phenomena there will be if different parameters are varied periodically and non-synchronously. Motivated by these questions, in this paper we analyze system (4) by periodically forcing d or/and K, to determine whether a "universal diagram" exists or not and to see the complex dynamics when both d and K are periodically forced. This work is organized as follows. In section 2 we give the method of investigating periodically forced systems. In section 3 we use numerical continuation [4] to generate bifurcation diagrams of the Poincaré map for the seasonally forced system with periodic d or K. Bifurcation diagrams when d and K are periodically forced are presented in section 4. We analyze these diagrams further in section 5, show that the seasonality mechanisms can give rise to multiple attractors and that catastrophic transitions can occur if suitable parameters are varied, and conclude that the forced predatorprey system (4) does not have a universal diagram when d or K is periodically forced. Moreover, we give some typical orbits in phase space and corresponding Poincaré sections to illustrate the quasiperiodic and chaotic solutions by MATLAB numerical simulations [23]. We end the paper with a brief conclusion. 2. The model and method of investigation. In predator-prey communities, the periodic presence of a superpredator exploiting the predator community can be taken into account by periodically varying the natural death rate of predators d. The periodicity in the amount of food available to the prey population leads to periodic carrying capacity K. Therefore, in this paper we will analyse the seasonally forced predator-prey system (4) with four different seasonality mechanisms. That is, 1. d = d 0 (1 + ε sin 2πt); 2. K = K 0 (1 + ε sin 2πt); 3. d = d 0 (1 + ε 2 cos 2πt) and K = K 0 (1 + ε 1 sin 2πt); 4. d = d 0 (1 + ε 2 sin 4πt) and K = K 0 (1 + ε 1 sin 2πt).
In order to study the dynamics of a periodically forced system of ODEs, one can transform the system into a higher dimensional autonomous system by augmenting the forced system with a decoupled pair of ODEs which have an oscillatory solution of the form required for the forcing [4]. As an example, we take the case when K is periodically forced. The higher dimensional autonomous system iṡ with initial conditions v(0) = 0 and w(0) = 1. The system (5c)-(5d) has an asymptotically stable periodic solution v = sin 2πt, w = cos 2πt.
In this way the equilibrium (x * , y * ) of the two dimensional unforced system corresponds to the periodic solution (x * , y * , sin 2πt, cos 2πt) of the four dimensional system with ε = 0. We can study whether the periodic solution can survive and whether the bifurcations of the equilibrium of the unforced system can be extended to the forced system as bifurcations of periodic solutions when ε = 0. The standard procedure for identifying and locating bifurcations of a periodic system such as (5) uses the Poincaré map P that transforms the continuous system into a discrete one by sampling the solution once in each forcing period. For our system the forcing period is one year, so the Poincaré map is The phase portrait of the Poincaré map P is composed of fixed points, regular and irregular invariant sets and all other orbits. In particular, we will consider • the fixed points of the kth iteration of the map, which correspond to the subharmonic periodic solutions of period k; • the closed and regular invariant sets of the map, which correspond to the quasiperiodic solutions (invariant tori); • the irregular invariant sets, which correspond to chaotic solutions (strange attractors). Fixed points of the map P bifurcate at parameter values where the fixed points change stability. For a periodically forced system of ODEs, the corresponding Poincaré return map undergoes bifurcations including Hopf (Neimark-Sacker), flip (period doubling) and tangent (fold) bifurcations. Following the convention in [10], we use h (k) , f (k) and t (k) to denote the above bifurcation curves respectively.
1. Hopf (Neimark-Sacker) curve h (k) . Along the curve h (k) , P has a period k fixed point with a pair of complex conjugate multipliers on the unit circle µ (k) 1,2 = exp(±iω). When the curve is crossed, the attracting (unstable) cycle of period k bifurcates into an attracting (unstable) quasiperiodic solution and an unstable (attracting) cycle of period k.

Flip bifurcation curve f (k) .
Along the curve f (k) , P has a period k fixed point with a multiplier µ (k) = −1. When the curve is crossed, a saddle (non-saddle) cycle of period k bifurcates into a non-saddle (saddle) cycle of period k and a saddle (nonsaddle) cycle of period 2k. 3. Tangent bifurcation curve t (k) .
3. Bifurcations of the system with periodic d or K. Based on the theoretical analysis in [22], regardless of the generalized Hopf bifurcation case the unforced system can have the following six different Hopf bifurcation cases.
Case 1. The unforced system undergoes two supercritical Hopf bifurcations when d is varied; it undergoes one supercritical Hopf bifurcation if K is varied.
Case 2. The unforced system undergoes one supercritical and one subcritical Hopf bifurcation when d is varied; it undergoes one supercritical Hopf bifurcation if K is varied.
Case 3. The unforced system undergoes one supercritical and one subcritical Hopf bifurcation when d is varied; it undergoes one subcritical Hopf bifurcation if K is varied.
Case 4. The unforced system undergoes one subcritical Hopf bifurcation when d or K is varied.
Case 5. The unforced system undergoes one supercritical Hopf bifurcation at the local minimum of F (x) = r m (1 − x K )(ax 2 + bx + 1) when d is varied; it undergoes one supercritical Hopf bifurcation if K is varied.
Case 6. The unforced system undergoes one supercritical Hopf bifurcation at the local maximum of F (x) when d is varied; it undergoes one supercritical Hopf bifurcation if K is varied.
We will explore and present bifurcation diagrams of the Poincaré map for the periodically forced system corresponding to the above six cases of the unforced system in this section. The correspondence of solutions for the original continuous system to those for the Poincaré map has been given above. Taking a = 4, b = 0, c = 1, m = 5π, r = 4.8π, K 0 = 0.89, d 0 = 1.1π, the unforced system undergoes two supercritical Hopf bifurcation at d 0 = d + and d 0 = d − (d + > d − ) when d 0 is varied; it undergoes a supercritical Hopf bifurcation at K 0 = K h when K 0 is varied. For these values, the unforced system oscillates on a stable limit cycle; the asymptotic period of the cycle (evaluated numerically) is T = 1.839 and T = 1.269 when d 0 approaches d + and d − respectively and it is T = 1.528 when K 0 approaches K h . The bifurcation diagrams of the forced system in the (ε, d 0 ) and (ε, K 0 ) planes are given in Fig. 1.
In Fig. 1a, the point H 1 and H 2 on the d 0 -axis correspond to the two supercritical Hopf bifurcations in the unforced system and are roots of the curve h   terminates at a point where AUTO reports an abnormal termination. The point T corresponds to the tangent bifurcation in the unforced system and curve t (1) originates from it.
Two period one fixed points (corresponding to two periodic solutions of the forced system) collide on t (1) then disappear while crossing t (1) to the above. One of the two fixed points is a saddle (if positive); the other loses its stability when it crosses to the region bounded by h 2 and ε = 0 and a stable closed invariant curve appears; in other words, the forced stable cycle of period one of system (4) with d periodically forced bifurcates into a stable torus; it becomes a saddle in the region bounded by f (1) and ε = 1.    Two pairs of period two fixed points appear when t 1 is crossed to the right, two saddle points, two attracting points. The two saddles disappear through the lower branch of f (1) on the left of the intersection of f (1) and t (2) 2 , through t (2) 2 on the right of the intersection. The two attracting points become saddles in the region bounded by f (2) and ε = 1 and disappear through the upper branch of f (1) . When t (2) 2 is crossed to the below, another pair of stable period two fixed points appears and they disappear through the lower branch of f (1) .
In Fig. 1b, the point H on the K 0 -axis corresponds to the supercritical Hopf bifurcation of the unforced system and it is the root of curve h  2 . A tangent bifurcation curve t (1) passes through A 1 . A torus bifurcation curve with period two h (2) starts from B 2 and terminates at a 1:1 resonance A 2 (Fig. 1d); One branch of tangent bifurcation curves with period two t (2) 1 passes through A 2 ; it starts from point T on K 0 -axis at which the stable limit cycle of the unforced system has period two and terminates at a point on f (1) . Another branch t (2) 2 originating from T terminates at another point on f (1) . In addition there is a closed curve f (2) .
The unstable period one fixed point becomes stable crossing h to the below; it becomes a saddle crossing to the region bounded by the closed curve f (1) ; the saddle becomes stable crossing f (1) to the right and it becomes an unstable point on the left of h (1) 2 . The saddle fixed point (if positive) and the non-saddle one collide on t (1) then disappear on the right of t (1) .
Two pairs of period two fixed points appear when t is crossed to the right, two saddles, two non-saddles. The two saddles disappear through the part of f (1) between t (2) 1 and t (2) 2 . The other two points are stable below h (2) , lose stability above h (2) ; they disappear through the part of f (1) on the right of the intersections of t Taking a = 4, b = 0, c = 1, m = 5π, K 0 = 0.92, d 0 = 1.2π, r = 8π, the unforced system undergoes a supercritical Hopf bifurcation at d 0 = d + and a subcritical one For these values, the unforced system oscillates on a stable limit cycle; the asymptotic period of the cycle (evaluated numerically) is T = 1.94 when d 0 approaches d + and it is T = 1.6 with K 0 approaching K h . The bifurcation diagrams of the forced system in the (ε, d 0 ) and (ε, K 0 ) planes are given in Fig. 2a and Fig. 2b respectively; enlargements of some subregions are shown in Figs. 3a, 3b and 3c.
In Fig. 2a, the point H 1 (H 2 ) on the d 0 -axis corresponds to the supercritical (subcritical) Hopf bifurcation in the unforced system and is the root of curve h 2 ) terminates at a 1:2 resonance B 1 (B 2 ). Curve f (1) passes through B 1 and B 2 . T corresponds to the tangent bifurcation in the unforced system and is the root of t (1) . t 2 intersect with f (1) . f (2) is period two flip bifurcation curve.
Two period one fixed points collide on t (1) then disappear while crossing t (1) to the above. One of the two fixed points is a saddle (if positive). The other point loses its stability when crossing to the below of h (1) 1 and a stable closed invariant curve appears; it becomes stable again when h (1) 2 is crossed to the below and an unstable closed invariant curve appears; it becomes a saddle in the region bounded by f (1) and ε = 1.
Two pairs of period two fixed points appear when t 1 is crossed to the right, two saddles, two attracting points. The two saddles disappear through f (1) (between the two intersections of t  with f (1) ) to the above. The two attracting points become saddles in the region surrounded by f (2) and ε = 1 then disappear through the upper part of f (1) . Besides, when f (1) (on the right of the intersection of t (2) 2 with f (1) ) is crossed to the above, two stable points appear then disappear above t The diagram in Fig. 2b is quite the same as that in Fig. 1b. We do not repeat its behaviors here.
For these values, the unforced system has an unstable limit cycle; the asymptotic period of the cycle (evaluated numerically) is T = 1.955 when d 0 approaches d − and it is T = 1.945 with K 0 approaching K h . The bifurcation diagrams of the forced system in the (ε, d 0 ) and (ε, K 0 ) planes are given in Fig. 2c and Fig. 2d respectively; enlargements of some subregions are shown in Fig. 3d and Fig. 3e.
In Fig. 2c, the point H 1 (H 2 ) on the d 0 -axis corresponds to the supercritical (subcritical) Hopf bifurcation in the unforced system and is the root of curve h 2 . Curve f (1) passes trough B. T 1 corresponds to the tangent bifurcation in the unforced system and is the root of t (1) . At T 2 the period of the unstable cycle is 2 and two branches of period two tangent bifurcation curve, t 2 ; a period two torus bifurcation curve h (2) originates from it.
Two period one fixed points collide on t (1) then disappear while crossing t (1) to the above. One of the two fixed points is a saddle (if positive). The other point loses its stability when crossing to the below of h  is crossed to the below and an unstable closed invariant curve appears; it becomes a saddle in the region bounded by f (1) and ε = 1.
Two pairs of period two fixed points appear when t      and f (1) . The other two non-saddles are attracting points in region 4, become unstable in region 3 and disappear when crossing f (1) from region 3 to 1 or from 4 to 2.
In Fig. 2d, the point H on the K 0 -axis corresponds to the subcritical Hopf bifurcation in the unforced system and is the root of curve h (1) . Point B is a 1:2 resonance on h (1) ; curve f (1) passes trough B. At T the period of the unstable cycle is 2 and two branches of period two tangent bifurcation curve, t  Two period one fixed points collide on t (1) then disappear crossing to the right of it. The positive one becomes attracting when crossing to the below of h (1) and an unstable closed invariant curve appears; it becomes a saddle in the region bounded by f (1) .
Two pairs of period two fixed points appear when t is crossed to the right, two saddles, two non-saddles. The two saddles only exist in the region bounded by t and f (1) . The other two non-saddles are attracting below h (2) , become unstable above h (2) and disappear in regions 1 and 2.
For these values, the unforced system has an unstable limit cycle; the asymptotic period of the cycle (evaluated numerically) is T = 1.85 if d 0 approaches d − and it is T = 1.924 when K 0 approaches K h . The bifurcation diagrams of the forced system in the (ε, d 0 ) and (ε, K 0 ) planes are given in Fig. 4a and Fig. 4b respectively.
In Fig. 4a the point H on the d 0 -axis corresponds to the subcritical Hopf bifurcation in the unforced system and is the root of curve h (1) . Point B is a 1:2 resonance on h (1) ; curve f (1) passes through B. At T the period of the unstable cycle is 2 and two branches of period two tangent bifurcation curve, t 2 , a period two torus bifurcation curve h (2) originates from it. There is a tangent curve t (1) above d 0 = 2 not shown in the figure.
Two period one fixed points collide on t (1) then disappear while crossing t (1) to the above. One of the two fixed points is a saddle (if positive). The other unstable point becomes stable when crossing to the below of h (1) and an unstable closed invariant curve appears; it becomes a saddle in the region bounded by f (1) and ε = 1.
Two pairs of period two fixed points appear when t is crossed to the below, two saddles, two unstable points. The two saddles only exist in the region bounded by t 2 , f (1) and ε = 1. The two unstable points become attracting below h (2) , then disappear through t     The diagram in Fig. 4b is quite similar to that in Fig. 2d except that there is a 1:3 resonance C on h (2) and a curve f (2) close to ε = 1. The behaviors are similar to those in Fig. 2d. So we do not repeat here.
Case 5. The unforced system (4) undergoes one supercritical Hopf bifurcation at the local minimum of F (x) whend = d − if d is varied; it undergoes one supercritical Hopf bifurcation if K is varied.
Taking a = 4, b = −1.5, c = 1, m = 5π, K 0 = 1.01, d 0 = 1.95π, r = 4π, the unforced system undergoes a supercritical Hopf bifurcation at d 0 = d − (K 0 = K h ) when d 0 (K 0 ) is varied. For these values, the unforced system oscillates on a stable limit cycle; the asymptotic period of the cycle (evaluated numerically) is T = 1.596 if d 0 approaches d − and it is T = 1.773 with K 0 approaching K h . The bifurcation diagrams of the forced system in the (ε, d 0 ) and (ε, K 0 ) planes are given in Fig. 5a and Fig. 5b respectively.
In Fig. 5a, point H and T 1 on d 0 -axis correspond to the supercritical Hopf bifurcation and the tangent bifurcation in the unforced system; they are the roots of h 2 . At point T 2 on d 0 -axis, the stable limit cycle of the unforced system has period two; from T 2 , two branches of tangent curves t 1 intersects with f (1) at P ; on t 2 . In Fig. 5d, F 1 and F 2 are the intersections of f (1) with ε = 1. Two period one fixed points collide on t (1) then disappear while crossing t (1) to the above. One of the two fixed points is a saddle (if positive). The other point loses its stability when crossing to the above of h is crossed to the right and becomes a saddle in the region bounded by f (1) and ε = 1.
Two pairs of period two fixed points appear when t is crossed to the right, two saddles, two non-saddles. The two saddles only exist in the region bounded by t and ε = 1. The other two non-saddles disappear through the part of f (1) P F 1 ; they are unstable in region 1 bounded by h 2 , f (1) and t (2) 1 , become saddles in region 2 bounded by f (2) and ε = 1; out of regions 1 and 2, they are attracting points. When the part of f (1) OF 2 is crossed to the above, two attracting points appear; they disappear above t 3 . In Fig. 5b, the point H on the d 0 -axis corresponds to the supercritical Hopf bifurcation in the unforced system and is the root of curve h    is crossed to the right, two saddles, two non-saddles. The two saddles disappear through the part of f (1) on the left of OP . The other two non-saddles disappear through the part of f (1) on the right of OP ; they are saddles in the closed region bounded by f (2) ; out of the closed region they are unstable (attracting) above (below) h Taking a = 4, b = 2, c = 1, m = 5π, K 0 = 0.75, d 0 = 0.75π, r = 7.5π, the unforced system undergoes a supercritical Hopf bifurcation at d 0 = d + (K 0 = K h ) when d 0 (K 0 ) is varied. For these values, the unforced system oscillates on a stable limit cycle; the asymptotic period of the cycle (evaluated numerically) is T = 1.912 if d 0 approaches d + and it is T = 1.798 with K 0 approaching K h . The bifurcation diagrams of the forced system in the (ε, d 0 ) and (ε, K 0 ) planes are given in Fig. 6a Fig. 6. Bifurcation diagrams of the forced system with a = 4, b = 2, c = 1, m = 5π, r = 7.5π. The two diagrams in (C) and (D) are partial enlargements of (B).
and Fig. 6b respectively. Some partial enlargements are presented in Fig. 6c and Fig. 6d. In Fig. 6a, point H and T 1 on d 0 -axis correspond to the supercritical Hopf bifurcation and the tangent bifurcation in the unforced system; they are roots of h (1) and t (1) respectively. h (1) terminates at a 1:2 resonance B; f (1) passes through B. At point T 2 on d 0 -axis, the stable limit cycle of the unforced system has period two; from T 2 , two branches of tangent curves t 1 intersects with f (1) . There are flip curves f (2) , f (4) , f (8) ,· · · in Fig. 6a.
Two period one fixed points collide on t (1) then disappear while crossing t (1) to the above. One of the two fixed points is a saddle (if positive). The other point loses its stability when crossing to the below of h (1) 1 and a stable closed invariant curve appears; it becomes a saddle in the region bounded by f (1) and ε = 1.
Two pairs of period two fixed points appear when t is crossed to the right, two saddles, two non-saddles. The two saddles disappear crossing f (1) from region 1 to 2. The other two non-saddles disappear when f (1) is crossed from region 2 to the above of it; in the region bounded by f (2) and ε = 1, they become saddles and two pairs of period 4 fixed points appear. When f (4) (f (8) ) is crossed to the region bounded by f (4) and ε = 1 (f (8) and ε = 1), period 8 (16) fixed points appear, period 4 (8) fixed points become saddles.
The diagram in Fig. 6b is very similar to the one in Fig. 1b. Thus we do not repeat the detailed description including behaviors of fixed points here.
In addition to the six cases analyzed above, there can be other more complicated cases. For example, it would be interesting to investigate the dynamics of the x y z  seasonally forced system when parameters are fixed near the codimension 2 or codimension 3 BT bifurcation. For the difficulty of detecting homoclinic bifurcations using this continuation technique, we will leave these cases for future studies.
4. Bifurcations of the system with periodic d and K. In this section, we consider the periodically forced system (4) with periodic and non-synchronous variations of d and K. To analyze bifurcation diagrams completely, we will not restrict the ranges of ε 1 and ε 2 in the following. That is, we will analyze the periodically forced systems mathematically without considering their biological meanings. For each seasonality mechanism we analyze two cases corresponding to two cases of the unforced system (with stable or unstable limit cycle).
Case 1. d = d 0 (1 + ε 2 cos 2πt) and K = K 0 (1 + ε 1 sin 2πt). The corresponding higher dimensional system that we will analyze is Taking a = 4, b = −1.5, c = 1, m = 5π, r = 4π, d 0 = 1.95π, K 0 = 1.01 the bifurcation diagram of the system (6) is presented in Fig. 7. For these values, the unforced system oscillates on a stable limit cycle. For the periodically forced system (6), its nontrivial periodic solutions only exist in the region bounded by the closed curve t (1) . The forced system has two nontrivial period one solutions, one of which is saddle-type. The other one is unstable in the closed region bounded by h 1 , and f 2 . In addition, in the region where it is unstable the system (6) can also have a stable quasiperiodic solution ((ε 1 , ε 2 ) = (0, 0)) or a chaotic solution. On the other hand, the forced system can have two period two solutions, one of which is saddle-type and only exists in the region bounded by t 2 and t (2) 3 . The other one is unstable in region 2, stable in region 3 and saddle-type in the regions bounded by f Taking a = 4, b = 0, c = 1, m = 5π, r = 2.15π, d 0 = 0.85π, K 0 = 0.92, the unforced system has an unstable limit cycle and the bifurcation diagram of system (6) is shown in Fig. 8. Near t (1) 1 (on the right) and t (1) 2 (on the left) there are two tanscritical curves (denoted by T 1 and T 2 ) not shown. The forced system has two nontrivial period one solutions between t and T 2 , out of which it is stable. In addition, in region 2 the forced system has two period two solutions; one is saddle-type; the other one is stable. In regions 3, 4 and 5, the system has only one period two solution; it is stable in region 3, saddle-type in regions 4 and 5. In a word, periodic solutions of period two only exist in regions 2, 3, 4 and 5.
forced system has two nontrivial period one solutions, one of which is saddle-type. The other one is saddle-type in regions 3 and 6 (inside f 2 ); it is stable in region 8, loses stability bifurcating into a stable quasiperiodic solution (or chaos in some subregion) while crossing h 2 (in regions 1, 2, 4, 5, and 7). In regions 2, 4, 5 and 7, the forced system has two period two solutions, one is saddle-type; the other one is unstable in regions 2 and 5, stable in regions 4 and 7. In regions 3, 6, 9 and 10, there is only one period two solution which is stable in regions 3 and 6, saddle-type in regions 9 and 10. In addition, in regions 9 and 10, there are also period four solutions and chaotic solution in some subregion.
Taking a = 4, b = 0, c = 1, m = 5π, r = 2.15π, d 0 = 0.85π, K 0 = 0.92, the unforced system has an unstable limit cycle and the bifurcation diagram of system (7) is presented in Fig. 10. Near t 2 and T 2 , out of which it is stable. On the other hand, in regions 1 and 2, the forced system has two period two solutions; one is stable; the other one is saddle-type. In regions 3 and 4, there is only one period two solution which is stable.
5. Discussion. In this section, we first analyse the (ε, d 0 ) bifurcation diagrams and the (ε, K 0 ) diagrams respectively, and compare (ε, d 0 ) bifurcation diagrams with (ε, K 0 ) diagrams to see if a "universal diagram" exists for the two different seasonality mechanisms. Then we will discuss about the cases when both d and K are periodically forced.
The bifurcation diagrams when d is periodically forced are shown in Figs. 1a, 2a, 2c, 4a, 5a and 6a. In Fig. 1a, the corresponding unforced system (ε = 0) undergoes two supercritical Hopf bifurcations when d 0 is varied. In this case the diagram is quite similar to the one in Fig. 2a ; even so we can not say the two diagrams are equivalent for the lower Hopf bifurcation in Fig. 2a is subcritical.
In Fig. 2a and Fig. 2c the two diagrams both correspond to the case when the unforced system undergoes one supercritical and one subcritical Hopf bifurcation. What makes the difference is that for the initial parameter values in Fig. 2a the unforced system oscillates on a stable limit cycle with asymptotic period T = 1.94, while in Fig. 2c the unforced system has an unstable limit cycle with asymptotic period T − = 1.955. The asymptotic period of the stable limit cycle in Fig. 2c is T + = 3.743 when d 0 approaches d + . In this case, we can see from the figure that curve h (1) 1 is not connected to f (1) . However this is not universal; if we take a larger r, for example r = 2.5, to make T + closer to 2, h 1 will terminate at a 1:2 resonance B 1 through which h (1) 1 is connected to f (1) . To compare the structure close to the lower Hopf bifurcation with the diagrams in Fig. 4a, we choose to illustrate the case in Fig. 2c here.
In Fig. 4a, there is a tangent bifurcation curve above d 0 = 2 not shown. We can see that the diagram is structurally equivalent to the diagram close to the lower Hopf bifurcation in Fig. 2c.
In Fig. 5a and Fig. 6a, the corresponding unforced system undergoes one supercritical Hopf bifurcation. The difference is that in Fig. 5a, the Hopf bifurcation occurs at the local minimum of F (x) and the stable limit cycle appears above h (1) , while in Fig. 6a it occurs at the local maximum of F (x) and the stable limit cycle appears below h (1) . Apparently the two diagrams are not structurally equivalent. In Fig. 5a, 1 and f (2) , and h (2) 2 connecting f (2) and f (1) . In addition, we get all resonances Fig. 5a including two different 1:2 resonances while in 6a there is only one 1:2 resonance.
In summary, when d is periodically forced, we obtain six bifurcation diagrams that are not all topologically equivalent to each other. We note that in Figs. 2c and 4a the diagrams close to the subcritical Hopf bifurcation are structurally equivalent. The six diagrams are also not equivalent to the universal diagram ( [16] , Fig. 3).
When K is periodically forced, the diagrams are shown in Figs. 1b, 2b, 2d, 4b, 5b and 6b. When the unforced system undergoes a supercritical Hopf bifurcation, diagrams in Figs. 1b, 2b, 5b and 6b are equivalent. When the unforced system undergoes a subcritical Hopf bifurcation the two diagrams in Figs. 2d and 4b are almost the same except that there is t (1) in Fig. 2d and there is a curve f (2) and a 1:3 resonance in Fig. 4b; we can also say they are equivalent in the sense that they both have unclosed f (1) infinitely close to ε = 1 with K 0 large enough.
In short, if K is periodically forced, we obtain two diagrams which are not topologically equivalent, one for supercritical Hopf bifurcation in the unforced system with closed f (1) , f (2) and all resonances, one for subcritical Hopf bifurcation in the unforced system. They are also not equivalent to the universal diagram for the system with Holling type II response function.
The bifurcation diagrams further confirm that seasonality can generate rather complex dynamics. We take the diagram in Fig. 5b as an example. Here we observe that seasonality can support the existence of multiple attractors. In region 1, the attractor is a stable cycle of period 1. Fix the value of ε on the right of A 2 and increase K 0 . A stable cycle of period 2 appears and coexists with the stable cycle of period 1 when t (2) 2 is crossed to the above (below h 2 . Secondly, even very small variations of a parameter can entail a radical change of behavior of the system as some of the bifurcations are catastrophic. In region 1 the system has one stable mode of behavior, a stable cycle of period 1. If K 0 is kept constant and ε is slowly increased, the stable cycle of period 1 becomes a saddle one when curve P B 1 is crossed. Then the system moves toward another attractor, a stable cycle of period 2. Now even if ε is slowly decreased, so that f (1) is crossed from the right, the stable mode of behavior remains the cycle of period 2. If ε is further reduced, we have another catastrophic transition bringing the system back to a period 1 cycle.
For the predator-prey system with response function of Holling type II, different seasonality mechanisms give rise to the same phenomenon. That is to say, the bifurcation diagrams are all equivalent to a universal diagram whichever parameter is periodically forced. In our system with generalized Holling type IV functional response, non-monotonic response function, such a universal diagram does not exist, that is, different periodic parameters have different effect on the system. Firstly the reason for the similarity of the six elementary mechanisms in [16] for system (1) with Holling type II response function lies partially in the similarity of the bifurcation diagrams of the unforced system. Early in 1972, May [14] showed that, for any positive constants c and m, if the predator death rate d satisfies d < mc, then system (1) with response function of Holling type II has a unique interior equilibrium in the positive cone. As K > 0 increases, there exists a unique supercritical Hopf bifurcation which occurs at and the stable limit cycle exists for all K > K h . There is also a transcritical bifurcation which occurs at We can see that the functions describing the bifurcations (8) and (9) have 1-to-1 correspondence between any two parameters; i.e., each parameter can be written as a unique function of the other parameters. As well, K h and K T are either strictly increasing or strictly decreasing as functions of each of the other parameters. Hence it is not surprising that the bifurcation diagrams in the (ε, q 0 ) plane for each of the six elementary mechanisms have similar structure. In our system with nonmonotonic generalized Holling type IV response function, such 1-to-1 relation does not hold any more, which leads to the fundamental differences between (ε, d 0 ) and (ε, K 0 ) diagrams. In other words, in the unforced system when d is varied we can have a tangent bifurcation and one or two Hopf bifurcations while there is only one Hopf bifurcation when K is varied. Secondly, when the unforced system oscillates on a stable limit cycle, the diagrams with K periodically forced have closed flip bifurcation curves f (1) and f (2) while diagrams with d periodically forced do not. Thirdly, when the unforced system has an unstable limit cycle, diagrams with different seasonality mechanisms are also not equivalent for the same reason as in the first and as f (1) may be infinitely close to ε = 1 or not. To our knowledge, the case when the system is forced with periodic and nonsynchronous variations of two parameters has not been investigated. We consider this case with two different seasonality mechanisms: d and K are varied periodically and non-synchronously with the same period; or the period of K is twice the period of d. Comparing Fig. 7 with Fig. 9 and Fig. 8 with Fig. 10, we can see that for the case when the unforced system has a stable limit cycle, the strong 1:2 resonances in Fig. 7 are different than those in Fig. 9. That is, the coefficients of the corresponding normal form are of opposite sign (see [11]). For the case when the unforced system has an unstable limit cycle the system (7) can not have chaotic solutions through cascade of period doublings while the system (6) can. What's more we can conclude that these differences originate from the two different seasonality mechanisms. When d and K are periodic with different and varied periods ω 1 and ω 2 , it would be interesting to study how the bifurcations are related to the varied periods. This project will be investigated in our future work.
In conclusion, a predator-prey system with generalized Holling type IV functional response is studied, and four elementary seasonality mechanisms are analysed by means of a continuation technique. The results are interpreted in detail and compared with each other. When the natural death rate of predators d varies periodically, we get six different bifurcation diagrams corresponding to different bifurcation cases of the unforced system. If the carrying capacity K is periodic, two different bifurcation diagrams are obtained, one corresponding to the supercritical Hopf bifurcation case, one to the subcritical case. We also give several bifurcation diagrams when both d and K are periodically varied. The results show that the periodically forced system has more complex dynamics than the unforced one even at very low degrees of seasonality. For example, the forced system can have different attractors such as periodic solutions of various periods, quasiperiodic solutions and chaotic solutions through torus destruction or cascade of period doublings, while the unforced system can only have attractors like attracting equilibria and limit cycles. In addition, the forced system can switch between the different attractors and catastrophic transitions can occur if suitable parameters are varied. Furthermore, the two elementary seasonality mechanisms when d or K is periodic lead to different bifurcation diagrams. In other words, a "universal diagram" for the periodically forced system does not exist, unlike the one in the seasonally forced system with Holling type II response function. Finally, we give some orbits in phase space and corresponding Poincaré sections of different types of solutions by MATLAB numerical simulations in Figs. 11-12. Quasi-periodic solution, chaos obtained from torus destruction and cascade of period doublings are shown in Fig. 11 , corresponding to the bifurcation case in Fig. 2a. An interesting special case is illustrated in Fig.  12, corresponding to the case in Fig. 4a. In our numerical study, we generate similar figures to illustrate behavior in other cases, but they are not shown to avoid redundancy.