`Discrete Dynamics in Nature and SocietyVolume 2011 (2011), Article ID 862494, 27 pageshttp://dx.doi.org/10.1155/2011/862494`
Research Article

## Codimension-Two Bifurcations of Fixed Points in a Class of Discrete Prey-Predator Systems

1Department of Applied Mathematics and Computer Science, Shahrekord University, P.O. Box 115, Saman Road, Shahrekord 88186-34141, Iran
2Department of Applied Mathematics and Computer Science, Ghent University, Krijgslaan 281-S9, 9000 Gent, Belgium

Received 9 December 2010; Revised 7 March 2011; Accepted 20 March 2011

Copyright © 2011 R. Khoshsiar Ghaziani et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The dynamic behaviour of a Lotka-Volterra system, described by a planar map, is analytically and numerically investigated. We derive analytical conditions for stability and bifurcation of the fixed points of the system and compute analytically the normal form coefficients for the codimension 1 bifurcation points (flip and Neimark-Sacker), and so establish sub- or supercriticality of these bifurcation points. Furthermore, by using numerical continuation methods, we compute bifurcation curves of fixed points and cycles with periods up to 16 under variation of one and two parameters, and compute all codimension 1 and codimension 2 bifurcations on the corresponding curves. For the bifurcation points, we compute the corresponding normal form coefficients. These quantities enable us to compute curves of codimension 1 bifurcations that branch off from the detected codimension 2 bifurcation points. These curves form stability boundaries of various types of cycles which emerge around codimension 1 and 2 bifurcation points. Numerical simulations confirm our results and reveal further complex dynamical behaviours.

#### 1. Introduction

The dynamic relationship between predators and their prey is one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance; see [16]. The prototype Lotka-Volterra predator-prey system has received a lot of attention from theoretical and mathematical biologists; see [26]. This model is described by the following system of ordinary differential equations: where and represent the densities of the prey and the predator, , , , and are the intrinsic growth rate of the prey, the capture rate, the death rate of the predator, and the conversion rate, respectively, and denote the intraspecific competition coefficients of the prey and the predator, and is the carrying capacity of the prey.

This and related models have been studied intensively in the previous decades and it has been noted, in particular, that they present a wealth of bifurcations of various codimensions. Organizing centers of codimension 3 have been investigated in great detail in [7] and references therein. For background on bifurcation theory, we refer to [8, 9].

Far-reaching generalizations of the model have been studied both analytically and numerically. We refer in particular to [10] where a five-parameter family of planar vector fields is studied that takes into account group defense strategy, competition between predators and a nonmonotonic response function. Though mathematical analysis is the prime tool also in this case, numerical methods are indispensable for a detailed study. In [10] not only the general-purpose languages Mathematica and Matlab were used, but also the dedicated bifurcation packages AUTO [11] and MatCont [12].

However, in recent years, many authors [46, 13] have suggested that discrete-time models are more appropriate than continuous ones, especially when the populations have nonoverlapping generations. Furthermore, discrete-time models often provide very effective approximations to continuous models which cannot be solved explicitly.

Though discrete models are by no means less complex than continuous ones, they have the computational (numerical) advantage that no differential equations have to be solved. In fact, the state-of-the-art in computational bifurcation study is in some respects more advanced for discrete systems than for continuous ones. So far, normal form coefficients for codimension 2 bifurcations of limit cycles are not computed by any available software. For discrete systems, such coefficients can be computed and they can even be used to start the computation of codimension 1 cycles rooted in such points; see [12, 14]. We will exploit this possibility.

In this paper we investigate a map in the plane for a pair of populations, first studied in [15, 16]. Danca et al. [15] investigated Neimark-Sacker bifurcations numerically. Murakami [16] discussed its branch points, flip bifurcations, and Neimark-Sacker bifurcations, establishing the sub- or supercritical character of the flip and Neimark-Sacker bifurcations by an explicit reduction to the center manifolds, obtaining a prediction for the invariant curve that branches off at the NS point. We undertake an analysis of the dynamics of the iterates of the map. This paper is organized as follows. In Section 2, we first discuss the map and the stability and bifurcations of its fixed points, proving the results in [16] by direct computations of the normal form coefficients. Next, in Section 3, we numerically compute curves of fixed points and bifurcation curves of the map and its iterates up to order 16 under variation of one and two parameters. We compute the critical normal form coefficients of all computed codim-1 and codim-2 bifurcations. These coefficients are powerful tools to compute stability boundaries of the map and its iterates and to switch to other bifurcation curves. In Section 4, we study more details on the bifurcation scenario of the system around a Neimark-Sacker bifurcation point. We conclude our work in Section 5 with a discussion of the obtained results.

