Hopf bifurcation analysis in a delayed diffusive predator-prey system with nonlocal competition and generalist predator

1.   Introduction

Dynamic systems are widely used in nature, such as in infectious disease models ^[1,2], vegetation-water models ^[3], and population models ^[4,5]. For example, Khan et al. ^[6] studied the dynamics of pine wilt disease with variable population sizes and showed control strategies for the possible elimination of the infection in the pine tree population. Khan et al. ^[7] obtained an efficient iterated homotopy perturbation transform method (IHPTM) to solve a mathematical model of HIV infection of CD4+ T cells. Khan et al. ^[8] proposed a new method to solve partial differential equations arising in the fields of science and engineering. Among the applications of dynamical systems, the predator-prey model is an important research topic ^[9,10,11,12]. Generally, predator-prey models assume that the prey is the only food source of the predator ^[13,14,15].

In ecosystems, predators are usually generalist predators. They feed on many types of species and can change their diet to another species when a focal prey population begins to run short ^{[16,17,18,19]}. In ^[17], the authors studied a predator-prey model with a generalist predator. They aimed to achieve biological control by generalists ^[17]. In ^[18], the authors studied the spatial pattern formation of a predator-prey model with generalist predator and the harvesting, refuge effect. To better understand the relationship between prey and generalist predators, Upadhyay and Agrawal ^[19] modified the Leslie-Gower model as

to describe the growth law of generalist predators. The term $dv^2$ represents the growth of predators, where $d$ is the mating rate of the predator. $e$ is another source of food for predators. They considered the following model:

$u(t)$ and $v(t)$ represent the prey and predator densities, respectively. The term $\frac{\omega u}{A+Bu+v}$ represents is the Beddington-DeAngelis functional response. $\tau$ is the gestation delay of the predator. Upadhyay and Agrawal studied the invariance, boundedness, and local and global stability and Hopf bifurcation.

In the real world, the spatial distribution of populations is often inhomogeneous, so it is more practical to use a reaction-diffusion model to describe the relationship between predators and prey. In addition, the reaction-diffusion predator-prey model shows more abundant dynamic properties such as spatial patterns and spatial nonhomogeneous periodic solutions. Considering this factor, Liu et al. ^[20] proposed the following model with the Crowley-Martin type functional response based on the model (1.1).

In the model (1.2), Liu et al. ^[20] used the Crowley-Martin-type functional response to reflect the impact of predators on prey. This considered the effect of interference among predators. $a$, $b$, and $c$ represent the capture rate, handling time, and magnitude of interference among predators, respectively. They mainly analyzed the instability and Hopf bifurcation induced by time delay ^[20]. Although they pointed out that time delay may cause spatially inhomogeneous periodic solutions, the numerical simulations did not show stably spatially inhomogeneous periodic solutions. This is because in the delayed reaction-diffusion predator-prey model, the spatial nonhomogeneous Hopf bifurcation curve is usually above the spatial homogeneous Hopf bifurcation curve. In addition, Turing instability cannot occur for the model (1.2). For model (1.2), it is unfortunate that the inhomogeneous distribution of prey and predator in space is not shown. This may be due to the lack of stable spatially inhomogeneous periodic solutions in the delayed reaction-diffusion predator-prey model.

In nature, competition within populations exists widely, and this competition is often nonlocal since resources are limited. In ^[21,22], the authors suggested that the internal competition of the population caused by the natural environment is related not only to the population density at the current location but also to the population density nearby. They measured this effect by weighting and integrating, and modified the $\frac{u}{K}$ as $\frac{1}{K} \int_{\Omega}G(x, y)u(y, t)dy$. $G(x, y)$ is a kernel function. Wu and Song studied Hopf and state-Hopf bifurcations in a diffusive predator-prey model with a nonlocal effect and delay ^[23]. Geng et al. studied a diffusive predator-prey model with nonlocal competition, including Hopf, Turing, double-Hopf, and Turing-Hopf bifurcations ^[24]. The works in ^[25,26] show that stable spatially inhomogeneous periodic solutions often exist in predator-prey models with nonlocal competition.

Inspired by the above work, we assume there is nonlocal competition in prey and modify the model (1.3) as follows.

