Modeling and analysis of a prey-predator system with prey habitat selection in an environment subject to stochastic disturbances

1.   Introduction

The study of population dynamics is an important branch of ecosystem models and has attracted the special interest of ecologists and applied mathematics at historical and contemporary stages ^[1,2]. Mathematical models play a key role in understanding changes in biological systems ^[3,4] and the effect of control measures ^[5]. In ecosystems, predation is an important way for species to interact and depend on each other. The pioneer work in mathematics that describes the interplay or dynamics between predator and prey belongs to the Lotka ^[6] and Volterra ^[7], which is well known as the Lotka-Voterra model. Subsequently, scholars have refined and enhanced the model in the application of different scenarios and proposed different kinds of uptake functions to characterize the predator's hunting ability ^[8,9,10,11]. The uptake function with a Holling type is representative and applicable to predation phenomena of different species, ranging from lower organisms such as algae and unicellular organisms to invertebrates and vertebrates, such as Holling-Ⅰ ^[12], Holling-Ⅱ ^[13], Holling-Ⅲ ^[14], Holling-(n+1) ^[15] and so on. The Holling-Ⅲ functional response showed that the predator's predation reached saturation when the number of prey reached a certain level ^[14], which is consistent with most of the actual situations. Therefore, it is of a great significance to study the dynamics of a predator-prey model with a Holling-Ⅲ functional response function.

Over the long period of evolution, prey have evolved with a series of anti-predator behaviors. The anti-predator behaviors of prey can be classified into two categories: (ⅰ) defensive counterattacks ^[16] and (ⅱ) morphological or behavioral changes, including camouflage ^[17], seeking refuge ^[18,19], and so on. The benefit of anti-predator behaviors is that they reduce the risk of predation. Red colobus monkeys show siege when they are threatened by chimpanzees ^[16], cuttlefish choose to camouflage by matching background features ^[17], northern pigtailed macaques select trees with abundant branches and elevated sleeping locations in order to minimize the threat of predation ^[18], and white-headed langurs opt for cliff edges or caves as their nocturnal habitats, thus seeking refuge from predators such as leopard cats ^[19]. Anti-predator behaviors are widespread in nature and are essential to consider in the modeling process. Tang and Xiao ^[20] proposed a model of the first type of anti-predation behaviors by adding anti-predation terms to reduce the rate of predator growth. It was demonstrated that anti-predator behaviors inhibit predator-prey oscillations and decreased the likelihood of the coexistence between prey and predator ^[21]. However, many animals lack strong defenses against attacks and they often exhibit the second type anti-predation behavior. An anti-predator model for habitat selection or vigilance was proposed by Ives and Dobson ^[22], which takes the form:

The existence of the anti-predator effect modulates the reproductive capacity of prey through the factor $vx$. For example, yellow-eyed juncos exhibit a heightened vigilance by dedicating more time to scanning for predators upon the release of a trained Harris hawk in their vicinity, thus resulting in a reduced time spent foraging for food ^[23]. It is hypothesized that the diminished feeding rates will ultimately either heighten the likelihood of starvation or hinder the ability to successfully rear offspring. Werner et al. discovered that when bass are present, bluegills will select not rich in resources but more hidden habitats, moreover, bluegills reared in the presence of bass attain only approximately 80$\%$ of the mass compared to those raised without the bass present ^[24]. Additionally, predators pay a price for the anti-predator effect, which is defined by the term $e^{-\varepsilon v}$ on the predator's hunting ability. Based on this consideration, a prey-predator model with a prey habitat selection is considered in this study.

In reality, population systems are typically subjected to a wide range of random disturbances. May has pointed out that the rate of growth, environmental carrying capacity, coefficient of competition, and other pertinent system parameters are all impacted by environmental noise to differing degrees ^[25]. The deterministic model has certain limitations, in some situations, frequently yields disappointing results. To more accurately characterize and predict the population ecosystems, it is necessary to investigate population models with stochastic disturbances. In natural system, the stochastic disturbances can be described by white noise ^[26], colored noise ^[27], Levy jumps ^[28], Markov switching ^[29], and others. Population systems are often subject to environmental fluctuations. Generally speaking, such fluctuations could be modeled by a colored noise. Moreover, if the colored noise is not strongly correlated, then it can be approximately modeled by a white noise, and the approximation works quite well ^[30]. Mao et al. observed that even the smallest white noise can suppress the population explosion and provided a classical method for the uniqueness of the solution ^[31]. Ji and Jiang ^[32] explored a stochastic predatory system that incorporated Beddington-DeAngelis uptake function and discussed the asymptotic property of the system model. Liu et al. ^[33] put forward a prey-predator system model which followed the Holling type-Ⅱ uptake function in an environment that involved stochastic disturbances. Zhang et al. ^[34] analyzed a prey-predator system model by considering the Holling type-Ⅱ uptake function and hyperbolic mortality in an environment subject to stochastic disturbances. Liu ^[35] analyzed the dynamics of a stochastic regime-switching predator-prey model with modified Leslie-Gower Holling-type Ⅱ schemes and prey harvesting. Li et al. ^[36] proposed a Melnikov-type method for chaos in a class of hybrid piecewise-smooth systems with impact and noise excitation under unilateral rigid constraint. In the present investigation, it is assumed that the environment is influenced by white noise in the modeling process.

The structure of this paper is organized as follows: in Section 2, a prey-predator system model with prey habitat selection and stochastic disturbances is formulated, which is followed by a presentation of some basic notations, definitions, and crucial lemmas used in this study; in Sections 3 and 4, the basic properties and dynamic behavior of the system model with or without stochastic disturbances are investigated; in Section 5, numerical simulations with discussions are presented to illustrate the main results; and finally, the work is summarized and further research directions are put forward.

2.   Mathematical model and basic knowledge

2.1.   Mathematical model

When prey habitat selection is considered, the prey-predator model can be described as follows:

where

● $x$ — prey's densities;

● $y$ — predator's densities;

● $r$ — prey's intrinsic growth rate;

● $K$ — prey's environmental carrying capacity;

● $v$ — intensity of anti-predator effect;

● $c$ — conversion efficiency;

● $b$ — interspecific competition coefficient;

● $d$ — natural mortality rate;

● $m$ — efficiency of anti-predator behavior;

● $e^{-m_{0}}$ — maximum predation rate.

Considering the changes in the natural environment and the universality of random disturbances, we introduce white noise into the Model (2.1) in order to analyze the effect of random disturbances on the Model (2.1). Let

where $\sigma_{i}$ represents the white noise intensity, which describes the strength of random disturbance in the environment. $B_{i}(t)$ stands for a Brownian motion, $i = 1, 2$. Then, the model subject to stochastic disturbances is represented as follows:

All the parameters in Model (2.1) and Model (2.2) are positive. In addition, it is assumed that $v\leq r$ and $c\geq \underline{c}\triangleq ade^{vm+m_0}$ in Model (2.1) for biological restriction.

2.2.   Basic knowledge

Let us denote the following:

Consider the following:

where ${\mathbf B}(t)$ is a standard Brownian motion defined on $\left(\Omega, \mathcal{F}, \{\mathcal F_{t}\}_{t\geq0}, {P}\right)$ equipped with $\{\mathcal{F}_{t}\}_{t\geq0}$, where ${\mathbf F}({\mathbf U}(t), t)\in\mathbb{R}^{l}\times[0, +\infty)$, ${\mathbf G}({\mathbf U}(t), t)_{l\times q}$ is a matrix. Define $\mathcal{L}$ as follows:

For ${V}({\mathbf U}(t), t)\in C^{2, 1}(\mathbb{R}^{l}\times[0, +\infty), \mathbb{R}_{+})$, there is

where $V_{t} = \frac{\partial V}{\partial t}$, $V_{{\mathbf U}} = \left(\frac{\partial V}{\partial U_{1}}, \frac{\partial V}{\partial U_{2}}, \cdots, \frac{\partial V}{\partial U_{l}}\right)$, $V_{{\mathbf UU}} = \left(\frac{\partial^{2}V}{\partial U_{i}\partial U_{j}}\right)_{l\times l}$. Then,

Definition 1 (Stochastically ultimate boundedness^[31]). The solution $(u(t), v(t))$ is defined as having the property of stochastically ultimate boundedness if for any $\varepsilon \in(0, 1)$, there exists a constant $\chi > 0$ such that

Definition 2 (Extinction and persistence^[26]). The species $u$ is defined as extinction when $\lim \limits_{t\rightarrow +\infty}u(t) = 0$, weakly persistence (in the mean) when $\langle u(t)\rangle^{*} > 0$, strong persistence (in the mean) when $\langle u(t)\rangle_{*} > 0$, persistence (in the mean) when $\langle u(t)\rangle > 0$, and weak persistence when $\lim \limits_{t\rightarrow +\infty}\sup u(t) > 0$.

Assume that ${\mathbf U}(t)\in E_{l}$ is a regular time-homogeneous Markov process characterized by the following stochastic differential equation:

The diffusion matrix of the process ${\mathbf U}(t)$ is denoted as follows:

