Hopf bifurcation analysis of a pine wilt disease model with both time delay and an alternative food source

1.   Introduction

Pine wilt disease is a forest disease that causes devastating harm to the ecological environment, especially pine forest resources. It is transmitted by the Bursaphelenchus xylophilus, which is spread by Monochamus alternatus. This disease is characterized by rapid onset, a long incubation period, and difficulties in management. Despite the considerable efforts that have been made, the scope of the epidemic continues to expand. Since it was first discovered in China in 1982, it has been listed as a key forest quarantine object in China^[1]. This type of pest and disease has resulted in significant losses to China's forestry resources ^[2]. In 2017, the total area of the disease exceeded 8 hectares and showed a continuous growth trend ^[3]. The pathogenesis of pine wilt disease is shown in Figure 1: when a Monochamus alternatus carrying Bursaphelenchus xylophilus feeds on or lays eggs in a healthy tree, the Bursaphelenchus xylophilus detaches from the Monochamus alternatus and enters the bite marks made by the Monochamus alternatus on the healthy tree, thus infecting the healthy tree with pine wilt disease. Trees infected with Bursaphelenchus xylophilus first turn fulvous, and then die after a period of time. Susceptible Monochamus alternatus beetles, when passing through these infected trees, become hosts for Bursaphelenchus xylophilus, which then parasitize them^[4]. As the cycle progresses, the spread of pine wilt disease accelerates. (Figure 1 is obtained by image reconstruction from https://www.vcg.com/).

