A mathematical study of effects of delays arising from the interaction of anti-drug antibody and therapeutic protein in the immune response system

1.   Introduction

In the field of biomathematics, the predator-prey model has been studied by many scholars. They explored the dynamical behavior among biological populations by establishing differential equations ^[1,2,3,4]. In 1925, the Lotka-Volterra model described the variation in population size between predator and prey in ^[5,6]. This model describes that population size changes are interacting. After that, in order to study the situation that predator populations have other food sources besides prey population, Leslie and Gower ^[7,8] proposed a Leslie-Gower type predator-prey model. This model mainly explores the fact that when the preferred food decreases, the number of predators also decreases. Then Aziz-Alaoui and Okiye added a constant to the denominator of the Leslie-Gower type functional response function in ^[9], calling it the modified Leslie-Gower type. This type limits the growth of predators due to severe shortages of preferred foods, despite the predators having other food sources. Besides, the impact of infectious diseases on populations in nature is also important. In order to study the dynamical behavior of populations in more complex situations, some factors influencing population size changes have been added to the Leslie-Gower type model. The authors of ^[10] describe the mechanism of disease transmission by using the Holling type Ⅱ. Some scholars have included disease factors in their studies of fractional differential equations as well, e.g., ^{[11,12,13,14]}.

However, in nature, biological populations are inevitably impacted by environment noise more or less. During the past decades, many investigators have focused on the study of stochastic biological models ^[15,16]. Among them, a predator-prey-parasite model with stochastic perturbations has been studied by Majumder et al. ^[17]. Parasitic infections divide prey populations into susceptible and infected populations, and infected populations lose fertility and do not heal again. Both susceptible and infected prey populations are preyed upon by predators, and the predators will not be infected by the disease. They constructed the following model:

where $x(0)\geq0$, $y(0)\geq0$, $z(0)\geq0$.

In addition, the fear of predators can also influence the birth rates and offspring survival of prey populations. Zanette et al. ^[18] verified this idea through experiments. It is noted in ^[19] that the mental state of juvenile prey can be mediated by predator-induced fear and this fear may have an impact on their survival rates as adults. Thus, many scholars have realized that the fear costs of predators can directly or indirectly affect the prey population. So it should be included in the predator-prey system. Based on this view, scholars have studied the fear effect of biological populations ^[20,21]. Qi and Meng ^[22] used $\left(\theta+\frac{K(1-\theta)}{K+z(t)}\right)$ to represent the fear function to measure the cost of fear. $\theta\in[0, 1]$ represents the cost of minimum fear and $K$ represents the level of prey fear of predator populations. Let $\kappa(\theta, K, z) = \theta+\frac{K(1-\theta)}{K+z}$; we have $\frac{\partial\kappa(\theta, K, z)}{\partial z} < 0$. From this, it is clear that the larger the predator population, the stronger the inhibitory effect on the growth of the prey population. Therefore, introducing a fear factor into the prey-predator system can help us to explore the variation of populations in different situations better.

In addition, natural species may be subjected to unexpected environmental disruptions such as epidemics, hurricanes, earthquakes and so on. Random perturbations described by Brownian can only characterize the continuous influence, but it does not describe sudden and drastic environmental changes very well. To explain this occurrence, Bao et al. included Lévy jumps into population models in ^[23,24]. Wu and Wang ^[25] considered the population dynamical behaviors of stochastic system with jumps. So we also want to introduce Lévy jumps into a stochastic predator-prey model.

Inspired by the articles above, we consider adding a fear factor to susceptible prey populations and considered more complex diseases, using the Holling type Ⅱ function to represent the spread of disease. Finally, we use Lévy noise to describe the situation when the population is subjected to drastic changes from the outside. We consider the following stochastic disease including predator-prey model with Lévy noise and fear effects:

with initial values $x(0) = x_{0}$, $y(0) = y_{0}$ and $z(0) = z_{0}$. In this model, $x(t)$, $y(t)$ and $z(t)$ denote the population densities of the susceptible prey, infected prey and predator population at time $t$ respectively. And the per capita maximum fertility rates of the prey and predator populations are written as $r$ and $\beta$ respectively; the intensity of interspecific competition in prey populations is denoted by $b$ and $c$; $\eta$ represents the disease transmission rate and $\frac{\eta x(t)y(t)}{b_{3}+x(t)}$ describes the spread of the disease; $a$ and $d$ are the consumption rates and $\gamma$ is the death rate of the infected prey population. $g$ represents the interspecific competition of predators and $m$ is the half-saturation constant of predators. $b_{1}$ and $b_{2}$ represent the prey saturation constants and predator disturbance respectively. There are no negative constants among any of the parameter values.

