Spatiotemporal dynamics of a predator-prey model with a gestation delay and nonlocal competition

1.   Introduction

Dendrolimus is a forest pest with large occurrence and wide damage area. Among them, Dendrolimus superans is mainly distributed in Northeast China, Russia Far East, Japan, and so on. Dendrolimus superans are mainly parasitic on coniferous trees such as Larix gmelinii Kuzen and Pinus tabuliformis Carriere. Their larvae eat a large number of needles, resulting in trees being unable to carry out normal photosynthesis, eventually leading to a wide range of trees being destroyed or killed. Dendrolimus superans is extremely harmful to agriculture and forestry. In addition, the outbreak of Dendrolimus superans leads to tree death, reduces vegetation coverage, changes in plant population structure, and damages forest health and ecosystem stability. The life cycle of Dendrolimus superans consists of four stages: Egg, larva, pupa, and adult. In the life cycle of Dendrolimus superans, their development typically encompasses several key stages: Adult, egg, larva, and pupa, with each stage transitioning closely to the next. Specifically, adult insects will actively approach trees as a crucial behavior for seeking suitable environments to lay eggs. After the adults deposit their eggs on pine needles, the eggs will eventually hatch into larvae. During the larval stage, insects crawl on pine needles; this behavior is likely in search of food resources to meet the demands of their rapid growth and development. Subsequently, the larvae enter the pupal stage, during which they undergo a series of complex physiological and morphological changes. Finally, they metamorphose into adults through a process of transformation, thus completing a full life cycle. This process not only showcases the morphological and behavioral characteristics of insects at different life stages but also reflects their important strategies for adapting to the environment and reproducing. Among these stages, the larval stage is the most damaging to forests. The spread mechanism of the Dendrolimus superans is illustrated in Figure 1.

There are three major methods for the prevention and control of Dendrolimus superans: Physical control, chemical control, and biological control. Among them, physical control is mainly aimed at the behavior patterns and characteristics of Dendrolimus superans. However, this method requires a lot of manpower and material resources, which is not suitable for large-scale prevention and control ^[1]. Although chemical control is effective for large-scale pests, the use of pesticides and other chemical reagents may cause environmental pollution. By comparison, biological control is an environmentally friendly and sustainable pest control technology. The use of natural biological control of pests reduces the use of chemical pesticides, maintains ecological balance, avoids pest resistance, and is harmless to human health. Therefore, it is suitable for long-term pest control. Dendrolimus superans mostly has two major natural enemies: Parasitic natural enemies and predatory natural enemies. Among the parasitic natural enemies, Trichogramma and Exorista civilis play an important role; in terms of predatory natural enemies, it is represented by Common Cuckoo, Cyanopica cyanus, and so on, as illustrated in Figure 2.

When exploring the interaction between Dendrolimus superans and Cyanopica cyanus, we frequently employ the predator-prey model for analysis. This model is crucial for understanding and studying the dynamic changes in the predator-prey relationship between these two species within the realm of mathematical biology. In other studies, scholars have continuously refined and optimized these models ^{[2,3,4,5,6,7]}. The flight distance of adult female Dendrolimus superans is about 1.5 km, which shows that there are some limitations in their dispersal ability. Therefore, in the study of predator-prey model, the research environment can be regarded as relatively isolated. In such a relatively closed natural environment, scholars are concerned about the various factors that affect the system. Some biological studies have shown that habitat complexity affects population size and growth trends, indicating the important role of habitat complexity in the construction of ecological communities ^[8,9,10]. In addition, the construction of mixed forests provides habitat for natural enemies, increases environmental complexity, creates an ecological environment that is not conducive to the growth and development of Dendrolimus superans, and can effectively control the population of Dendrolimus superans ^[11,12]. Therefore, it is of great significance to study the biological control of Dendrolimus superans and their impact on forest ecosystems by using mathematical models. Based on the Lotka-Volterra model^[13], Ma and Wang^[14] proposed a predator-prey model with time delay and habitat complexity, that is

