# Bifurcation Analysis of a Delayed Predator-Prey Model with Holling Type III Functional Response and Predator Harvesting

Authors: Uttam Das; T. K. Kar
Journal: Journal of Nonlinear Dynamics
 (2014)

---

## Abstract

This paper tries to highlight  a delayed prey-predator model with Holling type III functional response and harvesting to predator species. In this context, we have discussed local stability of the equilibria, and the occurrence of Hopf bifurcation of the system is examined by considering the harvesting effort as bifurcation parameter along with the influences of harvesting effort of the system when time delay is zero. Direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also studied by applying the normal form theory and the center manifold theorem. Lastly some numerical simulations are carried out to draw for the validity of the theoretical results.


---

## Body

## 1. Introduction and Model Description

Differential equation models for interactions between species are one of the classical applications of mathematics to biology, dating back to the first half of this century. The development and use of analytical techniques and the growth of computer power have progressively improved our understanding of these types of models.The study of population dynamics with harvesting is a subject of mathematical bioeconomics, which in turn related to the optimal management of renewable resources, Clark [1]. Generally the concept of optimal resource management is based on the standard cost benefit criterion which maximizes present values of net economic revenues. This criterion is relevant to both private and public management decisions, although the specification of costs and benefits are not necessarily the same in both cases.Regulation of exploitation of biological resources has become a problem of major concern nowadays in view of the dwindling resource stocks and the deteriorating environment. Exploitation reduces the biomass of the concerned species, exhibits oscillation, and even causes extinction of some other species. Rosenweig-MacAurtho model experiences oscillation under selective effort of Kar and Ghosh [2]. Legovic et al. [3] have concluded that harvesting the prey species at maximum level causes the extinction of the predator species in traditional prey-predator system. Kar and Ghosh [4] and Ghosh and Kar [5] show that harvesting the prey species at maximum sustainable yield (MSY) level never causes the extinction of the predator species in both ratio-dependent and Holling-Tanner prey-predator systems. It is shown that harvesting the prey species at MSY level may or may not drive the predator population to extinction if intraspecific competition is present among the predator species, Kar and Ghosh [4]. More recently, Ghosh and Kar [6] somehow managed an example to show that prey harvesting at MSY level may be sustainable in tritrophic food chain model. In conclusion, it is reasonable that harvesting the prey species to achieve MSY will cause the extinction of the predator if predation process follows the Holling type functional response. As MSY policy has a severe impact on other tropic levels in ecosystem, suitable management tools are required for smooth running of sustainable fishing. Creating marine protected areas, maximizing the economic rent, and taxation are some good examples to protect fish population from extinction, over exploitation, as well as to achieve desired stock levels for future generation. Kar and Matsuda [7] and Kar and Ghosh [8] have shown that marine protection can have positive impact from biological stock assessment but virtually reduces the economic rent of the fishermen. Recently, Kar and Das [9] have analyzed a prey-predator model with harvesting where they introduce a tax to regulate the fishery.It has been noticed over the years by many authors that time delays have come to play an important role in almost all branches of science, for example, ecology and biology. The time delay is considered in the population dynamics when the rate of change of the population is not only a function of the present population but also depends on the past population. Delay is frequently used in a predator-prey model to represent the biological process more accurately.Wangersky and Cunningham [10] have considered a predator-prey system with time delay and studied the effect of time delay on the system. Martin and Ruan [11] have studied the combined effects of prey harvesting and delay on the dynamics of predator-prey systems and focused on three very well-studied delay predator-prey models. For first and third models they were shown that the time delay could cause not only instability and oscillations but also the switching of stabilities, while the prey harvesting changes only the equilibrium values but not the properties of solutions. But in the second model, the time delay induces instability and bifurcation but there is no switching of stabilities; however, increasing the prey harvesting level will help the system to regain its stability. This indicates that the prey harvesting has a stabilizing effect on the dynamics of the model. Chen and Changming [12] have considered a delayed model with stage structured for prey and investigated nonnegative equilibria, existence of Hopf bifurcation, the stability and direction of Hopf bifurcation by applying the normal from theory, and the centre manifold argument. Zhang et al. [13] have considered a delayed predator-prey model with constant rate prey harvesting and Holling type III functional response. They studied local stability of equilibria and existence of Hopf bifurcation at interior equilibrium. Direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also studied by applying the normal form theory and the center manifold theorem. Also Kar [14], Kar and Matsuda [15], Kar and Pahari [16, 17], Kar and Ghorai [18], Kuang [19], and some other authors have discussed delayed predator-prey system.The purpose of the work is to illustrate the combined effects of harvesting and delay on the dynamics of predator-prey system.In this paper, we consider the following prey-predator system:(1)
d
x
d
t
=
r
x
(
1
-
x
k
)
-
m
x
2
y
a
+
x
2
,
d
y
d
t
=
-
d
y
+
m
β
x
2
(
t
-
τ
)
y
(
t
-
τ
)
a
+
x
2
(
t
-
τ
)
-
q
E
y
with initial conditions x
(
θ
)
=
ϕ
(
θ
), y
(
θ
)
=
ψ
(
θ
), ϕ
(
θ
)
≥
0, ψ
(
θ
)
≥
0, θ
∈
[
-
τ
,
0
], ϕ
(
0
)
>
0, and ψ
(
0
)
>
0. Here x and y denote the prey and predator populations, respectively, at time t. r and k, respectively, represent the intrinsic growth rate and carrying capacity of the prey. m is the maximum capture rate of the prey, d is the death rate of the predator, and β is the conversion factor. The term m
x
2
/
(
a
+
x
2
) denotes the functional response of the predator, which is termed as Holling type III response function (see [20]). We assume that the predator population is harvested according to catch-per-unit effort (CPUE) hypothesis [1] which describes that catch-per-unit effort is proportional to the stock level. Therefore, the harvest at time t is h
t
=
q
E
y, where q is the catchability coefficient and E is the harvesting effort. The time delay τ
≥
0 is a constant based on the assumption that the change rate of the predators depends on the number of prey and of predators present at some previous time. r, β, k, d, a, c, m, q, and E are positive constants.Forτ
=
0, system (1) becomes
(2)
d
x
d
t
=
r
x
(
1
-
x
k
)
-
m
x
2
y
a
+
x
2
,
d
y
d
t
=
-
d
y
+
m
β
x
2
y
a
+
x
2
-
q
E
y
.Our paper is organized in the following way: existence of equilibria of the system (2), the local stability criterion, and uniqueness of limit cycles at interior equilibrium are discussed in Section 2. Influence of harvesting effort E is discussed in Section 3. In presence of delay we studied stability and Hopf bifurcation in Section 4. Direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also studied by applying the normal form theory and the center manifold theorem in Section 5. Numerical simulations are given in Section 6. A brief concluding remark is given in Section 7.


