Bifurcation Analysis of a Hyphantria cunea-Chouioia cunea Yang Model with Nonlocal Competition for Biological Control

Xiaoxue Xia¹, Yujuan Gao², Yuting Ding^1*

¹Department of Mathematics, Northeast Forestry University, Harbin, China.
²Department of Economics and Management, Northeast Forestry University, Harbin, China.

DOI: 10.4236/ijmnta.2025.142002                

Abstract

Hyphantria cunea (H. cunea for short) causes serious forest pests and diseases, whereas its parasitic natural enemy, Chouioia cunea Yang (C. cunea for short) can effectively control the amount of H. cunea. Considering the eclosion time of pupae, nonlocal competition and biological control, we establish a biological control model with time delay and nonlocal competition. In this paper, we analyze the conditions for the existence of the Hopf bifurcation and the Turing bifurcation, then derive the normal form of Hopf bifurcation by using multiple time scales method. Finally, a group of real data is used for simulations. Comparing the models with and without nonlocal competition, we conclude that nonlocal competition can induce stable spatially inhomogeneous periodic solutions through Hopf bifurcation, while the model without nonlocal competition can induce spatially homogeneous periodic solutions through Hopf bifurcation. Comparing the models with and without biological control, we find that releasing limited amounts of C. cunea can quickly and effectively control the growth of H. cunea population. We also give the corresponding biological explanations and suggestions.

Keywords

H. cunea-C. Cunea Model, Nonlocal Competition, Hopf Bifurcation, Turing Bifurcation, Time Delay, Biological Control

1. Introduction

Hyphanthia cunea (also called Fall webworm, H. cunea for short) is a kind of quarantine pest that poses a danger to forestry and agriculture. H. cunea originated in North America, and it was accidentally introduced to central Europe and eastern Asia in the early 1940s [1]. Now it has spread to more than 30 countries, including South Africa, Argentina, China, Belarus, etc., and has been classified as an invasive and quarantine insect by organizations such as APPPC, COSA VE and EAEU EPPO (from Global Data Base: https://gd.eppo.int/taxon/HYPHCU/distribution). H. cunea has more than 600 recorded host plant species, which are widely distributed in several countries such as the United States, Japan, Korea, and China [2]. Due to its wide distribution across various host plants, H. cunea causes substantial economic losses globally. With the onset of global warming, H. cunea is increasingly spreading toward higher latitudes and elevations. Since it was first discovered in 1979 in Dandong City, Liaoning Province, China [3], the pest has spread to the Yangtze River basin, involving more than 10 provinces (cities and autonomous regions).

Both larvae and adults of the H. cunea will gnaw on leaves and even trunks. The species is highly reproductive, spreads rapidly, and causes significant damage to forests [4] [5]. H. cunea can threaten the survival of local species, diminish biodiversity and endanger ecological security. Figure 1 shows the condition of trees infected by the H. cunea. The leaves of the trees become transparent after being eaten by H. cunea (as in Figure 1(a) and Figure 1(b)), and the trunks of the trees turn yellow and wilt after being eaten (as in Figure 1(c) and Figure 1(d), Figure 1(c) is a partially enlarged image of Figure 1(d)) (The four images in Figure 1 are from: https://baike.baidu.com).

Figure 1. Leaves and trunks infected with the H. cunea.

The infestation of H. cunea mainly has the following control methods:

(1) Physical control method: This method primarily involves manual interventions such as cutting net screens, capturing adults, and digging pupae manually.

(2) Chemical pharmaceutical control method: This method mainly involves spraying susceptible trees with pesticides.

(3) Biological control method: This method involves releasing natural enemies of H. cunea, with the most effective and widely used being C. cunea.

The first two methods are effective, however, they either require a lot of manpower and material resources or can potentially damage the ecosystem. Biological control is the most environmentally friendly and economical method. In order to explore the biological control method of the H. cunea, Yang et al. [3] investigated the parasitic natural enemies of the H. cunea continuously from 1984 to 1986, and finally found a new natural enemy in the pupae of the H. cunea in Shanxi province. C. cunea is a dominant parasitic wasp with high parasitism and high fecundity. It is sensitive to lepidopteran pests such as H. cunea. C. cunea can insert its ovipositor into the pupae of insect pests, where it develops and grows, consuming all the nutrients and eventually causing the demise of H. cunea. H. cunea progresses through four life stages: egg, larva, pupa, and adult. The female wasps can lay eggs on the same day after eclosion under suitable temperature conditions. Female and male wasps mate within the host pupae and then emerge to seek a new host pupa for egg-laying. The release rate of C. cunea is determined by the density of H. cunea [6]. Figure 2 shows the reproductive mechanism of the coexistence of the H. cunea and the C. cunea in nature, with the red arrow indicating the artificial release of C. cunea. The left circle of Figure 2 shows the reproduction mechanism of the H. cunea, and the right circle shows the reproduction mechanism of the C. cunea. The intersection of the two circles is the H. cunea pupa.

Figure 2. Reproductive mechanism of H. cunea and C. cunea.

In recent years, some scholars have used ecological niche models (ENMs) to combine H. cunea occurrence records with environmental data. By creating relevant models to predict their ecological demand and potential geographic distribution. These experiments show that incorporating global data can substantially increase the forecasting precision of ENMs. This will provide valuable insights for the control of wedge nematode disease in China [7]. Several studies have investigated how temperature and season affect the behaviour of H. cunea populations and have created population dynamics models for its different life stages [8]. In addition, single population models with time delay were created [2]. Therefore, it is valuable to describe the dynamics of H. cunea and C. cunea by applying mathematical models.

In some sense, the richer dynamic properties of the model may be induced by time delay. Martin and Ruan [9] [10] proposed three types of time delay, these terms manifest within the prey-specific growth period, the predator response component of the predator equation, and the interaction component of the predator equation. Jiang et al. [11] [12] developed a predator-prey model incorporating time delay and explored a variety of dynamic properties associated with it. Time delay is crucial for the emergence of oscillatory and periodic behaviours. Considering that C. cunea can inhibit the growth rate of H. cunea, but this inhibition is not instantaneous. Therefore, the inclusion of time delay produces richer dynamic properties. The research on current time delay models has been a hot topic in mathematics, physics, biology and other fields. Furthermore, the introduction of nonlocal competition and its impact on model stability has been a subject of intense discussion in the literature over the past few decades. The works of Britton et al. [13] [14] have made significant contributions to the advancement of nonlocal competition. In their respective papers, they considered nonlocal competition and incorporated spatial convolution integrals into single species models. Britton et al. suggested that resource utilization in a given area is influenced by two factors: the density of the local population and the weighted average populations in surrounding areas. Song et al. [15]-[18] incorporated nonlocal competition into a predator-prey model with time delay and discovered that the interplay between nonlocal competition and time delay causes instability in the equilibrium.

The structure of the paper is as follows: Section 2 describes the process of establishing the mathematical model. Section 3 investigates the existence and stability of the positive equilibria and the existence of Hopf and Turing bifurcations near the equilibria. In Section 4, we utilize the multiple time scales method to derive the normal form of Hopf bifurcation for the model incorporating nonlocal competition. Section 5 presents numerical simulations that demonstrate our main results. Section 6 provides the conclusion.

2. Mathematical Modeling

In this section, we show the mathematical modeling process and explain it. Firstly, as noted in Ref. [19], the Lotka-Volterra and Nicholson-Bailey models were two primary frameworks traditionally employed to study predator-prey and parasitic-host interrelations. The Lotka-Volterra model results in neutrally stable cycles, while the Nicholson-Bailey system leads to diverging cycles and eventual extinction. Specifically, the Lotka-Volterra model can be expressed as follows:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{K}\right)-{\alpha }^{\prime }uv\mathrm{,}\\ \frac{\text{d}v}{\text{d}t}={\alpha }^{\prime }{\beta }^{\prime }uv-m^{\prime }v\mathrm{,}\end{array}$ (1)

here, $u$ and $v$ represent the prey and predator populations, respectively. The parameters $r$ and $K$ denote the highest per capita growth rate and environmental capacity of prey, respectively. ${\alpha }^{\prime }$ represents the capture rate of prey. ${\beta }^{\prime }$ represents the efficiency of converting ingested prey into new predators. $m^{\prime }$ represents the per capita mortality rate or the per capita nutritional needs for sustaining and replacing predators. Taking into account both biological and physical control, we establish the H. cunea-C. cunea model as follows:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=r\left(1-\frac{u}{K}\right)u-\alpha \beta uv-hu\mathrm{,}\\ \frac{\text{d}v}{\text{d}t}=\delta \alpha \beta uv-mv+cu\mathrm{,}\end{array}$ (2)

where $r$ and $K$ have the same meanings as in model (1). $u$ and $v$ represent the population density of H. cunea and C. cunea, respectively. $\alpha$ represents the probability of the H. cunea pupa being parasitized by the C. cunea. $\beta$ represents the proportionality coefficient. $h$ represents the mortality rate of H. cunea caused by physical control methods. $m$ represents the natural mortality rate of C. cunea. $\delta$ represents the probability of successful eclosion of the parasitized pupae. $c$ represents the release rate of C. cunea. This paper intends to study the method of reducing the density of H. cunea in more environmentally friendly ways. Therefore, the model (2) does not include the death rate of H. cunea due to pesticides. Reducing the population of H. cunea can be achieved by regularly releasing a quantified amount of C. cunea.

We also take into account the time delay, as the eclosion of H. cunea pupae is not instantaneous, regardless of parasitization. Since the consumption of environmental resources is affected by the flow of H. cunea, thus, the competition is not local. A plethora of researchers incorporated nonlocal competition by substituting the original $u/K$ with $\widehat{u}/K$ , where $\widehat{u}={\displaystyle {\int }_{\Omega }G\left(x\mathrm{,}y\right)u\left(y\mathrm{,}t\right)\text{d}y}$, which reflects intraspecific competition within prey populations. Note that $\Omega =\left(\mathrm{0,}l\pi \right)$ with $l$ is a positive parameter that measures the length of the habitat, and $\Omega$ represents the habitat volume of H. cunea. We assume that $G\left(x\mathrm{,}y\right)=1\backslash l\pi$ as $\left|\Omega \right|=l\pi$ , and we study the model (3) with Neumann boundary condition to facilitate the study of nonlocal competition. To sum up, we develop a H. cunea-C. cunea model with both nonlocal competition and time delay as follows:

$\{\begin{array}{l}\frac{\partial u\left(x,t\right)}{\partial t}=d_1\Delta u\left(x,t\right)+r\left(1-1/K{\displaystyle {\int }_{\Omega }G}\left(x,y\right)u\left(y,t\right)\text{d}y\right)u\left(x,t\right)\\ -\alpha \beta u\left(x,t-\tau \right)v\left(x,t-\tau \right)-hu\left(x,t\right),\mathrm{<mi>\; </mi>}x\in \Omega ,t>0,\\ \frac{\partial v\left(x,t\right)}{\partial t}=d_2\Delta v\left(x,t\right)+\delta \alpha \beta u\left(x,t-\tau \right)v\left(x,t-\tau \right)-mv\left(x,t\right)\\ +cu\left(x,t\right),\mathrm{<mi>\; </mi>}x\in \Omega ,t>0,\\ \frac{\partial u\left(x,t\right)}{\partial \overrightarrow{n}}=\frac{\partial v\left(x,t\right)}{\partial \overrightarrow{n}}=0,\mathrm{<mi>\; </mi>}x\in \partial \Omega ,t>0,\\ u\left(x,t\right)=u_0\left(x,t\right)\ge 0,\mathrm{<mi>\; </mi>}v\left(x,t\right)=v_0\left(x,t\right)\ge 0,\mathrm{<mi>\; </mi>}x\in \overline{\Omega },\mathrm{<mi>\; </mi>}t\in \left[-\tau ,0\right],\end{array}$ (3)

where $d_1>0$ , $d_2>0$ represent the diffusion coefficients of H. cunea and C. cunea, respectively. The other parameters $r$ , $K$ , $\alpha$ , $\beta$ , $h$ , $m$ , $\delta$ and $c$ represent the same meanings as the parameters in model (2), and all of them are non-negative.

Next, we will analyze the dynamic properties of model (3).

3. Existence and Stability of Equilibria and Existence of Bifurcations

In this section, we will consider the existence and stability of the equilibria for model (3), as well as the existence of Hopf and Turing bifurcations near the equilibria.

We know that model (3) has a trivial equilibrium $E_0=\left(\mathrm{0,0}\right)$ , and we make the following hypothesis:

(H₁) $r>h,$

if (H₁) holds, then the model (3) has a positive equilibrium $E_{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ , where

$\begin{array}{l}{u}^{\mathrm{<mo>*</mo>}}=\frac{m{v}^{\mathrm{<mo>*</mo>}}}{\delta \alpha \beta v^{\mathrm{<mo>*</mo>}}+c}\mathrm{,}\\ v^{\mathrm{<mo>*</mo>}}=\frac{I+\sqrt{I^2+4k^2\delta {\alpha }^2{\beta }^{2}c\left(r-h\right)}}{2k\delta {\alpha }^2{\beta }^2}\mathrm{,}\end{array}$

with $I=kr\delta \alpha \beta -kh\delta \alpha \beta +kc\alpha \beta +rm$ , it is obvious that $u^{*}>0$ , and $v^{\mathrm{<mo>*</mo>}}>0$ when $r>h$ .

Remark 1. The death rates from physical insect control methods are generally low, such as traps, barriers, electric shocks, etc., which are often non-lethal. The objective of physical control typically aims to diminish insect populations rather than inducing immediate widespread death.

3.1. Existence of Hopf Bifurcation

Next, we analyze the stability of the equilibria for the model (3), and model (3) can be shown as follows:

$\frac{\partial }{\partial t}\left(\begin{array}{c}u\left(x\mathrm{,}t\right)\\ v\left(x\mathrm{,}t\right)\end{array}\right)=D\left(\begin{array}{c}\Delta u\left(x\mathrm{,}t\right)\\ \Delta v\left(x\mathrm{,}t\right)\end{array}\right)+J_1\left(\begin{array}{c}u\left(x\mathrm{,}t\right)\\ v\left(x\mathrm{,}t\right)\end{array}\right)+J_2\left(\begin{array}{c}u\left(x\mathrm{,}t-\tau \right)\\ v\left(x\mathrm{,}t-\tau \right)\end{array}\right)+J_3\left(\begin{array}{c}\widehat{u}\left(x\mathrm{,}t\right)\\ \widehat{v}\left(x\mathrm{,}t\right)\end{array}\right)\mathrm{,}$ (4)

where

$\begin{array}{l}D=\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right),J_1=\left(\begin{array}{cc}{a}_1& 0\\ c& -m\end{array}\right),\\ J_2=\left(\begin{array}{cc}{b}_1& b_2\\ b_3& b_4\end{array}\right),J_3=\left(\begin{array}{cc}{c}_1& 0\\ 0& 0\end{array}\right),\end{array}$

with $a_1=r-h-r{u}^{\mathrm{<mo>*</mo>}}/K$ , $b_1=-\alpha \beta v^{\mathrm{<mo>*</mo>}}$ , $b_2=-\alpha \beta u^{\mathrm{<mo>*</mo>}}$ , $b_3=\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}$ , $b_4=\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}$ , $c_1=-r{u}^{\mathrm{<mo>*</mo>}}/K$ . For trivial equilibrium $E_0$ , the characteristic equation corresponding to model (3) at $E_0$ has two roots: ${\lambda }_1=2r-2h>0$ , ${\lambda }_2=-2m<0$ . Therefore, the trivial equilibrium $E_0$ is a saddle point.

