# The Study of a Predator-Prey Model with Fear Effect Based on State-Dependent Harvesting Strategy

Authors: Y. Tian; H. M. Li
Journal: Complexity
 (2022)

---

## Abstract

In presence of predator population, the prey population may significantly change their behavior. Fear for predator population enhances the survival probability of prey population, and it can greatly reduce the reproduction of prey population. In this study, we propose a predator-prey fishery model introducing the cost of fear into prey reproduction with Holling type-II functional response and prey-dependent harvesting and investigate the global dynamics of the proposed model. For the system without harvest, it is shown that the level of fear may alter the stability of the positive equilibrium, and an expression of fear critical level is characterized. For the harvest system, the existence of the semitrivial order-1 periodic solution and positive order-q (q≥1) periodic solution is discussed by the construction of a Poincaré map on the phase set, and the threshold conditions are given, which can not only transform state-dependent harvesting into a cycle one but also provide a possibility to determine the harvest frequency. In addition, to ensure a certain robustness of the adopted harvest policy, the threshold condition for the stability of the order-q periodic solution is given. Meanwhile, to achieve a good economic profit, an optimization problem is formulated and the optimum harvest level is obtained. Mathematical findings have been validated in numerical simulation by MATLAB. Different effects of different harvest levels and different fear levels have been demonstrated by depicting figures in numerical simulation using MATLAB.


---

## Body

## 1. Introduction

Prey-predator interaction is a crucial topic in theoretical ecology and evolutionary biology. The history of the study about the prey-predator interactions dates back long. The pioneering work to describe the prey-predator interactions in mathematics belongs to the Lotka-Volterra model [1, 2]. Subsequently the model was improved by adding logistic growth term for the prey and variety of population-dependent response functions [3–15]. A prototype model that captures the prey-predator interaction takes the form(1)dxdt=bx−dx−cx2−ypx,y,dydt=eϕpx,yy−my,where xt and yt represent the densities of prey and predator population, respectively, b, d, (d<b) and c represent the birth rate, natural death rate, and density-dependent decay rate due to the intraspecies competition, respectively, px,y represents the functional response, e is the efficiency of conversion, m is natural mortality of predator, and ϕ is a monotonically increasing function.Due to prey-predator interactions, predators always have an impact, direct, indirect, or both, on prey population. In model (1), the termypx,y models the direct impact of predator on prey by catching and killing behavior. Meanwhile fear of predation risk can be regarded as the indirect impact of predator on prey, and some theoretical ecologists and biologists have realised that a prey-predator model should involve not only direct killing but also the fear [16, 17]. The fieldwork of Zanette et al. [18] on song sparrows observed the impact of fear and found a reduction in reproduction by 40% in the number of offspring due to the fear of predation. Based on this phenomenon, Wang et al. [19] incorporated a predator-dependent fear factor into the birth rate of prey in model (1) (i.e., replace b by by=b/1+ky) with linear and Holling type-II functional response to explore the effect of fear on population dynamics. The results show that high level of fear could stabilize the system. Das and Samanta [20] investigated the impact of fear in exponential form on a stochastic prey-predator system when the predator is provided additional food. Sahoo and Samanta [21] investigated a two prey-one predator model by including the cost of fear into prey reproduction and switching mechanism in predation. Das et al. [22] developed and explored a predator-prey model incorporating the cost of perceived fear into the birth and death rates of prey species with Holling type-II functional response. Sarkar and Khajanchi [23] and Kumar and Kumari [24] incorporated a form of fear factor into the birth rate of prey by assuming a nonzero minimum cost of fear. The impact of fear has also been investigated on prey-predator systems with prey refuge [25–27], Allee effect [26], hunting cooperation [28], and additional food resource for predator [20, 29].The study of resource management including fisheries, forestry, and wildlife management has great importance. It is necessary to harvest the population but harvesting should be regulated in such a way that ecological sustainability as well as conservation of the species can be implemented in a long run. Besides, it is always hoped that the sustained ability can be achieved at a high level of productivity and good economic profit. In the past decade, scholars considered different kinds of harvest on the dynamics of the predator-prey system such as continuous harvesting [30–33] and intermittent harvesting [34–37]. Compared to fixed time harvest strategy, the state-dependent harvesting strategy takes the existing resources of species into full consideration and can maintain the sustainability of species in certain level. State-dependent harvested system can be described by the impulsive semidynamical system [38–45]. Recently, Lai et al. [46] proposed and studied a Lotka-Volterra predator-prey system incorporating both continuous harvesting and fear effect; that is,(2)dxdt=bx1+ky−dx−cx2−pxy−qExa1E+a2x,dydt=epxy−my,where k is the level of fear, E is the fishing effort used to harvest, q is the catchability coefficient, and a1 and a2 are constants. The harvest in (2) is continuous and in Michaelis-Menten type. However, in reality, the harvest of species should consider the aspect of ecological sustainability as well as conservation. Thus, in most cases, species are caught intermittently, not continuously.To the best of our knowledge, to this day, still no scholars investigated the dynamic behavior of the predator-prey system incorporating both fear effect and intermittent harvesting, which motivated us to study a predator-prey model incorporating fear effect based on state-dependent harvesting strategy. The aim of this study is to check the influence of fear level on the stability of the positive steady state of the system without harvest. Meanwhile, for the harvest system, it mainly discusses the existence of the order-q (q>0) periodic solution, since it provides a possibility to transform the state-dependent harvesting into a cycle one. Meanwhile, in order to make a maximum economic profit in the harvest process, the optimal control problem is discussed. The organization of this study is as follows. In the next section, we introduce the mathematical model for predator-prey system with fear effect based on state-dependent harvesting strategy. In the same section, we present some preliminaries used in the discussion of the system dynamics. Section 3 is dedicated to the existence and stability of semitrivial order-1 and positive order-1 periodic solution. We also study the existence of order-2 and order-3 periodic solution. In Section 4, we demonstrate different effects of different harvest levels and different fear levels by depicting figures in numerical simulation using MATLAB. The paper concludes in Section 5, in which we briefly summarize the biological indications of our analytical findings.


## 2. Model Formulation and Preliminaries

### 2.1. Model Formulation

In presence of predator population, the prey population may significantly change their behavior. Fear for predator population enhances the survival probability of prey population, and it can greatly reduce the reproduction of prey population [23]. In this study, we consider a predator-prey model introducing the cost of fear into prey reproduction with Holling type-II functional response and a saturation function ϕ in equation (1); that is,(3)dxdt=bx1+ky−dx−cx2−pxy1+h1x,dydt=epxy1+h1+h2px−my.Where the variables, model parameters, and their units/dimensions are given in Table 1. To achieve the commercial purpose of the fishery, it is necessary to harvest the population in such a way that ecological sustainability as well as conservation of the species can be implemented in a long run. The harvest can be continuous or intermittent. In this work, a state-dependent harvest strategy is considered. Let l be the harvest level of prey population; that is, when the density of prey population reaches level l, the harvest is implemented, resulting in a portion of prey and predator being caught. Let E denote the harvest effort, which is dependent on the harvest level l, and let q1 and q2 be the catchability coefficients of prey and predator populations. In addition, to avoid the extinction of predator, it is necessary to release a quantity of predator pups, denoted by τ, which is also dependent on level l. Based on this consideration, the model with state-dependent harvesting takes the following form:(4)dxdt=bx1+ky−dx−cx2−pxy1+h1xdydt=epxy1+h1+ph2x−my,x≠l,Δx=−q1ExΔy=−q2Ey+τ,x=l.Table 1
The description of the model parameters and variables and their units/dimensions.
SymbolDescriptionUnits/dimensionsxDensity of prey populationMassyDensity of predator populationMassbBirth rate of prey population not affected by predators1/timedNatural death rate of prey population1/timecDecay rate due to intraspecies competition1/mass.1/timepRate of predation1/mass.1/timeh1Handling time1/masseConversion rate of prey biomass to predator biomassDimensionlessh2Conversion time1/massmDeath rate of predator1/timekLevel of fear1/masslHarvesting level of prey populationMassEHarvesting effortDimensionlessq1Catchability coefficient of preyDimensionlessq2Catchability coefficient of predatorDimensionlessτQuantity of predator pups releasedMassDenoteK≜b−d/c. Then K is the carrying capacity of prey population in absence of predator. System (4) is considered in the domain S=x,y|0≤x≤K,y≥0 for ecological practices. The purpose of this paper is to analyze the dynamics of system (4). Besides, it is always hoped that the harvest can be achieved at a good economic profit, and this requires determining an optimal harvest level l. Next, some preliminaries are listed for the analysis of the harvest model (4).


### 2.2. Preliminaries

Let us consider a general planar system:(5)dxdt=Px,y,dydt=Qx,y,ifχx,y≠0,Δx=αx,y,Δy=βx,y,ifχx,y=0,where x,y∈Ω⊂ℝ2, and χx,y=0 describes the states at which the harvest is implemented; α and β describe the effects of the harvest strategy. Px,y and Qx,y are arbitrarily derivative with respect to x,y∈Ω; χ, α, and β are linearly dependent on x and y; that is, χx, χy, αx, αy, βx, and βy are constant.The dynamic system constituted by the solution mapping defined by system (5) is called an impulsive semicontinuous dynamic system, denoted asΩ,π;I,MIMP, where π=π1,π2: Ω×R⟶Ω, MIMP≜x,y|χx,y=0, and(6)I:MIMP⟶NPHA=IMIMP≜x′,y′|x′=x+αx,y,y′=y+βx,y,x,y∈MIMP.Letzt=xt,yt be the solution of system (5) with initial value z0=z0. Denote γz,z0=zt|t≥t0withzt0=z0, also denoted as γz in short. Denote zk=ztk+∈γz, where tk∈∏≜tk|k=1,2,… with ztk∈MIMP.Definition 1 (priodic solution [47–49]).
The solutionz=z^tof system (5) is said to be periodic if there exists positive integerm≥1such thatz^m=z^0. Denotek≜minm∈ℕ,z^m=z^0; then orbitγz^is said to be an order-kperiodic orbit of system (5).Definition 2 (orbitally stable [47–49]).
An orbitγz^is said to be orbitally stable if, for anyε>0, there is a neighborhoodVofz^so that, for allzinV, there is a reparameterization of time (a smooth, monotonic function)t^tsuch thatzt−z^t^t<εfor allt≥t0.Definition 3 (asymptotic orbital stability [47–49]).
γz^is said to be asymptotically orbitally stable if it is orbitally stable and additionallyVmay be chosen so that, for allz∈V, there exists a constantτzsuch thatzt−z^t−τz⟶0ast⟶∞.Definition 4 (Poincaré map).
LetNω≜S∈NPHA|∃TS>0such thatπS,TS∈MIMP. Define the Poincaré mapϕN:Nω⟶NPHAas follows:ϕNS=S+≜IπS,TS=πS,TS+απS,TS,βπS,TS.Remark 1.
If there exists a pointL∈Nωandq>0such thatϕNqL=LandϕNjL≠L(j<q), that is,Lis a fixed point ofϕNq, then system (4) admits an order-qperiodic solution.Lemma 1 (analogue of Poincaré criterion [47–49]).
The order-qT-periodic solutionzt=ξt,ηtof system (5) is orbitally asymptotically stable and enjoying the property of the asymptotic phase if the multiplierμ2satisfies the conditionμ2q<1, where(7)μ2q=∏j=1qΔjexp∫0T∂P∂xξt,ηt+∂Q∂yξt,ηtdt,with(8)Δj=P+∂β/∂y∂ϕ/∂x−∂β/∂x∂χ/∂y+∂χ/∂xP∂χ/∂x+Q∂χ/∂y+Q+∂α/∂x∂χ/∂y−∂α/∂y∂ϕ/∂x+∂χ/∂yP∂χ/∂x+Q∂χ/∂y,and P+=Pξτj+,ητj+,Q+=Qξτj+,ητj+, andP,Q,∂α/∂x,∂α/∂y,∂β/∂x,∂β/∂y,∂χ/∂x,∂χ/∂yare calculated at the pointξτj,ητj.


## 2.1. Model Formulation