where $u$ and $v$ describe the population of prey and predator, respectively. $r$ is the intrinsic growth rate of prey. $K$ is the capacity of environmental for prey. $c$ is the attack rate of predator. $\beta$ stands for the strength of habitat complexity. $h$ is the handling time of predator. $e$ is the conversion efficiency. $\tau$ is the gestation delay of predator. $d$ is the death rate of predator. $\alpha$ is a positive real number. All parameters in model (1.1) are positive.

The relationship between Dendrolimus superans and its natural enemies is typically characterized using the predator-prey model. The functional response, a key component in predator-prey models, is frequently employed to describe the predation ability of predators. Grounded in experimental findings, Holling^[15] developed three functional responses, which are described below:

In a specific area, we assume that the predation rate will reach saturation when the population of natural enemies is sufficiently large. This means that the predation rate will not increase further even if the number of natural enemies continues to grow. Conversely, when the population of natural enemies is low and begins to increase, the predation rate rises more rapidly than a linear function. This nonlinear increase is due to the fact that predators become more efficient at capturing prey as their numbers increase. A natural selection for such a predation rate is the Holling type-II functional response, which captures this nonlinear relationship between predator population and predation rate. In other words, the value of parameter $\alpha = 1$ in model (1.1).

In a relatively closed natural environment, the spatial distribution of predators and prey significantly influences population dynamics, resource allocation, and predation efficiency. The nonlocal competition is in the interaction of predators and prey in different spatial locations, so that the system can predict the population dynamics and competition results more accurately. Peng and Zhang^[16] investigated a predator-prey model incorporating collective behavior and nonlocal competition among prey. The findings indicated that nonlocal competition had a destabilizing effect on the predator-prey system.

In addition to geographical distribution, ecological conditions, and other factors, time delay also plays an important role in the maintenance of ecological balance. In the predator-prey model, gestation delay is a key biological factor affecting the system, which is mainly reflected in the time difference between the occurrence of predation behavior and the birth of offspring. This physiological trait can significantly change the functional response pattern of predators to prey: During the breeding interval, the predator's predation frequency may decrease significantly or decrease to zero. This time delay mechanism directly affects the predator's growth rate by adjusting the number of effective predation per unit time and affects the dynamic balance of the ecosystem. The researchers in ^[17,18,19] developed a type of predator-prey model incorporating gestation delay and examined the dynamic characteristics of the systems.

Inspired by the above, we incorporate the nonlocal competition and gestation delay into the model (1.1). In the ecological environment, spatial distribution is often inhomogeneous, and spatial diffusion often occurs within populations. Therefore, when we study predator-prey models, reaction-diffusion equations may be more realistic. Thus, the resulting system is given by:

where $u$ and $v$ denote the density of Dendrolimus superans (prey) and Cyanopica cyanus (predator), respectively, the parameters $r$, $K$, $h$, $c$, and $d$ have the same meanings as in model (1.1). $d_1$ stands for the diffusion coefficients of Dendrolimus superans. $d_2$ stands for the diffusion coefficients of Cyanopica cyanus. $g$ denotes the efficiency with which energy is transferred from the species Dendrolimus superans to Cyanopica cyanus. $a$ is the attack rate of predator on prey. $\tau$ is the gestation delay of Cyanopica cyanus. $r$, $K$, $a$, $c$, $h$, $g$, $d$, $d_1$, $d_2$ are a positive constant. We select $\Omega \in (0, l\pi)$ where $l > 0$. $\widetilde{u} = \frac{1}{K}\int_{\Omega }{G(x, y)u(y, t)}dy$ represents the nonlocal competition, and the kernel function is given by $G(x, y) = \frac{1}{\left| \Omega \right|} = \frac{1}{l\pi }$, which is based on the assumption that the competition intensity among prey individuals in the habitat is uniform, meaning that the competition between any two preys is identical. The Neumann boundary condition is employed in this study, suggesting that the habitat is closed and effectively preventing any prey or predator from entering or leaving, thereby maintaining a self-contained ecosystem. We aim to investigate the dynamics of a predator-prey model, with a particular focus on the stability and dynamic properties of the system as the time delay parameter varies and serves as the bifurcation parameter.

