EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 16, No. 1, 2023, 71-83
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

A Discrete Predator-prey Model with Allee and Refuge
Effect

Sinan Kapçak

College of Engineering and Technology, American University of the Middle East, Kuwait

Abstract. We consider a predator-prey model, where some prey are completely free from preda-
tion within a temporal or spacial refuge and the predator population is subject to Allee effect. We
study the effect of the presence of refuge and Allee effect on the stability and bifurcation of the
system. We investigate the existence and stability of the model as well as the stability region. We
also obtain the invariant manifolds of the system.

2020 Mathematics Subject Classifications: 37N25, 37D10, 39A28, 39A30, 39A60

Key Words and Phrases: Allee Effect, Refuge Effect, Invariant Manifolds

1. Introduction

In this paper, we investigate a predator-prey model, where some prey are completely
free from predation within a temporal or spacial refuge and predator are subject to Allee
effect. The most common type of spacial refuge, that we investigate here, takes the form
where a constant proportion of the prey population is protected. Some of studies have
investigated the influence of prey refuge and concluded that the refuge used by the prey
has a stabilizing effect on the predator-prey interaction and also that the prey species can
be prevented well from extinction. Some studies on refuge effect can be found in [3, 12, 13].

The Allee effect is characterized by a positive correlation between population size and
the mean individual fitness of a population. When the population density is low, the
population will experience a reduced overall growth rate, and may increase the risk of
extinction. If there are more individuals, population grows more rapidly, and the aggrega-
tion can improve the survival rate of individuals. Recently, many studies have been done
on predator-prey models with Allee effect. Livadiotis, Assas, Elaydi, Kwessi, Dennis [11]
investigated the impact of the Allee effect on the global dynamics of Beddington model.
Some other predator-prey studies with/without Allee effect were conducted in [5–7]/[1, 9].

DOI: https://doi.org/10.29020/nybg.ejpam.v16i1.4661

Email address: sinan.kapcak@aum.edu.kw (S. Kapçak)

https://www.ejpam.com 71 © 2023 EJPAM All rights reserved.



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 72

Maynard Smith [17] proposed the following predator-prey model:

Xt+1 = RXt −
R− 1

XE
X2

t − CXtYt,

Yt+1 =
r

XE
XtYt,

(1)

where Xt and Yt represent the prey and predator population in year t, respectively. XE is
the equilibrium value of Xt in the absence of predator. R and r are maximum reproductive
rates of prey and predator, respectively.

In [14], Murakami investigates the stability and bifurcation, including Neimark-Sacker
bifurcation, of model (1).

In this paper, we study model (1) with prey and predator subject to a refuge and
Allee effect, respectively. Regarding the refuge effect, we assume the presence of spacial
heterogeneity. The refuge model is due to Hassel [3]. In the absence of predator, the
growth of the prey population is logistic. The following model is proposed:

Xt+1 = R(1− b)Xt −
R− 1

XE
(1− b)2X2

t +RbXt −
(R− 1)b2

XE
X2

t − CbXtYt,

Yt+1 =
rb

XE
XtYt

Yt
A+ Yt

.

(2)

In equation (2), the expression (1 − b)Xt is the number of the prey in a protective
refuge at time t, 0 < b < 1, which has logistic growth represented by the expression
R(1− b)Xt − R−1

XE
(1− b)2X2

t . The rest of the prey, bXt, is affected by the predator. The

expression RbXt− (R−1)b2

XE
X2

t −CbXtYt gives this interaction. We provide support for the
stabilizing power of physical refuge. We show that with refuge present, the host and para-
sitoid can, apparently, persist together indefinitely. The expression Yt

A+Yt
is a mate-finding

Allee effect on predator.

Simplifying equation (2), we obtain the following system:

Xt+1 = RXt −
R− 1

XE
(2b2 − 2b+ 1)X2

t − CbXtYt,

Yt+1 =
rb

XE
XtYt

Yt
A+ Yt

.

(3)

