# Complex Dynamics in a Singular Delayed Bioeconomic Model with and without Stochastic Fluctuation

Authors: Xinyou Meng; Qingling Zhang
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

A singular delayed biological economic predator-prey system with and without stochastic fluctuation is proposed. The conditions of singularity induced bifurcation are gained, and a state feedback controller is designed to eliminate such bifurcation. Furthermore, saddle-node bifurcation is also showed. Next, the local stability of the positive equilibrium and the existence of Hopf bifurcation are obtained by analyzing the distribution of roots of the corresponding characteristic equation, and the hybrid control strategy is used to control the occurrence of Hopf bifurcation. In addition, some explicit formulas determining the spectral densities of the populations and harvest effort are given when the system is considered with stochastic fluctuation. Finally, numerical simulations are illustrated to verify the theoretical results.


---

## Body

## 1. Introduction

The dynamic relationship between predator and prey has long been and will still be one of the dominant themes in both biology and mathematical biology because of its universal existence and importance. In the description of dynamics interactions, a crucial element of all models is the classical definition of functional response. Lots of predator-prey models with Holling type [1], Leslie-Gower type [2], and Beddington-DeAngelis type [3, 4], and so forth have been investigated extensively by scholars. However, some predator-prey models in which prey population exhibits herd behavior, such as plankton-phytoplankton model [5], may appear in realistic world. Since the use of the square root properly accounts for the assumption that the interactions occur along the boundary of the population, Ajraldi et al. [6] proposed the following predator-prey system in which interaction terms use the square root of prey population rather than simply prey population:(1)dRtdt=rRt1-RtN-aRtSt,dStdt=-m~St+ae~RtSt,where R and S denote the prey and predator, respectively. The prey population exhibits a highly socialized behavior and lives in herds as the form aR(t)S(t); that is, the weaker individuals are being kept at the center of their herd for defensive purpose. Braza [7] also investigated the dynamics of system (1) and showed that the prey exhibits strong herd structure and the predator interacts with the prey along the outer corridor of the herd of prey.Recently, Yuan et al. [8] considered a predator-prey system as follows:(2)dXdt=rX1-XN-αXY1+thαX,dYdt=-sY2+cαXY1+thαX,where -sY2 represents the quadratic mortality for predator population. Predator-prey systems with such functional response have attracted little attention (see [6–9]).It is well-known that 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. Delay differential equations often show much more complicated dynamics than ordinary differential equations because a time delay can cause a stable equilibrium to become unstable and cause the population to fluctuate. Many authors have been devoted to investigating the time delay effect on the dynamics of system and obtained some results (see [10–15]). Considering the fact that there always exists a time delay in the conversion of the biomass of prey to that of predator in system (2), Xu and Yuan [9] introduced a time delay into system (2) and obtained the local stability and the existence of Hopf bifurcation of this system. However, bifurcate oscillation is harmful in some engineering applications, which has enormous potential in many technological disciplines such as power networks protection. Bifurcation control, which refers to the aim of designing a controller to suppress or reduce some existing bifurcation dynamics of a given system, can be useful. Various disciplines are attracted to bifurcation control and various methods of bifurcation control can be found [16–18].In addition, Gordon [19] proposed the economic theory of a common-property resource, which focuses on the effect of the harvest effort on the ecosystem from an economic perspective. If E(t) and Y(t) represent the harvest effort and the density of harvested population, respectively, then the total revenue TR=ωE(t)Y(t) and the total cost TC=cE(t), where ω represents unit price of harvested population and c represents the cost of harvest effort. Thus, an algebraic equation, which considers the economic interest v of the harvest effort on the harvested population, is established as follows: (3)EtωYt-c=v. Based on the Gordon [19] theory and theory of singular system, Zhang et al. [20] first proposed a class of singular biological economic systems. Some results on those systems, such as the stabilities, bifurcations, and chaos, can be found in [20–24].Based on the previous models, we establish the following predator-prey system consisting of two differential equations and an algebraic equation as follos:(4)dXds=rX1-XN-αXY1+thαX-EX,dYds=-d~Y2+b~αXs-τY1+thαXs-τ,0=Ep~X-c~-m~,where X, Y, and E are the prey, predator, and the harvest effort at the time s, respectively. τ denotes the time delay which means the growth rate of predator species depends on the number of the prey species τ units of time earlier. The parameters 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. -d~Y2 represents the quadratic mortality for predator population. p~ and c~ are unit price of harvested population and the cost of harvest effort. m~ is harvest interest of the harvest effort on the harvested prey.It is important to make some simplifying assumptions to discern the basic dynamics and to make the analysis more tractable. In order to simplify system (4), we use the following dimensionless transformations: x=X/N, y=αY/rN, e=E/r,  and t=rs. Then system (4) can be rewritten as follows:(5)dxdt=x1-x-xy1+ax-ex,dydt=-dy2+bxt-τy1+axt-τ,0=epx-c-m,where a=tnαN, b=αNb~/r, c=c~, d=Nd~/r, m=m~/r, and p=p~N.In the papers [6–9], the author used the simplifying assumption that a=0; that is, the average handling time is zero. In line with the work in [6–9], we also assume that a=0. Thus, system (5) takes the form:(6)dxdt=x1-x-xy-ex,dydt=-dy2+bxt-τy,0=epx-c-m. The initial conditions of system (6) are(7)xθ=ψθ≥0,yθ=ηθ≥0,eθ=ϕθ≥0,x0>0,y0>0,e0>0,θ∈-τ,0.In reality, the environmental fluctuation is one of the important components for ecological systems. Many natural phenomena do not follow the deterministic law and usually oscillate randomly around some average values. The deterministic approach has limitations on mathematical modeling, which makes the accurate prediction for the future dynamics of system very difficult. Stochastic differential equation models play a significant role in various dynamic analysis, because they can provide some additional degree of realism compared to their deterministic counterpart [25]. Recently, some results of the stochastic modeling of ecological populations are presented [26–30].In order to study the effects of the environmental fluctuations on a delayed singular prey-predator bioeconomic model, the following stochastic model corresponding to the delayed system (4) in the fluctuating environment is given:(8)dXds=rX1-XN-αXY1+thαX-EX+ξ~1t,dYds=-d~Y2+b~αXs-τY1+thαXs-τ+ξ~2t,0=Ep~X-c~-m~+ξ~3t. The perturbed terms ξ~i, i=1,2,3, are mutually independent Gaussian white noises characterized by 〈ξ~i〉 and 〈ξ~i(t)ξ~j(t1)〉=δijδ(t-t1), i,j=1,2,3. The symbol 〈·〉 is the ensemble average due to the effect of the fluctuating environment, δij is the Kronecker delta and represent the spectral density of the white noise, and δ is the Dirac delta function with t and t1 being the distinct times.By taking the same translations and same notations, we can obtain the following form:(9)dxdt=x1-x-xy-ex+ξ1t,dydt=-dy2+bxt-τy+ξ2t,0=epx-c-m+ξ3t, where ξ1(t)=ξ~1(t)/N, ξ2(t)=ξ~2(t)/rN, and ξ3(t)=ξ~3(t)/r. In addition, the initial conditions of system (9) are also similar to that of system (6).The rest of paper is organized as follows. In the next section, some conditions for the existence of the positive equilibrium, saddle-node bifurcation, singularity induced bifurcation, and Hopf bifurcation are obtained and bifurcation controls in deterministic model are also showed. Some explicit formulas determining the spectral densities of the populations and harvest effort in the singular bioeconomic system with the fluctuation environment are gotten in Section3. To support our theoretical predictions, some numerical simulations are included in Section 4. A brief discussion is also given in the last section.


## 2. Dynamics of the Deterministic Model

### 2.1. Existence of the Positive Equilibrium