The structure of this paper is as follows: In Section 2, we analyze the stability of positive constant steady states and the existence of Hopf bifurcation. In Section 3, we investigate the normal form of the Hopf bifurcation. In Section 4, we present numerical simulations. Finally, conclusions are drawn in Section 5.

2.   Stability analysis

System (1.2) has a boundary equilibria $E_0 = (0, 0)$ and a non-trivial equilibrium $E_1 = (u_*, v_*)$, where

If the following assumption holds:

then the system (1.2) must have a positive constant steady equilibrium $E_1 = (u_*, v_*)$.

By defining $U(x, t) = (u(x, t), v(x, t))^{\mathrm{T}}$, the linearized system for Eq (1.2) can be re-expressed as a differential equation at the equilibrium $E = (u_{0}, v_{0})\mathrm{, with}\; (u_{0}, v_{0}) = (0, 0)\; \mathrm{or}\; (u^{*}, v^{*}):$

where

with $a_{1} = r(1-\frac{{u_0}}{K})-\frac{ac{{v}_{0}}}{{{(1+ach{{u}_{0}})}^{2}}}, a_{2} = -\frac{ac{{u}_{0}}}{1+ach{{u}_{0}}}, b_{1} = \frac{gac{{v}_{0}}}{{{(1+ach{{u}_{0}})}^{2}}}, b_{2} = \frac{gac{{u}_{0}}}{1+ach{{u}_{0}}}, c_{1} = -\frac{ru_{0}}{K}.$

Especially, when $(u_0, v_0) = (u_*, v_*)$, $a_{1} = r(1-\frac{u_{*}}{K})\frac{hd}{g} > 0, a_{2} = -\frac{d}{g} < 0, b_{1} = r(g-hd)(1-\frac{u_{*}}{K}) > 0, \; b_{2} = d > 0, c_{1} = -\frac{ru_{*}}{K} < 0.$

Hence, the characteristic equation of (1.2) at $E = (u_{0}, v_{0})$ is given as follows:

where

with

2.1.   The case for $\tau=0$

When $\tau = 0$, Eq (2.1) for equilibrium $E_0 = (0, 0)$ becomes

For the case of $n = 0$ in Eq (2.2), the product of the two eigenvalues $-rd$ is negative. Therefore, the equilibrium $E_0 = (0, 0)$ is always unstable.

When $\tau = 0$, Eq (2.1) for equilibrium $E_1 = (u_*, v_*)$ becomes

Subsequently, we examine the conditions under which habitat complexity guarantees the stability of a positive constant steady state in the system (1.2).

$A_{n}+B_{n} = ({\frac{n}{l}})^{2}(d_{1}+d_{2})-a_{1}$. Due to $d_1+d_2 > 0$, there is $A_{n}+B_{n} > A_1+B_1 > A_0+B_0 = -(a_1+c_1)$. Assuming $A_0+B_0 = -(a_1+c_1) > 0$, there is $c < \frac{hd+g}{ahK(g-hd)}$. So we have $A_{1}+B_{1} = \frac{d_{1}+d_{2}}{l^{2}}-a_{1} > 0\Rightarrow \frac{d_{1}+d_{2}}{l^{2}}-\frac{rhd}{g}+\frac{rhd^{2}}{gk(g-hd)ac} > 0.$ Thus, if $rhdl^{2}-(d_{1}+d_{2})g > 0,$ that is $c < \frac{rhd^{2}l^{2}}{Ka(g-hd)[rhdl^{2}-(d_{1}+d_{2})g]}.$ Otherwise, if $rhdl^{2}-(d_{1}+d_{2})g < 0,$ that is $c > 0.$ Define:

Due to $a_{1} > 0, a_{2} < 0, b_{1} > 0, b_{2} > 0, \; C_{0}+D_{0} = -a_{2}b_{1} > 0$ holds. When $n = 0$, $-a_{2}b_{1} > 0$, $\frac{a_{1}}{2d_{1}}\geq0$ and $\triangle_1 = {{({{a}_{1}}{{d}_{2}})}^{2}}+4{{a}_{2}}{{b}_{1}}{{d}_{1}}{{d}_{2}} < 0$ $\Rightarrow c < \frac{dk}{(g-hd)a}{{(1-\frac{4{{d}_{1}}(g-hd)}{rh})}^{-1}}\triangleq{{c}^{*}}$, the function $C_{n}+D_{n} = (\frac{n}{l})^{4}d_{1}d_{2}-(\frac{n}{l})^{2}a_{1}d_{2}-a_{2}b_{1} > 0$ holds for any $n\in \mathbb{N}$.

In summary, if $c < \min\{\frac{hd+g}{ahK(g-hd)}, c_{0}, c_{*}\}$ holds, $A_{0}-B_{0} > 0, A_{n}-B_{n} > 0, C_{n}+D_{n} > 0$, for $n\in \mathbb{N}$ holds. Thus, we make the following hypothesis

When $(H_1)$ holds, the roots of characteristic equation (2.3) have negative real parts.

Theorem 1. For the model (1.2) with $\tau = 0$, the stability findings for equilibria are detailed below:

1) Equilibrium $E_0 = (0, 0)$ is always unstable.

2) Equilibrium $E_1 = (u_*, v_*)$ is always locally asymptotically stable when $(H_0)$ and $(H_1)$ hold.

3) Equilibrium $E_1 = (u_*, v_*)$ is always unstable when $(H_0)$ or $(H_1)$ does not hold.

2.2.   The case for $\tau\neq0$

Next, we investigate the existence of bifurcating periodic solutions near the positive constant steady state $E_1 = (u_*, v_*)$ for $\tau\neq0.$

When $\tau\neq0$, the characteristic equation describing the system (2.1) is presented as:

We may assume that $\lambda = \pm i\omega$ ($\omega > 0$) are a pair of purely imaginary roots of Eq (2.4). By substituting $\lambda = \pm i\omega$ into Eq (2.4) and separating the real and imaginary components, we derive:

that is,

with $n = 0, 1, 2\cdots.$ Then, for $n = 0, 1, 2, 3\cdots$ and $j = 0, 1, 2, 3\cdots$, we get

From (2.5), let $z = \omega^2$, we obtain:

Calculating the transversality condition, we get:

Denote

and

with $\Delta_n = (A_n^2 - 2D_n - B_n^2)^2-4(D_n^2 - C_n^2)$.

We define the following set

$G_1$ satis $(H_2)$ and $G_2$ satis $(H_3)$, respectively.

When $(H_2)$ holds, Eq (2.7) has the unique positive root ${{z}_{{{n}_{k}}}}$ for some $n = {{n}_{{{k}_{j}}}}\in G_1.$ Therefore, we can solve for ${{\omega}_{{{n}_{k}}}}$ and $\tau_{n_k}^{(j)}$, where

$\left.h^{\prime}(z_{n_k}) > 0, \text{so we have}\left. \mathrm{Re}(\frac{d\lambda}{d\tau_1})^{-1}\right|_{\tau = \tau_{n_k}^{(j)}} > 0.\right.$

When $(H_3)$ holds, Eq (2.7) has two positive root ${{z}_{{{n}_{t}, {i}}}}$ for some $n = {{n}_{{{k}_{t}}}}\in G_2.$ Therefore, we can solve for ${{\omega}_{{{n}_{t}, {i}}}}$ and $\tau_{n_t, i}^{(j)}, i = 1, 2$, where

Due to $\omega_{n_{t, 1}} < \omega_{n_{t, 2}},$ we have $h^{\prime}(z_{n_{t}, 1}) < 0$ and $h^{\prime}(z_{n_{t}, 2}) > 0.$ Thus $\mathrm{Re}(\frac{d\lambda}{d\tau})^{-1}|_{\tau = \tau_{n_{t}, 1}^{(j)}} < 0$, $\operatorname{Re}({\frac{d\lambda}{d\tau}})^{-1}|_{\tau = \tau_{n_{t}, 2}^{(j)}} > 0$ hold.

In summary, we can get the following conclusions.

Theorem 2. Under the framework of model (1.2) with $\tau \geq 0$, when $(H_0)$ and $(H_1)$ are valid and $\tau_{n}^{(j)}$ is defined by Eq (2.6), the subsequent conclusions hold.

1) Under the conditions $G_1 = \varnothing$ and $G_2 = \varnothing, \;$Eq (2.7) has no positive roots, and the positive constant steady state $E_1$ retains local asymptotic stability for all $\tau\geq0$.

2) Under the conditions $G_1\neq\varnothing$ and $G_2 = \varnothing, \;$Eq (2.7) admits a unique positive root $z_{n_k}$ for some $n_k\in G_1,$ the positive constant steady state $E_1$ of system (1.2) exhibits asymptotic stability for $0 < \tau_{1} < \tau_{c}^{*}$ and unstable for $\tau_{1} > \tau_{c}^{*}$, where $\tau_{c}^* = \min\{\tau_{n_k}^{(0)} \mid n_k\in G_1\}$. Additionally, system (1.2) undergoes $n_k\mathrm{-mode}$ Hopf bifurcation near $E_1$ when $\tau = \tau_{n_k}^{(j)}\; (n_k\in G_1, j = 0, 1, 2, 3\cdots).$

3) Under the conditions $G_1 = \varnothing$ and $G_2\neq\varnothing, \;$Eq (2.7) admits two positive roots $z_{n_t, i}\; (i = 1, 2)$ for some $n_t\in G_2,$ the positive constant steady state $E_1$ of system (1.2) exhibits asymptotic stability for $0 < \tau_{1} < \tau_{c}^{*}$, but a stability switch may occur for $\tau_{1} > \tau_{c}^{*}$, where $\tau_{c}^* = \min\{\tau_{n_t}^{(0)} \mid n_t\in G_1\}$. Furthermore, $n_t\mathrm{-mode}$ Hopf bifurcation emerges near $E_1$ at $\tau = \tau_{n_t}^{(j)}\; (n_t\in G_1, j = 0, 1, 2, 3\cdots).$

4) Under the conditions $G_1\neq\varnothing$ and $G_2\neq\varnothing, \;$Eq (2.7) admits positive roots $z_{n_k}$ and $z_{n_t, i}\; (i = 1, 2)$ and letting $n_p$ correspond to the smallest critical time delay $\tau_{c}^*$, where $n_p\in G_1$ or $n_p\in G_2$, the positive constant steady state $E_1$ of model (1.2) remains asymptotically stable for $0 < \tau_{1} < \tau_{c}^{*}$, but a stability switch may arise for $\tau_{1} > \tau_{c}^{*}$, where $\tau_{c}^* = \min\{\tau_{n_p}^{(0)} \mid n_p \in G_1\; or\; n_p \in G_2\}$. Furthermore, $n_p\mathrm{-mode}$ Hopf bifurcation occurs near $E_1$ at $\tau = \tau_{n_p}^{(j)}\; (n_p\in G_1\; or\; n_p\in G_2, j = 0, 1, 2, 3\cdots).$

3.   Normal form of Hopf bifurcation

In this section, we use the multiple time scales method to derive the normal form of the Hopf bifurcation for system (1.2).