$B_{i}(t) \ (i = 1, 2, 3)$ represents Brownian motion, and each value is independent of each other. In addition, there is a complete probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_{t}\}_{t\geq0}, P)$ with a filtration $\{\mathcal{F}_{t}\}_{t\geq0}$ that meets the normal requirements, and $B_{i}(t) \ (i = 1, 2, 3)$ is defined on this probability space. $\sigma_{i}^{2} \ \ (i = 1, 2, 3)$ represents the noise's level of intensity. The left limits of $x(t)$, $y(t)$ and $z(t)$ are represented by $x(t^{-})$, $y(t^{-})$ and $z(t^{-})$ respectively. $N({\rm d}t, {\rm d}u)$ is a Poisson counting measure which is defined on $\lambda({\rm d}u)$. The characteristic measure $\lambda$ on the measure subset $\Gamma$ of $[0, +\infty]$ such that $\lambda(\Gamma) < \infty$. $N({\rm d}t, {\rm d}u)$ is defined on ${\bf{R_{+}}}\times({\bf{R}}-\{0\})$, ${\bf{R}}_{+}: = (0, \infty)$. Besides, $\widetilde{N}(dt, du) = N(dt, du)-\lambda(du)dt$ is the corresponding martingale measure. $\gamma_{i}(u) \ (i = 1, 2, 3)$ measures the effect of a Lévy jump on prey and predator populations, $\gamma_{i}(u) > -1$ $(i = 1, 2, 3)$ for $u\in\Gamma$. In addition, it is important to note that Lévy jumps have a facilitating effect on the ecosystem when $\gamma_{i} > 0 \ (i = 1, 2, 3)$, such as an ocean red tide. When $\gamma_{i} < 0 \ (i = 1, 2, 3)$, Lévy jumps have a negative effect on the ecosystem, such as tsunamis and earthquakes. See ^[26,27] for specific examples.

In this article, we need to assume that the coefficients satisfy the following assumption:

Assumption 1. There exists a positive constant $K$ which gives:

It means that the intensity of Lévy noise is not very large. Here are some inequalities we will frequently use:

Lemma 1. (^[28]) Denote by $\Sigma(t)$ a local martingle vanishing at $t = 0$. Define

where $\langle\Sigma\rangle(t) = \langle\Sigma, \Sigma\rangle(t)$ stands for the Meyer's angle bracket process. If $\limsup\limits_{t\rightarrow +\infty}\overline{\Sigma}(t) < +\infty$, then $\lim\limits_{t\rightarrow +\infty}t^{-1}\Sigma(t) = 0$, a.s.

The rest of the research for this paper is as follows. In Section 2, we present the existence and uniqueness of the positive solution of System (1.1). In Section 3, we first studied the conditions for population extinction. Then we considered the existence and extinction of susceptible prey populations and predator populations under conditions where the disease is effectively prevented and infected prey populations are extinct. In Section 4, it proves the stochastic ultimate boundedness of System (1.1) and the existence of ergodic stationary distribution of the System (1.1) when Lévy noise does not exist. In Section 5, we select suitable parameters and use numerical simulations to prove our conclusions. Lastly, we briefly summarize the work of article.

2.   Existence and uniqueness of a global positive solution

In order to study the dynamical behavior of the system, we first verify that there is a globally unique positive solution for System (1.1). First, we give Lemma 2 to show that the positive solution of the system exists locally and uniquely, and then prove that the solution exists globally with Theorem 1.

Lemma 2. For any given initial value $(x_{0}, y_{0}, z_{0})\in{{\bf{R}}_{+}^{3}}$, there exists a unique local positive solution $(x(t), y(t), z(t))$ to System (1.1), as defined on the interval $t\in[0, \tau_{e})$, where $\tau_{e}$ is the explosion time.

Proof. Consider the equation

with initial values $N_{1}(0) = \ln{x_{0}}$, $N_{2}(0) = \ln{y_{0}}$ and $N_{3}(0) = \ln{z_{0}}$ on $t\geq0$. It is easy to see that the above equation satisfies the local Lipschitz condition. Therefore, System (1.1) has a unique local solution $(N_{1}(t), N_{2}(t), N_{3}(t))$ for $t\in[0, \tau_{e})$, and $\tau_{e}$ is the explosion time. Because $x(t) = e^{N_{1}(t)}$, $y(t) = e^{N_{2}(t)}$ and $z(t) = e^{N_{3}(t)}$, by Itô's formula, we get that $(x(t), y(t), z(t))$ is the unique local positive solution to System (1.1) with the initial value $(x_{0}, y_{0}, z_{0})$.

Theorem 1. For any given initial value $(x_{0}, y_{0}, z_{0})\in{{\bf{R}}_{+}^{3}}$, a unique global positive solution $(x(t), y(t), z(t))$ exists in System (1.1) for $\forall t\in[0, +\infty)$.

Proof. In order to prove that the solution is global, we need to prove $\tau_{e} = \infty$ a.s. Let $m_{0} > 0$ be sufficiently large to make $x_{0}, y_{0}, z_{0}\in[\frac{1}{m_{0}}, m_{0}]$ for each integer $m\geq m_{0}$; then, define the stopping time

Let $\inf\emptyset = \infty$ ($\emptyset$ refers to the empty set). Obviously, as $m$ tends to infinity, $\tau_{m}$ is increasing. Let $\tau_{\infty} = \lim\limits_{m\rightarrow \infty}\tau_{m}$ and $\tau_{\infty}\leq\tau_{e}$, a.s. If we can prove $\tau_{\infty} = \infty \ \ a.s.$ then we have $\tau_{e} = \infty \ \ a.s.$, for all $t\geq{0}$.

We can proof by contradiction. If not, there exists a pair constants $T\geq0$ and $\varepsilon\in(0, 1)$ such that $P(\tau_{\infty}\leq{T}) > \varepsilon$. Thus $\exists \ m_{1}\geq{m_{0}}$, we have

Define a $C^{2}$-function $V$: ${\bf{R}}_{+}^{3}\rightarrow {{\bf{R_{+}}}}$

Because $s-1-\ln{s} > 0$, for all $s > 0$, we have $V(x, y, z) > 0$. Applying Itô's formula, the following equation yields

$\mathcal{L}V:{\bf{R}}_{+}^{3}\rightarrow {{\bf{R_{+}}}}$ is given as follows and using Assumption 1, we obtain

where

is a positive constant.

We have that $z\leq{2(z-1-\ln{z})+\ln{4}}\leq{2V(x, y, z)+\ln{4}}$. We can write:

where $M_{1} = \max\left\{M+(\beta+g)\ln4, 2(\beta+g)\right\}$.

Combining (2.1) and (2.3), we can obtain that

The integration is taken at both ends of the inequality (2.4) from $0$ to $\tau_{m}\wedge{T}$, followed by the expectation, yielding

By Gronwall's inequality, we can get

Let $\Omega_{m} = {\tau_{m}\leq{T}}$; we have $P(\Omega_{m})\geq\varepsilon$. So for $\forall\omega\in\Omega_{m}$, at least one value of $x(\tau_{m}, \omega)$, $y(\tau_{m}, \omega)$ or $z(\tau_{m}, \omega)$ equals either $m$ or $\frac{1}{m}$. Note that $V(x(\tau_{m}), y(\tau_{m}), z(\tau_{m}))$ is no less than $(m-1-\ln{m})\wedge(\frac{1}{m}-1-\ln\frac{1}{m})$. Consequently,

where $1_{\Omega_{m}(\omega)}$ is the indicator function of $\omega_{m}$. Then, let $m\rightarrow \infty$; we deduce that

this is contradictory. Hence, we can get $\tau_{\infty} = \infty$.

The result is confirmed.

3.   Long time behavior of System (1.1)

In this section, we consider the long time behavior of System (1.1). The conditions when the susceptible prey population, the infected prey population and the predator population are all extinct are first considered. Then we explore the existence and extinction of susceptible prey populations and predator populations in the context of the effective prevention of infectious disease.

3.1.   Extinction

