According to the lack of biological resources on the earth, more and more people increasingly realized the importance of the modelling and research of biological system. The predator-prey is one of the most popular models that many researchers [1-8] have studied and acquired some valuable characters of dynamic behavior. For example, the stability of equilibrium, Hopf bifurcation, flip bifurcation, limit cycle and other relevant conducts. At the same time, the development and utilization of biological resources and artificial arrest has been researched commonly in the fields of fishery, wildlife and forestry management by some experts [9-11]. Most of them choose differential equations and difference equations to research biological models. It is well known that economic profit become more and more important and take a fundamental gradually situation in social development. In recent years, biological economic systems have been researched by many authors [12-16], who describe the system by differential-algebraic equations or differential-difference-algebraic equations.

Basic analysis model which applied by differential-algebraic equations and differential-difference-algebraic equations are familiar at present. However, there still exist some disadvantages in many systems such as harvesting function. In this paper, the main research is the stability and Hopf bifurcation of a biological-algebraic biological economic system, which is changed in some details and meaningful.

Our basic model is based on the following ratio-dependent predator-prey system with harvest

$ \begin{equation} \label{eq:1} \left\{\begin{array}{l} \dot{u}=u(r_{1}-\epsilon v), \\ \dot{v}=v(r_{2}-\theta \frac{v}{u})-\alpha vE^{*}, \end{array} \right. \end{equation} $

(1.1)

where $u$ and $v$ represent the predator density and prey density at time $t$, respectively, $\epsilon$, $\theta$ and $\alpha$ are all positive constants, and $r_{1}$ and $r_{2}$ stand for the densities of predator and prey populations, and $E$ represents harvesting effort. $\alpha E v$ denotes that the harvests for predator population are proportional to their densities at time $t$.

In 1954, Gordon [17] studied the effect of the harvest effort on ecosystem form an economic perspective and proposed the following economic principle:

$ \begin{aligned} {\rm Net~ Economic~ Revenue}(NER)={\rm Total~ Revenue}(TR)-{\rm Total~ Cost}(TC). \end{aligned} $

Associated with the system(1.1), an algebraic equation which considers the economic profit $m$ of the harvest effort on predator can be established as follows

$ \begin{aligned} E(t)(pv-q)=m, \end{aligned} $

where $E(t)$ represents the harvest effort, $p$ denotes harvesting reward per unit harvesting effort for unit weight, $c$ represents harvesting cost per unit harvesting effort. Combining the economic theory of fishery resources, we can establish a differential algebraic biological economic system

$ \begin{equation} \left\{\begin{array}{l} \dot{u}=u(r_{1}-\epsilon v), \\ \dot{v}=v(r_{2}-\theta \frac{v}{u})-\alpha vE^{*}, \\ 0=E^{*}(pv-q)-m. \end{array} \right. \end{equation} $

(1.2)

Nevertheless, the capture effect to predator is not always shown in the liner in nature based on many factors that can affect the predation such as the ability of search, illness and death. Therefore, the harvesting function of the system(1.2) has been modified as follows

$ \begin{equation} \left\{\begin{array}{l} \dot{u}=u(r_{1}-\epsilon v), \\ \dot{v}=v(r_{2}-\theta \frac{v}{u})-\alpha E^{*}\frac{v}{1+\gamma v}, \\ 0=E^{*}(p\frac{v}{1+\gamma v}-q)-m. \end{array} \right. \end{equation} $

(1.3)

To simplify the system (1.2), we use these dimensionless variables

$ \begin{aligned} x=\frac{\epsilon}{\theta}u, ~~ y=\epsilon v, ~~ E=\alpha E^{*}, ~~ \beta=\frac{\gamma}{\epsilon}, ~~ c=\frac{\epsilon q}{p}, ~~ \mu=\frac{\alpha\epsilon m}{p} \end{aligned} $

and then obtain the following system

$ \begin{equation} \left\{\begin{array}{l} \dot{x}=x(r_{1}-y), \\ \dot{y}=y(r_{2}-\frac{y}{x})-E\frac{y}{1+\beta y}, \\ 0=E(\frac{y}{1+\beta y}-q)-\mu. \end{array} \right. \end{equation} $

(1.4)

For simplicity, let

$ \begin{eqnarray*} &&f(Z, E, \mu)=\left(\begin{array}{cc} f_{1}(Z, E, \mu) \\ f_{2}(Z, E, \mu) \end{array}\right)=\left(\begin{array}{cc} x(r_{1}-y) \\ y(r_{2}-\frac{y}{x})-E\frac{y}{1+\beta y} \end{array}\right), \\ &&g(Z, E, \mu)=E(\frac{y}{1+\beta y}-q)-\mu, \end{eqnarray*} $

where $Z=(x, y)^{T}$, $\mu$ is a bifurcation parameter, which will be defined in the follows.

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

Next, the paper will be organized as follows. In Section 2, the stability of the positive equilibrium point is discussed by corresponding characteristic equation of system (2.2). In Section 3, provide Hopf bifurcation analysis of the system (1.4). In Section 4, use numerical simulations to illustrate the effectiveness of result. Then, give a brief conclusion in Section 5.

It is obvious that there exists an equilibrium in $R^{3}_{+}$ if only if this point $\chi_{0}:=(x_{0}, y_{0}, E_{0})^{T}$ is a real solution of the equations

$ \left\{ \begin{array}{*{35}{l}} x({{r}_{1}}-y)=0, \\ y({{r}_{2}}-\frac{y}{x})-E\frac{y}{1+\beta y}=0, \\ E(\frac{y}{1+\beta y}-q)-\mu =0. \\ \end{array} \right. $

(2.1)

By the calculation, we get

$ \begin{aligned} \chi_{0}=(x_{0}, y_{0}, E_{0})=(\frac{r_{1}G_{0}}{r_{2}G_{0}-\mu}, r_{1}, \frac{\mu(1+\beta r_{1})}{G_{0}}), \end{aligned} $

where

$ \begin{aligned} G=G(y)=y-c(1+\beta y), ~~~ G_{0}=G(y=r_{1})=r_{1}-c(1+\beta r_{1}). \end{aligned} $

According to this analysis procedure, this essay only concentrate on the interior equilibrium of the system (1.4). Based on the ecology meaningful of the interior equilibrium, the predator and the harvest effort to predator are all exist that it is the key point to the study. Thus, a simple assumption that the inequality $0<\mu<r_{2}G_{0}$ holds in this paper. Following, we use the linear transformation $\chi^{T}=QM^{T}$, where

$ \begin{aligned} M=(X, Y, \bar{E}), ~~~Q=\left(\begin{array}{ccc} 1&0&0 \\ 0&1&0 \\ 0&-\displaystyle\frac{E_{0}}{y_{0}(1+\beta y_{0})-c(1+\beta y_{0})^{2}}&1 \end{array}\right). \end{aligned} $

Then we obtain $D_{\chi}g(\chi)Q=(0, 0, \frac{y_{0}}{1+\beta y_{0}}-c)$, $M=(x, y, \frac{\mu}{G_{0}^{2}}y+E)$, where

$ \begin{aligned} \bar{E}=\frac{\mu}{G_{0}^{2}}y+E=\frac{E_{0}}{y_{0}(1+\beta y_{0})-c(1+\beta y_{0})^{2}}y+E. \end{aligned} $

Next, let $E=\bar{E}-\frac{\mu}{G_{0}^{2}}y$.Thus we transform the system (1.4) into

$ \begin{equation} \left\{\begin{array}{l} \dot{X}=X(r_{1}-Y), \\ \dot{Y}=Y[r_{2}-\frac{Y}{X}+(\frac{\mu}{G_{0}^{2}}Y-\bar{E})\frac{1}{1+\beta Y}], \\ 0=(\bar{E}-\frac{\mu}{G_{0}^{2}}Y)(\frac{Y}{1+\beta Y}-c)-\mu. \end{array} \right. \end{equation} $

(2.2)

From Section 1, we obtain

$ \begin{eqnarray*} &&f(M, \mu)=\left(\begin{array}{cc} f_{1}(M, \mu) \\ f_{2}(M, \mu) \end{array}\right)=\left(\begin{array}{cc} X(r_{1}-Y) \\ Y[r_{2}-\frac{Y}{X}+(\frac{\mu}{G_{0}^{2}}Y-\bar{E})\frac{1}{1+\beta Y}] \end{array}\right), \\ &&g(M, \mu)=(\bar{E}-\frac{\mu}{G_{0}^{2}}Y)(\frac{Y}{1+\beta Y}-c)-\mu. \end{eqnarray*} $

For the system (2.2), we consider the local parametric $\psi$, which defined as follows

$ \begin{aligned} (X, Y, \bar{E}^{T})=\psi(\bar{Z})=M^{T}_{0}+U_{0}\bar{Z}+V_{0}h(\bar{Z}), \end{aligned} $

where

$ \begin{aligned} U_{0}=\left(\begin{array}{ccc} 1&0 \\ 0&1 \\ 0&0 \end{array}\right), ~~~V_{0}=\left(\begin{array}{ccc} 0 \\ 0 \\ 1 \end{array}\right), ~~~\bar{Z}=(y_{1}, y_{2})^{T}, ~~~M_{0}=(X_{0}, Y_{0}, \bar{E_{0}}), \end{aligned} $

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

$ \begin{equation} \left\{\begin{array}{l} \dot{y_{1}}=f_{1}(\mu, \psi(\mu, M)), \\ \dot{y_{2}}=f_{2}(\mu, \psi(\mu, M)). \end{array} \right. \end{equation} $

(2.3)

More details about the definition can be found in the literature [18]. Based on the system (2.3), we can get Jacobian matrix $E(M_{0})$, which takes the form of

$ \begin{eqnarray*} E(M_{0})&=&\left(\begin{array}{cc} D_{M}f_{1}(M_{0 }) \\ D_{M}f_{2}(M_{0 }) \end{array}\right)\left(\begin{array}{cc} D_{M}G_{1}(M_{0 }) \\ U^{T}_{0} \end{array}\right)^{-1}\left(\begin{array}{cc} 0 \\ I_{1} \end{array}\right)\\ &=&\left(\begin{array}{cc} D_{X}f_{1}(M_{0 })&D_{Y}f_{1}(M_{0 }) \\ D_{X}f_{2}(M_{0 })&D_{Y}f_{2}(M_{0 }) \end{array}\right)\\ &=&\left(\begin{array}{cc} 0&-\frac{r_{1}G_{0}}{r_{2}G_{0}-\mu} \\ (r_{2}-\frac{\mu}{G_{0}})^{2}&-r_{2}\frac{\mu}{G_{0}^{2}}(2G_{0}+c) \end{array}\right). \end{eqnarray*} $

Then, the following theorem summarizes the stability of the positive equilibrium point of system (1.4).

Theorem 2.1  For the system (2.2)

(ⅰ) If $(r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}})^{2}\geqslant4r_{1}\frac{r_{2}G_{0}-\mu}{G_{0}}$ and $\mu<\min{\left\{\frac{r_{2}G_{0}^{2}}{2G_{0}+c}, r_{2}G_{0}\right\}}, $ the positive equilibrium point of system (1.4) is asymptotically stable; Otherwise when $\frac{r_{2}G_{0}^{2}}{2G_{0}+c}<\mu<r_{2}G_{0}, $ the positive equilibrium point of system (1.4) is unstable.

(ⅱ) If $(r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}})^{2}<4r_{1}\frac{r_{2}G_{0}-\mu}{G_{0}}$ and $\mu<\min{\left\{\frac{r_{2}G_{0}^{2}}{2G_{0}+c}, r_{2}G_{0}\right\}}, $ the positive equilibrium point of system (1.4) is a sink; Otherwise when $\frac{r_{2}G_{0}^{2}}{2G_{0}+c}<\mu<r_{2}G_{0}, $ the positive equilibrium point of system (1.4) is a source.

Proof  First, the characteristic equation of the matrix $E(M_{0})$ can be written as

$ \begin{equation} \lambda^{2}+(r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}})\lambda+r_{1}\frac{r_{2}G_{0}-\mu}{G_{0}}=0. \end{equation} $

(2.4)

Now donate $\Delta$ by

$ \begin{aligned} \Delta=(r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}})^{2}-4r_{1}\frac{r_{2}G_{0}-\mu}{G_{0}}. \end{aligned} $

