
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
[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(1−XN)−α√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(1−X(s−τ1)N)−α√XY1+thα√X−EX,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=αYr√N,e=Er,t=rs. |
Thus system (1.2) can be rewritten as following:
{dxdt=x(1−x(t−τ1))−√xy1+a√x−ex,dydt=−dy2+b√x(t−τ2)y1+a√x(t−τ2),0=e(px−c)−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(1−x(t−τ1))−√xy−ex,dydt=−dy2+b√x(t−τ2)y,0=e(px−c)−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(1−x)−√xy−ex=0, | (2.3a) |
−dy2+b√xy=0, | (2.3b) |
e(px−c)−m=0. | (2.3c) |
From Eq (2.3b) and (2.3c), we can obtain y=bd√x and e=mpx−c, respectively. Substituting the above values into Eq (2.3a), we can obtain that x satisfies the following equation
x2+(bd−1−cp)x+cp+mp−bcdp=0. | (2.4) |
It is obvious that system Eq (2.4) has positive roots if the condition Δ=(bd−1−cp)2−4(cp+mp−bcdp)≥0. If Δ=0, then Eq (2.4) has a unique positive root, defined x∗=cd+pd−bp2dp. That is, system (2.2) has only one positive equilibrium P∗(x∗,y∗,e∗) if the assumption
(H1): p>c+bpd
holds, here y∗=bd√x∗ and e∗=mpx∗−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∗i(x∗i,y∗i,e∗i)i=1,2, here y∗i=bd√x∗i,e∗i=mpx∗i−c(i=1,2) and x∗1=cd+pd−bp+√Δ2dp>x∗,x∗2=cd+pd−bp−√Δ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≤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.
For simplicity, let
f(X)=(f1(x,y,e)f2(x,y,e))=(x(1−x(t−τ1))−√xy−ex−dy2+b√x(t−τ2)y),g(X)=e(px−c)−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=(100010−pe∗px∗−c01). |
Then, we get DYg(P∗)Q=(0,0,px∗−c) and u(t)=x(t),v(t)=y(t),E(t)=pe∗px∗−cx(t)+e(t). Substituting the latter into system (2.2), we can obtain
{du(t)dt=u(t)(1−u(t−τ1))−√u(t)v(t)−E(t)u(t)+pe∗u2(t)pu∗−c,dv(t)dt=−dv2(t)+b√u(t−τ2)v(t),0=(E(t)−pe∗u(t)pu∗−c)(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=(u∗v∗E∗), h(Z(t))=(0,0,h3(y1(t),y2(t))T:R2→R 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)))+pe∗pu∗−c(u∗+y1(t))2,dy2(t)dt=−d(v∗+y2(t))2+b√u∗+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=(1−b2d−u∗−E∗+2pe∗u∗pu∗−c)y1(t)−√u∗y2(t)−u∗y1(t−τ1),dy2(t)dt=−b√u∗y2(t)+b22dy1(t−τ2). | (2.7) |
The corresponding characteristic equation of the linearized system of the above system is:
λ2+(b√u∗−A)λ−Ab√u∗+u∗(λ+b√u∗)e−λτ1+b2√u∗2de−λτ2=0, | (2.8) |
where A=1−b2d−u∗+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ˊe's theorem.
Lemma 2.1. Consider the exponential polynomial
P(λ,e−λτ1,⋯,e−λτm)=λn+p01λn−1+⋯+p0n−1λ+p0n+[p11λn−1+⋯+p1n−1λ+p1n]e−λτ1+⋯+[pm1λn−1+⋯+pmn−1λ+pmn]e−λτm, |
where τi≥0(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∗+b√u∗−A)λ+b√u∗(u∗+b2d−A)=0. | (2.9) |
Thus, if u∗+b2d−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+mc(pu∗−c)2<min{2u∗+bd,2u∗+b2d+b√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.
Case one: τ1=0 and τ2≠0.
In this case, Eq (2.8) is translated into the following form
λ2+A11λ+A10+B10e−λτ2=0, | (2.10) |
where A11=b√u∗+u∗−A,A10=b√u∗(u∗−A),B10=b2√u∗2d.
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:
−ω2−A11ωi+A10+B10(cosωτ2−isinωτ2)=0. |
Separating the real and imaginary parts, we obtain
{B10cosωτ2=ω2−A10,B10sinωτ2=−A11ω, | (2.11) |
which yields that
ω4+(A211−2A10)ω2+A210−B210=0. | (2.12) |
Further, if denote v=ω2, we can obtain Φ1(ν)=ν2+(A211−2A10)ν+A210−B210. It is obvious that Eq (2.12) dose not have positive root when A210−B210>0 and A211−2A10≥0 or Δ1≜(A211−2A10)2−4(A211−2A10)<0.
If the following assumption
(H3): A210−B210≤0
holds, then Eq (2.12) has a positive root ω+2, here ω+2=√2A10−A211+Δ12.
If the following assumption
(H4): A210−B210>0, A211−2A10<0 and Δ1>0
holds, then Eq (2.12) has two positive roots ω+2 and ω−2, ω±2=√2A10−A211±Δ12.
From Eq (2.11), if (H3) holds and denote
τ+2j=1ω+2{arccos(ω+2)2−A10B10+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++A211−2A10)}=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)2−A10B10+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⋯τ−2k−1<τ+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 N∈N, 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: τ1≠0 and τ2=0.
Eq (2.8) takes the form
λ2+A21λ++A20+(B21λ+B20)e−λτ1=0, | (2.13) |
where A21=b√u∗−A,A20=b√u∗2d−Ab√u∗,B21=u∗,B20=bu∗√u∗.
Let ωi(ω>0) be the root of Eq (2.13). We obtain
{B21ωsinωτ1+B20cosωτ1=ω2−A20,B21ωcosωτ1−B20sinωτ1=−A21ω, | (2.14) |
which gives
ω4+(A221−B221−2A20)ω2+A220−B220=0. | (2.15) |
Further, if denote v=ω2, we can obtain Φ2(ν)=ν2+ϱν+A220−B220, here ϱ=A221−B221−2A20. It is obvious that Eq (2.15) dose not have positive root when A221−B221−2A20>0 and A220−B220>0 or Δ2≜(A221−B221−2A20)2−4(A220−B220)<0. We give two assumptions as following:
(H5): A220−B220≤0;
(H6): A221−B221−2A20<0, A220−B220>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(B20−A21B21)(ω±1)2−A20B20B220+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+A221−2A20(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 τ1≥0.
(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 N∈N, 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 τ1≠0 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 τ∗1∈I1s, then we can obtain
{b2√u∗2dcosωτ2=ω2+Ab√u∗−bu∗√u∗cosωτ∗1−ωu∗sinωτ∗1,b2√u∗2dsinωτ2=(b√u∗−A)ω−bu∗√u∗sinωτ∗1+ωu∗cosωτ∗1, | (2.16) |
which yields
ω4+(A2+b2u∗2)ω2+A2b2u∗−2Au∗(ω2+b2u∗)cosωτ∗1−2u∗(ω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+Ab√u∗−bu∗√u∗cosωiτ∗1−ωiu∗sinωiτ∗1)b2√u∗+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 PRQR−PIQI>0 hold, then the following transversality condition
dReλ(τ2)dτ2|τ2=τ(i)2j>0 |
is satisfied, where
PR=b√u∗−A+u∗[(1−bτ∗1√u∗)cosωiτ∗1−ωiτ∗1sinωiτ∗1],PI=2ωi−u∗[(1−bτ∗1√u∗)sinωiτ∗1+ωiτ∗1cosωiτ∗1],QR=ω2i+Ab√u∗−bu∗√u∗cosωiτ∗1−ωiu∗sinωiτ∗1,QI=ωi(b√u∗−A)−bu∗√u∗sinωiτ∗1+ωiu∗cosωiτ∗1. |
Proof. Taking the derivative of λ with respect to τ2 in Eq (2.8), we obtain
(dλdτ2)−1=2λ+b√u∗−A+u∗(1−τ∗1λ−bτ∗1√u∗)e−λτ∗1b2√u∗2dλ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λ+b√u∗−A+u∗(1−τ∗1λ−bτ∗1√u∗)e−λτ∗1b2√u∗2dλe−λτ2−τ2λ]λ=ω∗i}=4d2ω∗u∗b4(Q2R+Q2I)sign{PRQR−PIQI}. |
If PRQR−PIQI>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 PRQR−PIQI>0 hold. Then system (2.2) is locally asymptotically stable when τ2∈[0,τ2∗) and τ1∈I1s, 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 τ1∗∈I1s, 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∗,τ1∗∈I1s. Without loss of generality, we assume that τ1∗<τ2∗. Let x1(t)=y1(t)−y∗1,x2(t)=y2(t)−y∗2,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))T∈R2, and Lμ:C→R2,F:R×C→R2 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)=(A−√u∗0−b√u∗),C(τ2)=(−u∗000),D(τ2)=(00b2√u∗2d0), |
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(θ))T∈C([−1,0],R2),
a11=b4du∗+mp(pu∗−c)2−mp(pu∗+c)(pu∗−c)2,a12=−12√u∗,a13=−1,a111=−3b(8du∗)2−2mp2(pu∗−c)3+4mp2(pu∗+2c)(pu∗−c)4,a112=−14u∗√u∗,b11=−b24du∗,b12=b2√u∗,b22=−2d,b111=3b28d(u∗)2,b112=−b4u∗√u∗. |
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μ(ϕ)=∫0−1dη(θ,μ)ϕ(θ),ϕ∈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),∫0−1dη(θ,μ)ϕ(θ), θ=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)ds, s∈(0,1],∫0−1ψ(−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ω∗τ2∗q(θ),A∗(0)q∗(s)=−iω∗τ2∗q∗(s). |
By some simple computation, we can obtain
q(θ)=Veiω∗τ2∗θ,q∗(s)=DV∗e−iω∗τ2∗s, |
where
V=(1,ρ)T,ρ=1√u∗(A−u∗e−iω∗τ1∗−ie−iω∗τ2∗),V∗=(1,ρ∗)T,ρ∗=2d√u∗2dω∗i+b2e−iω∗τ2∗,ˉD=2d2d+2dρ¯ρ∗+b2τ2∗e−iω∗τ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)=e−iω∗τ2∗z+eiω∗τ2∗ˉz+W(1)(−1),x2t(−1)=ρe−iω∗τ2∗z+¯ρ∗eiω∗τ2∗ˉz+W(2)(−1),x1t(−τ1∗τ2∗)=e−iω∗τ1∗z+eiω∗τ1∗ˉz+W(1)(−τ1∗τ2∗),x2t(−τ1∗τ2∗)=ρe−iω∗τ1∗z+¯ρ∗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ρ+a13e−iω∗τ1∗+b11¯ρ∗e−iω∗τ2∗+b12¯ρ∗ρe−iω∗τ2∗+b22¯ρ∗ρ2],g11=ˉD[2a11+a12(ρ+ˉρ)+a13(eiω∗τ1∗+e−iω∗τ1∗)+2b11¯ρ∗+b12¯ρ∗(ρeiω∗τ2∗+ˉρe−iω∗τ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)e−iω∗τ1∗]+3a111+a112(2+ˉρ)+b11¯ρ∗[2W(1)11(−1)e−iω∗τ2∗+W(1)20(−1)eiω∗τ2∗]+b12¯ρ∗[W(2)11(0)e−iω∗τ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¯ρ∗e−iω∗τ2∗+b11¯ρ∗(2ρ+e−2iω∗τ2∗)}. | (3.4) |
Since W20(θ) and W11(θ) appears in g21, we need to compute them as follows.
W20(θ)=ig20q(0)ω∗τ2∗eiω∗τ2∗θ+iˉg02ˉq(0)3ω∗τ2∗e−iω∗τ2∗θ+E1e2iω∗τ2∗θ,W11(θ)=−ig11q(0)ω∗τ2∗eiω∗τ2∗θ+iˉg11ˉq(0)ω∗τ2∗e−iω∗τ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+u∗e−2iω∗τ2∗−√u∗02iω∗+b√u∗−b22de−2iω∗τ2∗)−1(K11K21), |
and
E2=−(A−u∗−√u∗b2√u∗2db√u∗)−1(K12K22), |
here
K11=a11+a12ρ+a13e−iω∗τ1∗,K21=2a11+a12(ρ+ˉρ)+a13(eiω∗τ1∗+e−iω∗τ1∗),K12=b11e−iω∗τ2∗+b12ρe−iω∗τ2∗+b22ρ2,K22=2b11+b12(ρeiω∗τ2∗+ˉρe−iω∗τ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∗(g20g11−2|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)=(1−b2d−u∗−E∗+2pe∗u∗pu∗−c)y1(t)−√u∗y2(t)−u∗y1(t−τ1),.y2(t)=−b√u∗y2(t)+b22dy1(t−τ2)+k1y2(t)+k2y32(t). | (4.2) |
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
(H8)k1>A10A−u∗and(A10−k1(u∗−A))−(A11−k1)2<B10ek1−A112τ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+(A11−k1)λ+(A10−k1(u∗−A))+B10e−λτ2=0, | (4.3) |
where A11,A10 and B10 are the same to that in Eq (2.11)
Let g1(λ)=λ2+(A11−k1)λ+(A10−k1(u∗−A)) 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>A10A−u∗ and (A10−k1(u∗−A))−(A11−k1)2<B10ek1−A112τ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(1−x)−√xy−mx(p1−T)x−c,.y(t)=−dy2+b√xy. | (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((p1−T)x−c)+exT=e(p1x−c). |
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(p1x−c)e−δt+λ1(x(1−x)−√xy−mx(p1−T)x−c)+λ2(−dy2+b√xy). |
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,∂H∂T>0,0,∂H∂T<0. | (5.4) |
It is assumed that the optimal solution does not occur at T=Tmax 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
dλ1dt=−∂H∂xanddλ2dt=−∂H∂y. | (5.5) |
In order to obtain singular optimal equilibrium solution, we use steady state equations
{(1−x)−√xyx−x(p1−T)x−c=0,−dy+b√x=0. | (5.6) |
Thus, Eq (5.5) along with steady state Eq (5.6) gives
dλ1dt=ep1e−δt+λ2by2√x. | (5.7) |
From Eq (5.7) we get
λ2(t)=2ep1e−δtby√x. | (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(2dy−b√x)by√x. | (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(2dy−b√x)δby√x. |
From Eqs (5.8) and (5.10) we get
2ep1by√x−2ep1(2dy−b√x)δby√x=0. | (5.11) |
That is, 2dy−b√x−δ=0.
According to the interior equilibrium (x∗,y∗), we can get T=b2cd+p1(b2(d−b)−2dδ2)b2(d−b)−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): p−c−bpd=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∈I1s=(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.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.
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).
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 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:
0=qEsC−10+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
![]() |
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 |