In recent decades, there has been a spate of interest in bioeconomic analysis of exploitation of renewable resources like fisheries, exploitation of natural resources has become a matter of concern throughout the world. Therefore, it has become imperative to ensure scientific management of exploitation of biological resources. To insure the long-term benefits of humanity, there is a wide-range of interest in analysis and modelling of biological systems especially on predator-prey systems with or without delay. The inclusion of delays in these has illustrated more complicated and richer dynamics in terms of stability, bifurcation, periodic solutions and so on [1-10].

In this paper, the basic model we consider is based on the following coupled delayed-differential equations

$ \left\{ \begin{array}{*{35}{l}} \dot{u}=a-u-4\frac{uv(t-\tau )}{1+{{u}^{2}}}, \\ \dot{v}=\sigma b(u-\frac{uv(t-\tau )}{1+{{u}^{2}}}), \\ \end{array} \right. $

(1)

where $a$, $b$ represent growth rate of the prey and predator, and $\tau$ denotes the delay time for the prey density, $u$ and $v$ can be interpreted as the densities of prey and predator prey populations at time $t$.

It is well know that the harvesting has a strong impact on the dynamics of a model. The aim is to determine how much we can harvest, and there are basically several types of harvesting reported in the our usual literature :

(ⅰ) Constant harvesting where a constant number of individuals are harvested per unit of time[12, 13].

(ⅱ) Proportional harvesting $h(x)=qEx$ that means the number of individuals harvested per unit of time is proportional to the current population.

It was noticed that the proportionate harvesting embodies several unrealistic features like random search for prey, equal likelihood of being captured for every prey species, unbounded linear increase of $h(x)$ with $x$ for fixed $ < $ $E$ and unbounded linear increase of $h(x)$ with $E$ for fixed $x$. These restrictive features are largely removed in the nonlinear harvesting $H(x, E)=\frac{qEx}{m_{1}E+m_{2}x}$ [14-16], where $q$ is the catchability coefficient, $E$ is the effort applied to harvest individuals which is measured in terms of number of vessels being used to harvest the individuals population and $m_{1}$, $m_{2}$ are suitable positive constants. The functional $H(x, E)=\frac{qEx}{m_{1}E+m_{2}x}$ is more realistic in the sense that the above unrealistic features are largely removed. It may be noted that $H(x, E)\rightarrow\frac{qE}{m_{2}}$ as $x\rightarrow\infty$ and $H(x, E)\rightarrow\frac{qx}{m_{1}}$ as $E\rightarrow\infty$. This shows that the nonlinear harvesting function exhibits saturation effects with respect to both the stock abundance and the effort-level. Also the parameter $m_{1}$ is proportional to the ratio of the stock-level to the harvesting rate(catch-rate)at higher levels of effort and $m_{2}$ is proportional to the ratio of the effort-level to the harvesting rate at higher stock-levels.

In order to utilize the harvest rate that leads to the largest possible value for the total discounted net revenue which depends on the population level, we assume joint harvesting of prey where we use a more realistic form of the catch-rate function by Clark [15], we consider the following system

$ \left\{ \begin{array}{*{35}{l}} \dot{u}=a-u-4\frac{uv(t-\tau )}{1+{{u}^{2}}}-\frac{qEu}{{{m}_{1}}E+{{m}_{2}}u}, \\ \dot{v}=\sigma b(u-\frac{uv(t-\tau )}{1+{{u}^{2}}}). \\ \end{array} \right. $

(2)

In daily life, economic profit is a very important factor for governments, merchants and even every citizen, so it is necessary to research biological systems, which can be described by differential-algebraic equations or differential-difference-algebraic equations. In 1954, Gorden[11]studied the effect of the harvest effort on ecosystem from an economic perspective and proposed the following economic principle:

$ \text{Net}\ \text{Economic}\ \text{Revenue}(NER)=\text{Total}\ \text{Revenue}(TR)\text{-Total}\ \text{Cost}(TC). $

Associated with the system $(2)$, an algebraic equation which consider the economic profit $m$ of the harvest effort on prey can be established as follows:

$\frac{qE}{{{m}_{1}}E+{{m}_{2}}u}(pu(t)-c)=m. $

And then, we obtain a predator-prey biological economic model which takes the form of

$ \left\{ \begin{array}{*{35}{l}} \dot{u}=a-u-4\frac{uv(t-\tau )}{1+{{u}^{2}}}-\frac{qEu}{{{m}_{1}}E+{{m}_{2}}u}, \\ \dot{v}=\sigma b(u-\frac{uv(t-\tau )}{1+{{u}^{2}}}), \\ 0=\frac{qE}{{{m}_{1}}E+{{m}_{2}}u}(pu(t)-c)-m. \\ \end{array} \right. $

(3)

For convenience, let

$ \begin{array}{l} f(X(t),E(t)) = \left( {\begin{array}{*{20}{c}} {{f_1}(X(t),\;\;E(t))}\\ {{f_2}(X(t),\;\;E(t))} \end{array}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \left( {\begin{array}{*{20}{c}} {a - u - 4\frac{{uv(t - \tau )}}{{1 + {u^2}}} - \frac{{qEu}}{{{m_1}E + {m_2}u}}}\\ {\sigma b(u - \frac{{uv(t - \tau )}}{{1 + {u^2}}})} \end{array}} \right),\\ g(X(t),E(t)) = \frac{{qE}}{{{m_1}E + {m_2}u}}(pu(t) - c) - m, \end{array} $

where$X(t)=(u(t), v(t))^{T}$, $\tau$ is a bifurcation parameter, which will be difined in what follows.

In this paper, we mainly discuss the effects of the economic profit on the dynamics of system $(3)$ in region

$ R_{+}^{3}=Y(t)=((u(t),v(t),E(t))|u(t)\ge 0,v(t)\ge 0,E(t)\ge 0). $

The organization of this paper is as follows:regarding $\tau$ bifurcation parameter, we study the stability of the equilibrium point of the system $(3)$ and Hopf bifurcation of the positive equilibrium depending on $\tau$ where we show that positive equilibrium loses its stability and the system $(3)$ exhibits Hopf bifurcation in the second section. Then, based on the new normal form of the differential-algebraic system introduced by Chen et al.[17] and the normal form approach theory and center manifold theory introduced by Hassard et al.[18], we derive the formula for determining the properties of Hopf bifurcation of the system in the third section. Numerical simulations aimed at justifying the theoretical analysis will be reported in section 4. Finally, this paper ends with a discussion.

For system $(3)$, we can see that there an equilibrium in $R^{3}_{+}$ if and equations

$ \left\{ \begin{array}{*{35}{l}} a-u-4\frac{uv(t-\tau )}{1+{{u}^{2}}}-\frac{qEu}{{{m}_{1}}E+{{m}_{2}}u}=0, \\ \sigma b(u-\frac{uv(t-\tau )}{1+{{u}^{2}}})=0, \\ \frac{qE}{{{m}_{1}}E+{{m}_{2}}u}(pu(t)-c)-m=0. \\ \end{array} \right. $

(4)

Through a simple calculation, we obtain

$ \begin{array}{*{35}{l}} {{X}_{0}}=(\frac{({{m}_{2}}a-5{{m}_{1}}{{E}_{0}}-q{{E}_{0}})+\sqrt{{{({{m}_{2}}a-5{{m}_{1}}{{E}_{0}}-q{{E}_{0}})}^{2}}+20{{m}_{1}}{{m}_{2}}a{{E}_{0}}}}{10{{m}_{2}}}, \\ \ \ \ \ \ \ \ \ 1+u_{0}^{2},\frac{m{{m}_{2}}{{u}_{0}}}{q(p{{u}_{0}}-c)-m{{m}_{1}}}). \\ \end{array} $

In this paper we only concentrate on the interior equilibrium of the system $(2)$, since the biological meaning of the interior equilibrium implies that the prey, the predator and the harvest effort on prey all exist, which are relevant to our study. Thus, we assume that $m_{2}a-5m_{1}E_{0}-qE_{0}>0$, $q(pu_{0}-c)-mm_{1}>0$. In order to analyze the local stability of the positive equilibrium point for the system $(3)$, we first use the linear transformation$X^{T}=QN^{T}$, where

$ N = {(x,y,\bar E)^T},Q = \left( {\begin{array}{*{20}{c}} {\;\;\;\;\;\;\;\;\;\;1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;0}\\ {\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1\;\;\;0}\\ { - \frac{{p{m_2}{E^2} + {m_2}cE}}{{{m_2}{u_0}(p{u_0} - c)}}\;\;\;0\;\;\;1} \end{array}} \right). $

Then we have $D_{X}g(X_{0})Q=(0, 0, \frac{m_{2}u_{0}q(pu_{0}-c)}{(m_{1}E+m_{2}u)^{2}})$, $x=u$, $y=v$, $\bar{E}=\frac{pm_{2}E^{2}_{0}+m_{2}cE_{0}}{m_{_{2}u_{0}(pu_{0}-c)}}x+E$, for which system $(2)$ is transformed into

$ \left\{ \begin{array}{*{35}{l}} \dot{x}=a-x-4\frac{xy(t-\tau )}{1+{{x}^{2}}}-\frac{q(\bar{E}-\frac{p{{m}_{2}}E_{0}^{2}+{{m}_{2}}c{{E}_{0}}}{{{m}_{2}}{{u}_{0}}(p{{u}_{0}}-c)}x)x}{{{m}_{1}}(\bar{E}-\frac{p{{m}_{2}}E_{0}^{2}+{{m}_{2}}c{{E}_{0}}}{{{m}_{2}}{{u}_{0}}(p{{u}_{0}}-c)}x)+{{m}_{2}}x}, \\ \dot{y}=\sigma b(x-\frac{xy(t-\tau )}{1+{{x}^{2}}}), \\ 0=\frac{q(\bar{E}-\frac{p{{m}_{2}}E_{0}^{2}+{{m}_{2}}c{{E}_{0}}}{{{m}_{2}}{{u}_{0}}(p{{u}_{0}}-c)}x)}{{{m}_{1}}(\bar{E}-\frac{p{{m}_{2}}E_{0}^{2}+{{m}_{2}}c{{E}_{0}}}{{{m}_{2}}{{u}_{0}}(p{{u}_{0}}-c)}x)+{{m}_{2}}x}(px(t)-c)-m. \\ \end{array} \right. $

(5)

Now we derive the formula for determining the properties of the positive equilibrium point of the system $(5)$. First we consider the local parametric $\psi$ the third equation of system $(4)$ as the literature, which defined as follows:

$ {{[x(t),y(t),\bar{E}(t)]}^{T}}=\psi (Z(t))=N_{0}^{T}+{{U}_{0}}Z(t)+{{V}_{0}}h(Z(t)),g(\psi (Z(t)))=0, $

where

$ \begin{array}{*{35}{l}} {{U}_{0}}=\left( \begin{array}{*{35}{l}} 1\ \ \ 0 \\ 0\ \ \ 1 \\ 0\ \ \ 0 \\ \end{array} \right),{{V}_{0}}=\left( \begin{array}{*{35}{l}} 1 \\ 0 \\ 0 \\ \end{array} \right), \\ Z\left( t \right)={{\left( {{y}_{\text{1}}}\left( t \right),{{y}_{\text{2}}}\left( t \right) \right)}^{T}}, \\ {{N}_{\text{0}}}=\left( {{u}_{\text{0}}},{{v}_{\text{0}}},{{{\bar{E}}}_{0}} \right),h\left( Z\left( t \right) \right)=h\left( {{y}_{1}}\left( t \right),{{y}_{2}}\left( t \right) \right), \\ \end{array} $

$R^{2} \rightarrow R$ is a smooth mapping. Then we can obtain the parametric system of the $(4)$ as follows:

$ \left\{ \begin{array}{*{35}{l}} \overset{.}{\mathop{{{y}_{1}}}}\,=-{{a}_{1}}{{y}_{1}}(t)-{{a}_{2}}{{y}_{2}}(t-\tau )+{{a}_{3}}y_{1}^{2}(t)+{{a}_{4}}{{y}_{1}}(t){{y}_{2}}(t-\tau ), \\ \overset{.}{\mathop{{{y}_{2}}}}\,={{b}_{1}}{{y}_{1}}(t)-{{b}_{2}}{{y}_{2}}(t-\tau )+{{b}_{3}}y_{1}^{2}(t)+{{b}_{4}}{{y}_{1}}(t){{y}_{2}}(t-\tau ), \\ \end{array} \right. $

(6)

where

$ \begin{array}{l} {a_1} = \frac{{4{y_0} - 4x_0^2{y_0}}}{{{{(1 + x_0^2)}^2}}} - \frac{{mc}}{{{{(p{x_0} - c)}^2}}} + 1,{a_2} = \frac{{4{x_0}}}{{1 + x_0^2}},\\ {a_3} = \frac{{12{x_0}{y_0} - 4x_0^3{y_0}}}{{{{(1 + x_0^2)}^3}}} - \frac{{mc}}{{{{(p{x_0} - c)}^3}}},{a_4} = - \frac{{4(1 - x_0^2)}}{{{{(1 + x_0^2)}^2}}},\\ {b_1} = \frac{{\sigma b(1 + x_0^2 + x_0^2{y_0} - {y_0})}}{{{{(1 + x_0^2)}^2}}},{b_2} = \frac{{\sigma b{x_0}}}{{1 + x_0^2}},\\ {b_3} = - \frac{{\sigma b(3{x_0}{y_0} - x_0^3{y_0})}}{{{{(1 + x_0^2)}^3}}},{b_4} = - \frac{{\sigma b((1 - x_0^2))}}{{{{(1 + x_0^2)}^2}}}, \end{array} $

so we can get the linearized system of parametric system $(6)$ as follows:

$ \left\{ \begin{array}{*{35}{l}} \overset{.}{\mathop{{{y}_{1}}}}\,=-{{a}_{1}}{{y}_{1}}(t)-{{a}_{2}}{{y}_{2}}(t-\tau ), \\ \overset{.}{\mathop{{{y}_{2}}}}\,={{b}_{1}}{{y}_{1}}(t)-{{b}_{2}}{{y}_{2}}(t-\tau ). \\ \end{array} \right. $

(7)

The associated characteristic equation of the system $(7)$ is

$ \text{det}\left( \begin{matrix} \lambda +{{a}_{1}}\ \ \ \ \ {{a}_{2}}{{e}^{-\lambda \tau }} \\ \ \ \ -{{b}_{1}}\ \ \ \ \lambda +{{b}_{2}}{{e}^{-\lambda \tau }} \\ \end{matrix} \right)=0. $

This characteristic equation determines the local stability of the equilibrium solution

$ {{\lambda }^{2}}+({{a}_{1}}+{{b}_{2}}{{e}^{-\lambda \tau }})\lambda +({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}}){{e}^{-\lambda \tau }}=0. $

(8)

Case 1  When there is no time delay, i.e. $\tau=0$ in Eq. $(8)$, it becomes

$ {{\lambda }^{2}}+({{a}_{1}}+{{b}_{2}})\lambda +({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})=0. $

The associate eigenvalues are $\lambda_{1, 2}=\frac{-(a_{1}+b_{2})\pm\sqrt{(a_{1}+b_{2})^{2}-4(a_{1}b_{2}+a_{2}b_{1})})}{2}$ so that one has the following Lemma.

Lemma 1  If $a_{1}+b_{2}>0$, then the equilibrium point of the system $(2)$ with $\tau=0$ is asymptotically stable.

Case 2  Suppose now that $\tau\neq0$ in Eq. $(8)$. We will investigate location of the roots of the transcendental equation. First, we examine when this equation has pure imaginary roots $\lambda=\pm i\omega$ with $\omega$ real number and $\omega>0$. This is given by the following lemma.

Lemma 2  The characteristic equation $(8)$ associated with Eq. $(8)$ has one pure imaginary root.

Proof  Let $\lambda=\pm i\omega$ be a root of characteristic equation $(8)$ where $\omega>0$, then we have

$ -{{\omega }^{2}}+i\omega ({{a}_{1}}+{{b}_{2}}(\cos \omega \tau -i\sin \omega \tau ))+({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})(\cos \omega \tau -i\sin \omega \tau )=0. $

Separating real and imaginary parts, we have the following two equation

$ ({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})\cos \omega \tau +{{b}_{2}}\omega \sin \omega \tau ={{\omega }^{2}}. $

(9)

$ {{b}_{2}}\omega \cos \omega \tau -({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})\sin \omega \tau =-{{a}_{1}}\omega . $

(10)

By taking square of both sides of $(9)$ and $(10)$ and then adding them up, one obtains the following equation

$ {{\omega }^{4}}+(a_{1}^{2}-b_{2}^{2}){{\omega }^{2}}-{{({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})}^{2}}=0. $

(11)

Solving now this for $\omega^{2}$ leads to

$ \omega =\sqrt{\frac{b_{2}^{2}-a_{1}^{2}+\sqrt{{{(a_{1}^{2}-b_{2}^{2})}^{2}}+4{{({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})}^{2}})}}{2}}. $

(12)

Which is a unique positive root of $(8)$. From $(9)$ and $(10)$, we also obtain a sequence of the critical values of $\tau$ defined by

$ {{\tau }_{k}}=\frac{1}{\omega }\{{{\cos }^{-1}}(\frac{{{\omega }^{2}}{{a}_{2}}{{b}_{1}}}{{{({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})}^{2}}+b_{2}^{2}{{\omega }^{2}}})\}+\frac{2k\pi }{\omega },(k=0,1,2,3,...), $

(13)

this completes the proof.Notice that it may be seen easily that the purely imaginary root $\omega$ is simple. Let $\lambda(\tau)=\alpha(\tau)+i\omega(\tau)$ denote the roots of eq. $ (8)$ near $\tau=\tau_{k}$ satisfying conditions $\alpha(\tau_{k})=0$ and $\omega(\tau_{k})=\omega$, then we have the following transversality condition.

Lemma 3  The following transversality condition

$ \frac{d\text{Re}\{\lambda ({{\tau }_{k}})\}}{d\tau }>0\ (k=0,1,2,3,...) $

hold.

Proof  Differentiating eq. $(8)$ with respect to $\tau$, we get

$ \frac{d\lambda }{d\tau }=\frac{({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})\lambda +{{b}_{2}}{{\lambda }^{2}}{{e}^{-\lambda \tau }}}{2\lambda +{{a}_{1}}+{{b}_{2}}{{e}^{-\lambda \tau }}-{{b}_{2}}\lambda \tau {{e}^{-\lambda \tau }}-({{a}_{1}}{{b}_{2}}+{{a}_{2}}{{b}_{1}})\tau {{e}^{-\lambda \tau }}}. $

(14)

First substituting $\lambda=i\omega$ into it and then flipping it over and finally taking its real part, one obtains

$ \begin{array}{*{20}{l}} {{\rm{sign}}\{ {\rm{Re}}(\frac{{d\lambda }}{{d\tau }})\} {|_{\lambda = i\omega }} = {\rm{sign}}\{ {\rm{Re}}{{(\frac{{d\lambda }}{{d\tau }})}^{ - 1}}\} {|_{\lambda = i\omega }}}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \frac{{{{({a_1}{b_2} + {a_2}{b_1})}^2}(2{w^2} + a_1^2 + b_2^2) + ({{({a_1}{b_2})}^2} + 2b_2^2{\omega ^2} + b_2^4){\omega ^2}}}{{{{({\omega ^2}b_2^2 + {{({a_1}{b_2} + {a_2}{b_1})}^2})}^2}}}.} \end{array} $

(15)

Therefore ${\rm{sign}}\{ Re(\frac{{d\lambda }}{{d\tau }})\} {|_{\lambda = i\omega }} >0$. This completes the proof. Summarizing the above remarks and combining Lemmas, we have the following results on the distribution of roots of eq.$ (8)$.

Theorem 1  For the system $(3)$, the following statements are true

(ⅰ) The equilibrium point $(u_{0}, v_{0}, E_{0})$ is asymptotically stable for $\tau=0$ if $a+b>0$.

(ⅱ) $(u_{0}, v_{0}, E_{0})$ is asymptotically stable for $\tau < \tau_{0}$ and unstable $\tau>\tau_{0}$, where $\tau_{0}=\frac{1}{\omega}\{\cos^{-1}(\frac{\omega^{2}a_{2}b_{1}}{(a_{1}b_{2}+a_{2}b_{1})^{2}+b_{2}^{2}\omega^{2}})\}$. Furthermore, the system $(2)$ undergoes a Hopf bifurcation at $(u_{0}, v_{0}, E_{0})$ when $\tau=\tau_{0}$.

In this section, we investigate the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions based on the new normal form of the differential-algebraic system introduced by Chen et al. [17] and the normal form approach theory and center manifold theory introduced by Hassard et al.[18].

In the following part, we assume that the system $(3)$ undergoes Hopf bifurcations at the positive equilibrium $Y_{0}$ for $\tau=\tau_{k}$, that is, the system $(5)$ undergoes Hopf bifurcations at the positive equilibrium $N_{0}$ for $\tau=\tau_{k}$, and we let $i\omega$ is the corresponding purely imaginary root of the characteristic equation at the positive equilibrium $N_{0}$. In order to investigate the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions system $(3)$, we consider the parametric system $(6)$ of the system $(5)$. First by the transformation $\bar{y_{1}}=y_{1}$, $\bar{y_{2}}=y_{2}$, $t=\frac{t}{\tau}$, $\tau=\tau_{k}+\mu$, $\bar{Z}=(\bar{y_{1}}, \bar{y_{2}})$, for simplicity, we continue to use $Z$ said $\bar{Z}$, then the parametric system $(5)$ of the system $(4)$ is equivalent to the following Functional Differential Equation system in $C=C([-1,0], \Re^{2}) $,

$ \bar{Z}(t)={{L}_{\mu }}(Z(t))+f(\mu ,Z(t)), $

(16)

Where $\bar{Z}=(\bar{y_{1}}, \bar{y_{2}})^{T}$, and $L_{\mu}: C\to R$, $f : R\to R$ are given, respectively, by

$ {L_\mu }(\phi ) = ({\tau _k} + \mu )\left( {\begin{array}{*{20}{c}} { - {a_1}\;\;0}\\ {\;\;{b_1}\;\;\;0} \end{array}} \right){\phi ^T}(0) + ({\tau _k} + \mu )\left( {\begin{array}{*{20}{c}} {0\;\; - {a_2}}\\ {0\;\; - {b_2}} \end{array}} \right){\phi ^T}( - 1) $

and

$ f(\mu ,\phi )=({{\tau }_{k}}+\mu )\left( \begin{matrix} {{f}_{11}} \\ {{f}_{22}} \\ \end{matrix} \right), $

where

$ \begin{align} {{f}_{11}}={{a}_{3}}\phi _{1}^{2}(0)+{{a}_{4}}{{\phi }_{1}}(0){{\phi }_{2}}(-1), \\ {{f}_{22}}={{b}_{3}}\phi _{1}^{2}(0)+{{b}_{4}}{{\phi }_{1}}(0){{\phi }_{2}}(-1), \\ \end{align} $

and $\phi=(\phi_{1}, \phi_{2})$. By the Riesz representation theorem, there exists a matrix function whose components are bounded Variation function $\eta(\theta, \mu)$ in $\theta\in[-1,0]$ such that

$ {{L}_{\mu }}\phi =\int_{-1}^{0}{d}\eta (\theta ,\mu )\phi (\theta ),\phi \in C. $

In fact, we can choose

$ \eta (\theta ,\mu )=({{\tau }_{k}}+\mu )\left( \begin{matrix} -{{a}_{1}}\ \ 0 \\ \ \ {{b}_{1}}\ \ \ 0 \\ \end{matrix} \right)\delta (\theta )+({{\tau }_{k}}+\mu )\left( \begin{matrix} 0\ \ -{{a}_{2}} \\ 0\ \ -{{b}_{2}} \\ \end{matrix} \right)\delta (\theta +1), $

where

$ \delta (\theta )=\left\{ \begin{array}{*{35}{l}} 0,\quad \ \ \theta \ne 0, \\ 1,\quad \ \ \theta =0, \\ \end{array} \right. $

For $\phi^{1}([-1,0], \Re^{2})$, define

$ A(\mu )\phi (\theta )=\left\{ \begin{array}{*{35}{l}} \frac{d\phi (\theta )}{d\theta },\quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -1\le \theta < 0, \\ \int_{-1}^{0}{d}\eta (\theta ,\mu )\phi (\theta ),\quad \ \ \theta =0. \\ \end{array} \right. $

Then system $(15)$ is equivalent to

$ \dot{Z}(t)=A(\mu ){{Z}_{t}}+R(\mu ){{Z}_{t}}. $

(17)

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

$ {{A}^{*}}\psi (s)\left\{ \begin{align} -\frac{d\psi (s)}{ds},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\le \text{s} < 1, \\ \int_{-1}^{0}{d}{{\eta }^{T}}(s,0)\phi (-s),\ \ \ \ \ \ \ \ \ s=0, \\ \end{align} \right. $

(18)

and a linear inner product is given by

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

(19)

where $\eta(\theta)=\eta(\theta, 0)$. it is easy to verify that $A(0)$ and ${{A}^{*}}\left( 0 \right)$ are a pair of adjoint operators.

From the discussions in Section $2$, we know that $\pm i\omega$ are eigenvalues of $A(0)$. Thus, they are also eigenvalues of ${{A}^{*}}\left( 0 \right)$, Next we calculate the eigenvector $q(\theta)$ of $A$ belonging to $i\omega$ and eigenvector $q^{*}(s)$ of ${{A}^{*}}\left( 0 \right)$ belonging to $-i\omega$. Then it is not difficult to show that

$ q(\theta )={{(1,\beta )}^{T}}{{e}^{i\omega {{\tau }_{k}}\theta }},{{q}^{*}}(s)=G({{\beta }^{*}},1){{e}^{i\omega {{\tau }_{k}}s}}, $

where

$ \begin{array}{l} \beta = - \frac{{{a_1} + i\omega }}{{{a_2}{e^{ - i\omega }}}},{\beta ^*} = \frac{{i\omega - {b_2}{e^{ - i\omega }}}}{{{a_2}{e^{ - i\omega }}}},\\ \bar G = {\{ \overline {{\beta ^*}} + \beta + \beta {\tau _k}({a_2}{\beta ^*} + {b_2}){e^{ - i\omega {\tau _k}}}\} ^{ - 1}}. \end{array} $

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

Next, we study the stability of bifurcated periodic solutions. Using the same notations as in Hassard et al.[18]. We first compute the coordinates to describe the centre manifold $C_{0}$ at $\mu=0$. Define

$ \dot{z}(t)=\langle {{q}^{*}},{{Z}_{t}}\rangle ,W(t,\theta )={{Z}_{t}}-2\text{Re}\{z(t)q(\theta )\}. $

(20)

On the center manifold $C_{0}$, we have

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

(21)

In fact, $z$ and $\bar{z}$ are local coordinates for center manifold $C_{0}$ in the direction of $q$ and $q^{*}$. Note that $W$ is real if $Z_{t}$ is real. We consider only real solutions. For the solution $Z_{t}\in C_{0}$, since $\mu=0$ and Eq.$(15)$, we have

$ \begin{array}{*{35}{l}} \dot{z}=i\omega {{\tau }_{k}}z+\langle {{q}^{*}}(\theta ),f(0,W(z,\bar{z},\theta )+2Re[z(t)q(\theta )])\rangle \\ \ \ =i\omega {{\tau }_{k}}z+{{{\bar{q}}}^{*}}(0)f(0,W(z,\bar{z},0)+2Re[z(t)q(\theta )]), \\ \end{array} $

(22)

rewrite it as

$ \dot{z}=i\omega {{\tau }_{k}}z+g(z,\bar{z}), $

(23)

where

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

(24)

From $(16)$ and $(22)$, we have

$ \dot{W}=\overset{.}{\mathop{{{Z}_{t}}}}\,-\dot{z}-\dot{\bar{z}}\bar{q}=\left\{ \begin{align} AW-2\text{Re}\{{{{\bar{q}}}^{*}}(0)f(z,\bar{z})\}q(\theta ),\ \ \ \ \ \ \ \ \ -1\le \theta < 0, \\ AW-2\text{Re}\{{{{\bar{q}}}^{*}}(0)f(z,\bar{z})\}q(\theta )+f,\ \ \ \theta =0. \\ \end{align} \right. $

(25)

Rewrite $(25)$ as

$ \dot{W}=AW+H(z,\bar{z},\theta ), $

(26)

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 . $

(27)

Substituting the corresponding series into $(25)$ and comparing the coefficient, we obtain

$ \left( A-2i{{\omega }_{0}}{{\tau }_{0}} \right){{W}_{20}}=-{{H}_{20}}(\theta ),\ \ A{{W}_{11}}(\theta )=-{{H}_{11}}(\theta ). $

(28)

Notice that $q(\theta)=(1, \beta)^{T}e^{i\omega\tau_{k}\theta}, q^{*}(0)=G(\beta^{*}, 1)$, and $(20)$ we obtain

$ \begin{array}{*{20}{l}} {{y_{1t}}(0)\;\;\; = \;\;z + \bar z + W_{20}^{(1)}(0)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(1)}(0)(\theta )z\bar z + W_{02}^{(1)}(0)(\theta )\frac{{{{\bar z}^2}}}{2} + o(\mid z,\bar z{\mid ^3}),}\\ {{y_{2t}}(0)\;\;\; = \;\;\beta z + \overline {\beta z} + W_{20}^{(2)}(0)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(2)}(0)(\theta )z\bar z + W_{02}^{(2)}(0)(\theta )\frac{{{{\bar z}^2}}}{2} + o(\mid z,\bar z{\mid ^3}),}\\ {{y_{1t}}( - 1) = z{e^{ - i\omega \theta }} + \bar z{e^{i\omega \theta }} + W_{20}^{(1)}( - 1)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(1)}( - 1)(\theta )z\bar z + W_{02}^{(1)}( - 1)(\theta )\frac{{{{\bar z}^2}}}{2}}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\; + o(\mid z,\bar z{\mid ^3}),}\\ {{y_{2t}}( - 1) = \beta z{e^{ - i\omega \theta }} + \overline {\beta z} {e^{i\omega \theta }} + W_{20}^{(2)}( - 1)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(2)}( - 1)(\theta )z\bar z + W_{02}^{(2)}( - 1)(\theta )\frac{{{{\bar z}^2}}}{2}}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\; + o(\mid z,\bar z{\mid ^3}).} \end{array} $

According to $(22)$ and $(23)$, we know

$ g(z,\bar{z})={{\bar{q}}^{*}}(0){{f}_{0}}(z,\bar{z})=\bar{G}{{\tau }_{k}}(\overline{{{\beta }^{*}}},1)\left( \begin{matrix} {{f}_{11}} \\ {{f}_{22}} \\ \end{matrix} \right), $

(29)

where

$ \begin{array}{*{35}{l}} f_{11}^{0}={{a}_{3}}y_{1t}^{2}(0)+{{a}_{4}}{{y}_{1t}}(0){{y}_{2t}}(-1), \\ f_{22}^{0}={{b}_{3}}y_{1t}^{2}(0)+{{b}_{4}}{{y}_{1t}}(0){{y}_{2t}}(-1). \\ \end{array} $

By $(21)$, it follows that

$ \begin{array}{l} g(z, \bar z) = \bar G{\tau _k}\{ {a_3}\overline {{\beta ^*}} {(z + \bar z + W_{20}^{(1)}(0)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(1)}(0)(\theta )z\bar z + W_{02}^{(1)}(0)(\theta )\frac{{{{\bar z}^2}}}{2} + o(\mid z, \bar z{\mid ^3}))^2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + {a_4}\overline {{\beta ^*}} (z + \bar z + W_{20}^{(1)}(0)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(1)}(0)(\theta )z\bar z + W_{02}^{(1)}(0)(\theta )\frac{{{{\bar z}^2}}}{2} + o(\mid z, \bar z{\mid ^3}))\\ \;\;\;\;\;\;\;\;\;\;\;\;\; \times (\beta z{e^{ - i\omega \theta }} + \overline {\beta z} {e^{i\omega \theta }} + W_{20}^{(2)}( - 1)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(2)}( - 1)(\theta )z\bar z + W_{02}^{(2)}( - 1)(\theta )\frac{{{{\bar z}^2}}}{2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + o(\mid z, \bar z{\mid ^3})) + {b_3}(z + \bar z + W_{20}^{(1)}(0)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(1)}(0)(\theta )z\bar z + W_{02}^{(1)}(0)(\theta )\frac{{{{\bar z}^2}}}{2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + o(\mid z, \bar z{\mid ^3}){)^2} + {b_4}\overline {{\beta ^*}} (z + \bar z + W_{20}^{(1)}(0)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(1)}(0)(\theta )z\bar z + W_{02}^{(1)}(0)(\theta )\frac{{{{\bar z}^2}}}{2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + o(\mid z, \bar z{\mid ^3})) \times (\beta z{e^{ - i\omega \theta }} + \overline {\beta z} {e^{i\omega \theta }} + W_{20}^{(2)}( - 1)(\theta )\frac{{{z^2}}}{2} + W_{11}^{(2)}( - 1)(\theta )z\bar z\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + W_{02}^{(2)}( - 1)(\theta )\frac{{{{\bar z}^2}}}{2} + o(\mid z, \bar z{\mid ^3}))\} . \end{array} $

That is

$ \begin{array}{l} g(z, \bar z) = \bar G{\tau _k}\{ {z^2}[{a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} \beta {e^{-i\omega \theta }} + {b_3} + {b_4}\beta {e^{-i\omega \theta }}]\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ + }}z\bar z[2{a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} (\beta {e^{-i\omega \theta }} + \bar \beta {e^{i\omega \theta }}) + 2{b_3} + {b_4}(\beta {e^{-i\omega \theta }} + \bar \beta {e^{i\omega \theta }})]\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + {{\bar z}^2}[{a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} \beta {e^{i\omega \theta }} + {b_3} + {b_4}\beta {e^{i\omega \theta }}]\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + {z^2}\bar z[({a_3}\overline {{\beta ^*}} + \frac{1}{2}{a_4}\overline {{\beta ^*}} \beta {e^{i\omega \theta }} + {b_3} + \frac{1}{2}{b_4}\bar \beta {e^{i\omega \theta }})W_{20}^{(1)}(0)]\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + ({a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} \beta {e^{ - i\omega \theta }} + 2{b_3} + {b_4}\beta {e^{ - i\omega \theta }})W_{11}^{(1)}(0)\\ \;\;\;\;\;\;\;\;\;\;\;\;\; + (\frac{1}{2}{a_4}\overline {{\beta ^*}} + {b_4})W_{20}^{(2)}( - 1) + ({a_4}\overline {{\beta ^*}} + {b_4})W_{11}^{(2)}( - 1)\} . \end{array} $

