Global dynamics of a modified Leslie-Gower predator-prey model with Beddington-DeAngelis functional response and prey-taxis

1.   Introduction

In the ecosystem, the interaction of predator and prey is well-known and essential. It has become a hot topic to study the dynamic behavior of a predator-prey model. A famous predator-prey system established by Leslie and Gower in 1960 is of the form

where $u(t)$ and $v(t)$ are the population density of prey and predator at time $t$, respectively. Here $\lambda$ and $a$ are the intrinsic growth rate and the strength of competition among individuals of preys. The term $\frac{ev}{u}$ is called the Leslie-Gower term which means the loss in the predator population only due to rarity of its favorite food, where the parameter $e$ is a measure of the amount of food provided by the prey transformed into the birth of predator. The environment carrying capacity for the predator $\frac{u}{e}$ is not constant but proportional to the number of the prey. And they found that the predator could switch over to other preys even though its growth would be limited by the shortage of its favorite food. Hence, a positive constant d should be added into the denominator of the Leslie-Gower term, which is called a modified Leslie-Gower term $\frac{ev}{u+d}$.

The functional response function $\phi(u, v)$ represents the consumption of prey. It is particularly significant to select an appropriate response function to describe the relationship between the predator and the prey. As is known to all, the functional response can be classified into two types: prey-dependent and predator-dependent. The earliest functional response function $(\phi(u) = u)$ was proposed by Lotka and Volterra ^[1]. In the following research process, many scholars proposed several different response functions according to different predators and preys, among which Holling Ⅱ type $(\phi(u) = \frac{qu}{\alpha+bu})$ has been studied by a large number of researchers ^[2].

Recent accumulating evidence shows that predator-dependent is more realistic than prey-dependent in depicting the consuming of the prey. The classic example is Beddington-DeAngelis (abbreviated as B-D) functional response proposed by Beddington ^[3] and DeAngelis ^[4], which has the following form

where $q$ is the consumption rate; $\alpha, b, c$ mean the saturation constant, the saturation constant for an alternative prey and the predator interference, respectively. Compared with Holling-Ⅱ functional response, B-D functional response has an extra term $cv$ in the denominator modeling mutual interference among predators, which can exhibit more plentiful, more complicated and more acceptable dynamics ^[5,6,7]. In ^[8], Yu considered B-D functional response into the system (1.1). For this case, (1.1) becomes

with an initial condition $u_{0}(t) = u_0, v_{0}(t) = v_0$. He discussed the structure of nonnegative equilibria to (1.2) and their local stability. In addition, he applied the fluctuation lemma and Lyapunov direct method to get the global asymptotic stability of a positive equilibrium.

In the real world, the distribution of population density in a fixed bounded domain is inhomogeneous which makes that the population in high density area will spread to its low density area. Hence, establishing and studying various reaction-diffusion systems have been an effective way for researchers to further explore and predict biological evolution ^[9]. Through choosing appropriate scale transformation, (1.2) can be rewritten as

In fact, on top of random diffusion of the predator and the prey, it has been recognized that the spatial-temporal variations of the predator moves along the gradient direction of the prey. This kind of movement which is not random but directed is called prey-taxis, various types of predator-prey models with prey-taxis have received great attention among mathematical ecologist ^[10,11]. In detail, Wu, Shi and Wu ^[12] in 2016 was the first one that established the global boundedness for such model in higher dimension space with small $\chi > 0$, specific models are as follows:

Later on, Jin and Wang ^[13] analyzed the global boundedness for such model in two-dimensional domain and also prove the global stability of prey-taxis system. Wang and Wang ^[14] concerned the reaction-diffusion systems modeling the population dynamics of two predators and one prey with nonlinear prey-taxis in 2018. Wu and Ni ^[15] in 2021 proved the global existence and boundedness of solutions of a diffusive prey-predator model with prey-taxis and trophic interactions of three levels. And prey-taxis and predator-taxis also have an essential impact on pattern formation ^[16,17,18].

Coupled with the factors mentioned above, a modified Leslie-Gower predator-prey model with Beddington-DeAngelis functional response and prey-taxis can be formulated as:

Here $u(x, t)$ and $v(x, t)$ represent the densities of prey and predator at place $x$ and time $t$, $\Omega$ is a bounded domain in $\mathbb{R}^n$ with smooth boundary $\partial\Omega$ and $n$ is the outward unit normal vector on $\partial\Omega$. $\chi$ is called prey-taxis coefficient, and prey-taxis is called attractive (repulsive) if $\chi > 0$ $(\chi < 0)$. The parameters $d_1$ and $d_2$ are the diffusion rates of the prey and predator respectively. And we assume that all parameters are positive and have the same meaning as above. Our first main result is the following:

