Dynamical behavior of a stochastic predator-prey model with Holling-type III functional response and infectious predator

1.   Introduction

In ecological systems, the interaction between predator and prey is the most important to maintain ecosystem balance. The predator-prey model plays an important role in studying the relationship between two populations. Since Lotka^[1] and Volterra^[2] proposed the classical Lotka-Volterra predator-prey model, the predator-prey models have attracted much attention. Among them, many authors ^[3,4,5,6] proposed the predator-prey model with epidemic and considered the impact of epidemic on population.

As a matter of fact, the real ecosystem is inevitably affected by environmental noise. From the biological and mathematical point of view, the stochastic predator-prey model can predict future dynamics more accurately than the deterministic model. Therefore, many stochastic predator-prey models have been proposed. For example, Shi et al.^[7] considered a stochastic Holling-Type II predator-prey model with stage structure and refuge for prey. Liu et al.^[8] studied dynamics of stochastic predator-prey models with distributed delay and stage structure for prey. Ma et al.^[9] proposed a stochastic one-predator-two-prey time-delay system with jumps. Liu and Jiang^[10] studied the influence of the fear factor on the dynamics of a stochastic predator-prey model. Qi and Meng ^[11] proposed a stochastic predator-prey system with prey refuge and fear effect. These pieces of literature consider the influence of linear disturbance on the model.

Inspired by the above literatures, in this paper, we consider the epidemic disease and nonlinear perturbations into the model to accurately predict the future dynamics. The prey population is denoted by $X(t)$ at time $t$. In the presence of disease, the predator population $N(t)$ is divided into two classes, namely the susceptible predator $Y_{S}(t)$ and the infected predator $Y_{I}(t)$ at time $t$. The random perturbation may be dependent on the square of the state variables $X(t)$, $Y_{S}(t)$, and $Y_{I}(t)$, respectively. The Holling type III response functions $\frac{cX^{2}Y_{S}}{d+X^{2}}$ and $\frac{nX^{2}Y_{I}}{d+X^{2}}$ represent the functional response of the predator to the prey. Therefore, we propose the following stochastic predator-prey model:

where $r$ and $K$ are respectively the intrinsic growth rate and the environmental carrying capacity for prey. $d$ is the half-saturation constant, $\mu$ relates to the predators natural mortality rate, $\nu$ is the death rate of the predator due to disease. $c, n$ is the maximum value which per capita reduction rate can $X(t)$ attain. $h, k$ has a similar meaning to $c, n$. $\beta$ is a disease standard incidence disease-induced mortality rate of infected predators. $B_{i}(t)(i = 1, 2, 3)$ are mutually independent standard Brownian motions, and $\sigma_{il}(i = 1, 2, 3;l = 1, 2)$ are nonnegative and referred as their intensities of stochastic noises which are used to describe the volatility of perturbation.

This paper is organized as follows. In Section 2, we investigate the existence and uniqueness of the global positive solution for the stochastic predator-prey model. In Section 3, we establish sufficient conditions for the existence of an ergodic stationary distribution of the positive solutions to the model. In Section 4, we prove the extinction of the predator and prey populations under certain parametric restrictions. In Section 5, we give a summary of the main results and a series of numerical simulations to illustrate the theoretical results. Finally, concluding remarks are presented in Section 6.

2.   Existence and uniqueness of the global positive solution

In order to prove the existence and uniqueness of the solution, we first introduce some preliminaries that will be used in the rest of the paper.

We set $(\Omega, \mathcal{F}, \{\mathcal{F}_{t}\}_{t\geq 0}, \mathbb{P})$ to stand for a complete probability space with a filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets).

Generally speaking, an $n$-dimensional stochastic differential equation is given by:

where $F(t, X)$ represents a function in $[0, +\infty)\times\mathbb{R}^{n}$, and $G(t, X)$ is a $n\times m$ matrix, $F(t, X)$ and $G(t, X)$ satisfy the locally Lipschitz conditions in $X$. $B(t)$ is $m$-dimensional standard Brownian motion defined on the complete probability space. $C^{2, 1}(\mathbb{R}^{n}\times[0, +\infty), \mathbb{R}_{+})$ is a family of all nonnegative functions $V(x, t)$ which are defined on $\mathbb{R}^{n}\times[0, +\infty)$, such that this family of functions are continuously twice differentiable on $X$ and continuously once differentiable on $t$. The differential operator $\mathcal{L}$ for the stochastic differential Eq $(2.1)$

Applying $\mathcal{L}$ to a function $V(X, t)\in C^{2, 1}(\mathbb{R}^{n}\times[0, +\infty), \mathbb{R}_{+})$, we get

where

by Itô formula, when $X(t)\in \mathbb{R}^{n}$, we have

Next, by using the Lyapunov function method^[12], we shall show that the system $(1.1)$ has a unique local positive solution, then we show that this solution is global. And the main results are as follows.

Theorem 1. For any initial value $(X(0), Y_{S}(0), Y_{I}(0))\in\mathbb{R}_{+}^{3}$, then system $(1.1)$ has a unique positive solution $(X(t), Y_{S}(t), Y_{I}(t))$ for all $t\geq0$ almost surely, and the solution remains in $\mathbb{R}_{+}^{3}$ with probability 1.

Proof. It is obvious that the coefficients of the system $(1.1)$ satisfy the local Lipschitz condition, then for any given initial value $(X(0), Y_{S}(0), Y_{I}(0))\in \mathbb{R}_{+}^{3}$, there is a unique local solution $(X(t), Y_{S}(t), Y_{I}(t))$ for $t\in[0, \tau_{e})$, where $\tau_{e}$ is the explosion time ^[13]. Then $\tau_{e} = +\infty$ demonstrates that the solution of the system $(1.1)$ is global. At first, we prove that $X(t), Y_{S}(t)$ and $Y_{I}(t)$ do not explode to infinity at a finite time. Let $k_{0}\geq1$ be sufficiently large constant so that $X(0), Y_{S}(0)$ and $Y_{I}(0)$ lie within the interval $[\frac{1}{k_{0}}, k_{0}]$. For each integer $k\geq k_{0}$, we define the stopping time as follows:

It is easy to see that $\tau_{k}$ is increasing as $k\rightarrow +\infty$. We set $\tau_{\infty} = {\lim_{k\rightarrow +\infty}\tau_{k}}$, whence, $\tau_{\infty}\leq\tau_{e}$. If we can show that $\tau_{\infty} = +\infty$ a.s., then we can obtain that $\tau_{e} = +\infty$ a.s., and $(X(t), Y_{S}(t), Y_{I}(t))\in \mathbb{R}_{+}^{3}$ a.s.. If this statement is false, then there exists a pair of constants $T > 0$ and $\varepsilon\in(0, 1)$ such that $P\{\tau_{\infty}\leq T\} > \varepsilon$. Thus there exists an integer $k_{1}\geq k_{0}$ such that

Consider the $C^{2}$-function $V:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$ as follows:

where $0 < \alpha < 1$. Applying Itô formula leads to

where

where $K_{0}$ is a positive constant.

Therefore, we can have

Integrating $(2.6)$ from $0$ to $T\wedge\tau_{k}$, set $V(T\wedge\tau_{k}) = V(X(T\wedge\tau_{k}), Y_{S}(T\wedge\tau_{k}), Y_{I}(T\wedge\tau_{k}))$, and taking expectation on both sides yields

Set $\Omega_{k} = \{\tau_{k}\leq T\}$. By $(2.4)$, we have $P(\Omega_{k})\geq\varepsilon$ for $k\geq k_{1}$, we obtain

where $1\Omega_{k}$ is the indicator function of $\Omega_{k}$. For every $\omega\in\Omega_{k}$, there exists at least one of $X(\tau_{k}, \omega), Y_{S}(\tau_{k}, \omega)$ and $Y_{I}(\tau_{k}, \omega)$, which equals either $k$ or $\frac{1}{k}$. Thus we get

