A diffusive predator-prey model with generalist predator and time delay

1.   Introduction

Predator-prey model mainly describes the interaction between two populations with predation relationship. Since predator-prey relationship exists widely in nature, many scholars have studied the predator-prey models ^[1,2,3,4]. Considering the influence of different factors on the population, a variety of predator-prey models have been established ^[5,6]. Among these predator-prey models, Leslie-Gower predator-prey model is one of the classical model ^[7], with the following form

$u(t)$ and $v(t)$ stand for prey and predator's densities. $r_1$ and $K_1$ stand for the growth rate and the carrying capacity of the prey. $\varphi(u, v)$ is the functional response. The predator also follows the logistic growth law, where $r_2$ and $u/\beta$ stand for the growth rate of predator and the carrying capacity of the predator.

Another classical predator-prey model is Gauss predator-prey model ^[8], with the form

$c$ and $d$ are the conversion rate and death rate.

In predator-prey model, predators are mainly divided into specialist predators and generalist predators. Specialist predators feed almost exclusively on one specie of prey and require more specific environmental conditions. But, the generalist predators feed on many types of species, and can change its diet to another species when its a focal prey population begin to run short ^[9,10,11]. In ^[10], the authors studied a diffusive predator-prey model with generalist predator. They aimed to formalize the conditions in which spatial biological control can be achieved by generalists ^[10]. In ^[11], the authors studied the spatiotemporal dynamics and bifurcations of a diffusive predator-prey model with generalist predator and the combined the effect of linear prey harvesting and constant proportion of prey refuge. According to ^[10], the predator-prey model with generalist predators is of the following form

$K_2$ stands for the carrying capacity of the predator in absence of focal prey.

Predator-prey models with different functional responses can show different dynamic behaviors. In ^[10], the authors used the Type Ⅱ functional response to reflect the effect of predator to the prey. Holling Type Ⅱ functional response is a kind of prey-dependent functional response. Predator-dependent response function is also important. Such as Beddington-DeAngelis type ^[12], with the following form

where $B$, $C$, $A_1$ and $A_2$ stand for the maximum predator attack rate, the half-saturation constant, the effect of handling time and the magnitude of interference among predators. Some works all suggest that Beddington-DeAngelis functional response can enrich the dynamics of predator-prey model ^[13,14,15].

Moreover, time delay exists widely in the population model. The delayed predator-prey model has attracted wide attention from scholars ^[16,17,18]. These results show that time delay can enrich the dynamics of predator-prey model. Considering the generalist predator and discrete time delay $\tau$ in the resource limitation of the prey, we propose the following predator-prey model.

All the parameters in the model are positive. $u(x, t)$ and $v(x, t)$ stand for the densities of prey and predator at the location $x$ and time $t$, respectively. $E$ is conversion rate of prey. The boundary condition is Neumann boundary condition. The aim of this article is to study the dynamics of model (1.4) from the point of view of stability and bifurcation. Whether time delay can induce some new dynamic phenomena?

The organization of this paper is as follows. In Sec. 2, the existence of coexisting equilibrium of the model is given. In Sec. 3, stability of equilibria and the existence of Hopf bifurcation is considered. In Sec. 4, the property of Hopf bifurcation is analyzed. In Sec. 5, some numerical simulations are carried. In Sec. 6, a brief conclusion is given.

2.   Equilibrium analysis

Denote $\tilde{u} = u/K_1$, $\tilde{v} = u/K_2$ and $\tilde{t} = tr_2$, system (1.4) is changed to (after dropping tildes):

where

We assume $\Omega = (0, l\pi)$, where $l > 0$.

Solving the following equations,

We can obtained that $(0, 0)$, $(1, 0)$ and $(0, 1)$ are three boundary equilibria. And the coexisting equilibrium $(u_*, v_*)$ satisfying $v_* = \frac{(1-u_*) (1+b u_*)}{a-c+c u_*}$. Obviously, $v_* > 0$ implies $\mbox{max}\{0, \frac{c-a}{c}\} < u_* < 1$. In addition, form (2.3), we can easily obtained that