By comparing the coefficients with $(23)$, it follows that

$ \begin{array}{l} {g_{20}} = 2\bar G{\tau _k}[{a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} \beta {e^{-i\omega \theta }} + {b_3} + {b_4}\beta {e^{-i\omega \theta }}], \\ {g_{11}} = \bar G{\tau _k}[2{a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} (\beta {e^{-i\omega \theta }} + \bar \beta {e^{i\omega \theta }}) + 2{b_3} + {b_4}(\beta {e^{-i\omega \theta }} + \bar \beta {e^{i\omega \theta }})], \\ {g_{02}} = 2\bar G{\tau _k}[{a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} \beta {e^{i\omega \theta }} + {b_3} + {b_4}\beta {e^{i\omega \theta }}], \\ {g_{21}} = 2\bar G{\tau _k}\{ ({a_3}\overline {{\beta ^*}} + \frac{1}{2}{a_4}\overline {{\beta ^*}} \beta {e^{i\omega \theta }} + {b_3} + \frac{1}{2}{b_4}\bar \beta {e^{i\omega \theta }})W_{20}^{(1)}(0)\\ \;\;\;\;\;\;\; + ({a_3}\overline {{\beta ^*}} + {a_4}\overline {{\beta ^*}} \beta {e^{ - i\omega \theta }} + 2{b_3} + {b_4}\beta {e^{ - i\omega \theta }})W_{11}^{(1)}(0)\\ \;\;\;\;\;\;\; + (\frac{1}{2}{a_4}\overline {{\beta ^*}} + {b_4})W_{20}^{(2)}( - 1) + ({a_4}\overline {{\beta ^*}} + {b_4})W_{11}^{(2)}( - 1)\} . \end{array} $

