Extinction and stationary distribution of stochastic predator-prey model with group defense behavior

1.   Introduction

In ecosystems, predation is not the only relationship between predator and prey. Many preys can perceive threats from predators, and thus make a variety of different response behaviors ^[1], which has an impact on the population dynamics of predators and prey. Studies have shown that the impact of this indirect behavior can be as great as the direct capture and killing their prey by predators ^[2,3,4,5,6]. For example, sparrowhawks are more likely to attack individual redshanks or small groups of redshanks during hunting, while their predation success rate is significantly lower for larger groups of redshanks ^[7]. Swarming locusts release large amounts of the volatile compound phenylacetonitrile, while dispersing locusts hardly synthesize phenylacetonitrile. When locust is threatened by natural enemies, phenylacetonitrile, as an olfactory warning compound, can further synthesize highly toxic substance hydrocyanic acid to achieve the purpose of defense against natural enemies ^[8]. As a result, locusts can reproduce recklessly, which causes locust plagues and results in serious economic losses ^[9,10]. This behavior of prey is called group defense, which reduces the defense cost for individuals through 'number security' and thus increases survival rates. Freedman and Wolkowicz first investigated predator-prey interactions under the effects of prey group defense behavior and concluded that group defense behavior of sufficiently abundant prey would lead to extinction of predator populations ^[11]. From then on, more and more scholars have studied the more complex interspecific relationships of predator-prey systems under the effects of group defense ^{[12,13,14,15,16,17]}. Xiao and Ruan analyzed the dynamical complexities of the following predator-prey model with group defense behavior in ^[18]:

where $x(t)$ and $y(t)$ represent the population densities of prey and predator at time $t$, respectively. All parameter values of model (1.1) are non-negative. $r$ denotes the intrinsic growth rate of prey population; $K$ is the environmental capacity; $\alpha$ is the predation intensity; $\mu$ is the conversion efficiency of biomass; $D$ is the natural death rate of predator. The response function $f(x) = \alpha x{e^{ - \beta x}}$ represents the capture rate of predators on preys with group defense behavior, which satisfies

and

where $\beta \in \left[{0, 1} \right)$ reflects the group defense level of prey population. The predation rate of predator reaches the maximum value at $x = \frac{1}{\beta }$, and decreases if the group defense level of prey increases when $x > \frac{1}{\beta }$. The larger $\beta$ indicates that prey population can produce group defense behavior under low density. From ^[18], we know model (1.1) has complicated dynamics, when $\mu \alpha > e\beta D$, if ${x_1} < K < {1 / \beta }$, model (1.1) has three equilibria: two hyperbolic saddles $\left({0, 0} \right)$ and $\left({K, 0} \right)$, and a globally asymptotically stable equilibrium $\left({{x_1}, {y_1}} \right)$ in the interior of the first quadrant; if $K > {1 / \beta } + {x_2}$, model (1.1) has four equilibria: two hyperbolic saddles $\left({0, 0} \right)$ and $\left({{x_2}, {y_2}} \right)$, a hyperbolic stable node $\left({K, 0} \right)$ and an unstable equilibrium $\left({{x_1}, {y_1}} \right)$, and in this case, model (1.1) has no closed orbits, where $x_1$ and $x_2$ are roots of $\mu \alpha x{e^{ - \beta x}} - D = 0$. Moreover, model (1.1) may have a unique limit cycle or a homoclinic loop for different parameter values.

However, considering that in nature, various systems are inevitably affected by some uncertain environmental noises (weather change, human activities, etc.), which cannot be ignored, and it is more realistic and appropriate to use stochastic dynamical systems to describe these phenomena. Recently, many scholars considered introducing the stochastic pertubations into dynamical models ^{[19,20,21,22,23]} and have achieved meaningful results, such as the prediction and control of infectious disease ^{[24,25,26,27]}, pest management ^[28] and various predator-prey systems ^{[29,30,31,32]}. Based on ^[18] and above biological background, considering that small environmental fluctuations in nature mainly affect the parameters of the model, which can be characterized as white noise ^[33], the following stochastic perturbation predator-prey model with group defense is established:

Let ${B_1}(t)$ and ${B_2}(t)$ be independent standard Brownian motions defined on a complete probability space $\left({\Omega, {\cal F}, {{\left\{ {{{\cal F}_t}} \right\}}_{t \ge 0}}, {\cal P}} \right)$, ${\sigma _i} > 0\left({i = 1, 2} \right)$ is the intensity of environmental noise.

By constructing suitable Lyapunov functions and applying $Itô$ formula, scholars have studied the dynamics of the stochastic predator-prey models with hunting cooperation ^[3], prey refuge and fear effect ^[4,34], Beddington-DeAngelis functional response ^[35], Holling type Ⅱ-Ⅲ functional response ^[36], Allee effect ^[37], habitat complexity and prey aggregation ^[38], etc. As mentioned above, the dynamic properties of the deterministic predator-prey model (1.1) with group defense behavior are very complicated, so far, there is no literature on the stochastic predator-prey model (1.2) with group defense behavior. Therefore, one purpose of this paper is to study the effects of stochastic noise on the dynamics of the deterministic model with group defense behavior. We give how does the globally asymptotically stable equilibrium point or limit cycle of a deterministic model change under the effects of white noise. The other purpose of this paper is to analyze the impact of group defense behavior and other key factors on the extinction or persistence of the population.

The rest of this paper is arranged as follows. In Section 2, Lyapunov function analysis method is used for proving the existence and uniqueness of positive solution of model (1.2). In Section 3, the extinction and persistence of the population in model (1.2) is discussed, and proves that there is a unique ergodic stationary distribution. Numerical simulations are given in Section 4 and Section 5 to support our conclusions.

2.   Existence and uniqueness of the positive solution

In this paper, we define the state space as $\mathbb{R}_+^2 = \{(x_1, x_2)\in \mathbb{R}^2| x_1 > 0, x_2 > 0\}$.

Theorem 2.1. For any initial value $(x(0), y(0))\in \mathbb{R}_+^2$, there is a unique solution $(x(t), y(t))$ of model (1.2) on $t\ge0$ and the solution will remain in $\mathbb{R}_+^2$ with probability one, that is to say, $(x(0), y(0))\in \mathbb{R}_+^2$ for all $t\ge0$ almost surely (a.s.).

Proof. Since the coefficients of model (1.2) satisfy the local Lipschitz condition, for any given initial value $(x(0), y(0))\in \mathbb{R}_+^2$, there exists a unique local positive solution $(x(t), y(t))$ on $t \in \left[{0, {\tau _e}} \right)$, where $\tau _e$ denotes the explosion time. In order to prove that solution $(x(t), y(t))$ is global, it is sufficient to prove ${\tau _e} = \infty$.

Take a sufficiently large non-negative number ${n_0}$, such that $x\left(0 \right) \in \left[{\frac{1}{{{n_0}}}, {n_0}} \right]$, $y\left(0 \right) \in \left[{\frac{1}{{{n_0}}}, {n_0}} \right]$. For each integer $n \ge {n_0}$, define the stopping time as

Obviously, ${\tau _n}$ is increasing as $n \to \infty$. Set ${\tau _\infty } = \mathop {\lim }\limits_{n \to \infty } {\tau _n}$, hence ${\tau _\infty } \le {\tau _e}$ a.s. If we can prove ${\tau _\infty } = \infty$ a.s., then ${\tau _e} = \infty$ a.s.

If the statement is false, then there are a pair of constants $T > 0$ and $\varepsilon \in \left({0, 1} \right)$ such that

Thus, there is an integer ${n_1} \ge {n_0}$ such that

Define a ${C^2}$-function $V: \mathbb{R}_+^2 \to \mathbb{R}_+$ by

Making use of Itô formula, we obtain

where

$\tilde{P}$ is a positive constant. So we have

Integrating both sides of (2.2) from 0 to ${\tau _n} \wedge T = \min \left\{ {{\tau _n}, T} \right\}$ and then taking the expectation, we obtain

therefore

Let ${\Omega _n} = \left\{ {{\tau _n} \le T} \right\}$ for $n > {n_1}$. From (2.1), we obtain $P\left({{\Omega _n}} \right) > \varepsilon$. For each $\omega \in {\Omega _n}$, $x\left({{\tau _n}, \omega } \right)$ or $y\left({{\tau _n}, \omega } \right)$ equals either $n$ or $\frac{1}{n}$. So $V\left({x\left({{\tau _n}, \omega } \right), y\left({{\tau _n}, \omega } \right)} \right)$ is no less than either

Thus

holds. According to (2.3), we have

where ${I_{{\Omega _n}}}$ denotes the indicator function of ${\Omega _n}$. Letting $n \to \infty$, then we obtain

which leads to a contradiction, so we must have ${\tau _\infty } = \infty$ a.s. This completes the proof.

3.   Extinction and stationary distribution

3.1.   Extinction

In this section, we first discuss the conditions for the extinction of population. The following definition is given:

Definition 3.1. If $\mathop {\lim }\limits_{t \to + \infty } x(t) = 0$ a.s. holds, the prey population will be extinct with probability one; If $\mathop {\lim }\limits_{t \to + \infty } y(t) = 0$ a.s. holds, the predator population will be extinct with probability one.

Theorem 3.1. For model (1.2), if ${\sigma _1} > \sqrt {2r}$ holds, then the prey population $x(t)$ will be extinct with probability one.

Proof. Applying Itô's formula to the first equation of model (1.2), we have

Integrating both sides from 0 to $t$ and dividing by $t$, we obtain

Making use of the strong law of large numbers for local martingales ^[39] leads to $\mathop {\lim }\limits_{t \to \infty } { {{{\sigma _1}{B_1}(t)} \over t}} = 0{\rm{ }}$ a.s. Therefore

That is to say $\mathop {\lim }\limits_{t \to + \infty } x(t) = 0$ a.s., according to definition 3.1, the prey population will be extinct with probability one. This completes the proof.

Remark 3.1. Theorem 3.1 gives a sufficient condition for the extinction of prey population. When ${\sigma _1} > \sqrt {2r}$, the prey population will be extinct, and it is not difficult to draw from model (1.2) that if the prey population is extinct, it will lead to the extinction of predator population. When ${\sigma _1} \ne 0$, from the above proof process, it can be seen that even small effects of environmental noise will reduce the density of prey population in mean. Once the environmental noise reaches a certain level, the prey population will be extinct.

Theorem 3.2. For model (1.2), if ${\sigma _2} > \sqrt {2\left({{ {{\mu \alpha } \over \beta }} - D} \right)}$ holds, then the predator population $y(t)$ will be extinct with probability one.

Proof. Applying Itô's formula to the second equation of model (1.2), we have

Integrating both sides from 0 to $t$ and dividing by $t$, we obtain

So

Then we have $\mathop {\lim }\limits_{t \to \infty } y(t) = 0$. This completes the proof.

Remark 3.2. Note that when ${\sigma _2} \ne 0$, that is, in the presence of environmental noise, once the stochastic noise ${\sigma _2}$ is larger than the critical threshold ($\sqrt {2\left({{ {{\mu \alpha } \over \beta }} - D} \right)}$), then the predator population will tend to be extinct. From the expression of this threshold, it is not difficult to see that the lower group defense level $\beta$ of the prey is conductive to the survival of the predator population.

3.2.   Stationary distribution

Whether the predator and prey populations will continue to survive is also an important issue in the study of predator-prey model. Therefore, in this section, we will prove that model (1.2) has a unique ergodic stationary distribution through the equivalent conditions of Hasminskii theorem, indicating that both predators and prey populations will survive ^[40].

Theorem 3.3. If $R = {\left({D + { {{\sigma _2^2} \over 2}}} \right)^{ - 1}}\beta \left({r - { {{\sigma _1^2} \over 2}} - { {{D{{\left({{ {r \over K}} + { {{\mu \alpha } \over D}}r} \right)}^2}} \over {4\mu \alpha }}}} \right) > 1$, then model (1.2) has a unique stationary distribution $\pi(\cdot)$ with ergodicity.

Proof. Model (1.2) can be written as the following form

then we get its diffusion matrix

and there exists a positive number $G = \mathop {\min }\limits_{(x, y) \in \bar U} \left\{ {\sigma _1^2{x^2}, \sigma _2^2{y^2}} \right\}$, such that $\sum\limits_{i, j = 1}^2 {{a_{ij}}} {\xi _i}{\xi _j} = \sigma _1^2{x^2}{\xi _1}^2 + \sigma _2^2{y^2}{\xi _2}^2 \ge G{\left| \xi \right|^2}$, $\left({x, y} \right) \in \bar U$, $\xi = \left(\xi_{1}, \xi_{2}\right) \in \mathbb{R}^{2}$.

Construct a ${C^2}$-function $V: \mathbb{R}_+^2 \to \mathbb{R}_+$ by

where

It is easy to prove that $\bar V\left({x, y} \right)$ has a unique minimum point $\left({{x_0}, {y_0}} \right)$. Define a ${C^2}$-function $V: \mathbb{R}_+^2 \to \mathbb{R}_+$ by

where ${V_1} = - P\left[{\alpha \ln x + \ln y - \frac{{\beta \alpha }}{D}\left({\mu x + y} \right)} \right]$, ${V_2} = \frac{1}{2}{\left({\mu x + y} \right)^2} - \bar V\left({{x_0}, {y_0}} \right)$. Applying Itô's formula to (3.1), we have

Therefore

Choosing sufficiently small ${\varepsilon _1}$ and ${\varepsilon _2}$ such that

Then we consider the bounded open set

Clearly, $S^c = \mathbb{R}_+^2\backslash S = {S_1} \cup {S_2} \cup {S_3} \cup {S_4}$, where

Case 1. When $\left({x, y} \right) \in {S_1}$, and $xy < {\varepsilon _1}y < {\varepsilon _1}\left({1 + {y^2}} \right)$, we have

By (3.2) and the definition of $P$, we can obtain

Case 2. When $\left({x, y} \right) \in {S_2}$, and $xy < {\varepsilon _2}x < {\varepsilon _2}\left({1 + {x^3}} \right)$, we have

By (3.3) and the definition of $P$ again, we can obtain $LV \le - {{{P^*}\left({R - 1} \right)} \over 4} \le - 1$, for $\left({x, y} \right) \in {S_2}$.

In addition, by Young inequality, we have $xy \le { {2 \over 5}}{x^{{ {5 \over 2}}}} + { {3 \over 5}}{y^{{ {5 \over 3}}}}$, for $\forall x, y > 0$. Hence,

Case 3. When $\left({x, y} \right) \in {S_3}$, we have by (3.2) and (3.5)

where

By (3.4) and (3.6), we have $LV \le -1$, for $\forall (x, y) \in S_3$.

Case 4. When $\left({x, y} \right) \in {S_4}$, similar to the case 3, by (3.4) we have,

In conclusion, $LV \le -1$, $\forall (x, y) \in S^c$. According to the equivalent conditions of Hasminskii theorem in ^[40], we know that the model (1.2) is ergodic and has a unique stationary distribution.

4.   The effects of environmental noise on model (1.2)

In this section, we will discuss the effects of environmental noise on extinction and persistence of populations, and further investigate how the environmental noise affects the dynamics of deterministic system through numerical simulations.

