Bifurcation analysis of a special delayed predator-prey model with herd behavior and prey harvesting

1.   Introduction

The dynamic relationship between predator and prey has long been and will continue to be one of the dominant themes in both biology and mathematical biology due to its universal existence. In the description of the dynamics interactions, a crucial element of all models is the classic definition of functional response. Many kinds of predator-prey models with Holling type ^[1,2,3], Leslie-Gower type ^[4,5], and Beddington-DeAngelis type ^[6,7], etc. have been investigated extensively by scholars. However, as argued by the authors in ^[8], it is more appropriate to model the response function of prey which exhibits herd behavior in terms of the square root of the prey population in realistic world. For example, this may be entirely fitting for herbivores on a large savanna and their large predators. Since the prey population exhibits a highly socialized behavior and lives in herds, that is, the weaker individuals are kept at the center of their herd for defensive purpose, the interaction terms of predator-prey system use the square root of the prey population rather than the standard mass-action term, which properly accounts for the assumption that the interactions occur along the boundary of the population. Predator-prey systems with such functional response have attracted much attention (see ^{[8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]}). Recently, Yuan et al. ^[10] considered a predator-prey system as following:

where $-sY^2$ represents the quadratic mortality for predator population. They investigated the spatial dynamics of system (1.1) with Neumann boundary conditions and obtained Turing bifurcation. Considering the fact that there always exists a time delay in the conversion of the biomass of prey to that of predator in system (1.1), Xu and Yuan ^[11] studied the local stability and the existence of Hopf bifurcation of such system. Kooi and Venturino ^[14] studied ecoepidemic predator-prey model with prey herd behavior. Liu et al. ^[18] studied stationary distribution and extinction of a stochastic predator-prey model with herd behavior. Some authors (Meng and Wang ^[19], Yang ^[20,21], Souna et al. ^[22]) obtained dynamical behaviors of such delayed diffusive predator-prey model with herd schooling behavior, such as Hopf bifurcation, the steady-state bifurcation, diffusion driven instability and Turing-Hopf bifurcation. Bentout et al. ^[23] considered prey escaping from prey herd on three species fractional predator-prey interaction model, and showed the considered system undergoes an interesting behavior.

Time delays of one type or another have been incorporated into mathematical models of population dynamics due to maturation time, capturing time or other reasons. In general, delay differential equations exhibit much more complex dynamics than ordinary differential equations since a time delay could cause the stable equilibrium to become unstable and even cause the population to fluctuate. Many authors have devoted to investigating time delay effecting on the dynamics of biological systems, and obtained some results (see ^{[11,24,25,26,27,28,29,30]}).

In 1954, Gordon ^[31] proposed the economic theory of a common-property resource, which studies the effect of the harvest effort on the ecosystem from an economic perspective. In that reference, an algebraic equation is proposed to investigate the economic interest of the yield of the harvest effort, which takes form as follows: $\text{Net Economic Revenue (NER)} = \text{Total Revenue (TR)-Total Cost (TC)}$. Let $E(t)$ and $Y(t)$ represent the harvest effort and the density of harvested population, respectively. $TR = \omega E(t)Y(t)$ and $TC = c E(t)$, where $\omega$ represents unit price of harvested population, $c$ represents the cost of harvest effort. Associated with the above system, an algebraic equation, which considers the economic interest $v$ of the harvest effort on the harvested population, is established as follows:

In daily life, since economic profit is a very important factor for governments, merchants and even every citizen, it is necessary to research singular economic systems, which can be described by differential-algebraic equations or difference-algebraic equations. Luenberger ^[32] studied the dynamic invest and produce in the economic system by applying the descriptor system. As far as biological systems are concerned, the related research results are few. Recently, Zhang et al. ^[33] first proposed a class of singular biological economic systems, which are established by several differential equations (or several difference equations) and an algebraic equation, and gained some results on those systems, such as the stabilities, bifurcations and chaos. Meng and Wu ^[34] investigated the existence of Hopf bifurcation in a differential-algebraic predator-prey model with nonlinear harvesting. More results can be found in the references ^{[35,36,37,38,39]}.

Based on the above discussions, a class of singular biological economic predator-prey system with square root functional response and two delays have received surprisingly little attention in the literature. Thus, based on the previous system (1.1) and Gordon theory ^[31], we establish the following predator-prey system consisting of two differential equations and an algebraic equation as follows:

where $X, Y$ and $E$ denote the prey, predator and the harvest effort at the time $s$, respectively. $\tau_1$ denotes the feedback time delay of prey species to the growth of species itself, and $\tau_2$ is the time delay which means the growth rate of predator species depends on the number of the prey species $\tau_2$ units of time earlier. The parameter $r$ and $N$ are the growth rate of the prey and its carrying capacity. The parameter $\alpha$ is the search efficiency of $Y$ for $X$, $\tilde{b}$ is biomass conversion or consumption rate, and $t_h$ is Y's average handling time of $X$. $-\tilde{d}Y^2$ represents the quadratic mortality for predator population. $\tilde{p}$ and $\tilde{c}$ represent unit price of harvested population and the cost of harvest effort. $\tilde{m}$ denotes harvest interest of the harvest effort on the harvested population. So we will study the dynamics of such system with harvesting for the prey by selecting different values of two delays in this paper.