## 2. The Steady States and Their Stability

We now study the existence and nature of the steady states. Particularly we are interested in the interior equilibrium of the system. To begin with, we list all possible steady states of the system (2) as follows: P
0
(
0,0
), P
1
(
k
,
0
), and the interior equilibrium P
*
(
x
*
,
y
*
), where
(3)
x
*
=
(
q
a
E
+
a
d
m
β
-
d
-
q
E
)
1
/
2
,
y
*
=
r
m
k
x
*
(
k
-
x
*
)
(
a
+
(
x
*
)
2
)
.
If m
β
>
d
+
q
E holds, then the interior equilibrium P
*
(
x
*
,
y
*
) exists.For the stability analysis of the equilibria, let the variational matrix of the system (2) at an arbitrary point (
x
,
y
) be M
(
x
,
y
).AtP
0
(
0,0
), the eigenvalues of the variational matrix M
(
0,0
) are r and -
(
d
+
q
E
). Therefore, P
0
(
0,0
) is unstable.AtP
1
(
k
,
0
), the eigenvalues of the variational matrix M
(
k
,
0
) are -
r and -
(
d
+
q
E
). Therefore, P
1
(
k
,
0
) is asymptotically stable.The characteristic equation of the variational matrixM
(
x
*
,
y
*
) is given by
(4)
γ
2
+
γ
M
1
+
M
2
=
0
,
where
(5)
M
1
=
2
a
m
x
*
y
*
(
a
+
(
x
*
)
2
)
2
+
2
x
*
r
k
-
r
,
M
2
=
2
m
2
a
β
(
x
*
)
3
y
*
(
a
+
(
x
*
)
2
)
3
>
0
.
Routh-Hurwitz criterion states that all roots of the characteristic equation (4) have negative real parts if M
1
>
0.Therefore, interior equilibriumP
*
(
x
*
,
y
*
) will be asymptotically stable if M
1
>
0 and unstable if M
1
<
0.Theorem 1.
IfP
*
(
x
*
,
y
*
) exists with M
1
>
0, then P
*
(
x
*
,
y
*
) is locally asymptotically stable.


### 2.1. Uniqueness of Limit Cycles

It is known that, for prey-predator systems, existence and stability of a limit cycle are related to the existence and stability of a positive equilibrium. If the limit cycles do not exist, in this case the equilibrium is globally asymptotically stable. On the other hand if the positive equilibrium exists and is unstable, there must occur at least one limit cycle.Let us consider system (2) in the form
(6)
d
x
d
t
=
x
g
(
x
)
-
y
h
(
x
)
,
x
(
0
)
>
0
,
d
y
d
t
=
y
(
-
d
-
q
E
+
q
(
x
)
)
,
y
(
0
)
>
0
,
where g
(
x
)
=
r
(
1
-
(
x
/
k
)
), h
(
x
)
=
m
x
2
/
(
a
+
x
2
), and q
(
x
)
=
m
β
x
2
/
(
a
+
x
2
).Now we consider the following Theorem 4.2 of Kuang and Freedman [21] regarding uniqueness of limit cycles of the above system.Theorem 2.
Suppose for system (6)
(7)
d
d
x
(
x
g
′
(
x
)
+
g
(
x
)
-
x
g
(
x
)
(
h
′
(
x
)
/
h
(
x
)
)
-
d
-
q
E
+
q
(
x
)
)
≤
0
in 0
≤
x
<
x
* and x
*
<
x
≤
k. Then the above system has exactly one limit cycles which is globally asymptotically stable with respect to the set
(8)
{
(
x
,
y
)
:
x
>
0
,
y
>
0
}
∖
{
P
*
(
x
*
,
y
*
)
}
.Following Theorem2, we may state that when M
1
<
0, then system (6) has unique globally stable limit cycle. Thus we see that when the system is unstable, there exists unique globally stable limit cycle.


### 2.2. Bifurcation Analysis for the ParameterE

