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

Abstract: In this paper, we propose a predator-prey system with square root functional response, two delays and prey harvesting, in which an algebraic equation stands for the economic interest of the yield of the harvest effort. Firstly, the existence of the positive equilibrium is discussed. Then, by taking two delays as bifurcation parameters, the local stability of the positive equilibrium and the existence of Hopf bifurcation are obtained. Next, some explicit formulas determining the properties of Hopf bifurcation are analyzed based on the normal form method and center manifold theory. Furthermore, the control of Hopf bifurcation is discussed in theory. What’s more, the optimal tax policy of system is given. Finally, simulations are given to check the theoretical results.


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: Net Economic Revenue (NER) = Total Revenue (TR)-Total Cost (TC). Let E(t) and Y(t) represent the harvest effort and the density of harvested population, respectively. T R = ωE(t)Y(t) and TC = cE(t), where ω 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. τ 1 denotes the feedback time delay of prey species to the growth of species itself, and τ 2 is the time delay which means the growth rate of predator species depends on the number of the prey species τ 2 units of time earlier. The parameter r and N are the growth rate of the prey and its carrying capacity. The parameter α is the search efficiency of Y for X,b is biomass conversion or consumption rate, and t h is Y's average handling time of X. −dY 2 represents the quadratic mortality for predator population.p andc represent unit price of harvested population and the cost of harvest effort.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.

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: 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: (2. 2) 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 τ 1 and τ 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).

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 = b d √ x and e = 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 x * and e * = m px * −c . When ∆ > 0, Eq (2.4) has two positive roots, defined x * 1 and x * 2 ; that is, system (2.2) has two positive equilibria P 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 ≤ 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.

Stability and existence of Hopf bifurcation
For simplicity, let 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 Y g(P * )Q = (0, 0, px * − c) and u(t) = x(t), v(t) = y(t), E(t) = pe * px * −c x(t) + e(t). Substituting the latter into system (2.2), we can obtain (2.5) In order to derive the formula determining the properties of the positive equilibrium of the system (2.2), we consider the local parametric Ψ of the third equation of system (2.2) as in [40,41], which is given as follows: , h(Z(t)) = (0, 0, h 3 (y 1 (t), y 2 (t)) T : R 2 → R is a smooth mapping.
Thus, we obtain the following parametric system of the system (2.5): because g(Ψ(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: 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é's theorem.
It is obvious that λ = 0 is not a root of Eq (2.8). When there is no delay, that is, τ 1 = τ 2 = 0, Eq (2.8) is rewritten as follows: Thus, if u * + 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 √ u * 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 (2.9) have always negative real parts, that is, the positive equilibrium of system (2.2) with τ 1 = τ 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 τ 1 and τ 2 .
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 ±iω + 2 (resp. ±iω − 2 ) are a pair of purely imaginary root of Eq (2.10) with . Differentiating the both sides of Eq (2.10) with respect to τ 2 , then straightforward but tedious calculations lead to the following conclusion.
Proof. Taking the derivative of λ with respect to τ 1 in Eq (2.13), we obtain For simplify, we define ω ± 1 as ω and τ ± 1 j as τ 1 , and we can obtain This completes the proof.
We can summarize the preceding results as the following theorem. 2) is locally asymptotically stable for any τ 1 ≥ 0.
Lemma 2.6. Suppose (H7) and P R Q R − P I Q I > 0 hold, then the following transversality condition Proof. Taking the derivative of λ with respect to τ 2 in Eq (2.8), we obtain For simplify, defining ω i as ω * and τ (i) 2 j as τ 2 , we can obtain sign P R Q R − P I Q I .
If P R Q R − P I Q I > 0 holds, then this ends the proof.

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 τ 2 for fixed value τ 1 * ∈ 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.
(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).

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 g(t,k) = k 1 y 2 (t) + k 2 y 3 2 (t), (4.1) herek = (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 τ 1 = 0 and τ 2 0. The other cases can be investigated similarly, so those are omitted.
In order to give the corresponding result, we give the following assumption Theorem 4.1. If (H8) is true, then the control system (4.2) becomes stable when τ = τ 20 .
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.

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 τ 1 = τ 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 e −δt P(t, x, y, e, T )dt, (5.2) where δ 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 λ i = λ 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 max or T = 0. If ∂H ∂T = 0, then we have λ 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 λ 2 (t) = µ 2 (t)e −δt . (5.9) Taking the derivative of λ 2 with respect to t in (5.9), we have that The shadow prices µ i (t) = λ i (t)e δt (i = 1, 2) should remain constant when lim t→∞ λ i = 0 for i = 1, 2.
Thus the solution of Eq (5.10) is From Eqs (5.8) and (5.10) we get . At the same time, we can get the optimal equilibrium solution.

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 − 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 τ 20 = 2.955 when system (2.2) is considered under the case τ 1 = 0 and τ 2 0. From Theorem 2.1, we obtain that the corresponding waveform shown in Figure 2. That is, when τ 2 = 2.950 < τ 20 = 2.955, the system is locally asymptotically stable, but Hopf bifurcation occurs once τ 2 = 2.960 > τ 20 = 2.955. Similarly, we have that τ 10 = 1.901 when system (2.2) is considered under the case τ 1 0 τ 2 = 0 (see Figure 3).  Regarding τ 2 as the parameter with τ 1 = 1.0 ∈ I 1s = (0, 1.901), we can obtain that τ 2 * = 1.54. From Theorem 2.3, system (2.2) is locally asymptotically stable at the positive equilibrium for τ 2 ∈ (0, τ 2 * ) and unstable for τ 2 > τ 2 * (see Figure 4). Lastly, by the algorithm (3.5) derived in Section 3, we have C 1 (0) = −6.947 − 24.011i, then µ 2 = 1452.01, β 2 = −65.83. Thus, the Hopf bifurcation is supercritical, and the bifurcating periodic solutions are asymptotically stable, which is illustrated in Figure 5.  2) is stable with initial values "0.4, 0.25" when time delay τ 1 belongs to its stable interval and time delay τ 2 is less than its critical value τ 2 * ; (b) Hopf bifurcation occurs when time delay τ 1 belongs to its stable interval and τ 2 is greater than its critical value τ 2 * . 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 ).  Figure 7. The simulations show that species of system (2.2) is stable (resp. fluctuates) when the predator mortality is given by the linear form.

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 τ 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 τ 1 and τ 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 τ 2 is less than that in the literature [11] when τ 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 τ 2 (τ 2 * = 1.54) is less than that (τ 20 = 2.955) of system (2.2) without time delay τ 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 Figue 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 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 = sV with s being a proportionality constant andV 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.