Theorem 2. For any given initial value $(x_{0}, y_{0}, z_{0})\in{{\bf{R}}_{+}^{3}}$, the solution $(x(t), y(t), z(t))$ of system (1.1) has the following properties if Assumption 1 holds:

then the predator and prey populations will be extinctive almost surely, that is

Proof. First of all, we consider the prey population. We have ${\rm d}\{e^{t}\ln{x(t)}\} = e^{t}{\rm d}\ln{x(t)}+e^{t}\ln x(t)$ and the fundamental inequality $\ln{x}\leq{x-1}$ for all $x > 0$. Calculating by Itô's formula, we get

Integrating both sides of the above inequality simultaneously, we have

Then, we have

Denote $\Sigma_{1}(t) = \int_{0}^{t}\int_{\Gamma}\ln(1+\gamma_{1}(u))\widetilde{N}({\rm d}s, {\rm d}u)$; in light of Assumption 1,

where $F$ is a positive number. So we have $\int_{0}^{t}\frac{F}{(1+s)^{2}}{\rm d}s = \frac{t}{t+1} < \infty$; then, it follows from Lemma 2 that

Taking the limit on both sides of the inequality (3.2) and bringing in (3.3), we get

When $r(\theta+K(1-\theta))-\frac{\sigma_{1}^{2}}{2}-\int_{\Gamma}[\gamma_{1}(u)-\ln(1+\gamma_{1}(u))]\lambda({\rm d}u) < 0$, the susceptible prey population will be extinct.

Similarly, for the infected prey and predator populations, we have

When $\eta-\frac{\sigma_{2}^{2}}{2}-\int_{\Gamma}[\gamma_{2}(u)-\ln(1+\gamma_{2}(u))]\lambda({\rm d}u) < 0$ and $\beta-\frac{\sigma_{3}^{2}}{2}-\int_{\Gamma}[\gamma_{3}(u)-\ln(1+\gamma_{3}(u))]\lambda({\rm d}u) < 0$, the infected prey and predator population will be extinct.

The result is confirmed.

3.2.   Behavior of System (1.1) when y(t) is extinct

In this section, we expect that infectious diseases transmitted among prey populations will be effectively prevented, susceptible prey populations will no longer be infected and infected prey populations will gradually die out. Considering the extinction of the infected population under the condition that $\eta-\frac{\sigma_{2}^{2}}{2}-\int_{\Gamma}[\gamma_{2}(u)-\ln(1+\gamma_{2}(u))]\lambda({\rm d}u) < 0$, i.e., after the infectious disease is cured, we have $\lim\limits_{t\rightarrow \infty}y(t) = 0$. For $\forall \varepsilon > 0$, there exist $t_{1}$ and a set $\Omega_{\varepsilon}$ such that $P(\Omega_{\varepsilon})\geq1-\varepsilon$ and $\frac{\eta x(t)y(t)}{b_{3}+x(t)} < \varepsilon$; for $t_{1}\leq t$ and $\omega\in\Omega_{\varepsilon}$. Therefore, we next focused on the changes in susceptible and predator populations when infected prey populations perished.

Definition 1. (^[29]) The susceptible prey populations and predator populations are said to be persistent in mean if

Lemma 3. (^[30]) Let $b(t)\in{\bf{C}}(\Omega\times[0, +\infty), R_{+})$

1) If there exist two positive constants $T$ and $v_{0}$ such that $\ln b(t)\leq vt-v_{0}\int_{0}^{t}b(s){\rm d}s+\sigma B(t)$ for $\forall t\geq T$, where $\sigma > 0$, then

2) If there exist two positive constants $T$ and $v_{0}$ such that $\ln b(t)\geq vt-v_{0}\int_{0}^{t}b(s){\rm d}s+\sigma B(t)$ for $\forall t\geq T$, where $\sigma > 0$, then

Theorem 3. Suppose that $(x(t), y(t), z(t))$ denotes the positive solution to System (1.1) with the initial positive value $(x_{0}, y_{0}, z_{0}) > 0$; when infected prey populations tend to become extinct, that is, $\lim\limits_{t\rightarrow \infty}y(t) = 0$, we have the following