The characteristic equation corresponding to model (3) at $E_{\mathrm{<mo>*</mo>}}$ is as follows:

${\lambda }^2+\left(A_n-F{\text{e}}^{-\lambda \tau }\right)\lambda +B_n+T_n{\text{e}}^{-\lambda \tau }=0,n=0,1,2,3,\cdots ,$ (5)

where $F=b_1+b_4$ , when $n=0$ , $A_0=m-a_1-c_1$ , $B_0=-a_{1}m-c_{1}m$ , $T_0=a_1{b}_4-b_{1}m-b_{2}c+c_1{b}_4$ . When $n=\mathrm{1,2,3,}\cdots$ , $A_n=\frac{n^2}{l^2}{d}_1+\frac{n^2}{l^2}{d}_2+m-a_1$ , $B_n=\frac{n^4}{l^4}{d}_1{d}_2+\frac{n^2}{l^2}\left(m{d}_1-a_1{d}_2\right)-a_{1}m$ , $T_n=-\frac{n^2}{l^2}\left(d_1{b}_1+d_2{b}_1\right)+a_1{b}_4-b_{1}m-b_{2}c$ .

When $\tau =0$ , characteristic Eq. (5) becomes

${\lambda }^2+\left(A_n-F\right)\lambda +B_n+T_n=\mathrm{0,<mi>\; </mi>}n=\mathrm{0,1,2,3,}\cdots \mathrm{.}$

Next, we make the following hypothesis:

(H₂) $m-a_1-c_1-b_1-b_4>0,\mathrm{<mi>\; </mi><mi>\; </mi>}{a}_1{b}_4+c_1{b}_4-a_{1}m-c_{1}m-b_{1}m-b_{2}c>0,$

clearly, $A_n-F>A_0-F=m-a_1-c_1-b_1-b_4>0$ , and $B_n+T_n>B_0+T_0=a_1{b}_4+c_1{b}_4-a_{1}m-c_{1}m-b_{1}m-b_{2}c>0$ . This is because both $A_n-F$ and $B_n+T_n$ are monotonically increasing for $n>0$ . To sum up, if (H₂) holds, $A_n-F>0$ and $B_n+T_n>0$ hold for any $n=\mathrm{1,2,3,}\cdots$ .

Based on the above analysis, we give the following theorem:

Theorem 1. If hypothesis (H₁) holds, $E_{\mathrm{<mo>*</mo>}}$ is a positive equilibrium, and $E_0$ is a saddle point. If hypotheses (H₁) and (H₂) hold, the positive equilibrium $E_{\mathrm{<mo>*</mo>}}$ of model (3) for $\tau =0$ is locally asymptotically stable for any integer $n\in \mathbb{N}$ .

