Global dynamics and pattern formation for predator-prey system with density-dependent motion

1.   Introduction and main results

The dynamical relationship between predators and their prey is one of the dominant themes in ecology. The origin and theory of predator-prey model is due to pioneer work of Lotka and Volterra ^[1,2]. There have been a lot of studies on the dynamics of this particular type of model through developing various modifications of mathematical models of prey-predator interactions (e.g., ^{[3,4,5,6,7,8]}).

The non-random foraging strategies in the predator-prey dynamics, prey-taxis allows predators to move towards regions of higher prey density and to search more actively for prey. Such a prey-taxis model was derived by Kareiva and Odell in ^[9], and they studied predator aggregation in high prey density areas. It can be described as:

where $u(x, t)$ and $v(x, t)$ denote the population density of predators and preys at position $x$ and time $t$ respectively, and $D$ is a positive constant standing for the diffusion rate of the prey. The terms $\nabla \cdot (d(v) \nabla u)$ and $-\nabla \cdot (u \chi(v) \nabla v)$ account for the diffusion of predators with coefficient $d(v)$ and the prey-taxis with coefficient $\chi(v)$ respectively. $H_1(u, v)$ and $H_2(u, v)$ denote the predator-prey interactions. There mainly are three kinds of typical interspecific interactions: predator-prey, competition and mutualism, which can be represented as

where the functions $f(u)$ and $g(v)$ stand for the intra-specific interactions of predators and prey respectively. The parameters $c_1$ and $c_2$ are positive constants representing the coefficients of inter-specific interactions of predators and prey, and $F(v)$ is the so-called functional response function.

In particular, if $\chi(v) = -d'(v),$ the system (1.1) can be written as

the diffusion term $\Delta (d(v)u)$ with $d'(v) < 0$ is called the "density-suppressed motility" (see ^{[10,11,12,13,14,15]}), which can also characterize the incessant tumbling of cells at high concentration, resulting in a vanishing macroscopic motility. Here $d(v)$ is called the motility function, $d'(v) < 0$ means that the predator reduce its motility when encountering the prey. When $H_1(u, v) = 0$, $H_2(u, v) = u-v$ and $d(v) = c_0 v^{-k}$ decays algebraically in $v$, the solution may exist globally in two or higher dimensions. For example, Yoon and Kim in ^[16] proved that system (1.2) has a unique global bounded classical solution for any $k > 0$ under a smallness assumption on $c_0$ in any dimensions. The only global existence result without smallness assumptions was recently given by Ahn and Yoon ^[17]. Under the assumptions that there exist positive constants $\gamma_1, \gamma_2, \gamma_3$ such that $\gamma_1\leq d(v)\leq \gamma_2$ and $|d'(v)|\leq \gamma_3,$ Tao and Winkler in ^[18] proved the existence of global classical solutions in the 2-dimensional case and global weak solutions in 3-dimensions. While if $d(v)$ decays exponentially, the solution may blow-up in two dimensions with a critical mass, see ^[19,20,21] and so on. For $H_1(u, v)\neq 0$ with logistic growth on $f(u) = \mu u(1-u)$, there also are many interesting results. The global existence and asymptotic behavior of solutions was first established by Jin, Kim and Wang in ^[12] under certain conditions on $d(v)$ in two dimensions, which has been developed by many authors, please refer to ^{[19,22,23,24,25,26,27,28]} and so on.

In the above mentioned references, the authors assumed $d'(v) < 0$. While under normal circumstances, the predators will increase their motilities and keep chasing (such as wolves and sheep) when encountering the prey until they succeed. Therefore, it is meaningful for us to consider the case $d'(v)\geq 0$, and we shall consider general case for $d(v)$ without monotonicity assumptions in this paper.

In this paper, we consider the following density-dependent predator-prey system:

where $a_1, a_2 > 0$ represent the intrinsic growth rates of species, $b_1, b_2 > 0$ are the death rates due to intra-specific competition and $\alpha > 0$ denotes the intrinsic predation rate. $F(v)$ is the so-called functional response function accounting for the intake rate of predators as a function of prey density, and $d(v)$ is the motility function as we mentioned above. The most common types $F(v)$ in the literature are

● $F(v) = v$ (Lotka-Volterra type or Holling type Ⅰ);

● $F(v) = \frac{v}{\lambda +v}$ (Holling type Ⅱ);

● $F(v) = \frac{v^m}{\lambda^m+ v^m}$ (Holling type Ⅲ) with constants $\lambda > 0$ and $m > 1.$

We make the following assumptions throughout the whole paper:

(H1) $d(v) \in C^3([0, \infty))$ and $d(v) > 0$ on $[0, \infty).$

(H2) $F(v) \in C^1([0, \infty)), F(0) = 0, F(v) > 0$ in $(0, \infty)$ and $F'(v) > 0.$

In (1.3), the predator is called the specialist predator if $a_1 < 0$, Jin and Wang investigate the global boundedness, asymptotic stability and pattern formation of system (1.3) ^[29]. They study the dynamic behaviors of the predator and the prey under the condition $d'(v) < 0$. In this paper, we assume $a_1 > 0$ (the corresponding predator is called the generalist predator) and make no assumptions on the monotonicity of $d(v).$ Usually, a generalist species is able to thrive in a wide variety of environmental conditions and can make use of a variety of different resources (for example, a heterotroph with a varied diet). A specialist species can thrive only in a narrow range of environmental conditions or has a limited diet. Most organisms do not all fit neatly into either group. Some species are highly specialized (the most extreme case being monophagous, eating one specific type of food), others less so, and some can tolerate many different environments. In other words, there is a continuum from highly specialized to broadly generalist species. We mainly focus on exploring the global dynamics and spatial-temporal patterns for generalist predators with density-dependent motion for more general motility functions $d(v).$ We mention here Nakashima and Yamada in ^[30] studied the existence of positive solutions for boundary value problems of nonlinear elliptic systems which arise in the study of the Lotka-Volterra prey-predator models with cross-diffusion in the special case $d(v) = 1+\alpha v.$

We first derive the global boundedness and existence results for the classical solutions to the system (1.3).

Theorem 1.1 (Global boundedness). Let $\Omega\subset \mathbb{R}^2$ be a bounded domain with smooth boundary and the hypotheses (H1)-(H2) hold. Assume $(u_0, v_0)\in [W^{1, p}(\Omega)]^2$ with $p > 2$ and $u_0, v_0\gneqq 0$. Then the problem (1.3) has a unique global classical solution $(u, v) \in [C^0(\bar{\Omega}\times[0, \infty))\cap C^{2, 1}(\bar{\Omega}\times(0, \infty))]^2$ satisfying $u, v > 0$ for all $t > 0$. Furthermore there exists a constant $C > 0$ independent of $t$ such that

For the global stability, except for the hypotheses (H1)-(H2), we also need the following hypothesis on the compound function

(H3) The function $\phi(v)$ is continuously differentiable on $(0, \infty),$ $\phi(0) = \mathop{\lim}\limits_{v\to 0} \phi (v) > 0$ and $\phi'(v)\leq 0$ for any $v\geq 0.$

Remark 1.1. We remark that the hypothesis (H3) is not stringent, and can be satisfied by many forms by imposing some conditions on the parameters if needed. For example, if $F(v)$ is of Holling type Ⅰ or Holling type Ⅱ with $a_2\leq b_2 \lambda$, then (H3) is automatically satisfied. In general, if (H3) is violated, pattern formations such as periodic orbits or non-constant steady state may arise (see ^[40]).

Another relevant question is whether the interacting predator-prey population will arrive at the coexistence, exclusion or extinction eventually, which is always an important topic in population dynamics. One can easily compute that the system (1.3) has four possible steady states:

where $(u_*, v_*)$ satisfies

By constructing suitable Lyapunov functionals, we can obtain the following global stability of the coexistence steady state $(u_*, v_*)$ and semi-trivial steady state $\left(\frac{a_1}{b_1}, 0\right)$. In the context, we define

Theorem 1.2 (Global stability). Let $\Omega\subset \mathbb{R}^2$ be a bounded domain with smooth boundary and the hypotheses (H1)-(H2) hold. Assume $(u_0, v_0)\in [W^{1, p}(\Omega)]^2$ with $p > 2$ and $u_0, v_0\gneqq 0$ and $(u, v)$ is the solution of (1.3) obtained in Theorem 1.1.

(1) If (H3) holds and the parameters satisfy

then

and $(u_*, v_*)$ satisfies (1.5), where $K$ is defined in (1.6). Moreover, there exist some positive constants $\sigma, T_0$ and $C$ independent of $t$ such that

(2) If the parameters satisfy

then

where $\xi_1$ and $K_0$ are defined by $\xi_1 = \frac{1}{\alpha}+\frac{b_1b_2}{K_0^2\alpha^2}$ and $K_0 = \mathop{\max}\limits_{\bar\Omega}\frac{F(v)}{v}$ respectively. Moreover, there exist some positive constants $T_1$ and $C$ independent of $t$ such that

Remark 1.2. For the semi-trivial steady state $\left(0, \frac{a_2}{b_2}\right)$, in the special case of $F(v) = v$ (Holling type Ⅰ), we have showed that it is linear unstable in section 5. For more general $F$, the globally stability for $\left(0, \frac{a_2}{b_2}\right)$ is nontrivial and has to be left open in the current paper.

In the proof of global existence, the method used in ^[31] was based on a priori estimates for the energy functional $\int_\Omega u\ln u dx+ \int_\Omega |\nabla v|^2 dx$ to attain the $L^2$ estimates of the solutions. However, such a method of a priori estimates is only applicable for the case where the motility function $d(v)$ is constant. Therefore, the method in ^[31] is not adaptable to the model (1.3). In this paper, we first derive the $L^2$ estimates for $|\nabla v|$ and then directly establish the $L^2$ estimates for $u$ by the Gagliardo-Nirenberg inequality and regularity lemmas, in which we also need to prove the boundedness of $\int_t^{t+\tau}\int_\Omega u^2 dxds$ and $\int_t^{t+\tau}\int_\Omega |\Delta v|^2 dxds$. In the proof, we do not use the property of the self-adjoint realisation of $-\Delta+\delta$ (see ^[29]). Finally, we derive the boundedness for $u$ by the Moser iteration technique.

If $d(v)$ is constant, the system (1.3) has been studied from many aspects as we mentioned above. If $d(v)$ is non-constant as considered in this paper, we find that the system (1.3) can generate pattern formation as presented in section 5 under the condition $d'(v) > 0$. The pattern formation is obviously different from the ones in ^[29], in which the authors studied the case $d'(v) < 0$. In our case $a_1 > 0,$ the corresponding predator is the generalist predator, the pattern formation may not occur when $d'(v)\leq 0$ by linear instability analysis and numerical simulations.

The paper is organized as follows. In section 2, we present the local existence theorem with some preliminary results, and we derive the boundedness of $\|v(\cdot, t)\|_{L^\infty}$ and $\|\nabla v(\cdot, t)\|_{L^2}.$ In section 3, we derive the boundedness of $\|u(\cdot, t)\|_{L^\infty}$ by the technique of Moser iteration. In section 4, we construct suitable Lyapunov functionals and use the LaSalle invariance principle to prove the global stability and convergence rate stated in Theorem 1.2. In section 5, we further explore time-periodic patterns by linear instability analysis and numerical simulations.

2.   Local existence and preliminaries

In the sequel, we shall use $C$ or $C_i$ to denote a positive generic constant which may vary in the context. Without confusion, the integration variables $x$ and $t$ will be omitted, for instance $\int_0^a \int_{\Omega} f(x, t)dxdt$ will be abbreviated as $\int_0^a \int_{\Omega} f(x, t)$. Often $\|f\|_{L^p(\Omega)}$ will be written as $\|f\|_{L^p}$. The existence and uniqueness of local solutions to (1.3), which can be readily proved by the Amann theorem ^[32,33] or the well-established fixed pointed argument together with the parabolic regularity theory ^[12].

Lemma 2.1 (Local existence). Let the assumptions in Theorem 1.1 hold. Then there exists a constant $T_{max}\in (0, \infty]$ such that the problem (1.3) admits a unique classical solution

satisfying $u, v > 0$ for all $t > 0$. Moreover,

Proof. Denote $z = (u, v)$. Then the system (1.3) can be written as

where

Since the given initial value of $(u_0, v_0)$ are nonnegative and satisfy $0\leq (u_0, v_0)\in [W^{1, p}(\Omega)]^2$ with $p > 2$, and hence the matrix $P(z)$ is positive definite at $t = 0.$ This means that the system (1.3) is uniformly parabolic. Then the application of [33,Theorem 7.3] yields a $T_{max} > 0$ such that the system (1.3) possesses a unique solution $(u, v) \in [C^0(\bar \Omega \times [0, T_{max}))\cap C^{2, 1}(\bar \Omega \times (0, T_{max}))]^2$.

Next, we prove the positivity of $u$ and $v.$ To this end, we rewrite the first equation of (1.3) as

Applying the strong maximum principle to (2.2) with the Neumann boundary condition deduces that $u > 0$ for all $(x, t)\in \Omega \times (0, T_{max})$ due to the fact $u_0 \gneqq 0$. In a similar way, we can prove $v > 0$ for any $(x, t) \in \Omega \times (0, T_{max})$ by using the second equation in (1.3). In addition, since $P(z)$ is an upper triangular matrix, the blow-up criterion (2.1) follows from [35,Theorem 5.2] directly. Then the proof of Lemma 2.1 is completed.

Lemma 2.2. Let the assumptions in Theorem 1.1 hold. Then the solution of (1.3) satisfies

for all $t > 0,$ where $K$ is defined in (1.6), and

Proof. The proof is the similar to [31,Lemma 2.2], but for readers' convenience, we list the proof here. Since $u, v$ and $F(v)$ are nonnegative, we can derive from the second equation of (1.3) that

Let $\bar v(t)$ be the solution of the following ODE problem

Then we obtain from (2.6) that $\bar v(t)\leq K: = \max\{\frac{a_2}{b_2}, \|v_0\|_{L^\infty}\}.$ It is obvious that $\bar v(t)$ is one of the super-solution of the following PDE problem

and therefore, we have

where we have used the strong maximum principle to derive $V > 0.$ Combining (2.5), (2.7) with (2.8) and employing the comparison principle, one has

This proves (2.3).

In addition, since $v(a_2-b_2v) < 0$ for $v > K$, we can further deduce from (2.6) that

which together with (2.9) gives (2.4).

Lemma 2.3. Let the assumptions in Theorem 1.1 hold. Then the solution of (1.3) satisfies

Moreover, one has

where

Proof. Multiplying the second equation of (1.3) by $\alpha$ and adding the result into the first equation of (1.3), then integrating the result over $\Omega,$ we obtain

which implies that

Applying the Grönwall inequality to (2.13) yields

it follows that (2.10) is valid. Then integrating (2.13) over $(t, t+\tau)$, and using (2.14) to obtain

which implies (2.11).

By Lemma 2.3, we establish the estimates for $\|\nabla v\|_{L^2}$ and $\int_t^{t+\tau}\int_\Omega |\Delta v|^2 dxds.$

Lemma 2.4. Let the assumptions in Theorem 1.1 hold. Then there exists a constant $C > 0$ independent of $t$ such that

and

where $\tau$ and $\tilde T_{max}$ are defined in (2.12).

Proof. Multiplying the second equation of (1.3) by $-\Delta v$, integrating the result in $\Omega,$ using the assumption (H2) and the boundedness of $v$ (see (2.3)), we have

which implies

where $K$ is defined in (1.6).

Applying the Gagliardo-Nirenberg inequality, [29,Lemma 2.5] and noting the fact $\|v\|_{L^2}\leq K |\Omega|^{\frac{1}{2}}$ yields

Substituting (2.18) into (2.17), one has

which together with (2.11) yields (2.15).

Then integrating (2.19) over $(t, t+\tau)$, we obtain (2.16).

3.   Boundedness of solutions

In this section, we prove the boundedness of $\|u\|_{L^\infty}$ by the technique of Moser iteration. To the end, we first prove the boundedness of $\|u\|_{L^2}.$

Lemma 3.1. Let the assumptions in Theorem 1.1 hold. Then there exists a constant $C > 0$ independent of $t$ such that

Proof. Multiplying the first equation in (1.3) by $2u$ and integrating the result with respect to $x$ in $\Omega,$ one has

By the assumptions in (H1), (H2), and (2.3), we know that $d(v)\in C^3$ and $0 < v\leq K,$ where $K$ is defined in (1.6). Therefore, one has $d(v)\geq C_1,$ $0 < F(v)\leq F(K)$ and $|d'(v)|\leq C_2,$ then it follows from (3.2) that

Applying the Galiardo-Nirenberg inequality, one has

and

where we have used the boundedness of $\|\nabla v\|_{L^2}$ (see (2.15)) and [29,Lemma 2.5]. Combining (3.4) with (3.5) and applying the Young inequality yields

Then substituting (3.6) into (3.3), we obtain

For any $t\in(0, T_{max})$, by (2.10), there exists a nonnegative $t_0\in ((t-\tau)_+, t)$ with $\tau = \min\{1, \frac{1}{2}T_{max}\}$ such that

and by (2.16), one has

Then integrating (3.7) on $(t_0, t)$, and applying (3.8), (3.9) and the fact $t\leq t_0+\tau\leq t_0+1$, we have

which indicates (3.1) and thus completes the proof.

To derive the boundedness for $\|u\|_{L^\infty},$ we need the following regularity lemma.

Lemma 3.2. (^[37]) Assume that $\Omega \subset \mathbb{R}^n$ is a bounded domain with smooth boundary. Suppose that $y(x, t) \in C^{2, 1}(\bar\Omega \times (0, T_{max}))$ is the solution of

where $\phi(x, t) \in L^\infty ((0, T_{max}); L^p(\Omega))$. Then there exists a constant $C > 0$ such that

with

We will prove the boundedness of $u$ through the Moser iteration procedure.

Lemma 3.3. Let the assumptions in Theorem 1.1 hold. Then there exists a positive constant $C$ independent of $t$ such that

Proof. Multiplying the first equation of the system (1.3) by $u^{p-1}$ with $p\geq 2,$ and integrating the resulting equation by parts, we derive

Noting the fact $d(v)\geq C_1,$ $0 < F(v)\leq F(K)$ and $|d'(v)|\leq C_2,$ and using the Young inequality, we obtain from (3.11) that

which gives

for all $t\in(0, T_{max})$ and $p\geq 2.$ By Lemma 3.1, we have $\|u(\cdot, t)\|_{L^2}\leq C_3$, and thus we derive $\|\nabla v(\cdot, t)\|_{L^4}\leq C_4$ from Lemma 3.2. Then applying the Gagliardo-Nirenberg inequality and the Hölder inequality yields

and

Substituting (3.13) and (3.14) into (3.12), we obtain

by the Grönwall inequality, we derive

Let $p = 4$ in the above inequality and apply Lemma 3.2 again. We obtain that there exists a constant $C_7 > 0$ such that

Then using the Moser iteration procedure (see ^[34]), one derives (3.10) and thus proves Lemma 3.3.

Proof of Theorem 1.1. Theorem 1.1 is a consequence of Lemma 2.1, Lemma 2.2 and Lemma 3.3.

4.   Globally stability

In this section, we will investigate the asymptotical behavior of solutions solving system (1.3) and prove Theorem 1.2 based on Lyapunoval functional method along with Barbălat's lemma.

Let us first recall the basic Barbălat's lemma.

Lemma 4.1. (Barbălat's Lemma ^[36]) Suppose that $h: [1, \infty) \to \mathbb{R}$ is a uniformly continuous function such that $\lim_{t \to \infty} \int_1^t h(s)ds$ exists, then $\lim_{t \to \infty}h(t) = 0$.

Now we give the uniform estimates of the global solution.

Lemma 4.2. Let $(u, v)$ be the unique global bounded classical solution of (1.3) given by Theorem 1.1. Then for any given $0 < \alpha < 1,$ there exists a constant $C(\alpha) > 0$ such that

Proof. This proof is based on the standard regularity for parabolic equations. For readers' convenience, we sketch the proof here. Due to the boundedness of $(u, v),$ applying the interior $L^p$ estimate ^[39] to (1.3), we derive that

Using the Sobolev embedding theorem and (4.2), we derive

Applying (4.3) and the Schauder estimate ^[38] to the second equation of (1.3), we obtain

which implies

Rewrite the first equation in (1.3) as

where

Due to (4.3) and (4.4), we see that

Applying the Schauder estimate to (4.6) gives $\|u\|_{C^{2+\alpha, 1+\frac{\alpha}{2}}(\bar\Omega \times [i+1, i+3])} \leq C_6$ for all $i \geq 0.$ Thus

Then (4.1) follows from (4.5) and (4.7). This completes the proof of Lemma 4.2.

Next we shall prove the global stability of the coexistence steady state $(u_*, v_*)$ by constructing the following Lyapunov functional:

where $(u_*, v_*)$ satisfies (1.5).

Lemma 4.3. Let $\Omega\subset \mathbb{R}^2$ be a bounded domain with smooth boundary and the hypotheses (H1)-(H3) hold. Assume $(u_0, v_0)\in [W^{1, p}(\Omega)]^2$ with $p > 2$ and $u_0, v_0\gneqq 0$. If $(u, v)$ is the solution of (1.3) obtained in Theorem 1.1, and the parameters satisfy

then

where $(u_*, v_*)$ satisfies (1.5), and $K$ is defined in (1.6).

Proof. First by a similar argument as [31,Lemma 4.3], we deduce that

We compute the derivative of $V(t)$ by using the system (1.3) to derive

Then $I_1$ can be rewritten as

where $X^T$ denotes the transpose of $X: = (\nabla u, \nabla v)$, and

It is easy to check that the matrix $A$ is nonnegative definite if and only if

and thus

Next, we compute $I_2.$ By (1.5), we see that

we obtain from (H2) and (H3) that there exists a constant $T_1 > 0$ such that

Therefore, there exists a constant $C_1 > 0$ such that

where $\phi(v) = \frac{v(a_2-b_2v)}{F(v)}$ is defined in (1.4), $\xi_1$ and $\xi_2$ are lying between $v$ and $v_*.$ Combining (4.10), (4.11) with (4.12) gives

with $E(t) = \int_\Omega [(u-u_*)^2 + (v-v_*)^2]dx$.

Since $V(t) \geq 0,$ we have

It follows from the regularity of $u, v$ that $E(t)$ is uniformly continuous in $[1, \infty).$ An application of Lemma 4.1 yields

By Lemma 4.2, we derive that $u(\cdot, t)$ and $v(\cdot, t)$ are bounded for $t > 1$ in the space $W^{1, \infty}(\Omega).$ Applying the Gagliardo-Nirenberg inequality

to $u-u_*$ and $v-v_*$, respectively, we can obtain (4.9) from (4.14). This completes the proof of Lemma 4.3.

Lemma 4.4. Suppose that the conditions of Lemma 4.3 hold. Then there exist some positive constants $\sigma, T_0$ and $C$ independent of $t$ such that

Proof. By Lemma 4.3, we have

then applying L'Hôpital's rule to derive

Therefore, we obtain from the continuity that there exists a constant $t_1 > 0$ such that

In addition, by Lemma 4.1 in ^[31], we find that there exist some constants $C_1, C_2 > 0$ and $t_2 > 0$ such that

Combining (4.16), (4.17), and the definition of $V(t)$ in (4.8), we conclude that there exist some positive constants $C_3, C_4$ and $t_3$ such that for all $t > t_3$

This, together with (4.13) yields a constant $C_5 > 0$ such that

therefore, we derive that there exist some constants $C_6 > 0$ and $C_7 > 0$ such that

To finish the proof, we need the $L^\infty$ estimates for $u-u_*$ and $v-v_*.$ By Lemma 4.2 and the Gagliardo-Nirenberg inequality, we have

by a similar argument, one has

which, combined with (4.18) and (4.19), gives (4.15). Therefore, we complete the proof of Lemma 4.4.

Lemma 4.5. Let $\Omega\subset \mathbb{R}^2$ be a bounded domain with smooth boundary and the hypotheses (H1)-(H2) hold. Assume $(u_0, v_0)\in [W^{1, p}(\Omega)]^2$ with $p > 2$ and $u_0, v_0\gneqq 0$. If $(u, v)$ is the solution of (1.3) obtained in Theorem 1.1, and the parameters satisfies

then

where $K$ is defined in (1.6), $\xi_1$ and $K_0$ are defined by $\xi_1 = \frac{1}{\alpha}+\frac{b_1b_2}{K_0^2\alpha^2}$ and $K_0 = \mathop{\max}\limits_{\bar\Omega}\frac{F(v)}{v}$ respectively.

Next we prove the global stability of the generalist predator only steady state $(\frac{a_1}{b_1}, 0)$ by constructing the following Lyapunov functional:

Proof. It is obvious that $V_1(t)\geq 0$ from the definition. Next we compute the derivative of $V_1(t)$ by using the system (1.3) and (4.20) to obtain

Then $J_1$ and $J_2$ can be rewritten as

where $X_1^T$ and $Y_1^T$ denotes the transpose of $X_1: = \left(\frac{\nabla u}{u}, \nabla v\right)$ and $Y_1: = \left(u-\frac{a_1}{b_1}, v\right)$, and

and

By the definition of $\xi_1$ and $K_0$ defined in Lemma 4.5, and (4.20), we derive that

and

It follows that

with $E_1(t) = \int_\Omega\left(u-\frac{a_1}{b_1}\right)^2dx+\int_\Omega v^2dx.$

By a similar argument as Lemma 4.3, we obtain (4.21) and thus completes the proof of Lemma 4.5.

By a similar argument as Lemma 4.4, we can derive the following decay rate estimates.

Lemma 4.6. Suppose that the conditions in Lemma 4.5 hold. Then there exist some positive constants $T_1$ and $C$ independent of $t$ such that

Proof. Since

then combining this with (4.23), we obtain

Solving this ordinary differential inequality, for sufficiently large $t$, we arrive at

Using the same argument as in the proof of Lemma 4.4, we readily get (4.24) and thus complete the proof of Lemma 4.6.

5.   Applications and spatio-temporal patterns

In this section, we shall study the possible pattern formation generated by the system (1.3). As a typical example, we consider the case Holling type Ⅰ of $F(v) = v.$

5.1.   Linear instability analysis

To begin with, we first consider the linear stability of the system (1.3). We linearise the system (1.3) at an equilibrium $(u_s, v_s)$ and write the linearised system as

where $T$ denotes the transpose and

as well as

Let $W_k(x)$ be the eigenfunction of the following eigenvalue problem:

the allowable wavenumbers $k$ are discrete in a bounded domain, for instance, if $\Omega = (0, l),$ then $k = \frac{n \pi}{l}, \, \, n = 0, 1, 2, \cdots$. Since the system (5.1) is linear, the solution of $\Phi(x, t)$ has the form of

where $\rho$ is the temporal eigenvalue and $c_k (k\geq 0)$ are determined by the Fourier expansion of the initial conditions in terms of $W_k(x).$ Substituting (5.2) into (5.1) yields

which indicates that $\rho$ is the eigenvalue of the matrix $M_k$ with

Then the eigenvalue $\rho(k^2)$ satisfies

where

and

Next we proceed to consider the stability of equilibria $(0, 0)$, $\left(\frac{a_1}{b_1}, 0\right),$ $\left(0, \frac{a_2}{b_2}\right)$ and $\left(\frac{a_2\alpha+a_1b_2}{\alpha+b_1b_2}, \frac{a_2 b_1-a_1}{\alpha+b_1b_2}\right)$ in the presence of spatial structure.

● For the steady state $\left(0, 0\right)$, it can be easily check that $\rho$ satisfies

At least one of the roots is positive, therefore, $(0, 0)$ is always linearly unstable.

● For the steady state $\left(\frac{a_1}{b_1}, 0\right)$, it can be easily check that $\rho$ satisfies

Therefore if $a_2 = u_s,$ $\left(\frac{a_1}{b_1}, 0\right)$ is marginally stable. If $a_2 > u_s$, then there exists a $k$ such that one of the root is negative and the other is positive, and $\left(\frac{a_1}{b_1}, 0\right)$ is linearly unstable.

● For the steady state $\left(0, \frac{a_2}{b_2}\right)$, $\rho$ satisfies

Therefore there exists a $k$ such that one of the root is negative and the other is positive, and $\left(0, \frac{a_2}{b_2}\right)$ is linearly unstable.

● For the coexistence steady state $(u_s, v_s) = \left(\frac{a_2\alpha+a_1b_2}{\alpha+b_1b_2}, \frac{a_2 b_1-a_1}{\alpha+b_1b_2}\right)$, $\rho$ satisfies

with $M_1 = -k^2d(v_s)-b_1u_s, M_2 = -k^2 u_sd'(v_s)+\alpha u_s, M_3 = -v_s$ and $M_4 = -k^2-b_2v_s,$ where we have used the identity

Since $M_1+M_4 < 0,$ $(u_s, v_s)$ is linearly unstable iff

One root is negative and the other is positive. That is, we need the condition

which indicates that

We conclude that a steady-state bifurcation may occur if

and there are allowable wavenumbers $k$ such that

where $k_1^{\pm} = \frac{u_s v_sd'(v_s)-b_2 v_sd(v_s)-b_1u_s)\pm \sqrt{(b_2 v_sd(v_s)+b_1u_s-u_s v_sd'(v_s))^2-4(b_1b_2+\alpha) d(v_s) u_s v_s}}{2d(v_s)}.$

We remark here that when $d'(v)\leq 0$, the inequality (5.3) can not be hold for any $k$. Therefore, the coexistence steady state $(u_s, v_s)$ is linearly stable if $d'(v)\leq 0.$

5.2.   Spatio-temporal patterns

In this subsection, we take some examples to present the periodic patterns.

According to the condition (5.3) and the linear stability analysis in subsection 5.1, we fix the value of the parameters in all simulations as follows:

Then we obtain the possible steady states are $(u_s, v_s) = \left(\frac{3}{2}, \frac{1}{4}\right)$ or $(1, 0)$ or $(0, 1)$.

The numerical simulations of patterns are then shown in the following Figures 1–3.

Case 1. $(u_s, v_s) = (\frac{3}{2}, \frac{1}{4}).$ In this case, $k = \frac{n \pi}{10}, n = 0, 1, 2, \cdots,$ the condition (5.3) turns into

it can be checked that the above inequality is valid iff $1\leq n\leq 8$.

The numerical spatial-temporal patterns generated by the model (1.3) under the condition (5.4) are plotted in Figure 1(a), where we observe the spatially inhomogeneous temporal periodic patterns arising from the vicinity of equilibrium $(\frac{3}{2}, \frac{1}{4})$. The time distributions of predators and prey at a fixed space position plotted in Figure 1(b) show that predators and prey are periodic. We also plot the spatial distributions of predators and prey at a fixed time in Figure 1(b).