Theorem 1.1. Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with smooth boundary. Suppose that $(u_0, v_0)\in [W^{1, \infty}]^2$ with$u_0, v_0\geq 0(\not\equiv 0)$. Then the problem $(1.5)$ has a unique nonnegative global classical solution $(u, v)\in[C(\bar{\Omega}\times[0, \infty))\cap C^{2, 1}(\bar{\Omega}\times(0, \infty))]^2$ satisfying

where $C > 0$ is a constant independent of t, and in particular $0 < u\leq K_0$, where

It's easy to see that the system (1.5) admits four non-negative solutions:

$({\rm{i}})$ the trivial solution $E_0 = (0, 0)$;

$({\rm{ii}})$ the semi-trivial solutions $E_1 = (0, \dfrac{d}{e})$ and $E_2 = (\dfrac{\lambda}{a}, 0)$;

$({\rm{iii}})$ there exists a unique positive constant solution $E^{*} = :(u^*, v^*)$ when

holds, where

The following important property on positive constant equilibrium $E^*$ can be presented.

Theorem 1.2 (global stability). If condition $(H_0)$ and the following conditions are satisfied

then the positive constant equilibrium $E^{*}$ is globally asymptotically stable. Furthermore, it follows that

where the convergence is exponential.

Herein, we briefly outline the plan of this paper: Section 2 proves some estimates and the local existence of the global classical solutions; Section 3 addresses the boundedness and global existence of solutions; Section 4 analyzes the global stability of co-existence steady state.

2.   Local existence and preliminaries

In what follows, we shall abbreviate $\int_{\Omega}f dx$ as $\int_{\Omega}f$ and $\|f\|_{L^2(\Omega)}$ as $\|f\|_{L^2}$ for simplicity and use $c_i(i = 1, 2, 3 \cdots)$ to denote a generic constant which may vary in the context. We first state the existence of local-in-time classical solution of the system (1.5) by using the abstract theory (cf. ^[19]).

Lemma 2.1 (Local existence). Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with smooth boundary. Assume $(u_0, v_0)\in [W^{1, \infty}(\Omega)]^2$ with $u_0, v_0\geq 0(\not\equiv 0)$.Then there exists a positive constant $T_{\max}\in(0, \infty]$ (the maximal existence time) such that the problem $(1.5)$ has a unique classical solution $(u, v)\in[C(\bar{\Omega}\times [0, T_{\max}))\cap C^{2, 1}(\bar{\Omega}\times(0, T_{\max}))]^2$ satisfying $u, v\geq0$ for all $t > 0$. Moreover, we have

Proof. The local-in-time existence and uniqueness of classical solution to the problem (1.5) follow from Amann's theorem ^[20]. The specific proof steps can refer to [21,Lemma 2.1].

Lemma 2.2. Let $\Omega\subset\mathbb{R}^n$ be a bounded domain with sufficient smooth boundary. Under the conditions Theorem 1.1, the solution $(u, v)$ of the system $(1.5)$ satisfies

where $K_0$ is defined by $(1.7)$, and it further follows that

Proof. The proof procedure refers to Lemma 2.2 in ^[13].

Lemma 2.3 (see ^[21,25]). Let $\Omega\subset\mathbb{R}^n$ be a bounded domain with sufficient smooth boundary. Let $T\in(0, \infty]$ and suppose that $u\in C^0(\overline{\Omega}\times[0, T))\cap C^{2, 1}(\overline{\Omega}\times(0, T))$ is a solution of

where $g_0(u, v) = u+\lambda u-au^2-\dfrac{quv}{1+bu+cv}$ and $g_0(u, v)\in L^{\infty}((0, T); L^p(\Omega))$, then there exists a constant $C_1$ such that

Lemma 2.4. Let $(u, v)$ be the solution of the system $(1.5)$, then there exist two positive constants $M$ and $C_2$ such that

and

where $\gamma = \dfrac{e}{K_0+d}$, $\tau = :\min\{1, \frac{1}{2}T_{\max}\}$ and $\widetilde{T}_{\max} = :\begin{cases}T_{\max}-\tau, & {{if}}\; \; \; T_{\max} < \infty, \\\infty, & {{if}} \; \; \; T_{\max} = \infty.\end{cases}$

