2019 Guangzhou Workshop

Survival analysis of a stochastic predator–prey model with prey refuge and fear effect

Abstract
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

In this paper, we propose and investigate a stochastic Holling type-II predator–prey model with prey refuge and fear effect. We first prove the existence and uniqueness of the global positive solution. Then we perform the survival analysis of the model, including the existence of a unique ergodic stationary distribution and the extinction of the model. Numerical simulations are carried out to validate our analytical results. Our findings indicate that the white noise is adverse to the growth of predator and prey populations, and the increase of fear effect will lead to the decrease of predator density, but with no obvious effect on prey density.

Keywords:

Stochastic predator–prey model fear effect prey refuge stationary distribution extinction

MSC SUBJECT CLASSIFICATIONS:

60H10 92D25

1. Introduction

Population models, as an important part of ecology, have been extensively studied and explored for their rich dynamic behaviour with the aim to provide a theoretical guidance for the protection, development and utilization of biological resources [2,5]. Among the most important population models, predator–prey models play an important role in understanding the interaction of different species in volatile natural environments and have been widely investigated [1,4,9,12,30,33].

As is well known, it is common for predators to hunt preys in nature. However, not all preys are captured by the predators because the preys usually have refuges where they can evade the predation [10,18,27]. This indicates that a prey refuge can protect and prevent the extinction of preys, and thus, it plays an important role in the interaction between predator and prey populations. Recently, the research of the predator–prey models with prey refuges has attracted many mathematicians' attention and some interesting results have been obtained [19,23,25,32].

Notice that in most previous predator–prey models, only the direct killing of preys in the presence of predators is considered due to its obviousness in nature. But in reality, the preys may change their behaviour when the predator they faced is powerful [13,26,28]. Some studies have shown that all preys respond to the risk of predation and show all sorts of anti-predator responses, including habitat changes, foraging, vigilance and different physiological changes. Such anti-predator behaviour can reduce the long-term cost of reproduction and increase the probability of adult survival [4]. In addition, frightened preys often forage less, so their birth rate decline and they use some survival mechanisms like starvation [3].

Considering the effect of prey refuge and fear effect, in a recent paper [37], Zhang et al. proposed a predator–prey model which takes the following form: (1) $\begin{aligned} \frac{d x}{d t} & = \frac{α x}{1 + K y} - b x^{2} - \frac{β (1 - m) x y}{1 + a (1 - m) x}, \\ \frac{d y}{d t} & = - γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x}, \end{aligned}$ (1) where $x (t)$ and $y (t)$ denote, respectively, the densities of prey and predator population at time t, α is the intrinsic growth rate of prey, b is the strength of intraspecific competition of prey, γ is the death rate of the predator, c is the conversion coefficient, K refers to the level of fear and $β / a$ is the maximum number of prey eaten per predator per unit of time. All these parameters are assumed to be positive. The term $β x / (1 + a x)$ represents the Holling type-II functional response, and $(1 - m) x$ is the quantity of preys available to the predators, where $m \in [0, 1)$ is the protection rate of the prey refuge for prey. For this model, the authors performed a detailed stability and bifurcation analysis. We refer the readers to Ref. [37] for more details.

In fact, a real ecosystem is inevitably affected by environmental noise. Therefore, the shifting environmental effects can not be neglected, which indicates that deterministic prey–predator models have limitations in predicting dynamics accurately, while stochastic models can make it [38]. Many researchers introduced random disturbances into deterministic systems to reveal the influence of environmental noise [14,16,17,22,29,31,34–36,39–41]. Ripa et al. [24] examined the impact of environmental noise on populations and presented a general theory of environmental noise in ecological food webs. Mao et al. [21] proved that the explosion of population dynamics can be suppressed by environmental Brownian noise and obtained that even a sufficiently small noise perturbation can suppress explosions in population dynamics. Hence, it is of great significance to further incorporate environmental randomness into the deterministic model (1). Applying the technique used in Ref. [6] to include stochastic effects, we can obtain the stochastic version of model (1) as follows (see the appendix for specific proof): (2) $\begin{aligned} d x & = [\frac{α x}{1 + K y} - b x^{2} - \frac{β (1 - m) x y}{1 + a (1 - m) x}] d t + σ_{1} x d B_{1} (t), \\ d y & = [- γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x}] d t + σ_{2} y d B_{2} (t), \end{aligned}$ (2) where $B_{1} (t)$ and $B_{2} (t)$ are two independent standard Brownian motions defined in a complete probability space $(Ω, F, {F_{t}}_{t \geq 0}, P)$ with a filtration ${F_{t}}_{t \geq 0}$ satisfying the usual conditions, and $σ_{i}^{2}, i = 1, 2$ are the intensities of the white noises. In this paper, we will devote ourselves to the investigation on the dynamics of the stochastic model (2).

The organization of this paper is as follows. In the next section, we present some necessary notations and preliminary results of model (2). In Section 3, we prove the existence of a unique ergodic stationary distribution of the model, and in Section 4, we perform the extinction analysis of the model. Some numerical simulations are carried out in Section 5 to illustrate our theoretical results. We close the paper with conclusion in Section 6.

2. Preliminaries

In this section, we present some notations and basic results about model (2), which are the basis for further investigation on the dynamics of model (2).

Consider the following n-dimensional stochastic differential equation: (3) $d x (t) = f (x (t), t) d t + g (x (t), t) d B (t) f o r t \geq t_{0} .$ (3) Denote by $C^{2, 1} (R^{n} \times [t_{0}, \infty); R_{+})$ the family of all nonnegative functions $V (x, t)$ defined on $R^{n} \times [t_{0}, \infty)$ such that they are continuously twice differentiable in x and once in t. We introduce the differential operator defined in Ref. [20] $L = \frac{\partial}{\partial t} + \sum_{i = 1}^{n} f_{i} (x, t) \frac{\partial}{\partial x_{i}} + \frac{1}{2} \sum_{i, j = 1}^{n} [g^{T} (x, t) g (x, t)]_{i j} \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} .$ If L acts on a function $V \in C^{2, 1} (R^{n} \times [t_{0}, \infty); R_{+}),$ then $L V (x, t) = V_{t} (x, t) + V_{x} (x, t) f (x, t) + \frac{1}{2} t r a c e [g^{T} (x, t) V_{x x} (x, t) g (x, t)],$ where $V_{t} = \partial V / \partial t, V_{x} = (\partial V / \partial x_{1}, \dots, \partial V / \partial x_{n})$ and $V_{x x} = (\partial^{2} V / \partial x_{i} \partial x_{j})_{n \times n}$ . From Itô's formula, if $x (t) \in R^{n}$ , then $d V (x (t), t) = L V (x (t), t) d t + V_{x} (x (t), t) g (x (t), t) d B (t) .$ The following theorem is about the existence and uniqueness of the global positive solution of model (2).

Theorem 2.1

For any given initial value $(x (0), y (0)) \in R_{+}^{2}$ , there exists a unique solution $(x (t), y (t))$ of model (2) on $t \geq 0$ and the solution will remain in $R_{+}^{2}$ with probability one, that is to say, $(x (t), y (t)) \in R_{+}^{2}$ for all $t \geq 0$ almost surely (a.s.).

Proof.