There are three main types of prevention and control of pine wilt disease: physical, chemical, and biological. According to the study, the main transmission carrier of Bursaphelenchus xylophilus is Monochamus alternatus, so the prevention and control of pine wilt disease is mainly to control Monochamus alternatus. Physical control is primarily through the habits of the Bursaphelenchus xylophilus and is suitable for use on a small infestation scale. This control requires a lot of human and material resources and is too inefficient for controlling Bursaphelenchus xylophilus. Chemical control is mainly by spraying the canopy at ground level, but since the period of susceptibility of the Monochamus alternatus coincides with the rainy season, spraying has less impact on the Bursaphelenchus xylophilus. Although many improvements have been made to the chemical method, this method still has a certain impact on the ecological environment of the forest system and causes environmental pollution. Considering the protection of the environment and the healthy development of the ecology, biological control has become a hot topic of research. The natural enemies of Monochamus alternatus can be divided into two categories: parasitic and predatory natural enemies. Parasitic natural enemies mainly include Scleroderma guani Xiao et Wu ^[5] and Dastarcus helophoroides (Fairmaire); predatory natural enemies mainly include birds such as Dryocopus martius and Pica pica, as shown in Figure 2. (Figure 2 is obtained by image reconstruction from https://www.vcg.com/).

In recent years, some scholars have used mathematical models to study the spread of pine wilt disease ^[4,6]. In ^[6], Khan et al. constructed and analyzed the transmission dynamics model of pine wilt disease based on Caputo fractional derivatives. The results show that the fractional derivative can provide more efficient and flexible characterization for model dynamics, which has important theoretical value for optimizing the biological control strategy of pine wilt disease. In ^[4], the authors proposes a reaction-diffusion predator-prey model with nonlocal effect and memory diffusion for the prevention and control of pine wilt disease transmitted by Monochamus alternatus. The nonlocal effect factor and memory diffusion factor are discussed in detail. However, the additional food source factor is not discussed. Based on the inspiration of ^[4], in this paper, the effects of the pregnancy delay factor and additional food source factor on the dynamics of the system are analyzed. Therefore, it is interesting to study the dynamic relationship between Monochamus alternatus and natural enemies by a mathematical model in the prevention and control of pine wilt disease.

With regard to the predator-prey model, some scholars have studied the dynamic properties of the system through ordinary differential equations or partial differential equations ^[4,7,8,9]. When describing the dynamic relationship between Monochamus alternatus and natural enemies, there are also scholars who have adopted the predator-prey model, which is interesting for studying the population dynamics of two populations with a predator-prey relationship ^[10,11]. In^[10], Li and Ding established a predator-prey model with time delay for the control of Bursaphelenchus xylophilus based on the mutually beneficial symbiosis and parasitism among Bursaphelenchus xylophilus, Monochamus alternatus, and Dastarcus helophoroides. In ^[11], in order to control pine wilt disease, Hou et al. established a predator-prey model with prey-taxis and nonlocal intraspecific competition of prey, and studied the spatial form of the pine wilt disease transmission system. In this regard, this paper also wants to continue to study the dynamic influence of the time delay factor on the pine wilt disease model. The functional response function is the core and foundation of the predator-prey model. It represents the predator's predation ability. The functional response function is affected by many factors, such as the structural features of the habitat, the predatory capability of the predator, the evasive capacity of the prey, and other relevant factors^[12]. In ^[13], Holling proposed three types of functional response functions, which are as follows:

where $u$ denotes the density of the Monochamus alternatus, $m$ denotes the predation rate of the natural enemies, and $a$ denotes the semi-saturation rate of the natural enemies.

Many mathematical ecologists add the Holling II functional response function to the predator-prey model to reflect the dynamic complexity of the species interaction system, and show how to study the stability and Hopf bifurcation behavior in the interacting species system based on the Holling II functional response function ^{[14,15,16,17,18]}. In ^[19], Wang and Yu employed the Holling II functional response function to study the stability and Hopf bifurcation behavior of the Bazykin predator-prey ecosystem, and carried out theoretical and numerical studies. In ^[11], Hou et al. established a reaction-diffusion equation system including the Holling II response function, prey-taxis, and nonlocal intraspecific competition of prey to study the spatial pattern formation mechanism of pine wilt disease transmission. These studies have shown that the Holling II functional response function can enhance the dynamic effect of the predator-prey model, so in this paper, Holling II is used as the functional response function. Then we give the following model:

where $u$ and $v$ are the densities of Monochamus alternatus and natural enemies, $r$ represents the growth rate of Monochamus alternatus $u$ under the premise of avoiding natural enemies $v$, $k$ represents the maximum carrying capacity of the environment for Monochamus alternatus, $c_{0}$ represents the conversion rate of Monochamus alternatus eaten by natural enemies, and $D$ represents the mortality rate of the natural enemies $v$, where $m$, $k$, $a$, and $c_{0}$ are positive values.

The organizational structure of this paper is as follows. In Section 2, we establish a model with the Holling II functional response function, Beverton-Holt-like alternative food source, and pregnancy delay, and the model is dimensionless. In Section 3, we consider the existence of a constant steady state solution in the model. In Section 4, we study the stability of the constant steady state solution and the conditions of Hopf bifurcation in the model. In Section 5, we study the properties of Hopf bifurcation. In Section 6, we give some numerical simulations to confirm the correctness of our theoretical research. In Section 7, we finally provide a brief summary to conclude this paper.

2.   Model formulation

From the perspective of biology and practical applications, time delay factors play a crucial role in the predator-prey model. Incorporating the time delay effect into the predator-prey model can make the model closer to the ecological reality, and then reveal more abundant and complex dynamic behavior characteristics. Many scholars confirmed through rigorous research that time delay factors had a profound impact on the stability of population density^[20,21]. In^[22], the authors proposed a time-delayed susceptible-asymptomatic-infected-removed (SAIR) model considering the temporary immune characteristics, analyzed the impact of the COVID-19 vaccine on the epidemic dynamics, and discussed the impact of delay on system stability and Hopf bifurcation. Yang and Ding^[8] studied the effects of pests on plants by constructing a model containing a delay differential equation under the influence of temperature. Wang and Yang^[23] employed the gestation time delay of predators as the bifurcation parameter. They investigated the existence of Hopf bifurcation and ascertained the direction of the Hopf bifurcation and the stability of the resulting periodic solutions by analyzing the distribution of eigenvalues. These studies have shown that time delay can enrich the dynamic effect of the predator-prey model, and the time delay factor will induce the Hopf bifurcation and then produce the periodic solution phenomenon. In this paper, the periodic change of the density of Monochamus alternatus reflects the periodic change of the outbreak of pine wilt disease. The periodic outbreak of pine wilt disease has an impact on the ecological environment, so it is necessary to study the time delay factor.

Furthermore, many researchers have found that providing alternative food for predators can lead to a weaker aggression of predators and an increase in the growth rate for predators^[24]. In recent years, many well-known scholars have also added the factor of providing alternative food for predators to the predator-prey model, and have achieved good results^{[25,26,27,28]}. In ^[29], van Baalen et al. assumed that the alternative food had a fixed density, and analyzed the influence of the strategy of providing alternative food for predators on the dynamics of the predator-prey system. Therefore, in order to better reflect the complexity and dynamics of the predator-prey model, the present model takes into account this factor of providing alternative food for predators. Further, by the paper^[24], we use the Beverton-Holt-like alternative food source model, where the specific mathematical expression is

where $\varepsilon$ is the maximum per capita reproduction rate $v$, and $\lambda$ refers to the intensity coefficient that the density of natural enemies $v$ depends on. Therefore, the model in this paper is updated to a predator-prey model considering a Beverton-Holt-like alternative food source.

Similar to chemistry, the direct interrelationships of populations in biological population dynamics are often reflected by reaction-diffusion equations^[30,31,32]. At present, many scholars use the reaction-diffusion model to study the predator-prey model ^{[33,34,35,36]}. Therefore, the reaction-diffusion equation for pine wilt disease based on the Holling II functional response function, Beverton-Holt-like alternative food source, and pregnancy delay studied in this paper is

where $u(x, t)$ and $v(x, t)$ are the densities of Monochamus alternatus and natural enemies. $D_{1}$ and $D_{2}$ represent the diffusion coefficients of the Monochamus alternatus $u$ and the natural enemies $v$, respectively, and $\tau$ represents the pregnancy delay of the natural enemies. $D_{1}$, $D_{2}$, $r$, $k$, $m$, $a$, $c_{0}$, $D$, $\varepsilon$, and $\lambda$ are all positive numbers. $\Omega = (0, l\pi)$, $\nu$ is the outer normal vector of the boundary $\partial \Omega$, and $\Delta$ is denoted as a Laplace operator. $\frac{\varepsilon v}{1+\lambda v}$ refers to the influence of generalized predators on the model by preying on animals other than the Monochamus alternatus $u$ considered in this paper. The functional response function used in this paper is Holling type Ⅱ, and the term considering providing alternative food for predators is Beverton-Holt-like. The structure of the system is shown in Figure 3. As far as we know, the reaction-diffusion equation based on the Holling II functional response function, Beverton-Holt-like alternative food source, and delayed pregnancy has not been used to study pine wilt disease. The main purpose of this paper is to study the Hopf bifurcation phenomenon induced by diffusion-driven, delay-induced, and alternative food factors.

In order to simplify the calculation, we make (2.2) dimensionless. Denote $\tilde{u} = \frac{u}{k}$, $\tilde{v} = mv/r$, $\tilde{t} = tr$, and then system (2.2) is changed to (in order to better represent the equation, we have removed the tilde):

where $d_1 = \frac{D_{1}}r$, $d_2 = \frac{D_{2}}r$, $d = \frac{D}r$, ${s} = ak$, $c = \frac{c_{0}mk}r$, $e = \frac{r}\varepsilon$, $f = \frac{\lambda r^2}{\varepsilon m}$. We assume $\Omega = (0, l\pi)$, where $l > 0$.

3.   Equilibrium analysis

In this section, we prove the existence of positive solutions in the system. Solving the following equations,

we obtain that (0, 0) and (1, 0) are obviously equilibrium points, and there exists an equilibrium point ($u_*$, $v_*$), where $v_* = (1-u_*)\left(1+su_*\right)$, and $u_*$ is the root of ${H}(u) = 0$,

where

Then we make the following assumptions:

${\bf (H_0)}\; \; \; \; \; \; \; \; \; -d\left(e+f\right)+1 < 0, \; \; \mathrm{\; }and-e\left(d-c+ds\right)+1+s > 0.$

Theorem 1. If the parameters satisfy assumptions ${\bf (H_0)}$, $u_*\in\left(0, 1\right)$ is the root of ${H}(u)$. Then the system (2.3) has a positive equilibrium solution ($u_*$, $v_*$), where $v_* = (1-u_*)\left(1+su_*\right)$.

Proof. Because ${H}(u)$ is a continuous function, and ${H}(0) = -d\left(e+f\right)+1 < 0$, ${H}(1) = -e\left(d-c+ds\right)+1+s > 0$. Then by the intermediate value theorem, there is at least one point $u_*\in\left(0, 1\right)$ such that ${H}(u) = 0$. Since $u_*\in\left(0, 1\right)$, then $v_* = (1-u_*)\left(1+su_*\right) > 0$. □

4.   Stability analysis

We use a method similar to reference ^[37] to analyze the stability of the system. The linearized system (2.3) at $(u_*, v_*)$ is

where $U(x, t) = (u(x, t), v(x, t))^T$,

and

Denote $\mathbb{N}_{1} \triangleq\{0, 1, 2, 3\cdots\cdots\}$. The characteristic equation is

where

4.1.   The case of $\tau=0$

When $\tau = 0$, the characteristic equation (4.2) is in the following form:

where

and since (4.3) is a quadratic equation with one variable, then the eigenvalues satisfy

Then, we make the following hypotheses:

Theorem 2. Suppose $d_1$ = $d_2$ = 0, $\tau$ = 0, and they satisfy hypothesis (4.6). Then the positive equilibrium solution $(u_*, v_*)$ is locally asymptotically stable.

Proof. Hypotheses (4.4) and (4.6) imply that $tr_n < 0$ and $\Delta_n > 0$. Then the two roots of (4.3), i.e., the eigenvalues, have negative real parts, so the positive equilibrium solution $(u_*, v_*)$ is locally asymptotically stable. □

Divide the parameters into the following three cases:

Denote

where $\Delta_k$ refers to the value of $n$ in (4.4) belonging to $\mathbb{S}_1$ in $\Delta_n$.

Theorem 3. Suppose (4.6) holds and $\tau = 0$.

1) In Case 1 (or Case 2), the positive equilibrium solution $(u_*, v_*)$ in the differential system (2.3) is locally asymptotically stable;