If $\Delta\geq0$ and $\mu<\min{\left\{\frac{r_{2}G_{0}^{2}}{2G_{0}+c}, r_{2}G_{0}\right\}}, $ the eq. (2.4) has two negative real roots; when $\Delta\geq0$ and $\frac{r_{2}G_{0}^{2}}{2G_{0}+c}<\mu<r_{2}G_{0}, $ the eq. (2.4) has two positive real roots. We can obtain the part (ⅰ) of the theorem by the Routh-Hurwitz criteria, the part (ⅱ) can be similar proofed. Thus, we complete the proof of Theorem 2.1.

Remark 1  The local stability of $\chi_{0}$ is equivalent to the local stability of $M_{0}$.

Remark 2  When the roots of eq. (2.4) exist zero real parts, the system (1.4) will occur bifurcation, which will be discussed in the Section 3.

In this section, we discuss the Hopf bifurcation from the equilibrium point $\chi_{0}$ by choosing $\mu$ as the bifurcation parameter. Based on the Hopf bifurcation theorem in [19], we need find some sufficient conditions.

According to the definition of $\Delta$, we obtain

$ \begin{aligned} J_{\pm}=\frac{(2r_{2}G_{0}+r_{2}c-2r_{1}G_{0})G_{0}^{2}}{(2G_{0}+c)^{2}}\pm\sqrt{\frac{(2r_{2}G_{0}+r_{2}c-2r_{1}G_{0})^{2}G_{0}^{4}}{(2G_{0}+c)^{4}}+B}=A\pm\sqrt{A^{2}+B}, \end{aligned} $

