Pattern formation in a ratio-dependent predator-prey model with cross diffusion

1.   Introduction

This paper investigates the following ratio-dependent predator-prey model with cross-diffusion

where the habitat $\Omega$ is a bounded domain in ${\mathbb{R}}^n$ with smooth boundary $\partial\Omega$. The functions $u(x, t)$ and $v(x, t)$ denote the densities of prey and predator at space location $x$ and time $t$, respectively. The coefficients $b, c, k, l$ and $m$ are all positive constants; $\beta$ and $\rho$ are non-negative constants. The cross diffusion term $\Delta\left(\frac{\beta v}{1+\rho u}\right)$ indicates that the population pressure of predator species diminishes where the density of prey species is high. The prey and predator interact according to the ratio-dependent functional response which is given by the term $\frac{uv}{mu+v}$ when the predators hunt seriously.

As is known, kinds of biological populations correspond to various peculiar and interesting ecological phenomena. Among them, most of such phenomena appearing in ecology can be described by the reaction-diffusion models. These models typically consist of the diffusion term and the reaction term. Along with the innovative work ^[1] by Shigesada et al., many mathematicians and ecologists have been dedicated to investigate the population models with cross-diffusion from different mathematical perspectives. Especially for the bounded domain, the global existence in time, the existence and stability of positive steady states have been extensively studied in recent two decades (see ^{[2,3,4,5,6,7,8,9,10]} and the references therein). For the prey-predator relationship as Eq (1.1), the cross-diffusion term of this fractional type is consistent with the typical pattern of observed interaction between prey and predator. Biomathematical models with such cross-diffusion terms have attracted tremendous attention of scholars (see ^{[11,12,13,14]} and the references therein). Considering the saturation of predators for preys, the Holling-Ⅱ functional response $\frac{uv}{ mu+1}$ is widely used to describe the interactions of preys and predators. In ^[12], Wang and Li proved the global bifurcation branch of the positive solution for the prey-predator system with Holling-Ⅱ functional response and cross diffusion. They also investigated the limiting behavior of the steady states when the cross-diffusion approaches infinity. The model was further studied in ^[13], where the qualitative properties of positive solutions are obtained by applying the Leray-Schauder degree theory and bifurcation argument. Under Neumann boundary condition, Cao et al. ^[14] have considered the Turing instability in the prey-predator system with cross diffusion when the diffusion coefficients of prey are negative. However, in some situations, especially when predators have to search, share or compete for food, a more suitable prey-predator interaction should be the ratio-dependent functional response $\frac{uv}{mu+v}$. For a more detailed biological description of the ratio-dependent functional response, one can refer to ^{[15,16,17,18,19]} and the references therein. As far as we know, there are few studies concerned on the model with both cross-diffusion and ratio-dependent functional response. Recently in virtue of a priori estimates and bifurcation theory, Kumari and Mohan ^[20] proved the existence and the global bifurcation set of positive steady states to the Eq (1.1) with the bifurcation parameter $l$.

In the current work, we focus on the existence and stability of the nonconstant positive steady state for the ratio-dependent prey-predator model (1.1) with cross-diffusion, which are closely related to the pattern formation between predator and prey in ecology. In Section 2, it can be proved that there exists the nonconstant steady states of the nonlinear predator-prey model (1.1) bifurcating from the semi-trivial solution with the bifurcation parameter $\rho$. In Section 3, we study the stability of such bifurcating steady states by virtue of the spectral analysis. This paper ends with a discussion section containing both numerical simulations and ecological interpretations of our results.

2.   Existence of the positive steady states bifurcating from the semi-trivial steady state

At the beginning of this section, we introduce some symbols and known results. For $p > n$, define the Banach spaces

Let

then

By substituting Eq (2.2) into Eq (1.1), one can see that the $v-$ equation in Eq (1.1) becomes

where $g_u = \frac{V\rho\beta}{(1+\rho u+\beta)^2}$ and $g_{_V} = \frac{1+\rho u}{1+\rho u+\beta}$. Hence $(u, V)$ satisfies the following evolutional problem,

The monotonicity of $g = g(u, V, \rho)$ in $u$ and $V$ guarantees that a one-to-one correspondence is formed by Eq (2.1) between each solution $(u, v)$ of Eq (1.1) and each solution $(u, V)$ of Eq (2.3). Therefore, the positive nonconstant steady state $(u, v)$ of Eq (1.1) exists if and only if the positive nonconstant steady state $(u, V)$ of Eq (2.3) does. Obviously, the steady state $(u(x, t), V(x, t))$ of Eq (2.3) satisfies the following system,

For each $h(x)\in C^1(\bar{\Omega})$, let $\lambda_1(h)$ be the smallest eigenvalue of the following elliptic problem

then it is known that $\lambda_1(h)$ is simple and real. Moreover, $\lambda_1(h)$ is strictly increasing in $h(x)$. If $h(x)\equiv0$, then $\lambda_1(0)$ can be simply denoted by $\lambda_1$ with the corresponding positive eigenfunction $\Phi_{1}(x)$ normalized by ${||\Phi_{1}||_{L^2(\Omega)} = 1}$. It is noted that when the constant $a > \lambda_1$, the following Dirichlet problem of elliptic equation

has a unique positive solution denoted by $u(x) = \theta_{a}(x)$. Hence the stationary problem (2.4) has a semi-trivial solution $(u, V) = (\theta_k(x), 0)$ if $k > \lambda_1$.

For convenience of our later use, we will show the nonexistence results and some priori estimates of the positive solution for Eq (2.4), which have been proved in ^[20].

Lemma 2.1. (i) If $k\leqslant \lambda_1, l\leqslant\lambda_1- \frac{c}{m}$, then (2.4) has no positive solution.

(ii) Assume that $(u, V)$ and $(u, v)$ are positive steady states of Eqs (2.3) and (1.1), respectively, then for each $x\in\Omega$, it holds that

In the following, we try to obtain the positive solutions of Eq (2.4) emanating from the semi-trivial solution $(\theta_k(x), 0)$ by regarding $\rho$ as the bifurcation parameter.

Define the operator $F:\mathbf{X}\times\mathbb{R}^{+}\rightarrow \mathbf{Y}$ by

Denote the Fr$\acute{e}$chet derivative of $F$ with respect to $u$ and $V$ by $D_{(u, V)}F(u, V, \rho)$ as follows,

with $g_{_V} = g_{_V}(\theta_k, 0, \rho)$. For the sake of applying the local bifurcation theorem proposed by Crandall and Rabinowtiz ^[21], we firstly show that the kernel space of $D_{(u, V)}F(\theta_k, 0, \rho)$ is nontrivial. If $(\phi(x), \psi(x))\in{\rm Ker}D_{(u, V)}F(\theta_k, 0, \rho)$, it follows from Eq (2.7) that

To obtain the nontrivial solution of the system (2.8), we need to find some $\rho$ such that $S(k, l, \rho) = 0$ holds for any $k > \lambda_1$ and $l > \lambda_1-\frac{c}{m}$, where $S(k, l, \rho)$ is defined below.

Lemma 2.2. Assume that $\lambda_1(h)$ be the smallest eigenvalue of (2.5) and $\lambda_1 = \lambda_1(0)$. Let

and

then the set $\Gamma$ can be expressed by

where $\rho = \rho_1(k, l)$ is a positive continuous function.

Proof. Note that

and the mapping $q(x)\rightarrow \lambda_1(q):C(\bar{\Omega})\rightarrow {\mathbb{R}}$ is smooth and strictly increasing. Thus, it follows that

With the similar proof of Lemma 6 in ^[20], we can prove that there exists a continuous function $l(k, \rho)$ satisfying $S(k, l, \rho) = 0$ and $\lambda_1-\frac{c}{m} < l < (1+\beta)\lambda_1-\frac{c}{m}$ for any $k > \lambda_1, \; \rho > 0$. Then for any fixed $k_0 > \lambda_1$ and $l_0\in(\lambda_1-\frac{c}{m}, (1+\beta)\lambda_1-\frac{c}{m})$, there exists a unique $\rho_0 > 0$ such that $S(k_0, l_0, \rho_0) = 0.$ In virtue of the Implicit Function Theorem and Eq (2.9), it follows that there exists a small $\varepsilon > 0$ and a unique function $\rho = \rho_1(k, l)$ with $(k, l)\in (k_0-\varepsilon, k_0+\varepsilon)\times(l_0-\varepsilon, l_0+\varepsilon)$ such that $S(k, l, \rho_1(k, l)) = 0$. As $k_0$ and $l_0$ are arbitrary, it can be obtained that there exists a smooth function $\rho = \rho_1(k, l)$ such that $S(k, l, \rho_1(k, l)) = 0$ for $(k, l)\in (\lambda_1, +\infty)\times(\lambda_1-\frac{c}{m}, (1+\beta)\lambda_1-\frac{c}{m})$, which completes the proof.

For ease of notation, we simplify $\rho_1(k, l)$ as $\rho_1$. From Lemma 2.2, there exists a positive function denoted by $\psi^*(x)$ which solves the following eigenvalue problem

with $\int_\Omega(\psi^*)^2 \mathrm{d}x = 1$. Now we give the main result in this section.

Theorem 2.3. Assume $k > \lambda_1$ and $\lambda_1 < l+\frac{c}{m} < (1+\beta)\lambda_1$ hold. Then there exist a small $\delta > 0$ and a smooth function $\rho(s)\in C(0, \delta]$ with $\rho(0) = \rho_1$ such that the stationary problem (2.4) has a positive nonconstant solution $(u(x, s), V(x, s))$ bifurcating from $(\theta_k(x), 0)$, which satisfies the following expression:

where

and $(u_1(\cdot, s), V_1(\cdot, s))\in C[(0, \delta), X]$ satisfying $u_1(x, 0) = V_1(x, 0) = 0$.

Proof. By Lemma 2.2 and Eq (2.10), it can be deduced that

Under the Dirichlet boundary condition, the operator $-\Delta-k+2\theta_k$ is proved to be invertible in ^[22]. Moreover, when $\phi$ is positive, so is $(-\Delta-k+2\theta_k)^{-1}\phi$. Next, we will certify codim[$\rm{Range}$ $D_{(u, V)}F(\theta_k, 0, \rho_1)] = 1$. Suppose that $(\xi, \eta)\in \rm{Range}$ $D_{(u, V)}F(\theta_k, 0, \rho_1)$, then there exists $(\phi, \psi)\in X$ such that

By Eq (2.10) and applying Fredholm alternative theorem, there exists a solution $\tilde{\psi}$ for the second equation of Eq (2.12) if and only if $\int\eta\psi^* \mathrm{d}x = 0$. Substituting $\tilde{\psi}$ into the first equation of Eq (2.12), it has a unique solution $\tilde{\phi}$ due to the invertibility of the operator $-\Delta-k+2\theta_k$. Then, we can make a conclusion that codim[$\rm{Range}$ $D_{(u, V)}F(\theta_k, 0, \rho_1)] = 1$.

For the purpose of applying the local bifurcation argument at $(\phi, \psi, \rho) = (\theta_k, 0, \rho_1)$, we need to certify that

where

By contradiction, suppose that

Then we can find some $(\hat{\phi}, \hat{\psi})\in X$ such that

Multiplying Eq (2.15) by $\psi^*$ and integrating by parts, then by using Eq (2.10) we can obtain that

which is impossible due to the fact that the right-hand side of Eq (2.16) is positive. Thus, this completes the proof as we have verified the transversality condition.

3.   Stability of the bifurcating steady states for the Eq (2.3)

Thanks to some abstract theories of stability on the basis of analytic semigroup theory (see ^[23]), the stability of the steady state for the Eq (2.3) can be obtained by proving the spectral stability of the steady state in $\mathbf{W}^{1, p}(\Omega)\times\mathbf{W}^{1, p}(\Omega)$ for $p > n$. In this section, we mainly investigate the distribution of spectrum. Firstly, we study the bifurcation direction which is useful for subsequently analyzing the stability of the positive steady state $(u(x, s), V(x, s))$.

In the proof of Theorem 2.3, we see that

Define the adjoint operator $D_{(u, V)}F^*(\theta_k, 0, \rho_1)$ of $D_{(u, V)}F(\theta_k, 0, \rho_1)$ as follows:

Since $(\Delta+k-2\theta_k)$ is invertible, it follows that

Hence by the spectrum decomposition theorem, Eqs (2.11) and (3.1) imply that $X$ and $Y$ have the following direct decomposition

with $\oplus$ a direct sum in $Y$.

Lemma 3.1. For each fixed $k > \lambda_1$, the bifurcating direction satisfies

Proof. The derivative $\frac{d\rho(s)}{ds}\mid_{s = 0}$ can be expressed in the above form by using the bifurcation formula I.6.3 in ^[24]. Hereafter, it remains to evaluate $\frac{d\rho(s)}{ds}\mid_{s = 0} > 0$. From Eq (2.6), it can be calculated that

Thus,

By the facts that $\phi^* < 0$ and $\psi^* > 0$, we have

Meanwhile, Eq (2.14) implies that

which shows that $\rho'(0) > 0$.

For the rest of this section, we will study the stability of the nonconstant steady state $(u(x, s), V(x, s), \rho(s))$ near $(\theta_k, 0, \rho_1)$. Linearizing the Eq (2.3) at the steady state $(u(x, s), V(x, s))$, we can obtain that

The corresponding eigenvalue problem of the linearized Eq (3.2) with the eigenvalue $\sigma$ is as follows,

Rewrite Eq (3.3) as follows

Introduce an operator $H: X\times{\mathbb{R}}^{+}\rightarrow Y$ by

In virtue of Eqs (3.4) and (3.5), the Eq (3.3) can be transformed into the following system

With the aim of studying the spectral stability of the steady state, we need to certify that there is no eigenvalue with nonnegative real part of the linearized operator $H$.

Theorem 3.2. Suppose that $k > \lambda_1$, $\lambda_1 < l+ \frac{c}{m} < (1+\beta)\lambda_1$ and $\rho-\rho_1 > 0$ small enough. When $0 < s\leqslant\bar{\delta}$ with $\bar{\delta} > 0$ small enough, then the bifurcating steady state $(u(x, s), V(x, s))$ is locally asymptotically stable.

Proof. According to Eq (3.5), we can deduce that

From Eqs (2.10) and (2.11), it follows that

The above equation shows that $0$ is an eigenvalue of the operator $H(\theta_k, 0, \rho_1)$ with the corresponding eigenfunction $(\phi^*, \psi^*)^{\top}$.

Then, we prove that $0$ is the principal eigenvalue of $H(\theta_k, 0, \rho_1)$. By contradiction, we assume that the operator $H(\theta_k, 0, \rho_1)$ has a positive eigenvalue $\sigma_1$ with the corresponding eigenfunction $(\phi_1, \psi_1)\in X$ which satisfies

i.e.,

Assume $\psi_1 = 0$. Then Eq (3.10) implies

which is a contradiction due to the fact that $(-\Delta-k+2\theta_k)^{-1}$ is a positive operator, and hence $\psi_1\neq0$. From the second equation of Eq (3.10), it can be shown that

However by Lemma 2.2 and Eq (2.10), the principal eigenvalue of the operator $\Delta +\left(l+\frac{c}{m}\right)\frac{1+\rho\theta_k}{1+\rho\theta_k+\beta}$ is zero, which contradicts with $\sigma_1 > 0$. Thus all eigenvalues of $H(\theta_k, 0, \rho_1)$ except zero are negative.