Submitting $v_* = \frac{(1-u_*) (1+b u_*)}{a-c+c u_*}$ into it, yields $h(u_*) = 0$, where

Theorem 2.1. If $c > a-1$, system $(2.1)$ has at least one coexisting equilibrium $(u_*, v_*)$, where $u_*$ is the root of $h(u_*) = 0$ in interval $(\mathit{\mbox{max}}\{0, \frac{c-a}{c}\}, 1)$ and $v_* = \frac{(1-u_*) (1+b u_*)}{a-c+c u_*}$.

Proof. By direct calculation, we have $h(1) = a^2 (1 + b + e) > 0$, $h(\frac{c-a}{c}) = -\frac{a^2 (a b-bc-c)^2}{c^3}\leq0$ and $h(0) = -a (1-a+c)$. If $c > a$, then $\frac{c-a}{c} = \mbox{max}\{0, \frac{c-a}{c}\}$ and $h(\frac{c-a}{c}) < 0$. Then $h(u_*) = 0$ has at least one root in interval $(\frac{c-a}{c}, 1)$. If $a-1 < c\leq a$, then $0 = \mbox{max}\{0, \frac{c-a}{c}\}$ and $h(0) < 0$. Then $h(u_*) = 0$ has at least one root in interval $(0, 1)$. This completes the proof.

3.   Stability analysis

Linearize system (2.1) at $(u_*, v_*)$

where

and

The characteristic equation of (3.1) is

where $I = \mbox{diag}\{1, 1\}$ and $M_n = -n^2/l^2\mbox{diag}\{d_1, d_2\}$, $n \in \mathbb{N}_0$. Then, we have

where

3.1.   The non-delay model

When $\tau = 0$, the characteristic Eq (2.1) is

where

and the eigenvalues are given by

We make the following hypothesis

Proposition 3.1. If hypothesis ${\bf(H_1)}$ (or ${\bf(H_2)}$) holds, then $a_1-u_* < 0$.

Proof. It is easy to obtain that $a_1-u_* = -\frac{u_* \phi(u_*)}{a (1+b u_*)}$, where $\phi(u) = b c u^2+2 b(a-c) u+a-a b+b c$. By direct calculation, we have $\phi(1) = a (1+b) > 0$, $\phi(\frac{c-a}{c}) = a \left(1+b-\frac{a b}{c}\right)$, $\phi(0) = a-a b+b c$ and $\phi^{'}(\frac{c-a}{c}) = 0$. From the proof of Theorem 2.1, we can obtain that $u_*\in (\frac{c-a}{c}, 1)$ under hypothesis ${\bf(H_1)}$. And $\phi(\frac{c-a}{c}) > 0$, implying that $\phi(u_*) > 0$. Hence $a_1-u_* < 0$ under hypothesis ${\bf(H_1)}$. Similarly, we can verify that $a_1-u_* < 0$ under hypothesis ${\bf(H_2)}$.

Theorem 3.1. Suppose ${\bf(H_1)}$ (or ${\bf(H_2)}$) holds. Then the equilibrium $E_*(u_*, v_*)$ is locally asymptotically stable.

Proof. By the Proposition 3.1, we know that $a_1-u_* < 0$ under hypothesis ${\bf (H_1)}$ (or ${\bf (H_1)}$). Then we have $tr_n < 0$ and $\Delta_n > 0$ for $n \in \mathbb{N}_0$. This implies that all eigenvalues of (3.5) have negative real parts. Then the equilibrium $E_*(u_*, v_*)$ is locally asymptotically stable.

3.2.   The delay model

Now, we study the stability of $E_*(u_*, v_*)$ when $\tau > 0$. Let $i\omega$ ($\omega > 0$) be a solution of Eq (3.4), we have

Then we have

This leads to

Denote $z = \omega^2$, then (3.9) can be changed into

and the roots of (3.10) are

Under ${\bf (H_1)}$ (or ${\bf (H_2)}$), we have

By direct computation,

Fix parameters $r, \; a, \; b, \; c, \; e, \; d_1, \; d_2, \; l$, define

For $n\in\mathcal{D}$, if $z^{+} > 0$, Eq (3.4) has a pair of purely imaginary roots $\pm i\omega^+_n$ at $\tau^{j, +}_n$, $j \in \mathbb{N}_0$; if $z^{-} > 0$, Eq (3.4) has a pair of purely imaginary roots $\pm i\omega^-_n$ at $\tau^{j, -}_n$, $j \in \mathbb{N}_0$, where

From (3.12), we have $\tau^{0, \pm}_n < \tau^{j, \pm}_n$ $(j\in \mathbb{N})$. For $k\in \mathcal{D}$, define the smallest $\tau$ so that the stability will change, $\tau_* = \mbox{min}\{\tau^{0, \pm}_k \; \mbox{or}\; \tau^{0, +}_k \mid k \in \mathcal{D} \}$.

Lemma 3.1. Suppose ${\bf{(H_1)}}$ (or ${\bf (H_2)}$) holds.If $(A_n^2-2 B_n- r^2 u_*^2)^2-4(B_n^2-C_n^2) > 0$, then $\mathit{\text{Re}}(\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, +}_n} > 0$, $\mathit{\text{Re}} (\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, -}_n} < 0$ for $\tau \in \mathcal{D}$ and $j \in \mathbb{N}_0$.

Proof. Differentiating two sides of (3.4) with respect $\tau$, we have

Then

where $\Lambda = \omega^4 b_2^2 +C^2_n \omega^2 > 0$. Therefore $\text{Re}(\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, +}_n} > 0$, $\text{Re} (\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, -}_n} < 0$.

Theorem 3.2. Suppose ${\bf{(H_1)}}$ (or ${\bf (H_2)}$) holds. For system $(2.1)$, the following statements are true.

(i) $E_*(u_*, v_*)$ is locally asymptotically stable for all $\tau \geq0$ when $\mathbb{D} = \emptyset$.

(ii) $E_*(u_*, v_*)$ is locally asymptotically stable for $\tau \in [0, \tau_*)$, and unstable for $\tau \in [\tau_*, \tau_*+\epsilon)$ with some $\epsilon$ when $\mathcal{D}\neq \emptyset$.

(iii) System $(2.1)$ undergoes a Hopf bifurcation at the equilibrium$E_*(u_*, v_*)$ when $\tau = \tau^{j, +}_n$ $(or \; \tau = \tau^{j, -}_n)$, $j\in \mathbb{N}_0$, $n \in \mathcal{D}$ when $\mathcal{D}\neq \emptyset$.

Using the same process, we can obtain the following theorem about the stability of boundary equilibria.

Theorem 3.3. For system $(2.1)$, the following statements are true.

(i) $(0, 0)$ is always unstable for all $\tau \geq0$;

(ii) $(1, 0)$ is always unstable when $\tau = 0$;

(iii) $(0, 1)$ is locally asymptotically stable for $a > 1+c$ and $\tau \geq0$; and unstable for $a < 1+c$ and $\tau \geq0$.

4.   Properties of Hopf bifurcation

Now, we will study the property of Hopf bifurcation by the method in ^[19,20]. For a critical value $\tau^{j, +}_n$ (or $\tau^{j, -}_n)$, we denote it as $\tilde{\tau}$. Let $\tilde{u}(x, t) = u(x, \tau t)-u_*$ and $\tilde{v}(x, t) = v(x, \tau t)-v_*$, then the system (2.1) is (drop the tildes)