$({\bf{A1}}).$ If $r < \frac{\sigma_{1}^{2}}{2}-\int_{\Gamma}[\gamma_{1}(u)-\ln(1+\gamma_{1}(u))]\lambda({\rm d}u)$, $\beta > \frac{\sigma_{3}^{2}}{2}-\int_{\Gamma}[\gamma_{3}(u)-\ln(1+\gamma_{3}(u))]\lambda({\rm d}u)$, then the predator is persistent in mean and the susceptible prey is extinct, that is

$({\bf{A2}}).$ If $\frac{1}{b}\left(\frac{r\theta}{2}-\frac{a}{b_{2}}-\frac{\sigma_{1}^{2}}{2}-\int_{\Gamma}[\gamma_{1}(u)-\ln(1+\gamma_{1}(u))]\lambda({\rm d}u)\right) > 0$ and $\beta < \frac{\sigma_{3}^{2}}{2}-\int_{\Gamma}[\gamma_{3}(u)-\ln(1+\gamma_{3}(u))]\lambda({\rm d}u)$, then the susceptible prey is persistent in mean and the predator is extinct, that is

$({\bf{A3}}).$ If $\frac{1}{b}\left(\frac{r\theta}{2}-\frac{a}{b_{2}}-\frac{\sigma_{1}^{2}}{2}-\int_{\Gamma}[\gamma_{1}(u)-\ln(1+\gamma_{1}(u))]\lambda({\rm d}u)\right) > 0$ and $\beta > \frac{\sigma_{3}^{2}}{2}-\int_{\Gamma}[\gamma_{3}(u)-\ln(1+\gamma_{3}(u))]\lambda({\rm d}u)$, then both the susceptible prey and predator are persistent in mean, that is

Proof. According to Theorem 2, when $\eta-\frac{\sigma_{2}^{2}}{2}-\int_{\Gamma}[\gamma_{2}(u)-\ln(1+\gamma_{2}(u))]\lambda({\rm d}u) < 0$, $\lim\limits_{t\rightarrow +\infty}y(t) = 0$. So we know that $\exists \varepsilon_{2} > 0$ for $\forall T_{2} > 0$ when $t > T_{2}$ such that $0 < y(t)\leq\varepsilon_{2}$.

${\bf{(A1)}}$ Similarly, when $\lim\limits_{t\rightarrow +\infty}x(t) = 0$, we know that $\exists \varepsilon_{1} > 0$ for $\forall T_{1} > 0$ when $t > T_{1}$ such that $0 < x(t)\leq\varepsilon_{1}$. We obtain

Integrating both sides of the above formula, we have

Similarly, we have

Applying Lemma 3 and Assumption 1 to (3.4) and (3.5) respectively, we have

For $\forall \varepsilon_{1}, \varepsilon_{2}$, we can obtain

${\bf{(A2)}}$ When $\beta-\frac{\sigma_{3}^{2}}{2}-\int_{\Gamma}[\gamma_{3}(u)-\ln(1+\gamma_{3}(u))]\lambda({\rm d}u) < 0$, $z(t)$ comes to extinct, $\lim\limits_{t\rightarrow \infty}z(t) = 0$. By Itô's formula, we have

and

Integrating both sides of (3.8) and (3.9) from $0$ to $t$ and letting $\eta \varepsilon_{2} = \frac{r\theta}{2}$, we obtain

and

Similar to (A1), by Lemma 3 and Assumption 1, we can get

${\bf{(A3)}}$ From (3.5), we can deduce that

The conclusion is confirmed.

4.   Stochastic ultimate boundedness and stationary distribution without Lévy noise

Moreover, when $\gamma_{i}(u) = 0$ ($i = 1, 2, 3$), this means that the population will not suffer drastic environmental changes. So, System (1.1) produces the following system:

4.1.   Stochastic ultimate boundedness

In this part, we solve the stochastically ultimately bounded problem for the system solution. Before the proof we need preparation.

Definition 2. (^[31]) The solution of System (1.1) is called stochastically ultimately bounded, for any $\varepsilon\in(0, 1)$ if there exists a constant $H = H(\varepsilon)$ such that for any initial value $W_{0} = (x_{0}, y_{0}, z_{0})$ in ${\bf{R}}_{+}^{3}$, the solution $W(t) = (x(t), y(t), z(t))$ of System (1.1) has the property that

Theorem 4. The solution of System (1.1) is stochastically ultimately bounded for any initial value $W_{0} = (x_{0}, y_{0}, z_{0})$ in ${\bf{R}}_{+}^{3}$.

