The predator-prey systems play an important and fundamental role among the relationships between the biological population. There exists a strong interest during the last few decades to study predator-prey system in population dynamics by mathematicians and biologists. It has turned out that many researches in mathematical models contribute to understanding the interactions of predator and prey in these systems. Then various dynamics behaviors are achieved, such as steady states, oscillations, bifurcations and chaos depending on the model parameters, see [1-5].

Recently, the models governed by the ordinary and partial differential equations involving delays attracted people's attention and were successfully applied to many problems in many research areas (such as biology, population dynamics, neural networks, feedback controlled mechanical systems and so on, see [6-7]). In delayed differential equations, period solutions can arise through the Hopf bifurcation. Several methods for analyzing the nature of Hopf bifurcation had been proposed by many authors. Thereinto, the center manifold theory is one of the rigorous mathematical tools to study bifurcations of delay differential equations, see [13-20].

From the view of biology, many species have a life history that take them through an immature and mature stage, and the predator in the first stage often has no ability to attack prey and reproduce. Many author have investigated the models with stage-structure and obtained varies dynamics behaviors, see [4-12]. In [4], a predator-prey model with stage-structure and Holling type Ⅱ response is considered as follow:

$\left\{ \begin{array}{l}\dot{x}(t)=x(t)\left(r-ax(t)-\frac{a_{1}y_{2}(t)}{1+mx(t)}\right), \\ \dot{y}_{1}(t)=\frac{a_{2}x(t)y_{2}(t)}{1+mx(t)}-r_{1}y_{1}(t)-Dy_{1}(t), \\ \dot{y}_{2}(t)=Dy_{1}(t)-r_{2}y_{2}(t), \end{array} \right.$

(1.1)

where $x(t), y_{1}(t), y_{2}(t)$ represent the densities of prey, immature and mature predator, respectively, $ax(t), r_{1}, r_{2}$ represent the death rate of prey, immature and mature predator, respectively. $D$ is the convert rate from immature predator to mature predator. $\frac{ax}{1+mx}$ is response function of the mature predators' capture. The stability of model (1.1) was studied by Xu et al. [4] in detail.

Consider the time-delay due to the pregnant of the mature predator for advanced animals influences the stability of the equilibrium in delayed differential equations, we establish a delayed Holling Type Ⅲ response predator-prey system as follows:

$\left\{ \begin{array}{l}\dot{x}(t)=x(t)\left(r-ax(t)-\frac{a_{1}x(t)y_{2}(t)}{1+mx^{2}(t)}\right), \\ \dot{y}_{1}(t)=\frac{a_{2}x^{2}(t-\tau)y_{2}(t-\tau)}{1+mx^{2}(t-\tau)}-r_{1}y_{1}(t)-Dy_{1}(t), \\ \dot{y}_{2}(t)=Dy_{1}(t)-r_{2}y_{2}(t). \end{array} \right.$

(1.2)

The rest of this paper is organized as follows: in the next section, by analyzing the corresponding characteristic equations, the local stability of each of the feasible equilibrium of system (1.2) is discussed and the existence of a Hopf bifurcation at the coexistence equilibrium is established. The stability and the direction of periodic solutions bifurcating from Hopf bifurcations are investigated by using the normal form theory [23] and the center manifold theorem [24] in Section 3. The numerical simulations are carried out to support the theoretical analysis of the research in Section 4. Finally, a brief discussion is given in Section 5.

In this section, we analysis the local stability of each of feasible equilibria of system (1.2) and the existence of Hopf bifurcations at the coexistence equilibrium. Obviously, system (1.2) always has a trivial equilibrium $E_{0}(0, 0, 0)$ and a predator-extinction equilibrium $E_{1}(\frac{r}{a}, 0, 0)$. If the following condition holds:

$H_{1}:a_{2}hr^{2}>r_{2}(a^{2}+mr^{2})(D+r_{1}), $

then system (1.2) has a unique coexistence equilibrium $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$, where

$x^{*}=\sqrt{\frac{r_{2}(D+r_{1})}{a_{2}D-mr_{2}(D+r_{1})}}, \;\; y_{1}^{*}=\frac{r_{2}}{D}y^{*}_{2}, \;\; y^{*}_{2}=\frac{(r-ax^{*})(1+m(x^{*})^{2})}{a_{1}x^{*}}.$

(2.1)

The characteristic equation of system (1.2) at the trivial equilibrium $E_{0}(0, 0, 0)$ is of the form

$(\lambda-r)(\lambda+D+r_{1})(\lambda+r_{2})=0.$

(2.2)

It is apparently to see that the trivial equilibrium $E_{0}(0, 0, 0)$ is always unstable.

The characteristic equation of system (1.2) at the trivial equilibrium $E_{1}(\frac{r}{a}, 0, 0)$ takes the form

$(\lambda+r)(\lambda^{2}+(r_{2}+D+r_{1})\lambda+r_{2}(D+r_{1})-\frac{a_{2}D(x_{0})^{2}}{1+m(x_{0})^{2}}e^{-\lambda\tau})=0.$

(2.3)

Clearly, (2.3) has a negative real root $\lambda=-r$, other roots are determined by the follow equation:

$\lambda^{2}+P_{1}\lambda+P_{0}+Q_{0}e^{-\lambda\tau}=0, P_{0}=r_{2}(D+r_{1}), P_{1}=D+r_{1}+r_{2}, Q_{0}=-\frac{a_{2}r^{2}D}{a^{2}+mr^{2}}.$

Denote

$f(\lambda)=\lambda^{2}+(r_{2}+h+r_{1})\lambda+r_{2}(D+r_{1})-\frac{a_{2}D(x_{0})^{2}}{1+m(x_{0})^{2}}e^{-\lambda\tau}, $

if $(H1)$ holds, it is easy to show that, for $\lambda$ real,

$f(0)=-\frac{a_{2}Dr^{2}-r_{2}(a^{2}+mr^{2})(D+r_{1})}{a^{2}+mr^{2}}<0, \lim\limits_{\lambda\rightarrow\infty}f(\lambda)=\infty.$

Hence, $f(\lambda)=0$ has at least one positive real root. Therefore, if $(H1)$ holds, the equilibrium $E_{1}(\frac{r}{a}, 0, 0)$ is unstable. If $a_{2}Dr^{2}<r_{2}(a^{2}+mr^{2})(D+r_{1})$, it is readily seen from (2.3) that $E_{1}(\frac{r}{a}, 0, 0)$ is locally asymptotically stable when $\tau=0$. In this case, it is easy to show that

$P^{2}_{1}-2P_{0}=(D+r_{1})^{2}+r_{2}^{2}>0, \\ P_{0}^{2}-Q^{2}_{0}=(r_{2}(r_{1}+D)+\frac{a_{2}Dr^{2}}{a^{2}+mr^{2}})(\frac{r_{2}(r_{1}+D)(a^{2}+mr^{2})-a_{2}Dr^{2}}{a^{2}+mr^{2}})>0. $

By Theorem $3.4.1$ in Kuang [30], we see that, $E_{1}(\frac{r}{a}, 0, 0)$ is asymptotically stable for all $\tau>0$.

The characteristic equation of system at the coexistence equilibrium $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is of the form

$\lambda^{3}+p_{2}\lambda^{2}+p_{1}\lambda+p_{0}+(q_{1}\lambda+q_{0})e^{-\lambda\tau}=0,$

(2.4)

where

$p_{0}=r_{2}(D+r_{1})(\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*})^{2})^{2}}+2ax^{*}-r), \\ p_{1}=r_{2}(D+r_{1})+(D+r_{1}+r_{2})(\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*})^{2})^{2}}+2ax^{*}-r), \\ p_{2}=D+r_{1}+r_{2}-r+\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*})^{2})^2}+2ax^{*}, \\ q_{0}=\frac{r_{2}(D+r_{1})(r-2ax^{*})}{1+m(x^{*})^{2}}, q_{1}=-\frac{r_{2}(D+r_{1})}{1+m(x^{*})^{2}}. $

