Global existence and stability of three species predator-prey system with prey-taxis

1.   Introduction

Predator-prey models were developed to describe the dynamics of interactions between prey and predator species. The relationship between prey and predator has been explored in recent years due to its importance in ecology. In addition to the differential operators in the predator-prey system, predators also move toward the higher prey density, which is so-called the prey-taxis, and it plays an important role in pest control and biological control ^[1,2,3,4]. The first predator-prey model with prey-taxis was derived by Kareiva and Odell ^[5] to describe the predator aggregation phenomenon:

where $u = u(x, t)$ and $w = w(x, t)$ denote the predator and prey densities, respectively, and the term $\nabla\cdot(d(w)\nabla u)$ denotes the diffusion of predators with diffusion coefficient $d(w)$. The term $-\nabla\cdot(u\chi(w)\nabla w)$ represents the prey-taxis with $\chi(w)$ as prey-taxis coefficient. The parameter $d > 0$ is the diffusion coefficient of prey. The typical form of the functions $F(u, w) = auf(w)+h(u)$ and $G(u, w) = g(w)-buf(w),$ where $f(w)$ represents the functional response, for numerous functional response functions (see ^[6,7,8]) and the parameters $a, b\in \mathbb{R}$ describe the inter-specific interactions between predator-preys. The intra-specific interactions of predators and prey are described by the functions $h(u)$ and $g(w)$, respectively. The results related to variants of the above prey-taxis system have been studied by many authors, as one can refer to ^{[9,10,11,12,13,14,15,16,17,18]}, and nonlinear prey-taxis ^{[19,20,21,22,23,24]}. Moreover, the predator-prey system with prey-taxis and liquid surroundings was considered in ^[25], and proved global existence and large time behavior of solutions by using $L^p$ estimates and Lyapunov functionals, respectively.

In this paper, we consider a PDE model of indirect interactions between two prey species and a shared predator with homogeneous Neumann boundary conditions:

Where, $\Omega \subset \mathbb{R}^n (n \geqslant 1)$ is a bounded domain with smooth boundary $\partial \Omega$ and $\frac{\partial}{\partial \nu}$ represents the derivative with respect to outer normal of $\partial \Omega$, $u$ is the native prey, $v$ denotes the invasive prey and $w$ is the predator species. The parameters $d_1, d_2$ denote the random diffusion rates of prey, and $d_3$ and $d_4$ denote the random diffusion rate of predators and the chemical concentration, respectively. $K_u$ and $K_v$ are the carrying capacities for these prey species. The constants $\alpha_1$and $\alpha_2$ are intrinsic growth rate parameters. The parameters $T_1$ and $T_2$ usually represent the handling time required for catching and consuming a unit of prey type $u$ and $v$, respectively. The constants $c_1$ and $c_2$ are capture rates per unit prey density while the predator is searching. In particular, $c_1$ is the capture rate of prey $u$ and $c_2$ is the capture rate of prey $v$. In addition, $d$ is an intrinsic growth (death) rate for the predator and $\eta$ is a self-limiting or crowding coefficient for the predator. $\kappa$ is the production rate of chemical signal per individual prey $u$ and $\xi$ is the decay rate of the chemical signal. The positive constants $\gamma_1$ and $\gamma_2$ denote the transformation rates of the predator.

The terms $-\nabla\cdot(\chi_1w\nabla u)$ and $-\nabla\cdot(\chi_2 w\nabla v)$ denote the tendency of predators moves towards the high density of prey. The parameters $\chi_1$ and $\chi_2$ are the prey-taxis coefficients. The functions $\frac{c_1u}{1+c_1T_1u+c_2T_2v}$ and $\frac{c_2v}{1+c_1T_1u+c_2T_2v}$ these represent Holling type Ⅱ functional responses for two preys that are consumed in a single habitat, so that handling one prey reduces the time available to capture the other.

Let $\tilde{u} = \frac{u}{K_u}, \tilde{v} = \frac{v}{K_v}, \tilde{w} = dw,$ $d = \frac{d_1}{d_3}, \eta = \frac{d_2}{d_1},$ $L = \sqrt{\frac{d_3}{\alpha_1}}, T = \frac{L^2}{d_3} = \frac{1}{\alpha_1},$ $\tilde{t} = \frac{t}{T}, \tilde{x} = \frac{x}{L}, \tilde{y} = \frac{y}{L}, r = \alpha_2 T a_1 = c_1Td,$ $a_2 = c_1T_1K_u, a_3 = c_2T_2K_v, a_4 = c_2dT, a_5 = c_1\gamma_1K_ud, a_6 = c_2\gamma_2K_vd.$ Then, substituting these parameters into system (1.2) and dropping the tilde notation, we get a nondimensional system as follows:

Let us recall some existing works on three species predator-prey systems with prey-taxis. Very recently, Haskel and Bell ^[26] proved the existence of positive classical solutions for the two-prey one-predator system with prey competitions and prey taxis. In addition, they also established the pattern formation by using bifurcation analysis. Further, they also studied the bifurcation analysis of two competing prey with one shared predator model by using the theories of Crandall-Rabinowitz and Hopf bifurcation in ^[27]. The steady-state bifurcation analysis of the two-prey one-predator model with two prey taxis was studied by Xu et al. ^[28]. Jin et al. ^[29] considered the three-species food chain model in a two-dimensional bounded domain, and they also proved the global existence of classical solutions and global stability of constant steady states. The traveling wave solutions for a nonlocal dispersal predator-prey system with one predator and two prey was studied in ^[30]. Amorim et al. ^[31] studied the boundedness and global well-posedness of the spatio-temporal evolution of two competitive prey, and one predator model with the intra-specific competition. The global existence and boundedness of classical solutions for the two-predators and one-prey with competition in a bounded domain with Neumann boundary conditions were proved by Min et al. in ^[32]. Recently, the global existence of weak solutions to the two-prey one-predator system with prey-taxis, and competition between prey was proved in any dimension in ^[33]. For the similar mathematical structure of (1.3), we refer to ^[34,35,36]. Throughout this paper, we assume the system parameters are positive. To the best of author's knowledge, there is no article that discusses the well-posedness of the considered system (1.3). The main purpose of this article is to discuss the global dynamics of the system (1.3) in any dimension $(n \geqslant 1)$. In particular, we first prove the global existence of a classical solution in all dimensions, and then we investigate the global stability of steady states. Our main result regarding the global existence of classical solutions with uniform-in-time bound is stated below.

Theorem 1.1. $($Global existence$)$ Let $\Omega\subset\mathbb{R}^n, n \geqslant 1$ be a bounded domain with smooth boundary and let $d, \eta > 0, h_1, h_2 > 1, k \geqslant 2, a_i > 0, i = 1, \ldots, 6, \mu > 0, K_0 = \max\{1, \|u_0\|_{L^\infty}\}$ and $K_1 = \max\{1, \|v_0\|_{L^\infty}\}$. For any $(u_0, v_0, w_0)\in [W^{1, p}(\Omega)]^3$ with $p > n$ and $u_0, v_0, w_0 \geqslant 0(\not\equiv 0),$ if $d > \max\{5kK_0, \frac{5kK_1}{\eta}\}$ and $\chi_1$ satisfies

and $\chi_2$ satisfies

then there exists a unique global classical solution $(u, v, w)\in [C^0(\overline{\Omega}\times [0, \infty))\cap C^{2, 1}(\overline{\Omega}\times(0, \infty))]^3$ solving the problem (1.3). Moreover, the solution satisfies $u, v, w > 0$ for all $t > 0$ and

where $C > 0$ is a constant independent of $t.$

Next, we shall study the large time behaviour of the constant steady states $(u_s, v_s, w_s)$ of the system (1.3) solving the following system