In presence of predator population, the prey population may significantly change their behavior. Fear for predator population enhances the survival probability of prey population, and it can greatly reduce the reproduction of prey population [23]. In this study, we consider a predator-prey model introducing the cost of fear into prey reproduction with Holling type-II functional response and a saturation function ϕ in equation (1); that is,(3)dxdt=bx1+ky−dx−cx2−pxy1+h1x,dydt=epxy1+h1+h2px−my.Where the variables, model parameters, and their units/dimensions are given in Table 1. To achieve the commercial purpose of the fishery, it is necessary to harvest the population in such a way that ecological sustainability as well as conservation of the species can be implemented in a long run. The harvest can be continuous or intermittent. In this work, a state-dependent harvest strategy is considered. Let l be the harvest level of prey population; that is, when the density of prey population reaches level l, the harvest is implemented, resulting in a portion of prey and predator being caught. Let E denote the harvest effort, which is dependent on the harvest level l, and let q1 and q2 be the catchability coefficients of prey and predator populations. In addition, to avoid the extinction of predator, it is necessary to release a quantity of predator pups, denoted by τ, which is also dependent on level l. Based on this consideration, the model with state-dependent harvesting takes the following form:(4)dxdt=bx1+ky−dx−cx2−pxy1+h1xdydt=epxy1+h1+ph2x−my,x≠l,Δx=−q1ExΔy=−q2Ey+τ,x=l.Table 1
The description of the model parameters and variables and their units/dimensions.
SymbolDescriptionUnits/dimensionsxDensity of prey populationMassyDensity of predator populationMassbBirth rate of prey population not affected by predators1/timedNatural death rate of prey population1/timecDecay rate due to intraspecies competition1/mass.1/timepRate of predation1/mass.1/timeh1Handling time1/masseConversion rate of prey biomass to predator biomassDimensionlessh2Conversion time1/massmDeath rate of predator1/timekLevel of fear1/masslHarvesting level of prey populationMassEHarvesting effortDimensionlessq1Catchability coefficient of preyDimensionlessq2Catchability coefficient of predatorDimensionlessτQuantity of predator pups releasedMassDenoteK≜b−d/c. Then K is the carrying capacity of prey population in absence of predator. System (4) is considered in the domain S=x,y|0≤x≤K,y≥0 for ecological practices. The purpose of this paper is to analyze the dynamics of system (4). Besides, it is always hoped that the harvest can be achieved at a good economic profit, and this requires determining an optimal harvest level l. Next, some preliminaries are listed for the analysis of the harvest model (4).


## 2.2. Preliminaries

Let us consider a general planar system:(5)dxdt=Px,y,dydt=Qx,y,ifχx,y≠0,Δx=αx,y,Δy=βx,y,ifχx,y=0,where x,y∈Ω⊂ℝ2, and χx,y=0 describes the states at which the harvest is implemented; α and β describe the effects of the harvest strategy. Px,y and Qx,y are arbitrarily derivative with respect to x,y∈Ω; χ, α, and β are linearly dependent on x and y; that is, χx, χy, αx, αy, βx, and βy are constant.The dynamic system constituted by the solution mapping defined by system (5) is called an impulsive semicontinuous dynamic system, denoted asΩ,π;I,MIMP, where π=π1,π2: Ω×R⟶Ω, MIMP≜x,y|χx,y=0, and(6)I:MIMP⟶NPHA=IMIMP≜x′,y′|x′=x+αx,y,y′=y+βx,y,x,y∈MIMP.Letzt=xt,yt be the solution of system (5) with initial value z0=z0. Denote γz,z0=zt|t≥t0withzt0=z0, also denoted as γz in short. Denote zk=ztk+∈γz, where tk∈∏≜tk|k=1,2,… with ztk∈MIMP.Definition 1 (priodic solution [47–49]).
The solutionz=z^tof system (5) is said to be periodic if there exists positive integerm≥1such thatz^m=z^0. Denotek≜minm∈ℕ,z^m=z^0; then orbitγz^is said to be an order-kperiodic orbit of system (5).Definition 2 (orbitally stable [47–49]).
An orbitγz^is said to be orbitally stable if, for anyε>0, there is a neighborhoodVofz^so that, for allzinV, there is a reparameterization of time (a smooth, monotonic function)t^tsuch thatzt−z^t^t<εfor allt≥t0.Definition 3 (asymptotic orbital stability [47–49]).
γz^is said to be asymptotically orbitally stable if it is orbitally stable and additionallyVmay be chosen so that, for allz∈V, there exists a constantτzsuch thatzt−z^t−τz⟶0ast⟶∞.Definition 4 (Poincaré map).
LetNω≜S∈NPHA|∃TS>0such thatπS,TS∈MIMP. Define the Poincaré mapϕN:Nω⟶NPHAas follows:ϕNS=S+≜IπS,TS=πS,TS+απS,TS,βπS,TS.Remark 1.
If there exists a pointL∈Nωandq>0such thatϕNqL=LandϕNjL≠L(j<q), that is,Lis a fixed point ofϕNq, then system (4) admits an order-qperiodic solution.Lemma 1 (analogue of Poincaré criterion [47–49]).
The order-qT-periodic solutionzt=ξt,ηtof system (5) is orbitally asymptotically stable and enjoying the property of the asymptotic phase if the multiplierμ2satisfies the conditionμ2q<1, where(7)μ2q=∏j=1qΔjexp∫0T∂P∂xξt,ηt+∂Q∂yξt,ηtdt,with(8)Δj=P+∂β/∂y∂ϕ/∂x−∂β/∂x∂χ/∂y+∂χ/∂xP∂χ/∂x+Q∂χ/∂y+Q+∂α/∂x∂χ/∂y−∂α/∂y∂ϕ/∂x+∂χ/∂yP∂χ/∂x+Q∂χ/∂y,and P+=Pξτj+,ητj+,Q+=Qξτj+,ητj+, andP,Q,∂α/∂x,∂α/∂y,∂β/∂x,∂β/∂y,∂χ/∂x,∂χ/∂yare calculated at the pointξτj,ητj.


## 3. Main Works

Define(9)eM≜m+h1+h2pmKpK.TheneM is called the critical value of the conversion; that is, when e≤eM, the conversion is not enough to maintain the survival of predator and the predator population will go to extinction. Thus, in this work, it is reasonable to assume that e>eM.Besides, denote(10)x∗≜mep−mh1−pmh2,y∗≜p+kd+cx∗1+h1x∗2−4pkd+cx∗−b1+h1x∗−p−kd+cx∗1+h1x∗2pk,k∗≜−ph1h1d+cx∗−b+c1+h1x∗c1+h1x∗2c1+h1x∗+h1d+cx∗,T∗=x∗ph1y∗1+h1x∗2−c,D∗=bk1+ky∗2+p/1+h1x∗epx∗y∗1+h1x∗+h2px∗2>0.For system (3), the following result holds.Theorem 1.
There are three equilibria for system (3) whene>eM: two boundary saddlesO0,0andEK,0and one positive equilibriumE∗x∗,y∗. Moreover, one of the two following cases holds:(i)
E∗x∗,y∗is a stable focus or node in case ofk>k∗(ii)
E∗x∗,y∗is unstable in case ofk<k∗, and a unique stable limit cycle exists, denoted byΓLCProof.
The Jacobian matrixJE¯ of model (3) at the equilibrium E¯x¯,y¯ is(11)b1+ky¯−d−2cx¯−py¯1+h1x¯2−kbx¯1+ky¯2−px¯1+h1x¯epy¯1+h1x¯+h2px¯2epx¯1+h1x¯+h2px¯−m.
It is obvious thatO0,0 and EK,0 are saddles. At the equilibrium E∗x∗,y∗, the characteristic equation λ2−T∗λ+D∗=0. If y∗<yM≜a1+h1x∗2/ph1 holds, then T∗<0. By equation (6), y∗<yM if and only if k>k∗. Thus, the positive equilibrium E∗x∗,y∗ is locally asymptotically stable in case of k>k∗ and unstable in case of k<k∗. In this case, there exists a unique stable limit cycle for system (3).


### 3.1. Semitrivial Order-1 Periodic Solution forτ=0

Whenτ=0, there is yt≡0 for t>0 with y0=0. In this case, system (3) is reduced to the following system:(12)dxdt=b−dx1−xK,x≠l,Δx=−q1Ex,x=l.Settingx0=1−q1El, the solution of equation dx/dt=b−dx1−x/K with x0=x0 is(13)xt=K1−q1Elexpb−dtK−1−q1El+1−q1Elexpb−dt.Let(14)T≜1rlnK−1−q1El1−q1EK−l.Then there isxT=l and xT+=1−q1El=x0 by impulse effect. Thus, the following result holds.Theorem 2.
System (4) withτ=0has a semitrivial order-1 periodic solution forn−1T<t≤nT:(15)ξ¯t=K1−q1Elexpb−dt−n−1TK−1−q1El+1−q1Elexpb−dt−n−1T,η¯t=0,which is orbitally asymptotically stable when R0<1, where(16)R0≜1−q2E1+h1+h2plK−1−q1El1+1−q1Eh1+h2plK−lep/c1+Kh1+h2p1−q1EK−lK−1−q1Elm/b−d.Proof.
To discuss the stability ofξ¯t,η¯t, let us consider a small disturbance δ0. The trajectory starting from B11−q1El,δ0 is denoted by ξ^t,η^t. This disturbed trajectory first intersects the harvest set MIMP at point B2l,y1 when t=T+δt, and then it jumps to point B2+1−q1El,δ1. Thus, there is(17)ξ^0=1−q1El,η^0=δ0,ξ^T+δt=l,η^T+δt=y1.
Letδx=ξ^t−ξ¯t and δy=η^t−η¯t. Then δx0=ξ^0−ξ¯0=0, and δy0=η^0−η¯0=δ0. Setting δy1≜y1, for 0<t<T, the variables δx and δy can be expressed by the relation(18)δxtδyt=Φt0δ0+oδ02,where Φt=ϕij2×2 is the fundamental solution satisfying the variation equation.(19)Φ′t=b−d−2cξ¯t−kbξ¯t−pξ¯t1+h1ξ¯t0epξ¯t1+h1+ph2ξ¯t−mΦt,Φ0=I2.
According to the first-order Taylor expansion onη^t, there is δy1=η^T+δt≈δyT=ϕ22Tδ0, where(20)ϕ22T=1+h1+h2plK−1−q1El1+1−q1Eh1+h2plK−lep/c1+Kh1+h2p1−q1EK−lK−1−q1Elm/b−d.
By impulse effect, there isδ1=1−q2Eδy1=1−q2Eϕ22Tδ0. Thus, if inequality (10) holds, there is δ1<δ0. By the arbitrary of δ0, it concludes that the order-1 semitrivial periodic solution is orbitally asymptotically stable.Corollary 1.
The semitrivial order-1 periodic solutionξ¯t,η¯tis orbitally asymptotically stable if one of the two following cases is satisfied: (i)l≤x∗and (ii)l>x∗andE>E∗≜max0,ϕ22T−1/ϕ22Tq2.


### 3.2. Positive Order-K Periodic Solution forτ>0

Since the harvest may cause the extinction of predator whenτ=0, in order to keep the predator species from going extinct, it is necessary to reduce the harvest strength and release a certain quantity of predator pups.For0≤x≤K, define(21)yL˜x≜p+kd+cx1+h1x2+4pkb−d−cx1+h1x−p+kd+cx1+h1x2pk.LetN0 and M0 denote the intersection point between y=yL˜x and the phase set NPHA and the harvest set MIMP, respectively; G0 denotes the intersection point between y=τ and the phase set NPHA; in general τ≤τmax≜yN0=yL˜1−q1El. For a point S on NPHA with 0≤ys≤yN0, if the trajectory of system (4) starting from S1−q1El,ys intersects the harvest set MIMP, then it defines a function relationship between y and x for 1−q1El≤x≤l, denoted by y=yx,ys, which satisfies(22)dydx=epxy/1+h1+ph2x−mybx/1+ky−dx−cx2−pxy/1+h1x≐κx,y,y1−q1El,ys=ys.By equation (22), the function y=yx,ys can be expressed as follows:(23)yx,ys=ys+∫1−q1Elxκu,yu,ysdu.Definel¯≜maxl|yl,yN0≤yL˜l. When l≤l¯, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP, and denote the intersection point by N1; that is, N1=πN0,T0 for some T0>0. Define τf≜yL˜1−q1El−1−q2Eπ2N0,T0.


#### 3.2.1. Existence of Order-1 Periodic Solution