#### 2. The Model, the Fixed Points, Their Stability, and Bifurcations

In this paper, we study the map which is analogous to (1.1) for the case of predators and prey with nonoverlapping generations; see [16, 17]. It can also be seen as an approximate discretization of the continuous-time Lotka-Volterra model which is a simplification of (1.1) in which the intraspecific competition of the predator is ignored and the carrying capacity of the prey is 1. and represent the densities of the prey and predator, and , , , and are nonnegative parameters. Applying the forward Euler scheme to system (2.2) with the stepsize and assuming , we obtain the map (2.1) with nonnegative parameters , , and .

We note that in (2.1) one of the two parameters can be removed by a rescaling of . So, the system is in reality a two-parameter system and fully equivalent to the system studied in [16]. We will generally choose , as the unfolding parameters in the bifurcation study. In a way, it is the simplest possible discrete predator-prey model and, therefore, allows a reasonably complete analytical treatment as far as the fixed points of the map are concerned. However, we will see that even in this case the behaviour of cycles is very complicated and can only be studied by numerical methods.

We naturally start the bifurcation analysis of (2.1) with the calculation of the fixed points. These are the solutions to The origin is a fixed point of (2.1) but is of little interest. Two further nontrivial fixed points are and . We note that is biologically possible only if its coordinates are nonnegative, that is, . is biologically possible only if the following conditions hold: (i),(ii).

We start the local bifurcation analysis of the map (2.1) by linearization of around each of its fixed points. The Jacobian matrix is given by

The characteristic equation of is given by where and .

##### 2.1. Stability and Bifurcation of

Proposition 2.1. The fixed point is asymptotically stable for . It loses stability via branching for and there bifurcates to .

Proof. Eigenvalues of the Jacobian at are and 0. So, is stable if and loses stability at . In -space forms the curve with tangent vector . is represented in -space by the curve . When , these curves intersect at and the tangent vector in -space is , so it is clear that branches to for .

##### 2.2. Stability and Bifurcation of

The Jacobian matrix of (2.1) at is given by

Proposition 2.2. The fixed point is linearly asymptotically stable if and only if and . Moreover, it loses stability: (i)via branching for and there bifurcates to , (ii)via branching for and there bifurcates to if , (iii)via a supercritical flip for if .

Proof. The multipliers of are and . The fixed point is asymptotically stable if and only if and , that is, if and only if and . Boundary points of the stability region must satisfy one of three conditions: , , or .
In the first case, the conditions and are satisfied for nearby values , hence this is a real stability boundary. In Proposition 2.1, we proved that this is a branch point and the new branch consists of points.
In the second case, this is a stability boundary only if . The Jacobian (2.6) then has an eigenvalue +1 and it is checked easily that these boundary points are also points.
In the third case, this is a stability boundary only if . In this case, which means that loses stability via a period doubling point. For supercriticality of the period doubling point, it is sufficient to show that the corresponding critical normal form coefficient , derived by center manifold reduction is positive; see [9], Ch. 8 and [14]. Here, , and , are the second- and third- order multilinear forms, respectively, and and are the left and right eigenvectors of for the eigenvalue , respectively. These vectors are normalized by , , where is the standard scalar product in . We obtain The components of the multilinear form are given by where the state variable vector is for ease of notation generically denoted by instead of .
Let , then we have and find The third-order multilinear form is identically zero. The critical normal form coefficient is given by which is clearly positive. This completes the proof of supercriticality of the flip point at .

The stability region of , as obtained in Proposition 2.2, is shown in Figure 1.

Figure 1: Stability regions in -space. , are the stability regions of and , respectively. The stability region of is the union of and (which correspond to (i) and (iii) in Proposition 2.3, resp.), and the open interval that separates them (and corresponds to (ii) in Proposition 2.3).
##### 2.3. Stability and Bifurcation of

To study the stability of , we use the Jury's criteria; see [6, Section  A2.1]. Let be the characteristic polynomial of .

According to the Jury’s criteria, is asymptotically stable if the following conditions hold:

At , we have: We note that and are independent of .

Proposition 2.3. is linearly asymptotically stable if and only if one of the following mutually exclusive conditions holds: (i) and , (ii) and , (iii) and .

Proof. The criterion is easily seen to be equivalent to the condition or equivalently,
Next, the criterion is easily seen to be equivalent to
The criterion translates as