Lemma 1 (Ergodic stationary distribution ^[37]). For a given Markov process ${\mathbf U}(t)\in E_{l}$, it is called possessing a unique ergodic stationary distribution $\mu(\cdot)$ if $\exists I\subset E_{l}$ with boundary $\Gamma$ such that i) $\min(eig(\Lambda({\mathbf U})))$ is bounded in $I$ and some neighborhood of $I$; ii) For $u\in E_{l}\backslash I$, if a path originating from $u$ can reach $I$ in a finite average time $\tau$ and $\sup _{u \in S} E_{u} \tau < \infty$, $\forall S\subset E_{l}$, where $S$ is a compact subset.

Remark 1. To establish the condition i), it is sufficient to demonstrate that $\exists \theta > 0$ such that

To establish the condition ii), it is crucial to demonstrate the existence of a neighborhood $I$ along with $V\in C^{2}$ such that $\mathcal{L} V(u) < 0$, $\forall u\in E_{l} \backslash I$.

Lemma 2 (^[26]). For a given $u(t)\in C[\Omega\times[0, +\infty), R^{+}]$:

(1) If $\exists \rho_{_{0}}, T > 0$ with

for $t\geq T$, then

(2) If $\exists \rho_{_{0}}$, $T$ and $\rho$ with

for $t\geq T$, then $\langle u\rangle_{*}\geq \frac{\rho}{\rho_{0}}.$

Lemma 3 (^[37]). For $T_0 > 0$, $\alpha_1 > 0$ and $\alpha_2 > 0$, it has

3.   Dynamic properties of Model (2.1)

As the population densities are non-negative, the discussion on Model (2.1) is restricted in the region $\mathbb{R}^{2}/\mathbb{R}^2_{-}$. The equilibrium of Model (2.1) satisfies the following:

Obviously, Model (2.1) always has two equilibria: $O(0, 0)$ and $E_{B}(\overline{K}, 0)$, where $\overline{K}\triangleq (1-v/r)K$. $O(0, 0)$ is a saddle and constantly unstable.

Define the following:

Then, $x = 0$ and $y = \psi(x)$ are two $x$-isolines, $y = 0$ and $y = \phi(x)$ are two $y$-isolines.

Since

then it has $\psi'(x)\rightarrow -\infty$ as $x\rightarrow 0$. Denote $\overline{v}\triangleq r(1-{3\sqrt{3/a}}\overline{K}^{-2})$. For $v\geq \overline{v}$, there is $\psi'(x)\leq 0$ for $x\in (0, \overline{K}]$; for $0\leq v < \overline{v}$, it has $\psi'(\overline{K}/3) > 0$. Then, there exists $x_1\in (0, \overline{K}/3)$ and $x_2\in (\overline{K}/3, \overline{K})$ such that $\psi(x_1) = \psi(x_2) = 0$.

Similarly, it has $\phi(0) = -d/b$, $\phi'(x) > 0$ for $x\in [0, \overline{K}]$ and $\phi(x) = 0$ if and only if $x = \overline{x}\triangleq \sqrt{\frac{de^{vm+m_0}}{c-ade^{vm+m_0}}}$. Obviously, $\overline{x} < \overline{K}$ if and only if $c > \overline{c}\triangleq de^{vm+m_{0}}(a+\overline{K}^{-2})$.

Theorem 1. The boundary equilibrium $E_{B}(\overline{K}, 0)$ is globally asymptotically stable for $\underline{c}\leq c\leq \overline{c}$.

Proof. At $E_{B}(\overline{K}, 0)$, the Jacobian matrix is as follows:

and its characteristic roots are as follows:

If $c < \overline{c}$, then there is $\lambda_{2} < 0$, then $E_B(\overline{K}, 0)$ is a node and locally asymptotically stable. If $c = \overline{c}$, there is $\lambda_2 = 0$, then $E_B(\overline{K}, 0)$ is a saddle-node. Since for $c\leq \overline{c}$, there is $\phi(x) < 0$, i.e., $\text{d}y/\text{d}t < 0$ for $x\in [0, \overline{K})$, then $y(t)\rightarrow 0$ as $t\rightarrow \infty$. While $\text{d}x/\text{d}t > 0$ for $x\in [0, K)$ and $y = 0$, then $x(t)\rightarrow K$, i.e., $E_{B}(\overline{K}, 0)$ is globally asymptotically stable. □

For $c > \overline{c}$, there is $\overline{x} < \overline{K}$. In such case, $E_{B}(\overline{K}, 0)$ is unstable.

Theorem 2. For $c > \overline{c}$, Model (2.1) has at least one and no more than three interior equilibria. Moreover, if $\overline{v}\leq v < r$, then the interior equilibrium is unique. In addition, the interior equilibrium $\widehat{E}(\hat{x}, \hat{y})$ is locally asymptotically stable if and only if $\phi(\hat{x}) < \psi'(\hat{x}) < bc\hat{y}/(b\hat{y}+d)$.