By Theorem1, the dynamic behavior of system (3) varies with the model parameter k. Thus, the discussions will be divided according to parameter k and harvest level l.Case I:k>k∗Case I-1:l≤x∗Sincex≡1−q1El on NPHA, map ϕN in Definition 4 is only a function of y. Next, the Poincaré map ϕN will be characterized and its main property will be analyzed.For0≤ys≤yN0=yL˜1−q1El, by Definition 4, there is ϕNys=1−q1Eyl,ys+τ. Meanwhile, for ys>yN0, there exists a unique ys′∈0,yN0 and Ts>0 such that ys′=π21−q1El,ys,Ts. Then ϕNys=1−q2Eyl,π21−q1El,ys,Ts+τ. To sum up, there is(24)ϕNys=1−q2Eyl,ys+τ,0≤ys≤yN0,1−q2Eyl,π21−q1El,ys,Ts+τ,ys>yN0.Property 1.
For system (4), whenl≤x∗, the Poincaré map ϕN defined by equation (24) has the following properties:(i)
ϕNis continuous on0,+∞. Moreover,ϕNis increasing on0,yL˜1−q1Eland decreasing onyL˜1−q1El,+∞(ii)
ϕNis continuously differentiable on0,+∞, andϕNis concave on0,yL˜1−q1El(iii)
There exists a horizontal asymptoteϕN=τ; that is,ϕNys⟶τwhenys⟶+∞
DefinefNy≜ϕNy−y. The following result holds.Theorem 3.
There exists a unique positive order-1 periodic solution for system (4) when0<l≤x∗and0<τ≤τmax.Proof.
By Remark1, the existence of order-1 periodic solution is equivalent to the existence of a point L∈NPHA such that yL is a fixed point of ϕN. By Property 1 (i), fN is continuous on 0,+∞. Since fN0=ϕN0=τ>0 and fNy=ϕNy−y⟶−∞ as y⟶+∞, by the intermediary property of continuous function, there exists at least one yL>0 such that fNyL=0; that is, ϕNyL=yL. Thus, the trajectory of system (4) starting from L1−q1El,yL forms an order-1 periodic orbit.
Next, the location and uniqueness of the order-1 periodic orbit will be analyzed. By Property1 (i), ϕN achieves its maximum at y=yN0=yL˜1−q1El. It is obvious that fNτ>0.
Ifτ=τf, then ϕNyN0=yN0.
Ifτ<τf, then fNyN0<0, which means that yL∈τ,yN0, as shown in Figure 1(a). Since ϕNy is concave on 0,yN0, yL is unique.
Ifτ>τf, there is fNyN0>0; that is, ϕNyN0>yN0. Since ϕN is decreasing on yN0,+∞, ϕNϕNyN0<ϕNyN0; that is, fNϕNyN0<0. Besides, define yτ1≜maxy|ϕNy=yN0 and yτ2≜miny|ϕNy=yN0. Denote ymin≜minyτ1,ϕNy. Then there is yL∈yL˜1−q1El,ymin and yL is unique, as shown in Figure 1(b)).Figure 1
Illustration of the Poincaré mapϕN for different values of τ. The parameters are taken as b=0.7,d=0.2,c=0.005,p=0.1,h1=0.036,e=0.44,h2=1.44,m=0.2,q1=0.8,q2=0.6, and l=25%K. (a) τ=0.4286; (b) τ=2.5714.
(a)(b)Case I-2:l>x∗: in case of x∗<l≤l¯, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP. When l>l¯, the trajectory starting from point N0 does not intersect the harvest set MIMP. Denote yM2=maxy|ϕNy=yM0 and yM1=miny|ϕNy=yM0. Then the domain of ϕN is 0,yM1∪yM2,+∞. Define τM1≜yM1−1−q2EyL˜l and τM2≜yM2−1−q2EyL˜l.Theorem 4.
There exists a positive order-1 periodic solution for system (4) when (i)l≤l¯andτ∈0,τmaxor (ii)l>l¯andτ∈0,τM1∪τM2,τmax.Proof.
Whenx∗<l≤l¯, similar to the proof of Theorem 3, system (4) admits an order-1 periodic solution. For l>l¯, if yM1>1−q2EyL˜l, then, for τ∈0,τM1, there is ϕNyM1≤yM1. Combining with ϕNyG0>yG0, it can be concluded that system (4) admits an order-1 periodic solution. For τ∈τM2,τmax, there is ϕNyM2≥yM2 and ϕNyϕNyM2<ϕNyM2; thus, there exists yL∈yM2,ϕNyM2 such that ϕNyL=yL; that is, system (4) admits an order-1 periodic solution. □Case II:k<k∗: in this case, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP.Theorem 5.
There exists a positive order-1 periodic solution for system (4) when (i)l≤l¯andτ∈0,τmaxor (ii)l>l¯andτ∈0,τM1. Moreover, the order-1 periodic solution is unique whenl≤x∗.


#### 3.2.2. Stability of the Order-1 Periodic Solution

Letξt,ηt0≤t≤T be an order-1 periodic solution of system (4). Denote ξ0=ξ0, ξ1=ξT, η0=η0, and η1=ηT. Define(25)Θη1≜ln1−q2Eη1+τ1−q2Eη1b/1+kη1−d−cl−pη1/1+h1lb/1+k1−q2Eη1+τ−d−c1−q1El−p1−q2Eη1+τ/1+h11−q1El.Theorem 6.
The order-1T-periodic solutionξt,ηtis orbitally asymptotically stable if(26)∫0+Tph1ξtηt1+h1ξt2−cξtdt<Θη1.Proof.
In system (4), there is(27)Px,y=bx1+ky−dx−cx2−pxy1+h1x,Qx,y=epxy1+h1+ph2x−my,χx,y=x−l,αx,y=−q1Ex,βx,y=−q2Ey+τ.
Then(28)∂P∂x=b1+ky−d−2cx−py1+h1x2,∂Q∂y=epx1+h1+ph2x−m,∂χ∂x=1,∂χ∂y=0,∂α∂x=−q1E,∂α∂y=0,∂β∂x=0,∂β∂y=−q2E.
With a direct calculation, there is(29)P+I1+βy′χx′−βx′χy′+Q+I1+αx′χy′−αy′χx′PIχx′+QIχy′=1−q2EP1−q1El,1−q2EyL−+τPl,yL−,and(30)∫0+T∂P∂x+∂Q∂yξt,ηtdt=∫0+Tξdtξ+ηdtη+∫0+Tph1ξη1+h1ξ2−cξdt=lnξTξ0+lnηTη0+∫0+Tph1ξη1+h1ξ2−cξdt.
Thus, in case of (26), there is μ21<1; then, by Lemma 1, the order-1 T-periodic solution ξt,ηt is orbitally asymptotically stable.Theorem 7.
Forl≤x∗, ifτ≤τf, the order-1 periodic solution for system (4) is globally orbitally asymptotically stable.Proof.
By Theorem3, system (4) admits a unique order-1 periodic solution when l≤x∗. If τ≤τf, there exists a unique yL∈τ,yN0 such that ϕNyL=yL. Thus, for any y0, a sequence ykk=1,2,… is obtained under ϕN; that is, yk=ϕNyk−1. If y0<yL, then yk is a monotonically increasing sequence with yk<yL, so the limit is yL. Similarly, if y0∈yL,yN0, then yk is a monotonically bounded decreasing sequence, and the limit is yL. If y0>yN0, then y1=ϕNy0∈0,yN0; thus ykk=1,2,… is a monotonically bounded sequence with limit yL. To sum up, by the arbitrariness of y0, the order-1 periodic solution is globally attractive and so is globally orbitally asymptotically stable.


#### 3.2.3. Existence of Order-q (q≥2) Periodic Solution

Forl≤x∗, by Theorem 3, if τ≤τf, the order-1 periodic solution is orbitally asymptotically stable and globally attractive, which means that system (4) does not admit order-q (q≥2) periodic solution. For τ>τf, there exists unique yL2∈yN0,ϕNyN0 such that ϕNyL2=yL2. Let yL1∈0,yN0 such that ϕNyL1=yL2. Then ϕN2yL1=ϕNyL2=yL2. Meanwhile, let yN1∈0,yL1 and yN2∈yL2,+∞ such that ϕNyN1=ϕNyN2=yN0.Theorem 8.
Forl≤l¯andτ>τf, if (i)ϕN2yN0<yN0or (ii)ϕN2yN0≥yN0andμ21>1holds, system (4) admits an order-2 periodic solution.Proof.
Obviously,ϕN2yN1=ϕN2yN2=ϕNyN0. It can be easily checked that ϕN2y is increasing on 0,yN1 and yN0,yN2, and ϕN2y is decreasing on yN1,yN0 and yN2,+∞.(i)
ϕN2yN0<yN0. In this case, there is ϕNyN0>yN2; that is, ϕN2yN2>yN2; then ϕN2ϕNyN0<ϕNyN0. Besides, there is ϕN2yL1=yL2>yL1. Thus, there exist yP1∈yL1,yN0 and yP2∈yN2,ϕNyN0 such that ϕN2yP1=yP1 and ϕN2yP2=yP2. Moreover, there are ϕNyP1=yP2 and ϕNyP2=yP1.(ii)
ϕN2yN0≥yN0. In this case, there is ϕNyN0≤yN2; that is, ϕN2yN2≤yN2. For any y∈yN0,ϕNyN0, there is yN0≤ϕN2yN0≤ϕNy≤ϕNyN0. Next, it mainly discusses the property of ϕN on yN0,ϕNyN0. Let y0=yN0. Then y1=ϕNy0=ϕNyN0>yL2>yN0, y2=ϕNy1=ϕN2y0>y0, and y3=ϕNy2<ϕNy0=y1. Then, under ϕN, a sequence yk is obtained, where(31)y0<y2<y4<⋯<yL2<⋯<y5<y3<y1.
DenoteyP1=limk⟶∞y2k and yP2=limk⟶∞y2k+1. It is obvious that yP1≤yL2≤yP2. Since μ21>1, then yP1<yL2<yP2. Moreover, there is ϕNyP1=yP2 and ϕNyP2=yP1. It can be concluded that μ22<1, that is, the order-2 periodic solution is orbitally asymptotically stable and globally attractive.Theorem 9.
Forl≤l¯andτ>τf, system (4) admits an order-3 periodic solution if and only ifϕN2yN0≤yN1. Moreover, there is at least one order-3 periodic solution whenϕN2yN0=yN1, and there are at least two order-3 periodic solutions whenϕN2yN0<yN1.Proof.
“Necessity.” Proof by contradiction. Assume thatϕN2yN0>yN1.
SinceϕN2yN0<yL2, if ϕN2yN0≥yN0, system (4) admits a stable order-1 periodic solution or a stable order-2 periodic solution. Moreover, there is ϕNyN0≤yN2. If τ≤yN1, that is, ϕNτ≤yN0, there exist yR1∈0,yN1 and yR4∈yN2,+∞ such that ϕN2yR4=ϕN2yR4=yN0. It can be easily checked that ϕN3 is increasing on 0,yR1, yN1,yN0, and yN2,yR4, and ϕN3 is decreasing on yR1,yN1, yN0,yN2, and yR4,+∞. If τ>yN1, then R1 does not exist, ϕN3 is increasing on yN1,yN0 and yN2,yR4, and ϕN3 is decreasing on 0,yN1, yN0,yN2, and yR4,+∞. In any case, since ϕN3yN1=ϕN2yN0≥yN0>yN1, ϕN3yN0>ϕNyL2=yN0, and ϕN3yR4=ϕNyN0<yR4, it can be concluded that ϕN3y=y if and only if y=yL2; that is, the order-3 periodic solution does not exist.
IfyN1<ϕN2yN0<yN0, system (4) simultaneously admits an order-1 periodic solution and order-2 periodic solution. Moreover, there is ϕNyN0>yN2. If τ≤yN1, that is, ϕNτ≤yN0, there exist yR1∈0,yN1, yR2∈yL1,yN0, yR3∈yN0,yL2, and yR4∈yN2,+∞ such that ϕN2yR1=ϕN2yR2=ϕN2yR3=ϕN2yR4=yN0. It can be easily checked that ϕN3 is increasing on 0,yR1, yN1,yR2, yN0,yR3, and yN2,yR4, and ϕN3 is decreasing on yR1,yN1, yR2,yN0, yR3,yN2, and yR4,+∞. If τ>yN1, then R1 does not exist, ϕN3 is increasing on yN1,yR2, yN0,yR3, and yN2,yR4, and ϕN3 is decreasing on 0,yN1, yR2,yN0, yR3,yN2, and yR4,+∞. Since ϕN3yN1=ϕN2yN0≥yN0>yN1 and ϕN3yR4=ϕNyN0<yR4, it can be concluded that ϕN3y=y if and only if y=yL2; that is, the order-3 periodic solution does not exist.
“Sufficiency.” IfϕN2yN0≤yN1, then there is ϕN3yR1=ϕNyN0>0 and ϕN3yN1=ϕN2yN0≤yN1. If ϕN2yN0=yN1, then there exists an order-3 periodic solution since ϕNyN1≠yN1. If ϕN2yN0<yN1, then there exists at least one yf∈yR1,yN1 such that ϕN3yf=yf and ϕNyf>yf; that is, system (4) admits an order-3 periodic solution.
Next, the number of order-3 periodic solutions will be discussed:(i)
WhenϕN2yN0=yN1, there are ϕN3yN1=yN1, ϕN3yN0=yN0, and ϕN3ϕNyN0=ϕNyN0. Moreover, there are yN0=ϕNyN1 and yN1=ϕNϕNyN0; that is, system (4) admits at least one order-3 periodic solution.(ii)
WhenϕN2yN0<yN1, there are ϕN3yR1=ϕN3yR2=ϕN3yR3=ϕN3yR4=ϕNyN0>yR4>yR3>yR2>yR1, and ϕN3yN1<yN1, ϕN3yN0<yN0, ϕN3yN2<yN2, and ϕN3ϕNyN0=ϕN2τ<ϕNyN0; then there exist y1∈yR1,yN1, y2∈yN1,yR2, y3∈yR2,yN0, y4∈yN0,yR3, y5∈yN2,yR4, and y6∈yR4,ϕNyN0 such that ϕN3yi=yi (i=1,2,…,6). Moreover, there are ϕNy1=y3, ϕNy3=y6, ϕNy6=y1; ϕNy2=y4, ϕNy4=y5, and ϕNy5=y2; that is, system (4) admits at least two order-3 periodic solutions.


