Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
Contents lists available at SciVerse ScienceDirect
Nonlinear Analysis: Real World Applications
journal homepage: www.elsevier.com/locate/nonrwa
Resonance and bifurcation in a discrete-time predator–prey system with
Holling functional response
R. Khoshsiar Ghaziania
, W. Govaertsb
, C. Sonckb,∗
a
Department of Applied Mathematics, Shahrekord University, P.O. Box 115, Shahrekord, Iran
b
Department of Applied Mathematics and Computer Science, Ghent University, Krijgslaan 281-S9, B-9000 Gent, Belgium
a r t i c l e i n f o
Article history:
Received 3 August 2010
Accepted 15 November 2011
Keywords:
Iterated map
Normal form
Hyperbolic response
Stable cycle
a b s t r a c t
We perform a bifurcation analysis of a discrete predator–prey model with Holling
functional response. We summarize stability conditions for the three kinds of fixed points
of the map, further called F1, F2 and F3 and collect complete information on this in a
single scheme. In the case of F2 we also compute the critical normal form coefficient of
the flip bifurcation analytically. We further obtain new information about bifurcations of
the cycles with periods 2, 3, 4, 5, 8 and 16 of the system by numerical computation of the
corresponding curves of fixed points and codim-1 bifurcations, using the software package
MatContM. Numerical computation of the critical normal form coefficients of the codim-
2 bifurcations enables us to determine numerically the bifurcation scenario around these
points as well as possible branch switching to curves of codim-1 points. Using parameterdependent
normal forms, we compute codim-1 bifurcation curves that emanate at codim-2
bifurcation points in order to compute the stability boundaries of cycles with periods 4, 5,
8 and 16.
© 2011 Elsevier Ltd. All rights reserved.
1. Introduction
Interactions of different species may take many forms such as competition, predation, parasitism and mutualism. One
of the most important interactions is the predator–prey relationship. The dynamic relationship between predator and prey
is one of the dominant subjects in mathematical ecology due to its universal existence and importance. For more details
on different types of predator–prey systems we refer to [1] and the references cited therein. How predators respond to
changes in prey availability (functional response) is an issue of particular importance. A functional response specifies the
rate at which prey are consumed, per predator, as a function of the prey density. The type of functional response specified
can greatly affect model predictions, see [2,3].
Three general forms of functional response are commonly used in ecological models namely linear, hyperbolic, and
sigmoidal. The linear functional (Lotka–Volterra) response specifies a directly proportional relationship between the
consumption rate of an individual predator and the density of its prey. Holling [4] extended this to include a cap or limitation
(Holling’s type I), where there is an abrupt upper threshold representing predator satiation. The hyperbolic (respectively,
sigmoid) functional response, most commonly known as Holling’s type II (respectively type III) function, incorporate search
rate and predator handling time to produce a smooth asymptotic curve.
Discrete-time predator–prey models go back at least to [5] where the classical discrete-time Lotka–Volterra model is
introduced. They are further studied by many authors, see [6–15]. For example, Sacker and von Bremen [13] propose a
biological model for the genetic reproductive process and prove that for certain parameter values there exist stable invariant
∗ Corresponding author.
E-mail addresses: Khoshsiar@sci.sku.ac.ir (R. Khoshsiar Ghaziani), Willy.Govaerts@UGent.be (W. Govaerts), Charlotte.Sonck@UGent.be (C. Sonck).
1468-1218/$ – see front matter © 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/j.nonrwa.2011.11.009
1452 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
curves on which the map has a quasi-periodic behavior. Four typical discrete-time ecosystem models under the effects of
periodic forcing have been studied in [14]. In a simple discrete-time predator–prey model with Holling’s type I functional
response, chaotic dynamics can occur [16].
In this paper we consider the general case of a discrete-time predator–prey system with Holling type II response. Our
model is equivalent to model (30) in [12] and therefore to a special case of a model in [17] which itself is a discretization of
a system studied in [18].
In [19,20] the same model is studied and we follow their notation. The main differences with these papers are that (a)
we focus on results that have an immediate ecological meaning, (b) that we consider not only the map but also its iterates
so that periodic orbits come into the picture much more prominently, and (c) that we rely heavily on advanced numerical
methods, namely numerical continuation to obtain results that cannot be obtained analytically.
Point (a) means, for example, that we are interested mainly in stability and instability of fixed points of the map and its
iterates. Contrary to what is done in [19,20], the linear nature of the fixed points is of little interest to us. This leads to a
simple classification of the stability and instability regions of all fixed points of the map that can be summarized in a global
picture (Fig. 2).
For (b) we remark that the only numerical illustrations in [19,20] are orbit simulations. This does not allow to compute
stability boundaries numerically.
The present paper can also be seen as an extension of [21]; the map (2) reduces to the map (2.1) in [21] in the limit case
ϵ = 0 when the Holling functional response reduces to a linear response. Also, the bifurcation diagram in Fig. 2 reduces to
Fig. 1 in [21] in the case where ϵ tends to zero. However, the bifurcation behavior is not identical. We comment on this in
Section 2.3. The numerical computations in Sections 3.2–3.5 are done with respect to ϵ as a free parameter, so they cannot
be compared to any results in [21].
In Section 2 we introduce the model and discuss the stability and bifurcations of its fixed points. We derive analytically
the stability regions of all types of fixed points and their bifurcation behaviors. Moreover, we compute analytically the critical
normal form coefficients in the case of the period doubling bifurcation to prove supercriticality. These results correct those
in [19] and are consistent with those of [20]. In the degenerate case where the functional response reduces to a mass-action
law, they were already obtained in [12,11].
In Section 3 we numerically compute curves of codim-1 bifurcations and the critical normal form coefficients of codim-2
bifurcation points, using the Matlab toolbox MatContM [22,23]. These tools enable us to compute stability boundaries of
different cycles. In particular, we determine the bifurcation scenario of the map near an R4 resonance point, which involves
stable and unstable 4-cycles as well as 8-cycles and 16-cycles. As an example application we compute a region where a
stable 8-cycle coexists with a stable 4-cycle. As another application we compute a branch of neutral saddle 3-cycles. We
furthermore compute an Arnol’d tongue of period 5 and two related curves of codim-1 bifurcations which together delineate
a parameter region where stable period-5 cycles exist.
In Section 4 we summarize our results.
2. Holling type II predator–prey system, existence and stability of its fixed points
The Lotka–Volterra predator–prey system (see e.g. [24,25,15,26]) is a fundamental population model. More realistic
predator–prey systems were introduced by Holling [4] using the three kinds of functional responses for different species to
model the phenomena of predation. We first mention the continuous-time predator–prey model studied in [27] and later
studied in a discretized version in [19]:
˙x(t) = α0x(1 − x) − α
mxy
1 + ϵx
,
˙y(t) =

