Local bifurcation of a Ronsenzwing-MacArthur predator prey model with two prey-taxis

1.   Introduction

In addition to random diffusion of the predator and the prey, the model comprises also a prey-taxis term, which means that the spatial-temporal variations of the predator's velocity are directed by prey gradient. A chemotaxis model includes the responses of predators to the distribution of resources, that is, the foraging behavior of predators that move actively toward the higher prey density due to the prey defenses. The models with a prey-taxis term may undergo rich dynamics and generate different spatial patterns from other models without the prey-taxis.

The predator-prey systems with prey-taxis have been widely investigated from different point of view in recent years. The existence and uniqueness of weak solutions to the two-species predator-prey model with one prey-taxis has been proved in ^[1]. The existence and uniqueness of weak solutions to an $n \times m$ reaction-diffusion-taxis system has been extended in ^[3]. The global existence of classical solutions to a three-species predator-prey model with two prey-taxis including Holling Ⅱ functional response has been investigate in ^[4]. Since the pattern formation of the attraction-repulsion Keller-Segel system has been studied in ^[7], many progress on extended models has been developed in ^[8,17,9,13]. The global existence and uniqueness of classical solutions to a predator-prey model with nonlinear prey-taxis has been shown in ^[12,6,18]. Some results about the global bifurcation of solutions for a predator-prey model with one prey-taxis is obtained in ^[16].

In this paper, we consider the following general Ronsenzwing-MacArthur model with two prey-taxes under Neumann boundary conditions:

where $\Omega\subset \mathbb{R}^N$ is a bounded domain with smooth boundary; $u$ represent the density of the predator, and $v, w$ express the densities of two preys; $c$ is the death rate of the predator; $\Phi(v) > 0, \Psi(w) > 0$ represent the functional response of the prey; $f(v), g(w)$ are the growth function of the prey respectively, which satisfy

$(1)$ $f(0) = f(k_1) = 0, g(0) = g(k_2) = 0$.

$(2)$

The terms $\alpha u \nabla v$ and $\beta u\nabla w$ are directed toward the increasing population density of $v$ and $w$, respectively. In this way, the predators move in the direction of higher concentration of the prey species, where ${\alpha}, {\beta}$ indicate their prey-tactic sensitivity.

Finally, we assume that all of the functions $f, g, \Phi, \Psi$ are of $C^1$ class functions on $\mathbb{R}_+.$ It follows from a standard approach (eg., cf. ^{[2,11,12,14,16]}) that the system (1.1), under these conditions, is well-posed for non-negative initial data $0\le (u_0, v_0, w_0)\in W^{2, p}(\Omega)^3$ with $p > N.$

The steady state solutions of the system (1.1) satisfy

Using condition (2) and the positivity of both $\Phi$ and $\Psi,$ we find that a constant vector $(u, v, w)$ with $u\ge 0$ and $v > 0, w > 0$ is a solution of (1.2) iff

Particularly, we see that $(u, v, w): = (0, k_1, k_2)$ is a constant non-negative steady state solution of (1.2) (for any $c$). Moreover, by suitably choosing the pairs $\{f, \Phi\}$ and $\{g, \Psi\}$ we can make that the first two eqnations in (1.2) admit strictly positive solutions $(u^*, v^*, w^*)$ satisfying $u^* > 0,$ $0 < v^* < k_1, 0 < w^* < k_2.$ Finally, by choosing $c = {\Phi(v^*)}+ {\Psi(w^*)}$ we find that the second equation (1.2) is satisfied.

Below we let $(u^*, v^*, w^*)$ be a constant positive steady state solution of (1.2). Take the prey taxis coefficient $\alpha$ as the main parameter, we analyze the solutions bifurcating from $(u^*, v^*, w^*)$ for the system (1.2). The case of taking the prey taxis coefficient $\beta$ as the parameter is similar. Our main bifurcation result will be given in the next §2.

We would like to point out that there are other forms of the Ronsenzwing-MacArthur model which include the conversion rate. It remains interesting to extend our method for such models.

2.   Main Results

Let $\mathbb{N}$ be the set of all positive integers, and $\mathbb{N}_0: = \mathbb{N}\bigcup\{0\}.$ Let $p > N,$ and set $X: = \{u \in W^{2, p}({\Omega}):\dfrac{ \partial u}{ \partial\bf{n}} = 0\}$, $Y = L^p({\Omega}).$