2) If $\mathbb{S}_1 = \varnothing$, in Case 3 we have that the positive equilibrium solution $(u_*, v_*)$ in the differential system (2.3) is locally asymptotically stable.

Proof. From the hypothesis (4.6), we know that $tr_0 < 0$ and $\Delta_0 > 0$, and then we can easily see that for $n \in \mathbb{N}_0$, there are $tr_0 < 0$. When the parameters satisfy Case 1 (or Case 2), there are $\Delta_n > 0$ for $n \in \mathbb{N}_0$, and then it can be seen that the eigenvalues of the equation (4.3) have negative real parts, which shows that statement $(1)$ holds. Since when the parameters satisfy Case 3, for $n \in \mathbb{N}_0$, there are $\Delta_n > 0$, then statement $(2)$ holds. □

4.2.   The case of $\tau\neq 0$

We now assume that $\tau > 0$, and suppose that (4.6) and one of the conditions (1) or (2) in Theorem 3 hold. Then we assume that $i w$ ($\omega > 0$) is the solution of Eq (4.2). We obtain

Then we have

This leads to

Denote $z = w^2$, and Eq (4.9) is

Since we satisfy condition (1) or (2) of Theorem 3, then we have

By calculating, we arrive at

Define

Lemma 4.1. Suppose the condition (1) or (2) of Theorem 3 and hypothesis (4.6) hold.