In Figure 1, the parts of the 3 curves , , and where and are positive are depicted (resp., for , , and ). The curves and intersect solely at and the curves and at . It follows from the figure that the 3 inequalities are fullfilled when (i), (ii), or (iii) holds. This can also easily be shown algebraic. In Figure 1, the stability region of is thus the union of and (which correspond to (i) and (iii) in this proposition, resp.), and the open interval that separates them (and corresponds to (ii)).

Proposition 2.4. loses stability: (i)via a flip point when and , (ii)via a Neimark-Sacker point when and , (iii)via a branch point when and where it bifurcates to , (iv)via a branch-flip (BPPD) point when and , (v)via a resonance 1 : 2 point when and .

Proof. By Proposition 2.3 the stability boundary of consists of parts of three curves, namely, (1)Curve 1: . (2)Curve 2: . (3)Curve 3: . The points of Curve 1 which are on the stability boundary of satisfy , that is, they have an eigenvalue −1 and, thus, are period doubling points. The points of Curve 2 which are on the stability boundary satisfy , that is, they have two eigenvalues with product 1 and, thus, are Neimark-Sacker points. The points of Curve 3 which are on the stability boundary satisfy , that is, they have an eigenvalue 1. It can be checked easily that then branches to .
Combining this with Proposition 2.2 we find that the interior points of the boundary parts of Curves 1, 2, and 3 form the sets described in parts (i), (ii), and (iii) of Proposition 2.3, respectively.
At the intersection of Curves 1 and 3 (when and ), the point is a branch-flip point (BPPD) with eigenvalues −1 and 1. The intersection point of Curves 1 and 2 (when and ) is a resonance 1 : 2 point with double eigenvalue −1.

Proposition 2.5. The flip point in Proposition 2.4, part (i), is subcritical.

Proof. For the subcriticality of the flip bifurcation, it is sufficient to show that the normal form coefficient in (2.7) is always negative. By using the same procedure as we used in the proof of Proposition 2.2, we obtain To simplify computations, we do not normalize and , since rescaling does not change the sign of provided that is positive (it can be proved easily that this is the case). is computed as
Let , then we have and find The third-order multilinear form is zero. The critical normal form coefficient , which is clearly negative for . This completes the proof of subcriticality of the flip point.

Proposition 2.6. The Neimark-Sacker point in Proposition 2.4, part (ii), is supercritical.

Proof. To prove the supercriticality of the Neimark-Sacker point, it is sufficient to show that the real part of the corresponding critical normal form coefficient , is negative; see [14]. Here, is the argument of the critical multiplier, , and , are the second- and third-order multilinear forms, respectively. and are the left and right eigenvectors of These vectors are normalized by , where is the Hermitian inner product in . In the Neimark-Sacker point, (, so ). We get with characteristic polynomial . It follows that and , so . The eigenvectors are where is not normalized since rescaling of does not change the sign of (this can be proved easily) and we can simplify the computations by not normalizing. The calculation of the components of gives us and it follows that The remaining calculations are done with the help of Maple. The calculation of gives us where For and , we obtain The third-order multilinear form is identically zero, so we get The real part of equals (exact up to a positive factor due to the nonnormalization of ), which is clearly negative. This completes the proof of supercriticality of the Neimark-Sacker point in Proposition 2.4, part (ii).

The stability regions (i) and (iii) of obtained in Proposition 2.3 are depicted as and , respectively, in Figure 1. The region (ii) is the open interval on the common boundary of and . We have a complete description of the stability region of for all parameter combinations.

The biological interpretation of Figure 1 is as follows. The situation with no prey and no predators exists for all parameter values but is stable only in , that is, for . The situation with a fixed number of prey but no predators exists for all but is stable only in . Coexistence of fixed nonzero numbers of prey and predators is possible whenever and but is stable only in the union of and .

#### 3. Numerical Bifurcation Analysis of and

In this section, we perform a numerical bifurcation analysis by using the MATLAB package Cl_MatContM; see [12, 14]. The bifurcation analysis is based on continuation methods, whereby we trace solution manifolds of fixed points while some of the parameters of the map vary; see [18]. We note that in a two-parameter setting the boundaries of stability domains of cycles are usually curves of codimension 1 bifurcations, which necessarily have to be computed by numerical continuation methods.

To validate the model, we compute a number of stable cycles for parameter values in the regions where we claim they exist; these stable cycles were found by orbit convergence, independently of the continuation methods. They are represented in the Figures 2, 6, 11, and 14. Similarly, a stable closed invariant curve was computed by orbit convergence and is represented in Figure 4. The computation of an Arnold tongue in Section 3.5 by a continuation method also validates our computations.

