A diffusive predator-prey system with prey refuge and predator cannibalism

1.   Introduction

Predator-prey systems as one of the most important relationships between two populations have attracted the widespread attention and been extensively studied in both ecology and mathematical ecology. Based on ODE systems and PDE systems, various mathematical models have been built to understand and investigate predator-prey interaction. We refer the reader to the references ^{[1,2,3,4,5,6,7]} and references therein.

Cannibalism, defined more specifically as the killing and at least partial consumption of conspecifics, is widespread in nature ^[8,9]. It has been observed that cannibalism exists in different types of animals, such as, insects, fishes, zooplankton, isopods and amphibians. For example, in aquatic ecosystems, Shevtsova et. al. ^[11] have showed that adult Dreissena can feed on many small zooplankton species including rotifers, polyarthra vulgaris, protozoans, and cyclopoid copepopids. Chakraborty and Chattopadhyay ^[12] pointed out that the phenomenon of sexual cannibalism is very common in many families of spiders and scorpions. For more examples of cannibalism, please see references ^[13,14]. Cannibalism leads to a trophic structure and feedback loops within a population, and then it has a strong impact on population structure and dynamics. It is well explored in mathematical literatures that cannibalism can have either a stabilizing or a destabilizing effect on predator-prey systems ^{[10,12,15,16,17,18,19]}.

In order to preserve biodiversity and avoid species extinction, an effective strategy is to establish a refuge or a protection zone. In predator-prey interactions, prey species can exhibit spatial refugia which afford the prey some degree of protection from predation ^[20]. For example, Huffaker and Kennett ^[21] showed that cyclamen mites can use strawberry plants as physical barriers to avoid predation by Typhlodromus mites. Previous studies have shown that refugia have a stabilizing effect on prey-predator systems with different functional responses ^{[22,23,24,25]}. In the case of spatial distribution patterns and dispersal mechanisms, Du and Shi first in ^[26] investigated dynamics of a reaction-diffusion predator-prey system with a protection zone for the prey. In ^{[27,28,29,30,31,32,33]}, authors also studied the effect of a prey refuge or a protection zone in the diffusive predator-prey system.

Motivated by the existing studies and the above considerations, in this study, we consider the following diffusive predator-prey system with prey refuge and predator cannibalism

Here $\Omega$ is a bounded domain in ${\mathbb{R}}^n (n\geq1)$ with smooth boundary $\partial \Omega$ and $c, e_1, e_2\in[0, 1)$. All the variables and parameters of system (1.1) and their biological significance are listed in Table 1. When the spatial distribution is homogeneous and $d_u = d_v = c = 0$, system (1.1) reduces to an ODE system

In ^[18], Kohlmeier and Ebenhöh established the existence and stability of steady states of (1.2) and proved that cannibalism can have a stabilizing effect. Chakraborty and Chattopadhyay ^[12] showed that the paradox of enrichment does not hold for a higher cannibalism rate among predators for system (1.2). In ^[34], Prasad and Prasad gave the existence and stability of equilibria and analysed the existence of bifurcations for system (1.2) with provision of additional food.

There is increasing recognition that the understanding of patterns and mechanisms of spatial dispersal is a significant issue in the study of predator-prey system. Spatial heterogeneity can make predator-prey system exhibit more complex dynamic properties. Considering the effect of spatial diffusion coefficient on the dynamical properties of system (1.1) is the first research topic in the present paper. In view of the widespread existence of cannibalism, it is an interesting problem is to explore how cannibalism affects predator-prey systems. In addition, from the perspective of protecting biodiversity, we also discuss the effects of prey refuge.

The rest of the paper is organized as follows. In Section 2, we establish the global existence, dissipation and persistence of positive solutions of system (1.1). In Sections 3 and 4, we investigate the local and global stability of constant steady states, Turing instability, and nonexistence and existence of nonconstant steady state solutions. In Section 5, we consider the influence of prey refuge and predator cannibalism on predator and prey biomass density by using a systematic sensitivity analysis. In the discussion section, we summary our findings and state some biologically motivated mathematical questions for future study. Throughout this paper, numerical simulations under reasonable parameter values from literatures are presented to illustrate or complement our mathematical findings.

2.   Global existence, dissipation and persistence

This section is devoted to investigating global existence, dissipation and persistence of positive solutions of system (1.1).

Theorem 2.1. System (1.1) has a unique global solution $(u(x, t), v(x, t))$ such that $u(x, t) > 0$ and $v(x, t) > 0$ for $(x, t)\in \bar{ \Omega}\times(0, \infty)$.

Proof. It is clear that (1.1) is a mixed quasi-monotone system for the domain $\{u\geq0, v\geq0\}$. Let $(\underline{u}(x, t), \underline{v}(x, t)) = (0, 0)$ and $(\bar{u}(x, t), \bar{v}(x, t)) = (\bar{u}(t), \bar{v}(t))$, where $(\bar{u}(t), \bar{v}(t))$ satisfies

It follows from the existence and uniqueness theorem of solutions of ordinary differential equations that $(\bar{u}(t), \bar{v}(t))$ is global existence and $\bar{u}(t) > 0, \bar{v}(t) > 0$ for $t\geq0$. Note that

then from comparison principle of the parabolic equations, it is easy to verify that $u(x, t)\leq\bar{u}(t)$. Similarly, by $v_0(x)\leq \bar{v}_0$, we have $v(x, t)\leq\bar{v}(t)$. Then $(\bar{u}(x, t), \bar{v}(x, t))$ and $(\underline{u}(x, t), \underline{v}(x, t))$ are the coupled ordered upper and lower solutions of system (1.1). This means that there is a unique global solution $(u(x, t), v(x, t))$ satisfying

Moreover, by the strong maximum principle we see that $u(x, t) > 0$ and $v(x, t) > 0$ for $(x, t)\in \bar{ \Omega}\times(0, \infty)$.

Theorem 2.2. If $(u, v)$ is any solution of system (1.1), then

Proof. It is clear that

It follows from comparison principle of parabolic equations that the first inequality of (2.1) holds. This means that for any $\epsilon > 0$ there exists $T_1 > 0$ such that $u(x, t)\leq K+\epsilon$ for all $x\in \bar{\Omega}$ and $t\geq T_1$. Then

with boundary value $\partial v/ \partial n = 0, x\in\partial\Omega, t > T_1$ and initial value $v(x, T_1) > 0, x\in\bar{\Omega}$. Let $z_1(t)$ be a solution of

with $z_1(T_1) = \max\limits_{\bar{\Omega}}v(\cdot, T_1) > 0$. Note that

From the comparison principle, we conclude that the second inequality of (2.1) holds. This completes the proof.

Remark 2.1. It follows from Theorem 2.2 that

is a global attractor of (1.1) in $\mathbb{R}_+^2$ for any $\epsilon > 0$.

Theorem 2.3. If

then system (1.1) is {persistent}, that is,

Proof. It follows from the first equation of (1.1) that

From comparison principle of parabolic equations, the first inequality of (3.2) holds. Then for any $\epsilon > 0$ there is $T_2 > 0$ such that $u(x, t)\geq K(r\eta-a(1-c))/(r\eta)-\epsilon: = A$ for all $x\in \bar{\Omega}$ and $t\geq T_2$. By the second equation of (1.1), we have

Note that if $z_2(t)$ is a solution of

with $z_2(T_2) = \min\limits_{\bar{\Omega}}v(\cdot, T_2) > 0$, then

since (2.2) holds. This proves that the second inequality of (2.3) holds. The proof is completed.

3.   Analysis of constant steady states

In this section, we investigate the existence and stability of constant steady states of system (1.1). The constant steady states of (1.1) are listed below: the extinct steady state $E_0:(0, 0)$; the predator-extinction steady state $E_1:(K, 0)$; the coexistence steady state $E_2:(\bar{u}, \bar{v})$. To establish the stability of the above constant steady states of (1.1), we first make some notations. It is well-known that the operator $- \Delta$ in $\Omega$ with the homogeneous Neumann boundary condition has eigenvalues

where ${\mathbb{N}}_0: = {\mathbb{N}} \bigcup \{0\}$. Let $S(\mu_i)$ be the subspace generated by the eigenfunctions $\phi_{ij}$ corresponding to $\mu_i$, $m(\mu_i)$ be the multiplicity of $\mu_i$, and $\{\phi_{ij}\}_{j = 1}^{m(\mu_i)}$ be an orthonormal basis of $S(\mu_i)$. Define $X_{ij} = \{c\phi_{ij}:c\in {\mathbb{R}}^2\}$, $X_i = \bigoplus_{j = 1}^{m(\mu_i)}X_{ij}$, and

satisfying $X = \bigoplus_{i = 0}^{ \infty}X_i$. We linearize the system (1.1) about a constant steady state $(\hat{u}, \hat{v})$ and obtain

with domain $X_{H} = \left\{(\varphi, \psi)\in [C^1(\Omega\times {\mathbb{R}}^+)]^2: \partial{ \varphi}/ \partial\nu = \partial{\psi}/ \partial\nu = 0\right\}$, where

and

$(\hat{u}, \hat{v})$ is locally asymptotically stable if all eigenvalues of the operator $H_{(\hat{u}, \hat{v})}$ have negative real part, and it is unstable if at least one eigenvalue has positive real part. In the following subsections, we will discuss the existence, local stability and global stability of $E_0$, $E_1$ and $E_2$.

3.1.   The extinct steady state and predator-extinction steady state

This subsection focuses on the existence and stability of the extinct steady state $E_0$ and the predator-extinction steady state $E_1$. It is clear that $E_0$ and $E_1$ always exist.

Theorem 3.1. $E_0$ is always unstable with respect to (1.1).

Proof. It follows from (3.1) and (3.3) that the corresponding $k$-th characteristic equation for the linearized system of (1.1) at $E_0$ is

Note that two eigenvalues are $r$ and $-m$ when $k = 0$. This means that $E_0$ is unstable.

Theorem 3.2. If $m > e_1aK(1-c)/(h+K(1-c))$, then $E_1$ is globally asymptotically stable with respect to (1.1).

Proof. From (3.3), we have

It follows from (3.1) that the corresponding $k$-th characteristic equation for the linearized system of (1.1) at $E_1$ is

It is clear that $\lambda < 0$ for any $k\in {\mathbb{N}}_0$ if $m > e_1a(1-c)/(h+K(1-c))$, which implies that $E_1$ is locally asymptotically stable.

From (2.1), we conclude that if $m > e_1aK(1-c)/(h+K(1-c))$, then

This means that $v \rightarrow 0$ uniformly on $\bar{\Omega}$ as $t \rightarrow \infty$. For any $\epsilon > 0$ there exists $T > 0$ such that $v(x, t)\leq \epsilon$ for all $x\in\bar{ \Omega}$ and $t\geq T$. From the first equation of (1.1), we have

Note that if $z(t)$ is a solution of

with $z(T) = \min\limits_{\bar{\Omega}}u(\cdot, T) > 0$, then $\lim_{t\to \infty}z(t) = K$ since $\epsilon$ is sufficiently small. By using comparison principle of parabolic equations, we obtain $\liminf_{t\to\infty}\min_{\bar{\Omega}}u(\cdot, t)\geq K$ since $\epsilon$ is sufficiently small. Combining with the first inequality of (3.4) gives $u \rightarrow K$ uniformly on $\bar{\Omega}$ as $t\to\infty$, which means that $K$ is globally attractive. Hence, $E_1$ is globally asymptotically stable.

3.2.   The coexistence steady state

The interior coexistence steady state $(\bar{u}, \bar{v})$ can be obtained by solving

Let

A direct calculation gives

Note that $\bar{u} > K$ and $\bar{v} < 0$ when $m > e_1aK(1-c)/(h+K(1-c))$, $\bar{u} = K$ and $\bar{v} = 0$ when $m = e_1aK(1-c)/(h+K(1-c))$, $0 < \bar{u} < K$ and $\bar{v} > 0$ when $m < e_1aK(1-c)/(h+K(1-c))$. Therefore, we conclude that if

then system (1.1) has a unique coexistence positive constant steady state $E_2$.

We now establish the local stability and global stability of $E_2$. Let

We first give a relatively strong local stability criterion for $E_2$.

Theorem 3.3. If (3.7) and $\bar{u}\geq \alpha$ hold, then $E_2$ is locally asymptotically stable with respect to (1.1).

Proof. It follows from (3.3) that

where

The corresponding $k$-th characteristic equation for the linearized system of (1.1) at $E_2$ is

where

Note that if $T_k < 0$ and $D_k > 0$ for all $k\in {\mathbb{N}}_0$, then $E_2$ is locally asymptotically stable. From (3.5), we have

if $\bar{u}\geq \alpha$. Combining with (3.8) gives $T_k < 0$ and $D_k > 0$ for all $k\in {\mathbb{N}}_0$.

Let

Theorem 3.4. Assume that (3.7) and $\bar{u} < \alpha$ hold. If

then $E_2$ is locally asymptotically stable with respect to (1.1).

Proof. It follows from (3.9) and (3.10) that

for all $k\in {\mathbb{N}}_0$ if (3.11a) holds. A direct calculation gives $\bar{a}_{11}\bar{a}_{22}-\bar{a}_{12}\bar{a}_{21} > 0$. From the second equality of (3.10), we have the following two cases: (1) $-\bar{a}_{11}/\bar{a}_{22}\leq d_u/d_v$. It is clear that $D_k > 0$ for all $k\in {\mathbb{N}}_0$ and $\mu_k\in \Lambda$ since $d_ud_v > 0$ and $\bar{a}_{11}\bar{a}_{22}-\bar{a}_{12}\bar{a}_{21} > 0$; (2) $-\bar{a}_{11}/\bar{a}_{22} > d_u/d_v$. It is not difficult to show that if