1) For $n\in \mathbb{M}_1$, Eq (4.2) has a pair of purely imaginary roots $\pm i\omega_{n}$ at $\tau^{j}_n, \; j \in \mathbb{N}_1.$

2) For $n\in \mathbb{M}_2\cup\mathbb{M}_4$, Eq (4.2) has no purely imaginary root.

3) For $n\in \mathbb{M}_3$, Eq (4.2) has two pairs of purely imaginary roots $\pm i\omega_{n}^{\pm} \mathit{\text{at}} \tau_{n}^{j, \pm}, j \in \mathbb{N}_{1}$.

Proof. The roots of (4.10) are

Then

where $j\in\mathbb{N}_{1}.$ When $n\in \mathbb{M}_1$, since $Q_n < 0$, then obviously $z^{+}$ is a positive real root and $z^{-}$ is a negative real root. Then conclusion (1) holds. When $n\in \mathbb{M}_2$, Eq (4.10) has no roots; and when $n\in \mathbb{M}_4$, $z^{+}$ and $z^{-}$ are negative roots. Then, because of (4.11), conclusion (2) holds. When $n\in \mathbb{M}_3$, $z^{+}$ and $z^{-}$ are positive real roots, and by (4.11), conclusion (3) holds. □

Lemma 4.2. Assume (4.6) holds, and the parameters satisfy condition (1) or (2) of Theorem 3. Then when $\tau = \tau_{n}^{j, +}$, $Re\left(\frac{d\tau}{d\lambda}\right) > 0$, and when $\tau = \tau_{n}^{j, -}$, $Re\left(\frac{d\tau}{d\lambda}\right) < 0$ for $n \in \mathbb{M}_1\cup\mathbb{M}_3$ and $j \in \mathbb{N}_1$.

Proof. We adopt a method similar to that in references^[11,12,37], and the derivation of parameter $\lambda$ in Eq (4.2), and we obtain

Then

Therefore $[\text{Re}(\frac{d \lambda} {d \tau})^{-1}]_{\tau = \tau^{j, +}_n} > 0$, $[\text{Re}(\frac{d \lambda} {d \tau})^{-1}]_{\tau = \tau^{j, -}_n} < 0$, so $\text{Re} (\frac{d \tau} {d \lambda})_{\tau = \tau^{j, +}_n} > 0$ and $\text{Re} (\frac{d \tau} {d \lambda})_{\tau = \tau^{j, -}_n} < 0$. □

It is obvious from Eqs (4.11) and (4.12) that there must be $\tau^{0, \pm}_n < \tau^{j, \pm}_n$ $(j\in \mathbb{N}_1)$. For $n\in \mathbb{M}_1\cup\mathbb{M}_3$, define $\tau_{c} = \mbox{min}\{\tau^{0, \pm}_n or \; \tau^{0, +}_n \mid n \in \mathbb{M}_1\cup\mathbb{M}_3\}$. We adopt a standardized method similar to that in references ^[12,38], and it can be seen from the above that the following theorem holds.

Theorem 4. For system (2.3), suppose condition (1) or (2) of Theorem 3 and assumption (4.6) hold, and then we have the following conclusions.

1) If $\mathbb{M}_1\cup\mathbb{M}_3 = \varnothing$, then $(u_*, v_*)$ is locally asymptotically stable when $\tau \geq0$.

2) If $\mathbb{M}_1\cup\mathbb{M}_3\neq \varnothing$, then $(u_*, v_*)$ is locally asymptotically stable when $\tau \in [0, \tau_{c})$ and unstable for $\tau > \tau_{c}$.

3) A Hopf bifurcation occurs when $\tau = \tau^{j, +}_n (or \; \tau = \tau^{j, \pm}_n)$, $j\in \mathbb{N}_1$, $n \in \mathbb{M}_1\cup\mathbb{M}_3$.

5.   Properties of a Hopf bifurcation

In this section, we will use the multiple time scales method to derive the normal form of the Hopf bifurcation of model (2.3) with gestational time delay according to references ^[39,40].

When the critical value $\tilde{\tau} = \tau_{c}$, the characteristic equation (4.2) has a pair of pure imaginary roots $\lambda = \pm i \omega$. Then system (2.3) will undergo a Hopf bifurcation. We consider $\tau$ as a bifurcation parameter, $\tau = \tilde{\tau} + \varepsilon\mu$, where $\tilde{\tau}$ is the Hopf critical point, $\mu$ is the perturbation parameter, and $\varepsilon$ is the scale parameter. From the multiple time scales method, we can derive the normal form of the Hopf bifurcation for system (2.3). We use the Taylor expansion method to expand system (2.3) to the third order at $E = (u^{*}, v^{*})$, and take $\tilde{u}(x, t) = u(x, \tau t)-u_*$ and $\tilde{v}(x, t) = v(x, \tau t)-v_*$ into system (2.3). For convenience, we still use $u\left(x, t\right)$ and $v\left(x, t\right)$ instead of $\tilde{u}\left(x, t\right)$ and $\tilde{v}\left(x, t\right)$, so we get:

where $f_1 = \frac{v^*}{\left(1+su^*\right)^2}$, $f_2 = \frac{u^*}{1+su^*}$, $f_{3} = \frac{sv^{*}}{\left(1+su^{*}\right)^{3}}$, $f_4 = \frac{1}{\left(1+su^*\right)^2}$, $f_{5} = \frac{s^{2}v^*}{\left(1+su^*\right)^{4}}$, $f_6 = \frac{s}{\left(1+su^*\right)^3}$, $\beta_{1} = \frac{cv^{*}}{(1+su^{*})^{2}}$, $\beta_{2} = \frac{cu^{*}}{1+su^{*}}$, $\beta_{3} = \frac{csv^{*}}{\left(1+su^{*}\right)^{3}}$, $\beta_4 = \frac{c}{(1+su^*)^2}$, $\beta_5 = \frac{cs^2v^*}{(1+su^*)^4}$, $\beta_{6} = \frac{cs}{\left(1+su^*\right)^{3}}$, $\gamma_{1} = \frac{e}{(e+fv^*)^{2}}$, $\gamma_{2} = \frac{ef}{\left(e+fv^*\right)^{3}}$, $\gamma_{3} = \frac{ef^{2}}{\left(e+fv^{*}\right)^{4}}$.