The rest of paper is organized as follows. In the next section, with the help of new norm formal method ^[40,41], some conditions of the existence of the positive equilibrium, local stability and the existence of Hopf bifurcation are obtained. By using the normal form method and center manifold theory due to Hassard et al.^[42], some explicit formulas determining the direction and stability of periodic solutions bifurcating from Hopf bifurcations are showed in the Section 3. The control of Hopf bifurcations is given in the Section 4. The optimal tax policy of system is investigated baased on Pontryagin's maximum principle in the Section 5. To support our theoretical predictions, some numerical simulations are included in the Section 6. A brief discussion is also given in the last section.

2.   Stability analysis and existence of Hopf bifurcation

It is important to make some simplifying assumptions to discern the basic dynamics and to make the analysis more tractable. In order to simplify the system (1.2), we use the following dimensionless transformations:

Thus system (1.2) can be rewritten as following:

where $a = t_h \alpha \sqrt{N}, b = \frac{\alpha \sqrt{N}\tilde{b}}{r}, c = \tilde{c}, d = \frac{\sqrt{N}\tilde{d}}{r}, m = \frac{\tilde{m}}{r}, p = \tilde{p}N$.

In the papers ^[8,9,10,11], the authors used the simplifying assumption that $a = 0$, which implies that the average handling time is zero. In line with the work in ^[8,9,10,11], we also assume that $a = 0$. Thus, system (2.1) takes the form:

The initial conditions of system (2.2) have the form

In the next section, we will investigate the local stability of the positive equilibrium and the existence of Hopf bifurcation.

From the viewpoint of biological interpretation of system, we only consider the positive equilibrium. In this section, by choosing $\tau_1$ and $\tau_2$ as the bifurcation parameters and analyzing the corresponding characteristic equation of linearized system, we investigate the local stability of the positive equilibrium and the effects of the time delays on the dynamics of system (2.2).

2.1.   Existence of the positive equilibrium

There exists a positive equilibrium $P(x, y, e)$ of system (2.2), where the values $x, y$ and $e$ satisfy the following equations:

From Eq (2.3b) and (2.3c), we can obtain $y = \frac{b}{d}\sqrt{x}$ and $e = \frac{m}{px-c}$, respectively. Substituting the above values into Eq (2.3a), we can obtain that $x$ satisfies the following equation

It is obvious that system Eq (2.4) has positive roots if the condition $\Delta = (\frac{b}{d}-1-\frac{c}{p})^2-4(\frac{c}{p}+\frac{m}{p}-\frac{bc}{dp}) \geq 0$. If $\Delta = 0$, then Eq (2.4) has a unique positive root, defined $x^* = \frac{cd+pd-bp}{2dp}$. That is, system (2.2) has only one positive equilibrium $P^*(x^*, y^*, e^*)$ if the assumption

(H1): $p > c+\frac{bp}{d}$

holds, here $y^* = \frac{b}{d}\sqrt{x^*}$ and $e^* = \frac{m}{px^*-c}$. When $\Delta > 0$, Eq (2.4) has two positive roots, defined $x_1^*$ and $x_2^*$; that is, system (2.2) has two positive equilibria $P_i^*(x_i^*, y_i^*, e_i^*)\, i = 1, 2$, here $y_i^* = \frac{b}{d}\sqrt{x_i^*}, e_i^* = \frac{m}{px_i^*-c}\, (i = 1, 2)$ and $x_1^* = \frac{cd+pd-bp+\sqrt{\Delta}}{2dp} > x^*, x_2^* = \frac{cd+pd-bp-\sqrt{\Delta}}{2dp} < x^*$.

Remark 1: We can obtain some exact expressions of the positive equilibriums when $a = 0$ of system (2.1), which will be used to discuss the dynamical behaviours of system. If $a \leq 0$ of system (2.1), then it is difficult to investigate the properties of such system.

Remark 2: In the following, we only investigate the stability and bifurcation of system (2.2) at the unique positive equilibrium $P^*(x^*, y^*, e^*)$. Of course, we also study the dynamics of system (2.2) at the equilibrium $P_i^*(x_i^*, y_i^*, e_i^*)\, (i = 1, 2)$ by the same way.

2.2.   Stability and existence of Hopf bifurcation

For simplicity, let

where $X = (x, y, e)^T$.

In order to analyze the local stability of the positive equilibrium of the system (2.2), we first use the linear transformation $X(t) = QN(t)$, where

Then, we get $D_Yg(P^*)Q = (0, 0, px^*-c)$ and $u(t) = x(t), v(t) = y(t), E(t) = \frac{pe^*}{px^*-c}x(t)+e(t)$. Substituting the latter into system (2.2), we can obtain

In order to derive the formula determining the properties of the positive equilibrium of the system (2.2), we consider the local parametric $\Psi$ of the third equation of system (2.2) as in ^[40,41], which is given as follows:

where $u_0 = \begin{pmatrix} 1 & 0\\ 0 & 1\\ 0 & 0 \end{pmatrix}, \, v_0 = \begin{pmatrix} 0\\ 0 \\ 1 \end{pmatrix}, \, Z(t) = \begin{pmatrix} y_1(t)\\ y_2(t) \end{pmatrix}, \, N_0 = \begin{pmatrix} u^* \\ v^* \\ E^* \end{pmatrix}$, $h(Z(t)) = (0, 0, h_3(y_1(t), y_2(t))^T: R^2 \rightarrow R$ is a smooth mapping.

Thus, we obtain the following parametric system of the system (2.5):

because $g(\Psi(Z(t))) = 0$. Now, we give the linearized system of parametric system (2.6) at $(0, 0)$ as follows:

The corresponding characteristic equation of the linearized system of the above system is:

where $A = 1-\frac{b}{2d}-u^*+\frac{mc}{(pu^*-c)^2}$.

In order to investigate the distribution of roots of the transcendental equation (2.8), we introduce the following result which was proved by Ruan and Wei ^[24] using Rouch$\acute{\text{e}}$'s theorem.

Lemma 2.1. Consider the exponential polynomial

where $\tau_{i}\geq0 \, (i = 1, 2, \cdots, m)$, $p_{j}^{i}\, (i = 1, 2, \cdots, m;j = 1, 2, \cdots, n)$ are constants. As $(\tau_{1}, \tau_{2}, \cdots, \tau_{m})$ vary, the sum of the order of the zeros of $P(\lambda, e^{-\lambda\tau_{1}}, \cdots, e^{- \lambda\tau_{m}})$ on the open right half plane can change only if a zero appears on or crosses the imaginary axis.

It is obvious that $\lambda = 0$ is not a root of Eq (2.8). When there is no delay, that is, $\tau_1 = \tau_2 = 0$, Eq (2.8) is rewritten as follows:

Thus, if $u^*+\frac{b}{2d}-A < 0$, Eq (2.9) has at least one positive root, which implies that system (2.2) without time delay is unstable. If the condition

(H2): $1+\frac{mc}{(pu^*-c)^2} < \min\big\{2u^*+\frac{b}{d}, 2u^*+\frac{b}{2d}+b\sqrt{u^*}\big\}$

holds, the two roots of Eq (2.9) have always negative real parts. Thus, we have the following result.

Lemma 2.2. Assume that (H1) and (H2) hold, then the two roots of Eq ${\rm(2.9)}$ have always negative real parts, that is, the positive equilibrium of system ${\rm(2.2)}$ with $\tau_1 = \tau_2 = 0$ is locally asymptotically stable.

In what follows, we will discuss in detail the local stability of system around the positive equilibrium and the existence of Hopf bifurcations occurring at the positive equilibrium by selecting different values of $\tau_1$ and $\tau_2$.

Case one: $\tau_1 = 0$ and $\tau_2\neq 0$.

In this case, Eq (2.8) is translated into the following form

where $A_{11} = b\sqrt{u^*}+u^*-A, A_{10} = b\sqrt{u^*}(u^*-A), B_{10} = \frac{b^2\sqrt{u^*}}{2d}$.

Now, we can obtain the same result as Lemma 2.2 when $\tau_2 = 0$. For $\tau_2 > 0$, let $\omega {\rm i} \, (\omega > 0)$ be a root of Eq (2.10), then $\omega$ satisfies the following equation:

Separating the real and imaginary parts, we obtain

which yields that

Further, if denote $v = \omega^2$, we can obtain $\Phi_1(\nu) = \nu^2+(A_{11}^2-2A_{10})\nu+A_{10}^2-B_{10}^2$. It is obvious that Eq (2.12) dose not have positive root when $A_{10}^2-B_{10}^2 > 0$ and $A_{11}^2-2A_{10} \geq 0$ or $\Delta_1 \triangleq (A_{11}^2-2A_{10})^2-4(A_{11}^2-2A_{10}) < 0$.

If the following assumption

(H3): $A_{10}^2-B_{10}^2 \leq 0$

holds, then Eq (2.12) has a positive root $\omega_{2}^{+}$, here $\omega_{2}^{+} = \sqrt{\frac{2A_{10}-A_{11}^2+\Delta_1}{2}}$.

If the following assumption

(H4): $A_{10}^2-B_{10}^2 > 0$, $A_{11}^2-2A_{10} < 0$ and $\Delta_1 > 0$

holds, then Eq (2.12) has two positive roots $\omega_{2}^{+}$ and $\omega_{2}^{-}$, $\omega_{2}^{\pm} = \sqrt{\frac{2A_{10}-A_{11}^2\pm \Delta_1}{2}}$.

From Eq (2.11), if (H3) holds and denote

then $\pm {\rm i}\omega_{2}^{+}$ are a pair of purely imaginary root of Eq (2.10) with $\tau_2 = \tau_{2j}^{+}$.

Define $\tau _{20} = {\rm min}\{\tau_{2j}^{+}, j = 0, 1, 2, \cdots\}$. Further, if denote $\lambda(\tau) = \alpha(\tau)+{\rm i}\omega(\tau)$ be the root of (2.10), then it satisfies $\alpha(\tau_{20}) = 0, \omega(\tau_{20}) = \omega_{20}$. Next, we will investigate whether the transversality condition is satisfied. Differentiating the two sides of Eq (2.10) with respect to $\tau_2$, we get

For simplify, we define $\omega_{2}^{+}$ as $\omega$ and $\tau_{2j}^{+}$ as $\tau_2$, and obtain

Thus, we have the following conclusion.

Lemma 2.3. Assume that (H1) holds, then the transversality condition is satisfied. That is,

From Eq (2.11), if (H4) holds and we denote

then $\pm {\rm i}\omega_{2}^{+}$ (resp. $\pm i\omega_{2}^{-}$) are a pair of purely imaginary root of Eq (2.10) with $\tau_2 = \tau_{2j}^{+}$ (resp. $\tau_2 = \tau_{2j}^{-}$). Further, we can check that there exists a positive integral $k > 0$ such that $\tau_{2_{k+1}}^+ > \tau_{2_k}^-$ and $0 < \tau_{2_0}^+ < \tau_{2_0}^- < \tau_{2_1}^+ < \tau_{2_1}^-\cdots\tau_{2_{k-1}}^- < \tau_{2_k}^+ < \tau_{2_k}^- < \tau_{2_{k+1}}^+$.

Differentiating the both sides of Eq (2.10) with respect to $\tau_2$, then straightforward but tedious calculations lead to the following conclusion.

Lemma 2.4. Assume (H1) and $\Phi_1'(\nu)\neq 0$ hold, then the transversality condition

are satisfied, and $(\frac{{\rm d}\lambda}{{\rm d}\tau_2})^{-1}$ and $\Phi_1'(\nu)$ have the same sign.

Summarizing the above lemmas and applying Lemma 2.4 to Eq (2.10), we have the following theorem.

Theorem 2.1. Suppose that (H1) and (H2) are true.

(1) If (H3) holds, system $\rm(2.2)$ is locally asymptotically stable for $\tau_2 \in [0, \tau_{20})$ and unstable for $\tau_2 > \tau_{20}$. That is, system $\rm(2.2)$ undergoes a Hopf bifurcation at the positive equilibrium $P^*(x^*, y^*, e^*)$ when $\tau_2 = \tau_{2j}^+, \, j = 0, 1, 2, \cdots$.

(2) If (H4) and $\Phi_1'(\nu) \neq 0$, then there exists $N\in\mathbb{N}$, when $\tau_2 \in (0, \tau_{2_0}^+)\cup (\tau_{2_0}^-, \tau_{2_1}^+)\cup \cdots (\tau_{2_N}^-, \tau_{2_{N+1}}^+)$, all the roots of Eq ${\rm(2.10)}$ have negative real part, thus system $\rm(2.2)$ is locally asymptotically stable; when $\tau_2 \in (\tau_{2_0}^+, \tau_{2_0}^-)\cup (\tau_{2_1}^+, \tau_{2_1}^-)\cup \cdots \cup (\tau_{2_N}^+, \tau_{2_N}^-)\cup(\tau_{2_{N+1}}^+, +\infty)$, Eq ${\rm(2.10)}$ has at least one root with positive real part, thus system $\rm(2.2)$ is unstable. Hopf bifurcation occurs near the positive equilibrium $P^*$ when $\tau_2 = \tau_{2j}^{\pm}, \, j = 0, 1, 2, \cdots$.

Case two: $\tau_1 \neq 0$ and $\tau_2 = 0$.

Eq (2.8) takes the form

where $A_{21} = b\sqrt{u^*}-A, A_{20} = \frac{b\sqrt{u^*}}{2d}-Ab\sqrt{u^*}, B_{21} = u^*, B_{20} = bu^*\sqrt{u^*}$.

Let $\omega {\rm i} \, (\omega > 0)$ be the root of Eq (2.13). We obtain

which gives

Further, if denote $v = \omega^2$, we can obtain $\Phi_2(\nu) = \nu^2+\varrho \nu+A_{20}^2-B_{20}^2$, here $\varrho = A_{21}^2-B_{21}^2-2A_{20}$. It is obvious that Eq (2.15) dose not have positive root when $A_{21}^2-B_{21}^2-2A_{20} > 0$ and $A_{20}^2-B_{20}^2 > 0$ or $\Delta_2 \triangleq (A_{21}^2-B_{21}^2-2A_{20})^2-4(A_{20}^2-B_{20}^2) < 0$. We give two assumptions as following:

(H5): $A_{20}^2-B_{20}^2 \leq 0$;

(H6): $A_{21}^2-B_{21}^2-2A_{20} < 0$, $A_{20}^2-B_{20}^2 > 0$ and $\Delta_2 > 0$.

If the assumption (H5) holds, then Eq (2.15) has only one positive root, defined $\omega_1^+$. If the assumption (H6) holds, then Eq (2.15) has two positive roots $\omega_1^{\pm} = \sqrt{\frac{1}{2}[-\varrho\pm \Delta_2]}$, and $\omega_1^{+} > \omega_1^{-} > 0$.

From Eq (2.14), if denote

then $\pm {\rm i}\omega_{1}^{+}$ (resp. $\pm {\rm i}\omega_{1}^{-}$) are a pair of purely imaginary root of Eq (2.13) with $\tau_1 = \tau_{1_j}^+$ (resp. $\tau_1 = \tau_{1_j}^-$).

When (H6) holds, we define $\tau _{10} = \min\{\tau_{1_j}^{\pm}, \, j = 0, 1, 2, \cdots\}$. Let $\lambda(\tau) = \xi(\tau)+{\rm i}\omega(\tau)$ be the root of Eq (2.13) near $\tau = \tau_{10}$ satisfying $\xi(\tau_{10}) = 0$, $\omega(\tau_{10}) = \omega_{10}$. Further, if (H2) holds for $\tau_1 = 0$, $k$ switches from stability to instability to stability occur when the parameters are such that $0 < \tau_{1_0}^+ < \tau_{1_0}^- < \tau_{1_1}^+ < \cdots < \tau_{1_k}^+ < \tau_{1_k}^- < \tau_{1_{k+1}}^+$.

The following transversality condition is true.

Lemma 2.5. Suppose (H6) and $\Phi_2'(\omega^2) \neq 0$ hold, then the following transversality condition

are satisfied, $(\frac{{\rm d}\lambda}{{\rm d}\tau_1})^{-1}$ and $\Phi_2'(\nu)$ have the same sign.

Proof. Taking the derivative of $\lambda$ with respect to $\tau_1$ in Eq (2.13), we obtain

For simplify, we define $\omega_1^{\pm}$ as $\omega$ and $\tau_{1j}^{\pm}$ as $\tau_1$, and we can obtain

This completes the proof.

We can summarize the preceding results as the following theorem.

Theorem 2.2. For the system ${\rm(2.2)}$, if (H1) and (H2) hold, the following results are true.

(1) If Eq ${\rm(2.15)}$ does not have positive root, then system ${\rm(2.2)}$ is locally asymptotically stable for any $\tau_1 \geq 0$.

(2) If (H5) and $\Phi_2(\omega) \neq 0$ hold, then system ${\rm(2.2)}$ is locally asymptotically stable for $\tau_1 \in [0, \tau_{10})$ and unstable for $\tau_1 > \tau_{10}$. System ${\rm(2.2)}$ undergoes a Hopf bifurcation at the equilibrium $P^*$ when $\tau_1 = \tau_{1_j}^+ \, (j = 0, 1, 2, \cdots)$.

(3) If (H6) and $\Phi_2(\omega) \neq 0$, then there exists $N\in\mathbb{N}$, when $\tau_1 \in (0, \tau_{1_0}^+)\cup (\tau_{1_0}^-, \tau_{1_1}^+)\cup \cdots (\tau_{1_N}^-, \tau_{1_{N+1}}^+)$, all the roots of Eq ${\rm(2.15)}$ have negative real parts, thus system ${\rm(2.2)}$ is locally asymptotically stable; when $\tau_1 \in (\tau_{1_0}^+, \tau_{1_0}^-)\cup (\tau_{1_1}^+, \tau_{1_1}^-)\cup \cdots \cup (\tau_{1_N}^+, \tau_{1_N}^-)\cup(\tau_{1_{N+1}}^+, +\infty)$, Eq ${\rm(2.15)}$ has at least one root with positive real part, thus system ${\rm(2.2)}$ is unstable. Hopf bifurcation occurs near the positive equilibrium $P^*$ when $\tau_1 = \tau_{1_j}^{\pm}\, (j = 0, 1, 2, \cdots)$.

Remark 3 In case of $\tau_1\neq 0$ and $\tau_2 = 0$, whether (H6) or (H7) holds, we define the stable interval of time delay $\tau_1$ as $I_{1s} = [0, \tau_{10})$ (resp. $I_{1s} = [0, \tau_{1_0})\cup (\tau_{1_0}^-, \tau_{1_1}^+)\cup \cdots (\tau_{1_N}^-, \tau_{1_{N+1}}^+\;)$.

Case three: $\tau_1 > 0$ and $\tau_2 > 0$.

This case states that $\tau_2$ is regarded as a parameter and $\tau_1$ belongs to its stable interval. Let $\omega {\rm i} \, (\omega > 0)$ be a root of Eq (2.8) for any $\tau_1^*\in I_{1s}$, then we can obtain

which yields

Without loss of generality, suppose that

$\textbf{(H7)} \; \quad {\rm Eq}\; (2.17) \, \text{has at least finite positive roots}.$

If (H7) holds, the roots of Eq (2.17) defined by $\omega_1, \omega_2, \cdots, \omega_k$. For every fixed $\omega_i \, (i = 1, 2, \cdots, k)$, there exists a sequence

then $\pm {\rm i}\omega_{i}$ are a pair of purely imaginary root of (2.8) with $\tau_2 = \tau_{2j}^{(i)}$.

Define

Let $\lambda(\tau) = \xi(\tau)+{\rm i}\omega(\tau)$ be the root of Eq (2.8) near $\tau = \tau_{2*}$ satisfying $\xi(\tau_{2*}) = 0$, $\omega(\tau_{2*}) = \omega_{*}$. Further, we obtain the following lemma.

Lemma 2.6. Suppose (H7) and $P_RQ_R-P_IQ_I > 0$ hold, then the following transversality condition

is satisfied, where

Proof. Taking the derivative of $\lambda$ with respect to $\tau_2$ in Eq (2.8), we obtain

For simplify, defining $\omega_i$ as $\omega_*$ and $\tau_{2j}^{(i)}$ as $\tau_2$, we can obtain

If $P_RQ_R-P_IQ_I > 0$ holds, then this ends the proof.

Therefore, by the general Hopf bifurcation theorem for FDEs in Hale ^[25], we have the following results on stability and bifurcation of system (2.2).

Theorem 2.3. For system ${\rm(2.2)}$, suppose (H7) and $P_RQ_R-P_IQ_I > 0$ hold. Then system ${\rm(2.2)}$ is locally asymptotically stable when $\tau_2\in[0, \tau_{2*})$ and $\tau_1 \in I_{1s}$, and system ${\rm(2.2)}$ undergoes a Hopf bifurcation at the positive equilibrium $P^*$ when $\tau_2 = \tau_{2*}$. That is, system ${\rm(2.2)}$ has a branch of periodic solutions bifurcating from Hopf bifurcation near $\tau_2 = \tau_{2*}$.

3.   Direction and stability of the Hopf bifurcation

In this section, we shall study the direction of the Hopf bifurcation and stability of bifurcating periodic solution of system (2.2). The approach employed here is the normal form method and center manifold theorem due to Hassard et al. ^[42]. More precisely, we will compute the reduced system on the center manifold with the pair of conjugate complex, purely imaginary solution of the characteristic equation (2.8). By taking time delay $\tau_2$ for fixed value $\tau_{1*} \in I_{1s}$, we can determine the Hopf bifurcation direction to answer the question of whether the bifurcation branch of periodic solution exists locally for supercritical bifurcation or subcritical bifurcation.

Throughout this section, it is considered system (2.2) undergoes Hopf bifurcation at the positive equilibrium $P^*$ when $\tau_2 = \tau_{2*}, \tau_{1*} \in I_{1s}$. Without loss of generality, we assume that $\tau_{1*} < \tau_{2*}$. Let $x_1(t) = y_1(t)-y_1^{*}, x_2(t) = y_2(t)-y_2^{*}, t = t/\tau_2, \tau_2 = \tau_{2*}+\mu$ and $\bar{x}_i(t) = x_i(\tau_2t)$ and dropping the bars for simplification of notations, then system (2.2) is transformed into an functional differential equation in $C = C([-1, 0], R^2)$ as

where $x(t) = (x_1(t), x_2(t))^T\in R^2$, and $L_{\mu}:C \rightarrow R^2, F: R \times C \rightarrow R^2$ are given, respectively, by

and

where

and

here $\phi(\theta) = (\phi_1(\theta), \phi_2(\theta))^T\in C([-1, 0], \mathbb{R}^2)$,

System (3.1) turns to the linear problem

by the Riesz representation theorem. There exists $2\times 2$ matrix-valued function

such that

In fact, we can choose

Then Eq (3.2) is satisfied.

Next, for $\phi \in C([-1, 0], \mathbb{R}^2)$, we define the operator $A(\mu)$ as

and

Since $\tfrac{{\rm d}u_t}{{\rm d}\theta } = \tfrac{{\rm d}u_t}{{\rm d}t}$, then system (3.1) is equivalent to the following operator equation

where $u_t = u(t+\theta)$ for $\theta \in [-1, 0]$.

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

and a bilinear form

where $\eta(\theta) = \eta(\theta, 0)$. Then $A(0)$ and $A^{*}(0)$ are adjoint operators. From the above analysis, we obtain that $\pm {\rm i}\omega_{*}\tau_{2*}$ are the eigenvalues of $A(0)$ and therefore they are also the eigenvalues of $A^{*}(0)$. Let $q(\theta)$ be the eigenvector of $A(0)$ corresponding to ${\rm i}\omega_{*}\tau_{2*}$ and $q^{*}(s)$ be the eigenvector of $A^{*}(0)$ corresponding to $-{\rm i}\omega_{*}\tau_{2*}$, we have

By some simple computation, we can obtain

where

Then $< q^{*}, q > = 1$, and $< q^{*}, \bar{q} > = 0$.

In the remainder of this section, by using the same notations as in Hassard et al. ^[42], we compute the coordinates to describe the center manifold $C_0$ at $\mu = 0$ and obtain the following expressions:

From $g(z, \bar{z}) = \bar{q}^*(0)F_0(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}$, we obtain the coefficients determining the important quantities of the periodic solution:

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

where $E_1 = (E_1^{(1)}, E_1^{(2)})^T$ and $E_2 = (E_2^{(1)}, E_2^{(2)})^T$ are both two dimensional vectors and can be determined by

and

here

Furthermore, we can see that each $g_{ij}$ in (3.4) is determined by parameters and delays of system (2.2). Thus, we can compute the following quantities:

which determine the properties of bifurcation periodic solutions in the center manifold at the critical value $\tau_{2*}$. By the results of Hassard et al. ^[42], we have the following results.

Theorem 3.1. In ${\rm(3.5)}$, the following results hold.

(i) The sign of $\mu_{2}$ determines the direction of the Hopf bifurcation: If $\mu_{2} > 0$ (resp. $\mu_{2} < 0$), then the Hopf bifurcation is supercritical (resp. subcritical) and the bifurcation periodic solutions exist for $\tau_2 > \tau_{2*}$ ($\tau_2 < \tau_{2*}$).

(ii) The sign of $\beta_{2}$ determines the stability of the bifurcating periodic solutions: The bifurcation periodic solutions are stable (resp. unstable) if $\beta_{2} < 0$ (resp. $\beta_{2} > 0$).

(iii) The sign of $T_{2}$ determines the period of the bifurcating periodic solutions: The period increases (resp. decreases) if $T_{2} > 0$ (resp. $T_{2} < 0$).

4.   The control of Hopf bifurcation

The objective is to extend the range of parameter such that within this largest possible range system can remain its stable dynamical behavior, and yet the period-doubling bifurcations are being delayed and even being eliminated completely, or being stabilized to a desired unstable periodic orbit which is embedded in the chaotic attractor of the system.

The nonlinear controller for the population of the predator is taken as the form

here $\bar{k} = (k_{1}, k_{2})$ and $k_{1}, k_{2}$ are feedback gains.

Introducing the nonlinear controller into system (2.7), we can obtain

In what follows, we will only discuss one case $\tau_{1} = 0$ and $\tau_{2}\neq 0$. The other cases can be investigated similarly, so those are omitted.

In order to give the corresponding result, we give the following assumption

$\textbf{(H8)} \; k_{1} > \frac{A_{10}}{A-u^{*}} \, \text{and}\; \; (A_{10}-k_{1}(u^{*}-A))-(A_{11}-k_{1})^{2} < B_{10}e^{\frac{k_{1}-A_{11}}{2}\tau_{20}}.$

Theorem 4.1. If (H8) is true, then the control system (4.2) becomes stable when $\tau = \tau_{20}$.

Proof. When $\tau_{2} = \tau_{20}$, the associated characteristic transcendental equation of system (4.2) around $(0, 0)$ is

where $A_{11}, A_{10}$ and $B_{10}$ are the same to that in Eq (2.11)

Let $g_{1}(\lambda) = \lambda^{2}+(A_{11}-k_{1})\lambda+(A_{10}-k_{1}(u^{*}-A))$ and $g_{2}(\lambda) = -B_{10}e^{-\lambda \tau_{2}}$. In order to investigate the solution of system (4.3), we only need to discuss the solution of $g_{1}(\lambda) = g_{2}(\lambda)$. By a series of analysis, if $k_{1} > \frac{A_{10}}{A-u^{*}}$ and $(A_{10}-k_{1}(u^{*}-A))-(A_{11}-k_{1})^{2} < B_{10}e^{\frac{k_{1}-A_{11}}{2}\tau_{20}}$ hold, then we obtain that the whole roots of system (4.3) have the negative real part.

Remark 4 According to the theory of Hopf bifurcation, the nonlinear sate feedback method can be implemented to eliminate those unfavorable phenomena and stabilize the proposed marine ecosystems. Therefore, biological population can be controlled at states.

5.   Optimal tax policy

In this section, we will discuss the optimal tax policy of system (5.1). Our aim is to find the optimal tax policy that will maximize the benefits of fish population without extinction. We noticed that as the tax rate increased, the harvest decreased. Therefore, taxation plays a decisive role. When $\tau_{1} = \tau_{2} = 0$, system (2.2) becomes

However, we not only consider the net economic revenue of the harvested, but also consider social economic revenue of the regulatory agency based on Gupta et al. ^[43]. Therefore, we have

The present value $J(T)$ of a continuous time-stream of revenues is given by

where $\delta$ is the continuous annual discount rate which is fixed by harvesting agencies.

The control problem over an infinite time horizon is given by

We take advantage of the maximum principle to get the optimal solution of this problem. The associated Hamiltonian function is

where $\lambda_{i} = \lambda_{i}(t) (i = 1, 2)$ are adjoint variables corresponding to the variables $x, y$ respectively. The optimal control $T$ which maximizes $H$ must satisfy the following conditions:

It is assumed that the optimal solution does not occur at $T = T_{\text{max}}$ or $T = 0$. If $\frac{\partial H }{\partial T} = 0$, then we have $\lambda_{1}(t) = 0$. Based on Pontryagin's maximum principle ^[44], the adjoint variables must satisfy the following adjoint equations

In order to obtain singular optimal equilibrium solution, we use steady state equations

Thus, Eq (5.5) along with steady state Eq (5.6) gives

From Eq (5.7) we get

In order to obtain an optimal equilibrium, based on the method of reference ^[45], we rewrite Eq (5.8) as

Taking the derivative of $\lambda_{2}$ with respect to $t$ in (5.9), we have that

The shadow prices $\mu_{i}(t) = \lambda_{i}(t)e^{\delta t}(i = 1, 2)$ should remain constant when $\lim _{t\rightarrow \infty}\lambda_{i} = 0$ for $i = 1, 2$. Thus the solution of Eq (5.10) is

From Eqs (5.8) and (5.10) we get

That is, $2dy-b\sqrt{x}-\delta = 0$.

According to the interior equilibrium $(x^{*}, y^{*})$, we can get $T = \frac{b^{2}cd+p_{1}(b^{2}(d-b)-2d\delta^{2})}{b^{2}(d-b)-2d\delta^{2}}$. At the same time, we can get the optimal equilibrium solution.

6.   Numerical simulations

We will give some numerical results of system (2.2) to support the analytic results in this section. We consider the values of parameters as follows $b = 0.01, c = 0.02, d = 0.75, p = 2500, m = 604.33$. It is easy to check that the condition (H1): $p-c-\frac{bp}{d} = 2458.3 > 0$ is satisfied. Thus, system (2.2) without any time delay is locally asymptotically stable (see Figure 1).

We have that $\tau_{20} = 2.955$ when system (2.2) is considered under the case $\tau_1 = 0$ and $\tau_2 \neq 0$. From Theorem 2.1, we obtain that the corresponding waveform shown in Figure 2. That is, when $\tau_2 = 2.950 < \tau_{20} = 2.955$, the system is locally asymptotically stable, but Hopf bifurcation occurs once $\tau_2 = 2.960 > \tau_{20} = 2.955$. Similarly, we have that $\tau_{10} = 1.901$ when system (2.2) is considered under the case $\tau_1 \neq 0 \; \tau_2 = 0$ (see Figure 3).

Regarding $\tau_2$ as the parameter with $\tau_1 = 1.0 \in I_{1s} = (0, 1.901)$, we can obtain that $\tau_{2*} = 1.54$. From Theorem 2.3, system (2.2) is locally asymptotically stable at the positive equilibrium for $\tau_2\in(0, \tau_{2*})$ and unstable for $\tau_2 > \tau_{2*}$ (see Figure 4). Lastly, by the algorithm (3.5) derived in Section 3, we have $C_1(0) = -6.947-24.011 {\rm i}$, then $\mu_2 = 1452.01, \beta_2 = -65.83$. Thus, the Hopf bifurcation is supercritical, and the bifurcating periodic solutions are asymptotically stable, which is illustrated in Figure 5.

In addition, we show the effect of the herd behavior of prey species on the dynamics of system (2.2). Frome the Figure 6, we can see that the dynamics of system (2.2) is sensitive to the initial values of prey population. That is, when the number of prey species is small, the herd behavior can cause the fluctuation of population in system (2.2), even extinction. Inversely, much of prey species with herd behavior can keep this system be stable.

At last, Yuan et al. ^[10] pointed out that the Turing pattern is induced by quadratic mortality, which implies that if the predator mortality is given by the linear form, the spatial predator-prey system can not show the Turing structures. We also find similar results when the death rate of predator population in the system (2.2) takes the standard mass-action term (see Figure 7).

7.   Discussion and conclusions

After Ajraldi et al. ^[8] proposed a predator-prey system with square root functional response, many researchers have devoted to studying such systems (see ^{[8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]}). Braza ^[9] investigated the dynamics of such system and showed the prey exhibits strong herd structure and the predator interacts with the prey along the outer corridor of the herd of prey. Yuan et al. ^[10] also considered the quadratic mortality for predator population based on such system and proposed a spatial predator-prey model with Neumann boundary conditions. Xu and Yuan ^[11] introduced a time delay which means the growth rate of predator species to depend on the number of the prey species $\tau$ units of time earlier and studied the local stability and the existence of Hopf bifurcation. Some authors (Meng and Wang ^[19], Yang ^[20,21], Souna et al. ^[22]) obtained dynamical behaviors of such delayed diffusive predator-prey model with herd schooling behavior. Based on the above work, we not only introduce a new time delay denoting the feedback time delay of the prey species, but also add a algebraic equation defining the economic interest of the yield of the harvest effort, which yields system (1.2).

Firstly, we obtained the conditions for the existence of the positive equilibrium. In order to study the dynamics of system (1.2), we introduced the new normal form of differential-algebraic systems due to the work of literatures ^[40,41]. Then, we analyzed the local stability of the positive equilibrium and the existence of Hopf bifurcation by taking time delays $\tau_1$ and $\tau_2$ as bifurcation parameters. Next we gave the explicit formulas determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions by using the normal form theory and the center manifold theorem introduced by Hassard et al. ^[42]. From simulations, we found that the critical value of time delay $\tau_2$ is less than that in the literature ^[11] when $\tau_1 = 0$. From the biological point of view, the predator species has to shorten its time interval to survive when the prey species is predated by natural or man-made factors. Further, if we considered the time delay of prey species belonging to its stable interval, then we found that the critical value of time delay $\tau_2 \, (\tau_{2*} = 1.54)$ is less than that ($\tau_{20} = 2.955$) of system (2.2) without time delay $\tau_1$. From the view of population dynamics, in order to suit the real ecosystem, the predator species needs much less time to turn prey captured into its own newborns. In addition, in order to eliminate the bifurcation, the nonlinear state feedback controller is designed. Of course, the optimal tax policy is also discussed in theory. Unfortunately, we are unable to show the corresponding results in numerical simulations.

Of course, the herd behavior of prey population may not only be good for themselves but also for other population, such as the predator population and harvest effort (see Figue 6). Yuan et al. ^[10] pointed out that the spatial predator-prey system can not show the Turing structures if the predator mortality is given by the linear form. Similarly, we find that the prey species is stable or turns up periodic fluctuation and the predator population dies out as soon as fast when the predator mortality of system (2.2) takes the linear form (see Figure 7).

Aa pointed out in the remark, we only discuss the dynamics of system (2.1) under the case $a = 0$. We may investigate some other dynamics of system (2.1) with $a \neq 0$. For example, we may show some bifurcations and chaos only by numerical simulation. In addition, if the economic interest $m$ is taken as the parameter, then the singularity induced bifurcation may occur in such singular biological economic predator-prey system. In order to make the system more realistic, we can consider more factors into system. The linear harvesting is considered in system (1.2), but the following nonlinear harvesting function can also be considered for fish. Thus, the third equation of system (1.2) is replaced by the following form:

where $E = s\hat{V}$ with $s$ being a proportionality constant and $\hat{V}$ is number of same vessels per unit area, which is used to harvest the population. $X$ is the density of fish at time $t$, $C_{0}$ describes the degree of competition between the vessels, $q$ is the catchability co-efficient and $b$ stands for the proportionality constant for handling time. We leave this work in the future.

Acknowledgments

The author thanks the two reviewers for their comments. This work is supported by the National Natural Science Foundation of China(Grant Nos. 11661050, 62033002, 61833006 and 62071112), Ning Revitalization Talents Program (XLYC1807198), National Natural Science Foundation of Gansu Province (Grant No. 20JR10RA156) and the HongLiu First-class Disciplines Development Program of Lanzhou University of Technology.

Conflict of interest

All authors declare no conflicts of interest in this paper.