Figure 2: A stable 8-cycle of (2.2) for , , .
##### 3.1. Numerical Bifurcation of

By continuation of starting from , , , in the stable region of with free, we see that is stable when . It loses stability via a supercritical period doubling point (PD, the corresponding normal form coefficient is ) when , and via a branch point (BP) when crosses 1.4. The output of Run 1 is given by:label = BP, x = (0.000000 0.000000 1.000000) label = PD, x = (0.666667 0.000000 3.000000) normal form coefficient of PD = 9.000000e+000

The first two entries of x are the coordinate values of the fixed point , and the last entry of x is the value of the free parameter at the corresponding bifurcation point. We note that the normal form coefficient of the PD point is 9, confirming (2.11). We note furthermore that the detected bifurcation points in this run are in accordance with the statement of Proposition 2.2.

Beyond the PD point the dynamics of (2.1) is a stable 2-cycle. Cl_MatContM allows to switch to the continuation of this 2-cycle. It loses stability at a supercritical PD point: label = PD, x = (0.849938 0.000000 3.449490) normal form coefficient of PD = 4.062566e+002

A stable 4-cycle is born when . It loses stability at a supercritical PD point: label = PD, x = (0.884050 0.000000 3.544090) normal form coefficient of PD = 1.617268e+004

Thus, when a stable 8-cycle emerges. A stable 8-cycle is given by , where where , , and . This 8-cycle is represented in Figure 2. In this situation there are no predators and the number of prey repeats itself every 8 time spans. Switchings at PD points of the second and fourth iterates are given in Figure 3.

Figure 3: Curves of fixed points of the first, second and fourth iterates.
Figure 4: A stable closed invariant curve for , , and .

We note that for the map (2.1) is a logistic map. In fact, in this case the predator becomes extinct and the prey undergoes the period-doubling bifurcation to chaos when further increasing the parameter ; see [19].

Continuation of starting from the same parameter values as in Run 1, with as free parameter, leads to: label = BP, x = (0.500000 0.000000 2.000000)

The appearance of a branch point is consistent with Proposition 2.2 part (ii) which states that bifurcates to when .

##### 3.2. Numerical Bifurcation of

We now consider which is in the stable region for the parameter values , , and (stability follows from Proposition 2.3 part (i)). We do a numerical continuation of with control parameter , and call this Run 2. The output of Run 2 is given by: label = PD, x = (0.619048 1.666667 1.615385) normal form coefficient of PD = −5.905026e−001 label = NS, x = (0.357143 6.250000 2.800000) normal form coefficient of NS = −6.971545e−002

is stable when . It loses stability via a supercritical Neimark-Sacker (NS) point when , which is consistent with Proposition 2.4 part (ii) (). It also loses stability through a subcritical PD point when , which is consistent with Propositions 2.5 and 2.4 part (i) since .

The dynamics of the system prior to the PD point consists of an unstable 2-cycle that coexists with a stable fixed point. Beyond the NS point the dynamics of the system consists of a stable closed invariant which coexists with unstable fixed points of (2.1). For , a stable closed invariant curve is created around the unstable fixed point (see Figure 4). Now, we compute the period doubling curve, with and free, by starting from the PD point detected in Run 2. We call this Run 3. label = LPPD, x = (0.666667 0.00000 3 1.500000) Normal form coefficient for LPPD:[a/e, be] = −3.333334e−001, 4.541342e−006 label = R2, x = (0.444444 20.000000 9.000000 2.250000) Normal form coefficient for R2:[c, d] = 1.197634e−001, −7.185806e−001

Two codim-2 bifurcation points are detected on the flip curve, namely, a fold-flip LPPD and a resonance 2 bifurcation R2. We note that the LPPD point is in reality a branch-flip (BPPD) point, which by the software is detected as an LPPD point since BPPD points are ungeneric on a curve of PD points. This curve is shown in Figure 5 (left curve). Now, we compute the NS curve, with and free parameters, by starting from the NS point of Run 2. We call this Run 4. label = R4, x = (0.400000 10.000000 5.000000 2.500000 −0.000000) Normal form coefficient of R4: A = −1.240347e+000 + 6.201737e−001 i label = R3, x = (0.428571 15.000000 7.000000 2.333333 −0.500000) Normal form coefficient of R3: Re(c_1) = −5.000000e−001 label = R2, x = (0.444444 20.000000 9.000000 2.250000 −1.000000) Normal form coefficient of R2: [c, d] = 1.197646e−001, −7.185574e−001