If we solve the above system, we will find the following steady states

where

and $(0, 0, 0)$ is the extinction steady state, $(1, 0, 0)$ is the prey $u$ only steady state, $(0, 1, 0)$ is the prey $v$ only steady state. $E_*^1$ and $E_*^2$ denote the semi-coexistence steady state. Finally, $E_*$ denotes the coexistence steady state. Next, we shall explore the following question: which of the above seven homogeneous steady states will be asymptotically stable? As we know that the global stability of the cross-diffusion system is difficult and many techniques are not available, we try to use the Lyapunov functionals to prove the global stability of the homogeneous steady states under some conditions. Next, we state our stability results as in the following theorem:

Theorem 1.2. $($Global stability$)$ Assume the conditions in Theorem 1.1 hold. Let $(u, v, w)$ be the solution of (1.3) obtained in Theorem 1.1 and let $K_0 = \max\{1, \|u_0\|_{L^\infty}\}, K_1 = \max\{1, \|v_0\|_{L^\infty}\}, \Gamma_1 = \frac{a_5+(a_3a_5-a_2a_6)v_*}{a_1(1+a_2u_*+a_3v_*)}$ and $\Gamma_2 = \frac{a_5+(a_3a_5-a_2a_6)u_*}{a_1(1+a_2u_*+a_3v_*)}.$ Then the following results hold true.

● If $\mu > \frac{a_5+a_6}{1+a_2+a_3}, a_4 \leqslant \frac{a_1r(a_5+a_3a_5-a_2a_6)}{a_3a_5-a_6-a_2a_6},$ then the steady state $(1, 1, 0)$ is globally asymptotically stable.

● If $\mu < \frac{a_5+a_6}{1+a_2+a_3}, a_4 > \frac{a_1r(a_5+a_3a_5-a_2a_6)}{a_3a_5-a_6-a_2a_6},$ the steady state $E_*$ defined by (1.6) is globally asymptotically stable provided

where $\|u\|_{L^\infty}$ and $\|v\|_{L^\infty}$ depends on $d, \eta, a_1, a_2, a_3$ but independent of $\chi_1, \chi_2.$

The paper is organized as follows: In section 2, we first present some preliminary results, and then we state and prove the local existence of solutions of (1.3). Section 3 deals with the existence of globally bounded classical solutions as stated in Theorem 1.1. In Section 4, we establish the global stability results as stated in Theorem 1.2 by using the Lyapunov functionals and LaSalle's invariance principle.

2.   Preliminaries and local existence

In what follows, we shall abbreviate $\int_\Omega f \mathrm{d} x$ as $\int_\Omega f$ for simplicity. In this section, we first prove the local existence of classical solutions to the system (1.3) in any dimension $\Omega\subset\mathbb{R}^n, n \geqslant 1$ using the Amann's approach (cf. ^[37,38]). We use the following Gagliardo-Nirenberg interpolation inequality in the sequel.

Lemma 2.1. $($^[39]$)$ There exists a constant $C_4 > 0$, such that for all $u\in W^{1, q}(\Omega),$

where $p, q \geqslant 1$ which satisfies $p(n-q) < nq, m\in (0, p)$ with $a = \frac{\frac{n}{m}-\frac{n}{p}}{\frac{n}{m}+1-\frac{n}{q}}\in (0, 1).$

Lemma 2.2. $($^[39]$)$ There exists a constant $C_5 > 0$, such that for all $u\in W_{1, q}(\Omega)$, we have

where $p > 1$ and $q > 0.$

Lemma 2.3. $($^[40]$)$ Let $T > 0, \tau\in(0, T), \sigma \geqslant 0, a > 0, b \geqslant 0,$ and suppose that $f:[0, T)\rightarrow [0, \infty)$ is absolutely continuous, and satisfies

where $h \geqslant 0, h(t)\in L_{loc}^1([0, T))$ and