In the sequel, we let $(u^*, v^*, w^*)$ be a fixed constant positive steady state solution of (1.2) such that $u^* > 0,$ $0 < v^* < k_1, 0 < w^* < k_2.$.

2.1.   Prelimilary

By linearizing (1.2) around $(u^*, v^*, w^*),$ we have an eigenvalue problem as follows:

The following result says that the eigenvalue problem (2.1) can be reduced into a sequence of matrix eigenvalue problems.

Lemma 2.1. Let $\{ {\lambda}_n\}$ be the sequence of eigenvalues of $- {\Delta}$ with $Neumann$ boundary conditions satisfying $0 = {\lambda}_0 < {\lambda}_1\leq {\lambda}_2\leq \cdots$ For each $n,$ let $y_n(x)$ be the corresponding eigenfunction for $\lambda_n.$ Define

Then

1. Let $\mu$ be a complex number. Then $\mu$ is an eigenvalue of (2.1) iff there exists some $n\in\mathbb{N}_0$ such that $\mu$ is an eigenvalue of $A_n.$

2. $(u^*, v^*, w^*)$ is locally asymptotically stable with respect to (1.1) iff for every $n\in \mathbb{N}$, all eigenvalues of $A_n$ have negative real part.

3. $(u^*, v^*, w^*)$ is unstable with respect to (1.1) iff there exists an $n\in\mathbb{N}$, such that $A_n$ has at least one eigenvalue with nonnegative real part.

Proof. The proof of the first assertion can be done by direct computations. Here we omit the details. The rest two assertions follow from the principle of the linearized stability ^[5,10].

To continue, we use a direct calculation to find that the characteristic polynomial for each $A_n$ has the form

where

where

We impose the following conditions for the given constant equilibrium $(u^*, v^*, w^*):$

Under (2.5) we see that $a_2({\alpha}, {\lambda}_n) > 0$ for any $n\in \mathbb{N}$. Hence, as a result of applying the $Routh-Hurwitz$ criterion (see ^[8]), we have:

Corollary 2.2. The following assertions are true.

1. $(u^*, v^*, w^*)$ is locally asymptotically stable for (1.1) iff for each $n\in\mathbb{N}$ there holds that $(\mathit{\text{S}} 1)\; a_0({\alpha}, {\lambda}_n) > 0, \mathit{\text{and}} \, \, a_2({\alpha}, {\lambda}_n)a_1({\alpha}, {\lambda}_n)-a_0({\alpha}, {\lambda}_n) > 0.$

2. $(u^*, v^*, w^*)$ is unstable for (1.1) iff there exists an $n\in\mathbb{N}$ such that $(\mathit{\text{S}} 2)\; a_0({\alpha}, {\lambda}_n)\leq 0, \mathit{\text{or}} \, \, a_2({\alpha}, {\lambda}_n)a_1({\alpha}, {\lambda}_n)\leq a_0({\alpha}, {\lambda}_n).$

We now investigate the boundary between the stability and instability regimes

Let

be the steady state bifurcation curve, and

be the $Hopf$ bifurcation curve (see ^[14]). We study the steady state bifurcation form the given constant equilibrium $(u^*, v^*, w^*)$.

Since $a_0({\alpha}, p)$ is linear for ${\alpha}$, we solve ${\alpha}$ from the equation $a_0({\alpha}, p) = 0$ and obtain that ${\alpha} = {\alpha}_S(p)$ is given by

Correspondingly, we let ${\alpha}_H(p)$ be the solution of $T({\alpha}, p) = 0$, i.e., ${\alpha}_H(p)$ is the graph of function about $H$.

It is observed that the function ${\alpha}_S(p)$ has the following properties:

Lemma 2.3. Assume (2.5). Let ${\alpha}_S(p)$ be defined by (2.6). If $p^* > 0$ is a critical point of ${\alpha}_S(p),$ then $p^*$ is a local maximum point. Moreover, $\lim \limits_{p\rightarrow \infty} {\alpha}_S(p) = -\infty.$

Proof. Differentiating (2.6), we obtain

Assume $p > 0$ to be a critical point of ${\alpha}_S,$ i.e., ${\alpha}_S'(p) = 0.$ Then we have that

and

We find by (2.5) that all trems in the above bracket are positive, and $u^* {\Phi(v^*)} > 0.$ This implies that ${ {\alpha}_S}''(p) < 0$ and thus the critical point $p > 0$ is a local maximum point.