where $d_1$ and $d_2$ represent the diffusion coefficients of prey and predator, respectively. The term $\frac{a u}{(1+bu)(1+cv)}$ is the Crowley-Martin-type functional response. $\Omega$ is the prey and predator's living region. Just for the convenience of calculation, we choose $\Omega = (0, l\pi)$, where $l > 0$. The boundary condition is the Newman boundary condition, which means that the living region is closed and no prey and predator enter or leave the region. $\int_{\Omega}G(x, y)u(y, t)dy$ represents the nonlocal competition effect, and $G(x, y)$ is the kernel function. We assume that the competition strength among prey individuals in the habitat is the same; we choose $G(x, y) = \frac{1}{l\pi}$ as in previous works ^{[23,24,25,26]}.

To our knowledge, there is no work about the predator-prey model (1.3) with nonlocal competition in prey and generalist predators from the point of Hopf bifurcation. The aim of this paper is to study the combined effect of time delay and nonlocal competition on model (1.3). Compared with the model (1.2), do new dynamics appear, such as stably spatially inhomogeneous periodic solutions and Turing instability?

This paper is organized as follows: the stability and existence of Hopf bifurcation are studied in Section 2, the property of Hopf bifurcation is analyzed in Section 3, and numerical simulations are given in Section 4. Finally, a short conclusion is presented.

2.   Stability analysis

$(0, 0)$ and $(K, 0)$ are boundary equilibria of system (1.3). The existence of positive equilibria of system (1.3) was studied in ^[20], that is,

Lemma 2.1. ^[20]If the following condition

holds, then the model $(1.3)$ has a coexisting equilibrium point $(u_*, v_*)$, where

Next, we just denote $E_*(u_*, v_*)$ as a coexisting equilibrium. Linearize system (1.3) at $E_*(u_*, v_*)$

where

and

The characteristic equation is

where

2.1.   Nondelay model

When $\tau = 0$, the characteristic equation (2.3) is

where

Make the following hypothesis:

Theorem 2.2. For system $(1.3)$, assume $\tau = 0$ and ${\bf (H_0)}$ holds. Then, $E_*(u_*, v_*)$ is locally asymptotically stable under ${\bf (H_1)}$ and is Turing unstable under ${\bf (H_2)}$.

Proof. If ${\bf (H_1)}$ holds, then we can determine that the characteristic roots of (2.5) all have negative real parts. Then, $E_*(u_*, v_*)$ is locally asymptotically stable. If ${\bf (H_1)}$ holds, then the characteristic roots of (2.5) with $k\in \mathbb{N}$ have at least one positive real part, but with $n = 0$ all having a negative real part. This implies that $E_*(u_*, v_*)$ is Turing unstable.

Lemma 2.3. When $\tau = 0$, the Turing instability of $E_*(u_*, v_*)$ cannot occur for the system $(1.2)$, which lacks the nonlocal competition term.

Proof. When $\tau = 0$, for the system (1.3) without nonlocal competition, the characteristic equation (2.3) is

where

Assume $E_*(u_*, v_*)$ is locally asymptotically stable in the absence of spatial diffusion, which means that $a_1+\hat{a} < 0$. However, $T_n > 0$ and $D_n > 0$ when $a_1+\hat{a} < 0$ holds. Then, the Turing instability of $E_*(u_*, v_*)$ cannot occur.

2.2.   Delay model

Let $\mbox{i}\omega$ ($\omega > 0$) be a solution of Eq (2.3); then,

We can obtain

It leads to

Let $z = \omega^2$; then, (2.8) becomes

and the roots of (2.9) are

where $P_n = A_n^2-2 B_n-b_1^2$, $Q_n = B_n+C_n$, and $R_n = B_n-C_n$. If ${\bf (H_0)}$ and ${\bf (H_1)}$ hold, $Q_n > 0\; (n \in \mathbb{N}_0)$. By direct calculation, we have

Define

and

We have the following lemma.

Lemma 2.4. Assuming that ${\bf (H_0)}$ and ${\bf (H_1)}$ hold, the following results hold:

(1) Eq $(2.3)$ has a pair of purely imaginary roots$\pm i\omega^+_n$ at $\tau^{j, +}_n$ for $j \in \mathbb{N}_0$ and $n\in\mathbb{W}_1$.