Suppose that the characteristic equation (2.4) has a pair of pure imaginary roots $\lambda = \pm\; i\omega$ for $\tau = \tau_c$ Then, system (1.2) undergoes $n$-mode Hopf bifurcation near the equilibrium point $E_1$. Then, we derive the normal form of Hopf bifurcation of system (1.2) by the multiple time scales method. We choose time delay $\tau$ as a bifurcation parameter, where $\tau = \tau_c+ \varepsilon\mu$. The parameter $\varepsilon$ is a dimensionless scaling factor, $\mu$ is a perturbation parameter, and $\tau_c$ is given in Eq (2.8).

We let $u(x, t) \longrightarrow u(x, t)-u_*$, $v(x, t) \longrightarrow v(x, t)-v_*$, $\alpha = ac$, then model (1.2) can be rewritten as

where

Let $h = (h_{11}, h_{12})^T$ be the the eigenvector of the linear operator of system (3.1) corresponding to the eigenvalue $\lambda = i\omega$, and let $h^* = (h^*_{11}, h^*_{12})^T$ be the normalized eigenvector of the adjoint operator of the linear operator of system (3.1) corresponding to the eigenvalue $-i\omega$ satisfying the inner product $< h^{*}$, $h > = \overline{h^{*}}^{T}\cdot h = 1.$ By a simple calculation, we get

with $s = (\overline{h_{11}}h_{21}+h_{22}\overline{h_{12}})^{-1}$, and

Suppose the solution of Eq (3.1) is

where

The derivation with respect to $t$ is

where the differential operator $D_{i} = \frac{\partial}{\partial T_{i}}, i = 0, 1, 2, \cdots.$

We obtain

Denote

with $j = 1, 2, 3, \cdots.$

We take perturbations as $\tau = \tau_c+ \varepsilon\mu$ to deal with the delayed terms. By expand $U(x, t-\tau)$ at $U(x, T_0-\tau_c, T_1, T_2, ···),$ respectively, that is,

By substituting Eqs (3.3)–(3.5) into Eq (3.1), we obtain the expression for $\varepsilon$.

The solution of Eq (3.6) can be written as follows:

where $h_{11}$ and $h_{12}$ are given in Eq (3.2), and $c.c.$ means the complex conjugate of the preceding terms.

For the $\varepsilon^2$, we have

Substituting Eq (3.7) into the right side of Eq (3.8), we denote the coefficient vector of $e^{i\omega T_{0}}$ as $m_1$, from solvability condition $< h^{*}, (m_{1}, \cos\frac{n}{l}x) > = 0,$ we obtain

with

Assume the solution of Eq (3.8) to be as follows:

Denote

Substituting solutions Eqs (3.7) and (3.10) into Eq (3.8), we obtain

where

For the $\varepsilon^3$, we obtain

Substituting Eqs (3.7) and (3.10) into the right side of Eq (3.11), we denote the coefficient vector of $e^{i\omega T_{0}}$ as $m_2$, by solvability condition $< h^{*}, (m_{2}, \cos\frac{n}{l}x) > = 0,$ we obtain

where

with

According to the above analysis, the normal form of Hopf bifurcation for system (1.2) reduced on the center manifold is

where $M$ and $X$ are given by Eqs (3.9) and (3.13), respectively.

Let $G = re^{i\theta}$ and substitute it into Eq (3.14), and we obtain the Hopf bifurcation normal form in polar coordinates:

Therefore, we arrive at the following theorem.

Theorem 3. For system (3.15), if $\frac{\mathrm{Re}(M)\mu}{\mathrm{Re}(X)} < 0$ holds, model (1.2) exists periodic solutions near equilibrium $E_1 = (u_*, v_*)$.

1) If $\mathrm{Re}(M)\mu < 0$ the bifurcating periodic solutions reduced on the center manifold are unstable, and the direction of bifurcation is forward (backward) for $\mu > 0 \; (\mu < 0)$.

2) If $\mathrm{Re}(M)\mu > 0$ the bifurcating periodic solutions reduced on the center manifold are stable, and the direction of bifurcation is forward (backward) for $\mu > 0 \; (\mu < 0)$.