where

$ \begin{eqnarray*} &&A=\frac{(2r_{2}G_{0}+r_{2}c-2r_{1}G_{0})G_{0}^{2}}{(2G_{0}+c)^{2}}, \\ &&B=(4r_{1}r_{2}-r_{2}^{2})\cdot\frac{G_{0}^{4}}{(2G_{0}+c)^{2}}, \end{eqnarray*} $

here, we assume that $A^{2}+B\geqslant0$ in this paper.

Thus, for eq. (2.4), if $B>0$ and $0<\mu<\min{\left\{r_{2}G_{0}, J_{+}\right\}}$. Eq. (2.4) has one pair of imaginary roots.When $B>0, A>0, J_{-}<r_{2}G_{0}$ and $J_{-}<\mu<\min{\left\{r_{2}G_{0}, J_{+}\right\}}$, eq. (2.4) has one pair of imaginary roots.

In the case of meet the above conditions, we can get the roots as follows:

$ \begin{eqnarray*} \lambda_{1, 2}&=&-\frac{1}{2}(r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}})\pm\sqrt{r_{1}\frac{r_{2}G_{0}-\mu}{G_{0}}-\frac{1}{4} (r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}})^{2}}\\ &=&\alpha(\mu)\pm i\omega(\mu), \end{eqnarray*} $

where

$ \begin{eqnarray*} &&\alpha(\mu)=-\frac{1}{2}(r_{2}-\mu\frac{2G_{0}+c}{G_{0}^{2}}), \\ &&i\omega(\mu)=\sqrt{r_{1}\frac{r_{2}G_{0}-\mu}{G_{0}}-\alpha^{2}(\mu)}. \end{eqnarray*} $

By calculating, we obtain

$ \begin{align} & {{\mu }_{0}}=\frac{{{r}_{2}}G_{0}^{2}}{2{{G}_{0}}+c}, ~~~{\alpha }'({{\mu }_{0}})=\frac{2{{G}_{0}}+c}{G_{0}^{2}}>0, \\ & {{\omega }_{0}}=\omega ({{\mu }_{0}})=\sqrt{{{r}_{1}}{{r}_{2}}(\frac{1}{{{G}_{0}}}-\frac{{{G}_{0}}}{2{{G}_{0}}+c})}. \\ \end{align} $

(3.1)

Eq. (3.1) indicates that eq. (2.2) occurs Hopf bifurcation at $\mu_{0}$.

In order to calculate the Hopf bifurcation, we need to lead the normal form of the system (2.2) as follows

$ \begin{eqnarray} \dot{y_{1}}&=&\alpha(\mu)y_{1}-\omega(\mu)y_{2}+\frac{1}{2}a_{11}^{1}y_{1}^{2}+a_{12}^{1}y_{1}y_{2}+\frac{1}{2}a_{22}^{1}y_{2}^{2}+\frac{1}{6}a_{111}^{1}y_{1}^{3}\nonumber\\ &&+\frac{1}{2}a_{112}^{1}y^{2}_{1}y_{2}+\frac{1}{2}a_{122}^{1}y_{1}y_{2}^{2}+\frac{1}{6}a_{222}^{1}y^{3}_{2}+o(|\bar{Z}|^{4}), \nonumber\\ \dot{y_{2}}&=&\omega(\mu)y_{1}+\alpha(\mu)y_{2}+\frac{1}{2}a_{11}^{2}y_{1}^{2}+a_{12}^{2}y_{1}y_{2}+\frac{1}{2}a_{22}^{2}y_{2}^{2}+\frac{1}{6}a_{111}^{2}y_{1}^{3}\nonumber\\ &&+\frac{1}{2}a_{112}^{2}y^{2}_{1}y_{2}+\frac{1}{2}a_{122}^{2}y_{1}y_{2}^{2}+\frac{1}{6}a_{222}^{2}y^{3}_{2}+o(|\bar{Z}|^{4}). \end{eqnarray} $

(3.2)

From eq. (2.3), we have

$ \begin{eqnarray*} &&f_{1}(M, \mu)=X(r_{1}-Y), \\ &&f_{2}(M, \mu)=Y[r_{2}-\frac{Y}{X}+(\frac{\mu}{G_{0}^{2}}Y-\bar{E})\frac{1}{1+\beta Y}], \\ &&g(M, \mu)=(\bar{E}-\frac{\mu}{G_{0}^{2}}Y)(\frac{Y}{1+\beta Y}-c)-\mu. \end{eqnarray*} $

Then, we can easily obtain

$ \begin{eqnarray*} &&D_{N}f_{1}(M, \mu)=(r_{1}-Y, -X, 0), \\ &&D_{N}f_{2}(M, \mu)=(\frac{Y^{2}}{X^{2}}, r_{2}-\frac{2Y}{X}+F, \frac{-Y}{1+\beta Y}), \\ &&D_{N}g(M, \mu)=(0, \frac{\mu c}{G_{0}^{2}}-F, \frac{Y}{1+\beta Y}-c), \end{eqnarray*} $

where

$ F=\frac{\mu}{G_{0}^{2}}\cdot\frac{2Y+\beta Y^{2}}{(1+\beta Y)^{2}}-\frac{\bar{E}}{(1+\beta Y)^{2}} $

and

$ \begin{eqnarray*}D_{\psi}(\bar{Z}, \mu)&=&\left(\begin{array}{cc} D_{N}g(M, \mu) \\ U^{T}_{0} \end{array}\right)^{-1}\left(\begin{array}{ccc} 0&0 \\ 1&0 \\ 0&1 \end{array}\right)\\ &=&\left(\begin{array}{ccc} 0&1&0 \\ 1&0&1 \\ \frac{1+\beta Y}{G_{0}}&0&\frac{1+\beta Y}{Y-c(1+\beta Y)}(F-\frac{\mu c}{G_{0}^{2}}) \end{array}\right)\left(\begin{array}{ccc} 0&0 \\ 1&0 \\ 0&1 \end{array}\right)\\ &=&\left(\begin{array}{ccc} 1&0 \\ 0&1 \\ 0&\frac{1+\beta Y}{Y-c(1+\beta Y)}(F-\frac{\mu c}{G_{0}^{2}}) \end{array}\right).\end{eqnarray*} $

Then, we get

$ \begin{eqnarray*}&&f_{1y_{1}}(M, \mu)=D_{N}f_{1}(M, \mu)D_{y_{1}}\psi(\bar{Z}, \mu)=r_{1}-Y, \\ &&f_{1y_{2}}(M, \mu)=D_{N}f_{1}(M, \mu)D_{y_{2}}\psi(\bar{Z}, \mu)=-X, \\ &&f_{2y_{1}}(M, \mu)=D_{N}f_{2}(M, \mu)D_{y_{1}}\psi(\bar{Z}, \mu)=\frac{Y^{2}}{X^{2}}, \\ &&f_{2y_{2}}(M, \mu)=D_{N}f_{2}(M, \mu)D_{y_{2}}\psi(\bar{Z}, \mu)=\frac{-Y}{Y-c(1+\beta Y)}(F-\frac{\mu c}{G_{0}^{2}})+r_{2}-\frac{2Y}{X}+F. \end{eqnarray*} $

Thus, we have

$ \begin{eqnarray*} &&D_{N}f_{1y_{1}}(M, \mu)=(0, -1, 0), \\ &&D_{N}f_{1y_{2}}(M, \mu)=(-1, 0, 0), \\ &&D_{N}f_{2y_{1}}(M, \mu)=(-\frac{2Y^{2}}{X^{3}}, \frac{2Y}{X^{2}}, 0), \\ &&D_{N}f_{2y_{2}}(M, \mu)=(\frac{2Y}{X^{2}}, \cdots, \frac{c}{[Y-c(1+\beta Y)](1+\beta Y)}). \end{eqnarray*} $

Then, we obtain

$ \begin{eqnarray*} &&f_{1y_{1}y_{1}}(M, \mu)=D_{N}f_{1y_{1}}(M, \mu)D_{y_{1}}\psi(\bar{Z}, \mu)=0, \\ &&f_{1y_{1}y_{2}}(M, \mu)=D_{N}f_{1y_{1}}(M, \mu)D_{y_{2}}\psi(\bar{Z}, \mu)=-1, \\ &&f_{1y_{2}y_{2}}(M, \mu)=D_{N}f_{1y_{2}}(M, \mu)D_{y_{2}}\psi(\bar{Z}, \mu)=0, \\ &&f_{2y_{1}y_{1}}(M, \mu)=D_{N}f_{2y_{1}}(M, \mu)D_{y_{1}}\psi(\bar{Z}, \mu)=-\frac{2Y^{2}}{X^{3}}, \\ &&f_{2y_{1}y_{2}}(M, \mu)=D_{N}f_{2y_{1}}(M, \mu)D_{y_{2}}\psi(\bar{Z}, \mu)=\frac{2Y}{X^{2}}, \\ &&f_{2y_{2}y_{2}}(M, \mu)=D_{N}f_{2y_{2}}(M, \mu)D_{y_{2}}\psi(\bar{Z}, \mu)=\frac{2c}{G_{0}^{2}}(F+\beta F-\frac{\mu c}{G_{0}^{2}}-\frac{\mu}{G_{0}^{2}})-\frac{2}{X}. \end{eqnarray*} $

Substituting $M_{0}, \mu_{0}$ into above, we have

$ \begin{eqnarray*} &&f_{1y_{1}}(M_{0}, \mu_{0})=0, ~~~f_{1y_{2}}(M_{0}, \mu_{0})=-X_{0}, ~~~f_{2y_{1}}(M_{0}, \mu_{0})=\frac{Y_{0}^{2}}{X_{0}^{2}}, \\ &&f_{1y_{1}y_{2}}(M_{0}, \mu_{0})=-1, ~~~f_{2y_{1}y_{1}}(M_{0}, \mu_{0})=-\frac{2Y_{0}^{2}}{X_{0}^{3}}, ~~~f_{2y_{1}y_{2}}(M_{0}, \mu_{0})=\frac{2Y_{0}}{X_{0}^{2}}, \\ &&f_{2y_{2}y_{2}}(M_{0}, \mu_{0})=\frac{2c\mu}{G_{0}^{4}}(\beta c-1)+\frac{2}{X_{0}}. \end{eqnarray*} $

Now, we get

$ \begin{eqnarray*} &&D_{N}f_{2y_{1}y_{1}}(M, \mu)=(\frac{6Y^{2}}{X^{4}}, -\frac{4Y}{X^{3}}, 0), \\ &&D_{N}f_{2y_{1}y_{2}}(M, \mu)=(-\frac{4Y}{X^{3}}, \frac{2}{X^{2}}, 0), \\ &&D_{N}f_{2y_{2}y_{1}}(M, \mu)={\frac{2}{X^{2}}, \frac{G_{0}}{1+\beta Y}\cdot(\frac{2\mu}{G_{0}^{2}}-F(1-c\beta)), -\frac{c(1+\beta)}{G_{0}^{2}(1+\beta Y)^{2}})}. \end{eqnarray*} $

Finally, we obtain

$ \begin{eqnarray*} &&f_{2y_{1}y_{1}y_{1}}=D_{N}f_{2y_{1}y_{1}}(M, \mu)D_{y_{1}\psi}(\bar{Z}, \mu)=\frac{6Y^{2}}{X^{4}}, \\ &&f_{2y_{1}y_{1}y_{2}}=D_{N}f_{2y_{1}y_{1}}(M, \mu)D_{y_{2}\psi}(\bar{Z}, \mu)=-\frac{4Y}{X^{3}}, \\ &&f_{2y_{2}y_{2}y_{1}}=D_{N}f_{2y_{2}y_{2}}(M, \mu)D_{y_{1}\psi}(\bar{Z}, \mu)=\frac{2}{X^{2}}, \\ &&f_{2y_{2}y_{2}y_{2}}=D_{N}f_{2y_{2}y_{2}}(M, \mu)D_{y_{2}\psi}(\bar{Z}, \mu)=\frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu G_{0}}{1+\beta Y_{0}}-\mu c). \end{eqnarray*} $

Thus, we have eq. (3.3)

$ \begin{align} & \overset{.}{\mathop{{{y}_{1}}}}\,=-{{X}_{0}}{{y}_{2}}-{{y}_{1}}{{y}_{2}}, \\ & \overset{.}{\mathop{{{y}_{2}}}}\,=\frac{Y_{0}^{2}}{X_{0}^{2}}{{y}_{1}}-\frac{Y_{0}^{2}}{X_{0}^{3}}y_{1}^{2}+\frac{2{{Y}_{0}}}{X_{0}^{2}}{{y}_{1}}{{y}_{2}}+[\frac{c\mu }{G_{0}^{4}}(\beta c-1)+\frac{1}{{{X}_{0}}}]y_{2}^{2} \\ & \ \ \ \ \ \ +\frac{Y_{0}^{2}}{X_{0}^{4}}y_{1}^{3}-\frac{2Y}{{{X}^{3}}}y_{1}^{2}{{y}_{2}}+\frac{1}{{{X}^{2}}}{{y}_{1}}y_{2}^{2} \\ & \ \ \ \ \ \ \ +\frac{1}{6}\cdot \frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu {{G}_{0}}}{1+\beta {{Y}_{0}}}-\mu c)y_{2}^{3}+o(|\bar{Z}{{|}^{4}}). \\ \end{align} $

(3.3)

Comparing with the normal form (3.2), we chosse the nonsingular matrix

$ \begin{aligned} N=\left(\begin{array}{cc} X_{0}\sqrt{X_{0}}&0 \\ 0&Y_{0} \end{array}\right), \end{aligned} $

then we use the linear transformation $H=N\bar{Z}$, noticing $\omega_{0}=\frac{Y_{0}}{\sqrt{X_{0}}}$, we derive the normal form as follows

$ \begin{align} & \overset{.}{\mathop{{{u}_{1}}}}\,=-{{\omega }_{0}}{{u}_{2}}-{{Y}_{0}}{{u}_{1}}{{u}_{2}}, \\ & \overset{.}{\mathop{{{u}_{2}}}}\,={{\omega }_{0}}{{u}_{1}}-{{Y}_{0}}u_{1}^{2}+\frac{2{{Y}_{0}}}{\sqrt{{{X}_{0}}}}{{u}_{1}}{{u}_{2}}+[\frac{c\mu }{G_{0}^{4}}(\beta c-1)+\frac{1}{{{X}_{0}}}]{{Y}_{0}}u_{2}^{2} \\ & \ \ \ \ \ \ +{{Y}_{0}}\sqrt{{{X}_{0}}}u_{1}^{3}-2{{Y}_{0}}u_{1}^{2}{{u}_{2}}+\frac{{{Y}_{0}}}{\sqrt{{{X}_{0}}}}{{u}_{1}}u_{2}^{2} \\ & \ \ \ \ \ \ \ \ +\frac{1}{6}\cdot \frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu {{G}_{0}}}{1+\beta {{Y}_{0}}}-\mu c)Y_{0}^{2}u_{2}^{3}+o(|\bar{Z}{{|}^{4}}), \\ \end{align} $

(3.4)

where $H=(u_{1}, u_{2})^{T}$. Then,

$ \begin{eqnarray*}&&a_{11}^{1}=a_{22}^{1}=0, \\ &&a_{12}^{1}=-Y_{0}~~~~~~~a_{11}^{2}=-Y_{0}, ~~~~~~~a_{12}^{2}=\frac{2Y_{0}}{\sqrt{X_{0}}}, \\ &&a_{22}^{2}=[\frac{c\mu}{G_{0}^{4}}(\beta c-1)+\frac{1}{X_{0}}]Y_{0}, \\ &&a_{111}^{1}=a_{122}^{1}=0, \\ &&a_{111}^{2}=Y_{0}\sqrt{X_{0}}, ~~~~~~a_{112}^{2}=-2Y_{0}, ~~~~~~a_{122}^{2}=\frac{Y_{0}}{\sqrt{X_{0}}}, \\ &&a_{222}^{2}=\frac{1}{6}\cdot\frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu G_{0}}{1+\beta Y_{0}}-\mu c)Y_{0}^{2}.\end{eqnarray*} $

According to the Hopf bifurcation theorem in [19], now we only need to calculate the value of $a$

$ \begin{eqnarray*} 16a&=&[a_{12}^{1}(a_{11}^{1}+a_{22}^{1})-a_{12}^{2}(a_{11}^{2}+a_{22}^{2})-a_{11}^{1}a_{11}^{2}+a_{22}^{1}a_{22}^{2}]/\omega+(a_{111}^{1}+a_{122}^{1}+a_{112}^{2}+a_{222}^{2})\\ &=&(-a_{12}^{2}a_{11}^{2}-a_{12}^{2}a_{22}^{2})/\omega+a_{112}^{2}+a_{222}^{2}\\ &=&-\frac{2Y_{0}}{X_{0}}+2Y_{0}\frac{c\mu}{G_{0}^{4}}(1-\beta +\frac{1}{6}\cdot\frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu G_{0}}{1+\beta Y_{0}}-\mu c)Y_{0}^{2}.\end{eqnarray*} $

