On the Dynamics of a Stochastic Ratio-Dependent Predator-Prey System with Infection for the Prey

Jiying Ma¹, Haimiao Ren²

¹College of Science, University of Shanghai for Science and Technology, Shanghai, China.
²Department of Mathematics, University of Shanghai for Science and Technology, Shanghai, China.

DOI: 10.4236/ojapps.2021.104032                

Abstract

In this paper, we investigate the dynamics of a stochastic predator-prey model with ratio-dependent functional response and disease in the prey. Firstly, we prove the existence and uniqueness of the positive solution for the stochastic model by using conventional methods. Then we obtain the threshold  for the infected prey population, that is, the disease will tend to extinction if < 1, and it will exist in the long time if  > 1. Finally, the sufficient condition on the existence of a unique ergodic stationary distribution is obtained, which indicates that all the populations are permanent in the time mean sense. Numerical simulations are conducted to verify our analysis results.

Keywords

Stochastic Predator-Prey Model, Ratio-Dependent, Stationary Distribution, Extinction

1. Introduction

The research of eco-epidemiology involving ecological and epidemiological models is a significant field in mathematical biology. To our knowledge, Anderson and May were the first to study the spread and persistence of infectious diseases by formulating an eco-epidemiological prey-predator model [1]. Recently, a large number of researchers have devoted to the study of eco-epidemiological models (see [2] - [7] ). For example, Chakraborty et al. [2] have investigated the positivity and boundedness of the solutions for a predator-prey model with disease in prey population. Mondala et al. studied the local and global dynamical behavior of a predator-prey eco-epidemiological model with disease in predator [5].

In predator-prey system, functional response plays an important role in the population dynamics. Holling types functional response functions, namely Holling types I, II, III and IV, have been extensively used and investigated [8] [9]. In recent decades, Beddington-DeAngelis and Crowley-Martin type functional response are also widely chosen to model the predation [10] [11]. Li et al. [10] analyzed a stochastic predator-prey model with disease in the predator and Beddington-DeAngelis functional response. They showed that the stochastic system has a similar property to the corresponding deterministic system when the white noise is small enough. In many cases where the predator has to seek for the prey, the per capita predator growth rate should be a function of the ratio of prey to predator abundance in predator-prey model. Thus, the predator-prey models with ratio-dependent functional responses have been proposed and mathematically studied [12] [13] [14] [15]. Based on the literatures, we propose an eco-epidemiological model with infection in the prey and ratio-dependent functional responses as follows

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=rS\left(1-\frac{S}{k}\right)-\frac{\alpha SP}{mP+S+I}-bSI,\\ \frac{\text{d}I}{\text{d}t}=bSI-d_{1}I-\frac{\beta IP}{mP+S+I},\\ \frac{\text{d}P}{\text{d}t}=\frac{c\alpha SP}{mP+S+I}+\frac{c\beta IP}{mP+S+I}-d_{2}P,\end{array}$ (1)

where $S\left(t\right)$, $I\left(t\right)$ and $P\left(t\right)$ denote the densities of the susceptible prey, infected prey and predator respectively. Here, the susceptible prey is subject to the logistic growth, r is the intrinsic growth rate, and $\frac{r}{k}$ denotes the interspecific competition rate. The transmission of the disease in the prey is governed by the bilinear incidence rate $bSI$, where b represent the incidence rate of infected prey to susceptible prey, $S\left(t\right)$ and $I\left(t\right)$ denote the densities of the susceptible prey, infected prey. Moreover, the parameters $\alpha$ and $\beta$ represent the capturing rates of predator to the susceptible and infected prey, respectively; m is the so-called half saturation constant; $d_1$ and $d_2$ are the natural death rates of the infected prey and predator. All coefficients mentioned are positive.

As a matter of a fact, most realistic ecosystems are affected by environmental noise (see [16] [17] [18] [19]. Motivated by the method in [20], we introduce to system (1) Gaussian white noise which are directly proportional to S(t), I(t) and P(t), and obtain

$\{\begin{array}{l}\text{d}S=\left[rS\left(1-\frac{S}{k}\right)-\frac{\alpha SP}{mP+S+I}-bSI\right]\text{d}t+{\sigma }_{1}S\text{d}{B}_1,\\ \text{d}I=\left[bSI-d_{1}I-\frac{\beta IP}{mP+S+I}\right]\text{d}t+{\sigma }_{2}I\text{d}{B}_2,\\ \text{d}P=\left[\frac{c\alpha SP}{mP+S+I}+\frac{c\beta IP}{mP+S+I}-d_{2}P\right]\text{d}t+{\sigma }_{3}P\text{d}{B}_3,\end{array}$ (2)

where ${\left\{B_i\left(t\right)\right\}}_{t\ge 0}\left(i=1,2,3\right)$ are mutually independent standard Brownian motions, and ${\sigma }_i\left(i=1,2,3\right)$ denote the intensities of the white noise.

Throughout this article, let $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}{\left\{F_t\right\}}_{t\ge 0}\mathrm{,}\mathbb{P}\right)$ be a complete probability space with a filtration ${\left\{F_t\right\}}_{t\ge 0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while $\left\{F_0\right\}$ contains all $\mathbb{P}$ -null sets) and let $B_i\left(t\right)$ be defined on the complete probability space, $i=1,2,3$. Denote ${ℝ}_{+}^d=\left\{x=\left(x_1,\cdots ,x_d\right)\in {ℝ}^d:x_i>0,1\le i\le d\right\}.$

In this article, in order to better study the spread of infectious diseases among interacting populations, it is more practical to establish a more accurate random ecological infectious disease model. We will concentrate on the dynamics of the stochastic model (2). The rest of the article is organized as follows. In Section 2, the existence and uniqueness of the positive solution is proved for system (2). In Section 3, we analyze the extinction and persistence of the infected prey. In Section 4, we obtain the conditions on the existence of stationary distribution for model (2). In Section 5, numerical simulations are conducted to support the theoretical results. A conclusion is given in the last section.

2. Existence and Uniqueness

To begin with, we recall some basic notations in stochastic differential equation. let $X\left(t\right)$ be a regular time-homogeneous Markov process in ${ℝ}^d$ described by the stochastic differential equation

$\text{d}X\left(t\right)=f\left(X\left(t\right)\right)\text{d}t+g\left(X\left(t\right)\right)\text{d}B\left(t\right).$ (3)

The diffusion matrix of the process $X\left(t\right)$ is defined as $A\left(x\right)=\left(a_{ij}\left(x\right)\right)$, $a_{ij}\left(x\right)=g^i\left(x\right)g^j\left(x\right)$. Furthermore, the differential operator L is defined by

$LV\left(x\right)={\displaystyle \underset{i=1}{\overset{d}{\sum }}}{f}_i\left(x\right)\frac{\partial V\left(x\right)}{\partial x_i}+\frac{1}{2}{\displaystyle \underset{i,j=1}{\overset{d}{\sum }}}{\left[g^{\text{T}}\left(x\right)g\left(x\right)\right]}_{ij}\frac{{\partial }^{2}V\left(x\right)}{\partial x_i\partial x_j},$

where $V\in {\mathcal{C}}^2\left({ℝ}^d\mathrm{,}{ℝ}_{+}\right)$.

To investigate the dynamical behavior of the model, the first concern is whether the solution is global and positive. In this section, we show that there exists a unique global positive solution of system (2) by constructing an appropriate Lyapunov function.

