A stochastic predator–prey model with two competitive preys and Ornstein–Uhlenbeck process
Qun LiuSchool of Mathematics and Statistics, Key Laboratory of Applied Statistics of MOE, Northeast Normal University, Changchun, Jilin Province, People's Republic of ChinaCorrespondence[email protected] View further author information
Article: 2193211
|
Received 06 May 2022, Accepted 18 Nov 2022, Published online: 22 Mar 2023
?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.
Abstract
In this paper, a stochastic predator–prey model with two competitive preys and Ornstein–Uhlenbeck process is formulated and analysed, which is used to obtain a better understanding of the population dynamics. At first, we validate that the stochastic system has a unique global solution with any initial value. Then we analyse the stochastic dynamics of the model in detail, including pth moment boundedness, asymptotic pathwise estimation in turn. After that, we obtain sufficient conditions for the existence of a stationary distribution of the system by adopting stochastic Lyapunov function methods. In addition, under some mild conditions, we derive the specific form of covariance matrix in the probability density near the quasi-positive equilibrium of the stochastic system. Finally, numerical illustrative examples are depicted to confirm our theoretical findings.
Predator–prey interaction with multiple prey is common in nature, for example, Pelican birds in Salton sea eat Tilapia Oreochromis mossambicus with Avian botulism [Citation3]. In order to study the evolution law of predator–prey interaction with multiple preys, a lot of predator–prey models with multiple preys were established and analysed during the past few decades [Citation2, Citation14, Citation16, Citation23, Citation29]. The well-known predator–prey model with competitive interaction between two preys can be expressed as follows:
(1) (1) where and are all positive constants, and represent the population densities of two competitive species and stands for the population density of predator species. The biological significance of the parameters is as follows: and denote the intrinsic growth rates of the two competitive preys and represents the natural death rate of the predator population. The coefficients , and denote the intra-specific competitive coefficients of two prey species , and the predator species , respectively. and represent the capturing rate of the predator at time interval, and are the inter-specific competitive coefficients of two prey species , and , stand for the rate of conversion of nutrients of the predator species. All the parameters in the model (Equation1(1) (1) ) are assumed to be positive constants.
On the other hand, in the natural world, the growth of population is always disturbed by environmental noises [Citation4, Citation5], such as random changes in epidemics, droughts, floods, temperatures, etc. These changes sometimes are so strong that they can change the behaviours of population significantly [Citation32]. May [Citation24] pointed out that due to environmental fluctuations, the birth rates, death rates, competition coefficients, transmission coefficients and other parameters involved in the system should exhibit random fluctuation to some extent. Therefore, the parameters , and should be regarded as random variables involved in randomness. Accordingly, many scholars have introduced white noise in the deterministic population systems to study how the noise affects the dynamics of stochastic population models [Citation9, Citation10, Citation12, Citation13, Citation15, Citation19, Citation21, Citation25] and this way is common.
However, except using the common way to simulate the random disturbances, we can also assume that the parameters involved in the system satisfy the well-known Ornstein–Uhlenbeck process. The Ornstein–Uhlenbeck process is an Itô process which can be described by the following stochastic differential equation:
where θ, μ, σ are all positive constants and θ is a mean reversion constant that denotes the strength of the process attracted by the mean [Citation26], μ is the mean value, σ represents the degree of volatility around the mean value μ caused by shocks, denotes the standard Brownian motion.
By the work of Allen [Citation1], the stochastic mean-reverting process should incorporate the influences of environmental perturbations into the parameters of a system. So in this paper, by considering the second way of introducing randomness in model (Equation1(1) (1) ), we assume that , and are affected by three Ornstein–Uhlenbeck processes, i.e.
(2) (2) (3) (3) (4) (4) where , , are all positive constants, and denotes the speed of reversion; are the intensities of the volatility of process ; represent the mean reversion levels of , i = 1, 2, 3. , and are mutually independent standard Brownian motions defined on a complete probability space with a filtration satisfying the usual conditions [Citation22].
By performing stochastic integral operation to (Equation2(2) (2) )–(Equation4(4) (4) ), respectively, we have
where are the initial values of Ornstein–Uhlenbeck processes , i = 1, 2, 3. It is easy to see that the expected values and the variances of over an interval are as follows
According to the property of Brownian motion, we obtain that follows the distribution and follows the normal distribution . Apparently, the Ornstein–Uhlenbeck processes follow the distribution as t tends to zero, i = 1, 2, 3, which can be seen as a both mathematically and biologically reasonable method to simulate the random impacts of key parameters in a population system. Keeping in mind of the above facts and combining with (Equation2(2) (2) )–(Equation4(4) (4) ), we get the stochastic form corresponding to system (Equation1(1) (1) ) which is presented as follows:
(5) (5) Here the parameters involved in system (Equation5(5) (5) ) are the same as those in system (Equation1(1) (1) ).
The rest of the paper is arranged as follows. In Section 2, the existence and uniqueness of the global solution of the stochastic system (Equation5(5) (5) ) is proved which is necessary to study the dynamics of a population model. In Section 3, the boundedness of the pth moment of solutions to system (Equation5(5) (5) ) is shown. In Section 4, the asymptotic pathwise estimation of the solutions of system (Equation5(5) (5) ) is investigated. In Section 5, sufficient criteria for the existence of a stationary distribution of positive solutions to system (Equation5(5) (5) ) are derived. In Section 6, we gain the exact form of covariance matrix in the probability density around the origin point. In Section 7, through the numerical simulations of several examples, we validate the effectiveness of our theoretical findings. Finally, we present a brief conclusion to end the paper.
Next, we present the following notations which will be used in the whole paper:
Let be the set of all nonnegative functions defined on which means it contains all continuously twice differentiable in x. We use the notation to denote the indicator function of the set . If is a matrix, we use the notation to represent its transpose. If is an invertible matrix, its inverse matrix is denoted by . If G is a square matrix, we use the notation to denote its determinant. Let be the kth sequential principal submatrix of the matrix . In addition, for any two d-dimensional symmetric matrices and , we define
In this sense, if and .
Throughout this paper, the following hypothesis always holds:
Remark 1.1
The Assumption (A) is used to ensure the positive equilibrium of system (Equation1(1) (1) ) exists.
Now we will present basic theories on stochastic differential equations (for more details, we refer the reader to [Citation22]).
Let be a d-dimensional stochastic process on given by the following stochastic differential equation:
where , and is a d-dimensional Brownian motion defined on the complete probability space . We define the differential operator L as follows:
If L acts on a function , then
where and . If , then by means of Itô's formula, we gain
2. Existence and uniqueness of the global solution
To investigate the dynamical behaviour of system (Equation5(5) (5) ), we should first ensure whether the solution is global or not. To this end, we present the following theorem which is related to the existence and uniqueness of the global solution to system (Equation5(5) (5) ).
Theorem 2.1
Given initial value , there is a unique global solution to system (Equation5(5) (5) ) which will remain in almost surely (a.s.).
Proof.
It is noticed that all the coefficients of system (Equation5(5) (5) ) are locally Lipschitz continuous, then for any initial value , there exists a unique local solution of system (Equation5(5) (5) ) on the interval a.s., where represents the explosion time [Citation22]. Next, let be large enough such that , , , , and all lie within the interval . For each integer , define the stopping time as below [Citation22]
It is clear that is increasing as . Denote by , then a.s. If a.s. is true, then a.s. and a.s. In other words, to complete the proof, a.s. needs to be proved. If this assertion is not correct, then there are two constants T>0 and such that
Hence, there exists an integer such that
(6) (6) Define a -function by
where and are positive constants to be determined later. The nonnegativity of this function can be derived from
Applying Itô's formula [Citation22] to V, we obtain
(7) (7) where is defined by
where in the second inequality, we have used the Young inequality
Choose , such that and , then we have
Here is a positive constant which is independent of the variables , , , , and . The rest of the proof is similar to that of Zhang and Yuan [Citation31] and so it is omitted here. The proof is confirmed.
3. The pth moment boundedness
In this section, the boundedness of the pth moment of solutions to system (Equation5(5) (5) ) is shown by establishing some Lyapunov functions.
Theorem 3.1
Given initial value and any integer , there exists a positive constant such that the solution of system (Equation5(5) (5) ) has the property
Proof.
For any integer , define
where
Applying Itô's formula [Citation22] to , we have
Let
then A is positive definite. Let be the minimum eigenvalue of the matrix A, i.e., then we get
In addition, we have
Hence (8) (8) Similarly, applying Itô's formula [Citation22] to , we obtain (9) (9) where
From (Equation8(8) (8) ) and (Equation9(9) (9) ), it follows that (10) (10) where in the second inequality, we have used the Young inequality
and
Hence, for any , by means of (Equation10(10) (10) ), we obtain (11) (11) where
Furthermore, according to (Equation11(11) (11) ), we have (12) (12) By integrating (Equation12(12) (12) ) from 0 to t and then taking expectations on both sides, one has
Letting leads to that
The proof is confirmed.
Remark 3.1
In view of Theorem 3.1, it is easy to obtain that there exists a positive constant such that
4. Asymptotic pathwise estimation
In this section, we will verify pathwise properties of the solution of system (Equation5(5) (5) ).
Theorem 4.1
Given initial value , the solution of system (Equation5(5) (5) ) has the property
In addition
Proof.
According to Theorem 3.1, one can see that there is a positive constant such that
which together with the continuity of , and implies that there is a positive constant such that
(13) (13) For sufficiently small , , in view of (Equation13(13) (13) ), we get
where
So
Let be arbitrary. By Chebyshev's inequality [Citation22], we obtain
By Borel–Cantelli Lemma [Citation22], one obtains that for almost all ,
(14) (14) holds for all but finitely many k. As a result, there exists a for almost all , for which (Equation14(14) (14) ) holds whenever . Thus, for almost all , if and ,
Therefore
Letting gives that
and hence
Since the integer p is arbitrary and greater than 0, it can be seen that
This completes the proof.
5. Existence of a stationary distribution
In this section, sufficient conditions which are used to ensure the existence of a stationary distribution of the stochastic system (Equation5(5) (5) ) are presented.
Let be a homogeneous Markov process defined in which is described by the following stochastic differential equation:
(15) (15)
Assume that the vectors , are continuous functions with respect to X and satisfy the following conditions:
There is a constant B such that
(16) (16) There exists a nonnegative -function in such that outside some compact set. Then system (Equation15(15) (15) ) admits at least one stationary distribution on .
Remark 5.1
We can use the global existence of solutions of (Equation15(15) (15) ) to replace the conditions (Equation16(16) (16) ) in Lemma 5.1 and the conclusion can be found in Remark 5 of Xu et al. [Citation30].
Theorem 5.1
Let Assumption (A) hold. If the following conditions hold
and is a sufficiently small constant satisfying
then system (Equation5(5) (5) ) has at least one stationary distribution on , where and is the positive equilibrium of system (Equation1(1) (1) ).
Proof.
By Theorem 2.1, we have proved that there is a global solution to system (Equation5(5) (5) ). Hence, in order to complete the proof of Theorem 5.1, we only need to prove the condition . That is to say, we need to construct a nonnegative -function and a compact set such that for all .
Define a nonnegative -function as follows:
where is the positive equilibrium of system (Equation1(1) (1) ) which satisfies the following algebraic equations
(17) (17) For convenience, denote by
Applying Itô's formula [Citation22] to and using (Equation17(17) (17) ), we have (18) (18) where in the above inequality, we have used the Young inequality
and is a sufficiently small constant to ensure the coefficient in front of is negative.
Similarly, we obtain (19) (19) (20) (20) where , are sufficiently small constants to ensure the coefficients in front of , are negative, respectively.
Moreover, we have (21) (21) (22) (22) and (23) (23) Then according to (Equation18(18) (18) )–(Equation23(23) (23) ), we have
where in the second inequality we have used the Young inequality
and
such that
Let
and note that if the condition
holds, then the ellipsoid
lies entirely in . Thus, we can take to be a neighbourhood of the ellipsoid with , where and denotes the compact closure of . So for , , which implies that the condition in Lemma 5.1 holds. According to Lemma 5.1, we obtain that there is a solution of system (Equation5(5) (5) ) which is a stationary Markov process. This completes the proof.
6. Local density function for system (5)
In this section, we will get the specific form of covariance matrix in the probability density around the quasi-positive equilibrium of the stochastic system (Equation5(5) (5) ).
6.1. Equivalent transformation of system (5)
Firstly, we define a quasi-positive equilibrium involved in stochasticity by the equations
(24) (24) By solving Equation (Equation24(24) (24) ), we obtain that if Assumption (A) holds, then
where is the positive equilibrium of system (Equation1(1) (1) ).
Let . In view of Itô's integral and system (Equation5(5) (5) ), the corresponding linearised system of (Equation5(5) (5) ) around can be expressed as follows
(25) (25) where , , , , , , , , , , , .
6.2. Local density function for system (5)
In this subsection, we will obtain the specific form of covariance matrix in the probability density around the quasi-positive equilibrium of the stochastic system (Equation5(5) (5) ). To this end, we need to introduce a definition and a lemma.
The characteristic polynomial of the square matrix is defined as , then is called a Hurwitz matrix if and only if has all negative real-part eigenvalues, i.e.
where the complementary definition , j>n.
For the three-dimensional real algebraic equation , where , and
If is a three-dimensional symmetric matrix,
then is a positive definite matrix, where
Here is called the standard matrix.
Theorem 6.1
Let be a solution of system (Equation25(25) (25) ) with the initial value . If Assumption (A) holds, then there is a local density function of a multivariate normal distribution around the origin point , which takes the form
where denotes the covariance matrix of which is positive definite and takes the following form
If , , , then
If , , , then
If , , , then
If , , , then
If , , , then
If , , , then
If , , , then
If , , , then
here
and
Proof.
Let , ,
With these notations, system (Equation25(25) (25) ) can be written in the following equivalent form:
(26) (26) On the basis of Gardiner [Citation6], system (Equation26(26) (26) ) has a unique density function , which can be determined by the following six-dimensional Fokker–Planck equation
(27) (27) Next, we will give the accurate expression of the density function by solving Equation (Equation27(27) (27) ). Note that under a stationary case, then (Equation27(27) (27) ) becomes
Additionally, because the diffusion matrix H of system (Equation26(26) (26) ) is a constant matrix, then by the work of Gardiner [Citation6], we get
(28) (28) Because of the mutual independence of Brownian motions , and , in view of the finite independent superposition principle, system (Equation28(28) (28) ) can be equivalent to the sum of solutions of the following three sub-equations
where , , , are their corresponding solutions. One can easily get that and .
Let
In order to prove the matrix is a Hurwitz matrix, in view of Definition 6.1, we only need to show that all the eigenvalues of have negative real parts. To this end, define the characteristic equation of by
where
By direct computation, one has
By means of Lemma 6.1, we obtain that the matrix has all negative real-part eigenvalues.
Next, we will accomplish the proof of Theorem 6.1 in three steps.
Step 1. Consider the algebraic equation
(29) (29) where .
Define , where the ordering matrix is given by
Direct calculation leads to that
Define , where the elimination matrix takes the form
then we obtain
where . In view of the value of , we will divide the following discussion into two parts.
Case 1. If , consider the standard transformation matrix
such that
where
Moreover, define , where the standard transformation matrix takes the form
Direct calculation gives that
where
In addition, an equivalent algebraic equation of (6.6) can be described by
Let
where , then we have
(30) (30) By solving Equation (Equation30(30) (30) ), we obtain
where .
Noticeably, the matrix is positive semi-definite and its submatrix
is positive definite. Thus, the matrix is also positive semi-definite and there is a positive constant such that
Case 2. If , consider the standardized transformation matrix
such that
where
Then continue to do a matrix similarity transformation to , choose the matrix
such that
where
Furthermore, an equivalent algebraic equation of (6.6) can be described by
Denote by
then
(31) (31) By solving Equation (Equation31(31) (31) ), we obtain
It is obvious that the matrix is positive semi-definite and its submatrix
is also positive semi-definite. So the matrix is positive semi-definite and there is a positive constant such that
Step 2. Consider the algebraic equation
(32) (32) where .
Let , where the ordering matrix takes the form
Then
Let , where the elimination matrix takes the form
then we get
where . By means of the value of , we will divide the following discussion into two parts.
Case 1. If , consider the standardized transformation matrix
such that
where
Then continue to do a matrix similarity transformation to , choose the matrix
such that
where
Moreover, an equivalent algebraic equation of (Equation32(32) (32) ) can be described by
Let
where , then we have
(33) (33) By solving Equation (Equation33(33) (33) ), we have
where .
It is easy to see that the matrix is positive semi-definite and its submatrix
is positive definite. Accordingly, the matrix is also positive semi-definite and there is a positive constant such that
Case 2. If , consider the standardized transformation matrix
such that
where
Then continue to do a matrix similarity transformation to , choose the matrix
such that
where
Furthermore, an equivalent algebraic equation of (6.6) can be described by
Let
then
(34) (34) By solving Equation (Equation34(34) (34) ), we obtain
Obviously, the matrix is positive semi-definite and its submatrix
is also positive semi-definite. So the matrix is positive semi-definite and there is a positive constant such that
Step 3. Consider the algebraic equation
(29) (29) where .
Let , where the ordering matrix takes the form
By a simple computation, we have
Choose the elimination matrix
such that
where . According to the value of , we will divide the following discussion into two parts.
Case 1. If , consider the standardized transformation matrix
such that
where
Then continue to do a matrix similarity transformation to , choose the matrix
such that
where
Besides, an equivalent algebraic equation of (6.6) takes the following form
Let
where , we derive
(35) (35) By solving Equation (Equation35(35) (35) ), we have
where .
It is easy to see that the matrix is positive semi-definite and its submatrix
is positive definite. Therefore, the matrix is also positive semi-definite and there is a positive constant such that
Case 2. If , consider the standardized transformation matrix
such that
where
Then continue to do a matrix similarity transformation to , choose the matrix
such that
where
Furthermore, an equivalent algebraic equation of (6.6) can be described by
Denote by
then
(36) (36) By solving Equation (Equation36(36) (36) ), we obtain
It is clear that the matrix is positive semi-definite and its submatrix
is also positive semi-definite. Hence the matrix is positive semi-definite and there is a positive constant such that
Now we prove the matrix in Equation (Equation28(28) (28) ) is positive definite. If , , , then the covariance matrix
is positive definite. Similarly, we can show that under the other seven cases, the covariance matrix Σ is also positive definite. This completes the proof.
7. Numerical simulations
In this section, we will present some numerical simulations to show the effectiveness of our theoretical findings. By employing the Milstein higher-order method introduced in [Citation8], we obtain the corresponding discretization equations of system (Equation5(5) (5) ) which are given by:
(37) (37) where is the corresponding value of the jth iteration of the Equation (Equation37(37) (37) ). is the time increment which is positive, are the noise intensities, are mutually independent Gaussian random variables which follow the distribution for . We choose real parameter values from the previous references, and all the values of biological parameters are given in Table .
Table 1. List of parameter values of system (Equation5(5) (5) ).
To show the existence of a stationary distribution numerically, we choose and the other parameter values are given in Table . By simple computation, we obtain
and
That is to say, the conditions of Theorem 5.1 hold. According to Theorem 5.1, we get that system (Equation5(5) (5) ) admits a stationary distribution which indicates that all the species are coexistent a.s. in the long time. See Figure .
Figure 1. The left column displays the trajectory of the species , and of system (Equation5(5) (5) ) and their corresponding deterministic model (Equation1(1) (1) ) with noise intensities . The right column shows the frequency histograms and marginal density functions of , and in system (Equation5(5) (5) ).
Example 7.2
To illustrate the existence of the probability density near the origin point , we choose and the other parameter values are presented in Table . Direct calculation gives that and
In other words, the conditions of Theorem 6.1 hold. Accordingly, system (Equation5(5) (5) ) has a Gaussian probability density near the origin point . Hence, in view of Theorem 6.1, we calculate the specific form of the covariance matrix Σ,
and the corresponding probability density takes the form
Figure shows this.
Figure 2. Numerical simulations for: (i) the frequency histogram fitting density curves of , , , , and of system (Equation5(5) (5) ) with 200,000 iteration points, respectively. (ii) The marginal probability densities of , , , , and of system (Equation5(5) (5) ). All of the parameter values are the same as in Figure .
8. Conclusion
Predator–prey system is an important component of a real ecosystem. In fact, predator–prey interaction with multiple preys is common and more complicated in the nature. In this work, we concentrate on a stochastic predator–prey model with two competitive preys. In general, existing literature always simulated environmental fluctuations with white noise or coloured noise. While in this paper, we employ another way to introduce stochastic perturbations to simulate random interferences in the environment, i.e. the Ornstein–Uhlenbeck process. To ensure the mathematical and biological rationality of the stochastic system (Equation5(5) (5) ), we first show that the solution of the system is global. Then we study the dynamics of system (Equation5(5) (5) ) in detail. More precisely, analytical findings of this paper are presented as below:
Given initial value and any integer , there exists a positive constant such that the solution of system (Equation5(5) (5) ) has the property
Given initial value , the solution of system (Equation5(5) (5) ) has the property
In addition
Let Assumption (A) hold. If the following conditions hold
and is a sufficiently small constant satisfying
then system (Equation5(5) (5) ) has at least one stationary distribution on , where and is the positive equilibrium of system (Equation1(1) (1) ).
Let be a solution of system (Equation25(25) (25) ) with the initial value . If Assumption (A) holds, then there is a local density function of a multivariate normal distribution near the origin point , which takes the form
where the exact form of covariance matrix Σ can be obtained in the proof process of Theorem 6.1.
However, some other topics of interest are worthy of further consideration. For example, in this paper, we only suppose that the parameters satisfy the Ornstein–Uhlenbeck process, as a matter of fact, if we suppose that other parameters involved in system (Equation5(5) (5) ) satisfy the Ornstein–Uhlenbeck process, what kind of dynamic behaviour should system (Equation5(5) (5) ) exhibit? Secondly, as far as we know, most of key parameters in population systems are usually affected by coloured noise [Citation27]. Hence, the corresponding generalised stochastic system (Equation5(5) (5) ) with discrete Markov switching should be studied. Finally, in this paper, we do not consider the effects of spatial diffusion factor, if we take into account the factor, i.e. we formulate a stochastic partial differential equation model, then the system (Equation5(5) (5) ) will show more abundant dynamic behaviours. These matters will be the subject of our future work.
Acknowledgments
The author read and approved the final manuscript.
Availability of data and materials
Data sharing is not applicable to this article as no datasets are generated or analysed during the current study.
Disclosure statement
No potential conflict of interest was reported by the author(s).
Additional information
Funding
The author thanks the National Natural Science Foundation of P.R. China [grant number 12001090] and the Jilin Provincial Science and Technology Development Plan Project [grant number YDZJ202201ZYTS633].
References
E. Allen, Environmental variability and mean-reverting processes, Discrete Contin. Dyn. Syst. Ser. B 21 (2016), pp. 2073–2089.
H.J. Alsakaji, S. Kundu, and F.A. Rihan, Delay differential model of one-predator two-prey system with Monod-Haldane and Holling type II functional responses, Appl. Math. Comput. 397 (2021), p. 125919.
B. Han, B. Zhou, D. Jiang, T. Hayat, and A. Alsaedi, Stationary solution, extinction and density function for a high-dimensional stochastic SEI epidemic model with general distributed delay, Appl. Math. Comput. 405 (2021), p. 126236.
C. Ji, D. Jiang, and N. Shi, Analysis of a predator-prey model with modified Leslie-Gower and Holling-type II schemes with stochastic perturbation, J. Math. Anal. Appl. 359 (2009), pp. 482–498.
D. Jiang, N. Shi, and X. Li, Global stability and stochastic permanence of a non-autonomous logistic equation with random perturbation, J. Math. Anal. Appl 340 (2008), pp. 588–597.
Y. Li and H. Gao, Existence, uniqueness and global asymptotic stability of positive solutions of a predator-prey system with Holling II functional response with random perturbation, Nonlinear Anal.68 (2008), pp. 1694–1705.
Y. Li and L. Zhang, The stochastic interactions between predator and prey under Markovian switching: Competitive interaction between multiple prey, J Nonlinear Sci. Appl. 10 (2017), pp. 5622–5645.
Q. Liu and Q. Chen, A note on the stationary distribution of a three-species food web stochastic model with generalist predator, Appl. Math. Lett. 114 (2021), p. 106929.
M. Liu and P.S. Mandal, Dynamical behavior of a one-prey two-predator model with random perturbations, Commun. Nonlinear Sci. Numer. Simul. 28 (2015), pp. 123–137.
J. Lv and K. Wang, Asymptotic properties of a stochastic predator-prey system with Holling II functional response, Commun. Nonlinear Sci. Numer. Simul 16 (2011), pp. 4037–4048.
P.S. Mandal, Characterization of positive solution to stochastic competitor-competitor-cooperative model, Electron. J. Differ. Equations 2013 (2013), pp. 1–13.
K.S. Mathur, A prey-dependent consumption two-prey one predator eco-epidemic model concerning biological and chemical controls at different pulses, J. Franklin Inst. 353 (2016), pp. 3897–3919.
H. Wang, D. Jiang, T. Hayat, A. Alsaedi, and B. Ahmad, Stationary distribution of stochastic NP ecological model under regime switching, Physica. A. 549 (2020), p. 124064.
S. Wang, Z. Wang, C. Xu, and G. Jiao, Sensitivity analysis and stationary probability distributions of a stochastic two-prey one-predator model, Appl. Math. Lett. 116 (2021), p. 106996.
D. Xu, Y. Huang, and Z. Yang, Existence theorems for periodic Markov process and stochastic functional differential equations, Discrete Contin. Dyn. Syst. 24 (2009), pp. 1005–1023.
X. Zhang and R. Yuan, A stochastic chemostat model with mean-reverting Ornstein-Uhlenbeck process and Monod-Haldane response function, Appl. Math. Comput. 394 (2021), p. 125833.
S. Zhang, S. Yuan, and T. Zhang, A predator-prey model with different response functions to juvenile and adult prey in deterministic and stochastic environments, Appl. Math. Comput. 413 (2022), p. 126598.