Analysis of a stochastic Leslie-Gower predator-prey system with Beddington-DeAngelis and Ornstein–Uhlenbeck process

1.   Introduction

Predator-prey relationships are common in nature and have been the basis for the development of numerous biomathematical models by many researchers. A well-known predator-prey model was proposed by Leslie ^[1,2]. In contrast to the Lotka-Volterra model, the Leslie-Gower model provides a superior representation of the interactions between predators and prey. This model employs a more intricate function to depict this relationship, which more accurately mirrors the real-world situation. Furthermore, the Leslie-Gower model takes into consideration competition among the population of prey which is not mentioned in the Lotka-Volterra model. In order to study the dynamical behavior of populations under certain cases, the Beddington-DeAngelis functional response was introduced to the original model by Yu et al. ^[3]. This model is represented as

The natural growth rates of prey and predators are represented by $r_1$ and $r_2$, respectively. The intraspecific competition coefficient of prey is denoted by $b$. The amount of food provided by prey for the birth of predators is measured by $\beta$. The maximum of the mean reduction rate of prey is denoted by $\alpha$. The degree of protection provided by the environment to predators is measured by $k$. All of the above parameters are positive.

However, since any ecosystem inevitably suffers from environmental noise perturbation, the use of stochastic differential equations can describe the dynamic behavior of populations more accurately, and also benefit to explore the dynamic response and stability of the system under the influence of noise^[4]. There exist some methods to simulate the parameters of a varying environment. The first approach assumes that the parameters can be adequately modeled by a linear function of white noise (see example in ^[5,6,7]). The second approach assumes that the parameters satisfy a mean-reverting stochastic process, i.e., each parameter satisfies a specific stochastic differential equation.

The first method is considered to have limitations^[8]. In the following, mathematical methods will be employed to demonstrate this unreasonableness. It is assumed that a certain parameter of the population satisfies the following equation

where $\bar{r}$, which can be obtained through direct calculation, represents the average value of $r$ over a long term. $W(t)$ is the independent standard Brownian motion defined on a complete probability space $\left\{\Omega, \mathcal{F}, \left\{\mathcal{F}_t\right\}_{t \geq 0}, \mathbb{P}\right\}$ with the $\sigma-$filtration $\left\{\mathcal{F}_t\right\}_{t \geq 0}$ satisfying the usual conditions ^[4], and $\alpha > 0$ denotes the noise density of $W(t)$. We assume that for any time interval $[0, t]$, $\langle {r(t)} \rangle$ is the time average of $r(t)$. According to direct calculation, we obtain the result

where $\mathbb{N}\left(\cdot, \cdot\right)$ is the one dimension Gaussian distribution. It is evident that the average state $\langle {r(t)} \rangle$ on $[0, t]$ has a variance of $\dfrac{\alpha^2}{t}$, which approaches infinity as $t\to 0^+$. This leads to an unreasonable outcome where the stochastic fluctuation of certain parameter of the population $r(t)$ becomes very massive for a small time interval.

Another method is to assume that the parameters follow a mean-reverting stochastic process, i.e. each parameter obeys a certain stochastic differential equation

where $c > 0$ is the speed of reversion and $\sigma > 0$ is the intensity of volatility, respectively. $W(t)$ is the independent standard Brownian motion, which is the same as above. As stated by Mao^[9], $r(t)$ has a unique explicit solution in the form below

It is clearly indicated by the above equation that $r(t)$ follows the Gaussian distribution $\mathbb{N}(\mathbb{E}[r(t)], \mathbb{VAR}[r(t)])$ on $[0, t]$. It can be easily derived that $\mathbb{E}[r(t)] = \bar{r}+[r(0)-\bar{r}]e^{-c t}$ and $\mathbb{VAR}[r(t)] = \dfrac{\sigma^2}{2c}(1-e^{-2c t})$. Furthermore, it is evident that $\mathop {\lim }\limits_{t \to \infty } \mathbb{E}[r(t)] = \bar r$ and $\mathop {\lim }\limits_{t \to \infty } \mathbb{VAR}[r(t)] = 0$. Thus, for certain time intervals, the fluctuation of $r(t)$ will be sufficiently small, which is in line with the continuous perturbation property of environmental noise.