Then,

Lemma 2.4. $($^[39]$)$ Assume that $m\in\{0, 1\}, p\in[1, \infty]$ and $q\in (1, \infty).$ Then there exists some positive constant $C_1,$ such that

for any $\phi\in D((A+1)^\theta)$, where $\theta\in (0, 1)$ satisfies

If in addition $q \geqslant p,$ then there exist constants $C_2 > 0$ and $\gamma > 0$, such that for any $\phi\in L^p(\Omega),$

where the semigroup $\{e^{-t(A+1)}\}_{t \geqslant 0}$ maps $L^p(\Omega)$ into $D((A+1)^\theta).$ Moreover, for any $p\in(0, \infty)$ and $\epsilon > 0,$ there exist constants $C_3 > 0$ and $\gamma > 0$, such that

this is valid for all $\mathbb{R}^n-$ valued $\phi\in L^p(\Omega).$

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

which satisfies $(u, v, w) > 0$ for all $t > 0.$ Further,

Proof. Denote $\psi = (u, v, w).$ Then, the problem (1.3) can be written as

where

Since the eigenvalues of $A(\psi)$ are positive, the system (1.3) is normally parabolic(cf. ^[37,38]). Then, the application of Theorem 7.3 and Corollary 9.3 in ^[37] yields a $T_{max} > 0$, such that system (1.3) admits a unique solution $(u, v, w)\in [C^0(\overline\times[0, T_{max}))\cap C^{2, 1}(\overline{\Omega}\times (0, T_{max}))]^3.$ Next, we show the nonnegativity of $(u, v, w)$ by using the maximum principle. To do so, we need to rewrite the third equation of the system (1.3) as follows:

where $p_1(x, t) = \chi_1\nabla u+\chi_2\nabla v$ and $p_2(x, t) = \chi_1\Delta u+\chi_2\Delta v-\frac{a_5u+a_6v}{1+a_2u+a_3v}.$ Hence, we apply the maximum principle for parabolic equation with Neumann boundary condition to (2.12) and we get $w \geqslant 0$ for all $(x, t)\in \Omega\times (0, T_{max}).$ In addition, we also obtain $w > 0$ by strong maximum principle since the initial function $w_0\not\equiv0.$ In the same way, we can obtain that $u, v > 0$ for all $(x, t)\in \Omega\times (0, T_{max}).$ Because $A(\psi)$ is lower triangular, (2.9) follows from Theorem 5.2 in ^[41].

Lemma 2.5. The solution $(u, v, w)$ of the system (1.3) satisfies

where $K_0, K_1, K_2, \delta = \max\{a_4a_5\|u_0\|_{L^1}+a_1a_6\|v_0\|_{L^1}+a_1a_4\|w_0\|_{L^1}, \frac{c_2}{c_1}\}$ are positive constants independent of $t$.

Proof. The proof is similar to Lemma 2.2 in ^[21] but for reader's convenience, we provide the proof here. We have already proved that the solution $(u, v, w)$ of the system (1.3) is non-negative. Using this fact, we have

Let $u^*(t)$ be a solution to the following ODE problem

The solution of the above ODE satisfies $u^*(t) \leqslant K_0 = \max\{\|u_0\|_{L^\infty}, 1\}$, and in addition, $u^*(t)$ is a super solution of the following PDE problem

Hence, we have $U(x, t) \leqslant u^*(t)$ for all $(x, t)\in \overline{\Omega}\times (0, \infty).$ Using the strong maximum principle to the problem (2.18), we obtain $0 < U(x, t) \leqslant u^*(t)$ for all $(x, t)\in \overline{\Omega}\times (0, \infty)$. From (2.16)–(2.18), and using the comparison principle, we conclude that

which yields (2.13). Similarly, we can also prove (2.14).

Multiplying the first, second and third equations of (1.3) by $a_4a_5$, $a_1a_6$ and $a_1a_4$, respectively, and adding the resulting equations and then integrating it over $\Omega$, we get