The characteristic equation (4) has two purely imaginary roots if and only if M
1
=
0 and M
2
>
0 for some values of E, say E
*. Therefore, there is only one value of E
* at which we have a Hopf bifurcation. Since x
* and y
* contain E, it is not easy to find E explicitly. So we are trying to find E
=
E
* numerically by solving the system of equations with suitable values of the parameters:
(9)
r
(
1
-
x
*
k
)
-
m
x
*
y
*
(
a
+
(
x
*
)
2
)
=
0
,
-
d
+
m
β
(
x
*
)
2
a
+
(
x
*
)
2
-
q
E
=
0
,
2
a
m
x
*
y
*
(
a
+
(
x
*
)
2
)
2
+
2
x
*
r
k
-
r
=
0
.
Therefore the system undergoes a Hopf bifurcation at E
=
E
*. The system is stable for E
>
E
* and unstable E
<
E
*.


## 2.1. Uniqueness of Limit Cycles

It is known that, for prey-predator systems, existence and stability of a limit cycle are related to the existence and stability of a positive equilibrium. If the limit cycles do not exist, in this case the equilibrium is globally asymptotically stable. On the other hand if the positive equilibrium exists and is unstable, there must occur at least one limit cycle.Let us consider system (2) in the form
(6)
d
x
d
t
=
x
g
(
x
)
-
y
h
(
x
)
,
x
(
0
)
>
0
,
d
y
d
t
=
y
(
-
d
-
q
E
+
q
(
x
)
)
,
y
(
0
)
>
0
,
where g
(
x
)
=
r
(
1
-
(
x
/
k
)
), h
(
x
)
=
m
x
2
/
(
a
+
x
2
), and q
(
x
)
=
m
β
x
2
/
(
a
+
x
2
).Now we consider the following Theorem 4.2 of Kuang and Freedman [21] regarding uniqueness of limit cycles of the above system.Theorem 2.
Suppose for system (6)
(7)
d
d
x
(
x
g
′
(
x
)
+
g
(
x
)
-
x
g
(
x
)
(
h
′
(
x
)
/
h
(
x
)
)
-
d
-
q
E
+
q
(
x
)
)
≤
0
in 0
≤
x
<
x
* and x
*
<
x
≤
k. Then the above system has exactly one limit cycles which is globally asymptotically stable with respect to the set
(8)
{
(
x
,
y
)
:
x
>
0
,
y
>
0
}
∖
{
P
*
(
x
*
,
y
*
)
}
.Following Theorem2, we may state that when M
1
<
0, then system (6) has unique globally stable limit cycle. Thus we see that when the system is unstable, there exists unique globally stable limit cycle.


## 2.2. Bifurcation Analysis for the ParameterE

The characteristic equation (4) has two purely imaginary roots if and only if M
1
=
0 and M
2
>
0 for some values of E, say E
*. Therefore, there is only one value of E
* at which we have a Hopf bifurcation. Since x
* and y
* contain E, it is not easy to find E explicitly. So we are trying to find E
=
E
* numerically by solving the system of equations with suitable values of the parameters:
(9)
r
(
1
-
x
*
k
)
-
m
x
*
y
*
(
a
+
(
x
*
)
2
)
=
0
,
-
d
+
m
β
(
x
*
)
2
a
+
(
x
*
)
2
-
q
E
=
0
,
2
a
m
x
*
y
*
(
a
+
(
x
*
)
2
)
2
+
2
x
*
r
k
-
r
=
0
.
Therefore the system undergoes a Hopf bifurcation at E
=
E
*. The system is stable for E
>
E
* and unstable E
<
E
*.


## 3. Influence of the Harvesting EffortE

Let the interior equilibrium of the system (2) be P
*
(
x
*
,
y
*
), where
(10)
x
*
=
(
q
a
E
+
a
d
m
β
-
d
-
q
E
)
1
/
2
,
y
*
=
r
m
k
x
*
(
k
-
x
*
)
(
a
+
(
x
*
)
2
)
.
Now,
(11)
d
x
*
d
E
=
β
m
a
q
2
(
β
m
-
d
-
E
q
)
3
/
2
(
q
a
E
+
a
d
)
1
/
2
>
0
.
Therefore x
* is strictly increasing function of E.Also(12)
d
y
*
d
E
=
r
k
m
[
k
-
2
x
*
-
k
a
(
x
*
)
2
]
d
x
*
d
E
.
Since d
x
*
/
d
E
>
0, therefore
(13)
d
y
*
d
E
>
0
if
k
>
2
x
*
+
k
a
(
x
*
)
2
,
d
y
*
d
E
<
0
if
k
<
2
x
*
+
k
a
(
x
*
)
2
.
Thus as the harvesting effort E increases, the predator population increases when E is smaller than the threshold value. But if the harvesting effort E gradually increases above the threshold value, that is, as the harvesting effort E becomes large enough, then the increasing amount of the harvesting effort E can decreases predator population.To construct Figure1, we take r
=
0.8, a
=
1, β
=
0.8, m
=
0.9, d
=
0.3, k
=
18, and q
=
0.5 in appropriate units. From the above analysis we see that, as the harvesting effort E  increases, the prey and predator species increase. But if the harvesting effort E increases larger than a threshold value, that is, as the harvesting effort E becomes large enough, then increasing amount of the harvesting effort E can increase the prey species but decrease predator species and go to extinction of the predator species when E is large.Figure 1
Variation ofx
* and y
* with the change of the harvesting effort E.


## 4. Stability and Hopf Bifurcation in the Presence of Delay