where

Combining $(2.7)$–$(2.9)$, we can have

When $k\rightarrow +\infty$, we obtain

which is a contradiction. Thus, we must have $\tau_{e} = +\infty$ a.s. Consequently, $X(t)$, $Y_{S}(t)$, and $Y_{I}(t)$ are positive and global. Then the proof is complete.

3.   Stationary distribution

In this section, we give a sufficient condition for the existence of a stationary distribution of the positive solution of the system $(1.1)$.

Let $X(t)$ be a regular time-homogeneous Markov process in $\mathbb{R}^{d}$, and $X(t)$ is described by the following stochastic differential equation

The diffusion matrix of the process $X(t)$ is defined as follows

Lemma 1. ^[15] If there exists a bounded open domain $D\subset E_{d}$ with regular boundary $\Gamma$, having the following properties

(i) There is a normal number, such that$\sum\limits_{i, j = 1}^{d}a_{ij}(x)\xi_{i}\xi_{j}\geq M|\xi|^{2}, x\in D, \xi\in \mathbb{R}^{d};$

(ii) There exists a non-negative $C^{2}$ function $V$ such that $\mathcal{L}V$ is negative to any $x\in E_{d}\backslash D$;

then the Markov process $X(t)$ has a unique ergodic stationary distribution $\pi(\cdot)$.

Theorem 2. Let $(X(t), Y_{S}(t), Y_{I}(t))$ be a solution of the model $(1.1)$ for any given initial value $(X(0), Y_{S}(0), Y_{I}(0))\in\mathbb{R}_{+}^{3}$, if

such that

then the system $(1.1)$ has an ergodic stationary distribution $\pi(\cdot)$.

Proof. In order to prove Theorem 2, it suffices to verify conditions (i) and (ii) in Lemma 1 hold. Now we verify the condition (ii). For convenience defining the notations

Define $C^{2}$ functions $V_{1}:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$:

By Itô formula to $V_{1}(X, Y_{S}, Y_{I})$ and system $(1.1)$, we obtain

where

then

According to the expressions of $\rho_{1}$ and $\rho_{2}$, we can obtain that $H(K) = 0$, $H'(K) = 0$ and $H''(X) < 0$, thus, for all $X\in(0, +\infty)$, $H(X)\leq H(K) = 0$. Then

where

Define $C^{2}$ function $V_{2}:\mathbb{R}_{+}^{2}\rightarrow \mathbb{R}$,

By Itô formula to $V_{2}(Y_{S}, Y_{I})$ and system $(1.1)$, we get

We set $f(x) = -\frac{\eta(1-\theta)}{4}x^{\theta+2}+x^{2}, x\in[0, +\infty)$, then $f'(x) = -\frac{\eta(1-\theta)}{4}(\theta+2)x^{\theta+1}+2x$, as $x_{0} = \Big[\frac{8}{\eta(1-\theta)(\theta+2)}\Big]^{\frac{1}{\theta}}$, $f'(x_{0}) = 0$, and $f''(x_{0}) = -\frac{\eta(1-\theta)}{4}(1+\theta)(\theta+2)x_{0}^{\theta}+2 = -2\theta < 0$. Thus $f(x)\leq f(x_{0}) = \Big[\frac{8}{\eta(1-\theta)(\theta+2)}\Big]^{\frac{2}{\theta}}\frac{\theta}{\theta+2}$. Combining $(3.2)$ and $(3.3)$, we can have

Define $C^{2}$ functions $V_{3}:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$:

where $0 < \theta < 1$. By Itô formula to $V_{3}(X, Y_{S}, Y_{I})$ and system $(1.1)$, we obtain

where

Define a $C^{2}$ function $Q:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$

where $M > 0$ is sufficiently large, such that