Note that, in equation (3), when b = 1 and A = 0, there is no prey in a protective refuge
and no Allee effect. Hence we obtain equation (1) for this particular case. When b = 0,
all the prey population is protected and has logistic growth while the predator population
extincts in the next generation.



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 73

We apply yt = CbYt and xt =
rb
XE

Xt with α = R, β = 2b2−2b+1
rb , and γ = ACb. Hence

we obtain the following system:

xt+1 = xt(α− β(α− 1)xt)− xtyt,

yt+1 =
xty

2
t

γ + yt
,

(4)

where α, β > 0 and γ > 0.

2. Existence of the Fixed Points

In this section, we analyse the existence of fixed points of discrete system (4). Firstly,
we focus on the following isocline equations:

x = x(α− β(α− 1)x)− xy,

y =
xy2

γ + y
.

(5)

2.1. Extinction and Exclusion Fixed Points

In equation (5), if x = 0, we have the extinction fixed point P0 = (0, 0) for any
parameter values. If x ̸= 0 and y = 0, we obtain the exclusion fixed point P1 = ( 1β , 0).

2.2. Coexistence Fixed Points

If x ̸= 0 and y ̸= 0, the isoclines become

y = −β(α− 1)x+ α− 1,

y =
γ

x− 1
.

(6)

The intersection points of the isoclines in equation (6) give us the candidates for coexistence
fixed points which can be found by solving the following quadratic equation:

β(α− 1)x2 − (α− 1)(β + 1)x+ (α+ γ − 1) = 0. (7)

Note that in order to have a positive fixed point, by the second equation of system (6),
it is necessary that x-component of the fixed point is bigger than 1. Therefore, we have
a coexistence fixed point if and only if x∗ is a real root of equation (7) and x∗ > 1. The
possible roots of the equation (7) are the following:

x∗1,2 =
(α− 1)(β + 1)±

√
∆

2β(α− 1)
, (8)

where ∆ = (α− 1)2(β − 1)2 − 4β(α− 1)γ. In order to have positive fixed points, we must
have ∆ ≥ 0 and x∗i > 1 for some i ∈ {1, 2}. Under the condition ∆ > 0, we obtain two



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 74

different scenario: (a) there exist two coexistence fixed points if α > 1 and β < 1. (b)
there exists one coexistence fixed point if α < 1. For the case when ∆ = 0, the condition
that a coexistence fixed point exists is β < 1. We obtain the following theorem for the
condition of existence and the corresponding fixed points of system (4).

(a) Fixed points: P0, P1 (b) Fixed points: P0, P1, P
+
1

(c) Fixed points: P0, P1, P
+
0 (d) Fixed points: P0, P1, P

+
2 , P+

3

Figure 1: Four possible cases for existence of the fixed points. The dashed line is the isocline x = x(α− β(α−
1)x)− xy whereas the solid graph is the isocline y = xy2

γ+y
.

Theorem 1. The following table gives the existence condition for the fixed points of
system (4).



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 75

Fixed Point Condition for Existence

P0 = (0, 0) Always exists

P1 =
(

1
β , 0
)

Always exists

P+
0 =

(
β+1
2β , 2γβ

1−β

)
β < 1 and (β − 1)2 = 4βγ

α−1

P+
1 =

(
(α−1)(β+1)−

√
∆

2β(α−1) , −(α−1)(β−1)+
√
∆

2

)
α < 1

P+
2 =

(
(α−1)(β+1)+

√
∆

2β(α−1) , −(α−1)(β−1)−
√
∆

2

)
β < 1, α > 1, and (β − 1)2 > 4βγ

α−1

P+
3 =

(
(α−1)(β+1)−

√
∆

2β(α−1) , −(α−1)(β−1)+
√
∆

2

)
β < 1, α > 1, and (β − 1)2 > 4βγ

α−1

Figure 2: Existence regions for system (4) in α-β plane.

Figure 2 represents the existence regions for fixed points in α-β plane. By Theorem 1,
extinction fixed point P0 and exclusion fixed point P1 exist for any parameter values. P+