Since the coefficients of model (2) satisfy the local Lipschitz condition, then there exists a unique local solution $(x (t), y (t))$ of model (2) on $t \in [0, τ_{e})$ , a.s., where $τ_{e}$ denotes the explosion time. To show that this solution is global, we only need to prove $τ_{e} = \infty$ , a.s. Let $k_{0} > 0$ be sufficiently large such that both $x (0)$ and $y (0)$ lie within the interval $[1 / k_{0}, k_{0}]$ . For each integer $k \geq k_{0}$ , define the stopping time $τ_{k} = inf {t \in [0, τ_{e}) : min {x (t), y (t)} \leq \frac{1}{k}, o r max {x (t), y (t)} \geq k},$ where $τ_{k}$ is increasing as $k \to \infty$ . Set $τ_{\infty} = lim_{k \to \infty} τ_{k}$ , which implies $τ_{\infty} \leq τ_{e}$ , a.s. Hence to complete the proof, we only need to show that $τ_{\infty} = \infty$ , a.s. We prove this by contradiction. If $τ_{\infty} < \infty$ , then there exist a pair of constants $ϵ \in (0, 1)$ and T>0 such that $P (τ_{\infty} < T) > ϵ .$ Therefore, there exists $k_{1} \geq k_{0}$ such that for all $k \geq k_{1}$ , $P (Ω_{k}) \geq ϵ,$ where $Ω_{k} = {τ_{k} \leq T}$ . Then define a $C^{2}$ -function V: $R_{+}^{2} \to R_{+}$ by $V = x - \frac{γ}{2 c β (1 - m)} - \frac{γ}{2 c β (1 - m)} \ln \frac{2 c β (1 - m) x}{γ} + \frac{1}{c} (y - 1 - \ln y) .$ From Itô's formula, (4) $d V (x, y) = L V (x, y) d t + [σ_{1} x - \frac{γ σ_{1}}{2 c β (1 - m)}] d B_{1} (t) + (\frac{σ_{2} y}{c} - \frac{σ_{2}}{c}) d B_{2} (t),$ (4) where $\begin{aligned} L V (x, y) & = \frac{α x}{1 + K y} - b x^{2} - \frac{β (1 - m) x y}{1 + a (1 - m) x} - \frac{α γ}{2 c β (1 - m) (1 + K y)} + \frac{b γ x}{2 c β (1 - m)} \\ + \frac{γ}{2 c} \frac{y}{1 + a (1 - m) x} + \frac{γ σ_{1}^{2}}{4 c β (1 - m)} - \frac{γ}{c} y + \frac{β (1 - m) x y}{1 + a (1 - m) x} + \frac{γ}{c} \\ - \frac{β (1 - m) x}{1 + a (1 - m) x} + \frac{σ_{2}^{2}}{2 c} \\ \leq α x - b x^{2} + \frac{b γ x}{2 c β (1 - m)} + \frac{γ}{2 c} y + \frac{γ σ_{1}^{2}}{4 c β (1 - m)} - \frac{γ}{c} y + \frac{γ}{c} + \frac{σ_{2}^{2}}{2 c} \\ \leq - b x^{2} + [α + \frac{b γ}{2 c β (1 - m)}] x - \frac{γ}{2 c} y + \frac{γ σ_{1}^{2}}{4 c β (1 - m)} + \frac{γ}{c} + \frac{σ_{2}^{2}}{2 c} \\ \leq \frac{γ σ_{1}^{2}}{4 c β (1 - m)} + \frac{γ}{c} + \frac{σ_{2}^{2}}{2 c} + \tilde{J} := J . \end{aligned}$ Here, $\tilde{J} = sup_{(x, y) \in (0, \infty)} {- b x^{2} + [α + \frac{b γ}{2 c β (1 - m)}] x - \frac{γ}{2 c} y} .$ Integrating and taking the expectation of both sides of (4) from 0 to $τ_{k} \land T$ , we obtain (5) $E V (x (τ_{k} \land T), y (τ_{k} \land T)) \leq V (x (0), y (0)) + J E (τ_{k} \land T) \leq V (x (0), y (0)) + J T .$ (5) Noting that for every $ω \in Ω_{k}$ , $x (τ_{k}, ω)$ or $y (τ_{k}, ω)$ equals either k or $1 / k$ . Hence $\begin{aligned} V (x (τ_{k}), y (τ_{k})) & \geq [k - \frac{γ}{2 c β (1 - m)} - \frac{γ}{2 c β (1 - m)} \ln \frac{2 c β (1 - m) k}{γ}] \\ \land [\frac{1}{k} - \frac{γ}{2 c β (1 - m)} - \frac{γ}{2 c β (1 - m)} \ln \frac{2 c β (1 - m)}{γ k}] \\ \land [\frac{1}{c} (k - 1 - \ln k)] \land [\frac{1}{c} (\frac{1}{k} - 1 + \ln k)] . \end{aligned}$ From (5), we have $\begin{aligned} V (x (0), y (0)) + J T & \geq E (I_{Ω_{k}} V (x (τ_{k}), y (τ_{k}))) \\ \geq ϵ {[k - \frac{γ}{2 c β (1 - m)} - \frac{γ}{2 c β (1 - m)} \ln \frac{2 c β (1 - m) k}{γ}] \\ \land [\frac{1}{k} - \frac{γ}{2 c β (1 - m)} - \frac{γ}{2 c β (1 - m)} \ln \frac{2 c β (1 - m)}{γ k}] \\ \land [\frac{1}{c} (k - 1 - \ln k)] \land [\frac{1}{c} (\frac{1}{k} - 1 + \ln k)]} . \end{aligned}$ This leads to a contradiction as we let $k \to \infty$ , $\infty > V (x (0), y (0)) + J T = \infty .$ So we have that $τ_{\infty} = \infty$ , a.s.

Next, we present a lemma which gives a condition for the existence of a unique ergodic stationary distribution of system (3). Let $X (t)$ be a homogeneous Markov process in $R^{n}$ described by the following stochastic differential equation: $d X (t) = b (X) d t + \sum_{r = 1}^{k} g_{r} (X) d B_{r} (t),$ and the diffusion matrix is $A (x) = (a_{i j} (x)), a_{i j} (x) = \sum_{r = 1}^{k} g_{r}^{i} (x) g_{r}^{j} (x) .$

Lemma 2.2

[7]

The Markov process $X (t)$ has a unique ergodic stationary distribution $μ (\cdot)$ , if there exists a bounded domain $D \subset R^{n}$ with regular boundary Γ and

there is a positive number M such that $\sum_{i, j = 1}^{n} a_{i j} (x) ξ_{i} ξ_{j} \geq M {|ξ|}^{2}, x \in D, ξ \in R^{n} .$
There exists a nonnegative $C^{2} -$ function V such that LV is negative for any $R^{n} ∖ D$ . Then $P \{lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} f (X (t)) d t = \int_{R^{n}} f (x) μ (d x)\} = 1$ for all $x \in R^{n}$ , where $f (\cdot)$ is a function integrable with respect to the measure μ.

3. Stationary distribution

In this section, we establish sufficient conditions for the existence of a unique ergodic stationary distribution. First, we give the following assumptions:

(H1)	$γ > \frac{1}{2} max {σ_{1}^{2}, σ_{2}^{2}}$ ;
(H2)	$λ = \frac{2 α c β (1 - m)}{b + 2 a α (1 - m)} - \frac{2 α c β (1 - m)}{3 b [1 + a (1 - m) (2 α / b)]^{2}} - \frac{4 c_{2} α^{2}}{3 b} - \frac{c_{2} α σ_{1}^{2}}{2 b} - γ - \frac{σ_{2}^{2}}{2} > 0$ , where $c_{2}$ is a positive constant satisfying that (6) $c_{2} > \frac{3 c β a (1 - m)^{2}}{b {[1 + a (1 - m) \frac{2 α}{b}]}^{3}} + \frac{c β (1 - m)}{α {[1 + a (1 - m) \frac{2 α}{b}]}^{2}} .$ (6)

Theorem 3.1

Let $(x (t), y (t))$ be the solution of model (2) with any given initial value $(x (0), y (0)) \in R_{+}^{2}$ . If assumptions ( $H_{1}$ ) and $(H_{2})$ are satisfied, then model (2) admits a unique stationary distribution $μ (\cdot)$ and it has the ergodic property.

Proof.