Next, we explore the existence of Hopf bifurcation near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ for $\tau >0$ . Let $\lambda =\pm \text{i}\omega$ $\left(\omega >0\right)$ are the characteristic roots of Eq. (5). Separating the real part and imaginary part, we obtain

$\{\begin{array}{l}-{\omega }^2+B_n=F\omega \sin \omega \tau -T_n\cos \omega \tau \mathrm{,}\\ A_n\omega =F\omega \cos \omega \tau +T_n\sin \omega \tau \mathrm{,}\end{array}$ (6)

that is,

${\omega }^4+\left(-2B_n+A_n^2-F^2\right){\omega }^2+B_n^2-T_n^2=\mathrm{0,<mi>\; </mi>}n=\mathrm{0,1,2,3,}\cdots \mathrm{,}$ (7)

letting $z={\omega }^2$ , then we get

$h\left(z\right)=z^2+\left(-2B_n+A_n^2-F^2\right)z+B_n^2-T_n^2=0,\mathrm{<mi>\; </mi>}n=0,1,2,3,\cdots .$ (8)

By Eq. (6), we can get

$\{\begin{array}{l}\cos \omega \tau =\frac{-B_n{T}_n+T_n{\omega }^2+A_{n}F{\omega }^2}{F^2{\omega }^2+T_n^2}=E_n,\\ \sin \omega \tau =\frac{-{\omega }^{3}F+B_{n}F\omega +T_n{A}_n\omega }{F^2{\omega }^2+T_n^2}=G_n,\end{array}$

then for $j=\mathrm{0,1,2,3,}\cdots$ and $n=\mathrm{0,1,2,3,}\cdots$ , we get

${\tau }_n^{\left(j\right)}=\{\begin{array}{l}\frac{\arccos E_n+2j\pi }{\omega },G_n>0,\\ \frac{\pi +\arccos \left(-E_n\right)+2j\pi }{\omega },G_n<0,\end{array}$ (9)

it is clear that ${\tau }_n^{\left(j\right)}$ monotonically increases with $j$ . By calculation, the transversality condition

${\mathrm{Re}\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)|}_{\tau ={\tau }_n^{\left(j\right)}}=\frac{2{\omega }^2-2B_n+A_n^2-F^2}{F^2{\omega }^2+T_n^2}=\frac{h^{\prime }\left(z\right)}{F^2{\omega }^2+T_n^2}$

according to Eq. (8), it is clearly that $h^{\prime }\left(z\right)\ne 0$ for any $n\ge 0$ . Thus

${\text{Re}\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)|}_{\tau \mathrm{<mo>=</mo>}{\tau }_n^{\left(j\right)}}\ne 0$ holds, then we denote

${\tau }_c=\min \left\{{\tau }_n^{\left(j\right)}\mathrm{<mo>|</mo>}n\in \mathbb{N}\right\}\mathrm{,}$ (10)

where ${\tau }_n^{\left(j\right)}$ is given in Eq. (9).

Denote that $S_1$ and $S_2$ are sets of integers, which may be finite or infinite,

$S_1=\left\{n|B_n^2-T_n^2<0,n\in \mathbb{N}\right\},$

$S_2=\left\{n|B_n^2-T_n^2>0,\mathrm{<mi>\; </mi>}-2B_n+A_n^2-F^2<0,\mathrm{<mi>\; </mi>}\Delta >0,\mathrm{<mi>\; </mi>}n\in \mathbb{N}\right\},$

where $\Delta ={\left(A_n-F\right)}^2-4\left(B_n+T_n\right)$ . Combining the analysis above, we derive the following theorem.

Theorem 2. Assuming that hypotheses (H₁) and (H₂) are satisfied, and ${\tau }_n^{\left(j\right)}$ is given in Eq. (9). Thus, we have following conclusions:

(1) If $S_1=\varnothing$ and $S_2=\varnothing$ , then equilibrium $E_{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is locally asymptotically stable for $\tau \ge 0$ .

(2) If $S_1\ne \varnothing$ and $S_2=\varnothing$ , Eq. (8) has a unique positive root $z_{n_l\mathrm{,1}}$ for some $n_l\in S_1$ , then the equilibrium $E_{\mathrm{<mo>*</mo>}}$ of model (3) is locally asymptotically stable for $0<\tau <{\tau }_c$ and becomes unstable for $\tau >{\tau }_c$ , where ${\tau }_c=\min \left\{{\tau }_{n_l}^{\left(0\right)}\mathrm{<mo>|</mo>}{n}_l\in S_1\right\}$ . Additionally, model (3) undergoes $n_l$ -mode Hopf bifurcation near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ for $\tau ={\tau }_{n_l}^{\left(j\right)}$ with $n_l\in S_1$ and $j=\mathrm{0,1,2,3,}\cdots$ .

(3) If $S_1=\varnothing$ and $S_2\ne \varnothing$ , Eq. (8) has two positive roots $z_{n_s\mathrm{,}i}$ $\left(i=1,2\right)$ for some $n_s\in S_2$ , then the equilibrium $E_{\mathrm{<mo>*</mo>}}$ of the model (3) is locally asymptotically stable for $0<\tau <{\tau }_c$ , and there may be a switch of stability for $\tau >{\tau }_c$ in the model (3), where ${\tau }_c=\min \left\{{\tau }_{n_{s,1}}^{\left(0\right)}|n_s\in S_2\right\}$ . Additionally, model (3) undergoes $n_s$ -mode Hopf bifurcation near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ for $\tau ={\tau }_{n_s}^{\left(j\right)}$ with $n_s\in S_2$ and $j=\mathrm{0,1,2,3,}\cdots$ .

(4) If $S_1\ne \varnothing$ and $S_2\ne \varnothing$ , Eq. (8) has positive roots $z_{n_l\mathrm{,1}}$ and $z_{n_s\mathrm{,}i}$ $\left(i=\mathrm{1,2}\right)$ for some $n_l\in S_1$ , $n_s\in S_2$ , respectively. Then the equilibrium $E_{\mathrm{<mo>*</mo>}}$ of the model (3) is locally asymptotically stable for $0<\tau <{\tau }_c$ , and there may be a switch of stability for $\tau \mathrm{<mo>></mo>}{\tau }_c$ in the model (3), where ${\tau }_c=\min \left\{{\tau }_{n_T}^{\left(0\right)}|n_T\in S_1\text{or}{n}_T\in S_2\right\}$ . Additionally, model (3) undergoes $n_T$ -mode Hopf bifurcation near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ for $\tau ={\tau }_{n_T}^{\left(j\right)}$ with $n_T\in S_1$ , or $n_T\in S_2$ and $j=\mathrm{0,1,2,3,}\cdots$ .

3.2. Existence of Turing Bifurcation

In this subsection, we analyze the conditions for the existence of Turing bifurcation. Suppose $\lambda =0$ is the root of characteristic Eq. (5), then we get expression:

$B_n+T_n=\frac{n^4}{l^4}{d}_1{d}_2+\frac{n^2}{l^2}{d}_1\left(m-b_4\right)+J=\mathrm{0,}$ (11)

where $J=a_1{b}_4-c{b}_2=\delta {\alpha }^2{\beta }^2{u}^{*}{v}^{*}+c\alpha \beta u^{*}>0$ due to $\delta \mathrm{,}\alpha \mathrm{,}\beta \mathrm{,}c$ and $u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}$ are positive.

Letting $k=\frac{n^2}{l^2}$ , by simple calculation from Eq. (11), we get

$d_1=\frac{-J}{k^2{d}_2+k\left(m-b_4\right)}\mathrm{,}$

$d_2=\frac{-J-k{d}_1\left(m-b_4\right)}{k^2{d}_1}\mathrm{.}$

According to its biological significance that $d_1>0$ and $d_2>0$ , which means $k^2{d}_2+k\left(m-b_4\right)<0$ and $-J-k{d}_1\left(m-b_4\right)>0$ , then we have

$\sqrt{k}=\frac{n}{l}<\sqrt{\frac{b_4-m}{d_2}},$

$\sqrt{k}=\frac{n}{l}<\sqrt{\frac{-J}{d_1\left(m-b_4\right)}}.$

Then we making the following hypothesis:

(H₃) $b_4-m>0.$

We denote that $S_3=\left\{n|0<n<\min \left\{l\sqrt{\frac{b_4-m}{d_2}},\sqrt{\frac{-J}{d_1\left(m-b_4\right)}}\right\},n\in \mathbb{N}\right\}$ , and we

have the following theorem.

Theorem 3. If (H₃) is true and a positive integer $n_p\in S_3$ can be found, then the model (3) undergoes the $n_p$ -mode Turing bifurcation. If $S_3=\varnothing$ , then model (3) does not undergo Turing bifurcation.

4. Stability and Direction of Hopf Bifurcation

In this section, we derive the normal form of the Hopf bifurcation for the model (3) by using the multiple time scales method as referenced in ref. [20]-[22]. When $\tau ={\tau }_c$ , the characteristic equation Eq. (5) has a pair of pure imaginary roots $\lambda =\pm \text{i}\omega$ , at which model (3) undergoes $n_l$ , $n_s$ or $n_T$ -mode Hopf bifurcation at equilibrium $E_{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ . In the subsequent derivation, we simplify notation by using $n$ instead of $n_l$ , $n_s$ or $n_T$ . We consider the time delay $\tau$ serves as the bifurcation parameter in this paper, where $\tau ={\tau }_c+\epsilon \mu$ , with $\tau ={\tau }_c$ being the critical value for Hopf bifurcation. The parameter $\epsilon$ is a dimensionless scaling factor, and $\mu$ is a perturbation parameter. We also introduce the transformations: $t\mapsto t/\tau$ , $u\mapsto u-u^{\mathrm{<mo>*</mo>}}$ , $v\mapsto v-v^{\mathrm{<mo>*</mo>}}$ , then model (3) is transformed into

$\{\begin{array}{l}\frac{\partial u}{\partial t}=\tau [d_1\Delta u+r\left(1-\frac{\widehat{u}+u^{\mathrm{<mo>*</mo>}}}{K}\right)u-\alpha \beta u\left(t-1\right)v\left(t-1\right)\\ \begin{array}{c}\\ \end{array}-\alpha \beta v^{\mathrm{<mo>*</mo>}}u\left(t-1\right)-\alpha \beta u^{\mathrm{<mo>*</mo>}}v\left(t-1\right)-hu]\mathrm{,}\\ \frac{\partial v}{\partial t}=\tau [d_2\Delta v-mv+\delta \alpha \beta u\left(t-1\right)v\left(t-1\right)+\delta \alpha \beta v^{\mathrm{<mo>*</mo>}}u\left(t-1\right)\\ +\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}v\left(t-1\right)+cu]\mathrm{.}\end{array}$ (12)

Let $g$ be the eigenvector of the linear operator corresponding to the eigenvalue $\lambda =\text{i}\omega \tau$ of the linearized system of Eq. (12) for equilibrium $E^{\mathrm{<mo>*</mo>}}$ , and let $g^{\mathrm{<mo>*</mo>}}$ be the normalized eigenvector of the adjoint operator of the linear operator corresponding to the eigenvalue $\lambda =-\text{i}\omega \tau$ , satisfying

$\langle g^{\mathrm{<mo>*</mo>}}\mathrm{,}g\rangle ={\overline{g^{\mathrm{<mo>*</mo>}}}}^{\text{T}}g=1$. By calculation, it can be found

$g={\left(g_{11}\mathrm{,}{g}_{12}\right)}^{\text{T}}={\left(\mathrm{1,}\frac{\text{i}\omega +n^2{d}_1-a_1-b_1{\text{e}}^{-\text{i}\omega \tau }-c_1+\frac{r{u}^{\mathrm{<mo>*</mo>}}}{K}{x}_n}{b_2{\text{e}}^{-\text{i}\omega \tau }}\right)}^{\text{T}}\mathrm{,}$

$g^{\mathrm{<mo>*</mo>}}=d^{\mathrm{<mo>*</mo>}}{\left(g_{21}\mathrm{,}{g}_{22}\right)}^{\text{T}}=d^{\mathrm{<mo>*</mo>}}{\left(\frac{b_2{\text{e}}^{-\text{i}\omega \tau }}{-\text{i}\omega +n^2{d}_2-a_1-b_1{\text{e}}^{\text{i}\omega \tau }-c_1}\mathrm{,1}\right)}^{\text{T}}\mathrm{,}$ (13)

where $d^{\mathrm{<mo>*</mo>}}={\left(g_{21}\overline{g_{11}}+\overline{g_{22}}{g}_{12}\right)}^{-1}$ , and $x_n=\{\begin{array}{ll}\mathrm{1,}\hfill & \text{if}n=0\hfill \\ \mathrm{0,}\hfill & \text{if}n\in {\mathbb{N}}_{+}\hfill \end{array}$ .

According to the multiple scales method, we can assume the solution of Eq. (12) as follows:

$u\left(x,t\right)={\displaystyle \sum _{k=0}^{\infty }}{\epsilon }^k{u}_k\left(x,t_0,t_1,t_2,\cdots \right),$

$v\left(x,t\right)={\displaystyle \sum _{k=0}^{\infty }}{\epsilon }^k{v}_k\left(x,t_0,t_1,t_2,\cdots \right),$ (14)

where $t_i={\epsilon }^{i}t\mathrm{,}i=\mathrm{1,2,3,}\cdots$ , specially,

$u_{j,{\tau }_c}\left(x,t-1\right)=u_j\left(x,t_0-1,t_1,t_2,\cdots \right),$

$v_{j,{\tau }_c}\left(x,t-1\right)=v_j\left(x,t_0-1,t_1,t_2,\cdots \right),j=1,2,3,\cdots .$ (15)

By substituting (14) and (15) into (12), we compare the coefficients of $\epsilon \mathrm{,}{\epsilon }^2$ and ${\epsilon }^3$ on either side of the equation, resulting in the following expressions:

$\{\begin{array}{l}{D}_0{u}_1-{\tau }_c{d}_1\Delta u_1-{\tau }_{c}r{u}_1+\frac{r{\tau }_c{u}^{*}{u}_1}{K}+{\tau }_{c}h{u}_1+{\tau }_c\alpha \beta v^{*}{u}_{1{\tau }_c}+{\tau }_c\alpha \beta u^{*}{v}_{1{\tau }_c}=0,\\ D_0{v}_1-{\tau }_c{d}_2\Delta v_1+{\tau }_{c}m{v}_1-{\tau }_{c}c{u}_1-{\tau }_c\delta \alpha \beta v^{*}{u}_{1{\tau }_c}-{\tau }_c\delta \alpha \beta u^{*}{v}_{1{\tau }_c}=0.\end{array}$ (16)

$\{\begin{array}{l}{D}_0{u}_2-{\tau }_c{d}_1\Delta u_2-{\tau }_{c}r{u}_2+\frac{r{\tau }_c{u}^{*}{u}_2}{K}+{\tau }_{c}h{u}_2+{\tau }_c\alpha \beta v^{*}{u}_{2{\tau }_c}+{\tau }_c\alpha \beta u^{*}{v}_{2{\tau }_c}\\ =\mu d_1\Delta u_1+\mu r{u}_1-\frac{\mu r{u}^{*}{u}_1}{K}-\mu h{u}_1-\mu \alpha \beta v^{*}{u}_{1{\tau }_c}-\mu \alpha \beta u_{*}{v}_{1{\tau }_c}-\frac{{\tau }_{c}r{u}_1^2}{K}\\ -{\tau }_c\alpha \beta u_{1{\tau }_c}{v}_{1{\tau }_c}+{\tau }_c\alpha \beta v^{*}{D}_1{u}_{1{\tau }_c}+{\tau }_c\alpha \beta u^{*}{D}_1{v}_{1{\tau }_c}-D_1{u}_1,\\ D_0{v}_2-{\tau }_c{d}_2\Delta v_2+{\tau }_{c}m{v}_2-{\tau }_{c}c{u}_2-{\tau }_c\delta \alpha \beta v^{*}{u}_{2{\tau }_c}-{\tau }_c\delta \alpha \beta u^{*}{v}_{2{\tau }_c}\\ =\mu d_2\Delta v_1-\mu m{v}_1+\mu \delta \alpha \beta v^{*}{u}_{1{\tau }_c}+\mu \delta \alpha \beta v^{*}{v}_{1{\tau }_c}+{\tau }_c\delta \alpha \beta u_{1{\tau }_c}{v}_{1{\tau }_c}+\mu c{u}_1\\ -{\tau }_c\delta \alpha \beta v^{*}D-1u_{1{\tau }_c}-{\tau }_c\delta \alpha \beta u^{*}{D}_1{v}_{1{\tau }_c}-D_1{v}_1.\end{array}$ (17)

$\{\begin{array}{l}{D}_0{u}_3-{\tau }_c{d}_1\Delta u_3-{\tau }_{c}r{u}_3+\frac{r{\tau }_c{u}^{*}{u}_3}{K}+{\tau }_{c}h{u}_3+{\tau }_c\alpha \beta v^{*}{u}_{3{\tau }_c}+{\tau }_c\alpha \beta u^{*}{v}_{3{\tau }_c}\\ =-\mu \alpha \beta v^{*}{u}_{2{\tau }_c}+\mu \alpha \beta v^{*}{D}_1{u}_{1{\tau }_c}+\mu \alpha \beta u^{*}{D}_1{v}_{1{\tau }_c}-\frac{2{\tau }_{c}r{u}_1{u}_2}{K}\\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}-{\tau }_c\alpha \beta u_{1{\tau }_c}{v}_{2{\tau }_c}-{\tau }_c\alpha \beta u_{2{\tau }_c}{v}_{1{\tau }_c}+{\tau }_c\alpha \beta D_2{u}_{1{\tau }_c}+{\tau }_c\alpha \beta v^{*}{D}_1{u}_{1{\tau }_c}\\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}+{\tau }_c\alpha \beta u^{*}{D}_2{v}_{1{\tau }_c}+{\tau }_c\alpha \beta D_1{v}_{1{\tau }_c}-D_2{u}_1-D_1{u}_2+o\left(\mu \right),\\ D_0{v}_3-{\tau }_c{d}_2\Delta v_3+{\tau }_{c}m{v}_3-{\tau }_{c}c{u}_3-{\tau }_c\delta \alpha \beta v^{*}{u}_{3{\tau }_c}-{\tau }_c\delta \alpha \beta u^{*}{v}_{3{\tau }_c}\\ =\mu \delta \alpha \beta u_{1{\tau }_c}{v}_{1{\tau }_c}+\mu c{u}_2-\mu \delta \alpha \beta v^{*}{D}_1{u}_{1{\tau }_c}-\mu \delta \alpha \beta u^{*}{D}_1{v}_{1{\tau }_c}\\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}+{\tau }_c\delta \alpha \beta u_{1{\tau }_c}{v}_{2{\tau }_c}+{\tau }_c\delta \alpha \beta u_{2{\tau }_c}{v}_{1{\tau }_c}-{\tau }_c\delta \alpha \beta v^{*}{D}_2{u}_{1{\tau }_c}-{\tau }_c\delta \alpha \beta v^{*}{D}_1{u}_{1{\tau }_c}\\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}-{\tau }_c\delta \alpha \beta u^{*}{D}_2{v}_{1{\tau }_c}-{\tau }_c\delta \alpha \beta u^{*}{D}_1{v}_{1{\tau }_c}-D_2{v}_1-D_1{v}_2+o\left(\mu \right).\end{array}$ (18)

Referring to the notations in Refs. [20] [21], then Eq. (16) has the solution with following form:

$\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right)=Gg{\text{e}}^{\text{i}\omega {\tau }_c{T}_0}+\overline{G}\overline{g}{\text{e}}^{-\text{i}\omega {\tau }_c{T}_0}\mathrm{.}$ (19)

Therefore, by substituting Eq. (19) into the right part of Eq. (17), then we define the coefficient vector of ${\text{e}}^{\text{i}\omega {\tau }_c{T}_0}$ as $m_1$ . Therefore, according to the solvability condition, we have $\langle g^{\mathrm{<mo>*</mo>}}\mathrm{,}\left(m_1\mathrm{,}\cos nx\right)\rangle =0$ . Then we can solve $\frac{\partial G}{\partial T_1}$ , namely,

$\frac{\partial G}{\partial T_1}=M\mu G\mathrm{,}$ (20)

where

$M=\frac{\overline{g_{21}}{\varrho }_1+\overline{g_{22}}{\varrho }_2}{\overline{g_{21}}{\varrho }_3+\overline{g_{22}}{\varrho }_4}\mathrm{,}$

with

$\{\begin{array}{l}{\varrho }_1=\left(-d_1\frac{n^2}{l^2}+r-\frac{r{u}^{\mathrm{<mo>*</mo>}}}{K}-h-\alpha \beta v^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}\right)g_{11}-\alpha \beta u^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}{g}_{12}\mathrm{,}\\ {\varrho }_2=\left(-d_2\frac{n^2}{l^2}-m+\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}\right)g_{12}+\left(\delta \alpha \beta v^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}+c\right)g_{11}\mathrm{,}\\ {\varrho }_3={\tau }_c\alpha \beta v^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}{g}_{11}-{\tau }_c\alpha \beta u^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}{g}_{12}-g_{11}\mathrm{,}\\ {\varrho }_4={\tau }_c\delta \alpha \beta v^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}{g}_{12}+{\tau }_c\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}{\text{e}}^{-\text{i}\omega {\tau }_c}{g}_{12}+g_{12}\mathrm{.}\end{array}$