Many scholars have introduced standard white noise into the Leslie-Gower model, for example, a stochastic Leslie-Gower model with standard white noise has been studied by Zhao et al.^[10], and they demonstrated some good mathematical properties with biological significance. However, little research has been done on the Ornstein-Uhlenbeck process. Inspired by the above work, we also assume that the growth rates of both prey and predators are influenced by two mean-reverting Ornstein-Uhlenbeck processes

where $B_1(t)$ and $B_2(t)$ are two independent Brownian motions, $c_1 > 0$ and $c_2 > 0$ are the speed of reversion, and $\sigma_1 > 0$ and $\sigma_2 > 0$ are the intensities of volatility. Then we rewrite the system (1.1) as

The theoretical methods and techniques for dynamical analysis are well-developed. However, it should be noted that there are many essential differences between the methods that analyze the models driven by white noise and those driven by Ornstein-Uhlenbeck. Moreover, the introduction of Beddington-DeAngelis also increases the complexity of the models. For example, system (1.1) with stochastic fluctuation has a unique solution that is global and positive^[11]. However, the solution to system (1.6) is not necessarily positive due to the properties of the Ornstein-Uhlenbeck process. We attempt to develop some suitable methods and theories to obtain some dynamical properties of system (1.6), which are analogous to those of system (1.1) with linear white noise.

Currently, in order to study the dynamical properties of stochastic predator-prey models, many scholars have adopted the Ornstein-Uhlenbeck process as the driving force of stochastic systems. For example, Zhang et al. studied a three species predator-prey model driven by the Ornstein-Uhlenbeck process and demonstrated many important dynamical properties^[12]. Chen et al. studied a Leslie-Gower model driven by the Ornstein-Uhlenbeck process with a modified Holling-II functional response and demonstrated good dynamic properties^[13]. In addition, there are many applications of Ornstein-Uhlenbeck processes to the study of other stochastic systems, for example, Song et al. studied a stochastic SVEIS model with an Ornstein-Uhlenbeck process^[14], Wen et al. studied an SIB cholera model with saturated response rate and the Ornstein-Uhlenbeck process^[15], and Liu studied a stochastic HLIV model with virus production and Ornstein-Uhlenbeck process^[16].

We studied the asymptotic properties of the solution with respect to system (1.6). In subsection 2.1, the existence and uniqueness of the global solution of system (1.6) are proven. Additionally, the ultimate boundedness of system (1.6) is given in subsection 2.2. The existence of the stationary distribution is then shown in subsection 2.3. The extinction to the populations is discussed in subsection 2.4. Finally, our conclusions are verified by numerical simulations in subsection 2.5.

2.   Asymptotic property of the solution

2.1.   Existence and uniqueness of the global solution

For convenience, we need to define two necessary sets: $\mathbb{S}_n = (-n, n)\times(-n, n)\times(-n, n)\times(-n, n)$ and $\mathbb{R}_+^n = \{(x_1, \cdots, x_n)\in \mathbb{R}^n|x_m > 0, 0\le m\le n\}$, where $||\cdot||$ is the Euclidean norm. Next, we prove the existence and uniqueness of the global solution of system (1.6).

Theorem 2.1. For any initial value $(x(0), y(0), {r_1}(0), {r_2}(0))\in\mathbb{R}_+^2\times\mathbb{R}^2$, there exists a unique solution $(x(t), y(t), r_1(t), r_2(t))$ of system (1.6) on $t > 0$, and it will remain in $\mathbb{R}_+^2\times\mathbb{R}^2$ with probability one.

Proof. It is easy to verify that Eq (1.6) satisfies the linear growth condition and the local Lipschitz condition. So there exists a locally unique solution $(x(t), y(t), r_1(t), r_2(t))$ defined on $t\in [0, \tau_e)$ (see ^[9]). $\tau_e$ is the explosion time, so it suffices to prove that $\tau_e = \infty$. We let $n$ be sufficiently large so that both $\ln x(0)$ and $\ln y(0)$, $r_1(0)$ and $r_2(0)$ are in the interval $[-n, n]$, defining the stopping time as

Obviously $\tau_n$ is monotonically increasing as $n\to \infty$. Set $\tau_\infty = \mathop {\lim }\limits_{n \to + \infty } {\tau _n}$; here, $\tau_\infty \le \tau_e$. So, it is only necessary to prove that $\tau_\infty = \infty$. If not, there exists constants $T > 0$ and $\varepsilon \in (0, 1)$ such that $\mathbb{P}\{ {\tau _\infty } \le T\} > \varepsilon$. Therefore, there exists an integer $n_1 > n_0$ such that

