# A Stochastic Predator-Prey System with Stage Structure for Predator

Authors: Shufen Zhao; Minghui Song
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The authors introduce stochasticity into a predator-prey system with Beddington-DeAngelis functional response and stage structure for predator. We present the global existence and positivity of the solution and give sufficient conditions for the global stability in probability of the system. Numerical simulations are introduced to support the main theoretical results.


---

## Body

## 1. Introduction

The classical predator-prey model with Beddington-DeAngelis type functional response can be denoted as(1)
x
˙
(
t
)
=
x
[
b
1
-
a
11
x
-
a
12
y
1
+
m
x
+
n
y
]
,
y
˙
=
y
[
-
b
2
+
a
21
x
1
+
m
x
+
n
y
-
a
22
y
]
,
where x
(
t
) and y
(
t
) represent predator and prey densities at time t, respectively. b
i, a
i
j, m, and n are positive constants, i, j
=
1
,
2. For biological representation of each coefficient in (1) we refer the reader to [1, 2]. In model (1), it is assumed that all individuals of a single species have largely similar capabilities to hunt or to reproduce. But the life cycle of most animals consists of at least two stages, immature and mature, and the individuals in the first stage often can neither hunt or reproduce, being raised by their mature parents and there are recognizable morphological and behavioral differences that may exist between these stages. In [3], the authors studied the global properties of a predator-prey model with nonlinear functional response and stage structure for the predator, and the condition of the existence and the global stability of the positive steady states were established. However, May [4] pointed out that due to environment noise, the birth rate, carrying capacity, competition coefficients, and other parameters involved in the system exhibit random fluctuation to a greater or lesser extent. So, a lot of authors introduced stochastic noise into deterministic models to reveal the effect of environmental variability on the population dynamics in mathematical ecology [5–7]. In Liu and Wang [5], the authors investigated the global stability of stage-structured predator-prey models with Beddington-DeAngelis type functional response and with stage structure for the prey. The authors [5] also pointed out that there are some technical obstacles that cannot be overcome at present to investigate the stage structure on predator model. So, in this paper we are going to do some work on this problem. The following model,
(2)
d
x
d
t
=
x
(
γ
-
q
y
2
1
+
m
x
+
n
y
2
-
β
x
)
,
d
y
1
d
t
=
α
y
2
-
d
1
y
1
-
a
1
y
1
2
-
b
y
1
,
d
y
2
d
t
=
b
y
1
-
d
2
y
2
-
a
2
y
2
2
+
p
y
2
x
1
+
m
x
+
n
y
2
,is derived under the following assumptions.(H1)
The immature predator populationy
1: the birth rate into the immature population is proportional to the existing mature predator population with probability α
>
0; the death rate is proportional to the existing immature predator population with proportionality d
1
>
0; overcrowding rate of the immature predator population is a
1
>
0; the transformation rate from the immature predator to mature predator is proportional to the existing immature predator population with proportionality b
>
0.
(H2)
The mature predator populationy
2: d
2
>
0 and a
2
>
0 are the death rate and the overcrowding rate of the mature predator population, respectively, and only the mature predator population feeds the prey. It seems reasonable that a number of mammals, who are immature predators, are raised by their parents. p
/
q is the rate of conversion of nutrients into the reproduction of the predator.
(H3)
The prey populationx: the growth of the species is of the Lotka-Volterra nature and γ
>
0 is the birth rate; β
>
0 is the overcrowding rate. q
>
0 is the capturing rate of the predator.System (2) is greatly different from the model investigated in [3] for we comprehend that the effect of the response function will diminish the death rate of the predator and the predator does not only feed on the prey.Suppose that(
x
*
,
y
1
*
,
y
2
*
) is a positive equilibrium of (2). If we take the environmental noise into account, we can replace the birth rate of prey population and death rate of predator population by an average value plus a random fluctuation, respectively,
(3)
γ
+
σ
(
x
-
x
*
)
B
˙
(
t
)
,
d
1
+
σ
1
(
y
1
-
y
1
*
)
B
˙
1
(
t
)
,
d
2
+
σ
2
(
y
2
-
y
2
*
)
B
˙
2
(
t
)
,
where σ
2, σ
1
2, σ
2
2 represent the intensities of the noise and B
˙
(
t
), B
˙
1
(
t
), B
˙
2
(
t
) are standard white noise; namely, B
(
t
), B
1
(
t
), B
2
(
t
) are standard Brownian motion defined on a complete probability space (
Ω
,
ℱ
,
P
) with a filtration {
ℱ
t
}
t
∈
R
+ satisfying the usual condition (i.e., it is right continuous and increasing while ℱ
0 contains all P-null sets). So the corresponding stochastic system of (2) is
(4)
d
x
=
x
(
γ
-
q
y
2
1
+
m
x
+
n
y
2
-
β
x
)
d
t
+
σ
x
(
x
-
x
*
)
d
B
(
t
)
,
d
y
1
=
(
α
y
2
-
d
1
y
1
-
a
1
y
1
2
-
b
y
1
)
d
t
+
σ
1
y
1
(
y
1
-
y
1
*
)
d
B
1
(
t
)
,
d
y
2
=
(
b
y
1
-
d
2
y
2
-
a
2
y
2
2
+
p
x
y
2
1
+
m
x
+
n
y
2
)
d
t
-
σ
2
y
2
(
y
2
-
y
2
*
)
d
B
2
(
t
)
,
with the initial condition (
x
(
0
)
,
y
1
(
0
)
,
y
2
(
0
)
)
∈
R
+
3 where R
+
3
=
{
x
∈
R
3
∣
x
,
y
1
,
y
2
>
0
,
i
=
1,2
,
3
}.The paper is organized as follows. In Section2, we prove the existence, uniqueness, and the positivity of the solution to (4). In Section 3, we established the condition for the global stability of the positive equilibrium. We work out two simulation figures to illustrate our main results in Section 4. Section 5 gives the conclusions and future directions.