Since $W_{20}(\theta)$ and $W_{11}(\theta)$ appear in $g_{21}$, we still need to compute them.

From $(16)$ and $(24)$, 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 ). $

(30)

Comparing the coefficients of $(24)$ with $(25)$ gives that

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

(31)

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

(32)

From $(27)$ and the definition of A, we get

$ {{\dot{W}}_{20}}(\theta )=2i{{\omega }_{0}}{{\tau }_{0}}{{W}_{20}}(\theta )+{{g}_{20}}q(\theta )+{{\bar{g}}_{02}}\bar{q}(\theta ). $

Noting that $q(\theta)=q(0)e^{i\omega_{0}\tau_{0}\theta}$, we have

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

(33)

Similarly, from the definition of A, we have

$ \begin{array}{l} {{\dot W}_{11}}(\theta ) = {g_{11}}(\theta )q(\theta ) + {{\bar g}_{11}}(\theta )\bar q(\theta ), \\ {W_{11}}(\theta ) = - \frac{{i{g_{11}}}}{{{\omega _0}{\tau _0}}}q(0){e^{i{\omega _0}{\tau _0}}} + \frac{{i\overline {{g_{11}}} }}{{{\omega _0}{\tau _0}}}\bar q(0){e^{ - i{\omega _0}{\tau _0}}} + {M_2}. \end{array} $

