Next Article in Journal
Equivalence of Competitive Equilibria, Fuzzy Cores, and Fuzzy Bargaining Sets in Finite Production Economies
Next Article in Special Issue
Modeling and Stability Analysis of Within-Host IAV/SARS-CoV-2 Coinfection with Antibody Immunity
Previous Article in Journal
Soliton Solutions of Klein–Fock–Gordon Equation Using Sardar Subequation Method
Previous Article in Special Issue
Leveraging Geographically Distributed Data for Influenza and SARS-CoV-2 Non-Parametric Forecasting
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Predator–Prey Model with Beddington–DeAngelis Functional Response and Multiple Delays in Deterministic and Stochastic Environments

1
School of Science, Guilin University of Technology, Guilin 541004, China
2
School of Teacher Education, Qujing Normal University, Qujing 655011, China
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(18), 3378; https://doi.org/10.3390/math10183378
Submission received: 21 August 2022 / Revised: 12 September 2022 / Accepted: 14 September 2022 / Published: 17 September 2022

Abstract

:
In view of prey’s delayed fear due to predators, delayed predator gestation, and the significance of intra-specific competition between predators when their populations are sufficiently large, a prey–predator population model with a density-dependent functional response is established in a deterministic environment. We research the existence and asymptotic stability of the equilibrium statuses. Then, taking into consideration environmental disturbances, we extend the deterministic model to a stochastic model and research the existence and stationary distributions of stochastic solutions. Finally, we perform some numerical simulations to verify the theoretical results. Numerical examples indicate that fear, delays and environmental disturbance play crucial roles in the system stability of the equilibrium status.

1. Introduction

In natural ecosystems, the classical Lotka–Volterra model has been introduced to describe interactions between species [1]. One of the most important species interaction is predation. By incorporating logistic growth for a prey population and the predator population’s functional response (predation connecting the predator–prey interactions), the Lotka–Volterra model was improved to a prey–predator population model. The functional response formulates the feeding rate per predator on the prey population. This is a very key portion for researchers to investigate the kinetics of the model in a more realistic way. Various functional responses have been proposed to formulate realistic behavior among predator and prey species [2,3,4,5]. A Holling-type functional response describes well the predation speed within a regular range. However, in the real world, due to mutual disturbances among predators, an increase in predator population density may indicate a reduction in predation speed. Intra-specific competition within predator populations cannot be ignored when the predator population is significantly large. This is a crucial factor to modify the Holling-type response to a predator-dependent functional response [6,7,8]. By comparing 19 predator–prey systems, the authors of [6] found that when prey is abundant, the Beddington–DeAngelis functional response can preferably describe the predator population’s feeding behaviors.
Functional response formulates the direct feeding of predators, which is easy to find in nature, so a large number of dynamic studies have been performed, and many nice results have been obtained [1,2,3,4,5,6,7,8]. However, in the presence of a mighty predator, the prey species may change its behaviors more powerfully than due to the direct impact of killing, which has been verified by both theoretical studies and empirical data. Experimental facts show that almost all prey population react to the risk of predation with a variety of anti-predation behaviors resulting in a series of physiological changes such as vigilance, foraging and changing habitats. These anti-predation responses may reduce the long-term cost of reproduction and increase the probability of adult survival [9]. For example, an experiment on songbirds was executed in 2011 by the authors of [10]. They observed that the birthrate and offspring survival rate of songbirds were both affected by fear of predators, and the number of their offspring declined by 40%. Incorporating the cost of fear in the birthrate of prey, Wang et al. [11] proposed a prey–predator model and found that a large amount of fear could destabilize the system, accompanied by oscillations. They suggested that the death rate of prey is also influenced by fear of the predator. Further, the findings of Elliot et al. [12] also revealed that the intrinsic growth rate of prey declines in response to fear of the predator, with a reduction from 4.2 to 3.6.
Consequently, in the study of predator–prey dynamics, the fear induced by predators should be incorporated.
In reality, when prey perceives risk from a predator by vocal or chemical cues, they take some time to estimate the risk and decide whether it is necessary to present some anti-predation behavior or not, instead of having an instantaneous response. That is, there is always a time lag for the response of prey to the predation risk [13,14,15]. Actually, delay occurs in almost all biological activities, such as gestation, digestion, maturation and so on. For example, in prey–predator system, the predator’s consumption of prey cannot be instantaneous, and there must be a time delay for this conversion, which is called gestation delay [16]. Many results demonstrate that models with time delay present some rich dynamic behaviors, such as periodic orbits, stability switching, chaos and some bifurcations [17,18]. In view of the complexity of the system, multiple delays often appear in prey–predator systems. The authors of [16] researched a predator–prey system with two time lags of fear and gestation of predator populations. Their results indicated that these two delays jointly brought some crucial influence on the population dynamics.
In addition, predator–prey population systems are exposed to a variety of environmental perturbations, such as drought, disease, flood, earthquake, etc. These environmental components are very important in the study of system dynamics [19,20]. The disturbances can have large influences on the density of populations directly or the parameter values of models indirectly. Therefore, it is necessary to make use of a stochastic system to precisely analyze and predict changes to prey–predator dynamics [21,22]. In stochastic models, not only the mean tendency but also the variance are taken into account. Usually, for fixed initial data, similar results are generated by the deterministic model, but stochastic prediction values are to be given by the stochastic models. In reality, researchers propose many stochastic models to formulate the effects of white noise on prey–predator dynamics [23,24,25].
To the best of our knowledge, an interaction between predator and prey species together with fear in prey, Beddington–DeAngelis functional response, two time delays (delayed fear and delayed gestation) and environmental white noise has never been studied. Hence, we are encouraged to study its kinetics.
For dynamic systems, the existence of a stable equilibrium status is very important for deterministic models, so the main subject of our research is to study the existence and dynamics of the equilibrium status of above predator–prey model and then analyze the potential effects of fear, delay and stochastic parameters theoretically and numerically.
Our study is structured as follows: A deterministic population model with fear of predators and delays is presented, and the existence and global asymptotic stability of its equilibrium status are analyzed in Section 2. Then, the deterministic model is extended to a stochastic environment, and the stationary distribution of its solutions is studied in Section 3. The theoretical results are indicated by numerical examples, and the concrete influences of fear, delays and stochastic environmental disturbances are explored in Section 4. Finally, in Section 5, a summary discussion is given to conclude our research.

2. Deterministic Scenario

2.1. Model Formulation

The popular predator–prey system with Beddington–DeAngelis-type density-dependent functional response is as follows:
d X d T = X A 1 B 1 X P Y C 1 + X + C 2 Y , d Y d T = Y A 2 B 2 Y + Q P X C 1 + X + C 2 Y .
Taking the effect of fear of predators into account, the growth rate of prey will be changed. Then, we consider the modified growth rate of prey in the form A 1 1 + f Y , which is a monotonic decreasing function about both fear f and species Y, where f is a fear parameter. Consequently, the above model becomes
d X d T = X A 1 1 + f Y B 1 X P Y C 1 + X + C 2 Y , d Y d T = Y A 2 B 2 Y + Q P X C 1 + X + C 2 Y ,
where X is the prey population and Y is the predator population at time T. All parameter values in this model are positive for biological justification, and their biological meanings see Table 1.
For the convenience of our analysis in system dynamics, using the transformation as X = A 1 x B 1 , Y = A 1 P y B 1 and A 1 T = t , we can then get the system listed below:
d x d t = x 1 x 1 + f y c y 1 + a 1 x + a 2 y , d y d t = y c x 1 + a 1 x + a 2 y m e y ,
where c = P Q C 1 B 1 , a 1 = A 1 C 1 B 1 , a 2 = C 2 A 1 Q C 1 B 1 , m = A 2 A 1 , e = B 2 Q B 1 .
Incorporating the time lag of fear of prey induced by the predator and the time lag of the predator’s gestation, we get the delayed version of System (1) as follows:
d x d t = x 1 x 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y , d y d t = y c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y .
The initial history functions are
x ( θ ) = ϕ 1 ( θ ) , y ( θ ) = ϕ 2 ( θ ) , τ θ 0 ,
where ϕ = ( ϕ 1 , ϕ 2 ) T C ( [ τ , 0 ] , R + 2 ) , τ = max { τ 1 , τ 2 } , and τ 1 and τ 2 denote the fear and gestation delay, respectively; C ( [ τ , 0 ] , R + 2 ) is a Banach space of continuous functions with norm | | ϕ | | = sup θ [ τ , 0 ] { | ϕ ( θ ) | } , | ϕ ( θ ) | = ϕ 1 2 ( θ ) + ϕ 2 2 ( θ ) . We assume ϕ ( θ ) R + 2 = { ( x , y ) | x > 0 , y > 0 } for any θ [ τ , 0 ] so as to satisfy the biological justification.

