Analysis of a stochastic Leslie-Gower three-species food chain system with Holling-II functional response and Ornstein-Uhlenbeck process

1.   Introduction

The interactions between predators and their prey have always been and remain the main topics in both ecological and mathematical studies. Several researchers have established mathematical models that depend on the interactions between predators and prey. A famous model of predator-prey interactions was developed by Leslie ^[1,2]. In fact, in most ecosystems, there are many possible food chains that make up a food web, including a number of food chain systems consisting of three species ^[3]. Aziz-Alaoui therefore created a Leslie-Gower system for the tritrophic population to make the original model more realistic ^[4]. To improve the performance of the three-species system, Nindjin and Aziz-Alaoui introduced the Holling II functional response to the model ^[5]. Because biological population systems are inevitably subject to environmental disturbances, we need to introduce random noise to study them. Indeed, many studies have demonstrated that environmental noise causes a continuous variation of some of the key parameters of ecological dynamics, such as growth rates, around some average value ^[6,7,8,9].

To simulate the parameters of a changing environment, there are two common methods. One method is to assume that the parameters can be well represented by a linear function of white noise (see example ^[10]). Gaussian white noise can also simulate complex stochastic dynamical behavior in nonlinear systems (see example ^[11]). Similarly, we can assume that the growth rate or mortality rate of a population is affected by the linear white noise. This can be written as follows:

where $\bar{a}$, which can be obtained by direct calculation, represents the average value of $a$ over the long term. $B(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 ^[12], and $\alpha$ denotes the noise density of $B(t)$.

From a biological standpoint, this linear white noise is somewhat flawed. This flaw has been mildly criticized by Zhang et al. ^[13] and Cai et al. ^[14]. Next, we will explain this unreasonable part from both biological and mathematical perspectives. We assume that for any time interval $[0, t]$, $\langle {a(t)} \rangle$ is the time average of $a(t)$. According to direct calculation, we obtain:

where $\mathbb{N}\left(\cdot, \cdot\right)$ is a one-dimensional Gaussian distribution. It is evident that the average state $\langle {a(t)} \rangle$ on $[0, t]$ has a variance of $\frac{\alpha^2}{t}$, which approaches infinity as $t\to 0^+$. This leads to an unreasonable outcome where the stochastic fluctuation of the growth rate or mortality rate $a(t)$ becomes very large 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 (SDE). We suppose that the variable $\ln a$ is affected by an Ornstein-Uhlenbeck process ^[15,16], and it is described by the following stochastic differential equation:

where $\alpha$ is the speed of reversion and $\sigma$ is the intensity of volatility, respectively. $B(t)$ is the independent standard Brownian motion, which is the same as above. As stated by Mao ^[17], $\ln a(t)$ has a unique explicit solution in the form below:

which shows that the variable $\ln a(t) \sim \mathbb{N}\left(\ln \bar{a}+(\ln a(0)-\ln \bar{a}) e^{-\alpha t}, \sigma^2(1-\right. \; \left.\left.e^{-2 \alpha t}\right) /(2 \alpha)\right)$. It is noted that

In reality, many biological systems may have continuously changing environments due to the interaction of numerous variables. Consequently, mean-reverting processes may be more appropriate for modeling parameters. In most biological systems, both population growth and mortality rates are limited and tend to oscillate around an average value over time. Therefore, we use mean reversion processes to model stochastic perturbations to prevent unrealistic values of growth or mortality rates ^[15].

In this paper, we are interested in a stochastic three-species food chain Leslie-Gower model with the standard white noise and the Holling-II functional response as

All parameters $a_1, a_2, a_3, b_1, v_0, v_1, v_2, v_3, d_0, d_2, d_3$ are positive. Toxin-producing phytoplankton (TPP) has a stabilizing role in aquatic systems and can be used as a biological control mechanism ^[18]. We consider a realistic three-species food chain model consisting of TPP-zooplankton-fish populations. Since the top predator z is assumed to be a sexually reproducing species, the growth of z has two phases: A linear phase and a quadratic phase ^[19]. In this specific case, the TPP population (prey) of size x serves as the only food for the predatory zooplankton population of size y. This zooplankton population, in turn, serves as a favorite food for the generalist vertebrate predator fish population of size z. The predatory population dies out exponentially in the absence of its prey. And the loss to a predatory population is proportional to the reciprocal of the per capita availability of its most favorite food. $a_1$ and $a_3$ are the intrinsic growth rates of prey population TPP and top predator fish populations, respectively, and $a_2$ is the intrinsic mortality rate of zooplankton in the absence of the only food x. $b_1$ measures the strength of intra-specific competition among the individuals of the prey species x(TPP). The parameters $v_0, v_1, v_2$, and $v_3$ are the maximum values that the per capita growth rate can attain. For example, the release of the TPP toxin reduces zooplankton growth and leads to massive zooplankton mortality, which is the biological significance of $v_1$. $d_0$ quantifies the extent to which the environment provides protection to the prey TPP and can be thought of as a refuge or a measure of the effectiveness of the prey in evading a predator's attack. $d_2$ describes the value of zooplankton when the average mortality rate of y reaches $\frac{v_2}{2}$. The mortality of the top predator fish population due to the scarcity of it is favorite food, zooplankton, is measured by $d_3$. The third term on the right-hand side in the first equation of the system (1.4) is obtained by considering the probable effect of the density of the TPP's population on zooplankton's attack rate; the third term on the right-hand side in the second equation of the system (1.4) represents the per capita functional response of the vertebrate predator fish population ^[20]; the second term on the right-hand side in the third equation of the system (1.4) depicts the loss in the top predator (fish population).

Based on the above literature, we want to make system (1.4) more realistic. Therefore, we use the Orenstein-Uhlenbeck process to simulate the stochastic dynamic process. We let

where $B_i(t)$ are three independent Brownian motions, $\alpha_i > 0$ are the speed of reversion, and $\beta_i > 0$ are the intensity of volatility. We let $r_i = \ln a_i$, i = 1, 2, 3, and then we modify system (1.4) accordingly.

We begin by demonstrating that system (1.5) has a unique global solution in Section 3. Subsequently, the ultimate boundedness of the system (1.5) is given in Section 4. The existence of the stationary distribution is then shown in Section 5. Additionally, the extinction of populations is discussed in Section 6. Lastly, in Section 7, we carry out some computer numerical simulations to verify the above theoretical results.

2.   Preliminaries

We define two necessary sets: $\mathbb{G}_n = (-n, n)\times(-n, n)\times(-n, n)$ and $\mathbb{R}_+^n = \{(x_1, \cdots, x_n)\in \mathbb{R}^n|x_k > 0, 0\le k\le n\}$, $||\cdot||$ is the Euclidean norm. Considering a stochastic differential equation in the following form

The existence of a stationary solution to system (1.5) with any initial value is proved by the lemma of Khaminskii ^[21].

Lemma 1. (Khasminskii). Let the vectors $\psi(s, x), \sigma_1(s, x), \ldots, \sigma_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:

Further, there exists a non-negative function $V_*(x)$ such that

where $\mathbb{H}$ is a compact subset defined on $\mathbb{R}^m$. Then the Markov process (2.1) has at least one stationary solution $X(t)$, which has a stationary distribution $\omega(\cdot)$ on $\mathbb{R}^m$.

3.   Existence and uniqueness of the global solution

Theorem 1. For any initial value, $(x_0, y_0, z_0, {r_i}_0)\in\mathbb{R}_+^3\times\mathbb{R}^3$, there exists a unique solution $(x(t), y(t), z(t), r_i(t))$ of the system (1.5) on $t > 0$, where $i = 1, 2, 3$, and it will remain in $\mathbb{R}_+^3\times\mathbb{R}^3$ with probability one.

Proof. It is easy to verify that Eq (1.5) satisfies the linear growth condition and the local Lipschitz condition. So there exists a locally unique solution $(x(t), y(t), z(t), r_i(t))$ defined on $t\in [0, \tau_e)$ (see ^[17]). Where $i = 1, 2, 3$ and $\tau_e$ is the explosion time, we only need to verify $\tau_e = \infty$. We let $n$ be sufficiently large so that $x_0$, $y_0$ and $z_0$ are in the interval $[-n, n]$, defining the stopping time as

where $i = 1, 2, 3$. Obviously, $\tau_n$ is monotonically increasing as $n\to \infty$. Set $\tau_\infty = \mathop {\lim }\limits_{n \to + \infty } {\tau _n}$, obviously, $\tau_\infty \le \tau_e$. So, it is only necessary to prove that $\tau_\infty = \infty$. If not, then there exists a constant $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, z, r_i)$ is constructed as follows:

Applying Ito's formula, we have

where

After the routine calculation, it implies:

Integrating both sides of inequality (3.3) from $0$ to ${\tau _n} \wedge T$ and taking expectations, we get

When $\tau_n\le T$, let $\Omega_n = {\tau_n\le T}$. Then $\mathbb{P}(\Omega_n) > \varepsilon$ holds inequality (3.1). Notice that for any $\omega\in\Omega_n$, there exists $i$ such that $x(\tau_n, \omega), y(\tau_n, \omega), z(\tau_n, \omega) = -n \text{ or } n$.

According to Eq (3.4), it can be deduced that

as $n \to \infty$, we have

and Theorem 1 is proved. □

4.   Ultimate boundness

Ecosystems can only support a limited number of individuals due to finite resources. When a population exceeds this limit, the environment cannot provide enough sources for its survival. This means that population density has an upper bound and will stabilize over time. To show this theoretically, we need to define stochastically ultimate boundedness ^[22] and apply it to the system (1.5).

Definition 1. System (1.5) 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, z_0, {r_i}_0)$, where $i = 1, 2, 3$, and the solution of system (1.5) has the property that

Lemma 2. 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, z_0, {r_i}_0)$, where $i = 1, 2, 3$, such that the solution of system (1.5) has the property that

Proof. Define the Lyapunov function $W_0:\mathbb{R}_+^3\times\mathbb{R}^3\to\mathbb{R}$

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

where $\lambda = \theta \min\{\alpha_1, \alpha_2, \alpha_3\}$. Note that

where

Thus, we can obtain

Then we have

where $i = 1, 2, 3$, and then the result (4.2) holds by setting $M(\theta) = {3^{\frac{\theta }{2}}}\theta\frac{\kappa(\theta)}{\lambda}$. □

Theorem 2. The equation of system (1.5) is stochastically ultimately bounded.

Proof. By Lemma 2, there exists $M$ such that

Now, for any $\epsilon > 0$, let $\chi = \sqrt{3}\frac{\kappa(0.5)^2}{4\epsilon^2\lambda^2}$. Then, using Chebyshev's inequality, we can conclude that

Hence,

□

5.   Existence a stationary distribution

In biology, one of the primary goals is to comprehend the long-term behavior of systems. In this part of the paper, we demonstrate that system (1.5) possesses a stationary distribution that can be used to make long-term predictions for the system under the influence of stochastic perturbations. We also derive the sufficient conditions that guarantee the existence of such a stationary distribution for the system (1.5).

Definition 2. Define $M_0$ to be a natural number that satisfies the following conditions

where

where $i = 1, 2, 3$.

Theorem 3. For any initial value $(x(0), y(0), z(0), {r_i}(0))$, where $i = 1, 2, 3$, the system (1.5) has a stationary distribution with the definition of $M_0$.

Proof. We divide the relevant proof into the following two steps.

Step 1. Define the Lyapunouv Function $V_1:\mathbb{R}_+^3\times\mathbb{R}^3\to\mathbb{R}$

Applying Ito's formula to the definition of $M_0$, we obtain

Noting that the function $V_1$ tends to $\infty$ as $x$, $y$, or $z$ approaches the boundary of $\mathbb{R}_+$ or as $||(x, y, z, r_i)||\to\infty$. Thus, there exists a point $(x^0, y^0, z^0, r_i^0)$ in the interior of $\mathbb{R}_+^3\times\mathbb{R}^3$, at which the value of the function will be minimized. A non-negative function $V_2(x, y, z, r_i)$ can be constructed as

According to Ito's formula, we have

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

and we define

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

We verify $\mathcal{L}V_2(x, y, z, r_i)\le -1$ for any $(x, y, z, r_i)\in \left(\mathbb{R}_{+}^3 \times \mathbb{R}^3\right) \backslash \mathbb{H}_\varepsilon$. Noting that $\left(\mathbb{R}_{+}^3 \times \mathbb{R}^3\right) \backslash \mathbb{H}_\varepsilon = \bigcup_{k = 1}^{9} \mathbb{H}_{k, \varepsilon}^c$, where

Case 1. If $(x, y, z, r_i)$ is located in the set defined by $\mathbb{H}_{1, \varepsilon}^c$, then one can obtain the corresponding results by combining Eqs (5.2) and (5.3).

Case 2. If $(x, y, z, r_i)$ is located in the set defined by $\mathbb{H}_{2, \varepsilon}^c$ and $\mathbb{H}_{3, \varepsilon}^c$, consequently, from (5.2) and (5.3), we can obtain the relevant result.