Let $h = (h_{11}, h_{12})^{T}$ be the eigenvector of the eigenvalue $\lambda = \mathrm{i}\omega$ corresponding to the linear matrix of the linearized system of (5.1), $h^{*}$ is the eigenvector of the eigenvalue $\lambda = \mathrm{-i}\omega$ corresponding to the adjoint matrix of the linear matrix of the linearized system of (5.1), and we have $< h^*, h > = \overline{h^*}^T\cdot h = 1$. By calculation, we obtain

Denote

The solution of the equation is

where

The division of the derivation of $t$ is

where $D_i$ represents the differential operator, $D_i = \frac{\partial}{\partial T_i}, i = 0, 1, 2, \cdots$.

Denote

Available from (5.4),

In order to deal with the delay term, we take the perturbation as $\tau = \tilde{\tau} + \varepsilon\mu$. We expand $u(x, t-\tau), v(x, t-\tau)$ in $u(x, T_{0}-\tilde{\tau}, T_{1}, T_{2}, \cdots)$ and $v(x, T_{0}-\tilde{\tau}, T_{1}, T_{2}, \cdots)$. What we can get is

where $u_{j, 1} = u_{j}(x, T_{0}-1, T_{1}, T_{2}, \cdots), \ v_{j, 2} = v_{j}(x, T_{0}-1, T_{1}, T_{2}, \cdots), \ j = 1, 2, 3, \cdots.$

Bring (5.4)–(5.9) into (5.1), the following equation can be obtained by comparing the coefficients in front of $\varepsilon$.

Then the solution of (5.10) is

where $h_{11}$ and $h_{12}$ are given in Eq (5.2), and $G = G(T_1, T_2, \dots)$.

For the $\varepsilon^{2}$, we have

Substituting Eq (5.11) into Eq (5.12), by the solvability condition, we obtain

with

with

Suppose the solution of Eq (5.12) is as follows:

Denote

Substituting the solutions to Eq (5.11) and Eq (5.16) into the right side of Eq (5.12), we obtain

where

For the $\varepsilon^{3}$ term, we have

Substituting the solutions to Eq (5.11) and Eq (5.16) into the right side of Eq (5.18), we obtain the coefficient vector of term $e^{i{\omega}T_{0}}$, denoted as $m_{2}$, and let $< h^{*}, m_{2} > = 0$. We obtain

where

with

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

making $G\rightarrow G/\varepsilon$, and thus, Eq (5.22) becomes:

where $M_{1}$ and $\chi$ are given by Eq (5.13) and Eq (5.19), respectively.

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

Theorem 5. For system (5.24), if $\frac{\mathit{\text{Re}}(M_{1})\mu}{\mathit{\text{Re}}(\chi)} < 0$ holds, system (2.3) has periodic solutions near equilibrium $E = (u_{*}, v_{*})$.

1) If $\mathit{\text{Re}}(M_{1})\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 $\mathit{\text{Re}}(M_{1})\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)$.

6.   Numerical simulations

In this section, we conduct analysis based on relevant literature and make self-assumptions, determine the appropriate parameter values for the model, carry out numerical simulations to verify the correctness of the theoretical analysis, and provide theoretical support for the prevention and control of forest pests and diseases.

6.1.   Selection of parameter values

In this subsection, we perform data analysis and make self-assumptions to determine parameter values for simulations.

1) The coefficient $m$ and the semi-saturation rate $a$.

According to the results from reference ^[4], the predator's predation rate $m$ is 0.63. In addition, we may choose the semi-saturation rate of the predator $a$ to be 0.4.

2) Diffusion coefficients $D_{1}$, $D_{2}$.

Compared with the predatory natural enemy Picoides, the movement speed of Monochamus alternatus is relatively slow. Therefore, we choose the diffusion coefficient $D_{1}$ of Monochamus alternatus to be 0.06, and the diffusion coefficient $D_{2}$ of predatory natural enemies to be 0.1.

3) Birth rate of prey $r$.