By the condition (B.2) in Lemma 2.2, we need to construct a $C^{2}$ -continuous function $V : R_{+}^{2} \to R_{+}$ and a bounded closed $U_{ϵ}$ such that $L V \leq - 1$ for $(x, y) \in R_{+}^{2} ∖ U_{ϵ}$ . To this end, let $σ_{max} = max {σ_{1}, σ_{2}}$ and choose a constant $θ > 0$ such that $ρ = γ - (1 / 2) σ_{max}^{2} (θ + 1) > 0$ , and define a $C^{2}$ -function $\tilde{V} : R_{+}^{2} \to R$ as follows: $\tilde{V} = M [- c_{1} x + c_{2} (x - \frac{α}{b} \ln x) + \frac{c_{2} β (1 - m) α}{b γ} y - \ln y] + \frac{1}{θ + 2} {(x + \frac{y}{c})}^{θ + 2},$ where $c_{1} = (c β (1 - m)) / (3 α [1 + a (1 - m) (2 α / b)]^{2}) + (2 / 3) c_{2}$ and $c_{2}$ is defined in (6) (obviously, condition (6) implies that $c_{2} > c_{1}$ ); $M = (2 / λ) max {2, R_{1}}$ and $R_{1}$ is to be determined later. We can easily check that $\underset{ϵ \to 0, (x, y) \in R_{+}^{2} ∖ U_{ϵ}}{lim inf} \tilde{V} (x, y) = \infty,$ where $U_{ϵ} := (ϵ, 1 / ϵ) \times (ϵ, 1 / ϵ)$ , and $ϵ > 0$ is a sufficiently small number. Notice that $\tilde{V} (x, y)$ has a global minimum point $(x_{0}, y_{0})$ in the interior of $R_{+}^{2}$ . Then, we can define a nonnegative $C^{2}$ -function $V : R_{+}^{2} \to R_{+}$ by (7) $V (x, y) = \tilde{V} (x, y) - \tilde{V} (x_{0}, y_{0}) = M V_{1} (x, y) + V_{2} (x, y),$ (7) where $\begin{aligned} V_{1} (x, y) & = - c_{1} x + c_{2} (x - \frac{α}{b} \ln x) + \frac{c_{2} β (1 - m) α}{b γ} y - \ln y, \\ V_{2} (x, y) & = \frac{1}{θ + 2} {(x + \frac{y}{c})}^{θ + 2} - \tilde{V} (x_{0}, y_{0}) . \end{aligned}$ Applying Itô's formula to $V_{1} (x, y)$ , we have $\begin{aligned} L V_{1} & = - \frac{c_{1} α x}{1 + K y} + c_{1} b x^{2} + \frac{c_{1} β (1 - m) x y}{1 + a (1 - m) x} + \frac{c_{2} α x}{1 + K y} - c_{2} b x^{2} - \frac{c_{2} β (1 - m) x y}{1 + a (1 - m) x} \\ - \frac{c_{2} α^{2}}{b (1 + K y)} + c_{2} α x + \frac{c_{2} α β (1 - m) y}{b [1 + a (1 - m) x]} + \frac{c_{2} α σ_{1}^{2}}{2 b} - \frac{c_{2} β (1 - m) α}{b} y \\ + \frac{c_{2} c β^{2} (1 - m)^{2} α x y}{b γ [1 + a (1 - m) x]} + γ - \frac{c β (1 - m) x}{1 + a (1 - m) x} + \frac{σ_{2}^{2}}{2} \\ \leq \frac{α (c_{2} - c_{1}) x}{1 + K y} + c_{1} b x^{2} + c_{1} β (1 - m) x y - c_{2} b x^{2} + c_{2} α x + \frac{c_{2} α β (1 - m)}{b} y + \frac{c_{2} α σ_{1}^{2}}{2 b} \\ - \frac{c_{2} α β (1 - m)}{b} y + \frac{c_{2} c β^{2} (1 - m)^{2} α x y}{b γ} + γ - \frac{c β (1 - m) x}{1 + a (1 - m) x} + \frac{σ_{2}^{2}}{2} \\ \leq c_{2} α x - c_{1} α x + c_{1} b x^{2} - c_{2} b x^{2} + c_{2} α x - \frac{c β (1 - m) x}{1 + a (1 - m) x} \\ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y + \frac{c_{2} α σ_{1}^{2}}{2 b} + γ + \frac{σ_{2}^{2}}{2} \\ = c_{2} α x (2 - \frac{b}{α} x) + c_{1} α [x (\frac{b}{α} x - 1) - \frac{2 α}{b}] - \frac{c β (1 - m) x}{1 + a (1 - m) x} \\ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y + \frac{2 c_{1} α^{2}}{b} + \frac{c_{2} α σ_{1}^{2}}{2 b} + γ + \frac{σ_{2}^{2}}{2} \\ = - \frac{c β (1 - m) (2 α / b)}{1 + a (1 - m) (2 α / b)} + \frac{2 c_{1} α^{2}}{b} + \frac{c_{2} α σ_{1}^{2}}{2 b} + γ + \frac{σ_{2}^{2}}{2} + {- \frac{c β (1 - m) x}{1 + a (1 - m) x} \\ + \frac{c β (1 - m) (2 α / b)}{1 + a (1 - m) (2 α / b)} + c_{2} α x (2 - \frac{b}{α} x) + c_{1} α [x (\frac{b}{α} x - 1) - \frac{2 α}{b}]} \\ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y \\ = - \frac{c β (1 - m) (2 α / b)}{1 + a (1 - m) (2 α / b)} + \frac{2 c_{1} α^{2}}{b} + \frac{c_{2} α σ_{1}^{2}}{2 b} + γ + \frac{σ_{2}^{2}}{2} \\ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y + F (x), \end{aligned}$ where $\begin{aligned} F (x) & = - \frac{c β (1 - m) x}{1 + a (1 - m) x} + \frac{c β (1 - m) (2 α / b)}{1 + a (1 - m) (2 α / b)} + c_{2} α x (2 - \frac{b}{α} x) \\ + c_{1} α [x (\frac{b}{α} x - 1) - \frac{2 α}{b}] . \end{aligned}$ We can compute that $F^{'} (x) = - \frac{c β (1 - m)}{[1 + a (1 - m) x]^{2}} + 2 c_{2} α - 2 c_{2} b x + 2 c_{1} b x - c_{1} α$ and $F^{''} (x) = \frac{2 c β a (1 - m)^{2}}{[1 + a (1 - m) x]^{3}} - 2 c_{2} b + 2 c_{1} b .$ Therefore, $F^{'} (x) |_{x = 2 α / b} = - \frac{c β (1 - m)}{{[1 + a (1 - m) (2 α / b)]}^{2}} - 2 c_{2} α + 3 c_{1} α = 0,$ due to the equality $c_{1} = (c β (1 - m)) / (3 α [1 + a (1 - m) (2 α / b)]^{2}) + (2 / 3) c_{2}$ . Moreover, we can compute that $F^{''} (x) |_{x = 2 α / b} = \frac{2 c β a (1 - m)^{2}}{{[1 + a (1 - m) (2 α / b)]}^{3}} - 2 c_{2} b + 2 c_{1} b < 0,$ where we have used inequality (6). Hence, we have that for all $x \in R_{+}$ , $F (x) \leq F (\frac{2 α}{b}) = 0.$ Consequently, (8) $\begin{aligned} L V_{1} & \leq - \frac{c β (1 - m) (2 α / b)}{1 + a (1 - m) (2 α / b)} + \frac{2 c_{1} α^{2}}{b} + \frac{c_{2} α σ_{1}^{2}}{2 b} + γ + \frac{σ_{2}^{2}}{2} \\ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y \\ = - \frac{c β (1 - m) (2 α / b)}{1 + a (1 - m) (2 α / b)} + \frac{2 α c β (1 - m)}{3 b {[1 + a (1 - m) (2 α / b)]}^{2}} \\ + \frac{4 c_{2} α^{2}}{3 b} + \frac{c_{2} α σ_{1}^{2}}{2 b} + γ + \frac{σ_{2}^{2}}{2} \\ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y \\ := - λ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y, \end{aligned}$ (8) where λ is a positive number defined in assumption ( $H_{2}$ ).

Similarly, we can compute that (9) $\begin{aligned} L V_{2} & = {(x + \frac{y}{c})}^{θ + 1} [\frac{α x}{1 + K y} - b x^{2} - \frac{β (1 - m) x y}{1 + a (1 - m) x} - \frac{γ}{c} y + \frac{β (1 - m) x y}{1 + a (1 - m) x}] \\ + \frac{1}{2} (θ + 1) {(x + \frac{y}{c})}^{θ} [σ_{1}^{2} x^{2} + σ_{2}^{2} {(\frac{y}{c})}^{2}] \\ \leq {(x + \frac{y}{c})}^{θ + 1} [(α + γ) x - b x^{2} - γ x - \frac{γ}{c} y] \\ + \frac{1}{2} (θ + 1) {(x + \frac{y}{c})}^{θ} [σ_{1}^{2} x^{2} + σ_{2}^{2} {(\frac{y}{c})}^{2}] \\ \leq {(x + \frac{y}{c})}^{θ + 1} [G - γ (x + \frac{y}{c})] + \frac{1}{2} σ_{max}^{2} (θ + 1) {(x + \frac{y}{c})}^{θ + 2} \\ = G {(x + \frac{y}{c})}^{θ + 1} - γ {(x + \frac{y}{c})}^{θ + 2} + \frac{1}{2} σ_{max}^{2} (θ + 1) {(x + \frac{y}{c})}^{θ + 2} \\ = G {(x + \frac{y}{c})}^{θ + 1} - [γ - \frac{1}{2} σ_{max}^{2} (θ + 1)] {(x + \frac{y}{c})}^{θ + 2} \\ \leq F - \frac{ρ}{2} {(x + \frac{y}{c})}^{θ + 2} \\ \leq F - \frac{ρ}{2} [x^{θ + 2} + {(\frac{y}{c})}^{θ + 2}], \end{aligned}$ (9) where $\begin{aligned} G & = sup_{(x, y) \in R_{+}^{2}} {(α + γ) x - b x^{2}} < \infty, \\ F & = sup_{(x, y) \in R_{+}^{2}} {G {(x + \frac{y}{c})}^{θ + 1} - \frac{ρ}{2} {(x + \frac{y}{c})}^{θ + 2}} < \infty . \end{aligned}$ From (7)–(9), we have that (10) $\begin{aligned} L V & = L (V_{1} + V_{2}) \leq M {- λ + [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y} \\ - \frac{ρ}{2} [x^{θ + 2} + {(\frac{y}{c})}^{θ + 2}] + F . \end{aligned}$ (10) Next, we prove $L V \leq - 1$ on $R_{+}^{2} ∖ U_{ϵ}$ . Take $ϵ > 0$ sufficiently small such that the following inequalities hold: (11) $\begin{aligned} 0 < ϵ < \frac{λ b γ}{4 [c_{1} β (1 - m) b γ + c_{2} c β^{2} (1 - m)^{2} α]}, \end{aligned}$ (11) (12) $\begin{aligned} 0 < ϵ < \frac{ρ b γ}{4 M c^{θ + 2} [c_{1} β (1 - m) b γ + c_{2} c β^{2} (1 - m)^{2} α]}, \end{aligned}$ (12) (13) $\begin{aligned} 0 < ϵ < \frac{ρ b γ}{4 M [c_{1} β (1 - m) b γ + c_{2} c β^{2} (1 - m)^{2} α]}, \end{aligned}$ (13) (14) $\begin{aligned} - M λ - \frac{ρ}{4 ϵ^{θ + 2}} + R_{2} \leq - 1, \end{aligned}$ (14) (15) $\begin{aligned} - M λ - \frac{ρ}{4 c^{θ + 2} ϵ^{θ + 2}} + R_{3} \leq - 1, \end{aligned}$ (15) where $R_{2}$ and $R_{3}$ are constants given by (17) and (18). Denote $\begin{aligned} U_{ϵ}^{1} = \{(x, y) \in R_{+}^{2} : 0 < x < ϵ\}, U_{ϵ}^{2} = \{(x, y) \in R_{+}^{2} : 0 < y < ϵ\}, \\ U_{ϵ}^{3} = \{(x, y) \in R_{+}^{2} : x > \frac{1}{ϵ}\}, U_{ϵ}^{4} = \{(x, y) \in R_{+}^{2} : y > \frac{1}{ϵ}\} . \end{aligned}$ Obviously, $R_{+}^{2} ∖ U_{ϵ} = U_{ϵ}^{1} ⋃ U_{ϵ}^{2} ⋃ U_{ϵ}^{3} ⋃ U_{ϵ}^{4} .$ Thus we only need to validate $L V \leq - 1$ in each domain $U_{ϵ}^{i}$ , where i = 1, 2, 3, 4.

Case 1. On the domain $U_{ϵ}^{1}$ , we have $x y \leq ϵ (1 + y^{θ + 2})$ . It then follows from (10) that $\begin{aligned} L V & \leq - \frac{M λ}{4} + [- \frac{M λ}{4} + M ϵ (c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ})] \\ + [M ϵ (c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}) - \frac{ρ}{4 c^{θ + 2}}] y^{θ + 2} \\ - \frac{ρ}{4} x^{θ + 2} + (- \frac{M λ}{2} + R_{1}), \end{aligned}$ where (16) $R_{1} = sup_{(x, y) \in R_{+}^{2}} {- \frac{ρ}{4} (x^{θ + 2} + \frac{y^{θ + 2}}{c^{θ + 2}})} + F = F .$ (16) By $M = (2 / λ) max {2, R_{1}}$ , one can see that $M λ / 4 \geq 1$ . Hence we have that for all $(x, y) \in U_{ϵ}^{1}$ , $L V (x, y) \leq - \frac{M λ}{4} - \frac{ρ}{4} x^{θ + 2} \leq - \frac{M λ}{4} \leq - 1,$ where (11) and (12) have been used.