Proof. First, define a function $V:{\bf{R}}_{+}^{3}\rightarrow {\bf{R}}_{+}$

By Itô's formula, we can get

So we have

and

Rectifying (4.2) and (4.3) yields

Define the function

we can obtain

From the above equation we can obtain that there exists a positive number $G$ such that $\mathcal{L}W\leq Ge^{t}$. Therefore integrating both sides of (4.5) from $0$ to $t$ gives

Then, the expectations are taken at both ends of (4.6), so the following results can be obtained

Therefore, we have

Using Chebyshev inequality, we obtain

Next, taking the upper limit of (4.7) gives

where $\varepsilon\in(0, 1)$ and $H = \sqrt[2]{\frac{{2G_{1}}}{{\varepsilon}}}$.

The conclusion is confirmed.

4.2.   Existence of ergodic stationary distribution of System (4.1).

Let $X(t)$ be a homogeneous Markov process in $E_{h}$ ($E_{h}$ denotes Euclidean $h$-space) described by the stochastic equation

The diffusion matrix is given as follows:

If there exists a bounded domain $U\subset{E_{h}}$ with a regular boundary, then the following lemma holds:

Lemma 4. (^[32]) The Markov process $X(t)$ has a unique stationary distribution $\varsigma(\cdot)$ if it satisfies the following conditions:

(A.1):Suppose a positive number $M$ makes $\sum_{i, j = 1}^{d}a_{ij}(x)\zeta_{i}\zeta_{j}\geq{M|\zeta|^{2}}$, $x\in{U}$, $\xi\in{\bf{R}}^{d}$.

(A.2):There exists a $C^{2}-function$ such that $\mathcal{L}V$ is negative for $\forall{{\bf{R}}_{+}^{3}}\setminus{U}$. Then we have

for $\forall x\in{E_{h}}$, where $f(\cdot)$ is a function integrable with respect to the measure $\varsigma$.

To verify the condition (A.2), it is necessary to prove that there exist a neighborhood $U$ and a nonnegative function $V(x, y, z)$ such that $\mathcal{L}V$ is negative for any $E_{h}\backslash{U}$.

Lemma 5. (^[33]) For $\forall s > 0$, the following inequality holds:

Theorem 5. According to Lemma 4, for any initial value $(x_{0}, y_{0}, z_{0})$, there exists an ergodic stationary distribution for System (4.1) if the following conditions hold.

Proof. Now we prove the condition (A.2). According to the inequality $b(x-\frac{r}{b})^{2} > 0$, we derive

Define $V_{11}(x, y) = -\ln y+Ax$, where $A$ is a positive constant which will be determined later; we have

Then define $V_{12} = \frac{b}{2r^{2}}x+\frac{b}{r(b+1)}(-\ln x)$ and use Lemma 5; we can obtain

Next, we define $V_{13}(x, y) = V_{11}(x, y)+2\sqrt{\frac{A\eta r^{2}}{b}}V_{12}(x, y)+\frac{1}{\gamma}\left(c+\frac{2bn}{r(b+1)}\sqrt{\frac{A\eta r^{2}}{b}}\right)y$. Combining (4.9) and (4.10), we have the following inequality

Choose

So we have

where $B = \frac{b^{2}\eta \left(r-\frac{\sigma_{1}^{2}}{2}-\frac{a}{b_{2}}\right)^{2}}{r(b_{3}b+1)(b+1)^{2}}-(\gamma+\frac{d}{b_{2}}+\frac{\sigma_{2}^{2}}{2}) > 0$.

In the next step, we make

and $i\in(0, 1)$ is a positive number. By Itô's formula, we have

Let us choose $i$ sufficiently small such that

Finally, we define the function $V_{15} = \mathcal{R}V_{13}+V_{14}+3\mathcal{R}x+(m+x+y)z^{2}$, where $\mathcal{R}$ is a positive constant satisfying $-\mathcal{R}B+f_{1}^{u}+f_{2}^{u}+f_{3}^{u} = -2$, in which $f_{1}^{u}, f_{2}^{u} \ and \ f_{3}^{u}$ are bounded functions on $[0, \infty)$ and $f_{i}^{u} = \sup_{t\geq0}f_{i}(t) \ (i = 1, 2, 3)$; it will be determined later.