Suppose that the solution of (17) satisfies following form:

$\begin{array}{l}{u}_2={\displaystyle \sum _{k=0}^{\infty }}\left({\eta }_{0k}G\overline{G}+{\eta }_{1k}{\text{e}}^{2i\omega {\tau }_c{T}_0}{G}^2+\overline{{\eta }_{1k}}{\text{e}}^{-2\text{i}\omega {\tau }_c{T}_0}\overline{G^2}\right)\cos \frac{kx}{l}\mathrm{,}\\ v_2={\displaystyle \sum _{k=0}^{\infty }}\left({\xi }_{0k}G\overline{G}+{\xi }_{1k}{\text{e}}^{2\text{i}\omega {\tau }_c{T}_0}{G}^2+\overline{{\xi }_{1k}}{\text{e}}^{-2\text{i}\omega {\tau }_c{T}_0}\overline{G^2}\right)\cos \frac{kx}{l}\mathrm{,}\end{array}$ (21)

we substitute (21) into the left end of (17), it can be obtained by the method of undetermined coefficients, where

${\eta }_{0k}=\frac{{\zeta }_1{\zeta }_2-{\zeta }_3{\zeta }_4}{{\zeta }_6{\zeta }_2}\mathrm{,<mi>\; </mi>}{\eta }_{1k}=\frac{Q_1{Q}_2-Q_3{Q}_4}{Q_5{Q}_2-Q_3{Q}_4}\mathrm{,}$

