Global analysis and Hopf-bifurcation in a cross-diffusion prey-predator system with fear effect and predator cannibalism

1.   Introduction

In recent years, many scholars have established the prey-predator models to investigate the dynamic relations, multi-species systems^[1,2,3,4] have been studied. All these works are conducive to the analysis of the system, and to achieve the purpose of stabilizing the ecosystem through the modeling research between ecosystems. The problem of ecological species still deserves researchers' attentions. On the other hand, different functional response items play an important role in portraying different situations. For example, Holling type-Ⅱ function response is used in different systems ^[5,6] for scholars to investigate different systems. In nature, in order to escape from feeding behavior to avoid species extinction and be better to fit the actual environment, we always consider an effective strategy which is to establish a protection zone in the system. Many scholars are attracted by this interesting topic, and some researchers have established and investigated the predator-prey systems incorporating a prey refuge ^[7,8]. Thus, adding the prey refuge is significant. We consider a term $ \frac{1}{1+a(1-m)u} $ which can reflect the prey refuge in this article. Applying the idea of prey refuge is a very convincing hypothesis for the prey-predator system and it also adds a new feature in the predator-prey models with prey dependence.

In a recent paper ^[9], the author investigated a deterministic prey-predator system only with refuge, and proposed the following system,

However, we should notice that after the establishment of prey refuges, there may be competition among prey due to food shortage and other reasons, so we need to consider the factor of intraspecies competition. On the other hand, when the predator kills the prey, the anti-predator defence behavior may be altered on account of the fear of predators. There is now a popular saying: the prey will change its own behavior and physiological characteristics due to the fear of the predator as soon as the predator appears in front of the prey. In addition, this fear effect affects the normal birth rate of the prey. To a certain extent, it exceeds the impact of direct predation and even in some cases it becomes more influential than the direct predation. Hence it is important to investigate the fear effect. Many scholars have researched the fear effect in different stochastic models or ordinary differential equation(ODE) systems ^[10,11,12]. Authors in ^[13] believed that increasing the cost of fear effect or decreasing the strength of refuge can lead to the increase of predators. We multiply the factor $ \frac{1}{1+kv} $, which aims to better study the fear effect on the population, where $ k $ is the fear factor. In ^[14], the authors proposed a model with prey refuge and fear effect

Moreover, another factor that can not be ignored is predator cannibalism. The predator cannibalism has a strong impact on the dynamics systems^[15,16], it may affect the original stability of the system. Yet these behaviors are always omitted in the study of prey-predator systems. In ^[17], N.A. Schellhorn and D.A. Andow found a phenomenon that the effect of cannibalism was equal to or stronger than that of interspecific predation for both species. Authors in ^[18] explored a deterministic system as follows

Hence, to account for the above factors, we modify the system (1.3) and propose the following ordinary differential equation (ODE) system:

where $ r $ denotes the birth rate, $ \beta $ describes the capture rate. $ b $ represents the ability of competition between the prey population. $ c $ is the efficiency of food conversion. $ d $ is the loss rate of predator. $ 1+a(1-m)u $ represents the amount of preys which is available to the predators, and $ m \in[0, 1) $. $ \eta $ is the cannibalism rate between predators. $ k $ is the fear factor. All of the coefficients in this article are positive.

From the system (1.4), we know that:

(1) when $ k = 0 $, $ \eta = 0 $, then system (1.4) is similar to the system (1.1),

(2) when $ \eta = 0 $, which means there is no predator cannibalism in our system, then the system (1.4) becomes system (1.2),

(3) when $ k = 0 $, which means we do not consider the fear effect in our system, and the function response terms from $ \frac{\beta(1-m)uv}{1+a(1-m)u} $ to $ \frac{auv}{h+u+\eta v} $, then the system (1.4) becomes the system (1.3).

To sum up, we've improved on previous models. With the development of mathematical theory researches ^{[19,20,21,22,23,24]}. In order to describe the relationship of attraction and repulsion among species, a number of predator-prey systems with cross-diffusion was established ^[25,26]. Many articles have shown that these diffusion coefficients may lead to a Turing instability ^[27,28,29]. Hence, we consider the effects of diffusion coefficients and give the following reaction-diffusion system,

$ \Omega $ is a connected open region with a smooth boundary $ \partial\Omega $. $ \nu $ is the outward unit normal vector. $ d_{u} $ and $ d_{v} $ are self-diffusion coefficients. $ d_{12} $ and $ d_{21} $ are cross-diffusion coefficients. The other parameters have the same meaning as above. In ^[30,31,32], researchers investigated Hopf bifurcations of their systems to explain the complex spatiotemporal dynamics. Thus, we need to analyse the dynamics of system (1.5). We establish sufficient conditions where Hopf bifurcation occurs, as well as the stability. Authors in ^[33] have well combined theoretical results with experiments, we will also give some corresponding numerical simulation results.

We mainly analyze from the following aspects:

(i) the existence and persistence of the nonnegative solutions,

(ii) the local stability of system (1.5) at steady states,

(iii) Hopf bifurcation induced by the predator cannibalism in the deterministic system (1.4) and the reaction-diffusion system (1.5),

(iv) Turing instability of system (1.5).

The paper is structured as follows. The dynamics of the system, including existence of the solutions and the persistence of the system (1.5) are studied, and the proofs are given in Section 2. Sufficient conditions for the stability and Hopf bifurcation of the deterministic system are given in Section 3. Then, in Section 4, we give some conclusions for the reaction diffusion system: stability analysis, Turing instability. In Section 5, we discuss some well-posedness of the nonconstant steady state. Finally, numerical simulations and conclusion are given.

2.   Persistence of the system (1.5)

Utilizing the maximum principle and comparison principle for parabolic equations, we have:

Theorem 2.1. For system (1.5), it admits a unique global solution $ (u(x, t), v(x, t)) $ for all $ x\in \Omega, t > 0. $

Proof. Let $ \left(\bar{u}(x, t), \bar{v}(x, t)\right) $ = $ \left(\bar{u}(t), \bar{v}(t)\right) $, where $ \left(\bar{u}(t), \bar{v}(t)\right) $ satisfies

It is not difficult to find that $ (\bar{u}(t), \bar{v}(t)) $ is global existence according to the ordinary differential equations' existence and uniqueness theorem. Then we have

There exists $ T_{0} > 0 $ such that $ u(x, t)\leq \bar{u}(t) $ according to comparison theorem, similarly, for $ v(x, t) $, we have

By using the comparison theorem, there exists $ T_{0}\in(0, \infty) $ such that $ u(x, t)\leq\bar{u}(t) $, similarly, $ v(x, t) $gets the same conclusion. The solution $ (u(x, t), v(x, t)) $ satisfies,

for all $ x\in\Omega $, $ t\geq0 $.

Theorem 2.2. The non-negative solution $ (u(x, t), v(x, t)) $ in system (1.5) is bounded^[34],

if $ 0 < m < 1-\frac{bd}{(c\beta-ad)r} $ and $ c\beta-ad > 0 $ hold.

Proof. According to the comparison theorem of parabolic equation problems, it is clear that the first equation in system (1.5) satisfies $ \frac{\partial u}{\partial t}-d_{u}\Delta u-d_{12}\Delta uv \leq ru-bu^{2} $. Thus for a $ T > 0 $, $ t > T, $ $ u(x, t)\leq\frac{b}{r}+\epsilon $, where $ \epsilon $ is an arbitrary positive constant, then,

Define $ z_{1}(t) $ be a solution with respect to the equation

If

If

For any arbitrariness of $ \epsilon > 0 $, one has

This finishes the proof.

Proposition 2.1. Theorem 2.2 shows that, $ \Gamma: \left[0, \frac{r}{b}\right]\times\left[0, \frac{(c\beta-ad)(1-m)r-bd}{\eta b}\right] $ is the global attractor for the solutions of system (1.5), and any non-negative solution $ (u(x, t), v(x, t)) $ of system (1.5) is attracted to $ \Gamma $ for a sufficiently large $ t $.

Assumption 2.1. (i) $ 0 < m < 1-\frac{bd-1}{(c\beta-ad)r}, $

(ii) $ \eta br > \beta(1-m)\left[(c\beta-ad)(1-m)r-bd\right]\left[1+k(c\beta-ad)(1-m)r-bd\right], $

(iii) $ d < (c\beta-ad)(1-m)\left(\frac{r}{b+kb\bar{v}}-\frac{\beta(1-m)}{b}\bar{v}\right). $

Theorem 2.3. Under Assumption 2.1, the system (1.5) is persistent, and one has,

where, $ \underline{A} = \frac{r\eta}{\eta b+k(c\beta-ad)(1-m)r-bdk}-\frac{r\beta(1-m)^{2}(c\beta-ad)-\beta bd(1-m)}{\eta b^{2}} > 0. $

Proof. For the first equation in system (1.5), $ u(x, t) $ satisfies

One has $ u(x, t)\geq z_{2}(x, t) = \underline{A}-\epsilon, $ where $ u(x, t) $ is the upper solution.

An application of comparison of parabolic equation, the first equation of (1.5) holds. Obviously, $ u(x, t)\geq \underline{A} $ for all $ x\in \Omega, $ when there exists a $ T _{1} > 0 $ and $ t\geq T _{1}. $ Then we obtain

Note, $ z_{3}(t) $ satisfies

then $ {\lim\limits_{t\rightarrow +\infty}z_{3}(t) = -\frac{d}{\eta}+\frac{c\beta(1-m)\underline{A}}{(1+a(1-m)\underline{A})\eta}} > 0. $ Which means

Combining Theorems 2.2 and 2.3, we obtain that system (1.5) is persistent.

Then the proof is completed.

3.   The system (1.5) without diffusion: stability and Hopf bifurcation analysis

We discuss related properties of various constant equilibrium points. There are three types constant equilibrium points: extinction point $ E_{0}(0, 0), $ semi-trivia equilibrium point $ E_{1}(\frac{r}{b}, 0), $ and the positive equilibrium point $ E_{2}(u_{*}, v_{*}). $

Definition 3.1. $ U $ is an open region in $ \texttt{R}^{n} $, for any given $ \varepsilon > 0 $, and the equilibrium $ x_{*}\in U $, there exists a $ \delta > 0 $ and $ t_{0} > 0 $ such that $ \|x(t)-x_{*}\| < \varepsilon $ for $ t\geq t_{0} $. We can say the equilibrium point $ x_{*} $ is stable. We can find a $ T > 0 $ such that when $ t > T $, $ \lim_{t\rightarrow \infty}\|x(t)-x_{*}\| = 0 $, then the equilibrium point $ x_{*} $ is attractive. The equilibrium point $ x_{*} $ is said to be asymptotically stable if it is stable and attractive.

3.1.   Stability of $ E_{0} $ and $ E_{1} $

We first consider $ E_{0}(0, 0) $ and $ E_{1}(\frac{r}{b}, 0). $

(i) $ E_{0}(0, 0) $ exists and is always unstable;

(ii) $ E_{1}(\frac{r}{b}, 0) $ is always exists, and $ E_{1} $ is locally asymptotically stable when $ c\beta-ad\leq 0 $ or $ c\beta-ad > 0 $ and $ 1-\frac{bd}{(c\beta-ad)r} < m < 1 $. When $ c\beta-ad > 0 $ and $ 0\leq m < 1-\frac{bd}{(c\beta-ad)r}, $ $ E_{1} $ is unstable.

From (2.5) we know, when $ 1-\frac{bd}{(c\beta-ad)r} < m < 1 $,

Notice that there exists a $ T_{2} > 0, $ for any $ \epsilon > 0 $, such that $ v(x, t)\leq\epsilon $ for all $ x\in\Omega $ and $ t\geq T_{2} $. Then we have,

Let $ z_{4}^{'}(t) = z_{4}\left[\frac{r}{(1+k\epsilon)}-bz_{4}-\frac{\beta(1-m)\epsilon}{1+a(1-m)z_{4}}\right], $ then $ {\lim\limits_{t\rightarrow +\infty} z_{4}(t) = \frac{r}{b}} $ since $ \epsilon $ is sufficiently small.

Applying the comparison theorem, we get the conclusion that $ {\liminf\limits_{t\rightarrow +\infty}\min u(x, t)\geq\frac{r}{b}}. $ Combining with (2.3), $ E_{1} $ is global asymptotically stable. Now, we discuss $ E_{2}({u_{*}, v_{*}}) $ of the system (1.5). Firstly, in order to study the following parts we make some notations,

where

and

3.2.   Analysis of $ E_{2}(u_{*}, v_{*}) $

$ E_{2} $ satisfies the follow algebraic equations

According to equation (3.3), one has

Putting (3.4) into (3.2), yields

where

$ A = -abk\eta(c\beta-ad)(1-m)^{2} < 0, $

$ B = \left[abkd\eta-ab\eta^{2}-bk\eta^{2}-bk\eta (c\beta-ad)\right](1-m)-\beta k(c\beta-ad))^{2}(1-m)^{3}, $

$ C = r\eta^{2}a(1-m)-b\eta^{2}+bdk\eta-\beta\eta(1-m)^{2}(c\beta-ad)+2\beta kd(1-m)^{2}(c\beta-ad), $

$ D = r\eta^{2}+\beta d(\eta-kd)(1-m). $

From the Descarte's rule of sign, $ A < 0 $ always exists and the Equation (3.5) has one unique positive solution $ u_{*}, $ if any case holds

The Equation (3.5) has two positive solutions or no positive solution if

The Equation (3.5) has three positive solutions or one positive solution if

Theorem 3.1. (Ⅰ): Any condition in (3.6) holds, there exists one positive constant steady state for system (1.5).

(Ⅱ): Any condition in (3.7) holds, there exists two positive constant steady states or no positive steady state for system (1.5).

(Ⅲ): If the condition (3.8) holds, there exists three positive constant steady states or one positive steady state for system (1.5).

3.3.   The stability at $ E_{2}(u_{*}, v_{*}) $

We only consider the condition that $ B < 0, C > 0, D > 0. $ First, we investigate the stability of ODE system.

The Jacobian matrix at $ E_{2}(u_{*}, v_{*}) $ can be expressed as

where

For convenience, we have utilized the following notations

The characteristic polynomial can be expressed as

Define

$ b_{1} = \frac{a_{3}}{a_{1}}-\frac{a_{2}^{2}}{3a_{1}^{2}}, $ $ b_{2} = \frac{a_{4}}{a_{1}}+\frac{2a_{2}^{3}}{27a_{1}^{3}}-\frac{a_{2}a_{3}}{3a_{1}^{2}}, $ $ a_{1} = -a^{2}b\eta(1-m)^{2} $,

$ a_{2} = a(1-m)\left[(c\beta-ad)(\beta+m-1)-2b\eta\right] $, $ a_{3} = a\beta -(b(1-m)^{2}-\eta(1-m)(c\beta-2ad)-b\eta) $, $ a_{4} = d\eta. $

Lemma 3.1. $ E_{2} $ is stable for the ODE system (1.4) if $ \tilde{u}^{3}+b_{1}\tilde{u}+b_{2} = 0 $ and $ \texttt{Det}(J) > 0 $ hold, where

Proof. $ \widetilde{a}_{11}+\widetilde{a}_{22} = \widetilde{a}_{1}u_{*}^{3}+\widetilde{a}_{2}u_{*}^{2}+\widetilde{a}_{3}u_{*}+\widetilde{a}_{4} $, using Cardanor's formula, we have $ \widetilde{a}_{11}+\widetilde{a}_{22} = u_{*}^{3}+b_{1}u_{*}+b_{2}. $

If $ \tilde{u}^{3}+b_{1}\tilde{u}+b_{2} = 0, $ then $ \widetilde{a}_{11}+\widetilde{a}_{22} < 0 $ when $ \tilde{u} < u_{*} < \frac{r}{b}. $

Moreover, when $ \texttt{Det}(J) = \widetilde{a}_{11}\widetilde{a}_{22}-\widetilde{a}_{12}\widetilde{a}_{21} > 0, $ the characteristic polynomial has two strictly negative real parts.

For the sake of convenience, we introduce the following notations. $ q_{1} = -a^{2}b\eta(1-m)^{2}, $ $ q_{2} = a\beta(1-m)^{3}(c\beta-ad)-2ab\eta(1-m), $ $ q_{3} = -b\eta-ad\beta(1-m)^{2}. $ For the proof, we give the following assumptions.

Assumption 3.1. $ q_{2}^{2}-4q_{1}q_{3}\leq0 $ or $ 2b\eta\geq \beta(1-m)^{2}(c\beta-ad). $

Assumption 3.2. Under the condition $ q_{2}^{2}-4q_{1}q_{3} > 0, $ $ 2b\eta < \beta(1-m)^{2}(c\beta-ad), $ and $ \underline{A} < u_{*} < x_{1}, $ $ x_{2} < u_{*} < \overline{u}, $

where $ x_{1} = \frac{-q_{2}-\sqrt{q_{2}^{2}-4q_{1}q_{3}}}{2q_{1}} $, $ x_{2} = \frac{-q_{2}+\sqrt{q_{2}^{2}-4q_{1}q_{3}}}{2q_{1}}, $ and $ \underline{A} $ has the same meaning as before.

3.4.   Hopf bifurcation of ordinary differential system at $ E_{2}(u_{*}, v_{*}) $

We discuss the Hopf bifurcation points at the $ E_{2}(u_{*}, v_{*}), $ we make a translation which is to translate $ E_{2} $ to the origin by translating $ \hat{u} = u-u_{*} $, $ \hat{v} = v-v_{*}. $ Let $ \hat{u} $ and $ \hat{v} $ by $ u $, $ v $, respectively.

Applying Taylar series expansion theorem about the $ (u_{*}, v_{*}) = (0, 0) $, and the system (3.9) is rewritten as

$ J $ has the same meaning as $ J_{\omega} $ and

where $ a_{1} = -b+\frac{a\beta(1-m)^{2}v_{*}}{[a(1-m)u_{*}+1]^{3}}, $ $ a_{2} = -\frac{kr}{(1+kv_{*})^{2}}-\frac{\beta(1-m)}{[a(1-m)u_{*}+1]^{2}}, $

$ a_{3} = \frac{k^{2}ru_{*}}{(1+kv_{*})^{3}}, $ $ a_{4} = \frac{a\beta(1-m)^{2}}{[a(1-m)u_{*}+1]^{3}}, $ $ a_{5} = \frac{k^{2}r}{(1+kv_{*})^{3}}, $ $ a_{6} = -\frac{a\beta(1-m)^{3}v_{*}}{[a(1-m)u_{*}+1]^{4}}, $

$ a_{7} = -\frac{k^{3}ru_{*}}{(1+kv_{*})^{4}}, $ $ b_{1} = -\frac{ac\beta(1-m)^{2}v_{*}+a^{2}\eta(1-m)^{2}v_{*}^{2}}{[a(1-m)u_{*}+1]^{3}}, $ $ b_{2} = \frac{c\beta(1-m)+2a\eta(1-m)v_{*}}{[a(1-m)u_{*}+1]^{2}}, $

$ b_{3} = -\frac{\eta}{a(1-m)u_{*}+1}, $ $ b_{4} = -\frac{ac\beta(1-m)^{2}+2a^{2}\eta(1-m)^{2}v_{*}}{[a(1-m)u_{*}+1]^{3}}, $ $ b_{5} = \frac{a\eta(1-m)}{[a(1-m)u_{*}+1]^{2}}, $

$ b_{6} = \frac{a^{3}c\beta(1-m)^{3}v_{*}+a^{2}\eta(1-m)^{3}v_{*}^{2}}{[a(1-m)u_{*}+1]^{4}}, $ $ b_{7} = 0. $

The characteristic equation of $ J $ is

the roots of (3.11) are $ \kappa_{1} = r(p)+iw(p), $ $ \kappa_{2} = r(p)-iw(p), $ where $ r(p) = \frac{p_{0}-p}{2} = \frac{a_{11}+a_{22}}{2} $, $ w(p) = \sqrt{\texttt{Det}(J)-\left(\frac{1}{2}\texttt{Tr}(J)\right)^{2}}. $ When $ p = p_{0}, $ both of the two roots are a pair of imaginary. As we all known, when $ p > p_{0}, $ $ E_{2}(u_{*}, v_{*}) $ is asymptotically stable. When $ p < p_{0}, $ $ E_{2}(u_{*}, v_{*}) $ is unstable. Therefore, $ p = p_{0} $ is a bifurcation point. When $ \mid p-p_{0}\mid $ is small enough, and $ \kappa = i\sqrt{\texttt{Det}(J)} $. Putting $ \kappa = x+iy $ into (3.11), where $ x $ expresses the real part, and $ y $ expresses the imaginary part. We have

and when $ p = p_{0} $, we can obtain

Next, we set a matrix

clearly,

where $ \tilde{N} = -\frac{a_{21}}{w(p)} $, $ \tilde{M} = \frac{a_{22}-a_{11}}{2w(p)}. $

By the transformation,

then, the Equation (3.10) yields,

where

and

Using the polar coordinate form, (3.14) becomes

Then, when $ p = p_{0} $, according to the Taylor expansion that

We analyse the sign of the coefficient $ a(p_{0}) $,

with

Therefore, we get

