The dynamics of a delayed predator-prey model with square root functional response and stage structure

1.   Introduction

There are various biological species in nature, in which many interactions exist. Since the pioneering progress of the work of Lotka ^[1] and Volterra ^[2] in the 1920s on the interaction between predators and prey, the dynamics of predator-prey models have been widely studied by scholars ^[3,4,5,6].

The functional responses are the functions of the density of prey due to the predation rate. In the famous Lokta-Volterra model, the number of prey captured by predators per unit of time increases with the increase of the prey quantity, without considering the impact of digestion saturation on the system, which has some shortcomings. In 1965, Holling ^[7] developed three types of functional responses based on different species. The Holling-type functional responses can reflect the phenomenon of predator density gradually saturating as the number of prey increases, which sparked extensive research ^[8,9,10]. With further research, scholars have put forward some functional responses that are more in line with reality ^{[11,12,13,14]}. When exploring the herd behavior between populations, Ajraldi et al. ^[15] found that cattle herds gather together for defensive purposes. Usually, weaker individuals occupy the interior of the herd, while stronger individuals are outside. Due to the discovery that the number of cattle is directly proportional to the square root of density, they considered that using the square root of the prey density to simulate the functional response of prey exhibiting this herd behavior was appropriate. Braza ^[16] considered a predator-prey model with square root functional response, in which prey exhibit strong herd structure. Salman et al. ^[17] presented a discrete predator-prey model with square root functional response and studied the stability and bifurcation of the system. Suleman et al. ^[18] structured a slow-fast predator-prey model with herd behavior, where the presence and stability of fixed points were investigated. He and Li ^[19] introduced the square root response function into a Leslie-Gower predator-prey model and discussed the impact of square root functional response on the dynamic behaviors of the system. Furthermore, researchers also studied other systems with square root response functions ^{[20,21,22,23,24,25]}.

In ^[22], Mortuja et al. proposed a predator-prey model with nonlinear prey harvesting and square root functional response:

where $X\left(t \right)$ and $Y\left(t \right)$ represent the densities of the prey and predator population at the time $t$. $r$ denotes the growth rate of the prey and $p$ is the environmental carrying capacity. $\alpha$ represents the search efficiency of the predator for prey and $\beta$ represents the natural mortality rate of the predator in the absence of prey. ${t_h}$ is the average processing time of the predator for the prey, and $m$ is the consumption rate. In the nonlinear harvesting section, $q$ is the capture capability coefficient, and $E$ is the harvest effort. ${m_1}$ and ${m_2}$ are appropriate constants. They explored the stability of the system and discussed Hopf bifurcation and saddle-node bifurcation.

According to the characteristics and growth laws of organisms, young biological individuals cannot compete and reproduce for a long time. Therefore, it is more practical to divide the population into different stage structures. The predator-prey models with stage structure have been paid much attention by scholars. The predators or prey are divided into immature and mature stages in these models. Dai et al. ^[26] considered a Holling type-III predator-prey model with stage structure for prey and researched the asymptotic stability of the positive equilibrium point. Meng et al. ^[27] obtained global Hopf bifurcation in a predator-prey model with a stage structure for prey. A predator-prey model with a time delay and stage structure for predators was analyzed by Zhou ^[28]. He found that the combination of the time delay and stage structure can affect the dynamic behavior of the system. Zhao and Zeng ^[29] dealt with a predator-prey model with stage structure for predator and ratio-dependent functional response. They extended the existence of stationary distribution and obtained sufficient conditions for the extinction of predator populations.

In addition, interactions between populations inevitably cause time delays. For example, the reproductive cycle of a population, the time it takes for predators to absorb and convert their prey into their energy, and so on, are phenomena that exist in a time delay. From a dynamic perspective, differential equations with a time delay typically exhibit more complex dynamic behavior than ordinary differential equations. Introducing a time delay into predator-prey models has practical significance. Zhang et al. ^[30] came up with a predator-prey model with a discrete time delay, and the results showed that predators would become extinct with increasing the time delay. Yang and Jin ^[31] proposed and analyzed a diffusive model with a time delay. They studied the stability of the system and Hopf bifurcation. Some scholars have also studied predator-prey models with time delays ^{[32,33,34,35,36]}.