To simplify the proof, we omit all bracketed $t$. For any $t\le\tau_n$, a non-negative Lyapunov function $V_0(x, y, r_1, r_2):\mathbb{R}_+^2\times\mathbb{R}^2\to\mathbb{R}_+$ is constructed as

Applying Ito's formula, we have

where

Following our calculations, we have

Integrating both sides of inequality (2.3) from $0$ to ${\tau _n} \wedge T$ and computing expectations, we arrive at

When $\tau_n\le T$, we define $\Omega_n = {\tau_n\le T}$. From inequality (2.1), it follows that $\mathbb{P}(\Omega_n) > \varepsilon$. Observe that for any $\omega\in\Omega_n$, there exists an $n$ such that $\ln x(\tau_n, \omega), \ln y(\tau_n, \omega) = -n \text{ or } n$ and $r_1(\tau_n, \omega), r_2(\tau_n, \omega) = -n \text{ or } n$. By combining Eq (2.4), we derive

such that $n \to \infty$ leads to

The proof of Theorem 2.1 is complete. □

2.2.   Ultimate boundness

Since ecosystems have finite resources, population density cannot grow infinitely and will eventually converge to a certain value over time, and it is necessary to prove theoretically that system (1.6) is ultimately bounded. First, let us give the definition of stochastically ultimate boundedness.

Definition 2.1. ^[17] System (1.6) is said to be stochastically ultimately bounded: if, for any $\epsilon\in(0, 1)$, there is a positive constant $\chi = \chi(\omega)$ such that for any initial value $(x(0), y(0), {r_1}(0), {r_2}(0))$, the solution of system (1.6) has the property that

Next, we present a useful lemma, from which the stochastically ultimate boundedness will follow directly.

Lemma 2.1. Let $\theta\in (0, 1)$, then there is a positive constant $M = M(\theta)$ which is independent of the initial value $(x(0), y(0), {r_i}(0))$, where $i = 1, 2$, such that the solution of system (1.6) has the property that

Proof. Define a non-negative function

Applying the generalized Ito's formula and mathematical expectation to $e^{\lambda_\theta t}V_1$, we obtain

where $\lambda_\theta = \theta \min\{c_1, c_2\}$. Note that

where

Combining with Eq (2.8), we obtain

Then we have

By setting $M(\theta) = {2^{\frac{\theta }{2}}}\theta\dfrac{\psi(\theta)}{\lambda_\theta}$, the result (2.7) is obtained. □

Theorem 2.2. System (1.6) is stochastically ultimately bounded.

Proof. By Lemma 2.1, there exists $M > 0$ such that

Now, for any $\epsilon > 0$, let $\chi = \sqrt{2}\dfrac{\psi(0.5)^2}{4\epsilon^2\lambda_\theta^2}$. Then, by Chebyshev's inequality

Hence,

□

2.3.   Existence of the stationary distribution

In biology, a major goal is to study the behavior of systems over long periods of time. In this part we show that there is a stationary distribution for system (1.6) which might to make long term predictions for system (1.6) under stochastic perturbations. The sufficient conditions for the existence of a stationary distribution of system (1.6) is established in this part. First, let us give a useful lemma.

Lemma 2.2. (Khasminskii^[18]). Consider the stochastic system

Let the vectors $\xi(s, x), \upsilon_1(s, x), \ldots, \upsilon_l(s, x)$ $\left(s \in\left[t_0, T\right], x \in \mathbb{R}^m\right)$ be continuous functions of $(s, x)$ such that, for some constants $M$, the following conditions hold in the entire domain of definition:

Moreover, there exists a non-negative function $W(x)$ such that

where $\mathbb{H}$ is a compact subset defined on $\mathbb{R}^m$ and $LW(x)$ is the diffusion operator of the It$\hat{o}$ process with respect to the non-negative functions $W(x)$^[9]. Then, the Markov process (2.11) has at least one stationary solution $X(t)$, which has a stationary distribution on $\mathbb{R}^m$.