From the viewpoint of biological interpretation, we only consider the positive equilibrium. There exists a positive equilibriumP(x,y,e) of system (6), where the values x,y, and e satisfy the following equations: (10a)x1-x-xy-ex=0,(10b)-dy2+bxy=0,(10c)epx-c-m=0.From (10b) and (10c), we can obtain y=b/dx and e=m/px-c, respectively. Substituting the above values into (10a), we can obtain that x satisfies the following equation:(11)x2+bd-1-cpx+cp+mp-bcdp=0.By simple computation, we obtain that(12)msn=bcd-c+p4bd-1-cp2. Suppose that msn>0 and(1)
m>msn: there is no positive equilibrium;(2)
m=msn: there is only one positive equilibrium P∗(x∗,y∗,e∗) if the assumption
(H1): p>c+bp/dholds, wherey∗=b/dx∗, e∗=m/px∗-c, and x∗=cd+pd-bp/2dp;(3)
m<msn: there are two positive equilibriums P∗+(x∗+,y∗+,e∗+) and P∗-(x∗-,y∗-,e∗-), where y∗±=b/dx∗±, e∗±=m/px∗±-c, x∗±=cd+pd-bp±Δ/2dp, and Δ=(b/d-1-c/p)2-4(c/p+m/p-bc/dp).Remark 1.
In the following, we only investigate that system (6) has a unique positive equilibrium P∗(x∗,y∗,e∗). Of course, we can discuss the stability and bifurcation at the equilibrium P∗+(x∗+,y∗+,e∗+) (resp., P∗-(x∗-,y∗-,e∗-)) by the same way.


### 2.2.  Singularity Induced Bifurcation and Control

In this part, our main objective is to study the occurrence of singularity induced bifurcation and the effects of economic profit on system (6). System (6) without time delay takes the following form:(13)dxdt=x1-x-xy-ex,dydt=-dy2+bxy,0=epx-c-m. If the parameter m is taken as bifurcation parameter, the singularity induced bifurcation, which is first introduced by Venkatasubramanian et al. [31], may appear around the interior equilibrium. The singularity induced bifurcation does not occur in a normal ordinary differential equation system, which has been characterized by a singular system. Roughly speaking, a singularity induced bifurcation refers to a stability change of the singular system, which leads to an impulse phenomenon of the singular system and may even result in the collapse of the system (see [32, 33]). Thus, we give some results on the singularity induced bifurcation of singular system.Singular system, renamed as differential-algebraic system, can be described by the following form:(14)x˙=fx,y,μ,f:Rm+n+q⟶Rn,0=gx,y,μ,g:Rm+n+q⟶Rm,where x∈X⊂Rn and y∈Y⊂Rm denote differential parts and algebraic parts of state variable, respectively. μ∈Λ⊂Rq is the parametric variable, and f and g are vector functions with appropriate dimensions.Lemma 2.
We will consider system (14) with a one-dimensional parameter space. If we suppose the following conditions are satisfied at (0,0,μ0), where D is a differential operator, Δ=DEg is the matrix of partial derivatives of the components of g with respect to E, and Dyg is the matrix of partial derivatives of the components of g with respect to y:SI1:f(x,y,μ)=0,g(x,y,μ)=0, Dyg has an algebraically simple zero eigenvalue and trace[Dyfadj⁡(Dyg)(Dxg,Dyg)] is nonzero,SI2:DxfDyfDxgDyg is nonsingular,SI3:DxfDyfDμfDxgDygDμgDxΔDyΔDμΔ is nonsingular, and rankDxfDyfDμfDxgDygDpgDxΔDyΔDμΔ=rank⁡f(x,y,μ)+rank⁡g(x,y,μ),then there exists a smooth curve of equilibrium in Rn+m+1, which passes through (0,0,μ0) and is transversal to the singular surface at (0,0,μ0). When μ increases through μ0, one eigenvalue of this system moves from C- (the open complex left half plane) to C+ (the open complex right half plane) if M/N>0 (resp., from C+ to C- if M/N<0) along the real axis by diverging through the infinity. The other (n-1) eigenvalues remain bounded and stay away from the origin. The constants M and N can be computed: (15)M=-trace⁡Dyfadj⁡DygDxg,Dyg0,0,μ0,N=DμΔ-DxΔDyΔDxfDyfDxgDyg-1DμfDμg0,0,μ0.For simplicity, let(16)fXt,Et,μ=f1f2=x1-x-xy-ex-dy2+bxy,gXt,Et,μ=epx-c-m,where X(t)=(x(t),y(t))T, E(t)=e(t), μ=m. Let D be a differential operator: Δ(X(t),E(t),μ)=DEg(X(t),E(t),μ)=px(t)-c. Applying Lemma 2 to system (13), we can obtain the following result.Theorem 3.
If (H1) holds, then system (13) has a singularity induced bifurcation at the interior equilibrium P∗(x∗,y∗,e∗) when the bifurcation parameter m increases through zero. That is, a stability switch occurs as parameter m passes through zero.Proof.
It is obvious thatDEgP∗(x∗,y∗,e∗)=0 and m=0, which implies that DEg has a simple zero eigenvalue. Here P∗(x∗,y∗,e∗)≐(p/c,b/dp/c,1-b/d-p/c).
By simple computation, we have that(17)trace⁡DEfadj⁡DEgDXg,DEgP∗=-px∗e∗≠0. Furthermore, it can be calculated that (18)DXfDEfDXgDEgP∗=1-2x-y2x-e-x-xby2x-2dy+bx0pe0px-cP∗=-pdx∗y∗e∗≠0,which follows that DXfDEfDXgDEg is nonsingular at P∗.
By simple computing,(19)DXfDEfDμfDXgDEgDμgDXΔDEΔDμΔP∗=-x+y2x-x-x0by2x-dy00pe0px-c0p000P∗=pdx∗y∗≠0;it is easy to show that rank⁡(f(X,E,μ))=2 and rank⁡(g(X,E,μ))=1. Hence, rankDXfDEfDμfDXgDEgDμgDXΔDEΔDμΔ=rank⁡(f(X,E,μ))+rank⁡(g(X,E,μ))+1.
Thus, all conditions of Lemma2 are satisfied. System (13) has a singularity induced bifurcation at the positive equilibrium P∗(x∗,y∗,e∗) when the bifurcation parameter m=0.
Furthermore, the constantsM and N can be computed: (20)M=-trace⁡DEfadj⁡DEgDXg,DEgP∗x∗,y∗,e∗=px∗e∗≠0,N=DμΔ-DXΔ,DEΔDXfDEfDXgDEg-1DμfDμgP∗x∗,y∗,e∗=-1e∗≠0,which follows that (21)MN=-px∗e∗2<0. If parameter m increases through zero, then one eigenvalue of system (13) will move from C+ to C- along the real axis by diverging into the infinity. Thus, the stability of the positive equilibrium P∗(x∗,y∗,e∗) changes from unstable to stable.
The Jacobian of system (13) evaluated at P∗ takes the form: (22)J=∂f,g∂X,E=1-2x-y2x-e-x-xby2x-dy0pe0px-c. According to the leading matrix Ξ in system (13) and JP∗, here Ξ=100010000, the characteristic equation of system (13) at P∗ is det(λΞ-JP∗)=0. That is, the characteristic equation can be expressed: (23)det⁡λΞ-JP∗=px∗e∗λ+bx∗=0. It is obvious that the rest eigenvalue (denoted by λ2) has negative real part. Furthermore, it follows from Theorem 1.1 [34] that there is only one eigenvalue diverging to infinity, but the rest eigenvalue is continuous, nonzero, and cannot jump from on half open complex plane to another one as m increases through 0. Table 1 shows the change in the signs of the real parts of eigenvalues (λ1 and λ2) due to the variation of economic interest of harvest effort. According to Table 1, it can be concluded that system (13) is stable at P∗ as m<0 and unstable as m>0. Consequently, a stability switch occurs as m increases through 0.Table 1
Signs of real parts of eigenvalues of system (13) at P∗.
Parameter value
R
e
λ
1
R
e
λ
2
m
<
0
−
−
m
>
0
+
−Remark 4.
According to [35], system (13) along the positive equilibrium locus yields index one matrix pencil (Ξ,J) for m≠0 and index two matrix pencil (Ξ,J) for m=0. This shows that there exists one index jump at a bifurcation point m=0. It leads to the rapid expansion of the population from the point of view of biology. If this phenomenon lasts for a long time, the population will be beyond the carrying capacity of the environment and the prey-predator system will be out of balance, which is disastrous. With the purpose of economic interest of harvest effort at an ideal level as well as maintaining the sustainable development of the biological resource, some related measures should be taken to eliminate the impulse phenomenon and stabilize system (13) when the economic interest is positive. Thus, a state feedback controller is designed to stabilize system (13) at the positive equilibrium P∗.By using Theorem 3-1.2 [36], a state feedback controller u(t)=k(e(t)-e∗), where k is a feedback gain and e∗ is the component of the positive equilibrium P∗, can be applied to stabilize system (13) at P∗. Thus, a controlled system is as follows:(24)dxdt=x1-x-xy-ex,dydt=-dy2+bxy,0=epx-c+ke-e∗-m,where all variables and parameters have the same interpretations as that in system (13). The feedback gain can be given in the following result.Theorem 5.
Suppose the feedback gaink satisfies one of the following cases:(i)
ifbp-2cd>0, then k<2p(cd-bc-pd)/bp-2cd;(ii)
ifbp-2cd<0, then k>2p(cd-bc-pd)/bp-2cd;the controlled system (24) is stable at the positive equilibrium P∗.Proof.
The Jacobian of system (24) evaluated at the interior equilibrium P∗ has the form (25)J~P∗=1-2x∗-b2d-e∗-x∗-x∗b22d-bx∗0pe∗0k. According to the leading matrix Ξ in system (13) and J~P∗, the characteristic equation of the controlled system (24) at P∗ is (26)λ2+bx∗+px∗e∗k-1+2x∗+b2d+e∗λ+bx∗px∗e∗k-1+2x∗+b2d+e∗=0. By using the Routh-Hurwitz criteria [37], the sufficient and necessary condition for the stability of the controlled system (24) at P∗ is that the feedback gain k satisfies the following: if bp-2cd>0, then k<2p(cd-bc-pd)/bp-2cd; otherwise, k>2p(cd-bc-pd)/bp-2cd. This ends the proof.Remark 6.
After introducing the feedback controller into system (13), the controlled system (24) can be stabilized at the positive equilibrium. The elimination of the singularity induced bifurcation means the prey-predator system restores ecological balance. Both sustainable development of the prey-predator system and the ideal economic interest of harvesting can be obtained by enhancing or reducing the harvest effort on the prey.