When $\tau=0$, (2.4) becomes $\lambda^{3}+p_{2}\lambda^{2}+(p_{1}+q_{1})\lambda+p_{0}+q_{0}=0.$ Clearly, by calculation we derive that

$p_{0}+q_{0}=r_{2}(D+r_{1})(\frac{3(r-2ax^{*})}{1+m(x^{*})^{2}}+2ax^{*}-r), \\ p_{1}+q_{1}=(D+r_{1}+r_{2})(\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*})^{2})^{2}}+2ax^{*}-r)+r_{2}(D+r_{1})(\frac{m(x^{*})^{2}}{1+m(x^{*})^{2}}). $

Hence, by the Routh-Hurwitz theorem, the equilibrium $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is locally asymptotically stable when $\tau=0$ if the following holds: $p_2(p_{1}+q_{1})-(p_{0}+q_{0})>0$ and $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is unstable if the inequality in $(H2)$ is reversed.

If $i\omega (w>0)$ is a solution of (2.4), separating real and imaginary parts, we have

$-\omega^{3}+p_{1}\omega=q_{0}\sin\omega\tau-q_{1}\omega\cos\omega\tau , p_{2}\omega^{2}-p_{0}=q_{0}\cos\omega\tau+q_{1}\omega\sin\omega\tau .$

(2.5)

Squaring and adding the two equations of (2.5), it follows that

$\omega^{6}+(p^{2}_{2}-2p_{1})\omega^{4}+(p^{2}_{1}-2p_{0}p_{2}-q^{2}_{1})\omega^{2}+p^{2}_{0}-q^{2}_{0}=0 ,$

(2.6)

it is easy to show that

$\begin{aligned} p^{2}_{2}-2p_{1}&=(D+r_{1})^{2}+r^{2}_{2}+(\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*})^{2})^{2}}+2ax^{*}-r)^{2}, \\ p^{2}_{1}-2p_{0}p_{2}-q^{2}_{1}&=r^{2}_{2}(D+r_{1})^{2}(1-\frac{1}{(1+m(x^{*}))^{2}})+r^{2}_{2}(\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*}))^{2}}+2ax^{*}-r)^{2}\\ &\quad +(D+r_{1})^{2}(\frac{2a_{1}x^{*}y_{2}^{*}}{(1+m(x^{*}))^{2}}+2ax^{*}-r)^{2}, \\ p_{0}-q_{0}&=r_{2}(D+r_{1})(\frac{r-ax^{*}}{1+m(x^{*})^{2}}+2ax^{*}-r) . \end{aligned}$

Hence, if $p_{0}>q_{0}$, (2.6) has no positive real roots. Accordingly, if $(H2)$ and $p_{0}>q_{0}$ holds, then the equilibrium $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is locally asymptotically stable for all $\tau\geq0$. If $p_{0}<q_{0}, $ then (2.6) has a unique positive root $\omega_{0}$, that is, (2.5) has a pair of purely imaginary roots of the form $\pm\omega_{0}$. Denote

$\tau_{0n}=\frac{1}{\omega_{0}}\arccos\frac{q_{0}(p_{2}\omega^{2}_{0}-p_{0})+q_{1}\omega_{0}(\omega^{3}_{0}-p_{1}\omega_{0})}{q^{2}_{0}+q^{2}_{1}\omega^{2}_{0}}+\frac{2n\pi}{\omega_{0}}, n=0, 1, 2, \cdots .$

(2.7)

Noting that if $(H2)$ holds, $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is locally stable when $\tau=0$, by the general theory on characteristic equations of delay differential equations from [21] (Theorem 3.4.1), $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ remains stable for $\tau<\tau_{0}$, where $\tau_{0}=\tau_{00}$.

We claim that$ \frac{d(Re\lambda)}{d\tau}\Big|_{\tau=\tau_{0}}>0, $ this will show that there exists at least one eigenvalue with positive real part for $\tau>\tau_{0}$. Moreover, the conditions for the existence of a Hopf bifurcation [23] are then satisfied yielding a periodic solution. To this end, differentiating (2.4) with respect $\tau$, it follows that

$(3\lambda^{2}+2p_{2}\lambda+p_{1})\frac{d\lambda}{d\tau}+q_{1}e^{-\lambda\tau}\frac{d\lambda}{d\tau}-\tau(q_{1}\lambda+q_{0})e^{-\lambda\tau}\frac{d\lambda}{d\tau}=\lambda(q_{1}\lambda+q_{0})e^{-\lambda\tau} ,$

which yield

$\left(\frac{d\lambda}{d\tau}\right)^{-1}=\frac{3\lambda^{2}+2p_{2}\lambda+p_{1}}{-\lambda(\lambda^{3}+p_{2}\lambda^{2}+p_{1}\lambda+p_{0})}+\frac{q_{1}}{\lambda(q_{1}\lambda+q_{0})}-\frac{\tau}{\lambda}.$

Hence, it can be calculated that

${\rm sgn}\left\{\frac{d(Re\lambda)}{d\tau}\right\}_{\lambda=i\omega_{0}}={\rm sgn}\left\{Re\left(\frac{d\lambda}{d\tau}\right)^{-1}\right\}_{\lambda=i\omega_{0}}\\ ={\rm sgn}\left\{-\frac{(p_{1}-3\omega^{2}_{0})(\omega^{2}_{0}-p_{1})+2p_{2}(p_{0}-p_{2}\omega^{2}_{0})}{(\omega^{3}_{0}-p_{1}\omega_{0})^{2}+(p_{0}-p_{2}\omega^{2}_{0})^{2}}-\frac{q^{2}_{1}}{q^{2}_{0}+q^{2}_{1}\omega^{2}_{0}}\right\}.$

(2.8)

We derive from (2.5) that $(\omega^{3}_{0}-p_{1}\omega_{0})^{2}+(p_{0}-p_{2}\omega^{2}_{0})^{2}=q^{2}_{0}+q^{2}_{1}\omega^{2}_{0}.$ Hence, it follows that