Denote $\tau = \tilde{\tau}+\varepsilon,$ and $U = (u(x, t), v(x, t))^T$. In the phase space $\mathbb{C}_1: = C([-1, 0], X)$, (4.1) can be rewritten as

where $L_{\varepsilon}(\varphi)$ and $F(\varphi, \varepsilon)$ are

with

respectively, for $\phi = (\phi_1, \phi_2)^T \in \mathbb{C}_1$.

Consider the linear equation

We know that $\Lambda_n: = \{i \omega_n \tilde{\tau}, -i \omega_n \tilde{\tau}\}$ are characteristic roots of

By Riesz representation theorem, there exists a $2\times2$ matrix function $\eta^n(\sigma, \tilde{\tau})$, $(-1\le \sigma \le0)$, whose elements are of bounded variation functions such that

for $\phi \in C([-1, 0], \mathbb{R}^2)$.

Choose

where

Define the bilinear paring

for $\varphi \in C([-1, 0], \mathbb{R}^2)$, $\psi \in C([0, 1], \mathbb{R}^2)$. $A(\tilde{\tau})$ has a pair of simple purely imaginary eigenvalues $\pm i \omega_n \tilde{\tau}$, and they are also eigenvalues of $A^*$.

Define $p_1(\theta) = (1, \xi)^Te^{i\omega_n \tilde{\tau} \sigma}\; \; (\sigma \in[-1, 0]), \; \; q_1(r) = (1, \eta)e^{-i\omega_n \tilde{\tau} r}\; \; (r \in[0, 1]),$ where

Let $\Phi = (\Phi_1, \Phi_2)$ and $\Psi^* = (\Psi^*_1, \Psi^*_2)^T$ with

for $\theta \in[-1, 0]$, and

for $r \in[0, 1]$. Then we can compute by (4.9)

Define $(\Psi^*, \Phi) = (\Psi^* _j, \Phi _k) = \left(\begin{array}{c} D^*_1\; \; D^*_2\\ D^*_3\; \; D^*_4 \end{array} \right)$ and construct a new basis $\Psi$ for $P^*$ by

Then $(\Psi, \Phi) = I_2$. In addition, define $f_n: = (\beta^1_n, \beta^2_n)$, where

We also define

and

for $u = (u_1, u_2)$, $v = (v_1, v_2)$, $u, v\in X$ and $< \varphi, f_0 > = (<\varphi, f^1_0 > , < \varphi, f^2_0 > )^T$.

Rewrite Eq (4.1) as the following abstract form

where

The solution is

where

and

Then

Let $z = x_1-i x_2$, and notice that $p_1 = \Phi_1+i\Phi_2$. Then we have

and

Equation (4.13) is

where

From ^[19], $z$ satisfies

where

Let

then

and

with $\alpha_1 = \frac{a r v_* (b+b c v_*)}{(1+b u_*+c v_*)^3}$, $\alpha_2 = -\frac{a r (1+c v_*+b (u_*+2 c u_* v_*))}{(1+b u_*+c v_*)^3}$, $\alpha_3 = \frac{a r u_* (c+b c u_*)}{(1+b u_*+c v_*)^3}$, $\alpha_4 = -\frac{a b^2 r v_* (1+c v_*)}{(1+b u_*+c v_*)^4}$, $\alpha_5 = \frac{a b r \left(1-c^2 v_*^2+b (u_*+2 c u_* v_*)\right)}{(1+b u_*+c v_*)^4}$, $\alpha_6 = \frac{a c r \left(1-b^2 u_*^2+c v_*+2 b c u_* v_*\right)}{(1+b u_*+c v_*)^4}$, $\alpha_7 = -\frac{6 a c^2 r u_* (1+b u_*)}{(1+b u_*+c v_*)^4}$, $\beta_1 = -\frac{ e v_* (b+b c v_*)}{(1+b u_*+c v_*)^3}$, $\beta_2 = \frac{e (1+c v_*+b (u_*+2 c u_* v_*))}{(1+b u_*+c v_*)^3}$, $\beta_3 = -1-\frac{e u_* (c+b c u_*)}{(1+b u_*+c v_*)^3}$, $\beta_4 = \frac{b^2 e v_* (1+c v_*)}{(1+b u_*+c v_*)^4}$, $\beta_5 = -\frac{b e \left(1-c^2 v_*^2+b (u_*+2 c u_* v_*)\right)}{(1+b u_*+c v_*)^4}$, $\beta_6 = -\frac{c e \left(1-b^2 u_*^2+c v_*+2 b c u_* v_*\right)}{(1+b u_*+c v_*)^4}$, $\beta_7 = \frac{c^2 e u_* (1+b u_*)}{(1+b u_*+c v_*)^4}$.