Using Cauchy-Schwartz inequality, one has

which implies

Using Young's inequality $(-a^2-b^2 \leqslant -2ab, a, b > 0)$ yields

Substituting (2.22) into (2.20), we have that

Set $y(t) = \int_\Omega a_4a_5 u+ a_1a_6 v+ a_1a_4 w$ and choose $c_1 = \min\{1, r, \mu\}$ then (2.23) can be written as

where $c_2 = (a_4a_5+a_1a_6)|\Omega|$ and which yields (2.15) with the help of Gronwall's inequality. We further have from (2.17) that $\lim_{t \to \infty} \sup u(x, t) \leqslant 1.$ The proof of Lemma 2.5 is complete.

3.   Global existence and boundedness

In this subsection, we prove the global existence and boundedness of solutions. In order to prove the global existence, we first derive a uniform bound for $w$ in $L^{n+1}$ by using a weight function argument and the proof is inspired from ^[23], which also concerns the predator-prey taxis with a single prey population. It is worth mentioning that the method was initially developed in ^[42].

Lemma 3.1. Assume that $\chi_1$ and $\chi_2$ satisfy (1.4) and (1.5), respectively, and let $(u, v, w)$ be the solution of (1.3). Then, there exists a positive constant $c_0 > 0$, such that

Proof. Let us define the constants and weight functions

Now, by using the weight function and the first and second equation of (1.3), we obtain

where we used the boundedness of the functional response, $C > 0$. The above inequality can be written as

By using Young's inequality, we get

and

Again,

Substituting (3.7)–(3.10) into (3.6), one has

Now, we multiply the third equation of (1.3) by $\varphi_2(v)$, we have

By using Young's inequality, we arrive at

and

Substituting (3.13)–(3.16) into (3.12), we derive that

Now, adding (3.11) and (3.17), the resulting inequality becomes

Now, we do some computation to show that the terms involving $|\nabla u|^2$ and $|\nabla v|^2$ on the right-hand side of above inequality are dominated by $\int_\Omega w^k\varphi_1^{\prime\prime}|\nabla u|^2$ and $\int_\Omega w^k\varphi_2^{\prime\prime}|\nabla v|^2$, respectively. For $s \geqslant 0,$ define

and

Now, combining (3.19) and (3.20), one has

Similarly, we combine (2.1) and (2.2), we have that

Substituting (2.5) and (3.24) into (3.18), one has

where $c_1 = C+\frac{2\beta_1^2K_0^2}{k}$ and $c_2 = C+\frac{2r\beta_2^2K_1^2}{k}.$ Using Lemma 2.1, Lemma 2.2, and (3.3), we get the estimate

where $a = \frac{\frac{kn}{2}-\frac{n}{2}}{\frac{kn}{2}+1-\frac{n}{2}}\in (0, 1).$ Now using the fact (3.3) and from (3.26), one can obtain

Similarly, we get

From (3.27) and (3.28), we have

for all $t\in (0, T_{max})$ and where $\frac{1}{a} > 1.$ Set $y(t): = \frac{1}{k}\int_\Omega w^k(\varphi_1(u)+\varphi_2(v)).$ By using the inequality $x^p+y^p \geqslant n^{1-p}(x+y)^p, p \geqslant 1,$ we get

where $C_8: = \min\bigg\{\frac{(k-1)}{k^2C_6^{\frac{1}{a}}}, \frac{(k-1)}{k^2C_7^{\frac{1}{a}}}\bigg\} n^{1-\frac{1}{a}}$ and $\frac{1}{a} > 1.$ By using Lemma 2.3 and the fact that (3.3), one has

for all $t\in (0, T_{max})$, where $C(u_0, v_0, C_8, \tau) > 0.$ The proof of Lemma 3.1 is complete.

