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

A Host-parasitoid Dynamics with Allee and Refuge
effects

Burcin Kulahcioglu1, Unal Ufuktepe2,∗, Antonios Kalampakas2

1 Faculty of Engineering and Architecture, Beykoz University, Istanbul, Turkey
2 College of Engineering and Technology, American University of the Middle East, Kuwait

Abstract. We examine a discrete-time host–parasitoid model incorporating simultaneously an
Allee and a refuge effect on the host. We investigate existence of a positive fixed point, local
asymptotic stability, global stability of the fixed points and bifurcations. Numerical examples are
given for verification of the theoretical results and we compare the model with existing data from
the literature.

2020 Mathematics Subject Classifications: 92D25, 39A30, 92D30

Key Words and Phrases: Discrete dynamical systems, Host-parasitoid, Nicholson-Bailey, Bi-
furcation, Allee effect, Refuge effect

1. Introduction

Discrete dynamic systems (DDS) in general have been studied extensively and particu-
larly in the recent period in view of the Coronavirus disease pandemic [16],[21],[1],[17],[22].
In addition, it has been shown that population dynamics implementing DDS can produce
efficient computational results for host-parasitoid interactions with numerical simulations
[18],[9],[20], see also [11]. Parasitoids are insect species whose larvae developed by feeding
on the bodies of other arthropods, usually killing them. Larvae emerge from the host and
developed into free-living adults. The adults then lay their eggs in a subsequent generation
of hosts. The host-parasitoid dynamics are increasingly being studied for habitat man-
agement [13], controlling invasive pests [5] and [14], identifying movement among habitat
patches [4] and internal Host-Parasitoid variations [19]. The general framework for the
study of discrete-type host-parasitoid model is the following:

Ht+1 = rHtf(Ht, Pt),
Pt+1 = βHt (1− f(Ht, Pt)) ,

(1)

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v16i2.4707

Email addresses: burcinkulahcioglu@beykoz.edu.tr (B. Kulahcioglu),
unal.ufuktepe@aum.edu.kw (U. Ufuktepe), antonios.kalampakas@aum.edu.kw (A. Kalampakas)

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



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 737

whereHt and Pt are population sizes of host and parasitoids, respectively, in the generation
t. In the host population, r is the reproduction rate of hosts, and the function f(Ht, Pt)
is the fraction of host escaping parasitism. In the parasitoid equation, β is the average
number of eggs (larvae) released by parasitoid on a single host.

An ecologically interesting model developed to describe population dynamics of a cou-
pled host-parasitoid system is the Nicholson–Bailey model [12] which adds further assump-
tions to the general case. Since most parasitoid larvae require a specific life-stage of the
host, the parasitoid and host generations are lined to one another [7]. Hence the total
number of encounters with hosts by parasitoids is in direct proportion to host density and
the encounter number is distributed randomly among the available hosts. Using Poisson
distribution model 1 becomes:

Ht+1 = rHte
−aPt

Pt+1 = βHt

(
1− e−aPt

)
,

(2)

where e−aPt stands for probability of not to be infested by parasitoid and 1− e−aPt is the
probability of being infested by parasitoid at time t.

Initial studies on host-parasitoid interactions employed Equations (2), known as Nicholsan-
Bailey host-parasitoid model. Several other host-parasitoid models have since been devel-
oped as versions of Nicholsan-Bailey (see e.g. [15]-[2]).

In [11] the 2-year oscillations in abundance of Xestia moths were investigated by testing
a hypothesized Xestia-Ophion interaction. Using statistical modelling of time-series data
the authors provided evidence demonstrating the validity of an underlying host-parasitoid
hypothesis (see also [15]). The governing model is

Nt+2 = Ntg(Nt)f(Pt),
Pt+1 = Nt(1− f(Pt)),

(3)

where Nt and Pt are the population sizes of the host (Xestia) and parasitoid (Ophion),
respectively. The function g(N) denotes the population growth rate of the host and the
function f(P ) is the fraction of hosts that are not parasitised. This model assumes a
2-year host life cycle as opposed to single year life cycle of the parasitoid which explains
the apparent 2-year population oscillation. For this reason there are two independent host
populations one for even and one for odd years. Although there is no direct interaction
between the two host cohorts, they are coupled through the the parasitoid population.
When we checked their open source data with Minitab (see Figure 12., Figure 13., and
Figure 14.) there is only one fixed point which is (0,0) and the the fraction of the host
escaping parasitism is periodic and piece-vise exponential because of the seasonal behavior
of the parasitoid.

Unlike the predator in the general prey-predator models, the parasitoid in the host-
parasitoid models has a production rate closely defined by interactions between host and
parasitoid. Hence, generally it is sufficient to use a density dependent population model
for the host, which is also thought to stabilize the system [2]. In [10], 3 different model



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 738

types were proposed according to the ordering of parasitism and density dependence. In
this approach, when parasitism acts first and density dependence takes places only on the
survivors from parasitism (i.e. Htf(Pt)) the model type is:

Ht+1 = Htg(Ht, f(Pt))f(Pt),
Pt+1 = Ht (1− f(Pt)) .

(4)

Choosing f(Pt) = exp(−bPt), g(Ht, f(Pt)) = λ/(1 + kHtexp(−bPt)) and adding β multi-
plier to the second equation, we obtain:

Ht+1 =
λHt

1 + kHte−bPt
e−bPt ,

Pt+1 = βHt(1− e−bPt).
(5)

This model (5) was formulated and analyzed in [6]. In this setup, the host population in
the absence of the parasitoid is modeled by Beverton-Holt equation λH/(1 + kH). The
parameter β represents the average number of egg (larvae) released by parasitoid on a
single host and all parameters are positive. Allee and the refuge effects were separately
added to this model in [8].

In this paper, we incorporate both effects simultaneously and study the dynamics of
the resulting model. Our discrete-time host–parasitoid model is presented in section 2.
In section 3 we determine the conditions under which a positive fixed point exists and is
unique. In Section 4 we study local asymptotic stability, global fixed point stability and
types of bifurcations. Numerical examples for verification of the theoretical results are
given in Section 5 where the data of [11] are checked with Minitab in comparison to our
model.

2. Host–parasitoid model with Allee and refuge effects

In this section we will introduce a discrete-time host–parasitoid model incorporating
simultaneously an Allee and a refuge effect on the host. We start from the system of
equations (5). It has four parameters, namely, λ, b, k and β. We can simplify it by
substituting Nt = βHt and yt = bPt.

Nt+1 =
λNt

1 + k
βNte−yt

e−yt ,

yt+1 = bNt (1− e−yt) .

(6)

Let xt = bNt and c = k/(βb), then the model (6) can be written as follows:

xt+1 =
λxte

−yt

1 + cxte−yt
,

yt+1 = xt (1− e−yt) ,
(7)



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 739

where xt and yt denotes the population sizes of host and parasitoid, respectively at time
t. We add a constant proportion refuge effect and the mate limitation Allee effect to the
model (7). We obtain:

xt+1 =
(1−Ψ)λxt
1 + cxt

xt
s+ xt

+
Ψλxte

−yt

1 + cxte−yt

xt
s+ xt

,

yt+1 = Ψxt (1− e−yt) ,

(8)

where 0 < Ψ ≤ 1 is the proportion of host available to parasitoid, (1−Ψ) is proportion of
protective refuge, and s > 0 is the Allee effect constant. In the following sections we will
study the dynamics of (8).

3. Fixed Points

In this section we find the fixed points of the system, and identify conditions for the
existence and uniqueness of a positive fixed point.

x =
(1−Ψ)λx

1 + cx

x

s+ x
+

Ψλxe−y

1 + cxe−y

x

s+ x
,

y = Ψx (1− e−y) .

(9)

(i) F0 = (0, 0) is extinction fixed point for all values of parameters.

(ii) F1 = (A+
√
A2−4cs
2c , 0) and F2 = (A−

√
A2−4cs
2c , 0) are axial fixed points whereA = λ−cs−1

for (1 +
√
cs)

2 ≤ λ (extinction of the parasite).

(iii) Let F3 = (x∗, y∗) (the host and the parasite survive) be the positive fixed point,
where x∗ = y∗

Ψ(1−e−y∗ )
by the second equation of (9). First, we show the existence of F3.

If we substitute x = y/(Ψ(1− e−y)) in the first equation of (9), we obtain

1

λ
= f(y) + g(y), (10)

where

f(y) =
(1−Ψ)Ψey (ey − 1) y

(Ψ (ey − 1) s+ eyy) (ey(Ψ + cy)−Ψ)
, (11)

and

g(y) =
Ψ2 (ey − 1) y

(Ψ (ey − 1) + cy) (Ψ (ey − 1) s+ eyy)
. (12)

Since

g′(y) = − Ψ2eyy (ey − y − 1)

(Ψ (ey − 1) + cy) (Ψ (ey − 1) s+ eyy)2
−



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 740

Ψ3 (ey − 1) (1 + ey(y − 1))

(Ψ (ey − 1) + cy)2 (Ψ (ey − 1) s+ eyy)
< 0, (13)

the function g(y) is a decreasing function of y.

And we have

lim
y→0

g(y) =
Ψ2

(c+Ψ)(1 + Ψs)
(14)

Now, we continue with the graphical properties of f(y).

f ′(y) =
(1−Ψ)Ψey (ey − y − 1)

(
Ψ2 (ey − 1)2 s− ce2yy2

)
(Ψ (ey − 1) s+ eyy)2 (Ψ− ey(Ψ + cy))2

. (15)

If s < c
Ψ2 , then f ′(y) < 0, which means f(y) is decreasing function. And also

lim
y→0

f(y) =
(1−Ψ)Ψ

(c+Ψ)(1 + Ψs)
. (16)

Let h(y) = f(y) + g(y), h(y) is a decreasing function of y if s < c
Ψ2 . And

lim
y→0

h(y) =
(1−Ψ)Ψ

(c+Ψ)(1 + Ψs)
+

Ψ2

(c+Ψ)(1 + Ψs)
=

Ψ

(c+Ψ)(1 + Ψs)
(17)

The positive fixed point exists if λ = 1
h(y) where h(y) ̸= 0. 1

h(y) is an increasing function
of y and

lim
y→0

1

h(y)
=

(c+Ψ)(1 + Ψs)

Ψ
. (18)

As a result (See Figure 1) a positive fixed point exists and it is unique if s < c
Ψ2 and

λ > (c+Ψ)(1+Ψs)
Ψ . Otherwise, neither the existence nor the uniqueness is guaranteed (Figure

2 and Figure 3).

4. Local Asymptotic Stability and Global Behaviors

In this section we investigate local asymptotic stability, global stability and bifurca-
tions. The corresponding Jacobian matrix of (8) is as follows:

JA=


λx((ey+cx)2(x+s(2+cx))−Ψ(−1+ey)(ey(2s+x+csx)+cx(s−cx2)))

(s+x)2(1+cx)2(ey+cx)2
− λΨeyx2

(s+x)(ey+cx)2

Ψ(1− e−y) Ψe−yx

.

(i) JA(F0) is 0 matrix, so F0 is stable for all parameter values by Trace-determinant
Theorem (see Figure 4. , and Figure 9.).



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 741

0 2 4 6 8 10
0.02

0.04

0.06

0.08

0.10

0.12

y

1

Λ

hHyL

Figure 1: The positive fixed point of the model (8) exists if the two functions intersect. If s < c
Ψ2 and

λ > (c+Ψ)(1+Ψs)
Ψ

, they intersect once (unique positive fixed point). The figure depicts this case for λ = 10,
c = 1, s = 3, and Ψ = 0.5.

0 10 20 30 40

0.26

0.27

0.28

0.29

0.30

0.31

0.32

0.33

y

1

Λ

hHyL

Figure 2: The case s > c
Ψ2 for λ = 3.8, c = 0.01, s = 5, and Ψ = 0.5.Multiple intersections, multiple positive

fixed points.

0 2 4 6 8 10
0.0

0.1

0.2

0.3

0.4

0.5

y

1

Λ

hHyL

Figure 3: The case λ < (c+Ψ)(1+Ψs)
Ψ

for λ = 2, c = 1, s = 3, and Ψ = 0.5. The positive fixed point of the
model (8) exists if the two functions intersect.No intersection, no positive fixed point.

(ii) The Jacobian matrix evaluated at F1 is



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 742

Figure 4: The determinations of eigenvalues in all the regions in the Det-Trace plane [3].

JA(F1) =

 1−
√
A2−4cs

λ

Ψ(1−λ−cs−
√
A2−4cs)

2λc

0
Ψ(A+

√
A2−4cs)
2c

.

We have λ1 = 1−
√
A2−4cs

λ and λ2 =
Ψ(A+

√
A2−4cs)
2c . Hence, |λ1,2| < 1 if s < c

Ψ2 and

(1 +
√
cs)

2
< λ < (c+Ψ)(1+Ψs)

Ψ . Under these conditions F1 is stable (see Figure 4. ,
and Figure 11.).

(iii) The Jacobian matrix at F2 is

JA(F2) =

 1 +
√
A2−4cs

λ

Ψ(1−λ−cs+
√
A2−4cs)

2λc

0
Ψ(A−

√
A2−4cs)
2c

.

We have λ1 = 1 +
√
A2−4cs

λ and λ2 =
Ψ(A−

√
A2−4cs=det()
2c .

Case 1: If A2 = 4cs ⇒ λ1 = 1 and λ2 = ΨA
2c = det(JA) so we have non-hyperbolic

case and it is saddle-node (Fold) Bifurcation if 0 < λ2 < 1 then by the Center Man-
ifold theory[3] F2 is oscillatory: It is oscillatory source If ΨA > 2c. It is oscillatory
saddle if ΨA < 2c.

Case 2: If A2 − 4cs > 0 then λ1 > 1 and F2 is a spiral source. (see Figure 5. ,Figure
6, and Figure 10.)

Numerical simulations to investigate the behavior of the fixed point F3 will be implemented
in Section 5.

For the model (8), F0 = (0, 0) is globally asymptotically stable if λ < 1. We have

xt+1 =
(1−Ψ)λxt
1 + cxt

xt
s+ xt

+
Ψλxt

eyt + cxt

xt
s+ xt



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 743

Figure 5: Classification of the fixed points and bifurcations in all the regions in the Det-Trace plane.

By comparison we obtain

xt+1 ≤
λxt

1 + cxt

xt
s+ xt

<
λxt

1 + cxt
< λxt < xt

if λ < 1. Then
lim
t→∞

xt = 0,

which also implies
lim
t→∞

yt = 0.

Hence, the fixed point F0 is globally attracting if λ < 1, and it is locally asymptotically
stable for all values of parameters. Thus, it is globally asymptotically stable if λ < 1.

The global behavior of the other fixed points is omitted, since local stability of (0, 0)
does not allow any other fixed points to behave stable, globally. Since F0 is locally asymp-
totically stable for all values of parameters, at least with a very small initial value, the
system goes to extinction even with a large growth parameter λ. As a result, the other
points cannot be globally attractive.

5. Numerical Simulations

In this section we verify the theoretical results of our model by numerical simulations
and compare it by using the data of [11] by using Minitab. In order to investigate the
impact of λ to the model, we fix c = 1, Ψ = 0.5, and s = 0.5 and consider cases, which
satisfy the conditions s < c

Ψ2 and λ > (c+Ψ)(1+Ψs)
Ψ .

Assume that λ = 4. Then the positive fixed point is F3 = (2.16476, 0.160471) with
corresponding eigenvalues (0.834265, 0.608611). As a result, F3 is a sink. In Figure 6, we



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 744

0.5 1.0 1.5 2.0
x

0.05

0.10

0.15

0.20

y

8x, y<

Figure 6: Phase portrait for the model (8) for Ψ = 0.5, c = 1, s = 0.5 and λ = 4.

50 100 150 200

2.05

2.10

2.15

2.20

2.25

Figure 7: a) Time series of x for the model (8) for Ψ = 0.5, c = 1, s = 0.5 and λ = 4. b) Time series of (x, y)
for the model (8) for λ = 2, s = 3, c = 1, Ψ = 0.5.

give the phase portrait and we give the times series diagrams in Figure 7. If we assign a
large growth rate λ = 10, then the fixed point F3 = (5.31214, 2.41985) and eigenvalues are
0.312199±0.669214i. F3 is a spiral sink. We observe that the positive fixed point remains
stable for a large interval of growth parameter λ. In Figure 8, we give the bifurcation
diagram.

Finally, we give the basin of attraction for the model (8) in Figures 9-11. In the
analytical results, we show that the fixed point F0 = (0, 0) is stable for all values of
parameters. In the numerical simulations, we show that the positive fixed point F3 is
stable for Ψ = 0.5, c = 1, s = 0.5 and λ = 4. We present the basin of attraction to
show the set of points which are eventually iterated to either F0 or F3 under this set of
parameters. On the other hand, if we let s = 0, that means there is no Allee Effect. By
keeping the other parameters the same, only the positive fixed point is stable. For the
basin of attraction code, we refer to [18].



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 745

Figure 8: The bifurcation diagram of the model (8) for s = 0.5.

Figure 9: Initial point (2,1) λ = 2, s = 3, c = 1, Ψ = 0.5.

6. Conclusion

In this paper, the complex dynamics of a nonlinear discrete-time host-parasitoid sys-
tem with Allee and refuge effects are analyzed. We examined fixed point stability and
bifurcations and investigated global and local fixed point stability. By the study of nu-
merical simulations with the data of [11] using Minitab and our Mathematica codes we
find that the behavior of the system is simultaneously periodic and chaotic. A future
research direction is to create a specific model for this type of data with hybrid models.



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 746

Figure 10: Initial point (2,1) λ = 7, s = 2, c = 1, Ψ = 0.4.

Figure 11: Initial point (2,0) λ = 7, s = 2, c = 1, Ψ = 0.4



B. Kulahcioglu, U. Ufuktepe, A. Kalampakas / Eur. J. Pure Appl. Math, 16 (2) (2023), 736-750 747

Figure 12: Numerical simulation of the stability of the (0,0) fixed point.

Figure 13: Interactions of Ophion and Xest with respect to years.



REFERENCES 748

Figure 14: Time series of Xest.

References

[1] Joel Alba-Pérez and Jorge E Maćıas-Dı́az. Analysis of structure-preserving discrete
models for predator-prey systems with anomalous diffusion. Mathematics, 7(12):1172,
2019.

[2] John R Beddington. Mutual interference between parasites or predators and its effect
on searching efficiency. The Journal of Animal Ecology, pages 331–340, 1975.

[3] Saber N Elaydi. Discrete chaos: with applications in science and engineering. Chap-
man and Hall/CRC, 2007.

[4] Edward W Evans. Dispersal in host–parasitoid interactions: Crop colonization by
pests and specialist enemies. Insects, 9(4):134, 2018.

[5] James R Hepler, Kacie Athey, David Enicks, Paul K Abram, Tara D Gariepy, Elijah J
Talamas, and Elizabeth Beers. Hidden host mortality from an introduced parasitoid:
Conventional and molecular evaluation of non-target risk. Insects, 11(11):822, 2020.

[6] SOPHIA R-J Jang and Jui-Ling Yu. A discrete-time host-parasitoid model. In Pro-
ceedings of the Conference on Differential and Difference Equations and Applications.
Hindawi Publishing Corporation, pages 451–455, 2006.

[7] Abdul Qadeer Khan and Muhammad Naeem Qureshi. Dynamics of a modified
nicholson-bailey host-parasitoid model. Advances in Difference Equations, 2015(1):1–
15, 2015.

[8] Burcin Kulahcioglu and Unal Ufuktepe. A density dependent host-parasitoid model
with allee and refuge effects. In Computational Science and Its Applications–ICCSA



REFERENCES 749

2016: 16th International Conference, Beijing, China, July 4-7, 2016, Proceedings,
Part III 16, pages 228–239. Springer, 2016.

[9] Xiaorong Ma, Qamar Din, Muhammad Rafaqat, Nasir Javaid, and Yongliang Feng. A
density-dependent host-parasitoid model with stability, bifurcation and chaos control.
Mathematics, 8(4):536, 2020.

[10] RM May, MP Hassell, RM Anderson, and DW Tonkyn. Density dependence in host-
parasitoid models. The Journal of Animal Ecology, pages 855–865, 1981.

[11] Marko Mutanen, Otso Ovaskainen, Gergely Várkonyi, Juhani Itämies, Sean WJ
Prosser, Paul DN Hebert, and Ilkka Hanski. Dynamics of a host–parasitoid inter-
action clarified by modelling and dna sequencing. Ecology Letters, 23(5):851–859,
2020.

[12] Alexander J Nicholson and Victor A Bailey. The balance of animal populations.—part
i. In Proceedings of the zoological society of London, volume 105, pages 551–598. Wiley
Online Library, 1935.

[13] Ainara Peñalver-Cruz, Bruno Jaloux, and Blas Lavandero. The host-plant origin
affects the morphological traits and the reproductive behavior of the aphid parasitoid
aphelinus mali. Agronomy, 12(1):101, 2022.

[14] Serge Quilici and Pascal Rousse. Location of host and host habitat by fruit fly
parasitoids. Insects, 3(4):1220–1235, 2012.

[15] M Rost, G Várkonyi, and I Hanski. Patterns of 2-year population cycles in spatially
extended host–parasitoid systems. Theoretical Population Biology, 59(3):223–233,
2001.

[16] Bellie Sivakumar and Bhadran Deepthi. Complexity of covid-19 dynamics. Entropy,
24(1):50, 2022.

[17] Agus Suryanto, Isnani Darti, Hasan S. Panigoro, and Adem Kilicman. A fractional-
order predator–prey model with ratio-dependent functional response and linear har-
vesting. Mathematics, 7(11):1100, 2019.

[18] Unal Ufuktepe and Sinan Kapcak. Applications of discrete dynamical systems with
mathematica. RIMS Kyoto Proceedings, 1909:207–216, 2014.

[19] Saskya Van Nouhuys, Suvi Niemikapee, and Ilkka Hanski. Variation in a host–
parasitoid interaction across independent populations. Insects, 3(4):1236–1256, 2012.

[20] Gergely Várkonyi, Ilkka Hanski, Martin Rost, and Juhani Itämies. Host-parasitoid
dynamics in periodic boreal moths. Oikos, 98(3):421–430, 2002.

[21] Sandra Vaz and Delfim FM Torres. A discrete-time compartmental epidemiological
model for covid-19 with a case study for portugal. Axioms, 10(4):314, 2021.



REFERENCES 750

[22] Jianming Zhang, Lijun Zhang, and Yuzhen Bai. Stability and bifurcation analysis on
a predator–prey system with the weak allee effect. Mathematics, 7(5):432, 2019.