Take parameters as follows: $r = 0.5$, $\mu = 0.4$, $\alpha = 0.4$, $K = 4$, $D = 0.15$, $\beta = 0.38$. According to Theorem 2.1 in ^[18], since $e\beta D < \mu \alpha$ and ${x_1} < K < 1/\beta$, deterministic model (1.1) has a globally asymptotically stable equilibrium $\left({{x_1}, {y_1}} \right)$, which means that both populations are persistent (see the dotted line in Figure 1(a) and the red curve in Figure 1(b)). Now consider the effects of small environmental noise on stochastic model (1.2). Taking ${\sigma _1} = {\sigma _2} = 0.01$, then $R = {\left({D + \frac{{\sigma _2^2}}{2}} \right)^{ - 1}}\beta \left({r - \frac{{\sigma _1^2}}{2} - \frac{{D{{\left({\frac{r}{K} + \frac{{\mu \alpha }}{D}r} \right)}^2}}}{{4\mu \alpha }}} \right) \approx 1.0089 > 1$, it is clear from Theorem 3.3 that there is a unique ergodic stationary distribution for model (1.2), which implies that both populations are still persistent(see the solid line in Figure 1(a) and the blue curve in Figure 1(b)), and the solution of stochastic model (1.2) fluctuates near $\left({{x_1}, {y_1}} \right)$. In addition, we can also get the density function distribution of $x(t)$, $y(t)$, see Figure 1(c)(d). Moreover, we simulate the effects of environmental noises ${\sigma _1}$ and ${\sigma _2}$ on the threshold $R$. It is not difficult to find that $R > 1$ in Theorem 3.3 is a strong sufficient condition, which can only be satisfied when ${\sigma _1}$ and ${\sigma _2}$ are small, see Figure 2.

However, when the environmental noise $\sigma_1$ is small and $\sigma_2$ is large (${\sigma _1} = 0.01$, ${\sigma _2} = 1$), $\beta = 0.25$, by calculations we have ${\sigma _2} > \sqrt {2\left({{ {{\mu \alpha } \over \beta }} - D} \right)} \approx 0.9899$. From Theorem 3.2, we conclude that the predator population $y(t)$ goes extinct (see Figure 3(a)). We also find that the solution of stochastic model (1.2) fluctuates around the hyperbolic saddle $\left({K, 0} \right)$ of deterministic model (1.1), which means that environmental noise causes the solution of a deterministic model (1.1) to transform from tending to a globally asymptotically stable equilibrium point $\left({{x_1}, {y_1}} \right)$ to fluctuating near another equilibrium point $\left({K, 0} \right)$. In this case, predator population is extinct, and the prey fluctuates around the environmental capacity $K$.

Let environmental noise $\sigma_1$ continue to increase and take ${\sigma _1} = 1.3$, ${\sigma _2} = 1$. By calculations, we get ${\sigma _1} > \sqrt {2r} = 1$. Theorem 3.1 and Remark 3.1 show that both populations are extinct, see Figure 3(b). That is to say the environmental noise $\sigma_1$ is unfavorable for the persistence of prey and predator population.

In fact, the conditions in Theorem 3.1, Theorem 3.2 and 3.3 are only sufficient conditions. When the condition is not satisfied, the extinction or persistence of the population cannot be determined. Take Theorem 3.1 as an example, when ${\sigma _1} = 1$, ${\sigma _2} = 0.01$, although the condition ${\sigma _1} > \sqrt {2r}$ in Theorem 3.1 cannot hold, the prey population is still extinct, see Figure 4.

Take $r = 0.5$, $\mu = 0.4$, $\alpha = 0.4$, $K = 6$, $D = 0.15$, $\beta = 0.25$. Since $e\beta D < \mu \alpha$ and ${x_1}+1/\beta < K < x_2$, by the results in ^[18] and numerical simulation, we know deterministic model (1.1) has a stable limit cycle in the interior of the first quadrant. Now we investigate how $\sigma_1$ or $\sigma_2$ affects the dynamics of deterministic model (1.1) in this case, see Figure 5(a)–(d). When ${\sigma _1}$ and ${\sigma _2}$ are small (${\sigma _1} = {\sigma _2} = 0.01$ and ${\sigma _1} = {\sigma _2} = 0.05$), the solutions of model (1.2) still fluctuate around the stable limit cycle of model (1.1), see Figure 5(a), (b); As ${\sigma _1}$ and ${\sigma _2}$ increase gradually, environmental noises make the dynamics of deterministic model(1.1) very complex, see Figure 5(c) ((${\sigma _1} = {\sigma _2} = 0.1$)); If ${\sigma _1}$ and ${\sigma _2}$ continue to increase, according to Theorem 3.1 and Theorem 3.2, the population may become extinct. For example, take ${\sigma _1} = 0.05$, ${\sigma _2} = 0.65$, the solutions of model (1.2) tends to the boundary equilibrium point $(K, 0)$ of model (1.1) and predator population will be extinct, see Figure 5(d).

