Predator-prey systems with defense switching and density-suppressed dispersal strategy

1.   Introduction

Defense switching means that prey species pay more attention on guarding against the relatively more abundant population predator ^[1], certain fish species in Lake Tanganyika against two phenotypes (dextral and sinistral) of cichlid Perissodus microlepis ^[2] is a typical example. The dextral and sinistral phenotypes attack the prey fishes from the left-side and right-side, respectively. Pretend that the population of dextral individuals is more abundant, then prey fishes tend to be more defensive against the attacks from left-side, which leads to greater hunting success for sinistral individuals (relatively rare population). Based on a simple Lotka-Volterra equations of a two-predator one-prey system, Saleem et al. proposed a defensive switching model ^[1]

with the predatory rates functions

where $u: = u(x, t)$ and $v: = v(x, t)$ denote the density of two predators while $w: = w(x, t)$ is the prey density. The parameters $\alpha_j\; (j = 1, 2)$ account the death rate, $\beta_j > 0\; (j = 1, 2)$ are the predation coefficients and $\sigma > 0$ accounts for the growth rate of the prey species. Moreover, the predatory rates $f_1$ and $f_2$, also called "defensive switching functions", possess a characteristic property that the rate of the prey attacked by a predator will decrease if this predator population becomes much more abundant than the population of another predator. Specifically, when a predator population becomes large, the prey species guards against it more vigilantly and switches to another predator, which is in relatively small population, to keep its individual from being hunted too much. Such prey behaviors result in less successful hunting for abundant population predator species and more successful hunting for relatively rare one ^[1].

In light of the defensive switching model Eqs (1.1) and (1.2) in ^[1], Pang and Wang ^[3] considered the following reaction-diffusion system by introducing the random movements of species and the intra-specific interaction between the prey in Eq (1.1)

where constants $d_j > 0\; (j = 1, 2, 3)$ account for diffusion coefficients and the positive parameter $K$ represents the environmental carrying capacity to the prey species. We note that the interaction mechanism between the predators and prey in the defensive switching model is substantially different from that in the ratio-dependent predator-prey system ^[4] though they seem to have some similar structures.

In ^[1], Saleem et al. proved that the co-existence steady state is globally asymptotically stable except the case where the two predators have the same mortality rates; Otherwise, the system has a periodic solution. Pang and Wang ^[3] proved the co-existence steady states is globally asymptotically stable no matter $\alpha_1$ and $\alpha_2$ are equal or not and hence no pattern formation will arise from the system (1.3). If the term $d_1\Delta u$ is replaced by a cross diffusion term $\Delta (d_1u+\frac{ku}{\epsilon+v^2})\; (k, \epsilon > 0)$, however, they showed that the cross diffusion between predators can drive stationary patterns by using Leray-Schauder degree theory. Recently, the system (1.3) with prey-taxis (the cross diffusion between the predators and prey) was considered by Wang and Guo ^[5]. They established the existence of globally bounded solutions and global stability of constant steady states for small prey-tactic coefficient and numerically demonstrated that strong prey-taxis can induce the pattern formation if the two predators have different mortality rates. Subsequently, the existence of nonconstant steady states was obtained for some range of repulsive prey-taxis coefficient with more general functional response functions in ^[6]. From the above results, we find that the dispersal strategies (like cross-diffusion or prey-taxis) play important roles in determining the population distribution profiles.

In this paper, we shall consider a different dispersal strategy - density-suppressed diffusion, which was first used in ^[7] to describe the directed movement of predators in the predator-prey systems to fit the experimental observations. This type of diffusion assumes that the predator's diffusion decreasingly depends on the density distribution of the prey. As shown in ^[8], the density-suppressed diffusion can explain the heterogenous population distribution observed in the field experiment of ^[7] while random diffusion can not. Hence density-suppressed diffusion is a dispersal strategy employed in the predator-prey system. By taking into account the density-suppressed diffusion, the defensive switching model reads

where $d_j(w) > 0$ and $d'_j(w) < 0\; (j = 1, 2)$. Note that the property $d'_j(w) < 0$ means that predators will decrease their random diffusion rates at higher density of the prey species in order for predation. If we expand the density-suppressed diffusion