Remark 2.1. According to Xu et al.^[19] the condition (2.12) in Lemma 2 can be replaced by the existence of the globally unique solution to system (1.6).

Before we begin the proof, we define some notation and make some reasonable assumptions.

Definition 2.2. Define $N_0$ to be a natural number that satisfies the following conditions

where

Theorem 2.3. For any initial value $(x(0), y(0), {r_1}(0), {r_2}(0))$, if $\bar{r}_1 + \bar{r}_2 > 0$, then system (1.6) has at least a stationary distribution with the definition of $N_0$.

Proof. By Theorem 2.1, it is easy to know that there is a globally unique solution to system (1.6), so the description of $\mathbb{R}^m$ in the Lemma 2.2 should be changed to $\mathbb{R}_+^2\times\mathbb{R}^2$. Therefore, it is only necessary to verify that the following conditions hold:

We divide the relevant proof into the following two steps.

Step 1 Define the function

Applying Ito's formula, by the definition of $N_0$ we obtain

It should be noted that the function $V_2$ tends towards $\infty$ as $x$ or $y$ approaches the boundary of $\mathbb{R}_+$ or as $||(x, y, r_1, r_2)||\to\infty$. Consequently, there exists a point $(x^0, y^0, r_1^0, r_2^0)$ in the interior of $\mathbb{R}_+^2\times\mathbb{R}^2$, at which the value of function $V_2$ is minimized. A non-negative function $V_3(x, y, r_1, r_2)$ may be constructed as

Combining with Ito's formula, we have

Step 2 Considering a closed set $\mathbb{H}_\varepsilon$ in the form

we define

Let $\epsilon\in(0, 1)$ be a sufficiently small number such that the following inequalities hold:

After that, we will verify $LV_2(x, y, r_1, r_2)\le -1$ for any $(x, y, r_1, r_2)\in \left(\mathbb{R}_{+}^2 \times \mathbb{R}^2\right) \backslash \mathbb{H}_\epsilon$. Noting that $\left(\mathbb{R}_{+}^2 \times \mathbb{R}^2\right) \backslash \mathbb{H}_\epsilon = \bigcup_{k = 1}^6 \mathbb{H}_{k, \epsilon}^c$,

Below we will prove that when $(x, y, {r_1}, {r_2})$ belongs to the complement of ${\mathbb{H}_\epsilon}$ in $(\mathbb{R}_ + ^2 \times {\mathbb{R}^2})$, the value of $L{V_3}$ is less than or equal to $-1$. The proof can be discussed in five cases.

Case $\mathbf{1}$ If $(x, y, r_1, r_2)$ is located in the set defined by $\mathbb{H}_{1, \epsilon}^c$, then one can obtain the corresponding results by combining equation (2.15) and (2.16).

Case $\mathbf{2}$ If $(x, y, r_1, r_2)$ is located in the set defined by $\mathbb{H}_{2, \epsilon}^c$, it follows that similar conclusions could be calculated from (2.15) and (2.16).

Case $\mathbf{3}$ If $(x, y, r_1, r_2)$ is located in the set defined by $\mathbb{H}_{3, \epsilon}^c$, consequently, from (2.15) and (2.16), we can obtain the relevant result.

Case $\mathbf{4}$ If $(x, y, r_1, r_2)$ lie within the set demarcated by $\mathbb{H}_{4, \epsilon}^c$, the relevant conclusions can be deduced through (2.15) and (2.16).

Case $\mathbf{5}$ In the event that $(x, y, r_1, r_2)$ is situated within the set defined by $\mathbb{H}_{5, \epsilon}^c$, the associated findings could be calculated by (2.15) and (2.18).

Case $\mathbf{6}$ If the point $(x, y, r_1, r_2)$ belongs to the complement of the set $\mathbb{H}_{6, \epsilon}^c$, then we can use equations (2.15) and (2.17) to compute the relevant results.

We summarize the above cases and, by Eqs (2.15)–(2.17), it can be concluded that there exists a sufficiently small constant $\epsilon$ such that $L{V_3}(x, y, {r_1}, {r_2})\le -1$ for any $(x, y, r_1, r_2) \in \left(\mathbb{R}_{+}^2 \times \mathbb{R}^2\right) \backslash \mathbb{H}_\epsilon$, where $\epsilon$ satisfies $L{V_3}(x, y, {r_1}, {r_2}) \le -1$ for any $(x, y, r_1, r_2)\in \left(\mathbb{R}_{+}^2 \times \mathbb{R}^2\right) \backslash \mathbb{H}_\epsilon$ where $\epsilon$ satisfies that