Hence,

with

and

$\chi_{20} = \frac{1}{2} e^{-i \tau \omega _n} \left(-r+e^{i \tau \omega _n} (\alpha_1+\xi (\alpha_2+\alpha_3 \xi))\right)$,

$\chi_{11} = -\frac{1}{4} e^{-i \tau \omega _n} \left(\left(1+e^{2 i \tau \omega _n}\right) r-e^{i \tau \omega _n} (2 \alpha_1+2 \alpha_3 \bar{\xi} \xi +\alpha_2 (\bar{\xi} +\xi))\right)$,

$\chi_{1} = W_{11}^1(0) \left(-e^{-i \tau \omega _n} r+2 \alpha_1+\alpha_2 \xi \right)+W_{11}^2(0) (\alpha_2+2 \alpha_3 \xi)+\frac{1}{2} W_{20}^1(0) \left(-e^{i \tau \omega _n} r+2 \alpha_1 +\alpha_2 \bar{\xi} \right)$

$\; \; \; \; \; \; \; \; \; +\frac{1}{2} W_{20}^2(0) (\alpha_2+2 \alpha_3 \bar{\xi})-r W_{11}^1(-1)-\frac{r W_{20}^1(-1)}{2}$,

$\chi_{2} = \frac{1}{4} (3 \alpha_4+\alpha_5 (\bar{\xi} +2 \xi)+\xi (2 \alpha_6 \bar{\xi} +\alpha_6 \xi +3 \alpha_7 \bar{\xi} \xi))$,

$\varsigma_{20} = \frac{1}{2} (\beta_1+\xi (\beta_2+\beta_3 \xi))$, $\varsigma_{11} = \frac{1}{4} (2 \beta_1+2 \beta_3 \bar{\xi} \xi +\beta_2 (\bar{\xi} +\xi))$,

$\varsigma_{1} = W_{11}^1(0) (2 \beta_1+\beta_2 \xi)+W_{11}^2(0) (\beta_2+2 \beta_3 \xi)+W_{20}^1(0) \left(\beta_1+\frac{\beta_2 \bar{\xi} }{2}\right)+\frac{1}{2} W_{20}^2(0) (\beta_2+2 \beta_3 \bar{\xi})$,

$\varsigma_{2} = \frac{1}{4} (3 \beta_4+\beta_5 (\bar{\xi} +2 \xi)+\xi (2 \beta_6 \bar{\xi} +\beta_6 \xi +3 \beta_7 \bar{\xi} \xi))$.

Denote

Notice that

and we have

Then by (4.16), (4.18) and (4.23), we have $g_{20} = g_{11} = g_{02} = 0$, for $n = 1, 2, 3, \cdots$. If $n = 0$, we have:

And for $n\in\mathbb{N}_0$, $g_{21} = \tilde{\tau}(\gamma_1 \kappa_1 +\gamma_2 \kappa_2)$.

From ^[19], we have

and

where

Hence, we have

that is

Then

Therefore,

and

where

By the definition of $A_{\tilde{\tau}}$ and (4.25), we have

That is

where

By the definition of $A_{\tilde{\tau}}$ and (4.25), we have that for $-1\le \theta < 0$