4.   Numerical simulations

In this section, we perform numerical simulations to verify the correctness of the theoretical analysis. In summary, we choose $d = 0.3, \; g = 0.6, \; a = 0.6, \; d_1 = 0.01, \; d_2 = 0.2, \; c = 0.55$, $r = 0.41$, $K = 5$, $h = 0.41$, and $l = 1$.

Based on the data above, $(H_0)$ and $(H_1)$ always hold, and system (1.2) has one unstable boundary equilibrium $E_0 = (0, 0)$, and one nontrivial equilibrium $E_1 = (u_*, v_*) = (1.90585, 0.96710)$. What matters most to us is the stability of the nontrivial equilibrium $E_1$.

According to Theorem 1, then $E_1$ is local asymptotically stable. This means that if system (1.2) is without delay, although Cyanopica cyanus and Dendrolimus superans can coexist at this time, natural predators are capable of suppressing the reproduction of the pests.

When $\tau\neq0$, by a simple calculation, $(H_2)$ only holds for $n = 1$ there does not exist any $n$ such that $(H_3)$ holds. Thus, Eq (2.9) has a unique positive root $\omega = 0.1616$, which corresponds to the critical delay $\tau_{1}^{(0)} = 2.5897.$ From the definition of $\tau_{c}^*,$ we obtain $\tau_{c}^* = 2.5897$, theorem 2.2 supports the conclusion that the positive constant steady state $E_1$ achieves local asymptotic stability when $\tau\in [0, 2.5897)$, and instability when $\tau\in (2.5897, +\infty)$. We select $\tau = 1.5\in [0, 2.5897)$ to conduct the simulation. See Figure 3.

We choose $\tau = 2.8 > \tau^*_{c} = 2.5879,$ from Eqs (3.9) and (3.12), we obtain Re$(M) > 0$, Re$(X) < 0$. Thus, according to Theorem 3, system (1.2) will generate inhomogeneous periodic solutions near the positive constant steady state $E_1$ of model (1.2), and bifurcating periodic solutions are stable and forward (see Figure 4).

By analyzing Figures 3 and 4, we can draw the following conclusions.

1) When the gestation delay of the predator is less than the critical value, the predator can quickly control the number of pests and reach a stable state by preying on pests.

2) When the gestation delay of the predator is slightly more than the critical value, the predator has a certain control effect on the number of pests. This will cause the number of pest populations to erupt periodically, and the number of natural enemy populations will fluctuate periodically.

3) Comparing Figures 3 and 4, when the gestation delay of insectivorous birds is long, so it is necessary to artificially intervene in the population of pests to prevent the outbreak of insect pests. By releasing insectivorous birds to increase their number, we can indirectly reduce the increase in the number of pests and the outbreak of pests due to the prolonged pregnancy of predators, so that the population of pests can return to a controllable stable state.

5.   Conclusions

In this paper, aiming at the population density control problem of the pine caterpillar, considering the Holling-II type functional response, we developed a pest control model with gestation delay and nonlocal competition. For our analysis, we focused on two major aspects: first, the existence and stability of the positive constant steady state, and second, the conditions under which Hopf bifurcations emerge near this steady state. In the part of numerical simulation, we selected a set of suitable parameters for numerical simulation. It provided an explanation for effectively controlling the population density of the pine caterpillar and maintained the stability between natural enemies and the pine caterpillar species, so as to realize effective environmental protection and true green pest control. Therefore, gestation delay has a significant impact on the stability of the population. When the gestation delay is less than the critical delay, the predator can effectively control the pest population. When the gestation delay exceeds the critical delay, the pest population experiences periodic outbreaks.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was funded by the Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2024A001), and the College Students Innovations Special Project funded by Northeast Forestry University (No. DCLXY-2025011).

Conflict of interest

The authors declare there are no conflicts of interest.