mx
1 + ϵx
− β

y,
(1)
where α0, m, α, β and ϵ > 0 are parameters and x(t), y(t) represent the densities of the prey and the predator, respectively.
α0 is the intrinsic growth rate of the prey, m is a mass-action law constant, ϵ is a limitation parameter of the growth of the
predator population for increasing prey density, ax(1 − x) is a logistic function and
mxy
1+ϵx
is the Holling type II functional
response, β and α denote the death rate of the predator and conversion, respectively.
In the case of a predator–prey system with non-overlapping generations this can be replaced by a discrete system
F :

x
y

→



ax(1 − x) −
bxy
1 + ϵx
dxy
1 + ϵx


 , (2)
where a, b, d are nonnegative parameters and ϵ > 0. This model was studied in [19,20]. We note that the parameter b can
be absorbed by rescaling y. Therefore this is, from a bifurcation point of view, a three-parameter problem and we will indeed
see that b does not appear in any bifurcation equations.
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1453
The bifurcation analysis of (2) naturally starts with fixed points. The fixed points of (2) are the solutions (x∗
, y∗
) to
ax∗
(1 − x∗
) −
bx∗
y∗
1 + ϵx∗
= x∗
,
dx∗
y∗
1 + ϵx∗
= y∗
.
The origin F1 = (0, 0) is always a fixed point of (2). Two further fixed points of the system are given by F2 = (a−1
a
, 0) which
is biologically possible for a ≥ 1 and
F3 =

1
d − ϵ
,
d
d − ϵ

a
b

1 −
1
d − ϵ

−
1
b

. (3)
We note that F3 is biologically possible if its coordinates are nonnegative, i.e.,
a > 1, d ≥ ϵ +
a
a − 1
. (4)
We start the local bifurcation analysis of the map (2) by linearization of F around each of its fixed points. The Jacobian matrix
J(x, y) is given by:
J(x, y) =




a(1 − 2x) −
by
(1 + ϵx)2
−
bx
1 + ϵx
dy
(1 + ϵx)2
dx
1 + ϵx



 . (5)
The characteristic equation of J(x, y) is given by
λ2
− tr(J)λ + det(J) = 0, (6)
where tr(J) = a(1 − 2x) − by
(1+ϵx)2 + dx
1+ϵx
and det(J) = adx(1−2x)
1+ϵx
.
2.1. Stability of F1
Proposition 1. The fixed point F1 is asymptotically stable for 0 ≤ a < 1. It loses stability via branching for a = 1 and there
bifurcates to F2.
Proof. Eigenvalues of the Jacobian at F1 are a and 0. So F1 is stable if a < 1 and loses stability at a = 1. It remains to show
that F1 bifurcates to F2 at a = 1. To do this we consider the matrix (Fx − I|Fa), evaluated at F1:

a − 1 0 0
0 −1 0

. (7)
When a = 1, this matrix is clearly rank deficient. We choose vectors φ1 and φ2 which form a basis for the null space of J(F1)
and a vector ψ that spans the null space of J(F1)T
. A possible choice is:
φ1 =

1
√
2
, 0,
1
√
2
T
, φ2 = (1, 0, 0)T
, ψ = (1, 0)T
.
Now we consider the algebraic branching equation (ABE), see [22],
c11α2
+ 2c12αβ + c22β2
= 0, (8)
where cjk = ⟨ψ, F0
YY φjφk⟩ for j, k = 1, 2. Here the 2 × 3 × 3 tensor F0
YY of second derivatives of F with respect to the three
variables x, y, a taken at (0, 0) is given by:
F0
YY (:, :, 1) =

−2a −b 1
0 d 0

, (9)
F0
YY (:, :, 2) =

−b 0 0
d 0 0

, (10)
F0
YY (:, :, 3) =

1 0 0
0 0 0

. (11)
We obtain c11 = −2a+2
2
, c12 = −2a+1√
2
, c22 = −2a. For a = 1, we obtain c11 = 0. So the discriminant of (8) is ∆ =
c2
12 − c11c22 = c2
12 = 1
2
> 0.
If Y(s) is any branch of fixed points of (2) with Y(s0) = (0, 0, 1), then its derivative Ys(s0) can be written as Ys(s0) =
αφ1 + βφ2 where α and β are scaled roots of (8). We find
Y1s =

1
√
2
, 0,
1
√
2
T
, Y2s = (0, 0, 1)T
.
1454 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
The branch X(s) = (0, 0, s) of F1-equilibria has unit tangent vector Xs = (0, 0, 1)T
, i.e. Y2s. Differentiation of the branch of
F2-equilibria w.r.t. the parameter a leads to the vector ( 1
a2 , 0, 1)T
. When a = 1, the scaled tangent vector to this branch is

1√
2
, 0, 1√
2
T
, i.e. Y1s.
2.2. Stability of F2
The Jacobian matrix of (2) at F2 is given by
J(F2) =




2 − a
b(1 − a)
a + ϵ(a − 1)
0
d(a − 1)
a + ϵ(a − 1)



 . (12)