Lemma 3.2. Let $(u, v, w)$ be a solution of the system (1.3). Then, there exists a positive constant $c > 0$, such that

Proof. To obtain the $L^\infty-$ bound of $w$, we use the semigroup estimates. In order to do this, we first obtain that for any $\tau\in(0, T_{max}),$ there exists a constant $C > 0$, such that

Let $\tau\in (0, T_{max})$ be given such that $\tau < 1$ and also let $q: = n+1$ and $\theta\in \bigg(\frac{1}{2}(1+\frac{n}{q}), 1\bigg).$ To begin with, we rewrite the first equation of (1.3) as follows:

with $g(x, t): = u(1-u)- \frac{a_1uw}{1+a_2u+a_3v}+u.$ Then, by Lemma 3.1 and the fact that $0 < u \leqslant K_0, 0 < v \leqslant K_1$ (see Lemma 2.5) and $\frac{a_1u}{1+a_2u+a_3v} \leqslant \tilde{K}$, $\tilde{K} > 0$, we have

We apply the variation-of-constants formula to (3.34) and obtain

where $A_d = -d\Delta$. Then using (2.6), (2.7) and the estimate (3.35) in (3.36), one can derive

for all $t\in (\tau, T_{max})$, where $\Gamma(1-\theta) > 0.$ From the last inequality (3.37), we get the desired estimate

where $C_1$ is a generic constant which may vary line to line. Next, we obtain the bound for $\|v(\cdot, t)\|_{W^{1, \infty}}.$ To this end, we rewrite the second equation of (1.3) as follows:

with $h(x, t): = rv(1-v)- \frac{a_4vw}{1+a_2u+a_3v}+v.$ Then, by Lemma 3.1 and the fact that $0 < u \leqslant K_0, 0 < v \leqslant K_1$ (see Lemma 2.5) and $\frac{a_4v}{1+a_2u+a_3v} \leqslant \overline{K}$, $\overline{K} > 0$, we have

We apply the variation-of-constants formula to (3.34) and obtain

Then, using (2.6), (2.7) and the estimate (3.40) in (3.41), we find

for all $t\in (\tau, T_{max})$ where $\Gamma(1-\theta) > 0.$ From the last inequality (3.42), we obtain

Next we derive the $L^\infty$- bound of $w(\cdot, t)$. We rewrite the third equation of (1.3) as follows:

Then, applying the variation-of-constants formula to (3.44), one has

where $f(u, v, w) = \frac{a_5uw+a_6vw}{1+a_2u+a_3v}+(1-\mu)w.$ Now, we take the $L^\infty$-norm on both sides of the above equation, we have

First, we estimate the term $I_1$ as follows:

for all $t\in (\tau, T_{max})$ and $\theta\in(\frac{n}{2q}, 1)$ and $\mu > 0.$ In order to estimate the term $I_2$, we set $m = 0, q = n+1, p = \infty$ for (2.8) and we use (3.33) and (3.1), one has

where $\Gamma(\frac{1}{2}-\theta-\epsilon)$ is a Gamma function which is positive since $\frac{1}{2}-\theta-\epsilon > 0$ and $\mu, C_5 > 0.$

Next, we obtain the bound for $I_3.$ As in the estimate of $I_2,$ we set $m = 0, q = n+1, p = \infty$ for (2.8) and we use (3.33) and (3.1), one has

where $\Gamma(\frac{1}{2}-\theta-\epsilon)$ is a Gamma function which is positive since $\frac{1}{2}-\theta-\epsilon > 0$ and $\mu, C_8 > 0.$ Finally, we obtain the bound for $I_4$. To this end, we use (2.6) and (2.7) and let $m = 1, p = (n, \infty]$ and $q = n+1$. Hence, we can choose $\theta\in (\frac{1}{2}(1-\frac{n}{p}+\frac{n}{q}), 1)$. Then one has

Using the fact that $0 < u \leqslant K_0, 0 < v \leqslant K_1$ and(3.1), we can get