(34)

In what follows we shall seek appropriate $M_{1}$ and $M_{2}$ in $(32)$ and $(33)$.From $(25)$ and $(29)$, we have

$ H_{20}(\theta)=-g_{20}q(\theta)-\bar{g}_{02}\bar{q}(\theta)+2\tau_{k}A_{1}, $

(35)

$ H_{11}(\theta)=-g_{11}q(\theta)-\bar{g}_{11}\bar{q}(\theta)2\tau_{k}A_{2}, $

(36)

Where

$ \begin{array}{l} {A_1} = \left( {\begin{array}{*{20}{c}} {A_1^{(1)}}\\ {A_1^{(2)}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{a_3} + {a_4}\beta {e^{ - i\omega \theta }}}\\ {{b_3} + {b_4}\beta {e^{ - i\omega \theta }}} \end{array}} \right),\\ {A_2} = \left( {\begin{array}{*{20}{c}} {A_2^{(1)}}\\ {A_2^{(2)}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {2{a_3} + {a_4}(\beta {e^{ - i\omega \theta }} + \bar \beta {e^{i\omega \theta }})}\\ {2{b_3} + {b_4}(\beta {e^{ - i\omega \theta }} + \bar \beta {e^{i\omega \theta }})} \end{array}} \right). \end{array} $

Substituting $(32)-(35)$ into $(27)$ and noticing

$ \begin{array}{l} (i\omega {\tau _k}I - \int_{ - 1}^0 {{e^{i\omega {\tau _k}\theta }}} d\eta (\theta ))q(0) = 0,\\ ( - i\omega {\tau _k}I - \int_{ - 1}^0 {{e^{ - i\omega {\tau _k}\theta }}} d\eta (\theta ))q(0) = 0, \end{array} $

We obtain

$ \left( \begin{matrix} {{a}_{1}}+2i\omega \ \ \ \ \ \ \ \ \ \ \ {{a}_{2}}{{e}^{-2i\omega {{\tau }_{k}}}} \\ \ \ -{{b}_{1}}\ \ \ \ \ \ 2i\omega +{{b}_{2}}{{e}^{-2i\omega {{\tau }_{k}}}} \\ \end{matrix} \right){{M}_{1}}=2{{A}_{1}}, $

(37)

$ \left( \begin{matrix} {{a}_{1}}\ \ \ {{a}_{2}} \\ -{{b}_{1}}\ \ \ {{b}_{2}} \\ \end{matrix} \right){{M}_{2}}=2{{A}_{2}}. $

(38)

It is easy to obtain $M_{1}$ and $M_{1}$ from $(36)$ and $(37)$, that is

$ \begin{array}{l} M_1^{(1)} = \frac{{(4i\omega + 2{b_2}{e^{ - 2i\omega {\tau _k}}})A_1^{(1)} + 2{a_2}{e^{ - 2i\omega {\tau _k}}}A_1^{(2)}}}{{2i\omega ({a_1} + {b_2}{e^{ - 2i\omega {\tau _k}}}) + ({a_1}{b_2} + {a_2}{b_1}){e^{ - 2i\omega {\tau _k}}} - 4{\omega ^2}}}, \\ M_1^{(2)} = \frac{{2{b_1}A_1^{(1)} + 2({a_1} + 2i\omega )A_1^{(2)}}}{{2i\omega ({a_1} + {b_2}{e^{ - 2i\omega {\tau _k}}}) + ({a_1}{b_2} + {a_2}{b_1}){e^{ - 2i\omega {\tau _k}}} - 4{\omega ^2}}}, \\ M_2^{(1)} = \frac{{2{b_2}A_2^{(1)} - 2{a_2}A_2^{(2)}}}{{{a_1}{b_2} + {a_2}{b_1}}}, M_2^{(2)} = \frac{{2{b_1}A_2^{(1)} + 2{a_1}A_2^{(2)}}}{{{a_1}{b_2} + {a_2}{b_1}}}. \end{array} $

Therefore, we can compute the following values

$ \begin{array}{l} {C_1}(0) = \frac{i}{{2\omega {\tau _k}}}({g_{11}}{g_{20}} - 2|{g_{11}}{|^2} - \frac{1}{3}|{g_{02}}{|^2}) + \frac{{{g_{21}}}}{2}, \\ {\mu _2} = - \frac{{Re\{ {C_1}(0)\} }}{{Re\{ \lambda '({\tau _{{\tau _k}}})\} }}, {\beta _2} = 2Re\{ {C_1}(0)\}, \\ {T_2} = - \frac{{Im\{ {C_1}(0)\} + {\mu _2}Im\{ \lambda '({\tau _{{\tau _k}}})\} }}{{\omega {\tau _{{\tau _k}}}}}. \end{array} $

Theorem 2   $\mu_{2}$ determines the direction of Hopf bifurcation: If $\mu_{2}>0$, then the Hopf bifurcation is supercritical and bifurcating periodic solutions exist for $\tau>\tau_{0}$; and if $\mu_{2} <0$, then the Hopf bifurcation is subcritical and bifurcating periodic solutions exist for $\tau <\tau_{0}$. $\beta_{2}$ determines the stability of the bifurcating periodic solutions: Bifurcating periodic are stable if $\beta_{2} <0$; unstable if $\beta_{2}>0$. $T_{2}$ determines the period of the bifurcating solution: The period increase if $T_{2}>0$, decreases if $T_{2} <0$.

In this section, we present some numerical results of system $(4)$ at different values of $\tau$. From section $(3)$, we have determined the direction of a hopf bifurcation and the stability of the bifurcating periodic solutions. We consider the following system:

$ \left\{ \begin{array}{l}\dot{u}=4-u-4\frac{uv(t-\tau)}{1+u^{2}}-\frac{Eu}{0.5E+1.5u}, \\ \dot{v}=16(u-\frac{uv(t-\tau)}{1+u^{2}}), \\ 0=\frac{E}{0.5E+1.5u}(4u(t)-2)-0.01, \end{array} \right. $

(39)

which has an only positive equilibrium $X_{0}=(0.7987, 1.6379, 0.0101)$. By algorithms in section 2, we obtain $\tau_{0}=0.143$, $w=6.8630$. So by Theorem 1, the equilibrium point $E^{*}$ is asymptotically stable when $\tau\in[0, \tau_{0})=[0, 0.143)$ and unstable when $\tau>0.143$ and also Hopf bifurcation occurs at $\tau=\tau_{0}=0.143$ as it illustrated by computer simulations.

Now we determine the direction of a Hopf bifurcation with $\tau_{0}=0.143$ and the other properties of bifurcating periodic solutions based on the theory of Hassard et al.[18], as it is discussed before. By means of software Matlab7.0, we can obtain the following values $c_{1}(0)=-0.7922-0.4937i$, $\lambda'(\tau_{0})=0.0137 + 0.0011i$, it follows that $\mu_{2}=0.0109>0$, $\beta_{2}=-1.5844<0$, $T_{2}=0.4958>0$, from which and Theorem 2 we can conclude that the Hopf bifurcation of system $(3)$ occurring at $\tau_{0}=0.143$ is supercritical and the bifurcating periodic solution exist for $\tau>\tau_{0}$ and the bifurcating periodic solution is stable. So by Theorem 2, the positive equilibrium point $X_{0}$ of the system $(3)$ is locally asymptotically stable when $\tau=0.141<\tau_{0}$ as is illustrated by computer simulations in Fig. 1. And periodic solutions occur from $X_{0}$ when $\tau=0.145>\tau_{0}$ as is illustrated by computer simulation in Fig. 2. Here, we choose the initial conditions $x(0)=0.9, y(0)=1.5, E(0)=0.01$ in our simulations.

Nowadays, biological resources in the predator-prey system are mostly harvested and sold with the purpose of achieving the economic profit, and economic profit is a very important factor for governments, merchants and even every citizen, so it is necessary to research biological economic systems, which motivates the introduction of harvesting in the predator-prey system, in this paper, we provided a new and efficacious method for the qualitative analysis of the Hopf bifurcation of a differential-algebraic biological economic system with time delay, via numerical simulations we can conclude that the stability properties of the system could switch with parameter τ that is incorporated on the time delay on prey density in the differential-algebraic biological economic system.Form an economic perspective, the persistence and sustainable development of the predator-prey ecosystem will be very important, so with the purpose of maintaining the sustainable development of the biological resources in practice and application, more precise mathematic or physical architectures of the differential-algebraic biological economic system may be proposed, demonstrating a mature strategy rather than a concept, and dynamic property of differential-algebraic economic system should be analysed in practice or from experimental point of views in future works.