Proof. Since there are $\psi(\overline{x}) > 0$ and $\psi(\overline{K}) = 0$ for $c > \overline{c}$, $\phi(\overline{x}) = 0$, $\phi(\overline{K}) > 0$, according to the mediocrity theorem, there exists at least one $x^*\in (\overline{x}, \overline{K})$ such that $\psi(x^*) = \phi(x^*)$, i.e., Model (2.1) has at least one interior equilibrium $E^{*}(x^*, y^*)$, where $y^* = \psi(x^*)$, as illustrated in Figure 1.

For $\overline{v}\leq v < r$, $y = \psi(x)$ is monotonically decreasing on $(\overline{x}, \overline{K})$ and $y = \phi(x)$ is monotonically increasing on $(\overline{x}, \overline{K})$; therefore, $y = \psi(x)$ and $y = \phi(x)$ have a unique intersection point, i.e., the interior equilibrium $E(x^*, \psi(x^*))$ is unique (Figure 1(d)).

For $0 < v < \overline{v}$, $y = \phi(x)$ is monotonically increasing and tends to saturation $y = y_{\sup}\triangleq \frac{ce^{-vm-m_{0}}-ad}{ab}$. Alternatively, $y = \psi(x)$ is monotonically decreasing on $(0, x_1)$, monotonically increasing on $(x_1, x_2)$, and monotonically decreasing on $(x_2, \overline{K})$, which forms a S-type on $(0, \overline{K})$. Thus, $y = \psi(x)$ and $y = \phi(x)$ have no more than three intersections, i.e., Model (2.1) has no more than three interior equilibria (Figure 1(b)).

For the given interior equilibrium $\widehat{E}(\hat{x}, \hat{y})$, the Jacobian matrix is as follows:

Since

if $\phi(\hat{x}) < \psi'(\hat{x}) < bc\hat{y}/(b\hat{y}+d)$, then there are $\mathrm{Trace}(J_{\widehat{E}(\hat{x}, \hat{y})}) < 0$ and $\mathrm{Det}(J_{\widehat{E}(\hat{x}, \hat{y})})$ (i.e., the interior equilibrium $\widehat{E}(\hat{x}, \hat{y})$ is locally asymptotically stable). □

4.   Dynamical properties of Model (2.2)

In this segment, we will examine various characteristics of the stochastic prey-predator model.

4.1.   Global positivity

Theorem 3. For $\forall (x_{0}, y_{0})\in \mathbb{R}^{2}_{+}$, Model (2.2) possesses a unique, global solution $(x(t), y(t))$ for $t > 0$, which keeps positive with a probability with one.

Proof. Denote $x(t) = e^{p(t)}$, $y(t) = e^{q(t)}$. Consider the following:

It can be verified that Model (4.1) fulfills the Local Lipschitz Condition; then, for given $p_{0} = \ln x_{0}, q_{0} = \ln y_{0}$, it guarantees the existence of a unique, locally positive solution $(p(t), q(t))$ on the interval $[0, \tau_{e})$, where $\tau_{e}$ is the time of explosion. The It$\mathrm{\hat{o}}$ formula implies that $(x(t), y(t)) = (e^{p(t)}, e^{q(t)})$ is exactly the solution of Model (2.2).

Next, it is sufficient to demonstrate that $\tau_{e} = \infty$. Let $n_{0}$ be an integer satisfy $(x_{0}, y_{0})\subset\left[\frac{1}{n_{0}}, n_{0}\right]$. For $n\geq n_{0}$, define the following:

For an empty set $\phi$, define $\inf\phi = \infty$. It is evident that $\tau_n$ increases with $n\rightarrow \infty$. Define $\tau_{\infty} = \lim \limits_{n\rightarrow \infty}\tau_{n}$. Then, $\tau_{\infty}\leq\tau_{e}$. It is only necessary to prove that $\tau_{\infty} = \infty$.

By contradiction, it is assumed that $\tau_{\infty} < \infty$. Then, there exist $\varepsilon_{\infty}\in(0, 1)$ and $T_{\infty} > 0$ satisfying $P\{\tau_{\infty}\leq T_{\infty}\} > \varepsilon_{\infty}$. Thus, $\exists n_{\infty}\geq n_{0}$ such that

Define the following:

Clearly, $V(x, y)$ is nonnegative because $s-1-\ln s\geq0$ when $s > 0$. Then,

where

Thus, it can be concluded that

By taking expectations after integrating over the interval $(0, \tau_{n}\wedge T_{\infty})$, we have the following:

Denote $\Omega_{n} = \{w\in\Omega|\tau_{n} = \tau_{n}(w)\leq T_{\infty}\}$. Then, $P(\Omega_{n}) > \varepsilon$ by Eq (4.2). It can be obtained that $x(\tau_{n}, w)$ and $y(\tau_{n}, w)$ equal to either $\frac{1}{n}$ or $n$ for all $w\in\Omega_{n}$. From Eq (4.3), we obtain the following:

which implies

Therefore, we can obtain the following:

□

4.2.   Stationary distribution

Define the following:

where $\delta\in(0, \frac{r}{Ka})$.

Theorem 4. For $\forall (x_{0}, y_{0})\in \mathbb{R}^{2}_{+}$, Model (2.2) possesses a unique, ergodic stationary distribution for $t\geq0$ if $\lambda > 0$.

Proof. Since

its diffusion matrix is given by the following

Let $\theta\triangleq \min\{\sigma_{1}^{2}x^{2}, \sigma_{2}^{2}y^{2}\}$. Then, it has

so Lemma 1 condition $(i)$ holds.

To show Lemma 1 condition (ⅱ) holds, set the following:

where $H_{1}(x) = -M_{1}(\ln x+\ln y-x-y)$, $H_{2}(y) = \frac{1}{2}(cx+y)^{2}$, $M_{1} = \frac{2}{\lambda}\max\{2, \sup\limits_{(x, y)\in R_{+}^{2}}\{-\frac{b_{1}}{2}x^{3}+b_{2}x^{2}-\frac{b_{3}}{2}y^{3}+b_{4}y^{2}\}\} > 0$, and $b_{i} > 0$ (for $i = 1, 2, 3, 4$), will be specified lately. $H(x, y)$ tends to infinity when $(x, y)$ tends to the boundary of $\mathbb{R}^{2}_{+}$, so $H(x, y)$ has a lower bound at $(x', y')\in \mathbb{R}^{2}_{+}$. Define

For $\forall\delta\in(0, \frac{r}{ka})$, there are $\frac{xy}{1+ax^{2}}\leq xy$, $\frac{-cx^{2}}{1+ax^{2}}\leq-2\sqrt{\delta c}x+\delta(1+ax^{2})$ and $\frac{cx^{2 }y}{1+ax^{2}}\leq\frac{cxy}{2\sqrt{a}}$.

By the $\mathrm{It\hat{o}}$ formula, it has the following:

where $\lambda = r-v-d-\frac{\sigma_{1}^{2}}{2} -\frac{\sigma_{2}^{2}}{2}-\delta e^{-(vm+m_{0})}-\frac{\left[\frac{r}{K}+r-v-2e^{-(vm+m_{0})}\sqrt{\delta c}\right]^{2}}{4\left[\frac{r}{K}-\delta ae^{-(vm+m_{0})}\right]}-\frac{(b-d)^{2}}{4b}$.

Similarly,

Combining (4.4) and (4.5), we have the following:

where $b_{1} = \frac{rc^{2}}{K}$, $b_{2} = c^{2}\left(r+\frac{\sigma_{1}^{2}}{2}\right)$, $b_{3} = b$, $b_{4} = \frac{\sigma_{2}^{2}}{2}$, and $b_{5} = cr+M_{1}\left[e^{-(vm+m_{0})}+\frac{ce^{-(vm+m_{0})}}{2\sqrt{a}}\right]$.

Denote $U = [\varepsilon, 1/\varepsilon]\times[\varepsilon, 1/\varepsilon]$ for the given

where $M_{2} = \sup\limits_{(x, y)\in R_{+}^{2}}\left\{-\frac{b_{1}}{2}x^{3}+b_{2}x^{2}-\frac{b_{3}}{2}y^{3} +b_{4}y^{2}+\frac{b_{5}}{2}\left(x^{2}+y^{2}\right)\right\}$. Then, $\mathbb{R}^{2}_{+}\setminus U$ can be divided into four regions, $\mathbb{R}^{2}_{+}\setminus U = \Omega_{1}\bigcup \Omega_{2}\bigcup \Omega_{3}\bigcup \Omega_{4}$, in which

Next, we will prove that $\mathcal{L}W\leq-1$ on each $\Omega_{i}$ ($i = 1, 2, 3, 4$), respectively.

1) For $\forall (x, y)\in \Omega_{1}$, there is $xy\leq\varepsilon y\leq\varepsilon(1+y^{3})$. By (4.6) and (4.7), we obtain the following:

2) For $\forall(x, y)\in \Omega_{2}$, there is $xy\leq\varepsilon x\leq\varepsilon(1+x^{3})$. By (4.6) and (4.7), we obtain the following:

3) For $\forall (x, y)\in \Omega_{3}$, by (4.6)–(4.8), we obtain the following:

4) For $\forall (x, y)\in \Omega_{4}$, by (4.6)–(4.8), we obtain the following:

To sum up, condition (ⅱ) in Lemma 1 holds. Thus the Model (2.2) possesses a unique, ergodic stationary distribution. □

4.3.   Stochastically ultimate boundedness

In this section, we will go into more detail on whether the solution is always bounded.

Theorem 5. Model (2.2)'s solutions are stochastically ultimate bounded.

Proof. Let $V_{1}(x, y) = {x}^{\frac{1}{2}}+{y}^{\frac{1}{2}}.$ Then,

where

where $P_{0}$ represents a positive constant, expressed as

Then,

so

Thus,

and

Therefore, for any normal number $\varepsilon$, set $\chi = \frac{P_{0}^{2}}{\varepsilon^{2}}$. Utilizing Chebyshev's inequality^[31,37], there is the following:

□

4.4.   Extinction and persistence

The objective of this section is to explore the extinction and persistence for the populations.

Lemma 4. For Model (2.2), there are $\lim \limits_{t\rightarrow +\infty}\sup\{t^{-1}\ln x(t)\}\leq 0$ and $\lim \limits_{t\rightarrow +\infty}\sup\{t^{-1}\ln y(t)\}\leq 0$.

Proof. Since

then

Let $M(t) = \sigma_{1} \int_{0}^{t}e^{s} \mathrm{d}B_{1}(s)$. $M(t)$ can be regarded as a localized harness that exhibits the following quadratic variance function:

Take $\mu > 1$ and $\gamma > 1$. For $k\in \mathbb{N}$, denote $\alpha = e^{-\gamma k}$, $\beta = \mu e^{\gamma k}\ln k$, $T = \gamma k$. Then, according to the Exponential Martingale Inequality, it has the following:

Since $\sum \frac{1}{k^{\mu}} < \infty$, then by the Borel-Cantelli lemma^[31,37], there exists $\Omega\in \mathcal{F}$ with $P(\Omega) = 1$ and the integer-valued random variable $k_{0}(w)$ such that for any $\omega\in\Omega$ and $k\geq k_{0}(\omega)$, there is the following:

Hence,

where $Q(x(s)) = r-\frac{rx(s)}{K}-v-\frac{\sigma_{1}^{2}}{2}+\frac{\sigma_{1}^{2}}{2}e^{s-\gamma k}+\ln x(s)$. For any $s$ that satisfies $0\leq s\leq\gamma k$, $x(s) > 0$, a $k$-independent positive constant $\Phi_0$ can be found such that $Q(x(s))\leq \Phi_0$, then, $e^{t}\ln x(t)-\ln x_{0}\leq \Phi_0(e^{t}-1)+\mu e^{\gamma k}\ln k$ and $\ln x(t)\leq \Phi_0(1-e^{-t})+\mu e^{\gamma k-t}\ln k+e^{-t}\ln x_{0}$. Then,

Let $t\rightarrow \infty$; we obtain $\lim \limits_{t\rightarrow \infty} sup\frac{\ln x(t)}{t}\leq0$. In a similar way, one can establish $\lim \limits_{t\rightarrow \infty} sup\frac{\ln y(t)}{t}\leq0$. □

Denote $\eta\triangleq\lim \limits_{t\rightarrow +\infty}\sup x(t)$ and define the following:

Theorem 6. For Model (2.2)'s solution with a given $(x_{0}, y_{0})\in \mathbb{R}_{+}^{2}$, when $\Delta_{1} < 0$, both species $x$ and $y$ will eventually go extinct; when $\Delta_{1} > 0$, species $x$ can keep weakly persistent; when $\Delta_{1} > 0$ and $\Delta_{2} < 0$, species $x$ can keep persistent $\langle x(t)\rangle = \Delta_{1}$ and species $y$ will eventually go extinct; and when $\Delta_{1} > 0$ and $\Delta_{3} > 0$, both species $x$ and $y$ can keep weakly persistent.

Proof. 1) Since

then

Subsequently, it has

then

Therefore, $\lim \limits_{t\rightarrow +\infty}x(t) = 0$.

When $\Delta_{2} < 0$, it has the following:

Then,

Thus,

i.e.,

Therefore, $\lim \limits_{t\rightarrow \infty}y(t) = 0$.

When $\Delta_{2} > 0$, then for $\forall \varepsilon_{1}$ and $\varepsilon_{2} > 0$, $\exists T_{1}$ such that

hold for $t\geq T_{1}$. Thus,

By the arbitrariness of $\varepsilon_{1}$ and $\varepsilon_{2}$ and Lemma 2, there is $\lim \limits_{t\rightarrow \infty}y(t) = 0$.

2) By contradiction. Assume $\Sigma = \{\lim \limits_{t\rightarrow \infty} \sup x(t) = 0 \}$ with $P(\Sigma) > 0$. Thus,