Proposition 2. The fixed point F2 is asymptotically stable iff a ∈]1, 3[ and d < ϵ + a
a−1
. Moreover, it loses stability:
(i) via branching for a = 1 and there bifurcates to F1.
(ii) via branching for d = ϵ + a
a−1
and there bifurcates to F3 if 1 < a < 3.
(iii) via a supercritical flip for a = 3 if d < ϵ + 3
2
.
Proof. The eigenvalues of J(F2) are λ1 = 2 − a and λ2 = d(a−1)
a+ϵ(a−1)
. The fixed point F2 is asymptotically stable iff |λ1| < 1
and |λ2| < 1, i.e. iff a ∈]1, 3[ and d < ϵ + a
a−1
.
Boundary points of the stability region must satisfy one of three conditions: a = 1, d = ϵ + a
a−1
, or a = 3.
In the first case the conditions d < ϵ + a
a−1
and a < 3 are satisfied for nearby values a > 1, hence this is a real stability
boundary. In Proposition 1 we proved that this is a branch point and the new branch consists of F1 points.
In the second case this is a stability boundary only if 1 < a < 3. The Jacobian (12) then has an eigenvalue +1 and it is
checked easily that these boundary points are also F3 points.
In the third case, this is a stability boundary only if d < ϵ + 3
2
. In this case λ1 = −1 which means that F2 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 bPD,
bPD =
1
6

p, C(q, q, q) + 3B(q, (I − A)−1
B(q, q))

, (13)
derived by center manifold reduction is positive, see [28, Ch. 8], and [22]. Here A = J(F2), and B(., .), C(., ., .) are the second
and third order multilinear forms respectively, and p and q are the left and right eigenvectors of A for the eigenvalue −1,
respectively. These vectors are normalized by ⟨p, q⟩ = 1, ⟨q, q⟩ = 1, where ⟨., .⟩ is the standard scalar product in R2
. We
obtain:
q =

q1
q2

=

1
0

, (14)
and
p =

p1
p2

=


1
2b
3 + 2(d + ϵ)

 . (15)
The components of the multilinear form B(q, q) are given by:
[B(q, q)]1 =
2−
j,k=1
∂2

ax(1 − x) − bxy
1+ϵx

∂xj∂xk
qjqk = −2a = −6, (16)
[B(q, q)]2 =
2−
j,k=1
∂2
 dxy
1+ϵx

∂xj∂xk
qjqk = −
2dyϵ
(1 + ϵx)3
= 0, (17)
where the state variable vector is for ease of notation generically denoted by (x1, x2)T
instead of (x, y)T
.
Let ζ = (I − A)−1
B(q, q), then we have ζ =

−3
0

and find
[B (q, ζ)]1 = −2a(−3) = 18, [B (q, ζ)]2 =
6dyϵ
(1 + ϵx)3
= 0. (18)
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1455
Fig. 1. SF2 is the stability region of F2 for ϵ = 0.5.
The third order multilinear form C(q, q, q) is given by
[C(q, q, q)]1 =
2−
j,k,l=1
∂3

ax(1 − x) − bxy
1+ϵx

∂xj∂xk∂xl
qjqkql =
−6byϵ2
(1 + ϵx)4
= 0, (19)
[C(q, q, q)]2 =
2−
j,k,l=1
∂3
 dxy
1+ϵx

∂xj∂xk∂xl
qjqkql =
6dyϵ2
(1 + ϵx)4
= 0. (20)
The critical normal form coefficient bPD is given by
bPD =
1
6
pT

54
0

= 9, (21)
which is clearly positive. This completes the proof of supercriticality of the flip point at F2.
The stability region SF2 of F2, as obtained in Proposition 2, is shown in Fig. 1 for ϵ = 0.5. We note that for d < ϵ + 3
2
and
a close to but larger than 3 the map can behave in a stable alternating way. We note that the supercriticality of the flip
bifurcation in Proposition 2 was also obtained in [20] though our proof is different.
2.3. Stability of F3
To study the stability of F3 we use the Jury’s criteria, see [29, Section A2.1]. Let F(λ) = λ2
− tr(J(F3))λ + det(J(F3)) be
the characteristic equation of J(F3). Hence we have F(λ) = (λ − λ1)(λ − λ2) where λ1, λ2 are the eigenvalues of J(F3).
According to the Jury’s criteria F3 is asymptotically stable if the following conditions hold:
F(−1) = 1 + tr(J(F3)) + det(J(F3)) > 0,
F(1) = 1 − tr(J(F3)) + det(J(F3)) > 0,
1 − det(J(F3)) > 0.
(22)
At F3 we have:
J(F3) =




a

1 −
2
d − ϵ

−
d − ϵ
d

a

1 −
1
d − ϵ

− 1

−b
d
d − ϵ
b

a

1 −
1
d − ϵ

− 1

1



 . (23)
We note that:
tr(J(F3)) = 1 + a

1 −
2
d − ϵ

−
d − ϵ
d

a

1 −
1
d − ϵ

− 1

,
and
det(J(F3)) = a

1 −
2
d − ϵ

,
are independent of b.
1456 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
Proposition 3. F3 is asympotically stable iff one of the following mutually exclusive conditions holds:
(i) ϵ + 3
2
< d < 1
8

4ϵ + 9 +
√
16ϵ2 + 56ϵ + 81

and d−ϵ
d−ϵ−1
< a < (d−ϵ)(ϵ−3d)
2d(d−ϵ−2)+(d−ϵ)(−d+ϵ+1)
,
(ii) d = 1
8

4ϵ + 9 +
√
16ϵ2 + 56ϵ + 81

and d−ϵ
d−ϵ−1
< a < d−ϵ
d−ϵ−2
= (d−ϵ)(ϵ−3d)
2d(d−ϵ−2)+(d−ϵ)(−d+ϵ+1)
,
(iii) d > 1
8

4ϵ + 9 +
√
16ϵ2 + 56ϵ + 81

and d−ϵ
d−ϵ−1
< a < d−ϵ
d−ϵ−2
.
Proof. The criterion F(1) > 0 is easily seen to be equivalent to the condition
a >
d − ϵ
d − ϵ − 1
, (24)
or equivalently,
d > ϵ +
a
a − 1
, (25)
i.e. a slightly stronger version of the second condition in (4).
Next, the criterion det(J(F3)) < 1 is easily seen to be equivalent to
a <
d − ϵ
d − ϵ − 2
for all d > ϵ + 2. (26)
The criterion F(−1) > 0 translates as
a <
(d − ϵ)(3d − ϵ)
(d − ϵ)(d − ϵ − 1) − 2d(d − ϵ − 2)
(27)
for all d, ϵ that satisfy
(d − ϵ)(d − ϵ − 1) − 2d(d − ϵ − 2) > 0. (28)
The latter equation is easily seen to be equivalent to
d <
3 +
√
9 + 4ϵ2 + 4ϵ
2
. (29)
It is also easy to see that for a given ϵ > 0 the only values of a, d for which the inequalities in (24) and (27) are both equalities
is found for d = ϵ + 3
2
, a = 3.
Similarly, for a given ϵ > 0 there is a unique pair a, d for which the inequalities in (26) and (27) are both equalities. For
this point we have
d =
1
8

4ϵ + 9 +

16ϵ2 + 56ϵ + 81

, a =
−4ϵ +
√
16ϵ2 + 56ϵ + 81 + 9
−4ϵ +
√
16ϵ2 + 56ϵ + 81 − 7
. (30)
The stability regions (i) and (iii) of F3 obtained in Proposition 3 are depicted as S
(i)
F3 and S
(ii)
F3 , respectively, in Fig. 2 for
ϵ = 0.5. The region (ii) is the open interval on the common boundary of S
(i)
F3 and S
(ii)
F3 . We note that these regions are
qualitatively similar for all values ϵ > 0 and so we have a complete description of the stability region of F3 for all parameter
combinations.
2.4. Notes
Proposition 4 in [19] is somewhat similar to our Proposition 3 but is incomplete and needs several corrections. We give
a few examples.
• Consider fixed values d = 3.5, ϵ = 1. For these we have
d >
1
8

4ϵ + 9 +

16ϵ2 + 56ϵ + 81

≈ 3.17116.
By Proposition 3 F3 is stable when 1.666667 < a < 5.
The stability condition in Proposition 4, part (i) in [19] includes the condition
a >
(3d − ϵ)(d − ϵ)
(d − ϵ)(2d − ϵ) − 2d(d − ϵ − 2)
≈ 2.065217,
so this condition is not necessary.
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1457
Fig. 2. Stability regions in (d, a)-space for ϵ = 0.5. SF1, SF2 are the stability regions of F1, F2 respectively. The stability region SF3 of F3 is the union
of S
(i)
F3 and S
(ii)
F3 (which correspond to (i) and (iii) in Proposition 3, respectively), and the open interval that separates them (and corresponds to (ii)
in Proposition 3). a = (d−ϵ)(3d−ϵ)
(d−ϵ)(d−ϵ−1)−2d(d−ϵ−2)
indicates the PD curve, a = d−ϵ
d−ϵ−2
determines the NS curve, a = d−ϵ
d−ϵ−1
indicates the BP curve and
a = ϵ
2
+ 1
8
√
9 + 16ϵ2 + 56ϵ + 81 indicates the line that separates the regions S
(i)
F3 and S
(ii)
F3 .
• For the parameter values d = 3.7, ϵ = 2,
1
8

4ϵ + 9 +

16ϵ2 + 56ϵ + 81

≈ 4.1289
and so the first condition in part (i) of Proposition 3 is satisfied. The second condition leads to (approximately)
2.42857 < a < 4.53665.
But the condition a < d−ϵ
d−ϵ−2
= −5.666667 in Proposition 4, part (i) of [19] can never be satisfied for a > 0. Hence that
condition is not a necessary one.
• For d = 3.5, ϵ = 1, it is claimed in part (iii) of Proposition 4 in [19] that F3 is non-hyperbolic if
a =
(3d − ϵ)(d − ϵ)
(d − ϵ)(2d − ϵ) − 2d(d − ϵ − 2)
≈ 2.06522.
Part (i) of Proposition 3 in fact proves that this is a stable fixed point.
• For d = 3.5, ϵ = 1, it is claimed in part (iv) of Proposition 4, in [19] that F3 is a saddle if a < (3d−ϵ)(d−ϵ)
(d−ϵ)(2d−ϵ)−2d(d−ϵ−2)
≈
2.06522. Part (i) of Proposition 3 in fact proves that this is a stable fixed point for values of a slightly below 2.06522.
We further note that in the case ϵ = 0 our model (2) reduces to model (7) in [12] under a rescaling of the x-variable and
introduction of the parameters c, r with a = r + 1, b = c, d = (r + 1)c/r. By straightforward calculations one finds that
the new equations for the PD, BP and NS curve then must be replaced by c = 3r/(r + 4), c = 1, and c = 2, respectively.
Hence we exactly reproduce the results summarized in Fig. 1 of [12]. We note, however, that in [12] the subcriticality of the
PD bifurcations is proved (for ϵ = 0).
Model (30) in [12] uses the parameters c, r, γ . It is equivalent to (2) under a rescaling of the x− variable and by setting
a = r + 1, b = c/γ , ϵ = (r + 1)/(rγ ), d = c(r + 1)/(rγ ). In [12] the stability region of F3 is depicted in (r, c)-space in the
case γ = 1 and it is mentioned that subcritical flip points were found in this case.
The case ϵ = 0 also reduces to model (1) in [11] (in terms of the parameters α, β where a = α, b = 1, d = 1/β). We
note, however, that in [11] the supercriticality of the NS bifurcations is proved (for ϵ = 0).
2.5. Bifurcations of F3
Proposition 4. Suppose that the conditions (4) hold. Then F3 loses stability:
(i) via a flip point when a = (d−ϵ)(3d−ϵ)
(d−ϵ)(d−ϵ−1)−2d(d−ϵ−2)
and ϵ + 3
2
< d < 1
8

4ϵ + 9 +
√
16ϵ2 + 56ϵ + 81

,
(ii) via a Neimark–Sacker bifurcation when a = d−ϵ
d−ϵ−2
and d > 1
8

4ϵ + 9 +
√
16ϵ2 + 56ϵ + 81

,
(iii) via branching when a = d−ϵ
d−ϵ−1
and d > ϵ + 3
2
, where it bifurcates to F2,
(iv) via a fold-flip (LPPD) point when d = ϵ + 3
2
and a = 3,
(v) via a resonance 1:2 point when d = 1
8

4ϵ + 9 +
√
16ϵ2 + 56ϵ + 81

and a =
−4ϵ+
√
16ϵ2+56ϵ+81+9
−4ϵ+
√
16ϵ2+56ϵ+81−7
.
1458 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
Proof. By Proposition 3 (see also Fig. 2) the stability boundary of F3 consists of parts of three curves, namely
1. Curve 1: a = (d−ϵ)(3d−ϵ)
(d−ϵ)(d−ϵ−1)−2d(d−ϵ−2)
,
2. Curve 2: a = d−ϵ
d−ϵ−2
,
3. Curve 3: a = d−ϵ
d−ϵ−2
.
The points of Curve 1 which are on the stability boundary of F3 satisfy F(−1) = 0, i.e. they have an eigenvalue −1. The
points of Curve 2 which are on the stability boundary satisfy det(J(F3)) = 1, i.e. they have two eigenvalues with product 1.
The points of Curve 3 which are on the stability boundary satisfy F(1) = 0, i.e. they have an eigenvalue +1.
Combining this with Proposition 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 the Proposition, respectively.
Next, Curves 1 and 3 have an intersection point which has eigenvalues +1 and −1. This is the LPPD point in part (iv) of
the Proposition (but note that it is degenerate, in fact a BPPD point).
Finally, Curves 1 and 2 also have an intersection point. In this point both eigenvalues are −1. This is the 1 : 2(R2) in part
(v) of the the Proposition.
We remark that our numerical evidence indicates that the flip and Neimark–Sacker bifurcations in Proposition 4 are
sub- and supercritical, respectively in all cases considered in this paper. This is based on the numerical computation of the
normal form coefficients of these bifurcations (see [28, Ch. 8] and [22]). However, a detailed study in [20, Theorems 4.1 and
5.1] suggests that other cases might be possible for certain combinations of parameter values. Resonances are also possible,
cf. the discussion at the end of Section 5 in [20].
We further note that the bifurcation diagram in Fig. 2 reduces to Fig. 1 in [21] in the case where ϵ tends to zero. For this
limit case it was proved in [21] that the flip and Neimark–Sacker bifurcations of F3 are sub- and supercritical, respectively.
We note that for ϵ = 0.5 and (d, a) values close to but to the right of the right boundary of S
(ii)
F3 in Fig. 2 there exist stable
invariant curves on which the map behaves in a quasi-periodic way.
We will not investigate this further but concentrate on the numerical study of cycles and their stability since this leads
to finding parameter regions where stable periodic behavior with high periods can be expected.
3. Numerical bifurcation analysis of F2 and F3
In this section we perform a numerical bifurcation analysis by using the MATLAB package MatContM, see [22,23]. The
bifurcation analysis is based on continuation methods, tracing out the solution manifolds of fixed points while some of the
parameters of the map vary, see [30]. Test runs for these computations will be made available via [23].
3.1. Numerical analysis of F2
We continue F2 = (0.54545454, 0) starting with a = 2.2, b = 3, ϵ = 0.2, d = 0.5 with a free. We see that F2 is stable
when 1 < a < 3. It loses stability via a supercritical period doubling point when a = 3, and via a branch point when a crosses
1. These results are consistent with Proposition 2 since d < ϵ + 3
2
≤ d+ a
a−1
for all a ∈]1, 3[. 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 F2, and the last entry of x is the value of the
free parameter a at the corresponding bifurcation point. We note that the normal form coefficient of the PD point is 9,
confirming (21). The curve computed in Run 1 is presented in Fig. 3. Beyond the PD point the dynamics of (2) is a stable
2-cycle. MatContM allows to switch to the continuation of this 2-cycle. It loses stability at a supercritical PD point for
a = 3.449490. A stable 4-cycle is born when a > 3.449490. An instance is given by C4 =

X4
1 , X4
2 , X4
3 , X4
4

where X4
1 =
(0.87867696, 0), X4
2 = (0.37483245, 0), X4
3 = (0.82394508, 0), X4
4 = (0.51004805, 0). The corresponding parameter
values are a = 3.516128, b = 3, ϵ = 0.2, d = 0.5. We note that the y− coordinate equals zero in all four points of the
cycle. The 4-cycle loses stability via a supercritical PD point for a = 3.544090. Thus, when a > 3.544090 a stable 8-cycle
emerges. This 8-cycle loses stability at another supercritical PD point for a = 3.564407. In fact a cascade of period doublings
appears if we further increase a. We note that for F2 the map (2) is a logistic map, which is well know to have chaotic behavior
through a cascade of period doubling points, see [31].
Continuation of F2 starting from the same parameter values as in Run 1, with d as free parameter, leads to the discovery
of a branch point for d = 2.033333. This is consistent with Proposition 2 part (ii) which states that F2 bifurcates to F3 when
d = ϵ + a
a−1
= 0.2 + 2.2
1.2
= 2.033333.
3.2. Numerical analysis of F3
We now consider F3 = (0.4, 0.681333) which is in the stable region for the parameter values a = 4.1, b = 3, ϵ = 1, d =
3.5 (stability follows from Proposition 3 part (i)). We do a numerical continuation of F3 with control parameter ϵ, and call
this Run 2:
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1459
Fig. 3. Continuation of F2 in (x, a)-space.
Fig. 4. The stable closed invariant curve for a = 4.1, b = 3, ϵ = 0.84, d = 3.5.
label = NS, x = ( 0.378049 0.683638 0.854839 )
normal form coefficient of NS = -9.079782e+000
label = PD, x = ( 0.603958 0.439521 1.844256 )
normal form coefficient of PD = -7.976577e+000
label = BP, x = ( 0.756098 0.000000 2.177419)
F3 is stable when 0.854839 < ϵ < 1.844256. It loses stability via a supercritical Neimark–Sacker (NS) point when ϵ =
0.854839, which is consistent with Proposition 4 part (ii) ( d−ϵ
d−ϵ−2
= 4.1 = a). It loses stability through a subcritical PD point
when ϵ = 1.844256, which is consistent with Proposition 4 part (i) since
(d−ϵ)(3d−ϵ)
(d−ϵ)(d−ϵ−1)−2d(d−ϵ−2)
= 4.1 = a.
The dynamics of the system prior to the PD point consists of an unstable 2-cycle that coexists with a stable fixed point.
For ϵ < 0.854839, a stable closed invariant curve with quasi-periodic behavior is created around the unstable fixed point
F3. Such a curve is shown in Fig. 4. Now we compute the period doubling curve, with a and ϵ free, by starting from the PD
point detected in Run 2. We call this Run 3.
label = LPPD, x = ( 0.666667 0.000000 3.000000 2.000000 )
Normal form coefficient for LPPD:[a/e, be]= -1.200000e-001, -3.467480e-006,
label = R2, x = ( 0.464286 3.520833 14.000000 1.346154 )
Normal form coefficient for R2:[c, d]= 3.479148e+001, -1.346125e+002
Two codim-2 bifurcation points are detected on the flip curve, namely a fold-flip LPPD and a resonance 2 bifurcation R2,
see Fig. 5 (right curve). We note that the LPPD point is degenerate in the sense that it is really a BPPD point. The parameter
1460 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
Fig. 5. Flip and Neimark–Sacker bifurcation curves starting from points in Run 2.
values indeed satisfy the equation a = (d − ϵ)/(d − ϵ − 1) of the branch point curve, cf. Proposition 4. In general, a BPPD
point is ungeneric on a PD curve and therefore is not detected by the software. It follows that the normal form coefficients
given for the pretended LPPD point are meaningless.
Now we compute the NS curve, with a and ϵ free parameters, by starting from the NS point of Run 2. We call this Run 4.
label = R4, x = ( 0.428571 1.500000 7.000000 1.166667 -0.000000 )
Normal form coefficient of R4: A = -1.260354e+000 + 7.511202e-001 i
label = R3, x = ( 0.452381 2.506944 10.500000 1.289474 -0.500000 )
Normal form coefficient of R3: Re(c_1) = -5.439941e-001
label = R2, x = ( 0.464286 3.520833 13.999998 1.346154 -1.000000 )
Normal form coefficient of R2: [c, d] = 3.479148e+001, -1.346125e+002
label = CH, x = ( 0.122849 0.023342 1.325730 -4.640055 0.810610 )
Normal form coefficient of CH = 9.285903e+001
The computed curve of NS points is also shown in Fig. 5 (left curve). We note that the PD and NS curves intersect in an R2
point. The codim-2 bifurcations that are computed along the Neimark–Sacker curve are a resonance 1:2 (R2), resonance 1:3
(R3), resonance 1:4 (R4) and a Chenciner bifurcation (CH)(not in the physically relevant region). 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.
3.3. Orbits of period 4, 8 and 16
In this subsection we will describe parameter regions where predator and prey can coexist in a stable way and reproduce
their densities every fourth, eighth or sixteenth years. The regions with 4-cycle and 8-cycle stability overlap, so that
bistability occurs. The normal form coefficient A of the R4 point in Run 4 satisfies |A| > 1, hence two cycles of period 4
of the map are born. A stable 4-cycle for a = 6.401280, b = 3, ϵ = 1.093117 and d = 3.5 is given by: C4 = {X1, X2, X3, X4}
where X1 = (0.347307, 1.272391). In order to compute the stability region of this 4-cycle, we compute two fold curves of
the fourth iterate rooted at the R4 point. These curves exist since |A| > 1, see [28] and switching from an R4 point to the
fold curves of the fourth iterate is supported by MatContM. The stable fixed points of the fourth iterate exist in the wedge
between the two fold curves. The output of this continuation, Run 5, is given below and the fold curves (denoted by LP4
) are
shown in Fig. 7.
label = LPPD, x = ( 0.341617 1.808057 7.596591 1.084004 )
Normal form coefficient for LPPD:[a/e, be]= 5.904006e-001, -4.511886e+001,
label = LPPD, x = ( 0.085157 0.508424 3.947668 -0.035122 )
Normal form coefficient for LPPD:[a/e, be]= -1.863197e-001, -2.336863e+003,
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.
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1461
Fig. 6. Curve of fixed points of the fourth iterate starting from the 4-cycle C4.
Fig. 7. NS bifurcation curve, two branches of LP4
cycles rooted at the R4 point, a PD4
curve rooted at the LPPD point, two branches of LP8
cycles which are
born tangentially at the GPD points. These curves form stability boundaries of 4-cycles (Ω4
S ∪ Ω4,8
S ) and 8-cycles (Ω8
S ) and bound the bistability region of
4-cycles and 8-cycles (Ω4,8
S ).
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 = GPD, x = ( 0.341591 1.802206 7.570354 1.087021 )
Normal form coefficient of GPD = 8.180449e+005
label = R2, x = ( 0.317814 1.734061 7.225039 1.068919 )
Normal form coefficient for R2:[c, d]= -2.430015e+002, -1.078458e+003
label = GPD, x = ( 0.311884 1.760666 7.295843 1.056324 )
Normal form coefficient of GPD = -3.64862e+004
This curve is depicted in Fig. 7 and indicated by PD4
curve. We further compute a curve of fixed points of the fourth iterate
starting from the 4-cycle C4 with control parameter a. We call this Run 7. The curve is presented in Fig. 6.
The 4-cycle is stable in the wedge between the two LP4
curves, and loses stability when crossing the PD4
curve.
When C4 loses stability at the supercritical PD point corresponding to a = 7.284657, a stable 8-cycle is born which
coexists with an unstable 4-cycle until the second PD point (a = 7.483037) is reached. A stable 8-cycle is given by
C8 =

X8
1 , X8
2 , X8
3 , X8
4 , X8
5 , X8
6 , X8
7 , X8
8

where X8
1 = (0.510422, 1.384081) for a = 7.411918, b = 3, ϵ = 1.093117 and
d = 3.5. We can compute the stability boundaries of C8, by computing two fold curves of the eighth iterate by switching at
the GPD points in Run 6. Again, this is supported by MatcontM. These fold curves emanate tangentially to the fold curve of
the fourth iterate in Run 6. These curves are presented in Fig. 7 and indicated by LP8
curves. The region where C8 is stable is
bounded by the two fold curves of the eighth iterate and the lower part of the flip curve of the fourth iterate (shaded region
indicated by Ω8
S ).
1462 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
Fig. 8. A stable 16-cycle near the R4 point of Run 9 for a = 5.549119, b = 3, ϵ = 0.930419, d = 3.5.
Fig. 9. Two fold curves (LP16
) emanate from an R4 point.
Further we continue the 4-cycle (C4) with control parameter ϵ. The output of this continuation, Run 8, is:
label = NS, x = ( 0.328281 1.215657 1.057947 )
normal form coefficient of NS = -4.089749e+002
label = LP, x = ( 0.355258 1.300254 1.100154 )
normal form coefficient of LP =-7.154929e+000
The 4-cycle remains stable when 1.057947 < ϵ < 1.100154. Now we compute a NS-curve starting from the computed
NS point in Run 8, given below. This curve is depicted in Fig. 7 and indicated by NS4
curve. We call this Run 9.
label = R3, x = ( 0.331234 1.257772 6.541014 1.066957 -0.500000 )
Normal form coefficient of R3: Re(c_1) = -4.885966e+000
label = R2, x = ( 0.332363 1.439768 7.225039 1.068919 -1.000000 )
Normal form coefficient of R2: [c, d] = -1.046630e+002, -4.653504e+002
label = R4, x = ( 0.294522 0.915737 5.526230 0.928728 0.000000 )
Normal form coefficient of R4: A = -4.055960e+000 + -8.600805e-001 i
We note that we have bistability of the 4-cycle (C4) and 8-cycle (C8) in the region indicated by Ω4,8
S in Fig. 7. We note that
the NS curve of the fourth iterate in Run 9 and the PD curve of the fourth iterate of Run 6 intersect in an R2 point. Now we
consider the R4 point computed in Run 9. Since |A| > 1 (A is the corresponding normal form coefficient of the R4 point), two
cycles of period 16 of the map are born. A stable 16-cycle is given in Fig. 8. In order to compute the stability region of this
16-cycle, we compute two fold curves of the sixteenth iterate rooted at the R4 point. These curves exist since |A| > 1. The
stable fixed points of the sixteenth iterate exist in the wedge between the two fold curves. The output of this continuation,
Run 10 is given below. The fold curves are shown in Fig. 9.
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1463
Fig. 10. Curve of Neutral Saddles of the third iterate.
label = LPPD, x = ( 0.294382 0.903595 5.468882 0.911109 )
Normal form coefficient for LPPD:[a/e, be]= 4.600977e+000, -1.127121e+002,
label = R1, x = ( 0.221519 0.406043 4.347547 0.363383 )
normal form coefficient of R1 = -1
label = LPPD, x = ( 0.222150 0.409925 4.354848 0.368810 )
Normal form coefficient for LPPD:[a/e, be]= -5.087744e-001, -1.179672e+009,
label = R1, x = ( 0.181745 0.167415 3.975582 -0.000054 )
normal form coefficient of R1 = -1
label = LPPD, x = ( 0.302337 0.997939 5.668515 0.954589 )
Normal form coefficient for LPPD:[a/e, be]= 3.260385e+001, -1.518085e+001,
label = LPPD, x = ( 0.324821 1.218070 6.191127 1.018418 )
Normal form coefficient for LPPD:[a/e, be]= -8.188768e-001, -3.916020e+005,
label = CP, x = ( 0.338557 1.721717 7.495970 1.027620 )
Normal form coefficient of CP s= -7.702949e+003
We note in particular the existence of a cusp point of 16-cycle.
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 [28, Fig. 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 curve of Neutral Saddles of fixed points of the third iterate emanates [28, Ch 9]. We compute this curve by
branch switching at the R3 point. This curve is presented in Fig. 10. We note that, however interesting, it is not a bifurcation
curve. In general, continuation of the neutral saddle curve from a 1:3 resonance may help to find a Neimark–Sacker curve
of period 3 cycles. The closed curve destroys this hope.
3.5. Computation of Arnol’d 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 Fig. 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 measures 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 Arnol’d tongue as a limit point of cycles. From this, the boundary curves can be computed by a
continuation in two free parameters. In Fig. 11 we present an Arnol’d tongue rooted in a weak 2: 5 resonant Neimark–Sacker
1464 R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465
Fig. 11. An Arnol’d tongue rooted in a weak 2 : 5 resonant Neimark–Sacker point.
Fig. 12. A stability region of 5-cycles bounded by LP5
, PD5
and NS5
curves.
point. Its computation was started from a stable 5-cycle with x = 0.460832, y = 3.136574, a = 12.7, b = 3, ϵ = 1.327634,
and d = 3.5. The boundary LP5
curve further contains a R1 and LPPD bifurcation points.
To further compute the stability region of the 5-cycle inside the tongue we compute the NS5
and PD5
curves rooted in
these codimension 2 points. These curves intersect in an R2 point and are depicted in Fig. 12. The stability region is indicated
as Ω5
S .
From the ecological point of view, this means that we have described a parameter region where predator and prey can
coexist in a stable way and reproduce their densities every fifth year.
4. Concluding remarks
We investigated the dynamical behavior of a discrete-time predator–prey model with Holling type II functional response.
In Section 2, we focused on the stability and possible bifurcations of three types of fixed points of the model denoted
F1, F2 and F3 respectively. We established the stability condition and branching behavior of F1 in Proposition 1. Conditions
under which F2 may bifurcate to a flip or a branch point, are derived in Proposition 2. We proved supercriticality of
the flip bifurcations of F2 by computing the corresponding normal form coefficient. Proposition 3 provides the necessary
and sufficient conditions under which F3 is stable. All possible bifurcations of F3 are given in Proposition 4. In Section 3,
we computed curves of fixed points and codim 1 bifurcations of cycles. In particular, we computed curves of flip and
Neimark–Sacker bifurcations of the fourth iterate and fold curves of the fourth, eighths and sixteenth iterates. We computed
two curves of folds of the eighth iterate that are born tangentially at GPD points on the flip curve of the fourth iterate. These
curves bound the stability region of an 8-cycle that is born when a fixed point of the fourth iterate crosses a supercritical flip
R. Khoshsiar Ghaziani et al. / Nonlinear Analysis: Real World Applications 13 (2012) 1451–1465 1465
point. We note the bistability of the 4- and 8-cycles in Fig. 7. Furthermore, curves of fold points of the sixteenth iterate are
computed which bound the stability region of a 16-cycle that appears near a resonance 4 point of the fourth iterate. Finally,
we described a parameter region inside an Arnol’d tongue where stable 5-cycles of the map occur.
Acknowledgments
We thank H.G.E. Meijer (Twente, NL) and two referees for several detailed remarks that led to a substantial improvement
in the content as well as the presentation of the paper. The first author also would like to thank Shahrekod University for
the financial support through a research grant.
References
[1] E.M. Elabbasy, S.H. Saker, Dynamics of a class of non-autonomous systems of two non-interacting preys with common predator, J. Appl. Math. Comput.
17 (2005) 195–215.
[2] K.G. Magnusson, O.K. Palsson, Predator–prey interactions of cod and capelin in Icelandic waters, in: ICES Marine Science Symposium, vol. 193, 1991,
pp. 153–170.
[3] H. Gao, H. Wei, W. Sun, X. Zhai, Functions used in biological models and their influence on simulations, Indian J. Mar. Sci. 29 (2000) 230–237.
[4] C.S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Can. 45 (1965)
1–60.
[5] J. Maynard Smith, Mathematical Ideas in Biology, Cambridge University Press, 1968.
[6] L. Edelstein-Keshet, Mathematical Models in Biology, SIAM, Philadelphia, PA, 2005.
[7] S.R.J. Jang, Allee effects in a discrete-time host–parasitoid model, J. Difference Equ. Appl. 12 (2) (2006) 165–181.
[8] M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, 2001.
[9] X. Liu, D. Xiao, Complex dynamic behaviors of a discrete-time predator–prey system, Chaos Solitons Fractals 32 (2007) 80–94.
[10] R.E. Mickens, A nonstandard finite-difference scheme for the Lotka–Volterra system, Appl. Numer. Math. 45 (2003) 309–314.
[11] K. Murakami, Stability and bifurcation in a discrete-time predator–prey model, J. Difference Equ. Appl. 10 (2007) 911–925.
[12] M.G. Neubert, M. Kot, The subcritical collapse of predator populations in discrete-time predator–prey models, Math. Biosci. 110 (1992) 45–66.
[13] R.J. Sacker, H.F. Von Bremen, A new approach to cycling in a 2-locus 2-allele genetic model, J. Difference Equ. Appl. 9 (5) (2003) 441–448.
[14] D. Summers, C. Justian, H. Brian, Chaos in periodically forced discrete-time ecosystem models, Chaos Solitons Fractals 11 (2000) 2331–2342.
[15] Y. Takeuchi, Global Dynamical Properties of Lotka–Volterra Systems, World Scientific, Singapore, 1996.
[16] M. Danca, S. Codreanu, B. Bako, Detailed analysis of a nonlinear prey–predator model, J. Biol. Phys. 23 (1997) 11–20.
[17] K.P. Hadeler, I. Gerstmann, The discrete Rosenzweig model, Math. Biosci. 98 (1990) 49–72.
[18] M. Rosenzweig, Paradox of enrichment: destabilization of exploitation ecosystems in ecological time, Science 171 (1971) 385–387.
[19] H.N. Agiza, E.M. Elabbasy, H. EL-Metwally, A.A. Elsadany, Chaotic dynamics of a discrete prey–predator model with Holling type II, Nonlinear Anal.
RWA 10 (2009) 116–129.
[20] S. Li, W. Zhang, Bifurcations of a discrete prey–predator model with Holling type II functional response, Discrete Contin. Dyn. Syst. 14 (2010) 159–176.
[21] R. Khoshsiar Ghaziani, W. Govaerts, C. Sonck, Codimension-two bifurcations of fixed points in a class of discrete prey–predator systems, Discrete Dyn.
Nat. Soc. 2011 (2011) doi:10.1155/2011/862494, Article ID 862494, 27 pages.
[22] W. Govaerts, R. Khoshsiar Ghaziani, Yu.A. Kuznetsov, H.G.E. Meijer, Numerical methods for two-parameter local bifurcation analysis of maps, SIAM J.
Sci. Comput. 29 (6) (2007) 2644–2667.
[23] W. Govaerts, Yu.A. Kuznetsov, MatCont: a Matlab software project for the numerical continuation and bifurcation study of continuous and discrete
parameterized dynamical systems. www.sourceforge.net.
[24] K. Gopalsamy, Stability and Oscillations in Delay Differential Equations of Population Dynamics, Kluwer Academic, Dordrecht, Norwell, MA, 1992.
[25] A.J. Lotka, Elements of Mathematical Biology, Dover, New York, 1926.
[26] V. Volterra, Opere Matematiche: Memorie e Note, vol. V, Acc. Naz. dei Lincei, Roma, Cremon, 1926.
[27] S.B. Hsu, T.W. Hwang, Global stability for a class of predator–prey systems, SIAM J. Appl. Math. 55 (1995) 763–783.
[28] Yu.A. Kuznetsov, Elements of Applied Bifurcation Theory, third ed., Springer-Verlag, New York, 2004.
[29] J.D. Murray, Mathematical Biology, second corrected ed., Springer, Berlin, Heidelberg, New York, 1993.
[30] E.L. Allgower, K. Georg, Numerical Continuation Methods: An Introduction, Springer-Verlag, 1990.
[31] R.L. Kraft, Chaos, cantor sets, and hyperbolicity for the logistic maps, Amer. Math. Monthly 106 (1999) 400–408.