In recent years, more and more attention were given to the reaction-diffusion system of a predator-prey model with prey-taxis. For example, for the existence and uniqueness of weak solutions [1, 11], the global existence and uniqueness of classical solutions [2, 17, 18], pattern formation induced by the prey-taxis [12], global bifurcation for the predator-prey model with prey-taxis [13], boundedness or blow up in a chemotaxis system [14-16].

In this paper, we study the following reaction-diffusion system of a predator-prey model with Beddington-DeAngelis functional response and prey-taxis

$ \left\{\begin{array}{ll} u_t-d_1\Delta u+\nabla(u\chi(u)\nabla v)=-nu-hu^2+e\frac{muv}{au+bv+c}, ~~x \in \Omega, ~t>0, \\ v_t-d_2\Delta v=rv(1-\frac{v}{K})-\frac{muv}{au+bv+c}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x \in \Omega, ~t>0, \\ \frac{\partial u}{\partial \nu}=\frac{\partial v}{\partial \nu}=0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x \in \partial\Omega, t>0, \\ u(x, 0)=u_0(x)\geq 0, ~v(x, 0)=v_0(x)\geq 0, ~~~~~~~~~~~~~~~~~~~~~~x \in \Omega, \end{array}\right. $

(1.1)

where $\Omega$ is a bounded domain in $\mathbb{R}^N$ with smooth boundary $\partial\Omega$, initial data $u_0(x), v_0(x)\in C^{2+\alpha}(\bar{\Omega})$ compatible on $\partial \Omega$, and $\nu$ is the normal outer vector on $\partial \Omega$, $u$ and $v$ represent the densities of the predator and prey, respectively, $d_1, d_2, n, h, e, m, a, b, c, r, K$ are positive constants that stand for diffusion coefficients, death rate of $u$, intra-specific competition of $u$, conversion rate, consumption rate, predator interference, prey saturation constant, another saturation constant, intrinsic growth and carrying capacity of $v$, respectively.

The Beddington-DeAngelis model (1.1) is similar to the well known Holling type Ⅱ model with an extra term $au$ in the denominator which models the mutual interference among predators. In 2008, Ainseba et al. [1] proposed a Holling type Ⅱ model with prey-taxis and established the existence of weak solution by Schauder fixed point theorem and the uniqueness via duality technique. In 2010, Tao [2] gave the global existence and uniqueness of classical solution to Ainseba's model by contraction mapping principle together with $L^p$ estimates and Schauder estimates of parabolic equations. In 2015, He and Zheng [3] proved further more that the global classical solution is globally bounded.

There were also many works published for model (1.1). For the ODE system corresponding to (1.1), Cantrell and Cosner [4] presented some qualitative analysis of solutions from the view point of permanence and the existence of a global asymptotic stable positive equilibrium; Hwang [5] demonstrated that the local asymptotic stability of the positive equilibrium implies its global asymptotic stability. Chen and Wang [6] presented the qualitative analysis of system (1.1) from the view point of local asymptotic stability of the positive constant steady state and the existence and nonexistence of a nonconstant positive steady state. Haque [7] investigated the the influence of intra-specific competition among predators in the original Beddington-DeAngelis predator-prey model and offered a detailed mathematical analysis of the model. Yan and Zhang [8] studied model (1.1) without prey-taxis and obtained that the diffusion can destabilize the positive constant steady state of the system.

However, the emergency of the prey-taxis makes it more difficult to deal with the original problems. It is known that the global existence and bounedeness of solutions in (1.1) without prey-taxis can be easily obtained by using energy estimates and bootstrap arguments. In this paper, however, we will prove that the global classical solutions of (1.1) are moreover globally bounded by using the $L^p-L^q$ estimates for the Neumann heat semigroup and $L^p$ estimates with Moser's iteration of parabolic equations.

Throughout this paper, we assume that $\chi(u) \in C^1([0, +\infty)), ~\chi(u)\equiv 0$ for $u>u_m$, and $\chi'(u)$ is Lipschitz continuous, i.e., $|\chi'(u_1)-\chi'(u_2)|\leq L |u_1-u_2|$ for any $u_1, u_2\in [0, +\infty)$, where $u_m$ and $L$ are two positive constants. The assumption of $\chi$ is a regularity requirement for our qualitative analysis, and the assumption that $\chi(u)\equiv 0$ for $u>u_m$ has a clear biological interpretation [1]. Our main result is stated as follows.

Theorem 1  Under the assumptions for $\chi$ and initial data described above, the unique nonnegative classical solution of (1.1) is globally bounded.

The paper is organized as follows. We introduce some known results as preliminaries in Section 2. In Section 3, we give the proof of Theorem 1.

First we introduce the well-known classical $L^p-L^q$ estimates for the Neumann heat semigroup on bounded domains.

Lemma 1(see Lemma 1.3 in [9])  Suppose $(e^{t\Delta})_{t>0}$ is the Neumann heat semigroup in $\Omega$, and let $\lambda_1$ denote the first nonzero eigenvalue of $-\Delta$ in $\Omega$ under Neumann boundary conditions. Then there exist $C_1, C_2>0$ only depending on $\Omega$ such that the following estimates hold

(ⅰ) If $1\leq q\leq p\leq +\infty$, then

$ \|\nabla e^{t\Delta}\omega\|_{L^p(\Omega)}\leq C_1(1+t^{-\frac{1}{2}-\frac{N}{2}(\frac{1}{q}-\frac{1}{p})})e^{-\lambda_1t}\|\omega\|_{L^q(\Omega)}, ~t>0 $

for all $\omega\in L^q(\Omega)$;

(ⅱ) If $2\leq q\leq p\leq +\infty$, then

$ \|\nabla e^{t\Delta}\omega\|_{L^p(\Omega)}\leq C_2(1+t^{-\frac{N}{2}(\frac{1}{q}-\frac{1}{p})})e^{-\lambda_1t}\|\nabla\omega\|_{L^q(\Omega)}, ~t>0 $

for all $\omega\in W^{1, q}(\Omega)$.

One can obtain the boundedness of $v$ based on the comparison principle of ODEs.

Lemma 2  Let $(u, v)$ be a solution of (1.1). Then $u\geq 0$ and $0\leq v\leq K_1=\max\{K, \max_{\bar{\Omega}}v_0(x)\}$.

In this section, we give proof of Theorem 1, which is motivated by Tao and Winkler [10].

Proof of Theorem 1

Step 1  Boundedness of $\|u\|_{L^1(\Omega)}$.

Integrating the sum of the first equation and $e$ times of the second equation in (1.1) on $\Omega$ by parts, we have

$ \begin{array}{ll} \displaystyle\frac{d}{dt}\int_\Omega u+\frac{d}{dt}\int_\Omega ev &=\displaystyle -n\int_\Omega u-h\int_\Omega u^2+re\int_\Omega v-\frac{re}{K}\int_\Omega v^2\\ &\leq \displaystyle-(n+h)\int_\Omega u-r\int_\Omega ev+\left(reK+\frac{h}{4}\right)|\Omega|, \end{array} $

(3.1)

where we use that $\displaystyle\int_\Omega v\leq \frac{1}{2K}\displaystyle\int_\Omega v^2+\frac{K}{2}|\Omega|$ and $-\displaystyle\int_\Omega u^2\leq -\displaystyle\int_\Omega u+\frac{|\Omega|}{4}$. Define $y(t)=\displaystyle\int_\Omega u+\displaystyle\int_\Omega ev$ for $t>0$. Then $y'(t)+C_3y(t)\leq C_4$ for all $t>0$ by (3.1) with $C_3=\min\{n+h, r\}$ and $C_4=(reK+\frac{h}{4})|\Omega|$. This yields $y(t)\leq C_5=y(0)+\frac{C_4}{C_3}$ for all $t>0$ by the Gronwall inequality.

Step 2  Boundedness of $\|u\|_{L^p(\Omega)}$ with $p\geq2$.

Multiplying the first equation of (1.1) by $u^{p-1}$ and integrate on $\Omega$ by parts, combining Lemma 2.2, we have

$ \begin{align} & \ \ \ \ \ \ \frac{1}{p}\frac{d}{dt}\int\limits_{\Omega }{{{u}^{p}}}+{{d}_{1}}(p-1)\int\limits_{\Omega }{{{u}^{p-2}}}|\nabla u{{|}^{2}} \\ & =-n\int\limits_{\Omega }{{{u}^{p}}}+em\int\limits_{\Omega }{\frac{{{u}^{p}}v}{au+bv+c}}-h\int\limits_{\Omega }{{{u}^{p+1}}}+(p-1)\int\limits_{\Omega }{\chi }(u){{u}^{p-1}}\nabla u\cdot \nabla v \\ & \le \left( \frac{em{{K}_{1}}}{b{{K}_{1}}+c}-n \right)\int\limits_{\Omega }{{{u}^{p}}}+\frac{{{d}_{1}}(p-1)}{2}\int\limits_{\Omega }{{{u}^{p-2}}}|\nabla u{{|}^{2}}+\frac{p-1}{2{{d}_{1}}}\int\limits_{\Omega }{\chi }{{(u)}^{2}}{{u}^{p}}|\nabla v{{|}^{2}}. \\ \end{align} $

Together with $\chi(u)\leq M$ due to $\chi\in C^1$ and $\chi\equiv 0$ for $u\geq u_m$. This yields

$ \frac{1}{p}\frac{d}{dt}\int_\Omega u^p+\frac{d_1(p-1)}{2}\int_\Omega u^{p-2}|\nabla u|^2\leq (\frac{emK_1}{bK_1+c}-n)\int_\Omega u^p+ \frac{(p-1)M^2u_m^p}{2d_1}\int_\Omega|\nabla v|^2. $

(3.2)

Multiply the second equation of (1.1) by $-\Delta v$, and integrate on $\Omega$ by parts to get

$ \begin{array}{ll} \displaystyle\frac{d}{dt}\int_\Omega |\nabla v|^2+2d_2\int_\Omega |\Delta v|^2 &=\displaystyle2r\int_\Omega |\nabla v|^2-\frac{4r}{K}\int_\Omega v|\nabla v|^2+2m\int_\Omega \frac{uv}{au+bv+c}\Delta v\\ &\leq \displaystyle2r\int_\Omega|\nabla v|^2+\frac{2mK_1}{bK_1+c}\int_\Omega u|\Delta v|\\ &\leq \displaystyle2r\int_\Omega|\nabla v|^2+\frac{\epsilon}{2}\int_\Omega |\Delta v|^2+\frac{2m^2K_1^2}{\epsilon(bK_1+c)^2}\int_\Omega u^2 \end{array} $

by the Young inequality. Choosing $\epsilon=2d_2$, we have

$ \frac{d}{dt}\int_\Omega|\nabla v|^2+d_2\int_\Omega |\Delta v|^2\leq 2r \int_\Omega |\nabla v|^2+\frac{m^2K_1^2}{d_2(bK_1+c)^2}\int_\Omega u^2. $

(3.3)

Noting $\frac{d_1(p-1)}{2}\displaystyle\int_\Omega u^{p-2}|\nabla u|^2=\frac{2d_1(p-1)}{p^2}\int_\Omega|\nabla u^{\frac{p}{2}}|^2, ~p>2$. From (3.2) and (3.3), combining Young's inequality, we obtain

$ \frac{1}{p}\frac{d}{dt}\int_\Omega u^p+\frac{d}{dt}\int_\Omega|\nabla v|^2+\frac{2d_1(p-1)}{p^2}\int_\Omega|\nabla u^{\frac{p}{2}}|^2+d_2\int_\Omega |\Delta v|^2~~~~~~\\ \le \left(\frac{emK_1}{bK_1+c}+\epsilon_1-n\right)\int_\Omega u^p+\left(2r+\frac{(p-1)M^2u_m^p}{2d_1}\right)\int_\Omega|\nabla v|^2+C_6 $

(3.4)

with $C_6>0$ depending on $\epsilon_1$. By the Sobolev interpolation inequality and Lemma 2.2, we have for any $\epsilon_2>0$ that

$ \int_\Omega|\nabla v|^2\leq \epsilon_2\int_\Omega |\Delta v|^2+C_7\int_\Omega|v|^2\leq \epsilon_2\int_\Omega |\Delta v|^2+C_8, $

(3.5)

where $C_7, C_8>0$ depending on $\epsilon_2$. Applying the Gagliardo-Nirenberg inequality yields

$ \int_\Omega u^p=\int_\Omega|u^{\frac{p}{2}}|^2\leq C_9\|\nabla u^{\frac{p}{2}}\|_2^{2\theta}\|u^{\frac{p}{2}}\|_{\frac{2}{p}}^{2(1-\theta)} +C_9\|u^{\frac{p}{2}}\|_{\frac{2}{p}}^2 $

with $0<\theta=\frac{Np-N}{Np-N+2}<1$ and $C_9>0$. By Young's inequality,

$ \int_\Omega u^p\leq \epsilon_3\|\nabla u^{\frac{p}{2}}\|_2^{2} +C_{10}\|u^{\frac{p}{2}}\|_{\frac{2}{p}}^2=\epsilon_3\|\nabla u^{\frac{p}{2}}\|_2^{2}+C_{10}\|u\|_1^p $

for any $\epsilon_3>0$, with $C_{10}>0$ depending on $\epsilon_3$. By Step 1, we know $\|u\|_1\leq C_5$. So

$ \int_\Omega u^p\leq \epsilon_3\|\nabla u^{\frac{p}{2}}\|_2^{2}+C_{11} $

(3.6)

with $C_{11}>0$. Now fix $\epsilon_2, \epsilon_3$ such that $(2r+\frac{(p-1)M^2u_m^p}{2d_1})\epsilon_2=\frac{d_2}{2}$ and $(\frac{emK_1}{bK_1+c}+\epsilon_1)\epsilon_3=\frac{2d_1(p-1)}{p^2}$. From (3.4)-(3.6), we have

$ \frac{1}{p}\frac{d}{dt}\int_\Omega u^p+\frac{d}{dt}\int_\Omega|\nabla v|^2\leq -n\int_\Omega u^p-\left(2r+\frac{(p-1)M^2u_m^p}{2d_1}\right)\int_\Omega|\nabla v|^2+C_{12} $

with $C_{12}>0$. Define $z(t)=\frac{1}{p}\displaystyle\int_\Omega u^p+\int_\Omega |\nabla v|^2, ~t>0$. Then $z'(t)\leq -C_{13}z(t)+C_{12}$ for all $t>0$ with $C_{13}=\min\{pn, 2r+\frac{(p-1)M^2u_m^p}{2d_1}\}$. By the Gronwall inequality, we have $z(t)\leq C_{14}=z(0)+\frac{C_{12}}{C_{13}}$ for all $t>0$.

Step 3  Boundedness of $\|\nabla v\|_{L^p(\Omega)}$ with $p\geq2$.

Define $f(u, v)=rv(1-\frac{v}{K})-\frac{muv}{au+bv+c}$. From Lemma 2.2 and Step 2, there exists a constant $C_{15}$ such that $\sup\limits_{t>0}\|f\|_{L^p(\Omega)}\leq C_{15}<+\infty.$ The variation of constants formula for $v$ yields

$ v(\cdot, t)=e^{d_2t\Delta}v_0+\int_0^te^{d_2(t-s)\Delta}f(u(s), v(s))ds, ~t>0. $

By Lemma 2.1, we conclude that

$ \begin{array}{ll} \displaystyle\|\nabla v\|_{L^p(\Omega)} &\leq\displaystyle \|\nabla e^{d_2t\Delta}v_0\|_{L^p(\Omega)}+\int_0^t\|\nabla e^{d_2(t-s)\Delta}f(u(s), v(s))\|_{L^p(\Omega)}ds\\ &\leq \displaystyle 2C_2e^{-\lambda_1t}\|\nabla v_0\|_{L^p(\Omega)}+\\ & C_1\int_0^t(1+d_2^{-\frac{1}{2}}(t-s)^{-\frac{1}{2}})e^{-\lambda_1d_2(t-s)}\|f(s)\|_{L^p(\Omega)}ds\\ &\leq \displaystyle 2C_2e^{-\lambda_1t}\|\nabla v_0\|_{L^p(\Omega)}+C_1C_{15}\int_0^t(1+d_2^{-\frac{1}{2}}s^{-\frac{1}{2}})e^{-\lambda_1d_2s}ds\\ &\leq \displaystyle 2C_2\|\nabla v_0\|_{L^p(\Omega)}+C_1C_{15}\left(\frac{1}{\lambda_1d_2}+d_2^{-\frac{1}{2}}(2+\frac{1}{\lambda_1d_2})\right)\\ & ~\textrm{for all} ~t>0. \end{array} $

Step 4  Global boundedness.

On the basis of Steps 2 and 3, using Lemma A.1 in [10], we can obtain the global boundedness of solutions to (1.1) by the standard Moser iterative technique.

Remark  we used the assumption that $\chi'(u)$ is Lipschitz continuous, which is a necessary condition for existence of the global solutions (see [2]).

On the other hand, the intra-specific competition term $hu^2$ makes our estimates easier, which is a "good" term. This also coincides with Haque's [7] result that competition among the predator population is beneficial for both populations co-existence.