2.2. The Positivity and Permanence of Solutions

We begin with a comparison theorem and the definition of permanence.
Lemma 1
([26]). For equation
d x ( t ) d t = p x ( t τ ) q x ( t ) r x 2 ( t ) ,
where p , q , r , τ > 0 and x ( t ) > 0 for t [ τ , 0 ] , we then have
lim t x ( t ) = p q r , i f p > q , lim t x ( t ) = 0 , i f p < q .
Definition 1 (Permanence).
For System (2) with positive initial values ϕ i ( θ ) > 0 , θ [ τ , 0 ] ( i = 1 , 2 ) , we assume that there exist two positive constants N and M satisfying 0 < N M and
min lim inf t x ( t ) , lim inf t y ( t ) N , max lim sup t x ( t ) , lim sup t y ( t ) M ,
then, (2) is said to be permanent.
To assure System (2) is reasonable in terms of biological justification, we begin with the proof of positivity and boundedness of the solutions of (2). By integrating both sides of the former equality of (2), we have
x ( t ) = x ( 0 ) exp 0 t ( 1 x ) 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y d s > 0 .
Similarly, we deduce that
y ( t ) = y ( 0 ) exp 0 t c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y d s > 0 .
By the second equation of (2), then
d y d t c a 1 y ( t τ 2 ) m y e y 2 .
By Lemma 1, then
y c a 1 d e T if c a 1 d > 0 .
On the other hand, by use of the first equation of (2), we have
d x ( t ) d t x ( t ) 1 x ( t ) .
For the above inequality, by the same integrating, we have
x ( t ) 1 .
That is to say, all solutions of Model (2) are bounded.
Further, from the former equality of Model (2), we get
d x d t x 1 x 1 + f T c a 2 .
Then
x ( t ) ( 1 1 + f T c a 2 ) ( 1 + f T ) N if 1 1 + f T > c a 2 .
Using (2) again, then
d y d t y N 1 + a 1 + a 2 T m e y .
Thus
y N 1 + a 1 + a 2 T m e S if N 1 + a 1 + a 2 T m > 0 .
We summarize the above analysis that System (2) is permanent.

2.3. Existence of Equilibrium Status

We plan to investigate the existence of a coexisting equilibrium status of (2) in this subsection. Clearly, the equilibrium status of (2) satisfies the following equalities:
1 x 1 + f y c y 1 + a 1 x + a 2 y = 0 , c x 1 + a 1 x + a 2 y = m + e y .
By the second equality of (3), the equilibrium point x * should meet
x * = ( m + e y * ) ( 1 + a 2 y * ) c a 1 ( m + e y * ) ,
which is positive under the condition c a 1 ( m + e y * ) > 0 . Putting x = x * into the first equality, then y * satisfies the following equality:
e f y 3 + ( m f + e ) y 2 + m y ( d + e y ) ( 1 + a 2 y ) c a 1 ( m + e y ) ( 1 ( d + e y ) ( 1 + a 2 y ) c a 1 ( m + e y ) ) = 0 .
Eliminate the denominators; then, we have
a 1 2 e 3 f y 5 + D 1 y 4 + D 2 y 3 + D 3 y 2 + D 4 y + D 5 = 0 ,
where
D 1 = a 1 2 e 2 ( e + m f ) 2 a 1 e ( c a 1 m ) e f + e 2 a 2 2 ) , D 2 = a 1 2 e 2 m + ( c a 1 m ) 2 e f 2 a 1 e ( c a 1 m ) ( e + m f ) + a 2 e ( 2 e + a 1 e + 2 a 2 m ) , D 3 = ( c a 1 m ) 2 ( e + m f ) 2 a 1 e ( c a 1 m ) m + e 2 + a 1 e 2 + 2 m e a 1 a 2 + 4 m e a 2 + m 2 2 a 2 2 c e a 2 , D 4 = ( c a 1 m ) 2 m + 2 m e + 2 m e a 1 + m 2 a 1 a 2 + 2 m 2 a 2 c m a 2 c e , D 5 = m ( m ( c a 1 m ) ) .
By verification, if max a 2 , m , 2 a 1 e m e + m f < c a 1 m < min a 1 m , a 1 m ( e + f ) 2 e , a 1 ( e + m f ) 2 f , then D i > 0 ( i = 1 , 2 , 3 , 4 ) and D 5 < 0 . By Descartes’ rule of signs, we conclude that there exists a unique positive solution y * if the above condition holds. Hence, System (2) has a unique positive equilibrium E * ( x * , y * ) .

2.4. Global Asymptotic Stability

