In recent years, the increasingly serious problem of environmental degradation and resource shortage, made the analysis and modeling of biological systems more terested. The predator-prey system played a crucial role among the relationships between the biological population, and it naturally attracted much attention both for mathematicians and biologists, especially on predator-prey systems with or without time delay. As we know, delay differential equation models exhibit much more complicated dynamics than differential equation models without delay, see [1-12]. A lot of researchers studied the dynamics of predator-prey models with harvesting and obta ed many dynamic behaviors, such as stability of equilibrium, Hopf bifurcation, periodic solution, Bogdanov-Takens bifurcation, Neimark-Sacker bifurcation, and so on, see [10-15].

In [16], Lucas studied the dynamic properties of the following Leslie-Gower predator-prey system

$\begin{equation} \left\{\begin{array}{l} \dot{x}=x(a-y), \\ \dot{y}=y(d-k\frac{y}{x}), \end{array}\right. \end{equation}$

(1.1)

where $x$ and $y$ denote prey and predator population densities at time $t$, respectively, $a$, $d$, and $k$ are positive constants that represent the prey intrinsic growth rate, predator mortality rate, and the maximum value of the per capita reduction rate of $x$ due to $y$, respectively.

At present, economic profit is a very important factor for merchants, governments and even every citizen, so it is necessary to research biological economic systems, which are often described by differential-algebraic equations or differential difference-algebraic equations.

In 1954, Gordon [13] studied the effect of the harvest effort on ecosystem from an economic perspective and proposed the following economic theory:

Net Economic Revenue (NER) = Total Revenue (TR) -Total Cost (TC).

This provides theoretical evidence for the establishment of differential-algebraic equation.

Based on the economic theory as mentioned above and system (1.1), Liu and Fu [12] considered the following Leslie-Gower predator-prey system

$\begin{equation} \left\{\begin{array}{l} \dot{x}=x(a-y-E), \\ \dot{y}=y(d-k\frac{y}{x}), \\ 0=E(px-c)-v. \end{array}\right. \end{equation}$

(1.2)

They investigated the Hopf bifurcation of the above system without considering the effect of time delay and the harvesting of predator.

As is known to all, delay differential equation models exhibit much more complicated dynamics than ord ary differential equation models, see [1-12], as was po ted by Kuang [17] that any model of species dynamics without delays is an approximation at best. When we considered the model with non-selective harvesting, namely at the same time there are also the harvesting of predator and harvesting of the prey the model, it will be more l e with the actual situation of the predator-prey systems.

Motivated by the above discussion, this paper, by choosing the time delay as a bifurcation parameter and consider the predator-prey systems with non-selective harvesting, we vestigate a modified Leslie-Gower predator-prey systems with non-selective harvesting and time delay described by the following system

$\begin{equation} \left\{\begin{array}{l} \frac{dx}{dt}=x(t)(r_{1}-by(t-\tau)-E(t)), \\ \frac{dy}{dt}=y(t)(r_{2}-k\frac{y(t)}{x(t)}-E(t)), \\ 0=E(t)(p_{1}x(t)+p_{2}y(t)-c_{1}-c_{2})-m, \end{array}\right. \end{equation}$

(1.3)

where $p_{1}>0$ and $p_{2}>0$ are harvesting reward per unit harvesting effort for unit prey and predator, respectively; $c_{1}$ and $c_{2}$ are harvesting cost per unit harvesting effort for prey and predator, respectively; $m$ is the economic profit per unit harvesting effort.

In this paper, we ma ly discuss the effects of economic profit on the dynamics of system (1.3) the region $R_{+}^{3}=\{(x, y, E)|x>0, y>0, E>0)\}$.

For convenience, we let

$\begin{eqnarray*} & & f(X(t), E(t))=\left(\begin{array}{l} f_{1}(X(t), E(t))\\ f_{2}(X(t), E(t))\\ \end{array} \right) =\left( \begin{array}{l} x(t)(r_{1}-by(t-\tau)-E(t))\\ y(t)(r_{2}-k\frac{y(t)}{x(t)}-E(t))\\ \end{array} \right), \\ & & g(X(t), E(t)) =E(t)(p_{1}x(t)+p_{2}y(t)-c_{1}-c_{2})-m, \end{eqnarray*}$

where $X_{t}=(x, y)^{T}$.

The rest of the paper is arranged as follows: Section 2, the local stability of the positive equilibrium po ts are vestigated by the corresponding characteristic equation of system (1.3). In Section 3, by using the normal form and Hopf bifurcation theorem, we study the Hopf bifurcation of the nonnegative equilibrium depending on the parameter where we show that the positive equilibrium loses its stability and system (1.3) exhibits Hopf bifurcation the second section. In Section 4, the theoretical result is supplied by a numerical example. F ally, this paper ends with a brief discussion.

In this section, we discuss the local stability of a positive equilibrium for system (1.3). Now, we try to f d all possible positive equilibrium po ts of system (1.3). A po t $Y_{0}=(x_{0}, y_{0}, E_{0})$ is an equilibrium po t of system (1.3) if and only if $Y_{0}$ satisfy the following equations

$\begin{equation} \left\{\begin{array}{l} x(t)(r_{1}-by(t-\tau)-E(t))=0, \\ y(t)(r_{2}-k\frac{y(t)}{x(t)}-E(t))=0, \\ E(t)(p_{1}x(t)+p_{2}y(t)-c_{1}-c_{2})-m=0. \end{array}\right. \end{equation}$

(2.1)

From (2.1), we can easy get $E_{0}$ satisfy

$\begin{equation} E^{3}+\gamma_{1}E_{2}+\gamma_{2}E+\gamma_{3}=0, \end{equation}$

(2.2)

where

$\begin{eqnarray*} & & \gamma_{1}=\frac{bc_{1}+bc_{2}-kp_{1}-p_{2}r_{1}-p_{2}r_{2}}{p_{2}}, \\ & & \gamma_{2}=\frac{kp_{1}r_{1}+r_{1}p_{2}r_{2}+bm-bc_{1}c_{2}-bc_{2}r_{2}}{p_{2}}, ~~ \gamma_{3}=\frac{-bmr_{2}}{p_{2}}<0.\end{eqnarray*}$

Based on the root and coefficient relationship of equation and $\gamma_{3} < 0$, we can f d at least one positive root $E_{0}$, so system (1.3) has at least one positive equilibrium po t $Y_{0}=(x_{0}, y_{0}, E_{0})=(\frac{k(r_{1}-E_{0})}{b(r_{2}-E_{0})}, \frac{r_{1}-E_{0}}{b}, E_{0})$, where $r_{1}>E_{0}$, $r_{2}>E_{0}.$

Now, we derive the formula for determining the properties of the positive equilibrium point of system (1.3). As in [13], first we consider the local parametric $\psi$ of the third equation of system (1.3), which is defined as follows

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

where

$\begin{eqnarray*} & & U_{0}= \left( \begin{array}{lcr} 1&0 \\ 0&1 \\ 0&0 \\ \end{array} \right), ~~V_{0}= \left( \begin{array}{lcr} 0 \\ 0 \\ 1 \\ \end{array} \right), \\ & & Z(t)=(y_{1}, y_{2})^{T}, ~~N_{0}=(x_{0}, y_{0}, \bar{E}_{0}), \\ & & h(Z(t))=h(h_{1}(y_{1}(t), y_{2}(t)), h_{2}(y_{1}(t), y_{2}(t)), h_{3}(y_{1}(t), y_{2}(t))): R^{2}\longrightarrow R\end{eqnarray*}$

is a smoothing mapping, that is

$$x(t)=x_{0}+y_{1}(t), y(t)=y_{0}+y_{2}(t), \bar{E}(t)=\bar{E_{0}}+h_{3}(y_{1}(t), y_{2}(t)).$$

Then we can obtain the parametric system of system (1.3) as follows

$\begin{equation} \left\{\begin{array}{l} \frac{dy_{1}(t)}{dt}=(y_{1}(t)+x_{0})(r_{1}-b(y_{2}(t)+y_{0})-(\bar{E_{0}}+h_{3}(y_{1}(t), y_{2}(t))), \\ \frac{dy_{2}(t)}{dt}=(y_{2}(t)+y_{0})(r_{2}-k\frac{(y_{2}(t)+y_{0})}{(y_{1}(t)+x_{0})}-(\bar{E_{0}}+h_{3}(y_{1}(t), y_{2}(t))). \end{array}\right. \end{equation}$

(2.3)

Noticing that $g(\psi(Z(t)))=0$, so we can get the l earized system of parametric system (2.3) at (0, 0) as follows

$\begin{equation} \left\{\begin{array}{l} \frac{dy_{1}(t)}{dt}=-bx_{0}y_{2}(t-\tau), \\ \frac{dy_{2}(t)}{dt}=\frac{ky_{0}^{2}}{x_{0}^{2}}y_{1}(t)-\frac{ky_{0}}{x_{0}}y_{2}(t). \end{array}\right. \end{equation}$

(2.4)

From (2.4), we can obtain the characteristic equation of the linearized system of parametric system (2.2) at (0, 0) as follows

$\begin{equation} \lambda^{2}+\frac{ky_{0}}{x_{0}}\lambda+b\frac{ky_{0}^{2}}{x_{0}}e^{-\lambda\tau}=0. \end{equation}$

(2.5)

By eq. (2.5), when $\tau=0$, it is obvious that $\frac{ky_{0}}{x_{0}}>0$ and $\frac{kby_{0}^{2}}{x_{0}}>0$, then, two roots of eq. (2.5) has always negative teal parts, i.e., the positive equilibrium po t of system (1.3) is locally asymptotically stable.

Now, based on the above discussion, we study the local stability around the positive equilibrium po t for system (1.3) and the existence of Hopf bifurcation occurring at the positive equilibrium po t when $\tau>0$.

If $i\omega$ is a root of eq. (2.5), and substituting $i\omega$ ($\omega$ is a positive real number) into eq. (2.5), and separating the real and imaginary parts, two transcendental equations can be obtained as follows

$\begin{eqnarray} & & \omega^{2}=\frac{kby_{0}^{2}}{x_{0}}\cos(\omega\tau), \\ \end{eqnarray}$

(2.6)

$\begin{eqnarray} & & \frac{ky_{0}}{x_{0}}\omega=\frac{kby_{0}^{2}}{x_{0}}\sin(\omega\tau). \end{eqnarray}$

(2.7)

Since $\sin(\omega\tau)^{2}+\cos(\omega\tau)^{2}=1$ and adding (2.6) and (2.7), we obtain

$\begin{equation} \omega^{4}+(\frac{ky_{0}}{x_{0}})^{2}\omega^{2}-(\frac{kby_{0}^{2}}{x_{0}})^{2}=0. \end{equation}$

(2.8)

From (2.8), $(\frac{ky_{0}}{x_{0}})^{2}>0$ and $-(\frac{kby_{0}^{2}}{x_{0}})^{2} < 0$, we can easy find that eq. (2.5) has a unique positive root $\omega_{0}$, that is

$\begin{equation} \omega_{0}=\sqrt{\frac{-(\frac{ky_{0}}{x_{0}})^{2}+\sqrt{(\frac{ky_{0}}{x_{0}})^{2}+4(\frac{kby_{0}^{2}}{x_{0}})^{2}}}{2}}. \end{equation}$

(2.9)

Substituting $\omega_{0}$ into (2.6) and solving for $\tau$, we get

$\begin{equation} \tau_{n}=\frac{1}{\omega_{0}}\arccos(\frac{x_{0}\omega_{0}^{2}}{kby_{0}^{2}})+\frac{2n\pi}{\omega_{0}}, n=0, 1, 2, \cdots. \end{equation}$

(2.10)

Thus when $\tau=\tau_{n}$, the characteristic equation (2.5) has a pair of purely imaginary roots ${i\omega_{0}}$.

Lemma 2.1  Denote by $\lambda_{n}(\tau)=\eta_{n}(\tau)+i\omega_{n}(\tau)$ the root of (2.5) such that $\eta_{n}(\tau_{n})=0$, $\omega_{n}(\tau_{n})=\omega_{0}$, $n=0, 1, 2, \cdots.$ Then the following transversality condition $\eta'_{n}(\tau_{n})$ is satisfied.

Proof Differentiating eq. (2.5) with respect to $\tau$, we obtain

$$(\frac{d\lambda}{d\tau})^{-1}=\frac{[2\lambda x_{0}+ky_{0}-kby_{0}^{2}\tau e^{-\lambda\tau}]}{kby_{0}^{2}\tau e^{-\lambda\tau}}=-\frac{1}{\lambda^{2}}-\frac{\tau}{\lambda}-\frac{1}{\lambda^{2}+\frac{ky_{0}}{x_{0}}\lambda}.$$

Noting that

$${\rm sign}\{{\rm Re}(\frac{d\lambda}{dt})\}_{\lambda=i\omega}={\rm sign}\{{\rm Re}(\frac{d\lambda}{dt})^{-1}\}_{\lambda=i\omega}={\rm sign}\{\frac{2\omega_{0}^{2}+(\frac{ky_{0}}{x_{0}})^{2}}{\omega_{0}^{2}(\omega_{0}^{2}+\frac{ky_{0}}{x_{0}}^{2})}\}>0.$$

The proof is completed.

From the above analysis and [17, 18], we have the following results.

Theorem 2.1  (i) For system (1.3), its positive equilibrium point $Y_{0}$ is locally asymptotically stable for $\tau\in[0, \tau_{0})$ and unstable for $\tau>\tau_{0}$.

(ii) System (1.3) undergoes Hopf bifurcation at the positive equilibrium point $Y_{0}$ for $\tau=\tau_{n}, n=0, 1, 2, \cdots.$

In this section, we investigate the direction of bifurcation and the stability of bifurcation periodic orbits from the positive equilibrium point $Y_{0}$ of system (1.3) at $\tau=\tau_{0}$ by using the normal form approach theory and center manifold theory introduced by Hassard [15].

Now, we re-scare the time by

$$t=\frac{t}{\tau}, \bar{y_{1}}=y_{1}, \bar{y_{2}}=y_{2}, \tau=\tau_{0}+\mu, \bar{Z}=(\bar{y_{1}}, \bar{y_{2}})$$

for simplicity, we continue to use $Z$ said $\bar{Z}$, then the parametric system (2.3) of system (1.3) is equivalent to the following functional differential equation (FDE) system in $C=C([-1, 0], \mathbb{R}^{2})$,

$\begin{equation} \dot{Z}(T)=L_{\mu}Z(T)+f(\mu, Z_{t}), \end{equation}$

(3.1)

where $Z(T)=(y_{1}(t), y_{2}(t))^{T}$, and $L_{\mu}: C\rightarrow\mathbb{R}, f: \mathbb{R}\times C\rightarrow\mathbb{R}$ are given, respectively, by

$$L_{\mu}(\phi)=(\tau_{0}+\mu) \left[ \begin{array}{lcr} 0&0 \\ a_{21}&a_{22} \\ \end{array} \right]\phi^{T}(0)+(\tau_{0}+\mu) \left[ \begin{array}{lcr} 0&a_{12} \\ 0&0 \\ \end{array} \right]\phi^{T}(-1), $$

where

$$ a_{12}=-bx_{0}, ~~ a_{21}=\frac{ky^{2}_{0}}{x^{2}_{0}}, ~~ a_{22}=-\frac{ky_{0}}{x_{0}}, $$

and $f(\mu, \phi)=(\tau_{0}+\mu) \left(\begin{array}{lcr} f_{11} \\ f_{22} \\ \end{array} \right)$, where

$\begin{eqnarray*} & & f_{11}=-b\phi_{1}(0)\phi_{2}(-1), \\ & & f_{22}=-\frac{ky^{2}_{0}}{x^{3}_{0}}\phi^{2}_{1}(0)+\frac{2ky_{0}}{x_{0}}\phi_{1}(0)\phi_{2}(0)-\frac{k}{x_{0}}\phi_{2}^{2}(0)+\cdots, \end{eqnarray*}$

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

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

In fact, choosing

$$\eta(\theta, \mu)=(\tau_{0}+\mu) \left[ \begin{array}{lcr} 0&0 \\ a_{21}&a_{22} \\ \end{array} \right]\delta(\theta)+(\tau_{0}+\mu) \left[ \begin{array}{lcr} 0&a_{12} \\ 0&0 \\ \end{array} \right]\delta(\theta+1), $$

where $\delta (\theta )=\left\{ \begin{array}{*{35}{l}} 0,&\theta \ne 0, \\ 1,&\theta =0. \\ \end{array} \right.$$\text{For}\ \phi \in {{C}^{1}}([-1,0],{{\mathbb{R}}^{2}})$, define

$\begin{equation} A(\mu)\phi(\theta)=\left\{\begin{array}{l} \frac{d\phi(\theta)}{d\theta}, ~~~~~~~~~~~~~~~~~~~~~~-1\leq\theta<0, \\ \displaystyle\int^{0}_{-1}d\eta(\theta, \mu)\phi(\theta), ~~\theta=0. \end{array}\right. \end{equation}$

(3.2)

Then system (3.1) can be rewritten as

$\begin{equation} \dot{Z}=A(\mu)Z_{t}+R(\mu)Z_{t}. \end{equation}$

(3.3)

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

$\begin{equation} A^{*}\psi(s)=\left\{\begin{array}{l} -\frac{d\psi(s)}{ds}, ~~~~~~~~~~~~~~~~~~~~~~~0<s\leq1, \\ \displaystyle\int^{0}_{-1}d\eta^{T}(s, 0)\psi(-s), ~~s=0, \end{array}\right. \end{equation}$

(3.4)

where $\eta^{T}$ is the transpose of the matrix $\eta$.

For $\phi\in C^{1}([-1, 0], \mathbb{R}^{2})$ and for $\psi\in C^{1}([0, 1], (\mathbb{R}^{2})^{*})$, in order to normalize the eigenvectors of operator $A$ and adjoint operator $A^{*}$, we define a bilinear inner product

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

(3.5)

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

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

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

where

$$\beta=-\frac{i\omega}{bx_{0}e^{-i\omega\tau_{0}}}, ~~ \beta^{*}=\frac{i\omega-\frac{ky_{0}}{x_{0}}}{bx_{0}e^{-i\omega\tau_{0}}}, ~~ \bar{G}=(\beta+\bar{\beta^{*}}-bx_{0}\tau_{0}\beta\bar{\beta^{*}}e^{-i\omega\tau_{0}})^{-1}.$$

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

In the reminder of this section, we use the same notations as those in [15]. We first compute the coordinates to describe the center manifold $C_{0}$ at $\mu=0$. Define

$\begin{equation} \dot{z}(t)=\langle q^{*}, Z_{t}\rangle, W(t, \theta)=Z_{t}-2{\rm Re}{z(t)q(\theta)}. \end{equation}$

(3.6)

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

$\begin{equation} 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. \end{equation}$

(3.7)

In fact, $z$ and $\bar{z}$ are local coordinates for center manifold $C_{0}$ in the direction of $q$ and $\bar{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 (3.3), we have

$\begin{eqnarray}\dot{z} & = & i\omega\tau_{0}z+\langle q^{*}(\theta), f(0, W(z, \bar{z}, \theta)+2{\rm Re}[z(t)q(\theta)])\rangle, \nonumber\\ & = & i\omega\tau_{0}z+\bar{q}^{*}f(0, W(z, \bar{z}, 0)+2{\rm Re}[z(t)q(\theta)]), \end{eqnarray}$

(3.8)

rewrite it as

$\begin{equation} \dot{z}=i\omega\tau_{0}z+\mathrm{g}(z, \dot{z}), \end{equation}$

(3.9)

where

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

(3.10)

From (3.3) and (3.8), we have

$\begin{equation} \dot{W}=\dot{Z}_{t}-\dot{z}q-\dot{\bar{z}}q =\left\{\begin{array}{l} AW-2{\rm Re}\{\bar{q}^{*}f(z, \bar{z})q(\theta))\}, ~~~~~~~-1\leq\theta<0, \\ AW-2{\rm Re}\{\bar{q}^{*}f(z, \bar{z}q(\theta))\}+f, ~~\theta=0. \end{array}\right. \end{equation}$

(3.11)

Rewrite (3.11) as

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

(3.12)

where

$\begin{equation} 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. \end{equation}$

(3.13)

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

$\begin{equation} (A-2i\omega\tau_{0})W_{20}(\theta)=-H_{20}(\theta), AW_{11}(\theta)=-H_{11}(\theta). \end{equation}$

(3.14)

Notice that

$$q(\theta)=(1, \beta)^{T}e^{i\omega\tau_{0}\theta}, ~~ q^{*}(0)=G(\beta^{*}, 1), $$

and (3.6) we obtain

$\begin{eqnarray*} & & y_{1t}(0)=z+\bar{z}+W_{20}^{1}(0)\frac{z^{2}}{2}+W_{11}^{1}(0)z\bar{z}+W_{02}^{1}(0)\frac{\bar{z}^{2}}{2}+\cdots, \\ & & y_{2t}(0)=\beta z+\bar{\beta}\bar{z}+W_{20}^{2}(0)\frac{z^{2}}{2}+W_{11}^{2}(0)z\bar{z}+W_{02}^{2}(0)\frac{\bar{z}^{2}}{2}+\cdots, \\ & & y_{1t}(-1)=\beta ze^{-i\omega\tau_{0}}+\bar{\beta}\bar{z}e^{i\omega\tau_{0}}+W_{20}^{1}(-1)\frac{z^{2}}{2}+W_{11}^{1}(-1)z\bar{z}+W_{02}^{1}(-1)\frac{\bar{z}^{2}}{2}+\cdots, \\ & & y_{2t}(-1)=\beta ze^{-i\omega\tau_{0}}+\bar{\beta}\bar{z}e^{i\omega\tau_{0}}+W_{20}^{2}(-1)\frac{z^{2}}{2}+W_{11}^{2}(-1)z\bar{z}+W_{02}^{2}(-1)\frac{\bar{z}^{2}}{2}+\cdots. \end{eqnarray*}$

According to (3.8) and (3.9), we know that

$\begin{equation} \mathrm{g}(z, \bar{z})=\bar{q}^{*}(0)f_{0}(z, \bar{z})=\bar{G}\tau_{0}(\bar{\beta}^{*}, 1)\left(\begin{array}{lcr} f^{0}_{11} \\ f^{0}_{22} \\ \end{array} \right), \end{equation}$

(3.15)

where

$\begin{eqnarray*} & & f^{0}_{11}= -by_{1t}(0)y_{2t}(-1), \\ & & f^{0}_{22}=-\frac{ky_{0}^{2}}{x_{0}^{3}}y_{1t}(0)^{2}+\frac{2ky_{0}}{x_{0}}y_{1t}(0)y_{2t}(0)-\frac{k}{x_{0}}y_{2t}(0)^{2}+\cdots.\end{eqnarray*}$

By (3.7), it follows that

$\begin{eqnarray*}\mathrm{g}(z, \bar{z}) & = & \bar{G}\tau_{0}\{-b\bar{\beta}^{*}[z+\bar{z}+W_{20}^{1}(0)\frac{z^{2}}{2}+W_{11}^{1}(0)z\bar{z}+W_{02}^{1}(0)\frac{\bar{z}^{2}}{2}]\\ & & \times[\beta ze^{-i\omega\tau_{0}}+\bar{\beta}\bar{z}e^{i\omega\tau_{0}}+W_{20}^{2}(-1)\frac{z^{2}}{2}+W_{11}^{2}(-1)z\bar{z}+W_{02}^{2}(-1)\frac{\bar{z}^{2}}{2}]\\ & & -\frac{ky_{0}^{2}}{x_{0}^{3}}[z+\bar{z}+W_{20}^{1}(0)\frac{z^{2}}{2}+W_{11}^{1}(0)z\bar{z}+W_{02}^{1}(0)\frac{\bar{z}^{2}}{2}]^{2}\\ & & +\frac{2ky_{0}}{x_{0}}[z+\bar{z}+W_{20}^{1}(0)\frac{z^{2}}{2}+W_{11}^{1}(0)z\bar{z}+W_{02}^{1}(0)\frac{\bar{z}^{2}}{2}]\\ & & \times[\beta z+\bar{\beta}\bar{z}+W_{20}^{2}(0)\frac{z^{2}}{2}+W_{11}^{2}(0)z\bar{z}+W_{02}^{2}(0)\frac{\bar{z}^{2}}{2}]\\ & & -\frac{k}{x_{0}}[\beta z+\bar{\beta}\bar{z}+W_{20}^{2}(0)\frac{z^{2}}{2}+W_{11}^{2}(0)z\bar{z}+W_{02}^{2}(0)\frac{\bar{z}^{2}}{2}]^{2}+\cdots\}. \end{eqnarray*}$

That is,

$\begin{eqnarray*}\mathrm{g}(z, \bar{z}) & = & \bar{G}\tau_{0}\{z^{2}[-b\beta\bar{\beta}^{*}e^{-i\omega\tau_{0}}-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\beta-\frac{k}{x_{0}}\beta^{2}]\\ & & +z\bar{z}[-b\bar{\beta}^{*}(\bar{\beta}e^{i\omega\tau_{0}}+\beta e^{-i\omega\tau_{0}})-2\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}(\bar{\beta}+\beta)-2\frac{k}{x_{0}}\beta\bar{\beta}]\\ & & +\bar{z}^{2}[-b\bar{\beta}^{*}\bar{\beta}e^{i\omega\tau_{0}}-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\bar{\beta}-\frac{k}{x_{0}}\bar{\beta^{2}}]\\ & & +z^{2}\bar{z}[(-b\beta\bar{\beta}^{*}e^{i\omega\tau_{0}}-2\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\beta)W_{11}^{1}(0)+(\frac{2ky_{0}}{x_{0}}-2\frac{k}{y_{0}})W_{11}^{2}(0)\\ & & +(\frac{-b\beta\bar{\beta}^{*}e^{i\omega\tau_{0}}}{2}-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{\frac{2ky_{0}}{x_{0}}\bar{\beta}}{2})W_{20}^{1}(0)+(\frac{\frac{2ky_{0}}{x_{0}}}{2}-\frac{k}{x_{0}}\bar{\beta})W_{20}^{2}(0)\\ & & -b\bar{\beta}^{*}W_{11}^{2}(-1)+\frac{-b\bar{\beta}^{*}}{2}W_{20}^{2}(-1)]+\cdots\}.\end{eqnarray*}$

Comparing the coefficients with (3.10), it follows that

$\begin{eqnarray*} \mathrm{g}_{20} & = & 2\bar{G}\tau_{0}[-b\beta\bar{\beta}^{*}e^{-i\omega\tau_{0}}-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\beta-\frac{k}{x_{0}}\beta^{2}], \\ \mathrm{g}_{11} & = & \bar{G}\tau_{0}[-b\bar{\beta}^{*}(\bar{\beta}e^{i\omega\tau_{0}}+\beta e^{-i\omega\tau_{0}})-2\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}(\bar{\beta}+\beta)-2\frac{k}{x_{0}}\beta\bar{\beta}], \\ \mathrm{g}_{02} & = & [-b\bar{\beta}^{*}\bar{\beta}e^{i\omega\tau_{0}}-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\bar{\beta}-\frac{k}{x_{0}}\bar{\beta^{2}}], \\ \mathrm{g}_{21} & = & [(-b\beta\bar{\beta}^{*}e^{i\omega\tau_{0}}-2\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\beta)W_{11}^{1}(0)+(\frac{-b\beta\bar{\beta}^{*}e^{i\omega\tau_{0}}}{2}-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{\frac{2ky_{0}}{x_{0}}\bar{\beta}}{2})W_{20}^{1}(0)\\ & & +(\frac{2ky_{0}}{x_{0}}-2\frac{k}{y_{0}})W_{11}^{2}(0)+(\frac{\frac{2ky_{0}}{x_{0}}}{2}-\frac{k}{x_{0}}\bar{\beta})W_{20}^{2}(0)-b\bar{\beta}^{*}W_{11}^{2}(-1)+\frac{-b\bar{\beta}^{*}}{2}W_{20}^{2}(-1)]. \end{eqnarray*}$

Now we compute $W_{20}(\theta)$ and $W_{11}(\theta).$ From (3.11) and (3.15), we have that for $\theta\in[-1, 0), $

$\begin{equation} H(z, \bar{z}, \theta)=-2{\rm Re}\{\bar{q}^{*}(0)f(z, \bar{z})q(\theta)\}=-\mathrm{g}(z, \bar{z})q(\theta)-\bar{\mathrm{g}}(z, \bar{z})\bar{q}(\theta). \end{equation}$

(3.16)

Comparing the coefficients with (3.13), we can obtain that

$\begin{equation} H_{20}(\theta)=-\mathrm{g}_{20}q(\theta)-\bar{\mathrm{g}}_{02}\bar{q}(\theta), H_{11}(\theta)=-\mathrm{g}_{11}q(\theta)-\bar{\mathrm{g}}_{11}\bar{q}(\theta). \end{equation}$

(3.17)

Substituting the above equalities into (3.14), it follows that

$\begin{equation} \left\{\begin{array}{l} \dot{W}_{20}(\theta)=2i\omega W_{20}(\theta)+\mathrm{g}_{20}q(\theta)+\bar{\mathrm{g}}_{02}\bar{q}(\theta), \\ \dot{W}_{11}(\theta)=\mathrm{g}_{11}q(\theta)+\bar{\mathrm{g}}_{11}\bar{q}(\theta). \end{array}\right. \end{equation}$

(3.18)

Solving (3.18), we have

$\begin{equation} \left\{\begin{array}{l} W_{20}(\theta)=\frac{i\mathrm{g}_{20}}{\tau_{0}\omega}q(0)e^{i\omega\tau_{0}\theta}+\frac{i\bar{\mathrm{g}}_{02}}{3\omega\tau_{0}}\bar{q}(0)e^{-i\omega\tau_{0}\theta}+Ee^{2i\omega\tau_{0}}, \\[0.5em] W_{11}(\theta)=-\frac{i\mathrm{g}_{11}}{\tau_{0}\omega}q(0)e^{i\omega\tau_{0}\theta}+\frac{i\bar{\mathrm{g}}_{11}}{\omega\tau_{0}}\bar{q}(0)e^{-i\omega\tau_{0}\theta}+F. \end{array}\right. \end{equation}$

(3.19)

In what follows, we seek appropriate $E$ and $F$ in (3.19). From (3.11) and (3.15), we have

$\begin{eqnarray} & & H_{20}(0)=-\mathrm{g}_{20}q(0)-\bar{\mathrm{g}}_{02}\bar{q}(0)+2\tau_{0} \left( \begin{array}{lcr} d_{1} \\ d_{2} \\ \end{array} \right), \\ \end{eqnarray}$

(3.20)

$\begin{eqnarray} & & H_{11}(0)=-\mathrm{g}_{11}q(0)-\bar{\mathrm{g}}_{11}\bar{q}(0)+2\tau_{0} \left( \begin{array}{lcr} d_{1} \\ d_{2} \\ \end{array} \right), \end{eqnarray}$

(3.21)

where

$\begin{eqnarray*} & & d_{1}=-b\beta^{*}e^{-i\omega\tau_{0}}, d_{2}=-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}\beta-\frac{k}{x_{0}}\beta^{2}, \\ & & d_{3}=-bRe(\beta^{*}e^{i\omega\tau_{0}}), d_{4}=-\frac{ky_{0}^{2}}{x_{0}^{3}}+\frac{2ky_{0}}{x_{0}}Re(\beta)-\frac{k}{x_{0}}\beta\bar{\beta}. \end{eqnarray*}$

Substituting (3.19)--(3.21) into (3.14) and noting that

$\begin{eqnarray*} & & (i\omega\tau_{0}I-\int^{0}_{-1}e^{i\omega\tau_{0}\theta}d\eta(\theta))q(0)=0, \\ & & (-i\omega\tau_{0}I-\int^{0}_{-1}e^{-i\omega\tau_{0}\theta}d\eta(\theta))q(0)=0.\end{eqnarray*}$

We obtain

$\begin{eqnarray} \left(\begin{array}{lcr} 2i\omega&bx_{0}e^{-2i\omega_{0}\tau_{0}} \\ -\frac{ky_{0}^{2}}{x_{0}^{2}}&2i\omega+\frac{ky_{0}}{x_{0}} \\ \end{array} \right) E =2 \left( \begin{array}{lcr} d_{1} \\ d_{2} \\ \end{array} \right), \\ \end{eqnarray}$

(3.22)

$\begin{eqnarray} \left(\begin{array}{lcr} 0&bx_{0}\\ -\frac{ky_{0}^{2}}{x_{0}^{2}}&\frac{ky_{0}}{x_{0}} \\ \end{array} \right) F =2 \left( \begin{array}{lcr} d_{3} \\ d_{4} \\ \end{array} \right). \end{eqnarray}$

(3.23)

It is easy to obtain $E$ and $F$ from (3.22) and (3.23), that is

$\begin{eqnarray*} & & E^{(1)}=\frac{(2i\omega+\frac{ky_{0}}{x_{0}})d_{1}-bx_{0}e^{-2i\omega_{0}\tau_{0}}}{\frac{2ky_{0}i\omega}{x_{0}}+\frac{bky_{0}^{2}}{x_{0}}e^{-2i\omega_{0}\tau_{0}}-4\omega^{2}}, \\ & & E^{(2)=}\frac{\frac{ky_{0}^{2}}{x_{0}^{2}}d_{3}+2i\omega d_{4}}{\frac{2ky_{0}i\omega}{x_{0}}+\frac{bky_{0}^{2}}{x_{0}}e^{-2i\omega_{0}\tau_{0}}-4\omega^{2}}, \\ & & F^{(1)}=\frac{\frac{ky_{0}}{x_{0}}d_{3}-bx_{0}d_{4}}{\frac{bky_{0}^{2}}{x_{0}}}, F^{(2)}=\frac{d_{3}}{bx_{0}}.\end{eqnarray*}$

Therefore we can compute the following quantities

$\begin{eqnarray*} & & \mathcal{C}_{1}(0)=\frac{i}{2\tau_{0}\omega_{0}}(\mathrm{g}_{20}\mathrm{g}_{11}-2|\mathrm{g}_{11}|^{2}-\frac{1}{3}|\mathrm{g}_{02}|^{2})+\frac{\mathrm{g}_{21}}{2}, \\ & & \mu_{2}=-\frac{{\rm Re}\mathcal{C}_{1}(0)}{{\rm Re}\lambda'(\tau_{0})}, t_{2}=-\frac{{\rm Im}\mathcal{C}_{1}(0)+\mu_{2}{\rm Im}\lambda'(\tau_{0})}{\tau_{0}\omega_{0}}, \beta_{2}=2{\rm Re}\mathcal{C}_{1}(0), \end{eqnarray*}$

which determine the direction of Hopf bifurcation and stability of bifurcated periodic solutions of system (1.3) at the critical value $\tau_{0}$.

Theorem 3.1 (i) The direction of Hopf bifurcation is determined by the sign of $\mu_{2}$: the Hopf bifurcation is supercritical if $\mu_{2}>0$ and the Hopf bifurcation is subcritical if $\mu_{2} < 0$.

(ii) The stability of bifurcated periodic solution is determined by $\beta{2}$: the periodic solution are stable if $\beta_{2}>0$ and unstable if $\beta_{2} < 0$.

(iii) The period of bifurcation periodic solution is determined by $t_{2}$: the period increase if $t_{2}>0$, decrease if $t_{2} < 0$.

As an example we consider the differential-algebraic predator-prey system (1.3) with the parameters $r_{1}=1.6$, $r_{2}=1.3$, $b=k=m=0.5$, $p_{1}=7$, $p_{2}=6$, $c_{1}=5$, $c_{2}=3$, that is,

$\begin{equation} \left\{\begin{array}{l} \frac{dx}{dt}=x(t)(1.6-0.5y(t-\tau)-E(t)), \\ \frac{dy}{dt}=y(t)(1.3-0.5\frac{y(t)}{x(t)}-E(t)), \\ 0=E(t)(7x(t)+6y(t)-5-3)-0.5.\\ \end{array}\right. \end{equation}$

(4.1)

And by the discussions in Section 2 and Section 3, we determine the stability of the positive equilibrium point and Hopf bifurcation. Here, for convenience, we only discuss one of the positive equilibrium point $Y_{0}$ of system (4.1), and others positive equilibrium points of system (4.1) can be similar studied. we can easily get $Y_{0}=(2.0053, 3.1480, 0.0256)$, and by computing, we get $\omega_{0}=0.9942$, $\tau_{0}=0.6473$. So by Theorem 2.1, the equilibrium point $Y_{0}$ is asymptotically stable when $\tau\in[0, \tau_{0})=[0, 0.6473)$ and unstable when $\tau>0.6473$.

When $\tau=0$, we can easily show that the positive equilibrium point

$$Y_{0}=(2.0053, 3.1480, 0.0256)$$

is asymptotically stable.

By the theory of Hassard [15], as it is discussed in former section, we also determine the direction of Hopf bifurcation and the other properties of bifurcating periodic solution. By computing, we can obtain the following values $\mathcal{C}_{1}(0)=0.5303-0.4428i$, $\lambda'(\tau_{0})=1.6352+ 1.1431i$, it follows that $\mu_{2}= -0.3243 < 0$, $\beta_{2}=1.0607>0$, $t_{2}=1.2643>0$, from which and Theorem 3.1 we conclude that the Hopf bifurcation of system(4.1) occurring at $\tau_{0}=0.6473$ is subcritical and the bifurcating periodic solution exists when $\tau$ cross $\tau_{0}$ to the left and the bifurcating periodic solution is unstable.

By Theorem 3.1, the positive equilibrium po t $Y_{0}$ of system (4.1) is locally asymptotically stable when $\tau=0.62 < \tau_{0}$ as is illustrated by computer simulation Fig. 1. And periodic solutions occur from $Y_{0}$ when $\tau=0.682>\tau_{0}$ as is illustrated by computer simulation Fig. 2.

Nowadays, economic profit is a very important factor for governments, merchants, and even citizen, and the harvested biological resources in the predator-prey systems are usually sold as commodities in the market in order to achieve the economic interest. So modelling and qualitative analysis for bio-economic system are necessary.

Compared with most other researches on dynamics of predator-prey population, see [1, 5, 12, 18], the main contribution of this paper lies in the following aspect. The predator-prey system we consider incorporate delay and non-selective harvesting, which could make our model more realistic and the analysis result in this paper is more scientific. So our paper provide a new ideal and a efficacious method for the qualitative analysis of the Hopf bifurcation of the differential-algebraic biological economic system. In addition, stage structure, diffusion effects, disease effects may be incorporated into our bio-economic system, which would make the bio-economic system exhibit much more complicated dynamics.