we find that the density-suppressed diffusion indeed intrinsically includes both random diffusion and advection (prey-taxis) components. The differences from the spatial models considered in ^[3,5,6] is that here both diffusion and advection coefficient are not constant but functions of the prey density. It was shown that the density-suppressed diffusion in the predator-prey systems may generate spatially non-homogeneous patterns that the random diffusion can not do ^[8] and may bring substantially different dynamics ^[9].

Beyond the predator-prey systems, the density-suppressed diffusion has already been commonly used in the modeling of other biological processes such as the chemotaxis ^[10,11], bacterial movement ^[12,13] and so on. Since the possible degeneracy caused by the density-suppressed diffusion brings considerable challenges for analysis, the studies of these biological models with density-suppressed diffusion have been increasingly attracting attentions and produced many interesting analytical results ^{[14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]} alongside rich numerical simulations demonstrating complex dynamics and patterns ^[29,30,31].

The purpose of this paper is to study the global dynamics of the defensive switching model with density-suppressed diffusion including global existence and asymptotic behavior of solutions as well as pattern formations. Specifically we consider the following problem

where $\Omega\subset{\mathbb{R}}^2$ is a bounded domain with smooth boundary. The parameters $\alpha_j\; (j = 1, 2), \; \beta_j\; (j = 1, 2, 3), \; \sigma$ and $K$ are all positive constants. Note that we have assumed $d_3 = 1$ without loss of generality and consider more general constant $\beta_3$ compared to Eq (1.4). We suppose the density-suppressed diffusion coefficients $d_j(w)$ satisfy the following conditions:

$(H_0)$ $d_j(w)\in C^2([0, \infty))$ with $d_j(w) > 0$ and $d_j'(w)\leq 0$ for all $w\geq 0$, $\ j = 1, 2$.

We further suppose the initial data satisfy

Then our first result on the global existence of solutions is stated in the following theorem.

Theorem 1.1 (Global boundedness). Let $\Omega\subset{\mathbb{R}}^2$ be a bounded domain with smooth boundary. Assume the assumption $(H_0)$ holds and the initial data $(u_0, v_0, w_0)$ satisfy Eq (1.6). Then there exists a uniquely determined triple $(u, v, w)$ of nonnegative functions which solves Eq (1.5) classically in $\Omega\times(0, \infty)$ and satisfies

where $M_0 > 0$ is a constant independent of $t$. Particularly, one has

Remark 1.2. Note that in works ^[5,6] considering the defensive switching model with prey-taxis, the smallness of prey-taxis coefficient was required to ensure the global boundedness of solutions. Here we obtain the existence of global classical solutions with uniform-in-time boundedness without any smallness assumptions on the parameters.

Next, we shall show the global stability of constant steady states. In fact it is straightforward to find that a positive constant steady state $(u_*, v_*, w_*)$ exists if and only if $K > \frac{\alpha_1}{\beta_1}+\frac{\alpha_2}{\beta_2}$, where $u_*, v_*, w_* > 0$, satisfying

which can be solved to obtain

Then we can show that $(u_*, v_*, w_*)$ is globally asymptotically stable under certain conditions, as stated in the following theorem.

Theorem 1.3 (Global stabilization). Suppose the conditions in Theorem 1.1 hold and let $(u, v, w)$ be the solution obtained in Theorem 1.1. If $K > \frac{\alpha_1}{\beta_1}+\frac{\alpha_2}{\beta_2}$ and

then the co-existence steady state $(u_*, v_*, w_*)$ is globally asymptotically stable.

The rest of this paper is organized as follows. In Section 2, we present the local existence theorem of solutions and establish some preliminary results. Then in Section 3, we prove Theorem 1.1. We prove Theorem 1.3 in Section 4 and explore the pattern formation in Section 5.

2.   Local existence and preliminaries

In the sequel, the integral $\int_\Omega f(x)dx$ and $\|f\|_{L^p(\Omega)}$ will be abbreviated as $\int_\Omega f$ and $\|f\|_{L^p}$, respectively. The generic constants $c_j$ or $K_j$ for $j = 1, 2, \cdots$, are independent of $t$ and will vary in the context. Below we present the local existence result of Eq (1.5), which can be proved in a similar way as ^[8] by applying Amann's theorem ^[36,37], and we omit the details for brevity.

Lemma 2.1 (Local existence). Let $\Omega\subset{\mathbb{R}}^2$ be a bounded domain with smooth boundary. Assume that the initial data $(u_0, v_0, w_0)$ satisfies (1.6) and suppose the hypothesis $(H_0)$ holds. Then there exists $T_{\max}\in(0, \infty]$ such that the system (1.5) admits a unique classical solution $(u, v, w)\in [C(\bar{\Omega}\times[0, T_{\max}))\cap C^{2, 1}(\bar{\Omega}\times(0, T_{\max}))]^3$ satisfying $u, v, w > 0$ for all $t > 0$. Moreover, if $T_{\max} < \infty$, then

Using the similar argument as in [38,Lemma 2.2], we obtain the global boundedness of $w$ immediately.

Lemma 2.2. Let $(u, v, w)$ be the classical solution of Eq (1.5) obtained in Lemma 2.1. Then it holds that

Furthermore, one has

Lemma 2.3. Let $(u, v, w)$ be the classical solution of (1.5) obtained in Lemma 2.1. Assume there is a constant $c_1 > 0$ such that

then one has

with

Proof. From the third equation of (1.5), we have

with $g(u, v, w): = w-\frac{\beta_3 uvw}{u+v}+\sigma w\left(1-\frac{w}{K}\right)$. Since $0 < w\leq K_*$ (see Eq (2.1)) and $u, v > 0$, we have

which, together with Eq (2.3), gives

With Eq (2.6), we use the results in [39,Lemma 1] to obtain Eq (2.4) with Eq (2.5) directly.

Now we will show some basic boundedness properties of the solution $(u, v, w)$ obtained in Lemma 2.1.

Lemma 2.4. Let $(u, v, w)$ be the classical solution of (1.5) obtained in Lemma 2.1. Then it holds that

and

where the constants $K_1, \; K_2 > 0$ are independent of $t$ and $\tau: = \min\{1, T_{\max}/2\}.$

Proof. Multiplying the first and second equations of (1.5) by $\beta_3$, and multiplying the third equation of (1.5) by $(\beta_1+\beta_2)$, then adding the resulting equations, one obtains

Thus, integrating Eq (2.9) with respect to $x$ over $\Omega$, we have

Using Young's inequality, one has

which, substituted into Eq (2.10), gives

where $\gamma_0: = \min\{\alpha_1, \alpha_2, 1\}$ and $c_1: = \frac{(\sigma+1)^2(\beta_1+\beta_2)K|\Omega|}{4\sigma}$. Then, applying the Grönwall's inequality to Eq (2.11) gives

which yields Eq (2.7).

Next, we shall show Eq (2.8) holds based on some ideas in ^[8]. Under the homogeneous Neumann boundary conditions, we define a shifted Laplacian operator $\mathcal{B} = -\Delta +\gamma_1$ with $\gamma_1: = \min\Big\{\frac{\alpha_1}{2d_1(0)}, \frac{\alpha_2}{2d_2(0)}\Big\} > 0$. Then $\mathcal{B}$ is a sectorial operator in $L^p(\Omega)$ for all $p\in (1, \infty)$ ^[32], and one can easily show that its inverse $\mathcal{B}^{-1}$ satisfies

and

for some constant $c_2 > 0$. Using the assumptions in ($H_0$) and Eq (2.1), we can derive that there exist two positive constants $\delta_1$ and $\delta_2$ independent of $t$ such for all $j = 1, 2$ one has

and

Using the definition of $\mathcal{B}$, we can rewrite Eq ($2.9$) as

Then applying the facts $\gamma_1: = \min\Big\{\frac{\alpha_1}{2d_1(0)}, \frac{\alpha_2}{2d_2(0)}\Big\} > 0$ and Eq (2.15) as well as $0 < w\leq K_*$ (see Eq (2.1)), we can derive that

where $c_3: = (\gamma_1+\sigma)(\beta_1+\beta_2)K_* > 0$.

We multiply Eq (2.16) by ${\mathcal{B}}^{-1}[\beta_3(u+ v)+(\beta_1+\beta_2)w]\geq0$ and integrate the result with respect to $x$ over $\Omega$ to obtain

which combined with Eq (2.15) enables us to find a positive constant $\gamma_2: = \min\{d_1(K_*), d_2(K_*), 1\}$ such that

On one hand, using Eq (2.12), the Hölder inequality and Young's inequality, one has

On the other hand, using Eq (2.13), we get

Defining $z(t): = \int_\Omega|{\mathcal{B}}^{-\frac{1}{2}}[\beta_3(u+ v)+(\beta_1+\beta_2)w]|^2$, and substituting Eqs (2.18) and (2.19) into Eq (2.17), we obtain

Then applying the Grönwall's inequality to Eq (2.20) alongside Eq (2.13), one has

Integrating Eq (2.20) over $(t, t+\tau)$ for $\tau = \min\{1, T_{\max}/2\}$ and applying Eq (2.21), we get

which gives Eq (2.8). Then we complete the proof of Lemma 2.4.

Lemma 2.5. Let $(u, v, w)$ be the classical solution of (1.5) and $\tau = \min\{1, T_{\max}/2\}$. Then it holds that

and

as well as

where the positive constants $K_3$, $K_4$ and $K_5$ are independent of $t.$

Proof. Multiplying the third equation of (1.5) by $w$, integrating the resulting equation by parts and using $0 < w\leq K_*$ in Eq (2.1), we obtain

Integrating Eq (2.25) over $(t, t+\tau)$ and using the facts $0 < \tau\leq 1$ and Eq (2.1), one has Eq (2.22) directly.

Next, we multiply the third equation of (1.5) by $-\Delta w$, use the integration by parts to the resulting equation, and apply Young's inequality as well as Eq (2.1) to get

which entails

where $c_1: = \frac{3\sigma^2K_*^4|\Omega|}{K^2} +3\sigma^2K_*^2|\Omega|$. Applying a Gagliardo-Nirenberg type inequality derived in [16,Lemma 2.5] and noting the fact $\|w\|_{L^2}\leq K_*|\Omega|^{\frac{1}{2}}$, we have

which substituted into Eq (2.26) gives

with $c_3: = 3\beta_3^2K_*^2+c_1+\frac{(2+c_2)c_2K_*^2|\Omega|}{2}$. For any $t\in(0, T_{\max})$, by Eq (2.22) there exists a $t_0: = t_0(t)\in((t-\tau)_+, t)$ such that $t_0\geq 0$ and

Moreover, from Eq (2.8), we know that

Multiplying Eq (2.27) by $e^{t}$ and integrating the results with respect to $t$ over $(t_0, t)$, and using the facts Eqs (2.28) and (2.29) and $t_0 < t\leq t_0+\tau\leq t_0+1$, one has

Hence Eq (2.23) holds. On the other hand, integrating Eq (2.27) along with Eqs (2.8) and (2.23), it holds that

which gives Eq (2.24). Then we complete the proof of Lemma 2.5.

3.   Global existence: Proof of Theorem 1.1

In this section, we shall show the existence of global classical solution based on some ideas in ^[8]. To this end, we first establish the a priori $L^2$-estimates of $u$ and $v$.

Lemma 3.1 ($L^2$-Boundedness). Suppose that the assumptions in Lemma 2.1 hold and let $(u, v, w)$ be the classical solution of (1.5) which is defined on its maximal existence time interval $[0, T_{\max})$. Then the solution of (1.5) satisfies

where the constant $K_6 > 0$ is independent of $t$.

Proof. Multiplying the first equation of (1.5) by $u$, integrating the results with respect to $x$ over $\Omega$ by parts and using the fact $u, v > 0$ for all $t\in(0, T_{\max})$, we have

which along with Eqs (2.14) and (2.15) indicates

Applying Young's inequality and the Hölder inequality, one gets

which updates Eq (3.2) as

Using the Gagliardo-Nirenberg inequality, the following inequality [16,Lemma 2.5]

and Eq (2.23), we obtain

where $c_3: = \frac{c_2^2K_4^2}{d_1(K_*)}$ and $c_4: = \frac{c_2^2K_4^4}{2d_1(K_*)}+\frac{d_1(K_*)}{2}+c_2K_4^2$.

Then substituting Eq (3.4) into Eq (3.3), we obtain

where $c_5: = 2\beta_1K_*+c_4+c_3$. For any $t\in(0, T_{\max})$, by Eq (2.8) there exists a $t_0: = t_0(t)\in((t-\tau)_+, t)$ such that $t_0\geq 0$ and

On the other hand, from Eq (2.24), we know that

Then integrating Eq (3.5) with respect to $t$ over $(t_0, t)$, and using the Eqs (3.6) and (3.7) and $t\leq t_0+\tau\leq t_0+1$, one has

Similarly, we can show that

Then the combination of Eqs (3.8) and (3.9) gives Eq (3.1). The proof of Lemma 3.1 is completed.

Lemma 3.2 ($L^{\infty}$-Boundedness). Suppose that the assumptions in Lemma 2.1 hold and let $(u, v, w)$ be the classical solution of (1.5) defined on its maximal existence time interval $[0, T_{\max})$. Then it holds that

where the constant $K_7 > 0$ is independent of $t.$

Proof. Multiplying the first equation of (1.5) by $u^{p-1}$ with $p > 2$ and integrating the resulting equation over $\Omega$ by parts, one has

Using Eqs (2.14) and (2.15) and $0 < \frac{vw}{u+v}\leq w\leq K_*$ which is satisfied by the fact $u, v > 0$ for all $t\in(0, T_{\max})$, from Eq (3.11) we get

which together with the fact

gives

Noting the fact $\|u\|_{L^2}\leq K_6$ in Eq (3.1), and using Lemma 2.3, we can find a constant $c_1 > 0$ such that

Then we can use the Hölder inequality, Gagliardo-Nirenberg inequality and Young's inequality as well as Eq (3.13) to obtain

where we have used the fact that

Moreover, using the Gagliardo-Nirenberg inequality and Eq (3.15) again, we obtain

Then substituting Eqs (3.14) and (3.16) into Eq (3.12) yields

where

Then applying the Grönwall's inequality to Eq (3.17), we get

Taking $p = 4$ in Eq (3.18), one enables $\|u(\cdot, t)\|_{L^4}\leq c_5,$ which together with Lemma 2.3 gives

where the constant $c_6 > 0$ is independent of $t$ and $p$. Then using the Moser iteration process ^[33], there exists a constant $c_7 > 0$ independent of $t$ such that

Similarly, we can find a constant $c_8 > 0$ such that

Then the combination of Eqs (3.19) and (3.20) gives Eq (3.10). The proof of Lemma 3.2 is completed.

Proof of Theorem 1.1. From Lemma 3.2, we know there exists a constant $c_1 > 0$ independent of $t$ such that

which along with the extensibility criterion in Lemma 2.1 and Lemma 2.2 gives Theorem 1.1.

4.   Global stabilization: Proof of Theorem 1.3

Under the condition $K > \frac{\alpha_1}{\beta_1}+\frac{\alpha_2}{\beta_2}$, the system (1.5) has only one positive constant steady states $(u_*, v_*, w_*)$ which is defined in Eq (1.8). In what follows, we shall show that the co-existence steady state $(u_*, v_*, w_*)$ is globally asymptotically stable based on the following energy functional:

where

Then we have the following results.

Lemma 4.1. Let $(u, v, w)$ be the solution of system (1.5) obtained in Theorem 1.1. If $K > \frac{\alpha_1}{\beta_1}+\frac{\alpha_2}{\beta_2}$ and

then the co-existence steady states $(u_*, v_*, w_*)$ is globally asymptotically stable.

Proof. Letting $\psi(f) = f-f_*\ln f$ and applying Taylor's expansion, for all positive $f_*$ and $f$, one has

where $\xi$ is between $f$ and $f_*$. Then choosing $f = u$ and $f_* = u_*$, we obtain from Eq (4.3) that there exists $\xi_1$ is between $u$ and $u_*$ such that

and $\mathcal{J}_u(t) = 0$ iff $u = u_*$. Using the similar way, one has $\mathcal{J}_v(t), \; \mathcal{J}_w(t)\geq0$ and "=" holds iff $v = v_*, \; w = w_*$. Hence according to the definition of $\mathcal{E}(t)$ in Eq (4.1), we obtain that $\mathcal{E}(u, v, w) > 0$ for all $(u, v, w)\ne (u_*, v_*, w_*)$ and moreover $\mathcal{E}(u_*, v_*, w_*) = 0$.

Next, we shall show that $\frac{d}{dt}\mathcal{E}(t)\leq 0$ under the conditions Eq (4.2) and $K > \frac{\alpha_1}{\beta_1}+\frac{\alpha_2}{\beta_2}$. In fact, using the equations in Eq (1.5), one has

which can be rewritten as

where

After some calculations, we can check that the matrix $B$ is positive definite provided

which, together with the fact $\frac{\beta_3u_*v_*}{u_*+v_*} = \sigma(1-\frac{w_*}{K})$ (see Eq (1.7)), gives

where the function

If $\|w_0\|_{L^\infty}\leq K$, one can easily check that Eq (4.5) holds provided Eq (4.2) is true. Now, we consider the cases $\|w_0\|_{L^\infty}\geq K$. In fact if Eq (4.2) holds, then there exists a small $\varepsilon_0 > 0$ such that

From the assumptions $(H_0)$, we know that $F(w)$ belongs to $C^1([0, \infty))$, which combined with the fact Eq (2.2), implies that

and thus for the chosen $\varepsilon_0$ in Eq (4.6), there exists a $T_* > 0$ such that

Then the combination of Eqs (4.6) and (4.7) gives

Hence the matrix $B$ is positive and there is a constant $c_1 > 0$ such that

On the other hand, using the fact $\frac{\alpha_1}{\beta_1} = \frac{ v_*w_*}{u_*+v_*}$ (see Eq (1.7)), one gets

which together with the fact

gives

Similarly, we can use the fact $\frac{\alpha_2}{\beta_2} = \frac{u_*w_*}{u_*+v_*}$ (see Eq (1.7)) to obtain

At last, using the fact $\sigma = \frac{\beta_3u_*v_*}{u_*+v_*}+\frac{\sigma w_*}{K}$ (see Eq (1.7)), we have

Combining Eqs (4.9) and (4.10) with Eq (4.11), we have

where

Then we can update Eq (4.12) as

Substituting Eqs (4.8) and (4.13) into Eq (4.4) gives for all $t\geq T_*$

which indicates $\frac{d\mathcal{E}(t)}{dt}\leq 0$ for all $u, v, w$. Moreover, one can check that if $\frac{d\mathcal{E}(t)}{dt} = 0$, then $\nabla u = \nabla v = 0$ and $w = w_*$, which implies that

with $\tilde{u}_*$ and $\tilde{v}_*$ are some positive constants. Since $(u, v, w) = (\tilde{u}_*, \tilde{v}_*, w_*)$ is a solution of (1.5), then one has

which gives

Thus $\frac{d\mathcal{E}(t)}{dt} = 0$ if and only if $(u, v, w) = (u_*, v_*, w_*)$. By using the LaSalle's invariant principle (cf. [41,Theorem 5.24]), we get that the co-existence steady state $(u_*, v_*, w_*)$ is globally asymptotically stable.

Proof of Theorem 1.3. Theorem 1.3 is a consequence of Lemma 4.1.

5.   Spatio-temporal patterns

In this section, we assume $K > \frac{\alpha_1}{\beta_1}+\frac{\alpha_2}{\beta_2}$ and will study the effect of density-suppressed dispersals on the dynamics of the system (1.5).

5.1.   Linear instability analysis

Denote $d_j: = d_j(w_*), d_j^{'}: = d_j^{'}(w_*)\; (j = 1, 2)$ and linearize (1.5) at $(u_*, v_*, w_*)$ to obtain

where $\mathcal{T}$ denotes the transpose and

with

Then by the standard linear stability principle, the linear stability of $(u_*, v_*, w_*)$ is determined by the eigenvalues of the matrix $(-\mu_kA+B)$ where the sequence $\{\mu_k\}_{k = 0}^{\infty}$ with $0 = \mu_0 < \mu_1\leq \mu_2\leq \mu_3\leq \dots$ denotes the eigenvalues of $-\Delta$ under the Neumann boundary condition. The characteristic equation of $(-\mu_kA+B)$ is

where

Hence, by direct calculations, one has

where $\gamma_1, \gamma_2, \gamma_3$ are positive constants and $A_1A_2-A_3 > 0$ are independent of $d'_j\; (j = 1, 2)$ and $\mu_k$ (see details in the Appendix for the precise definitions of $\gamma_i$ and $A_i$, $i = 1, 2, 3$), and

as well as

Then using the Routh-Hurwitz criterion, we have the results stated as below.

Proposition 5.1. The $(u_*, v_*, w_*)$ is linearly stable if one of the following conditions holds:

1) $d_1(w), d_2(w)$ are constants with $\sigma > 0$;

2) $d'_j(w) < 0$ with $\sigma > 0$ if $\alpha_1 = \alpha_2$ or $\sigma\geq\frac{K\beta_1\beta_2f(\beta)|\alpha_1-\alpha_2|}{(\alpha_1\beta_2+\alpha_2\beta_1)^2}$ if $\alpha_1\not = \alpha_2$, where

and $j = 1, 2$.

Proof. In the case that $d_j(w)\; (j = 1, 2)$ are constants, one has

On the other hand, considering the case $d'_j(w) < 0$. If $\alpha_1 = \alpha_2$, we get

As for $\alpha_1\not = \alpha_2$ with

we can calculate to get

All cases discussed above indicate $D_1(\mu_k)D_2(\mu_k, |d_j'|)-D_3(\mu_k, |d_j'|) > 0$ for all $k\in \mathbb{N}$. Then by the Routh-Hurwitz criterion ^[34], we know $(u_*, v_*, w_*)$ is linearly stable.

Hence, we are left to study the possible patterns bifurcating from the co-existence steady state $(u_*, v_*, w_*)$ in the following range of parameters:

$(H_1)$ $\alpha_1\ne \alpha_2$, $0 < \sigma < \frac{K\beta_1\beta_2f(\beta)|\alpha_1-\alpha_2|}{(\alpha_1\beta_2+\alpha_2\beta_1)^2}$ and $d_i'(w) < 0$, where $i = 1$ if $\alpha_1 < \alpha_2$ or $i = 2$ if $\alpha_1 > \alpha_2$ and $f(\beta)$ is defined in (5.4).

Proposition 5.2. Let $\sigma, K, \alpha_1, \alpha_2, \beta_1, \beta_2$ and $\beta_3$ be fixed positive parameters, and suppose the assumption $(H_1)$ holds. Then $(u_*, v_*, w_*)$ is linearly unstable provided that $|d_1^{'}(w_*)|$ is large enough if $\alpha_1 < \alpha_2$ (or $|d_2'(w_*)|$ is large enough if $\alpha_2 < \alpha_1$) and there is some $k_0\geq 1$ such that

Proof. Noting that $H(\mu_k, |d_j'|)$ can be rewritten as

Since $(H_1)$ holds, then for some fixed $\mu_{k_0} > 0$ satisfying Eq (5.5), we have

when $|d_1'(w_*)|$ is large enough if $\alpha_1 < \alpha_2$ (or $|d_2'(w_*)|$ is large enough if $\alpha_2 < \alpha_1$) and $(u_*, v_*, w_*)$ is linearly unstable by the Routh-Hurwitz criterion ^[40].

Remark 5.3. The above results imply that the density-suppressed diffusion can induce the instability of $(u_*, v_*, w_*)$ and trigger the pattern formation.

5.2.   Numerical simulations

In this subsection, we aim to give some numerical simulations to complement the previous analysis in Section 5.1. Without loss of generality, we assume $\alpha_1 < \alpha_2$, and define the set

as the bifurcation curve ^[35], where the linearized system (5.1) admits an eigenvalue with zero real part and the curve $\mathcal{H}$ is the graph of the function

where the constants $\gamma_j, A_j > 0\; (j = 1, 2, 3)$ are given in the Appendix. Then the characteristic Eq (4.3) has a pair of purely imaginary eigenvalues if $|d_1'(w_*)| = d_{\mathcal{H}}(\mu_{k_0})$ for some $k_0\in\mathbb{N}$ and $\mu_{k_0}$ satisfying Eq (5.5) and the Hopf bifurcation may emerge.

We first fix parameter values as follows:

Then by the assumption $(H_1)$, one has

Hence, we fix $\sigma = \frac{1}{10}$ in what follows without loss of generality. By Eq (1.8) we have

In addition, we choose the motility function $d_1(w)$ as