According to (4.3), we can obtain

where

With the condition that $\frac{c^{2}}{\gamma}-\frac{\eta}{2\gamma}\left(c+\frac{2b\eta}{r(b+1)}\sqrt{\frac{A\eta r^{2}}{b}}\right) > 0$, $2b-\frac{\eta}{2\gamma}\left(c+\frac{2b\eta}{r(b+1)}\sqrt{\frac{A\eta r^{2}}{b}}\right) > 0$. Denote

Then we have

Therefore, we can deduce that

for $\forall(x, y, z)\in {{\bf{R}}_{+}^{3}\setminus{U}}$, which implies that the condition (A.2) in Lemma 4 is satisfied.

The next step will be to prove the condition (A.1) and the following is the diffusion matrix of the System (4.1):

Choose $\widetilde{M} = \min_{x, y, z\in{U\subset{\bf{R}}_{3}^{+}}}\{\sigma_{1}^{2}x^{2}, \sigma_{2}^{2}y^{2}, \sigma_{3}^{2}z^{2}\}$ such that

Therefore we can conclude that the condition (A.1) in Lemma 4 holds. Further, from Lemma 4, we can infer that System (4.1) is ergodic and has a unique stationary distribution.

5.   Illustrative examples

In this section, to verify the conditions obtained by theorems, we take the determined initial values $x_{0} = 2.9$, $y_{0} = 1.4$ and $z_{0} = 0.5$ for the numerical simulation of System (1.1). In addition, let $\sigma_{1} = \sigma_{2} = \sigma_{3} = 0.6$ and $\gamma_{1} = \gamma_{2} = \gamma_{3} = 0.06$. The figures for the numerical simulations are as follows, where the left figure shows the numerical simulation of the stochastic model with white noise and Lévy noise, and the right figure shows the numerical simulation of the deterministic model.

Example 1.

As shown on the left in Figure 1, in order to verify the case of extinction of both the predator and prey populations in Theorem 2, we chose the appropriate parameters $r = 0.1, \theta = 0.09, K = 0.35, b = 0.8, a = 0.27, b_{1} = 0.16, b_{2} = 0.5,$ $\eta = 0.015, b_{3} = 0.08, c = 0.08, d = 0.004, \gamma = 0.06, \beta = 0.013, g = 0.01, m = 0.01$, $\sigma_{1} = \sigma_{2} = \sigma_{3} = 0.6$ and $\Gamma_{1} = \Gamma_{2} = \Gamma_{3} = 0.06$; then, we have $0.1\times(0.09+0.35\times0.01)-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx-0.172 < 0$, $0.015-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx-0.003 < 0$ and $0.13-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx-0.168 < 0$. In this case, the populations of $x(t), y(t)$ and $z(t)$ all tend to become extinct, in accordance with the conclusion obtained in Theorem 2. Compared to the graph on the right, the curve of the stochastic model converges to zero with sharp fluctuations, while the curve of the deterministic model is smooth and takes longer to converge to zero. This shows that random factors accelerate the rate of population extinction when the conditions of Theorem 2 hold.

Example 2.

When $0.015-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx-0.003 < 0$, the infected prey population tends to become extinct. This is illustrated in Figures 2–4.

In order to verify the conditions of (A1) in Theorem 3, the numerical simulation we made by selecting suitable parameters is shown in Figure 2. Figure 2 shows that when $r = 0.1, \theta = 0.09, K = 0.35, b = 0.8, a = 0.27, b_{1} = 0.16, b_{2} = 0.5, \eta = 0.015, b_{3} = 0.08,$ $c = 0.08, d = 0.004, \gamma = 0.06, \beta = 0.38, g = 0.01, m = 0.01$, $\sigma_{1} = \sigma_{2} = \sigma_{3} = 0.6$ and $\Gamma_{1} = \Gamma_{2} = \Gamma_{3} = 0.06$; then, we have $0.1-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx-0.08 < 0$, $0.015-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx-0.003 < 0$ and $0.38-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u)\approx0.199 > 0$. At this time, susceptible and infected prey populations $x(t)$ and $y(t)$ tend to become extinct and predator populations $y(t)$ tend to persist. Compared to the deterministic model, the left curve with the effect of white noise and Levy noise changes more dramatically.