${\rm sgn}\left\{\frac{d({\rm Re}\lambda)}{d\tau}\right\}_{\lambda=i\omega_{0}}={\rm sgn}\left\{\frac{3\omega^{4}_{0}+2(p^{2}_{2}-2p_{1})\omega^{2}_{0}+p^{2}_{1}-2p_{0}p_{2}-q^{2}_{1}}{q^{2}_{0}+q^{2}_{1}\omega^{2}_{0}}\right\}>0.$

Therefore, the transversal condition holds and a Hopf bifurcation occurs at $\omega=\omega_{0}, \tau=\tau_{0}$.

We have the following theorem on stability and Hopf bifurcation.

Theorem 2.1 For system (1.2), we have the following:

(ⅰ) The equilibrium $E_{0}(0, 0, 0)$ is always unstable;

(ⅱ) The predator-extinction equilibrium $E_{1}(\frac{r}{a}, 0, 0)$ is locally asymptotically stable if $a_{2}Dr^{2}<r_{2}(D+r_{1})(a^{2}+mr^{2})$; $E_{1}(\frac{r}{a}, 0, 0)$ is unstable if $a_{2}Dr^{2}<r_{2}(D+r_{1})(a^{2}+mr^{2})$;

(ⅲ) Let $(H1)$ and $p_2(p_{1}+q_{1})>p_{0}+q_{0}$ holds. If $p_{0}>q_{0}$, then the coexistence equilibrium $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is locally asymptotically stable for all $\tau\geq0$; if $p_{0}<q_{0}$, then there exist a positive number $\tau_{0}$ such that $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ is locally asymptotically stable if $0<\tau<\tau_{0}$ and is unstable if $\tau>\tau_{0}$. Futher, system (1.2) undergoes a Hopf bifurcation at $E^{*}(x^{*}, y^{*}_{1}, y^{*}_{2})$ when $\tau=\tau_{0}$.

In the following part, we will investigate the direction of these Hopf bifurcations and stability of bifurcated periodic solutions arising through Hopf bifurcations based on the normal form theory and center manifold theorem introduced by Hassard et al. [24]. We have achieved the conditions under which a family of periodic solutions bifurcated from the positive equilibrium of system (1.2) when the delay crosses through the critical value $\tau_{0n}$.

Which determine the direction of Hopf bifurcation and stability of bifurcated periodic solutions of system (1.2) at the critical value $\tau_{0n}$:

$g_{20}=2\tau_{0}\bar{G}\hat{g}_{20}, \; g_{11}=2\tau_{0}\bar{G}\hat{g}_{11}, \; g_{02}=2\tau_{0}\bar{G}\hat{g}_{02}, \; g_{21}=2\tau_{0}\bar{G}\hat{g}_{21}, $

where

