Analysis of a stochastic fear effect predator-prey system with Crowley-Martin functional response and the Ornstein-Uhlenbeck process

1.   Introduction

A predator-prey model is an interdependent and restrictive survival model for different populations in nature. It is significant to ecology, and this model has served as the foundation for numerous research projects. The classic predator-prey model was first proposed by Lotka ^[1] and Volterra^[2]. Many functional responses, such as Beddington-DeAngelis ^[3,4], Leslie-Gower ^[5], Crowley-Martin, etc., have been introduced into predator feeding models in order to to study the relationships between populations better. The functional responses, i.e., how a predator consumes the prey species, are a central component of the theory on a consumer's resource interactions and have a significant effect on the dynamical properties^[6]. Functional responses can be classified into two main categories: prey-dependent^[7,8] and both prey-dependent and predator-dependent ^[9]. First, Crowley and Martin^[10] proposed the predator-prey model with Crowley-Martin-type functional response. The Crowley-Martin functional response, incorporating both prey and predator abundances, provides a more realistic perspective from an ecological standpoint ^[11]. In 2014, Meng et al.^[12] studied a predator-prey system with Crowley-Martin functional response and stage structure for prey. The local stability of equilibria (two boundary equilibria and a positive equilibrium) has been analyzed. However, functional responses just reflect what may happen regarding direct killing. With the development of ecology, some scholars have found that predators not only affect the number of prey through direct killing, but also change the behavior and physiology of prey, thus reducing the number of prey^[13,14]. Actually, when the prey is faced with a predator, they are able to sense a crisis and develop a fear of the predator, known as the fear effect. The fear felt by the prey produces a chain reaction in various ecosystems, thus changing the stability of the system. Therefore, it is necessary to study this fear effect. Wang et al.^[15] added the fear factor to a prey's birth rate $\mu \left({x, y} \right) = \frac{r}{{1 + ky}}$ and the result obtained show that the fear effect can interplay with maturation delay between juvenile prey and adult prey in determining the long-term population dynamics. Das et al.^[16] investigated the effect of the fear function with exponential form (i.e., $\mu \left({x, y} \right) = r{e^{-ky}}$) on prey when a predator is provided with additional food under environmental perturbations. The results show that the higher the fear coefficient, the higher the bait population increases, and a lower predator population tends to extinction. Kumar et al.^[17] combined the fear effect $\mu \left({x, y} \right) = \frac{r}{{1 + ky}}$ with Hassell-Varley functional response to analyze the stability and bifurcation of a predator-prey system, and studied the dynamics of the system in the presence of fear of predation risk. Li and Tian ^[18] proposed a predator-bait model with exponential fear effect (i.e., $\mu \left({x, y} \right) = r{e^{-ky}}$) and Hassell-Varley functional response. The application of a feedback control strategy demonstrated that for predators forming a tight group, appropriately increasing the fear level could stabilize the system. Sarkara et al.^[19] analyzed a prey-predator system introducing the cost of fear function $f\left({\alpha, \eta, P} \right)$ into prey reproduction with Holling type-II functional response. In 2023, a deterministic predator-prey model with a fear effect and Crowley-Martin functional response was proposed and discussed by Zhang^[20] et al., as follow:

where $x\left(t\right)$ and $y\left(t\right)$ reflect the prey population density and the predator population at time $t$. For other parameters, see Table 1. All parameters in the model are assumed to be positive constants, where $\frac{r}{{1 + fy\left(t \right)}}$ represents the fear effect of the prey, and $\frac{{\beta x}}{{\left({1 + ax\left(t \right)}\right)\left({1 + by\left(t \right)}\right)}}$ represents the Crowley-Martin functional response term, which indicates that the interference between predators not only exists when predators handle the prey, but also when they look for the prey.

In nature, however, the development of populations can be disturbed by a variety of uncertain environmental factors, and deterministic population models do not take into account the influence of these random factors and are therefore not suitable for describing reality. For this reason, for the study of populations, it is necessary to add a stochastic disturbance term to the deterministic model. The nature of the model will also change, and it is important to study the nature of the stochastic model.

Environmental noise naturally affects population systems in nature. Many parameters in ecological dynamics should fluctuate around some average values. Mao et al.^[21] demonstrated that even small amounts of environment noise can have a large impact on species populations, which means that stochastic population models can provide additional authenticity compared to deterministic population models. Therefore, it is generally assumed that environmental noise primarily affects the fundamental parameters of the model for studying the dynamic properties of ecosystems in different environments ^[22]. Based on the fact that population death rates are easily affected by environmental fluctuations, we assume $r$ and $m$ in the stochastic predator-prey model are two random variables $r\left(t \right)$ and $m\left(t \right)$.

There are two common methods for simulating small disturbances in the environment in accordance with the current literature. The most common method to describe environmental disturbances is to introduce white noise into the deterministic model^{[23,24,25,26]}. Another method is to incorporate the mean-reverting Ornstein-Uhlenbeck process to simulate environmental perturbations, which has been demonstrated to be a practical and biologically meaningful method. Let us consider the first method, which introduces white noise into the known deterministic model. Assume that Gaussian linear white noise perturbs the intrinsic growth rate and the mortality rate in the model.

where $\bar r$ and $\bar m$ signify the long-term average levels of $r\left(t \right)$ and $m\left(t \right)$, respectively. $B_i\left(t \right)$, i = 1, 2, represents two independent standard Brownian motions defined on a complete probability space $\left\{ \varOmega, \mathcal {F}, \left\{ \mathcal{F}_t \right\} _{t\ge 0}, {\mathbb{P}} \right\}$ with a filtration $\left\{ \mathcal{F}_t \right\} _{t\ge 0}$ adhering to the usual conditions. Additionally, $\alpha_i > 0$, i = 1, 2, indicates the noise density of $B_i\left(t \right)$. We assume that for any time interval $\left[ 0, t \right]$, we have

where $\left < r\left(t \right) \right >$ and $\left < m\left(t \right) \right >$ are the time average of $r\left(t \right)$ and $m\left(t \right)$. ${\mathbb{N}}\left(\cdot, \cdot \right)$ denotes one-dimensional normal distribution. From the above formula, it is easy to see that the variance of $\left < r\left(t \right) \right >$ and $\left < m\left(t \right) \right >$ are $\frac{{{\alpha ^2}}}{t}$, which tends to infinity as $t \to {0^ + }$. This means that the mean value of the perturbation parameter will be more variable in a short amount of time. The use of Gaussian linear white noise to simulate small disturbances in the environment is unreasonable.

According to the above literature, in Section 2, we add environmental perturbations to the predator-prey model with fear effect (1.1) and introduce some lemmas and assumptions to assist the proof below. Section 3 shows several dynamical properties of the system (2.3), including the existence and uniqueness of global solutions, ultimate boundedness, and the existence of the stationary distribution. Additionally, we obtained sufficient conditions for species extinction. In Section 4, numerical simulations are conducted to verify the theoretical results. Finally, we present some conclusions in Section 5.

2.   Model formulation and preliminaries

2.1.   Model formulation

Now, the Ornstein-Uhlenbeck process is introduced into the deterministic model. In other words, each parameter adheres to a certain stochastic differential equation (SDE). When we directly disturb the contact rate $m$ using the Ornstein-Uhlenbeck process, it may yield negative values due to the inherent properties of this type of process, meaning that non-negativity cannot be guaranteed. Motivated by Allen's work ^[27], we propose that the variable $lnm$ is influenced by the Ornstein-Uhlenbeck process, which can be described by the following stochastic differential equation. According to the literature, it will be determined by the following formula:

where $\beta _i > 0$ and $\sigma _i > 0$ (i = 1, 2) represent the speed of reversion and the volatility intensity, respectively. As stated by Mao^[28], by performing random integration operations, we can obtain the following unique solution:

which shows that the variables

According to the above discussion, note that

and

In addition, ${\mathbb{VAR}}\langle {\ln m\left(t \right)} \rangle = \frac{{\sigma _2^2}}{3}t + O\left({{t^2}} \right)$. This indicates that the variation of $\ln m\left(t \right)$ will be adequately minor within a small interval, which is consistent with the facts. Therefore, this method is justifiable for modeling the random effects of key parameters from both biological and mathematical viewpoints^[29].

As a result, introducing the Ornstein-Uhlenbeck process to perturb parameters $r$ and $m$ is more appropriate than Gaussian linear white noise for reflecting the actual situation. Based on the above analysis, by combining model (1.1) and model (2.1), we let $g = lnm$, and we can obtain a stochastic model of the following form:

In this paper, the model establishes a stochastic fear effect predator-prey model with the Crowley-Martin functional response and the Ornstein-Uhlenbeck process, which provides greater stability in environmental variability and the ability of organisms to respond to external changes. Second, the Ornstein-Uhlenbeck process is able to better model environmental perturbations compared to Brownian motion since species do not grow rapidly over a short period of time. The population density of Brownian motion may grow rapidly, which is not consistent with the facts. The mean-reverting Ornstein-Uhlenbeck process is employed to represent minor environmental fluctuations, which offers a more realistic approach compared to assuming that population parameters follow a linear distribution in Gaussian white noise. Numerous scholars have utilized the Ornstein-Uhlenbeck process to study the dynamic properties of stochastic predator-prey models. For instance, Liu ^[30] analyzed a stochastic predator-prey model with two competitive preys and the Ornstein-Uhlenbeck process. Zhou et al.^[31] formulated and analyzed a stochastic nonautonomous population model with Allee effects and two mean-reverting Ornstein-Uhlenbeck processes. Additionally, Liu and Jiang ^[32] explored a stochastic logistic model by incorporating diffusion with two Ornstein-Uhlenbeck processes, which is a stochastic nonautonomous system.

2.2.   Preliminaries

For the convenience of the proof, we first define two sets ${\mathbb{G}}_n = \left(-n, n \right) \times \left(-n, n \right)$ and ${\mathbb{R}}_{+}^{n} = \left\{ \left(x_1, ...x_n \right) \in {\mathbb{R}}^n|x_k > 0, 1\le k\le n \right\}$. Then, we consider the following form of an $m$-dimensional stochastic differential equation (SDE) to introduce Lemma 2.1.

According to Khasminskii ^[33], we will give the existence of stable solutions of system (2.3) through the following lemma.

Lemma 2.1. Let the vectors $\zeta \left(s, x \right), \nu _1\left(s, x \right), \nu _2\left(s, x \right), \cdot \cdot \cdot, \nu _m\left(s, x \right) \left(s\notag\right.\left.\in \left[ t_0, T \right], x\in {\mathbb{R}}^m \right)$ be continuous functions of $\left(s, x \right)$, such that for some constants M, the following conditions hold in the entire domain of definition:

Moreover, ${\mathbb{E}}$ is a compact subset defined on ${\mathbb{R}}^m$. So there exists a non-negative function $V^*\left(x \right)$ such that

where ${\mathbb{E}}$ is a compact subset defined on ${\mathbb{R}}^m$. Then, the Markov process (2.4) has at least one stationary solution $D\left(x \right)$, which has a stationary distribution $\iota \left(\cdot \right)$ on ${\mathbb{R}}^m$.

Remark 1. Based on Remark 5 of Xu et al^[34], conditions (2.6) and (2.7) in Lemma 2.1 can be replaced by the global existence of solutions of system (2.3).

Definition 2.2.1. The solution of system (2.3) is said to be stochastically and ultimately bounded, if for any $\varepsilon \in \left({0, 1} \right)$, there is a positive number $\psi = \psi \left(\omega \right)$ such that for any initial value $x\left(0 \right), y\left(0 \right), r\left(0 \right), g\left(0 \right)\in {\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$, the solution of system (2.3) satisfies

Definition 2.2.2. Define $\varPi _1$ to be a natural number that satisfies the following conditions $M_1$,

where

Lemma 2.2. [Strong law of large numbers]: Let $M = {\left\{ {{M_t}} \right\}_{t \ge 0}}$ be a real-valued continuous local martingale vanishing at $t = 0$. Then $\mathop {\lim }\limits_{t \to \infty } {\langle {M, M} \rangle _t} = \infty$, $\mathop {\lim }\limits_{t \to \infty } \frac{{{M_t}}}{{{{\langle {M, M} \rangle }_t}}} = 0$ a.s., and $\mathop {\lim }\limits_{t \to \infty } \sup \frac{{{{\langle {M, M} \rangle }_t}}}{{{M_t}}} < \infty$, $\mathop {\lim }\limits_{t \to \infty } \frac{{{M_t}}}{t} = 0$ a.s. More generally, if $A = {\left\{ {{A_t}} \right\}_{t \ge 0}}$ is a continuous adapted increasing process such that $\mathop {\lim }\limits_{t \to \infty } {A_t} = \infty$ and $\int_0^\infty {\frac{{d{{\langle {M, M} \rangle }_t}}}{{{{\left({1 + {A_t}} \right)}^2}}}} < \infty$, then, $\mathop {\lim }\limits_{t \to \infty } \frac{{{M_t}}}{{{A_t}}} = 0$ a.s.

Remark 2. If ${\bar \varphi _1} > 0, {\bar \varphi _2} < 0$, the population $x\left(t \right), y\left(t \right)$ are weakly persistent, and it is clear from Theorem 3.4 that the survival and extinction of the population is only related to the mean ${\bar \varphi _1}, {\bar \varphi _2}$.

3.   Results

First, we prove the relevant properties of the solution of system (2.3) through the following theorem. It is important to note that in nature, because $x$ and $y$ represent the population size of the species of the system (2.3), they cannot take on negative values. Therefore, it is necessary to demonstrate both the existence of global solutions $\left(x\left(t \right), y\left(t \right), r\left(t \right), g\left(t \right) \right)$ for system (2.3) and the positivity properties of $x\left(t \right), y\left(t \right)$.

3.1.   Existence and uniqueness of the global solution

Theorem 3.1. For any initial value condition $\left(x\left(0 \right), y\left(0 \right), r\left(0 \right), g\left(0 \right) \right) \in {\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$, system (2.3) has a unique solution $\left(x\left(t \right), y\left(t \right), r\left(t \right), g\left(t \right) \right)$ on $t > 0$, and it will remain in ${\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$ with a probability of one.

Proof. For any initial value $\left({x\left(0 \right), y\left(0 \right), r\left(0 \right), g\left(0 \right)} \right)\in {\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$, it is can be easily demonstrated that the coefficients of the equations in the system (2.3) satisfy the local Lipschitz conditions. Therefore, there exists a unique local solution $\left(x\left(t \right), y\left(t \right), r\left(t \right), g\left(t \right) \right) \in {\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$ on $\left[ 0, \tau _e \right)$, where $\tau _e$ represents the explosion time^[21].

To prove that the model has a global positive solution, we need only show that $\tau _e = \infty$ a.s. By defining a necessary set ${\mathbb{G}}_n = \left(-n, n \right) \times \left(-n, n \right)$, we can always find a sufficiently large integer $n_0$ such that $\left(\ln x\left(0 \right), \ln y\left(0 \right), r\left(0 \right), g\left(0 \right) \right) \in {\mathbb{G}}_{n_0}$. For any integer $n\ge n_0$, we define a stopping time set $\tau _n$ as follows:

Clearly, $\tau _n$ monotonically increased as $n$ increased. For convenience, let ${\tau _\infty } = \mathop {\lim }\limits_{n \to \infty } {\tau _n}$ and $inf\varnothing = \infty$, which implies $\tau _{\infty}\le \tau _e$. To prove Theorem 3.1, it suffices to verify $\tau _{\infty} = \infty\ $ a.s. Consider the contradiction; That is, $\tau _{\infty} < \infty\ $ a.s. This implies that there exist constants $T > 0$ and $\varepsilon \in \left(0, 1 \right)$ such that ${\mathbb{P}}\left\{ \tau _{\infty}\le T \right\} > \varepsilon$. Hence, there exists a positive number $n_1 > n_0$, such that

For any $t\le \tau _n$, using the inequality $x-1\ge \ln x\ \left(x > 0 \right)$, we can construct a non-negative $C^2$-function $V\left(x\left(t \right), y\left(t \right), r\left(t \right), g\left(t \right) \right)$ as follows:

Applying $It\hat{o}$'s formula, we have

By defining

we can obtain that

Following our calculations, we obtain

By integrating the inequality from 0 to $\tau _n\land T$ and subsequently taking the expectation of both sides of the inequality (3.5), we get

For all $n_1 > n_0$, let $\Omega_n = \left\{ \tau _n\le T \right\}$, and then we have ${\mathbb{P}}\left(\Omega_n \right) \ge \varepsilon$. Note that for any $\omega \in \Omega_n$, $lnx$, $lny$, $r$ or $g$ equal either $-n$ or $n$, so there is

where $I_{\Omega_n\left(\omega \right)}$ represents the characteristic function. As $n\rightarrow \infty$, we have

which leads to a contradiction. Therefore, we have $\tau _{\infty} = \infty$. This concludes the proof of Theorem 3.1. □

3.2.   Ultimate boundedness

Due to the limited resources in ecosystems, population density cannot increase indefinitely and will eventually stabilize at a certain value over time. It is essential to theoretically demonstrate that system (2.3) is ultimately bounded. First, we will define stochastically ultimate boundedness^[35] in Definition 2.2.1.

Lemma 3.1. For any initial value $x\left(0 \right), y\left(0 \right), r\left(0 \right), g\left(0 \right)\in {\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$, the solution $\left({x\left(t \right), y\left(t \right), r\left(t \right), g\left(t \right)} \right)$ of system (2.3) has the property

Let $q \in \left({0, 1} \right)$, where $M\left(q \right)$ is a positive constant independent of the initial value $x\left(0 \right), y\left(0 \right), r\left(0 \right), g\left(0 \right)$.

Proof. We define a non-negative $C^2$-function $V_1\left(x\left(t \right), y\left(t \right), r\left(t \right), g\left(t \right) \right) :{\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2\to \mathbb{R}$ by

Applying the generalized $It\hat o$ formula, we obtain

To simplify the notation, the subsequent proof replaces $x\left(t \right)$, $y\left(t \right)$, $r\left(t \right)$, and $g\left(t \right)$ with $x$, $y$, $r$, and $g$, respectively, so we have

Taking the mathematical expectation of $e^{\eta t}V_1$, we obtain

where $\eta = q\min \left\{ \beta _1, \beta _2\right\}$. Noting that

Combining formula (3.10) and formula (3.11), we obtain

and then we have

By setting $M\left(q \right) = 2^{\frac{q}{2}}\frac{q\kappa\left(q \right)}{\eta}$, the result (3.8) holds. □

Theorem 3.2. The solutions of system (2.3) are stochastic and ultimately bounded.

Proof. According to Lemma 3.1, $M\left(q \right)$ exists such that $\mathop {\lim }\limits_{t \to \infty } \sup {\mathbb{E}}\sqrt {|\left({x, y} \right)|} \le M\left(q \right)$. Now, Chebyshev's inequality is applied. For any $\varepsilon > 0$, let $\psi = \sqrt 2 \frac{{\kappa{{\left({0.5} \right)}^2}}}{{4{\varepsilon ^2}{\eta ^2}}}$. Then, we can get

Then, $\mathop {\lim }\limits_{t \to \infty } \sup P\left({\left| {\left({x, y} \right)} \right| > \psi } \right) \le \frac{M}{{\frac{M}{\varepsilon }}} = \varepsilon$. This completes the proof of Theorem 3.2. □

3.3.   Existence of a stationary distribution

In the field of biology, a major objective is to analyze the behavior of systems over long periods. This section seeks to establish sufficient conditions for the existence of a stationary distribution in system (2.3). This will demonstrate the persistence of each species in system (2.3). By Theorem 3.1, it is easy to know that there is a globally unique solution to system (2.3), so the description of ${\mathbb{R}}^{m}$ in Lemma 2.1 should be changed to ${\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$.

Theorem 3.3. For any initial value $\left(x\left(0 \right), y\left(0 \right) r\left(0 \right), g\left(0 \right) \right) \in {\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$, system (2.3) has a stationary distribution with the definition 2 of $M_1$ on ${\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$.

Proof. First, a $C^2$-function $V_2\left(x, y, r, g \right) :{\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2\rightarrow {\mathbb{R}}$ is defined by

By applying $It\hat o$'s formula to $V_2$ and using the definition of $M_1$, we obtain

Then,

According to the expression of function $V_2\left(x, y, r, g \right)$, it can be clearly seen that when $x$ and $y$ tend to infinity, function $V_2\left(x, y, r, g \right)$ will also become infinite. Thus, we can obtain a point $\left(x^0, y^0, r^0, g^0 \right)$ inside ${\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$, where the value of $V_2\left(x, y, r, g \right)$ at that point will reach a minimum value. Combining the above discussion and the requirements of the constructed function in Lemma 2.1, we can construct a non-negative $C^2$-function $V_3\left(x, y, r, g \right)$, whose specific expression is as follows:

According to the application principle of $It\hat o$'s formula, it is known that for the function $V_2\left(x, y, r, g \right)$ under study, adding a constant $V_2\left(x^0, y^0, r^0, g^0 \right)$ to the end does not affect the expression of the result. Therefore, $V_2\left(x, y, r, g \right)$ and $V_3\left(x, y, r, g \right)$ have the same operator, that is to say,

Subsequently, a closed set ${\mathbb{E}}_{\varepsilon}$ is constructed as follows:

Let $\varepsilon$ be a numerical value within the range of 0 and 1. It is of such minute magnitude that all three subsequent inequalities can be fulfilled.

where

In order to simplify the proof, we partitioned the complement of the closed set into 6 distinct regions, namely $\left({\mathbb{R}}_{+}^2 \times {\mathbb{R}}^2\right) \backslash {\mathbb{E}}_{\varepsilon} = U_{i = 1}^{6} {\mathbb{{\mathbb{E}}}}_{i}$, where

Now, we only need to demonstrate that $\mathcal{L}V_3\left(x, y, r, g \right) \le -1$ for all values of $\left(x, y, r, g \right) \in \left({\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2 \right) \backslash {\mathbb{E}}_{\varepsilon}$. Considering the partition of the above complement set, we prove this through the following six cases.

Case 1. If $\left(x, y, r, g \right) \in {\mathbb{E}}_{1, \varepsilon}$, then one can derive the corresponding results by combining Eqs (3.14) and (3.15), and we obtain

Case 2. If $\left(x, y, r, g\right) \in {\mathbb{E}}_{2, \varepsilon}$, it follows that similar conclusions can be calculated from (3.14) and (3.15), and we obtain

Case 3. If $\left(x, y, r, g\right) \in {\mathbb{E}}_{3, \varepsilon}$, consequently, from (3.14) and (3.15), we obtain

Case 4. If $\left(x, y, r, g\right) \in {\mathbb{E}}_{4, \varepsilon}$, from (3.14) and (3.15), we obtain

Case 5. If $\left(x, y, r, g\right) \in {\mathbb{E}}_{5, \varepsilon}$, from (3.14) and (3.16), we obtain

Case 6. If $\left(x, y, r, g\right) \in {\mathbb{E}}_{6, \varepsilon}$, from (3.14) and (3.17), we obtain

In summary of the aforementioned cases, it can be concluded that there exists a sufficiently small constant $\varepsilon > 0$ such that $\mathcal{L}V_3\left(x, y, r, g \right) \le -1$ for any $\left(x, y, r, g \right) \in \left({\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2 \right) \backslash {\mathbb{E}}_{\varepsilon}$.

Where ${{\Pi _2} \le-1},$

In addition, if $\varPi _2 > 1$, then $\varepsilon$ not only needs to satisfy the aforementioned condition, but also needs to fulfill the following additional conditions:

This confirms Condition (2.7) in Lemma 2.1. Therefore, system (2.3) possesses a stationary distribution on ${\mathbb{R}}_{+}^{2}\times {\mathbb{R}}^2$. This completes the proof of Theorem 3.3. □

3.4.   Extinction

Theorem 3.4. We define

when ${{\bar \varphi }_1} < 0$ and ${{\bar \varphi }_2} > 0$, and then $x\left(t\right)$ and $y\left(t\right)$ are extinct.

Proof. Based on Eq (1.2) and the definition of the Ornstein-Uhlenbeck process, we obtain

Equation (3.20) clearly indicates that $r\left(t \right)$ and $m\left(t \right)$ follow the Gaussian distribution ${\mathbb{N}}\left({{\mathbb{E}}\left[ {r\left(t\right)}\right], {\mathbb{VAR}}\left[ {r\left(t \right)} \right]} \right)$, ${\mathbb{N}}\left({{\mathbb{E}}\left[ {m\left(t\right)}\right], {\mathbb{VAR}}\left[ {m\left(t \right)} \right]} \right)$ on $\left[ {0, t} \right]$. It is easy to infer that

Therefore, ${\beta _i}\int_0^t {{e^{-{\alpha _i}\left({t-s} \right)}}d{B_i}\left(s \right)} \sim {\mathbb{N}}\left({0, \frac{{\beta _i^2}}{{2{\alpha _i}}}\left({1-{e^{-2{\alpha _i}t}}} \right)} \right)$, $i = 1, 2.$ Then, it is equivalent to $\frac{{\beta _i}}{{\sqrt {2{\alpha _i}} }}\sqrt {1-{e^{-2{\alpha _i}t}}} \frac{{d{B_i}\left(t \right)}}{{dt}}.$ According to Chen et al.^[36], we let ${\gamma_i }\left(t \right) = \frac{{\beta _i}}{{\sqrt {2{\alpha _i}} }}\sqrt {1-{e^{-2{\alpha _i}t}}} \frac{{d{B_i}\left(t \right)}}{{dt}}$, where ${{B_i}\left(t \right)}$ represents a standard Brownian motion. Equation (3.20) can be written as:

Subsequently, we modify system (2.3) accordingly:

Applying $It\hat{o}$'s formula to $lnx\left(t\right)$ and $lny\left(t\right)$, we can get

and from Eq (3.23), we obtain

Then taking the upper limit, we get

According to Lemma 2.2, $\int_0^t {{\gamma _2}\left(s \right)d} {B_2}\left(s \right)$ follows from the strong law of large numbers theorem for martingales. Therefore, $\mathop {\lim }\limits_{t \to \infty } \frac{1}{t}\int_0^t {{\gamma _2}\left(s \right)d} {B_2}\left(s \right) = 0$. Based on the definition of the Ornstein-Uhlenbeck process, if ${{\bar \varphi }_1} + {\varepsilon _1} < 0$, then $\mathop {\lim }\limits_{t \to \infty } x\left(t \right) = 0$. Similarly, when ${{\bar \varphi }_2} > 0$, then $\mathop {\lim }\limits_{t \to \infty } y\left(t \right) = 0$. Theorem 3.4 is proved. □

4.   Computer simulations

Thus far, we have rigorously demonstrated several of the dynamic properties of the system. By constructing Lyapunov functions, we demonstrated the global existence and uniqueness of the solutions. We also derived estimates on the upper bounds of the moments of the solutions. Furthermore, we proved the existence of stationary distributions of the solutions. Next, we will verify our conclusions through numerical simulations.

First, we discretize system (2.3) using the Milstein scheme of higher-order discretization. The discretization equation is as follows:

Immediately afterward, it is necessary to introduce the biological significance of the parameters related to the process being modeled. $\Delta t > 0$ represents the time increment, and ${\xi _i}$, ${\psi _i}$ are two independent stochastic variables that adhere to the standard Gaussian distribution. In addition, $\left({{x^i}, {y^i}, {r^i}, {g^i}} \right)$ denotes the value corresponding to the $i$th iteration of the discretization Eq (4.1), where $i = 1, 2, \cdot \cdot \cdot$. Computer simulations can then be performed to gain insights into the dynamics of the biological system (see Table 1).

Based on the biological significance of the above parameters, we set reasonable values from Jørgensen^[37] in Table 2.

Example 1. We chose the combination $(\mathscr{C}_1)-(\mathscr{C}_3)$ in Table 2 as the biological parameter values for the system (2.3) and made the total number of iterations ${T_{\max }}$ = 2000. We can obtain Figure 1.

Remark 3. From Figure 1, the first line graph represents the trend of population $x$, $y$ whose growth rate is disturbed by the OU process, which together with Theorem 3.1 shows that the corresponding populations $x\left(t \right)$, $y\left(t \right)$ have been fluctuating and growing around the mean value under the influence of various stochastic disturbances. Figure 1 represents $r$, $g$ of the OU process disturbance, showing that the OU process disturbance makes the growth rate fluctuate randomly. Therefore, the growth rate is not a determinant of the population density, and other factors such as environmental disturbances also have a great influence on the survival of the population. It can be seen that different combinations of coefficients have different solutions, all of which exist and are unique. Therefore, the conclusion of Theorem 3.1 can be verified.

On the basis of Theorem 3.1, we next study the effect of environmental perturbations on the population. We chose the parameter combination $(\mathscr{C}_1)$ and varied the value of the parameter ${\sigma _i}$, $i = 1, 2.$ We assumed that ${\sigma _i}$ was 0.02, 0.08 and 0.12, respectively. From Figure 2, we can see that the environmental perturbation has a great influence on the survival of the population. In order to verify the existence and uniqueness of the solution of the system (2.3) more clearly, we performed 100 simulations of Theorem 3.1. All 100 paths are shown as gray lines in Figure 3, and the green solid line represents the average of the 100 simulated paths. It can be seen that different combinations of coefficients have different existence and unique solutions, and thus, the conclusion of Theorem 3.1 can be verified.

Example 2. We chose the combination of $(\mathscr{C}_1)-(\mathscr{C}_3)$ as the biological parameter values for the system (2.3), with a total iteration count of ${T_{\max }}$ = 1500. We can obtain Figure 4.

Remark 4. From Figure 4, we can see the expected values of the three coefficient combinations $(\mathscr{C}_1)-(\mathscr{C}_3)$. These are less than the upper limit $K\left(q \right)$ and this upper limit is not infinite. As time t increases, the probability $P$ gradually stabilizes and becomes greater than a constant. So $\mathop {\lim }\limits_{t \to \infty } \sup P\left\{ {\sqrt {{x^2} + {y^2}} \le \psi } \right\} \ge 1-\varepsilon$. This indicates that Theorem 3.2 holds. From a biological point of view, since environmental resources are finite, no population of any organism can grow indefinitely, which is consistent with reality.

Example 3. We chose the $(\mathscr{C}_1)-(\mathscr{C}_3)$ combination as the biological parameter of the system (2.3). Let the iteration count ${T_{\max }}$ = 2000. Figure 5 shows that system (2.3) has a stationary distribution, and Theorem 3.3 is proven.

Remark 5. From Figure 5, it can be seen that population $x\left(t \right)$ and population $y\left(t \right)$ approximately obey the normal distribution, which indicates that the growth of the population will finally reach a stable state under the influence of the random environmental disturbances received.

To further verify the conclusions of Theorems 3.2 and 3.3, we follow Theorem 3.1 to select 100 paths for simulation. In Figure 6, the left figure shows the simulation results of Theorem 3.2, and the right figure shows the simulation results of Theorem 3.3. Obviously, both theorems are valid.

Example 4. We chose the combination as the biological parameter value (2.3) for the system. We made the total number of iterations ${T_{\max }}$ = 1000. We can obtain Figure 7. From the figure, it can be seen that system (2.3) has extinction, and Theorem 3.4 holds.

5.   Conclusions

In this paper, the model establishes a stochastic predation model with a fear effect, Crowley-Martin functional response, and the Ornstein-Uhlenbeck process. Then, we have proved the existence and uniqueness of solutions, the ultimate boundedness, the existence of stationary distributions, and extinction. Under conditions of not too much environmental noise, the population is able to keep fluctuating around the mean value, and with limited environmental resources, there is a limit to the growth of the population. The model has a stationary distribution when the parameters satisfy certain conditions. In the case of $r-\frac{{\beta _1^2}}{{4{\alpha _1}}} < 0$ and $m-{\rm{ }}\frac{{\beta _2^2}}{{4{\alpha _2}}} > 0$, the population will be extinct.

Simulations of system (2.3) show that due to the finite nature of environmental resources, the intrinsic growth rate also fluctuates around the mean level. The $q$-order moment of the model is less than a certain value. From the biological point of view, due to the limited resources in the natural environment, any biological population will not grow indefinitely, and the model maintains the healthy growth, in line with the biological law. The population $x\left(t \right)$, $y\left(t \right)$ approximately obeys a normal distribution, which in biological sense indicates the persistence of predators and bait; if $r-\frac{{\beta _1^2}}{{4{\alpha _1}}} < 0$ and $m-{\rm{ }}\frac{{\beta _2^2}}{{4{\alpha _2}}} > 0$, the population will be extinct. Unfavorable stochastic environmental disturbances also accelerate the extinction of the population.

Compared to previous work, the inclusion of the Ornstein-Uhlenbeck process makes the existence and uniqueness of the global positive solution in the traditional model into the existence and uniqueness of the global solution, and the proofs of the remaining parts of the Lyapunov function are changed in an innovative way. The Ornstein-Uhlenbeck process is used to more accurately represent the stochastic properties of the environment, exhibiting a relatively stable variation pattern. At present, the significance of the Ornstein-Uhlenbeck process has not been widely discussed in population dynamics models. Previous studies mainly focus on the study of Crowley-Martin's periodicity and balance. Therefore, there is some value in studying other properties of the system (2.3).

The models can be applied to develop ecological and species conservation strategies. The Ornstein-Uhlenbeck process can help ecologists understand how to cope with environmental fluctuations, and it can also be applied to biodiversity studies to explore how predation and fear effects affect interrelationships among different species and the service functions of ecosystems. In addition, the analysis of population extinction has deepened our understanding of the biological background of the model, reflecting the profound impact of extinction events on species evolution, ecosystem balance, and biodiversity. It also provides theoretical support for us to develop conservation strategies in practical applications.

In fact, there is still room for improvement in our model. There is a time-delay in the interaction between populations in many natural ecosystems. However, we did not consider the effect of time-delay on the model. The effects of ecological changes in nature are often delayed, so it is necessary to add a time-delay term to the model. This will be improved in the future. Second, there are often some drastic environmental changes in nature, such as volcanic eruptions, earthquakes, etc., which also affect the population. Therefore, in future studies, to make our model more complete, we will add Lévy jumps to the model to simulate drastic environmental changes in nature.

Use of Generative-AI tools declaration

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

Author contributions

Jingwen Cui and Xiaohui Ai: conceptualization, writing-original draft; Hao Liu: software, visualization. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This study was supported by the National Natural Science Foundation of China under Grant (No. 11401085), the Central University Basic Research Grant (2572021DJ04), the Heilongjiang Postdoctoral Grant (LBH-Q21059), and the Innovation and Entrepreneurship Training Program for University Students (S202410225108).

The authors are grateful to the anonymous reviewers for their helpful comments and suggestions.

Conflict of interest

The authors declare that there is no conflict of interest.