Case 2. On the domain $U_{ϵ}^{2}$ , we have $x y \leq ϵ (1 + x^{θ + 2})$ . It then follows from (10) that $\begin{aligned} L V & \leq - \frac{M λ}{4} + [- \frac{M λ}{4} + M ϵ (c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ})] \\ + [M ϵ (c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}) - \frac{ρ}{4}] x^{θ + 2} \\ - \frac{ρ}{4 c^{θ + 2}} y^{θ + 2} + (- \frac{M λ}{2} + R_{1}) . \end{aligned}$ In view of (11) and (13), we get $L V (x, y) \leq - \frac{M λ}{4} - \frac{ρ}{4 c^{θ + 2}} y^{θ + 2} \leq - \frac{M λ}{4} \leq - 1.$

Case 3. On $U_{ϵ}^{3}$ , we have from (14) that for all $(x, y) \in U_{ϵ}^{3}$ , $\begin{aligned} L V & \leq - M λ - \frac{ρ}{4} x^{θ + 2} + M [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y - \frac{ρ}{4} x^{θ + 2} \\ - \frac{ρ}{2 c^{θ + 2}} y^{θ + 2} + F \\ \leq - M λ - \frac{ρ}{4} x^{θ + 2} + R_{2} \leq - M λ - \frac{ρ}{4 ϵ^{θ + 2}} + R_{2} \leq - 1, \end{aligned}$ where (17) $R_{2} = sup_{(x, y) \in R_{+}^{2}} {M [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y - \frac{ρ}{4} x^{θ + 2} - \frac{ρ}{2 c^{θ + 2}} y^{θ + 2} + F} .$ (17)

Case 4. On $U_{ϵ}^{4}$ , we have from (15) that for all $(x, y) \in U_{ϵ}^{4}$ , $\begin{aligned} L V & \leq - M λ - \frac{ρ}{4 c^{θ + 2}} y^{θ + 2} + M [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y - \frac{ρ}{2} x^{θ + 2} \\ - \frac{ρ}{4 c^{θ + 2}} y^{θ + 2} + F \\ \leq - M λ - \frac{ρ}{4 c^{θ + 2}} \frac{1}{ϵ^{θ + 2}} + R_{3} \leq - 1, \end{aligned}$ where (18) $R_{3} = sup_{(x, y) \in R_{+}^{2}} {M [c_{1} β (1 - m) + \frac{c_{2} c β^{2} (1 - m)^{2} α}{b γ}] x y - \frac{ρ}{2} x^{θ + 2} - \frac{ρ}{4 c^{θ + 2}} y^{θ + 2} + F} .$ (18) To sum up, we can take a sufficiently small ϵ such that $L V (x, y) \leq - 1 f o r a l l (x, y) \in R_{+}^{2} ∖ U_{ϵ},$ meaning that condition (B.2) in Lemma 2.2 holds.

In addition, we find $M_{0} = min_{(x, y) \in U_{ϵ}} {σ_{1}^{2} x^{2}, σ_{2}^{2} y^{2}}$ such that $\sum_{i, j = 1}^{2} a_{i j} (x) ξ_{i} ξ_{j} = σ_{1}^{2} x^{2} ξ_{1}^{2} + σ_{2}^{2} y^{2} ξ_{2}^{2} \geq M_{0} {∥ξ∥}^{2}, x \in U_{ϵ}, ξ = (ξ_{1}, ξ_{2}) \in R^{2},$ which implies that condition (B.1) in Lemma 2.2 also holds. Therefore, thanks to Lemma 2.2, model (2) has a unique stationary distribution $μ (\cdot)$ and it is ergodic. This completes the proof.

4. Extinction

In this section, we derive sufficient criteria for the extinction of model (2) in two scenarios, one is both the predator and the prey go extinct, the other is the predator goes to extinction but the prey persists.

Theorem 4.1

Let $(x (t), y (t))$ be the solution of model (2) with any positive initial value $(x (0), y (0)) \in R_{+}^{2}$ . We have

if $α < σ_{1}^{2} / 2$ , then both $x (t)$ and $y (t)$ will be extinct, a.s., that is, $lim_{t \to \infty} x (t) = 0, lim_{t \to \infty} y (t) = 0, a . s .$
if $α > σ_{1}^{2} / 2$ , $c β (1 - m) ((α - (σ_{1}^{2} / 2)) / (b)) - γ - (σ_{2}^{2} / 2) < 0$ , then $y (t)$ will go to extinction, a.s., i.e. $lim_{t \to \infty} y (t) = 0, a . s .$ Moreover, the distribution of $x (t)$ converges weakly to the measure which has the following density: $π (x) = C σ_{1}^{- 2} x^{(2 α / σ_{1}^{- 2} 2)} e^{- (2 b / σ_{1}^{2}) x}, x \in (0, \infty),$ where $C = (σ_{1}^{2} (σ_{1}^{2} / 2 b)^{1 - (2 α / σ_{1}^{2})}) / (Γ ((2 α / σ_{1}^{2}) - 1))$ is a constant such that $\int_{0}^{\infty} π (x) d x = 1.$

Proof.

$(i)$ Making use of Itô's formula on $\ln x (t)$ , we get $\begin{aligned} d \ln x (t) & = (\frac{α}{1 + K y (t)} - b x (t) - \frac{β (1 - m) y (t)}{1 + a (1 - m) x (t)} - \frac{σ_{1}^{2}}{2}) d t + σ_{1} d B_{1} (t) \\ \leq (α - \frac{σ_{1}^{2}}{2}) d t + σ_{1} d B_{1} (t), \end{aligned}$ from which we have $\frac{\ln x (t)}{t} \leq \frac{\ln x (0)}{t} + α - \frac{σ_{1}^{2}}{2} + \frac{σ_{1} B_{1} (t)}{t} .$ Noticing that $lim_{t \to \infty} B_{1} (t) / t = 0$ , a.s. and $α < σ_{1}^{2} / 2$ , we have $\underset{t \to \infty}{lim sup} \frac{\ln x (t)}{t} \leq α - \frac{σ_{1}^{2}}{2} < 0,$ which means that (19) $lim_{t \to \infty} x (t) = 0, a . s .$ (19) Similarly, by using Itô's formula to $\ln y (t)$ , we obtain $\begin{aligned} d \ln y (t) & = (- γ + \frac{c β (1 - m) x (t)}{1 + a (1 - m) x (t)} - \frac{σ_{2}^{2}}{2}) d t + σ_{2} d B_{2} (t) \\ \leq (c β (1 - m) x (t) - γ - \frac{σ_{2}^{2}}{2}) d t + σ_{2} d B_{2} (t) . \end{aligned}$ Then, (20) $\frac{\ln y (t)}{t} \leq \frac{\ln y (0)}{t} + \frac{c β (1 - m)}{t} \int_{0}^{t} x (s) d s - γ - \frac{σ_{2}^{2}}{2} + \frac{σ_{2} B_{2} (t)}{t} .$ (20) It follows from (19) and $lim_{t \to \infty} B_{2} (t) / t = 0$ , a.s., we get $\underset{t \to \infty}{lim sup} \frac{\ln y (t)}{t} \leq - γ - \frac{σ_{2}^{2}}{2} < 0,$ which implies that $lim_{t \to \infty} y (t) = 0,$ a.s.

$(i i)$ For any positive initial value $(x (0), y (0)) \in R_{+}^{2}$ , we have $d x (t) \leq x (t) (α - b x (t)) d t + σ_{1} x (t) d B_{1} (t),$ then by Lemma A.1 in Ref. [11] and the comparison theorem, we obtain (21) $\underset{t \to \infty}{lim sup} \frac{1}{t} \int_{0}^{t} x (s) d s \leq \frac{α - (σ_{1}^{2} / 2)}{b} a . s .$ (21) Combining (20) and (21), it follows that if $α > σ_{1}^{2} / 2$ , $c β (1 - m) ((α - (σ_{1}^{2} / 2)) / b) - γ - (σ_{2}^{2} / 2) < 0$ , then we have that $\underset{t \to \infty}{lim sup} \frac{\ln y (t)}{t} \leq c β (1 - m) (\frac{α - (σ_{1}^{2} / 2)}{b}) - γ - \frac{σ_{2}^{2}}{2} < 0,$ that is, $lim_{t \to \infty} y (t) = 0$ , a.s.

On the other hand, consider the following one-dimensional stochastic differential equation: (22) $d X (t) = X (t) (α - b X (t)) d t + σ_{1} X (t) d B_{1} (t),$ (22) with the initial value $X (0) = x (0) > 0$ . We know from Ref. [15] that model (22) has the ergodic property and the invariant density is given by $π (x) = C σ_{1}^{- 2} x^{(2 α / σ_{1}^{2}) - 2} e^{- (2 b / σ_{1}^{2}) x}, x \in (0, \infty),$ where $C = (σ_{1}^{2} (σ_{1}^{2} / 2 b)^{1 - (2 α / σ_{1}^{2})}) / (Γ ((2 α / σ_{1}^{2}) - 1))$ is a constant satisfying that $\int_{0}^{\infty} π (x) d x = 1.$

Since $lim_{t \to \infty} y (t) = 0$ , a.s., for any small ϵ, there are T and a set $Ω_{ϵ} \subset Ω$ such that $P (Ω_{ϵ}) > 1 - ϵ$ , $(β (1 - m) x y) / (1 + a (1 - m) x) \leq β (1 - m) x y \leq ϵ x$ and $α / (1 + K y) \geq α / (1 + K ϵ)$ for $t \geq T$ and $ω \in Ω_{ϵ}$ . Then, we have $\begin{aligned} [x (t) (\frac{α}{1 + K ϵ} - b x (t)) - ϵ x (t)] d t + σ_{1} x (t) d B_{1} (t) \\ \leq d x (t) \leq x (t) (α - b x (t)) d t + σ_{1} x (t) d B_{1} (t), \end{aligned}$ that is to say, when $ϵ \to 0$ , we obtain that the distribution of the process $x (t)$ converges weakly to the measure with the density of $π (x)$ . This completes the proof.

5. Numerical simulations

In this section, we will perform some numerical simulations using MATLAB R2016b to illustrate the effect of white noise, fear effect and a prey refuge on the dynamics of model (2). The numerical scheme obtained through Milstein's higher order method [8] applied on stochastic model (2) under consideration is given by $\begin{aligned} x_{j + 1} & = x_{j} + x_{j} [(\frac{α}{1 + K y_{j}} - b x_{j} - \frac{β (1 - m) y_{j}}{1 + a (1 - m) x_{j}}) Δ t \\ + σ_{1} ζ_{j} \sqrt{Δ t} + \frac{1}{2} σ_{1}^{2} Δ t (ζ_{j}^{2} - 1)], \\ y_{j + 1} & = y_{j} + y_{j} [(\frac{c β (1 - m) x_{j}}{1 + a (1 - m) x_{j}} - γ) Δ t + σ_{2} ξ_{j} \sqrt{Δ t} + \frac{1}{2} σ_{2}^{2} Δ t (ξ_{j}^{2} - 1)], \end{aligned}$ where $ζ_{j}$ and $ξ_{j}$ are two independent Gaussian random variables $N (0, 1)$ for $j = 1, 2, \dots, n$ . All numerical simulations reported here are carried out with the choice of time stepping $Δ t = 0.01$ . In this section, we always take the following parameter values: (23) $α = 0.6, b = 0.3, β = 0.3, c = 0.8, a = 0.3, γ = 0.1,$ (23) and assume that the initial value of model (2) is $(x (0), y (0)) = (0.6, 0.5)$ .

First, in order to investigate the impact of white noise on the dynamics of model (2), we choose all parameter values in (23) and $m = 0.1, K = 0.3, σ_{1} = σ_{2} = 0.01$ . By calculating, we get from (6) that $c_{2} > 0.148$ . Now we take $c_{2} = 0.1481$ , then we compute $λ = 0.0118 > 0$ . Thus Assumptions $(H_{1})$ and $(H_{2})$ are satisfied. By Theorem 3.1, there exists a unique ergodic stationary distribution for the solutions of model (2). Figure 1 confirms this.

Figure 1. (a) and (c): The asymptotic behaviour of the solutions to stochastic model (2) around the positive equilibrium of model (1) with initial value $(x (0), y (0)) = (0.6, 0.5)$ ; (b) and (d): The density function diagrams of $x (t)$ and $y (t)$ , respectively. The parameters are taken as (23) and m = 0.1, K = 0.3, $σ_{1} = σ_{2} = 0.01$ .

Display full size

To obtain deep insights of the influences of white noise on population dynamics, we keep the model parameter values the same as in (23) but let $σ_{1} = 1.1$ , which implies that condition $α < σ_{1}^{2} / 2$ holds. Then we get that both prey and predator populations are extinction by Theorem 4.1, shown as in Figure 2. Now we take $σ_{1} = 0.1, σ_{2} = 0.9$ and the other parameters remain unchanged, we can easily check that the condition $α > σ_{1}^{2} / 2, c β (1 - m) ((α - (σ_{1}^{2} / 2) / b) - γ - \frac{σ_{2}^{2}}{2} < 0$ are satisfied. By Theorem 4.1, we know that the prey population is persistent and the predator population will be extinct (see Figure 3). These show that the random disturbance will change the dynamics when the environment noise is large enough.

Figure 2. Numerical simulation for model (1) and model (2) with initial value $(x (0), y (0)) = (0.6, 0.5)$ . The parameters are taken as (23) and m = 0.1, K = 0.3, $σ_{1} = 1.1$ , $σ_{2} = 0.01$ .

Display full size

Figure 3. Numerical simulation for model (1) and model (2) with initial value $(x (0), y (0)) = (0.6, 0.5)$ . The parameters are taken as (23) and m = 0.1, K = 0.3, $σ_{1} = 0.1$ , $σ_{2} = 0.9$ .

Display full size

Next, in order to illustrate the influence of fear effect to model (2) through numerical simulation, we choose different values of K, say K = 0, K = 0.1 and K = 1, $σ_{1} = σ_{2} = 0.01$ . For the remaining parameter values, we keep them the same as in (23). We can compute that these parameters satisfy Assumptions $(H_{1})$ and $(H_{2})$ and therefore model (2) has a stationary distribution, which means that both prey and predator populations are persistent. Moreover, from Figure 4 we find that the increase of fear effect will reduce the density of predator but with no obvious effect on the density of prey.

Figure 4. Numerical simulation for model (1) and model (2) with initial value $(x (0), y (0)) = (0.6, 0.5)$ and different K, respectively. The parameters are taken as (23) and m = 0.1, $σ_{1} = σ_{2} = 0.01$ .

Display full size

Finally, we numerically illustrate the impact of a prey refuge to model (2) and choose different values of m, say m = 0, m = 0.2, m = 0.5 and m = 0.7, $σ_{1} = σ_{2} = 0.01$ . For the remaining parameter values, we keep them the same as in (23). We can easily compute that these parameters satisfy assumptions $(H_{1})$ and $(H_{2})$ . So we obtain that there is a unique ergodic stationary distribution for model (2), which means that both prey and predator populations are persistent. From Figure 5, we can see that with the increase of the amount of refuges m, the density of the prey increases while the predator density decreases. This means that the increase of the amount of refuges m is favourable to prey population but it is adverse to the persistence of predator population. Moreover, we can also see that large K may lead to the predator density decreasing fast. When we choose $σ_{1} = 1.1, σ_{2} = 0.01$ and the other parameters are the same as in Figure 5, we simply calculate that condition $α < σ_{1}^{2} / 2$ holds. Then we obtain that both prey and predator populations are extinction by Theorem 4.1, shown in Figure 6. The results again show that the large noise is very destructive to the persistence of the population.

Figure 5. The solutions of model (2) with the initial value $(x (0), y (0)) = (0.6, 0.5)$ and different K,m, respectively. The parameters are taken as (23) and $σ_{1} = σ_{2} = 0.01$ .

Display full size

Figure 6. The solutions of model (2) with the initial value $(x (0), y (0)) = (0.6, 0.5)$ and different K,m, respectively. The parameters are taken as (23) and $σ_{1} = 1.1, σ_{2} = 0.01$ .

Display full size

6. Conclusion

This paper focuses on a stochastic Holling type-II prey–predator model with a prey refuge and fear effect. Mathematically, the sufficient criteria for the existence of a unique ergodic stationary distribution and the extinction of the model have been obtained. Ecologically, we get the following conclusions:

Random disturbance may change the dynamic behaviour of the population, especially when the noise is large, it may lead to the extinction of the prey and predator populations.
By simulations, we find that the increase of fear effect K will reduce the density of predator, but the effect on the density of prey is not obvious. In addition, one can see that the increase of the amount of refuges m is favourable to prey population but it is adverse to the persistence of predator population.

There are some interesting themes worthy of further research. On the one hand, we can incorporate some other environmental noises into model (2), such as coloured noise, Poisson noise and so on. On the other hand, to make the model (2) more realistic, we can further consider the factors such as the influence of impulsive perturbations and delay. We leave these for future investigations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Research is supported by the National Natural Science Foundation of China (Grant Nos. 11671260 and 12071293).

References

J.R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, J. Anim. Ecol. 44 (1975), pp. 331–340. [Crossref], [Web of Science ®], [Google Scholar]
F. Brauer and C. Castillo-Chvez, Mathematical Models in Population Biology and Epidemiology, Springer, Berlin, 2001. [Crossref], [Google Scholar]
S. Creel and D. Christianson, Relationships between direct predation and risk effects, Trends Ecol. Evol. 23 (2008), pp. 194–201. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
W. Cresswell, Predation in bird populations, J. Ornithol. 152 (2011), pp. 251–263. [Crossref], [Web of Science ®], [Google Scholar]
D.L. DeAngelis, R.A. Goldstein, and R.V. O'Neill, A model for tropic interaction, Ecology, 56 (1975), pp. 881–892. [Crossref], [Web of Science ®], [Google Scholar]
R. Durrett, Stochastic Calculus: A Practical Introduction, Probability Stochastics, 1st ed., CRC Press, Boca Raton, FL, 1996. [Google Scholar]
R.Z. Has'minskii, Stochastic Stability of Differential Equations, Sijthoff Noordhoff, Alphen aan den Rijn, The Netherlands, 1980. [Crossref], [Google Scholar]
D. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001), pp. 525–546. [Crossref], [Web of Science ®], [Google Scholar]
C.S. Holling, The functional response of predators to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Can. 97 (1965), pp. 5–60. [Crossref], [Google Scholar]
Y. Huang, F. Chen, and L. Zhong, Stability analysis of a prey–predator model with Holling type III response function incorporating a prey refuge, Appl. Math. Comput. 182 (2006), pp. 672–683. [Crossref], [Web of Science ®], [Google Scholar]
C. Ji and D. Jiang, Qualitative analysis of a stochastic ratio-dependent predator–prey system, J. Comput. Appl. Math. 235 (2011), pp. 1326–1341. [Crossref], [Web of Science ®], [Google Scholar]
D. Jia, T. Zhang, and S. Yuan, Pattern dynamics of a diffusive toxin producing phytoplankton–zooplankton model with three-dimensional patch, Int. J. Bifurcat. Chaos 29 (2019), pp. 1930011. [Crossref], [Web of Science ®], [Google Scholar]
S.L. Lima, Nonlethal effects in the ecology of predator–prey interactions, Bioscience 48 (1998), pp. 25–34. [Crossref], [Web of Science ®], [Google Scholar]
M. Liu and M. Deng, Permanence and extinction of a stochastic hybrid model for tumor growth, Appl. Math. Lett. 94 (2019), pp. 66–72. [Crossref], [Web of Science ®], [Google Scholar]
Q. Liu, D. Jiang, T. Hayat, and A. Alsaedi, Dynamics of a stochastic predator-prey model with stage structure for predator and Holling type II functional response, J. Nonlinear Sci. 28 (2018), pp. 1151–1187. [Crossref], [Web of Science ®], [Google Scholar]
C. Liu, L. Wang, Q. Zhang, and Y. Yan, Dynamical analysis in a bioeconomic phytoplankton zooplankton system with double time delays and environmental stochasticity, Physica A 15 (2017), pp. 682–698. [Crossref], [Google Scholar]
M. Liu and Y. Zhu, Stability of a budworm growth model with random perturbations, Appl. Math. Lett. 79 (2018), pp. 13–19. [Crossref], [Web of Science ®], [Google Scholar]
Z. Ma, F. Chen, C. Wu, and W. Chen, Dynamic behaviors of a Lotka-Volterra predator–prey model incorporating a prey refuge and predator mutual interference, Appl. Math. Comput. 219 (2013), pp. 7945–7953. [Crossref], [Web of Science ®], [Google Scholar]
Z. Ma, S. Wang, W. Li, and Z. Li, The effect of prey refuge in a patchy predator–prey system, Math. Biosci. 243 (2013), pp. 126–130. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
X. Mao, Stochastic Differential Equations and Applications, Woodhead Publishing, Cambridge, 2007. [Google Scholar]
X. Mao, G. Marion, and E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stochast. Proc. Appl. 97 (2002), pp. 95–110. [Crossref], [Google Scholar]
X. Meng, F. Li, and S. Gao, Global analysis and numerical simulations of a novel stochastic eco-epidemiological model with time delay, Appl. Math. Comput. 339 (2018), pp. 701–726. [Crossref], [Web of Science ®], [Google Scholar]
D. Mukherjee, The effect of prey refuges on a three species food chain model, Differ. Equ. Dyn. Syst.22 (2014), pp. 413–426. [Crossref], [Google Scholar]
J. Ripa, P. Lundberg, and V. Kaitala, A general theory of environmental noise in ecological food webs, Am. Nat. 151 (1998), pp. 256–263. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
S. Sarwardi, P.K. Mandal, and S. Ray, Analysis of a competitive prey–predator system with a prey refuge, Biosystems 110 (2012), pp. 133–148. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
S. Sasmal, Population dynamics with multiple Allee effects induced by fear factors induced by fear factors – a mathematical study on prey–predator, Appl. Math. Model. 64 (2018), pp. 1–14. [Crossref], [Web of Science ®], [Google Scholar]
A. Sih, Prey refuges and predator–prey stability, Theor. Popul. Biol. 31 (1987), pp. 1–12. [Crossref], [Web of Science ®], [Google Scholar]
X. Wang, L. Zanette, and X. Zou, Modelling the fear effect in predator–prey interactions, J. Math. Biol. 73 (2016), pp. 1–38. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
C. Xu and S. Yuan, Competition in the chemostat: a stochastic multi-species model and its asymptotic behavior, Math. Biosci. 280 (2016), pp. 1–9. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
C. Xu, S. Yuan, and T. Zhang, Global dynamics of a predator–prey model with defence mechanism for prey, Appl. Math. Lett. 62 (2016), pp. 42–48. [Crossref], [Web of Science ®], [Google Scholar]
C. Xu, S. Yuan, and T. Zhang, Average break-even concentration in a simple chemostat model with telegraph noise, Nonlinear Anal. Hybrid. Syst. 29 (2018), pp. 373–382. [Crossref], [Web of Science ®], [Google Scholar]
S. Yan, D. Jia, T. Zhang, and S. Yuan, Pattern dynamic in a diffusive predator–prey model with hunting cooperations, Chaos Solitons Fract. 130 (2020), pp. 109428. [Crossref], [Web of Science ®], [Google Scholar]
F. Yi, J. Wei, and J. Shi, Bifurcation and spatiotemporal patterns in a homogeneous diffusive predator–prey system, J. Differ. Equ. 246 (2009), pp. 1944–1977. [Crossref], [Web of Science ®], [Google Scholar]
X. Yu and S. Yuan, Asymptotic properties of a stochastic chemostat model with two distributed delays and nonlinear perturbation, Discrete Contin. Dyn. B 25 (2020), pp. 2373–2390. [Web of Science ®], [Google Scholar]
X. Yu, S. Yuan, and T. Zhang, Asymptotic properties of stochastic nutrient-plankton food chain models with nutrient recycling, Nonlinear Anal. Hybrid. Syst. 34 (2019), pp. 209–225. [Crossref], [Web of Science ®], [Google Scholar]
S. Yuan, D. Wu, G. Lan, and H. Wang, Noise-induced transitions in a nonsmooth predator–prey model with stoichiometric constraints, Bull. Math. Biol. 82 (2020), pp. 167. [Crossref], [Web of Science ®], [Google Scholar]
H. Zhang, Y. Cai, S. Fu, and W. Wang, Impact of the fear effect in a prey–predator model incorporating a prey refuge, Appl. Math. Comput. 356 (2019), pp. 328–337. [Crossref], [Web of Science ®], [Google Scholar]
S. Zhang and X. Meng, Dynamics analysis and numerical simulations of a stochastic non-autonomous predator–prey system with impulsive effects, Nonlinear Anal. Hybrid. Syst. 26 (2017), pp. 19–37. [Crossref], [Web of Science ®], [Google Scholar]
Y. Zhao, S. Yuan, and J. Ma, Survival and stationary distribution analysis of a stochastic competitive model of three species in a polluted environment, Bull. Math. Biol. 77 (2015), pp. 1285–1326. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
S. Zhao, S. Yuan, and H. Wang, Threshold behavior in a stochastic algal growth model with stoichiometric constraints and seasonal variation, J. Differ. Equ. 268 (2020), pp. 5113–5139. [Crossref], [Web of Science ®], [Google Scholar]
W. Zuo and D. Jiang, Stationary distribution and periodic solution for stochastic predator–prey systems with nonlinear predator harvesting, Commun. Nonlinear Sci. 36 (2016), pp. 65–80. [Crossref], [Web of Science ®], [Google Scholar]

Appendix 1: Construction of the stochastic model (2)

By the methods used in Ref. [6] about stochastic effects, we first consider a discrete time Markov chain. For $x (t)$ and $y (t)$ in model (2), given $Δ t > 0$ , we present a process $X^{(Δ t)} (t) = (x^{(Δ t)} (t), y^{(Δ t)} (t))^{T}, t = 0, Δ t, 2 Δ t, \dots$ with initial value $X^{(Δ t)} (0) = (x^{(Δ t)} (0), y^{(Δ t)} (0))^{T} \in R_{+}^{2}$ . Let normal distribution random variable sequence ${R_{i}^{(Δ t)} (l)}_{l = 0}^{\infty}$ satisfies (A1) $E [R_{i}^{(Δ t)} (l)] = 0, E [R_{i}^{(Δ t)} (l)]^{2} = σ_{i}^{2} Δ t, E [R_{i}^{(Δ t)} (l)]^{4} = o (Δ t),$ (A1) where $i = 1, 2,$ $l = 0, 1, 2, \dots$ and $σ_{i}^{2}$ express the intensities of stochastic effects. On each period $[l Δ t, (l + 1) Δ t)$ , we assume that $X^{(Δ t)} (t)$ grows according to model (2) and is also affected by the random $(x^{(Δ t)} R_{1}^{(Δ t)} (l), y^{(Δ t)} R_{2}^{(Δ t)} (l))^{T}$ . Specifically, for $l = 0, 1, \dots$ , we get $\begin{aligned} x^{(Δ t)} ((l + 1) Δ t) & = x^{(Δ t)} (l Δ t) + x^{(Δ t)} (l Δ t) R_{1}^{(Δ t)} (l) \\ + [x^{(Δ t)} (l Δ t) (\frac{α}{1 + K y^{(Δ t)} (l Δ t)} - b x^{(Δ t)} (l Δ t) \\ - \frac{β (1 - m) y^{(Δ t)} (l Δ t)}{1 + a (1 - m) x^{(Δ t)} (l Δ t)})] Δ t, \\ y^{(Δ t)} ((l + 1) Δ t) & = y^{(Δ t)} (l Δ t) + y^{(Δ t)} (l Δ t) R_{2}^{(Δ t)} (l) \\ + [- γ y^{(Δ t)} (l Δ t) + \frac{c β (1 - m) x^{(Δ t)} (l Δ t) y^{(Δ t)} (l Δ t)}{1 + a (1 - m) x^{(Δ t)} (l Δ t)}] Δ t . \end{aligned}$ Next we demonstrate that $X^{(Δ t)} (t)$ converges to a diffusion process as $Δ t \to 0$ . And we need to determine the drift coefficient and diffusion coefficient of the diffusion process. Let $P^{(Δ t)} (\bar{x}, d z)$ denote the transition probabilities of the homogeneous Markov chain ${X^{(Δ t)} (l Δ t)}$ , that is $P^{(Δ t)} (\bar{x}, B) = P r o b {X^{(Δ t)} ((l + 1) Δ t) \in B | X^{(Δ t)} = \bar{x}}$ for all $\bar{x} = (x, y) \in R_{+}^{2}$ and all Borel sets $B \subset R_{+}^{2}$ . From (A1) $\begin{aligned} \frac{1}{Δ t} \int (z_{1} - x) P^{(Δ t)} (\bar{x}, d z) & = x (\frac{α}{1 + k y} - b x - \frac{β (1 - m) y}{1 + a (1 - m) x}) + \frac{x}{Δ t} E R_{1}^{(Δ t)} (0) \\ = x (\frac{α}{1 + k y} - b x - \frac{β (1 - m) y}{1 + a (1 - m) x}), \\ \frac{1}{Δ t} \int (z_{2} - y) P^{(Δ t)} (\bar{x}, d z) & = - γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x} + \frac{y}{Δ t} E R_{2}^{(Δ t)} (0) \\ = - γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x} . \end{aligned}$ To determine the diffusion coefficients, we consider the following moments: $g_{i j}^{(Δ t)} (x) = \frac{1}{Δ t} \int (z_{i} - x_{i}) (z_{j} - x_{j}) P^{(Δ t)} (x, d z), i, j = 1, 2.$ By (A1) $\begin{aligned} |g_{11}^{(Δ t)} (x) - σ_{1}^{2} x^{2}| & = |\frac{1}{Δ t} E {Δ t [x (\frac{α}{1 + K y} - b x - \frac{β (1 - m) y}{1 + a (1 - m) x})] + R_{1}^{(Δ t)} (0) x}^{2} - σ_{1}^{2} x^{2}| \\ = | Δ t {[x (\frac{α}{1 + K y} - b x - \frac{β (1 - m) y}{1 + a (1 - m) x})]}^{2} \\ + 2 x^{2} (\frac{α}{1 + K y} - b x - \frac{β (1 - m) y}{1 + a (1 - m) x}) E R_{1}^{(Δ t)} (0) \\ + \frac{1}{Δ t} x^{2} E {[R_{1}^{(Δ t)} (0)]}^{2} - σ_{1}^{2} x^{2} | \\ = Δ t {[x (\frac{α}{1 + K y} - b x - \frac{β (1 - m) y}{1 + a (1 - m) x})]}^{2}, \\ |g_{22}^{(Δ t)} (x) - σ_{2}^{2} y^{2}| & = |\frac{1}{Δ t} E {Δ t [- γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x}] + R_{2}^{(Δ t)} (0) y}^{2} - σ_{2}^{2} y^{2}| \\ = | Δ t {[- γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x}]}^{2} + 2 y^{2} [- γ + \frac{c β (1 - m) x}{1 + a (1 - m) x}] E R_{2}^{(Δ t)} (0) \\ + \frac{1}{Δ t} y^{2} E {[R_{2}^{(Δ t)} (0)]}^{2} - σ_{2}^{2} y^{2} | \\ = Δ t {[- γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x}]}^{2} . \end{aligned}$ Therefore, $lim_{Δ t \to 0^{+}} sup_{| | x | | \leq L} |g_{11}^{(Δ t)} (x) - σ_{1}^{2} x^{2}| = 0, lim_{Δ t \to 0^{+}} sup_{| | x | | \leq L} |g_{22}^{(Δ t)} (x) - σ_{2}^{2} y^{2}| = 0$ for all $0 < L < \infty$ . Similarly, $lim_{Δ t \to 0^{+}} sup_{| | x | | \leq L} |g_{12}^{(Δ t)} (x)| = 0.$ In addition, from (A1), we obtain that for all $0 < L < \infty$ $lim_{Δ t \to 0^{+}} sup_{| | x | | \leq L} \frac{1}{Δ t} \int | | z - x | |^{3} P^{(Δ t)} (x, d z) = 0.$ Finally, extend the definition of $X^{(Δ t)} (t)$ to all $t \geq 0$ by setting $X^{(Δ t)} (t) = X^{(Δ t)} (l Δ t)$ for all $t \in [l Δ t, (l + 1) Δ t)$ . On the basis of Theorem 7.1, Lemma 8.2 in Ref. [6], we can verify that as $t \to 0$ , $X^{(Δ t)} (t)$ converges weakly to the solution of the stochastic system as follows: $\begin{aligned} d x & = [\frac{α x}{1 + K y} - b x^{2} - \frac{β (1 - m) x y}{1 + a (1 - m) x}] d t + σ_{1} x d B_{1} (t), \\ d y & = [- γ y + \frac{c β (1 - m) x y}{1 + a (1 - m) x}] d t + σ_{2} y d B_{2} (t), \end{aligned}$ where $B_{1} (t)$ and $B_{2} (t)$ are two independent Brownian motions with the intensities represented by two positive constants $σ_{1}$ and $σ_{2}$ .

Journal of Biological Dynamics

Survival analysis of a stochastic predator–prey model with prey refuge and fear effect

2019 Guangzhou Workshop

Survival analysis of a stochastic predator–prey model with prey refuge and fear effect

Abstract
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

1. Introduction

2. Preliminaries

[7]

3. Stationary distribution

4. Extinction

5. Numerical simulations

Published online:

Published online:

Published online:

Published online:

Published online:

Published online:

6. Conclusion

Disclosure statement

Funding

References

Related research

Information for

Open access

Opportunities

Help and information

Journal of Biological Dynamics

Survival analysis of a stochastic predator–prey model with prey refuge and fear effect

2019 Guangzhou Workshop

Survival analysis of a stochastic predator–prey model with prey refuge and fear effect

AbstractFormulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

1. Introduction

2. Preliminaries

[7]

3. Stationary distribution

4. Extinction

5. Numerical simulations

Published online:

Published online:

Published online:

Published online:

Published online:

Published online:

6. Conclusion

Disclosure statement

Additional information

Funding

References

Related research

Information for

Open access

Opportunities

Help and information

Keep up to date

Abstract
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.