Proof. By means of $u(x, t)\leq K_0$, the second equation of the system (1.5) becomes

Integrating this equation over $\Omega$, it follows that

Then applying the Cauchy-Schwarz inequality, we have $\gamma\int_{\Omega}v^2\geq\frac{\gamma}{|\Omega|}(\int_{\Omega}v)^2$ which implies

It is obvious that $v(\cdot, x)$ satisfies (2.4) by the ODE methods. Then integrating (2.7) over $(t, t+\tau)$ and using (2.4), we can obtain (2.5).

Lemma 2.5. Let $(u, v)$ be the solution of the system $(1.5)$, then there exist two positive constants $C_3$ and $C_4$ independent of $t$ such that

and

where $\tau$ and $\widetilde{T}_{\max}$ are defined by Lemma 2.4.

Proof. We multiply the first equation of the system (1.5) by $-\Delta u$, and integrate the result by parts to have

then

where $A = \frac{q^2K_0^2}{d_{1}(1+bK_0)^2}$. Multiplying (2.7) by $\frac{A}{\gamma}$ and adding the result to (2.11), which yields

Adding $\int_{\Omega}|\nabla u|^2+\frac{A}{\gamma}\int_{\Omega}v$ to the both sides of this equation, we can get

By the sobolev interpolation inequality and Lemma 2.2, we have for any $\varepsilon > 0$ and a constant $m^* = :C_\varepsilon K_0^2|\Omega|$ that

which updates (2.13) to

Let $y(t) = :\int_{\Omega}|\nabla u|^2+\frac{A}{\gamma}\int_{\Omega}v$, we obtain that

By the Gronwall inequality, there exists $T > 0$ such that $t > T$, we have $y(t)\leq y(T)+m^{**}(T_{\max}-T) = :M^*$. Hence, $\|\nabla u\|_{L^2} = (\int_{\Omega}|\nabla u|^2)^{\frac{1}{2}} < (M^*)^{\frac{1}{2}} = C_3$. On the other hand, by integrating (2.13) over $(t, t+\tau)$ and further calculating, we have $\int_{t}^{t+\tau}\int_{\Omega} |\Delta u|^2\leq\frac{1}{d_1}[(2\lambda+1)C_3\tau+\frac{2A}{\gamma}M] = C_4.$

3.   Boundedness of solutions

In this section, we shall use some related estimates derived in the previous section to further show the boundedness and existence of global classical solutions for the system (1.5). Motivated by Jin, Kim and Wang ^[25], we shall prove the following Gronwall-type inequality

which yields the uniform-in-time boundedness of $\|v(\cdot, t)\|_{L^2}$. Based on the parabolic regularity, we can get $\|v(\cdot, t)\|_{L^{\infty}}$ is uniformly bounded, which along with Lemma 2.1 extends a local solution to a global one.

Lemma 3.1 ($L^2$-estimate). Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with smooth boundary. If $(u, v)$ is a solutionof the system $(1.5)$, then there exists a constant $C_5 > 0$ such that

Proof. Multiplying the inequality (2.6) by $v$ and integrating the results over $\Omega$, we have

And applying the $H\ddot{o}lder$ inequality and Young's inequality, we have

and

It follows that

where $c_1 = \frac{32}{27\gamma^2}|\Omega|$. According to Gagliardo-Nirenberg inequality, we can get

Using the fact $\|\nabla u\|_{L^2}\leq c_4$ in Lemma 2.5, we derive from

Substituting (3.4) and (3.5) into (3.3) gives

where

then

where $c_6 = \frac{(\chi^2c_2c_5)^2}{d_2^3}, \; c_7 = \frac{(\chi^2c_2c_5+d_2^2)^2}{2d_2^3}$. It follows that

Furthermore, we can get the following estimate of the second term to the right of the inequality

Finally, letting $c_8 = :\frac{4c_7^3}{27\gamma^2}|\Omega|+c_1$, one has from (3.7) that

Noting (2.5) and (2.10), the rest of this proof is completed by using the same proof method as [25,Theorem 3.1].

Lemma 3.2 ($L^\infty$-estimate). Suppose that the conditions in Lemma 3.1 hold, then the solution of the system $(1.5)$ satisfies

where the constant $C_6 > 0$ independent of t.

Proof. Using $v^{p-1}$ with ${p\geq 2}$ as a test function for the equation (2.6) and integrating the results over $\Omega$, we have