### 2.3. Saddle-Node Bifurcation

Ifm is taken as bifurcation parameter, the following result is true.Theorem 7.
If the assumption (H1) holds, then system (13) undergoes saddle-node bifurcation at m=msn.Proof.
The characteristic polynomial of system (13) at the positive equilibrium P∗(x∗,y∗,e∗) is given: (27)Fλ=λ2+l1λ+l0,where l1=bx∗+2x∗+b/2d-1-mc/(px∗-c)2 and l0=bx∗(2x∗+b/d-1-mc/(px∗-c)2). Clearly, it has a simple zero eigenvalue with left eigenvector usn=(-b,1) and right eigenvector vsn=(x∗,b/2d)T. And (28)usnDμfRP∗=usnDμf-DEfDEg-1DμgRP∗=bx∗px∗-c,usnDX2fRvsn,vsnP∗=usn∑i=1i=2eivsnTDXDXfiTvsnP∗=b22d+2bx∗+by∗4x∗b-1+2bmpcx∗px∗-c3,where ei(i=1,2) are unit vectors.
According to [31], system (13) undergoes saddle-node bifurcation when m=msn holds. This completes the proof.


### 2.4. Hopf Bifurcation and Control

From the discussion above, we know that if (H1) holds, there is only one positive equilibrium P∗(x∗,y∗,e∗); if m<msn, there are two positive equilibriums P∗+(x∗+,y∗+,e∗+) and P∗-(x∗-,y∗-,e∗-). In this subsection, we only investigate the dynamical behavior of system (6) at the positive equilibrium P∗+(x∗+,y∗+,e∗+) but neglect similar results on the other positive equilibria.For simplicity, let(29)fW=f1x,y,ef2x,y,e=x1-x-xy-ex-dy2+bxt-τy,gW=epx-c-m,where W=(x,y,e)T. In order to analyze the local stability of the positive equilibrium of system (6), we first use the linear transformation WT(t)=QNT(t), where (30)Nt=ut,vt,Et,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 (6), we can obtain(31)dutdt=ut1-u-utvt-Etut+pe∗+u2tpx∗+-c,dvtdt=-dv2t+but-τvt,0=Et-pe∗+utpx∗+-cput-c-m. Further, in order to derive the formula determining the properties of the positive equilibrium of system (31), we consider the local parametric Ψ of the third equation of system (31) as the literatures [38, 39], which is given: (32)NTt=ΨZt=N0T+u0Zt+v0hZt,gΨZt=0,where(33)u0=100100,v0=001,Zt=y1ty2t,N0=u∗v∗E∗=x∗+y∗+px∗+e∗+px∗+-c+e∗+,hZt=00h3y1t,y2t:R2⟶Ris a smooth mapping. The corresponding characteristic equation of the linearized system of the above system is(34)λ2+bu∗-Aλ-Abu∗+b22du∗e-λτ=0,where A=1-b/2d-u∗+mc/(pu∗-c)2.In order to investigate the distribution of roots of the transcendental equation (34), we introduce the following result proved by Ruan and Wei [10].Lemma 8.
Consider the exponential polynomial:(35)Pλ,e-λτ1,…,e-λτm=λn+p10λn-1+⋯+pn-10λ+pn0+p11λn-1+⋯+pn-11λ+pn1e-λτ1+⋯+p1mλn-1+⋯+pn-1mλ+pnme-λτm,where τi≥0(i=1,2,…,m) and pji (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 (34). When τ=0, (34) is rewritten as follows:(36)λ2+bu∗-Aλ+bu∗b2d-A=0.Thus, if b/2d-A<0, (36) has at least one positive root, which implies that system (6) without any time delay is unstable. If the condition(H2): 1+mc/(pu∗-c)2<min⁡{b/d+2u∗,b/2d+2u∗+bu∗}holds, two roots of (36) have always negative real parts. Thus, we have the following result.Lemma 9.
Assume that (H1) and (H2) hold, then the two roots of (36) have always negative real parts; that is, system (6) with τ=0 is locally asymptotically stable.Letωi(ω>0) be the root of (34), and we have the following equations by substituting it into (34) and separating the real and imaginary parts:(37)b2u∗2dcos⁡ωτ=ω2+Abu∗,b2u∗2dsin⁡ωτ=bu∗-Aω,which gives(38)ω4+a1ω2+a0=0,where a1=A2+b2u∗>0 and a0=b2u∗(A2-b2/4d2).Letω2=θ; (38) takes the form(39)Φθ=θ2+a1θ+a0=0.It is obvious that (39) has no real roots when Δ=a12-4a0<0. When a0>0, (39) has two negative roots; when Δ>0 and a0<0, (39) has one positive root. Thus, we have the following results.Lemma 10.
For (38), the following results are true:(i)
whena0>0, (38) does not have any positive real root;(ii)
whenΔ>0 and a0<0, (38) has only one positive root, defined by ω+=-a1+Δ/2.If condition (ii) of Lemma10 holds, and we denote (40)τj+=1ω+arccos⁡2dAbu∗+ω+2b2u∗+2jπ,j=0,1,2,…,then ±iω+ are a pair of purely imaginary roots of (34) with τ=τj+. Next, we will check whether the following transversality conditions are satisfied.Lemma 11.
Suppose that (H1) holds, the transversality condition (41)dRe⁡λτdττ=τj+>0is satisfied, and (dλ/dτ)-1 and Φ′(θ) have the same sign.Proof.
Differentiating the two sides of (34) with respect to τ, we get (42)dλdτ-1=2λ+bu∗-Aλb2u∗/2de-λτ-τλ.
For simplify, we defineω0 as ω and τj+ as τ, and we can obtain (43)sign⁡dRe⁡λdττ=τj+=sign⁡Re⁡dλdτ-1τ=τj+=sign⁡2ω2+bu∗-A2+2Abu∗bu∗-A2ω2+ω2+Abu∗2=1bu∗-A2ω2+ω2+Abu∗2sign⁡Φ′θ. This completes the proof.Summarizing the above lemmas and applying Lemma11 to (6), we have the following theorem.Theorem 12.
Suppose that (H1) holds, the following results are true.(1)
Whena0>0, system (6) is locally asymptotically stable for any τ≥0.(2)
WhenΔ>0 and a0<0, system (6) is locally asymptotically stable for τ∈[0,τ0) and unstable when τ>τ0 and system (6) undergoes a Hopf bifurcation at the equilibrium P∗+ when τ=τj(j=0,1,2,…).It has been realized that bifurcation control can be useful when bifurcation oscillation is harmful in some engineering applications, which has great potential in many technological disciplines. The aim of bifurcation control is to delay the onset of an inherent bifurcation, change the parameter value of an existing bifurcation point, stabilize a bifurcated solution or branch, and so forth. In order to control the Hopf bifurcation to achieve some desirable behaviors, a new hybrid control strategy was proposed in [40]. Recently, some scholars succeed in applying this method in various fields [40–43]. However, as far as we known, there is little attention on Hopf bifurcation of the singular biological economic system by use of this method. Extending this method on singular system is first time.By choosing the hybrid control strategy, the controlled model with delay is as follows:(44)dxdt=x1-x-xy-ex,dydt=α-dy2+bxt-τy+1-αy-y∗+,0=epx-c-m,where α(0<α<1) is the hybrid control parameter. The controlled system (44) reduces to the original system (6) if α=1. By selecting appropriate control parameter α, the Hopf bifurcation can be delayed (or advanced) or even eliminated completely without changing the equilibrium of the original system.The controlled system (44) has the same equilibrium as system (6). Let w1(t)=x(t)-x∗+,w2(t)=y(t)-y∗+, and w3(t)=e(t)-e∗+. The corresponding characteristic equation of the linear system of the controlled system (44) has the form (45)px∗+-cλ2+b11λ+b10+b00e-λτ=0,where b11=pe∗+x∗+/px∗+-c-1+2x∗++b/2d+e∗++bx∗+α-1+α,  b10=(-bx∗+α+1-α)(1-2x∗+-b/2d-e∗+-pe∗+x∗+/px∗+-c), and  b00=-αb2x∗+/2d<0. Therefore, in order to study the local stability of the equilibrium of system (44), we will study the following equation:(46)λ2+b11λ+b10+b00e-λτ=0.Whenτ=0, (46) is given:(47)λ2+b11λ+b10+b00=0.It is obvious that the positive equilibrium P∗ is locally asymptotically stable when τ=0 if the following condition(H3) b11>0 and b10+b00>0holds.Letωi(ω>0) be a root of (46). We have that(48)b00cos⁡ωτ=ω2+b10,b00sin⁡ωτ=-b11ω,which yields to(49)ω4+b112-2b10ω2+b102-b002=0.By simple computation,b112-2b10=(pe∗+x∗+/px∗+-c-1+2x∗++b/2d+e∗+)2>0 is obtained. If b102-b002>0, then (49) does not have any positive real root. If b102-b002≤0, then (49) has only one positive real root ω~+=1/2[2b10-b112+(b112-2b10)2-4(b102-b002)]. From (48), we can get (50)τ~j+=1ω~+arccos⁡ω~+2-b10b00+2jπ,j=0,1,2,…,and then ±iω~+ are a pair of purely imaginary roots of (46) with τ=τ~j+. Define τ~0=min⁡{τ~j+,j=0,1,2,…}. Further, if λ(τ)=α(τ)+iω(τ) is denoted by the root of (46), then it satisfies α(τ~0)=0, ω(τ~0)=ω~0.Next, we will check whether the following transversality condition is satisfied. Differentiating the two sides of (46) with respect to τ and applying the implicit theorem, we get (51)dλdτ-1=2λ+b11λb00e-λτ-τλ.For simplify, we defineω~+ as ω and τ~j+ as τ, and we can obtain (52)sign⁡dRe⁡λdττ=τ~j+=sign⁡Re⁡dλdτ-1τ=τ~j+=sign⁡Re⁡2λ+b11λb00e-λτλ=ω~+i=sign⁡2ω2+b112-2b10b112ω2+ω2+b102.Since b112-2b10>0, then the transversality condition (53)dRe⁡λτdττ=τ~j+>0is satisfied.Summarizing the above results, we have the following theorem.Theorem 13.
If (H3) holds, then we have the following results on controlled system (44).(1)
Ifb11>0 and b10+b00>0 or b102-b002>0, then the positive equilibrium of controlled system (44) is locally asymptotically stable for any τ≥0.(2)
Ifb102-b002≤0, the positive equilibrium of controlled system (44) is locally asymptotically stable for τ∈[0,τ~0) and unstable when τ>τ~0, and controlled system (44) undergoes a Hopf bifurcation at the equilibrium P∗ when τ=τ~0.Remark 14.
For the controlled system (44), it has been proven in Theorem 13 that we can delay the onset of Hopf bifurcation without changing the original equilibrium by choosing appropriate parameter.


## 2.1. Existence of the Positive Equilibrium

From the viewpoint of biological interpretation, we only consider the positive equilibrium. There exists a positive equilibriumP(x,y,e) of system (6), where the values x,y, and e satisfy the following equations: (10a)x1-x-xy-ex=0,(10b)-dy2+bxy=0,(10c)epx-c-m=0.From (10b) and (10c), we can obtain y=b/dx and e=m/px-c, respectively. Substituting the above values into (10a), we can obtain that x satisfies the following equation:(11)x2+bd-1-cpx+cp+mp-bcdp=0.By simple computation, we obtain that(12)msn=bcd-c+p4bd-1-cp2. Suppose that msn>0 and(1)
m>msn: there is no positive equilibrium;(2)
m=msn: there is only one positive equilibrium P∗(x∗,y∗,e∗) if the assumption
(H1): p>c+bp/dholds, wherey∗=b/dx∗, e∗=m/px∗-c, and x∗=cd+pd-bp/2dp;(3)
m<msn: there are two positive equilibriums P∗+(x∗+,y∗+,e∗+) and P∗-(x∗-,y∗-,e∗-), where y∗±=b/dx∗±, e∗±=m/px∗±-c, x∗±=cd+pd-bp±Δ/2dp, and Δ=(b/d-1-c/p)2-4(c/p+m/p-bc/dp).Remark 1.
In the following, we only investigate that system (6) has a unique positive equilibrium P∗(x∗,y∗,e∗). Of course, we can discuss the stability and bifurcation at the equilibrium P∗+(x∗+,y∗+,e∗+) (resp., P∗-(x∗-,y∗-,e∗-)) by the same way.


## 2.2.  Singularity Induced Bifurcation and Control

In this part, our main objective is to study the occurrence of singularity induced bifurcation and the effects of economic profit on system (6). System (6) without time delay takes the following form:(13)dxdt=x1-x-xy-ex,dydt=-dy2+bxy,0=epx-c-m. If the parameter m is taken as bifurcation parameter, the singularity induced bifurcation, which is first introduced by Venkatasubramanian et al. [31], may appear around the interior equilibrium. The singularity induced bifurcation does not occur in a normal ordinary differential equation system, which has been characterized by a singular system. Roughly speaking, a singularity induced bifurcation refers to a stability change of the singular system, which leads to an impulse phenomenon of the singular system and may even result in the collapse of the system (see [32, 33]). Thus, we give some results on the singularity induced bifurcation of singular system.Singular system, renamed as differential-algebraic system, can be described by the following form:(14)x˙=fx,y,μ,f:Rm+n+q⟶Rn,0=gx,y,μ,g:Rm+n+q⟶Rm,where x∈X⊂Rn and y∈Y⊂Rm denote differential parts and algebraic parts of state variable, respectively. μ∈Λ⊂Rq is the parametric variable, and f and g are vector functions with appropriate dimensions.Lemma 2.
We will consider system (14) with a one-dimensional parameter space. If we suppose the following conditions are satisfied at (0,0,μ0), where D is a differential operator, Δ=DEg is the matrix of partial derivatives of the components of g with respect to E, and Dyg is the matrix of partial derivatives of the components of g with respect to y:SI1:f(x,y,μ)=0,g(x,y,μ)=0, Dyg has an algebraically simple zero eigenvalue and trace[Dyfadj⁡(Dyg)(Dxg,Dyg)] is nonzero,SI2:DxfDyfDxgDyg is nonsingular,SI3:DxfDyfDμfDxgDygDμgDxΔDyΔDμΔ is nonsingular, and rankDxfDyfDμfDxgDygDpgDxΔDyΔDμΔ=rank⁡f(x,y,μ)+rank⁡g(x,y,μ),then there exists a smooth curve of equilibrium in Rn+m+1, which passes through (0,0,μ0) and is transversal to the singular surface at (0,0,μ0). When μ increases through μ0, one eigenvalue of this system moves from C- (the open complex left half plane) to C+ (the open complex right half plane) if M/N>0 (resp., from C+ to C- if M/N<0) along the real axis by diverging through the infinity. The other (n-1) eigenvalues remain bounded and stay away from the origin. The constants M and N can be computed: (15)M=-trace⁡Dyfadj⁡DygDxg,Dyg0,0,μ0,N=DμΔ-DxΔDyΔDxfDyfDxgDyg-1DμfDμg0,0,μ0.For simplicity, let(16)fXt,Et,μ=f1f2=x1-x-xy-ex-dy2+bxy,gXt,Et,μ=epx-c-m,where X(t)=(x(t),y(t))T, E(t)=e(t), μ=m. Let D be a differential operator: Δ(X(t),E(t),μ)=DEg(X(t),E(t),μ)=px(t)-c. Applying Lemma 2 to system (13), we can obtain the following result.Theorem 3.
If (H1) holds, then system (13) has a singularity induced bifurcation at the interior equilibrium P∗(x∗,y∗,e∗) when the bifurcation parameter m increases through zero. That is, a stability switch occurs as parameter m passes through zero.Proof.
It is obvious thatDEgP∗(x∗,y∗,e∗)=0 and m=0, which implies that DEg has a simple zero eigenvalue. Here P∗(x∗,y∗,e∗)≐(p/c,b/dp/c,1-b/d-p/c).
By simple computation, we have that(17)trace⁡DEfadj⁡DEgDXg,DEgP∗=-px∗e∗≠0. Furthermore, it can be calculated that (18)DXfDEfDXgDEgP∗=1-2x-y2x-e-x-xby2x-2dy+bx0pe0px-cP∗=-pdx∗y∗e∗≠0,which follows that DXfDEfDXgDEg is nonsingular at P∗.
By simple computing,(19)DXfDEfDμfDXgDEgDμgDXΔDEΔDμΔP∗=-x+y2x-x-x0by2x-dy00pe0px-c0p000P∗=pdx∗y∗≠0;it is easy to show that rank⁡(f(X,E,μ))=2 and rank⁡(g(X,E,μ))=1. Hence, rankDXfDEfDμfDXgDEgDμgDXΔDEΔDμΔ=rank⁡(f(X,E,μ))+rank⁡(g(X,E,μ))+1.
Thus, all conditions of Lemma2 are satisfied. System (13) has a singularity induced bifurcation at the positive equilibrium P∗(x∗,y∗,e∗) when the bifurcation parameter m=0.
Furthermore, the constantsM and N can be computed: (20)M=-trace⁡DEfadj⁡DEgDXg,DEgP∗x∗,y∗,e∗=px∗e∗≠0,N=DμΔ-DXΔ,DEΔDXfDEfDXgDEg-1DμfDμgP∗x∗,y∗,e∗=-1e∗≠0,which follows that (21)MN=-px∗e∗2<0. If parameter m increases through zero, then one eigenvalue of system (13) will move from C+ to C- along the real axis by diverging into the infinity. Thus, the stability of the positive equilibrium P∗(x∗,y∗,e∗) changes from unstable to stable.
The Jacobian of system (13) evaluated at P∗ takes the form: (22)J=∂f,g∂X,E=1-2x-y2x-e-x-xby2x-dy0pe0px-c. According to the leading matrix Ξ in system (13) and JP∗, here Ξ=100010000, the characteristic equation of system (13) at P∗ is det(λΞ-JP∗)=0. That is, the characteristic equation can be expressed: (23)det⁡λΞ-JP∗=px∗e∗λ+bx∗=0. It is obvious that the rest eigenvalue (denoted by λ2) has negative real part. Furthermore, it follows from Theorem 1.1 [34] that there is only one eigenvalue diverging to infinity, but the rest eigenvalue is continuous, nonzero, and cannot jump from on half open complex plane to another one as m increases through 0. Table 1 shows the change in the signs of the real parts of eigenvalues (λ1 and λ2) due to the variation of economic interest of harvest effort. According to Table 1, it can be concluded that system (13) is stable at P∗ as m<0 and unstable as m>0. Consequently, a stability switch occurs as m increases through 0.Table 1
Signs of real parts of eigenvalues of system (13) at P∗.
Parameter value
R
e
λ
1
R
e
λ
2
m
<
0
−
−
m
>
0
+
−Remark 4.
According to [35], system (13) along the positive equilibrium locus yields index one matrix pencil (Ξ,J) for m≠0 and index two matrix pencil (Ξ,J) for m=0. This shows that there exists one index jump at a bifurcation point m=0. It leads to the rapid expansion of the population from the point of view of biology. If this phenomenon lasts for a long time, the population will be beyond the carrying capacity of the environment and the prey-predator system will be out of balance, which is disastrous. With the purpose of economic interest of harvest effort at an ideal level as well as maintaining the sustainable development of the biological resource, some related measures should be taken to eliminate the impulse phenomenon and stabilize system (13) when the economic interest is positive. Thus, a state feedback controller is designed to stabilize system (13) at the positive equilibrium P∗.By using Theorem 3-1.2 [36], a state feedback controller u(t)=k(e(t)-e∗), where k is a feedback gain and e∗ is the component of the positive equilibrium P∗, can be applied to stabilize system (13) at P∗. Thus, a controlled system is as follows:(24)dxdt=x1-x-xy-ex,dydt=-dy2+bxy,0=epx-c+ke-e∗-m,where all variables and parameters have the same interpretations as that in system (13). The feedback gain can be given in the following result.Theorem 5.
Suppose the feedback gaink satisfies one of the following cases:(i)
ifbp-2cd>0, then k<2p(cd-bc-pd)/bp-2cd;(ii)
ifbp-2cd<0, then k>2p(cd-bc-pd)/bp-2cd;the controlled system (24) is stable at the positive equilibrium P∗.Proof.
The Jacobian of system (24) evaluated at the interior equilibrium P∗ has the form (25)J~P∗=1-2x∗-b2d-e∗-x∗-x∗b22d-bx∗0pe∗0k. According to the leading matrix Ξ in system (13) and J~P∗, the characteristic equation of the controlled system (24) at P∗ is (26)λ2+bx∗+px∗e∗k-1+2x∗+b2d+e∗λ+bx∗px∗e∗k-1+2x∗+b2d+e∗=0. By using the Routh-Hurwitz criteria [37], the sufficient and necessary condition for the stability of the controlled system (24) at P∗ is that the feedback gain k satisfies the following: if bp-2cd>0, then k<2p(cd-bc-pd)/bp-2cd; otherwise, k>2p(cd-bc-pd)/bp-2cd. This ends the proof.Remark 6.
After introducing the feedback controller into system (13), the controlled system (24) can be stabilized at the positive equilibrium. The elimination of the singularity induced bifurcation means the prey-predator system restores ecological balance. Both sustainable development of the prey-predator system and the ideal economic interest of harvesting can be obtained by enhancing or reducing the harvest effort on the prey.


## 2.3. Saddle-Node Bifurcation

Ifm is taken as bifurcation parameter, the following result is true.Theorem 7.
If the assumption (H1) holds, then system (13) undergoes saddle-node bifurcation at m=msn.Proof.
The characteristic polynomial of system (13) at the positive equilibrium P∗(x∗,y∗,e∗) is given: (27)Fλ=λ2+l1λ+l0,where l1=bx∗+2x∗+b/2d-1-mc/(px∗-c)2 and l0=bx∗(2x∗+b/d-1-mc/(px∗-c)2). Clearly, it has a simple zero eigenvalue with left eigenvector usn=(-b,1) and right eigenvector vsn=(x∗,b/2d)T. And (28)usnDμfRP∗=usnDμf-DEfDEg-1DμgRP∗=bx∗px∗-c,usnDX2fRvsn,vsnP∗=usn∑i=1i=2eivsnTDXDXfiTvsnP∗=b22d+2bx∗+by∗4x∗b-1+2bmpcx∗px∗-c3,where ei(i=1,2) are unit vectors.
According to [31], system (13) undergoes saddle-node bifurcation when m=msn holds. This completes the proof.


## 2.4. Hopf Bifurcation and Control

From the discussion above, we know that if (H1) holds, there is only one positive equilibrium P∗(x∗,y∗,e∗); if m<msn, there are two positive equilibriums P∗+(x∗+,y∗+,e∗+) and P∗-(x∗-,y∗-,e∗-). In this subsection, we only investigate the dynamical behavior of system (6) at the positive equilibrium P∗+(x∗+,y∗+,e∗+) but neglect similar results on the other positive equilibria.For simplicity, let(29)fW=f1x,y,ef2x,y,e=x1-x-xy-ex-dy2+bxt-τy,gW=epx-c-m,where W=(x,y,e)T. In order to analyze the local stability of the positive equilibrium of system (6), we first use the linear transformation WT(t)=QNT(t), where (30)Nt=ut,vt,Et,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 (6), we can obtain(31)dutdt=ut1-u-utvt-Etut+pe∗+u2tpx∗+-c,dvtdt=-dv2t+but-τvt,0=Et-pe∗+utpx∗+-cput-c-m. Further, in order to derive the formula determining the properties of the positive equilibrium of system (31), we consider the local parametric Ψ of the third equation of system (31) as the literatures [38, 39], which is given: (32)NTt=ΨZt=N0T+u0Zt+v0hZt,gΨZt=0,where(33)u0=100100,v0=001,Zt=y1ty2t,N0=u∗v∗E∗=x∗+y∗+px∗+e∗+px∗+-c+e∗+,hZt=00h3y1t,y2t:R2⟶Ris a smooth mapping. The corresponding characteristic equation of the linearized system of the above system is(34)λ2+bu∗-Aλ-Abu∗+b22du∗e-λτ=0,where A=1-b/2d-u∗+mc/(pu∗-c)2.In order to investigate the distribution of roots of the transcendental equation (34), we introduce the following result proved by Ruan and Wei [10].Lemma 8.
Consider the exponential polynomial:(35)Pλ,e-λτ1,…,e-λτm=λn+p10λn-1+⋯+pn-10λ+pn0+p11λn-1+⋯+pn-11λ+pn1e-λτ1+⋯+p1mλn-1+⋯+pn-1mλ+pnme-λτm,where τi≥0(i=1,2,…,m) and pji (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 (34). When τ=0, (34) is rewritten as follows:(36)λ2+bu∗-Aλ+bu∗b2d-A=0.Thus, if b/2d-A<0, (36) has at least one positive root, which implies that system (6) without any time delay is unstable. If the condition(H2): 1+mc/(pu∗-c)2<min⁡{b/d+2u∗,b/2d+2u∗+bu∗}holds, two roots of (36) have always negative real parts. Thus, we have the following result.Lemma 9.
Assume that (H1) and (H2) hold, then the two roots of (36) have always negative real parts; that is, system (6) with τ=0 is locally asymptotically stable.Letωi(ω>0) be the root of (34), and we have the following equations by substituting it into (34) and separating the real and imaginary parts:(37)b2u∗2dcos⁡ωτ=ω2+Abu∗,b2u∗2dsin⁡ωτ=bu∗-Aω,which gives(38)ω4+a1ω2+a0=0,where a1=A2+b2u∗>0 and a0=b2u∗(A2-b2/4d2).Letω2=θ; (38) takes the form(39)Φθ=θ2+a1θ+a0=0.It is obvious that (39) has no real roots when Δ=a12-4a0<0. When a0>0, (39) has two negative roots; when Δ>0 and a0<0, (39) has one positive root. Thus, we have the following results.Lemma 10.
For (38), the following results are true:(i)
whena0>0, (38) does not have any positive real root;(ii)
whenΔ>0 and a0<0, (38) has only one positive root, defined by ω+=-a1+Δ/2.If condition (ii) of Lemma10 holds, and we denote (40)τj+=1ω+arccos⁡2dAbu∗+ω+2b2u∗+2jπ,j=0,1,2,…,then ±iω+ are a pair of purely imaginary roots of (34) with τ=τj+. Next, we will check whether the following transversality conditions are satisfied.Lemma 11.
Suppose that (H1) holds, the transversality condition (41)dRe⁡λτdττ=τj+>0is satisfied, and (dλ/dτ)-1 and Φ′(θ) have the same sign.Proof.
Differentiating the two sides of (34) with respect to τ, we get (42)dλdτ-1=2λ+bu∗-Aλb2u∗/2de-λτ-τλ.
For simplify, we defineω0 as ω and τj+ as τ, and we can obtain (43)sign⁡dRe⁡λdττ=τj+=sign⁡Re⁡dλdτ-1τ=τj+=sign⁡2ω2+bu∗-A2+2Abu∗bu∗-A2ω2+ω2+Abu∗2=1bu∗-A2ω2+ω2+Abu∗2sign⁡Φ′θ. This completes the proof.Summarizing the above lemmas and applying Lemma11 to (6), we have the following theorem.Theorem 12.
Suppose that (H1) holds, the following results are true.(1)
Whena0>0, system (6) is locally asymptotically stable for any τ≥0.(2)
WhenΔ>0 and a0<0, system (6) is locally asymptotically stable for τ∈[0,τ0) and unstable when τ>τ0 and system (6) undergoes a Hopf bifurcation at the equilibrium P∗+ when τ=τj(j=0,1,2,…).It has been realized that bifurcation control can be useful when bifurcation oscillation is harmful in some engineering applications, which has great potential in many technological disciplines. The aim of bifurcation control is to delay the onset of an inherent bifurcation, change the parameter value of an existing bifurcation point, stabilize a bifurcated solution or branch, and so forth. In order to control the Hopf bifurcation to achieve some desirable behaviors, a new hybrid control strategy was proposed in [40]. Recently, some scholars succeed in applying this method in various fields [40–43]. However, as far as we known, there is little attention on Hopf bifurcation of the singular biological economic system by use of this method. Extending this method on singular system is first time.By choosing the hybrid control strategy, the controlled model with delay is as follows:(44)dxdt=x1-x-xy-ex,dydt=α-dy2+bxt-τy+1-αy-y∗+,0=epx-c-m,where α(0<α<1) is the hybrid control parameter. The controlled system (44) reduces to the original system (6) if α=1. By selecting appropriate control parameter α, the Hopf bifurcation can be delayed (or advanced) or even eliminated completely without changing the equilibrium of the original system.The controlled system (44) has the same equilibrium as system (6). Let w1(t)=x(t)-x∗+,w2(t)=y(t)-y∗+, and w3(t)=e(t)-e∗+. The corresponding characteristic equation of the linear system of the controlled system (44) has the form (45)px∗+-cλ2+b11λ+b10+b00e-λτ=0,where b11=pe∗+x∗+/px∗+-c-1+2x∗++b/2d+e∗++bx∗+α-1+α,  b10=(-bx∗+α+1-α)(1-2x∗+-b/2d-e∗+-pe∗+x∗+/px∗+-c), and  b00=-αb2x∗+/2d<0. Therefore, in order to study the local stability of the equilibrium of system (44), we will study the following equation:(46)λ2+b11λ+b10+b00e-λτ=0.Whenτ=0, (46) is given:(47)λ2+b11λ+b10+b00=0.It is obvious that the positive equilibrium P∗ is locally asymptotically stable when τ=0 if the following condition(H3) b11>0 and b10+b00>0holds.Letωi(ω>0) be a root of (46). We have that(48)b00cos⁡ωτ=ω2+b10,b00sin⁡ωτ=-b11ω,which yields to(49)ω4+b112-2b10ω2+b102-b002=0.By simple computation,b112-2b10=(pe∗+x∗+/px∗+-c-1+2x∗++b/2d+e∗+)2>0 is obtained. If b102-b002>0, then (49) does not have any positive real root. If b102-b002≤0, then (49) has only one positive real root ω~+=1/2[2b10-b112+(b112-2b10)2-4(b102-b002)]. From (48), we can get (50)τ~j+=1ω~+arccos⁡ω~+2-b10b00+2jπ,j=0,1,2,…,and then ±iω~+ are a pair of purely imaginary roots of (46) with τ=τ~j+. Define τ~0=min⁡{τ~j+,j=0,1,2,…}. Further, if λ(τ)=α(τ)+iω(τ) is denoted by the root of (46), then it satisfies α(τ~0)=0, ω(τ~0)=ω~0.Next, we will check whether the following transversality condition is satisfied. Differentiating the two sides of (46) with respect to τ and applying the implicit theorem, we get (51)dλdτ-1=2λ+b11λb00e-λτ-τλ.For simplify, we defineω~+ as ω and τ~j+ as τ, and we can obtain (52)sign⁡dRe⁡λdττ=τ~j+=sign⁡Re⁡dλdτ-1τ=τ~j+=sign⁡Re⁡2λ+b11λb00e-λτλ=ω~+i=sign⁡2ω2+b112-2b10b112ω2+ω2+b102.Since b112-2b10>0, then the transversality condition (53)dRe⁡λτdττ=τ~j+>0is satisfied.Summarizing the above results, we have the following theorem.Theorem 13.
If (H3) holds, then we have the following results on controlled system (44).(1)
Ifb11>0 and b10+b00>0 or b102-b002>0, then the positive equilibrium of controlled system (44) is locally asymptotically stable for any τ≥0.(2)
Ifb102-b002≤0, the positive equilibrium of controlled system (44) is locally asymptotically stable for τ∈[0,τ~0) and unstable when τ>τ~0, and controlled system (44) undergoes a Hopf bifurcation at the equilibrium P∗ when τ=τ~0.Remark 14.
For the controlled system (44), it has been proven in Theorem 13 that we can delay the onset of Hopf bifurcation without changing the original equilibrium by choosing appropriate parameter.


## 3. The Model within the Fluctuating Environment

In this section, we will investigate the dynamics of system (9) in a fluctuation environment. Here, the method can be found in [25]. A linearized system of system (9) at the positive equilibrium P∗+ is given:(54)x˙t=a11xt+a12yt+a13et+ξ1t,y˙t=a21xt-τ+a22yt+ξ2t,0=a31xt+a33et+ξ3t,where (55)a11=1-b2d-2x∗+-mpx∗+-c,a12=-x∗+,a13=-x∗+,a21=b22d,a22=-bx∗+,a31=mppx∗+-c,a33=px∗+-c.The following system can be given by taking the Fourier transform of equations in (54):(56)ξ~s=BsX~s,where ξ~(s)=(ξ~1(s),ξ~2(s),ξ~3(s)), X~(s)=(x~(s),y~(s),e~(s)), and (57)Bs=is-a11-a12-a13-a21e-isτis-a220-a310-a33≐B11B12B13B21B22B23B31B32B33. Thus, a new equation is given by(58)X~s=B-1sξ~s,where (59)A-1s=k11k12k13k21k22k23k31k32k33≐1det⁡BsB11∗B21∗B31∗B12∗B22∗B32∗B13∗B23∗B33∗,and the symbol Bij∗ denotes the algebraic cofactor of Bij, i,j=1,2,3.If the functionΥ(t) has zero mean value, then the fluctuation intensity of the components is the frequency band s and s+ds is SΥ(s)ds, where the spectral density SΥ(s) is formally defined by (60)SΥsds=limT→∞⁡ξ~s2T. Thus, (61)SΥsds=limT→∞⁡ξ~s2T=limT→∞⁡1T∬T/2T/2ξtξt′¯eist-t′dtdt′.Therefore,(62)Sxs=∑j=13k1j2Sξjs,Sys=∑j=13k2j2Sξjs,Ses=∑j=13k3j2Sξjs,where the equation 〈ξi〉=0 and 〈ξi(t)ξj(t1)〉=δijδ(t-t1), i,j=1,2,3, are used to obtain (62). Furthermore, from Sξj(s)=1, j=1,2,3, the fluctuation intensities of x(t),y(t), and e(t) are written as follows (63)σx2=12π∫-∞+∞Sxsds=12π∑j=13∫-∞+∞k1j2Sξjsds=12π∑j=13∫-∞+∞k1j2ds,σy2=12π∫-∞+∞Sysds=12π∑j=13∫-∞+∞k2j2Sξjsds=12π∑j=13∫-∞+∞k2j2ds,σe2=12π∫-∞+∞Sesds=12π∑j=13∫-∞+∞k3j2Sξjsds=12π∑j=13∫-∞+∞k3j2ds.By simple computation, the fluctuation intensities ofx(t),y(t), and e(t) are given by (64)σx2=12π∑j=13∫-∞+∞Bj1∗det⁡Bs2ds,σy2=12π∑j=13∫-∞+∞Bj2∗det⁡Bs2ds,σe2=12π∑j=13∫-∞+∞Bj3∗det⁡Bs2ds.It should be pointed out that whenτ>0 the explicit values of the spectral densities of the populations and harvest efforts are difficult to obtain as the evaluation of the integrals is not easy at all. However, at a given value of time delay τ, the fluctuation intensities can be determined by numerical integration.


## 4. Numerical Simulations

With the help ofMATLAB, some numerical results of system (6) are provided to substantiate the analytic results in this section.The values of parameters area=1, b=0.6, c=1, d=1, p=50, and m=1. It is easy to check that the condition (H1): p-c-bp/d=19>0 is satisfied, and two roots x∗+=0.33689 and x∗-=0.08311 can be obtained by solving (11). So, two equilibriums P∗+(0.33689,0.34825,0.06311) and P∗-(0.08311,0.17298,0) of system (6) are obtained, respectively. Since there is zero component of P∗-, then P∗- contradicts the biological interpretations of interior equilibrium. Hence, the positive equilibrium P∗+ is considered in the following part, which is denoted by P∗ for implicity.


### 4.1. Singularity Induced Bifurcation and State Feedback Control

By selecting harvest effortm as the parameter, singularity induced bifurcation at the interior equilibrium P∗ is obtained. By using Theorem 3, it can be shown that system (13) has a singularity induced bifurcation at the interior equilibrium P∗, and a stability switch occurs as m increases through 0. Further, a state feedback controller u(t) can be applied to stabilize system (13) at P∗. According to Theorem 5, if the feedback gain k satisfies k>-17.24, then the controlled system (24) is stable at P∗ and the singularity induced bifurcation of system (13) is also eliminated. The dynamical responses of system (13) and the controlled system (24) can be shown in Figures 1 and 2.Figure 1
The dynamical responses of system (13) without feedback control.
(a)
(b)Figure 2
The dynamical responses of harvest effort of system (24) with feedback control.


### 4.2. Hopf Bifurcation and Hybrid Control

Since the assumption (H2) is true, the positive equilibrium P∗ of system (6) without any time delay is locally asymptotically stable (see Figure 3).Figure 3
The positive equilibriumP∗ of system (6) with τ=0 is locally asymptotically stable.Based on the given parameter values, the critical value of time delayτ0=2.955 can be obtained by solving the corresponding expression. From Theorem 12, the corresponding waveforms are shown in Figures 4 and 5. That is, when τ=2.950<τ0=2.955, the positive equilibrium P∗ of the system is locally asymptotically stable (see Figure 4), and Hopf bifurcations occur once τ=2.960>τ0=2.955 (see Figure 5).Figure 4
The positive equilibriumP∗ of system (6) with τ=2.950<τ0=2.955 is locally asymptotically stable.Figure 5
Hopf bifurcation occurs at the positive equilibrium whenτ=2.960>τ0=2.955.Now we choose appropriate value of parameterα to control singular biological economic system (6). By choosing α=0.85, which has the same equilibrium point as that of the original system, we find that the periodic solutions of the controlled system (44) are eliminated (Figure 6). Of course, the onset of the Hopf bifurcation may be advanced or delayed if choosing appropriate parameter values.Figure 6
The simulations show that period of Hopf bifurcation changes with the controlled parameterα.


### 4.3. System with the Fluctuating Environment

For numerical simulation of the stochastic version corresponding to the delayed model, the same set of data is used to simulate for the stochastic model by increasing the magnitude ofτ step by step. Figure 7 exhibits the period orbits of the stochastic model for τ=1.5. The fluctuation intensity for x(t) is σx2=0.43012. Figure 8 depicts the oscillation orbits in the case of the singular stochastic model with time delay τ=2.960. The fluctuation intensities for x(t) are σx2=0.3158. Compared with the oscillation observed in a deterministic environment, the amplitudes of oscillation for x(t),y(t), and E(t) are enhanced. By comparing the frequencies of the oscillations between Figures 4-5 and Figures 7-8, it is clear to see that the environmental fluctuation plays a crucial role in determining the magnitude of oscillation (as the magnitude of delay parameters is the same in both cases). Figure 8 has shown that the intensity of fluctuation for the population and harvest effort increase gradually from their steady state values as the delay parameter τ increases.Figure 7
Oscillation solutions of stochastic system (9) in the case of τ=1.5.Figure 8
Oscillation solutions of stochastic system (9) in the case of τ=2.960.


## 4.1. Singularity Induced Bifurcation and State Feedback Control

By selecting harvest effortm as the parameter, singularity induced bifurcation at the interior equilibrium P∗ is obtained. By using Theorem 3, it can be shown that system (13) has a singularity induced bifurcation at the interior equilibrium P∗, and a stability switch occurs as m increases through 0. Further, a state feedback controller u(t) can be applied to stabilize system (13) at P∗. According to Theorem 5, if the feedback gain k satisfies k>-17.24, then the controlled system (24) is stable at P∗ and the singularity induced bifurcation of system (13) is also eliminated. The dynamical responses of system (13) and the controlled system (24) can be shown in Figures 1 and 2.Figure 1
The dynamical responses of system (13) without feedback control.
(a)
(b)Figure 2
The dynamical responses of harvest effort of system (24) with feedback control.


## 4.2. Hopf Bifurcation and Hybrid Control

Since the assumption (H2) is true, the positive equilibrium P∗ of system (6) without any time delay is locally asymptotically stable (see Figure 3).Figure 3
The positive equilibriumP∗ of system (6) with τ=0 is locally asymptotically stable.Based on the given parameter values, the critical value of time delayτ0=2.955 can be obtained by solving the corresponding expression. From Theorem 12, the corresponding waveforms are shown in Figures 4 and 5. That is, when τ=2.950<τ0=2.955, the positive equilibrium P∗ of the system is locally asymptotically stable (see Figure 4), and Hopf bifurcations occur once τ=2.960>τ0=2.955 (see Figure 5).Figure 4
The positive equilibriumP∗ of system (6) with τ=2.950<τ0=2.955 is locally asymptotically stable.Figure 5
Hopf bifurcation occurs at the positive equilibrium whenτ=2.960>τ0=2.955.Now we choose appropriate value of parameterα to control singular biological economic system (6). By choosing α=0.85, which has the same equilibrium point as that of the original system, we find that the periodic solutions of the controlled system (44) are eliminated (Figure 6). Of course, the onset of the Hopf bifurcation may be advanced or delayed if choosing appropriate parameter values.Figure 6
The simulations show that period of Hopf bifurcation changes with the controlled parameterα.


## 4.3. System with the Fluctuating Environment

For numerical simulation of the stochastic version corresponding to the delayed model, the same set of data is used to simulate for the stochastic model by increasing the magnitude ofτ step by step. Figure 7 exhibits the period orbits of the stochastic model for τ=1.5. The fluctuation intensity for x(t) is σx2=0.43012. Figure 8 depicts the oscillation orbits in the case of the singular stochastic model with time delay τ=2.960. The fluctuation intensities for x(t) are σx2=0.3158. Compared with the oscillation observed in a deterministic environment, the amplitudes of oscillation for x(t),y(t), and E(t) are enhanced. By comparing the frequencies of the oscillations between Figures 4-5 and Figures 7-8, it is clear to see that the environmental fluctuation plays a crucial role in determining the magnitude of oscillation (as the magnitude of delay parameters is the same in both cases). Figure 8 has shown that the intensity of fluctuation for the population and harvest effort increase gradually from their steady state values as the delay parameter τ increases.Figure 7
Oscillation solutions of stochastic system (9) in the case of τ=1.5.Figure 8
Oscillation solutions of stochastic system (9) in the case of τ=2.960.


## 5. Conclusion and Discussion

After Ajraldi et al. [6] proposed a predator-prey system (1) with square root functional response, Xu and Yuan [9] introduced a time delay and studied the effect of time delay on system (2). In addition, Zhang et al. [30] studied the dynamics of a predator-prey system with stochastic fluctuation. Based on the above work, we propose a singular delayed biological economic system with herd behavior and within the deterministic environment or fluctuating environment, which yields system (4) and system (8).For the deterministic model, we obtain the conditions for the existence of the positive equilibrium. Singularity induced bifurcation which only occurs in the singular system is also obtained. Thus, a state feedback controller is designed to eliminate the bifurcation phenomenon. The condition of saddle-node bifurcation is also given. In order to study time delay on the effect of system (6), we introduce the new normal form of differential-algebraic systems due to the work of literatures [38, 39] and analyze the local stability of the positive equilibrium and the existence of Hopf bifurcation by taking time delay τ as bifurcation parameter. From simulations, we found that the critical value of time delay τ is less than that in the literature [9]. 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. Since bifurcation oscillation is harmful in some field, the hybrid control strategy is introduced into system (6). Thus, the Hopf bifurcation may advance, delay, and even eliminate by selecting the proper value of parameter. In fact, Hopf bifurcation of the controlled system (44) with hybrid control strategy can be eliminated when á = 0.85. In addition, for the stochastic model, we find that the populations and harvest effort oscillate when the system is effected by the stochastic fluctuation (see Figures 7 and 8).However, for controlled system (24), the explicit formulas determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions by using the center manifold theorem are not investigated. In addition, noise is internal to system (8) since both the birth and the death are randomly varying due to the outer environment, season change, climate, and weather. For simplicity, the function defined the Gaussian white noise t as the derivative of the Wiener process W(t). We leave this work in the future.



---

*Source: 302494-2015-10-27.xml*