Therefore, $\lim \limits_{t\rightarrow \infty} \sup{t^{-1}\ln x(t)} = r-v-\frac{\sigma_{1}^{2}}{2} = \frac{r}{K}\Delta_{1} > 0$. That means that $\lim \limits_{t\rightarrow \infty} \sup t^{-1}{\ln x(t)} > 0$ holds for any $w\in \Sigma$. There is an inclusion relation $\Sigma\subseteq \{w:\lim \limits_{t\rightarrow \infty} \sup t^{-1}{\ln x(t)} > 0\}$, then $P\{w:\lim \limits_{t\rightarrow \infty} \sup t^{-1}{\ln x(t)} > 0\}\geq P(\Sigma) > 0$, which contradicts to $\lim \limits_{t\rightarrow \infty} \sup t^{-1}{\ln x(t)}\leq0$. Therefore, species $x$ is weakly persistent.

3) It is known that $y(t)$ is extinct as $\lim \limits_{t\rightarrow \infty}y(t) = 0$ when $\Delta_{2} < 0$. From $\lim \limits_{t\rightarrow +\infty}\frac{\ln x_{0}}{t} = 0$, it is known that for $\forall\varepsilon_{3}, \varepsilon_{4} > 0$, $\exists T_{2}$ such that $-\varepsilon_{3}\leq-e^{-(vm+m_{0})}\frac{xy}{1+ax^{2}}\leq\varepsilon_{3}$ and $-\varepsilon_{4}\leq\frac{\ln x_{0}}{t}\leq\varepsilon_{4}$ hold for $t\geq T_{2}$. Obviously,

Based on Lemma 2, it can be deduced that $\langle x(t)\rangle^{*}\leq \frac{K}{r}\left(r-v+\varepsilon_{3}+\varepsilon_{4}-\frac{\sigma_{1}^{2}}{2}\right)$ and $\langle x(t)\rangle_{*}\geq \frac{K}{r}\left(r-v-\varepsilon_{3}-\varepsilon_{4}-\frac{\sigma_{1}^{2}}{2}\right)$. Therefore, $\langle x(t)\rangle = \frac{K}{r}\left(r-v-\frac{\sigma_{1}^{2}}{2}\right) = \Delta_{1} > 0$ (i.e., species $x$ keeps in persistence (in the mean)).

4) Similarly, by using the converse method, one can show that $y$ is also weakly persistent. □

5.   Numerical simulations

This section presents numerical simulations to illustrate the theoretical consequences presented in the study.

5.1.   Verification of the main results

For Model (2.1) with the model parameters $r = 0.4$, $K = 100$, $v = 0.36$, $m = 0.5$, $m_0 = 0.1$, $a = 15$, $c = 0.8$, $b = 0.02$ and $d = 0.01$, there exist three interior equilibria: $E^*_1(0.26, 0.46)$, $E^*_2(2.5, 1.5)$ and $E^*_3(7.5, 1.5)$, where $E^*_1$ is a locally asymptotically stable focus, $E^*_2$ is a saddle and unstable, and $E^*_3$ is a locally asymptotically stable nodel, as illustrated in Figure 2.

For Model (2.2), let consider the following discrete form ^[38]:

where $\varepsilon_{1j}$, $\varepsilon_{2j}$ follow $\mathcal{N}(0, 1)$, $j = 1, 2, \cdots, n$. Take $\sigma_{1} = 0.02$, $\sigma_{2} = 0.01$, $\delta = 0.0001$. For the above model parameters, there is $\lambda = -0.00286 < 0$, and the condition in Theorem 4 does not hold. When $v$ increases (e.g. $v = 0.38$), Model (2.1) has a unique interior equilibrium $E^*(0.193, 0.213)$, which is a globally asymptotically stable focus, as presented in Figure 3(a). In such case, the condition in Theorem 4 does not hold due to $\lambda = -0.00136$. For $a = 8$, Model (2.1) has a globally asymptotically stable focus $E^*(0.155, 0.2)$, as presented in Figure 3(b). In such case, there is $\lambda = 0.000144 > 0$, and the condition in Theorem 4 holds. As depicted in Figure 4, for Model (2.2), there exists an ergodic stationary distribution. It is observed that $(x(t), y(t))$ is stable at $(0.155, 0.2)$ without white noise. $x(t)$ varies around 0.155 and $y(t)$ fluctuates around 0.2 when the white noise is presented.

5.2.   Impact of the noise level and parameters

For Model (2.1) with the model parameters $r = 0.6$, $K = 100$, $a = 2$, $v = 0.3$, $m = 0.5$, $m_{0} = 0.01$, $c = 0.8$, $b = 0.01$ and $d = 0.05$, there is a unique interior equilibrium $E^*(0.1914, 0.208)$, which is locally asymptotically stable.