where $A = \sqrt{\bar{a}_{12}\bar{a}_{21}(\bar{a}_{12}\bar{a}_{21}-\bar{a}_{11} \bar{a}_{22})}$, then $(d_u\bar{a}_{22}+d_v\bar{a}_{11})^2-4d_ud_v(\bar{a}_{11}\bar{a}_{22}-\bar{a}_{12}\bar{a}_{21}) < 0$, which implies that $D_k > 0$ for all $k\in {\mathbb{N}}_0$. Hence, if (3.12) holds, then we have $D_k > 0$ for all $k\in {\mathbb{N}}_0$ since

The proof is completed.

Remark 3.1. The local stability of $E_2$ is independent of diffusion coefficient $d_u, d_v$ when $\bar{u}\geq \alpha$ in Theorem 3.3, and related to diffusion coefficient $d_u, d_v$ when $\bar{u} < \alpha$ in Theorem 3.4.

Let

where

By direct calculation, we conclude that $B_2\geq B_1$ when $h\geq K(1-c)$; $B_1 > B_2 > B_3$ when $h < K(1-c)$ and $\eta > A_1$; $B_1 > B_3\geq B_2$ when $h < K(1-c)$ and $A_2 < \eta\leq A_1$.

We next investigate the global stability of $E_2$ by using the upper and lower solutions method.

Theorem 3.5. $E_2$ is globally asymptotically stable with respect to (1.1) if $(m, \eta, h)\in \Delta_1, \Delta_2, \Delta_3, \Delta_4$.

Proof. Note that if $(m, \eta, h)\in \Delta_1, \Delta_2, \Delta_3, \Delta_4$, then (3.7) and (2.2) hold. It follows from Theorem 2.2 that

From Theorem 2.3, we have

For any $0 < \epsilon < \underline{v}_1$ there exists a $T > 0$ such that $v\geq \underline{v}_1- \epsilon$ and $v\leq \bar{v}_1+ \epsilon$ for all $(t, x)\in[T, \infty)\times\bar{ \Omega}$. Then

and

which imply that

Let

A direct calculation gives

Hence,

We construct four sequences $\{\underline{u}_i\}$, $\{\underline{v}_i\}$, $\{\bar{u}_i\}$ and $\{\bar{v}_i\}$ by

It follows from (3.12) and (3.13) that

Then we have

and $0 < \underline{\psi}\leq\bar{\psi}, 0 < \underline{\phi}\leq\bar{\phi}$. By (3.13), we get

and then

We now prove $\underline{\psi} = \bar{\psi}$ and $\underline{\phi} = \bar{\phi}$. Two equations subtraction in (3.16b) gives

which means that if $\bar{\phi} = \underline{\phi}$, then $\bar{\psi} = \underline{\psi}$, and vice versa. Substituting (3.16b) into (3.16a), we obtain

(3.17) minus (3.18) gives

If $\bar{\psi}\neq\underline{\psi}$, then

This shows that if $(m, \eta, h)\in \Delta_1, \Delta_2$, then $\bar{\psi}+\underline{\psi} < 0$, which is a contradiction. (3.17) plus (3.18) gives

This proves that if $(m, \eta, h)\in \Delta_3, \Delta_4$, then $\bar{\psi}\underline{\psi} < 0$, which is a contradiction. The above results suggest that $\bar{\psi} = \underline{\psi}$ and $\bar{\phi} = \underline{\phi}$ if $(m, \eta, h)\in \Delta_1, \Delta_2, \Delta_3, \Delta_4$. Combining with (3.5), we have $\bar{\psi} = \underline{\psi} = \bar{u}$ and $\bar{\phi} = \underline{\phi} = \bar{v}$. From (3.14) and (3.15), we obtain $\lim_{t \rightarrow \infty}(u(x, t), v(x, t)) = (\bar{u}, \bar{v})$ uniformly on $\bar{ \Omega}$. The proof is complete.

3.3.   Turing instability

It has proved that diffusion could destabilize an otherwise stable steady state of the reaction-diffusion system and lead to nonuniform spatial patterns. This kind of instability, essentially originated in landmark work of Turing ^[35], is usually called Turing instability or diffusion-driven instability.

We assume that $\bar{u} < \alpha$, (3.11a) and

hold. Then the quadratic equation $d_ud_v \omega^2-(d_u\bar{a}_{22}+d_v\bar{a}_{11}) \omega+\bar{a}_{11}\bar{a}_{22}-\bar{a}_{12}\bar{a}_{21} = 0$ has two real positive roots

Theorem 3.6. Assume that $\bar{u} < \alpha$, (3.11a) and (3.19) hold. Then we have the following conclusions:

(ⅰ) if $\Lambda\cap(\omega_1(d_u, d_v), \omega_2(d_u, d_v)) = \emptyset$, then $E_2$ is locally asymptotically stable with respect to (1.1);

(ⅱ) if $\Lambda\cap(\omega_1(d_u, d_v), \omega_2(d_u, d_v))\neq\emptyset$, then the positive constant steady state $E_2$ of system (1.1) is Turing unstable;

(ⅲ) for a fixed $d_v > 0$, there exists $d_* > 0$ such that $E_2$ is Turing unstable when $0 < d_u < d_*$;

(ⅳ) there exists $d^* > 0$ such that $E_2$ is locally asymptotically stable when $d_v > d^*$ and $d_u > \bar{a}_{11}/\mu_1$.

Proof. Obviously, (ⅰ) and (ⅱ) hold. Note that

for a fixed $d_v > 0$ and

for a fixed $d_u > 0$. This implies that (ⅲ) and (ⅳ) hold.

3.4.   Simulations

In this subsection, we do some numerical simulations to illustrate our analysis of steady states for system (1.1). This has been showed that at some stage in the life cycle, 90% of some zooplankton's food is obtained by cannibalism ^[12]. This also means that cannibalism is widespread in aquatic systems. Therefore, here the set of parameter values we use is derived from the phytoplankton-zooplankton system. The values of all parameters are listed in Table 2.

In mathematical theory, the total extinction of predator and prey will never occur since $E_0$ is unstable (see Theorem 3.1). However, this can happen in real nature when the predator and prey density become very small. Figure 1 and Figure 2 show solutions of (1.1) converge to constant steady states $E_1$ or $E_2$ for different parameter value $m$ while other parameters are from Table 2. For the case of $m = 0.18$, one can see that the extinction of predator with prey reaching its carrying capacity ($E_1$) is a possible outcome of system (1.1) (see Theorem 3.2 and Figure 1). For $m = 0.06$, predator and prey can coexist together at a positive constant steady state $E_2$ (see Figure 2). In Figure 3, Turing instability may arise from system (1.1) if (ⅱ) or (ⅲ) in Theorem 3.6 holds. Turing instability destroys the spatial symmetry and causes the pattern formation which is stationary in time and oscillatory in space ^[6,39].

4.   Nonconstant positive steady state solutions

As an indication of dynamical complexity of predator-prey systems, it is important to investigate the existence of nonconstant positive steady state solutions, also called stationary patterns, in the spatially inhomogeneous case. In this section, we explore the nonexistence and existence of nonconstant positive steady state solutions of (1.1), which satisfy

4.1.   A priori estimates and nonexistence of nonconstant solutions

To establish the existence and nonexistence of nonconstant positive steady state solutions, we need to derive some a priori estimates for positive solutions of (4.1). We introduce the following maximum principle.

Lemma 4.1. (Maximum principle ^[5,40]) Assume that $f\in C(\Omega)$ and $c_j\in C(\Omega)$ with $j = 1, 2, \cdots, n$.

(ⅰ) If $\omega\in C^2(\Omega)\cap C^1(\bar \Omega)$ satisfies

and $\omega(x_0) = \max \limits_{x\in\bar \Omega} \omega(x)$, then $f(x_0)\geq0;$

(ⅱ) If $\omega\in C^2(\Omega)\cap C^1(\bar \Omega)$ satisfies

and $\omega(x_0) = \min \limits_{x\in\bar \Omega} \omega(x)$, then $f(x_0)\leq0.$

We first have a priori upper bound estimates for any positive solution of (4.1).

Lemma 4.2. Assume that $(u(x), v(x))$ is a positive solution of (4.1). If (3.7) holds, then

Proof. Let $u(x_1) = \max\limits_{\bar{ \Omega}}u(x)$, $v(x_2) = \max\limits_{\bar{ \Omega}}v(x)$. From (4.1), we have

By Lemma 4.1, we obtain $ru(x_1)(1-u(x_1)/K)\geq0$, which means that $\max\limits_{\bar{ \Omega}}u(x)\leq K$. It follows from (4.1) that

Lemma 4.1 shows that

Then $\max\limits_{\bar{ \Omega}}v(x)\leq (K(e_1a-m)(1-c)-mh)/(\eta(a(1-e_2)+m))$.

We now establish the nonexistence of nonconstant positive solutions of (4.1) if the diffusion coefficients $d_u$ and $d_v$ are large.

Theorem 4.1. If (3.7) holds, then there is a $\hat{d} > 0$ such that system (4.1) has no nonconstant positive solution when $d_u, d_v\geq \hat{d}$.

Proof. Let $(u, v)$ be a positive solution of system (4.1), and denote $\tilde{u} = |\Omega|^{-1}\int_{\Omega}udx, \tilde{v} = |\Omega|^{-1}\int_{\Omega}vdx.$ Then $\int_{ \Omega}(u- \tilde{u})dx = \int_{ \Omega}(v- \tilde{v})dx = 0.$ Multiplying the first equation of system (4.1) by $u- \tilde{u}$, and integrating over $\Omega$, we obtain