Case 2. $(u_s, v_s) = (1, 0).$ In this case, from the above linear analysis, we know that $\rho$ satisfies

therefore, then one root is positive and the other is negative for the above equation iff $n = 0, 1, 2, 3$.

Similarly, we observe the spatially inhomogeneous temporal periodic patterns arising from the vicinity of equilibrium $(1, 0)$. The time distributions of predators and prey at a fixed space position plotted in Figure 2(b) show that predators and prey are periodic. We also plot the spatial distributions of predators and prey at a fixed time in Figure 2(b).

Case 3. $(u_s, v_s) = (0, 1).$ In this case, we know that $\rho$ satisfies

therefore, then one root is positive and the other is negative for the above equation iff $n = 0$.

We also can observe the spatially inhomogeneous temporal periodic patterns arising from the vicinity of equilibrium $(0, 1)$. The time distributions of predators and prey at a fixed space position plotted in Figure 3(b) show that predators and prey are periodic. We also plot the spatial distributions of predators and prey at a fixed time in Figure 3(b).

In summary, we conclude that the time-periodic patterns have been obtained when $d(v)$ is monotone increasing and satisfies some conditions. While from the linear stability analysis, we know that the coexistence steady state and the semi-trivial steady states are both linearly stable. These results indicate that the motility function $d(v)$ can trigger pattern formation and is a factor inducing the spatial heterogeneity of populations. In addition, if $a_1 < 0,$ ($u$ represents the specialist predator) Jin and Wang in ^[29] have proved that the semi-trivial steady state $\left(0, \frac{a_2}{b_2}\right)$ is linearly stable if $\alpha \frac{a_2}{b_2}\leq -a_1.$ While our results indicate that when $u$ is a generalist predator ($a_1 > 0$), the semi-trivial steady state $\left(0, \frac{a_2}{b_2}\right)$ is linearly unstable, because the generalist predator $u$ can gain food from other preys.

Acknowledgments

The authors are grateful to the referee for his/her valuable comments, which greatly improved the exposition of our paper. The research of T. Feng was partially supported by the Special Basic Cooperative Research Programs of Yunnan Provincial Undergraduate Universities Association (Grant No. 2019F001-078, 202101BA070001-132). The research of L. Wu was supported by the 2020 Hong Kong Scholars Program (Project ID P0031250 and Primary Work Programme YZ3S).

Conflict of interest

The authors declare there is no conflict of interest.