In order to verify the conditions of (A2) in Theorem 3, the numerical simulation we made by selecting suitable parameters is shown in Figure 3. We chose $r = 21.45$ and $\beta = 0.01$ and the other parameters take the same value. Then we have $21.45-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u) > 0$ and $0.01-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u) < 0$. At this point, as seen in Figure 3, susceptible prey populations tend to persist and predator populations tend to become extinct. Whereas in the right picture, $x(t)$ is gradually leveling off and $z(t)$ is not extinct, but slowly rising.

In order to verify the conditions of (A3) in Theorem 3, the numerical simulation we made by selecting suitable parameters is shown in Figure 4. We set $r = 21.45$ and $\beta = 0.19$ and the other parameters take the same value. Then we have $21.45-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u) > 0$ and $0.19-\frac{0.6^{2}}{2}-\int_{\Gamma}[0.06-\ln(1+0.06)]\lambda({\rm d}u) > 0$. From Figure 3, we can see that both the susceptible prey populations and predator populations tend to be persistent. This satisfies the condition given by (A3) in Theorem 3. z(t) grows more rapidly in the deterministic model than in the model with white noise versus Lévy noise.

Example 3.

To verify the condition of Theorem 5 that System (1.1) has an ergodic stationary distribution, we chose the following parameter values: $r = 18, \theta = 0.09, K = 0.35, b = 0.8, a = 0.27, b_{1} = 0.16, b_{2} = 4.6, \eta = 0.78,$ $b_{3} = 0.18, c = 0.7, d = 0.004, \gamma = 0.19, \beta = 0.65, g = 0.012$ and $m = 0.012$. From the conditions obtained in Theorem 5, it follows that $\frac{\eta}{2\gamma}\left(c+\frac{2b\eta}{r(b+1)}\sqrt{\frac{A\eta r^{2}}{b}}\right)\approx1.863 < \frac{c^{2}}{\gamma}\approx2.579$, $\frac{\eta}{2\gamma}\left(c+\frac{2b\eta}{r(b+1)}\sqrt{\frac{A\eta r^{2}}{b}}\right)\approx1.863 < 3b = 2.4$, $\frac{b^{2}\eta \left(r-\frac{\sigma_{1}^{2}}{2}-\frac{a}{b_{2}}\right)^{2}}{r(b_{3}b+1)(b+1)^{2}}\approx2.360 > \gamma+\frac{d}{b_{2}}+\frac{\sigma_{2}^{2}}{2} = 0.37$, $\frac{2d}{b_{2}}+2\gamma+\sigma_{2}^{2}\approx0.741 < \eta = 0.78$ and $\frac{\sigma_{3}^{2}}{2} = 0.18 < 0.65$. As shown in Figure 5, the sample paths are concentrated within regions of circularity or ellipticity, which indicates that the system is stochastically stable.

6.   Conclusions

This paper discusses the dynamical properties of predator-prey models with fear effects and disease transmission in prey population. Meanwhile, susceptible prey populations, infected prey populations and predator populations are affected by white noise and Lévy noise. By constructing the appropriate Lyapunov equation, we proved the uniqueness of the global positive solution of System (1.1). From Theorem 2, it is clear that under certain conditions, Lévy noise can lead to population extinction. Furthermore, the fear effect can also lead to population size change. When the fear effect is too large, it is more likely to result in population extinction. We also explored the existence and extinction of susceptible prey populations and predator populations under conditions when infectious diseases prevalent among prey populations are effectively prevented and infected prey populations die out in Theorem 3. Then, we studied the stochastic ultimate boundedness of System (4.1) and the ergodic stationary distribution under certain conditions without the influence of Lévy noise. Finally, the numerical simulations were performed at the end to further illustrate the validity of the theoretical results. Adding the influence of environmental factors to the model made it more consistent with the predator-prey relationship in the ecosystem. This shows that stochastic factors have an effect on the behavior of population dynamics and in some cases can accelerate the extinction of populations. The fear effect also affects the population size change; when the fear effect is stronger, the population is more likely to become extinct.

Acknowledgments

This research was supported by the Natural Science Foundation of Heilongjiang Province (No. LH2022A002, No. LH2021A018), National Natural Science Foundation of China (No. 12071115) and Fundamental Research Funds for the Central Universities (No. 2572020BC09).

Conflict of interest

The authors declare that there is no conflict of interest.