We will investigate the dynamics of delay system (1).Let(
x
*
,
y
*
) be the only interior equilibrium of the system (1), and let X
=
x
-
x
*, Y
=
y
-
y
* be the perturbed variables. After removing the nonlinear terms, we obtain the linear variational system, by using equilibrium conditions as
(14)
d
X
d
t
=
a
1
X
+
a
2
Y
,
d
Y
d
t
=
a
3
X
(
t
-
τ
)
+
a
4
Y
(
t
-
τ
)
+
a
5
Y
,
where
(15)
a
1
=
r
-
2
x
*
r
k
-
2
a
m
x
*
y
*
(
a
+
(
x
*
)
2
)
2
,
a
2
=
-
m
(
x
*
)
2
a
+
(
x
*
)
2
,
a
3
=
2
m
β
a
x
*
y
*
(
a
+
(
x
*
)
2
)
2
,
a
4
=
m
β
(
x
*
)
2
a
+
(
x
*
)
2
,
a
5
=
-
(
d
+
q
E
)
.
The characteristic equation of the linearized system (14) is the transcendental equation of the following form:
(16)
λ
2
+
P
λ
+
Q
+
(
R
λ
+
S
)
e
-
λ
τ
=
0
,
where P
=
-
(
a
1
+
a
5
), Q
=
a
1
a
5, R
=
-
a
4, and S
=
a
1
a
4
-
a
2
a
3.It is known that the steady state is asymptotically stable if all roots of the characteristic equation (16) have negative real parts. For a nonlinear delay equation, there are two types of stabilities: absolute stability (independent of the delay) and conditional stability (depending of the delay).Case 1 (τ
=
0).
Equation (16) becomes
(17)
λ
2
+
(
P
+
R
)
λ
+
Q
+
S
=
0
.
Assume that all roots of (16) have negative real parts, which is true if (
H
1
)
P
+
R
>
0,
(
H
2
)
Q
+
S
>
0.Therefore, the interior equilibrium (
x
*
,
y
*
) is locally asymptotically stable if (
H
1
) and (
H
2
) hold.Case 2 (τ
≠
0).
We want to determine if the real part of some root increases to reach zero and eventually becomes positive asτ varies. If i
ω is a root of (16), then by separating the real and imaginary parts, we get
(18)
S
cos
⁡
ω
τ
+
R
ω
sin
ω
τ
=
ω
2
-
Q
,
(19)
R
ω
cos
⁡
ω
τ
-
S
sin
ω
τ
=
-
P
ω
.
From (19), we obtain the fourth-order equation for ω as
(20)
ω
4
-
(
R
2
-
P
2
+
2
Q
)
ω
2
+
(
Q
2
-
S
2
)
=
0
.
From (19) we also have cos
⁡
ω
τ
=
(
S
(
ω
2
-
Q
)
-
ω
2
P
R
)
/
(
ω
2
R
2
+
S
2
) and sin
ω
τ
=
(
ω
R
(
ω
2
-
Q
)
+
ω
P
S
)
/
(
ω
2
R
2
+
S
2
).
From (20), it follows that, if (
H
3
)
(
P
2
-
R
2
-
2
Q
)
>
0 and (
Q
2
-
S
2
)
>
0,then (20) does not have positive roots. Therefore, characteristic equation (16) does not have purely imaginary roots. Since (
H
1
) and (
H
2
) ensure that all roots of (17) have negative real parts by Rouche’s theorem, it follows that the roots of (20) have negative real part too. This is summarized as follows.Theorem 3.
If(
H
1
)
–
(
H
3
) hold, then all roots of (16) have negative real parts for all τ
≥
0.On the other hand, if(
H
4
)
(
Q
2
-
S
2
)
<
0,then (20) has a unique positive root ω
0
2.Substitutingω
0
2 into (19) and solving for τ, we get
(21)
τ
n
=
1
ω
0
arctan
(
ω
0
P
S
+
ω
0
R
(
ω
0
2
-
Q
)
S
(
ω
0
2
-
Q
)
-
ω
0
2
P
R
)
+
2
n
π
ω
0
,
n
=
0
,
1
,
2,3
,
…
.
If (
H
5
)
(
R
2
-
P
2
+
2
Q
)
>
0, (
Q
2
-
S
2
)
>
0, and (
R
2
-
P
2
+
2
Q
)
2
>
4
(
Q
2
-
S
2
), then (20) has two positive roots ω
+
2 and ω
-
2. Substituting ω
±
2 into (19), we get
(22)
τ
k
±
=
1
ω
±
arctan
(
ω
±
P
S
+
ω
±
R
(
ω
±
2
-
Q
)
S
(
ω
±
2
-
Q
)
-
ω
±
2
P
R
)
+
2
k
π
ω
±
,
k
=
0,1
,
2
,
3
,
…
.Differentiating (16) with respect to τ, we obtain
(23)
d
λ
d
τ
=
(
R
λ
2
+
S
λ
)
e
-
λ
τ
(
2
λ
+
P
)
+
(
R
-
R
λ
τ
-
S
τ
)
e
-
λ
τ
,
(
d
λ
d
τ
)
-
1
=
R
λ
(
R
λ
+
S
)
-
(
2
λ
+
P
)
λ
(
λ
2
+
P
λ
+
Q
)
-
τ
λ
,
by using e
-
λ
τ
=
-
(
(
λ
2
+
P
λ
+
Q
)
/
(
R
λ
+
S
)
).Thus(24)
sign
⁡
{
d
d
τ
(
Re
λ
)
}
λ
=
i
ω
=
sign
⁡
{
Re
(
d
λ
d
τ
)
-
1
}
λ
=
i
ω
=
sign
⁡
{
2
ω
2
-
(
R
2
-
P
2
+
2
Q
)
R
2
ω
2
+
S
2
}
,
by using (
ω
2
-
Q
)
2
+
P
2
ω
2
=
S
2
+
R
2
ω
2.Theorem 4.
If(
H
1
), (
H
2
), and (
H
4
) hold, then equilibrium (
x
*
,
y
*
) is asymptotically stable for τ
<
τ
0 and unstable τ
>
τ
0. Further, as τ increases through τ
0, (
x
*
,
y
*
) bifurcates into small amplitude periodic solutions, where
(25)
τ
0
=
1
ω
0
arctan
(
ω
0
P
S
+
ω
0
R
(
ω
0
2
-
Q
)
S
(
ω
0
2
-
Q
)
-
ω
0
2
P
R
)
.Proof.
Forτ
=
0, (
x
*
,
y
*
) is asymptotically stable if conditions (
H
1
) and (
H
2
) hold. Hence, by Butler’s lemma, (
x
*
,
y
*
) remain stable for τ
<
τ
0. We now have to show that
(26)
d
d
τ
(
Re
λ
)
|
τ
=
τ
0
,
ω
=
ω
0
>
0
.
This will signify that there exists at least one eigenvalue with positive real part forτ
>
τ
0. Moreover, the conditions of Hopf bifurcation (see [22]) are then satisfied yielding the required periodic solution.
From (24), we have
(27)
sign
⁡
{
d
d
τ
(
Re
λ
)
}
λ
=
i
ω
0
=
sign
⁡
{
(
R
2
-
P
2
+
2
Q
)
2
-
4
(
Q
2
-
S
2
)
R
2
ω
0
2
+
S
2
}
.
Therefore,
(28)
d
d
τ
(
Re
λ
)
|
τ
=
τ
0
,
ω
=
ω
0
>
0
.
Therefore, the transversality condition holds, and hence, Hopf bifurcation occurs atω
=
ω
0, τ
=
τ
0. This completes the proof.Theorem 5.
Letτ
k
± be defined in (22). If (
H
1
), (
H
2
), and (
H
5
) hold, then there is a positive integer m such that there are m switches from stability to instability and to stability. In other words, when τ
∈
[
0
,
τ
0
+
)
,
(
τ
0
-
,
τ
1
+
)
,
…
,
(
τ
m
-
1
-
,
τ
m
+
), the equilibrium is stable, and when τ
∈
[
τ
0
+
,
τ
0
-
)
,
(
τ
1
+
,
τ
1
-
)
,
…
,
(
τ
m
-
1
+
,
τ
m
-
), the equilibrium is unstable. Therefore, there are bifurcations at (
x
*
,
y
*
) when τ
=
τ
k
±, k
=
0
,
1
,
2
,
3
,
….Proof.
If conditions(
H
1
), (
H
2
), and (
H
5
) hold, then to prove the theorem, we need only to verify the transversality conditions by Cushing [23] and Cushing and Saleem [24]
(29)
d
d
τ
(
Re
λ
)
|
τ
=
τ
k
+
,
ω
=
ω
+
>
0
,
d
d
τ
(
Re
λ
)
|
τ
=
τ
k
-
,
ω
=
ω
-
<
0
.
From (24), we have
(30)
sign
⁡
{
d
d
τ
(
Re
λ
)
}
λ
=
i
ω
+
=
sign
⁡
{
(
R
2
-
P
2
+
2
Q
)
2
-
4
(
Q
2
-
S
2
)
R
2
ω
+
2
+
S
2
}
.
Therefore,
(31)
d
d
τ
(
Re
λ
)
|
τ
=
τ
k
+
,
ω
=
ω
+
>
0
.
Again,
(32)
sign
⁡
{
d
d
τ
(
Re
λ
)
}
λ
=
i
ω
-
=
sign
⁡
{
-
(
R
2
-
P
2
+
2
Q
)
2
-
4
(
Q
2
-
S
2
)
R
2
ω
-
2
+
S
2
}
.
Therefore,
(33)
d
d
τ
(
Re
λ
)
|
τ
=
τ
k
-
,
ω
=
ω
-
<
0
.
Hence, the transversality conditions hold. This completes the proof.