And, for any $\Pi_2 > 1$, we set

According to the discussion above, condition (2.13) in Lemma 2 is verified. Therefore, system (1.6) has at least a stationary distribution. □

2.4.   Extinction

Theorem 2.4. For any initial value $(x(0), y(0), r_1(0), r_2(0))\in \mathbb{R}_+^2 \times \mathbb{R}^2$, the solution $(x(t), y(t), r_1(t), r_2(t))$ of system (1.6) has the property that

In particular, if$\; \bar{r}_1 < 0, \bar{r}_2 < 0$, then $x(t), y(t)$ are extinct.

Proof. Applying the $\mathrm{It\hat{o} }$ formula to $\ln x(t), \ln y(t)$, we can get

Integrating from $0$ to $t$, we have

Then, combining the strong law of large numbers^[20] and the definition of the Ornstein–Uhlenbeck process, we have

According to (2.25) and (2.26), we obtain

Then,

When $\bar{r}_1 < 0, \bar{r}_2 < 0$, implying $\lim\limits_{t\to \infty} x(t) = 0, \lim\limits_{t \to \infty} y(t) = 0$, then $x(t), y(t)$ are extinct. Theorem 2.4 is proved. □

2.5.   Numerical simulations

In this section, we will use computer simulations to verify our conclusions. Using the Milstein higher order method^[21], we obtain the discretization equation for system (1.6). The corresponding discretization equation is as

where $\Delta t > 0$ denotes the time increment, and $\eta_j$ and $\xi_j$ are two independent stochastic variables which follow the Gaussian distribution $\mathbb{N}(0, 1)$. Besides, $(x^j, y^j, r_1^j, r_2^j)$ is the corresponding value of the $jth$ iteration of the discretization Eq (2.29). We will use some different combinations of biological parameters in Table 1 to simulate.

From Tables 1 and 2, we choose the combination ($\mathcal{A}_1$) as the value of the biological parameters of system (1.6). Obviously, the existence and uniqueness of a solution of system (1.6) is shown (see Figure 1). We choose the combination ($\mathcal{A}_2$) as the parameters of system (1.6), the result of Figure 2 demonstrates the $\theta th$ order moments of the solutions of system (1.6) are bounded, and system (1.6) is ultimately bounded. Then, we choose the combination ($\mathcal{A}_3$) as the biological parameters of system (1.6), and the stationary distribution is explained in Figure 3. Finally, we simulated the extinction of system (1.6) using the parameter combination ($\mathcal{A}_4$, $\mathcal{A}_5$) (see Figure 4).

3.   Discussion

This paper first introduced a two-species Leslie-Gower model and further reviewed the seminal work of previous scholars using these models to simulate the stochastic effects of population systems. The use of Ornstein-Uhlenbeck enables a more realistic simulation of the stochastic properties of the environment with a relatively stable pattern of variation compared to the methods in ^[5,6,7], which are using standard white noise.

Therefore, we considered and studied a class of Leslie-Gower models that contain Ornstein-Uhlenbeck and Beddington-DeAngelis functional responses in system (1.6).

We have proven many important theoretical results of system (1.6), including the existence and uniqueness of solutions, the ultimate boundedness, the existence of stationary distributions, and the extinction of the populations. We verified the correctness of related conclusions using numerical simulation methods. These theoretical contributions serve to enrich the theory of population dynamics and establish a mathematical basis for the practical application of population dynamics changes.

In summary, this paper proposes and studies a Leslie-Gower population model with environmental fluctuations containing the Ornstein-Uhlenbeck mean-reverting process. This model can describe the stochastic changes in environmental conditions more accurately than previous population models and serves to enrich the theoretical study of the impact of random environmental factors on population dynamics.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was supported by National Natural Science Foundation of China (11401085), Fundamental Research Projects of Chinese Central Universities (2572021DJ04), Heilongjiang Province Postdoctoral Funding Program (LBH-Q21059), and College Students' Innovative Entrepreneurial Training Plan Program (202210225161).

Conflict of interest

The authors declare there is no conflicts of interest.