With respect to the long-time dynamic behavior of E * ( x * , y * ) , we have the below conclusion:
Theorem 1.
For System (2) with initial data ϕ i ( θ ) > 0 , θ [ τ , 0 ] ( i = 1 , 2 ) , we assume that
1 1 + f T > c a 1 y * 1 + a 1 x * + a 2 y * + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) , e > c ( 1 + a 2 y * ) 1 + a 1 x * + a 2 y * .
Then, the equilibrium status E * ( x * , y * ) of Model (2) is globally asymptotically stable, i.e., for any solution ( x ( t ) , y ( t ) ) of (2), we have lim t x ( t ) = x * , lim t y ( t ) = y * .
Proof. 
Define a functional listed below:
V 1 ( x ) = x x * x * ln x / x * .
Due to the monotonicity of V 1 ( x ) with respect to x, obviously V 1 ( x ) is positive on t [ 0 , ) . Differentiate V 1 ( x ) on variable t along System (2); then
d V 1 ( x ) d t = ( x x * ) 1 x 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y = ( x x * ) 1 x 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y 1 x * 1 + f y * + c y * 1 + a 1 x * + a 2 y * = ( x x * ) ( ( 1 + f y * ) ( x x * ) + f ( 1 x * ) ( y ( t τ 1 ) y * ) ( 1 + f y ( t τ 1 ) ) ( 1 + f y * ) c ( 1 + a 1 x * ) ( y y * ) c a 1 y * ( x x * ) ( 1 + a 1 x + a 2 y ) ( 1 + a 1 x * + a 2 y * ) ) 1 1 + f T c a 1 y * 1 + a 1 x * + a 2 y * ( x x * ) 2 .
Similarly, define
V 2 ( y ) = y y * y * ln y / y * .
Function V 2 ( y ) : R R + is positive. Differentiating V 2 ( y ) about variable t along (2), we get
d V 2 ( y ) d t = ( y y * ) c x 1 + a 1 x + a 2 y d e y = ( y y * ) c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) d e y c x * 1 + a 1 x * + a 2 y * + d + e y * = e ( y y * ) 2 + ( y y * ) c ( 1 + a 2 y * ) ( x ( t τ 2 ) x * ) c a 2 x * ( y ( t τ 2 ) y * ) ( 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) ) ( 1 + a 1 x * + a 2 y * ) e ( y y * ) 2 + c ( 1 + a 2 y * ) 1 + a 1 x * + a 2 y * ( x ( t τ 2 ) x * ) ( y y * ) e ( y y * ) 2 + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( ( x ( t τ 2 ) x * ) 2 + ( y y * ) 2 ) e c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( y y * ) 2 + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( ( x ( t τ 2 ) x * ) 2 .
Introduce V 3 ( t ) below to eliminate the terms of delay that appeared in the above inequality:
V 3 ( t ) = t τ 2 t c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( ( x ( s τ 2 ) x * ) 2 d s .
Then
d V 2 ( y ) d t + d V 3 ( t ) d t e c ( 1 + a 2 y * ) 1 + a 1 x * + a 2 y * ( y y * ) 2 + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( ( x x * ) 2 .
Let V ( t ) = V 1 ( x ) + V 2 ( y ) + V 3 ( t ) ; together with the previous condition, we then have
d V ( t ) d t = 1 1 + f T c a 1 y * 1 + a 1 x * + a 2 y * c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( x x * ) 2 e c ( 1 + a 2 y * ) 1 + a 1 x * + a 2 y * ( y y * ) 2 < 0 .
By functional differential equation theory [18], we conclude that the equilibrium status E * ( x * , y * ) of (2) is globally asymptotically stable. □
Remark 1.
From the biological point of view, if the population growth rate of prey affected by fear is large, such that 1 1 + f T > c a 1 y * 1 + a 1 x * + a 2 y * + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) , then the prey is persistent. On the other hand, compared to the conversion rate of the prey, if the interspecific competition rate of the predator is not too low, such that e > c ( 1 + a 2 y * ) 1 + a 1 x * + a 2 y * , then the predator population is inhibited from increasing excessively and remains persistent.
When both populations hold, then the prey and predator can coexist near the equilibrium state and the system has global stability, which is biologically reasonable.
Theorem 2.
For System (1) with initial data x ( 0 ) > 0 , y ( 0 ) > 0 , E * ( x * , y * ) is globally asymptotically stable provided that 1 1 + f T > c a 1 y * 1 + a 1 x * + a 2 y * .
Proof. 
We apply the same techniques presented in Theorem 1 to prove the conclusion. Take
V 1 ( x ) = x x * x * ln x / x * , V 2 ( y ) = y y * y * ln y / y * ,
which are the same as before. Compute the derivatives of V 1 ( x ) and V 2 ( y ) with respect to t along System (1); then
d V 1 ( x ) d t = ( x x * ) 1 x 1 + f y c y 1 + a 1 x + a 2 y 1 x * 1 + f y * + c y * 1 + a 1 x * + a 2 y * 1 1 + f T c a 1 y * 1 + a 1 x * + a 2 y * ( x x * ) 2 c ( 1 + a 1 x * ) ( 1 + a 1 x + a 2 y ) ( 1 + a 1 x * + a 2 y * ) ( x x * ) ( y y * ) ,
and
d V 2 ( y ) d t = ( y y * ) c x 1 + a 1 x + a 2 y d e y c x * 1 + a 1 x * + a 2 y * + d + e y * e ( y y * ) 2 + c ( 1 + a 1 y * ) ( 1 + a 1 x + a 2 y ) ( 1 + a 1 x * + a 2 y * ) ( x x * ) ( y y * ) .
Let V ( t ) = V 1 ( x ) + V 2 ( y ) . Together with the condition 1 1 + f T > c a 1 y * 1 + a 1 x * + a 2 y * , we can derive that
d V ( t ) d t = ( 1 1 + f T c a 1 y * 1 + a 1 x * + a 2 y * ) ( x x * ) 2 e ( y y * ) 2 < 0 .
Hence, the result holds. □
Remark 2.
By Theorem 1, it seems that delays have no effect on stability because they are absent from the restricted conditions; while applying the same techniques, Theorem 2 reveals that the assumption guaranteeing the global stability of the equilibrium status E * is lower than that in Theorem 1 with no delays ( τ 1 = τ 2 = 0 ). This implies that, in fact, the delays have some important influence on global system stability, as presented by simulations in Section 4.

3. Stochastic Scenario

3.1. Model Formulation

As stated before, predator–prey population systems are often influenced by environmental disturbances. In this section, through consideration of environmental white noise, we extend model (2) to the following stochastic model.
Theoretically, environmental fluctuation may disturb all parameters in Model (2), but the demographic randomness mainly influences the birth velocity of the prey population and the mortality velocity of the predator population [27]. Mathematical researchers often introduce white noise to model environmental perturbations. By the central limit theorem, the error terms of the influenced parameters follow a normal distribution. References [28,29] showed that the stochastic disturbance was related to the difference between existing populations and the equilibrium status. Therefore, perturbing the growth velocity of the prey population and the natural mortality velocity of the predator population, that is, adding a random disturbance coefficient to the original constants:
A 1 A 1 + σ 1 ( X X * ) , A 2 A 2 + σ 2 ( Y Y * ) .
Using the transformation in Section 2.1 together with the method in [9], we then obtain the following stochastic model:
d x d t = x 1 x 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y + σ 1 x ( x x * ) d w 1 ( t ) , d y d t = y c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y σ 2 y ( y y * ) d w 2 ( t ) ,
where w 1 ( t ) and w 2 ( t ) are two Brownian motions that are independent and defined on a complete probability space ( Ω , F , P ) , and σ i 2 ( i = 1 , 2 ) denotes the intensity of white noise.

3.2. Existence and Stochastic Permanence of Model (4)

We investigate the existence of stochastic solutions and stochastic permanence of Model (4) in this subsection. Suppose x ( t ) is a homogeneous Markovian process in R n described by the stochastic equation:
d x ( t ) = f ( t , x ( t ) ) d t + k = 1 n g k ( t , x ( t ) ) d w k ( t ) .
The diffusion matrix is A ( x ) = ( a i j ( x ) ) , a i j ( x ) = k = 1 n g k i ( x ) g k j ( x ) . For Lyapunov functional V ( x ( t ) ) , we define the Lyapunov operator L of V ( x ( t ) ) as follows:
L V ( x ( t ) ) = V x ( x ( t ) ) f ( t , x ( t ) ) + 1 2 t r a c e [ g T ( t , x ( t ) ) V x x ( X ( t ) ) g ( t , x ( t ) ) ] .
Before study of the existence and positiveness of solutions, we give a lemma.
Lemma 2.
Let x ( t ) be a stochastic process satisfying
d x ( t ) = ( b ( t ) a ( t ) x ( t ) ) d t + σ ( t ) x ( t ) d w ( t ) .
For System (6), there exists a unique positive solution x(t) with any initial value x ( 0 ) = x 0 > 0 that is global and is represented by
x ( t ) = e 0 t ( a ( u ) + σ 2 ( u ) 2 ) d u + σ ( u ) d w ( u ) x 0 + 0 t b ( s ) e 0 t ( a ( u ) + σ 2 ( u ) 2 ) d u σ ( u ) d w ( u ) d s ,
or
x ( t ) = x 0 e 0 t a ( s ) d s + 0 t e a ( s ) ( t s ) ( b ( s ) d s + σ ( s ) x ( s ) d w ( s ) ) .
Proof. 
First, we verify
x ( t ) = e 0 t ( a ( u ) + σ 2 ( u ) 2 ) d u + σ ( u ) d w ( u ) x 0 + 0 t b ( s ) e 0 t ( a ( u ) + σ 2 ( u ) 2 ) d u σ ( u ) d w ( u ) d s ,
is the solution of (6). Let ξ ( t ) = 0 t ( a ( u ) + σ 2 ( u ) 2 ) d u + σ ( u ) d w ( u ) . Then ξ ( t ) is a stochastic process, and
d ξ ( t ) = a ( t ) + σ 2 ( t ) 2 d t + σ ( t ) d w ( t ) .
Define V ( ξ ) = e ξ . Applying Itô’s formula to V ( ξ ) results in
d V ( ξ ) = e ξ a ( t ) σ 2 ( t ) 2 + 1 2 σ 2 ( t ) d t + e ξ σ ( t ) d w ( t ) .
Consequently, we have
d x ( t ) = x 0 d ( e ξ ) + d e ξ 0 t b ( s ) e ξ d s = x 0 d ( e ξ ) + d ( e ξ ) × 0 t b ( s ) e ξ d s + e ξ × d 0 t b ( s ) e ξ d s + d ( e ξ ) × d 0 t b ( s ) e ξ d s = x 0 + 0 t b ( s ) e ξ d s e ξ a ( t ) σ 2 ( t ) 2 + 1 2 σ 2 ( t ) d t + σ ( t ) d w ( t ) + e ξ b ( t ) e ξ d t + e ξ a ( t ) σ 2 ( t ) 2 + 1 2 σ 2 ( t ) d t + e ξ σ ( t ) d w ( t ) b ( t ) e ξ d t = a ( t ) x ( t ) d t + σ ( t ) x ( t ) d w ( t ) + b ( t ) d t = ( b ( t ) a ( t ) x ( t ) ) d t + σ ( t ) x ( t ) d w ( t ) .
Next, we illustrate that
x ( t ) = x 0 e 0 t a ( s ) d s + 0 t e a ( s ) ( t s ) ( b ( s ) d s + σ ( s ) x ( s ) d w ( s ) )
is the solution of (6). Let
ρ ( t ) = 0 t e a ( s ) ( t s ) ( b ( s ) d s + σ ( s ) x ( s ) d w ( s ) ) ;
then ρ ( t ) is a stochastic process, and
d ρ = b ( t ) d t + σ ( t ) x ( t ) d w ( t ) a ( t ) 0 t e a ( s ) ( t s ) ( b ( s ) d s + σ ( s ) x ( s ) d w ( s ) ) .
An easy computation gives
d x ( t ) = a ( t ) x 0 e 0 t a ( s ) d s + b ( t ) d t + σ ( t ) x ( t ) d w ( t ) a ( t ) 0 t e a ( s ) ( t s ) ( b ( s ) d s + σ ( s ) x ( s ) d w ( s ) ) = a ( t ) x 0 e 0 t a ( s ) d s + 0 t e a ( s ) ( t s ) ( b ( s ) d s + σ ( s ) x ( s ) d w ( s ) ) + b ( t ) d t + σ ( t ) x ( t ) d w ( t ) = ( b ( t ) a ( t ) x ( t ) ) d t + σ ( t ) x ( t ) d w ( t ) .
Hence, the result holds. □
If the coefficients of (6) are all constants, we can obtain the following corollary.
Corollary 1.
Suppose x ( t ) is a stochastic process satisfying
d x ( t ) = ( b a x ( t ) ) d t + σ x ( t ) d w ( t ) .
For System (7), there exists a unique positive solution x ( t ) with any initial value x ( 0 ) = x 0 > 0 that is global and is represented by
x ( t ) = e ( a + σ 2 2 ) t + σ w ( t ) x 0 + b 0 t e ( a + σ 2 2 ) s σ w ( s ) d s .
If b = 0 , then
x ( t ) = x 0 e ( a + σ 2 2 ) t + σ w ( t ) .
Remark 3.
In [30], there are some errors in the use of Itô’s formula to solve the stochastic equation. For example, by applying Itô’s formula, the authors stated that the solution of the following equation (page 1511, line-4):
d H 1 ( t ) = r g ζ d 1 + β ( 1 ζ ) β H 1 + ζ 2 σ 1 2 H 1 d t ζ σ 1 H 1 d w 1
is
H 1 ( t ) = r g ζ d 1 ζ 2 σ 1 2 + β ( 1 ζ ) β + H 1 ( 0 ) r g ζ d 1 ζ 2 σ 1 2 + β ( 1 ζ ) β e ( ζ 2 σ 1 2 g ζ + d 1 β ( 1 ζ ) β ) t ζ σ 1 w 1 ( t ) .
Easy verification implies that this is wrong. Actually, one can derive from Corollary 1 that the solution to the above equation is
x ( t ) = e ( g ζ d 1 ζ 2 σ 1 2 + β ( 1 ζ ) β ) ζ 2 σ 1 2 ( t ) 2 ζ σ 1 w ( t ) H 1 ( 0 ) + r 0 t e g ζ d 1 ζ 2 σ 1 2 + β ( 1 ζ ) β + ζ 2 σ 1 2 2 ζ σ 1 w ( s ) d s .
Now, by use of Corollary 1, we discuss the existence and uniqueness of a global positive solution of system (4).
Theorem 3.
For System (4) with initial data ϕ i ( θ ) > 0 , θ [ τ , 0 ] ( i = 1 , 2 ) , there exists a positive solution ( x ( t ) , y ( t ) ) a.s. that is positive and global.
Proof. 
It is not difficult to find that the parameters of Model (4) are local Lipschitz. By use of the theory of existence of stochastic processes, we assert that Model (4) has a unique positive solution on [ τ , τ e ] , where τ e denotes the explosion time of solutions. Thus, we only want to prove τ e = . For this purpose, we apply the comparison method of stochastic systems.
By the positivity of solutions of (4), we get from (4) that
d x ( t ) x ( t ) d t + σ 1 x ( t ) ( x ( t ) x * ) d w 1 ( t ) .
Let X 1 ( t ) be the solution of
d X 1 ( t ) = X 1 ( t ) d t + σ 1 M X 1 ( t ) d w 1 ( t ) , X 1 ( 0 ) = X 1 0 ,
where M is the super boundness of x ( t ) defined in Theorem 4. Then by use of Corollary 1, the solution of the previous equation is
X 1 ( t ) = X 1 0 e ( 1 σ 1 2 M 2 2 ) t + M w 1 ( t ) .
By comparison theory of stochastic differential equations, we have
x ( t ) X 1 ( t ) .
From the second equation of (4), then
d y ( t ) c a 2 m y ( t ) d t σ 2 y ( t ) ( y ( t ) y * ) d w 2 ( t ) .
Let Y 1 ( t ) be the solution of
d Y 1 ( t ) = c a 2 m Y 1 ( t ) + σ 2 y * Y 1 ( t ) d w 2 ( t ) , Y 1 ( 0 ) = Y 1 0 .
Again, using Corollary 1 gives
Y 1 ( t ) = Y 1 0 e ( c a 2 m σ 2 2 y * 2 2 ) t + σ 2 y * w 1 ( t ) .
Then
y ( t ) Y 1 ( t ) .
On the other hand, from the second equation of (4), we have
d y ( t ) m y ( t ) d t σ 2 y ( t ) ( y ( t ) y * ) d w 2 ( t ) .
Let Y 2 ( t ) be the solution of
d Y 2 ( t ) = m Y 2 ( t ) d t σ 2 M Y 2 ( t ) d w 2 ( t ) , Y 2 ( 0 ) = Y 2 0 ,
where M is the super boundness of y ( t ) defined in Theorem 4. By use of Corollary 1 again, we have the solution of the previous equation:
Y 2 ( t ) = Y 2 0 e ( m + σ 2 2 M 2 2 ) t σ 2 M w 2 ( t ) .
Therefore,
y ( t ) Y 2 ( t ) .
Similarly, we have
d x ( t ) ( 1 + c a Y 1 ( t ) ) x ( t ) d t σ 1 x * x ( t ) d w 1 ( t ) .
Let X 2 ( t ) be the solution of
d X 2 ( t ) = ( 1 + c a Y 1 ( t ) ) X 2 ( t ) d t σ 1 x * X 2 ( t ) d w 1 ( t ) , X 2 ( 0 ) = X 2 0 .
By verification, the formal solution is
X 2 ( t ) = X 2 0 e 0 t ( 1 + c a Y 1 ( s ) ) d s + σ 1 2 x * 2 2 t σ 1 x * w 1 ( t ) .
By Lemma 2, we have
x ( t ) X 2 ( t ) .
Summarizing the above conclusions, we obtain
X 2 ( t ) x ( t ) X 1 ( t ) , Y 2 ( t ) y ( t ) Y 1 ( t ) , for t [ 0 , τ e ) .
Clearly, the solutions X 1 ( t ) , X 2 ( t ) , Y 1 ( t ) , Y 2 ( t ) exist for all t > 0 , so we conclude that the explosion time τ e = and the solution is global. □
Theorem 4.
Let X ( t ) = ( x ( t ) , y ( t ) ) T be a solution of (4) with initial value ϕ i ( θ ) > 0 , θ [ τ , 0 ] ( i = 1 , 2 ) . Then X ( t ) is stochastically permanent; that is, there are two constants M > 0 and N > 0 such that the solution X ( t ) of System (4) satisfies
P { | X ( t ) | M } 1 ε , P { | X ( t ) | N } 1 ε ,
for any ε ( 0 , 1 ) , where P is the probability of an event.
Proof. 
Define V 1 ( x ) = e t x , V 2 ( y ) = e t y and V ( x , y ) = V 1 ( x ) + V 2 ( y ) . Applying Itô’s formula to V ( x , y ) along System (4), then
d V ( x , y ) = e t x + e t x ( 1 x 1 + f y c y 1 + a 1 x + a 2 y ) d t + e t σ 1 x ( x x * ) d w 1 ( t ) + e t y + e t y c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y d t + e t σ 2 y ( y y * d w 2 ( t ) e t x + e t x ( 1 x ) d t + e t σ 1 x ( x x * ) d w 1 ( t ) + ( c / a 1 m ) 2 4 e e t d t + e t σ 2 y ( y y * ) d w 2 ( t ) ( e t x + e t ( 2 x x 2 x ) d t + e t σ 1 x ( x x * ) d w 1 ( t ) + ( c / a 1 m ) 2 4 e e t d t + e t σ 2 y ( y y * ) d w 2 ( t ) ( e t x + e t ( 1 x ) d t + e t σ 1 x ( x x * ) d w 1 ( t ) + ( c / a 1 m ) 2 4 e e t d t + e t σ 2 y ( y y * ) d w 2 ( t ) = e t 1 + ( c / a 1 m ) 2 4 e + e t σ 1 x ( x x * ) d w 1 ( t ) + e t σ 2 y ( y y * ) d w 2 ( t ) .
Taking the mathematical expectation, then
E 0 t d V ( x , y ) 1 + ( c / a 1 m ) 2 4 e ( e t 1 ) .
That is,
E V ( x , y ) V ( 0 , 0 ) + 1 + ( c / a 1 m ) 2 4 e ( e t 1 ) .
Therefore, we have
E ( x + y ) 1 + ( c / a 1 m ) 2 4 e + V ( 0 , 0 ) 1 + ( c / a 1 m ) 2 4 e e t .
Note that
| X ( t ) | = ( x 2 + y 2 ) 1 2 x + y .
Consequently,
E | X ( t ) | 1 + ( c / a 1 m ) 2 4 e + V ( 0 , 0 ) 1 + ( c / a 1 m ) 2 4 e e t .
Let M be a large enough positive number such that 1 + ( c / a 1 m ) 2 4 e M = ε , which is a sufficiently small positive constant. Applying Chebyshev’s inequality, we have
P { | X ( t ) | M } E | X ( t ) | M 1 + ( c / a 1 m ) 2 4 e M = ε .
Taking the superior limit of both sides of the above inequality leads to
lim sup t P { | X ( t ) | M } E | X ( t ) | M ε ,
which implies
P { | X ( t ) | M } 1 ε .
On the other hand,
E V 1 ( x ) = x ( 0 ) + 0 t e s E x 1 + 1 x 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y d t x ( 0 ) + 0 t e s E x 1 + 1 x 1 + f y ( t τ 1 ) c a 2 d t .
By the first conclusion, we have
1 + 1 x 1 + f y ( t τ 1 ) 1 + 1 x 1 + f M .
Let f ( x ) = x 1 + 1 x 1 + f M c a 2 . To find the supremum value of f ( x ) , we compute
sup x 0 f ( x ) = a 2 c ( 1 + f M ) 4 a 2 2 ( 1 + f M ) .
By the definition of supremum, then
E V ( x ) x ( 0 ) + a 2 c ( 1 + f M ) 4 a 2 2 ( 1 + f M ) ( e t 1 ) .
Consequently, as t , we have
E ( x ) a 2 c ( 1 + f M ) 4 a 2 2 ( 1 + f M ) + x 0 a 2 c ( 1 + f M ) 4 a 2 2 ( 1 + f M ) e t a 2 c ( 1 + f M ) 4 a 2 2 ( 1 + f M ) ι .
Let ς be a large enough positive number such that 1 ι ς = ε is a sufficiently small positive constant. Applying Chebyshev’s inequality again, we have
P 1 x ς 1 ς E ( x ) 1 ς ι = ε ,
i.e.,
P x 1 ς ε ,
which implies
P x 1 ς 1 ε .
Applying Itô’s formula to V 2 ( y ) = e t y , then
d V 2 ( y ) = e t y + e t y c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y d t + e t σ 2 y ( y y * d w 2 ( t ) e t y 1 + c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y d t + e t σ 2 y ( y y * ) d w 2 ( t ) .
Thus,
1 + c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) m e y 1 + c / ς 1 + a 1 M + a 2 M m e y .
Let f ( y ) = y 1 + c / ς 1 + a 1 M + a 2 M m e y . By the same reasoning, we deduce that the supremum value of f ( y ) is as follows:
sup x 0 f ( y ) = 1 + c / ς 1 + a 1 M + a 2 M m 2 4 e ϖ .
Consequently,
E ( y ) ϖ + x ( 0 ) ϖ e t ϖ , t .
Take ς as before such that 1 ϖ ς = ε > 0 is small enough. Applying Chebyshev’s inequality, we have
P y 1 ς ε .
It is clear that | X ( t ) | = ( x 2 + y 2 ) 1 2 2 ς N ; then
P { X ( t ) N } ε .
This completes the proof. □
Remark 4.
By the second equation of (4), we have d y y ( c a 1 m e y ) + σ 2 y ( y y * ) d w 2 ( t ) . The numerical solution is
y = exp [ ( c / a 1 m + y * 2 σ 2 2 / 2 ) t σ 2 y * w 2 ( t ) ] y 0 + e 0 t exp [ ( c / a 1 m + y * 2 σ 2 2 / 2 ) s σ 2 y * w 2 ( s ) ] d s .
Thus, approximately, y c / a 1 m + y * 2 σ 2 2 / 2 e = T ˜ if y 0 e c / a 1 m + y * 2 σ 2 2 / 2 , which is used in the proof of the stationary distribution.

3.3. Stability in Distribution

For the proof of existence and uniqueness of a stable distribution ψ of (4) that is ergodic, we present a useful lemma.
Lemma 3
([25]). Suppose O R n is a region with a regular boundary. For System (5), if the below conditions hold:
(1) 
There exists a constant M > 0 such that i , j = 1 n a i , j ( x ) ξ i ξ j M | ξ | 2 , x O , ξ R n .
(2) 
There exists a non-negative function V ( x ) C 2 such that L V ( x ) is negative for any x R n / O .
Then the Markovian process X ( t ) is stable in distribution, where the distribution ψ ( · ) is unique and ergodic; Hence,
P lim t 1 t 0 t f ( X ( s ) ) d s = R n f ( x ) ψ ( d x ) = 1 ,
for any x R n , where function f ( · ) is integrable about the given measure ψ.
Theorem 5.
Suppose E * ( x * , y * ) is the unique positive equilibrium status of System (2) such that
1 1 + f T ˜ > c a 1 y * 1 + a 1 x * + a 2 y * + f ( 1 x * ) 2 ( 1 + f y * ) + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) + σ 1 2 x * 2 , e > c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) + f ( 1 x * ) 2 ( 1 + f y * ) + σ 2 2 y * 2 ,
Then for system (4) with initial value ϕ i ( θ ) > 0 , θ [ τ , 0 ] ( i = 1 , 2 ) , there exists a unique stationary distribution ψ , and it has an ergodic property.
Proof. 
By Theorem 3, we know that the region R + 2 is invariant for process ( x ( t ) , y ( t ) ) . By Lemma 3, we only need to find a bounded domain O R + 2 satisfying those conditions. We recall that ( x * , y * ) is the positive equilibrium status of (1) satisfying
1 x * 1 + f y * = c y * 1 + a 1 x * + a 2 y * , c x * 1 + a 1 x * + a 2 y * = m + e y * .
Next, we go to construct the required functional. Define V 1 ( x ) : D R + R + as follows:
V 1 ( x ) = 0 x x * v v + x * d v .
Obviously, it is positive for all x x * . Differentiating V 1 ( x ) along the solution of (4), we get
L V 1 ( x ) ( 4 ) = ( x x * ) 1 x 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y + σ 1 2 x * 2 ( x x * ) 2 = ( x x * ) ( 1 x ) 1 + f y ( t τ 1 ) c y 1 + a 1 x + a 2 y 1 x * 1 + f y * + c y * 1 + a 1 x * + a 2 y * + σ 1 2 x * 2 ( x x * ) 2 1 1 + f T ˜ c a 1 y * 1 + a 1 x * + a 2 y * σ 1 2 x * 2 ( x x * ) 2 + f ( 1 x * ) 1 + f y * ( y ( t τ 1 ) y * ) ( x x * ) 1 1 + f T ˜ c a 1 y * 1 + a 1 x * + a 2 y * σ 1 2 x * 2 f ( 1 x * ) 2 ( 1 + f y * ) ( x x * ) 2 + f ( 1 x * ) 2 ( 1 + f y * ) ( y ( t τ 1 ) y * ) 2 .
Define
V 2 ( t ) = f ( 1 x * ) 2 ( 1 + f y * ) t τ 1 t ( y ( s ) y * ) 2 d s .
Differentiating V 2 ( t ) on variable t gives
d V 2 ( t ) d t = f ( 1 x * ) 2 ( 1 + f y * ) ( y y * ) 2 f ( 1 x * ) 2 ( 1 + f y * ) ( y ( t τ 1 ) y * ) 2 .
Then
L ( V 1 ( x ) + V 2 ( t ) ) 1 1 + f T ˜ c a 1 y * 1 + a 1 x * + a 2 y * σ 1 2 x * 2 f ( 1 x * ) 2 ( 1 + f y * ) ( x x * ) 2 + f ( 1 x * ) 2 ( 1 + f y * ) ( y y * ) 2 .
Define
V 3 ( y ) = 0 y y * v v + y * d v .
Similarly, it is defined positively in t [ 0 , ) . Differentiating V 3 ( y ) along the solution of (4) again, we have
L V 3 ( y ) ( 4 ) = ( y y * ) c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) d e y + σ 2 2 y * 2 ( y y * ) 2 = ( y y * ) c x ( t τ 2 ) 1 + a 1 x ( t τ 2 ) + a 2 y ( t τ 2 ) d e y c x * 1 + a 1 x * + a 2 y * + d + e y * + σ 2 2 y * 2 ( y y * ) 2 e σ 2 2 y * 2 ( y y * ) 2 + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( ( x ( t τ 2 ) x * ) 2 + ( y y * ) 2 ) e σ 2 2 y * 2 c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( y y * ) 2 + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( x ( t τ 2 ) x * ) 2 .
In order to eliminate the delay term, we define
V 4 ( t ) = c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) t τ 2 t ( x ( s ) x * ) 2 d s .
Thus,
L ( V 3 ( y ) + V 4 ( t ) ) e σ 2 2 y * 2 c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( y y * ) 2 + c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) ( x x * ) 2 .
Let V ( t ) = V 1 ( x ) + V 2 ( t ) + V 3 ( y ) + V 4 ( t ) , then
L V ( t ) = 1 1 + f T ˜ c a 1 y * 1 + a 1 x * + a 2 y * f ( 1 x * ) 2 ( 1 + f y * ) c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) σ 1 2 x * 2 ( x x * ) 2 e c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) f ( 1 x * ) 2 ( 1 + f y * ) σ 2 2 y * 2 ( y y * ) 2 .
Under the given condition, obviously L V ( t ) < 0 , which implies that Condition (2) of Lemma 3 holds. On the other hand, it is clear that there exists a M > 0 such that
i , j = 1 2 a i j ( x , y ) ξ i ξ j = σ 1 2 x 2 ( x x * ) 2 ξ 1 2 + σ 2 2 y 2 ( y y * ) 2 ξ 2 2 M | ξ | 2 ,
for all ( x , y ) O , ξ R 2 . Hence, Condition (1) of Lemma 3 is also satisfied. Together with the positive invariant property, System (4) therefore has a stable stationary distribution ψ ( · ) , and it is ergodic. □
If τ 1 = τ 2 = 0 , then System (4) leads to (8) as follows.
d x d t = x 1 x 1 + f y c y 1 + a 1 x + a 2 y + σ 1 x ( x x * ) d w 1 ( t ) , d y d t = y c x 1 + a 1 x + a 2 y m e y + σ 2 y ( y y * ) d w 2 ( t ) .
For (8), by the same manner, we have the following finding.
Corollary 2.
For (8) with initial value ( x 0 , y 0 ) R + 2 , there is a stationary distribution ψ with an ergodic property provided that 1 1 + f T ˜ > c a 1 y * 1 + a 1 x * + a 2 y * + σ 1 2 x * 2 and e > σ 2 2 x * 2 .
If f = 0 , τ 1 = τ 2 = 0 , by use of the ergodic property of ψ ( · ) , we have the following proposition.
Proposition 1.
Suppose the conditions in Theorem 5 hold. Then
lim t 1 t 0 t x ( s ) d s = R + 2 z 1 ψ ( d z 1 , d z 2 ) , a . e . , lim t 1 t 0 t y ( s ) d s = R + 2 z 2 ψ ( d z 1 , d z 2 ) , a . e .
Proof. 
By Theorem 4, we know that E ( x ) M , E ( y ) M for all t [ 0 , ) , where M = 1 + ( c / a 1 m ) 2 4 e . By the ergodic property, for all n > 0 we get
lim t 1 t 0 t ( x ( s ) n ) d s = R + 2 ( z 1 n ) ψ ( d z 1 , d z 2 ) , a . e . ,
lim t 1 t 0 t ( y ( s ) n ) d s = R + 2 ( z 2 n ) ψ ( d z 1 , d z 2 ) , a . e .
By the dominated convergence theorem, we have
E [ lim t 1 t 0 t ( x ( s ) n ) d s ] = lim t 1 t 0 t E ( x ( s ) n ) d s ] M ,
E [ lim t 1 t 0 t ( y ( s ) n ) d s ] = lim t 1 t 0 t E ( y ( s ) n ) d s ] M .
Thus,
R + 2 ( z 1 n ) ψ ( d z 1 , d z 2 ) M , R + 2 ( z 2 n ) ψ ( d z 1 , d z 2 ) M .
Taking the limit as n leads to
R + 2 z 1 ψ ( d z 1 , d z 2 ) M , R + 2 z 2 ψ ( d z 1 , d z 2 ) M .
That is to say, z 1 and z 2 are integrable with respect to the measure ψ . From Thereom 5, we know that the measure ψ has an ergodic property; then
lim t 1 t 0 t x ( s ) d s = R + 2 z 1 ψ ( d z 1 , d z 2 ) , a . e . ,
lim t 1 t 0 t y ( s ) d s = R + 2 z 2 ψ ( d z 1 , d z 2 ) , a . e .
This completes the proof. □
Remark 5.
By Proposition 1, we find that the mean and variance of the stationary distribution ψ ( · ) can be easily computed. It is clear that the mean
x ¯ = lim t 1 t 0 t x ( s ) d s = R + 2 z 1 ψ ( d z 1 , d z 2 ) , a . e . ,
y ¯ = lim t 1 t 0 t y ( s ) d s = R + 2 z 2 ψ ( d z 1 , d z 2 ) , a . e .
and the covariance matrix of ψ ( · ) is
Π = a 11 a 12 a 21 a 22 ,
where
a 11 = lim t 1 t 0 t ( x ( s ) x ¯ ) 2 d s , a 12 = a 21 = lim t 1 t 0 t ( x ( s ) x ¯ ) ( y ( s ) y ¯ ) d s , a 22 = lim t 1 t 0 t ( y ( s ) y ¯ ) 2 d s .
The mean and variance of ψ ( · ) are useful for population system.