From Lemma 4.2, $\max\limits_{ \tilde{ \Omega}}v(x)\leq (K(e_1a-m)(1-c)-mh)/(\eta(a(1-e_2)+m)): = \delta$. Multiplying the second equation of system (4.1) by $v- \tilde{v}$, and integrating over $\Omega$, we have

Let

Hence, by the Poincaré inequality, we get

This means that if $\min\{d_u, d_v\} > \max\{C_1/\mu_1, C_2/\mu_2\}$, then $\nabla(u- \tilde{u}) = \nabla(v- \tilde{v}) = 0$ and $u\equiv \tilde{u}, v\equiv \tilde{v}$.

4.2.   Existence of nonconstant positive steady state solutions

In this part, we explore the existence of nonconstant positive solutions to (4.1) by using degree theory. To do this, we recall the following Harnack inequality.

Lemma 4.3. (Harnack inequality ^[5,41]) If $u\in C^2(\Omega)\cap C^1(\bar \Omega)$ is a positive solution of

where $b\in C(\Omega)\cap L^{\infty}(\Omega)$, then there exists a positive constant $L$ which depends only on $M$, satisfying $\|b\|_{\infty}\leq M$, such that

We now establish a prior lower bound estimates for positive solutions of system (4.1).

Lemma 4. If $(u(x), v(x))$ is a positive solution of (4.1) and (3.7) holds, then there exists a constant $\underline{C} > 0$ depending possibly on $d_u, d_v, \Omega, n$ and parameters of (4.1), such that

Proof. Let $u(x_3) = \min\limits_{\bar{ \Omega}}u(x)$. From (4.1) and Lemma 1, we have

Hence,

Let

There is a positive constant $M$ depending on $d_u, d_v, \Omega, n$ and parameters of (4.1) such that $\|b_1\|_{\infty}\leq M$, $\|b_2\|_{\infty}\leq M$ since (4.2) holds. By using Harnack inequality in Lemma 4.3, there exists a positive constant $L$ which depends only on $M$ such that

It only need to prove that there exists a $\bar{L} > 0$ such that

If it is not true, then there exists a sequence of positive solutions $\{(u_n(x), v_n(x))\}_{n = 1}^ \infty$ such that

From the standard regularity theorem for the elliptic equations, there exists a subsequence of $\{(u_n, v_n)\}_{n = 1}^ \infty$, which we still denote by $\{(u_n, v_n)\}_{n = 1}^ \infty$, and two nonnegative functions $\hat{u}, \hat{v}\in C^2(\bar{ \Omega})$ such that $u_n \rightarrow \hat{u}$ and $v_n \rightarrow \hat{v}$ in $C^2(\bar{ \Omega})$ as $n \rightarrow \infty$. By (4.2), (4.5) and (4.4), we have $0 < \hat{u}\leq K$ and either $\hat{u}\equiv 0, \hat{v}\not\equiv 0$ or $\hat{u}\not\equiv 0, \hat{v}\equiv0$. Note that $(u_n, v_n)$ is a positive solution of (4.1), then

If $\hat{u}\equiv 0, \hat{v}\not\equiv 0$, then

for sufficiently large $n$ since $u_n \rightarrow 0$ as $n \rightarrow \infty$. It is a contradiction to (4.6b) since $v_n > 0$. If $\hat{u}\not\equiv 0, \hat{v}\equiv0$, then from (4.6a), we obtain $\int_ \Omega \hat{u}(1-\hat{u}/K)dx = 0$. It follows from $0 < \hat{u}\leq K$ that $\hat{u}\equiv K$. Thus, we have

as $n \rightarrow \infty$ since (3.7) holds. This contradicts (4.6b).

Summarizing the discussion above, we conclude that (4.5) holds, which implies that (4.3) holds. This completes the proof.

We now investigate the existence of nonconstant positive solutions of system (4.1) by using the Leray-Schauder degree theory (^[42]) and the methods in ^[5,43]. Denote

where $\bar{C} = \max\{K, (K(e_1a-m)(1-c)-mh)/(\eta(a(1-e_2)+m))\}$ and $X$ can be found in (3.2). Note that if (3.7) holds, then (4.1) has a unique positive constant solution $E_2 = (\bar{u}, \bar{v})$. Let

with $U = (u, v)^T\in X$ and $(I-\triangle)^{-1}$ be the inverse of $I-\triangle$. Then system (4.1) can be rewritten as

where $I-\triangle$ satisfies the homogeneous Neumann boundary condition. Frechét derivative of system (4.8) with respect to $(u, v)$ at $(\bar{u}, \bar{v})$ is