for all $t\in (\tau, T_{max}).$ Hence, we have

for all $t\in (\tau, T_{max})$ and where $\Gamma(1-\theta)$ is a Gamma function and it is positive since $1-\theta > 0$ and $\nu > 0.$ Since $p > n$, Sobolev embedding theorem yields that

Substituting the estimates (3.47), (3.48), (3.49), (3.53) into (3.46) which yields (3.32). Hence, this completes the proof.

Proof of Theorem 1.1. From Lemma 2.5, we obtain $\|(u(\cdot, t), v(\cdot, t))\|_{L^\infty(\Omega)} \leqslant C$. Further, we also obtain the bound for $\|w(\cdot, t)\|_{L^\infty}$ from Lemma 3.2. By noticing these results now, we can conclude that

where $c$ is a positive constant. From the criterion (2.9), we obtain that $T_{max} = \infty$ and hence $\|u(\cdot, t)\|_{L^\infty}+\|v(\cdot, t)\|_{L^\infty}+\|w(\cdot, t)\|_{L^\infty} \leqslant c\ \mbox{for all} \ t\in (0, \infty).$ The proof of Theorem 1.1 is complete.

4.   Global stability of solutions

In this section, we shall prove the global stability of solutions of (1.3) by constructing some suitable Lyapunov functionals, and then we use the LaSalle's principle.

4.1.   Global stability of Prey only state

Lemma 4.1. Let $(u, v, w)$ be the solution of (1.3) and let $\Gamma_1 = \frac{(a_5(1+a_3-a_2a_6))}{a_1(1+a_2+a_3)}, \Gamma_2 = \frac{(a_6(1+a_2)-a_3a_5)}{a_4(1+a_2+a_3)}.$ Then, if $\mu \geqslant \frac{a_5+a_6}{1+a_2+a_3}$ and $a_4 < \frac{a_1r(a_5+a_3a_5-a_2a_6)}{a_3a_5-a_6-a_2a_6}$, it holds that

Proof. Let us define the energy functional

First we need to prove that $\mathcal{E}(t) \geqslant 0$ and $\mathcal{E}(t) = 0$ iff $(u, v, w) = (1, 1, 0).$ To this end, let $\psi(x) = x-g_*\ln x$ and by Taylor's formula, one has

where we choose $\delta$ in between $g_*$ and $g.$ Now putting $g = u$ and $g_* = 1$ in the last equation, we have

Similarly, we can show that

Therefore, we get $\mathcal{E}(u, v, w) = \mathcal{E}(1, 1, 0) = 0$ and $\mathcal{E}(u, v, w) \geqslant 0$ for $(u, v, w)\neq (1, 1, 0).$ Differentiating (4.2) with respect to $t$, and then substituting equations from (1.3), we have

By using the assumptions of $\Gamma_1$ and $\Gamma_2$ in (4.3), one has

which yields

for all $(u, v, w)$ and also the equality holds when $(u, v, w) = (1, 1, 0).$ At last, the LaSalle's invariance principle (cf. ^[43], Theorem 3) yields that the solutions $(u, v, w)$ converge to the constant steady state $(1, 1, 0)$ as time $t\rightarrow \infty$.

4.2.   Global stability of coexistence state

Lemma 4.2. Let $\Gamma_1 = \frac{a_5+(a_3a_5-a_2a_6)v_*}{a_1(1+a_2u_*+a_3v_*)}$ and $\Gamma_2 = \frac{a_6+(a_2a_6-a_3a_5)u_*}{a_4(1+a_2u_*+a_3v_*)}$ be positive constants. Let $(u, v, w)$ be the solution of (1.3). If $\mu < \frac{a_5+a_6}{1+a_2+a_3},$ $a_4 > \frac{a_1r(a_5+a_3a_5-a_2a_6)}{a_3a_5-a_6-a_2a_6}$ and

then it holds that