Different levels of white noise are set to illustrate how white noise affects Model (2.2)'s dynamics. Taking $\sigma_{1} = 0.2$, $\sigma_{2} = 0.1$, it has $\Delta_{1} > 0$, $\Delta_3 > 0$. Then, species $x$ and species $y$ can keep in a weak persistence, as shown in Figure 5(a), (b). This means that the white noise will have little impact on the populations. Let $\sigma_{2}\rightarrow 0.8$; it has $\Delta_{1} > 0$, $\Delta_{2} < 0$. Therefore, species $x$ can keep in persistence (in the mean), while species $y$ will eventually go extinct, as shown in Figure 5(c), (d). It is evident that a larger white noise induces the predator to be eventually extinct. Let $\sigma_{1}\rightarrow 0.8$; it has $\Delta_{1} < 0$. Then, species $x$ and species $y$ will eventually go to extinction, as shown in Figure 5(e), (f). Make it clear that the higher the noise level, the more significant the impact on the species of prey.

The predator's hunting ability is affected by a term $e^{-vm}$, which decreases when $m$ or $v$ increases. For the parameter $m$, we select a range of diverse values and take $\sigma_{1} = 0.2$, $\sigma_{2} = 0.1$. The solution is presented in Figure 6. It is evident that $x(t)$ greatly increases while $y(t)$ increases to be a lesser extent when $m$ rises from 0.5 to 5. There is little change in $x(t)$ when $m$ increases from 5 to 10, while $y(t)$ will change from weakly persistent to eventually extinct. As a result, a modest increase in $m$ benefits $y$, a large increase in $m$ leads to the extinction of $y$, and an increase in $m$ promotes the growth of $x$.

By varying the parameter $v$, the solutions are presented in Figure 7. Obviously, for a smaller $v$, species $x$ and species $y$ are both persistent; species $y$ becomes extinct when $v$ increases, and when the anti-predator level becomes excessively strong, species $x$ and species $y$ will eventually go extinct.

6.   Conclusions

The current work explored a Holling-Ⅲ prey-predator model by incorporating prey habitat selection in an environment with stochastic disturbances. For the model without stochastic disturbances, it was shown that $O(0, 0)$ is constantly unstable, while $E(\overline{K}, 0)$ is stable when $c\leq \overline{c}\triangleq de^{vm+m_{0}}(a+\overline{K}^{-2})$. Additionally, it was shown that Model (2.1) possesses not more than three interior equilibria, and the interior equilibrium $\widehat{E}(\hat{x}, \hat{y})$ is locally asymptotically stable if and only if $\phi(\hat{x}) < \psi'(\hat{x}) < bc\hat{y}/(b\hat{y}+d)$.

Additionally, we investigated the dynamics of Model (2.2). By devising an appropriate Lyapunov function, it obtained the conditions for an ergodic stationary distribution. Furthermore, it was demonstrated that Model (2.2)'s solutions are stochastically, eventually bounded. The results suggest that white noise will eventually cause the extinction of prey and predators when it has a large impact on the prey $x(t)$ ($\Delta_{1} < 0$). When the white noise has more impact on the predators and less disturbance to the prey, the prey can persist (in the mean), and the predators become extinct ($\Delta_{1} > 0$, $\Delta_{2} < 0$). The prey and the predator are weakly persistent when white noise has a small effect on them.

Compared to deterministic models, it is well known that sometimes less intense white noise contributes to the survival of prey; however, white noise with a high intensity is unfavorable to the survival of the population. When the hunting ability of the predator is affected by the $e^{-vm}$ term, an increase in $m$ favors the growth of $x$, a small increase in $m$ would benefit $y$, while a large increase in $m$ would lead to the demise of $y$. The population survival is significantly impacted by anti-predatory behavior. It indicated that $x(t)$ and $y(t)$ can persistently survive when the anti-predator level is small. The population $y(t)$ becomes extinct with an increase of the anti-predator level. The prey and predator eventually go extinct when the anti-predator level becomes excessively strong. Incorporating habitat selection and white noise interference into predator-prey models allows for more accurate estimates of changes in natural populations and a deeper understanding of the mechanisms of species interactions and ecological balance. The results suggest that populations may tend to go extinct when the white noise is high or the habitat selection behavior is strong. Therefore, we can carry out an early warning and take effective intervention and protection measures, such as establishing protected areas, improving the habitat environment, controlling the intensity of the white noise interference, and so on.

Use of AI tools declaration

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

Acknowledgement

The work was supported by the National Natural Science Foundation of China (No. 11401068).

Conflict of interest

The authors declare there is no conflicts of interest.