(2) Eq $(2.3)$ has two pairs of purely imaginary roots$\pm i\omega^\pm_n$ at $\tau^{j, \pm}_n$ for $j \in \mathbb{N}_0$ and $n\in\mathbb{W}_2$.

(3) Eq $(2.3)$ has no purely imaginary root for $n\in\mathbb{W}_3$.

Lemma 2.5. Assume ${\bf (H_0)}$ and ${\bf (H_1)}$ hold. 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 $n \in \mathbb{W}_1\cup \mathbb{W}_2$ and $j \in \mathbb{N}_0$.

Proof. By (2.3), we have

Then

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$.

Denote $\tau_* = min\{\tau^0_n|\; n\in\mathbb{W}_1\cup \mathbb{W}_2\}$. We have the following theorem:

Theorem 2.6. Assume that ${\bf (H_0)}$ and ${\bf (H_1)}$ hold. Then, the following statements are true for system $(1.3) :$

(1) $E_*(u_*, v_*)$ is locally asymptotically stable for $\tau > 0$ when $\mathbb{W}_1\cup \mathbb{W}_2 = \varnothing$.

(2) $E_*(u_*, v_*)$ is locally asymptotically stable for $\tau \in [0, \tau_*)$ when $\mathbb{W}_1\cup \mathbb{W}_2\neq \varnothing$.

(3) $E_*(u_*, v_*)$ is unstable for $\tau \in (\tau_*, \tau_*+\varepsilon)$ for some $\varepsilon > 0$ when $\mathbb{W}_1\cup \mathbb{W}_2\neq \varnothing$.

(4) Hopf bifurcation occurs at$(u_*, v_*)$ when $\tau = \tau^{j, +}_n \; (\tau = \tau^{j, -}_n)$, $j\in \mathbb{N}_0$, $n \in \mathbb{W}_1\cup \mathbb{W}_2$.

3.   Property of Hopf bifurcation

In ^[27,28], we studied the property of Hopf bifurcation. For fixed $j\in \mathbb{N}_0$ and $n\in \mathbb{W}_1\cup \mathbb{W}_2$, we denote $\tilde{\tau} = \tau^{j, \pm}_n$. Let $\bar{u}(x, t) = u(x, \tau t)-u_*$ and $\bar{v}(x, t) = v(x, \tau t)-v_*$. Drop the bar, and (1.3) can be written as

We rewrite system (4) as follows:

where

Define the real-valued Sobolev space as

the complexification of $X \; X_{\mathbb{C}}: = X \oplus \mbox{i}X = \left\{x_1+\mbox{i}x_2|\; x_1, x_2\in X\right\}.$ and the inner product

for $\tilde{u} = (u_1, u_2)^{T}$, $\tilde{v} = (v_1, v_2)^{T}$, $\tilde{u}$, $\tilde{v}\in X_{\mathbb{C}}$. The phase space $\mathscr{C}: = C([-1, 0], X)$ has the sup norm. Then, we can write $\phi_t \in \mathscr{C}$, $\phi_t(\theta) = \phi(t+\theta)$ or $-1\leq \theta \leq0$. Denote $\beta_n^{(1)}(x) = (\gamma_n(x), 0)^T$, $\beta_n^{(2)}(x) = (0, \gamma_n(x))^T$, and $\beta_n = \{\beta_n^{(1)}(x), \beta_n^{(2)}(x) \}$, where $\{\beta_n^{(i)}(x) \}$ is an an orthonormal basis of $X$. We define the subspace of $\mathscr{C}$ as

There exists a $2\times2$ matrix function $\eta^n(\sigma, \tilde{\tau}) \; -1\le \sigma \le0$ such that

for $\phi \in \mathscr{C}$. The bilinear form on $\mathscr{C}^* \times \mathscr{C}$ is defined by

for $\phi \in \mathscr{C}$, $\psi \in \mathscr{C}^*$. Defining $\tau = \tilde{\tau}+\mu$, the system undergoes a Hopf bifurcation at $(0, 0)$ when $\mu = 0$, with a pair of purely imaginary roots $\pm \mbox{i}\omega_{n_0}$. Let $A$ denote the infinitesimal generators of the semigroup, and $A^*$ be the formal adjoint of $A$ under the bilinear from (3.3). Define the following function:

Choose $\eta_{n_0}(0, \tilde{\tau}) = \tilde{\tau}[(-n_0^2/l^2)D+L_1+L_3 \delta(n_{n_0})]$, $\eta_{n_0}(-1, \tilde{\tau}) = -\tilde{\tau} L_2$, $\eta_{n_0}(\sigma, \tilde{\tau}) = 0$ for $-1 < \sigma < 0$. Let

be the eigenfunctions of $A(\tilde{\tau})$ and $A^*$ corresponding to $i \omega_{n_0} \tilde{\tau}$, respectively. We can choose $p(0) = (1, p_1)^T$, $q(0) = M(1, q_2)$, where

and

Then, (4) can be rewritten in abstract form as

where

respectively, for $\phi = (\phi_1, \phi_2)^T \in \mathscr{C}$ and $\hat{\phi}_1 = \frac{1}{l\pi}\int_{0}^{l\pi}\phi dx$. Then, the space $\mathscr{C}$ can be decomposed as $\mathscr{C} = P\oplus Q$, where

Then, system (3.6) can be rewritten as

and $\hat{U_t} = \frac{1}{l\pi}\int_{0}^{l\pi} U_t dx$, where

Then, we have

There exists a center manifold $\mathcal{C}_0$, and $\omega$ can be written as follows near $(0, 0)$:

Then, restrict the system to the center manifold as $\dot{z}(t) = \mbox{i}\omega_{n_0} \tilde{\tau} z(t)+g(z, \bar{z})$. Denoting

By direct computation, we have

where $I_2 = \int_{0}^{l\pi} \gamma^2_{n_0}(x) dx$, $I_3 = \int_{0}^{l\pi} \gamma^3_{n_0}(x) dx$, $I_4 = \int_{0}^{l\pi} \gamma^4_{n_0}(x) dx$, and

Now, we compute $W_{20}(\theta)$ and $W_{11}(\theta)$ for $\theta \in [-1, 0]$ to give $g_{21}$. By (3.7), we have

where

Comparing the coefficients of (3.8) with (3.9), we have

Then, we have

where $E_1 = \sum_{n = 0}^{\infty}E_1^{(n)}$, $E_2 = \sum_{n = 0}^{\infty}E_2^{(n)}$,

and $\hat{F}_{20} = 2(\varsigma_1, \varsigma_2)^T$, $\hat{F}_{11} = 2(\varrho_1, \varrho_2)^T$.

Thus, we can obtain

Theorem 3.1. For any critical value $\tau^j_n$ ($n\in \mathbb{S}, \; j \in \mathbb{N}_0$), we have the following results:

(1) When $\mu_2 > 0$ (resp. $< 0$), the Hopf bifurcation is forward (resp. backward).

(2) When $\beta_2 < 0$ (resp. $> 0$), the bifurcating periodic solutions on the center manifold are orbitally asymptotically stable(resp. unstable).

(3) When $T_2 > 0$ (resp. $< 0$), the period increases (resp. decreases).

4.   Numerical simulations

To verify our theoretical results, we give the following numerical simulations. Fix parameters

The bifurcation diagram of system $(1.3)$ with the parameter of interference magnitude among predators $c$ is given in Figure 1, where ${\bf (H_0)}$ and ${\bf (H_1)}$ hold. We can see that with an increase in parameter $c$, the Hopf bifurcation curves decrease, which implies that the stable region of positive equilibrium $(u_*, v_*)$ will decrease. This means that increasing the interference magnitude among predators is not conducive to the homogeneous distribution of prey and predators in space. This causes spatial oscillations in prey and predator densities. Since the inhomogeneous Hopf bifurcation curve $\tau_1^0$ is always under the homogeneous Hopf bifurcation curve $\tau_0^0$. This spatial oscillation of prey and predator density is inhomogeneous.

Fixing $c = 0.1$, we can determine that $(u_*, v_*)\approx (9.0000, 1.9048)$ is the unique positive equilibrium, and ${\bf (H_1)}$ holds. By direct calculation, we have $\tau_* = \tau_1^0\approx 0.0758 < \tau_0^0\approx 0.1002$, $\mu_2\approx 0.0105$, $\beta_2 \approx -0.0584$ and $T_2\approx 0.8602$. If we choose $\tau = 0.05 < \tau_*$, then $(u_*, v_*)$ is locally asymptotically stable for model (1.3) (shown in Figure 2) and model (1.2) (shown in Figure 3). This means that the nonlocal competition term has no effect on the stability of the model (1.3). Prey and predator will coexist in a spatially homogeneous form, and their densities will converge to the positive equilibrium $(u_*, v_*)$.

