# Influence of Gestation Delay and Predator’s Interference in Predator-Prey Interaction under Stochastic Environment

Authors: Debaldev Jana
Journal: International Journal of Stochastic Analysis
 (2014)

---

## Abstract

Previous experimental and theoretical studies suggest that predator’s interference in predator-prey relationship provides better descriptions of predator’s feeding over a range of predator-prey abundances. Also biological delays and environmental stochasticity play an important role to describe the system and its values. In this present study, I consider a Gaussian white-noise induced stochastic predator-prey model with the Beddington-DeAngelis functional response and gestation delay. Stochastic stability is measured by second order moment terms by calculating the nonequilibrium fluctuation of the nondelayed system and Fourier transform technique depicts the fluctuation of stochastic stability by introducing time lag. Different dynamical behaviors for both situations have been illustrated numerically also. The biological implications of my analytical and numerical findings are discussed critically.


---

## Body

## 1. Introduction

In ecology, predation describes a biological interaction where a predator (an organism that is hunting) feeds on its prey (the organism that is attacked). Predators may or may not kill their prey prior to feeding on them, but the act of predation often results in the death of its prey and the eventual absorption of the prey’s tissue through consumption. Predators can have profound impacts on the dynamics of their prey that depend on how predator consumption is affected by prey density (the predators functional response). Understanding the prey-predator relationship is the central goal in ecology and a very significant component of this is the predators rate of feeding upon prey. The simplification of the food web structure due to predator interference allows qualitative predictions concerning the response of a food web to an external perturbation [1–5]. Ginzburg and Akcakaya [6] demonstrated a positive relationship between the abundances of all trophic levels and concluded that the ratio-dependent models provide a better representation of predator-prey interactions than prey-dependent models. Predators functional response, defined as the amount of prey catch per predator per unit of time, is affected by the structure of prey habitat and predators hunting ability [7, 8]. Predator's functional response is defined as the amount of prey catch per predator per unit of time. After the pioneering studies of [9], that is, linear, only prey-dependent and ratio-dependent functional response, response functions of predators which depends on both prey and predator abundances, can provide better descriptions of predator feeding over a range of predator-prey abundances because of predator interference. More recent theoretical work has demonstrated that the mathematical form of the feeding rate can influence the distribution of predators through space [10], the stability of enriched predator-prey systems [11, 12] correlations between nutrient enrichment and the biomass of higher trophic levels [11], and the length of food chains [13]. So, a prey-predator model with Beddington-DeAngelis functional response is represented by
(1)
d
x
d
t
=
r
x
(
1
-
x
k
)
-
a
x
y
1
+
β
y
+
α
x
,
d
y
d
t
=
e
a
x
y
1
+
β
y
+
α
x
-
μ
y
,
where r, k, a, α, β, e, and μ are assumed to be positive here.Delay is frequently used in a predator-prey model to make the model biologically more realistic. Recently, researchers are using more than one delay to study the effect of past history of the system populations [14–18], as in reality time delays occur in almost every biological situation [19] and are assumed to be one of the reasons of regular fluctuations in population density [20]. Delay is frequently introduced in a biologically realistic predator-prey model. In a review paper of predator-prey models with discrete delay, Ruan [21] discussed different types of delays and the dynamics of the corresponding models. Reproduction of predator after consuming prey is not instantaneous but mediated by some time lag required for gestation. So, time lag between capturing prey and its corresponding contribution to predator's growth is an important factor to model formulation and its biological clarification [21]. So, the model system (1) with gestation delay τ is represented by
(2)
d
x
d
t
=
r
x
(
1
-
x
k
)
-
a
x
y
1
+
β
y
+
α
x
,
d
y
d
t
=
e
a
x
(
t
-
τ
)
y
1
+
β
y
(
t
-
τ
)
+
α
x
(
t
-
τ
)
-
μ
y
and this model is extensively studied by [2].Deterministic model in ecology does not usually incorporate environmental fluctuation; they are often justified by the implicit assumption that, in large populations, stochastic deviations are small enough to be ignored. Deterministic model will prove ecologically useful only if the dynamical patterns they reveal are still in evidence when stochastic effects are introduced. For terrestrial system, the environmental variability is large at both short and long time periods and could be expected to develop internal mechanisms to the system which would cope with short term variability and minimize the effects of long term variations; hence analysis of the system with white noise gives better results. Uncertain growth of populations is usually considered as an effect of environmental stochasticity. Reproduction of species depends on various factors, such as temperature, humidity, parasites and pathogens, and environmental pollution [22]. Since physical and biological environments of populations are not totally predictable, the growth of populations should be considered as a stochastic process rather than a deterministic one [23]. In spite of some shortcomings, Gaussian white noise has been proved extremely useful to model rapidly fluctuating phenomena [24, 25]. The basic mechanism and factors of population growth like the resources and vital rates birth, death, and so forth change nondeterministically due to random environment and they are the main parameters subject to coupling of the system with its environment [25–35]. Here we assume that fluctuations in the environment will manifest themselves mainly as fluctuations in the growth coefficients of the prey and predator since these are the main parameters subject to coupling of a prey-predator pair with its environment [34, 35]. So, growth terms of both prey and predator are perturbed by additive white noise [25, 33]. Therefore, model system (2) by introducing the environmental stochasticity in the form of Gaussian white noise is represented by
(3)
d
x
d
t
=
x
[
r
(
1
-
x
k
)
+
η
1
(
t
)
]
-
a
x
y
1
+
β
y
+
α
x
,
d
y
d
t
=
y
[
-
μ
+
η
2
(
t
)
]
+
e
a
x
(
t
-
τ
)
y
1
+
β
y
(
t
-
τ
)
+
α
x
(
t
-
τ
)
,
where the perturbed terms η
1
(
t
) and η
2
(
t
) are assumed to be the independent Gaussian white noise.In the next section, we perform the stochastic scenario of nondelayed system of the model system (1). We perform the stochastic scenario of delayed system (3) of the model system in Section 3. Numerical computation is performed in Section 4 and the paper ends with a brief discussion in Section 5.


## 2. Stochastic Scenario of Nondelayed System

I have so far considered an exact deterministic representation for the external force which derives the birth and death rates. In this section, I introduce the environmental stochasticity in the form of Gaussian white noise in the growth terms of both prey and predator equations. The model system (1), in this case, will be represented by the following stochastic differential equations:
(4)
d
x
d
t
=
x
[
r
(
1
-
x
k
)
+
η
1
(
t
)
]
-
a
x
y
1
+
β
y
+
α
x
,
d
y
d
t
=
y
[
-
μ
+
η
2
(
t
)
]
+
e
a
x
y
1
+
β
y
+
α
x
,
where the perturbed terms η
1
(
t
) and η
2
(
t
) are assumed to be the independent Gaussian white noise satisfying the conditions
(5)
〈
η
j
(
t
)
〉
=
0
,
〈
η
j
(
t
1
)
η
j
(
t
2
)
〉
=
ϵ
j
δ
(
t
1
-
t
2
)
for
j
=
1,2
.
Here ϵ
j
>
0
(
j
=
1,2
) are the intensities or strengths of the random perturbations, δ is the Dirac delta function defined by
(6)
δ
(
t
)
=
{
1
,
for
t
1
=
t
2
,
0
,
otherwise,
and 〈
·
〉 represents the ensemble average of the stochastic process.To study the behavior of the system (4) about the steady state E
*, we give the transformation x
′
=
ln
⁡
x, y
′
=
ln
⁡
y; x
=
u
+
x
*
,
y
=
v
+
y
* and then the system (4) reduces to the following Itô type stochastic differential equations in terms of deviation variables (
u
,
v
):
(7)
d
u
d
t
=
a
1
u
+
b
1
u
2
+
c
1
v
+
d
1
v
2
+
e
1
u
v
+
η
1
(
t
)
,
d
v
d
t
=
a
2
u
+
b
2
u
2
+
c
2
v
+
d
2
v
2
+
e
2
u
v
+
η
2
(
t
)
,
where
(8)
a
1
=
-
r
x
*
k
+
a
α
x
*
y
*
(
1
+
β
y
*
+
α
x
*
)
2
,
b
1
=
-
r
k
+
a
α
y
*
(
1
+
β
y
*
)
(
1
+
β
y
*
+
α
x
*
)
3
,
c
1
=
-
a
x
*
(
1
+
α
x
*
)
(
1
+
β
y
*
+
α
x
*
)
2
,
d
1
=
a
β
x
*
(
1
+
α
x
*
)
(
1
+
β
y
*
+
α
x
*
)
3
,
e
1
=
-
a
(
1
+
α
x
*
+
β
y
*
+
2
α
β
x
*
y
*
)
(
1
+
β
y
*
+
α
x
*
)
3
,
a
2
=
e
a
y
*
(
1
+
β
y
*
)
(
1
+
β
y
*
+
α
x
*
)
2
,
b
2
=
-
e
a
α
y
*
(
1
+
β
y
*
)
(
1
+
β
y
*
+
α
x
*
)
3
,
c
2
=
-
e
a
β
x
*
y
*
(
1
+
β
y
*
+
α
x
*
)
2
,
d
2
=
-
e
a
β
x
*
(
1
+
α
x
*
)
(
1
+
β
y
*
+
α
x
*
)
3
,
e
2
=
e
a
(
1
+
α
x
*
+
β
y
*
+
2
α
β
x
*
y
*
)
(
1
+
β
y
*
+
α
x
*
)
3
.
The solutions {
u
(
t
)
,
v
(
t
)
} of (7) subject to known initial values {
u
(
t
0
)
,
v
(
t
0
)
} determine the statistical behavior of the model system (4) near the steady state E
* at time t
>
t
0.


### 2.1. Statistical Linearization: Moment Equations

A method [24, 36, 37] for estimating the dynamical statistical properties of the solutions of nonlinear stochastic differential equations is presented in this section. The nonlinear stochastic differential equation is linearized within a small interval of the independent variable and statistical properties are expressed analytically within the interval. Long term behavior of the solution process is obtained by appropriately matching the approximate solutions at the boundaries between intervals. The method is applied to this model for which the exact time-dependent moments can be obtained by numerical methods. The calculations demonstrate that the method represents a significant improvement over the method of statistical linearization in time regimes far from equilibrium. The statistical linearization of the system (7) is represented by the following system of linear equations:
(9)
d
u
d
t
=
α
1
u
+
β
1
v
+
f
1
+
η
1
(
t
)
,
d
v
d
t
=
α
2
u
+
β
2
v
+
f
2
+
η
2
(
t
)
,
where the errors in the above linearization are given by
(10)
e
1
=
a
1
u
+
b
1
u
2
+
c
1
v
+
d
1
v
2
+
e
1
u
v
-
α
1
u
-
β
1
v
-
f
1
,
e
2
=
a
2
u
+
b
2
u
2
+
c
2
v
+
d
2
v
2
+
e
2
u
v
-
α
2
u
-
β
2
v
-
f
2
.
The unknown coefficients α
i
,
β
i, and f
i
(
i
=
1,2
) of (9) are to be determined from the minimization of the averages of the squares of errors in (10). We determine the unknown coefficients [24, 26, 27] by demanding that
(11)
∂
∂
α
i
〈
e
i
2
〉
=
∂
∂
β
i
〈
e
i
2
〉
=
∂
∂
f
i
〈
e
i
2
〉
=
0
,
i
=
1,2
.
Also, I use the following expressions [24]: I now express 〈
u
3
〉
,
〈
u
4
〉
,
〈
u
2
v
〉
,
〈
u
v
2
〉
,
〈
u
2
v
2
〉, and 〈
u
3
v
〉 in terms of the first two moments of each of the variables and the correlation coefficient using a bivariate Gaussian distribution. Since I am interested only in the first few moments, it is convenient to use the characteristic function
(12)
χ
(
ν
1
,
ν
2
)
=
exp
⁡
[
1
2
i
〈
u
〉
ν
1
+
i
〈
v
〉
ν
2
f
f
f
f
f
f
-
1
2
{
σ
1
2
ν
1
2
+
σ
2
2
ν
2
2
+
2
ρ
12
σ
1
σ
2
ν
1
ν
2
}
]
,
σ
1
2
=
〈
u
2
〉
-
〈
u
〉
2
,
σ
2
2
=
〈
v
2
〉
-
〈
v
〉
2
,
ρ
12
=
〈
u
v
〉
-
〈
u
〉
〈
v
〉
σ
1
σ
2
,
using
(13)
〈
u
n
v
m
〉
=
(
-
1
)
n
+
m
∂
n
+
m
∂
ν
1
n
∂
ν
2
m
[
χ
(
ν
1
,
ν
2
)
]
|
ν
1
=
ν
2
=
0
.
I get
(14)
〈
u
4
〉
=
3
〈
u
2
〉
2
-
2
〈
u
〉
4
,
〈
u
2
v
2
〉
=
〈
u
2
〉
〈
v
2
〉
+
2
〈
u
v
〉
2
-
2
〈
u
〉
2
〈
v
〉
2
,
〈
u
3
v
〉
=
3
〈
u
2
〉
〈
u
v
〉
-
2
〈
u
〉
3
〈
v
〉
,
〈
u
3
〉
=
3
〈
u
〉
〈
u
2
〉
-
2
〈
u
〉
3
,
〈
v
3
〉
=
3
〈
v
〉
〈
v
2
〉
-
2
〈
v
〉
3
,
〈
u
2
v
〉
=
2
〈
u
〉
〈
u
v
〉
-
2
〈
u
〉
2
〈
v
〉
+
〈
u
2
〉
〈
v
〉
,
〈
u
v
2
〉
=
2
〈
u
〉
〈
u
v
〉
-
2
〈
u
〉
〈
v
〉
2
+
〈
u
〉
〈
v
2
〉
.
Then expressions for α
i, β
i, and f
i
(
i
=
1,2
) are given by
(15)
α
i
=
a
i
+
2
b
i
〈
u
〉
+
e
i
〈
v
〉
,
β
i
=
c
i
+
2
d
i
〈
v
〉
+
e
i
〈
u
〉
,
f
i
=
b
i
(
〈
u
2
〉
-
2
〈
u
〉
2
)
+
d
i
(
〈
v
2
〉
-
2
〈
v
〉
2
)
+
e
i
(
〈
u
v
〉
-
2
〈
u
〉
〈
v
〉
)
.
The coefficients are the functions of the parameters involved with the model system and also of the different moments involving u and v. After some algebraic manipulations, I obtain the following system of first two moments:
(16)
d
〈
u
〉
d
t
=
a
1
〈
u
〉
+
b
1
〈
u
2
〉
+
c
1
〈
v
〉
+
d
1
〈
v
2
〉
+
e
1
〈
u
v
〉
,
d
〈
v
〉
d
t
=
a
2
〈
u
〉
+
b
2
〈
u
2
〉
+
c
2
〈
v
〉
+
d
2
〈
v
2
〉
+
e
2
〈
u
v
〉
,
d
〈
u
2
〉
d
t
=
2
[
a
1
〈
u
2
〉
+
b
1
〈
u
3
〉
+
c
1
〈
u
v
〉
f
f
f
f
+
d
1
〈
u
v
2
〉
+
e
1
〈
u
2
v
〉
]
+
2
ϵ
1
,
d
〈
v
2
〉
d
t
=
2
[
a
2
〈
u
v
〉
+
b
2
〈
u
2
v
〉
+
c
2
〈
v
2
〉
fffffff
+
d
2
〈
v
3
〉
+
e
2
〈
u
v
2
〉
]
+
2
ϵ
2
,
d
〈
u
v
〉
d
t
=
a
1
〈
u
v
〉
+
b
1
〈
u
2
v
〉
+
c
1
〈
v
2
〉
+
d
1
〈
v
3
〉
+
e
1
〈
u
v
2
〉
+
a
2
〈
u
2
〉
+
b
2
〈
u
3
〉
+
c
2
〈
u
v
〉
+
d
2
〈
u
v
2
〉
+
e
2
〈
u
2
v
〉
,
where I have used the relations
(17)
〈
u
η
1
〉
=
ϵ
1
,
〈
u
η
2
〉
=
〈
v
η
1
〉
=
0
,
〈
v
η
2
〉
=
ϵ
2
.
Let me now assume that the system size expansion is valid such that the correlations ϵ
i
(
i
=
1,2
) given by (17) decrease with the increase of the population size and they are assumed to be of the order of the inverse of the population size N [24, 27, 28]:
(18)
ϵ
i
∝
o
[
1
N
]
,
i
=
1,2
.
Therefore, using the expressions (14) and (17) and keeping the lowest order terms and replacing the averages 〈
u
〉 and 〈
v
〉 by their steady state values 〈
u
〉
=
〈
v
〉
=
0 [29], I get the following reduced equations for second order moments:
(19)
[
D
-
2
a
1
]
〈
u
2
〉
=
2
c
1
〈
u
v
〉
,
[
D
-
2
c
2
]
〈
v
2
〉
=
2
a
2
〈
u
v
〉
,
[
D
-
a
1
-
c
2
]
〈
u
v
〉
=
a
2
〈
u
2
〉
+
c
1
〈
v
2
〉
,
where D stands for the operator d
/
d
t.


### 2.2. Nonequilibrium Fluctuation and Stability Analysis

Eliminating〈
u
2
〉 and 〈
v
2
〉 from the equations of (19), I get the following third order linear ordinary differential equation in 〈
u
v
〉:
(20)
[
D
3
+
3
A
D
2
+
3
B
D
+
C
]
〈
u
v
〉
=
0
.
Let 〈
u
v
〉
=
e
m
t be a trial solution of (20) and the auxiliary equation is given by
(21)
m
3
+
3
A
m
2
+
3
B
m
+
C
=
0
,
where
(22)
A
=
-
(
a
1
+
c
2
)
,
B
=
2
3
{
(
a
1
+
c
2
)
2
+
2
(
a
1
c
2
-
a
2
c
1
)
}
,
C
=
-
4
(
a
1
+
c
2
)
(
a
1
c
2
-
a
2
c
1
)
.
Let H
=
A
2
-
B. Then the nature and structure of the roots of (21) will solely be determined by the quantities A and H. I discuss the following two cases.Case 1 (H
<
0).
In this case, roots of (21) are given by
(23)
m
1
=
-
A
,
m
2,3
=
-
A
±
i
3
H
0
,
where
H
0
=
-
H
(
>
0
)
.
The solutions of the linear system (20) are then given by
(24)
〈
u
v
〉
=
A
1
e
-
A
t
+
e
-
A
t
[
A
2
cos
⁡
(
3
H
0
t
)
f
f
f
f
f
f
f
f
f
f
f
+
A
3
sin
(
3
H
0
t
)
]
,
〈
v
2
〉
=
B
1
e
-
A
t
+
e
-
A
t
[
B
2
cos
⁡
(
3
H
0
t
)
ffffffffffffffffffffff
+
B
3
sin
(
3
H
0
t
)
]
+
P
1
e
2
a
1
t
,
〈
u
2
〉
=
C
1
e
-
A
t
+
e
-
A
t
[
C
2
cos
⁡
(
3
H
0
t
)
ffffffffffffffffffffff
+
C
3
sin
(
3
H
0
t
)
]
+
P
2
e
2
c
2
t
,
where A
i, B
i, C
i, (
i
=
1,2
,
3
), P
1, and P
2 are constants. Thus, each of 〈
u
v
〉
,
〈
u
2
〉
,
〈
v
2
〉, given by (24), converges with increasing time if A
>
0, that is, if a
1
+
c
2
<
0, depicting the stochastically stable situation of the system (4) in the sense of second order moments [38]. On the other hand, if a
1
+
c
2
>
0, then A
<
0 and the system is stochastically unstable.Case 2 (H
>
0).
In this case, roots of (21) are given by
(25)
m
1
=
-
A
,
m
2,3
=
-
A
±
3
H
.
The solutions of the linear system (20) are then given by
(26)
〈
u
v
〉
=
F
1
e
-
A
t
+
F
2
e
(
-
A
+
3
H
)
t
+
F
3
e
(
-
A
-
3
H
)
t
,
〈
v
2
〉
=
G
1
e
-
A
t
+
G
2
e
(
-
A
+
3
H
)
t
+
G
3
e
(
-
A
-
3
H
)
t
+
Q
1
e
2
a
1
t
,
〈
u
2
〉
=
K
1
e
-
A
t
+
K
2
e
(
-
A
+
3
H
)
t
+
K
3
e
(
-
A
-
3
H
)
t
+
Q
2
e
2
c
2
t
,
where F
i, G
i, K
i, (
i
=
1,2
,
3
), Q
1, and Q
2 are constants. In this case, we observe that each of 〈
u
v
〉
,
〈
u
2
〉
,
〈
v
2
〉, given by (26), converges with increasing time when A
>
0 (i.e., a
1
+
c
2
<
0) and 3
H
<
A hold simultaneously and, then system will said to be stochastically stable in the sense of second order moment. The system, on the other hand, is unstable when a
1
+
c
2
>
0 or 3
H
<
A or both of them hold. Thus, exchange of stability occurs when a
1
+
c
2
=
0 and 3
H
=
A.


## 2.1. Statistical Linearization: Moment Equations

A method [24, 36, 37] for estimating the dynamical statistical properties of the solutions of nonlinear stochastic differential equations is presented in this section. The nonlinear stochastic differential equation is linearized within a small interval of the independent variable and statistical properties are expressed analytically within the interval. Long term behavior of the solution process is obtained by appropriately matching the approximate solutions at the boundaries between intervals. The method is applied to this model for which the exact time-dependent moments can be obtained by numerical methods. The calculations demonstrate that the method represents a significant improvement over the method of statistical linearization in time regimes far from equilibrium. The statistical linearization of the system (7) is represented by the following system of linear equations:
(9)
d
u
d
t
=
α
1
u
+
β
1
v
+
f
1
+
η
1
(
t
)
,
d
v
d
t
=
α
2
u
+
β
2
v
+
f
2
+
η
2
(
t
)
,
where the errors in the above linearization are given by
(10)
e
1
=
a
1
u
+
b
1
u
2
+
c
1
v
+
d
1
v
2
+
e
1
u
v
-
α
1
u
-
β
1
v
-
f
1
,
e
2
=
a
2
u
+
b
2
u
2
+
c
2
v
+
d
2
v
2
+
e
2
u
v
-
α
2
u
-
β
2
v
-
f
2
.
The unknown coefficients α
i
,
β
i, and f
i
(
i
=
1,2
) of (9) are to be determined from the minimization of the averages of the squares of errors in (10). We determine the unknown coefficients [24, 26, 27] by demanding that
(11)
∂
∂
α
i
〈
e
i
2
〉
=
∂
∂
β
i
〈
e
i
2
〉
=
∂
∂
f
i
〈
e
i
2
〉
=
0
,
i
=
1,2
.
Also, I use the following expressions [24]: I now express 〈
u
3
〉
,
〈
u
4
〉
,
〈
u
2
v
〉
,
〈
u
v
2
〉
,
〈
u
2
v
2
〉, and 〈
u
3
v
〉 in terms of the first two moments of each of the variables and the correlation coefficient using a bivariate Gaussian distribution. Since I am interested only in the first few moments, it is convenient to use the characteristic function
(12)
χ
(
ν
1
,
ν
2
)
=
exp
⁡
[
1
2
i
〈
u
〉
ν
1
+
i
〈
v
〉
ν
2
f
f
f
f
f
f
-
1
2
{
σ
1
2
ν
1
2
+
σ
2
2
ν
2
2
+
2
ρ
12
σ
1
σ
2
ν
1
ν
2
}
]
,
σ
1
2
=
〈
u
2
〉
-
〈
u
〉
2
,
σ
2
2
=
〈
v
2
〉
-
〈
v
〉
2
,
ρ
12
=
〈
u
v
〉
-
〈
u
〉
〈
v
〉
σ
1
σ
2
,
using
(13)
〈
u
n
v
m
〉
=
(
-
1
)
n
+
m
∂
n
+
m
∂
ν
1
n
∂
ν
2
m
[
χ
(
ν
1
,
ν
2
)
]
|
ν
1
=
ν
2
=
0
.
I get
(14)
〈
u
4
〉
=
3
〈
u
2
〉
2
-
2
〈
u
〉
4
,
〈
u
2
v
2
〉
=
〈
u
2
〉
〈
v
2
〉
+
2
〈
u
v
〉
2
-
2
〈
u
〉
2
〈
v
〉
2
,
〈
u
3
v
〉
=
3
〈
u
2
〉
〈
u
v
〉
-
2
〈
u
〉
3
〈
v
〉
,
〈
u
3
〉
=
3
〈
u
〉
〈
u
2
〉
-
2
〈
u
〉
3
,
〈
v
3
〉
=
3
〈
v
〉
〈
v
2
〉
-
2
〈
v
〉
3
,
〈
u
2
v
〉
=
2
〈
u
〉
〈
u
v
〉
-
2
〈
u
〉
2
〈
v
〉
+
〈
u
2
〉
〈
v
〉
,
〈
u
v
2
〉
=
2
〈
u
〉
〈
u
v
〉
-
2
〈
u
〉
〈
v
〉
2
+
〈
u
〉
〈
v
2
〉
.
Then expressions for α
i, β
i, and f
i
(
i
=
1,2
) are given by
(15)
α
i
=
a
i
+
2
b
i
〈
u
〉
+
e
i
〈
v
〉
,
β
i
=
c
i
+
2
d
i
〈
v
〉
+
e
i
〈
u
〉
,
f
i
=
b
i
(
〈
u
2
〉
-
2
〈
u
〉
2
)
+
d
i
(
〈
v
2
〉
-
2
〈
v
〉
2
)
+
e
i
(
〈
u
v
〉
-
2
〈
u
〉
〈
v
〉
)
.
The coefficients are the functions of the parameters involved with the model system and also of the different moments involving u and v. After some algebraic manipulations, I obtain the following system of first two moments:
(16)
d
〈
u
〉
d
t
=
a
1
〈
u
〉
+
b
1
〈
u
2
〉
+
c
1
〈
v
〉
+
d
1
〈
v
2
〉
+
e
1
〈
u
v
〉
,
d
〈
v
〉
d
t
=
a
2
〈
u
〉
+
b
2
〈
u
2
〉
+
c
2
〈
v
〉
+
d
2
〈
v
2
〉
+
e
2
〈
u
v
〉
,
d
〈
u
2
〉
d
t
=
2
[
a
1
〈
u
2
〉
+
b
1
〈
u
3
〉
+
c
1
〈
u
v
〉
f
f
f
f
+
d
1
〈
u
v
2
〉
+
e
1
〈
u
2
v
〉
]
+
2
ϵ
1
,
d
〈
v
2
〉
d
t
=
2
[
a
2
〈
u
v
〉
+
b
2
〈
u
2
v
〉
+
c
2
〈
v
2
〉
fffffff
+
d
2
〈
v
3
〉
+
e
2
〈
u
v
2
〉
]
+
2
ϵ
2
,
d
〈
u
v
〉
d
t
=
a
1
〈
u
v
〉
+
b
1
〈
u
2
v
〉
+
c
1
〈
v
2
〉
+
d
1
〈
v
3
〉
+
e
1
〈
u
v
2
〉
+
a
2
〈
u
2
〉
+
b
2
〈
u
3
〉
+
c
2
〈
u
v
〉
+
d
2
〈
u
v
2
〉
+
e
2
〈
u
2
v
〉
,
where I have used the relations
(17)
〈
u
η
1
〉
=
ϵ
1
,
〈
u
η
2
〉
=
〈
v
η
1
〉
=
0
,
〈
v
η
2
〉
=
ϵ
2
.
Let me now assume that the system size expansion is valid such that the correlations ϵ
i
(
i
=
1,2
) given by (17) decrease with the increase of the population size and they are assumed to be of the order of the inverse of the population size N [24, 27, 28]:
(18)
ϵ
i
∝
o
[
1
N
]
,
i
=
1,2
.
Therefore, using the expressions (14) and (17) and keeping the lowest order terms and replacing the averages 〈
u
〉 and 〈
v
〉 by their steady state values 〈
u
〉
=
〈
v
〉
=
0 [29], I get the following reduced equations for second order moments:
(19)
[
D
-
2
a
1
]
〈
u
2
〉
=
2
c
1
〈
u
v
〉
,
[
D
-
2
c
2
]
〈
v
2
〉
=
2
a
2
〈
u
v
〉
,
[
D
-
a
1
-
c
2
]
〈
u
v
〉
=
a
2
〈
u
2
〉
+
c
1
〈
v
2
〉
,
where D stands for the operator d
/
d
t.


## 2.2. Nonequilibrium Fluctuation and Stability Analysis

Eliminating〈
u
2
〉 and 〈
v
2
〉 from the equations of (19), I get the following third order linear ordinary differential equation in 〈
u
v
〉:
(20)
[
D
3
+
3
A
D
2
+
3
B
D
+
C
]
〈
u
v
〉
=
0
.
Let 〈
u
v
〉
=
e
m
t be a trial solution of (20) and the auxiliary equation is given by
(21)
m
3
+
3
A
m
2
+
3
B
m
+
C
=
0
,
where
(22)
A
=
-
(
a
1
+
c
2
)
,
B
=
2
3
{
(
a
1
+
c
2
)
2
+
2
(
a
1
c
2
-
a
2
c
1
)
}
,
C
=
-
4
(
a
1
+
c
2
)
(
a
1
c
2
-
a
2
c
1
)
.
Let H
=
A
2
-
B. Then the nature and structure of the roots of (21) will solely be determined by the quantities A and H. I discuss the following two cases.Case 1 (H
<
0).
In this case, roots of (21) are given by
(23)
m
1
=
-
A
,
m
2,3
=
-
A
±
i
3
H
0
,
where
H
0
=
-
H
(
>
0
)
.
The solutions of the linear system (20) are then given by
(24)
〈
u
v
〉
=
A
1
e
-
A
t
+
e
-
A
t
[
A
2
cos
⁡
(
3
H
0
t
)
f
f
f
f
f
f
f
f
f
f
f
+
A
3
sin
(
3
H
0
t
)
]
,
〈
v
2
〉
=
B
1
e
-
A
t
+
e
-
A
t
[
B
2
cos
⁡
(
3
H
0
t
)
ffffffffffffffffffffff
+
B
3
sin
(
3
H
0
t
)
]
+
P
1
e
2
a
1
t
,
〈
u
2
〉
=
C
1
e
-
A
t
+
e
-
A
t
[
C
2
cos
⁡
(
3
H
0
t
)
ffffffffffffffffffffff
+
C
3
sin
(
3
H
0
t
)
]
+
P
2
e
2
c
2
t
,
where A
i, B
i, C
i, (
i
=
1,2
,
3
), P
1, and P
2 are constants. Thus, each of 〈
u
v
〉
,
〈
u
2
〉
,
〈
v
2
〉, given by (24), converges with increasing time if A
>
0, that is, if a
1
+
c
2
<
0, depicting the stochastically stable situation of the system (4) in the sense of second order moments [38]. On the other hand, if a
1
+
c
2
>
0, then A
<
0 and the system is stochastically unstable.Case 2 (H
>
0).
In this case, roots of (21) are given by
(25)
m
1
=
-
A
,
m
2,3
=
-
A
±
3
H
.
The solutions of the linear system (20) are then given by
(26)
〈
u
v
〉
=
F
1
e
-
A
t
+
F
2
e
(
-
A
+
3
H
)
t
+
F
3
e
(
-
A
-
3
H
)
t
,
〈
v
2
〉
=
G
1
e
-
A
t
+
G
2
e
(
-
A
+
3
H
)
t
+
G
3
e
(
-
A
-
3
H
)
t
+
Q
1
e
2
a
1
t
,
〈
u
2
〉
=
K
1
e
-
A
t
+
K
2
e
(
-
A
+
3
H
)
t
+
K
3
e
(
-
A
-
3
H
)
t
+
Q
2
e
2
c
2
t
,
where F
i, G
i, K
i, (
i
=
1,2
,
3
), Q
1, and Q
2 are constants. In this case, we observe that each of 〈
u
v
〉
,
〈
u
2
〉
,
〈
v
2
〉, given by (26), converges with increasing time when A
>
0 (i.e., a
1
+
c
2
<
0) and 3
H
<
A hold simultaneously and, then system will said to be stochastically stable in the sense of second order moment. The system, on the other hand, is unstable when a
1
+
c
2
>
0 or 3
H
<
A or both of them hold. Thus, exchange of stability occurs when a
1
+
c
2
=
0 and 3
H
=
A.


## 3. Stochastic Scenario of Delayed System

Consider(27)
d
x
d
t
=
x
[
r
(
1
-
x
k
)
+
η
1
(
t
)
]
-
a
x
y
1
+
β
y
+
α
x
,
d
y
d
t
=
y
[
-
μ
+
η
2
(
t
)
]
+
e
a
x
(
t
-
τ
)
y
1
+
β
y
(
t
-
τ
)
+
α
x
(
t
-
τ
)
.
Again using the transformations x
′
=
ln
⁡
x, y
′
=
ln
⁡
y, and x
=
u
+
x
*
,
y
=
v
+
y
* and assuming the delay to be very small, the system (27) (to a first approximation) can be written as
(28)
d
u
d
t
=
a
1
u
+
c
1
v
+
η
1
,
d
v
d
t
=
a
2
u
(
t
-
τ
)
+
c
2
v
(
t
-
τ
)
+
η
2
,
where a
i and c
i
(
i
=
1,2
) are given by (8).


### 3.1. Fourier Transforms: Spectral Density

Taking Fourier transform of both sides of each of the equations in (28) and following [20, 26], I obtain
(29)
η
-
1
(
s
)
=
i
s
u
-
(
s
)
-
a
1
u
-
(
s
)
-
c
1
v
-
(
s
)
,
η
-
2
(
s
)
=
i
s
v
-
(
s
)
-
a
2
u
-
(
s
)
e
-
i
s
τ
-
c
2
v
-
(
s
)
e
-
i
s
τ
,
where n
-
(
s
)
=
∫
-
∞
+
∞
‍
n
(
t
)
e
-
i
s
t
d
t.The system of equations (29) can be written in matrix form as
(30)
A
X
=
B
,
where
(31)
A
=
(
-
a
1
+
i
s
-
c
1
-
a
2
e
-
i
s
τ
-
c
2
e
-
i
s
τ
+
i
s
)
,
X
=
(
u
-
(
s
)
v
-
(
s
)
)
,
B
=
(
η
-
1
η
-
2
)
.
Now
(32)
M
=
det
⁡
A
=
{
(
a
1
c
2
-
a
2
c
1
)
cos
⁡
(
s
τ
)
-
s
c
2
sin
(
s
τ
)
-
s
2
}
-
i
{
(
a
1
c
2
-
a
2
c
1
)
sin
⁡
(
s
τ
)
+
s
c
2
cos
⁡
⁡
(
s
τ
)
+
a
1
s
}
.
I assume that A
-
1 exists. Then I have A
-
1
=
(
a
i
j
)
2
×
2, where
(33)
a
11
=
-
c
2
cos
⁡
(
s
τ
)
+
i
{
s
+
c
2
sin
(
s
τ
)
}
M
,
a
12
=
c
1
M
,
a
21
=
a
2
(
cos
⁡
(
s
τ
)
-
i
sin
(
s
τ
)
)
M
,
a
22
=
-
a
1
+
i
s
M
.
Then the solution of (30) can be written as
(34)
u
-
(
s
)
=
∑
j
=
1
2
a
1
j
η
j
,
v
-
(
s
)
=
∑
j
=
1
2
a
1
j
η
j
.
Now following [20, 26] and using (34), the spectral density of u is given by
(35)
S
u
(
ω
)
=
lim
⁡
T
→
∞
1
T
∫
∫
-
T
/
2
T
/
2
〈
u
(
t
)
u
(
t
′
)
〉
exp
⁡
⁡
{
i
ω
(
t
′
-
t
)
}
d
t
d
t
′
=
∑
j
=
1
2
|
a
1
j
|
2
S
η
j
(
ω
)
.
Similarly the spectral density of v is given by
(36)
S
v
=
∑
j
=
1
2
|
a
2
j
|
2
S
η
j
(
ω
)
.
Therefore the fluctuation intensity (variance) of u is given by
(37)
σ
u
2
=
1
2
π
∫
-
∞
+
∞
S
u
(
ω
)
d
ω
=
1
2
π
∫
-
∞
+
∞
‍
∑
j
=
1
2
|
a
1
j
|
2
S
η
j
(
ω
)
d
ω
=
1
2
π
∫
-
∞
+
∞
∑
j
=
1
2
|
a
1
j
|
2
d
ω
,
since S
η
j
(
ω
)
=
1.Similarly the fluctuation intensity ofv is given by
(38)
σ
v
2
=
1
2
π
∫
-
∞
+
∞
∑
j
=
1
2
|
a
2
j
|
2
d
ω
.Letτ be very small and the following two conditions hold:(i)
(
a
1
c
2
-
a
2
c
1
)
<
min
⁡
{
(
c
1
2
+
c
2
2
)
/
(
1
+
c
2
τ
)
,
(
1
+
c
2
τ
)
(
a
1
2
+
a
2
2
)
},
(ii)
4
(
1
+
c
2
τ
)
(
a
1
c
2
-
a
2
c
1
)
<
{
(
a
1
c
2
-
a
2
c
1
)
τ
-
a
1
-
c
2
},where (
a
1
c
2
-
a
2
c
1
)
=
(
e
a
x
*
y
*
/
(
1
+
β
y
*
+
α
x
*
)
2
)
[
r
β
x
*
/
K
+
a
/
(
1
+
β
y
*
+
α
x
*
)
]
>
0.Then after some calculation, I obtain the fluctuation intensities ofu and v as
(39)
σ
u
2
=
c
1
2
+
c
2
2
-
(
a
1
c
2
-
a
2
c
1
)
(
1
+
c
2
τ
)
2
ℏ
(
a
1
c
2
-
a
2
c
1
)
,
σ
v
2
=
(
1
+
c
2
τ
)
(
a
1
2
+
a
2
2
)
-
(
a
1
c
2
-
a
2
c
1
)
2
(
1
+
c
2
τ
)
ℏ
(
a
1
c
2
-
a
2
c
1
)
,
where ℏ
=
[
{
(
a
1
c
2
-
a
2
c
1
)
τ
-
a
1
-
c
2
}
2
-
4
(
1
+
c
2
τ
)
(
a
1
c
2
-
a
2
c
1
)
]
1
/
2.Let(40)
Ω
1
=
c
1
2
+
c
2
2
-
(
a
1
c
2
-
a
2
c
1
)
(
1
+
c
2
τ
)
2
(
a
1
c
2
-
a
2
c
1
)
,
Ω
2
=
(
1
+
c
2
τ
)
(
a
1
2
+
a
2
2
)
-
(
a
1
c
2
-
a
2
c
1
)
2
(
1
+
c
2
τ
)
(
a
1
c
2
-
a
2
c
1
)
.
Following the criteria of stability in the stochastic environment [30], it is seen that the deterministic stability criteria are not enough to determine the stability of the system in a rapidly fluctuating environment. ℏ
≪
Ω
i
,
(
i
=
1,2
), in rapidly fluctuating environment; then populations exhibit abnormally large fluctuations which rapidly lead to extinction. In the intermediate region where ℏ and Ω
i
,
(
i
=
1,2
), are commensurate, the populations are likely to undergo significant fluctuations, even though they persist for long time. Finally, if ℏ
≫
Ω
i
,
(
i
=
1,2
), population fluctuations are relatively small and the environment is effectively deterministic. These results are in good agreement with those of [25, 30–33].


## 3.1. Fourier Transforms: Spectral Density

Taking Fourier transform of both sides of each of the equations in (28) and following [20, 26], I obtain
(29)
η
-
1
(
s
)
=
i
s
u
-
(
s
)
-
a
1
u
-
(
s
)
-
c
1
v
-
(
s
)
,
η
-
2
(
s
)
=
i
s
v
-
(
s
)
-
a
2
u
-
(
s
)
e
-
i
s
τ
-
c
2
v
-
(
s
)
e
-
i
s
τ
,
where n
-
(
s
)
=
∫
-
∞
+
∞
‍
n
(
t
)
e
-
i
s
t
d
t.The system of equations (29) can be written in matrix form as
(30)
A
X
=
B
,
where
(31)
A
=
(
-
a
1
+
i
s
-
c
1
-
a
2
e
-
i
s
τ
-
c
2
e
-
i
s
τ
+
i
s
)
,
X
=
(
u
-
(
s
)
v
-
(
s
)
)
,
B
=
(
η
-
1
η
-
2
)
.
Now
(32)
M
=
det
⁡
A
=
{
(
a
1
c
2
-
a
2
c
1
)
cos
⁡
(
s
τ
)
-
s
c
2
sin
(
s
τ
)
-
s
2
}
-
i
{
(
a
1
c
2
-
a
2
c
1
)
sin
⁡
(
s
τ
)
+
s
c
2
cos
⁡
⁡
(
s
τ
)
+
a
1
s
}
.
I assume that A
-
1 exists. Then I have A
-
1
=
(
a
i
j
)
2
×
2, where
(33)
a
11
=
-
c
2
cos
⁡
(
s
τ
)
+
i
{
s
+
c
2
sin
(
s
τ
)
}
M
,
a
12
=
c
1
M
,
a
21
=
a
2
(
cos
⁡
(
s
τ
)
-
i
sin
(
s
τ
)
)
M
,
a
22
=
-
a
1
+
i
s
M
.
Then the solution of (30) can be written as
(34)
u
-
(
s
)
=
∑
j
=
1
2
a
1
j
η
j
,
v
-
(
s
)
=
∑
j
=
1
2
a
1
j
η
j
.
Now following [20, 26] and using (34), the spectral density of u is given by
(35)
S
u
(
ω
)
=
lim
⁡
T
→
∞
1
T
∫
∫
-
T
/
2
T
/
2
〈
u
(
t
)
u
(
t
′
)
〉
exp
⁡
⁡
{
i
ω
(
t
′
-
t
)
}
d
t
d
t
′
=
∑
j
=
1
2
|
a
1
j
|
2
S
η
j
(
ω
)
.
Similarly the spectral density of v is given by
(36)
S
v
=
∑
j
=
1
2
|
a
2
j
|
2
S
η
j
(
ω
)
.
Therefore the fluctuation intensity (variance) of u is given by
(37)
σ
u
2
=
1
2
π
∫
-
∞
+
∞
S
u
(
ω
)
d
ω
=
1
2
π
∫
-
∞
+
∞
‍
∑
j
=
1
2
|
a
1
j
|
2
S
η
j
(
ω
)
d
ω
=
1
2
π
∫
-
∞
+
∞
∑
j
=
1
2
|
a
1
j
|
2
d
ω
,
since S
η
j
(
ω
)
=
1.Similarly the fluctuation intensity ofv is given by
(38)
σ
v
2
=
1
2
π
∫
-
∞
+
∞
∑
j
=
1
2
|
a
2
j
|
2
d
ω
.Letτ be very small and the following two conditions hold:(i)
(
a
1
c
2
-
a
2
c
1
)
<
min
⁡
{
(
c
1
2
+
c
2
2
)
/
(
1
+
c
2
τ
)
,
(
1
+
c
2
τ
)
(
a
1
2
+
a
2
2
)
},
(ii)
4
(
1
+
c
2
τ
)
(
a
1
c
2
-
a
2
c
1
)
<
{
(
a
1
c
2
-
a
2
c
1
)
τ
-
a
1
-
c
2
},where (
a
1
c
2
-
a
2
c
1
)
=
(
e
a
x
*
y
*
/
(
1
+
β
y
*
+
α
x
*
)
2
)
[
r
β
x
*
/
K
+
a
/
(
1
+
β
y
*
+
α
x
*
)
]
>
0.Then after some calculation, I obtain the fluctuation intensities ofu and v as
(39)
σ
u
2
=
c
1
2
+
c
2
2
-
(
a
1
c
2
-
a
2
c
1
)
(
1
+
c
2
τ
)
2
ℏ
(
a
1
c
2
-
a
2
c
1
)
,
σ
v
2
=
(
1
+
c
2
τ
)
(
a
1
2
+
a
2
2
)
-
(
a
1
c
2
-
a
2
c
1
)
2
(
1
+
c
2
τ
)
ℏ
(
a
1
c
2
-
a
2
c
1
)
,
where ℏ
=
[
{
(
a
1
c
2
-
a
2
c
1
)
τ
-
a
1
-
c
2
}
2
-
4
(
1
+
c
2
τ
)
(
a
1
c
2
-
a
2
c
1
)
]
1
/
2.Let(40)
Ω
1
=
c
1
2
+
c
2
2
-
(
a
1
c
2
-
a
2
c
1
)
(
1
+
c
2
τ
)
2
(
a
1
c
2
-
a
2
c
1
)
,
Ω
2
=
(
1
+
c
2
τ
)
(
a
1
2
+
a
2
2
)
-
(
a
1
c
2
-
a
2
c
1
)
2
(
1
+
c
2
τ
)
(
a
1
c
2
-
a
2
c
1
)
.
Following the criteria of stability in the stochastic environment [30], it is seen that the deterministic stability criteria are not enough to determine the stability of the system in a rapidly fluctuating environment. ℏ
≪
Ω
i
,
(
i
=
1,2
), in rapidly fluctuating environment; then populations exhibit abnormally large fluctuations which rapidly lead to extinction. In the intermediate region where ℏ and Ω
i
,
(
i
=
1,2
), are commensurate, the populations are likely to undergo significant fluctuations, even though they persist for long time. Finally, if ℏ
≫
Ω
i
,
(
i
=
1,2
), population fluctuations are relatively small and the environment is effectively deterministic. These results are in good agreement with those of [25, 30–33].


## 4. Numerical Simulation

In the numerical experiments of the stochastic system, I consider the set of parameter values taken asr
=
0.4, k
=
5, a
=
0.8, α
=
0.2, β
=
1.5, e
=
0.6, and μ
=
0.2 and approximate the solutions of the system (4) (similarly for the Ito SDEs (7)) by Euler-Maruyama method. For c
=
0.45, we compute the values of A and H as 0.4222 (>0) and −0.1801 (<0), respectively. Thus, conditions in Case 1 are satisfied and the coexisting equilibrium point E
* is therefore stable in the sense of second order moment (Figure 1(a)). It means that the second moments of the system state vectors approach zero asymptotically. At a
=
0.1, the above quantities become A
=
-
0.1475 (<0) and H
=
-
0.4485 (<0) and the coexisting equilibrium point E
* is unstable (Figure 1(b)). The time evolutions and phase plane diagrams of the system (4) for the above two cases are shown in Figure 2. If I choose μ
=
0.2
,
β
=
1.5
,
a
=
0.8, then H
=
0.8878 (>0), A
=
1.7723 (>0), and 3
H = 1.632 (>0). Thus, all conditions of Case 2 are satisfied and the system (4) becomes stable in the sense of second order moments (Figure 3(a)). If I consider μ
=
0.2
,
β
=
0.5
,
a
=
0.1, keeping other parameter values unchanged, then H
=
0.4284 (>0), A
=
-
1.2979 (<0), and 3
H
=
1.1336 (>0) required for the stability are violated. Therefore, the second order moments diverge with increasing time (Figure 3(b)) and E
* becomes unstable. The corresponding time evolutions and phase plane diagrams are shown in Figure 4.Behavior of the model system (4) in the sense of second order moments: (a) stable behavior for a
=
0.8 and (b) unstable behavior for a
=
0.1. Other parameters are as in text.
(a)
H
<
0, A
>
0
(b)
H
<
0, A
<
0Figure 2
(a)-(c) and (b)-(d) depict, respectively, the time evolution and corresponding phase trajectory of the stochastic model system (4) for a
=
0.8 and a
=
0.1. (a-b) describes the stable behavior and (c-d) describes the unstable behavior of the system (4). The rest of the parameters are as in text.Behavior of the model system (4) in the sense of second order moments: (a) stable behavior for μ
=
0.2, β
=
1.5, and a
=
0.8 and (b) unstable behavior for μ
=
0.2, β
=
0.5, and a
=
0.1. Other parameters are as in text.
(a)
H
>
0
,
A
>
0
,
(
3
H
)
1
/
2
<
A
(b)
H
>
0
,
A
<
0
,
(
3
H
)
1
/
2
>
AFigure 4
(a)-(c) and (b)-(d) depict, respectively, the time evolution and corresponding phase trajectory of the stochastic model system (4) for μ
=
0.2, β
=
1.5, and a
=
0.8 and μ
=
0.2, β
=
0.5, and a
=
0.1. (a-b) describes the stable behavior and (c-d) describes the unstable behavior of the system (4). The rest of the parameters are as in text.


## 5. Discussion

Understanding the relationship between predator and prey is central goal in ecology. Predators rate of feeding upon prey is one of the most significant components in the study of predator-prey relationship. Predator's interference in prey-predator interaction is supposed to have significant role in the stability of their interaction. Also the time lag between prey capture and its corresponding positive feedback to predator's growth rate is a very significant part to realize the system. The purpose of this work is to observe the extent to which predator's interference and gestation delay drive the population dynamics of a predator-prey interaction under fluctuating environment. Results show that predator's interference plays a significant role in changing the stochastic stability of the system. To study the effect of environmental fluctuation on the time-delayed predator-prey system (27), I have superimposed Gaussian white noises on (2) and then studied nonequilibrium fluctuation and stability of the resulting stochastic model (27) by using Fourier transform technique. Following the criteria of stability in the stochastic environment [30], it is seen that the environmental noises have a destabilizing effect on the system when ℏ
≪
Ω
i
(
i
=
1,2
). Also the deterministic system and the noise-induced stochastic system may behave alike with respect to stability when ℏ
≪
Ω
i
(
i
=
1,2
). Further, in the intermediate region where ℏ and Ω
i
(
i
=
1,2
) are commensurate, the populations are likely to undergo significant fluctuations, even though they persist for long time. It is well known that natural populations of plants and animals neither increase indefinitely to blanket the world nor become extinct (except in some rare cases due to some rare reasons). Hence, in practice, we often want to keep the prey and predator population to an acceptable level in finite time. In order to accomplish this, we strongly suggest that, in realistic field situations (where effect of time-delay and environmental fluctuation can never be violated), the parameters of the system should be regulated in such a way that E
* is deterministically stable and ℏ
≪
Ω
i
(
i
=
1,2
).



---

*Source: 501836-2014-04-08.xml*