## 5. Stability and Direction of Hopf Bifurcation

Without loss of generality, denote the critical valuesτ
j by τ
~, and set τ
=
τ
~
+
μ. Then μ
=
0 is a Hopf bifurcation value of the system (1). Thus, we can work in the phase space C
=
C
(
[
-
τ
~
,
0
]
,
ℜ
2
). Let u
1
(
t
)
=
x
(
t
)
-
x
*, u
2
(
t
)
=
y
(
t
)
-
y
*.Then system (1) is transformed into
(34)
u
˙
1
(
t
)
=
a
1
u
1
(
t
)
+
a
2
u
2
(
t
)
+
∑
i
+
j
≥
2
1
i
!
j
!
f
i
j
(
1
)
u
1
i
(
t
)
u
2
j
(
t
)
,
u
˙
2
(
t
)
=
a
3
u
1
(
t
-
τ
)
+
a
4
u
2
(
t
-
τ
)
+
a
5
u
2
(
t
)
∑
i
+
j
+
l
≥
2
1
i
!
j
!
l
!
f
i
j
l
(
2
)
u
1
i
(
t
-
τ
)
×
u
2
j
(
t
-
τ
)
u
2
l
(
t
)
,
where
(35)
f
i
j
(
1
)
=
∂
i
+
j
f
(
1
)
∂
x
i
∂
y
i
|
(
x
*
,
y
*
)
,
f
i
j
l
(
2
)
=
∂
i
+
j
+
l
f
(
2
)
∂
x
i
(
t
-
τ
)
∂
y
j
(
t
-
τ
)
∂
y
l
|
(
x
*
,
y
*
,
y
*
)
,
i
,
j
,
l
≥
0
,
f
(
1
)
=
r
x
(
t
)
(
1
-
x
(
t
)
k
)
-
m
x
2
(
t
)
y
(
t
)
a
+
x
2
(
t
)
,
f
(
2
)
=
-
d
y
(
t
)
+
m
β
x
2
(
t
-
τ
)
y
(
t
-
τ
)
a
+
x
2
(
t
-
τ
)
-
q
E
y
(
t
)
.
We rewrite (34) as
(36)
u
˙
(
t
)
=
L
μ
(
u
t
)
+
f
(
μ
,
u
t
)
,
where u
(
t
)
=
(
u
1
(
t
)
,
u
2
(
t
)
)
T
∈
R
2, u
t
(
θ
)
∈
C is defined by u
t
(
θ
)
=
u
(
t
+
θ
) and L
μ
:
C
→
R, f
:
R
×
C
→
R are given by
(37)
L
μ
(
ϕ
)
=
(
a
1
a
2
0
a
5
)
ϕ
(
0
)
+
(
0
0
a
3
a
4
)
ϕ
(
-
τ
~
)
,
f
(
μ
,
ϕ
)
=
(
∑
i
+
j
≥
2
1
i
!
j
!
f
i
j
(
1
)
ϕ
1
i
(
0
)
ϕ
2
j
(
0
)
∑
i
+
j
+
l
≥
2
1
i
!
j
!
l
!
f
i
j
l
(
2
)
ϕ
1
i
(
-
τ
~
)
ϕ
2
j
(
-
τ
~
)
ϕ
2
l
(
0
)
)
,
respectively. By the Riesz representation theorem, there exists a two-dimensional matrix η
(
θ
,
μ
), θ
∈
[
-
τ
~
,
0
], whose elements are functions of bounded variation such that
(38)
L
μ
(
ϕ
)
=
∫
-
τ
~
0
d
η
(
θ
,
μ
)
ϕ
(
θ
)
,
ϕ
∈
C
(
[
-
τ
~
,
0
]
,
ℜ
2
)
.
In fact, we can choose
(39)
η
(
θ
,
μ
)
=
(
a
1
a
2
0
a
5
)
δ
(
θ
)
-
(
0
0
a
3
a
4
)
δ
(
θ
+
τ
~
)
,
where δ is the Dirac delta function.Forϕ
(
θ
)
∈
C
1
(
[
-
τ
~
,
0
]
,
ℜ
2
), define the operator A
(
μ
) as
(40)
A
(
μ
)
ϕ
(
θ
)
=
{
d
ϕ
(
θ
)
d
θ
,
θ
∈
[
-
τ
~
,
0
)
,
∫
-
τ
~
0
d
η
(
μ
,
s
)
ϕ
(
s
)
θ
=
0
,
R
(
μ
)
ϕ
(
θ
)
=
{
0
,
θ
∈
[
-
τ
~
,
0
)
,
f
(
μ
,
ϕ
)
θ
=
0
.
Then system (36) is equivalent to the operator equation of the form
(41)
u
˙
t
=
A
(
μ
)
u
t
+
R
(
μ
)
u
t
,
where u
(
t
)
=
(
x
(
t
)
,
y
(
t
)
)
T and u
t
=
u
(
t
+
θ
) for θ
∈
[
-
τ
~
,
0
].Forψ
(
θ
)
∈
C
1
(
[
0
,
τ
~
]
,
(
ℜ
2
)
*
), define
(42)
A
*
ψ
(
s
)
=
{
-
d
ψ
(
s
)
d
s
,
s
∈
(
0
,
τ
~
]
,
∫
-
τ
~
0
d
η
T
(
ξ
,
0
)
ψ
(
-
ξ
)
s
=
0
and define the bilinear form
(43)
〈
ψ
(
s
)
,
ϕ
(
θ
)
〉
=
ψ
-
(
0
)
ϕ
(
0
)
-
∫
-
τ
~
0
∫
ξ
=
0
θ
ψ
-
(
ξ
-
θ
)
d
η
(
θ
)
ϕ
(
ξ
)
d
ξ
,
where η
(
θ
)
=
η
(
θ
,
0
). Then A
(
0
) and A
* are adjoint operators.By the discussion in Section4, we know that ±
i
ω
0 are eigenvalues of A
(
0
), and therefore they are also eigenvalues of A
*. We first need to compute the eigenvector of A
(
0
) and A
* corresponding to i
ω
0 and -
i
ω
0, respectively.By direct computation, we can obtain thatq
(
θ
)
=
(
1
ρ
)
e
i
ω
0
θ and θ
∈
[
-
τ
~
,
0
] is the eigenvector of A
(
0
) corresponding to the eigenvalue i
ω
0, where
(44)
ρ
=
a
3
e
-
i
ω
0
τ
~
-
a
4
e
-
i
ω
0
τ
~
-
a
5
+
i
ω
0
.
Also we have q
*
(
s
)
=
D
(
1
σ
)
e
i
ω
0
s, and s
∈
[
0
,
τ
~
] is the eigenvector of A
* corresponding to the eigenvalue -
i
ω
0, where
(45)
σ
=
-
a
2
a
4
e
i
ω
0
τ
~
+
a
5
+
i
ω
0
.
Then we choose
(46)
D
=
1
1
+
ρ
-
σ
+
(
a
3
+
a
4
ρ
-
)
σ
τ
~
e
i
ω
0
τ
~
,
such that 〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1 and 〈
q
*
(
s
)
,
q
-
(
θ
)
〉
=
0.Using the same algorithms and similar computation process as in Hassard et al. [25], we obtain the following quantities:
(47)
g
20
=
D
-
{
f
20
(
1
)
+
2
f
11
1
ρ
+
σ
-
(
f
200
(
2
)
e
-
2
i
ω
0
τ
~
+
2
f
110
(
2
)
ρ
e
-
2
i
ω
0
τ
~
)
}
,
g
11
=
D
-
{
f
20
(
1
)
+
2
f
11
1
Re
{
ρ
}
+
σ
-
(
f
200
(
2
)
+
2
f
110
(
2
)
Re
{
ρ
}
)
}
,
g
02
=
D
-
{
f
20
(
1
)
+
2
f
11
1
ρ
-
+
σ
-
(
f
200
(
2
)
e
2
i
ω
0
τ
~
+
2
f
110
(
2
)
ρ
-
e
2
i
ω
0
τ
~
)
}
,
g
21
=
D
-
{
f
30
(
1
)
+
f
21
(
1
)
ρ
-
+
2
f
21
(
1
)
ρ
+
f
11
(
1
)
W
20
(
2
)
(
0
)
+
f
11
(
1
)
W
20
(
1
)
(
0
)
ρ
-
+
σ
-
(
f
300
(
2
)
e
-
i
ω
0
τ
~
+
f
210
(
2
)
(
2
ρ
e
-
i
ω
0
τ
~
+
ρ
-
e
-
i
ω
0
τ
~
)
)
+
f
110
(
2
)
(
W
20
(
1
)
(
-
τ
~
)
ρ
-
e
i
ω
0
τ
~
+
W
20
(
2
)
(
-
τ
~
)
e
i
ω
0
τ
~
+
2
W
11
(
1
)
(
-
τ
~
)
ρ
e
-
i
ω
0
τ
~
+
2
W
11
(
2
)
(
-
τ
~
)
e
-
i
w
0
τ
~
)
}
,
where
(48)
W
20
(
θ
)
=
i
g
20
ω
0
q
(
0
)
e
i
ω
0
θ
+
i
g
-
02
3
ω
0
q
-
(
0
)
e
-
i
ω
0
θ
+
E
1
e
2
i
ω
0
θ
,
W
11
(
θ
)
=
-
i
g
11
ω
0
q
(
0
)
e
i
ω
0
θ
+
i
g
-
11
ω
0
q
-
(
0
)
e
-
i
ω
0
θ
+
E
2
,
where E
1
=
(
E
1
(
1
)
,
E
1
(
2
)
)
∈
R
2 and E
2
=
(
E
2
(
1
)
,
E
2
(
2
)
)
∈
R
2 are constant vectors, where(49)
E
1
(
1
)
=
1
A
|
f
20
(
1
)
+
2
f
11
(
1
)
ρ
-
a
2
f
200
(
2
)
e
-
2
i
ω
0
τ
~
+
2
f
110
(
2
)
ρ
e
-
2
i
ω
0
τ
~
2
i
ω
0
-
a
5
-
a
4
e
-
2
i
ω
0
τ
~
|
,
E
1
(
2
)
=
1
A
|
2
i
ω
0
-
a
1
f
20
(
1
)
+
2
f
11
(
1
)
ρ
-
a
3
e
-
2
i
ω
0
τ
~
f
200
(
2
)
e
-
2
i
ω
0
τ
~
+
2
f
110
(
2
)
ρ
e
-
2
i
ω
0
τ
~
|
,
A
=
|
2
i
ω
0
-
a
1
-
a
2
-
a
3
e
-
2
i
ω
0
τ
~
2
i
ω
0
-
a
5
-
a
4
e
-
2
i
ω
0
τ
~
|
,
E
2
(
1
)
=
-
1
a
2
a
3
|
f
20
(
1
)
+
2
f
11
(
1
)
Re
{
ρ
}
-
a
2
f
200
(
2
)
+
2
f
110
(
2
)
Re
{
ρ
}
0
|
,
E
2
(
2
)
=
-
1
a
2
a
3
|
-
a
1
f
20
(
1
)
+
2
f
11
(
1
)
Re
{
ρ
}
-
a
3
f
200
(
2
)
+
2
f
110
(
2
)
Re
{
ρ
}
|
.
Therefore, we can compute the following values:
(50)
C
1
(
0
)
=
i
2
ω
0
(
g
20
g
11
-
2
|
g
11
|
2
-
1
3
|
g
02
|
2
)
+
1
2
g
21
,
μ
2
=
-
Re
{
C
1
(
0
)
}
Re
{
λ
′
(
τ
~
)
}
,
β
2
=
2
Re
{
C
1
(
0
)
}
,
T
2
=
-
Im
⁡
{
C
1
(
0
)
}
+
μ
2
Im
⁡
{
λ
′
(
τ
~
)
}
ω
0
.
Here μ
2 determines the direction of the Hopf bifurcation, β
2 determines the stability of the bifurcating periodic solutions, and T
2 determines the periods of the bifurcating periodic solutions.Based on the discussion above, we can obtain the following results.Theorem 6.
The sign ofμ
2 determines the direction of Hopf bifurcation: if μ
2
>
0  (
μ
2
<
0
), then the Hopf bifurcation is supercritical (subcritical); β
2 determines the stability of the bifurcating periodic solution; the bifurcating periodic solution is stable (unstable) if β
2
<
0
(
β
2
>
0
); T
2 determines the period of the bifurcating periodic solution; the period increases (decreases) if T
2
>
0
(
T
2
<
0
).


## 6. Numerical Simulation

As the problem is not a case study, the real-world data are not available for this model. We, therefore, take here some hypothetical data with the sole purpose of illustrating the results that we have established in the previous sections.(i) Let us consider the parameters of the system (2) as r
=
0.2, a
=
1, β
=
0.8, d
=
0.5, k
=
10, q
=
2, and m
=
0.9 in appropriate units. For this value of parameters we get the critical value of E as E
*
=
0.0949124. Thus it is easy to verify that for this set of parameters the system (2) is locally asymptotically stable around its interior equilibrium (
x
*
,
y
*
) for E
>
E
* and is unstable for E
<
E
*. Thus for E
=
E
*
=
0.0949124, the system (2) undergoes a Hopf bifurcation. Now for E
=
0.1, we have interior equilibrium (
5.916
,
0.552
) which is asymptotically stable (see Figures 2 and 3), but for E
=
0.05, the interior equilibrium (
2.236
,
0.463
) is unstable (see Figures 4 and 5). Thus taking E as control parameter, it is possible to drive the prey-predator system to require equilibrium and to prevent the cycle behaviour of the system.Figure 2
The system is asymptotically stable forE
=
0.1
>
0.0949124
=
E
*.Figure 3
Phase plane trajectories corresponding to different initial levels. The figure clearly indicates that the equilibrium(
5.916
,
0.552
) is asymptotically stable when E
=
0.1
>
0.0949124
=
E
*.Figure 4
The system is unstable forE
=
0.05
<
0.0949124
=
E
*.Figure 5
There is a stable limit cycle surrounding(
2.236
,
0.463
) for E
=
0.05
<
0.0949124
=
E
*.(ii) For the valuesr
=
1.5, a
=
1, β
=
0.2, m
=
0.9, d
=
0.02, k
=
5, q
=
2, and E
=
0.06, then the system (1) has a positive equilibrium (
1.87083
,
2.50892
). The conditions (
H
1
) and (
H
2
) are satisfied, so the positive equilibrium (
1.87083
,
2.50892
) is asymptotically stable for τ
=
0 (see Figure 6). The conditions (
H
1
), (
H
2
), and (
H
4
) are satisfied, so we compute that ω
0
=
0.227539 and τ
0
=
0.799696. Therefore, by Theorem 4, we see that the positive equilibrium (
x
*
,
y
*
) is stable when τ
<
τ
0 (see Figures 7 and 8) and unstable when τ
>
τ
0 (see Figures 9 and 10). Thus the system undergoes a Hopf bifurcation at τ
0. Thus when τ
0
=
0.799696, ω
0
=
0.227539, we compute λ
′
(
τ
0
)
=
-
0.030746
-
0.015498
i, C
1
(
0
)
=
-
0.441228
-
1.442407
i, μ
2
=
-
14.350865, β
2
=
-
0.882456, and T
2
=
5.361681. Since μ
2
<
0, β
2
<
0, and T
2
>
0, therefore, the Hopf bifurcation is subcritical, the bifurcating periodic solutions at τ
0 are stable, and period of bifurcating periodic solutions increases.Figure 6
Both the prey and predator populations converge to their equilibrium values.Figure 7
The equilibrium(
1.87083
,
2.50892
) is asymptotically stable for τ
=
0.1
<
0.799696
=
τ
0.Figure 8
Phase plane trajectories corresponding to different initial levels. The figure clearly indicates that the equilibrium(
1.87083
,
2.50892
) is asymptotically stable when τ
=
0.1.Figure 9
Bifurcating periodic solutions from(
1.87083
,
2.50892
) occur when τ
=
1
>
0.799696
=
τ
0.Figure 10
Bifurcating periodic solutions from(
1.87083
,
2.50892
) occur when τ
=
1
>
0.799696
=
τ
0.


## 7. Concluding Remarks

This paper deals with a delayed prey-predator model with Holling type III functional response and harvesting of predator species. Oscillatory behavior and existence of limit cycles in harvested predator-prey system are common in nature. In the absence of delay, we have proved that exactly one stable limit cycle occurs when positive equilibrium is unstable. Since harvesting is associated with economic interests, our analysis shows that the harvesting of predator species plays an important role in the shaping of the dynamical behavior of the system. From Figure1, we see that harvesting effort affects the interior equilibrium of system (2), and, as the harvesting effort increases, the interior equilibrium may be vanished, which means that the population of predator may become extinct. It is also observed from the obtained results that harvesting effort can cause an unstable equilibrium to become stable and even a simple Hopf bifurcation occurred when the parameter E passes through its critical value.Taking delay as bifurcation parameters, we see that a Hopf bifurcation occurs, whenever delay increases a critical value. Furthermore, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are determined by the normal form theory and center manifold theorem.



---

*Source: 543041-2014-03-03.xml*