where the constant $D > 0$ and $d_1'(w) = -{D}{e^{-D(w-\frac{4}{3})}}$. Then Eq (5.5) is equivalent to

Using Theorem 1.3 and Proposition 5.2, the co-existence steady state $(\frac{4}{39}, \frac{4}{117}, \frac{4}{3})$ may lose its stability when the parameter $D > 0$ satisfies

and

with $0 < \eta < \frac{73}{240}$.

In our numerical simulation, we set the initial value $(u_0, v_0, w_0)$ as a small random perturbation of $(\frac{4}{39}, \frac{4}{117}, \frac{4}{3})$ and hence $K_* = 2$. Moreover, we perform numerical simulations in one dimension and set $\Omega = (0, 16\pi)$. Choosing $D > 0$ satisfies (5.9) and (5.10) with $\eta = (k/16)^2 < \frac{73}{240}$, which indicates allowable unstable modes are $k = 1, 2, 3, 4, 5, 6, 7, 8$. Next, we will consider the following two cases for $d_2(w)$.

Case 1: $d_2(w) = 1$. In this case, $d_2' = d_2'(w_*) = 0$ and according to the definition of $d_{\mathcal{H}}(\eta)$ in Eq (5.6), we can calculate to obtain that

and Eq (5.9) can be updated as

with $D > 4$. Then after some calculations, one gets $d_{\mathcal{H}}(\frac{1}{256})\approx 647.963$, $d_{\mathcal{H}}(\frac{4}{256})\approx223.026$, $d_{\mathcal{H}}(\frac{9}{256})\approx159.956$, $d_{\mathcal{H}}(\frac{16}{256})\approx162.600$, $d_{\mathcal{H}}(\frac{25}{256})\approx204.933$, $d_{\mathcal{H}}(\frac{36}{256})\approx 305.214$, $d_{\mathcal{H}}(\frac{49}{256})\approx 548.498$ and $d_{\mathcal{H}}(\frac{64}{256})\approx 1450.71$. Hence, $d_{\mathcal{H}}(\frac{9}{256})\approx159.956$ is the critical value for possible pattern formation. In our numerical simulations shown in Figure 1, we choose $D = 180$ and find the spatio-temporal pattern. In particular the time evolutionary profiles of solutions are oscillatory, which implies the bifurcation might be Hopf bifurcation.

Case 2: ${d_2(w) = {e^{-50(w-\frac{4}{3})}}}$. For this case, $d_2' = d_2'(w_*) = -50$ and hence

and Eq (5.9) can be simplified as

One also can calculate to get that $d_{\mathcal{H}}(\frac{1}{256})\approx 722.767$, $d_{\mathcal{H}}(\frac{4}{256})\approx306.961$, $d_{\mathcal{H}}(\frac{9}{256})\approx260.876$, $d_{\mathcal{H}}(\frac{16}{256})\approx291.910$, $d_{\mathcal{H}}(\frac{25}{256})\approx381.793$, $d_{\mathcal{H}}(\frac{36}{256})\approx 567.953$, $d_{\mathcal{H}}(\frac{49}{256})\approx 997.112$ and $d_{\mathcal{H}}(\frac{64}{256})\approx 2546.860$. Therefore, $d_{\mathcal{H}}(\frac{9}{256})\approx260.876$ is the critical value for possible pattern formation. In our numerical simulations, we choose $D = 268$ and find the spatio-temporal pattern shown in Figure 2. Again oscillatory time evolutionary profiles indicate the solution bifurcating from the co-existence steady state may undergo a Hopf bifurcation.

Acknowledgments

The authors are grateful to four referees for their valuable comments, which greatly improved the exposition of our paper. The research of H. Y. Jin was supported by the NSF of China (No. 11871226), Guangdong Basic and Applied Basic Research Foundation (No. 2020A1515010140, No. 2022B1515020032), Guangzhou Science and Technology Program (No. 202002030363).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

Below we show the details for the parameters $\gamma_1, \gamma_2, \gamma_3 > 0$ and $A_1A_2-A_3 > 0$ present in Eq (5.3). By direct calculations, we obtain

As for

which implies $\gamma_2 > 0$ and is independent of $d'_j\; (j = 1, 2)$ and $\mu_k$. Now, we calculate the other parameters:

and

Hence, a direct calculation gives