Furthermore, $Q(X, Y_{S}, Y_{I})$ is continuous, and $(X^{*}, Y_{S}^{*}, Y_{I}^{*})$ is a minimum value point of $Q(X, Y_{S}, Y_{I})$ in $\mathbb{R}_{+}^{3}$. Therefore, a non-negative $C^{2}$ function is defined as follows

by Itô formula and combining $(3.4)$ and $(3.5)$, we get,

Structure compact set

where $\varepsilon_{i}(0 < \varepsilon_{i} < 1, i = 1, 2, 3)$ is a sufficiently small constant and satisfying

For the sake of discussion, we're going to break $\mathbb{R}_{+}^{3}\backslash D$ down into six areas:

Next, we will prove that $\mathcal{L}V(X, Y_{S}, Y_{I})\leq-1$ for any $(X, Y_{S}, Y_{I})\in D^{c} = \mathbb{R}_{+}^{3}\backslash D$, we can discuss under six cases:

Case 1. For any $(X, Y_{S}, Y_{I})\in D_{1}$, by $(3.6)$–$(3.8)$, we obtain

Case 2. For any $(X, Y_{S}, Y_{I})\in D_{2}$, by $(3.6)$, $(3.7)$ and $(3.9)$, we have

Case 3. For any $(X, Y_{S}, Y_{I})\in D_{3}$, by $(3.6)$, $(3.7)$ and $(3.10)$, we get

Set

We get

Case 4. For any $(X, Y_{S}, Y_{I})\in D_{4}$, by $(3.11)$ and $(3.14)$, we have

Case 5. For any $(X, Y_{S}, Y_{I})\in D_{5}$, by $(3.12)$ and $(3.14)$, we obtain

Case 6. For any $(X, Y_{S}, Y_{I})\in D_{6}$, by $(3.13)$ and $(3.14)$, we have

For sufficiently small positive numbers $\varepsilon_{i}(i = 1, 2, 3)$, we obtain

In order to verify conditions (i) in Lemma 1 hold we set

then, for all $(X, Y_{S}, Y_{I})\in D$, the diffusion matrix of model $(1.1)$ is given as follows:

where $\xi = (\xi_{1}, \xi_{2}, \xi_{3})\in \mathbb{R}^{3}$.

From the Lemma 1, the conclusion of the Theorem 2 holds. We complete the proof of Theorem 2.

4.   Extinction

In this section, we shall prove the extinction of the predator and the prey populations under certain assumptions.

Theorem 3. For any given initial value $(X(0), Y_{S}(0), Y_{I}(0))\in\mathbb{R}_{+}^{3}$, if $r-\frac{\sigma_{11}^{2}}{2} < 0$, then

That is to say that $(X(t), Y_{S}(t), Y_{I}(t))$ exponentially converges to $(0, 0, 0)$ a.s.

Proof. By Itô formula, and system $(1.1)$, we obtain

Integral $(4.1)$ from $0$ to $t$, and divide by $t$, we get

Using the strong law of large numbers and $(4.2)$, and taking the upper bound and the limit, we have

the upper formula indicates that

Therefore, for any $\varepsilon > 0$, there are a constant $T$ and a set $\Omega_{\varepsilon}\subset\Omega$ satisfying $\mathbb{P}(\Omega_{\varepsilon}) > 1-\varepsilon$, and $\frac{X^{2}}{d+X^{2}}\leq\frac{\varepsilon^{2}}{d+\varepsilon^{2}}$ for $t\geq T$ and $\omega\in\Omega_{\varepsilon}$. Therefore, by Itô formula, and system $(1.1)$, we obtain

Taking the superior limit on both side of $(4.5)$ and noting

Letting $\varepsilon\rightarrow0$ leads to

which implies that

We complete the proof of Theorem 3.

5.   Numerical simulations

To conform the analytical results above, we use Milsteins higher order method ^[16,17] to find the strong solutions of system $(1.1)$. The discrete equations of system $(1.1)$ are described by