4. Numerical Analysis

Some numerical examples are given to visualize the globally asymptotically stable equilibrium status, and the important impact of various parameters on system dynamics are to be explored in this section.
First, we validate our theoretical results by numerical examples. Fix a set of parameter values given below:
c = 1.5 , a 1 = 3 , a 2 = 0.3 , m = 0.3 , e = 0.8 .
We vary the fear value to see how it affects the equilibrium status. Take f = 0 , i.e., there is no fear effect on System (2); then, the equilibrium status is E 0 * ( 0.9669 , 0.0867 ) . For f = 5 and f = 50 , the corresponding equilibrium statuses are E * ( 0.9532 , 0.0850 ) and E ^ * ( 0.8597 , 0.0726 ) , respectively. That is, if the fear f induced by the predator is large, then both the predator and prey populations are smaller at equilibrium. That means the value of fear from the predator brings some negative influence to the system dynamics (see Figure 1).
Then, we investigate the global asymptotic stability of equilibrium statuses of E 0 * ( 0.9669 , 0.0867 ) and E * ( 0.9532 , 0.0850 ) . For the fixed values of parameters in (9) and E 0 * ( 0.9669 , 0.0867 ) , we compute the numerical solutions and depict the time-series curves (Figure 2), which are in accordance with our theoretical findings of Theorem 2 because 1 > D 1 = 0.0994 .
For f = 5 and E * ( 0.9532 , 0.0850 ) , since 1 1 + f T = 0.5714 > D 1 = 0.0985 , by use of Theorem 2 again, it should be still globally stable, which is verified by the numerical graph, see Figure 3. If τ 1 = τ 2 = 0.2 , then the equilibrium status also keeps stable, see Figure 4.
Now, we validate the stochastic results by performing numerical simulations in MATLAB. The stochastic terms are simulated through Milstein’s higher-order method [31]. We take the parameter values as in (9) and σ 1 = σ 2 = 1 and solve the numerical solutions of System (4) (see Figure 5).
This reveals that System (4) has a stationary distribution; that is, the equilibrium status E * ( x * , y * ) remains globally asymptotically stable. We summarize all simulation results in Table 2.
Further, we explore the effects of fear, delays and stochastic parameters on the system dynamics of the equilibrium state of (2) numerically. Set the parameter values of (2) as follows:
c = 0.73 , a 1 = 1.98 , a 2 = 0.01 , m = 0.0695 , e = 0.05 .
Let f = 5 , then, the equilibrium state of (2) is E * ( 0.1701 , 0.4612 ) . The numerical curves imply that it remains globally asymptotically stable (see Figure 6). When the fear is large enough, for example f = 1000 , then the equilibrium state is E ˜ * ( 0.1214 , 0.0381 ) , which is unstable (see Figure 7). From a biological point of view, a too high amount of fear can make the prey exhibit such anti-predation responses as foraging, vigilance and even starvation, reducing the population. Certainly, due to the shortage of prey, the predator population also declines, which leads to a coexistence equilibrium state with lower species densities.
Next, we investigate the influence of time delay and white noise on the system dynamics of E * . Let τ = 0.5 and solve (2) numerically; then we get Figure 8, which indicates the stability of E * fails due to the impact of delays. We alter the values of stochastic parameters to explore the effect of white noise. Let σ 1 = σ 2 = 0.02 ; Figure 9 shows that a small amount of white noise does not influence the numerical stability of E * , whereas with greater amounts of white noise, i.e., σ 1 = σ 2 = 0.5 , the simulation reveals that E * is unstable (see Figure 10). From a biological angle, if the delay and white noise are relatively small, then their effects can be ignored since this kind of negative role is not enough to affect the long-time behavior of species, but if the effects are large enough, then the biomass will change rapidly in a short period of time so that the negative roles of delay and white noise cannot be offset, which makes the system unstable.
We summarize the simulation results in Table 3.