### 3.3. Optimal Harvest Level Determination

To achieve the commercial purpose of the fishery, it is necessary to harvest the population, and it is always hoped that the sustained ability can be achieved at a good economic profit. For the harvest problem, it is necessary to determine the controlled valuesE and τ and harvest level l, and this in general involves the optimization theory [48, 49].Letl be the harvest level, which is a decision variable. Theorems 3 and 4 show that system (4) admits an order-1 periodic solution when l≤l¯ or l>l¯ with τ≤yM1−1−q2EyL˜l. Since the harvest effort and yield of released predator are dependent on the harvest level, then it is assumed that El and τl take the following forms:(32)El=E1+E2−E1l−l1l2−l1,τl=τ1+τ2−τ1l−l1l2−l1,where l1≤l≤l2, l1 and l2 are minimum and maximum of the harvest level, and E1 (E2) and τ1 (τ2) are the harvest effort and yield of released predator at the harvest level l1 (l2).Letc1 and c2 be the unit selling prices of prey and predator, let c3 be the unit cost of harvest, and let c4 be the unit cost in breeding predator. Then the benefits from harvest can be described as Fbenefitl=c1q1Ell+c2q2ElηTl−c3E2l−c4τl. The objective is to maximize the unit benefits; that is,(33)maxFbenefitlTlsuch thatl1≤l≤l¯,l¯<l≤l¯,τl≤yM1−1−q2EyL˜l.


## 3.1. Semitrivial Order-1 Periodic Solution forτ=0

Whenτ=0, there is yt≡0 for t>0 with y0=0. In this case, system (3) is reduced to the following system:(12)dxdt=b−dx1−xK,x≠l,Δx=−q1Ex,x=l.Settingx0=1−q1El, the solution of equation dx/dt=b−dx1−x/K with x0=x0 is(13)xt=K1−q1Elexpb−dtK−1−q1El+1−q1Elexpb−dt.Let(14)T≜1rlnK−1−q1El1−q1EK−l.Then there isxT=l and xT+=1−q1El=x0 by impulse effect. Thus, the following result holds.Theorem 2.
System (4) withτ=0has a semitrivial order-1 periodic solution forn−1T<t≤nT:(15)ξ¯t=K1−q1Elexpb−dt−n−1TK−1−q1El+1−q1Elexpb−dt−n−1T,η¯t=0,which is orbitally asymptotically stable when R0<1, where(16)R0≜1−q2E1+h1+h2plK−1−q1El1+1−q1Eh1+h2plK−lep/c1+Kh1+h2p1−q1EK−lK−1−q1Elm/b−d.Proof.
To discuss the stability ofξ¯t,η¯t, let us consider a small disturbance δ0. The trajectory starting from B11−q1El,δ0 is denoted by ξ^t,η^t. This disturbed trajectory first intersects the harvest set MIMP at point B2l,y1 when t=T+δt, and then it jumps to point B2+1−q1El,δ1. Thus, there is(17)ξ^0=1−q1El,η^0=δ0,ξ^T+δt=l,η^T+δt=y1.
Letδx=ξ^t−ξ¯t and δy=η^t−η¯t. Then δx0=ξ^0−ξ¯0=0, and δy0=η^0−η¯0=δ0. Setting δy1≜y1, for 0<t<T, the variables δx and δy can be expressed by the relation(18)δxtδyt=Φt0δ0+oδ02,where Φt=ϕij2×2 is the fundamental solution satisfying the variation equation.(19)Φ′t=b−d−2cξ¯t−kbξ¯t−pξ¯t1+h1ξ¯t0epξ¯t1+h1+ph2ξ¯t−mΦt,Φ0=I2.
According to the first-order Taylor expansion onη^t, there is δy1=η^T+δt≈δyT=ϕ22Tδ0, where(20)ϕ22T=1+h1+h2plK−1−q1El1+1−q1Eh1+h2plK−lep/c1+Kh1+h2p1−q1EK−lK−1−q1Elm/b−d.
By impulse effect, there isδ1=1−q2Eδy1=1−q2Eϕ22Tδ0. Thus, if inequality (10) holds, there is δ1<δ0. By the arbitrary of δ0, it concludes that the order-1 semitrivial periodic solution is orbitally asymptotically stable.Corollary 1.
The semitrivial order-1 periodic solutionξ¯t,η¯tis orbitally asymptotically stable if one of the two following cases is satisfied: (i)l≤x∗and (ii)l>x∗andE>E∗≜max0,ϕ22T−1/ϕ22Tq2.


## 3.2. Positive Order-K Periodic Solution forτ>0

Since the harvest may cause the extinction of predator whenτ=0, in order to keep the predator species from going extinct, it is necessary to reduce the harvest strength and release a certain quantity of predator pups.For0≤x≤K, define(21)yL˜x≜p+kd+cx1+h1x2+4pkb−d−cx1+h1x−p+kd+cx1+h1x2pk.LetN0 and M0 denote the intersection point between y=yL˜x and the phase set NPHA and the harvest set MIMP, respectively; G0 denotes the intersection point between y=τ and the phase set NPHA; in general τ≤τmax≜yN0=yL˜1−q1El. For a point S on NPHA with 0≤ys≤yN0, if the trajectory of system (4) starting from S1−q1El,ys intersects the harvest set MIMP, then it defines a function relationship between y and x for 1−q1El≤x≤l, denoted by y=yx,ys, which satisfies(22)dydx=epxy/1+h1+ph2x−mybx/1+ky−dx−cx2−pxy/1+h1x≐κx,y,y1−q1El,ys=ys.By equation (22), the function y=yx,ys can be expressed as follows:(23)yx,ys=ys+∫1−q1Elxκu,yu,ysdu.Definel¯≜maxl|yl,yN0≤yL˜l. When l≤l¯, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP, and denote the intersection point by N1; that is, N1=πN0,T0 for some T0>0. Define τf≜yL˜1−q1El−1−q2Eπ2N0,T0.


### 3.2.1. Existence of Order-1 Periodic Solution

By Theorem1, the dynamic behavior of system (3) varies with the model parameter k. Thus, the discussions will be divided according to parameter k and harvest level l.Case I:k>k∗Case I-1:l≤x∗Sincex≡1−q1El on NPHA, map ϕN in Definition 4 is only a function of y. Next, the Poincaré map ϕN will be characterized and its main property will be analyzed.For0≤ys≤yN0=yL˜1−q1El, by Definition 4, there is ϕNys=1−q1Eyl,ys+τ. Meanwhile, for ys>yN0, there exists a unique ys′∈0,yN0 and Ts>0 such that ys′=π21−q1El,ys,Ts. Then ϕNys=1−q2Eyl,π21−q1El,ys,Ts+τ. To sum up, there is(24)ϕNys=1−q2Eyl,ys+τ,0≤ys≤yN0,1−q2Eyl,π21−q1El,ys,Ts+τ,ys>yN0.Property 1.
For system (4), whenl≤x∗, the Poincaré map ϕN defined by equation (24) has the following properties:(i)
ϕNis continuous on0,+∞. Moreover,ϕNis increasing on0,yL˜1−q1Eland decreasing onyL˜1−q1El,+∞(ii)
ϕNis continuously differentiable on0,+∞, andϕNis concave on0,yL˜1−q1El(iii)
There exists a horizontal asymptoteϕN=τ; that is,ϕNys⟶τwhenys⟶+∞
DefinefNy≜ϕNy−y. The following result holds.Theorem 3.
There exists a unique positive order-1 periodic solution for system (4) when0<l≤x∗and0<τ≤τmax.Proof.
By Remark1, the existence of order-1 periodic solution is equivalent to the existence of a point L∈NPHA such that yL is a fixed point of ϕN. By Property 1 (i), fN is continuous on 0,+∞. Since fN0=ϕN0=τ>0 and fNy=ϕNy−y⟶−∞ as y⟶+∞, by the intermediary property of continuous function, there exists at least one yL>0 such that fNyL=0; that is, ϕNyL=yL. Thus, the trajectory of system (4) starting from L1−q1El,yL forms an order-1 periodic orbit.
Next, the location and uniqueness of the order-1 periodic orbit will be analyzed. By Property1 (i), ϕN achieves its maximum at y=yN0=yL˜1−q1El. It is obvious that fNτ>0.
Ifτ=τf, then ϕNyN0=yN0.
Ifτ<τf, then fNyN0<0, which means that yL∈τ,yN0, as shown in Figure 1(a). Since ϕNy is concave on 0,yN0, yL is unique.
Ifτ>τf, there is fNyN0>0; that is, ϕNyN0>yN0. Since ϕN is decreasing on yN0,+∞, ϕNϕNyN0<ϕNyN0; that is, fNϕNyN0<0. Besides, define yτ1≜maxy|ϕNy=yN0 and yτ2≜miny|ϕNy=yN0. Denote ymin≜minyτ1,ϕNy. Then there is yL∈yL˜1−q1El,ymin and yL is unique, as shown in Figure 1(b)).Figure 1
Illustration of the Poincaré mapϕN for different values of τ. The parameters are taken as b=0.7,d=0.2,c=0.005,p=0.1,h1=0.036,e=0.44,h2=1.44,m=0.2,q1=0.8,q2=0.6, and l=25%K. (a) τ=0.4286; (b) τ=2.5714.
(a)(b)Case I-2:l>x∗: in case of x∗<l≤l¯, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP. When l>l¯, the trajectory starting from point N0 does not intersect the harvest set MIMP. Denote yM2=maxy|ϕNy=yM0 and yM1=miny|ϕNy=yM0. Then the domain of ϕN is 0,yM1∪yM2,+∞. Define τM1≜yM1−1−q2EyL˜l and τM2≜yM2−1−q2EyL˜l.Theorem 4.
There exists a positive order-1 periodic solution for system (4) when (i)l≤l¯andτ∈0,τmaxor (ii)l>l¯andτ∈0,τM1∪τM2,τmax.Proof.
Whenx∗<l≤l¯, similar to the proof of Theorem 3, system (4) admits an order-1 periodic solution. For l>l¯, if yM1>1−q2EyL˜l, then, for τ∈0,τM1, there is ϕNyM1≤yM1. Combining with ϕNyG0>yG0, it can be concluded that system (4) admits an order-1 periodic solution. For τ∈τM2,τmax, there is ϕNyM2≥yM2 and ϕNyϕNyM2<ϕNyM2; thus, there exists yL∈yM2,ϕNyM2 such that ϕNyL=yL; that is, system (4) admits an order-1 periodic solution. □Case II:k<k∗: in this case, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP.Theorem 5.
There exists a positive order-1 periodic solution for system (4) when (i)l≤l¯andτ∈0,τmaxor (ii)l>l¯andτ∈0,τM1. Moreover, the order-1 periodic solution is unique whenl≤x∗.


### 3.2.2. Stability of the Order-1 Periodic Solution