Clearly, we see from (2.5) that $\lim\limits_{p\rightarrow\infty} {\alpha}_S(p) = -\infty,$ since $u^* > 0$ and ${\Phi(v^*)} > 0.$

2.2.   Steady State Bifurcation

In this section, we inverstigate the global steady state bifurcation from $(u^*, v^*, w^*)$ near $\alpha = \alpha_{S}$. According to Lemma 2.1, we have the following result.

Proposition 2.4. Let $\{ {\lambda}_n\}$ be the sequence of eigenvalues of $- {\Delta}$ with $Neumann$ boundary conditions, such that $0 = {\lambda} _0 < {\lambda} _1\leq {\lambda} _2\leq \cdots$, let $y_n(x)$ be the eigenfunction corresponding to ${\lambda}_n$ ($n\in\mathbb{N}$). Let ${\alpha}_S(p)$ be defined by (2.6). For $n\in \mathbb{N}$ we define

Then the eigenvalue problem (2.1) has an eigenvalue $\mu = 0$ if and only if ${\alpha} = {\alpha}^S_n$ for some $n\in \mathbb{N}$, and the corresponding eigenfunction is $V_ny_n$, where $V_n$ satisfies $A_nV_n = 0$ with $A_n$ defined as in (2.2).

We recall the following global bifurcation theorem (see ^[8,11]):

Lemma 2.5. Let $V$ be an open connected subset of ${\mathbb R} \times X$ and $(\lambda_0, u_0)\in V$, and let $F$ be a continuously differentiable mapping from $V$ into $Y$. Suppose that

(1) $F(\lambda, u_0) = 0$ for $(\lambda, u_0)\in V$,

(2) the partial derivative $D_{\lambda u}F(\lambda, u)$ exists and is continuous in $(\lambda, u)$ near $(\lambda_0, u_0)$,

(3) $D_uF(\lambda_0, u_0)$ is a $Fredholm$ operator with index $0$, and $\dim \mathcal{N}(D_uF(\lambda_0, u_0)) = 1$,

(4) $D_{\lambda}(D_uF(\lambda_0, u_0))[w_0]\notin \mathcal{R}(D_uF(\lambda_0, u_0))$, where $w_0\in X$ spans $\mathcal{N}(D_uF(\lambda_0, u_0))$.

Let $Z$ be any complement of $span\{w_0\}$ in $X$. Then there exists an open interval $I_1 = (-\epsilon, \epsilon)$ and continuous functions $\lambda:I_1\rightarrow R, \; \psi:I_1\rightarrow Z$, such that $\lambda(0) = \lambda_0, \; \psi(0) = 0$, and, if $\; u(s) = u_0+sw_0+s\psi(s)\;$for $s\in I_1$, then $F(\lambda(s), u(s)) = 0$. Moreover, $F^{-1}(\{0\})$ near $(\lambda_0, u_0)$ consists precisely of the cures $u = u_0$ and $\Gamma = \left\{(\lambda(s), u(s)):\; s\in I_1\right\}$. If in addition, $D_uF(\lambda, u)$ is a $Fredholm$ operator for all $(\lambda, u)\in V$, then the curve $\Gamma$ is contained in $\mathcal{C}$, which is a connected component of $\bar S$ where $S = \left\{(\lambda, u)\in V:\; F(\lambda, u) = 0, \; u\neq u_0\right\}$; and either $\mathcal{C}$ is not compact in $V$, or $\mathcal{C}$ contains a point $(\lambda_*, u_0)$ with $\lambda_*\neq \lambda_0$.

We obtain the result for the global bifurcation of steady state solutions in predator prey taxis system (1.1) as follows.

Theorem 2.6. Assume (2.5) as well as ${\beta} > 0.$ Let ${\alpha}^S_n$ be defined as (2.7). Moreover, we assume that the following conditions (A1)-(A2) hold true:

$(A1)\;$For some $j\in \mathbb{N}$, ${\lambda}_j$ is a simple eigenvalue of $- {\Delta}$ in ${\Omega}$ with $Neumann$ boundary conditions, and the corresponding eigenfunction is $y_j(x)$.

$(A2)\;$For any $n\in \mathbb{N}$, $H({\alpha}_j^S, {\lambda}_n)\ne 0$, and if $n\neq j$, then ${\alpha}_j^S\neq {\alpha}_n^S$.