0

exists only on the borderline of E2 and E3. This is the curve where ∆ = 0. Region E1 is
the region for existence of P+

1 . E2 is the region where ∆ > 0 and we have two distinct
positive fixed points P+

2,3.



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 76

3. Stability Analysis of System (4)

In this section, we analyse the stability of system (4). The Jacobian matrix of the
system is

J(x, y) =

(
α− 2β(α− 1)x− y −x

y2

γ+y
xy(2γ+y)
(γ+y)2

)
.

3.1. Stability of the Extinction Fixed Point

For the extinction fixed point P0 = (0, 0), the Jacobian matrix J is given by

J0 =

(
α 0
0 0

)
.

The eigenvalues of the matrix are λ1 = α and λ2 = 0. Hence the fixed point P0 is locally
asymptotically stable if α < 1.

3.2. Stability of the Exclusion Fixed Point

For the exclusion fixed point P1 =
(

1
β , 0
)
, we have the following Jacobian matrix:

J1 =

(
2− α − 1

β

0 0

)
.

The eigenvalues of the Jacobian matrix are λ1 = 2 − α and λ2 = 0. Hence the fixed
point P1 is locally asymptotically stable if 1 < α < 3.

3.3. Stability of the Coexistence Fixed Points

By Theorem 1, the coexistence fixed points exist if ∆ ≥ 0. We start with the following
case:

3.3.1. Stability of P+
1 , P+

2 , and P+
3 : if ∆ > 0

For the case when ∆ > 0, we investigate three possible fixed points P+
1 , P+

2 , P+
3 . Trace

and determinant of the Jacobian matrices J+
1 , J+

2 , J+
3 , respectively, at these fixed points

are given by

tr(J+
1,2,3) =

1

2

(
βγ + γ ±

√
∆

α+ γ − 1
− α(β + 1)±

√
∆+ 6

)
,

det(J+
1,2,3) =

1

2

(
βγ + γ ±

√
∆

α+ γ − 1
− α(β + 1)± 3

√
∆+ 4

)
,



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 77

where ∆ = (α− 1)2(β − 1)2 − 4β(α− 1)γ.

We obtained the stability conditions applying the Tr-Det Formula

|tr(J)| − 1 < det(J) < 1.

Theorem 2. For system (4),

i. P+
1 is unstable.

ii. P+
2 is unstable.

iii. P+
3 is locally asymptotically stable if

2 + 3
√
∆ < α(α+ β)−K < 6 + 2

√
∆,

where K = βγ+γ−
√
∆

α+γ−1 and ∆ = (α− 1)2(β − 1)2 − 4β(α− 1)γ.

Neimark-Sacker Bifurcation

Saddle-Node Bifurcation
Flip Bifurcation

0 2 4 6 8 10
α

0.2

0.4

0.6

0.8

1.0

1.2

1.4

β

Figure 3: (γ = 0.26) Stability region of coexistence fixed point P+
3

Figure 3 represents the stability region of the coexistence fixed point P+
3 in the param-

eter plane α-β. Types of bifurcations are also given. In Table 1 and 2, we give several
parameter values and corresponding phase diagrams. One can see in the table that the
system also exhibits Neimark-Sacker bifurcation (see [4, 8, 10, 15, 16]) for some values of
parameters. This is the case when det(J+

3 ) = 1 and |tr(J+
3 )| < 2. When the parameter

point is in the stability region and close to the Neimark-Sacker borderline, the orbit is
spiraling and just after the parameter point (α, β) crosses the N-S curve, orbits converge
to a closed invariant curve in the phase diagram.



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 78

Table 1: (γ = 0.26) Stability of the coexistence fixed point.

α-β Plane Phase Diagram

3.3.2. Stability of P+
0 : Non-hyperbolic case or the case ∆ = 0

There is only one positive fixed point when ∆ = 0 and β < 1 by Theorem 1: P+
0 =(

β+1
2β , 2γβ

1−β

)
. Applying this point to Jacobian matrix and eliminating γ using the equation