Letξt,ηt0≤t≤T be an order-1 periodic solution of system (4). Denote ξ0=ξ0, ξ1=ξT, η0=η0, and η1=ηT. Define(25)Θη1≜ln1−q2Eη1+τ1−q2Eη1b/1+kη1−d−cl−pη1/1+h1lb/1+k1−q2Eη1+τ−d−c1−q1El−p1−q2Eη1+τ/1+h11−q1El.Theorem 6.
The order-1T-periodic solutionξt,ηtis orbitally asymptotically stable if(26)∫0+Tph1ξtηt1+h1ξt2−cξtdt<Θη1.Proof.
In system (4), there is(27)Px,y=bx1+ky−dx−cx2−pxy1+h1x,Qx,y=epxy1+h1+ph2x−my,χx,y=x−l,αx,y=−q1Ex,βx,y=−q2Ey+τ.
Then(28)∂P∂x=b1+ky−d−2cx−py1+h1x2,∂Q∂y=epx1+h1+ph2x−m,∂χ∂x=1,∂χ∂y=0,∂α∂x=−q1E,∂α∂y=0,∂β∂x=0,∂β∂y=−q2E.
With a direct calculation, there is(29)P+I1+βy′χx′−βx′χy′+Q+I1+αx′χy′−αy′χx′PIχx′+QIχy′=1−q2EP1−q1El,1−q2EyL−+τPl,yL−,and(30)∫0+T∂P∂x+∂Q∂yξt,ηtdt=∫0+Tξdtξ+ηdtη+∫0+Tph1ξη1+h1ξ2−cξdt=lnξTξ0+lnηTη0+∫0+Tph1ξη1+h1ξ2−cξdt.
Thus, in case of (26), there is μ21<1; then, by Lemma 1, the order-1 T-periodic solution ξt,ηt is orbitally asymptotically stable.Theorem 7.
Forl≤x∗, ifτ≤τf, the order-1 periodic solution for system (4) is globally orbitally asymptotically stable.Proof.
By Theorem3, system (4) admits a unique order-1 periodic solution when l≤x∗. If τ≤τf, there exists a unique yL∈τ,yN0 such that ϕNyL=yL. Thus, for any y0, a sequence ykk=1,2,… is obtained under ϕN; that is, yk=ϕNyk−1. If y0<yL, then yk is a monotonically increasing sequence with yk<yL, so the limit is yL. Similarly, if y0∈yL,yN0, then yk is a monotonically bounded decreasing sequence, and the limit is yL. If y0>yN0, then y1=ϕNy0∈0,yN0; thus ykk=1,2,… is a monotonically bounded sequence with limit yL. To sum up, by the arbitrariness of y0, the order-1 periodic solution is globally attractive and so is globally orbitally asymptotically stable.


### 3.2.3. Existence of Order-q (q≥2) Periodic Solution

Forl≤x∗, by Theorem 3, if τ≤τf, the order-1 periodic solution is orbitally asymptotically stable and globally attractive, which means that system (4) does not admit order-q (q≥2) periodic solution. For τ>τf, there exists unique yL2∈yN0,ϕNyN0 such that ϕNyL2=yL2. Let yL1∈0,yN0 such that ϕNyL1=yL2. Then ϕN2yL1=ϕNyL2=yL2. Meanwhile, let yN1∈0,yL1 and yN2∈yL2,+∞ such that ϕNyN1=ϕNyN2=yN0.Theorem 8.
Forl≤l¯andτ>τf, if (i)ϕN2yN0<yN0or (ii)ϕN2yN0≥yN0andμ21>1holds, system (4) admits an order-2 periodic solution.Proof.
Obviously,ϕN2yN1=ϕN2yN2=ϕNyN0. It can be easily checked that ϕN2y is increasing on 0,yN1 and yN0,yN2, and ϕN2y is decreasing on yN1,yN0 and yN2,+∞.(i)
ϕN2yN0<yN0. In this case, there is ϕNyN0>yN2; that is, ϕN2yN2>yN2; then ϕN2ϕNyN0<ϕNyN0. Besides, there is ϕN2yL1=yL2>yL1. Thus, there exist yP1∈yL1,yN0 and yP2∈yN2,ϕNyN0 such that ϕN2yP1=yP1 and ϕN2yP2=yP2. Moreover, there are ϕNyP1=yP2 and ϕNyP2=yP1.(ii)
ϕN2yN0≥yN0. In this case, there is ϕNyN0≤yN2; that is, ϕN2yN2≤yN2. For any y∈yN0,ϕNyN0, there is yN0≤ϕN2yN0≤ϕNy≤ϕNyN0. Next, it mainly discusses the property of ϕN on yN0,ϕNyN0. Let y0=yN0. Then y1=ϕNy0=ϕNyN0>yL2>yN0, y2=ϕNy1=ϕN2y0>y0, and y3=ϕNy2<ϕNy0=y1. Then, under ϕN, a sequence yk is obtained, where(31)y0<y2<y4<⋯<yL2<⋯<y5<y3<y1.
DenoteyP1=limk⟶∞y2k and yP2=limk⟶∞y2k+1. It is obvious that yP1≤yL2≤yP2. Since μ21>1, then yP1<yL2<yP2. Moreover, there is ϕNyP1=yP2 and ϕNyP2=yP1. It can be concluded that μ22<1, that is, the order-2 periodic solution is orbitally asymptotically stable and globally attractive.Theorem 9.
Forl≤l¯andτ>τf, system (4) admits an order-3 periodic solution if and only ifϕN2yN0≤yN1. Moreover, there is at least one order-3 periodic solution whenϕN2yN0=yN1, and there are at least two order-3 periodic solutions whenϕN2yN0<yN1.Proof.
“Necessity.” Proof by contradiction. Assume thatϕN2yN0>yN1.
SinceϕN2yN0<yL2, if ϕN2yN0≥yN0, system (4) admits a stable order-1 periodic solution or a stable order-2 periodic solution. Moreover, there is ϕNyN0≤yN2. If τ≤yN1, that is, ϕNτ≤yN0, there exist yR1∈0,yN1 and yR4∈yN2,+∞ such that ϕN2yR4=ϕN2yR4=yN0. It can be easily checked that ϕN3 is increasing on 0,yR1, yN1,yN0, and yN2,yR4, and ϕN3 is decreasing on yR1,yN1, yN0,yN2, and yR4,+∞. If τ>yN1, then R1 does not exist, ϕN3 is increasing on yN1,yN0 and yN2,yR4, and ϕN3 is decreasing on 0,yN1, yN0,yN2, and yR4,+∞. In any case, since ϕN3yN1=ϕN2yN0≥yN0>yN1, ϕN3yN0>ϕNyL2=yN0, and ϕN3yR4=ϕNyN0<yR4, it can be concluded that ϕN3y=y if and only if y=yL2; that is, the order-3 periodic solution does not exist.
IfyN1<ϕN2yN0<yN0, system (4) simultaneously admits an order-1 periodic solution and order-2 periodic solution. Moreover, there is ϕNyN0>yN2. If τ≤yN1, that is, ϕNτ≤yN0, there exist yR1∈0,yN1, yR2∈yL1,yN0, yR3∈yN0,yL2, and yR4∈yN2,+∞ such that ϕN2yR1=ϕN2yR2=ϕN2yR3=ϕN2yR4=yN0. It can be easily checked that ϕN3 is increasing on 0,yR1, yN1,yR2, yN0,yR3, and yN2,yR4, and ϕN3 is decreasing on yR1,yN1, yR2,yN0, yR3,yN2, and yR4,+∞. If τ>yN1, then R1 does not exist, ϕN3 is increasing on yN1,yR2, yN0,yR3, and yN2,yR4, and ϕN3 is decreasing on 0,yN1, yR2,yN0, yR3,yN2, and yR4,+∞. Since ϕN3yN1=ϕN2yN0≥yN0>yN1 and ϕN3yR4=ϕNyN0<yR4, it can be concluded that ϕN3y=y if and only if y=yL2; that is, the order-3 periodic solution does not exist.
“Sufficiency.” IfϕN2yN0≤yN1, then there is ϕN3yR1=ϕNyN0>0 and ϕN3yN1=ϕN2yN0≤yN1. If ϕN2yN0=yN1, then there exists an order-3 periodic solution since ϕNyN1≠yN1. If ϕN2yN0<yN1, then there exists at least one yf∈yR1,yN1 such that ϕN3yf=yf and ϕNyf>yf; that is, system (4) admits an order-3 periodic solution.
Next, the number of order-3 periodic solutions will be discussed:(i)
WhenϕN2yN0=yN1, there are ϕN3yN1=yN1, ϕN3yN0=yN0, and ϕN3ϕNyN0=ϕNyN0. Moreover, there are yN0=ϕNyN1 and yN1=ϕNϕNyN0; that is, system (4) admits at least one order-3 periodic solution.(ii)
WhenϕN2yN0<yN1, there are ϕN3yR1=ϕN3yR2=ϕN3yR3=ϕN3yR4=ϕNyN0>yR4>yR3>yR2>yR1, and ϕN3yN1<yN1, ϕN3yN0<yN0, ϕN3yN2<yN2, and ϕN3ϕNyN0=ϕN2τ<ϕNyN0; then there exist y1∈yR1,yN1, y2∈yN1,yR2, y3∈yR2,yN0, y4∈yN0,yR3, y5∈yN2,yR4, and y6∈yR4,ϕNyN0 such that ϕN3yi=yi (i=1,2,…,6). Moreover, there are ϕNy1=y3, ϕNy3=y6, ϕNy6=y1; ϕNy2=y4, ϕNy4=y5, and ϕNy5=y2; that is, system (4) admits at least two order-3 periodic solutions.


## 3.2.1. Existence of Order-1 Periodic Solution

By Theorem1, the dynamic behavior of system (3) varies with the model parameter k. Thus, the discussions will be divided according to parameter k and harvest level l.Case I:k>k∗Case I-1:l≤x∗Sincex≡1−q1El on NPHA, map ϕN in Definition 4 is only a function of y. Next, the Poincaré map ϕN will be characterized and its main property will be analyzed.For0≤ys≤yN0=yL˜1−q1El, by Definition 4, there is ϕNys=1−q1Eyl,ys+τ. Meanwhile, for ys>yN0, there exists a unique ys′∈0,yN0 and Ts>0 such that ys′=π21−q1El,ys,Ts. Then ϕNys=1−q2Eyl,π21−q1El,ys,Ts+τ. To sum up, there is(24)ϕNys=1−q2Eyl,ys+τ,0≤ys≤yN0,1−q2Eyl,π21−q1El,ys,Ts+τ,ys>yN0.Property 1.
For system (4), whenl≤x∗, the Poincaré map ϕN defined by equation (24) has the following properties:(i)
ϕNis continuous on0,+∞. Moreover,ϕNis increasing on0,yL˜1−q1Eland decreasing onyL˜1−q1El,+∞(ii)
ϕNis continuously differentiable on0,+∞, andϕNis concave on0,yL˜1−q1El(iii)
There exists a horizontal asymptoteϕN=τ; that is,ϕNys⟶τwhenys⟶+∞
DefinefNy≜ϕNy−y. The following result holds.Theorem 3.
There exists a unique positive order-1 periodic solution for system (4) when0<l≤x∗and0<τ≤τmax.Proof.
By Remark1, the existence of order-1 periodic solution is equivalent to the existence of a point L∈NPHA such that yL is a fixed point of ϕN. By Property 1 (i), fN is continuous on 0,+∞. Since fN0=ϕN0=τ>0 and fNy=ϕNy−y⟶−∞ as y⟶+∞, by the intermediary property of continuous function, there exists at least one yL>0 such that fNyL=0; that is, ϕNyL=yL. Thus, the trajectory of system (4) starting from L1−q1El,yL forms an order-1 periodic orbit.
Next, the location and uniqueness of the order-1 periodic orbit will be analyzed. By Property1 (i), ϕN achieves its maximum at y=yN0=yL˜1−q1El. It is obvious that fNτ>0.
Ifτ=τf, then ϕNyN0=yN0.
Ifτ<τf, then fNyN0<0, which means that yL∈τ,yN0, as shown in Figure 1(a). Since ϕNy is concave on 0,yN0, yL is unique.
Ifτ>τf, there is fNyN0>0; that is, ϕNyN0>yN0. Since ϕN is decreasing on yN0,+∞, ϕNϕNyN0<ϕNyN0; that is, fNϕNyN0<0. Besides, define yτ1≜maxy|ϕNy=yN0 and yτ2≜miny|ϕNy=yN0. Denote ymin≜minyτ1,ϕNy. Then there is yL∈yL˜1−q1El,ymin and yL is unique, as shown in Figure 1(b)).Figure 1
Illustration of the Poincaré mapϕN for different values of τ. The parameters are taken as b=0.7,d=0.2,c=0.005,p=0.1,h1=0.036,e=0.44,h2=1.44,m=0.2,q1=0.8,q2=0.6, and l=25%K. (a) τ=0.4286; (b) τ=2.5714.
(a)(b)Case I-2:l>x∗: in case of x∗<l≤l¯, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP. When l>l¯, the trajectory starting from point N0 does not intersect the harvest set MIMP. Denote yM2=maxy|ϕNy=yM0 and yM1=miny|ϕNy=yM0. Then the domain of ϕN is 0,yM1∪yM2,+∞. Define τM1≜yM1−1−q2EyL˜l and τM2≜yM2−1−q2EyL˜l.Theorem 4.
There exists a positive order-1 periodic solution for system (4) when (i)l≤l¯andτ∈0,τmaxor (ii)l>l¯andτ∈0,τM1∪τM2,τmax.Proof.
Whenx∗<l≤l¯, similar to the proof of Theorem 3, system (4) admits an order-1 periodic solution. For l>l¯, if yM1>1−q2EyL˜l, then, for τ∈0,τM1, there is ϕNyM1≤yM1. Combining with ϕNyG0>yG0, it can be concluded that system (4) admits an order-1 periodic solution. For τ∈τM2,τmax, there is ϕNyM2≥yM2 and ϕNyϕNyM2<ϕNyM2; thus, there exists yL∈yM2,ϕNyM2 such that ϕNyL=yL; that is, system (4) admits an order-1 periodic solution. □Case II:k<k∗: in this case, the trajectory of system (4) starting from N0 will intersect the harvest set MIMP.Theorem 5.
There exists a positive order-1 periodic solution for system (4) when (i)l≤l¯andτ∈0,τmaxor (ii)l>l¯andτ∈0,τM1. Moreover, the order-1 periodic solution is unique whenl≤x∗.