Case 3. For any $(x, y, z, r_i)$, which is located in the set defined in $\mathbb{H}_{j, \varepsilon}^c$, according to Eqs (5.2) and (5.3), we obtain

Case 4. If $(x, y, z, r_i) \; \in \; \mathbb{H}_{7, \varepsilon}^c$, it follows from (5.2) and (5.4) that

Case 5. If $(x, y, z, r_i)$ lie within the set demarcated by $\mathbb{H}_{8, \varepsilon}^c$, the relevant conclusion is deduced through (5.2) and (5.5).

Case 6. In the event that $(x, y, z, r_i)$ is situated within the set defined by $\mathbb{H}_{9, \varepsilon}^c$, the associated findings can be calculated by (5.2) and (5.6).

Defining

In summary, we can obtain that there exists a sufficiently small constant $\varepsilon$ such that $L V_2\left(x, y, z, r_i\right) \leq-1$ for any $\left(x, y, z, r_i\right) \in\left(\mathbb{R}_{+}^3 \times \mathbb{R}^3\right) \backslash \mathbb{H}_\varepsilon$, where $\varepsilon$ satisfies that

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

□

6.   Extinction

Theorem 4. We define

When $\varphi_1 < 0, \varphi_2 < 0, \varphi_3 < 0$, then $x(t), y(t), z(t)$ are extinct.

Proof. From Eq (1.1) and the definition of the Ornstein-Uhlenbeck process, we have

Equation (6.1) clearly indicates that $a(t)$ follows the Gaussian distribution $\mathbb{N}\left(E\left[a_i(t)\right], \operatorname{VAR}\left[a_i(t)\right]\right)$ on $[0, \mathrm{t}]$. It is easy to derive that $E\left[a_i(t)\right] = \bar{a}_i+\left[a_i(0)-\bar{a}_i\right] e^{-\alpha_i t}$ and $\operatorname{VAR}\left[a_i(t)\right] = \frac{\beta_i^2}{2 \alpha_i}\left(1-e^{-2 \alpha_i t}\right)$. Therefore, we find that the term $\beta_i \int_0^t e^{-\alpha_i(t-s)} d B_i(s)$ follows the normal distribution $\mathbb{N}\left(0, \frac{\beta_i^2}{2 \alpha_i}\left(1-e^{-2 \alpha_i t}\right)\right)$. Then, it is equivalent to $\frac{\beta_i}{\sqrt{2 \alpha_i}} \sqrt{1-e^{-2 \alpha_i t}} \frac{\mathrm{dB}_i(t)}{\mathrm{d} t}$. Let $\sigma_i(t) = \frac{\beta_i}{\sqrt{2 \alpha_i}} \sqrt{1-e^{-2 \alpha_i t}}$, where $B_i(t)$ stands for a standard Brownian motion. Thus, (6.1) can be almost surely written as follows: ^[23]

And then we modify system (1.5) accordingly.

Applying Itô formula to $\ln x(t), \ln y(t), \ln z(t)$, we can get

Integrating from 0 to $t$, we have

From Eq (6.4), we obtain

Then, according to the strong law of large numbers and the definition of the Ornstein-Uhlenbeck process, if $\bar{\varphi}_1+\varepsilon_1 < 0$, we have $\lim _{t \rightarrow +\infty} x(t) = 0$. Similarly, when $\bar{\varphi}_2 < 0$ and $\bar{\varphi}_3 < 0, \lim _{t \rightarrow +\infty} y(t) = 0$ and $\lim _{t \rightarrow +\infty} z(t) = 0$. Theorem 4 is proved. □

So far, we have proved that there is a unique global solution for the system (1.5), the solution to the system (1.5) is stochastically ultimately bounded, the system (1.5) has a stationary distribution, and the system (1.5) has an extinction. Next, we will perform some computer simulations to verify our conclusions.

7.   Numerical simulations

In this section, we will use numerical simulations to verify our conclusions. Using the Milstein higher order method, we obtain the discretization equation for system (1.5). The corresponding discretization equation is as follows:

where $\Delta t > 0$ denotes the time increment, and $\eta_j$, $\xi_j$, and $\zeta_j$ are three independent stochastic variables that follow the standard Gaussian distribution $\mathbb{N}(0, 1)$. Besides, $(x^j, y^j, z^j, a_i^j)$ is the corresponding value of the $jth$ iteration of the discretization Eq (7.1), where $i = 1, 2, 3$ and $j = 1, 2, \cdots$. We will use some different combinations of biological parameters in Table 1.