5. Conclusions

For predator–prey systems for which the predator scale is sufficiently large, it is reasonable to use a density-dependent functional response to describe predation, and the fear of prey due to a powerful predator cannot be ignored; hence, we propose a delayed prey–predator model with Beddington–DeAngelis functional response and fear of a predator in a deterministic environment, and then analyze the global asymptotic stability of the equilibrium status of the model with or without time delays. Then, we extend the deterministic model to random environments and get a stochastic model. Applying Lyapunov functional theory, we investigate the distribution stability. Our main results are summarized in Table 4.
The theoretical results are indicated by above examples (see Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5). The important effects of multiple time lags (delayed fear and delayed gestation), fear of predators and environmental white noise are revealed by simulations (Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10). From Figure 1, we find that the parameter value of fear affects the value of the equilibrium status E * ( x * , y * ) . In other words, the equilibrium E * ( x * , y * ) is smaller if the value of fear of the predator is larger. Figure 2, Figure 3, Figure 4 and Figure 5 show the long-time behavior of E * ( x * , y * ) ; i.e., it is globally stable when the conditions in the main results are satisfied. Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 reveal that the parameter values have large impacts on the stability of System (5), which can change the dynamics of the equilibrium state from stable to unstable. We deduce that there may be a threshold for every parameter for system stability. In other words, we numerically observe that factors such as fear of prey, time delay of fear, time delay of gestation, and environmental disturbance can greatly influence the stability of the equilibrium status of System (2). In detail, the main contributions of this research are as follows:
  • Theorems 1 and 2 show the effects of fear on the global stability of Systems (1) and (2), respectively. Remark 2 reveals the effect of delay on system dynamics.
  • Lemma 2 gives the representation of the solution of the stochastic differentialEquation (6). By use of the comparison method together with Lemma 2, the existence of the solution of System (4) is proved, which is distinct from the Lyapunov functional methods applied in some existing research [9,15,32].
  • Theorem 3 shows theoretically the effects of environmental white noise and fear on the distribution stability.
  • The effects of fear, white noise and delays are validated by simulation examples visually and numerically in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10.