Figure 5: Flip and Neimark-Sacker bifurcation curves starting from points in Run 2. We note that the LPPD point is in reality a BPPD point.
Figure 6: A stable 4-cycle near an R4 point for , , and .

The computed curve of NS points is also shown in Figure 5 (right curve). The codim-2 bifurcations that are computed along the Neimark-Sacker curve are a resonance 1 : 2 (R2), resonance 1 : 3 (R3), and a resonance 1 : 4 (R4) point. In addition to the coordinates of the bifurcation point, parameter values and the real part of the Neimark-Sacker multiplier at the bifurcation point are output. We note that the PD curve crosses over the NS curve at a resonance 1 : 2 (R2) point, which is consistent with Proposition 2.4 part (v). It also hits the BP curve at a branch-flip point BPPD which is consistent with Proposition 2.4 part (iv).

##### 3.3. Orbits of Period 4, 8, 16 and 32

The normal form coefficient of the R4 point in Run 4 satisfies , hence two cycles of period 4 of the map are born. A stable 4-cycle for , , and is given by: , where . We present this cycle in Figure 6. In order to compute the stability region of this 4-cycle, we compute the twofold curves of the fourth iterate rooted at the R4 point. These curves exist since , see [9], and switching from an R4 point to the fold curves of the fourth iterate is supported by Cl_MatContM. The stable fixed points of the fourth iterate exist in the wedge between the twofold curves. The output of this continuation, Run 5, is given below. The curves are shown in Figure 7. label = LPPD, x = (0.077546 7.084151 3.948550 3.637018) Normal form coefficient for LPPD:[a/e, be] = −1.638977e−001, −1.231773e+001 label = LPPD, x = (0.540029 8.419070 5.511767 2.615029) Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002

Figure 7: Twofold curves (LP4) emanate from an R4 point on an NS curve.

We can further compute the stability boundaries of the 4-cycle. This region is bounded by the two just computed limit point curves and a period doubling curve of the fourth iterate rooted at the detected LPPD points on the branches of LP4 curves. Continuation of the flip curve of the fourth iterate emanated at the LPPD of Run 5 is given below. We call this Run 6.

label = LPPD, x = (0.862908 2.257241 3.948550 3.637018) Normal form coefficient for LPPD:[a/e, be] = −1.638977e−001, −9.218471e+001 label = LPPD, x = (0.540029 8.419070 5.511767 2.615029) Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002 label = GPD, x = (0.541277 8.238085 5.478361 2.609243) Normal form coefficient of GPD = 6.483721e−001 label = R2, x = (0.582936 7.330070 5.625162 2.664169) Normal form coefficient for R2:[c, d] = 6.130829e−002, −3.509615e+000

By superposing the flip curve on Figure 7, we obtain Figure 8. We further compute a curve of fixed points of the fourth iterate starting from the 4-cycle presented in Figure 6, with control parameter . The output of this continuation, Run 7, is given below. The curve is presented in Figure 9. label = LP, x = (0.418135 9.911757 2.502419) normal form coefficient of LP = −2.095861e−001 label = PD, x = (0.554614 8.335785 2.864961) normal form coefficient of PD = 1.909178e−001 label = NS, x = (0.546431 9.030507 2.608544) normal form coefficient of NS = −2.177981e+000

Figure 8: Twofold curves (LP4) of the fourth iterate emanate from an R4 point, and a flip curve of the fourth iterate (PD4) is rooted at an LPPD point.
Figure 9: Curve of fixed points of the fourth iterate starting from the 4-cycle .

The 4-cycle remains stable when . Now, we compute an NS-curve starting from the computed NS point in Run 8 and call this Run 9. This curve is superposed on Figure 8 and is depicted in Figure 10. label = R4, x = (0.594620 7.423051 4.615965 2.666797 0.000000)Normal form coefficient of R4: A = −4.453277e+000 + −1.491255e+000 i label = R3, x = (0.529139 9.928228 5.218594 2.614545 −0.500000)Normal form coefficient of R3: Re(c_1) = −6.065820e+000 label = R2, x = (0.513005 11.383897 5.625162 2.664169 −1.000000)Normal form coefficient of R2: [c, d] = 1.004097e−001, −5.657220e+000

Figure 10: Stability region of the 4-cycle . The boundaries consist of NS, LP4, PD4, and NS4 curves.
Figure 11: A stable 16-cycle around the R4 point of Run 9 for , , and .