Based on the real datasets ^{[4,18,24,25,26,27]}, we use the TPP-zooplankton-fish population food chain system as an example. After determining the range of values for each parameter of the system (1.5) (e.g., the values of the intensity of volatility $\beta_i$(i = 1, 2, 3) belong to [0, 5] ^[13,14,28]), we fixed some of the parameter values by making reasonable assumptions in combination in Table 2.

Example 1. In Table 1, we choose the combination $\mathcal{A}_1$-$\mathcal{A}_3$ as the value of biological parameters of the system (1.5). Clearly, it follows from Theorem 1 that the system (1.5) is stochastically ultimately bounded. By choosing the total number of iterations, $T_{max} = 2000$. Then, we obtain Figure 1.

Remark 1. It can be seen from Figure 1 that the growth rates are disturbed around the given mean value, which reflects the characteristics of the mean regression of the Ornstein-Uhlenbeck process. At the same time, it can be seen that different coefficient combinations have different solutions, and the solutions are all existing and unique. Thus, the conclusion of Theorem 1 can be verified.

Example 2. In Table 1, we choose the combination $\mathcal{A}_1$-$\mathcal{A}_3$ as the value of the biological parameters of the system (1.5). Obviously, it follows from Lemma 2 that the $\theta th$ order moments of the solutions of the system (1.5) are bounded, and the solutions of the system (1.5) are ultimately bounded. By choosing the total number of iterations, $T_{max} = 2000$. Then, we obtain Figure 2.

Remark 2. It can be seen from Figure 2 that the expected value of the above three coefficient combinations $(A_1, A_2, \, and \, A_3)$ is less than an upper bound. With the increase in time t, the probability P is gradually stable and greater than a constant, which means $\limsup_{{\rm{t}} \to \infty }P(|(x, y, z)| > \chi)\le\epsilon$. which verifies the conclusion of Theorem 2.

Example 3. Based on Table 1, we select the $\mathcal{A}_3$ combination as the biological parameters of the system (1.5). By Theorem 3, we know that the distribution of the system (1.5) exists. We set the number of iterations to $T = 2000$. From Figure 3, we can see the system (1.5) has a stationary distribution.

Example 4. In Table 1, we choose the combination $\mathcal{A}_4$–$\mathcal{A}_6$ as the value of the biological parameters of the system (1.5). As proven in Theorem 4, the system (1.5) possesses extinction. By choosing the total number of iterations, $T_{max} = 2000$. Then, we obtain Figure 4.

8.   Conclusions

This paper focuses on the mathematical properties of the three-species food chain stochastic Leslie-Gower model with the Ornstein-Uhlenbeck process. On the basis of previous research on food web models, we consider that the real dynamic environment may be affected by random factors such as seasonal changes and natural disasters, and then random noise is introduced. Compared to using standard white noise methods, simulations with Ornstein-Uhlenbeck processes are more realistic and reliable. Therefore, we demonstrate the feasibility of studying a more complex Leslie-Gower model with the Ornstein-Uhlenbeck process. The main contributions of this paper are as follows:

We investigate the good properties of stochastic mathematical models with the Ornstein-Uhlenbeck process, providing a new approach to studying the stochastic properties of the Leslie-Gower systems. We derive the conditions for the existence and uniqueness of the global solution under stochastic disturbance, as well as the conditions for the stochastic ultimate boundedness, the existence of a stationary distribution, and the extinction.

However, there are still many open questions that warrant further research. In this paper, we only study the effect of stochastic noise driven by the Ornstein-Uhlenbeck process on the three-species food chain Leslie-Gower model, while the study of Markov transformation and Lévy jump is still scarce, which is one of the future directions. Most existing studies assume that the system parameters are constant or time-invariant, or they rarely consider the case of periodic coefficients. However, in reality, the ecological environment that supports populations often exhibits cyclical patterns. Consequently, it would be more realistic to study stochastic population systems with periodic coefficients.

Author contributions

Ruyue Hu and Xiaohui Ai: Conceptualization, Writing-original draft; Yifan Wu: Software, Visualization, Writing-original draft; Chi Han: Software, Visualization. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

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

Acknowledgments

Thanks to all anonymous reviewers for their valuable comments and helpful suggestions that made the article better.

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(202210225154).

Conflict of interest

The authors declare there is no conflicts of interest.