∆ = 0, we obtain

J+
0 =

(
1
2(β − α(β + 1) + 3) −β+1

2β

− (α−1)(β−1)β
β+1

2
β+1

)
.

The eigenvalues of the Jacobian matrix are λ1 = 1 and λ2 =
−α(β+1)2+β2+2β+5

2(β+1) . The fixed

point P+
0 is non-hyperbolic. If |λ2| > 1, then it is unstable. If |λ2| < 1, we have to apply

the center manifold theory [2].

It is more convenient to make a change of variables in system (4) so we can have a
shift from the point P+

0 to (0, 0). Let u = x− β+1
2β and v = y− 2γβ

1−β . Then the new system



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 79

Table 2: (γ = 0.26) Stability of the coexistence fixed point.

α-β Plane Phase Diagram

is

ut+1 =

(
β + 1

2β
+ ut

)(
α− (α− 1)β

(
β + 1

2β
+ ut

))
−
(
β + 1

2β
+ ut

)(
2βγ

1− β
+ vt

)
− β + 1

2β
,

vt+1 =

(
β+1
2β + ut

)(
2βγ
1−β + vt

)2(
2βγ
1−β + vt

)
+ γ

− 2βγ

1− β
.

(9)



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 80

The Jacobian matrix of the system (9) at (0, 0) is

J̃+
0 =

(
A B
C D

)
=

(
1
2(β − α(β + 1) + 3) −β+1

2β

− (α−1)(β−1)β
β+1

2
β+1

)
.

Now we can write the equations in system (9) as

ut+1 = Aut +Bvt + f̃(ut, vt),

vt+1 = Cut +Dvt + g̃(ut, vt),
(10)

where f̃(ut, vt) and g̃(ut, vt) are the rest of the expressions.

Let us assume that the map h, which represents a center manifold, takes the form

h(u) = (β − αβ)u+ au2 + bu3 +O(u4), a, b ∈ R.

Now we have to compute the constants a and b. The function h must satisfy the center
manifold equation

h(Au+Bh(u) + f̃(u, h(u)))− Cu−Dh(u)− g̃(u, h(u)) = 0. (11)

The Taylor series expansions, at the point u = 0, are evaluated for the equation above.
Equating the coefficients of the series and using the equation ∆ = 0, after some manipu-
lations, we obtain

a = − 4(α− 1)β2

α(β + 1)2 − β2 − 3

b =
16β3

(
α3(β + 1)3 − α2(β + 1)(3β(β + 1) + 4) + α

(
3β
(
β2 + β + 1

)
+ 7
)
− β3 + β − 4

)
(α(β + 1)2 − β2 − 3)3

.

Thus on the center manifold v = h(u) we find the map P (u). Because of the lengthy
expressions we omit it here. Calculations show that P ′(0) = 1. Since 0 < β < 1 and
|λ2| < 1, we have α > 1 and α(β + 1)2 > β2 + 3. Therefore,

P ′′(0) =
4(α− 1)β(β + 1)

α(β + 1)2 − β2 − 3
> 0.

Hence, for the map P , the origin is semistable from the left.

Now, we are going to find the stable manifold, which exists when |λ2| < 1 or

−α(β + 1)2 + β2 + 2β + 5 < 2(β + 1).



S. Kapcak / Eur. J. Pure Appl. Math, 16 (1) (2023), 71-83 81

Since the stable manifold is tangent to the eigenvector at the point, let us take

h(u) =
2(β − 1)β

(β + 1)2
u+ au2 + bu3 +O(u4), a, b ∈ R.

This map must satisfy the centre manifold equation (11). We calculate map Q on the
stable manifold and found that

Q′(0) =
−α(β + 1)2 + β2 + 2β + 5

2(β + 1)
.

Because of the long output of the computations we omit them here.

Stable Manifold
Center M

anifold

1.6 1.7 1.8 1.9 2.0
x

0.45

0.50

0.55

0.60

0.65

0.70

0.75

y