The stability region of the 4-cycle is bounded by the LP4, PD4, and NS4 curves (see Figure 10).

Now, we consider the R4 point computed in Run 9. Since the corresponding normal form coefficient satisfies , two cycles of period 16 of the map are born. A stable 16-cycle for , , and is given by , where . We present this cycle in Figure 11.

In order to compute the stability region of this 16-cycle, we compute twofold curves of the sixteenth iterate rooted at the R4 point. These curves exist since . The stable fixed points of the sixteenth iterate exist in the wedge between the twofold curves. The output of this continuation, Run 10, is given below, and the fold curves are shown in Figure 12. label = LPPD, x = (0.557836 8.068783 4.651131 2.670940) Normal form coefficient for LPPD:[a/e, be] = 1.235611e+001, 3.799853e−001 First Lyapunov coefficient for second iterate = 3.799853e−001 label = LPPD, x = (0.519225 8.927562 4.765661 2.658381) Normal form coefficient for LPPD:[a/e, be] = 9.129886e+000, 5.166118e+000 First Lyapunov coefficient for second iterate = 5.166118e+000 label = LPPD, x = (0.469475 10.678780 5.153847 2.636682) Normal form coefficient for LPPD:[a/e, be] = 1.298255e−001, −4.388341e+000 label = LPPD, x = (0.580415 9.268192 5.370935 2.641181) Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 6.631508e+001 First Lyapunov coefficient for second iterate = 6.631508e+001 label = LPPD, x = (0.578148 9.853334 5.628124 2.676611) Normal form coefficient for LPPD:[a/e, be] = −9.656910e−001, −2.522009e+004 label = R1, x = (0.577592 9.941486 5.660840 2.682153) normal form coefficient of R1 = −1

Figure 12: Twofold curves (LP16) emanate from an R4 point on the NS4 curve.

We further compute a curve of fixed points of the sixteenth iterate starting from the 16-cycle presented in Figure 11, with control parameter . The output of this continuation, Run 11, is given below. The curve is presented in Figure 13. label = LP, x = (0.598114 7.473352 4.646130)normal form coefficient of LP = 2.965904e+000 label = PD, x = (0.595604 7.701225 4.700640) normal form coefficient of PD = 2.829176e+002

Figure 13: Curve of fixed points of the sixteenth iterate starting from the 16-cycle .
Figure 14: A stable 32-cycle for , , and .

The 16-cycle remains stable when . The dynamics of the system beyond the supercritical PD is a stable 32-cycle which coexists with the unstable fixed points of the sixteenth iterate of the system.

A stable 32-cycle for , , and is given by: , where . We present this cycle in Figure 14.

Now, we compute a PD-curve starting from the computed PD point in Run 10, and call this Run 12. This curve is superposed on Figure 12, and is depicted in Figure 15. The boundaries of the stability region of the 16-cycle consist of the LP16 and PD16 curves, as illustrated in Figure 15.label = LPPD, x = (0.600292 7.504485 4.651131 2.670940) Normal form coefficient for LPPD:[a/e, be] = 1.235612e+001, 1.579190e+003 First Lyapunov coefficient for second iterate = 1.579190e+003 label = LPPD, x = (0.587235 7.998201 4.765661 2.658381) Normal form coefficient for LPPD:[a/e, be] = 9.129867e+000, 9.328278e+002 First Lyapunov coefficient for second iterate = 9.328278e+002 label = GPD, x = (0.529363 9.922888 5.130266 2.635413) Normal form coefficient of GPD = 1.292549e+007 label = LPPD, x = (0.526233 10.046755 5.153847 2.636682) Normal form coefficient for LPPD:[a/e, be] = 1.298254e−001, −7.143714e+002 label = LPPD, x = (0.521191 10.840009 5.370935 2.641181) Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 1.217953e+002 First Lyapunov coefficient for second iterate = 1.217953e+002 label = GPD, x = (0.521992 10.863772 5.384512 2.642770) Normal form coefficient of GPD = −6.188850e+003

Figure 15: Stability region of the 16-cycle . The boundaries consist of LP16 and PD16 curves.
##### 3.4. Orbits of Period 3

Next, we consider the resonance 1 : 3 (R3) point in Run 4. Since its normal form coefficient is negative, the bifurcation picture near the R3 point is qualitatively the same as presented in [9, Figure 9.12]. In particular, there is a region near the R3 point where a stable invariant closed curve coexists with an unstable fixed point. For parameter values close to the R3 point, the map has a saddle cycle of period three.