Based on the above discussions and inspiration from the work of Mortuja et al. ^[22], in this paper, we consider introducing a stage structure into the system (1.1) and dividing the predator into two stages: an immature predator and a mature predator, and the prey is only preyed upon by mature predators. In addition, the time delay of the predator is also considered. The predator-prey model with stage structure and a time delay is as follows:

where $X\left(t \right)$ represents the population density of the prey, and ${Y_1}\left(t \right)$ and ${Y_2}\left(t \right)$ denote the population density of immature and mature predators, respectively. $D$ is the rate at which immature predators become mature predators. $\beta$ and $d$ respectively express the natural mortality rates of immature and mature predators. $\tau$ is taken as a gestation delay of the mature predator. The other parameters $r$, $p$, $\alpha$, ${t_h}$, $m$, $q$, $E$, ${m_1}$, and ${m_2}$ represent the same meaning as in system (1.1).

For convenience, we have scaled the variables and the parameters of the system (1.2). The variables are scaled as

and the other parameters are made dimensionless as follows:

In addition, to avoid heavy notation, we still use $t$, $\tau$, $D$ to represent $\bar t$, $\bar \tau$, $\bar D$. Through the above transformation, the dimensionless form of the system (1.2) is as follows:

with initial conditions

where $R_{ + 0}^3 = \left\{ {\left({{x_1}, {x_2}, {x_3}} \right):{x_i} \geqslant 0, i = 1, 2, 3} \right\}$.

According to the above definition, all parameters here are positive.

The rest of this article is organized as follows: In Section 2, we discuss positivity, boundedness, the stability of the equilibrium points, and the existence of Hopf bifurcation. Section 3 is devoted to the nature of Hopf bifurcation. In Section 4, some numerical simulations are conducted to verify our analytical results. The article ends with a simple summary.

2.   Local stability and Hopf bifurcation analysis

In this section, we will study the positivity and boundedness of the solution, the local stability of the equilibria of the system (1.3), and prove the existence of Hopf bifurcation.

2.1.   Positivity and boundedness of the solution

2.1.1.   Positivity of the solution

Considering the biological background, positivity means the survival of the population, so it is necessary to prove the positivity of the solution. The system (1.3) can be rewritten into the following matrix form:

where $X = {\left({x\left(t \right), {y_1}\left(t \right), {y_2}\left(t \right)} \right)^T} \in {R^3}$, and $J\left(X \right)$ can be written as:

Because $J\left(X \right)$ and $\frac{{\partial J}}{{\partial X}}$ are continuous on $R_ + ^3$, then $J:R_ + ^3 \to {R^3}$ is locally Lipschitz. Through the standard theory of the ordinary differential equation system, the unique solution of the system (1.3) under any initial condition $J\left(0 \right) = {J_0} = \left({x\left(0 \right), {y_1}\left(0 \right), {y_2}\left(0 \right)} \right) \in R_ + ^3$ can be obtained.

Then, the system (1.3) can be rewritten as:

where

Taking the second equation as an example, we can obtain

that is

Integrating both sides of the above equation from $0$ to $t$, we obtain

then

and we have

Therefore, we have

Similarly, it can be concluded that $x\left(t \right) \geqslant 0$, ${y_2}\left(t \right) \geqslant 0$. Therefore, for any $t > 0$, the solution $X\left(t \right) = \left({x\left(t \right), {y_1}\left(t \right), {y_2}\left(t \right)} \right)$ satisfying the initial condition $X\left(0 \right) = \left({x\left(0 \right), {y_1}\left(0 \right), {y_2}\left(0 \right)} \right) \in R_ + ^3$ remains nonnegative in the region $R_ + ^3$.

Next, we show that all solutions of system (1.3) under the initial condition are bounded.

2.1.2.   Boundedness of the solution

Boundedness is indicated as a result of limited resources. We know that all solutions of system (1.3) with the initial condition are positive. From the first equation of system (1.3), we obtain $\dot x\left(t \right) \leqslant x\left(t \right)\left({1 - x\left(t \right)} \right)$. Considering that $\dot s\left(t \right) = s\left(t \right)\left({1 - s\left(t \right)} \right)$, the solution that satisfies the critical condition is $s\left(t \right) = \frac{{s\left(0 \right) - 1}}{{s\left(0 \right) - 1 + {e^{ - t}}}}$. By comparing $x\left(t \right)$ and $s\left(t \right)$, it is obtained that $x\left(t \right) \leqslant s\left(t \right)$ for all $t \geqslant 0$. Then, all solutions of system (1.3) satisfy $0 < x\left(t \right) \leqslant 1$ for $\forall t \geqslant 0$.