Next, there are two cases should be discussed. That is $a>0$ and $a<0$. Based on the Hopf bifurcation theorem in [19], we obtain Theorem 3.1.

Theorem 3.1  For the system (2.2), there exist an $\varepsilon>0$ and two small enough neighborhoods $P_{1}$ and $P_{2}$ of $\chi_{0}(\mu)$, where $P_{1}\subset P_{2}$.

(ⅰ) If

$ \begin{aligned} 2Y_{0}\frac{c\mu}{G_{0}^{4}}(1-\beta c) +\frac{1}{6}\cdot\frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu G_{0}}{1+\beta Y_{0}}-\mu c)Y_{0}^{2}>\frac{2Y_{0}}{X_{0}}, \end{aligned} $

then

(1) when $\mu_{0}<\mu<\mu_{0}+\varepsilon$, $\chi_{0}(\mu)$ is unstable, and repels all the points in $P_{2}$;

(2) when $\mu_{0}-\varepsilon<\mu<\mu_{0}$, there exist at least one periodic solution in $\bar{P_{1}}$, which is the closure of $P_{1}$, one of them repel all the points in $\bar{P_{1}}\backslash \{\chi_{0}(\mu)\}$, and also have another periodic solution (may be the same that) repels all the points in $P_{2} \backslash \bar{P_{1}}$, and $\chi_{0}(\mu)$ is locally asymptotically stable.

(ⅱ) If

$ \begin{aligned} 2Y_{0}\frac{c\mu}{G_{0}^{4}}(1-\beta c) +\frac{1}{6}\cdot\frac{1-\beta c}{G_{0}^{2}}(\frac{2\mu G_{0}}{1+\beta Y_{0}}-\mu c)Y_{0}^{2}<\frac{2Y_{0}}{X_{0}}, \end{aligned} $

then

(1) when $\mu_{0}-\varepsilon<\mu<\mu_{0}, \chi_{0}(\mu)$ is locally asymptotically stable, and repels all the points in $P_{2}$;

(2) when $\mu_{0}<\mu<\mu_{0}+\varepsilon$, there exist at least one periodic solution in $\bar{P_{1}}$, one of them repel all the points in $\bar{P_{1}}\backslash\{\chi_{0}(\mu)\}$, and also have another periodic solution (may be the same that) repels all the points in $P_{2}\backslash \bar{P_{1}}$, and $\chi_{0}(\mu)$ is unstable.

Proof  The Theorem 3.1 can be similarly proved as the Hopf bifurcation theorem in [19], so we omit the process here.

In this section, we give a numerical example of the system (1.4) with the parameters $r_{1}=3, r_{2}=1, c=1, \beta=0.195$, then, the system (1.4) becomes

$ \left\{ \begin{array}{*{35}{l}} \dot{x}=x(3-y), \\ \dot{y}=y(1-\frac{y}{x})-E\frac{y}{1+0.195y}, \\ 0=E(\frac{y}{1+0.195y}-q)-\mu . \\ \end{array} \right. $

(4.1)

By simple computing, the only positive equilibrium point of above system is

$ \chi(\mu_{0})=(4.7578, 3, 0.5856), $

and the Hopf bifurcation value $\mu_{0}=\frac{r_{2}G_{0}^{2}}{2G_{0}+c}=\frac{2.0002225}{3.83}$.

Therefore, by Theorem 3.1, we can easily show that the positive equilibrium point $\chi_{0}(\mu)$ of the system (4.1) is locally asymptitically stable when $\mu=0.505<\mu_{0}$ as is illustrated by computer simulations in Fig. 1; periodic solutions occur from $\chi_{0}(\mu)$ when $\mu=0.5195<\mu_{0}$ as is illustrated in Fig. 2; the positive equilibrium point $\chi_{0}(\mu)$ of the system (4.1) is unstable when $\mu=0.535>\mu_{0}$ as is illustrated in Fig. 3.

Based on the above inference and calculation, we have found that economic effect will influence the stability of differential-algebraic biological economic system. For instance, according to those statistics and graphs, if people fix the economic index at a high level, over the bifurcation value of Hopf-bifurcation, the system will become unstable that means people have destroyed the economic balance even led to the extinction of ecologic species. Therefore, with an aim to realize the harmonious sustainable development co-existence between man and nature, we should not seek economic effect blindly and control it within a certain limit, such as less than bifurcation value.

In addition, we can make some improvements in our model. For example, we do not consider the influence of time delays and double harvesting that is, human harvesting will harvest predator and prey at the same time. So it is necessary for us to go on with our research in these aspects in the future.