From Poincare-Andronov Hopf bifurcation theorem^[38], $ r^{'}(p_{0})\mid_{p = p_{0}} = -\frac{1}{2} < 0, $ and we have the following conclusions:

(i): If $ a(p_{0}) < 0, $ the bifurcating periodic solution is stable and the direction of the Hopf bifurcation is subcritical;

(ii): If $ a(p_{0}) > 0, $ the bifurcating periodic solution is unstable and the direction of the Hopf bifurcation is supercritical.

4.   The system (1.5) with diffusion: stability analysis

Linearizing system (1.5) at $ (u_{*}, v_{*}) $ which yields

where

Define $ A_{i} = -\mu_{i}D_{w}+J. $ and the characteristic equation of $ A_{i} $ at $ E_{2} $ is

and $ \lambda $ is the eigenvalue of the matrix $ A_{i}. $

4.1.   The system (1.5) without cross-diffusion: stability analysis

Theorem 4.1. For $ E_{2}(u_{*}, v_{*}), $ if Assumptions 3.1 and 3.2 hold, the system (1.5) is locally asymptotically stable at $ E_{2} $ without cross-diffusion.

Proof. When we do not consider the cross-diffusion, the $ \texttt{Tr}(A_{i}) $ and $ \texttt{Det}(A_{i}) $ in the Equation (4.1) are

and

respectively. When all eigenvalues have negative real parts, we have the conclusion that $ E_{2}(u_{*}, v_{*}) $ is locally asymptotically stable. In other words, when $ \texttt{Tr}(A_{i}) < 0 $, and $ \texttt{Det}(A_{i}) > 0 $, $ E_{2}(u_{*}, v_{*}) $ is locally asymptotically stable. Then we analyze the numerator of $ \widetilde{a}_{11}, $ let

(Ⅰ): if $ q_{2}^{2}-4q_{1}q_{3}\leq0 $ or $ 2b\eta\geq\beta(1-m)^{2}(c\beta-ad), $ then $ \widetilde{a}_{11} < 0 $ always holds,

(Ⅱ): if $ q_{2}^{2}-4q_{1}q_{3} > 0, $

when $ 2b\eta < \beta(1-m)^{2}(c\beta-ad), $ we get one negative solution $ x_{1} = \frac{-q_{2}-\sqrt{q_{2}^{2}-4q_{1}q_{3}}}{2q_{1}} $, and one positive solution $ x_{2} = \frac{-q_{2}+\sqrt{q_{2}^{2}-4q_{1}q_{3}}}{2q_{1}} $. Combining with the image, if $ \underline{A} < u_{*} < x_{1} $ or $ x_{2} < u_{*} < \overline{u}, $ then we obtain $ \widetilde{a}_{11} < 0. $ Combining with (4.7) and (4.8), it is not difficult to see, $ \texttt{Tr}(A_{i}) < 0 $ and $ \texttt{Det}(A_{i}) > 0 $ hold, the system (1.5) is locally asymptotically stable.

Theorem 4.2. Under the condition $ x_{1} < u_{*} < x_{2}, $ if

and

$ E_{2}(u_{*}, v_{*}) $ is locally asymptotically stable, where\everymath{ }$ R_{0} = \frac{2b\eta+\eta(c\beta-ad)(1-m)}{\beta(c\beta-ad)(1-m)^{2}}. $

Proof. It follows from Equation (4.7) that

If (4.2) holds, $ \texttt{Tr}(A_{i}) < 0. $ Now, we analyze the sign of $ \texttt{Det}(A_{i}), $

Then we get $ (\mu_{i})_{\texttt{critic}} = \frac{d_{u}\widetilde{a}_{22}+d_{v}\widetilde{a}_{11}}{2d_{u}d_{v}}, $ $ \texttt{Det}(A_{i}) $ takes the minimum value at this point,

Note,

Let $ \tilde{\Delta} = 0 $, we have

Then, we have,

Therefore, $ \tilde{\Delta} = 0 $ has two positive real roots

We get the following conclusions:

(i): if $ \frac{d_{v}}{d_{u}}\leq-\frac{\widetilde{a}_{22}}{\widetilde{a}_{11}} $, it is clear that $ \texttt{Det}(A_{i}) > 0 $ for all $ i\in N_{0} $,

(ii): if $ \frac{d_{v}}{d_{u}} > -\frac{\widetilde{a}_{22}}{\widetilde{a}_{11}}, $ $ \left(\frac{d_{v}}{d_{u}}\right)_{1} < \frac{d_{v}}{d_{u}} < \left(\frac{d_{v}}{d_{u}}\right)_{2}, $ which means $ \tilde{\Delta} < 0 $, and $ \texttt{Det}(A_{i}) > 0 $ for all $ i\in N_{0} $.

The proof is completed.

Theorem 4.3. $ E_{2}(u_{*}, v_{*}) $ is global asymptotically stable when $ R_{1} > 0 $, $ R_{2} > 0, $

where $ R_{1} = b-\frac{\beta(1-m)v_{*}}{u_{*}}-\frac{r}{2}, $ $ R_{2} = \frac{\left(\beta+a\beta(1-m)u_{*}\right)(\eta+a\eta(1-m)u_{*})}{(c\beta+a\eta v_{*})[1+a(1-m)u_{*}][1+a(1-m)\frac{r}{b}]}-\frac{r}{2}. $

Proof. We first give the Liapunov function:

let $ \iota = \frac{\beta+a\beta(1-m)u_{*}}{c\beta+a\eta v_{*}} > 0, $ then

According the maximum principle in ^[36], let $ M $ be a positive constant, such that

where $ (u(x, t), v(x, t)) $ is the solution of (1.4). Moreover, following the Theorem in ^[35] that

On the basis of (4.6), that

We obtain

which means $ E_{2}(u_{*}, v_{*}) $ is global asymptotically stable.

4.2.   Turing instability

With the conclusions we have already drawn in Section 3.3, the ordinary differential equation system is stable if $ \tilde{a}_{11}+\tilde{a}_{22} < 0 $. As we all known, diffusion can make a stable system be unstable, we call this kind of instability as Turing instability^[37]. We assume $ \tilde{a}_{11}d_{v}+\tilde{a}_{22}d_{u} > 0, $ if $ \left(\frac{d_{v}}{d_{u}}\right)_{2} < \frac{d_{v}}{d_{u}}, $ which means $ \tilde{\Delta}(d_{u}, d_{u}) > 0. $ Then there exists positive real part in the following equation

where

under the condition (4.2), $ \texttt{Tr}(A_{i}) < 0 $ is always true. For a fixed $ d_{v} > 0, $

It implies the same conclusion when $ d_{u} > 0 $ is fixed. Then, Turing instability occurs.

4.3.   The system (1.5) with cross-diffusion: stability analysis

In this subsection, we take the cross-diffusion coefficients into our system, we can analyze the stability and give the following equations

Under the condition $ p_{0}-p < 0 $, it is easy to check that $ \texttt{Tr}(A_{i}) < 0 $. Now, let us consider the sign of $ \texttt{Det}(A_{i}). $

It is clear that $ d_{u}d_{v}+d_{12}d_{v}v_{*}+d_{u}d_{21}u_{*} > 0 $ always holds, and assume that $ \widetilde{a}_{11}\widetilde{a}_{22}-\widetilde{a}_{12}\widetilde{a}_{21} > 0 $. When $ d_{21} = 0, $ and $ d_{12}\neq0, $ then the coefficient of primary term in Equation (4.8) becomes

always holds, which means $ d_{12} $ can not change the stability of the system.

When $ d_{21}\neq0 $ and $ d_{12} = 0 $, then the coefficient of primary term in Equation (4.8) becomes

there exists a $ d_{21}^{*} $, when $ d_{21} > d_{21}^{*} $, then the system (1.5) becomes unstable. We obtain that $ d_{21} $ can arise a region of the system, in which Turing instability may appear. The cross-diffusion coefficient $ d_{12} $ is almost the complete opposite to the role of the $ d_{21} $. In the Figure 2, it is clear that as $ d_{21} $ increases, the system (1.5) becomes unstable. By comparing two columns in Figure 2, we also obtain that the ability of cannibalism can affect the stability of system (1.5). The smaller the cannibalism rate, the more likely instability is to occur. modes, and that's almost the complete opposite to the role of the cross-diffusion $ d_{21} $.

5.   Nonconstant positive steady states

We derive the properties of nonconstant positive steady state of system (1.5). A priori estimate and a standard approach based on the Leray-Schauder degree theory are used to obtain the corresponding conclusions.

5.1.   A priori estimate

Theorem 5.1. (Upper Bounds) Assume $ \tilde{M} $ is a fixed positive number, and $ \frac{d_{12}}{d_{u}}, \frac{d_{21}}{d_{v}} < \tilde{M} $, when $ 0 < m < 1-\frac{bd}{(c\beta-ad)r}, $ $ (u(x), v(x)) $ is a positive solution, one has

then the reaction-diffusion system has upper bounds.

Proof. Define $ \Upsilon_{1} = d_{u}u+d_{12}uv, $ $ \Upsilon_{2} = d_{v}v+d_{21}uv $, then we have

Define $ x_{1}\in\Omega $, then applying Harnack Inequality in ^[40] to the equation (2.2), that

One gets $ u(x_{1})\leq\frac{r}{b}, $ $ v(x_{1})\leq\frac{br-a(1-m)r^{2}}{b\beta(1-m)} $. Then we have

It follows from (2.3), that

Then $ u(x_{2})\leq\frac{\eta v(x_{2})}{c\beta(1-m)}, $ $ v(x_{2})\leq\frac{(c\beta-ad)(1-m)r-bd}{\eta b}. $

Theorem 5.2. (Lower Bounds) Let $ \hat{D} > 0 $, and $ \frac{d_{12}}{d_{u}}, \frac{d_{21}}{d_{v}}\leq \hat{D} $, then we can identify a constant $ C^{*} > 0, $ which satisfies $ C^{*}\leq(u(x), v(x))\leq \max\left(\frac{(c\beta-ad)(1-m)r-bd}{\eta b}, \frac{r}{b}\right). $

Proof. Let $ u_{m} = \min u(x) $, according to Maximum principle in ^[39], we have

then, $ \|C_{1}\|_{\infty}\leq C^{*}, $ $ \|C_{2}\|_{\infty}\leq C^{*}, $ other parameters are in system (1.5).

Applying the Maximum principle^[39], we have the following conclusion

where $ \tilde{L} $ is a positive constant. On the contrary, we suppose that there is a sequence of positive solution $ \{u_{i}(x), v_{i}(x)\}_{i = 1}^{\infty} $ satisfies

and $ (u_{i}, v_{i}) $ satisfies

By the regularity theorem for the elliptic equations, we can find some $ \{u_{i}(x), v_{i}(x)\}_{i = 1}^{\infty} $ and two nonnegative functions $ \tilde{u}, \tilde{v}\in C^{2}(\Omega) $, which can make $ (u_{i}(x), v_{i}(x))\rightarrow (\tilde{u}, \tilde{v}) $ in $ [C^{2}(\Omega)]^{2} $ as $ i\rightarrow \infty. $ By (5.4), we have $ \tilde{u}\equiv0 $ or $ \tilde{v}\equiv0 $. Integrating (5.5), we obtain

If $ \tilde{v}\neq0 $, and $ \tilde{u}\equiv0, $ then

when $ i\rightarrow \infty $, $ u_{i}\rightarrow0 $ and $ v_{i} > 0. $ It is a contradiction to the second equation of (5.6).

If $ \tilde{u}\neq0, \tilde{v}\equiv0, $ we have $ \int_{\Omega}\tilde{u}(r-b\tilde{u}) {\rm{d}}x = 0. $ It follows from $ 0 < \tilde{u}\leq\frac{r}{b}. $

Therefore, $ \frac{c\beta(1-m)u_{i}-\eta v_{i}}{1+a(1-m)u_{i}}-d\leq\frac{c\beta(1-m)\frac{r}{b}}{1+a(1-m)\frac{r}{b}}-d < 0, $ $ \frac{c\beta(1-m)r}{b+ar(1-m)}-d < 0 $ which is in contrast to the conclusion that $ \frac{c\beta(1-m)r}{b+ar(1-m)}-d > 0 $ in (2.6). Then there exist lower bounds in the reaction-diffusion system.

5.2.   Non-existence of non-constant positive steady states

We first give the steady-state problem with

Theorem 5.3. Denote $ d^{*} > 0 $, and $ d^{*} = \max\left\{\frac{M_{1}}{\mu_{1}}, \frac{M_{2}}{\mu_{1}}\right\}, $ $ \mu_{1} $ is the smallest positive eigenvalue of the operator $ -\Delta $, such that if $ \min\{d_{1}, d_{2}\}\geq d^{*} $, the system (5.7) has no positive non-constant steady-state. And let

Proof. Define $ (u(x), v(x)) $ be any positive solution of system (5.7), denote $ \bar{u} = \frac{1}{|\Omega|}\int_{\Omega}u(x) {\rm{d}}x. $ Let multiply the first equation of system (5.7) by $ u-\bar{u}, $ we arrive at

By similar arguments as (5.8), one yields

By the Poincare inequality, we obtain

That means if $ \min\{d_{u}, d_{v}\} > \max\left\{\frac{M_{1}}{\mu_{1}}, \frac{M_{2}}{\mu_{1}}\right\}, $ we obtain $ u = \bar{u} $ and $ v = \bar{v}, $ $ (u, v) $ is a constant solution has been proved.

5.3.   Non-constant positive steady states

We always require Assumption 2 or Assumption 3 holds. And the definition of $ \omega $ in this section is different from the previous section, $ d_{21} $ and $ d_{12} $ are enough small. Next, let us explore the existence of nonconstant positive solution of the system (5.7).

Define $ X = \left\{(u, v)^{T}\in [C^{1}(\bar{\Omega})]: \partial_{\nu}u = \partial_{\nu}v = 0\; \texttt{on} \; \partial_{\Omega}\right\}, $ and $ E_{2} = (u_{*}, v_{*}) $ is the positive solution of the system (5.7). $ \Theta = \left\{(u, v)\in X | (c^{-1}\leq u(x), v(x)\leq c, \; \texttt{for}\; \texttt{all} \; x\in \bar{\Omega}\right\}. $ Where $ c $ is a positive constant and $ c $ is guaranteed to exist by upper bound and lower bound. System (5.7) is given by

where

with

Denote $ (I-\Delta)^{-1} $ be the inverse of $ (I-\Delta) $. The system (5.7) can be rewritten as

and $ \omega $ is the positive solution of system (5.9), when it satisfies the equation (5.10).

Theorem 5.2 suggests that $ F(d_{u}, d_{v}, \omega) $ has a unique positive solution $ (u_{*}, v_{*}). $ Then

In fact, it is clear that the eigenvalue $ \lambda(1+\mu_{i}) $ satisfies the following matrix

where $ \textbf{A} = \tilde{\Theta}^{-1}G_{\omega}(u, v). $ Then,

where

If

then $ S(\mu) = 0 $ has two positive roots

Set $ S_{p} = \{\mu_{0}, \mu_{1}, \mu_{2}\ldots\} $ and $ B = B(d_{u}, d_{v}) = \{u\geq0, \; \mu_{-}(d_{u}, d_{v}) < \mu < \mu_{+}(d_{u}, d_{v})\}. $

Let $ m(\mu_{j}) $ be the multiplicity of $ m(\mu_{j}), $ we find that

Clearly,

Lemma 5.1. If $ S(\mu_{i})\neq0 $ for all $ i\geq0 $ and $ \mu_{i}\in S_{p}, $ then$ \rm index $$ (F(\cdot), \omega^{*}) = (-1)^{\sigma} $, where

On the basis of the following Lemma, in order to obtain the index of $ F((\cdot), \omega^{*}) $, we only need to determine the range of $ \mu_{i} $ for which $ S(d_{u}, d_{v}, \mu_{i})\neq0 $.

Theorem 5.4. $ d_{12} $, $ d_{21} $ are fixed and $ d_{21} < d_{21}^{*} $. Under the condition of $ d_{u} > 0 $ and $ a_{11}+a_{22} > 0 $, $ \mu_{-} $ and $ \mu_{+} $ are the two positive solutions, if there exist $ i, j\in N $ such that

and $ \sum\limits_{k = i+1}^{j}m(\mu_{k}) $ is odd, we can find a $ \hat{d}_{v} > 0, $ such that for all $ d_{v}\geq\hat{d}_{v}, $ the system (1.4) has at least one nonconstant positive solution.

Proof. By virtue of the Theorem 5.3, the system (5.10) has constant positive solution when $ d_{u}, d_{v} > {d}^{*} $. From (5.11) and $ \frac{\tilde{a}_{11}}{d_{u}+d_{12}v_{*}}\in (\mu_{q}, \mu_{q+1}), $ we get that when $ d_{v} $ is large enough, then

And there exists a sufficient large $ d_{0} $ such that $ S(\mu) > 0. $ We choose $ {\tilde{d}_{u}} > d^{*} $ such that $ \frac{\tilde{a}_{11}}{\tilde{d}_{u}+d_{12}v_{*}} < \mu_{1} $ and $ {\tilde{d}_{v}} > d^{*}, $ such that

Assume the conclusion is not true, there is no non-constant positive solution for some $ \tilde{d}\geq \hat{d}_{v} $ in the system (5.7). Fixing $ d_{v} = \tilde{d}, $ for $ t\in[0, 1] $, we define

where $ t\in[0, 1] $, and when $ t = 0 $,

when $ t = 1 $,

The solution $ \omega $ in the system (5.9) is equivalent to finding the fixes point of $ \Psi(\omega, 1) $ with $ t = 1. $ We obtain that

Let define $ (u_{*}, v_{*})\in\Theta $ for any solution of system (5.10) on $ \Omega $.

According Leray-Schauder degree, one has

which is a contradiction.

6.   Numerical simulations and conclusions

First, we give the following discretization equations:

the time increment $ \Delta t > 0 $, and space increment $ \Delta x > 0 $. The Discrete initial value $ u_{j}^{0} = \varphi_{1}(x_{j}) = 0 $, $ v_{j}^{0} = \varphi_{2}(x_{j}) = 0, $ and the bounded conditions are $ u_{-1}^{i} = u_{0}^{i}, $ $ u_{M}^{i} = u_{M+1}^{i}, $ $ v_{-1}^{i} = v_{0}^{i}, $ $ v_{M}^{i} = v_{M+1}^{i}, $ $ \Delta t = 0.1, $ and $ \Delta x = 0.5. $ And the relevant images are given. First of all, we set some parameters which satisfy the Descarte's rule of sign, and the Equation (3.2) has one unique solution $ u = u_{*} $, $ v = v_{*}. $

We choose $ d_{u} = 0.1 $, $ d_{v} = 0.1 $; The initial value $ (u_{0}, v_{0}) = \left(1.625+0.2\texttt{sin}\left(\frac{1}{2}\pi x\right); 0.45+0.2\texttt{sin}\left(\frac{1}{2}\pi x\right)\right) $.

In Figure 3a and b, we choose $ m = 0.5, $ the solution of system (1.5) tends to the $ E_{2} $, while $ m = 0.95 $ in Figure 3c and d, the solution tends to the $ E_{1}. $ At the same time, we can obtain that the variation of $ m $ does not make spatial inhomogeneous solution appear. We change the value of $ \eta $. $ \eta = 0.075 $ in Fig. 4a and b, the system exhibits temporal periodic patterns. $ \eta = 0.2 $ in Figure 4c and d, the system becomes stable. $ \eta = 0.8 $ in Figure 4e and f, the system keep stable, and the density of predator $ v(x, t) $ is decreasing. Combining Figure 4 and Figure 1, we are not difficult to find that the cannibalism rate can influence the stability of systems (1.4) and (1.5). To sum up, too small cannibalism rate can destabilize the system, but too large cannibalism rate can make the predator go extinct. Cannibalism has certain positive effects on the stability of the system which is different from our previous cognition. Increasing $ k = 0.2 $ further to $ k = 0.9 $, we can see in Figure 5 that the stability do not change, but the density of the population is decreasing, we know large fear effect can cause the prey and predator decrease but not be extinct. we keep $ k = 0.1 $, $ \eta = 0.2 $, and increase $ m = 0.2 $ to $ m = 0.9 $ in Figure 6, and we observe that the density of the prey is increasing, while the predator is decreasing. When the $ m $ is large enough, the predator may go extinct. From Figure 7a to Figure 7d, we can obtain that the diffusion ratio of predator to prey plays important roles, Turing instability occurs when $ \frac{d_{v}}{d_{u}} $ at a large value, it can break stability of the coexisting state. In Figure 8a, for prey, as $ m $ increases, the number of prey presents an increasing trend which means the prey refuge can protect the prey, the greater the number of prey refuge, the stronger protection the prey will accept. In Figure 8b, for predator, the number of predators present a state that the number of predator is decreasing. We can see obviously from the image that the adding a prey refuge will affect the normal predation behavior of the predator, and reduce the number of predators. It is possible that too much prey refuge could upset the ecological balance and affect the species. In Figure 9a and Figure 9b, prey refuge is beneficial to the prey, and the number of predators is reduced by the protection. In Figure 10, when $ \eta $ is small, the system may be unstable, and if we increase the value of $ m $, it can promote the stability of the system. In Figure 11, we give the bifurcation diagrams for prey and predator for spatial system, it is not difficult to find that the diffusion can extend the time of unstable.

7.   Conclusions

This is a diffusive predator-predator system considering predator cannibalism and taking the effects of predator refuge and fear effect as practical matters. To begin with, we modify the previous systems and establish a new reaction-diffusion system. We give a theoretical analysis for the existence of the positive solutions and stability of constant steady states from Theorem 2.1 to Theorem 2.3. We obtain the sufficient criteria for stability of each steady state. Considering the influence of diffusion factors, we study the Turing instability caused by diffusion, which also produces spatial inhomogeneous patterns when the diffusion coefficients satisfy the condition $ \left(\frac{d_{v}}{d_{u}}\right)_{1} < \frac{d_{v}}{d_{u}} < \left(\frac{d_{v}}{d_{u}}\right)_{2}. $

We discuss the Hopf bifurcation of ODE system, we investigate the Hopf point $ p_{0}, $ and get the solution that when $ p_{0} = p, $ the $ \texttt{Tr}(A_{i}) = 0, $ then Hopf bifurcation could occur. We also verify the results by numerical simulation. In addition, Leray-Schauder degree theory is applied to discuss the constant steady states.

The results in this paper show that:

(1) If the system is based on the assumption of spatially homogeneous and we do not consider the fear effect, predator cannibalism, we can obtain the same conclusion with the system (1.1). When we only consider the prey refuge and fear effect, the fear of predator can affect the density of prey and predator, which is different from the conclusion in ^[14], the decrease of the predator in our system may be due to the predator cannibalism. In non-spatial systems, small predator cannibalism can cause Hopf-bifurcation while a large value of predator cannibalism can make the system be stable. However, the density of prey and predator may decrease.

(2) Both diffusion and cannibalism can disturbance the stability of the spatiotemporal system. From the numerical results that small cannibalism can induce the periodic solution, and large cannibalism can cause the decrease of predator population and prey population; large diffusion rate can make stable system be unstable, and cause the occurrence of Turing instability.

From the images obtained from the numerical simulation, we can conclude that in real life, the purpose of protecting the prey can be achieved by appropriately establishing a prey refuge. For the stability of the ecosystem, we should not build prey refuges indefinitely which may lead to the extinction of predators. Putting a proper amount of predators and increasing the competition between predators are beneficial to the growth of prey and predators. However, if the amount exceeds a certain amount, it will cause the reduction of predators and increase of prey, which may eventually lead to the extinction of predators.

In traditional papers, some scholars considered prey-predator systems without diffusion-reaction, and others considered the systems without predator cannibalism or prey refuge. We improve the previous systems and take more factors into account, we can obtain from the conclusion that diffusion factors, prey refuge and cannibalism do have a huge impact on a prey-predator system. In the models^[42,43], the authors have considered the effect of delay for their system, which made the system more realistic. In ^[44], the authors explored the stochastic partial differential prey-predator system, and in ^[45], Hopf bifurcation has been studied. In the future research, we hope to consider a time delay in the diffusion. We will investigate more interesting and actual systems.

Acknowledgments

This work was supported by the Research Fund for the Taishan Scholar Project of Shandong Province of China, Shandong Provincial Natural Science Foundation of China (ZR2019MA003), the SDUST Innovation Fund for Graduate Students (YC20210217).

Conflict of interest

The authors declare there is no conflict of interest.