Furthermore, a Neutral Saddle bifurcation curve of fixed points of the third iterate emanates [9, Chapter  9]. We compute this curve by branch switching at the R3 point. This curve is presented in Figure 16.

Figure 16: Curve of Neutral Saddle bifurcations of the third iterate for .
##### 3.5. Computation of Arnold Tongues

It is well known that near a Neimark-Sacker curve there exists a dense array of resonance tongues, generalizing the isolated tongue of period 4 in Figure 7. The tongues locally form an open and dense set of the parameter plane. There are also quasiperiodic invariant circles in between that correspond to a set of positive measure in the parameter plane.

So far, no numerical methods have been implemented to specifically compute the boundaries of the resonance tongues that are rooted in weakly resonant Neimark-Sacker points (unlike the strong resonant 1 : 4 case). However, since they are limit point curves of fixed points of cycles with known periods, they can be computed relatively easily if the cycles inside the tongue are globally stable (which depends on the criticality of the Neimark-Sacker curve and the noncritical multipliers as well). It is sufficient to find a fixed point of cycles inside the tongue by orbit convergence and to continue it in one free parameter to find a point on the boundary of the Arnold tongue as a limit point of cycles. From this, the boundary curves can be computed by a continuation in two free parameters.

In Figure 17, we present an Arnold Tongue rooted in a weak 2 : 7 resonant Neimark-Sacker point. Its computation started from a stable 7-cycle with , , , , and . We note that the boundary curves contain further bifurcation points.

Figure 17: An Arnold tongue rooted in a weak 2 : 7 resonant Neimark-Sacker point.

#### 4. Numerical Simulation

To reveal the qualitative dynamical behaviours of (2.2) near the NS point in Run 1, we present a complete bifurcation sequence that is observed for different values of . We fix the parameters , and assume that is free.

Figure 18 shows that is a stable attractor for . We note that for the given parameters value, is a stable fixed point consistent with Proposition 2.3 part (iii). The behaviour of (2.2) at , so before the NS point at , is depicted in Figure 19.

Figure 18: An attracting fixed point for system (2.2) for , , and .
Figure 19: Phase portrait for the system (2.2) before the NS point (, , and ).

Figure 20 demonstrates the behaviour of the model after the NS bifurcation, for . From Figures 19 and 20 it turns out that the fixed point loses its stability through an NS bifurcation, when varies from 2.78 to 2.81. The dynamics of (2.2) beyond the NS point is a stable closed invariant curve which coexists with unstable fixed point .

Figure 20: Phase portrait for the system (2.2) for , , and .

As is increased further, however, the phase portrait starts to fold. We see that the circle, after being stretched, shrunk, and folded, creates new phenomena due to the breakdown of the closed curve; see Figure 21 for .

Figure 21: The breakdown of the closed invariant curve of the system (2.2) for , , and .

For increasing , we obtain the multiple invariant closed curves brought about the NS bifurcation point of iterates of (2.2). In these cases, higher bifurcations of the torus occurs as the system moves out of quasiperiodic region by increasing . The dynamics move from one closed curve to another periodically, but the dynamics in each closed curve, may be periodic or quasiperiodic. Figure 22 presents the set of closed curves around the NS bifurcation.

Figure 22: The existence of multiple closed curves of the system (2.2) for , , and .

Moreover, the closed curves may break and lead to multiple fractal tori on which the dynamics of (2.2) are chaotic. Figures 23 and 24 present strange attractors for (2.2) with and , respectively, which exhibit a fractal structure.

Figure 23: Chaotic attractor for the system (2.2) for , , and .
Figure 24: Chaotic attractor for the system (2.2) for , , and .

#### 5. Concluding Remarks

In this paper, we studied a planar map that models a predator-prey interaction with nonoverlapping generations. We derived analytically a complete description of the stability regions of the fixed points of the system, namely, , , and . We showed that the system undergoes branching, period doubling, and Neimark-Sacker bifurcations. In particular, we determined criticality of the flip, and Neimark-Sacker bifurcations for and analytically by direct computation of the normal form coefficients without explicit reduction to the center manifold. To support the analytical results and reveal the further complex behaviour of the system, we employed numerical continuation methods to compute curves of codimension 1 and 2 bifurcation points. In particular, we computed curves of fixed points of different cycles as well as curves of fold, flip and Neimark-Sacker bifurcations of the fourth iterate, and fold and flip bifurcations of the sixteenth iterate. These curves form stability boundaries of different cycles of the system. These tasks were possible by means of the computation of normal form coefficients and branch-switching algorithms. We further used numerical simulation methods to reveal chaotic behaviours and a strange attractor near the Neimark-Sacker bifurcation.