${\xi }_{0k}=\frac{{\zeta }_5{\zeta }_1+{\zeta }_3{\zeta }_6}{{\zeta }_2{\zeta }_6+{\zeta }_4{\zeta }_5}\mathrm{,<mi>\; </mi>}{\xi }_{1k}=\frac{Q_1{Q}_6-Q_5{Q}_3}{Q_5{Q}_2+Q_3{Q}_4{Q}_6}\mathrm{,}$ with

$\{\begin{array}{l}{Q}_1=-\frac{{\tau }_{c}r}{K}{g}_{11}^2-{\tau }_c\alpha \beta {\text{e}}^{-2\text{i}\omega {\tau }_c}{g}_{11}{g}_{12}\mathrm{,}\\ Q_2=2\text{i}\omega {\tau }_c+{\tau }_c{d}_2\frac{k^2}{l^2}+{\tau }_{c}m-{\tau }_c\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}{\text{e}}^{-2\text{i}\omega {\tau }_c}\mathrm{,}\\ Q_3={\tau }_c\delta \alpha \beta g_{11}{g}_{12}{\text{e}}^{-2\text{i}\omega {\tau }_c}\mathrm{,}\\ Q_4={\tau }_c\alpha \beta u^{\mathrm{<mo>*</mo>}}{\text{e}}^{-2\text{i}\omega {\tau }_c}\mathrm{,}\\ Q_5=2\text{i}\omega {\tau }_c+{\tau }_c{d}_1\frac{k^2}{l^2}+{\tau }_c\alpha \beta v^{\mathrm{<mo>*</mo>}}{\text{e}}^{-2\text{i}\omega {\tau }_c}-{\tau }_{c}r+\frac{{\tau }_{c}r{u}^{\mathrm{<mo>*</mo>}}}{K}+{\tau }_{c}h\mathrm{,}\\ Q_6=-{\tau }_c\delta \alpha \beta v^{\mathrm{<mo>*</mo>}}{\text{e}}^{-2\text{i}\omega {\tau }_c}-c{\tau }_c\mathrm{,}\end{array}$ and