where $\xi_{1j}, \xi_{2j}$, and $\xi_{3j}(j = 1, 2, \cdots)$ are independent Gaussian random variables $N(0, 1)$, $\sigma_{il}(i = 1, 2, 3;l = 1, 2)$ are intensities of white noises.

The parameter values are chosen as follows: $r = 0.5, K = \frac{5}{3}, \beta = 0.128, c = 0.52, n = 0.4, h = 0.32, k = 0.13, d = 0.5, \mu = 0.1, \delta = 0.046$, and the step size $\Delta t = 0.01.$

We give some numerical simulations to illustrate our theoretical results of Theorem 2. Figure 1 is the model of stochastic system $(1.1)$ with $\sigma_{11} = 0.02, \sigma_{12} = 0.01$, $\sigma_{21} = 0.02, \sigma_{22} = 0.045$, $\sigma_{31} = 0.02, \sigma_{32} = 0.3$, $\rho_{1} = 1.9404$. By computation, we get that $r+K\sigma_{11}\sigma_{12}-K^{2}\sigma_{12}^{2} = 0.5001 > 0$, $\rho_{1} = 1.9404 > \frac{1}{r+K\sigma_{11}\sigma_{12}-K^{2}\sigma_{12}^{2}}\Big[\frac{2d(h+k)K}{(d+K^{2})^{2}}+\frac{h+k}{d}\Big] = 1.9394$, $\rho = \frac{(h+k)K^{2}}{d+K^{2}}-\rho_{1}K^{3}\sigma_{12}^{2}-\rho_{1}K^{2}\sigma_{11}\sigma_{12}-\frac{\rho_{1}K}{2}\sigma_{11}^{2} -2\mu-\beta-\delta-\sigma_{21}^{2}-\sigma_{31}^{2} = 0.0041 > 0$. The conditions of the Theorem 2 are satisfied. According to Theorem 2, Figure 1 shows that there is a unique ergodic stationary distribution of the model $(1.1)$.

Figure 2 shows the stochastic system $(1.1)$ with $\sigma_{11} = 1.02, \sigma_{12} = 0.01$, $\sigma_{21} = 0.02, \sigma_{22} = 0.045$, $\sigma_{31} = 0.02, \sigma_{32} = 0.3$, we get that $\sigma_{11}^{2} = 1.0404 > 1 = 2r$, which satisfies the conditions of Theorem 3, this is consistent with our conclusion in Theorem 3, when the intensities of white noises are sufficiently large, the populations of the predator and the prey of the stochastic system $(1.1)$ are extinct.

6.   Conclusions

In this paper, we construct and analyze a stochastic predator-prey model with Holling-type III functional response and infectious predator. The effect of stochastic perturbations on the ergodic stationary distribution and the possible extinction of the predator and the prey have been studied in detail. By establishing a suitable Lyapunov function, the existence and uniqueness of the global positive solution of the system are proved. Then, by establishing a series of suitable Lyapunov functions, we established the conditions of the stationary distribution and the ergodic to the system $(1.1)$. In addition, we derived sufficient criteria for the extinction of the predator and the prey populations. That is, if the environment disturbance is big enough, then the entire predator-prey system $(1.1)$ can die out. Our results show that predator and prey populations have the ability to adapt to external environmental disturbances. If the intensity of environmental perturbation is small enough such that the conditions in Theorem 2 are satisfied, then the predator-prey system is permanent in a sense. Finally, by using Milstein's scheme, we carry out a series of numerical simulations to illustrate our results.

In the model, the combination of epidemic and white noise (Brownian motion) has a great impact on the dynamics, complexity and extinction of the predator and the prey populations. The existence of the ergodic stationary distribution of the positive solutions to the proposed model is a very important problem for the predator-prey system and affects the survival of the species in the environment. Other exciting research points deserve more investigation if we consider time-delays such as ^[18], which will be considered in future work.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No.11971127) and the Key Research Project of Higher school in Henan Province (No.20B110017; No.22B110006).

Conflict of Interest

The authors declare that they have no competing of interests regarding the publication of this paper.