Next, we define that $V\left(t \right) = x\left({t - \tau } \right) + \frac{{{y_1}\left(t \right)}}{b} + \frac{{{y_2}\left(t \right)}}{b}$. The derivative of $V\left(t \right)$ along the solution of system (1.3) is

and we choose $\mu = \min \left\{ {{d_1}, {d_2}} \right\}$ as a constant. Then,

Let $\frac{{{{\left({1 + \mu } \right)}^2}}}{4} = M$, and we obtain that

Then, for a sufficiently large $t$, there is $0 < V\left(t \right) \leqslant \frac{M}{\mu }$. Therefore, the solutions $x\left(t \right)$, ${y_1}\left(t \right)$, and ${y_2}\left(t \right)$ under the initial conditions are bounded.

2.2.   Existence of equilibrium points

The system (1.3) has nonnegative equilibrium points.

(i) System (1.3) always has a trivial equilibrium point ${E_0}\left({0, 0, 0} \right)$.

(ii) (a) If (H1) $1 - s > 0$, (H2) ${\left({1 - s} \right)^2} - 4\left({h - s} \right) > 0$, and (H3) $h > s$ hold, system (1.3) has two predator-free equilibrium points ${E_i}\left({{x_i}, 0, 0} \right)\;\left({i = 1, 2} \right)$, where

(b) If (H1), (H2), and (H4) $h < s$ hold, only a boundary equilibrium ${E_1}\left({{x_1}, 0, 0} \right)$ exists.

(c) When (H1) and (H5) ${\left({1 - s} \right)^2} = 4\left({h - s} \right)$ hold, system (1.3) has a predator-free equilibrium point ${E_3}\left({{x_3}, 0, 0} \right)$, where ${x_3} = \frac{{1 - s}}{2}$.

(iii) The unique positive equilibrium point ${E^*}\left({x', {{y'}_1}, {{y'}_2}} \right)$ exists in the system (1.3) if condition (H6) ${c_0} > {c_1} > 0$, where ${c_0} = 1 - x'$, ${c_1} = \frac{h}{{x' + s}} > 0$, and (H7) $bD - a{d_2}({d_1} + D) > {d_2}\left({{d_1} + D} \right)$ are met, where $x' = {\left[{\frac{{{d_2}\left({{d_1} + D} \right)}}{{bD - a{d_2}\left({{d_1} + D} \right)}}} \right]^2}$, ${y'_1} = \frac{{b\left({c{}_0 - {c_1}} \right)}}{{{d_1} + D}}x'$, ${y'_2} = \frac{D}{{{d_2}}}{y'_1}$.

2.3.   Stability of equilibrium points

In this part, we mainly discuss the stability of the trivial equilibrium point ${E_0}$ and the positive equilibrium point ${E^*}$.

2.3.1.   Trivial equilibrium point ${E_0}\left({0, 0, 0} \right)$

Theorem 2.1. If (H4) $h < s$ holds, the equilibrium point ${E_0}\left({0, 0, 0} \right)$ is unstable. The equilibrium point ${E_0}\left({0, 0, 0} \right)$ is locally asymptotically stable, when (H3) $h > s$ is satisfied.

Proof. The characteristic equation at the equilibrium point ${E_0}$ is

All of the three characteristic roots of the equation can be obtained: ${\lambda _1} = - {d_2} < 0$, ${\lambda _2} = - \left({{d_1} + D} \right) < 0$, ${\lambda _3} = 1 - \frac{h}{s}$. If condition (H4) $h < s$ holds, the eigenvalue ${\lambda _3}$ is a positive real root, and therefore, ${E_0}$ is unstable. If condition (H3) $h > s$ holds, ${\lambda _3}$ is a negative real root. All eigenvalues are negative, which leads to ${E_0}$ being locally asymptotically stable. The proof is complete.

2.3.2.   Positive equilibrium point ${E^*}\left({x', {{y'}_1}, {{y'}_2}} \right)$

Next, we analyze the stability of the unique positive equilibrium point ${E^*}$ of the system (1.3).

Let $\bar x\left(t \right) = x\left(t \right) - x'\left(t \right)$, ${\bar y_1}\left(t \right) = {y_1}\left(t \right) - {y'_1}\left(t \right)$, and ${\bar y_2}\left(t \right) = {y_2}\left(t \right) - {y'_2}\left(t \right)$. For convenience, we still use $x\left(t \right)$, ${y_1}\left(t \right)$, ${y_2}\left(t \right)$ to represent $\bar x\left(t \right)$, ${\bar y_1}\left(t \right)$, ${\bar y_2}\left(t \right)$. We get

where the coefficients and functions of the system (2.2) are listed in Appendix A.

Then, the system (1.3) is linearized to

Therefore, the characteristic equation of system (2.3) is written as

where

Next, we will discuss two different situations.

When $\tau = 0$, Eq (2.4) becomes

where ${l_{20}} = {l_{10}} + {h_{11}}$, ${l_{21}} = {l_{11}} + {h_{10}}$.

Therefore, by the Routh-Hurwitz discriminant method, we obtain that all of the roots of Eq (2.5) have the negative real part, if (H8) ${l_{12}} > 0$, ${l_{20}} > 0$, ${l_{12}}{l_{21}} - {l_{20}} > 0$ holds. This means that the positive equilibrium point ${E^*}$ is locally asymptotically stable.

When $\tau \ne 0$, let $i\omega \left({\omega > 0} \right)$ be a pure imaginary root of Eq (2.4). By separating the real and imaginary parts, it can be obtained that

By squaring and adding, we have

where

Let ${\omega ^2} = \upsilon$, and Eq (2.7) can be written as

Let

We can get $g\left(0 \right) = {e_{10}}$, $\mathop {\lim }\limits_{\upsilon \to + \infty } g\left(\upsilon \right) = + \infty$, and the derivative of Eq (2.9) is $g'\left(\upsilon \right) = 3{\upsilon ^2} +$$2{e_{12}}\upsilon + {e_{11}}$.

Therefore, the roots of Eq (2.8) can be studied by using the method of reference ^[37], along with the following lemma.

Lemma 2.1. For Eq (2.8),

(i) if (H9) ${e_{10}} \geqslant 0$, $\Delta = e_{12}^2 - 3{e_{11}} \leqslant 0$ holds, then Eq (2.8) does not have positive roots.

(ii) if (H10) ${e_{10}} \geqslant 0$, $\Delta = e_{12}^2 - 3{e_{11}} > 0$, $\upsilon = \frac{{ - {e_{11}} + \sqrt \Delta }}{3} > 0$, $g\left(\upsilon \right) \leqslant 0$ or (H11) ${e_{10}} < 0$ holds, then Eq (2.8) has positive roots.

Therefore, we suppose that Eq (2.8) has positive roots. Without loss of generality, we assume that it has three positive roots, denoted by ${\upsilon _1}$, ${\upsilon _2}$ and ${\upsilon _3}$, respectively. Then Eq (2.7) also has three positive roots, defined as ${\omega _k} = \sqrt {{\upsilon _k}}$, $k = 1, 2, 3$. Furthermore, the critical value corresponding to the time delay $\tau _k^{\left(j \right)}$ can be expressed as

Let $\lambda \left(\tau \right) = \varsigma \left(\tau \right) + i\omega \left(\tau \right)$ be the root of Eq (2.4) and $\varsigma \left({\tau _k^j} \right) = 0$, thus $\pm {\omega _k}$ is a pair of purely imaginary roots of Eq (2.4) with $\tau = \tau _k^{\left(j \right)}$. We define $\tau ' = \mathop {\min }\limits_{k = 1, 2, 3;j = 0, 1, 2, \cdots } \left\{ {\tau _k^{\left(j \right)}} \right\}$ and $\omega ' = {\omega _k}$.

Lemma 2.2. If (H12) $g'\left({{{\omega '}^2}} \right) \ne 0$ holds, then ${\left\{ {\frac{{d\left({\text{Re} \lambda } \right)}}{{d\tau '}}} \right\}_{\lambda = i\omega '}} \ne 0$.

Proof. Taking the derivative of Eq (2.4) on both sides with respect to $\tau$, we obtain

Further calculation reveals that

According to Eq (2.6), it can be calculated that

Since ${\left\{ {\frac{{d\left({\text{Re} \lambda } \right)}}{{d\tau '}}} \right\}_{\lambda = i\omega '}}$ and ${\left\{ {\text{Re} {{\left({\frac{{d\lambda }}{{d\tau '}}} \right)}^{ - 1}}} \right\}_{\lambda = i\omega '}}$ have the same symbol, then

So, if (H12) $g'\left({{{\omega '}^2}} \right) \ne 0$ is met, ${\left\{ {\frac{{d\left({\text{Re} \lambda } \right)}}{{d\tau '}}} \right\}_{\lambda = i\omega '}} \ne 0$ can be obtained.

The proof of the transversality condition is complete.

In summary, the following conclusions can be drawn.

Theorem 2.2. For system (1.3), suppose (H6), (H7), and (H8) hold, then the following results satisfy:

(i) If (H9) holds, the positive equilibrium point ${E^*}$ is locally asymptotically stable when $\tau ' \geqslant 0$.

(ii) If (H10) or (H11) and (H12) are satisfied, the positive equilibrium point ${E^*}$ is locally asymptotically stable when $\tau \in \left[{0, \tau '} \right)$ and unstable when $\tau > \tau '$. In addition, when $\tau = \tau '$, a Hopf bifurcation occurs at the positive equilibrium point ${E^*}$.

3.   Property analysis of Hopf bifurcation

In this section, we will apply the central manifold theorem and the normal form method in ^[38] to analyze the properties of Hopf bifurcation at the positive equilibrium point ${E^*}$.

Let $\tau = \tau ' + \xi$, $\xi \in R$, $t = s\tau$, $x\left({s\tau } \right) = \tilde x\left(s \right)$, ${y_1}\left({s\tau } \right) = {\tilde y_1}\left(s \right)$, and ${y_2}\left({s\tau } \right) = {\tilde y_2}\left(s \right)$, where $x = \tilde x$, ${y_1} = {\tilde y_1}$, ${y_2} = {\tilde y_2}$, and $t = s$. Then the system (1.3) can be described as the following functional differential equation in $C = C\left({\left[{- 1, 0} \right], {R^3}} \right)$

where $u\left(t \right) = {\left({x\left(t \right), {y_1}\left(t \right), {y_2}\left(t \right)} \right)^T} \in C$, ${u_t}\left(\theta \right) = u\left({t + \theta } \right) = {\left({x\left({t + \theta } \right), {y_1}\left({t + \theta } \right), {y_2}\left({t + \theta } \right)} \right)^T} \in C$, and ${L_\xi }:C \to {R^3}$, $F:R \times C \to {R^3}$ are written as

and

where the specific form of ${L_\xi }\left(\phi \right)$ and $F\left({\xi, \phi } \right)$ are given in Appendix B.

Therefore, according to the Riesz representation theorem, there exists a bounded variation function $\eta \left({\theta, \xi } \right)$ of $3 \times 3$ where $\theta \in \left[{- 1, 0} \right]$, such that

${L_\xi }\phi = \int_{ - 1}^0 {d\eta \left({\theta, \xi } \right)} \phi \left(\theta \right)$, for $\phi \in C$. (3.4)

We can choose

where $\delta \left(\theta \right) = \left\{ {\begin{array}{*{20}{c}} {0, \; \; \theta \ne 0, } \\ {1, \; \; \theta = 0.} \end{array}} \right.$

For $\phi \in C\left({\left[{- 1, 0} \right], {R^3}} \right)$, we define

and

Then, Eq (3.1) can be rewritten in the form of the following operator equation

For $\varphi \in C\left({\left[{- 1, 0} \right], {{\left({{R^3}} \right)}^*}} \right)$, where ${\left({{R^3}} \right)^*}$ is a three-dimensional row vector space, define the adjoint operator ${A^*}$ of $A\left(0 \right)$

For $\phi \in C\left({\left[{- 1, 0} \right], {R^3}} \right)$ and $\varphi \in C\left({\left[{- 1, 0} \right], {{\left({{R^3}} \right)}^*}} \right)$, we define bilinear inner product maps

where $\eta \left(\theta \right) = \eta \left({\theta, 0} \right)$, $A = A\left(0 \right)$, and ${A^*}$ are adjoint operators.

According to the analysis in the previous section, there is a pair of pure imaginary roots $\pm i\omega '\tau '$ for $A\left(0 \right)$ and ${A^*}$. We assume that $q\left(\theta \right) = {\left({1, {q_1}, {q_2}} \right)^T}{e^{i\omega '\tau '\theta }}$ is the eigenvector of the feature root $A\left(0 \right)$ corresponding to $i\omega '\tau '$, and ${q^*}\left(s \right) = \Omega \left({1, q_1^*, q_2^*} \right){e^{i\omega '\tau 's}}$ is the eigenvector of the feature root ${A^*}$ corresponding to $- i\omega '\tau '$. After some calculations, it can be concluded that

Thus, combined with Eq (3.10), there is

Hence, it is obtained that

so that $\left\langle {{q^*}\left(s \right), q\left(\theta \right)} \right\rangle = 1$, $\left\langle {{q^*}\left(s \right), \bar q\left(\theta \right)} \right\rangle = 0$.

Next, adopting the approach in ^[38] and the calculation process similar to ^[37], we obtain the coefficients that can be applied to determine the nature of the Hopf bifurcation.

The details of the derivations of coefficients ${g_{20}}$, ${g_{11}}$, ${g_{02}}$ and ${g_{21}}$ are given in Appendix C.

Then, the following formulas can be calculated to obtain:

Through the above analysis, the following conclusions can be drawn.

Theorem 3.1. For system (1.3),

(i) the sign of ${\rho _2}$ can determine the direction of the Hopf bifurcation, if ${\rho _2} > 0\left({{\rho _2} < 0} \right)$, then the Hopf bifurcation is supercritical (subcritical).

(ii) the sign of ${\beta _2}$ can determine the stability of the bifurcation periodic solutions, if ${\beta _2} < 0\left({{\beta _2} > 0} \right)$, then the bifurcation periodic solutions are stable (unstable).

(iii) the sign of ${T_2}$ determines the periodic change of the bifurcating periodic solutions, if ${T_2} > 0\left({{T_2} < 0} \right)$, it represents an increase (decrease) in the period of the bifurcating periodic solutions.

4.   Numerical simulations

In this section, we use MATLAB software to perform some numerical simulations to demonstrate the results of the analysis in the above sections.

In the system (1.3), we select a set of parameters $a = 0.4$, $b = 0.4$, $h = 0.2$, $s = 0.4$, ${d_1} = 0.1$, ${d_2} = 0.1$, and $D = 0.1$. We found the trivial equilibrium point ${E_0}\left({0, 0, 0} \right)$ is always unstable.

When the conditions (H6) and (H7) hold, system (1.3) has a unique positive equilibrium ${E^*}\left({0.3906, 0.2784, 0.2784} \right)$. For $\tau = 0$, (H8) is also met under this set of parameters. Then the positive equilibrium ${E^*}$ is locally asymptotically stable, which is presented in Figure 1.

For $\tau \ne 0$, we have $\omega ' \approx 0.0710$ and $\tau ' \approx 8.9781$, and the transversality condition is satisfied. According to Theorem 2.2, when $\tau = 5.75 < \tau '$, the positive equilibrium ${E^*}$ is locally asymptotically stable, and two populations will be stable eventually, which is shown in Figure 2. For $\tau = 9.73 > \tau '$, ${E^*}$ becomes unstable. Thus, the Hopf bifurcation occurs when $\tau ' \approx 8.9781$ and a family of periodic solutions bifurcate from the positive equilibrium ${E^*}$. From Figure 3, it can be seen that the two populations are unstable. It can further prove the instability of the system caused by a time delay, inducing species to fluctuate.

Furthermore, after the calculation of Eq (3.13), we can gain ${c_1}\left(0 \right) = - 43.2389 - 125.4193i$, ${\beta _2} = - 86.4778$, ${\rho _2} = 7.6070 \times {10^4}$, and ${T_2} = 428.9949$. From Theorem 3.1, it can be obtained that the Hopf bifurcation is supercritical, the bifurcating periodic solutions are stable, and its period increases. From a biological perspective, the stable bifurcating periodic solutions indicate the coexistence of prey and predator in an oscillatory mode.

5.   Conclusions

A delayed predator-prey model with square root functional response and stage structure was studied in this paper inspired by the work of Mortuja et al. ^[22], where predators were divided into two stages: immature predators and mature predators, and $\tau$ represented the gestation delay of the mature predator in prey and predator interaction.

First, we analyzed the positivity and boundedness of the solution, and studied the existence of three types of different equilibrium points, as well as the relevant stability conditions for the trivial equilibrium and positive equilibrium of the system (1.3). We found that the equilibrium ${E_0}$ was always unstable, and if (H8) ${l_{12}} > 0$, ${l_{20}} > 0$, ${l_{12}}{l_{21}} - {l_{20}} > 0$ were satisfied, then the positive equilibrium ${E^*}$ was locally asymptotically stable when $\tau = 0$(see Figure 1). On this basis, the sufficient conditions for the Hopf bifurcation of the system (1.3) at the positive equilibrium point ${E^*}$ when $\tau > 0$ were also studied, and the critical value of the time delay was calculated. When the time delay was less than the critical value $\tau '$, the positive equilibrium ${E^*}$ was locally asymptotically stable (see Figure 2). When the time delay increased and passed the critical value $\tau '$, the stability of positive equilibrium ${E^*}$ switched from locally asymptotically stable to unstable and generated a Hopf bifurcation, as well as a family of periodic solutions bifurcated from the positive equilibrium ${E^*}$ (see Figure 3). The research results indicate that the time delay $\tau$ plays an important role in the system (1.3) which can affect the stability of this system and the existence of the Hopf bifurcation. In addition, by using the central manifold theorem and normal form theory, explicit formulas were deduced to determine the direction of the Hopf bifurcation and the properties of the periodic solutions. In the end, numerical simulations validated our theoretical analysis. Mortuja et al. ^[22] observed that the harvesting rate can influence population coexistence and ecological balance. Compared with reference ^[22], we introduced time delay and stage structure into the system and have shown that gestation delay is a very important factor to understand the dynamics of the system. The stable bifurcating periodic solutions mean that the prey, immature predator, and mature predator can oscillate periodically to survive together. This is valuable from the view of ecology.

This paper investigated a predator-prey model with a stage structure for predators. What kind of dynamic behavior would occur if the prey had a stage structure or if both populations had a stage structure? We will leave this for future research.

Use of AI tools declaration

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

Acknowledgments

The authors would like to thank the editors and referees for their careful reading of the manuscript and valuable suggestions. This work was supported by the National Natural Science Foundation of China (Grant Nos. 12102148 and 12372012).

Conflict of interest

The authors declare that there are no conflicts of interest.

Appendix

Appendix A. The coefficients and functions of the system (2.2)

Appendix B. The specific form of ${L_\xi }\left(\phi \right)$ and $F\left({\xi, \phi } \right)$

${L_\xi }:C \to {R^3}$, $F:R \times C \to {R^3}$ are written as

and

where

Appendix C. The details of the derivations of coefficients ${g_{20}}$, ${g_{11}}$, ${g_{02}}$, and ${g_{21}}$

We need to calculate the coordinate representation of the center manifold ${C_0}$ at the origin. We define

on the central manifold ${C_0}$, and there is

We obtain

where

with

From Eq (C1), we get

therefore, it is obtained that

From the above formula, it can be calculated that

By comparing the coefficients of Eqs (C4) and (C5), we get

Through further calculation, we obtain the unknown coefficients ${W_{20}}\left(\theta \right)$ and ${W_{11}}\left(\theta \right)$ of ${g_{21}}$

where ${M_1} = {\left({M_1^{\left(1 \right)}, M_1^{\left(2 \right)}, M_1^{\left(3 \right)}} \right)^T} \in {R^3}$ and ${M_2} = {\left({M_2^{\left(1 \right)}, M_2^{\left(2 \right)}, M_2^{\left(3 \right)}} \right)^T} \in {R^3}$ are constant vectors, and given

with