Loading [MathJax]/jax/output/SVG/jax.js
Research article

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

  • Received: 05 February 2021 Accepted: 18 March 2021 Published: 26 March 2021
  • MSC : 92D25, 34D23, 34H05

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

    Citation: Xin-You Meng, Fan-Li Meng. Bifurcation analysis of a special delayed predator-prey model with herd behavior and prey harvesting[J]. AIMS Mathematics, 2021, 6(6): 5695-5719. doi: 10.3934/math.2021336

    Related Papers:

    [1] Heping Jiang . Complex dynamics induced by harvesting rate and delay in a diffusive Leslie-Gower predator-prey model. AIMS Mathematics, 2023, 8(9): 20718-20730. doi: 10.3934/math.20231056
    [2] Kwadwo Antwi-Fordjour, Rana D. Parshad, Hannah E. Thompson, Stephanie B. Westaway . Fear-driven extinction and (de)stabilization in a predator-prey model incorporating prey herd behavior and mutual interference. AIMS Mathematics, 2023, 8(2): 3353-3377. doi: 10.3934/math.2023173
    [3] Yingyan Zhao, Changjin Xu, Yiya Xu, Jinting Lin, Yicheng Pang, Zixin Liu, Jianwei Shen . Mathematical exploration on control of bifurcation for a 3D predator-prey model with delay. AIMS Mathematics, 2024, 9(11): 29883-29915. doi: 10.3934/math.20241445
    [4] San-Xing Wu, Xin-You Meng . Dynamics of a delayed predator-prey system with fear effect, herd behavior and disease in the susceptible prey. AIMS Mathematics, 2021, 6(4): 3654-3685. doi: 10.3934/math.2021218
    [5] Na Min, Hongyang Zhang, Xiaobin Gao, Pengyu Zeng . Impacts of hunting cooperation and prey harvesting in a Leslie-Gower prey-predator system with strong Allee effect. AIMS Mathematics, 2024, 9(12): 34618-34646. doi: 10.3934/math.20241649
    [6] Rongjie Yu, Hengguo Yu, Min Zhao . Steady states and spatiotemporal dynamics of a diffusive predator-prey system with predator harvesting. AIMS Mathematics, 2024, 9(9): 24058-24088. doi: 10.3934/math.20241170
    [7] Jie Liu, Qinglong Wang, Xuyang Cao, Ting Yu . Bifurcation and optimal harvesting analysis of a discrete-time predator–prey model with fear and prey refuge effects. AIMS Mathematics, 2024, 9(10): 26283-26306. doi: 10.3934/math.20241281
    [8] Jawdat Alebraheem . Asymptotic stability of deterministic and stochastic prey-predator models with prey herd immigration. AIMS Mathematics, 2025, 10(3): 4620-4640. doi: 10.3934/math.2025214
    [9] Chenxuan Nie, Dan Jin, Ruizhi Yang . Hopf bifurcation analysis in a delayed diffusive predator-prey system with nonlocal competition and generalist predator. AIMS Mathematics, 2022, 7(7): 13344-13360. doi: 10.3934/math.2022737
    [10] Kottakkaran Sooppy Nisar, G Ranjith Kumar, K Ramesh . The study on the complex nature of a predator-prey model with fractional-order derivatives incorporating refuge and nonlinear prey harvesting. AIMS Mathematics, 2024, 9(5): 13492-13507. doi: 10.3934/math.2024657
  • 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.



    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:

    {dXdt=rX(1XN)αXY1+thαX,dYdt=sY2+cαXY1+thαX, (1.1)

    where sY2 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. TR=ω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:

    E(t)(ωY(t)c)=v.

    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:

    {dXds=rX(1X(sτ1)N)αXY1+thαXEX,dYds=˜dY2+˜bαX(sτ2)Y1+thαX(sτ2),0=E(˜pX˜c)˜m, (1.2)

    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 th is Y's average handling time of X. ˜dY2 represents the quadratic mortality for predator population. ˜p and ˜c 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.

    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:

    x=XN,y=αYrN,e=Er,t=rs.

    Thus system (1.2) can be rewritten as following:

    {dxdt=x(1x(tτ1))xy1+axex,dydt=dy2+bx(tτ2)y1+ax(tτ2),0=e(pxc)m, (2.1)

    where a=thαN,b=αN˜br,c=˜c,d=N˜dr,m=˜mr,p=˜pN.

    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:

    {dxdt=x(1x(tτ1))xyex,dydt=dy2+bx(tτ2)y,0=e(pxc)m. (2.2)

    The initial conditions of system (2.2) have the form

    x(θ)=ψ(θ)0,y(θ)=η(θ)0,e(0)>0,θ[τ,0),τ=max{τ1,τ2}.

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

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

    x(1x)xyex=0, (2.3a)
    dy2+bxy=0, (2.3b)
    e(pxc)m=0. (2.3c)

    From Eq (2.3b) and (2.3c), we can obtain y=bdx and e=mpxc, respectively. Substituting the above values into Eq (2.3a), we can obtain that x satisfies the following equation

    x2+(bd1cp)x+cp+mpbcdp=0. (2.4)

    It is obvious that system Eq (2.4) has positive roots if the condition Δ=(bd1cp)24(cp+mpbcdp)0. If Δ=0, then Eq (2.4) has a unique positive root, defined x=cd+pdbp2dp. That is, system (2.2) has only one positive equilibrium P(x,y,e) if the assumption

    (H1): p>c+bpd

    holds, here y=bdx and e=mpxc. When Δ>0, Eq (2.4) has two positive roots, defined x1 and x2; that is, system (2.2) has two positive equilibria Pi(xi,yi,ei)i=1,2, here yi=bdxi,ei=mpxic(i=1,2) and x1=cd+pdbp+Δ2dp>x,x2=cd+pdbpΔ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 a0 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 Pi(xi,yi,ei)(i=1,2) by the same way.

    For simplicity, let

    f(X)=(f1(x,y,e)f2(x,y,e))=(x(1x(tτ1))xyexdy2+bx(tτ2)y),g(X)=e(pxc)m,

    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

    N(t)=(u(t)v(t)E(t)),Q=(100010pepxc01).

    Then, we get DYg(P)Q=(0,0,pxc) and u(t)=x(t),v(t)=y(t),E(t)=pepxcx(t)+e(t). Substituting the latter into system (2.2), we can obtain

    {du(t)dt=u(t)(1u(tτ1))u(t)v(t)E(t)u(t)+peu2(t)puc,dv(t)dt=dv2(t)+bu(tτ2)v(t),0=(E(t)peu(t)puc)(pu(t)c)m. (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:

    N(t)=Ψ(Z(t))=N0+u0Z(t)+v0h(Z(t)),g(Ψ(Z(t)))=0,

    where u0=(100100),v0=(001),Z(t)=(y1(t)y2(t)),N0=(uvE), h(Z(t))=(0,0,h3(y1(t),y2(t))T:R2R is a smooth mapping.

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

    {dy1(t)dt=(u+y1(t))[1(y1(tτ1)+u)]u+y1(t)(v+y2(t))(u+y1(t))(E+h3(y1(t),y2(t)))+pepuc(u+y1(t))2,dy2(t)dt=d(v+y2(t))2+bu+y1(tτ2)(v+y2(t)), (2.6)

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

    {dy1(t)dt=(1b2duE+2peupuc)y1(t)uy2(t)uy1(tτ1),dy2(t)dt=buy2(t)+b22dy1(tτ2). (2.7)

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

    λ2+(buA)λAbu+u(λ+bu)eλτ1+b2u2deλτ2=0, (2.8)

    where A=1b2du+mc(puc)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ˊe's theorem.

    Lemma 2.1. Consider the exponential polynomial

    P(λ,eλτ1,,eλτm)=λn+p01λn1++p0n1λ+p0n+[p11λn1++p1n1λ+p1n]eλτ1++[pm1λn1++pmn1λ+pmn]eλτm,

    where τi0(i=1,2,,m), pij(i=1,2,,m;j=1,2,,n) are constants. As (τ1,τ2,,τm) vary, the sum of the order of the zeros of P(λ,eλτ1,,eλτm) on the open right half plane can change only if a zero appears on or crosses the imaginary axis.

    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:

    λ2+(u+buA)λ+bu(u+b2dA)=0. (2.9)

    Thus, if u+b2dA<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+mc(puc)2<min{2u+bd,2u+b2d+bu}

    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.

    Case one: τ1=0 and τ20.

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

    λ2+A11λ+A10+B10eλτ2=0, (2.10)

    where A11=bu+uA,A10=bu(uA),B10=b2u2d.

    Now, we can obtain the same result as Lemma 2.2 when τ2=0. For τ2>0, let ωi(ω>0) be a root of Eq (2.10), then ω satisfies the following equation:

    ω2A11ωi+A10+B10(cosωτ2isinωτ2)=0.

    Separating the real and imaginary parts, we obtain

    {B10cosωτ2=ω2A10,B10sinωτ2=A11ω, (2.11)

    which yields that

    ω4+(A2112A10)ω2+A210B210=0. (2.12)

    Further, if denote v=ω2, we can obtain Φ1(ν)=ν2+(A2112A10)ν+A210B210. It is obvious that Eq (2.12) dose not have positive root when A210B210>0 and A2112A100 or Δ1(A2112A10)24(A2112A10)<0.

    If the following assumption

    (H3): A210B2100

    holds, then Eq (2.12) has a positive root ω+2, here ω+2=2A10A211+Δ12.

    If the following assumption

    (H4): A210B210>0, A2112A10<0 and Δ1>0

    holds, then Eq (2.12) has two positive roots ω+2 and ω2, ω±2=2A10A211±Δ12.

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

    τ+2j=1ω+2{arccos(ω+2)2A10B10+2jπ},j=0,1,2,,

    then ±iω+2 are a pair of purely imaginary root of Eq (2.10) with τ2=τ+2j.

    Define τ20=min{τ+2j,j=0,1,2,}. Further, if denote λ(τ)=α(τ)+iω(τ) be the root of (2.10), then it satisfies α(τ20)=0,ω(τ20)=ω20. Next, we will investigate whether the transversality condition is satisfied. Differentiating the two sides of Eq (2.10) with respect to τ2, we get

    (dλdτ2)1=2λ+A11λB10eλτ2τ2λ.

    For simplify, we define ω+2 as ω and τ+2j as τ2, and obtain

    sign{d(Reλ)dτ2|τ2=τ+2j}=sign{Re[dλ(τ2)dτ2]1τ2=τ+2j}=sign{Re(2λ+A11λB10eλτ2τ2λ)|λ=ω+2i}=sign{1B210(2ω2++A2112A10)}=1B210sign{Φ1(ν)}.

    Thus, we have the following conclusion.

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

    dReλ(τ2)dτ2|τ2=τ+2j>0.

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

    τ±2j=1ω±2{arccos(ω±2)2A10B10+2jπ},j=0,1,2,,

    then ±iω+2 (resp. ±iω2) are a pair of purely imaginary root of Eq (2.10) with τ2=τ+2j (resp. τ2=τ2j). Further, we can check that there exists a positive integral k>0 such that τ+2k+1>τ2k and 0<τ+20<τ20<τ+21<τ21τ2k1<τ+2k<τ2k<τ+2k+1.

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

    Lemma 2.4. Assume (H1) and Φ1(ν)0 hold, then the transversality condition

    Sign{dReλ(τ2)dτ2|τ2=τ+2j}=Sign{Φ1(ν)},Sign{dReλ(τ2)dτ2|τ2=τ2j}=Sign{Φ1(ν)}

    are satisfied, and (dλdτ2)1 and Φ1(ν) 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 (2.2) is locally asymptotically stable for τ2[0,τ20) and unstable for τ2>τ20. That is, system (2.2) undergoes a Hopf bifurcation at the positive equilibrium P(x,y,e) when τ2=τ+2j,j=0,1,2,.

    (2) If (H4) and Φ1(ν)0, then there exists NN, when τ2(0,τ+20)(τ20,τ+21)(τ2N,τ+2N+1), all the roots of Eq (2.10) have negative real part, thus system (2.2) is locally asymptotically stable; when τ2(τ+20,τ20)(τ+21,τ21)(τ+2N,τ2N)(τ+2N+1,+), Eq (2.10) has at least one root with positive real part, thus system (2.2) is unstable. Hopf bifurcation occurs near the positive equilibrium P when τ2=τ±2j,j=0,1,2,.

    Case two: τ10 and τ2=0.

    Eq (2.8) takes the form

    λ2+A21λ++A20+(B21λ+B20)eλτ1=0, (2.13)

    where A21=buA,A20=bu2dAbu,B21=u,B20=buu.

    Let ωi(ω>0) be the root of Eq (2.13). We obtain

    {B21ωsinωτ1+B20cosωτ1=ω2A20,B21ωcosωτ1B20sinωτ1=A21ω, (2.14)

    which gives

    ω4+(A221B2212A20)ω2+A220B220=0. (2.15)

    Further, if denote v=ω2, we can obtain Φ2(ν)=ν2+ϱν+A220B220, here ϱ=A221B2212A20. It is obvious that Eq (2.15) dose not have positive root when A221B2212A20>0 and A220B220>0 or Δ2(A221B2212A20)24(A220B220)<0. We give two assumptions as following:

    (H5): A220B2200;

    (H6): A221B2212A20<0, A220B220>0 and Δ2>0.

    If the assumption (H5) holds, then Eq (2.15) has only one positive root, defined ω+1. If the assumption (H6) holds, then Eq (2.15) has two positive roots ω±1=12[ϱ±Δ2], and ω+1>ω1>0.

    From Eq (2.14), if denote

    τ±1j=1ω±1{arccos(B20A21B21)(ω±1)2A20B20B220+B221(ω±1)2+2jπ},j=0,1,2,,

    then ±iω+1 (resp. ±iω1) are a pair of purely imaginary root of Eq (2.13) with τ1=τ+1j (resp. τ1=τ1j).

    When (H6) holds, we define τ10=min{τ±1j,j=0,1,2,}. Let λ(τ)=ξ(τ)+iω(τ) be the root of Eq (2.13) near τ=τ10 satisfying ξ(τ10)=0, ω(τ10)=ω10. Further, if (H2) holds for τ1=0, k switches from stability to instability to stability occur when the parameters are such that 0<τ+10<τ10<τ+11<<τ+1k<τ1k<τ+1k+1.

    The following transversality condition is true.

    Lemma 2.5. Suppose (H6) and Φ2(ω2)0 hold, then the following transversality condition

    dReλ(τ1)dτ1|τ1=τ+1j>0anddReλ(τ1)dτ1|τ1=τ1j<0

    are satisfied, (dλdτ1)1 and Φ2(ν) have the same sign.

    Proof. Taking the derivative of λ with respect to τ1 in Eq (2.13), we obtain

    (dλdτ1)1=(2λ+A21)eλτ1+B21λ(B21λ+B20)τ1λ.

    For simplify, we define ω±1 as ω and τ±1j as τ1, and we can obtain

    sign{d(Reλ)dτ1|τ1=τ±1j}=sign{Re(dλdτ1)1|τ1=τ±1j}=sign{Re[(2λ+A21)λ(λ2+A21λ+A20)]|λ=±iω±1+Re[B21λ(B21λ+B20)]λ=±iω±1}=sign{2ω2+A2212A20(A20ω2)2+(A21ω)2+B221B220+(B21ω)2}=1B220+B221ω2sign{Φ2(ω2)}.

    This completes the proof.

    We can summarize the preceding results as the following theorem.

    Theorem 2.2. For the system (2.2), if (H1) and (H2) hold, the following results are true.

    (1) If Eq (2.15) does not have positive root, then system (2.2) is locally asymptotically stable for any τ10.

    (2) If (H5) and Φ2(ω)0 hold, then system (2.2) is locally asymptotically stable for τ1[0,τ10) and unstable for τ1>τ10. System (2.2) undergoes a Hopf bifurcation at the equilibrium P when τ1=τ+1j(j=0,1,2,).

    (3) If (H6) and Φ2(ω)0, then there exists NN, when τ1(0,τ+10)(τ10,τ+11)(τ1N,τ+1N+1), all the roots of Eq (2.15) have negative real parts, thus system (2.2) is locally asymptotically stable; when τ1(τ+10,τ10)(τ+11,τ11)(τ+1N,τ1N)(τ+1N+1,+), Eq (2.15) has at least one root with positive real part, thus system (2.2) is unstable. Hopf bifurcation occurs near the positive equilibrium P when τ1=τ±1j(j=0,1,2,).

    Remark 3 In case of τ10 and τ2=0, whether (H6) or (H7) holds, we define the stable interval of time delay τ1 as I1s=[0,τ10) (resp. I1s=[0,τ10)(τ10,τ+11)(τ1N,τ+1N+1).

    Case three: τ1>0 and τ2>0.

    This case states that τ2 is regarded as a parameter and τ1 belongs to its stable interval. Let ωi(ω>0) be a root of Eq (2.8) for any τ1I1s, then we can obtain

    {b2u2dcosωτ2=ω2+Abubuucosωτ1ωusinωτ1,b2u2dsinωτ2=(buA)ωbuusinωτ1+ωucosωτ1, (2.16)

    which yields

    ω4+(A2+b2u2)ω2+A2b2u2Au(ω2+b2u)cosωτ12u(ω3+b2uω)sinωτ1=0. (2.17)

    Without loss of generality, suppose that

    (H7)Eq(2.17)has at least finite positive roots.

    If (H7) holds, the roots of Eq (2.17) defined by ω1,ω2,,ωk. For every fixed ωi(i=1,2,,k), there exists a sequence

    τ(i)2j=1ωi{arccos2d(ω2i+Abubuucosωiτ1ωiusinωiτ1)b2u+2jπ},i=1,2,,k,j=0,1,2,,

    then ±iωi are a pair of purely imaginary root of (2.8) with τ2=τ(i)2j.

    Define

    τ2=τ(i)20=mini{1,2,,k}{τ(i)2j},j=0,1,2,.

    Let λ(τ)=ξ(τ)+iω(τ) be the root of Eq (2.8) near τ=τ2 satisfying ξ(τ2)=0, ω(τ2)=ω. Further, we obtain the following lemma.

    Lemma 2.6. Suppose (H7) and PRQRPIQI>0 hold, then the following transversality condition

    dReλ(τ2)dτ2|τ2=τ(i)2j>0

    is satisfied, where

    PR=buA+u[(1bτ1u)cosωiτ1ωiτ1sinωiτ1],PI=2ωiu[(1bτ1u)sinωiτ1+ωiτ1cosωiτ1],QR=ω2i+Abubuucosωiτ1ωiusinωiτ1,QI=ωi(buA)buusinωiτ1+ωiucosωiτ1.

    Proof. Taking the derivative of λ with respect to τ2 in Eq (2.8), we obtain

    (dλdτ2)1=2λ+buA+u(1τ1λbτ1u)eλτ1b2u2dλeλτ2τ2λ.

    For simplify, defining ωi as ω and τ(i)2j as τ2, we can obtain

    sign{d(Reλ)dτ2|τ2=τ(i)2j}=sign{Re(dλdτ2)1|τ2=τ(i)2j}=sign{Re[2λ+buA+u(1τ1λbτ1u)eλτ1b2u2dλeλτ2τ2λ]λ=ωi}=4d2ωub4(Q2R+Q2I)sign{PRQRPIQI}.

    If PRQRPIQI>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 (2.2), suppose (H7) and PRQRPIQI>0 hold. Then system (2.2) is locally asymptotically stable when τ2[0,τ2) and τ1I1s, and system (2.2) undergoes a Hopf bifurcation at the positive equilibrium P when τ2=τ2. That is, system (2.2) has a branch of periodic solutions bifurcating from Hopf bifurcation near τ2=τ2.

    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 τ1I1s, 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 τ2=τ2,τ1I1s. Without loss of generality, we assume that τ1<τ2. Let x1(t)=y1(t)y1,x2(t)=y2(t)y2,t=t/τ2,τ2=τ2+μ and ˉxi(t)=xi(τ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],R2) as

    ˙x(t)=Lμ(xt)+F(μ,xt), (3.1)

    where x(t)=(x1(t),x2(t))TR2, and Lμ:CR2,F:R×CR2 are given, respectively, by

    Lμ(ϕ)=(τ2+μ)[B(τ2)ϕ(0)+C(τ2)ϕ(τ1τ2)+D(τ2)ϕ(1)],

    and

    F(μ,xt)=(τ2+μ)(F1,F2)T,

    where

    B(τ2)=(Au0bu),C(τ2)=(u000),D(τ2)=(00b2u2d0),

    and

    {F1=a11ϕ21(0)+a12ϕ1(0)ϕ2(0)+a13ϕ1(0)ϕ1(τ1τ2)+a111ϕ31(0)+a112ϕ21(0)ϕ2(0),F2=b11ϕ21(1)+b12ϕ1(1)ϕ2(0)+b22ϕ22(0)+b111ϕ31(1)+b112ϕ21(1)ϕ2(0),

    here ϕ(θ)=(ϕ1(θ),ϕ2(θ))TC([1,0],R2),

    a11=b4du+mp(puc)2mp(pu+c)(puc)2,a12=12u,a13=1,a111=3b(8du)22mp2(puc)3+4mp2(pu+2c)(puc)4,a112=14uu,b11=b24du,b12=b2u,b22=2d,b111=3b28d(u)2,b112=b4uu.

    System (3.1) turns to the linear problem

    ˙x(t)=Lμxt

    by the Riesz representation theorem. There exists 2×2 matrix-valued function

    η(,μ):[1,0]R2×2,

    such that

    Lμ(ϕ)=01dη(θ,μ)ϕ(θ),ϕC. (3.2)

    In fact, we can choose

    η(θ,μ)={(τ2+μ)(B+C+D),θ=0,(τ2+μ)(C+D),θ[τ1τ2,0),(τ2+μ)D,θ(1,τ1τ2),0,θ=1.

    Then Eq (3.2) is satisfied.

    Next, for ϕC([1,0],R2), we define the operator A(μ) as

    A(μ)ϕ(θ)={dϕdθθ[1,0),01dη(θ,μ)ϕ(θ), θ=0,

    and

    R(μ)ϕ(θ)={0θ[1,0),F(μ,ϕ), θ=0.

    Since dutdθ=dutdt, then system (3.1) is equivalent to the following operator equation

    ˙u(t)=A(μ)ut+R(μ)ut, (3.3)

    where ut=u(t+θ) for θ[1,0].

    For ψC([1,0],(R2)), we further define the adjoint A of A as

    A(μ)ψ(s)={dψ(s)dss(0,1],01ψ(s)dη(s,μ),s=0,

    and a bilinear form

    ψ(s),ϕ(θ)=ˉψT(0)ϕ(0)0θ=1θξ=0ˉψT(ξs)dη(θ)ϕ(ξ)dξ,

    where η(θ)=η(θ,0). Then A(0) and A(0) are adjoint operators. From the above analysis, we obtain that ±iωτ2 are the eigenvalues of A(0) and therefore they are also the eigenvalues of A(0). Let q(θ) be the eigenvector of A(0) corresponding to iωτ2 and q(s) be the eigenvector of A(0) corresponding to iωτ2, we have

    A(0)q(θ)=iωτ2q(θ),A(0)q(s)=iωτ2q(s).

    By some simple computation, we can obtain

    q(θ)=Veiωτ2θ,q(s)=DVeiωτ2s,

    where

    V=(1,ρ)T,ρ=1u(Aueiωτ1ieiωτ2),V=(1,ρ)T,ρ=2du2dωi+b2eiωτ2,ˉD=2d2d+2dρ¯ρ+b2τ2eiωτ2ρ¯ρ.

    Then <q,q>=1, and <q,ˉ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 C0 at μ=0 and obtain the following expressions:

    x1t(0)=z+ˉz+W(1)(0),x2t(0)=ρz+¯ρˉz+W(2)(0),x1t(1)=eiωτ2z+eiωτ2ˉz+W(1)(1),x2t(1)=ρeiωτ2z+¯ρeiωτ2ˉz+W(2)(1),x1t(τ1τ2)=eiωτ1z+eiωτ1ˉz+W(1)(τ1τ2),x2t(τ1τ2)=ρeiωτ1z+¯ρeiωτ1ˉz+W(2)(τ1τ2).

    From g(z,ˉz)=ˉq(0)F0(z,ˉz)=g20z22+g11zˉz+g02ˉz22+g21z2ˉz2, we obtain the coefficients determining the important quantities of the periodic solution:

    {g20=2ˉD[a11+a12ρ+a13eiωτ1+b11¯ρeiωτ2+b12¯ρρeiωτ2+b22¯ρρ2],g11=ˉD[2a11+a12(ρ+ˉρ)+a13(eiωτ1+eiωτ1)+2b11¯ρ+b12¯ρ(ρeiωτ2+ˉρeiωτ2)+2b22¯ρρˉρ],g02=2ˉD[a11+a12ˉρ+a13eiωτ1+b11¯ρeiωτ2+b12¯ρˉρeiωτ2++b22¯ρˉρ2],g21=2ˉD{a11(2W(1)11(0)+W(1)20(0))+a12[12W(2)11(0)W(2)20(0)+12ˉρW(1)20(0)+ρW(1)11(0)]+a13[W(1)11(τ1τ2)+12W(1)20(τ1τ2)+12W(1)20(0)eiωτ2+W(1)11(0)eiωτ1]+3a111+a112(2+ˉρ)+b11¯ρ[2W(1)11(1)eiωτ2+W(1)20(1)eiωτ2]+b12¯ρ[W(2)11(0)eiωτ2+12W(2)20(0)eiωτ2+ρW(1)11(1)+12ˉρW(1)20(1)]+b22¯ρ[2ρW(2)11(0)+ˉρW(2)20(0)]+3b11¯ρeiωτ2+b11¯ρ(2ρ+e2iωτ2)}. (3.4)

    Since W20(θ) and W11(θ) appears in g21, we need to compute them as follows.

    W20(θ)=ig20q(0)ωτ2eiωτ2θ+iˉg02ˉq(0)3ωτ2eiωτ2θ+E1e2iωτ2θ,W11(θ)=ig11q(0)ωτ2eiωτ2θ+iˉg11ˉq(0)ωτ2eiωτ2θ+E2,

    where E1=(E(1)1,E(2)1)T and E2=(E(1)2,E(2)2)T are both two dimensional vectors and can be determined by

    E1=2(2iωA+ue2iωτ2u02iω+bub22de2iωτ2)1(K11K21),

    and

    E2=(Auub2u2dbu)1(K12K22),

    here

    K11=a11+a12ρ+a13eiωτ1,K21=2a11+a12(ρ+ˉρ)+a13(eiωτ1+eiωτ1),K12=b11eiωτ2+b12ρeiωτ2+b22ρ2,K22=2b11+b12(ρeiωτ2+ˉρeiωτ2)+2b22ρˉρ.

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

    {C1(0)=i2ωτ2(g20g112|g11|2|g02|23)+g212,μ2=Re{C1(0)}Re{λ(τ2)},β2=2Re{C1(0)},T2=Im{C1(0)}+μ2Im{λ(τ2)}ωτ2, (3.5)

    which determine the properties of bifurcation periodic solutions in the center manifold at the critical value τ2. By the results of Hassard et al. [42], we have the following results.

    Theorem 3.1. In (3.5), the following results hold.

    (i) The sign of μ2 determines the direction of the Hopf bifurcation: If μ2>0 (resp. μ2<0), then the Hopf bifurcation is supercritical (resp. subcritical) and the bifurcation periodic solutions exist for τ2>τ2 (τ2<τ2).

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

    (iii) The sign of T2 determines the period of the bifurcating periodic solutions: The period increases (resp. decreases) if T2>0 (resp. T2<0).

    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)=k1y2(t)+k2y32(t), (4.1)

    here ˉk=(k1,k2) and k1,k2 are feedback gains.

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

    {.y1(t)=(1b2duE+2peupuc)y1(t)uy2(t)uy1(tτ1),.y2(t)=buy2(t)+b22dy1(tτ2)+k1y2(t)+k2y32(t). (4.2)

    In what follows, we will only discuss one case τ1=0 and τ20. The other cases can be investigated similarly, so those are omitted.

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

    (H8)k1>A10Auand(A10k1(uA))(A11k1)2<B10ek1A112τ20.

    Theorem 4.1. If (H8) is true, then the control system (4.2) becomes stable when τ=τ20.

    Proof. When τ2=τ20, the associated characteristic transcendental equation of system (4.2) around (0,0) is

    λ2+(A11k1)λ+(A10k1(uA))+B10eλτ2=0, (4.3)

    where A11,A10 and B10 are the same to that in Eq (2.11)

    Let g1(λ)=λ2+(A11k1)λ+(A10k1(uA)) and g2(λ)=B10eλτ2. In order to investigate the solution of system (4.3), we only need to discuss the solution of g1(λ)=g2(λ). By a series of analysis, if k1>A10Au and (A10k1(uA))(A11k1)2<B10ek1A112τ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.

    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

    {.x(t)=x(1x)xymx(p1T)xc,.y(t)=dy2+bxy. (5.1)

    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

    P(t,x,y,e,T)=e((p1T)xc)+exT=e(p1xc).

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

    J(T)=0eδtP(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

    maxJ(T)=0eδP(t,x,y,e,T)dt. (5.3)

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

    H(t,x,y,e,T)=e(p1xc)eδt+λ1(x(1x)xymx(p1T)xc)+λ2(dy2+bxy).

    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:

    ˉT={Tmax,HT>0,0,HT<0. (5.4)

    It is assumed that the optimal solution does not occur at T=Tmax or T=0. If HT=0, then we have λ1(t)=0. Based on Pontryagin's maximum principle [44], the adjoint variables must satisfy the following adjoint equations

    dλ1dt=Hxanddλ2dt=Hy. (5.5)

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

    {(1x)xyxx(p1T)xc=0,dy+bx=0. (5.6)

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

    dλ1dt=ep1eδt+λ2by2x. (5.7)

    From Eq (5.7) we get

    λ2(t)=2ep1eδtbyx. (5.8)

    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

    dμ2dtδμ2=2ep1(2dybx)byx. (5.10)

    The shadow prices μi(t)=λi(t)eδt(i=1,2) should remain constant when limtλi=0 for i=1,2. Thus the solution of Eq (5.10) is

    μ2=2ep1(2dybx)δbyx.

    From Eqs (5.8) and (5.10) we get

    2ep1byx2ep1(2dybx)δbyx=0. (5.11)

    That is, 2dybxδ=0.

    According to the interior equilibrium (x,y), we can get T=b2cd+p1(b2(db)2dδ2)b2(db)2dδ2. At the same time, we can get the optimal equilibrium solution.

    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): pcbpd=2458.3>0 is satisfied. Thus, system (2.2) without any time delay is locally asymptotically stable (see Figure 1).

    Figure 1.  Dynamical responses of differential-algebraic system model (2.2) with τ1=τ2=0.

    We have that τ20=2.955 when system (2.2) is considered under the case τ1=0 and τ20. 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 τ10τ2=0 (see Figure 3).

    Figure 2.  The trajectories of system (2.2) with initial values "0.4, 0.25" and different values of time delay τ2. (a) System (2.2) is locally asymptotically stable when τ2=2.950<τ20=2.955; (b) Hopf bifurcation occurs around the positive equilibrium when τ2=2.960>τ20=2.955.
    Figure 3.  The trajectories of system (2.2) with initial values "0.4, 0.25" and different values of time delay τ1. (a) System (2.2) is locally asymptotically stable when τ1=1.900<τ10=1.901; (b) Hopf bifurcation occurs around the positive equilibrium when τ1=1.902>τ10=1.900.

    Regarding τ2 as the parameter with τ1=1.0I1s=(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 C1(0)=6.94724.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.

    Figure 4.  The simulations show the phase portrait of system. (a) System (2.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.
    Figure 5.  This simulation shows that the periodic solutions bifurcating from Hopf bifurcation are stable according to the Theorem 4.1.

    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.

    Figure 6.  The simulations show that herd behavior of system (2.2) depends on the initial value of the prey species. (a) Prey population dies out and harvest effort fluctuates wildly when the initial values are "0.2, 0.5, 1.2"; (b) Prey and predator of system (2.2) with sub-figure(a); (c) System (2.2) is stable when the initial values are "0.8, 0.5, 1.2".

    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.

    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 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 a0. 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:

    0=qEsC10+E+bsX(˜pX˜c)v

    where E=sˆV with s being a proportionality constant and ˆ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, C0 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.

    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.

    All authors declare no conflicts of interest in this paper.



    [1] C. S. Holling, The functional response of invertebrate predators to prey density, Mem. Entomol. Soc. Can., 98 (1966), 5–86.
    [2] R. M. May, Stability and Complexity in Model Ecosystemsm, New Jersey: Princeton University Press, 1973.
    [3] H. I. Freedman, Deterministic Mathematical Models in Population Ecology, New York: Dekker, 1980.
    [4] P. H. Leslie, J. C. Gower, The properties of a stochastic model for the predator-prey type of interaction between two species, Biometrika, 47 (1960), 219–234. doi: 10.1093/biomet/47.3-4.219
    [5] P. H. Crowley, E. K. Martin, Functional responses and interference within and between year classes of a dragonfly population, J. N. Am. Benthol. Soc., 8 (1989), 211–221. doi: 10.2307/1467324
    [6] J. R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, J. Anim. Ecol., 44 (1975), 331–340. doi: 10.2307/3866
    [7] D. L. DeAngelis, R. A. Goldstein, R. V. O'Neill, A model for tropic interaction, Ecology, 56 (1975), 881–892. doi: 10.2307/1936298
    [8] V. Ajraldi, M. Pittavino, E. Venturino, Modeling herd behavior in population systems, Nonlinear Anal.: Real World Appl., 12 (2011), 2319–2338. doi: 10.1016/j.nonrwa.2011.02.002
    [9] P. A. Braza, Predator-prey dynamics with square root functional responses, Nonlinear Anal.: Real World Appl., 13 (2012), 1837–1843. doi: 10.1016/j.nonrwa.2011.12.014
    [10] S. L. Yuan, C. Q. Xu, T. H. Zhang, Spatial dynamics in a predator-prey model with herd behavior, Chaos, 23 (2013), 033102. doi: 10.1063/1.4812724
    [11] C. Q. Xu, S. L. Yuan, Stability and Hopf bifurcation in a delayed predator-prey system with herd behavior, Abstr. Appl. Anal., 2014 (2014), 568943.
    [12] E. Venturino, S. Petrovskii, Spatiotemporal behavior of a prey-predator system with a group defense for prey, Ecol. Complexity, 14 (2013), 37–47. doi: 10.1016/j.ecocom.2013.01.004
    [13] C. Q. Xu, S. L. Yuan, T. H. Zhang, Global dynamics of a predator-prey model with defense mechanism for prey, Appl. Math. Lett., 62 (2016), 42–48. doi: 10.1016/j.aml.2016.06.013
    [14] B. W. Kooi, E. Venturino, Ecoepidemic predator-prey model with feeding satiation, prey herd behavior and abandoned infected prey, Math. Biosci., 274 (2016), 58–72. doi: 10.1016/j.mbs.2016.02.003
    [15] M. Banerjee, B. W. Kooi, E. Venturino, An ecoepidemic model with prey herd behavior and predator feeding saturation response on both healthy and diseased prey, Math. Model. Nat. Phenom., 12 (2017), 133–161. doi: 10.1051/mmnp/201712208
    [16] X. S. Tang, Y. L. Song, Bifurcation analysis and Turing instability in a diffusive predator-prey model with herd behavior and hyperbolic mortality, Chaos Solitons Fractals, 81 (2015), 303–314. doi: 10.1016/j.chaos.2015.10.001
    [17] X. S. Tang, Y. L. Song, T. H. Zhang, Turing-Hopf bifurcation analysis of a predator-prey model with herd behavior and cross-diffusion, Nonlinear Dyn., 86 (2016), 73–89. doi: 10.1007/s11071-016-2873-3
    [18] Q. Liu, D. Q. Jiang, T. Hayat, A. Alsaedi, Stationary distribution and extinction of a stochastic predator-prey model with herd behavior, J. Franklin Inst., 355 (2018), 8177–8193. doi: 10.1016/j.jfranklin.2018.09.013
    [19] X. Y. Meng, J. G. Wang, Dynamical analysis of a delayed diffusive predator-prey model with schooling behavior and Allee effect, J. Biol. Dyn., 14 (2020), 826–848. doi: 10.1080/17513758.2020.1850892
    [20] W. B. Yang, Analysis on existence of bifurcation solutions for a predator-prey model with herd behavior, Appl. Math. Model., 53 (2018), 433–446. doi: 10.1016/j.apm.2017.09.020
    [21] W. B. Yang, Effect of cross-diffusion on the stationary problem of a predator-prey system with a protection zone, Comput. Math. Appl., 76 (2018), 2262–2271. doi: 10.1016/j.camwa.2018.08.025
    [22] F. Souna, A. Lakmeche, S. Djilali, Spatiotemporal patterns in a diffusive predator-prey model with protection zone and predator harvesting, Chaos Solitons Fractals, 140 (2020), 110180. doi: 10.1016/j.chaos.2020.110180
    [23] S. Bentout, S. Djilali, S. Kumar, Mathematical analysis of the influence of prey escaping from prey herd on three species fractional predator-prey interaction model, Phys. A, 572 (2021), 125840. doi: 10.1016/j.physa.2021.125840
    [24] S. G. Ruan, J. J. Wei, On the zeros of transcendental functionas with applications to stability of delay differential equations with two delays, Dyn. Contin. Discrete Impulsive Sys., 10 (2003), 863–874.
    [25] J. K. Hale, Theory of Functional Differential Equations, New York: Springer-Verlag, 1977.
    [26] Y. Kuang, Delay Differerntial Equations with Application in Population Dynamics, New York: Academic Press, 1993.
    [27] E. Beretta, Y. Kuang, Global analyses in some delayed ratio-dependent predator-prey systems, Nonlinear Anal.: Theory Methods Appl., 32 (1998), 381–408. doi: 10.1016/S0362-546X(97)00491-4
    [28] X. Y. Meng, J. Li, Dynamical behavior of a delayed prey-predator-scavenger system with fear effect and linear harvesting, Int. J. Biomath., 2021. Available from: https://doi.org/10.1142/S1793524521500248.
    [29] X. Y. Meng, N. N. Qin, H. F. Huo, Dynamics of a food chain model with two infected predators, Int. J. Bifurcation Chaos, 31 (2021), 2150019. doi: 10.1142/S021812742150019X
    [30] S. Djilali, B. Ghanbari, The influence of an infectious disease on a prey-predator model equipped with a fractional-order derivative, Adv. Differ. Equations, 2021 (2021), 20. Available from: https://doi.org/10.1186/s13662-020-03177-9.
    [31] H. S. Gordon, The economic theory of a common property resource: The fishery, J. Political Econ., 62 (1954), 124–142. doi: 10.1086/257497
    [32] D. G. Luenberger, A. Arbel, Singular dynamic Leontief systems, Econometrics, 45 (1997), 991–995.
    [33] Y. Zhang, Q. L. Zhang, Chaotic control based on descriptor bioeconomic systems, Control Decis., 22 (2007), 445–447+452.
    [34] X. Y. Meng, Y. Q. Wu, Bifurcation and control in a singular phytoplankton-zooplankton-fish model with nonlinear fish harvesting and taxation, Int. J. Bifurcation Chaos, 28 (2018), 1850042. doi: 10.1142/S0218127418500426
    [35] K. Chakraborty, M. Chakraborty, T. K. Kar, Bifurcation and control of a bioeconomic model of a prey-predator system with a time delay, Nonlinear Anal.: Hybrid Syst., 5 (2011), 613–625. doi: 10.1016/j.nahs.2011.05.004
    [36] Q. L. Zhang, C. Liu, X. Zhang, Complexity, Analysis and Control of Singular Biological Systems, London: Springer Press, 2012.
    [37] G. D. Zhang, Y. Shen, B. S. Chen, Hopf bifurcation of a predator-prey system with predator harvesting and two delays, Nonlinear Dyn., 73 (2013), 2119–2131. doi: 10.1007/s11071-013-0928-2
    [38] C. Liu, W. L. Ping, Q. L. Zhang, Y. Yan, Dynamical analysis in a bioeconomic phytoplankton zooplankton system with double time delays and environmental stochasticity, Phys. A, 482 (2017), 682–698. doi: 10.1016/j.physa.2017.04.104
    [39] B. Babaei, , M. Shafiee, Analysis and behavior control of a modified singular prey-predator model, Eur. J. Control, 49 (2019), 107–115. doi: 10.1016/j.ejcon.2019.01.001
    [40] W. M. Boothby, An Introduction to Differential Manifolds and Riemannian Geometry, New York: Academic Press, 1986.
    [41] B. S. Chen, X. X. Liao, Y. Q. Liu, Normal forms and bifurcations for the differential-algebraic systems, Acta Math. Appl. Sin., 23 (2000), 429–433.
    [42] B. D. Hassard, N. D. Kazarinoff, Y. H. Wan, Theory and Applications of Hopf Bifurcation, London: Cambridge University Press, 1981.
    [43] R. P. Gupta, M. Banerjee, P. Chandra, Period doubling cascades of prey-predator model with nonlinear harvesting and control of over exploitation through taxation, Commu. Nonlinear Sci. Numer. Simul., 19 (2014), 2382–2405. doi: 10.1016/j.cnsns.2013.10.033
    [44] M. Li, B. S. Chen, H. W. Ye, A bioeconomic differential algebraic predator-prey model with nonlinear prey harvesting, Appl. Math. Model., 42 (2017), 17–28. doi: 10.1016/j.apm.2016.09.029
    [45] P. D. N. Srinivasu, Bioeconomics of a renewable resource in presence of a predator, Nonlinear Anal.: Real. World Appl., 2 (2001), 497–506. doi: 10.1016/S1468-1218(01)00006-2
  • This article has been cited by:

    1. Anal Chatterjee, Samares Pal, A predator-prey model for the optimal control of fish harvesting through the imposition of a tax, 2023, 13, 2146-5703, 68, 10.11121/ijocta.2023.1218
    2. Xiuzhen Fan, Feng Zhou, Yan Li, Stationary pattern and Hopf bifurcation of a diffusive predator–prey model, 2021, 0003-6811, 1, 10.1080/00036811.2021.2021186
    3. Ruizhi Yang, Dan Jin, Wenlong Wang, A diffusive predator-prey model with generalist predator and time delay, 2022, 7, 2473-6988, 4574, 10.3934/math.2022255
    4. Md Golam Mortuja, Mithilesh Kumar Chaube, Santosh Kumar, Bifurcation analysis of a discrete type prey-predator model with Michaelis–Menten harvesting in predator, 2023, 0, 0932-0784, 10.1515/zna-2023-0022
    5. Salam Mohammed Ghazi Al-Mohanna, Yong-Hui Xia, Fear Effect on a Predator–Prey Model with Non-Differential Fractional Functional Response, 2023, 7, 2504-3110, 312, 10.3390/fractalfract7040312
    6. Shihua Ding, Rui Yang, Hopf and Turing–Hopf bifurcation analysis of a delayed predator–prey model with schooling behavior, 2023, 74, 0044-2275, 10.1007/s00033-023-02099-2
    7. Miao Peng, Rui Lin, Zhengdi Zhang, Lei Huang, The dynamics of a delayed predator-prey model with square root functional response and stage structure, 2024, 32, 2688-1594, 3275, 10.3934/era.2024150
    8. Liping Wu, Zhongyi Xiang, Dynamic analysis of a predator-prey impulse model with action threshold depending on the density of the predator and its rate of change, 2024, 9, 2473-6988, 10659, 10.3934/math.2024520
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2777) PDF downloads(202) Cited by(8)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog