# Bifurcation Analysis of a Discrete-Time Two-Species Model

Authors: A. Q. Khan
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

We study the local dynamics and bifurcation analysis of a discrete-time modified Nicholson–Bailey model in the closed first quadrantR
+
2. It is proved that model has two boundary equilibria: O
0
,
0
,
A
ζ
1
−
1
/
ζ
2
,
0, and a unique positive equilibrium B
r
e
r
/
e
r
−
1
,
r under certain parametric conditions. We study the local dynamics along their topological types by imposing method of Linearization. It is proved that fold bifurcation occurs about the boundary equilibria: O
0
,
0
,
A
ζ
1
−
1
/
ζ
2
,
0. It is also proved that model undergoes a Neimark–Sacker bifurcation in a small neighborhood of the unique positive equilibrium B
r
e
r
/
e
r
−
1
,
r and meanwhile stable invariant closed curve appears. From the viewpoint of biology, the stable closed curve corresponds to the period or quasi-periodic oscillations between host and parasitoid populations. Some simulations are presented to verify theoretical results. Finally, bifurcation diagrams and corresponding maximum Lyapunov exponents are presented for the under consideration model.


---

## Body

## 1. Introduction

The usual framework for the discrete-time host–parasite models is:(1)
X
n
+
1
=
b
X
n
f
X
n
,
Y
n
Y
n
+
1
=
c
X
n
1
−
f
X
n
,
Y
n
,whereX
n and Y
n represent the population size of the host and parasite in successive generations n and n
+
1 respectively. The parameter b is the host finite rate of increase in the absence of parasites, c is the biomass conversion constant and f is the function defining the fractional survival of hosts from parasitism. The simplest version of this model is that of Nicholson, and Nicholson and Bailey who explored in depth a model in which the proportion of hosts escaping parasitism is given by the zero term of the Poisson distribution [1–3]:(2)
f
X
n
,
Y
n
=
e
−
a
Y
n
,wherea
Y
n are the mean encounters per host. Thus, 1
−
e
−
a
Y
n is the probability of a host will be attacked. Using (2) in (1), one gets(3)
X
n
+
1
=
b
X
n
e
−
a
Y
n
Y
n
+
1
=
c
X
n
1
−
e
−
a
Y
n
.In 2014, Qureshi et al. [4] have investigated the asymptotic behavior of the following Nicholson–Bailey model:(4)
X
n
+
1
=
R
x
n
e
−
a
y
n
Y
n
+
1
=
x
n
1
−
e
−
a
y
n
,whereR, a and initial conditions x
0
,
y
0 are positive real numbers. Further in 2015, Khan and Qureshi [5] have investigated the dynamics of the following modified Nicholson–Bailey model:(5)
X
n
+
1
=
b
X
n
e
−
a
Y
n
1
+
d
X
n
Y
n
+
1
=
c
X
n
1
−
e
−
a
Y
n
,wherea
,
b
,
c
,
d and initial conditions x
0
,
y
0 are positive real numbers. Our aim in this paper is to explore the local dynamics along with topological classification and bifurcation analysis of the model (5). First, we make the following rescaling transformations:(6)
X
n
=
x
n
b
,
Y
n
=
y
n
a
,then system (5) becomes(7)
x
n
+
1
=
b
x
n
e
−
y
n
1
+
d
/
b
x
n
y
n
+
1
=
a
c
b
x
n
1
−
e
−
y
n
.For simplicity, we assume thatb
=
a
c, and then model (7) becomes:(8)
x
n
+
1
=
ζ
1
x
n
e
−
y
n
1
+
ζ
2
x
n
y
n
+
1
=
x
n
1
−
e
−
y
n
,whereζ
1
=
b
>
0 and ζ
2
=
d
/
b
>
0.The rest of the paper is organized as follows: Section2 deals with the study of existence of equilibria of the model (8). In Section 3, we study the local dynamics and existence of bifurcations about equilibria: O
0
,
0
,
A
ζ
1
−
1
/
ζ
2
,
0
,
B
r
e
r
/
e
r
−
1
,
r of the model. Section 4 deals with the study of Neimark–Sacker bifurcation about B
r
e
r
/
e
r
−
1
,
r of the model (8). Numerical simulations along with discussion are presented in the last Section.


## 2. Existence of Equilibria of the Discrete-Time Model (8)

In this Section, we study the existence of equilibria of the model (8) in R
+
2. The results about the existence of equilibria are summarized as follows:Lemma 1.
Discrete-time model (8) has at least two boundary equilibria and the unique positive equilibrium point in R
+
2.  More precisely,(i)
For all parametric valuesζ
1 and ζ
2, model (8) has a unique equilibrium point: O
0
,
0;
(ii)
Ifζ
1
>
1 then model has boundary equilibrium point: A
ζ
1
−
1
/
ζ
2
,
0;
(iii)
Suppose that(9)
z
=
G
y
^
=
e
y
^
+
ζ
2
y
^
e
2
y
^
e
y
^
−
1
,
and z
=
ζ
1 and when ζ
1
>
ζ
2
+
1, the curve z
=
G
y
^ intersect the line z
=
ζ
1 at r, say. If ζ
1
>
ζ
2
+
1 then there exist a unique r such that B
r
e
r
/
e
r
−
1
,
r has a unique positive equilibrium point of (8).Proof.
For finding number of equilibria of the model (8), we have to solve the following system of equations:(10)
x
^
=
ζ
1
x
^
e
−
y
^
1
+
ζ
2
x
^
y
^
=
x
^
1
−
e
−
y
^
.
(i)
Letx
^
=
0, then 1
s
t equation of system (10) satisfied identically and from 2
n
d equation we obtain y
^
=
0. So system (10) has always equilibrium O
0
,
0 for all parameter values ζ
1
,
ζ
2
>
0.
(ii)
Lety
^
=
0, then 2
n
d equation of (10) satisfied identically and from 1
s
t equation we obtain x
^
=
ζ
1
−
1
/
ζ
2. Hence system has boundary equilibrium A
ζ
1
−
1
/
ζ
2
,
0 if  ζ
1
>
1.
(iii)
Now we locate the unique positive equilibrium of (10) in R
+
2. For this, let x
^
≠
0, then (10) becomes
(11)
1
=
ζ
1
e
−
y
^
1
+
ζ
2
x
^
y
^
=
x
^
1
−
e
−
y
^
.
Now eliminatingx
^ from (11), one gets(12)
ζ
1
=
e
y
^
+
ζ
2
y
^
e
2
y
^
e
y
^
−
1
.
Denote,(13)
G
y
^
=
e
y
^
+
ζ
2
y
^
e
2
y
^
e
y
^
−
1
.
Then they
^-coordinates of positive equilibria of (8) are the corresponding y
^-coordinates of the point of intersection of z
=
G
y
^ and z
=
ζ
1 with y
^
>
0. By calculating derivative of F
y
^, one get(14)
G
′
y
^
=
e
y
^
+
ζ
2
e
2
y
^
e
y
^
−
1
2
y
^
+
1
−
y
^
e
y
^
e
2
y
^
−
1
2
>
0
,
∀
y
^
>
0
.
Moreover(15)
lim
y
^
→
0
+
G
y
^
=
lim
y
^
→
0
+
e
y
^
+
ζ
2
y
^
e
2
y
^
e
y
^
-
1
,
=
ζ
2
+
1
.
So, ifζ
1
≤
ζ
2
+
1, then there exists no intersection point of z
=
G
y
^ and z
=
ζ
1. This implies that model (8) has no positive equilibria if ζ
1
≤
ζ
2
+
1. And if ζ
1
>
ζ
2
+
1, then there exists a unique point of intersection ζ
1
,
r of z
=
G
y
^ and z
=
ζ
1 with r
>
0 (see Figure 1). Therefore, if ζ
1
>
ζ
2
+
1 then (8) has positive equilibrium point and the positive equilibrium point of (8) is unique. We denote it by B
r
e
r
/
e
r
−
1
,
r where r is the positive solution of (12).Figure 1
z
=
G
y
^
,
z
=
ζ
1. (a) ζ
1
<
ζ
2
+
1. (b) ζ
1
=
ζ
2
+
1. (c) ζ
1
>
ζ
2
+
1.
(a)
(b)
(c)


## 3. Local Dynamics and Existence of Bifurcations about Equilibria:O
0
,
0, A
ζ
1
−
1
/
ζ
2
,
0, B
r
e
r
/
e
r
−
1
,
r of the Model  (8)

In this Section, we will study the local dynamics of (8) about O
0
,
0
,
A
ζ
1
−
1
/
ζ
2
,
0, and B
r
e
r
/
e
r
−
1
,
r. The Jacobian matrix J
x
^
,
y
^ of (8) about equilibrium x
^
,
y
^ becomes(16)
J
x
^
,
y
^
=
ζ
1
e
−
y
^
1
+
ζ
2
x
^
2
−
ζ
1
x
^
e
−
y
^
1
+
ζ
2
x
^
1
−
e
−
y
^
x
^
e
−
y
^
.And its characteristic equation is(17)
κ
2
−
p
x
^
,
y
^
κ
+
q
x
^
,
y
^
=
0
,where(18)
p
x
^
,
y
^
=
ζ
1
e
−
y
^
1
+
ζ
2
x
^
2
+
x
^
e
−
y
^
,
q
x
^
,
y
^
=
ζ
1
x
^
e
−
2
y
^
1
+
ζ
2
x
^
2
+
ζ
1
x
^
e
−
y
^
1
−
e
−
y
^
1
+
ζ
2
x
^
.Lemma 2.
For equilibriumO, the following holds:(i)
O is a sink if ζ
1
<
1;
(ii)
O is never source;
(iii)
O is a saddle if ζ
1
>
1;
(iv)
O is nonhyperbolic if ζ
1
=
1.From Lemma2, we can see that one of the eigenvalues about the equilibrium O
0
,
0 is 1. So fold bifurcation may occurs when parameter vary in the small neighborhood of ζ
1
=
1.Lemma 3.
ForA
ζ
1
−
1
/
ζ
2
,
0, the following holds:(i)
A
ζ
1
−
1
/
ζ
2
,
0 is a sink if ζ
2
>
ζ
1
−
1;
(ii)
A
ζ
1
−
1
/
ζ
2
,
0 is never source;
(iii)
A
ζ
1
−
1
/
ζ
2
,
0 is a saddle if ζ
2
<
ζ
1
−
1;
(iv)
A
ζ
1
−
1
/
ζ
2
,
0 is nonhyperbolic if ζ
2
=
ζ
1
−
1.We can easily see that if condition (iv) of Lemma3 hold then one of the eigenvalues about equilibrium A
ζ
1
−
1
/
ζ
2
,
0 is 1. So fold bifurcation may occur when parameters vary in a small neighborhood of ζ
2
=
ζ
1
−
1. And we denote the parameters satisfying ζ
2
=
ζ
1
−
1 as(19)
F
A
=
ζ
1
,
ζ
2
:
ζ
2
=
ζ
1
−
1
,
ζ
1
,
ζ
2
>
0
.Hereafter, we will investigate the local dynamics of (8) about B
r
e
r
/
e
r
−
1
,
r by using Lemma 2.2 of [6]. The Jacobian matrix J
B
r
e
r
/
e
r
−
1
,
r of linearized system of (8) about B
r
e
r
/
e
r
−
1
,
r is(20)
κ
2
−
p
κ
+
q
=
0
,where(21)
p
=
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
+
r
e
r
−
1
,
q
=
r
e
r
−
1
+
ζ
2
r
e
r
+
r
.Moreover eigenvalues ofJ
B
r
e
r
/
e
r
−
1
,
r about B
r
e
r
/
e
r
−
1
,
r is given by(22)
κ
1
,
2
=
−
p
±
Δ
2
,where(23)
Δ
=
p
2
−
4
q
=
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
+
r
e
r
−
1
2
−
4
l
e
r
−
1
+
ζ
2
r
e
r
+
r
=
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
.Hereafter, we will give the topological classification of (8) about B
r
e
r
/
e
r
−
1
,
r according to the sign of Δ
=
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r.Lemma 4.
ForB
r
e
r
/
e
r
−
1
,
r, the following holds:
(i)B
r
e
r
/
e
r
−
1
,
r is Locally Asymptotically Node if
(24)
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
≥
0
,
0
<
ζ
2
<
e
r
−
1
2
e
r
r
−
r
+
1
e
r
−
1
;
(ii)B
r
e
r
/
e
r
−
1
,
r is Unstable Node if
(25)
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
≥
0
,
ζ
2
>
e
r
−
1
2
e
r
r
−
r
+
1
e
r
−
1
;
(iii)B
r
e
r
/
e
r
−
1
,
r is nonhyperbolic if
(26)
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
≥
0
,
ζ
2
=
e
r
−
1
2
e
r
r
−
r
+
1
e
r
−
1
.Lemma 5.
ForB
r
e
r
/
e
r
−
1
,
r, following statements holds:(i)B
r
e
r
/
e
r
−
1
,
r is Locally Asymptotically Focus if
(27)
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
<
0
,
0
<
ζ
2
<
e
r
1
−
r
−
1
r
e
r
r
−
1
;
(ii)B
r
e
r
/
e
r
−
1
,
r is Unstable Focus if
(28)
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
<
0
,
ζ
2
>
e
r
1
−
r
−
1
r
e
r
r
−
1
;
(iii)B
r
e
r
/
e
r
−
1
,
r is nonhyperbolic if
(29)
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
<
0
,
ζ
2
=
e
r
1
−
r
−
1
r
e
r
r
−
1
.If condition (iii) of Lemma5 holds then we obtain that eigenvalues of B
r
e
r
/
e
r
−
1
,
r are a pair of conjugate complex numbers with modulus one. So Neimark–Sacker bifurcation exists by the variation of parameter in a small neighborhood of ζ
2
=
e
r
1
−
r
−
1
/
r
e
r
r
−
1. For simplicity, we denote the parameters satisfying ζ
2
=
e
r
1
−
r
−
1
/
r
e
r
r
−
1 as(30)
N
B
=
ς
1
,
ς
2
:
e
r
−
1
e
r
−
1
+
ς
2
r
e
r
−
r
e
r
−
1
2
−
4
r
<
0
,
ς
2
=
e
r
1
−
r
−
1
r
e
r
r
−
1
,
0
<
r
<
r
∗
.