Theorem 2.1. For any initial value $\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)\in {ℝ}_{+}^3$ , there exists a unique positive solution $\left(S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right)$ of system (2) on $t\ge 0$ and the solution will remain in ${ℝ}_{+}^3$ with probability one, that is to say, $\left(S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right)\in {ℝ}_{+}^3$ for all $t\ge 0$ almost surely (a.s.).

Proof. Since the coefficients of system (2) are locally Lipschitz continuous, then for any initial value $\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)\in {ℝ}_{+}^3$ there is a unique local solution

$\left(S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right)$ on $t\in \left[\mathrm{0,}{\tau }_e\right)$, where ${\tau }_e$ is the explosion time. Now, let us

show that this solution is global, i.e., ${\tau }_e=\infty$ a.s.. Let $n_0>0$ be sufficiently large such that $S\left(0\right)$, $I\left(0\right)$ and $P\left(0\right)$ lie within the interval $\left[\frac{1}{n_0}\mathrm{,}{n}_0\right]$. For each integer $n\ge n_0$, define stopping-times

${\tau }_n=\inf \left\{t\in \left[0,{\tau }_e\right):\min \left\{S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right\}\le \frac{1}{n_0}\text{or}\max \left\{S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right\}\ge n\right\},$

set $\inf \varnothing =\infty$ ( $\varnothing$ represents the empty set). It is clear that ${\tau }_n$ is increasing as $n\to \infty$. Let ${\tau }_{\infty }={\lim }_{n\to \infty }{\tau }_n$, then ${\tau }_{\infty }\le {\tau }_e$ a.s.. In the following, it only needs to show that ${\tau }_{\infty }=\infty$ a.s.. If this statement is violated, then there exists a constant $T>0$ and $\epsilon \in \left(\mathrm{0,1}\right)$ such that $\mathbb{P}\left\{{\tau }_{\infty }\le T\right\}>\epsilon$. Consequently, there exists an integer $n_1\ge n_0$ such that

$\mathbb{P}\left\{{\tau }_n\le T\right\}\ge \epsilon \text{for}\text{all}n\ge n_1\mathrm{.}$ (4)

Define a $C^2$ -function $V\mathrm{<mo>\; :\; </mo>}{ℝ}_{+}^3\to {ℝ}_{+}$ by

$\tilde{V}\left(t\right)=\tilde{V}\left(S\mathrm{,}I\mathrm{,}P\right)=\left(S-a-a\ln \frac{S}{a}\right)+\left(I-1-\ln I\right)+\frac{1}{c}\left(P-1-\ln P\right)\mathrm{,}$

where a is a positive constant to be determined later. The nonnegativity of the function can be obtained from $u-1-\log u\ge 0$ for any $u>0$.

Let $n\ge n_0$ and $T>0$ be arbitrary. Applying Itô’s formula to $\tilde{V}$, we obtain that

$\begin{array}{l}\text{d}\tilde{V}=\left(1-\frac{a}{S}\right)\text{d}S+\frac{a{\sigma }_1^2}{2}\text{d}t+\left(1-\frac{1}{I}\right)\text{d}I+\frac{{\sigma }_2^2}{2}\text{d}t+\frac{1}{c}\left(1-\frac{1}{P}\right)\text{d}P+\frac{{\sigma }_3^2}{2c}\text{d}t\\ :=L\tilde{V}\text{d}t+\left(S-a\right){\sigma }_1\text{d}{B}_1+\left(I-1\right){\sigma }_2\text{d}{B}_2+\frac{\left(P-1\right){\sigma }_3\text{d}{B}_3}{c},\end{array}$

where

$\begin{array}{c}L\tilde{V}=\left(S-a\right)\left(r-\frac{r}{k}S-\frac{\alpha P}{mP+S+I}-bI\right)+\left(I-1\right)\left(bS-d_1-\frac{\beta P}{mP+S+I}\right)\\ +\frac{P-1}{c}\left(\frac{c\alpha S}{mP+S+I}+\frac{c\beta I}{mP+S+I}-d_2\right)+\frac{a{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2c}\\ =rS-\frac{r}{k}{S}^2-\frac{\alpha SP}{mP+S+I}-bSI-ra+\frac{ra}{k}S+\frac{a\alpha P}{mP+S+I}+abI+bSI-d_{1}I\\ -\frac{\beta PI}{mP+S+I}-bS+d_1+\frac{\beta P}{mP+S+I}+\frac{\alpha SP}{mP+S+I}+\frac{\beta PI}{mP+S+I}-\frac{d_2}{c}P\\ -\frac{\alpha S}{mP+S+I}-\frac{\beta I}{mP+S+I}+\frac{d_2}{c}+\frac{a{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2c}\end{array}$

$\begin{array}{c}\le rS+\frac{ra}{k}S-\frac{r}{k}{S}^2+\frac{a\alpha P}{mP+S+I}+\left(ab-d_1\right)I+\frac{\beta P}{mP+S+I}\\ -ra+d_1+\frac{d_2}{c}+\frac{a{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2c}.\end{array}$

Choose $a=\frac{d_1}{b}$ such that $ab-d_1=0$, then

$\begin{array}{c}L\tilde{V}\le rS+\frac{r}{k}\frac{d_1}{b}S-\frac{r}{k}{S}^2+\frac{\alpha }{m}\frac{d_1}{b}+\frac{\beta }{m}-\frac{r{d}_1}{b}+d_1+\frac{d_2}{c}+\frac{a{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2c}\\ \le \frac{k{\left(r+\frac{r{d}_1}{kb}\right)}^2}{4r}+\frac{\alpha }{m}\frac{d_1}{b}+\frac{\beta }{m}-\frac{r{d}_1}{b}+d_1+\frac{d_2}{c}+\frac{d_1}{b}\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2c}:=N,\end{array}$

where N is a positive parameter. Thus, we can obtain

$\text{d}\tilde{V}\le L\tilde{V}\text{d}t+\left(S-a\right){\sigma }_1\text{d}{B}_1+\left(I-1\right){\sigma }_2\text{d}{B}_2+\frac{\left(P-1\right){\sigma }_3\text{d}{B}_3}{c}.$ (5)

Interacting and taking the expectation of both sides of (5) yield

$E\left(V\left(S\left({\tau }_n\wedge T\right),I\left({\tau }_n\wedge T\right),P\left({\tau }_n\wedge T\right)\right)\right)\le V\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)+NT.$ (6)

Set ${\Omega }_n={\tau }_n\le T$ for $n\ge n_1$, then by (4), we have $\mathbb{P}\left({\Omega }_n\right)\ge \epsilon$. Noting that for every $\omega \in {\Omega }_n$, there is at least one of $S\left({\tau }_n\mathrm{,}\omega \right)\mathrm{,}I\left({\tau }_n\mathrm{,}\omega \right)\mathrm{,}P\left({\tau }_n\mathrm{,}\omega \right)$ that equals to either n or $\frac{1}{n}$, and

$V\left(S\left({\tau }_n\right)\mathrm{,}I\left({\tau }_n\right)\mathrm{,}P\left({\tau }_n\right)\right)\mathrm{,}$

is no less than

$n-1-\ln n\text{or}\frac{1}{n}-1-\ln \frac{1}{n}\mathrm{.}$

That is,

$V\left(S\left({\tau }_n\right)\mathrm{,}I\left({\tau }_n\right)\mathrm{,}P\left({\tau }_n\right)\right)\ge \left(n-1-\ln n\right)\wedge \left(\frac{1}{n}-1-\ln \frac{1}{n}\right)\mathrm{.}$

By (6), we can obtain

$\begin{array}{c}\infty >V\left(S\left(0\right),I\left(0\right),P\left(0\right)\right)+NT\ge E\left(1_{{\Omega }_n\left(\omega \right)}V\left(S\left({\tau }_n\right),I\left({\tau }_n\right),P\left({\tau }_n\right)\right)\right)\\ =\mathbb{P}\left({\Omega }_n\right)V\left(S\left({\tau }_n\right),I\left({\tau }_n\right),P\left({\tau }_n\right)\right)\ge \epsilon V\left(S\left({\tau }_n\right),I\left({\tau }_n\right),P\left({\tau }_n\right)\right)\\ \ge \epsilon \left[\left(n-1-\ln n\right)\wedge \left(\frac{1}{n}-1-\ln \frac{1}{n}\right)\right].\end{array}$

Taking $n\to \infty$ induces $\infty >+\infty$, which is a contradiction. Hence, we have ${\tau }_{\infty }=\infty$, a.s.. The conclusion is confirmed.

3. Extinction and Persistence

According to the theory in [21], the basic reproductive number $R_0$ is a threshold to control whether the disease will spread. If $R_0\le 1$, the disease disappear; If $R_0>1$, the infectious population will be persistence in the mean. It is easy to conclude the basic reproductive number $R_0=\frac{bk}{d_1}$ and system (1) has the following properties:

• If $R_0\le 1$, the disease-free equilibria $E_0\left(k,0,0\right)$ is globally asymptotically stable, the disease disappear;

• If $R_0>1$, the infectious prey population will be persistence in the mean.

In this section, we turn to establish sufficient criteria on the extinction and persistence of infected prey population for the stochastic system (2). Before giving our main results, we first recall the following lemma.

Lemma 3.1. ( [22] ) Let $X\left(t\right)\in C\left(\Omega \times \left[\mathrm{0,}\infty \right)\mathrm{,}{R}_{+}\right)$ .

1) If there exists $T>0$, ${\lambda }_0>0$, $\lambda$ and $n_i$ such that

$\ln X\left(t\right)\le \lambda t-{\lambda }_0{\displaystyle {\int }_0^t}X\left(s\right)\text{d}s+{\displaystyle \underset{i=1}{\overset{j}{\sum }}}{n}_{i}B\left(t\right)\text{a}\text{.s}.,\text{for}t\ge T,$

then

$\{\begin{array}{l}{\langle X\left(t\right)\rangle }^{*}\le \frac{\lambda }{{\lambda }_0}\text{a}\text{.s}.,\text{if}\lambda \ge 0;\\ \underset{t\to \infty }{lim}X\left(t\right)=0\text{a}\text{.s}.,\text{if}\lambda <0.\end{array}$

where ${\langle X\left(t\right)\rangle }^{\mathrm{<mo>\; *\; </mo>}}=\underset{t\to \infty }{\lim \sup }\langle X\left(t\right)\rangle$.

2) If there exists $T>0$, ${\lambda }_0>0$, $\lambda >0$ and $n_i$ such that

$\ln X\left(t\right)\ge \lambda t-{\lambda }_0{\displaystyle {\int }_0^t}X\left(s\right)\text{d}s+{\displaystyle \underset{i=1}{\overset{j}{\sum }}}{n}_{i}B\left(t\right)\text{a}\text{.s}\text{.},\text{for}t\ge T,$

then

${\langle X\left(t\right)\rangle }_{\mathrm{<mo>\; *\; </mo>}}\ge \frac{\lambda }{{\lambda }_0}\text{a}\text{.s}\mathrm{.,}$

where ${\langle X\left(t\right)\rangle }_{\mathrm{<mo>\; *\; </mo>}}=\underset{t\to \infty }{\lim \inf }\langle X\left(t\right)\rangle$.

Theorem 3.2. Let $\left(S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right)$ be the solution of system (2) with any initial value $\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)\in {ℝ}_{+}^3$. If $r>\frac{{\sigma }_1^2}{2}$, $\alpha +\beta <\frac{d_2+\frac{{\sigma }_3^2}{2}}{c}$ and $R_0^s=\frac{bk\left(r-\frac{{\sigma }_1^2}{2}\right)}{r\left(d_1+\frac{{\sigma }_2^2}{2}\right)}<1$, then