$\{\begin{array}{l}{\zeta }_1=-\frac{2r{\tau }_c}{K}{g}_{11}\overline{g_{11}}-{\tau }_c\alpha \beta g_{11}\overline{g_{12}}-{\tau }_c\alpha \beta g_{12}\overline{g_{11}},\\ {\zeta }_2={\tau }_c{d}_2\frac{k^2}{l^2}+{\tau }_{c}m-{\tau }_c\delta \alpha \beta u^{*},\\ {\zeta }_3={\tau }_c\delta \alpha \beta g_{11}\overline{g_{12}},\\ {\zeta }_4={\tau }_c\alpha \beta u^{*},\\ {\zeta }_5=-{\tau }_c\delta \alpha \beta v^{*}-c{\tau }_c,\\ {\zeta }_6={\tau }_c{d}_1\frac{k^2}{l^2}+{\tau }_c\alpha \beta v^{*}-{\tau }_{c}r+\frac{{\tau }_{c}r{u}^{*}}{K}+{\tau }_{c}h.\end{array}$

Next, substitute (19) and (21) into (18), then compare the coefficients, and denoting the coefficient vector of ${\text{e}}^{\text{i}\omega {\tau }_c{T}_0}$ as $m_2$ . Therefore, according to the solvability condition, we have $\langle g^{\mathrm{<mo>*</mo>}}\mathrm{,}\left(m_2\mathrm{,}\cos nx\right)\rangle =0$ . Specially, $\mu$ is a perturbation parameter, and ${\mu }^2$ has minimal influence of small unfolding parameters. Thus, we have

$\frac{\partial G}{\partial T_2}=X{G}^2\overline{G}\mathrm{,}$ (22)

where

$X=\frac{\left[\overline{g_{21}}{\iota }_1-\overline{g_{22}}{\iota }_2\right]}{\left[\overline{g_{21}}{\iota }_3+\overline{g_{22}}{\iota }_4\right]D_k}\mathrm{,}$

with

$\{\begin{array}{l}{\iota }_1=\frac{2r{\tau }_c}{K}{g}_{11}{\eta }_{00}{C}_0+\frac{2r{\tau }_c}{K}\overline{g_{11}}{\eta }_{10}{C}_0+{\tau }_c\alpha \beta g_{12}\left({\eta }_{00}{C}_0+{\eta }_{02n}{C}_{2n}\right)\\ +{\tau }_c\alpha \beta \overline{g_{12}}\left({\eta }_{10}{C}_0+{\eta }_{12n}{C}_{2n}\right)+{\tau }_c\alpha \beta g_{11}{\text{e}}^{-\text{i}\omega {\tau }_c}\left({\xi }_{00}{C}_0+{\xi }_{02n}{C}_{2n}\right)\\ +{\tau }_c\alpha \beta \overline{g_{11}}{\text{e}}^{-\text{i}\omega {\tau }_c}\left({\xi }_{10}{C}_0+{\xi }_{12n}{C}_{2n}\right)\mathrm{,}\\ {\iota }_2={\tau }_c\delta \alpha \beta g_{11}{\text{e}}^{-\text{i}\omega {\tau }_c}\left({\xi }_{00}{C}_0+{\xi }_{02n}{C}_{2n}\right)+{\tau }_c\delta \alpha \beta \overline{g_{11}}{\text{e}}^{-\text{i}\omega {\tau }_c}\left({\xi }_{10}{C}_0+{\xi }_{12n}{C}_{2n}\right)\\ +{\tau }_c\delta \alpha \beta g_{12}{\text{e}}^{-\text{i}\omega {\tau }_c}\left({\eta }_{00}{C}_0+{\eta }_{02n}{C}_{2n}\right)+{\tau }_c\delta \alpha \beta \overline{g_{12}}{\text{e}}^{-\text{i}\omega {\tau }_c}\left({\eta }_{10}{C}_0+{\eta }_{12n}{C}_{2n}\right)\mathrm{,}\\ {\iota }_3={\tau }_c\alpha \beta v^{\mathrm{<mo>*</mo>}}{g}_{11}{\text{e}}^{-\text{i}\omega {\tau }_c}+{\tau }_c\alpha \beta u^{\mathrm{<mo>*</mo>}}{g}_{12}{\text{e}}^{-\text{i}\omega {\tau }_c}-g_{11}\mathrm{,}\\ {\iota }_4=-{\tau }_c\delta \alpha \beta v^{\mathrm{<mo>*</mo>}}{g}_{11}{\text{e}}^{-\text{i}\omega {\tau }_c}-{\tau }_c\delta \alpha \beta u^{\mathrm{<mo>*</mo>}}{g}_{12}{\text{e}}^{-\text{i}\omega {\tau }_c}-g_{12}\mathrm{,}\end{array}$

and

$C_k={\displaystyle {\int }_0^{l\pi }{\cos }^2\frac{nx}{l}\text{d}x\cos \frac{kx}{l}\text{d}x}=\{\begin{array}{l}\frac{l\pi }{2}\mathrm{,}k=\mathrm{0,}\\ \frac{l\pi }{4}\mathrm{,<mi>\; </mi>}k=2n\ne \mathrm{0,}\\ \mathrm{0,}k\ne 2n\ne 0.\end{array}$

$D_k={\displaystyle {\int }_0^{l\pi }\cos \frac{nx}{l}\cos \frac{nx}{l}\text{d}x}=\{\begin{array}{l}l\pi \mathrm{,}k=\mathrm{0,}\\ \frac{l\pi }{2}\mathrm{,}k\ne 0.\end{array}$

We drive the normal form of Hopf bifurcation for model (3) reduce on the center manifold as follows:

$\dot{G}=M\mu G+X{G}^2\overline{G}\mathrm{,}$ (23)

where $M$ and $X$ are defined by (20) and (22), respectively. Substituting $G=r{\text{e}}^{\text{i}\theta }$ into the Eq. (23), and we derive the amplitude equation on the center manifold as follows:

$\{\begin{array}{l}\dot{r}=\text{Re}\left(M\right)\mu r+\text{Re}\left(X\right)r^3\mathrm{,}\\ \dot{\theta }=\text{Im}\left(M\right)\mu +\text{Im}\left(X\right)r^2\mathrm{,}\end{array}$

then we state the following theorem.

Theorem 4. If the condition $\frac{\text{Re}\left(M\right)\mu }{\text{Re}\left(X\right)}<0$ holds, the model (3) generates periodic solutions near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ . If $\text{Re}\left(M\right)\mu <0$ $\left(\text{Re}\left(M\right)\mu >0\right)$ , the bifurcating periodic solutions are unstable (stable). Additionally, the direction of bifurcating periodic solution is forward (backward) if $\mu >0$ ($\mu <0$ ).

5. Numerical Simulations

In this section, we first analyze the parameter ranges, then simulations are performed numerically.

5.1. Selection of Parameter Values

(1) Releasing rate $c$ and natural mortality $m$ of C. cunea.

In the process of biological control of H. cunea, the release rate of C. cunea depends on the density of H. cunea. According to Ref. [6], the ratio of the C. cunea to the H. cunea is 3:1 in areas with more pests. Thus, we first consider the case when $c=3$ , and we also consider what happens when $c$ takes different values. According to Ref. [23], the mortality rate varies across different life stages of insects and may range from 0.1 to 0.9. Consequently, we set the natural mortality rate of C. cunea as $m=0.3$ in our paper.

(2) Self-diffusion coefficients $d_1$ and $d_2$ .

Considering the distinct life stages, it's evident that the egg and pupa stages are relatively stationary. Consequently, the self-diffusion coefficient $d_1$ for H. cunea can be conservatively estimated at 0.0003 due to limited mobility during these stages. C. cunea lays eggs in the pupae of H. cunea, thus, it moves with the movement of the H. cunea pupae but slightly faster. Therefore, we set its diffusion coefficient to $d_2=0.0004$ .

(3) Net birth rate $r$ and physical control fatality rate $h$ of H. cunea.

First of all,

$r=r_1-r_2-r_3\mathrm{,}$

where $r_1$ is birth rate of H. cunea. $r_2$ is natural death rate of H. cunea, and $r_3$ is death rate of predatory. According to the life table in Ref. [24], it showed mortality records for H. cunea across all stages due to various causes, and the range of H. cunea birth rate is $r_1\in \left(\mathrm{0.824,0.955}\right)$ , thus, we select $r_1=0.9$ . The range of H. cunea natural mortality rate is $r_2\in \left(\mathrm{0.105,0.509}\right)$ , and we select $r_2=0.3$ . The range of death rate of predatory is $r_3\in \left(\mathrm{0.069,0.238}\right)$ , and we select $r_3=0.15$ . Therefore, $r=0.9-0.3-0.15=0.45$ . Moreover, we select $h=0.3<0.45$ according to the analysis of Remark 1.

(4) Parasitism rate $\alpha$ and survival rate $\delta$ of C. cunea.

In Ref. [25], the experiment of rearing C. cunea with Tussah (Silkworm) pupae was carried out by three-cut method and inoculating method, and the range of parasitism rate is $\alpha \in \left[\mathrm{0.49,0.95}\right]$ . The range of survival rate of H. cunea is $\delta \in \left[\mathrm{0.85,1}\right]$ . In our paper, we select $\alpha =0.78$ and $\delta =0.93$ , this is because selecting intermediate high values for $\alpha$ and $\delta$ can provide a more robust and conservative estimate in experimental design. This helps ensure the reliability of experimental results and offers higher feasibility and effectiveness in practical applications.

(5) Other dimensionless parameters $K$ and $\beta$ .

We may select $K=18$ and $\beta =0.65$ for the environmental capacity and proportionality coefficient of H. cunea.

To sum up the analysis, we select, $d_1=0.0003$ , $d_2=0.0004$ , $r=0.45$ , $K=18$ , $\beta =0.65$ , $\alpha =0.78$ , $h=0.3$ , $\delta =0.93$ , $m=0.25$ , $c=3$ .

5.2. Simulations

In this subsection, we will show the existence of various spatiotemporal patterns for model (2) and model (3).

5.2.1. Non-spatial Model (2)

The hypotheses (H₁) and (H₂) hold when the parameters $d_1=0.0003$ , $d_2=0.0004$ , $r=0.45$ , $K=18$ , $\beta =0.65$ , $\alpha =0.78$ , $h=0.3$ , $\delta =0.93$ , $m=0.25$ are fixed and $c\in \left(\mathrm{0,4}\right)$ . According to Theorem 1, the equilibrium $E_{\mathrm{<mo>*</mo>}}$ is locally asymptotically stable for $\tau =0$ and $n=0$ . Figure 3 shows the change in the value of $E_{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ when $c\in \left(\mathrm{0,4}\right)$ . Fig. 3 shows that as the release rate of C. cunea increases, and then the value of $u^{\mathrm{<mo>*</mo>}}$ gradually decreases.

Biological interpretation 1:

Figure 3 shows that the population of the H. cunea reaches highest value when $c=0$ . This is attributed to the lack of effective natural enemies, which led to the rapid growth of the population. Furthermore, it can be observed that the populations of H. cunea and C. cunea tend to stabilize when $c>3$ . This is because they can quickly find and parasitize most of the H. cunea pupae when the amount of C. cunea is large enough, thereby inhibiting the increase in the H. cunea population. At the moment, the continuous parasitic action of the C. cunea keeps the H. cunea population at a low and stable level, rather than continuing to decline to near extinction.

Figure 3. The change of the equilibrium $E_{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ as $c$ increase.

Case 1: $c=0$ :

Considering the case with biological control, that is, $c=3$ , the equilibrium $E_{\mathrm{<mo>*</mo>}}=\left(\mathrm{0.0235,0.2947}\right)$ . This set of parameters satisfies hypotheses (H₁) and (H₂) according to Theorem 1, and the positive equilibrium $E_{\mathrm{<mo>*</mo>}}$ of model (2) is locally asymptotically stable (see Figure 4). Figure 4(a) shows that the vector field distribution of equilibria for non-spatial model (2) when $c=3$ . It is evident that $E_0$ is a saddle point, and $E_{\mathrm{<mo>*</mo>}}$ is locally asymptotically stable. Figure 4(b) shows the locally asymptotic stability of the equilibrium $E_{\mathrm{<mo>*</mo>}}$ .

(a) (b)

Figure 4. The stability of equilibrium $E_{\mathrm{<mo>*</mo>}}$ of model (3) when $c=3$ .

Case 2: $c=0$ :

Next, we consider the case of biological control is absent, that is, $c=0$ , the equilibrium $E_{\mathrm{<mo>*</mo>}}=\left(\mathrm{0.5302,0.2697}\right)$ . This set of parameters makes $A_0-F=0$ , then model (2) may generate periodic solution near the positive equilibrium $E_{\mathrm{<mo>*</mo>}}$ . Figure 5(a) shows that the vector field distribution of equilibria for non-spatial model (2) when $c=0$ , and it can be seen that the $E_0$ is a saddle point and the model (2) generates a stable periodic solution near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ . Figure 5(b) shows periodic solution near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ .

(a) (b)

Figure 5. The vector field near the equilibria $E_0$ and $E_{\mathrm{<mo>*</mo>}}$ , and periodic solution near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ of model (2) when $c=0$ .

Biological interpretation 2:

(1) Figure 4 reflects the effectiveness of the strategy of controlling the H. cunea population through the artificial release of C. cunea. Initially, both the H. cunea and C. cunea populations show significant oscillations, then gradually stabilize. The H. cunea population decreases rapidly, while the C. cunea population increases and gradually stabilizes. This indicates that the artificial release of C. cunea can effectively control the H. cunea population.

(2) Figure 5 indicates that both H. cunea and C. cunea experience periodic natural eruptions when $c=0$ . This erupt periodically is the result of the complex interactions between the H. cunea and C. cunea, reflecting the population regulation mechanisms in the natural ecosystem.

5.2.2. Spatial Model (3)

Case 1: $c=3$ :

The set of parameters $d_1=0.0003$ , $d_2=0.0004$ , $r=0.45$ , $K=18$ , $\beta =0.65$ , $\alpha =0.78$ , $h=0.3$ , $\delta =0.93$ , $m=0.25$ , $c=3$ satisfies hypotheses (H₁) and (H₂), $S_1=\left\{\mathrm{1,2,3,}\cdots \mathrm{,32}\right\}$ , $S_2=\varnothing$ and $S_3=\varnothing$ , then according to Theorem 2 and Eq. (10), we get the critical value ${\tau }_c={\tau }_1^{\left(0\right)}=3.0578$ . Therefore, the positive equilibrium $E_{\mathrm{<mo>*</mo>}}$ is locally asymptotically stable for $\tau \in \left(\mathrm{0,3.0578}\right)$ and unstable for $\tau \in \left(\mathrm{3.0578,}+\infty \right)$ . Figure 6 shows the asymptotic stability of the equilibrium $E_{\mathrm{<mo>*</mo>}}$ when $\tau =1.5\in \left(\mathrm{0,3.0578}\right)$ .

Figure 6. For $c=3$ , when $\tau =1.5$ , $E_{\mathrm{<mo>*</mo>}}$ of model (3) is locally asymptotically stable for $c=3$ .

By calculating Eq. (20) and Eq. (22), we obtain

$M=0.2933-0.1031\text{i}\mathrm{,}X=-8.3874+11.7617\text{i}\mathrm{,}$

according to Theorem 4, the bifurcating periodic solution is stable and forward. The model (3) will generate inhomogeneous periodic solutions near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ of model (3) for $\tau =3.058>{\tau }_c=3.0578$ (see Figure 7).

Figure 7. For $c=3$ , the inhomogeneous stable and forward periodic solutions for model (3) near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ when $\tau =3.058>{\tau }_c=3.0578$ .

We also consider the dynamics for the model (3) without nonlocal competition. Figure 8 shows that model (3) only generates stable homogeneous periodic solutions near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ for $\tau =3.058>{\tau }_c=3.0578$ if there is no nonlocal competition.

Figure 8. For $c=3$ , the homogeneous stable and forward periodic solutions for model (3) without nonlocal competition near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ when $\tau =3.058>{\tau }_c=3.0578$ .

Case 2: $c=0$ :

Considering spatial model (3) without biological control, that is, $c=0$ . This set of parameters satisfies hypotheses (H₁) and (H₂). Similarly, we get the critical value ${\tau }_c={\tau }_1^{\left(0\right)}=0.0102$ , thus, the equilibrium $E_{\mathrm{<mo>*</mo>}}$ is locally asymptotically stable for $\tau \in \left(\mathrm{0,0.0102}\right)$ and becomes unstable for $\tau \in \left(\mathrm{0.0102,}+\infty \right)$ . By calculation, $S_1=\left\{\mathrm{1,2,3,}\cdots \mathrm{,27}\right\}$ , $S_2=\varnothing$ , $S_3=\varnothing$ and $M=0.0107+0.0771\text{i}$ , $X=-0.0730+0.1089\text{i}$ . According to Theorem 4, the bifurcating periodic solution is stable and forward, then model (3) with $c=0$ will generate inhomogeneous periodic solutions near $E_{\mathrm{<mo>*</mo>}}$ for $\tau =0.011>{\tau }_c=0.0102$ (see Figure 9).

Figure 9. For $c=0$ , the inhomogeneous stable and forward periodic solutions for model (3) near the equilibrium $E_{\mathrm{<mo>*</mo>}}$ when $\tau =0.011>{\tau }_c=0.0102$ .

Biological interpretation 3:

(1) Figure 6 indicates that releasing C. cunea within the time delay ${\tau }_c$ can effectively control the population of H. cunea. The H. cunea population is gradually controlled both temporally and spatially, and tends to stabilize.

(2) Figure 7 shows that releasing C. cunea after the critical value ${\tau }_c$ will lead to periodic eruptions and instability in the H. cunea population. This means the H. cunea population cannot be effectively controlled and continues to pose a threat to the environment.

(3) By comparing Figure 7 and Figure 8, we find that the model (3) with nonlocal competition produces richer dynamics, which means the population dynamics are more similar to the behavior of the real ecosystem. In this way, the dynamic changes of populations in nature can be better simulated and understood. The phenomenon in Figure 8 may be due to the fact that H. cunea and C. cunea gather at the boundary when they fly to the boundary. When the density of H. cunea decreases, C. cunea flies in the opposite direction. Eventually, it causes large periodic oscillations on both sides, with equal numbers flying in opposite directions.

(4) By comparing Figure 7 and Figure 9, we find that when $c=0$ , the stability interval of model (3) becomes smaller, and the amplitude of the inhomogeneous periodic solution becomes bigger. In summary, without biological control, the equilibrium $E_{\mathrm{<mo>*</mo>}}$ becomes unstable more rapidly and the model is prone to periodic eruptions. This highlights that the artificial release of C. cunea is a critical measure for maintaining ecological balance and effectively controlling the population of the H. cunea. It also reflects the inadequacy of natural populations regulation mechanisms.

Case 3: $c\in \left(\mathrm{0,10}\right)$ :

By calculation, the hypotheses (H₁) and (H₂) hold when $d_1=0.0003$ , $d_2=0.0004$ , $r=0.45$ , $K=18$ , $\beta =0.65$ , $\alpha =0.78$ , $h=0.3$ , $\delta =0.93$ , $m=0.25$ , and $c\in \left(\mathrm{0,10}\right)$ . When $c\le 3$ , we get the critical value ${\tau }_c$ at $n=1$ , and when $c>3$ , we get the critical value ${\tau }_c$ at $n=0$ . Figure 10 shows that ${\tau }_c$ increases gradually with the increases of $c$ .

Figure 10. The value of ${\tau }_c$ with $c$ gradually increasing.

Biological interpretation 4:

The nonlocal competition influences model (3) when $c\le 3$ but does not impact model (3) when $c>3$ . Therefore, nonlocal competition has little effect on the model when the population of C. cunea is large enough to cover every area of the forest, leaving no escape for H. cunea.

Suggestions:

(1) According to the Biological interpretation 1, the amount of releasing C. cunea should be controlled. Releasing too much C. cunea will put pressure on the environment and cause unnecessary economic losses.

(2) According to Biological interpretation 2 and 3, it is necessary to select the appropriate time and place for releasing C. cunea. Releasing C. cunea can effectively control the density of H. cunea before most pupae emerge. At the same time, the deployment of releasing C. cunea is optimized to ensure their migration to areas with high H. cunea concentrations, thereby improving parasitic efficiency on H. cunea pupae. The migration patterns of H. cunea are regularly monitored to ensure effective targeting and control of the H. cunea population in different spatial locations.

(3) According to Biological interpretation 4, the influence of nonlocal competition can be ignored in the process of analyzing the models when $c>3$ .

(4) Based on the above analysis, our findings demonstrate that the population of artificially released C. cunea only requires reaching threefold the population size of H. cunea to achieve rapid suppression of H. cunea growth. This optimal ratio suggests that excessive release of C. cunea is unnecessary in practical applications, thereby effectively conserving both labor resources and economic expenditures.

6. Conclusions

In this paper, we studied a reaction-diffusion model of H. cunea and C. cunea with both time delay and nonlocal competition for biological control. We analyzed the existence and stability of the equilibria and the existence of Hopf and Turing bifurcations near the equilibria. We used the multiple time scales method to derive the normal form of the Hopf bifurcation reduced on the center manifold. Finally, a suitable set of parameters was selected for conducting numerical simulations. Spatial and non-spatial models were equally important aspects of our research. In the non-spatial model, we compared and analyzed the phenomena resulting from the with and without biological control. In the spatial model, we compared models with and without nonlocal competition, as well as with and without biological control. We found that the nonlocal competition can induce stable inhomogeneous periodic solutions of Hopf bifurcating. In the model without biological control, the critical time delay of the model decreased and generated periodic solutions with large amplitudes.

Results indicated that nonlocal competition significantly altered the spatial distribution and oscillation patterns of populations within the ecosystem. This phenomenon emphasised the need for geographically and environmentally sensitive biological control strategies. Effective biological control measures, particularly through the timely release of parasitic natural enemies, provided sustainable management and protection of H. cunea populations and ecosystems. The theoretical framework of the model, encompassing bifurcation analysis and the stability of periodic solutions, offered valuable theoretical insights for analyzing other invasive species-natural enemy systems. Future studies could validate the applicability of similar systems by altering the functional response structure of the model or adjusting the dispersal mechanism.

Acknowledgements

This study was funded by Fundamental Research Funds for the Central Universities (Grant No. HFW230600022) and Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2024A001).

NOTES

*Corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ge, X., He, S., Zhu, C., Wang, T., Xu, Z. and Zong, S. (2018) Projecting the Current and Future Potential Global Distribution of Hyphantria cunea (Lepidoptera: Arctiidae) Using CLIMEX. Pest Management Science, 75, 160-169. https://doi.org/10.1002/ps.5083
[2]	Lu, H., Song, H. and Zhu, H. (2017) A Series of Population Models for Hyphantria cunea with Delay and Seasonality. Mathematical Biosciences, 292, 57-66. https://doi.org/10.1016/j.mbs.2017.07.010
[3]	Yang, Z., Wei, J. and Wang, X. (2006) Mass Rearing and Augmentative Releases of the Native Parasitoid Chouioia cunea for Biological Control of the Introduced Fall Webworm Hyphantria cunea in China. Biocontrol, 51, 401-418. https://doi.org/10.1007/s10526-006-9010-z
[4]	Moon, J., Won, S., Choub, V., Choi, S., Ajuna, H.B. and Ahn, Y.S. (2022) Biological Control of Fall Webworm Larva (Hyphantria cunea Drury) and Growth Promotion of Poplar Seedlings (Populus × canadensis Moench) with Bacillus Licheniformis PR2. Forest Ecology and Management, 525, Article ID: 120574. https://doi.org/10.1016/j.foreco.2022.120574
[5]	Yamanaka, T., Tatsuki, S. and Shimada, M. (2001) Flight Characteristics and Dispersal Patterns of Fall Webworm (Lepidoptera: Arctiidae) Males. Environmental Entomology, 30, 1150-1157. https://doi.org/10.1603/0046-225x-30.6.1150
[6]	Xin, B., Liu, P., Zhang, S., Yang, Z., Daane, K.M. and Zheng, Y. (2017) Research and Application of Chouioia cunea Yang (Hymenoptera: Eulophidae) in China. Biocontrol Science and Technology, 27, 301-310. https://doi.org/10.1080/09583157.2017.1285865
[7]	Wen, X., Fang, G., Chai, S., He, C., Sun, S., Zhao, G., et al. (2024) Can Ecological Niche Models Be Used to Accurately Predict the Distribution of Invasive Insects? A Case Study of Hyphantria cunea in China. Ecology and Evolution, 14, Article No. 11159. https://doi.org/10.1002/ece3.11159
[8]	Misra, A.K., Yadav, A. and Sajid, M. (2024) Modeling the Effects of Baculovirus to Control Insect Population in Agricultural Fields. Mathematics and Computers in Simulation, 218, 420-434. https://doi.org/10.1016/j.matcom.2023.12.002
[9]	Duan, D., Niu, B. and Wei, J. (2019) Hopf-Hopf Bifurcation and Chaotic Attractors in a Delayed Diffusive Predator-Prey Model with Fear Effect. Chaos, Solitons & Fractals, 123, 206-216. https://doi.org/10.1016/j.chaos.2019.04.012
[10]	Chen, S., Shi, J. and Wei, J. (2012) The Effect of Delay on a Diffusive Predator-Prey System with Holling Type-II Predator Functional Response. Communications on Pure and Applied Analysis, 12, 481-501. https://doi.org/10.3934/cpaa.2013.12.481
[11]	An, Q. and Jiang, W. (2018) Bifurcations and Spatiotemporal Patterns in a Ratio‐dependent Diffusive Holling‐tanner System with Time Delay. Mathematical Methods in the Applied Sciences, 42, 440-465. https://doi.org/10.1002/mma.5299
[12]	Sun, H., Hao, P. and Cao, J. (2024) Dynamics Analysis of a Diffusive Prey-Taxis System with Memory and Maturation Delays. Electronic Journal of Qualitative Theory of Differential Equations, 40, 1-20. https://doi.org/10.14232/ejqtde.2024.1.40
[13]	Britton, N.F. (1989) Aggregation and the Competitive Exclusion Principle. Journal of Theoretical Biology, 136, 57-66. https://doi.org/10.1016/s0022-5193(89)80189-4
[14]	Furter, J. and Grinfeld, M. (1989) Local vs. Non-Local Interactions in Population Dynamics. Journal of Mathematical Biology, 27, 65-80. https://doi.org/10.1007/bf00276081
[15]	Wu, S. and Song, Y. (2020) Spatiotemporal Dynamics of a Diffusive Predator-Prey Model with Nonlocal Effect and Delay. Communications in Nonlinear Science and Numerical Simulation, 89, Article ID: 105310. https://doi.org/10.1016/j.cnsns.2020.105310
[16]	Hou, Y. and Ding, Y. (2024) Bifurcation Analysis of Pine Wilt Disease Model with Both Memory-Based Diffusion and Nonlocal Effect. Nonlinear Dynamics, 112, 10723-10738. https://doi.org/10.1007/s11071-024-09598-5
[17]	Ma, Y. and Yang, R. (2024) Hopf-Hopf Bifurcation in a Predator-Prey Model with Nonlocal Competition and Refuge in Prey. Discrete and Continuous Dynamical Systems B, 29, 2582-2609. https://doi.org/10.3934/dcdsb.2023193
[18]	Hou, Y., Ding, Y. and Jiang, W. (2025) Spatial Pattern Formation in Pine Wilt Disease Model with Prey-Taxis and Nonlocal Competition. Discrete and Continuous Dynamical Systems B, 30, 3479-3504. https://doi.org/10.3934/dcdsb.2025030
[19]	Abrams, P.A. (2000) The Evolution of Predator-Prey Interactions: Theory and Evidence. Annual Review of Ecology and Systematics, 31, 79-105. https://doi.org/10.1146/annurev.ecolsys.31.1.79
[20]	Ding, Y., Liu, G. and Zheng, L. (2023) Equivalence of MTS and CMR Methods Associated with the Normal Form of Hopf Bifurcation for Delayed Reaction-Diffusion Equations. Communications in Nonlinear Science and Numerical Simulation, 117, Article ID: 106976. https://doi.org/10.1016/j.cnsns.2022.106976
[21]	Nayfeh, A.H. (2007) Order Reduction of Retarded Nonlinear Systems—The Method of Multiple Scales versus Center-Manifold Reduction. Nonlinear Dynamics, 51, 483-500. https://doi.org/10.1007/s11071-007-9237-y
[22]	Achouri, H. (2024) Normal Forms for Retarded Functional Differential Equations Associated with Zero-Double-Hopf Singularity with 1: 1 Resonance. Electronic Journal of Qualitative Theory of Differential Equations, 59, 1-23. https://doi.org/10.14232/ejqtde.2024.1.59
[23]	Zhao, Z. (2020) Concepts and Applications of Insect Demography. Pest Management Science, 47, 904-911.
[24]	Wei, J., Yang, Z. and Su, Z. (2003) Use of Life Table to Evaluate Control of the Fall Webworm by the Parasitoid Chouioia cunea (Hymenoptera: Eulophidae). Acta Entomologica Sinica, 46, 318-324.
[25]	Tian, X., Wang, H. and Jiang, F. (2002) Reproduction and Biological Characteristic of Chouioia cunea. Journal of Forestry Research, 13, 331-333. https://doi.org/10.1007/bf02860102