As

and

we have

That is

where

Similarly, from (4.26), we have

That is

Similarly, we have

where

Thus, we have

Theorem 4.1. When $\mu_2 > 0$ (or $\mu_2 < 0$), the bifurcating periodic solutions exists for $\tau > \tau^{j, \pm}_n$ (or $\tau < \tau^{j, \pm}_n$), andare orbitally asymptotically stable(or unstable) when $\beta_2 < 0$ (or $\beta_2 > 0$).

5.   Numerical simulations

According to the reference ^[21], we choose $r_1 = 19.3$, $K_1 = 400$, $r_2 = 8.8$, $K_2 = 5$, $B = 10.76$, $C = 60.6$, $A_1 = 0.00728$, $A_2 = 1$. By (2.2), we can obtain $a$, $b$, $c$, $e$ and $r$. Fix $d_1 = 0.1$, $d_2 = 0.2$, $l = 2$. We give the coexisting equilibrium (Figure 1), and bifurcation diagram (Figure 2) of system (2.1) with parameter $E$. It shows that the density of prey (predator) in the coexisting equilibrium decreases (increases) and the stable region increases with the increase of parameter $E$. This indicates that increasing parameter $E$ is benefit to prey and predator to reach the steady state (coexisting equilibrium).

When $E = 8.8$, we can obtain $(u_*, v_*)\approx(0.0037, \; 1.2370)$ is a unique coexisting equilibrium. Hypothesis ${\bf (H_1)}$ is satisfied. We have $\tau_* = \tau^0_0 \approx 0.2440$. Then $E_*(u_*, v_*)$ is local stable when $\tau \in [0, \tau_*)$ (shown in Figure 3). When $\tau = \tau_*$, Hopf bifurcation occurs. We can obtain

Then, when $\tau > 3.4595$, the local stable bifurcating periodic solutions exists (shown in Figure 4).

6.   Conclusions

In this paper, we propose a diffusive predator-prey system with generalist predator and time delay in the resource limitation of the prey. We obtained that system (2.1) has three boundary equilibria: $(0, 0)$ (predator and prey extinction equilibrium), $(1, 0)$ (predator extinction equilibrium) and $(0, 1)$ (prey extinction equilibrium). We mainly analyze the stability and Hopf bifurcation of coexisting equilibrium. By the theory of normal form and center manifold method, we give some parameters that determining the property of Hopf bifurcation: Bifurcation direction and the stability of the bifurcating periodic solution.

Since the predators are generalist type and have other food resource, they will not be extinct. This is in agreement with the Theorem 3.3. When the predator attack rate is large enough $a > 1+c$, all the prey will be caught by the predator. This will lead to the extinction of the prey, and the predator will reach a balanced state. It is also in agreement with the Theorem 3.3 that $(0, 1)$ is local asymptotically stable for $a > 1+c$. When the predator attack rate is not large enough $a < 1+c$, then the prey and predator will coexist.

The conversion rate $E$ can affect the the density of prey (predator) in the coexisting equilibrium. With the increase of conversion rate, the density of prey (predator) will decrease (increase) and the stable region will increase. In addition, time delay will also affect the stability of the equilibrium point when the parameters satisfying the condition $\mathcal{D}\neq \emptyset$. Specifically, when the time delay is small than the critical value $\tau_*$, the prey and the predator will coexist and tend to the coexisting equilibrium. But when time delay is larger than the critical value $\tau_*$, the prey and predator will exhibit oscillatory behavior. In addition, the spatial inhomogeneous periodic solutions may exist, but they are generally unstable.

Acknowledgments

This research is supported by Fundamental Research Funds for the Central Universities (No. DL13BB17), Heilongjiang Provincial Natural Science Foundation (No. A2018001), Postdoctoral Science Foundation of China (No. 2019M651237) and National Nature Science Foundation of China (No. 11601070).

Conflict of interest

The authors declare that they have no competing interests.