In the following, we investigate the distribution of spectrum for the linear operator $H(u(x, s), V(x, s), \rho(s))$ by means of perturbation arguments (see Corollary 1.13 in ^[25]). For $0 < s < \bar{\delta}$ with $\bar{\delta}$ small enough, then the linear operator $H(u(x, s), V(x, s), \rho(s))$ admits an eigenvalue $\sigma(s)$ perturbed from zero and the corresponding function $(\phi_1(x, s), \psi_1(x, s))\in X_R$, where

with $\sigma(0) = 0$ and $\phi_1(x, 0) = \psi_1(x, 0) = 0$. Analogously, there exists an eigenvalue $\mu(\rho)$ perturbed from zero with the corresponding continuous differential functions $(\phi_1(x, \rho), \psi_1(x, \rho))\in X_{R}$ satisfying

with $\mu(\rho_1) = 0$ and $\phi_1(x, \rho_1) = \psi_1(x, \rho_1) = 0$. Differentiating Eq (3.12) with respect to $\rho$ at $\rho = \rho_1$, it yields that

where $\mu'$ denotes $\frac{d\mu}{d\rho}$. From Eq (3.1), we have

It follows from Eq (3.7) that the left side of Eq (3.13) is equal to

Finally, combining the above results with Eq (3.7), it indicates that when $\rho < \rho_1$, the semi-trivial steady state $(\theta_k, 0)$ is stable and unstable when $\rho\geqslant \rho_1$. According to Theorem 1.16 in ^[25], we have

where $\dot{\sigma}(s) = \frac{d\sigma}{ds}$. Moreover by Lemma 3.1, it can be deduced that $\dot{\sigma}(0) < 0$. This implies that $\sigma(s) < 0$ for $0 < s < \bar{\delta}$ when $\bar{\delta}$ is sufficiently small, which also shows that the bifurcating steady state $(u(x, s), V(x, s))$ is locally asymptotically stable.

4.   Discussion

The study of spatial patterns in the distribution of organisms is an important issue in ecology. Over the past decades, a large number of papers have been published to gain a better understanding of classical prey-dependent models and ratio-dependent predator-prey models. Extinction of one or both populations in predator-prey systems has occupied the most of the predator-prey literature. Therefore, it makes sense to study that how diffusion affects the stability of the predator-prey coexistence equilibrium, i.e., the spatial patterns, in the ratio-dependent model.

The theoretical results in Sections 2 and 3 have shown that for some region in the parameter space of the ratio-dependent model, coexistence phenomenon can appear. In this section, we perform numerical simulations to demonstrate the spatial-temporal behaviors of the ratio-dependent predator-prey model (1.1) in one dimensional case. For the sake of convenience, we take the domain $\Omega = (0, 1)$. Besides that, we choose the mesh size of space and time to be $\Delta x = 0.01$ and $\Delta t = 0.01$. To study the effect of cross-diffusion on the dynamics of Eq (1.1), we fix $k = 30 > \lambda_1 = \pi^2, \; l = 10, \; m = 1.1, \; b = 10, \; c = 3, \; \rho = 2, \; \beta = 1$ and choose the initial data to be small perturbations of $(\theta_k(x), \; 0)$ in simulations. The numerical simulations support our theoretical results on the existence and stability of the bifurcating solutions (see Figure 1). Biologically, the results imply that the species $u, v$ can coexist if predators hunt and favor to stay in the region where density of prey species is high. We hope that the observations in this paper will help experimental ecologists to carry out some experimental setups that will ensure biodiversity.

Acknowledgments

We are very grateful to the anonymous referees for careful reading and valuable comments, which led to an improvement of the original manuscript. The first author was supported by National Natural Science Foundation of China (No. 61902238, No. 12201426) and Shanghai Science and Technology Innovation Plan Project under Grant (No. 22692194600). The second author was supported by National Natural Science Foundation of China (No. 12201426).

Conflict of interest

The authors declare there is no conflicts of interest.