## 2. Existence of the Global Positive Solution

Theorem 1.
For any initial value(
x
(
0
)
,
y
1
(
0
)
,
y
2
(
0
)
)
∈
R
+
3, system (4) has a unique global positive solution (
x
(
t
)
,
y
1
(
t
)
,
y
2
(
t
)
) on t
>
0 with probability one.Proof.
We see that the coefficients of the system are locally Lipschitz continuous, so, for any given initial valuesx
(
0
)
>
0, y
1
(
0
)
>
0, y
2
(
0
)
>
0, there is a unique maximal local solution (
x
(
t
)
,
y
1
(
t
)
,
y
2
(
t
)
) on t
∈
[
0
,
τ
e
), where τ
e is the explosion time [8]. To show this solution is global, we need to show that τ
e
=
∞. Define the stopping time by
(5)
τ
k
=
inf
⁡
⁡
{
t
∈
(
0
,
∞
)
,
y
i
(
t
)
=
1
k
or
x
(
t
)
=
1
k
or
x
(
t
)
=
k
}
or
y
i
(
t
)
=
k
,
i
=
1,2
,
or
x
(
t
)
=
1
k
or
x
(
t
)
=
k
}
,
where, throughout this paper, we set inf
⁡
⁡
∅
=
∞. Clearly, τ
k is increasing as k
→
∞. Set τ
∞
=
lim
⁡
k
→
∞
τ
k, whence τ
∞
≤
τ
e a.s. If we can show that τ
∞
=
∞ a.s., then τ
e
=
∞ a.s. Namely, to complete the proof, it is sufficient to show that τ
∞
=
∞ a.s. If this statement is false, then there is a pair of constants T
>
0 and ϵ
∈
(
0,1
), such that P
{
τ
∞
≤
T
}
>
ϵ. Hence, there is an integer k
1
≥
k
0 such that
(6)
P
{
τ
k
≤
T
}
≥
ϵ
∀
k
≥
k
1
,
where k
0 is satisfying 1
/
k
0
<
y
i
(
0
)
<
k
0, i
=
1,2, and 1
/
k
0
<
x
(
0
)
<
k
0.
Define aC
2-function V: R
+
3
→
R
+ by
(7)
V
(
x
,
y
1
,
y
2
)
=
(
x
-
1
-
0.5
ln
⁡
⁡
x
)
+
(
y
1
-
1
-
0.5
ln
⁡
⁡
y
1
)
+
(
y
2
-
1
-
0.5
ln
⁡
⁡
y
2
)
.
Using Itô’s formula, we get
(8)
d
V
=
1
2
(
x
-
1
/
2
-
x
-
1
)
·
d
x
+
1
2
(
y
1
-
1
/
2
-
y
1
-
1
)
·
d
y
1
+
1
2
(
y
2
-
1
/
2
-
y
2
-
1
)
·
d
y
2
+
1
2
(
-
1
4
x
-
3
/
2
+
1
2
x
-
2
)
·
(
d
x
)
2
+
1
2
(
-
1
4
y
1
-
3
/
2
+
1
2
y
1
-
2
)
·
(
d
y
1
)
2
+
1
2
(
-
1
4
y
2
-
3
/
2
+
1
2
y
2
-
2
)
·
(
d
y
2
)
2
,
(9)
L
V
=
1
2
(
x
1
/
2
-
1
)
(
γ
-
q
y
2
1
+
m
x
+
n
y
2
-
β
x
)
+
1
2
(
-
1
4
x
1
/
2
+
1
2
)
σ
2
(
x
-
x
*
)
2
+
1
2
(
y
1
-
1
/
2
-
y
1
-
1
)
(
α
y
2
-
d
1
y
1
-
a
1
y
1
2
-
b
y
1
)
+
1
2
(
-
1
4
y
1
1
/
2
+
1
2
)
σ
1
2
(
y
1
-
y
1
*
)
2
+
1
2
(
y
2
-
1
/
2
-
y
2
-
1
)
×
(
b
y
1
-
d
2
y
2
-
a
2
y
2
2
+
p
x
y
2
1
+
m
x
+
n
y
2
)
+
1
2
(
-
1
4
y
2
1
/
2
+
1
2
)
σ
2
2
(
y
2
-
y
2
*
)
2
=
1
2
(
x
1
/
2
-
1
)
(
γ
-
q
y
2
1
+
m
x
+
n
y
2
-
β
x
)
+
1
2
(
y
2
-
1
/
2
-
y
2
-
1
)
b
y
1
+
1
2
(
y
1
-
1
/
2
-
y
1
-
1
)
α
y
2
+
1
2
(
y
1
1
/
2
-
1
)
(
-
d
1
-
a
1
y
1
-
b
)
+
1
2
(
y
2
1
/
2
-
1
)
(
-
d
2
-
a
2
y
2
+
p
y
2
1
+
m
x
+
n
y
2
)
+
1
2
(
-
1
4
x
1
/
2
+
1
2
)
σ
2
(
x
-
x
*
)
2
+
1
2
(
-
1
4
y
1
1
/
2
+
1
2
)
σ
1
2
(
y
1
-
y
1
*
)
2
+
1
2
(
-
1
4
y
2
1
/
2
+
1
2
)
σ
2
2
(
y
2
-
y
2
*
)
2
.
Now, we pay attention to the term (
1
/
2
)
(
y
1
-
1
/
2
-
y
1
-
1
)
α
y
2. If y
1
<
1, then (
1
/
2
)
(
y
1
-
1
/
2
-
y
1
-
1
)
α
y
2
<
0, so, this term can be omitted from the right side of the inequality. If y
1
>
1, then (
1
/
2
)
(
y
1
-
1
/
2
-
y
1
-
1
)
α
y
2
<
(
1
/
2
)
(
y
1
1
/
2
-
1
)
α
y
2. The similar argument can be taken on (
1
/
2
)
(
y
2
-
1
/
2
-
y
2
-
1
)
b
y
1. So we get the following inequality:
(10)
L
V
≤
1
2
[
-
β
x
3
/
2
+
β
x
+
γ
x
1
/
2
+
q
n
-
γ
]
+
b
2
y
1
(
y
2
-
1
)
+
α
2
y
2
(
y
1
-
1
)
+
1
2
[
-
a
1
y
1
3
/
2
-
(
d
1
+
b
)
y
1
1
/
2
+
a
1
y
1
+
d
1
+
b
]
+
1
2
[
-
a
2
y
2
3
/
2
+
(
p
n
-
d
)
y
2
1
/
2
+
a
2
y
2
+
d
2
]
+
1
2
σ
2
(
-
1
4
x
5
/
2
+
1
2
x
2
+
1
2
x
*
x
3
/
2
-
x
*
x
-
1
4
x
2
*
2
x
1
/
2
+
1
2
x
2
*
2
)
+
1
2
σ
1
2
(
-
1
4
y
1
5
/
2
+
1
2
y
1
2
+
1
2
y
1
*
y
1
3
/
2
-
y
1
*
y
1
-
1
4
y
1
*
2
y
1
1
/
2
+
1
2
y
1
*
2
)
+
1
2
σ
2
2
(
-
1
4
y
2
5
/
2
+
1
2
y
2
2
+
1
2
y
2
*
y
2
3
/
2
-
y
2
*
y
2
-
1
4
y
2
*
2
y
2
1
/
2
+
1
2
y
2
*
2
)
≤
1
2
[
-
β
x
3
/
2
+
β
x
+
γ
x
1
/
2
+
q
n
]
+
1
2
[
+
(
a
1
-
b
+
α
2
)
y
1
+
d
1
+
b
[
b
2
y
1
2
-
a
1
y
1
3
/
2
-
(
d
1
+
b
)
y
1
1
/
2
]
b
2
y
1
2
-
a
1
y
1
3
/
2
-
(
d
1
+
b
)
y
1
1
/
2
+
(
a
1
-
b
+
α
2
)
y
1
+
d
1
+
b
[
b
2
y
1
2
-
a
1
y
1
3
/
2
-
(
d
1
+
b
)
y
1
1
/
2
]
+
1
2
[
+
(
a
2
-
α
+
b
2
)
y
2
+
d
2
]
+
(
a
2
-
α
+
b
2
)
y
2
+
d
2
]
α
2
y
2
2
-
a
2
y
2
3
/
2
+
(
p
n
-
d
)
y
2
1
/
2
+
(
a
2
-
α
+
b
2
)
y
2
+
d
2
]
+
1
2
σ
2
(
-
1
4
x
5
/
2
+
1
2
x
2
+
1
2
x
*
x
3
/
2
-
x
*
x
-
1
4
x
2
*
2
x
1
/
2
+
1
2
x
2
*
2
)
+
1
2
σ
1
2
(
-
1
4
y
1
5
/
2
+
1
2
y
1
2
+
1
2
y
1
*
y
1
3
/
2
-
y
1
*
y
1
-
1
4
y
1
*
2
y
1
1
/
2
+
1
2
y
1
*
2
)
+
1
2
σ
2
2
(
-
1
4
y
2
5
/
2
+
1
2
y
2
2
+
1
2
y
2
*
y
2
3
/
2
-
y
2
*
y
2
-
1
4
y
2
*
2
y
2
1
/
2
+
1
2
y
2
*
2
)
≤
K
,
where K is positive constant. Integrating both sides of (8) from 0 to τ
k
∧
T and then taking the expectations lead to
(11)
V
(
x
(
τ
k
∧
T
)
,
y
1
(
τ
k
∧
T
)
,
y
2
(
τ
k
∧
T
)
)
≤
V
(
x
(
0
)
,
x
1
(
0
)
,
y
2
(
0
)
)
+
K
T
.
Setting Ω
k
=
P
{
τ
k
≤
T
} for k
≥
k
1, then by the inequality (6) we have P
(
Ω
k
)
≥
ϵ; note that for every ω
∈
Ω
k, V
(
x
(
τ
k
,
ω
)
,
y
1
(
τ
k
,
ω
)
,
y
2
(
τ
k
,
ω
)
) is no less than min
⁡
⁡
{
k
-
1
-
0.5
ln
⁡
⁡
k, 1
/
k
-
1
-
0.5
ln
⁡
⁡
(
1
/
k
)
}.
It then follows from (11) that
(12)
E
[
I
Ω
k
V
(
x
(
τ
k
)
,
x
1
(
τ
k
)
,
y
2
(
τ
k
)
)
]
≤
V
(
x
(
0
)
,
y
1
(
0
)
,
y
2
(
0
)
)
+
K
T
,
where I
Ω
k is the indicator function of Ω
k and E
[
I
Ω
k
V
(
x
(
τ
k
)
,
y
1
(
τ
k
)
,
y
2
(
τ
k
)
)
] ≥ ϵ
min
⁡
⁡
{
k
-
1
-
0.5
ln
⁡
⁡
k, 1
/
k
-
1
-
0.5
ln
⁡
⁡
(
1
/
k
)
}. So,
(13)
ϵ
min
⁡
⁡
{
k
-
1
-
0.5
ln
⁡
⁡
k
,
1
k
-
1
-
0.5
ln
⁡
⁡
1
k
}
≤
V
(
x
,
y
1
(
0
)
(
0
)
,
y
2
(
0
)
)
+
K
T
.
Letting k
→
∞ leads to the contradiction
(14)
∞
>
V
(
x
(
0
)
,
y
1
(
0
)
,
y
2
(
0
)
)
+
K
T
≥
∞
.
This contradiction shows that τ
∞
=
∞, which completes the proof.