It is clear that $\zeta$ is an eigenvalue of $F_{U}(d_1, d_2, \bar{u}, \bar{v})$ on $X_i$ with $i\in {\mathbb{N}}_0$ if and only if $\zeta(1+\mu_i)$ is an eigenvalue of the matrix

Then

where

From Subsection 3.3, if $\bar{u} < \alpha$ and (3.19) hold, then $S(d_u, d_v, \mu) = 0$ has two positive roots

and

for a fixed $d_u > 0$. Let

Lemma 4.5. (^[5]) If $S(d_u, d_v, \mu_i)\neq 0$ for all $\mu_i\in\Lambda$, then ${\rm{index}}(F(d_u, d_v, \cdot), (\bar{u}, \bar{v})) = (-1)^{\sigma},$ where $\sigma = \sum\limits_{\mu_i\in \mathcal W(d_u, d_v) \cap \Lambda}m(\mu_i)$ when $\mathcal W(d_u, d_v)\cap \Lambda\neq\phi$ and $\sigma = 0$ when $\mathcal W(d_u, d_v)\cap \Lambda = \phi$. In particular, if $S(d_u, d_v, \mu) > 0$ for all $\mu\geq0$, then $\sigma = 0$.

Theorem 4.2. Assume that $d_u > 0$, $\bar{u} < \alpha$, and (3.19) hold. If $\bar{a}_{11}/d_u\in (\mu_q, \mu_{q+1})$ for some $q\in {\mathbb{N}}$ and $\sum_{i = 1}^q m(\mu_i)$ is odd, then there is a positive constant $\tilde{d}_v$ such that for any $d_v > \tilde{d}_v$, (4.1) has at least one nonconstant positive solution.

Proof. It follows from (4.9) and $\bar{a}_{11}/d_u\in (\mu_q, \mu_{q+1})$ that there exists a sufficient large $d_0$ such that for any $d_v > d_0$

From Theorem 4.1, system (4.1) has no nonconstant positive solution for any $d_u, d_v > \hat{d}$. We choose $\tilde{d}_u > \hat{d}$ such that $\bar{a}_{11}/ \tilde{d}_u < \mu_1$ and $\tilde{d}_v > \max\{\hat d, d_0\}$ such that

Assume that the conclusion of Theorem 4.2 is not true. Then there is some $d_v$ such that system (4.1) has no nonconstant positive solution for $d_v\geq \tilde{d}_v$. For $\kappa\in[0, 1]$, we let $D_\kappa = \;\mbox{diag}\;(\kappa d_u+(1-\kappa) \tilde d_u, \kappa d_v+(1-\kappa) \tilde d_v)$ and consider the following system

where $G(U)$ is defined in (4.7). Obviously, (4.12) is equivalent to

Note that $\Phi(U, 1) = F(d_u, d_v, U)$, $\Phi(U, 0) = F(\tilde d_u, \tilde d_v, U)$ and

The above arguments show that $\Phi(U, 1) = 0$ and $\Phi(U, 0) = 0$ have no nonconstant positive solution.

From (4.10) and (4.11), we have

which imply that

By using Lemmas 4.2 and 4.4, we obtain $(u, v)\in\Theta$ for any solution $(u, v)$ of system (4.1) on $\bar{ \Omega}$. Then $\Phi(U, \kappa)\neq0$ on $\partial \Theta$ for all $\kappa\in[0, 1]$. It follows from the homotony invariance of Leray-Schauder degree that

Note that $\Phi(U, 0) = 0$ and $\Phi(U, 1) = 0$ have only the constant solution $(\bar{u}, \bar{v})$ in $\Theta$ and hence,

which is a contradiction to (4.13). The proof is complete.

5.   Influence of prey refuge and predator cannibalism on biomass

The predator and prey biomass density in an ecosystem is an important index for avoiding population extinction and protecting biological diversity. In this section, we will investigate the influence of prey refuge and predator cannibalism in (1.1) on predator and prey biomass density. To facilitate the discussion below, we let $\Omega = [0, 40]$ and use the spatial average of $u(x, t)$ and $v(x, t)$ defined as

We consider the effect of predator cannibalism rate $\eta$. In Figure 4, we compare the (spatial averaged) coexistence constant or nonconstant steady states $(U, V)$ for different values $\eta$. From Figure 4 left panel, one can observe that prey biomass density is increasing gradually with the increase of $\eta$. This shows that predator cannibalism is beneficial to prey biomass density. From Figure 4 right panel, there exists a $\eta^*$ such that predator biomass density is increasing gradually when $0.125 < \eta < \eta^*$, and decreasing gradually when $\eta > \eta^*$. This confirms that appropriate predator cannibalism ($\eta = \eta*$) has a positive effect on predator biomass density, and then high predator cannibalism has a negative effect on predator biomass density.