Adding $\int_{\Omega}v^p$ to both sides of the above equation and using the $H\ddot{o}lder$ inequality and Young's inequality, we end up with

which implies

Multiplying the above inequality by $p$ and integrating the results over $\Omega$, we obtain

where $c_9 = \frac{2}{p+1}(\frac{2p}{(p+1)\gamma})^p|\Omega|$. By means of

the inequality (3.9) becomes

Noting the fact $\|v(\cdot, t)\|_{L^2} < C_5$ and $\|\nabla u(\cdot, t)\|_{L^4} < c_{10}$ in Lemma 3.1 and Lemma 2.3, then one has

Owing to Gagliardo-Nirenberg inequality, we have

Define $c_{12} = \frac{c_{10}^2c_{11}\chi^2p(p-1)}{2d_2}$, it follows that

which together with (3.10) gives

where

Through Gronwall's inequality and (3.11), we can derive

Then choosing $p = 4$ in (3.12) and using Lemma 2.3, we can find a constant $c_{14}$ independent of $p$ such that $\|\nabla u(\cdot, t)\|_{L^{\infty}} < c_{14}$. Then applying Moser iteration procedure (cf. ^[22]), one has (3.8). This completes the proof.

On account of Lemma 3.2 and Lemma 2.3, we can get the global boundedness of solutions to (1.5) by the Moser iteration procedure (cf. ^[19]). Next, we will show the following results on the global existence of solutions.

Lemma 3.3 (global existence).Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with smooth boundary. Assume $(u_0, v_0)\in[W^{1, \infty}(\Omega)]^2$ with $u_0, v_0\geq0(\not\equiv 0)$, then the system $(1.5)$ has a unique global classical solution

satisfying (1.6).

Proof. From Lemma 3.2 and Lemma 2.3, which together with the local existence results in Lemma 2.1 completesthe proof of this Lemma.

4.   Proof of Theorem 1.2

We are now in the position to derive the global stability of $E^* = (u^*, v^*)$.

Proof. Let $(u(x, t), v(x, t))$ be any solution of the system (1.5), we construct the following Lyapunov function

It is clear from the fact $W^*(\omega) = 0$ if $\omega = (u^*, v^*)$ and $W^*(\omega) > 0$ for all $\omega\neq(u^*, v^*)$. That is, $W^*(\omega)$ is a positive definite function. Furthermore, from definition of $W^*$ and results of Theorem 1.1, we have $W^*(\omega)\leq C_7$ for a constant $C_7 > 0$ independent of $t > 0$ (see ^[13,21]). Next, we take the derivative of $W^{*}$ with regard to $t$ along the trajectory of the system (1.5) and arrive at

where

and

where

The equation (4.1) can be further written as

It is obvious that $I_{22} < 0$ if and only if the matrix in the integrand of (4.2) is positive definite, which is equivalent to $k(u, v) > 0$ and $\rho(u, v) = k(u, v)m(u, v)-l^2(u, v) > 0$, where

By calculation and the condition, we can get

then $I_{21} < 0$ can be determined by the first fraction in (1.8) and we can see that $k(u, v) > 0$ and $\rho(u, v) > 0$ from the second fraction in (1.8). Here it is clearly that the coexistence state $(u^*, v^*)$ is globally asymptotically stable by the LaSalle's invariant principle and there exists a $t_0 > 0$ so that for all $t > t_0$ the following inequality holds:

for the specific procedures of the above equation, we can refer to the proof of [13,Lemma 4.3] and [21,Lemma 4.5] which further yields the exponential decay rate in $L^{\infty}-$norm from (1.9).

Remark 1. Theorem 1.2 discusses the global stability under the assumption that $\chi > 0$. If $\chi\leq0$, the lighter condition $q < \frac{acde}{(K_0+d)(be+c)}$ is needed to satisfy the global stability. That's to say, if there is no prey-taxis phenomenon $(\chi = 0)$ or the prey can gather to form a group that can resist foreign enemies $(\chi < 0)$, the co-existence steady state is globally asymptotically stable when the competition between predators and preys is weak. Once prey-taxis phenomenon occurs $(\chi > 0)$, the above state may require weaker competitiveness to maintain its global stability.

Acknowledgments

Partially supported by National Natural Science Foundation of China (No. 11571086), National Natural Science Foundation of Heilongjiang Province, China(No. LH2020A019) and Harbin Normal University Innovation Fund, China (No. HSDSSCX2021-16).

Conflict of interest

The authors declare there is no conflicts of interest.