## 3. Global Behavior

Supposez
=
z
(
t
) is the solution of the following n-dimensional stochastic differential equation:
(15)
d
z
(
t
)
=
f
(
z
(
t
)
,
t
)
d
t
+
g
(
z
(
t
)
,
t
)
d
B
(
t
)
,
and z
* is the equilibrium position of (15).From the stability theory of stochastic differential equations, we only need to find a Lyapunov functionV
(
z
) satisfying L
V
(
z
)
≤
0 and the identity holds if and only if z
=
z
* [9], where z
=
z
(
t
) is the solution of the n-dimensional stochastic differential equation (15) and d
V
(
x
(
t
)
,
t
)
=
L
V
d
t
+
V
x
(
x
(
t
)
,
t
)
g
(
t
)
d
B
(
t
).Theorem 2.
Ifa
1
-
(
σ
1
2
/
2
)
>
0, a
2
-
(
σ
2
2
/
2
)
>
0, and β
-
(
σ
2
/
2
)
-
(
q
m
y
2
*
/
(
1
+
m
x
*
+
n
y
2
*
)
)
>
0, then the positive equilibrium (
x
*
,
y
1
*
,
y
2
*
) of model (4) is globally asymptotically stable with probability one.Proof.
System (4) can be rewritten as
(16)
d
x
=
q
x
[
m
y
2
*
(
x
-
x
*
)
-
(
1
+
m
x
*
)
(
y
2
-
y
2
*
)
]
(
1
+
m
x
+
n
y
2
)
(
1
+
m
x
*
+
n
y
2
*
)
d
t
-
β
x
(
x
-
x
*
)
d
t
+
σ
x
(
x
-
x
*
)
d
B
(
t
)
,
d
y
1
=
α
y
1
*
[
y
1
(
y
2
-
y
2
*
)
-
y
2
(
y
1
-
y
1
*
)
]
d
t
-
a
1
y
1
(
y
1
-
y
1
*
)
d
t
+
σ
1
y
1
(
y
1
-
y
1
*
)
d
B
1
(
t
)
,
d
y
2
=
b
y
2
*
[
y
2
(
y
1
-
y
1
*
)
-
y
1
(
y
2
-
y
2
*
)
]
d
t
+
p
y
2
[
(
1
+
n
y
2
*
)
(
x
-
x
*
)
-
n
x
*
(
y
2
-
y
2
*
)
]
(
1
+
m
x
+
n
y
2
)
(
1
+
m
x
*
+
n
y
2
*
)
d
t
-
a
2
y
2
(
y
2
-
y
2
*
)
d
t
+
σ
2
y
2
(
y
2
-
y
2
*
)
d
B
2
(
t
)
.
Define
(17)
V
(
x
1
,
x
2
,
y
)
=
c
1
(
x
-
x
*
-
x
*
ln
⁡
⁡
(
x
x
*
)
)
+
c
2
(
y
1
-
y
1
*
-
y
1
*
ln
⁡
⁡
(
y
1
y
1
*
)
)
+
c
3
(
y
2
-
1
-
y
2
*
ln
⁡
⁡
y
2
y
2
*
)
,
where c
i
(
i
=
1,2
,
3
) are positive numbers to be determined. Applying Itô’s formula to system (16) gives
(18)
L
V
=
c
2
{
-
a
1
(
y
1
-
y
1
*
)
2
+
σ
1
2
2
(
y
1
-
y
1
*
)
2
}
(
y
1
-
y
1
*
)
α
y
1
*
[
(
y
2
-
y
2
*
)
-
y
2
y
1
(
y
1
-
y
1
*
)
]
-
a
1
(
y
1
-
y
1
*
)
2
+
σ
1
2
2
(
y
1
-
y
1
*
)
2
}
+
c
3
{
-
a
2
(
y
2
-
y
2
*
)
2
+
σ
2
2
2
(
y
2
-
y
2
*
)
2
}
(
y
2
-
y
2
*
)
b
y
2
*
[
(
y
1
-
y
1
*
)
-
y
1
y
2
(
y
2
-
y
2
*
)
]
-
a
2
(
y
2
-
y
2
*
)
2
+
σ
2
2
2
(
y
2
-
y
2
*
)
2
}
+
c
3
p
(
y
2
-
y
2
*
)
(
1
+
n
x
*
)
(
x
-
x
*
)
-
n
x
*
(
y
2
-
y
2
*
)
(
1
+
m
x
+
n
y
2
)
(
1
+
m
x
*
+
n
y
2
*
)
+
c
1
q
(
x
-
x
*
)
m
y
2
*
(
x
-
x
*
)
-
(
1
+
m
x
*
)
(
y
2
-
y
2
*
)
(
1
+
m
x
+
n
y
2
)
(
1
+
m
x
*
+
n
y
2
*
)
-
c
1
β
(
x
-
x
*
)
2
+
c
1
σ
2
2
(
x
-
x
*
)
2
.
Setc
2
=
y
1
*
/
α, c
3
=
y
2
*
/
b, c
1
=
c
3
p
(
1
+
n
x
2
*
)
/
q
(
1
+
m
y
*
)
=
(
y
2
*
/
b
)
(
p
/
q
)
(
(
1
+
m
x
*
)
/
(
1
+
n
y
*
)
). Then we have
(19)
L
V
=
{
-
y
2
y
1
(
y
1
-
y
1
*
)
2
+
2
(
y
1
-
y
1
*
)
(
y
2
-
y
2
*
)
-
y
1
y
2
(
y
2
-
y
2
*
)
2
}
-
c
2
(
a
1
-
σ
1
2
2
)
(
y
1
-
y
1
*
)
2
-
c
3
(
a
2
-
σ
2
2
2
)
(
y
2
-
y
2
*
)
2
-
c
3
p
n
x
*
(
y
2
-
y
2
)
2
(
1
+
m
x
+
n
y
2
)
(
1
+
m
x
*
+
n
y
2
*
)
-
c
1
(
β
-
σ
2
2
)
(
x
-
x
*
)
2
+
c
1
q
y
2
*
m
(
x
-
x
*
)
2
(
1
+
m
x
+
n
y
2
)
(
1
+
m
x
*
+
n
y
2
*
)
≤
-
[
y
2
y
1
(
y
1
-
y
1
*
)
-
y
1
y
2
(
y
2
-
y
2
*
)
]
2
-
y
1
*
α
(
a
1
-
σ
1
2
2
)
(
y
1
-
y
1
*
)
2
-
(
y
2
*
b
a
2
-
y
2
*
b
σ
2
2
2
)
(
y
2
-
y
2
*
)
2
-
y
2
*
b
p
q
1
+
n
x
*
1
+
m
y
*
×
(
β
-
σ
2
2
-
q
m
y
2
*
1
+
m
x
*
+
n
y
2
*
)
(
x
-
x
*
)
2
≤
-
y
1
*
α
(
a
1
-
σ
1
2
2
)
(
y
1
-
y
1
*
)
2
-
y
2
*
b
(
a
2
-
σ
2
2
2
)
×
(
y
2
-
y
2
*
)
2
-
y
2
*
b
p
q
1
+
m
x
*
1
+
n
y
*
×
(
β
-
σ
2
2
-
q
m
y
2
*
1
+
m
x
*
+
n
y
2
*
)
(
x
-
x
*
)
2
.
The condition in Theorem 2 implies L
V
≤
0, and the identity holds if and only if (
x
,
y
1
,
y
2
)
=
(
x
*
,
y
1
*
,
y
2
*
). By Theorem  2.1 in [9] and the description of that theorem, we get the conclusion.