$\underset{t\to \infty }{\lim \sup }\frac{\ln I\left(t\right)}{t}\le \left(d_1+\frac{{\sigma }_2^2}{2}\right)\left(R_0^s-1\right)<0\text{a}\text{.s}.,$

$\underset{t\to \infty }{\lim \sup }\frac{\ln P\left(t\right)}{t}\le -\left[d_2+\frac{{\sigma }_3^2}{2}-c\left(\alpha +\beta \right)\right]<0\text{a}\text{.s}\mathrm{.,}$

and

$\underset{t\to \infty }{lim}\frac{1}{t}{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s=\frac{k\left(r-\frac{{\sigma }_1^2}{2}\right)}{r}\text{a}\text{.s}\mathrm{..}$

Proof. By the Itô’s formula, we have

$\begin{array}{c}\text{d}\ln S=\left(r-\frac{r}{k}S-\frac{\alpha P}{mP+S+I}-bI-\frac{{\sigma }_1^2}{2}\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\\ \le \left(r-\frac{{\sigma }_1^2}{2}-\frac{r}{k}S\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\mathrm{,}\end{array}$ (7)

integrating Equation (7) from 0 to t and dividing it by t, we obtain

$\frac{1}{t}\ln \frac{S\left(t\right)}{S\left(0\right)}\le r-\frac{{\sigma }_1^2}{2}-\frac{r}{k}\frac{{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s}{t}+\frac{{\displaystyle {\int }_0^t}{\sigma }_1\text{d}{B}_1\left(\theta \right)}{t},$

applying Lemma 3.1, we have

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s\le \frac{k\left(r-\frac{{\sigma }_1^2}{2}\right)}{r}\text{a}\text{.s}\mathrm{..}$ (8)

Similarly, by the Itô’s formula, we have

$\begin{array}{c}\text{d}\ln I=\left(bS-d_1-\frac{\beta P}{mP+S+I}-\frac{{\sigma }_2^2}{2}\right)dt+{\sigma }_2\text{d}{B}_2\left(t\right)\\ \le \left[bS-\left(d_1+\frac{{\sigma }_2^2}{2}\right)\right]\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right),\end{array}$ (9)

integrate Equation (9) from 0 to t and divide it by t yields, we obtain

$\frac{\text{1}}{t}\ln \frac{I\left(t\right)}{I\left(0\right)}\le \frac{b{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s}{t}-\left(d_1+\frac{{\sigma }_2^2}{2}\right)+\frac{{\displaystyle {\int }_0^t}{\sigma }_1\text{d}{B}_2\left(\theta \right)}{t}\mathrm{,}$

so

$\underset{t\to \infty }{{\displaystyle \lim \sup }}\frac{\ln I\left(t\right)}{t}\le -\left(d_1+\frac{{\sigma }_2^2}{2}\right)+\frac{bk\left(r-\frac{{\sigma }_1^2}{2}\right)}{r}=\left(d_1+\frac{{\sigma }_2^2}{2}\right)\left(R_0^s-1\right)<0\text{a}\text{.s}\mathrm{..}$

By the Itô’s formula, we also have

$\begin{array}{c}\text{d}\ln P=\left(\frac{c\alpha S}{mP+S+I}+\frac{c\beta I}{mP+S+I}-d_2-\frac{{\sigma }_3^2}{2}\right)\text{d}t+{\sigma }_3\text{d}{B}_3\left(t\right)\\ \le \left(c\alpha +c\beta -d_2-\frac{{\sigma }_3^2}{2}\right)\text{d}t+{\sigma }_3\text{d}{B}_3\left(t\right)\mathrm{,}\end{array}$ (10)

integrating Equation (9) from 0 to t and dividing it by t, one can get

$\frac{\text{1}}{t}\ln \frac{P\left(t\right)}{P\left(0\right)}\le \left(c\alpha +c\beta -d_2-\frac{{\sigma }_3^2}{2}\right)+\frac{{\displaystyle {\int }_0^t}{\sigma }_3\text{d}{B}_3\left(\theta \right)}{t}\mathrm{,}$

so we can obtain that

$\underset{t\to \infty }{{\displaystyle \lim \sup }}\frac{\ln P\left(t\right)}{t}\le -\left[d_2+\frac{{\sigma }_3^2}{2}-c\left(\alpha +\beta \right)\right]<0\text{a}\text{.s}\mathrm{..}$

On the other hand,

$\text{d}\ln S=\left(r-\frac{r}{k}S-\frac{\alpha P}{mP+S+I}-bI-\frac{{\sigma }_1^2}{2}\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\mathrm{,}$

since $\underset{t\to \infty }{{\displaystyle \lim \sup }}\frac{\ln P\left(t\right)}{t}<0\text{a}\text{.s}\mathrm{.}$, there exists an arbitrarily small constant $\epsilon >0$ such that when $t>T$, we obtain $\frac{P}{mP+S+I}<\epsilon$, so

$\text{d}\ln S\ge \left(r-\frac{{\sigma }_1^2}{2}-\frac{r}{k}S-\alpha \epsilon \right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\mathrm{,}$ (11)

integrating Equation (11) from 0 to t and dividing it by t, we have

$\frac{\text{1}}{t}\ln \frac{S\left(t\right)}{S\left(0\right)}\ge r-\frac{{\sigma }_1^2}{2}-\frac{r}{k}\frac{{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s}{t}-\alpha \epsilon +\frac{{\displaystyle {\int }_0^t}{\sigma }_1\text{d}{B}_1\left(\theta \right)}{t}\mathrm{,}$

applying Lemma 3.1 and the arbitrariness of $\epsilon$, we get

$\underset{t\to \infty }{\lim \inf }\frac{1}{t}{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s\ge \frac{k\left(r-\frac{{\sigma }_1^2}{2}\right)}{r}\text{a}\text{.s}\mathrm{..}$ (12)

From (8) and (12), we obtain

$\underset{t\to \infty }{lim}\frac{1}{t}{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s=\frac{k\left(r-\frac{{\sigma }_1^2}{2}\right)}{r}\text{a}\text{.s}\mathrm{..}$

Theorem 3.3. Let $\left(S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)\right)$ be the solution of system (2) with any initial value $\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)\in {ℝ}_{+}^3$. If $r>\frac{{\sigma }_1^2}{2}$, $\alpha +\beta <\frac{d_2+\frac{{\sigma }_3^2}{2}}{c}$ and $R_0^s>1$, then

$\underset{t\to \infty }{\lim \inf }\frac{\ln I\left(t\right)}{t}\ge \left(d_1+\frac{{\sigma }_2^2}{2}\right)\left(R_0^s-1\right)>0\text{a}\text{.s}\mathrm{..}$

Proof. By the Itô’s formula, one we can obtain

$\text{d}\ln I=\left(bS-d_1-\frac{\beta P}{mP+S+I}-\frac{{\sigma }_2^2}{2}\right)\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right),$

according to Theorem 3.2, there exists an arbitrarily small constant $\epsilon >0$ such

that when $t>T$, we obtain $\frac{P}{mP+S+I}<\epsilon$, so

$\text{d}\ln I\ge \left[bS-\beta \epsilon -\left(d_1+\frac{{\sigma }_2^2}{2}\right)\right]\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right),$ (13)

integrate Equation (13) from 0 to t and divide it by t yields, we can obtain

$\begin{array}{c}\frac{1}{t}\ln \frac{I\left(t\right)}{I\left(0\right)}\ge \frac{b{\displaystyle {\int }_0^t}S\left(s\right)\text{d}s}{t}-\beta \epsilon -\left(d_1+\frac{{\sigma }_2^2}{2}\right)+\frac{{\displaystyle {\int }_0^t}{\sigma }_1\text{d}{B}_2\left(\theta \right)}{t}\\ =\frac{bk\left(r-\frac{{\sigma }_1^2}{2}\right)}{r}-\beta \epsilon -\left(d_1+\frac{{\sigma }_2^2}{2}\right)\\ \ge \left(d_1+\frac{{\sigma }_2^2}{2}\right)\left(R_0^s-1\right)>0\text{a}\text{.s}\mathrm{..}\end{array}$

According to Theorem 3.2 and 3.3, $R_0^s$ is the threshold for the infected prey population. The disease will go to extinction if $R_0^s<1$, and it will exist in the long time if $R_0^s>1$.

4. Stationary Distribution

Now we present the following lemma.

Lemma 4.1. ( [23] ) The Markov process $X\left(t\right)$ described by Equation (3) has a unique ergodic stationary distribution $\mu \left(\mathrm{.}\right)$ if there exists a bounded domain $D\subset {ℝ}^d$ with regular boundary $\Gamma$ [(B.1).] There is a positive number M such that ${\displaystyle {\sum }_{i,j=1}^d}{a}_{ij}\left(x\right){\xi }_i{\xi }_j\ge M{\left|\xi \right|}^2$ , $x\in D$ , $\xi \in {ℝ}^d$ . [(B.2).] There exists a nonnegative C²-function V such that LV is negative for any ${ℝ}^d\backslash D$ . Then

${\mathbb{P}}_x\left\{\underset{T\to \infty }{lim}\frac{1}{T}{\displaystyle {\int }_0^T}f\left(X\left(t\right)\right)\text{d}t={\displaystyle {\int }_{{ℝ}^d}}f\left(x\right)\mu \left(\text{d}x\right)\right\}=\mathrm{1,}$

for all $x\in {ℝ}^d$, where $f\left(\mathrm{.}\right)$ is a function integrable with respect to the measure $\mu$.

Theorem 4.2. Assume that ${\hat{R}}_0^s=\frac{r+d_2+\frac{{\sigma }_3^2}{2}-\frac{{\sigma }_1^2}{2}-\frac{\alpha }{m}-c\alpha -c\beta }{\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)}>1$ , and

environmental noises are small enough that ${\sigma }_2^2<2d_1$ and ${\sigma }_3^2<2d_2$ . Then for any initial value $\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)\in {ℝ}_{+}^3$ , there exists a unique stationary distribution for system (2) and it is ergodic.

Proof. In order to prove Theorem 4.2, first we need to verify Lemma 4.1. To verify B.2, we need to proof there exists a neighborhood $D\subset {ℝ}_{+}^3$ and a nonnegative C²-function V such that for any $\left(S\mathrm{,}I\mathrm{,}P\right)\subset {ℝ}_{+}^3\backslash D$, LV is negative.

Define a C²-function

$\bar{V}=M\left(-\log S-\frac{r}{bk}\log I+\frac{b}{d_1}I\right)+\frac{1}{2}{\left(S+I+\frac{P}{c}\right)}^2\mathrm{,}$

where

$M=\frac{bk}{r\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)}\max \left\{2+{\bar{f}}_1+f_2+f_3\mathrm{,2}+f_1+{\bar{f}}_2+f_3\mathrm{,1}+f_1+f_2+{\bar{f}}_3\right\}$

and functions $f_1$, $f_2$, $f_3$, ${\bar{f}}_1$, ${\bar{f}}_2$ and ${\bar{f}}_3$ will be determined later. There exists a unique minimum point $\left(\tilde{S}\mathrm{,}\tilde{I}\mathrm{,}\tilde{P}\right)$ of $\bar{V}$.

Define a nonnegative C²-Lyapunov function

$V=M\left(-\log S-\frac{r}{bk}\log I+\frac{b}{d_1}I+\log P\right)+\frac{1}{2}{\left(S+I+\frac{P}{c}\right)}^2-\bar{V}\left(\tilde{S}\mathrm{,}\tilde{I}\mathrm{,}\tilde{P}\right)\mathrm{,}$

denote

$V_1=M\left(-\log S-\frac{r}{bk}\log I+\frac{b}{d_1}I+\log P\right),V_2=\frac{1}{2}{\left(S+I+\frac{P}{c}\right)}^2.$

By Itô’s formula, we get

$\begin{array}{c}L{V}_1=-\frac{1}{S}\left[rS\left(1-\frac{S}{k}\right)-\frac{\alpha SP}{mP+S+I}-bSI\right]+\frac{{\sigma }_1^2}{2}\\ -\frac{r}{bk}\frac{1}{I}\left(bSI-d_{1}I-\frac{\beta IP}{mP+S+I}\right)+\frac{r{\sigma }_2^2}{2bk}\\ +\frac{b}{d_1}\left(bSI-d_{1}I-\frac{\beta IP}{mP+S+I}\right)\\ +\frac{c\alpha S}{mP+S+I}+\frac{c\beta I}{mP+S+I}-d_2-\frac{{\sigma }_3^2}{2}\end{array}$

$\begin{array}{c}\le -\left(r-\frac{{\sigma }_1^2}{2}-\frac{\alpha }{m}\right)+\frac{r}{k}S+bI-\frac{r}{k}S+\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\\ +\frac{b^2}{d_1}SI-bI+c\alpha +c\beta -d_2-\frac{{\sigma }_3^2}{2}\\ =-\left(r+d_2+\frac{{\sigma }_3^2}{2}-\frac{{\sigma }_1^2}{2}-\frac{\alpha }{m}-c\alpha -c\beta \right)+\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)+\frac{b^2}{d_1}SI,\end{array}$

and

$\begin{array}{c}L{V}_2=\left(S+I+\frac{P}{c}\right)\left(rS-\frac{r}{k}{S}^2-d_{1}I-\frac{d_2}{c}P\right)+\frac{1}{2}\left({\sigma }_1^2{S}^2+{\sigma }_2^2{I}^2+\frac{{\sigma }_3^2}{c^2}{P}^2\right)\\ =r{S}^2-\frac{r}{k}{S}^3-d_{1}SI-\frac{d_2}{c}SP+rSI-\frac{r}{k}I{S}^2-d_1{I}^2-\frac{d_2}{c}IP+\frac{r}{c}PS\\ -\frac{r}{ck}P{S}^2-\frac{d_1}{c}PI-\frac{d_2}{c^2}{P}^2+\frac{1}{2}\left({\sigma }_1^2{S}^2+{\sigma }_2^2{I}^2+\frac{{\sigma }_3^2}{c^2}{P}^2\right)\\ \le r{S}^2-\frac{r}{k}{S}^3+rSI-d_1{I}^2+\frac{r}{c}PS-\frac{d_2}{c^2}{P}^2+\frac{1}{2}\left({\sigma }_1^2{S}^2+{\sigma }_2^2{I}^2+\frac{{\sigma }_3^2}{c^2}{P}^2\right),\end{array}$

according to Young inequality, $xy\le frac25x^{\frac{5}{2}}+\frac{3}{5}{y}^{\frac{5}{3}}$, so

$\begin{array}{c}L{V}_2\le r{S}^2-\frac{r}{k}{S}^3+\left(\frac{2r}{5}+\frac{2r}{5c}\right)S^{\frac{5}{2}}+\frac{{\sigma }_1^2}{2}{S}^2-\left(d_1-\frac{{\sigma }_2^2}{2}\right)I^2+\frac{3r}{5}{I}^{\frac{5}{3}}\\ -\frac{1}{c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{P}^{\frac{5}{3}}\mathrm{,}\end{array}$

so

$\begin{array}{c}LV\le M\left[-\left(r+d_2+\frac{{\sigma }_3^2}{2}-\frac{{\sigma }_1^2}{2}-\frac{\alpha }{m}-c\alpha -c\beta \right)+\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)+\frac{b^2}{d_1}SI\right]\\ -\frac{r}{k}{S}^3+\left(\frac{2r}{5}+\frac{2r}{5c}\right)S^{\frac{5}{2}}+\left(\frac{{\sigma }_1^2}{2}+r\right)S^2-\left(d_1-\frac{{\sigma }_2^2}{2}\right)I^2\\ +\frac{3r}{5}{I}^{\frac{5}{3}}-\frac{1}{c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{P}^{\frac{5}{3}}\\ =-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+\frac{M{b}^2}{d_1}SI-\frac{r}{k}{S}^3+\left(\frac{2r}{5}+\frac{2r}{5c}\right)S^{\frac{5}{2}}\\ +\left(\frac{{\sigma }_1^2}{2}+r\right)S^2-\left(d_1-\frac{{\sigma }_2^2}{2}\right)I^2+\frac{3r}{5}{I}^{\frac{5}{3}}-\frac{1}{c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{P}^{\frac{5}{3}},\end{array}$

where

$\begin{array}{l}{\hat{R}}_0^s=\frac{r+d_2+\frac{{\sigma }_3^2}{2}-\frac{{\sigma }_1^2}{2}-\frac{\alpha }{m}-c\alpha -c\beta }{\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)},\\ f_1\left(S\right)=-\frac{r}{k}{S}^3+\left(\frac{2r}{5}+\frac{2r}{5c}\right)S^{\frac{5}{2}}+\left(\frac{{\sigma }_1^2}{2}+r\right)S^2,\end{array}$

$f_2\left(I\right)=-\left(d_1-\frac{{\sigma }_2^2}{2}\right)I^2+\frac{3r}{5}{I}^{\frac{5}{3}},f_3\left(P\right)=-\frac{1}{c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{P}^{\frac{5}{3}}.$

We aim to prove that $LV\le -1$, consider the bounded set D

$D=\left\{\left(S\mathrm{,}I\mathrm{,}P\right)\mathrm{<mo>\; :\; </mo>}\epsilon \le S\le \frac{1}{\epsilon }\mathrm{,}\epsilon \le I\le \frac{1}{\epsilon }\mathrm{,}\epsilon \le P\le \frac{1}{\epsilon }\right\}\mathrm{,}$

then ${ℝ}_{+}^3\backslash D=D_1^c\cup D_2^c\cup D_3^c\cup D_4^c\cup D_5^c\cup D_6^c$, with

$D_1^c=\left\{\left(S,I,P\right)\in {ℝ}_{+}^3:0<S<\epsilon \right\},D_2^c=\left\{\left(S,I,P\right)\in {ℝ}_{+}^3:0<I<\epsilon \right\},$

$D_3^c=\left\{\left(S,I,P\right)\in {ℝ}_{+}^3:0<P<\epsilon \right\},D_4^c=\left\{\left(S,I,P\right)\in {ℝ}_{+}^3:S>\frac{1}{\epsilon }\right\},$

$D_5^c=\left\{\left(S,I,P\right)\in {ℝ}_{+}^3:I>\frac{1}{\epsilon }\right\},D_6^c=\left\{\left(S,I,P\right)\in {ℝ}_{+}^3:P>\frac{1}{\epsilon }\right\},$

and $\epsilon$ is a sufficiently small positive constant satisfying the following conditions

$\begin{array}{l}-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+\frac{b^{2}M\epsilon }{d_1}+f_1\left(S\right)+{\bar{f}}_2\left(I\right)+f_3\left(P\right)\le -1,\\ \frac{b^{2}M\epsilon }{d_1}<\frac{1}{2}\left(d_1-\frac{{\sigma }_2^2}{2}\right),\end{array}$ (14)

$\begin{array}{l}-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+\frac{b^{2}M\epsilon }{d_1}+{\bar{f}}_1\left(S\right)+f_2\left(I\right)+f_3\left(P\right)\le -1,\\ \frac{b^{2}M\epsilon }{d_1}<\frac{r}{2k},\end{array}$ (15)

$\begin{array}{l}-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+f_1\left(S\right)+f_2\left(I\right)\\ -\frac{1}{c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{\epsilon }^{\frac{5}{3}}\le -1,\end{array}$ (16)

$-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+A-\frac{r}{2k{\epsilon }^3}\le -\mathrm{1,}$ (17)

$-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+B-\frac{1}{2{\epsilon }^2}\left(d_1-\frac{{\sigma }_2^2}{2}\right)\le -\mathrm{1,}$ (18)

$-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+C-\frac{1}{2{\epsilon }^2{c}^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)\le -\mathrm{1,}$ (19)

where

$A=\underset{\left(S,I,P\right)\in {ℝ}_{+}^3}{{\displaystyle \sup }}\left\{{\bar{f}}_1\left(S\right)+f_2\left(I\right)+f_3\left(P\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}\right\}<\infty ,$

$B=\underset{\left(S,I,P\right)\in {ℝ}_{+}^3}{{\displaystyle \sup }}\left\{f_1\left(S\right)+{\bar{f}}_2\left(I\right)+f_3\left(P\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}\right\}<\infty ,$

$C=\underset{\left(S,I,P\right)\in {ℝ}_{+}^3}{{\displaystyle \sup }}\left\{f_1\left(S\right)+f_2\left(I\right)+{\bar{f}}_3\left(P\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}\right\}<\infty .$

Case 1. If $\left(S\mathrm{,}I\mathrm{,}P\right)\in D_1^c$, $SI<\epsilon I<\epsilon \left(1+I^2\right)$, we have

$\begin{array}{l}LV\le -M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+\frac{b^{2}M\epsilon }{d_1}\\ +f_1\left(S\right)+\left[-\frac{d_1-\frac{{\sigma }_2^2}{2}}{2}+\frac{b^{2}M\epsilon }{d_1}\right]I^2+{\bar{f}}_2\left(I\right)+f_3\left(P\right)\le -1,\end{array}$

which follows from (14), where

${\bar{f}}_2\left(I\right)=-\frac{1}{2}\left(d_1-\frac{{\sigma }_2^2}{2}\right)+\frac{3r}{5}{I}^{\frac{5}{3}}\mathrm{.}$

Case 2. If $\left(S\mathrm{,}I\mathrm{,}P\right)\in D_2^c$, $SI<\epsilon S<\epsilon \left(1+S^3\right)$, we have

$\begin{array}{c}LV\le -M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+\frac{b^{2}M\epsilon }{d_1}+{\bar{f}}_1\left(S\right)\\ +\left(-\frac{r}{2k}+\frac{b^{2}M\epsilon }{d_1}\right)S^3+f_2\left(I\right)+f_3\left(P\right)\\ \le -1,\end{array}$

which follows from (15), where

${\bar{f}}_1\left(S\right)=-\frac{r}{2k}{S}^3+\left(\frac{2r}{5}+\frac{2r}{5c}\right)S^{\frac{5}{2}}+\left(\frac{{\sigma }_1^2}{2}+r\right)S^2\mathrm{.}$

Case 3. If $\left(S\mathrm{,}I\mathrm{,}P\right)\in D_3^c$, $SI\le \frac{2}{5}{S}^{\frac{5}{2}}+\frac{3}{5}{I}^{\frac{5}{3}}$, we have

$\begin{array}{c}LV\le -M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+f_1\left(S\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}\\ +f_2\left(I\right)+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}-\frac{1}{c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{\epsilon }^{\frac{5}{3}}\\ \le -1,\end{array}$

which follows from (16).

Case 4. If $\left(S\mathrm{,}I\mathrm{,}P\right)\in D_4^c$, $SI\le \frac{2}{5}{S}^{\frac{5}{2}}+\frac{3}{5}{I}^{\frac{5}{3}}$, $-\frac{r}{2k}{S}^3<-\frac{r}{2k}{\epsilon }^3$, we have

$\begin{array}{c}LV\le M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+{\bar{f}}_1\left(S\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}\\ +f_2\left(I\right)+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}+f_3\left(P\right)-\frac{r}{2k}{\epsilon }^3\\ =-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+A-\frac{r}{2k{\epsilon }^3}\\ \le -1,\end{array}$

which follows from (17).

Case 5. If $\left(S\mathrm{,}I\mathrm{,}P\right)\in D_5^c$, $SI\le \frac{2}{5}{S}^{\frac{5}{2}}+\frac{3}{5}{I}^{\frac{5}{3}}$, $-\frac{1}{2}\left(d_1-\frac{{\sigma }_2^2}{2}\right)I^2<-\frac{1}{2{\epsilon }^2}\left(d_1-\frac{{\sigma }_2^2}{2}\right)$, we have

$\begin{array}{c}LV\le -M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+f_1\left(S\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}\\ +{\bar{f}}_2\left(I\right)+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}+f_3\left(P\right)-\frac{r}{2k}{\epsilon }^3\\ =-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+B-\frac{1}{2{\epsilon }^2}\left(d_1-\frac{{\sigma }_2^2}{2}\right)\\ \le -1,\end{array}$

which follows from (18).

Case 6. If $\left(S\mathrm{,}I\mathrm{,}P\right)\in D_6^c$, $SI\le \frac{2}{5}{S}^{\frac{5}{2}}+\frac{3}{5}{I}^{\frac{5}{3}}$, $-\frac{1}{2c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2\mathrm{<mo>\; <\; </mo>}-\frac{1}{2c^2{\epsilon }^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)$, we have

$\begin{array}{c}LV\le -M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+f_1\left(S\right)+\frac{2}{5}\frac{b^{2}M}{d_1}{S}^{\frac{5}{2}}\\ +f_2\left(I\right)+\frac{3}{5}\frac{b^{2}M}{d_1}{I}^{\frac{5}{3}}+{\bar{f}}_3\left(P\right)-\frac{1}{2c^2{\epsilon }^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)\\ =-M\left[\frac{r}{bk}\left(d_1+\frac{\beta }{m}+\frac{{\sigma }_2^2}{2}\right)\left({\hat{R}}_0^s-1\right)\right]+C-\frac{1}{2{\epsilon }^2{c}^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)\\ \le -1,\end{array}$

which follows from (19), where ${\bar{f}}_3\left(P\right)=-\frac{1}{2c^2}\left(d_2-\frac{{\sigma }_3^2}{2}\right)P^2+\frac{3r}{5c}{P}^{\frac{5}{3}}$.

The proof of B.2 in Lemma 4.1 is completed. We get the conclusion that $LV\le 1$ on $D\subset {ℝ}_{+}^3$.

Direct computation shows that the diffusion matrix of system (2) is given by

$A=\left(\begin{array}{ccc}{\sigma }_1^2{S}^2& 0& 0\\ 0& {\sigma }_2^2{I}^2& 0\\ 0& 0& {\sigma }_3^2{P}^2\end{array}\right)\mathrm{.}$

It is clearly that the matrix A is positive definite for any compact subset of ${ℝ}_{+}^3$. The proof of B.1 in Lemma 4.1 is verified. According to Lemma 4.1, we know system (2) admits a stationary distribution.

5. Numerical Simulations

In this section, we conduct Example 1 - 3 numerical simulations to show the effect of noise on the dynamics of the system. Applying Milstein’s higher-order method [24] to system (2), we obtain the corresponding discretization equation as follows

$\{\begin{array}{l}{S}^{k+1}=S^k+\left(r{S}^k\left(1-\frac{S^k}{k}\right)-\frac{\alpha S^k{P}^k}{m{P}^k+S^k+I^k}-b{S}^k{I}^k\right)\Delta t+{\sigma }_1{S}^k\sqrt{\Delta t}{\xi }_k+\frac{{\sigma }_1^2}{2}{S}^k\left(\Delta t{\xi }_k^2-\Delta t\right),\\ I^{k+1}=I^k+\left(b{S}^k{I}^k-d_1{I}^k-\frac{\beta I^k{P}^k}{m{P}^k+S^k+I^k}\right)\Delta t+{\sigma }_2{I}^k\sqrt{\Delta t}{\eta }_k+\frac{{\sigma }_2^2}{2}{I}^k\left(\Delta t{\eta }_k^2-\Delta t\right),\\ P^{k+1}=P^k+\left(\frac{c\alpha S^k{P}^k}{m{P}^k+S^k+I^k}+\frac{c\beta I^k{P}^k}{m{P}^k+S^k+I^k}-d_2{P}^k\right)\Delta t+{\sigma }_3{P}^k\sqrt{\Delta t}{\varsigma }_k+\frac{{\sigma }_3^2}{2}{P}^k\left(\Delta t{\varsigma }_k^2-\Delta t\right),\end{array}$ (20)

where the time increment $\Delta t>0$, ${\xi }_k$, ${\eta }_k$ and ${\varsigma }_k$, $k=1,2,3,\cdots ,n$, are independent Gaussian random variables with normal distribution $N\left(\mathrm{0,1}\right)$, and ${\sigma }_i$, $1\le i\le 3$, are noise intensities. The unit of all time is days, the unit of the population $S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)$ are the individual.

Example 1. Choose the initial value $S\left(0\right)=5$, $I\left(0\right)=5$, $P\left(0\right)=2$, and $r=0.1$, $k=9$, $b=0.007$, $m=1$, $\alpha =0.02$, $\beta =0.02$, $c=0.8$, $d_1=0.1$, $d_2=0.05$.

By simple computations,

$R_0=0.63<1,\alpha +\beta =0.04<0.075=\frac{d_2+\frac{{\sigma }_3^2}{2}}{c},R_0^s=0.5925<1.$

The numerical simulation is shown in Figure 1, from which one can see that the susceptible prey is persistent, the infectious prey and predator will die out.

Example 2. Choose the initial value $S\left(0\right)=5$, $I\left(0\right)=5$, $P\left(0\right)=2$, and $r=0.1$, $k=9$, $b=0.015$, $m=1$, $\alpha =0.02$, $\beta =0.02$, $c=0.8$, $d_1=0.1$, $d_2=0.05$.

By simple computations,

$R_0=1.35>1,\alpha +\beta =0.04<0.075=\frac{d_2+\frac{{\sigma }_3^2}{2}}{c},R_0^s=1.2214>1.$

The conclusion of Theorem 3.3 holds, and the numerical simulation is shown in Figure 2. We note that the susceptible prey and the infectious prey will persist and the predator is going to die out.

Figure 1. Numerical simulation of the solution of system (1) and (2) for Example 1, respectively, where ${\sigma }_1=0.05,{\sigma }_2={\sigma }_3=0.1$.

Figure 2. Numerical simulation of the solution of system (1) and (2) for Example 2, respectively, ${\sigma }_1={\sigma }_2={\sigma }_3=0.1$.

Example 3. Choose the initial value $S\left(0\right)=5$, $I\left(0\right)=5$, $P\left(0\right)=5$, and $r=0.1$, $k=9$, $b=0.007$, $m=1$, $\alpha =0.02$, $\beta =0.02$, $c=0.8$, $d_1=0.01$, $d_2=0.01$, the (a) of Figure 3, ${\sigma }_1={\sigma }_2={\sigma }_3=0.02$, the (b) of Figure 3, ${\sigma }_1={\sigma }_2={\sigma }_3=0.01$. The red, black and blue dotted line represent the solution $S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)$ of model (2); the red, black and blue solid line represent the solution $S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}P\left(t\right)$ of model (1) for the same initial value $S\left(0\right)=5$, $I\left(0\right)=5$, $P\left(0\right)=5$.

Under this condition, by simple computations, (a) and (b) of Figure 3 satisfy

${\hat{R}}_0^s=1.2099>1,{\sigma }_2^2=0.0004<0.02=2d_1,{\sigma }_3^2=0.0004<0.02=2d_2.$

${\hat{R}}_0^s=1.2160>1,{\sigma }_2^2=0.0001<0.02=2d_1,{\sigma }_3^2=0.0001<0.02=2d_2.$

The numerical simulation is shown in Figure 3, which is consistent with our conclusion in Theorem 4.2. The difference between (a) and (b) of Figure 3 is the intensity of white noise. We can conclude that with the noise intensity decreases, the dynamics of stochastic system (2) is getting close to the deterministic system (1). Figure 4 shows the simulation of density functions, where system (2) has a unique stationary distribution.

Figure 3. Numerical simulation of the solution of system (1) and (2) for Example 3, respectively, picture(a) ${\sigma }_1={\sigma }_2={\sigma }_3=0.02$, picture(b) ${\sigma }_1={\sigma }_2={\sigma }_3=0.01$.

Figure 4. The density functions of $S\left(t\right)$, $I\left(t\right)$, and $P\left(t\right)$, respectively. Subfigures (a), (b), (c) ${\sigma }_1={\sigma }_2={\sigma }_3=0.02$. Subfigures (d), (e), (f) ${\sigma }_1={\sigma }_2={\sigma }_3=0.01$.

6. Conclusions

A stochastic predator-prey model with ratio-dependent functional response and disease in the prey population is investigated. We have shown that the existence of a unique positive global solution of the stochastic model (2). The threshold $R_0^s$ between stochastic persistence in the mean and extinction for infectious

prey population is given. If $r>\frac{{\sigma }_1^2}{2}$, $\alpha +\beta <\frac{d_2+\frac{{\sigma }_3^2}{2}}{c}$ and $R_0^s<1$, we obtained that the population of the infected prey and the predator will die out in time mean sense. If $r>\frac{{\sigma }_1^2}{2}$, $\alpha +\beta <\frac{d_2+\frac{{\sigma }_3^2}{2}}{c}$ and $R_0^s>1$, we note that the susceptible prey and the infectious prey are persistent and predator is going to die out. By constructing some suitable Lyapunov function, the existence of stationary distribution for both populations is established under certain parametric restrictions. If ${\hat{R}}_0^s>1$, and environmental noises are small enough that ${\sigma }_2^2<2d_1$ and ${\sigma }_3^2<2d_2$, then for any initial value $\left(S\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}P\left(0\right)\right)\in {ℝ}_{+}^3$, there exists a unique stationary distribution for system (2) and it is ergodic.