In light of our findings, in order to decrease the negative effects of delayed fear and white noise in order to maintain ecological balance, some measures may be undertaken, such as (i) construction of refuge zones to reduce the prey’s fear of predators, (ii) provision of additional food for predators to protect the prey while maintaining the predator population, and (iii) application of some ‘biological control’ strategies such as releasing predators to control prey, and so on.
Some interesting topics deserve further study. Firstly, by simulation examples, some oscillations occur, and there may exist parameter thresholds for stability. Then, it is unknown if the Hopf bifurcations occur at some critical value of fear f or time delay τ 1 and τ 2 . Secondly, due to sudden environmental perturbations such natural biological/ecological disasters and sudden changes in weather (hurricanes), which are usually modeled by environmental Lévy noise, it is necessary to introduce Lévy noise to the population model. Thirdly, there is uncertainty created by interactions between different species and shortage of food resources in real-world systems. A review of the literature shows that this uncertainty plays a key role in using mathematical models to predict population dynamics (see, for example [33,34] and references cited therein). Thus, incorporated uncertainty into our model to examine changes in system dynamics is a very interesting topic to be explored. We leave these subjects for future research.

Author Contributions

Y.S. contributed to the conception, model design, analysis, interpretation, drafting of the manuscript, revision and proofreading. W.K. contributed to simulation analysis, interpretation and proofreading. All authors have read and agreed to the published version of the manuscript.