## 3.2.2. Stability of the Order-1 Periodic Solution

Letξt,ηt0≤t≤T be an order-1 periodic solution of system (4). Denote ξ0=ξ0, ξ1=ξT, η0=η0, and η1=ηT. Define(25)Θη1≜ln1−q2Eη1+τ1−q2Eη1b/1+kη1−d−cl−pη1/1+h1lb/1+k1−q2Eη1+τ−d−c1−q1El−p1−q2Eη1+τ/1+h11−q1El.Theorem 6.
The order-1T-periodic solutionξt,ηtis orbitally asymptotically stable if(26)∫0+Tph1ξtηt1+h1ξt2−cξtdt<Θη1.Proof.
In system (4), there is(27)Px,y=bx1+ky−dx−cx2−pxy1+h1x,Qx,y=epxy1+h1+ph2x−my,χx,y=x−l,αx,y=−q1Ex,βx,y=−q2Ey+τ.
Then(28)∂P∂x=b1+ky−d−2cx−py1+h1x2,∂Q∂y=epx1+h1+ph2x−m,∂χ∂x=1,∂χ∂y=0,∂α∂x=−q1E,∂α∂y=0,∂β∂x=0,∂β∂y=−q2E.
With a direct calculation, there is(29)P+I1+βy′χx′−βx′χy′+Q+I1+αx′χy′−αy′χx′PIχx′+QIχy′=1−q2EP1−q1El,1−q2EyL−+τPl,yL−,and(30)∫0+T∂P∂x+∂Q∂yξt,ηtdt=∫0+Tξdtξ+ηdtη+∫0+Tph1ξη1+h1ξ2−cξdt=lnξTξ0+lnηTη0+∫0+Tph1ξη1+h1ξ2−cξdt.
Thus, in case of (26), there is μ21<1; then, by Lemma 1, the order-1 T-periodic solution ξt,ηt is orbitally asymptotically stable.Theorem 7.
Forl≤x∗, ifτ≤τf, the order-1 periodic solution for system (4) is globally orbitally asymptotically stable.Proof.
By Theorem3, system (4) admits a unique order-1 periodic solution when l≤x∗. If τ≤τf, there exists a unique yL∈τ,yN0 such that ϕNyL=yL. Thus, for any y0, a sequence ykk=1,2,… is obtained under ϕN; that is, yk=ϕNyk−1. If y0<yL, then yk is a monotonically increasing sequence with yk<yL, so the limit is yL. Similarly, if y0∈yL,yN0, then yk is a monotonically bounded decreasing sequence, and the limit is yL. If y0>yN0, then y1=ϕNy0∈0,yN0; thus ykk=1,2,… is a monotonically bounded sequence with limit yL. To sum up, by the arbitrariness of y0, the order-1 periodic solution is globally attractive and so is globally orbitally asymptotically stable.


## 3.2.3. Existence of Order-q (q≥2) Periodic Solution

Forl≤x∗, by Theorem 3, if τ≤τf, the order-1 periodic solution is orbitally asymptotically stable and globally attractive, which means that system (4) does not admit order-q (q≥2) periodic solution. For τ>τf, there exists unique yL2∈yN0,ϕNyN0 such that ϕNyL2=yL2. Let yL1∈0,yN0 such that ϕNyL1=yL2. Then ϕN2yL1=ϕNyL2=yL2. Meanwhile, let yN1∈0,yL1 and yN2∈yL2,+∞ such that ϕNyN1=ϕNyN2=yN0.Theorem 8.
Forl≤l¯andτ>τf, if (i)ϕN2yN0<yN0or (ii)ϕN2yN0≥yN0andμ21>1holds, system (4) admits an order-2 periodic solution.Proof.
Obviously,ϕN2yN1=ϕN2yN2=ϕNyN0. It can be easily checked that ϕN2y is increasing on 0,yN1 and yN0,yN2, and ϕN2y is decreasing on yN1,yN0 and yN2,+∞.(i)
ϕN2yN0<yN0. In this case, there is ϕNyN0>yN2; that is, ϕN2yN2>yN2; then ϕN2ϕNyN0<ϕNyN0. Besides, there is ϕN2yL1=yL2>yL1. Thus, there exist yP1∈yL1,yN0 and yP2∈yN2,ϕNyN0 such that ϕN2yP1=yP1 and ϕN2yP2=yP2. Moreover, there are ϕNyP1=yP2 and ϕNyP2=yP1.(ii)
ϕN2yN0≥yN0. In this case, there is ϕNyN0≤yN2; that is, ϕN2yN2≤yN2. For any y∈yN0,ϕNyN0, there is yN0≤ϕN2yN0≤ϕNy≤ϕNyN0. Next, it mainly discusses the property of ϕN on yN0,ϕNyN0. Let y0=yN0. Then y1=ϕNy0=ϕNyN0>yL2>yN0, y2=ϕNy1=ϕN2y0>y0, and y3=ϕNy2<ϕNy0=y1. Then, under ϕN, a sequence yk is obtained, where(31)y0<y2<y4<⋯<yL2<⋯<y5<y3<y1.
DenoteyP1=limk⟶∞y2k and yP2=limk⟶∞y2k+1. It is obvious that yP1≤yL2≤yP2. Since μ21>1, then yP1<yL2<yP2. Moreover, there is ϕNyP1=yP2 and ϕNyP2=yP1. It can be concluded that μ22<1, that is, the order-2 periodic solution is orbitally asymptotically stable and globally attractive.Theorem 9.
Forl≤l¯andτ>τf, system (4) admits an order-3 periodic solution if and only ifϕN2yN0≤yN1. Moreover, there is at least one order-3 periodic solution whenϕN2yN0=yN1, and there are at least two order-3 periodic solutions whenϕN2yN0<yN1.Proof.
“Necessity.” Proof by contradiction. Assume thatϕN2yN0>yN1.
SinceϕN2yN0<yL2, if ϕN2yN0≥yN0, system (4) admits a stable order-1 periodic solution or a stable order-2 periodic solution. Moreover, there is ϕNyN0≤yN2. If τ≤yN1, that is, ϕNτ≤yN0, there exist yR1∈0,yN1 and yR4∈yN2,+∞ such that ϕN2yR4=ϕN2yR4=yN0. It can be easily checked that ϕN3 is increasing on 0,yR1, yN1,yN0, and yN2,yR4, and ϕN3 is decreasing on yR1,yN1, yN0,yN2, and yR4,+∞. If τ>yN1, then R1 does not exist, ϕN3 is increasing on yN1,yN0 and yN2,yR4, and ϕN3 is decreasing on 0,yN1, yN0,yN2, and yR4,+∞. In any case, since ϕN3yN1=ϕN2yN0≥yN0>yN1, ϕN3yN0>ϕNyL2=yN0, and ϕN3yR4=ϕNyN0<yR4, it can be concluded that ϕN3y=y if and only if y=yL2; that is, the order-3 periodic solution does not exist.
IfyN1<ϕN2yN0<yN0, system (4) simultaneously admits an order-1 periodic solution and order-2 periodic solution. Moreover, there is ϕNyN0>yN2. If τ≤yN1, that is, ϕNτ≤yN0, there exist yR1∈0,yN1, yR2∈yL1,yN0, yR3∈yN0,yL2, and yR4∈yN2,+∞ such that ϕN2yR1=ϕN2yR2=ϕN2yR3=ϕN2yR4=yN0. It can be easily checked that ϕN3 is increasing on 0,yR1, yN1,yR2, yN0,yR3, and yN2,yR4, and ϕN3 is decreasing on yR1,yN1, yR2,yN0, yR3,yN2, and yR4,+∞. If τ>yN1, then R1 does not exist, ϕN3 is increasing on yN1,yR2, yN0,yR3, and yN2,yR4, and ϕN3 is decreasing on 0,yN1, yR2,yN0, yR3,yN2, and yR4,+∞. Since ϕN3yN1=ϕN2yN0≥yN0>yN1 and ϕN3yR4=ϕNyN0<yR4, it can be concluded that ϕN3y=y if and only if y=yL2; that is, the order-3 periodic solution does not exist.
“Sufficiency.” IfϕN2yN0≤yN1, then there is ϕN3yR1=ϕNyN0>0 and ϕN3yN1=ϕN2yN0≤yN1. If ϕN2yN0=yN1, then there exists an order-3 periodic solution since ϕNyN1≠yN1. If ϕN2yN0<yN1, then there exists at least one yf∈yR1,yN1 such that ϕN3yf=yf and ϕNyf>yf; that is, system (4) admits an order-3 periodic solution.
Next, the number of order-3 periodic solutions will be discussed:(i)
WhenϕN2yN0=yN1, there are ϕN3yN1=yN1, ϕN3yN0=yN0, and ϕN3ϕNyN0=ϕNyN0. Moreover, there are yN0=ϕNyN1 and yN1=ϕNϕNyN0; that is, system (4) admits at least one order-3 periodic solution.(ii)
WhenϕN2yN0<yN1, there are ϕN3yR1=ϕN3yR2=ϕN3yR3=ϕN3yR4=ϕNyN0>yR4>yR3>yR2>yR1, and ϕN3yN1<yN1, ϕN3yN0<yN0, ϕN3yN2<yN2, and ϕN3ϕNyN0=ϕN2τ<ϕNyN0; then there exist y1∈yR1,yN1, y2∈yN1,yR2, y3∈yR2,yN0, y4∈yN0,yR3, y5∈yN2,yR4, and y6∈yR4,ϕNyN0 such that ϕN3yi=yi (i=1,2,…,6). Moreover, there are ϕNy1=y3, ϕNy3=y6, ϕNy6=y1; ϕNy2=y4, ϕNy4=y5, and ϕNy5=y2; that is, system (4) admits at least two order-3 periodic solutions.


## 3.3. Optimal Harvest Level Determination

To achieve the commercial purpose of the fishery, it is necessary to harvest the population, and it is always hoped that the sustained ability can be achieved at a good economic profit. For the harvest problem, it is necessary to determine the controlled valuesE and τ and harvest level l, and this in general involves the optimization theory [48, 49].Letl be the harvest level, which is a decision variable. Theorems 3 and 4 show that system (4) admits an order-1 periodic solution when l≤l¯ or l>l¯ with τ≤yM1−1−q2EyL˜l. Since the harvest effort and yield of released predator are dependent on the harvest level, then it is assumed that El and τl take the following forms:(32)El=E1+E2−E1l−l1l2−l1,τl=τ1+τ2−τ1l−l1l2−l1,where l1≤l≤l2, l1 and l2 are minimum and maximum of the harvest level, and E1 (E2) and τ1 (τ2) are the harvest effort and yield of released predator at the harvest level l1 (l2).Letc1 and c2 be the unit selling prices of prey and predator, let c3 be the unit cost of harvest, and let c4 be the unit cost in breeding predator. Then the benefits from harvest can be described as Fbenefitl=c1q1Ell+c2q2ElηTl−c3E2l−c4τl. The objective is to maximize the unit benefits; that is,(33)maxFbenefitlTlsuch thatl1≤l≤l¯,l¯<l≤l¯,τl≤yM1−1−q2EyL˜l.