Figure 4: (α = 2.75 and β = 0.4) Phase Diagram of system (4) with the invariant manifolds.

Stable and centre manifolds are given in the phase diagram in Figure 4 for the case
∆ = 0. The dashed curve is the centre manifold for which the fixed point is semi-stable.
The solid curve is the stable manifold. An orbit approaching the fixed point is also shown.

4. Conclusions

In this paper, we investigated the stability of a predator-prey model with refuge and
Allee effect. We showed that the presence of a safe refuge, where a portion of the host
is in a safe refuge from predation has a stabilizing effect on the model. The conditions



REFERENCES 82

for the existence of the fixed points are found. We also obtained the invariant manifolds
for the positive fixed point P+

0 . Furthermore, we obtained the stability region for the
coexistence fixed point P+

3 . By numerical computations, we confirm our analytic results.
The Mathematica codes displaying phase diagram can be found in [18].

References

[1] R Asheghi. Bifurcations and dynamics of a discrete predator–prey system. Journal
of Biological Dynamics, 8(1):161–186, 2014.

[2] S Elaydi. Discrete Chaos: With Applications in Science and Engineering, Second
Edition. Chapman & Hall/CRC, 2008.

[3] M P Hassell. The Dynamics of Arthropod Predator-Prey Systems, volume 111.
Princeton University Press, 2020.

[4] ANW Hone, MV Irle, and GW Thurura. On the neimark–sacker bifurcation in a
discrete predator-prey system. Journal of Biological Dynamics, 4(6):594–606, 2010.

[5] S R-J Jang. Allee effects in a discrete-time host-parasitoid model. Journal of Differ-
ence Equations and Applications, 12(2):165–181, 2006.

[6] S R-J Jang. Discrete-time host–parasitoid models with allee effects: Density de-
pendence versus parasitism. Journal of Difference Equations and Applications,
17(04):525–539, 2011.

[7] SR-J Jang and S.L. Diamond. A host-parasitoid interaction with allee effects on the
host. Computers and Mathematics with Applications, 53:89–103, 2007.

[8] S Kapçak, S Elaydi, and Ü Ufuktepe. Stability of a predator–prey model with refuge
effect. Journal of Difference Equations and Applications, 22(7):989–1004, 2016.

[9] S Kapçak, Ü Ufuktepe, and S Elaydi. Stability and invariant manifolds of a general-
ized beddington host-parasitoid model. Journal of Biological Dynamics, 7(1):233–253,
2013.

[10] MRS Kulenovic, E Pilav, and E Silic. Naimark-sacker bifurcation of a certain second
order quadratic fractional difference equation. Journal of Mathematical and Compu-
tational Science, 4(6):1025–1043, 2014.

[11] G Livadiotis, L Assas, B Dennis, S Elaydi, and E Kwessi. A discrete-time host–
parasitoid model with an allee effect. Journal of Biological Dynamics, 9(1):34–51,
2015.

[12] R M May. Host-parasitoid systems in patchy environments: a phenomenological
model. The Journal of Animal Ecology, pages 833–844, 1978.



REFERENCES 83

[13] R M May, MP Hassell, RM Anderson, and DW Tonkyn. Density dependence in
host-parasitiod models. Journal of Animal Ecology, 50:855–865, 1981.

[14] K Murakami. Stability and bifurcation in a discrete-time predator–prey model. Jour-
nal of Difference Equations and Applications, 13(10):911–925, 2007.

[15] R J Sacker. Introduction to the 2009 re-publication of the ‘neimark–sacker bifurcation
theorem’. Journal of Difference Equations and Applications, 15(8-9):753–758, 2009.

[16] R John Sacker. On invariant surfaces and bifurcation of periodic solutions of ordinary
differential equations. New York University, 1964.

[17] J M Smith. Mathematical ideas in biology. CUP Archive, 1968.

[18] Ü Ufuktepe and S Kapçak. Applications of discrete dynamical systems with mathemat-
ica (study of mathematical software and its effective use for mathematics education).
RIMS Kyoto University, 1909:207–216, 2014.