Funding

The paper is supported by the NSF of China (no. 11861027) and partially by the Project of Education Department of Yunnan Province (2022J0807).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are very grateful to the anonymous referees and the editor for their helpful comments and suggestions, which have greatly improved the presentation of this paper.

Conflicts of Interest

The authors declare that they have no competing interest with regards to this paper.

References

  1. Volterra, V. Fluctuations in the abundance of a species considered mathematically. Nature 1926, 118, 558–560. [Google Scholar] [CrossRef]
  2. Holling, C.S. The functional response of predator to prey density and its role in mimicry and population regulation. Mem. Entomol. Soc. Can. 1965, 45, 1–60. [Google Scholar] [CrossRef]
  3. Shao, Y.F.; Li, Y. Dynamical analysis of a stage structured predator-prey system with impulsive diffusion and generic functional response. Appl. Math. Comput. 2013, 220, 472–481. [Google Scholar] [CrossRef]
  4. Xu, D.S.; Liu, M.; Xu, X.F. Analysis of a stochastic predator-prey system with modified Leslie-Gower and Holling-type IV schemes. Physica A 2020, 537, 122761. [Google Scholar] [CrossRef]
  5. Roy, J.; Barman, D.; Alam, S. Role of fear in a predator-prey system with ratio-dependent functional response in deterministic and stochastic environment. BioSystems 2020, 197, 104176. [Google Scholar] [CrossRef] [PubMed]
  6. Skalski, G.T.; Gilliam, J.F. Functional response with predator interference: Viable aiternative to the Holling type II model. Ecology 2001, 82, 3083–3092. [Google Scholar] [CrossRef]
  7. Beddington, J.R. Mutual interference between parasites of predators and its effect on searching efficiency. J. Anim. Eccol. 1975, 44, 331–341. [Google Scholar] [CrossRef]
  8. DeAngelis, D.L.; Goldsten, R.A.; Neill, R. A model for trophic interaction. Ecology 1975, 56, 881–892. [Google Scholar] [CrossRef]
  9. Xia, Y.X.; Yuan, S.L. Survival analysis of a stochastic predator-prey model with prey refuge and fear effect. J. Biol. Dynam. 2020, 14, 871–892. [Google Scholar] [CrossRef]
  10. Zanette, L.Y.; Allen, A.F.; Clinchy, M. Percieved predation risk reduces the number of spring songbirds produe per year. Science 2011, 334, 1398–1401. [Google Scholar] [CrossRef]
  11. Wang, X.; Zanette, L.; Zou, X. Modelling the fear effect in predator-prey interactions. J. Math. Biol. 2016, 73, 1179–1204. [Google Scholar] [CrossRef] [PubMed]
  12. Elliot, K.H.; Betini, G.S.; Norris, D.R. Fear creats an Alleeeffects: Experimental evidence from seasonal populations. Proc. Biol. Sci. 2017, 284, 20170878. [Google Scholar]
  13. Panday, P.; Samanta, S.; Pal, N.; Chattopadhyay, J. Delay induced multiple stability switch and chaos in a predator-prey model with fear effect. Math. Comput. Simul. 2020, 172, 134–158. [Google Scholar] [CrossRef]
  14. Sahoo, D.; Samanta, G.P. Comparison between two tritrophic food chain models with multiple delays and anti-predation effect. Int. J. Biomath. 2021, 14, 2150010. [Google Scholar] [CrossRef]
  15. Shao, Y. Global stability of a delayed predator-prey system with fear and Holling-type II functional response in deterministic and stochastic environments. Math. Comput. Simul. 2022, 200, 65–77. [Google Scholar] [CrossRef]
  16. Wang, Y.; Zou, X. On a predator-prey system with digestion delay and anti-predation strategy. J. Nonlinear Sci. 2020, 30, 1579–1605. [Google Scholar] [CrossRef]
  17. Xiao, Y.; Chen, L. Modelling and analysis of a predator-prey model with disease in the prey. Math. Biosci. 2001, 171, 59–82. [Google Scholar] [CrossRef]
  18. Kuang, Y. Delay Differential Equations with Applications in Population Dynamics; Academic Press: New York, NY, USA, 1993. [Google Scholar]
  19. Upadhyay, R.K.; Mukhopadhyay, A.; Iyengar, S.R. Influence of environmental noise on the dynamics of a realistic ecological model. Fluct. Noise Lett. 2007, 7, 61–77. [Google Scholar] [CrossRef]
  20. Liu, M.; Wang, K. Stochastic Lotka-Volterra systems with Lévy noise. J. Math. Anal. Appl. 2014, 410, 750–763. [Google Scholar] [CrossRef]
  21. Zhao, Y.; Yuan, S. Stability in distribution of a stochastic hybrid competitive Lotka-Volterra model with Lévy jumps. Chaos Solitons Fractals 2016, 85, 98–109. [Google Scholar] [CrossRef]
  22. Mao, X. Stochastic Differential Equations and Applications; Horwood Publishing: Chichester, UK, 1997. [Google Scholar]
  23. Liu, M.; Wang, K. Global asymptotic stability of a stochastic Lotka-Volterra model with infinite delays. Commun. Nonlinear Sci. Numer. Simul. 2012, 17, 3115–3123. [Google Scholar] [CrossRef]
  24. Liu, Q. The effects of time-dependent delays on global stability of stochastic Lotka-Volterra competitive model. Physica A 2015, 420, 108–115. [Google Scholar] [CrossRef]
  25. Has’minskii, R.Z. Stochastic Stability of Differential Equations; Sijthoff Noordhoff: Alphen annden Rijn, The Netherlands, 1980. [Google Scholar]
  26. Song, X.Y.; Chen, L.S. Optimal harvesting and stability with stage structure for a two species competitive system. Math. Biosci. 2001, 170, 173–186. [Google Scholar] [CrossRef]
  27. May, R. Stability in randomly fluctuating deterministic environments. Am. Nat. 1973, 107, 621–650. [Google Scholar] [CrossRef]
  28. Mao, X.; Yuan, C.; Zou, J. Sttochastic differential delay equations of population dynamics. J. Math. Anal. Appl. 2005, 304, 296–320. [Google Scholar] [CrossRef]
  29. Huang, Y.; Liu, Q.; Liu, Y.L. Global asymptotic stability of a general stochastic Lotka-Volterra system with delays. Appl. Math. Lett. 2013, 26, 175–178. [Google Scholar] [CrossRef]
  30. Das, A.; Samanta, G.P. Modelling the fear effect in a two-species predator-prey system under the influence of toxic substances. Rend. Circ. Mat. Palermo Ser. 2 2021, 70, 1501–1526. [Google Scholar] [CrossRef]
  31. Higham, D.J. An aigorithmic introduction to numerical simulations of stochastic differential equations. SIAM Rev. 2001, 43, 525–546. [Google Scholar] [CrossRef]
  32. Kong, W.L.; Shao, Y. The long time behavior of equilibrium status of a predator-prey system with delayed fear in deterministic and stochastic scenarios. J. Math. 2022, 2022, 3214358. [Google Scholar] [CrossRef]
  33. Waezizadeha, T.; Mehrpooya, A.; Rezaeizadehc, M.; Yarahmadiand, S. Mathematical models for the effects of hypertension and stress on kidney and their uncertainty. Math. Bios. 2018, 305, 77–95. [Google Scholar] [CrossRef]
  34. Brehme, M.; Koschmieder, S.; Montazeri, M.; Copland, M.; Oehler, V.G.; Radich, J.P.; Brmmendorf, T.H.; Schuppert, A. Combined population dynamics and entropy modelling supports patient stratification in Chronic cyeloid leukemia. Sci. Rep. 2016, 6, 24057. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Time-series curves of equilibrium statuses of System (2). The blue line is for E 0 * = E 0 * ( 0.9669 , 0.0867 ) ( f = 0 ) , the green line is for E * = E * ( 0.9532 , 0.0850 ) ( f = 5 ) , and the red line is for E ^ * = E ^ * ( 0.8597 , 0.0726 ) ( f = 50 ) .