#### Acknowledgments

This paper was greatly improved by the comments of four reviewers who read the paper from different viewpoints. The authors thank them all and thank the editorial staff of DDNS for the fast and efficient handling of the manuscript. R. K. Ghaziani also would like to thank Shahrekord University for the financial support through a research grant.

#### References

1. A. A. Berryman, “The origins and evolution of predator-prey theory,” Ecology, vol. 73, no. 5, pp. 1530–1535, 1992.
2. M. Fan and K. Wang, “Periodic solutions of a discrete time nonautonomous ratio-dependent predator-prey system,” Mathematical and Computer Modelling, vol. 35, no. 9-10, pp. 951–961, 2002.
3. H. Fang and J. D. Cao, “Global existence for positive periodic solutions to a class of predator-prey systems,” Journal of Biomathematics, vol. 15, no. 4, pp. 403–407, 2000.
4. H. I. Freedman, Deterministic Mathematical Models in Population Ecology, vol. 57 of Monographs and Textbooks in Pure and Applied Mathematics, Marcel Dekker, New York, NY, USA, 1980.
5. B. S. Goh, Management and Analysis of Biological Populations, Elsevier, Amsterdam, The Netherlands, 1980.
6. J. D. Murray, Mathematical Biology, vol. 19 of Biomathematics, Springer, Berlin, Germany, 2nd edition, 1993.
7. A. D. Bazykin, Nonlinear Dynamics of Interacting Populations, vol. 11 of World Scientific Series on Nonlinear Science. Series A: Monographs and Treatises, A. I. Khibnik and B. Krauskopf, Eds., World Scientific, River Edge, NJ, USA, 1998.
8. V. I. Arnold, Geometrical Methods in the Theory of Ordinary Differential Equations, vol. 250 of Grundlehren der Mathematischen Wissenschaften, Springer, New York, NY, USA, 1983.
9. Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, vol. 112 of Applied Mathematical Sciences, Springer, New York, NY, USA, 3rd edition, 2004.
10. H. W. Broer, K. Saleh, V. Naudot, and R. Roussarie, “Dynamics of a predator-prey model with non-monotonic response function,” Discrete and Continuous Dynamical Systems. Series A, vol. 18, no. 2-3, pp. 221–251, 2007.
11. E. J. Doedel, R. A. Champneys, T. F. Fairgrieve, Yu. A. Kuznetsov, B. Sandstede, and X. J. Wang, “AUTO2000: Continuation and Bifurcation Sotware for Ordinary Differential Equations (with HomCont), Users’ Guide,” Concordia University, Montreal, Canada 1997–2000, http://indy.cs.concordia.ca.
12. W. Govaerts and Yu. A. Kuznetsov, “Matcont: a Matlab software project for the numerical continuation and bifurcation study of continuous and discrete parameterized dynamical systems,” http://sourceforge.net/.
13. R. P. Agarwal, Difference Equations and Inequalities: Theory, Methods, and Applicationa, vol. 228 of Monographs and Textbooks in Pure and Applied Mathematics, Marcel Dekker, New York, NY, USA, 2nd edition, 2000.
14. W. Govaerts, R. K. Ghaziani, Yu. A. Kuznetsov, and H. G. E. Meijer, “Numerical methods for two-parameter local bifurcation analysis of maps,” SIAM Journal on Scientific Computing, vol. 29, no. 6, pp. 2644–2667, 2007.
15. M. Danca, S. Codreanu, and B. Bakó, “Detailed analysis of a nonlinear prey-predator model,” Journal of Biological Physics, vol. 23, no. 1, pp. 11–20, 1997.
16. K. Murakami, “Stability and bifurcation in a discrete-time predator-prey model,” Journal of Difference Equations and Applications, vol. 13, no. 10, pp. 911–925, 2007.
17. S. Li and W. Zhang, “Bifurcations of a discrete prey-predator model with Holling type II functional response,” Discrete and Continuous Dynamical Systems. Series B, vol. 14, no. 1, pp. 159–176, 2010.
18. E. L. Allgower and K. Georg, Numerical Continuation Methods: An Introduction, vol. 13 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1990.
19. R. L. Kraft, “Chaos, Cantor sets, and hyperbolicity for the logistic maps,” The American Mathematical Monthly, vol. 106, no. 5, pp. 400–408, 1999.