Acknowledgements

The bell for graduation sounded immediately, and three years passed unconsciously. Looking back on the three years, I find that there are sweet and bitter, tension, anxiety, joy, joy and tears. I am very grateful to my tutor, friends and family for their concern, care and tolerance to me. During my postgraduate period, my tutor Ma Jiying was strict in study and gave careful guidance. She also cared about my life and gave me meticulous care. I am very grateful to my mentor Ma Jiying. From I began writing my paper, horse the teacher give me a clear research direction are pointed out, in many times in the discussion class for my proposed amendments, encounter problems ma teacher and I discuss research, finally during the period of thesis writing finalized, teacher ma for all aspects of the draft, including writing format, paper typesetting, English grammar, logic of paper, proof of modified earnestly repeatedly. Every time we met, Ms. Ma urged me and my sister to understand the principle of each theorem and to derive it by ourselves. We should not only know how it is, but also know why it is so that we can successfully complete my graduation thesis. I have also learned from Mr. Ma a lot of hard research, serious study of the spirit of scientific research, here, I would like to express the most sincere thanks to my respected mentor.

Second, thank Wei Guoliang dean, professor held, Liu Xiping professor, the original three led, professor jammeh, an associate professor in the school and they taught me a very rich academic professional theory knowledge, really thank ginger teachers and Wang Ting counselor, would also like to thank the teacher elder sister of my horse sasha and yi qing, each encounter problems they were positive and I explore together, thanks to pool Lin Wei accompany and encouragement, I thank my roommate for the care and help.