Figure 1. Time-series curves of equilibrium statuses of System (2). The blue line is for E 0 * = E 0 * ( 0.9669 , 0.0867 ) ( f = 0 ) , the green line is for E * = E * ( 0.9532 , 0.0850 ) ( f = 5 ) , and the red line is for E ^ * = E ^ * ( 0.8597 , 0.0726 ) ( f = 50 ) .
Mathematics 10 03378 g001
Figure 2. The asymptotic stability of E 0 * ( x * , y * ) of model (2).
Figure 2. The asymptotic stability of E 0 * ( x * , y * ) of model (2).
Mathematics 10 03378 g002
Figure 3. Global stability of E * ( x * , y * ) of Model (2).
Figure 3. Global stability of E * ( x * , y * ) of Model (2).
Mathematics 10 03378 g003
Figure 4. Global stability of E * ( x * , y * ) of Model (2) with τ 1 = τ 2 = 0.2 .
Figure 4. Global stability of E * ( x * , y * ) of Model (2) with τ 1 = τ 2 = 0.2 .
Mathematics 10 03378 g004
Figure 5. The stationary distribution of (4) with σ 1 = σ 2 = 1 : (a) time-series curves of x ( t ) and y ( t ) ; (b) density distribution of x ( t ) ; and (c) density distribution of y ( t ) .
Figure 5. The stationary distribution of (4) with σ 1 = σ 2 = 1 : (a) time-series curves of x ( t ) and y ( t ) ; (b) density distribution of x ( t ) ; and (c) density distribution of y ( t ) .
Mathematics 10 03378 g005
Figure 6. Global stability of E * ( x * , y * ) of Model (2) with f = 5 : (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Figure 6. Global stability of E * ( x * , y * ) of Model (2) with f = 5 : (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Mathematics 10 03378 g006
Figure 7. Equilibrium status E ˜ * of Model (2) with f = 1000 is numerically unstable: (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Figure 7. Equilibrium status E ˜ * of Model (2) with f = 1000 is numerically unstable: (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Mathematics 10 03378 g007
Figure 8. The equilibrium status E * of (2) with τ = 0.5 is numerically unstable: (a) time-series curves of x ( t ) and y ( t ) , (b) is the phase diagram.
Figure 8. The equilibrium status E * of (2) with τ = 0.5 is numerically unstable: (a) time-series curves of x ( t ) and y ( t ) , (b) is the phase diagram.
Mathematics 10 03378 g008
Figure 9. The global stability of E * ( x * , y * ) of (4) with σ 1 = σ 2 = 0.02 : (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Figure 9. The global stability of E * ( x * , y * ) of (4) with σ 1 = σ 2 = 0.02 : (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Mathematics 10 03378 g009
Figure 10. The equilibrium status E * ( x * , y * ) of (4) with σ 1 = σ 2 = 0.5 is numerically unstable: (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Figure 10. The equilibrium status E * ( x * , y * ) of (4) with σ 1 = σ 2 = 0.5 is numerically unstable: (a) time-series curves of x ( t ) and y ( t ) and (b) phase diagram.
Mathematics 10 03378 g010
Table 1. Biological meanings of all parameters.
Table 1. Biological meanings of all parameters.
ParameterBiological Meaning
A 1 Intrinsic growth rate of prey
A 2 Natural mortality rate of predator
B 1 Intra-specific competition coefficient of prey
B 2 Density dependence rate of predator
PMaximal relative increase of predation
QConversion factor
C 1 Half-saturation constant
C 2 Magnitude of interference among predators
fFear effect of prey induced by predator
Table 2. Numerical simulations of (4) with parameters given in (9).
Table 2. Numerical simulations of (4) with parameters given in (9).
ParametersStability of EquilibriaSimulation Figures
FearTime DelaysStochastic EffectEquilibria
000 E 0 * yesFigure 2
500 E * yesFigure 3
50.20 E * yesFigure 4
50.21 E * yesFigure 5
where E 0 * = E 0 * ( 0.9669 , 0.0867 ) and E * = E * ( 0.9532 , 0.0850 ) .
Table 3. Numerical simulations of (4) with parameters given in (10).
Table 3. Numerical simulations of (4) with parameters given in (10).
ParametersStability of EquilibriaSimulation Figures
FearTime DelayStochastic EffectEquilibria
500 E * yesFigure 6
100000 E ˜ * noFigure 7
50.50 E * noFigure 8
500.02 E * yesFigure 9
500.5 E * noFigure 10
where E * = E * ( 0.1701 , 0.4612 ) and E ˜ * = E ˜ * ( 0.1214 , 0.0381 ) .
Table 4. Stability conditions of the equilibrium points E * ( x * , y * ) .
Table 4. Stability conditions of the equilibrium points E * ( x * , y * ) .
Parameters Stability Conditions
Fear ( f ) Stochastic Effect ( σ ) Time Delay ( τ )
nonono 1 > D 1
nonoyes 1 > D 1 + D 3
noyesno 1 > D 1 + σ 2 2 x * 2 , e > σ 2 2 y * 2
noyesyes 1 > D 1 + D 2 + D 3 + σ 2 2 x * 2 , e > D 1 + D 2 + σ 2 2 y * 2
yesnono 1 1 + f T > D 1
yesnoyes 1 1 + f T > D 1 + D 3
yesyesno 1 1 + f T ˜ > D 1 + σ 2 2 y * 2 , e > σ 2 2 y * 2
yesyesyes 1 1 + f T ˜ > D 1 + D 2 + D 3 + σ 2 2 y * 2 , e > D 1 + D 2 + σ 2 2 y * 2
where T = c / a 1 m e ,   T ˜ = c / a 1 m + y * 2 σ 2 2 / 2 e ,   D 1 = c a 1 y * 1 + a 1 x * + a 2 y * ,   D 2 = f ( 1 x * ) 2 ( 1 + f y * ) ,   D 3 = c ( 1 + a 2 y * ) 2 ( 1 + a 1 x * + a 2 y * ) .
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Shao, Y.; Kong, W. A Predator–Prey Model with Beddington–DeAngelis Functional Response and Multiple Delays in Deterministic and Stochastic Environments. Mathematics 2022, 10, 3378. https://doi.org/10.3390/math10183378

AMA Style

Shao Y, Kong W. A Predator–Prey Model with Beddington–DeAngelis Functional Response and Multiple Delays in Deterministic and Stochastic Environments. Mathematics. 2022; 10(18):3378. https://doi.org/10.3390/math10183378

Chicago/Turabian Style

Shao, Yuanfu, and Weili Kong. 2022. "A Predator–Prey Model with Beddington–DeAngelis Functional Response and Multiple Delays in Deterministic and Stochastic Environments" Mathematics 10, no. 18: 3378. https://doi.org/10.3390/math10183378

APA Style

Shao, Y., & Kong, W. (2022). A Predator–Prey Model with Beddington–DeAngelis Functional Response and Multiple Delays in Deterministic and Stochastic Environments. Mathematics, 10(18), 3378. https://doi.org/10.3390/math10183378

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop