STRANGE ATTRACTORS IN A PREDATOR–PREY SYSTEM WITH NON-MONOTONIC RESPONSE FUNCTION AND PERIODIC PERTURBATION

. A system of ordinary diﬀerential equations of a predator–prey type, depending on nine parameters, is studied. We have included in this model a nonmonotonic response function and time periodic perturbation. Using numerical continuation software, we have detected three codimension two bifurcations for the unperturbed system, namely cusp, Bogdanov-Takens and Bautin bifurcations. Furthermore, we concentrate on two regions in the pa- rameter space, the region where the Bogdanov-Takens and the region where Bautin bifurcations occur. As we turn on the time perturbation, we ﬁnd strange attractors in the neighborhood of invariant tori of the unperturbed system.

The model itself has been modified to make it more realistic, or to include other behaviours which have been observed in nature. Some of those variations include generalization of the response functions as were introduced by C.S. Holling [12,13]. Other types of variations are to include space-dynamics in the system (see for example in [14,21]) or delay (see [15,22]), or hyberbolic mortality rate (see [21]).
In [23] a predator-prey type of dynamical system with a nonmonotonic response function is considered. This response function includes the so-called group defence mechanism, i.e., the success rate of the predation decreases as the density of the prey becomes larger. In [4] a system similar to [23] with an addition of a competition factor among the predators is studied.
We study a predator-prey type of dynamical system similar to the system in [4] with a focus on describing the mechanism that creates strange attractors due to the presence of the time periodic perturbation. This type of perturbation was introduced by Rinaldy et al [18] in order to model the seasonal variation of the carrying capacity of the environment. (For example, the volume of water in a pond varies periodically between the rainy season and the dry season.) Strange attractors have received a lot of attention in the literature. The existence of a strange attractor in a predator-prey system might increase the unpredictability of the behaviour of the populations. As is known, two-dimensional systems of time-independent differential equations cannot be chaotic. Thus, adding a timedependent perturbation might change the situation completely.
The predator-prey type of system which we are studying has interesting dynamics. In [10,11] we have observed the occurrence of an interesting swallowtail bifurcation of a periodic solution. In the absence of the periodic perturbation, three out of four bifurcations from a list of codimension two bifurcations of equilibria in [16] are present in the system, i.e., the cusp bifurcation, Bogdanov-Takens bifurcation, and Bautin bifurcation. These bifurcations are connected to each other in a codimension three point as is realized in [4].
In this paper, we produce numerically constructed phase-portraits for these bifurcations from the original system of differential equations. This turns out to be quite challenging due to the fact that the values of some of the parameters are small. This numerically constructed phase-portrait helps in determining the area in the parameter space where we suspect that a strange attractor exists.
Outline of this paper. We will start with a formulation of the problem. The original system we are studying is a system of two ordinary differential equations which is periodically perturbed. For our purpose, it is more convenient to work with an autonomous system. For that reason we attach to the original system another system of two ordinary differential equations where the periodic perturbation term can be seen as one of the component of a solution to the attached system for a particular initial condition. To study the creation of the strange attractor, we consider the time-2π stroboscopic map of the four-dimensional system.
We proceed with the analysis by studying the bifurcations in the unperturbed system (the system with no periodic perturbation). In this paper we concentrate on describing the dynamics and bifurcations in two regions in parameter space, namely the region where we have Bogdanov-Takens bifurcation (region B) and a region where we have Bautin bifurcation (region A). The phase-portraits of the system in these regions are numerically produced to describe the dynamics of the system there.
In Section 4, we include the perturbation by letting the perturbation parameter become nonzero. Thus, we turn into the stroboscopic map and produce the phase portrait of the two dimensional map, also numerically. Two strange attractors are observed in the system, which correspond to the stable and unstable tori in the unperturbed system.

Problem formulation.
2.1. The model. Let x and y denote the density of prey and predator populations, respectively. Consider a system of two-ordinary differential equations, i.e.

471
The natural birth rate of the prey has been scaled to one, while 1 λ corresponds to the carrying capacity of the environment. The natural death rate of the predator is δ, and the competition rate between predators is µ. We have assumed also that the parameter that measures the rate of the conversion from captured prey to predation has been scaled to one.
The function P (x) is called the response function. In [12,13] Holling suggested three types of functional responses. The first type is a strictly increasing monotonic linear function when 0 ≤ x ≤ C for some C, and then constant for x ≥ C. The second and the third type of response functions are also increasing functions with respect to the density of the prey, but the limit of the function as the density become very large is constant. The difference between the the second and the third type is that the third type goes to zero faster than the second type as the density goes to zero. Zhu et al [23] generalized the response function to a non-monotonic response function (type IV).
Following [11] (see also [4,23]), we have used the so-called type IV functional response. For this type of functional response, there exists a critical value such that the predation is increasing when the density of the prey is lower and decreasing for larger density than that critical value. Thus, the prey has power over the predator when they are grouped together (the so-called group defence mechanism). The parameter that controls this effect is α, while β measures the saturation rate of the predator.
We also include periodic perturbation in the system, by setting: which brings system (1) into a nonautonomous system of ordinary differential equations. This type of perturbation was introduced in [18,11].

2.2.
The autonomous version of the model. For our purpose, it is more convenient to work with a system of autonomous differential equations. This is done by attaching to (1) a particular system of ordinary differential equations, i.e.: Solving this system with initial condition u(0) = 0 and v(0) = 1, we have: (u, v) = (sin(ωt), cos(ωt)) as an asymptotically stable solution. The system then can be represented as:ẋ Furthermore, by restricting ourselves to the initial condition: (x(0), y(0), u(0), v(0)) = (x • , y • , 0, 1) and projecting the solution to (x, y)−plane, we derive back system (1) with the periodic perturbation λ = λ 0 (1 + ε sin ωt). Moreover the stable nature of this solution prevents us from having a growing numerical error as we employ a numerical integration method to integrate (2). This extension to a four-dimensional system brings the original system to a higher dimensional autonomous system of ordinary differential equations. One of the advantages of this extension is that we can prove a classical result on the period of periodic solution of a system of differential equations which is defined by a periodic vector field, in a rather straightforward manner (see [11]). (2). Let ξ(t; ξ 0 ) = (x(t), y(t), u(t), v(t)) be the solution of System (2), with initial condition: ξ(0) = ξ 0 = (x 0 , y 0 , 0, 1). Using this, consider the time-2π solution map:

The Poincaré map for System
Furthermore, let us consider a projection mapping: defines a stroboscopic map (or also known as Poincaré map) where the bar denote the time-2π upshift. This map P ϕ is used to study the dynamics of System (2) for ε = 0. Clearly, the phase-portrait for System (2) can be obtained by applying an S 1 −action to the phase portrait of (4). In particular, for ε = 0, the phase-portrait of (4) is similar to the phase-portrait of (1). We proceed with the analysis by setting ε = 0 in (2), which is equivalent to analyzing system (1) without the periodic perturbation. Note that we have two equilibria of System (2) for ε = 0, i.e.: (0, 0) and ( 1 λ0 , 0). The equilibrium (0, 0) is saddle type with linear eigenvalues: 1 and −δ. The eigenvalue 1 correspond to the situation of non-existing predator, thus the prey will grow exponentially towards the second equilibrium: ( 1 λ0 , 0). The eigenvalue −δ corresponds to the situation of the absence of prey; the predator will go to extinction. The stability of the equilibrium (0, 0) is independent to the chosen bifurcation parameters α and β.
The stability of the second equilibrium ( 1 λ0 , 0), however, depends on both bifurcation parameters. The eigenvalues of the linearized vector field in the vicinity of the equilibrium ( 1 λ0 , 0) are: −1 and Interesting dynamics occurs when the equilibrium ( 1 λ0 , 0) is unstable. Since both of the equilibria found so far are unstable, then another equilibrium (or other invariant structure, such as periodic solution) can be found. For α = 0.002, and β = 0.25 we found another equilibrium: (x, y) := (1.816, 1.434).
Using a numerical continuation software AUTO [6], we will follow this equilibrium while varying the bifurcations parameters α and β. 3. The dynamics of System (1) without periodic perturbation. Our focus in this paper is on describing the dynamics and bifurcations of system (1) (or similarly System (2)). In particular, we are interested on the creation of a strange attractor (or repeller) in the system. Following [11], we fix a value for some of the parameters, i.e. δ = 1.1, λ 0 = 0.01, µ = 0.1, and ω = 1. We have chosen the parameter α and β for the bifurcation parameters. In this section, we consider the dynamics of the system (1) for ε = 0. Following the equilibrium (x, y) := (1.816, 1.434), which is found at the parameter values: α = 0.002 and β = 0.25, by varying α we found a fold bifurcation. The locus of this bifurcation can be continued by varying α and β to construct the fold curve F 1 (see Figure 1). On this fold curve we have found a codimension two bifurcation where the equilibrium has double zero eigenvalues. We have indicated the region on the (β, α) plane where this bifurcation occurs as B in Figure 1.
There are two branches of fold bifurcation curves which are labelled by F 1 and F 2 . The two curves coalesce at the cusp bifurcation point. The curve which is plotted using a thickened solid line corresponds to Hopf bifurcations for one of the equilibria in the interior of the phase space. This curve is divided into two parts which are labelled by H 1 and H 3 , while H 2 is a point between the two branches of Hopf bifurcation curves.
The curve plotted using the dashed line is homoclinic bifurcation curve. This curve is branching out from a codimension two bifurcation point BT which is located in the domain labelled by B, i.e.: This domain is magnified and the bifurcation diagrams in this area are plotted in Figure 2 (the diagram on the first row). We have divided the domain B into several parts: B 1 , B 2 , B 3 , B 4 , and B 5 , which are bounded by the bifurcation curves. The phase portrait of system (1) for ε = 0, on each of these domain is plotted in Figure 2. The transition from B 1 , to B 2 , then to B 5 , and then to B 3 , corresponds to the Bogdanov-Takens bifurcations. In Appendix A.1 we have shown that indeed the point BT corresponds to a Bogdanov-Takens bifurcation.
There is a part of the curve Hom which coincides with the F 2 -curve. This part is labelled by F 2 and it corresponds to a the occurrence of an orbit homoclinic to a degenerate equilibrium with one zero eigenvalue.
Let a region A be the area defined by (β, α) ∈ [0.02, 0.12] × [−0.001, 0.013]. In this region we have another codimension two point H 2 , at where an equilibrium of system (1) for ε = 0, undergoes the so-called: Bautin bifurcation or Degenerate Hopf bifurcation (see [16]). Branching out from the locus of this Bautin bifurcation there is the curve indicated by FLC in Figure 1. On this curve the periodic solution created via the Hopf bifurcation undergoes a fold bifurcation of limit cycles. See also Figure 3 for a more details on the bifurcation sets and the phase portrait on each of the regions in the neighborhood of the Bautin bifurcation H 2 -point. In Appendix A.2 we present the analytic confirmation of the Bautin bifurcation.
In the next section we will take the periodic perturbation into consideration. An immediate consequence of this periodic perturbation is that a non-degenerate equilibrium system (1) for ε = 0 turns into a periodic solution as ε = 0. Our interest, however, is to look for strange attractors in System (2), either for positive time or for negative time. 4. Strange attractors and chaos. The values for the parameters that we have chosen are: δ = 1.1, λ 0 = 0.01, µ = 0.1, ω = 1, α = 0.007, and β = 0.08. This choice corresponds to the region in parameter space (β, α), labeled by A 5 in Figure 3 (see also the diagram labelled by A 5 in Figure 4 for the corresponding phase-portrait).
We start integrating System (2) at an initial condition which is close to the domain of attraction of the attractor for some period of time to make sure that we have reached the attractor. Then we derive the corresponding variational equations in the neighborhood of this attractor, and integrate this equation for some period of time to compute the Lyapunov characteristic exponent for the attractor. The computation is stopped after observing that five nonzero digits of the exponents are no longer varying.
A negative time strange attractor. The invariant structures of System (2) for ε = 0 consist of two stable equilibria (plotted using dotted circles), an unstable periodic solution (plotted using dashed line), and a saddle type equilibrium (plotted using the thickened cross sign); see the diagram labelled by A 5 in Figure 4. Recall that the unstable periodic solution surrounds one of the stable equilibria.  B 1 , B 2 , B 3 , B 4 , and B 5 . The diagrams on the second and third rows are the phase portraits in each of these regions. The four phase portraits for (α, β) in B 1 , B 2 , B 5 , and B 3 , correspond to the Bogdanov-Takens bifurcations. The transition from phase portrait when (β, α) ∈ B 3 to (β, α) ∈ B 4 , or when (β, α) ∈ B 5 to (β, α) ∈ B 4 corresponds to fold bifurcation of equilibrium of system (1) for ε = 0; the latter with the creation of orbit homoclinic to the degenerate equilibrium.   Figure 3. We have plotted the magnification of the region A of Figure 1. We have indicated seven regions on that diagram, i.e. A j , j = 1, 2, . . . , 7. These regions are separated from each other by the bifurcation curves: F 2 -where a fold bifurcation occurs-, H 1,3 -where a Hopf bifurcation occurs-, Hom -where a homoclinic bifurcation occurs -, and FLC -where a fold bifurcation of limit cycles occurs.
By applying an S 1 −action to the periodic solution and the stable equilibrium, we have constructed an unstable invariant torus and a stable periodic solution, for (2) for ε = 0.
For ε = 0.07, we have performed a backward integration for System (2) starting at an initial condition in the neighborhood of the stable periodic solution. The solution then converges to a neighborhood of the invariant torus. By doing this we have found an evidence for the existence of a strange negative-time attractor or strange repeller. See Figure 5 for the plot of this strange repeller. A similar procedure can be done also for other situations, i.e. in Figure 4 for the diagrams labelled by A 1 , and A 6 .
The destruction of a stable torus. Let us now consider the situation in a region in Figure 3 which is labelled by A 1 . In Figure 6, we have plotted Poincaré section for System (2) for ε = 0 (the diagrams on the left hand side) and ε = 0.07 (the diagrams on the right hand side). On these diagrams we focus on describing the situation near the stable torus, which is created from applying an S 1 −action to the stable periodic solution for the ε = 0 case. We have observed that again, as ε become nonzero, the stable torus starts to break-up into a positive time strange-attractor.
Note that the strange negative time attractor which is described previously is related to the unstable torus. The size of this negative time attractor in phase space is smaller than the positive time attractor. Thus, in Figure 6, the Poincaré section of this negative time strange attractor looks like an ordinary torus. Furthermore, as  Figure 4. We have plotted seven diagrams which correspond to the phase portraits of system (1), for parameter value of (β, α) in region: A j , j = 1, 2, . . . , 7 and ε = 0. The topological changes of these phase portraits when the parameter moves from A 1 , to A 2 , to A 3 , and back, are in agreement with the scenario of Bautin (or degenerate Hopf) bifurcation.
β approaches the Hopf line (see region A 1 in Figure 3), the negative time strange attractor become unstable torus and then collides into the stable periodic solution.
Remark 1. On the positive Lyapunov exponent of the attractor. In Table 1 we have listed a positive Lyapunov exponent and Kaplan-Yorke dimension of the attractors found for several values of β and fixed α = 0.005. All of the computed exponents are rather small, i.e. of order 10 −5 . We suspect that this due to the fact that the values of the parameters of our interest in this paper are rather small.  Figure 5.  Note that the negative attractor is considered as a positive attractor for negative time dynamical system. Thus we multiply the vector field in System (2) by −1 and then integrate it in the neighborhood of its positive attractor.
Transient chaotic behaviour. It is evident in the System (2) for ε = 0 that we have found a fold bifurcatino of limit cycles, where a stable limit cycle collides with an unstable limit cycle. Exactly at the bifurcation point, a degenerate periodic solution which is attracting from one side and repelling from the other side is created. At ε = 0.07, in the neighborhood in the parameter space of the fold of limit cycle bifurcation, a transient behaviour is created.
To illustrate the above described behaviour, let us consider a solution which is started at an initial condition at (70, 30, 0, 1). In a rather short time (less then 50 units of integration time), the solution is attracted to the neighborhood of the torus which is a product of the degenerate periodic solution and an S 1 -action in (u, v)-direction. The solution stays in the neighborhood of the torus for a long time, i.e. for at least 800 units of integration time. Finally, the solution leaves the neighborhood of the degenerate torus, as it is captured into the domain of attraction of the stable periodic solution. In Figure 7 we have plotted the function log y(t) against log x(t) using the dotted line. On the same diagram, we have plotted the iteration of the Poincaré map for the corresponding solution using the * symbol. From the iteration of the Poincaré map, we have observed that there are three locations on the degenerate torus which are visited more often by the solution. For the diagram on the right hand side in Figure 7, we have plotted the function x(t) against t. In this diagram the transient behaviour is evident between the times 200 and 1000.
The biological relevance of this study. One can think of the mathematical analysis, which is done here and also in [11], as providing a map in the parameter space which indicates locations at where interesting behaviour can be expected. In real life applications however, most of the time one can only observed qualitatively (not quantitatively) similar behaviour as is observed here. This mathematical study then provide a plausible explanation for the origin of such observations in real life applications.
Consider Figure 4 for example. From the diagram labeled as A 1 , one can see that there is an asymptotically stable periodic solution. However, this solution is not the only stable invariant structure in the system; there is also a stable equilibrium. The implication of this is that in some part of phase-space the densities of the population of prey and predator will converge to limiting values. In our case, these limiting values can be quite small. On the other part of phase-space the densities of each population are attracted to periodic orbits. Note that the boundary of the domains of attraction of these two stable structures is an unstable periodic solution. This means that we can have two initial conditions which are close to each other, but the limiting behavior of the corresponding solution can be very different.
In particular, the chaotic transient behaviour we have observed here might be interesting from the application point of view. If we start from the initial condition with large density of predator and density of prey, after a short time the solution will be captured into non-periodic oscillations and stay there for a longer period of time. Thus, the chaotic behaviour occurs for a rather long period of time, but after that it disappears when the solution is captured into a stable attracting equilibrium. Similar behaviour is observed in [20] where the author considered a toy model for atmospheric regime transitions.
Furthermore, we have observed also in this paper that in some regions of the parameter space, we have chaotic oscillations which are caused by destruction of a torus of the unperturbed system. Some preliminary studies that we have done indicate that the size of these chaotic solution is rather small. Thus, in real life situation, this might not be observable.
One last remark is about the size of the parameter ε which is chosen to be small in this paper. In the case where the periodic variation models the seasonal variation of the volume of water in a pond, the variation can be rather small in compare to the volume of the pond. Thus, this study is relevant for that situation. We note also that for ε = 0 the system is two-dimensional and autonomous, so that no chaotic behaviour is possible. Thus, having the periodic variations in the carrying capacity is essential for the existence of the chaotic behaviour.
Appendix A. Analytical confirmation for some of the codimension two bifurcations . In this appendix we will follow closely the description of the Bautin and Bogdanov-Takens bifurcations in the book [16] to prove the existence of the bifurcation. Let us fix the parameters at the value of: λ = 0.01, δ = 1.1, and µ = 0.1, in (2). The corresponding bifurcation points and their related equilibria are computed using the numerical continuation software AUTO97 [6]. Linearizing (1) in the vicinity of (x 0 , y 0 ) gives double zero eigenvalues up to the order of 10 −6 . The linearized matrix reads These vectors are chosen as solutions of : , is a standard dot product in R 2 . Using ξ 0 and ξ 1 as bases, we can write: x y = w 1 ξ 0 + w 2 ξ 1 .