Proof. The coexistence steady state $(u_*, v_*, w_*)$ of (1.3) satisfies equations

Let us define the Lyapunov functional $\mathcal{E}(u, v, w)$ as

where

Next, we take the derivative of $\mathcal{E}(t)$ with respect to $t$ along the trajectory of the system (1.3) and we obtain

Using the definition of $\mathcal{F}_1(t)$, the first equation of (1.3) and the fact that $(1- u_*)-\frac{a_1 w_*}{1+a_2 u_*+a_3 v_*} = 0$, we estimate the term $I_1$ as follows:

Now, let us simplify the numerator in the integrand of the third integral on the R.H.S. of (4.13) as follows:

Substituting (4.14) into the last integral on the R.H.S of (4.13), we obtain

Again, inserting (4.16) into (4.13), we end up with

Similarly, we estimate the term $I_2.$ Using the definition of $\mathcal{F}_2(t)$, the second equation of (1.3) and the fact that $r(1- v_*)-\frac{a_4 w_*}{1+a_2 u_*+a_3 v_*} = 0,$ we estimate $I_2$ as follows:

Now, let us simplify the numerator in the integrand of the third integral on the R.H.S. of (4.18) as follows:

Substituting (4.19) into (4.18) and simplifying, we obtain

Finally, we estimate the term $I_3.$ Using the definition of $\mathcal{F}_3(t)$, the third equation of (1.3) and the fact that $-\mu +\frac{a_5 u_*}{1+a_2 u_*+a_3 v_*}+\frac{a_6 v_*}{1+a_2 u_*+a_3 v_*} = 0$, we find

Further, we simplify the numerator in the integrand of the third integral on the R.H.S. of (4.21), one has that

Now, we rewrite the fourth term in the R.H.S of (4.21) using the above estimate, we obtain

Substituting (4.22) into (4.21), we have

Furthermore, inserting (4.17), (4.20) and (4.23) into (4.12) and let $p(u, v) = \frac{1}{(1+a_2u+a_3v)(1+a_2 u_*+a_3 v_*)},$ we arrive at

Applying the Cauchy's inequality, we get

Inserting the last inequality (4.25) into (4.24), and letting $Z = (\frac{|\nabla u|}{u}, \frac{|\nabla v|}{v}, \frac{|\nabla w|}{w})$, we obtain

where

and the symmetric matrix is denoted by

The above matrix $B$ is positive definite if (4.8) holds. Therefore, we check that

and

Hence, there exists a positive constant $\alpha$, such that

Noting the facts that $1-u_*-\frac{a_1w_*}{1+a_2u_*+a_3v_*} = 0$ and $r(1-v_*)-\frac{a_4w_*}{1+a_2u_*+a_3v_*} = 0$, one has that

and

where we used the assumptions (4.6) and (4.7). Hence, we can conclude that

which yields that $\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{E}(t) \leqslant 0$ for all $u, v, w$ and the equality holds if $\nabla u = \nabla v = \nabla w = 0$. Therefore by applying LaSalle's invariance principle (cf. ^[43], Theorem 3) we can say that the solutions of (1.3) converges to the coexistence steady state $(u_*, v_*, w_*)$ as time $t\rightarrow \infty$.

Proof of Theorem 1.2. Theorem 1.2 is a consequence of Lemma 4.1 and Lemma 4.2.

Remark 4.1. We note that the condition (4.8) is a strong requirement, which implies that $\chi_1$ and $\chi_2$ have to be small enough.

Acknowledgments

The author is grateful to the referees for insightful comments which largely help improve the exposition of this paper. The author also thanks Prof. Zhi-An Wang from Hong Kong Polytechnic University for the fruitful discussions and valuable suggestions. This work was supported by the Postdoc Matching Fund Schemes of the Hong Kong Polytechnic University (Project ID P0036730 and a/c no. W18M).

Conflict of interest

The author declares that there are no conflicts of interest.