From the analysis of reference ^[41], we can get that the monthly egg production of Monochamus alternatus is 78.6, and the daily egg production of Monochamus alternatus is 2.62. The formula for the birth rate of Monochamus alternatus defined therein can be obtained from reference ^[41]. By a simple calculation, the birth rate of Monochamus alternatus is about 1. Therefore, we choose $r$ = 1.

4) Maximum average growth rate $\varepsilon$ and strength coefficient density $\lambda$.

For the maximum average growth rate of predators, we may choose $\varepsilon$ = 0.5; for the strength coefficient density of predators, we may choose $\lambda$ = 0.7.

5) Other parameters $D$, $k$, $c$.

For the environmental capacity of Monochamus alternatus, we may choose $k$ = 50, $c$ = 0.5, $D$ = 0.78 for simulations.

Due to the dimensionless characteristic of this paper, according to (2.3), the summary analysis is

6.2.   Simulations

In this section, we perform numerical simulations for model (2.3). The simulation results can provide reference and a theoretical basis for preventing and controlling the outbreak of pine wilt disease.

When the above parameters are satisfied, we have $-d(e+f)+1 = -2.2955 < 0$ and $-e(d-c+ds)+1+s = 19.74 > 0$, so we have satisfied ${\bf (H_0)}$. Then we have $E = (u_{*}, v_{*}) = (0.4569, \; 5.5059)$, $a_1+a_3+b_2 = -0.0277706 < 0$, and $a_1 (b_2+a_3)-a_2 b_1 = 0.0360561732433 > 0$. We have $(b_{2}+a_{3})d_{1}+a_{1}d_{2} = -0.000364161 < 0.$ So we have satisfied (4.6) and $\bf{Case\; 1}$. When $\tau$ = 0, ${\bf (H_0)}$, (4.6), and $\bf{Case\; 1}$ hold, according to Theorem 3, so we know that the constant steady state solution $E = (u_*, v_*) = (0.4569, \; 5.5059)$ is locally asymptotically stable, see Figure 4.

Biological interpretation 1: This means that if model (2.3) is without a time delay term, although natural enemies and Monochamus alternatus can coexist at this time, natural enemies can suppress the reproduction of Monochamus alternatus. Then we can know that pine wilt disease can be effectively treated.

When $\tau\neq 0$, by simple calculation, when $n\in \mathbb{K}_1, \mathbb{K}_1\triangleq \{0, 1, 2\}$ and $n\in \mathbb{M}_1$, the above parameters satisfy the conditions of Lemma 4.1, so there is a pair of pure imaginary roots $\omega_{0}$. Therefore, Eq (4.2) has a unique positive root $\omega_{0} = 0.1630394$, where $\omega_{0} = 0.1630394$ corresponds to the critical delay $\tau_{0}^{0} \approx 0.492222$. According to the definition of $\tau_{c}$, we choose $\tau_{c} = 0.492222$. According to Theorem 4, we determine that when $\tau_{1} = 0.35\in[0, 0.492222)$, the constant steady state solution $E = (u_*, v_*)$ is locally asymptotically stable, see Figure 5.

We choose $\tau_{2} = 0.55 > \tau_{c} = 0.492222$, and thus, according to Theorem 4, system (2.3) will generate homogeneous periodic solutions near the positive constant steady state solution $E$ of system (2.3), see Figure 6.

Biological interpretation 2:

1) Without the interference of other factors, when the predator's gestation period is lower than the critical value $\tau_{c} = 0.492222$, the large increase of predators leads to the rapid control of the population of Monochamus alternatus, which has made it reach a stable state.

2) When the gestation period is slightly higher than the critical value $\tau_{c}$, the natural enemies will exert a certain degree of control on the population of Monochamus alternatus, but they cannot effectively inhibit the reproduction of these pests. Therefore, the formation of pests has the characteristics of periodic outbreaks. In the periodic control of pests, the number of natural enemies also shows a cyclical growth. Therefore, when the natural enemies of Monochamus alternatus have the characteristics of long gestation, we need to take some manual intervention measures. We can increase the number of natural enemies by artificially providing them with an additional food source, indirectly reducing the negative impacts of their long gestation period and returning pest levels to a manageable and stable state.

Then, we are going to modify the magnitude of the time delay parameter within the stable interval $[0, 0.492222)$ of system (2.3). When $\tau_{3} = 0.25$, the time required for the solution to remain stable is less than the time required for $\tau_{1} = 0.35$, as shown in Figure 7.