$\begin{aligned} \hat{g}_{20}&=-(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho^{2}_{1}\bar{\rho}^{*}_{1}-\frac{2a_{1}x^{*}}{c^{2}_{2}}\rho_{1}\bar{\rho}^{*}_{1}+\frac{a_{2}c_{3}}{c^{3}_{2}}e^{-2i\omega_{0}\tau_{0}}\rho^{2}_{1}\bar{\rho}^{*}_{2}+\frac{2a_{2}x^{*}}{c^{2}_{2}}e^{-2i\omega_{0}\tau_{0}}\rho_{1}\bar{\rho}^{*}_{2}, \\ \hat{g}_{11}&=-(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho_{1} \bar{\rho}_{1}\bar{\rho}^{*}_{1}-\frac{a_{1}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1})\bar{\rho}^{*}_{1}\\&\quad +\frac{a_{2}c_{3}}{c^{3}_{2}}e^{-2i\omega_{0}\tau_{0}}\rho^{2}_{1} \bar{\rho}^{*}_{2}+\frac{a_{2}x^{*}}{c^{2}_{2}}e^{-2i\omega_{0}\tau_{0}}(\rho_{1}+\bar{\rho}_{1})\bar{\rho}^{*}_{2}, \\ \hat{g}_{02}&=-a\bar{\rho}^{2}_{1}\bar{\rho}^{*}_{1}-\frac{a_{1}c_{3}}{c^{3}_{2}}\bar{\rho}^{2}_{1}\bar{\rho}^{*}_{1}-\frac{2a_{1}x^{*}}{c^{2}_{2}}\bar{\rho}_{1}\bar{\rho}^{*}_{1} +\frac{a_{2}c_{3}}{c^{3}_{2}}e^{-2i\omega_{0}\tau_{0}}\bar{\rho}^{2}_{1}\bar{\rho}^{*}_{2}+\frac{2a_{2}x^{*}}{c^{2}_{2}}e^{-2i\omega_{0}\tau_{0}}\bar{\rho}_{1}\bar{\rho}^{*}_{2}, \\ \hat{g}_{21}&=\frac{6mx^{*}(c_{3}+c_{2}y^{*}_{2})\rho^{2}_{1}\bar{\rho}_{1}}{c^{4}_{2}}(a_{1}\bar{\rho}^{*}_{1}-a_{2}\bar{\rho}^{*}_{2}e^{-i\omega_{0}\tau_{0}}) \\&\quad +\frac{c_{2}-4mx^{*}}{c^{3}_{2}}(\rho^{2}_{1}+2\rho_{1}\bar{\rho}_{1})(a_{2}\bar{\rho}^{*}_{2}e^{-i\omega_{0}\tau_{0}}-a_{1}\bar{\rho}^{*}_{1})\\&\quad -\frac{a_{1}x^{*}\bar{\rho}^{*}_{1}}{c^{2}_{2}}(W^{(3)}_{20}(0)\bar{\rho}_{2}+W^{(3)}_{20}(0)+2W^{(1)}_{11}(0)\rho_{2}+2W^{(3)}_{11}(0)) \\&\quad +\frac{a_{2}x^{*}\bar{\rho}^{*}_{2}}{c^{2}_{2}}(W^{(3)}_{20}(-1)\bar{\rho}_{1}e^{i\omega_{0}\tau_{0}}\\ &\quad +W^{(3)}_{20}(-1)e^{i\omega_{0}\tau_{0}}+2W^{(1)}_{11}(-1)\rho_{1}e^{i\omega_{0}\tau_{0} }+2W^{(3)}_{11}(-1)e^{i\omega_{0}\tau_{0}})+\bar{\rho}_{1}\bar{\rho}^{*}_{1}(W^{(1)}_{20}(0)\\&\quad +2W^{1}_{11}(0))(-a-\frac{a_{1}c_{3}}{c^{3}_{2}})+\bar{\rho}_{1}\bar{\rho}^{*}_{2}(W^{(1)}_{20}(-1)+2W^{1}_{11}(-1))\frac{a_{2}c_{3}}{c^{3}_{2}} e^{i\omega_{0}\tau_{0}}, \end{aligned}$

where the terms $\bar{G}, W^{1}_{11}(0), W^{3}_{11}(0), W^{1}_{11}(-1), W^{3}_{11}(-1), W^{1}_{20}(0), W^{3}_{20}(0), W^{1}_{20}(-1), W^{3}_{20}(-1)$ are calculated in Appendix. Now, we exhibit these coefficients we can evaluate the following values

$ C_{1}(0)=\frac{i}{2\omega_{0}\tau_{0}}(g_{11}g_{20}-2|g_{11}|^{2}-\frac{1}{3}|g_{02}|^{2})+\frac{g_{21}}{2}, \\ \mu_{2}=-\frac{{\rm Re}C_{1}(0)}{{\rm Re}\lambda'(\tau_{0})}, \; \beta_{2}=2{\rm Re}\{C_{1}(0)\}, \; T_{2}=\frac{{\rm Im}C_{1}(0)+\mu_{2}{\rm Im}\lambda'(\tau_{0})}{\omega_{0}\tau_{0}}. $

Now using the quantities above, the properties of the Hopf bifurcation is determined by the following results.

Theorem 3.1 For system (1.2), we have the following:

(ⅰ) If $\mu_{2}>0 (\mu_{2}<0)$, then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for $\tau>\tau_{0} (\tau<\tau_{0})$.

(ⅱ) $\beta_{2}$ determines the stability of the bifurcating periodic solutions: the bifurcating solutions on the center manifold are stable (unstable) if $\beta_{2}<0 (\beta_{2}>0)$.

(ⅲ) And $T_{2}$ determines the period of the bifurcating periodic solutions: the period increase (decrease) if $T_{2}>0 (T_{2}<0)$. From the discussion in Section 2, we know that $\operatorname{Re}({\lambda }'({{\tau }_{0}}))>0$, therefore we have the following result.

Theorem 3.2 The direction of the Hopf bifurcation of system (1.2) at the origin when $\tau=\tau_{j} (j=0, 1, 2, \cdots)$ is supercritical (subcritical) and the bifurcating periodic solutions on the center manifold are stable (unstable) if Re$\{C_{1}(0)\}<0(>0)$; particularly, when $\tau=\tau_{0}$, the stability of the bifurcating periodic solutions is the same as that on the center manifold.

We now give an example to illustrate algorithm for determining the existence of Hopf bifurcation in Section 2 and the direction and stability of Hopf bifurcation in Section 3. Let $r=25, r_{1}=\frac{2}{19}, r_{2}=\frac{2}{19}, a=16, a_{1}=4, a_{2}=\frac{7}{2}, m=16, D=\frac{1}{2}$, it is easy to show that

$a_{2}Dr^{2}-r_{2}(a^{2}+mr^{2})(D+r_{1})=440.3206>0.$

Hence, system (1.2) has a unique coexistence equilibrium

$E^{*}(\sqrt{\frac{92}{1055}}, \frac{3325}{92}\sqrt{\frac{92}{1055}}-\frac{2128}{1055}, \frac{63175}{368}\sqrt{\frac{92}{1055}}-\frac{10108}{1055}).$

By calculation, we have

$ p_{0}-q_{0}\approx-0.3258<0, \;\;p_{2}(p_{1}+q_{1})-(p_{0}+q_{0})\approx1.6236>0, \\ \tau\approx2.3940, \;\; C_{1}(0)\approx-14.3078+2.4964i, \\ \mu_{2}\approx1.6612, \;\; \beta_{2}\approx-28.6155, \;\; T_{2}\approx24.9693. $

By Theorem 3.1, we obtain the result that the Hopf bifurcation of system (1.2) occurring at $\tau_{0}=2.3940$ is subcritical and the bifurcating periodic solution exist when $\tau$ cross $\tau_{0}$ to the left and the bifurcating periodic solution is stable, and the equilibrium $E^{*}$ is locally asymptotically stable if $0<\tau<\tau_{0}$ as is illustrated by computer simulations in Fig. 1 and is unstable if $\tau>\tau_{0}$ as is illustrated by computer simulations in Fig. 2. System (4.1) will show the complicated dynamical behaviors, here, we choose the initial condition $x^{0}=0.6, y^{0}_{1}=12, y^{0}_{2}=42$ in our simulations.

In this paper, we conclude the stability properties of this system based on the improvement of a stage structure predator-prey system proposed in [7]. By deriving the equation describing the flow on the center manifold, we determine the direction of the Hopf bifurcations and the stability of the bifurcating periodic solutions. It shows that under some conditions, the time delay due to the gestation of the mature predator may destabilize the coexistence equilibrium of the system and cause the population to fluctuate.

In addition, with the development of the society, the ordinary differential equations of the ecosystem is no longer satisfy mankind's significance. Correspondingly, the economic profit is becoming a important factor for governments, merchants and even every citizen. Therefore, it is necessary that we should study the differential-algebraic system centered on economic profit and commercial harvesting, and investigates the efforts of economic profit on the dynamics of the predator-prey system.

We denote $\bar{x}(t)=x(t)-x^{*}, \bar{y}_{1}(t)=y_{1}(t)-y^{*}_{1}, \bar{y}_{2}(t)=y_{2}(t)-y^{*}_{2}$. By the transformation, the system (1.2) can be taken the form of

$\dot{\bar{x}}=c\bar{x}-\frac{r_{2}(D+r_{1})}{D}{\bar{y}_{2}}-a\bar{x}^{2}-a_{1}h(\bar{x}, \bar{y}_{2}) , \\ \dot{\bar{y}}_{1}=c_{2}\bar{x}(t-\tau)+\frac{a_{1}}{a_{2}}\frac{r_{2}(D+r_{1})}{D}\bar{y}_{2}(t-\tau)-(D+r_{1})\bar{y}_{1}+a_{2}h(\bar{x}(t-\tau), \bar{y}_{2}(t-\tau)), \\ \dot{\bar{y}}_{2}=D\bar{y}_{1}-r_{2}\bar{y}_{2},$

(1)

where

$ h(x, y)=\frac{c_{1}xy(x+2x^{*})+x^{2}(c_{3}-2mx^{*}y^{*}_{2}x)}{c^{2}_{1}(c_{1}+m(x^{2}+2xx^{*}))} , \\ c_{1}=1+m(x^{*})^{2}, \;\; c_{2}=\frac{2a_{1}x^{*}y_{2}^{*}}{c^{2}_{1}}, \;\;c_{3}=y_{2}^{*}(1-3m(x^{*})^{2}), \;\; c=r-2ax^{*}-c_{2} . $

Then we imitate [25] to do the presentations as follow.

Let $t=s\tau, \bar{x}(s\tau)=\hat{x}(s), $ $\bar{y}_{1}(s\tau)=\hat{y}_{1}(s), $ $\bar{y}_{2}(s\tau)=\hat{y}_{2}(s), $ $\tau=\tau_{0}+\mu, \mu\in R$, $\tau_{0}$ is defined by (2.7), then system (1.2) can be transformed as an FDE in $C=C([-1,0], R^{3})$, we will denote $x=\hat{x}, y_{1}=\hat{y}_{1}, y_{2}=\hat{y}_{2}$,

$\dot{x}(t)=(\tau_{0}+\mu)\Big(cx(t)-\frac{r_{2}(D+r_{1})}{D}y_{2}(t)-ax^{2}(t)-a_{1}h(x, y_{2})\Big) , \\ \dot{y}_{1}(t)=(\tau_{0}+\mu)\Big(c_{2}x(t-1)-\frac{a_{1}}{a_{2}}\frac{r_{2}(D+r_{1})}{D}y_{2}(t-1)-(D+r_{1})y_{1}(t)+a_{2}h(x(t-1), y_{2}(t-1)\Big), \\ \dot{y}_{2}(t)=(\tau_{0}+\mu)\Big(Dy_{1}-r_{2}y_{2}\Big).$

(2)

For $\phi=(\phi_{1}, \phi_{2}, \phi_{3})^{T}\in C([-1,0], R^{3})$, system (1) is equivalent to the following operator equation

$\dot{U}(t)=L(\mu)U_{t}+F(\mu)U_{t},$

(3)

where $U=(x, y_{1}, y_{2}), U_{t}=U(t+\theta)$ for $\theta\in [-1,0]$.

$\begin{array}{l}F(\mu, \phi)=\left(\begin{array}{c} f_{1} \\ f_{2} \\ 0 \end{array}\right)=(\tau_{0}+\mu)\left(\begin{array}{c} -a\phi_{1}^{2}(0)-a_{1}h(\phi_{1}(0), \phi_{3}(0)) \\ a_{2}h(\phi_{1}(-1), \phi_{3}(-1)) \\ 0 \end{array}\right). \end{array}$

We define

$L_{\mu}\phi=(\tau_{0}+\mu)\left(\begin{array}{ccc} c 0 -\frac{r_{2}(D+r_{1})}{D} \\ 0 -(D+r_{1}) 0 \\ 0 D -r_{2} \end{array}\right)\phi(0)+(\tau_{0}+\mu)\left(\begin{array}{ccc} 0 0 0 \\ c_{2} 0 \frac{a_{1}}{a_{2}}\frac{r_{2}(D+r_{1})}{D} \\ 0 0 0 \end{array}\right)\phi(-1).$

By the Riesz representation theorem, there exists a $3\times3$ matrix functions $\eta(\theta, \mu):[-1,0]\rightarrow R^{3\times3}$ whose elements are of bounded variation such that

$L_{\mu}\phi=\displaystyle\int^{0}_{-1}[d\eta(\theta, \mu)]\phi(\theta) \;\mathrm{for}\; \phi\in C([-1,0], R^{3}),$

(4)

in fact, we chose $\eta(\theta, \mu)=B\delta(\theta)+K\delta(\theta+1), $ then $(3)$ is satisfied.

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

$A(\mu, \phi)= \begin{cases}\frac{d\phi(\theta)}{d\theta}, \quad \ \ -1\leq\theta<0, \\ \displaystyle\int^{0}_{-1}[d\eta(s, \mu)]\phi(s), \quad \ \ \theta=0, \end{cases}\\ R(\mu, \phi)= \begin{cases}0, \quad \ \ -1\leq\theta<0, \\ F(\mu, \theta), \quad \ \ \theta=0. \end{cases} $

Then system (2) is equivalent to the following operator equation:

$\dot{U}(t)=A(\mu)U_{t}+R(\mu)U_{t} .$

For $\psi\in C([-1,0], (R^{3})^{*})$, the adjoint operator $A^{*}$ of the $A$ is defined as

$A^{*}\psi(s)= \begin{cases}-\frac{d\psi(s)}{ds}, \quad \ \ -1\leq \mathrm{s}<0, \\ \displaystyle\int^{0}_{-1}[d\eta^{T}(t, 0)]\phi(-t), \quad \ \ s=0. \end{cases}$

For $\phi\in C([-1,0], R^{3})$ and $\psi\in C([-1,0], (R^{3})^{*})$, define a bilinear form

$\langle\phi, \psi\rangle=\bar{\psi}^{T}(0)\phi(0)-\int^{0}_{-1}\int^{\theta}_{\xi=0}\bar{\psi}^{T}(\xi-\theta)d\eta(\theta)\phi(\xi)d\xi, $

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

From the discussion in Section 2, we know that $q(\theta)$ and $q^{*}(s)$ are eigenvectors of $A$ and $A^{*}$ corresponding to $iw_{0}\tau_{0}$ and $-iw_{0}\tau_{0}$, respectively. We defined $q(\theta)=(\rho_{1}, \rho_{2}, 1)^{T}e^{i\omega_{0}\tau_{0}\theta}, $ then $A(0)q(\theta)=iw_{0}t_{0}q(\theta)$, it follows the above definition that

$\begin{array}{l}(\tau_{0})\left(\begin{array}{ccc} c-iw_{0}\tau_{0} 0 -\frac{r_{2}(D+r_{1})}{D} \\ 0 -(D+r_{1})-iw_{0}\tau_{0} 0 \\ 0 D -r_{2}-iw_{0}\tau_{0} \end{array}\right)q(0)=\left(\begin{array}{ccc} 0 \\ 0 \\ 0 \end{array}\right), \end{array}$

then, we can easily get $q(\theta)=(\rho_{1}, \rho_{2}, 1)^{T}e^{iw_{0}\tau_{0}\theta}=q(0)e^{iw_{0}\tau_{0}\theta}$, where $\rho_{1}=\frac{r_{2}(D+r_{1})}{D(c-i\omega_{0})}, $ $\rho_{2}=\frac{r_{2}+i\omega_{0}}{D}$. Similarly, by definition of $A^{*}$,

$\begin{array}{l}(\tau_{0})\left(\begin{array}{ccc} c+iw_{0}\tau_{0} 0 -\frac{r_{2}(D+r_{1})}{D} \\ 0 -(D+r_{1})+iw_{0}\tau_{0} 0 \\ 0 D -r_{2}+iw_{0}\tau_{0} \end{array}\right)q^{*}(0)=\left(\begin{array}{ccc} 0 \\ 0 \\ 0 \end{array}\right) \end{array}$

and $q^{*}(\theta)=\bar{G}(\rho^{*}_{1}, \rho^{*}_{2}, 1)^{T}e^{iw_{0}\tau_{0}\theta}=q^{*}(0)e^{iw_{0}\tau_{0}\theta}$, where $\rho^{*}_{1}=\frac{-c_{2}De^{-i\omega_{0}}}{(c+i\omega)(D+r_{1}-i\omega)}, $ $\rho^{*}_{2}=\frac{D}{D+r_{1}-i\omega_{0}}, $ in order to have $\langle q^{*}(s), q(\theta)\rangle=1$, we evaluate the value $\bar{G}$, By the definition of the bilinear inner product, we have

$\begin{array}{l} \quad \langle q^{*}(0), q(0)\rangle\\ =\bar{G}(\bar{\rho}^{*}_{1}, \bar{\rho}^{*}_{2}, 1)(\rho_{1}, \rho_{2}, 1)^{T}-\int^{0}_{-1}\displaystyle\int^{\theta}_{\xi=0}\bar{G}(\bar{\rho}^{*}_{1}, \bar{\rho}^{*}_{2}, 1)e^{-iw_{0}(\xi-\theta)}d\eta(\theta, 0)(\rho_{1}, \rho_{2}, 1)e^{iw_{0}\xi}d\xi \\ =\bar{G}\left(\rho_{1}\bar{\rho}^{*}_{1}+\rho_{2}\bar{\rho}^{*}_{2}+1-\displaystyle\int^{0}_{-1}(\rho^{*}_{1}, \rho^{*}_{2}, 1)e^{-iw_{0}\theta}\theta\eta(\theta)(\rho_{1}, \rho_{2}, 1)^{T}\right)\\ =\bar{G}\left(1+\rho_{1}\bar{\rho}^{*}_{1}+\rho_{2}\bar{\rho}^{*}_{2}+(c_{2}\rho_{1}+\frac{a_{1}r_{2}(D+r_{1})}{a_{2}D})\tau_{0}e^{-i\omega_{0}\tau_{0}}\bar{\rho}^{*}_{2}\right). \end{array}$

Thus, we can choose $\bar{G}$ as

$\bar{G}=\frac{1}{1+\rho_{1}\bar{\rho}^{*}_{1}+\rho_{2}\bar{\rho}^{*}_{2}+(c_{2}\rho_{1}+\frac{a_{1}r_{2}(D+r_{1})}{a_{2}D})\tau_{0}e^{-i\omega_{0}\tau_{0}}\bar{\rho}^{*}_{2}},$

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

According to the algorithms given in [22] and utilizing a computation process akin to that of Wei and Ruan [25], we compute the coordinates to describe the center manifold $C_{0}$ at $\mu=0$. We define

$z(t)=\langle q^{*}, x_{t}\rangle, \;\; W(t, \theta)=x_{t}(\theta)-2{\rm Re}{z(t)q(\theta)},$

(5)

thus, we have $ W(t, \theta)=W(z(t), \bar{z}(t), \theta), $ where

$W(z(t), \bar{z}(t), \theta)=W(z, \bar{z}) =\left(\begin{array}{ccc} W^{(1)}_{20}(\theta) \\ W^{(2)}_{20}(\theta) \\ W^{(3)}_{20}(\theta) \end{array}\right)\frac{z^{2}}{2}+\left(\begin{array}{ccc} W^{(1)}_{11}(\theta) \\ W^{(2)}_{11}(\theta) \\ W^{(3)}_{11}(\theta) \end{array}\right)z\bar{z}+\left(\begin{array}{ccc} W^{(1)}_{02}(\theta) \\ W^{(2)}_{02}(\theta) \\ W^{(3)}_{02}(\theta) \end{array}\right)\frac{\bar{z}^{2}}{2}+\cdots, $

in fact, $z$ and $\bar{z}$ are local coordinates for center manifold $C_{0}$ in the direction of $q^{*}$ and $\bar{q}^{*}$. Noting that $W$ is also real if $x_{t}$ is real, For solution $x_{t}\in C_{0}$ of (3.4),

$\dot{z}(t)=i\omega_{0}\tau_{0}z+\bar{q}^{*}(\theta)F(0, W(z(t), \bar{z}(t), \theta))+2{\rm Re}{zq(\theta)} :=i\omega_{0}\tau_{0}z+\bar{q}^{*}(\theta)F_{0}=i\omega_{0}\tau_{0}z+g(z, \bar{z}). $

By using (5), we have $x_{t}=(\phi_{1}, \phi_{2}, \phi_{3})=W(t, \theta)+zq(\theta)+\bar{z}\bar{q}(\theta)$ and $q(\theta)=(\rho_{1}, \rho_{2}, 1)^{T}e^{iw_{0}\tau_{0}\theta}, $ then we achieve

$\begin{equation*}\begin{array}{l}\phi_{1}(0)=\rho_{1}z+\bar{\rho}_{1}\bar{z}+W^{1}_{20}(0)\frac{z^{2}}{2}+W^{1}_{11}(0)z\bar{z}+W^{1}_{02}(0)\frac{\bar{z}^{2}}{2}, \\ \phi_{2}(0)=\rho_{2}z+\bar{\rho}_{2}\bar{z}+W^{2}_{20}\frac{z^{2}}{2}+W^{2}_{11}(0)z\bar{z}+W^{2}_{02}(0)\frac{\bar{z}^{2}}{2}, \\ \phi_{3}(0)=z+\bar{z}+W^{3}_{20}\frac{z^{2}}{2}+W^{3}_{11}(0)z\bar{z}+W^{3}_{02}(0)\frac{\bar{z}^{2}}{2}, \\ \phi_{1}(-1)=\rho_{1}ze^{-iw_{0}\tau_{0}\theta}+\bar{\rho}_{1}\bar{z}e^{iw_{0}\tau_{0}\theta}+W^{1}_{20}(-1)\frac{z^{2}}{2}+W^{1}_{11}(-1)z\bar{z}+W^{1}_{02}(-1)\frac{\bar{z}^{2}}{2}, \\ \phi_{2}(-1)=\rho_{1}ze^{-iw_{0}\tau_{0}\theta}+\bar{\rho}_{1}\bar{z}e^{iw_{0}\tau_{0}\theta}+W^{2}_{20}(-1)\frac{z^{2}}{2}+W^{2}_{11}(-1)z\bar{z}+W^{2}_{02}(-1)\frac{\bar{z}^{2}}{2}, \\ \phi_{3}(-1)=ze^{-iw_{0}\tau_{0}\theta}+\bar{z}e^{iw_{0}\tau_{0}\theta}+W^{3}_{20}(-1)\frac{z^{2}}{2}+W^{3}_{11}(-1)z\bar{z}+W^{3}_{02}(-1)\frac{\bar{z}^{2}}{2}. \end{array}\tag{6} \end{equation*}$

(6)

From the above presentation, we have

$g(z, \bar{z})=g_{20}\frac{z^{2}}{2}+g_{11}z\bar{z}+g_{02}\frac{\bar{z}^{2}}{2}+g_{21}\frac{z^{2}\bar{z}}{2}+\cdots=\bar{q}^{*}F_{0}(0, x_{t})=\bar{G}\tau_{0}(\bar{\rho}^{*}_{1}, \bar{\rho}^{*}_{2}, 1)(f_{1}, f_{2}, 0)^{T},$

(7)

where

$\begin{aligned}{l}f_{1}&=-a\phi_{1}^{2}(0)-a_{1}h(\phi_{1}(0), \phi_{3}(0))\\ &=-(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\phi^{2}_{1}(0)-\frac{2a_{1}x^{*}}{c^{2}_{2}}\phi_{1}(0)\phi_{3}(0)\\&\quad -\frac{a_{1}(c_{2}-4mx^{*})}{c^{3}_{2}}\phi^{2}_{1}(0)\phi_{3}(0)+O((|\phi_{1}(0)|+|\phi_{3}(0)|)^{3}), \\ f_{2}&=a_{2}h(\phi_{1}(-1), \phi_{3}(-1))\\ &=\frac{a_{2}c_{3}}{c^{3}_{2}}\phi^{2}_{1}(0)+\frac{2a_{2}x^{*}}{c^{2}_{2}}\phi_{1}(0)\phi_{3}(0)\\ &\quad +\frac{a_{2}(c_{2}-4mx^{*})}{c^{3}_{2}}\phi^{2}_{1}(0)\phi_{3}(0)+O((|\phi_{1}(-1)|+|\phi_{3}(-1)|)^{3}). \end{aligned}$

(8)

Thus, By comparing the coefficients, we obtain $g_{20}=2\tau_{0}\bar{G}\hat{g}_{20}, $ $g_{11}=2\tau_{0}\bar{G}\hat{g}_{11}, $ $g_{02}=2\tau_{0}\bar{G}\hat{g}_{02}, $ $g_{21}=2\tau_{0}\bar{G}\hat{g}_{21}, $ where

$ \begin{aligned}\hat{g}_{20}&=-(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho^{2}_{1}\bar{\rho}^{*}_{1}-\frac{2a_{1}x^{*}}{c^{2}_{2}}\rho_{1}\bar{\rho}^{*}_{1} +\frac{a_{2}c_{3}}{c^{3}_{2}}e^{-2i\omega_{0}\tau_{0}}\rho^{2}_{1}\bar{\rho}^{*}_{2} +\frac{2a_{2}x^{*}}{c^{2}_{2}}e^{-2i\omega_{0}\tau_{0}}\rho_{1}\bar{\rho}^{*}_{2}, \\ \hat{g}_{11}&=-(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho_{1} \bar{\rho}_{1}\bar{\rho}^{*}_{1}-\frac{a_{1}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1})\bar{\rho}^{*}_{1}+\frac{a_{2}c_{3}}{c^{3}_{2}} e^{-2i\omega_{0}\tau_{0}}\rho^{2}_{1} \bar{\rho}^{*}_{2}+\frac{a_{2}x^{*}}{c^{2}_{2}}e^{-2i\omega_{0}\tau_{0}}(\rho_{1}+\bar{\rho}_{1})\bar{\rho}^{*}_{2}, \\ \hat{g}_{02}&=-a\bar{\rho}^{2}_{1}\bar{\rho}^{*}_{1}-\frac{a_{1}c_{3}}{c^{3}_{2}}\bar{\rho}^{2}_{1}\bar{\rho}^{*}_{1} -\frac{2a_{1}x^{*}}{c^{2}_{2}}\bar{\rho}_{1}\bar{\rho}^{*}_{1} +\frac{a_{2}c_{3}}{c^{3}_{2}}e^{-2i\omega_{0}\tau_{0}}\bar{\rho}^{2}_{1}\bar{\rho}^{*}_{2} +\frac{2a_{2}x^{*}}{c^{2}_{2}}e^{-2i\omega_{0}\tau_{0}}\bar{\rho}_{1}\bar{\rho}^{*}_{2}, \\ \hat{g}_{21}&=\frac{6mx^{*}(c_{3}+c_{2}y^{*}_{2})\rho^{2}_{1}\bar{\rho}_{1}}{c^{4}_{2}}(a_{1} \bar{\rho}^{*}_{1}-a_{2}\bar{\rho}^{*}_{2}e^{-i\omega_{0}\tau_{0}}) +\frac{c_{2}-4mx^{*}}{c^{3}_{2}}(\rho^{2}_{1}+2\rho_{1}\bar{\rho}_{1})(a_{2}\bar{\rho}^{*}_{2}e^{-i\omega_{0}\tau_{0}}-a_{1}\bar{\rho}^{*}_{1})\\ &\quad -\frac{a_{1}x^{*}\bar{\rho}^{*}_{1}}{c^{2}_{2}}(W^{(3)}_{20}(0)\bar{\rho}_{2}+W^{(3)}_{20}(0)+2W^{(1)}_{11}(0)\rho_{2}+2W^{(3)}_{11}(0)) +\frac{a_{2}x^{*}\bar{\rho}^{*}_{2}}{c^{2}_{2}}(W^{(3)}_{20}(-1)\bar{\rho}_{1}e^{i\omega_{0}\tau_{0}}\\&\quad +W^{(3)}_{20}(-1)e^{i\omega_{0}\tau_{0}}+2W^{(1)}_{11}(-1)\rho_{1}e^{i\omega_{0}\tau_{0} }+2W^{(3)}_{11}(-1)e^{i\omega_{0}\tau_{0}})\\&\quad +\bar{\rho}_{1}\bar{\rho}^{*}_{1}(W^{(1)}_{20}(0)+2W^{1}_{11}(0))(-a-\frac{a_{1}c_{3}}{c^{3}_{2}})+\bar{\rho}_{1}\bar{\rho}^{*}_{2}(W^{(1)}_{20}(-1)+2W^{1}_{11}(-1))\frac{a_{2}c_{3}}{c^{3}_{2}} e^{i\omega_{0}\tau_{0}}. \end{aligned}$

Next, we calculate the $W^{i}_{20}(0), W^{i}_{11}(0), W^{j}_{20}(-1), W^{j}_{11}(-1) (i=1, 2, 3; j=1, 2)$ in $g_{21}$. From $(4)$ and $(5)$, it is easily to derive:

$\begin{array}{l}W^{'}=\begin{cases}AW-2Re{\bar{q}^{*}(0)\bar{F}q(\theta)}, \quad \ \ -1\leq\theta<0, \\ AW-2Re{\bar{q}^{*}(0)\bar{F}q(\theta)}, \quad \ \ \theta=0, \end{cases}:=AW+H(z, \bar{z}, \theta) , \end{array}$

where

$H(z, \bar{z}, \theta)=H_{20}(\theta)\frac{z^{2}}{2}+H_{11}(\theta)z\bar{z}+H_{02}(\theta)\frac{\bar{z}^{2}}{2}+\cdots.$

(9)

Substituting $(6)$ to $(5)$, we obtain

$\left(A-2i\omega_{0}\tau_{0}\right)W_{20}=-H_{20}(\theta),$

(10)

$AW_{11}(\theta)=-H_{11}(\theta) .$

(11)

According the previous definition, we know that for $\theta \in [-1, 0), $

$H(z, \bar{z}, \theta)=-\bar{q}^{*}(0)F_{0}q(\theta)-q^{*}(0)\bar{F}_{0}\bar{q}(\theta)=-g(z, \bar{z})q(\theta)-\bar{g}(z, \bar{z})\bar{q}(\theta) .$

(12)

Comparing the coefficients of $(6)$ with $(8)$ gives that

$H_{20}(\theta)=-g_{20}q(\theta)-\bar{g}_{02}\bar{q}(\theta) ,$

(13)

$H_{11}(\theta)=-g_{11}q(\theta)-\bar{g}_{11}\bar{q}(\theta) .$

(14)

It follow from (7) and (10) that

$\dot{W}_{20}(\theta)=2i\omega_{0}\tau_{0}W_{20}(\theta)+g_{20}q(\theta)+\bar{g}_{02}\bar{q}(\theta).$

By the previous definition, we denote $q(\theta)=q(0)e^{i\omega_{0}\tau_{0}\theta}$, and by solving $W_{20}$, we obtain

$W_{20}(\theta)=\frac{ig_{20}}{\omega_{0}\tau_{0}}q(0)e^{i\omega_{0}\tau_{0}}+\frac{ig_{20}}{3\omega_{0}\tau_{0}}\bar{q}(0) e^{-i\omega_{0}\tau_{0}}+M_{1}e^{2i\omega_{0}\tau_{0}\theta} .$

We can also obtain the follow result by utilizing the same method.

$\dot{W}_{11}(\theta)=g_{11}(\theta)q(\theta)+\bar{g}_{11}(\theta)\bar{q}(\theta) , W_{11}(\theta)=-\frac{ig_{11}}{\omega_{0}\tau_{0}}q(0)e^{i\omega_{0}\tau_{0}}+\frac{ig_{11}}{\omega_{0}\tau_{0}}\bar{q}(0)e^{-i\omega_{0}\tau_{0}}+M_{2} .$

Now, we will try to find $E_{1}$ and $E_{2}$. From the definition of A and (7), we obtain that

$\int^{0}_{-1}d\eta(\theta)W_{20}(\theta)=2iw_{0}W_{20}-H_{20}(0)$

and

$\int^{0}_{-1}d\eta(\theta)W_{11}(\theta)=-H_{11}(0), $

where $d\eta(\theta)=\eta(\theta, 0).$

Substituting (13) and (14) and noticing that

$\left(iw_{0}\tau_{0}I-\displaystyle\int^{0}_{-1}e^{iw_{0}\tau_{0}\theta}d\eta(\theta)\right)q(0)=0, \left(-iw_{0}\tau_{0}I-\displaystyle\int^{0}_{-1}e^{-iw_{0}\tau_{0}\theta}d\eta(\theta)\right)q(0)=0,$

(15)

we obtain

$\begin{array}{l}\left(-2iw_{0}\tau_{0}I+\int^{0}_{-1}e^{2iw_{0}\tau_{0}\theta}d\eta(\theta)\right)M_{1}=-2\tau_{0}\left(\begin{array}{ccc} -(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho^{2}_{1}-\frac{2a_{1}x^{*}}{c^{2}_{2}}\rho_{1} \\ a_{2}e^{-2i\omega_{0}\tau_{0}}(\frac{c_{3}}{c^{3}_{2}}\rho^{2}_{1}+\frac{2x^{*}}{c^{2}_{2}}\rho_{1}) \\ 0 \end{array}\right) \end{array}$

which is equivalent to

$\begin{array}{l}\left(\begin{array}{ccc} 2i\omega_{0}-c 0 -\frac{r_{2}(D+r_{1})}{D} \\ -c_{2}e^{-2i\omega} 2i\omega_{0}-(D+r_{1}) -\frac{a_{1}}{a_{2}}\frac{r_{2}(D+r_{1})}{D}e^{-2i\omega} \\ 0 -D 2i\omega_{0}-r_{2} \end{array}\right)M_{1}=\left(\begin{array}{ccc} -(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho^{2}_{1}-\frac{2a_{1}x^{*}}{c^{2}_{2}}\rho_{1} \\ a_{2}e^{-2i\omega_{0}\tau_{0}}(\frac{c_{3}}{c^{3}_{2}}\rho^{2}_{1}+\frac{2x^{*}}{c^{2}_{2}}\rho_{1}) \\ 0 \end{array}\right). \end{array}$

Similarly, from (13) and (15), $-AW_{11}(\theta)=H_{11}(\theta)$ and for $\theta=0, $

$ \begin{array}{l}\displaystyle\int^{0}_{-1}d\eta(\theta)W_{11}(\theta)=g_{11}q(0)-\bar{g}_{11}\bar{q}(0)-\tau_{0}\left(\begin{array}{ccc} -(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho_{1}\bar{\rho}_{1}-\frac{a_{1}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1}) \\ \frac{a_{2}c_{3}}{c^{2}_{2}}\rho_{1}\bar{\rho}_{1}+\frac{a_{2}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1}) \\ 0 \end{array}\right), \end{array}$

from which we can get

$\begin{array}{cc}\quad \frac{g_{11}}{iw_{0}} \int^{0}_{-1}d\eta(\theta)q(0)e^{iw_{0\tau_{0}}}-\frac{\bar{g}_{11}}{iw_{0}\tau_{0}} \int^{0}_{-1}d\eta(\theta)\bar{q}(0)e^{-iw_{0\tau_{0}}} +E_{2}\int^{0}_{-1}d\eta(\theta)\\ =g_{11}q(0)-\bar{g}_{11}\bar{q}(0)-\tau_{0}\left(\begin{array}{ccc} -2(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho_{1}\bar{\rho}_{1}-\frac{2a_{1}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1}) \\ \frac{2a_{2}c_{3}}{c^{2}_{2}}\rho_{1}\bar{\rho}_{1}+\frac{2a_{2}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1}) \\ 0 \end{array}\right), \end{array}$

which is equivalent to

$\begin{array}{l}\left(\begin{array}{ccc} -c 0 -\frac{r_{2}(D+r_{1})}{D} \\ -c_{2} -(D+r_{1}) -\frac{a_{1}}{a_{2}}\frac{r_{2}(D+r_{1})}{D} \\ 0 -D -r_{2} \end{array}\right)M_{2}=\left(\begin{array}{ccc} -(a+\frac{a_{1}c_{3}}{c^{3}_{2}})\rho_{1}\bar{\rho}_{1}-\frac{a_{1}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1}) \\ \frac{a_{2}c_{3}}{c^{2}_{2}}\rho_{1}\bar{\rho}_{1}+\frac{a_{2}x^{*}}{c^{2}_{2}}(\rho_{1}+\bar{\rho}_{1}) \\ 0 \end{array}\right) \end{array}.$

From the above calculation, we can get the vector $M_{1}$ and $M_{2}$ and put them in the equations (14) and (15), we can also obtain the terms

$W^{1}_{11}(0), W^{3}_{11}(0), W^{1}_{11}(-1), W^{3}_{11}(-1), W^{1}_{20}(0), W^{3}_{20}(0), W^{1}_{20}(-1), W^{3}_{20}(-1),$

(16)

and achieve the quantity $g_{21}$. For this propose, we express each $g_{ij}$ in terms of the parameters and delay. and then, we can evaluate the following values:

$ C_{1}(0)=\frac{i}{2\omega_{0}\tau_{0}}(g_{11}g_{20}-2|g_{11}|^{2}-\frac{1}{3}|g_{02}|^{2})+\frac{g_{21}}{2}, \\ \mu_{2}=-\frac{{\rm Re} C_{1}(0)}{{\rm Re}\lambda'(\tau_{0})}, \beta_{2}=2{\rm Re} \{C_{1}(0)\}, \\ T_{2}=-\frac{{\rm Im} C_{1}(0)+\mu_{2}{\rm Im} \lambda'(\tau_{0})}{\omega_{0}\tau_{0}}. $