When predator cannibalism is low ($\eta = 0.08$), Figure 5 shows that system (1.1) can produce Hopf bifurcation which destroys the temporal symmetry and induces periodic oscillations that are uniform in space and periodic in time for carrying capacity $K = 10$ of prey. But when predator cannibalism is high, predator and prey biomass density converge to a positive constant steady state (see Figure 2). These indicate that high predator cannibalism may stabilize the predator-prey system, and prevent the paradox of enrichment.

Prey refuge is an effective strategy for protecting prey population and avoiding over-predation. Figure 6 shows that prey refuge has a beneficial influence on prey biomass density, and a negative influence on predator biomass density. From the perspective of biodiversity conservation, prey refuge in point has a better effect for maintaining the persistence of predator-prey system (see Theorems 2.2 and 2.3). Excessive or low prey refuge is likely to destroy the balance of ecosystems.

6.   Discussion

In this paper, we analyze a diffusive predator-prey system (1.1) with prey refuge and predator cannibalism. We now roughly summarize our main results as below: (1) system (1.1) is dissipation and persistence (see Theorems 2.2 and 2.3); (2) the existence, local and global stability of constant steady states are established (see Theorems 3.2, 3.3, 3.4 and 3.5); Turing instability caused by diffusion is given (see Theorem 3.6); (3) the nonexistence and existence of nonconstant steady state solutions is investigated (see Theorems 4.1 and 4.2); (4) Studies show that appropriate predator cannibalism has a positive effect on predator-prey ecosystem (see Figure 4).

We do some theoretical analysis to explore threshold conditions for the regime shift from extinction to coexistence of predator and prey. Our results show that the total extinction of predator and prey will never occur, but this can happen ecologically even though the equilibrium at the origin $E_0$ is unstable. This is because that organisms are discrete and can be completely eliminated when the densities become very small. If $m > e_1aK(1-c)/(h+K(1-c))$, then predator is extinct and prey reaches its maximum environmental capacity. The above condition also shows that the possibility of predator extinction increases with the gradual increase of prey refuge ratio. This means that excessive prey refuge has a negative effect on predator-prey system, and is also not conducive to biodiversity conservation. Predator and prey can coexist together in three different forms: constant steady state: nonconstant steady state; periodic oscillations in time or space.

In previous studies, it has been widely believed that predator cannibalism has a negative effect on predator biomass density. However, our studies point out that appropriate predator cannibalism can not only increase prey biomass density, but also enhance predator biomass density under the right circumstance. From the ecological point of view, the reason why this can happen is that appropriate predator cannibalism can moderately reduce predator pressure of prey and enhance prey biomass density that leads to an increase in predator biomass density. On the other hand, it is worth noting that high predator cannibalism may stabilize the predator-prey system, and prevent the paradox of enrichment. This is because that high predator cannibalism increases intraspecific competition among predators, and then reduces the possibility of population oscillation. Results above indicate appropriate predator cannibalism has a positive effect on predator-prey ecosystem.

Spatial environmental parameters $d_u, d_v$ have an important influence on dynamical properties of system (1.1). If diffusion coefficients $d_u, d_v$ are sufficiently large, then predator and prey are evenly distributed in space. By contrast, when $d_u$ is very low for a fixed $d_v$, $E_2$ loses its stability and Turing instability occurs. This produces a steady state solution of spatial inhomogeneity called the pattern formation. This implies that spatial distribution patterns and dispersal mechanisms can make predator-prey system exhibit more complex dynamical properties. Our numerical simulations also show that diffusion coefficients $d_u, d_v$ have no significant effect on predator and prey oscillation. By comparing Figure 7 and Figure 5, when $d_u$ takes three different values: 0.001, 0.1 and 100 for a fixed $d_v = 0.1$, there is no obvious change in the period and amplitude with time for predator and prey biomass density. This indicates that diffusion coefficients do not have a fundamental impact on the paradox of enrichment in system (1.1).

This paper attempts to investigate dynamics of system (1.1) and the influence of prey refuge and predator cannibalism. It is important to understand the existence and stability of Hopf bifurcation when predator cannibalism rate $\eta$ changes, which are not discussed in this paper. In view of the important role of Allee effects or time delay in the predator-prey system, it will be of interest to further model dynamic properties of system (1.1) with Allee effects or time delay.

Acknowledgments

The research was partially supported by Heilongjiang Postdoctoral Funds for scientific research initiation (Q17148), NSFHLJ-JJ2019LH0178 and the State Scholarship Fund ensured by CSC (201806620040).

Conflict of interest

All authors declare no conflicts of interest in this paper.