Finally, I would like to express my heartfelt thanks to my parents and my boyfriend for their care and support and help in my study. Thank you all.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Anderson, R. and May, R. (1981) The Population Dynamics of Microparasites and Their Invertebrate Hosts. Philosophical Transactions of the Royal Society B Biological Sciences, 291, 451-524. https://doi.org/10.1098/rstb.1981.0005
[2]	Chakraborty, S., Kooi, B. and Chattopadhyay, J. (2015) Revealing the Role of Predator Interference in a Predator-Prey System with Disease in Prey Population. Ecological Complexity, 21, 100-111. https://doi.org/10.1016/j.ecocom.2014.11.005
[3]	Wang, L., Xu, R. and Feng, G. (2016) Modelling and Analysis of an Eco-Epidemiological Model with Time Delay and Stage Structure. Journal of Applied Mathematics and Computing, 50, 175-197. https://doi.org/10.1007/s12190-014-0865-3
[4]	Yang, R., Liu, M. and Zhang, C. (2018) A Diffusive Predator-Prey System with Additional Food and Intra-Specific Competition among Predators. International Journal of Biomathematics, 11, Article ID: 1850060. https://doi.org/10.1142/S1793524518500602
[5]	Mondala, A., Pal, A. and Samanta, G. (2019) On the Dynamics of Evolutionary Leslie-Gower Predator-Prey Eco-Epidemiological Model with Disease in Predator. Ecological Genetics and Genomics, 10, Article ID: 100034. https://doi.org/10.1016/j.egg.2018.11.002
[6]	Shaikh, A., Das, H. and Sarwardi, S. (2020) Dynamics of an Eco-Epidemiological System with Disease in Competitive Prey Species. Journal of Applied Mathematics and Computing, 62, 525-545. https://doi.org/10.1007/s12190-019-01295-6
[7]	Ghanbari, B. (2021) On the Modeling of an Eco-Epidemiological Model Using a New Fractional Operator. Results in Physics, 21, Article ID: 103799. https://doi.org/10.1016/j.rinp.2020.103799
[8]	Holling, C. (1965) The Functional Response of Predator to Prey Density and Its Role in Mimicry and Population Regulation. The Memoirs of the Entomological Society of Canada, 97, 5-60. https://doi.org/10.4039/entm9745fv
[9]	Huang, Y., Shi, W., Wei, C. and Zhang, S. (2021) A Stochastic Predator-Prey Model with Holling II Increasing Function in the Predator. Journal of Biological Dynamics, 15, 1-18. https://doi.org/10.1080/17513758.2020.1859146
[10]	Li, S. and Wang, X. (2015) Analysis of a Stochastic Predator-Prey Model with Disease in the Predator and Beddington-DeAngelis Functional Response. Advances in Difference Equations, 1, Article No. 224. https://doi.org/10.1186/s13662-015-0448-0
[11]	Melese, D. and Feyissa, S. (2021) Stability and Bifurcation Analysis of a Diffusive Modified Leslie-Gower Prey-Predator Model with Prey Infection and Beddington DeAngelis Functional Response. Heliyon, 7, e06193. https://doi.org/10.1016/j.heliyon.2021.e06193
[12]	Arditi, R. and Ginzburg, L. (1989) Coupling in Predator-Prey Dynamics: Ratio-Dependence. Journal of Theoretical Biology, 139, 311-326. https://doi.org/10.1016/S0022-5193(89)80211-5
[13]	Xiao, Y. and Chen, L. (2002) A Ratio-Dependent Predator-Prey Model with Disease in the Prey. Applied Mathematics and Computation, 131, 397-414. https://doi.org/10.1016/S0096-3003(01)00156-4
[14]	Shi, H. and Li, Y. (2015) Global Asymptotic Stability of a Diffusive Predator-Prey Model with Ratio-Dependent Functional Response. Applied Mathematics and Computation, 250, 71-77. https://doi.org/10.1016/j.amc.2014.10.116
[15]	Kim, K., Choi, W. and Ahn, I. (2020) Diffusive Epidemiological Predator-Prey Models with Ratio-Dependent Functional Responses and Nonlinear Incidence Rate. Journal of Dynamical and Control Systems, 1-15. https://doi.org/10.1007/s10883-020-09512-3
[16]	Mukherjee, D. (2003) Stability Analysis of a Stochastic Model for Prey-Predator System with Disease in the Prey. Nonlinear Anal: Model. Control, 8, 83-92. https://doi.org/10.15388/NA.2003.8.2.15186
[17]	Li, D., Cui, J., Liu, M. and Liu, S. (2015) The Evolutionary Dynamics of Stochastic Epidemic Model with Nonlinear Incidence Rate. Bulletin of Mathematical Biology, 77, 1705-1743. https://doi.org/10.1007/s11538-015-0101-9
[18]	Meng, X., Zhao, S. and Feng, T. (2016) Dynamic of a Novel Nonlinear Stochastic SIS Epidemic Model with Double Epidemic Hypothesis. Journal of Mathematical Analysis and Applications, 433, 227-242. https://doi.org/10.1016/j.jmaa.2015.07.056
[19]	Bai, C., Deng, M. and Du, B. (2016) Analysis of Stochastic Two-Prey One-Predator Model with Lvy Jumps. Physica A: Statistical Mechanics and its Applications, 445, 176-188. https://doi.org/10.1016/j.physa.2015.10.066
[20]	Liu, Q. and Jiang, D. (2019) Stationary Distribution and Extinction of a Stochastic One-Prey Two-Predator Model with Holling Type II Functional Response. Stochastic Analysis and Applications, 37, 321-345. https://doi.org/10.1080/07362994.2019.1566005
[21]	Liu, S., Ruan, S. and Zhang, X. (2017) Nonlinear Dynamics of Avian Influenza Epidemic Models. Mathematical Biosciences, 283, 118-135. https://doi.org/10.1016/j.mbs.2016.11.014
[22]	Liu, M. and Wang, K. (2011) Survival Analysis of a Stochastic Cooperation System in a Polluted Environment. Journal of Biological Systems, 19, 183-204. https://doi.org/10.1142/S0218339011003877
[23]	Khasminskii, R. (2012) Stochastic Stability of Differential Equations. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-23280-0
[24]	Higham, D. (2001) An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations. SIAM Review, 43, 525-546. https://doi.org/10.1137/S0036144500378302