If we choose $\tau_* < \tau = 0.09 < \tau_0^0$, then $(u_*, v_*)$ is unstable, and inhomogeneous periodic solutions exist for the model (1.3) (shown in Figure 4). This means that prey and predator will coexist in the form of spatially inhomogeneous oscillations, and their densities will distribute inhomogeneously in space. To compare our result with the work in ^[20], we give the numerical simulations of model $(1.2)$ under the same parameter $\tau = 0.09$ in Figure 5. We can see that $(u_*, v_*)$ is locally asymptotically stable $(1.2)$. This means that prey and predators will still coexist in a spatially homogeneous form and that their densities will converge to positive equilibrium $(u_*, v_*)$.

5.   Conclusions and discussion

Assuming there is nonlocal competition caused by limited resources in prey, we modified the model (1.2) to the model (1.3) with the nonlocal competition in prey, gestation delay in predator, and generalist predator. We study the local stability of the positive equilibrium and Hopf bifurcation at the positive equilibrium by using the time delay as a parameter. Through the center manifold and normal form method, we give some parameters that can determine the direction of Hopf bifurcation, the stability of bifurcating periodic solutions, and the period of periodic solutions.

For the time delay when $\tau = 0$, we determine that nonlocal competition can induce Turing instability of the coexisting equilibrium point $(u_*, v_*)$, but Turing instability cannot occur in the model (1.2) without nonlocal competition. Similar to the results in ^[20], time delays can affect the stability of the positive equilibrium. The positive equilibrium is stable when the delay is shorter than the critical value but unstable when the delay is longer than the critical value, and prey and predator densities will produce periodic oscillations.

Through the bifurcation diagram, we find that increasing the interference magnitude among predators is not conducive to the homogeneous distribution of prey and predators in space. This causes spatial oscillations in prey and predator densities. Compared with the work in ^[20], the numerical simulations show that the inhomogeneous Hopf bifurcation curve $\tau_1^0$ is always under the homogeneous Hopf bifurcation curve $\tau_0^0$ for the model (1.3) with nonlocal competition. When the time delay is longer than the critical value $\tau_1^0$ and shorter than $\tau_0^0$, the prey and predator will coexist in space, but their densities will oscillate periodically in the form of spatially inhomogeneous oscillations for the model (1.3). However, the prey and predator will still coexist in a spatially homogeneous form, and their densities will converge to the positive equilibrium for the model (1.2). The model (1.3) is more realistic than the model (1.2) in the real world since it is difficult for the population density to exist in a completely uniform way. Through our theoretical analysis and numerical simulation, we show that stable inhomogeneous periodic solutions exist when the delay crosses the critical value, which is different from the work in ^[20].

In the inhomogeneous periodic oscillations of prey and predator densities (Figure 4), we find that the amplitude near the boundary of the region is larger than that at the center of the region. This is because the prey needs to avoid the hunting of predators. They will escape from the interior of the area and move out to the boundary. However, because the living area is closed, they cannot cross the boundary. Thus, they can only move inward again. After reaching the interior of the area, they have to escape the hunting of predators again. This forms the periodic oscillation of this mode.

In applying the predator-prey model, we suggest that scholars consider the nonlocal competition of prey since the dynamics of the models are different in the absence and presence of nonlocal competition. We hope that our model can be applied to ecological environment protection and population control. In future work, we will study the predator-prey model with a nonlocal competition term, nonconstant kernel function, and the Allee effect.

Acknowledgments

We would like to thank you for following the instructions above very closely in advance. It will definitely save us lot of time and expedite the process of your paper's publication. This research is supported by the Fundamental Research Funds for the Central Universities (Grant No. 2572022BC01), Postdoctoral Program of Heilongjiang Province (No. LBH-Q21060), and College Students Innovations Special Project funded by Northeast Forestry University (No. 202210225157).

Conflict of interest

The authors declare that they have no competing interests.