Then there hold following assertions:

1. The system (1.2) has a unique one-parameter family $\Gamma_j = \{(\hat{U}_j(s), \hat{ {\alpha}}_j(s)):-\varepsilon < s < \varepsilon\}$ of nontrivial solutions near $(u, v, w, {\alpha}) = (u^*, v^*, w^*, {\alpha}_j^S)$. More precisely, there exists $\varepsilon > 0$ and $C^\infty$ function $s\mapsto (\hat{U}_j(s), \hat{ {\alpha}}_j(s))$ from $s\in(-\varepsilon, \varepsilon)$ to $X^3\times {\mathbb R}$ satisfying

and

such that $h_{1, j}(0) = h_{2, j}(0) = h_{3, j}(0) = 0$;

2. The set $\Gamma_j$ is a subset of a connected component $\mathcal{C}_j$ of $\bar{S}$, where

and either $\mathcal{C}_j$contains another point $(u^*, v^*, w^*, {\alpha}_k^S)$ with ${\alpha}_k^S\neq {\alpha}_j^S$ or $\mathcal{C}_j$ is unbounded.

Proof. We define a mapping $F:X^3\times {\mathbb R}\rightarrow Y_0\times Y^2$ by

We have that $F(u^*, v^*, w^*, {\alpha}_j^S) = 0$, and $F$ is continuously differentiable. We will prove our result by applying Lemma 2.5 to $F.$

To this end, we must check that $F$ satisfies all requirements of Lemma 2.5. It will be completed in several steps.

$(1)$ For each $U = (u, v, w)$ the derivative $F_U(u^*, v^*, w^*, {\alpha}_j^S)$ is a $Fredholm$ operator with index zero (see ^[15]), and the kernel space $N(F_U(u^*, v^*, w^*, {\alpha}_j^S))$ is one-dimensional.

It left to show $N(F_U(u^*, v^*, w^*, {\alpha}_j^S))\neq\{0\}.$ We note that

Here $A, B$ are given by (2.4). Let $(\phi, \psi, \varphi)(\neq0)\in F_U(u^*, v^*, w^*, {\alpha}_j^S)$, then from Lemma 2.1 there exists $j\in \mathbb{N}$, such that $0$ is an eigenvalue of $A_j$, and the corresponding eigenvector is

According to condition $(A1)$, the eigenvector is unique up to a constant multiple. Hence we have

that is $\dim(N(F_U(u^*, v^*, w^*, {\alpha}_j^S))) = 1$.

$(2)$ $F_{ {\alpha} U}(u^*, v^*, w^*, {\alpha}_j^S)[(a_j^*, b_j^*, c_j^*)y_j]\not \in R(F_U(u^*, v^*, w^*, {\alpha}_j^S))$.

Define

where $(\bar{a_j}, \bar{b_j}, \bar{c_j})$ is a non-zero eigenvector fo the eigenvalue $\mu = 0$ of $A_n^T$, $A_n^T$ the transpose of $A_n$ and $(\bar{a_j}, \bar{b_j}, \bar{c_j}): = (a_j^*, b_j^*, c_j^*)y_j.$ Furthermore, if $(h_1, h_2, h_3, r)\in R(F_U(u^*, v^*, w^*, {\alpha}_j^S))$, then there exists $(\phi_1, \psi_1, \varphi_1)\in X^3$, such that

Define

and its adjoint operator

Then

where $\langle\cdot, \cdot\rangle$ is the inner product in $[L^2({\Omega})]^3$. If

then

Since (2.9) defines a codimension-1 set in $Y_0\times Y^2\times {\mathbb R}$, we get that

Therefore, $R(F_U(u^*, v^*, w^*, {\alpha}_j^S))$ is given by (2.8).

Note that

Therefore,

Here we have used the condition $A > 0$ from (2.5). Thus $F_{ {\alpha} U}(u^*, v^*, w^*, {\alpha}_j^S)[(a_j^*, b_j^*, c_j^*)y_j]\not \in R(F_U(u^*, v^*, w^*, {\alpha}_j^S))$.