In summary, the environmental noise makes the dynamics of the deterministic system more complex, and the greater the intensity of the environmental noise, the less conducive to the survival of the population.

5.   Sensitivity analysis on threshold $R$

In order to identify the key factors affecting populations stability, a sensitivity analysis is performed on threshold $R$.

Sensitivity analysis is a method of quantitative uncertainty analysis for any complex model. Here we analyze the sensitivity of the parameter to $R$ by calculating the Partial Rank Correlation Coefficient (PRCC) value of the parameter, and determine the key parameters that affect the threshold. When the $PRCC$ value of the parameter is positive, the parameter is positively correlated with $R$. Conversely, the parameter is negatively correlated with $R$. We make the following provisions: if $\left| {PRCC} \right| > 0.4$, then there is a strong correlation between the input parameter and the output variable, i.e., the parameter has a strong influence on $R$; if $0.2 < \left| {PRCC} \right| < 0.4$, this parameter has a moderate influence on $R$. If $\left| {PRCC} \right| < 0.2$, then this parameter has a weak influence on $R$. For more details, see ^[41].

Take $\beta = 0.25$, ${\sigma _1} = {\sigma _2} = 0.01$, the sample size is 1500 and all parameters vary simultaneously. From Figure 6, the $PRCC$ value of $r$, $\beta$, $K$, $D$ are positive, it shows that parameter $r$, $\beta$, $K$, $D$ are positively correlated with the threshold $R$. The increase of these parameters increases $R$, which is conductive to the persistence of populations. And the $PRCC$ value of $\mu$, $\alpha$, $\sigma _1$, $\sigma _2$ are negative, so $\mu$, $\alpha$, $\sigma _1$, $\sigma _2$ are negatively correlated with the threshold $R$. The increase of these parameters decreases $R$, which is harmful to the populations. Among them, $\left| {PRCC} \right|_{r, \sigma_{1}} > 0.4$, so $r$, $\sigma_{1}$ have a strong correlation with $R$; $0.2 < \left| {PRCC} \right|_{\mu, \alpha, D} < 0.4$, so $\mu$, $\alpha$ and $D$ have a moderate influence on $R$; $\left| {PRCC} \right|_{\beta, K, \sigma_{2}} < 0.2$, so $\beta$, $K$ and $\sigma_{2}$ have a weak influence on $R$.

6.   Conclusions

In this paper, a stochastic predator-prey model with group defense behavior is established. The dynamics of stochastic model (1.2) are studied through theoretical analysis and numerical simulations.

According to Theorem 3.1 and 3.2, when the environmental noise has a great impact on prey population, both populations will be extinct; when the environmental noise has a great impact on predator population, the predator population will be extinct, which will destroy the stability between two populations. By numerical simulations, we show that the environmental noise makes the dynamics of the deterministic system more complex.

Under the influence of environmental noise, the appropriate group defense level of prey can help the survival of two populations, and maintain the stability of the relationship between two populations(see Figure 7(a), $\beta = 0.25$). From Theorem 2.2 in ^[18], if $e\beta D < \mu \alpha$ and $K > \frac{1}{\beta} + {x_2}$, then model (1.1) has a hyperbolic stable node $\left({K, 0} \right)$. So the larger group defense level is not conducive to the persistence of predator population, resulting in inundation of prey population (see Figure 7(b), $\beta = 0.39$). The threshold $R$, which plays a crucial role in maintaining stability of the relationship between two populations, is strongly influenced by the intrinsic growth rate of prey population $r$ and the intensity of environmental noise $\sigma_1$.

In this paper, we only consider the influence of group defense behavior and environmental noise on the dynamics of stochastic predator-prey model. However, there are other more complex influences in nature. For example, in the study of pest management with group defense behavior, the effects of releasing natural enemies, impulsive spraying pesticide and environmental noise on agricultural production should be considered simultaneously. We will enrich it in future work.

Acknowledgments

This work is supported by National Natural Science Foundation of China (12171004).

Conflict of interest

No potential conflict of interest was reported by the authors.