## 4. Numerical Simulations

In this section, we will use the Euler method and the Milstein method mentioned in [10] to substantiate the analytical findings. For system (4), consider the discretization equations
(20)
x
(
k
+
1
)
-
x
(
k
)
=
x
(
k
)
(
γ
-
q
y
2
(
k
)
1
+
m
x
(
k
)
+
n
y
2
(
k
)
-
β
x
(
k
)
)
+
σ
x
(
k
)
(
x
(
k
)
-
x
*
)
Δ
t
ξ
(
k
)
+
σ
2
2
(
x
(
k
)
-
x
*
)
[
(
ξ
(
k
)
)
2
-
1
]
,
y
1
(
k
+
1
)
-
y
1
(
k
)
=
[
α
y
2
(
k
)
-
d
1
y
1
(
k
)
-
a
1
(
y
1
(
k
)
)
2
-
b
y
1
(
k
)
]
Δ
t
-
σ
1
y
1
(
k
)
(
y
1
(
k
)
-
y
1
*
)
Δ
t
ζ
(
k
)
+
σ
1
2
2
(
y
1
(
k
)
-
y
1
*
)
[
(
ζ
(
k
)
)
2
-
1
]
,
y
2
(
k
+
1
)
-
y
2
(
k
)
=
[
+
p
y
2
(
k
)
x
(
k
)
1
+
m
x
(
k
)
+
n
y
2
(
k
)
]
b
y
1
(
k
)
-
d
2
y
2
(
k
)
-
a
2
(
y
2
(
k
)
)
2
+
p
y
2
(
k
)
x
(
k
)
1
+
m
x
(
k
)
+
n
y
2
(
k
)
]
Δ
t
-
σ
2
y
2
(
k
)
(
y
2
(
k
)
-
y
2
*
)
Δ
t
η
(
k
)
+
σ
2
2
2
(
y
2
(
k
)
-
y
2
*
)
[
(
η
(
k
)
)
2
-
1
]
,
where ξ
(
k
), ζ
(
k
), η
(
k
), k
=
1,2
,
…
,
n are the Gaussian random variables which follow N
(
0,1
).In Figure1, we show the dynamics of the deterministic model with parameters α
=
1.2, d
1
=
0.26, b
=
0.58, a
1
=
0.01, q
=
1.5, p
=
1, β
=
0.375, d
2
=
0.32, a
2
=
0.13, m
=
1, n
=
0.081, γ
=
5.9, σ
=
σ
1
=
σ
2
=
0; then x
*
=
12.37, y
1
*
=
13.1, y
2
*
=
10.6.Figure 1
The solution of (2) with initial values x
(
0
)
=
10, y
1
(
0
)
=
10, y
2
(
0
)
=
12.In Figure2, we choose the same parameter values as Figure 1 except that σ
=
0.4, σ
1
=
0.45, σ
2
=
0.375, which satisfy the condition in Theorem 2, so Figure 2 clearly supports the conclusion of Theorem 2.Figure 2
The solution of the stochastic model (4) with the same parameters as in Figure 1.In Figure3, we choose α
=
0.9, d
1
=
0.34, b
=
0.4, q
=
0.9, p
=
0.4, β
=
0.15, d
2
=
0.12, a
1
=
0.1, a
2
=
0.6, m
=
0.8, n
=
0.5, γ
=
0.75, σ
=
σ
1
=
σ
2
=
0; then x
*
=
3.066, y
1
=
1.157, y
2
=
1.1.; in Figure 4, we choose the same parameter values as Figure 1 except that σ
1
=
0.2; σ
2
=
0.42; σ
3
=
0.25. So the conditions of our theoretical results hold. Obviously, the numerical simulations are indeed confirming our analytical results.Figure 3
The solution of (2) with initial values x
(
0
)
=
3, y
1
(
0
)
=
3.6,y
2
(
0
)
=
3.2.Figure 4
The solution of the stochastic model (4) with the same parameters as in Figure 3.


## 5. Discussion

In this paper, a stochastic predator-prey model with stage structure for the predator has been proposed and investigated. We discuss the biological significance of the model and establish sufficient conditions for global asymptotic stability of the model. These results are important because from the biological point of view, a global stable positive equilibrium means that the community consisting of two species is a stable biotic community in which all species will coexist. To the best of our knowledge, the present paper is the first attempt to study system (4).Although we only consider the global stability of the positive equilibrium, some interesting questions deserved investigation, like the stage structure effect on the long term behavior of the system. In fact, in (4) we have supposed that the predator is not only feeding on prey; we can also discuss the case in which the predator feeds on prey only. We want to mention that we are unable to give the sufficient conditions under which system (4) or (2) has a positive equilibrium, for there are some technical obstacles that cannot be overcome at present stage. However, the values in Figure 1 show that the system (2) has the positive equilibrium position in some case, and we leave this for future work.



---

*Source: 518695-2014-04-30.xml*