## 4. Numerical Simulations and Optimization

In this section, we compute some numerical simulations regarding the existence and stability of the periodic solution for the predator-prey model (4). It is quite difficult to verify the mathematical model simulations with realistic parameter values. We take a hypothetical set of parameter values to illustrate our analytical findings. The model parameters are as follows:b=0.7, d=0.2, c=0.005, p=0.1, h1=0.036, h2=1.44, e=0.44, and m=0.2. The control parameters are as follows: E1=0.2, E2=1, and τ1=0.


### 4.1. Numerical Simulations

Sincee>eM, the boundary equilibrium E100,0 is unstable, and system (3) has a positive equilibrium E∗. From equation (6), it can be observed that the fear effect factor k only affects the value y∗. Figure 2 illustrates the dependence of isoline dx/dt=0 on x for different fear level k. As illustrated, the positive equilibrium becomes stable from unstable with increasing of k. By equation (6), there is k∗=0.0376. To verify the theoretical results obtained in the above section, the simulations are implemented by considering different combinations of k, q1, q2, τ2, and l.Figure 2
The dependence of isolinedx/dt=0 on x for different fear level k with parameters given in numerical section. The red asterisk represents the positive equilibrium. For k=0,0.01,0.02,0.03, the positive equilibrium E∗ is unstable and, for k=0.04, the positive equilibrium E∗ is locally asymptotically stable.Case I:k=0.04, q1=0.8, and q2=0.6.I-(1):τ2=0. By Theorem 2, system (4) admits an order-1 semitrivial periodic solution for any l≤K=100, which is expressed by equation (10). Moreover, for l≤l¯=34.45, by Corollary 1, the semitrivial order-1 periodic solution is orbitally asymptotically stable. Meanwhile, for l=50%K, by Theorem 2, there is R0=0.5673<1; that is, the order-1 semitrivial periodic solution is orbitally asymptotically stable, as presented in Figure 3.Figure 3
The time series evolution of prey populationxt, predator population yt, and the phase portrait diagram of model (4) for τ2=0 with parameters given in numerical section. The phase portrait diagram shows that the semitrival order-1 periodic solution is orbitally asymptotically stable and predator population tends to extinction.I-(2):τ2>0. Firstly, for τ2=2, there is l¯=34.45%K. For l=25%K≤x∗, Theorem 4 and Theorem 7 indicate that system (4) admits a globally asymptotically stable positive order-1 periodic solution, as illustrated in Figure 1. Meanwhile, for l>x∗, function ϕN of system (4) for l=30%K, l=50%K, and l=60%K is presented in Figure 4. It can be observed that system (4) admits an order-1 periodic solution for l=30%K and l=50%K since the inequality τ<yM1−1−q2EyL˜l in Theorem 4 holds. When l=60%K, the direction of the inequality has changed, and the trajectory of system (4) will tend to the positive equilibrium E∗25,4.93 after finite impulses, as shown in Figure 5. But this is not always true; it is dependent on the value of τ2. As τ2 goes up to 9, the inequality τ>yM2−1−q2EyL˜l in Theorem 4 holds; then system (4) admits a positive order-1 periodic solution, as shown in Figure 6(a). However, for l=80%K, the direction of the inequality is changed again, and the trajectory of system (4) will tend to the positive equilibrium E∗25,4.93 after finite impulses, as shown in Figure 6(b).Figure 4
Illustration of the Poincaré mapϕN of system (4) for τ2=2 and different harvest level l with parameters given in numerical section: (a) l=30%K; (b) l=50%K; (c) l=60%K. The dotted line represents v=y and the intersection point is the fixed point of the Poincaré map ϕN. For l=30%K and l=50%K, the Poincaré map ϕN has a unique fixed point; that is, system (4) admits a unique order-1 periodic solution. For l=60%K, the Poincaré map ϕN does not have fixed point; that is, system (4) does not admit order-1 periodic solution.
(a)(b)(c)Figure 5
The time series evolution of prey populationxt, predator population yt, and the phase portrait diagram of model (4) for l=60%K and τ2=2 with parameters given in numerical section. The phase portrait diagram shows that the trajectory of system (4) will eventually tend to the positive equilibrium E∗. In this case, the positive equilibrium E∗ is globally asymptotically stable.Figure 6
The phase portrait diagram of model (4) forτ2=9 and different harvest level l with parameters given in numerical section: (a) l=60%K; (b) l=80%K. For l=60%K, there is τ2>τM2, and system (4) admits an order-1 periodic solution; for l=80%K, there is τ2∈ (τM1,τM2), and the trajectory of system (4) will tend to the positive equilibrium E∗ after finite impulses.
(a)(b)Case II.k=0.04 and q2=0.II-(1):τ2=0 and q1=0.8. Notice from Figure 1 that a higher catching rate for predators (e.g., q2=0.6) will cause the predator species extinction for l=50%K. When the catch for the predator is very small or ignored, that is, q2=0, the order-1 semitrivial periodic solution is unstable by Corollary 1 (i.e., R0=1.02>1). The function fNy=ϕNy−y is presented in Figure 7 and it can be observed that fNy>0, and, for any initial condition, the trajectory of system (4) will eventually tend to the positive equilibrium E∗25,4.93.Figure 7
Illustration of the successor functionfNy=ϕNy−y for q2=0 with parameters given in numerical section. The result indicates that the function is always greater than zero, which means that the Poincaré map ϕN of system (4) does not have fixed point and the trajectory will tend to the positive equilibrium E∗.II-(2):τ2>0. For τ2=8 and q1=0.8, by Theorem 3, system (4) admits a unique positive order-1 periodic solution for l=25%K. To show the existence of order-2 periodic solution, the catching rate for prey q1 is selected as a key parameter to verify how does the dynamic behavior of the system change. Function ϕN for different catching rate for prey is presented in Figure 8. It can be observed that system (4) admits a unique globally asymptotically stable order-1 periodic solution for a higher catching rate for prey, for example, q1=0.8, as shown in Figure 8(a). As the catching rate for prey goes down, for example q1=0.5, condition (ii) ϕN2yN0≥yN0 and μ21=1.2612>1 in Theorem 8 holds; then system (4) admits a stable order-2 periodic solution, as shown in Figure 8(b). Meanwhile, for q1=0.2, condition (i) ϕN2yN0<yN0 in Theorem 8 holds, and then system (4) admits an order-2 periodic solution, as shown in Figure 8(c).Figure 8
Illustration of the Poincaré mapsϕN and ϕN2 of system (4) for l=25%K and τ2=8 and different values of q1 with parameters given in numerical section. For q1=0.8, system (4) admits a unique orbitally asymptotically stable order-1 periodic solution; for q1=0.5, there is ϕN24.92>4.92 and μ1>1; in this case, system (4) admits a stable order-2 periodic solution; (c) for q1=0.2, there is ϕN24.92<4.92, and system (4) admits a stable order-2 periodic solution.
(a)(b)(c)Case III:k=0.01 and q1=q2=0.6.The positive equilibriumE∗ becomes unstable when k=0.01, and system (3) admits a limit cycle ΓLC. Since the existence and stability of semitrivial order-1 periodic solution for τ=0 do not depend on the fear factor k, the results are the same as those in Case I-(1) and are omitted hereby. So it mainly discusses the dynamic behavior of system (4) for τ2>0. It is easily checked that l¯=32.7%K.Firstly, forl=25%K≤x∗, by Theorem 3, system (4) admits a positive order-1 periodic solution for any τ, as shown in Figure 9.Figure 9
Illustration of the Poincaré mapϕN of system (4) for l=25%K and different values of τ2 with parameters given in numerical section. The results indicate that ϕN has a unique fixed point for any τ2; that is, system (4) admits a unique order-1 periodic solution.Next, let us considerl=50%K>l¯. It is easily checked that τM1=1.2495 (i.e., τ2=2.1867); then, by Theorem 4, there exists an order-1 periodic solution for system (4) when τ2≤2.1867. Here it should be pointed out that the condition given in Theorem 4 is only a sufficient one; in fact, as long as τ2≤2.55, system (4) admits an order-1 periodic solution, as illustrated in Figure 10.Figure 10
Illustration of the successor functionfNy of system (4) for l=50%K and different values of τ2 with parameters given in numerical section. The results indicate that fNy has a zero point in case of τ2<2.55; that is, system (4) admits an order-1 periodic solution when τ2<2.55.With an increase ofτ2, the existence of order-1 periodic solution cannot be guaranteed. For example, for any τ2∈2.56,3.58, there does not exist order-1 periodic solution for system (4). System (4) admits an order-3 periodic solution for τ2=2.9 (Figure 11), an order-5 periodic solution for τ2=3 (Figure 12), an order-2 periodic solution for τ2=3.1 (Figure 13), an order-1 periodic solution for τ2=3.6 (Figure 14), and an order-4 periodic solution for τ2=4.2 (Figure 15).Figure 11
Illustration of functionsϕN3y and ϕNy of system (4) for l=50%K and τ2=2.9 with parameters given in numerical section. The time series evolution of prey population xt, predator population yy, and the phase portrait diagram demonstrate the order-3 periodic solution.Figure 12
Illustration of functionsϕN5y and ϕNy of system (4) for l=50%K and τ2=3 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-5 periodic solution.Figure 13
Illustration of functionsϕN2y and ϕNy of system (4) for l=50%K and τ2=3.1 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-2 periodic solution.Figure 14
Illustration of functionϕNy of system (4) for l=50%K and τ2=3.6 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-1 periodic solution.Figure 15
Illustration of functionsϕN4y and ϕNy of system (4) for l=50%K and τ2=4.2 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-4 periodic solution.


### 4.2. Optimization

To achieve a good economic profit, it is necessary to find a levell∗ at which the benefits from harvest are maximal. Assume that c1=100, c2=5c1=500, and c3=20%c2=100. Denote σ≜c3/c1. In order to sustain the harvest, the releasing yield of predator should not be too large, so, in this part, it is assumed that τ2=1.Fork=0.04, system (4) admits an order-1 periodic solution when l≤70%K; the dependencies of period T and benefit Fbenefit on the harvest level l are presented in Figure 16. It can be seen that period T increases as l increases. Meanwhile the benefit function Fbenefit climbs up and then declines as l increases. For different σ, Fbenefit achieves its maximum at different lσ∗. When the cost of harvest is ignored, that is, σ=0, there is l0∗=40%K. As σ goes up, Fbenefit goes down. When σ=8, Fbenefit achieves its maximum at l8∗=50%K.Figure 16
The dependence of periodT and benefit Fbenefit on the harvest level l for k=0.04 with parameters given in numerical section.For the case ofk=0.01, system (4) admits an order-1 periodic solution when l≤l2, and the dependencies of period T and benefit Fbenefit on the harvest level l are presented in Figure 17. The benefit funtion Fbenefit climbs up and then declines as l increases. For different σ, Fbenefit almost achieves its maximum at lσ∗=50%K.Figure 17
The dependence of periodT and benefit Fbenefit on the harvest level l for k=0.01 with parameters given in numerical section.


## 4.1. Numerical Simulations

Sincee>eM, the boundary equilibrium E100,0 is unstable, and system (3) has a positive equilibrium E∗. From equation (6), it can be observed that the fear effect factor k only affects the value y∗. Figure 2 illustrates the dependence of isoline dx/dt=0 on x for different fear level k. As illustrated, the positive equilibrium becomes stable from unstable with increasing of k. By equation (6), there is k∗=0.0376. To verify the theoretical results obtained in the above section, the simulations are implemented by considering different combinations of k, q1, q2, τ2, and l.Figure 2
The dependence of isolinedx/dt=0 on x for different fear level k with parameters given in numerical section. The red asterisk represents the positive equilibrium. For k=0,0.01,0.02,0.03, the positive equilibrium E∗ is unstable and, for k=0.04, the positive equilibrium E∗ is locally asymptotically stable.Case I:k=0.04, q1=0.8, and q2=0.6.I-(1):τ2=0. By Theorem 2, system (4) admits an order-1 semitrivial periodic solution for any l≤K=100, which is expressed by equation (10). Moreover, for l≤l¯=34.45, by Corollary 1, the semitrivial order-1 periodic solution is orbitally asymptotically stable. Meanwhile, for l=50%K, by Theorem 2, there is R0=0.5673<1; that is, the order-1 semitrivial periodic solution is orbitally asymptotically stable, as presented in Figure 3.Figure 3
The time series evolution of prey populationxt, predator population yt, and the phase portrait diagram of model (4) for τ2=0 with parameters given in numerical section. The phase portrait diagram shows that the semitrival order-1 periodic solution is orbitally asymptotically stable and predator population tends to extinction.I-(2):τ2>0. Firstly, for τ2=2, there is l¯=34.45%K. For l=25%K≤x∗, Theorem 4 and Theorem 7 indicate that system (4) admits a globally asymptotically stable positive order-1 periodic solution, as illustrated in Figure 1. Meanwhile, for l>x∗, function ϕN of system (4) for l=30%K, l=50%K, and l=60%K is presented in Figure 4. It can be observed that system (4) admits an order-1 periodic solution for l=30%K and l=50%K since the inequality τ<yM1−1−q2EyL˜l in Theorem 4 holds. When l=60%K, the direction of the inequality has changed, and the trajectory of system (4) will tend to the positive equilibrium E∗25,4.93 after finite impulses, as shown in Figure 5. But this is not always true; it is dependent on the value of τ2. As τ2 goes up to 9, the inequality τ>yM2−1−q2EyL˜l in Theorem 4 holds; then system (4) admits a positive order-1 periodic solution, as shown in Figure 6(a). However, for l=80%K, the direction of the inequality is changed again, and the trajectory of system (4) will tend to the positive equilibrium E∗25,4.93 after finite impulses, as shown in Figure 6(b).Figure 4
Illustration of the Poincaré mapϕN of system (4) for τ2=2 and different harvest level l with parameters given in numerical section: (a) l=30%K; (b) l=50%K; (c) l=60%K. The dotted line represents v=y and the intersection point is the fixed point of the Poincaré map ϕN. For l=30%K and l=50%K, the Poincaré map ϕN has a unique fixed point; that is, system (4) admits a unique order-1 periodic solution. For l=60%K, the Poincaré map ϕN does not have fixed point; that is, system (4) does not admit order-1 periodic solution.
(a)(b)(c)Figure 5
The time series evolution of prey populationxt, predator population yt, and the phase portrait diagram of model (4) for l=60%K and τ2=2 with parameters given in numerical section. The phase portrait diagram shows that the trajectory of system (4) will eventually tend to the positive equilibrium E∗. In this case, the positive equilibrium E∗ is globally asymptotically stable.Figure 6
The phase portrait diagram of model (4) forτ2=9 and different harvest level l with parameters given in numerical section: (a) l=60%K; (b) l=80%K. For l=60%K, there is τ2>τM2, and system (4) admits an order-1 periodic solution; for l=80%K, there is τ2∈ (τM1,τM2), and the trajectory of system (4) will tend to the positive equilibrium E∗ after finite impulses.
(a)(b)Case II.k=0.04 and q2=0.II-(1):τ2=0 and q1=0.8. Notice from Figure 1 that a higher catching rate for predators (e.g., q2=0.6) will cause the predator species extinction for l=50%K. When the catch for the predator is very small or ignored, that is, q2=0, the order-1 semitrivial periodic solution is unstable by Corollary 1 (i.e., R0=1.02>1). The function fNy=ϕNy−y is presented in Figure 7 and it can be observed that fNy>0, and, for any initial condition, the trajectory of system (4) will eventually tend to the positive equilibrium E∗25,4.93.Figure 7
Illustration of the successor functionfNy=ϕNy−y for q2=0 with parameters given in numerical section. The result indicates that the function is always greater than zero, which means that the Poincaré map ϕN of system (4) does not have fixed point and the trajectory will tend to the positive equilibrium E∗.II-(2):τ2>0. For τ2=8 and q1=0.8, by Theorem 3, system (4) admits a unique positive order-1 periodic solution for l=25%K. To show the existence of order-2 periodic solution, the catching rate for prey q1 is selected as a key parameter to verify how does the dynamic behavior of the system change. Function ϕN for different catching rate for prey is presented in Figure 8. It can be observed that system (4) admits a unique globally asymptotically stable order-1 periodic solution for a higher catching rate for prey, for example, q1=0.8, as shown in Figure 8(a). As the catching rate for prey goes down, for example q1=0.5, condition (ii) ϕN2yN0≥yN0 and μ21=1.2612>1 in Theorem 8 holds; then system (4) admits a stable order-2 periodic solution, as shown in Figure 8(b). Meanwhile, for q1=0.2, condition (i) ϕN2yN0<yN0 in Theorem 8 holds, and then system (4) admits an order-2 periodic solution, as shown in Figure 8(c).Figure 8
Illustration of the Poincaré mapsϕN and ϕN2 of system (4) for l=25%K and τ2=8 and different values of q1 with parameters given in numerical section. For q1=0.8, system (4) admits a unique orbitally asymptotically stable order-1 periodic solution; for q1=0.5, there is ϕN24.92>4.92 and μ1>1; in this case, system (4) admits a stable order-2 periodic solution; (c) for q1=0.2, there is ϕN24.92<4.92, and system (4) admits a stable order-2 periodic solution.
(a)(b)(c)Case III:k=0.01 and q1=q2=0.6.The positive equilibriumE∗ becomes unstable when k=0.01, and system (3) admits a limit cycle ΓLC. Since the existence and stability of semitrivial order-1 periodic solution for τ=0 do not depend on the fear factor k, the results are the same as those in Case I-(1) and are omitted hereby. So it mainly discusses the dynamic behavior of system (4) for τ2>0. It is easily checked that l¯=32.7%K.Firstly, forl=25%K≤x∗, by Theorem 3, system (4) admits a positive order-1 periodic solution for any τ, as shown in Figure 9.Figure 9
Illustration of the Poincaré mapϕN of system (4) for l=25%K and different values of τ2 with parameters given in numerical section. The results indicate that ϕN has a unique fixed point for any τ2; that is, system (4) admits a unique order-1 periodic solution.Next, let us considerl=50%K>l¯. It is easily checked that τM1=1.2495 (i.e., τ2=2.1867); then, by Theorem 4, there exists an order-1 periodic solution for system (4) when τ2≤2.1867. Here it should be pointed out that the condition given in Theorem 4 is only a sufficient one; in fact, as long as τ2≤2.55, system (4) admits an order-1 periodic solution, as illustrated in Figure 10.Figure 10
Illustration of the successor functionfNy of system (4) for l=50%K and different values of τ2 with parameters given in numerical section. The results indicate that fNy has a zero point in case of τ2<2.55; that is, system (4) admits an order-1 periodic solution when τ2<2.55.With an increase ofτ2, the existence of order-1 periodic solution cannot be guaranteed. For example, for any τ2∈2.56,3.58, there does not exist order-1 periodic solution for system (4). System (4) admits an order-3 periodic solution for τ2=2.9 (Figure 11), an order-5 periodic solution for τ2=3 (Figure 12), an order-2 periodic solution for τ2=3.1 (Figure 13), an order-1 periodic solution for τ2=3.6 (Figure 14), and an order-4 periodic solution for τ2=4.2 (Figure 15).Figure 11
Illustration of functionsϕN3y and ϕNy of system (4) for l=50%K and τ2=2.9 with parameters given in numerical section. The time series evolution of prey population xt, predator population yy, and the phase portrait diagram demonstrate the order-3 periodic solution.Figure 12
Illustration of functionsϕN5y and ϕNy of system (4) for l=50%K and τ2=3 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-5 periodic solution.Figure 13
Illustration of functionsϕN2y and ϕNy of system (4) for l=50%K and τ2=3.1 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-2 periodic solution.Figure 14
Illustration of functionϕNy of system (4) for l=50%K and τ2=3.6 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-1 periodic solution.Figure 15
Illustration of functionsϕN4y and ϕNy of system (4) for l=50%K and τ2=4.2 with parameters given in numerical section. The time series evolution of prey population xt, predator population yt, and the phase portrait diagram demonstrate the order-4 periodic solution.


## 4.2. Optimization

To achieve a good economic profit, it is necessary to find a levell∗ at which the benefits from harvest are maximal. Assume that c1=100, c2=5c1=500, and c3=20%c2=100. Denote σ≜c3/c1. In order to sustain the harvest, the releasing yield of predator should not be too large, so, in this part, it is assumed that τ2=1.Fork=0.04, system (4) admits an order-1 periodic solution when l≤70%K; the dependencies of period T and benefit Fbenefit on the harvest level l are presented in Figure 16. It can be seen that period T increases as l increases. Meanwhile the benefit function Fbenefit climbs up and then declines as l increases. For different σ, Fbenefit achieves its maximum at different lσ∗. When the cost of harvest is ignored, that is, σ=0, there is l0∗=40%K. As σ goes up, Fbenefit goes down. When σ=8, Fbenefit achieves its maximum at l8∗=50%K.Figure 16
The dependence of periodT and benefit Fbenefit on the harvest level l for k=0.04 with parameters given in numerical section.For the case ofk=0.01, system (4) admits an order-1 periodic solution when l≤l2, and the dependencies of period T and benefit Fbenefit on the harvest level l are presented in Figure 17. The benefit funtion Fbenefit climbs up and then declines as l increases. For different σ, Fbenefit almost achieves its maximum at lσ∗=50%K.Figure 17
The dependence of periodT and benefit Fbenefit on the harvest level l for k=0.01 with parameters given in numerical section.


## 5. Conclusion and Discussions

In this paper, we have discussed the dynamics of a harvested prey-predator model, where the prey is provided with fear effect. For the system without harvest (3), there exists a unique positive equilibrium. To verify the stability of the equilibrium, a critical level of fear factork∗ is characterized (i.e., equation (6)). When the impact of fear on prey is small, that is, 0≤k<k∗, the positive equilibrium is unstable and a limit cycle exists. As the impact of fear grows and exceeds k∗, the positive equilibrium becomes stable and the limit cycle disappears (Figure 2). In any case, predators coexist with prey and the system is persistent.For the system with harvest (4), if we do not consider the release of predator (i.e.,τ=0), system (4) admits a semitrivial order-1 periodic solution for any harvest level (Figure 3). Moreover, the semitrivial order-1 periodic solution is orbitally stability when the harvest level is not higher than the first equilibrium component (i.e., l≤x∗). Meanwhile, for the case of l>x∗, the semitrivial order-1 periodic solution is orbitally stable when a strong harvest intensity is implemented. This means that the system can be disrupted, and predators will go extinct if the harvest is not properly planned. To maintain the ecological health and avoid the extinction of predator populations, it is necessary to reduce the catch rate of predators (Figure 7) or release a certain quantity of predator pups. In the second case, that is, τ>0, system (4) admits a positive order-1 periodic solution when l≤l¯ (Figures 1, 4, and 9). Moreover, the order-1 periodic solution is orbitally asymptotically stable and globally attractive when l≤x∗ and τ≤τf. Meanwhile, for τ>τf, system (4) admits an order-2 periodic solution in case of ϕN2yN0≥yN0 and μ21>1 (Figure 8). In case of l>l¯, system (4) also admits a positive order-1 periodic solution for τ≤τM1 (Figures 4 and 10) or τ≥τM2 (Figure 6). But, for τ∈τM1,τM2, when k>k∗, the trajectory of system (4) will tend to the positive equilibrium E∗25,4.93 after finite impulses (Figure 5). Meanwhile, in case of k<k∗, the dynamic behavior of system (4) depends heavily on parameter τ. For different value of τ, system (4) may admit an order-kk=1,2,3,4,5 periodic solution (Figures 11–15).To achieve a good economic profit, optimization withτ2=1 is carried out and the results show that the benefits from harvest depend on the unit selling prices of prey and predator, as well as the unit cost of harvest. For given c1=100, c2=500, and c4=100, the benefit function first climbs up and then declines as l increases. For k=0.04, the economic profit Fbenefit achieves its maximum at l∗=40%K when σ=2. As σ goes up, Fbenefit goes down, and Fbenefit achieves its maximum at l∗=50%K when σ=8. Meanwhile, in case of k=0.01, Fbenefit almost achieves its maximum at l∗=50%K.



---

*Source: 9496599-2022-01-11.xml*