## 4. Bifurcation Analysis aboutB
r
e
r
/
e
r
−
1
,
r of the Model (8)

This Section deals with the study of Neimark–Sacker bifurcation of the model (8) about B
r
e
r
/
e
r
−
1
,
r. Consider parameter ζ
2 in a small neighborhood of ζ
2
∗, i.e., ζ
2
=
ζ
2
∗
+
ϵ, where ϵ
≪
1, then (8) becomes:(31)
x
n
+
1
=
ζ
1
x
n
e
−
y
n
1
+
ζ
2
∗
+
ϵ
x
n
y
n
+
1
=
x
n
1
−
e
−
y
n
.The characteristic equation ofJ
B
r
e
r
/
e
r
−
1
,
r about B
r
e
r
/
e
r
−
1
,
r of (31) is(32)
κ
2
−
q
1
ϵ
κ
+
q
2
ϵ
=
0
,where(33)
q
1
ϵ
=
e
r
−
1
e
r
−
1
+
ζ
2
∗
+
ϵ
r
e
r
+
r
e
r
−
1
,
q
2
ϵ
=
r
e
r
−
1
+
ζ
2
∗
+
ϵ
r
e
r
+
r
.The roots of characteristic equation ofJ
B
r
e
r
/
e
r
−
1
,
r about B
r
e
r
/
e
r
−
1
,
r are(34)
κ
1
,
2
=
q
1
ϵ
±
ι
4
q
2
ϵ
−
q
1
2
ϵ
2
=
e
r
−
1
2
e
r
−
1
+
ζ
2
∗
+
ϵ
r
e
r
+
r
2
e
r
−
1
±
ι
2
4
r
−
e
r
−
1
e
r
−
1
+
ζ
2
∗
+
ϵ
r
e
r
−
r
e
r
−
1
2
.
(35)
κ
1
,
2
=
q
2
ϵ
,
d
κ
1
,
2
d
ϵ
ϵ
=
0
=
−
1
2
e
r
r
−
1
2
<
0
.Additionally, we required that whenϵ
=
0, κ
1
,
2
m
≠
1, m
=
1
,
2
,
3
,
4, which corresponds to q
1
0
≠
−
2
,
0
,
1
,
2. Since q
1
0
2
−
4
q
2
0
<
0 and q
2
0
=
1. Thus q
1
0
2
<
4 and hence q
1
0
≠
±
2. So we only require that q
1
0
≠
0
,
1. By computation, we get(36)
r
e
r
−
1
≠
e
r
−
1
r
−
1
r
,
e
r
r
−
1
+
1
r
.Letu
n
=
x
n
−
x
∗
,
v
n
=
y
n
−
y
∗  then equilibrium B
r
e
r
/
e
r
−
1
,
r of system (8) transform into O
0
,
0. By calculating, we obtain(37)
u
n
+
1
=
ζ
1
u
n
+
x
∗
e
−
v
n
+
y
∗
1
+
ζ
2
∗
+
ϵ
u
n
+
x
∗
−
x
∗
v
n
+
1
=
u
n
+
x
∗
1
−
e
−
v
n
+
y
∗
−
y
∗
,wherex
∗
=
r
e
r
/
e
r
−
1
,
y
∗
=
r. Hereafter, when ϵ
=
0, normal form of system (37) is studied. Expanding (37) up to third order about u
n
,
v
n
=
0
,
0 by Taylor series, we get(38)
u
n
+
1
=
Λ
11
u
n
+
Λ
12
v
n
+
Λ
13
u
n
2
+
Λ
14
u
n
v
n
+
Λ
15
v
n
2
+
Λ
16
u
n
3
+
Λ
17
u
n
2
v
n
+
Λ
18
u
n
v
n
2
+
Λ
19
v
n
3
+
o
u
n
+
v
n
3
u
n
+
1
=
Λ
21
u
n
+
Λ
22
v
n
+
Λ
23
u
n
v
n
+
Λ
24
v
n
2
+
Λ
25
u
n
v
n
2
+
Λ
26
v
n
3
+
o
u
n
+
v
n
3
,where(39)
Λ
11
=
1
1
+
ζ
2
∗
x
∗
,
Λ
12
=
−
x
∗
,
Λ
13
=
−
ζ
2
∗
1
+
ζ
2
∗
x
∗
2
,
Λ
14
=
−
1
1
+
ζ
2
∗
x
∗
,
Λ
15
=
1
2
x
∗
,
Λ
16
=
ζ
2
∗
1
+
ζ
2
∗
x
∗
3
,
Λ
17
=
ζ
2
∗
1
+
ζ
2
∗
x
∗
2
,
Λ
18
=
1
2
1
+
ζ
2
∗
x
∗
,
Λ
19
=
−
1
6
x
∗
,
Λ
21
=
e
y
∗
−
1
e
y
∗
,
Λ
22
=
x
∗
e
y
∗
,
Λ
23
=
1
e
y
∗
,
Λ
24
=
−
1
2
x
∗
e
y
∗
,
Λ
25
=
−
1
2
1
e
y
∗
,
Λ
26
=
1
6
x
∗
e
y
∗
.Now, let(40)
η
=
e
r
−
1
2
e
r
−
1
+
ζ
2
∗
r
e
r
+
r
2
e
r
−
1
,
ζ
=
1
2
4
r
−
e
r
−
1
e
r
−
1
+
ζ
2
∗
r
e
r
−
r
e
r
−
1
2
,and invertible matrixT defined by(41)
T
=
Λ
12
0
η
−
Λ
11
−
ζ
.Using following translation(42)
u
n
v
n
=
Λ
12
0
η
−
Λ
11
−
ζ
X
^
n
Y
^
n
,(38) gives:(43)
X
n
+
1
^
Y
n
+
1
^
=
η
−
ζ
ζ
η
X
n
^
Y
n
^
Φ
^
X
n
^
,
Y
n
^
Ψ
^
X
n
^
,
Y
n
^
,where(44)
Φ
^
X
n
^
,
Y
n
^
=
Ω
11
X
n
^
2
+
Ω
12
X
n
^
Y
n
^
+
Ω
13
Y
n
^
2
+
Ω
14
X
n
^
3
+
Ω
15
X
n
^
2
Y
n
^
+
Ω
16
X
n
^
Y
n
^
2
+
Ω
17
Y
n
^
3
+
o
X
n
^
+
Y
n
^
3
,
Ψ
^
X
n
^
,
Y
n
^
=
Ω
21
X
n
^
2
+
Ω
22
X
n
^
Y
n
^
+
Ω
23
Y
n
^
2
+
Ω
24
X
n
^
3
+
Ω
25
X
n
^
2
Y
n
^
+
Ω
26
X
n
^
Y
n
^
2
+
Ω
27
Y
n
^
3
+
o
X
n
^
+
Y
n
^
3
,(45)
Ω
11
=
Λ
12
Λ
13
+
η
−
Λ
11
Λ
14
+
Λ
15
η
−
Λ
11
Λ
12
,
Ω
12
=
−
ζ
Λ
14
+
2
Λ
15
η
−
Λ
11
Λ
12
,
Ω
13
=
ζ
2
Λ
15
Λ
12
,
Ω
14
=
Λ
12
Λ
12
Λ
16
+
Λ
17
η
−
Λ
11
+
η
−
Λ
11
2
Λ
18
+
Λ
19
η
−
Λ
11
Λ
12
,
Ω
15
=
−
ζ
Λ
12
Λ
17
+
η
−
Λ
11
2
Λ
18
+
3
Λ
19
η
−
Λ
11
Λ
12
,
Ω
16
=
ζ
2
Λ
18
+
3
Λ
19
η
−
Λ
11
Λ
12
,
Ω
17
=
−
ζ
3
Λ
19
Λ
12
,
Ω
21
=
η
−
Λ
11
ζ
Λ
12
Λ
13
+
Λ
14
η
−
Λ
11
−
Λ
23
+
η
−
Λ
11
Λ
15
η
−
Λ
11
Λ
12
−
Λ
24
,
Ω
22
=
2
η
−
Λ
11
Λ
15
η
−
Λ
11
Λ
12
−
Λ
24
−
Λ
12
Λ
14
η
−
Λ
11
−
Λ
23
,
Ω
23
=
ζ
Λ
15
η
−
Λ
11
Λ
12
−
Λ
24
,
Ω
24
=
η
−
Λ
11
ζ
Λ
12
η
−
Λ
11
Λ
12
Λ
16
+
Λ
17
η
−
Λ
11
+
Λ
18
η
−
Λ
11
Λ
12
−
Λ
25
+
η
−
Λ
11
2
Λ
19
η
−
Λ
11
Λ
12
−
Λ
26
,
Ω
25
=
−
η
−
Λ
11
Λ
12
Λ
17
ζ
+
2
Λ
18
η
−
Λ
11
Λ
12
−
Λ
25
+
3
η
−
Λ
11
Λ
19
η
−
Λ
11
Λ
12
−
Λ
26
,
Ω
26
=
ζ
Λ
12
Λ
18
η
−
Λ
11
Λ
12
−
Λ
25
+
3
η
−
Λ
11
Λ
19
η
−
Λ
11
Λ
12
−
Λ
26
,
Ω
27
=
−
ζ
2
Λ
19
η
−
Λ
11
Λ
12
−
Λ
26
.In addition,(46)
Φ
^
X
n
^
X
n
^
0
,
0
=
2
Ω
11
,
Φ
^
X
n
^
Y
n
^
0
,
0
=
Ω
12
,
Φ
^
Y
n
^
Y
n
^
0
,
0
=
2
Ω
13
,
Φ
^
X
n
^
X
n
^
X
n
^
0
,
0
=
6
Ω
14
,
Φ
^
X
n
^
X
n
^
Y
n
^
0
,
0
=
2
Ω
15
,
Φ
^
X
n
^
Y
n
^
Y
n
^
0
,
0
=
2
Ω
16
,
Φ
^
Y
n
^
Y
n
^
Y
n
^
0
,
0
=
6
Ω
17
,
(47)
Ψ
^
X
n
^
X
n
^
0
,
0
=
2
Ω
21
,
Ψ
^
X
n
^
Y
n
^
0
,
0
=
Ω
22
,
Ψ
^
Y
n
^
Y
n
^
0
,
0
=
2
Ω
23
,
Ψ
^
X
n
^
X
n
^
X
n
^
0
,
0
=
6
Ω
24
,
Ψ
^
X
n
^
X
n
^
Y
n
^
0
,
0
=
2
Ω
25
,
Ψ
^
X
n
^
Y
n
^
Y
n
^
0
,
0
=
2
Ω
26
,
Ψ
^
Y
n
^
Y
n
^
Y
n
^
0
,
0
=
6
Ω
27
.In order for (43) to undergo Neimark–Sacker bifurcation, it is required that following discriminatory quantity, i.e., Ψ
≠
0 (see [6–13]).(48)
Ψ
=
−
R
e
1
−
2
κ
¯
κ
¯
2
1
−
κ
τ
11
τ
20
−
1
2
τ
11
2
−
τ
02
2
+
R
e
κ
¯
τ
21
,where(49)
τ
02
=
1
8
Φ
^
X
n
^
X
n
^
−
Φ
^
Y
n
^
Y
n
^
+
2
Ψ
^
X
n
^
Y
n
^
+
ι
Ψ
^
X
n
^
X
n
^
−
Ψ
^
Y
n
^
Y
n
^
+
2
Φ
^
X
n
^
Y
n
^
0
,
0
,
τ
11
=
1
4
Φ
^
X
n
^
X
n
^
+
Φ
^
Y
n
^
Y
n
^
+
ι
Ψ
X
n
^
X
n
^
+
Ψ
^
Y
n
^
Y
n
^
0
,
0
,
τ
20
=
1
8
Φ
^
X
n
^
X
n
^
−
Φ
^
Y
n
^
Y
n
^
+
2
Ψ
^
X
n
^
Y
n
^
+
ι
Ψ
^
X
n
^
X
n
^
−
Ψ
^
Y
n
^
Y
n
^
−
2
Φ
^
X
n
^
Y
n
^
0
,
0
,
τ
21
=
1
16
Φ
^
X
n
^
X
n
^
X
n
^
+
Φ
^
X
n
^
Y
n
^
Y
n
^
+
Ψ
^
X
n
^
X
n
^
Y
n
^
+
Ψ
^
Y
n
^
Y
n
^
Y
n
^
+
ι
Ψ
^
X
n
^
X
n
^
X
n
^
+
Ψ
^
X
n
^
Y
n
^
Y
n
^
−
Φ
^
X
n
^
X
n
^
Y
n
^
−
Φ
^
Y
n
^
Y
n
^
Y
n
^
0
,
0
.After calculating, we get(50)
τ
02
=
1
4
Ω
11
−
Ω
13
+
Ω
22
+
ι
Ω
21
−
Ω
23
+
Ω
12
,
τ
11
=
1
2
Ω
11
+
Ω
13
+
ι
Ω
21
+
Ω
23
,
τ
20
=
1
4
Ω
11
−
Ω
13
+
Ω
22
+
ι
Ω
21
−
Ω
23
−
Ω
12
,
τ
21
=
1
8
3
Ω
14
+
Ω
16
+
Ω
25
+
3
Ω
27
+
ι
3
Ω
24
+
Ω
26
−
Ω
15
−
3
Ω
17
.Based on this analysis and Neimark–Sacker bifurcation Theorem discussed in [12, 13], we arrive at the following Theorem:Theorem 1.
IfΨ
≠
0 then model (8) undergoes a Neimark–Sacker bifurcation about B
r
e
r
/
e
r
−
1
,
r as the parameters ζ
1
,
ζ
2 go through N
B.  Additionally, attracting (respectively repelling) invariant closed curve bifurcate from B
r
e
r
/
e
r
−
1
,
r if Ψ
<
0 (respectively Ψ
>
0).According to Neimark–Sacker bifurcation discussed in [12, 13], the bifurcation is called supercritical Neimark–Sacker bifurcation if the discriminatory quantity Ψ
<
0. In the following Section, numerical simulations guarantee that supercritical Neimark–Sacker bifurcation occurs for the model (8). Biologically, attracting closed curve indicates that both parasitoid and host populations will coexist under the periodic or quasi-periodic oscillations with long time.


## 5. Numerical Simulations and Discussion

This work deals with the study of local dynamics and bifurcation analysis of a discrete-time two-species model inR
+
2. We proved that the model has two boundary equilibria: O
,
A
ζ
1
−
1
/
ζ
2
,
0 and a unique positive equilibrium point B
r
e
r
/
e
r
−
1
,
r under certain parametric conditions. We studied the local dynamics along with topological classification about equilibria: O
0
,
0, A
ζ
1
−
1
/
ζ
2
,
0, B
r
e
r
/
e
r
−
1
,
r
, and conclusion is presented in Table 1. We proved that about A
ζ
1
−
1
/
ζ
2
,
0 there may exist a fold bifurcation when parameters of (8) are located in the set: F
A
=
ζ
1
,
ζ
2
:
ζ
2
=
ζ
1
−
1
,
ζ
1
,
ζ
2
>
0. We also proved that if  ζ
2
=
e
r
1
−
r
−
1
/
r
e
r
r
−
1  then eigenvalues J
B
r
e
r
/
e
r
−
1
,
r about B
r
e
r
/
e
r
−
1
,
r are pair of complex conjugate with modulus one and thus in particular supercritical Neimark–Sacker bifurcation occurs under the bifurcation curve:(51)
N
B
=
ζ
1
,
ζ
2
:
e
r
−
1
e
r
−
1
+
ζ
2
r
e
r
−
r
e
r
−
1
2
−
4
r
<
0
,
ζ
2
=
e
r
1
−
r
−
1
r
e
r
r
−
1
,
0
<
r
<
r
∗
.Table 1
Number of equilibria along their qualitative behavior of (8).
E.P
Corresponding behavior
O
0
,
0
Sink ifζ
1
<
1; never source; saddle if ζ
1
>
1; nonhyperbolic if ζ
1
=
1.
A
ζ
1
−
1
/
ζ
2
,
0
Sink ifζ
2
>
ζ
1
−
1; never source; saddle if ζ
2
<
ζ
1
−
1; nonhyperbolic if ζ
2
=
ζ
1
−
1.
B
r
e
r
/
e
r
−
1
,
r
Locally asymptotically node if
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r
≥
0
and
0
<
ζ
2
<
(
e
r
−
1
)
2
/
e
r
r
−
r
+
1
e
r
−
1
;
Unstable node if
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r
≥
0
and
ζ
2
>
(
e
r
−
1
)
2
/
e
r
r
−
r
+
1
e
r
−
1
;
Nonhyperbolic (for real eigenvalues) if
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r
≥
0
and
ζ
2
=
(
e
r
−
1
)
2
/
e
r
r
−
r
+
1
e
r
−
1
;
Locally asymptotically focus if
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r
<
0
and
0
<
ζ
2
<
e
r
1
−
r
−
1
/
r
e
r
r
−
1
;
Unstable focus if
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r
<
0
and
ζ
2
>
e
r
1
−
r
−
1
/
r
e
r
r
−
1
;
Nonhyperbolic (for complex eigenvalues) if
e
r
−
1
/
e
r
−
1
+
ζ
2
r
e
r
−
r
/
e
r
−
1
2
−
4
r
<
0
and
ζ
2
=
e
r
1
−
r
−
1
/
r
e
r
r
−
1
.Biologically, existence of stable closed curves implies that there exist the periodic or quasiperiodic oscillations between host and parasitoid populations. Finally, numerical simulations are provided to verify theoretical discussion. These numerical simulations presented in Figures2–5 agree with our theoretical discussion. Figure 2 shows that B
r
e
r
/
e
r
−
1
,
r of the model (8) is Locally Asymptotically Focus when ζ
2
<
ζ
2
∗, where ζ
2
∗
≈
0.5000046770 as presented in Figures 2(a)–2(i) by choosing ζ
1
=
2.9899. But when ζ
2 goes through the bifurcation value ζ
2
∗, equilibrium B
r
e
r
/
e
r
−
1
,
r of (8) is Unstable Focus. Meanwhile, an attracting closed invariant curve bifurcates from B
r
e
r
/
e
r
−
1
,
r of the model (8) as presented in Figures 3(a)–3(i). Moreover, bifurcation diagrams along with Maximum Lyapunov Exponent in this case, are plotted and drawn in Figure 4. Finally 3D bifurcation diagrams are also plotted and drawn in Figure 5.Figure 2
Phase portraits of the model (8). (a) ζ
2
=
0.01665 with (0.006, 0.04). (b) ζ
2
=
0.042365 with (0.006, 0.04). (c) ζ
2
=
0.235 with (0.006, 0.004). (d) ζ
2
=
0.477 with (0.06, 0.04). (e) ζ
2
=
0.4789 with (0.06, 0.04). (f) ζ
2
=
0.4789235 with (0.06, 0.04). (g) ζ
2
=
0.466 with (0.06, 0.04). (h) ζ
2
=
0.4332 with (0.06, 0.04). (i) ζ
2
=
0.45545 with (0.06, 0.04).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 3
Phase portraits of the model (8). (a) ζ
2
=
2.9 with (0.1, 0.02). (b) ζ
2
=
2.989 with (0.1, 0.2). (c) ζ
2
=
2.9899 with (0.1, 0.2). (d) ζ
2
=
3.6 with (0.001, 0.002). (e) ζ
2
=
3.754 with (0.001, 0.002). (f) ζ
2
=
4.1 with (0.01, 0.02). (g) ζ
2
=
4.231 with (1.001, 0.2). (h) ζ
2
=
4.5 with (0.01, 0.0000002). (i) ζ
2
=
5.9 with (0.01, 0.002).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 4
Bifurcation diagram and their corresponding maximum lyapunov exponent of the model (8) about B
r
e
r
/
e
r
−
1
,
r. (a)-(b) Bifurcation diagram of the model if ζ
2
∈
0.32
,
1.85 and initial condition 1.1
,
1.2. (c) Maximum lyapunov exponent corresponding to (a)-(b).
(a)
(b)
(c)Figure 5
3D bifurcation diagrams of the model (8).
(a)
(b)
(c)



---

*Source: 2954059-2020-01-09.xml*