Biological interpretation 3: If the time of the predator's gestation is kept within the critical time $\tau_{c} = 0.492222$ (years), the spread of pine wilt disease will be controlled and there will be no large-scale outbreaks. In particular, the time required to maintain local asymptotic stability is about 960 years when $\tau = 0.25$, and 1750 years when $\tau = 0.35$. Although the time required to maintain stability is not in line with the actual situation, it is obvious that the shorter the pregnancy time of the predator, the faster pine wilt disease can be controlled, and the less damage it will cause to forestry resources.

It is worth noting that although it is difficult to alter the gestation time of natural enemies of Monochamus alternatus, we can adopt additional alternative foods for natural enemies of Monochamus alternatus for the purpose of shortening the gestation period of natural enemies. Although we make the predator obtain additional alternative food, in the two populations with a predator-prey relationship, when the number of predator groups suddenly increases sharply, for the long term, it will lead to a decrease in the predation rate. So let $\varepsilon$ = 0.6, $m$ = 0.3, the remaining parameters remain unchanged, and the constant steady state solution $E_{1} = (0.6534, \; 4.8758)$, see Figure 8.

Biological interpretation 4: If additional alternative food for the natural enemies of Monochamus alternatus is added, it is indeed possible to control Monochamus alternatus faster and reduce the risk of outbreaks of pine wilt disease. If $\tau$ = 0.25, then the time required to maintain stability is 49 years. If $\tau$ = 0.35, then the time required to maintain stability is 60 years. However, it is worth noting that although we have provided additional alternative food for the natural enemies of Monochamus alternatus, it will indeed lead to a sharp increase in its number in the short term, but this does not mean that the long-term stable number of natural enemies of Monochamus alternatus will be improved. On the contrary, due to the limited food resources, the pressure of competition within the population will increase, which will lead to a decrease in its balanced number.

7.   Conclusions

In this paper, considering a Holling II-type functional response, we established a reaction-diffusion predator-prey model with both a Beverton-Holt-like alternative food source and pregnancy delay aimed at the control of pine wilt disease. We analyzed the existence and stability of the normal constant steady state solution, and the existence of a Hopf bifurcation near the normal constant steady state solution. We derived the normal form of the Hopf bifurcation for the system with a Beverton-Holt-like alternative food source and pregnancy delay based on the multiple time scales method. Through the data analysis of the relevant literature, we selected a set of appropriate numbers for numerical simulation. By analyzing the dynamic properties, we conclude that the combined effect of the Beverton-Holt-like alternative food source and the pregnancy delay can induce a stable homogeneous bifurcation periodic solution. The biological explanation of the prevention and control of pine wilt disease is also given. Through biological interpretation, this paper demonstrates the feasibility of maintaining the stable and controllable population density of natural enemies and Monochamus alternatus, so as to achieve effective protection of forestry resources and real green biological control.

Furthermore, this paper also demonstrates that the time delay of pregnancy is of great significance for the prevention and treatment of pine wilt disease, and the length of pregnancy time has a significant impact on the prevention and treatment of pine wilt disease. Therefore, we can take some intervention measures to reduce the negative effects caused by the long pregnancy of Monochamus alternatus' natural enemies. Specifically, the effect of introducing an alternative food source for natural enemies can be achieved by increasing the maximum average growth rate of natural enemies and reducing the predation rate of natural enemies. It is worth noting that by introducing an intervention strategy of an alternative food source, the system dynamics behavior can be significantly changed: when $\tau$ = 0.25, the system stability time is shortened from 960 years of the original model to 49 years, and when $\tau$ = 0.35, the stability time is reduced from 1750 to 60 years. This indicates that the intervention strategy through an alternative food source will significantly reduce the time required to control the population of Monochamus alternatus. At the same time, this will also reduce the possibility of an outbreak of the population of Monochamus alternatus and restore pine wilt disease to a controllable and stable state. Therefore, the mathematical model proposed in this paper has certain research value.

Author's contributions

The idea of this research was introduced by C. Wang and R. Yang. All authors contributed to the main results and numerical simulations.

Use of AI tools declaration

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

Acknowledgments

This research is supported by the Fundamental Research Funds for the Central Universities (No. 2572022DJ05), Postdoctoral Program of Heilongjiang Province (No. LBHQ21060), and the Northeast Forestry University College Student Innovation and Entrepreneurship Training Program Project (No. DCLXY-2025015).

Conflict of interest

The authors declare that they have no competing interests.