Finally, we can apply the argument in Lemma 2.3 of ^[15] to obtain that the operator $F_U(u^*, v^*, w^*, {\alpha}_j^S)$ is a $Fredholm$ operator with index zero for any $(u, v, w, {\lambda})\in X^3\times {\mathbb R}$. Hence, all conditions in Lemma 2.5 are satisfied, and we get that solutions bifurcating from $(u^*, v^*, w^*, {\alpha}_j^S)$ are on a connected component $\mathcal{C}_j$ of the set of nontrivial solutions of (1.2). We see that all solutions on $\mathcal{C}_j$ are positive. This is apparently true for solutions near the bifurcation point $(u^*, v^*, w^*, {\alpha}_j^S)$ by $u^* > 0, v^* > 0, w^* > 0$. From the equation of $v$ and $w$ in (1.2), we know that if $u$ is nonnegative then there exists $(u, v, w, {\lambda})\in \mathcal{C}_j$, such that $v(x) > 0, w(x) > 0$ in $\bar{ {\Omega}}$, and when $x\in\bar{ {\Omega}}$ we have $u(x) = 0$. But from the equation of $u$ in (1.2), which is linear about $u$, we see that this is a contradiction by the strong maximum principle. Thus all solution are positive on $\mathcal{C}_j$.

2.3.   Numerical Simulations

Some numerical simulations of (1.1) are shown in this section, here the functional response $\Phi(v) = \dfrac{m_1 v}{a_1+v}$, $\Psi(w) = \dfrac{m_2 v}{a_2+v}$ are taken as Holling type Ⅱ, and $f(v) = v(k_1-v)$ and $g(w) = w(k_2-w)$ are taken as Logistic growth.

For simplicity, we take $a_1 = a_2 = a,$ $k_1 = k_2 = k$ and $m_1 = m_2 = m.$ Then the Eqn. (1.3) for a solution $(u, v, w)$ with $u > 0$ becomes

Thus, we assume that

Then, under condition (2.11), we obtain solution of (2.10) and thus constant steady state $(u^*, v^*, w^*)$ of (1.2) as follows:

In the following, we fix

Then $m$ satisfies condition (2.11) iff

Moreover, we take $\Omega = (0, 30\pi)$ (one-dimensional space).

$1.$ We take $m_1 = 0.35, m_2 = 0.35$ so that $m: = 0.35$ satisfies condition (2.12). Then the constant steady state $(u^*, v^*, w^*)$ is locally asymptotically stable for $\alpha = \beta = 0$, see Figure 1, where the initial value is $(u_0, v_{0}, w_{0}) = (3.3+0.01\sin(10x); 1.4+0.02\sin(x); 1.4+0.02\sin(x))$. From Theorem 2.6, the constant steady state $(u^*, v^*, w^*)$ becomes unstable where a steady state bifurcates from it for $\alpha = 100, \beta = 1$, see Figure 2.

$2.$ We further find that the taxis coefficient $\alpha$ and $\beta$ play a stability role for (1.1). In fact, we take $m_1 = 0.43, m_2 = 0.43$ so that $m: = 0.43$ satisfies condition (2.12). Then there is a periodic solution bifurcating from the positive constant steady state $(u^*, v^*, w^*)$ for $\alpha = \beta = 0$, see Figure 3, where the initial value is $(u_0, v_{0}, w_{0}) = (3.9, 0.75, 0.75)$. However, the positive large taxis coefficients $\alpha = 1\times 10^5$ and $\beta = 1\times10^4$ hold back the oscillation and make $(u^*, v^*, w^*)$ more stable, see Figure 4, where the initial value is $(u_0, v_{0}, w_{0}) = (2, 1, 1)$.

We might expect such an ecosystem (1.1) to exhibit a rich dynamical interplay among the three species. In this paper, we show that this is indeed the case, see Figure 1–Figure 4.

3.   Conclusion

We investigate the steady state bifurcation analysis in a general Ronsenzwing-MacArthur predator prey model with two prey-taxis under Neumann boundary conditions. Comparing with the dynamic analysis of reaction diffusion predator prey systems without taxis, the bifurcation analysis become more complicated and the bifurcation results in this paper cover most of reaction diffusion taxis predator prey systems. The results show that the rich dynamics in predator prey systems with two prey taxis and we will study other properties introduced by two taxis term in the future.

4.   Acknowledgements

We deeply thank three anonymous reviewers for their kind helps and valuable comments. We thank also Prof. Sen-Zhong Huang for valuable discussions. The work partially is supported by grants from NSFC(No. 11601105 and No. 11471091), NSFH(No. A2017007), GBD1317047, and HUDF2019101.

5.   Conflict of interest

The authors declare there is no conflict of interest.