Ecological dynamics of age selective harvesting of fish population: Maximum sustainable yield and its control strategy

Chaos, Solitons and Fractals 93 (2016) 111–122

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Ecological dynamics of age selective harvesting of fish population: Maximum sustainable yield and its control strategy Debaldev Jana a, Rashmi Agrawal b, Ranjit Kumar Upadhyay b,∗, G.P. Samanta c a

Department of Mathematics, SRM University, Kattankulathur, Tamil Nadu-603 203, India Department of Applied Mathematics, Indian Institute of Technology (Indian School of Mines), Dhanbad - 826 004, Jharkhand, India c Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah-711 103, India b

a r t i c l e

i n f o

Article history: Received 18 July 2016 Revised 21 September 2016 Accepted 23 September 2016

Keywords: Age selective harvesting Maximum sustainable yield Bionomic equilibrium Optimal harvesting

a b s t r a c t Life history of ecological resource management and empirical studies are increasingly documenting the impact of selective harvesting process on the evolutionary stable strategy of both aquatic and terrestrial ecosystems. In the present study, the interaction between population and their independent and combined selective harvesting are framed by a multi-delayed prey-predator system. Depending upon the age selection strategy, system experiences stable coexistence to oscillatory mode and vice versa via Hopfbifurcation. Economic evolution of the system which is mainly featured by maximum sustainable yield (MSY), bionomic equilibrium and optimal harvesting vary largely with the commensurate age selections of both population because equilibrium population abundance becomes age-selection dependent. Our study indicates that balance between harvesting delays and harvesting intensities should be maintained for better ecosystem management. Numerical examples support the analytical findings. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Age or body size is generally considered to be one of the most important and significant traits of a species because it correlates with many aspects of its biology, from life history to ecology [7,15,36,50]. The views of fisheries management over the last few decades, have oscillated between alarm and trust in management progress. The predominant scientific policies for remedying the world fishing crisis aim at maximum sustainable yield (MSY) by selective harvesting which is by adjusting gear selectivity, fishing effort etc. Law and Grey [30] by a theoretical work on northeast Arctic cod pointed out the concept of fisheries-induced evolution and its effects on the yield of an exploited fish stock. In this study they show how fishing pressure creates a strong selection for individuals with early maturation, leading to a change in the age at maturation of the stock, towards earlier maturation. Correlated life-history traits in response to size-selective harvesting have been demonstrated in laboratory populations of the Atlantic silverside, Menidia menidia [10,52]. Different case studies on Hilsa fish at Hooghly estuary say a real scenario of ecological and economical status of India and Bangladesh. De and Datta [12] did an experiment of Indian shad,

∗

Corresponding author. Fax: +913262296559. E-mail address: [email protected] (R.K. Upadhyay).

http://dx.doi.org/10.1016/j.chaos.2016.09.021 0960-0779/© 2016 Elsevier Ltd. All rights reserved.

Tenualosa ilisha (Hamilton) commonly known as Hilsa, of the freshwater zone of Hooghly estuary by using the length- frequency method. They used different techniques e.g. Von Bertalanffy growth equation to estimate the length-weight relationship of Hilsa fish. It is well known that the rate of growth of fish population may vary from one environment to other or in the same environment from year to year due to changes in ecology including changes in food availability, density dependent growth factors etc [38,22,43,39,44,42]. The low yield of Hilsa in present day situation has arrived a question to find out the factors like indiscriminate killing of juveniles, establishment of Farakha Barrage, decreasing depth of estuary, pollution state [3]. Also, it is observed by the empirical data on Serampore-Uttarpara belt [3] that there is a very important biological as well as economic impact of lengthweight relationship of Hilsa fish on fishery management. Over fishing and over exploitation of juveniles play a very negative role in any economically profitable fishery which in the broad sense stress on the biological and economic equilibrium and that also emphasize on the evolutionary consequences of the population and its environment. Gear selection for particular size of fish [3], perfect boat choice [20], government policies to prevent juvenile fish and protect overfishing [3,20] etc are the prominent solutions and real application to prevent ecosystem management and biodiversity. Presently, management objectives emphasize not only the biological yield of fisheries (MSY) but also include economic and social considerations. Many coastal rural communities are solely de-

112

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

pendent on particular fisheries [3,22,38,39,42–44]. In most cases, fisheries management is conducted for regulating fishing mortality [20,46]. Methods to regulate fishing mortality include direct methods such as size and catch limits and indirect methods such as seasonal closures, area closures, gear restrictions and licensing [3]. In direct methods, the most important is catch limits (number of tons per year). Among the indirect methods, the amount of time that fishermen spends to harvest a species is reduced by seasonal closing of fishery which actually reduces the total harvest. But this method has some drawbacks, sometimes it may encourage fishermen to give more effort during the fishing season and virtually the total harvest is not reduced. For the purpose of economically profitable as well as biologically beneficial system, bionomic equilibrium gives us constraints to have profitable harvesting [31,17,18]. For the interest of getting MSY and bionomic equilibrium of a profitable fishery, optimization is a technique which maximizes the profit without any harmful effect on the system. Optimal harvesting policy gives us a strategy which helps us in both maximizing the profit by minimum level of effort as well as eradicating the risk of extinction of the harvested population [27,14]. Proper harvesting or age/ size selective harvesting is a very significant methodology that dampen fluctuations by making overcompensatory dynamics under compensatory [6,8,9,25,33]. Combined effect of harvesting and delay on prey-predator system are studied by several researchers [5,32,53,54]. Recently, Jana et al. [21] has been established a age-selective harvesting model by using the technique of Arino et al. [2]. For this consequence in the present work, we consider two delays: one in the fish prey harvesting term and the other in fish predator harvesting term. Thus, we are considering selective harvesting of both groups species. Biologically, it implies that we are restricting harvesting of species below some certain age so that they can grow up to some specific size or age and thus protecting juvenile fish population. According to new paradigm of ecosystem based fisheries management (EBFM), selective fishing is widely encouraged in contrary to nonselective fishing which has many adverse impacts [37,41]. Zhou et al. [56] proposed a balanced exploitation approach and commented that selective harvesting should be employed with reduced exploitation rates. The objective of this article is to study the effect of selective harvesting of either or both species on the dynamics of the system. It is also our intention to verify rationality of some of the above observations on EBFM. The organization of the paper is as follows: In Section 2, we formulate the mathematical model. Also, boundedness, equilibria analysis and existence of Hopf-bifurcation for different cases are derived in this Sect. In Section 3, we present the sustainable yield and MSY policy of the harvested population. Bionomic equilibrium of the model is presented in Section 4. In Section 5, we have presented optimal harvesting policy of the model system. Finally, a brief discussion has been drown from the manuscript in Section 6. 2. Formulation of mathematical model and its analysis The Rosenweig-McArthur type prey-predator system was studied by Kuang and Freedman [29] which shows an unique limit cycle in the positive quadrant. The case of instantaneous harvesting of both the population has been studied by Srinivasu et al. [49] through the following model:





x dx β xy = αx 1 − − − E1 x, dt k 1 + ax dy cβ xy = − γ y − E2 y, dt 1 + ax

y is the predator population density. The predation process obeys the Holling type II functional response with attack rate β , attack rate multiplied by handling time is a, food conversion efficiency c and food independent mortality rate γ . Also, it is assumed that both prey and predator are harvested instantaneously by the rates E1 and E2 respectively. It is well known that the rate of change of the population depends on three components: growth, death, and intraspecific competition (crowding or direct interference). The death and the intraspecific competition rates together known as the decline rate. We assume that the decline rate is instantaneous and that the death rate is given by a linear term whereas the intraspecific competition rate is given by the quadratic term. For this, we can write the classical logistic ODE model [2] as follows:

dN (t ) = bN (t ) − μN (t ) − δ N 2 (t ), dt

(2)

where b, μ and δ are the birth, death and death due to intraspecific competition respectively of the population N(t) and we choose α = b − μ (i.e. intrinsic growth rate of N(t)). In this equation, the decline rate is given by {μN (t ) + δ N 2 (t )} and the growth rate is given by bN(t). Both are assumed to depend on the present population size. We continue to assume that the decline rate is instantaneous, but we now assume that the growth rate is proportional to the number of individuals in the population (t − τ ) time units in the past, that manage to survive until time t [2]. To obtain an equation that describes how many individuals alive at time (t − τ ) are still alive at time t, we solve the following first order ODE for N(t) as a function of N (t − τ ):

dN (t ) = −μN (t ) − δ N 2 (t ). dt

(3)

Using the technique of separation of variables and integrating both sides from (t − τ ) to t, we obtain [2]

N (t ) =

μN (t − τ ) . μeμτ + δ (eμτ − 1 )N (t − τ )

(4)

Now, if the population does not posses the death due to intraspecific competition (i.e. δ = 0), then the survival rate of population N(t) is given by

N (t ) = e−μτ N (t − τ ).

(5)

Since we assume that growth in the population at time t is proportional to those individuals alive at time (t − τ ) who survive until time t, we replace x and y in the monomial with coefficient E1 and E2 respectively in (1) by (4) and (5). According to the assumption of selective harvesting, our model system (1) becomes





x E1 μx(t − τ1 ) dx β xy = αx 1 − − − , dt k 1 + ax μeμτ1 + δ (eμτ1 − 1 )x(t − τ1 ) dy cβ xy = − γ y − E2 e−γ τ2 y(t − τ2 ), dt 1 + ax

(6)

where α = b − μ (b and μ are linear birth and death rate of prey(x) population), δ = αk . τi (i = 1, 2 ) are the delays regarding age selection for harvesting of x and y population respectively. The initial conditions of (6) are given as

x ( θ ) = φ1 ( θ ) ≥ 0 , y ( θ ) = φ 2 ( θ ) ≥ 0 ,

θ ∈ [−τ , 0], φi (0 ) > 0 (i = 1, 2 ), where φ : [−τ , 0] → 2 are continuous and τ = max[τ1 , τ2 ] with norm

(1)

where x is the logistically growing prey population density with intrinsic growth rate α to its environmental carrying capacity k and

||φ|| = sup {|φ1 (θ )|, |φ2 (θ )|}, −τ ≤θ ≤0

such that φ = (φ1 , φ2 ).

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

113

(i)

(ii)

16

20

14 12 10

x

*

y*

15 8 10

6

5 10

10

4 10 10 5

5

5

5

Fig. 1. Equilibrium biomass: (i) for x∗ and (ii) for y∗ of the system (6) for different combinations of τ1 − τ2 . Other parameters are as in the text.

Remark. Considering the assumption under which this model has been derived, it only makes sense to consider initial data that is positive over the entire interval. Under this assumption, solutions of (6) remain positive. Boundedness, whether a situation has a clearly defined beginning or end. The boundedness may be interpreted as a natural restriction to growth as a consequence of limited resources. The right hand side of system (6) is completely continuous and locally Lipschitzian on C ([−τ , 0]; R2 ; τ = max(τ1 , τ2 )), the Banach space of continuous mapping from [−τ , 0] → R2 , the solution (x(t), y(t)) of (6) with initial conditions x(θ ) > 0, y(θ ) > 0, ∀θ ∈ [−τ , 0] exists and is unique on [0, ζ ), where 0 < ζ ≤ +∞. From the previous discussions, it is evident that x(t) > 0, y(t) > 0, ∀t ≥ 0. From system (6), it follows that



x dx ≤ αx 1 − dt k



⇒ lim sup x(t ) ≤ k. t→∞

If possible, let us assume that y(t) ≥ M, ∀ t ≥ T, where M and T are sufficiently large positive numbers. Therefore,

cβ ky dy ≤ −γM dt 1 + ak



⇒ y(t )≤y(0 ) exp

 =









cβ kt γ M (1 + ak ) cβ kt + 1− exp 1 + ak cβ k 1+ak







γ M (1 + ak ) cβ kt γ M (1 + ak ) y (0 ) − exp + < 0, cβ k 1 + ak cβ k

for sufficiently large values of t, a contradiction. This implies y(t) is also bounded. Hence the system (6) is ecologically well-posed. The system (6) has only one coexistence equilibrium point which is given by E∗ (x∗ , y∗ ), where

γ + E2 e−γ τ2 , cβ − (γ + E2 e−γ τ2 )a    x∗ E1 μ 1 + ax∗ y∗ = α 1− − . β k μeμτ1 + δ (eμτ1 − 1 )x∗ x∗ =

(7)

Feasibility or biologically positivity studies aim to objectively and rationally uncover the strengths and weaknesses of an existing proposed venture, opportunities and threats present in the environment, the resources require to carry through and ultimately the

prospects for success. Therefore, the coexistence equilibrium abundance of the system (6) is feasible if

(i ) cβ > aγ , cβ − aγ (ii ) E2 > , a

(iii ) τ2 > ( i v ) τ1 >

1

γ 1

μ

log

E2 a , cβ − aγ

 log

1 μ + δ x∗

 δ x∗ +

E 1 μk α ( k − x∗ )



.

We choose parameters as α = 1, k = 60, β = 0.1, a = 0.08, c = 0.7, γ = 0.33, μ = 0.01, E1 = 0.5, E2 = 0.2. First we plot the equilibrium biomass (x∗ in Fig. 1(i) and y∗ in Fig. 1(ii)) of the system (6) with different combinations of τ1 − τ2 . This figure depicts the total scenario of the dependence of the equilibrium biomass of the system (6) on the age selection (τ 1 and τ 2 ) of prey and predator population respectively. With increasing the length of prey age selection, predator population’s equilibrium biomass increases whereas increasing the length of predator age selection decreases the equilibrium biomass of prey population. Linearizing the system (6) at E∗ (x∗ , y∗ ) and its corresponding characteristic equation is given by

λ2 + Aλ + B + (C1 λ + D1 )e−λτ1 + (C2 λ + D2 )e−λτ2 + Ee−λ(τ1 +τ2 ) = 0,

(8)

where A = −(a11 + a22 ), B = (a11 a22 − a12 a21 ), C1 = −b11 , D1 = a22 ∗ b11 , C2 = −c22 , D2 = a11 c22 , E = b11 c22 , a11 = α (1 − 2kx ) − β y∗ , (1+ax∗ )2

∗

βx a12 = − 1+ ax∗ ,

E μ2 eμτ1 − {μeμτ1 +1δ (eμτ1 −1 )x∗ }2 ,

∗

cβ y a21 = (1+ , ax∗ )2

a22 =

c β x∗ 1+ax∗

−γ,

b11 =

c22 = −E2 e−γ τ2 . The system (6) is stable,

unstable and experiences Hopf bifurcation according as Re(λ) < 0, > 0 and = 0 respectively. In Fig. 2, We plot the stability region in τ 1 τ 2 -plane by using Eq. (8) in which black line represents the Hopf bifurcation line where system (6) switches its stability and in the stable and unstable region, system (6) converges to its coexistence equilibrium point E∗ and near about stable limit cycle of E∗ respectively. When τi = 0, i = 1, 2, the system (6) is stable at E∗ (Appendix A), then we move on with the system with two different ways to observe the dynamical behavior of the system (6). First we consider Eq. (8) with instantaneous prey harvesting, i.e., τ1 = 0 and then we obtain a critical value

114

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

3. Sustainable yield and MSY policy

Fig. 2. Stability regions of the system (6) are depicted in τ 1 τ 2 -plane where the system shows stable and unstable dynamics in ‘Stable Region’ and ‘Unstable Region’ respectively and stable and unstable regions are separated by the line which is denoted by ‘Hopf bifurcation line’ on which the system switches its stability through Hopf-bifurcation. Parameters are as in the text.

of predator age selection τ 2 as τ20 , here the system is stable if τ2 < τ20 , unstable if τ2 > τ20 and system experiences Hopf bifurcation when τ2 = τ20 . Now, we stay with τ 2 in its stable interval (τ2 ∈ [0, τ20 )), suppose parameters satisfy conditions of the Theorem 2.1, H3 and H5 (given in Appendix A), then the equilibrium E∗ (x∗ , y∗ ) is asymptotically stable when τ1 ∈ (0, τ10 ), unstable when τ1 > τ10 and a Hopf bifurcation occurs when τ1 = τ10 , where

τ1 0 =

1 M 1 cos−1 . ω0 N1

Detailed proof is given in Appendix A. For the total dynamical feature of the system (6) with respect to τ 1 when τ2 = 1.5, we obtain a complete graphical representation of the above case. Fig. 3(i) plot the functions T0 (τ 1 ) and Fig. 3 (ii) is the bifurcation diagram of the system (6) with respect to τ 1 when τ2 = 1.5 in the three-dimensional space (τ 1 , x, y). Fig. 3(iii) and (iv) are the time evolutions of system (6), if we choose τ1 = 0.27 < 0.338, then the system (6) is stable, i.e., all the solutions of system (6) with initial population strength (30, 5) converge to coexistence equilibrium E∗ and if we select τ1 = 0.4 > 0.338, then the system (6) is unstable, i.e., all the solutions of system (6) with initial population strength (30, 5) converge to the near about stable limit cycle of coexistence equilibrium E∗ . Therefore, the system (6) experiences a Hopf bifurcation when τ 1 crosses the critical value τ10 = 0.338. From Fig. 3(ii), graphically it is clear that amplitude of the oscillation increases with increasing delay parameter. In the second way, we consider Eq. (8) with instantaneous predator harvesting, i.e., τ2 = 0 and then we get a critical value of prey age selection τ 1 as τ˜10 , here the system is stable if τ1 < τ˜10 , unstable if τ1 > τ˜10 and system experiences Hopf bifurcation when τ1 = τ˜10 . Now we stay with τ 1 in its stable interval (τ1 ∈ [0, τ˜10 )) and regard τ 2 as a parameter (τ 2 ∈ (0, ∞)). Suppose parameters satisfy conditions of the Theorem 2.1, H3 and H5 (given in Appendix A). Then the equilibrium E∗ (x∗ , y∗ ) is locally asymptotically stable when τ2 ∈ (0, τ t20 ), unstable when τ2 > τ t20 and a Hopf bifurcation occurs when τ2 = τ t20 , where

τ t 20 =

 1  M 2 cos−1 . ω0 N2

Detailed proof is given in Appendix A.

The sustainable yield of the harvested population is the ecological yield that can be harvested without reducing the population drastically i.e. the population required to maintain ecosystem health at the same or increasing level over time. This yield usually varies over time with the needs of the ecosystem to maintain itself. This yield depends on many constraints [9,51,35,28]. As an example, NOAA’s [55] report shows that in fishery during oil spill there is drastic decrease in harvested fish population. During that time the sustainable yield is very less as the ecological yield is needed by the fish population to re-establish itself [55]. In case of overfishing, we can harvest huge amount of fish population of all age classes randomly, but this harvesting policy is short term. It can be cured by harvesting a certain amount of fish population after a particular age class of the fish population [48,45]. By this method we can overcome the extinction risk of the population. This is the main goal of maximum sustainable yield policy by which we can get satisfactory amount of yield without imposing harsh biological effect on the harvested population for a long period of time [55,45]. τ τ In this case, the total yield (Y (E11 , E22 )) at equilibrium E∗ (x∗ , y∗ ) is

Y (E1τ1 , E2τ2 ) = E1 x∗ + E2 y∗ E1 (γ + E2 e−γ τ2 ) cE2 = + cβ − (γ + E2 e−γ τ2 )a cβ − (γ + E2 e−γ τ2 )a    γ + E2 e−γ τ2 × α 1− k{cβ − (γ + E2 e−γ τ2 )a}



E1 μ(cβ − (γ + E2 e−γ τ2 )a ) μeμτ1 (cβ − (γ + E2 e−γ τ2 )a ) + δ (eμτ1 − 1 )(γ + E2 e−γ τ2 )    p E 1 μq pE1 cE2 = + α 1− − , (9) q q qk qμeμτ1 + δ (eμτ1 − 1 ) p −

where p = (γ + E2 e−γ τ2 ) and q = cβ − pa. Yield management is a variable pricing strategy, based on understanding, anticipating and influencing consumer behavior in order to maximize revenue or profits from a fixed, perishable resource [34]. In population ecology and economics, maximum sustainable yield (MSY) is theoretically the largest yield or catch that can be taken from the species stock over a long time period. By an observation, Rothschild et al. [45] point out that with respect to the management goals, the commonly used age-based production model (ABPM) which is devised by Shepherd [47] gives a problematic perception of stock abundance relative to maximum sustainable yield (MSY). Therefore, a management system needs to be designed that creates feasible goals and uses improved and simpler metrics to measure better performance [45]. Now, we depict graphically the total yield-effort perspective of our system (6). First, we fix the age selections of both prey and predator respectively, numerically we choose τ1 = 0.4, τ2 = 0.5. We τ τ plot sustainable yield Y (E11 , E22 ) with varying prey and predator harvesting effort in 0 < Ei < 0.5, (i = 1, 2 ) (Fig. 4(i)). By this graph we obtain two different observations which are shown in Fig. 4(ii) and (iii) separately. First, the graph shows that the yield of the system increases gradually with increasing prey harvesting effort (E1 ) for a constant value of predator harvesting effort (E2 ). This argument is verified by Fig. 4(iii), i.e., choose E2 = 0.2, we see that yield (Y) increases gradually with increasing E1 . Second, at any degree of prey harvesting effort (E1 ), we see that sustainable yield of the system initially increases up to a certain level of the strength of predator harvesting effort (E2 ) and after that it gradually decreases. Therefore, it is clear that the sustainable yield of the system (6) reaches at its maximum level (maximum sustainable yield,

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

115

Fig. 3. Total dynamical feature of the system (6) with respect to τ 1 when τ2 = 1.5. (i) is the plot of the functions T0 (τ 1 ). (ii) is the bifurcation diagram of the system (6) with respect to τ 1 when τ2 = 1.5 in the three-dimensional space (τ 1 , x, y). (iii) and (iv) are the time evolutions for τ1 = 0.27 < 0.338 (system is stable) and τ1 = 0.4 > 0.338 (system is unstable) respectively. A Hopf bifurcation exists at τ1 = 0.338. Other parameters are as in the text.

Fig. 4. Yield-effort scenario of system (6) with a particular age selection of both population τ1 = 0.4, τ2 = 0.5. (i) Variation of sustainable yield (Y (E1τ1 , E2τ2 )) w.r.t the variation of harvesting efforts Ei , (i = 1, 2 ), figure shows that yield attains it’s maximum value (maximum sustainable yield - MSY) w.r.t E2 which is pointed by black colored line but not by E1 . (ii) portraits of last figure at E1 = 0.2 which shows that at E2 = 0.413 sustainable yield of the system maximizes (MSYE2 ) and (iii) portraits of figure (i) at E2 = 0.2 which shows that there is no MSY with the variation of E1 . Other parameters are as in the text.

116

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

Fig. 5. Characteristic variation of economic rent or net revenue R at any time at a constant harvesting effort of both population (E1 = 0.4, E2 = 0.5) w.r.t the simultaneous variation of population abundance (x − y). Red, black, blue and green regions say that harvesting of both is not profitable, harvesting of prey is not profitable but harvesting of predator is profitable, harvesting of predator is not profitable but harvesting of prey is profitable and harvesting of both is profitable. Parameters are: c1 = 2.5, c2 = 3.5, p1 = 2.5, p2 = 3, E1 = 0.5 and E2 = 0.4. Other parameters are as in the text. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

MSY) at a critical value of E2 , here we denote it by black line in Fig. 4(i). This observation is also verified by Fig. 4(ii), here we select E1 = 0.2. Now we see that yield (Y) of the system (6) increases and again decreases with increasing E2 , and the slope of the graph switches. The yield maximizes at E2 = 0.413 and it reaches to maximum sustainable yield (MSYE2 ). 4. Bionomic equilibrium Economic equilibrium is a state where economic forces such as supply and demand are balanced and in the absence of external influences the values of economic variables will not change. The biological equilibrium and the economic equilibrium both vary with the age-selection strategy. Thus synchronization of biological and economic equilibrium which is better description of economically profitable and stable system is visualized by bionomic equilibrium. Bionomic equilibrium gives us constraints to have profitable harvesting i.e., in which condition we can harvest either one of the population or all population together [31,17,18], so that we can get an economically profitable as well as biologically beneficial system. Age-selective harvesting of the population of any system controls its bionomic equilibrium when the biological equilibrium depends upon the selective harvesting [15,23,48]. To discuss the bionomic equilibrium of the system (6), we consider the parameters such as:

c1 = cost for unit effort for prey species,

p2 = price for unit biomass of predator. The net economic rent or net revenue R at any time is given by



p 1 μx − c1 E1 μeμτ1 + δ (eμτ1 − 1 )x

E 1 μx = 0, μeμτ1 + δ (eμτ1 − 1 )x

(10)

where R1 and R2 represent the net revenues for the prey and predator respectively. The bionomic equilibrium (x∞ , y∞ , E1∞ , E2∞ )

(11)

cβ xy − γ y − E2 e−γ τ2 y = 0, 1 + ax



(12)



p 1 μx

μeμτ1 + δ (eμτ1 − 1 )x

− c1 E1 + ( p2 e−γ τ2 y − c2 )E2 = 0.

(13)

To find the bionomic equilibrium, we consider the following cases: p μx

Case A. If c1 > μeμτ1 +δ1(eμτ1 −1 )x and c2 > p2 e−γ τ2 y i.e., the costs are greater than the revenue for the prey and predator both, then harvesting of prey and predator is not practicable, then the prey and predator harvesting is stopped (Ei = 0, i = 1, 2) (red region in Fig. 5). p μx

Case B. If c1 > μeμτ1 +δ1(eμτ1 −1 )x i.e., the cost is greater than the revenue for the prey, then harvesting of prey is not practicable and the prey harvesting is stopped (E1 = 0). Then only predator harvesting remains operational (i.e., c2 < p2 e−γ τ2 y). Therefore E1 = 0 c and c2 < p2 e−γ τ2 y, we have y∞ = p2 eγ τ2 (black region in Fig. 5). cβ xy 1+ax

−γy =

Case C. If c2 > p2 e−γ τ2 y i.e., the cost is greater than the revenue for the predator, then harvesting of predator is not practicable and the predator harvesting is stopped (E2 = 0). Then only prey harp μx vesting remains operational (i.e., c1 < μeμτ1 +δ1(eμτ1 −1 )x ). Therefore p μx

+( p2 e−γ τ2 y − c2 )E2 = R1 + R2

1 + ax

k

2

p1 = price for unit biomass of prey,



  x β xy αx 1 − − −

The point (x∞ , E2∞ ) will be any point on the line E2 e−γ τ2 y in the first quadrant of the xE2 -plane.

c2 = cost for unit effort for predator species,

R=

is given by the following simultaneous equations

c μeμτ1

E2 = 0 and c1 < μeμτ1 +δ1(eμτ1 −1 )x , we have x∞ = p μ−δ1 (eμτ1 −1 )c 1 1 (blue region in Fig. 5). The point (y∞ , E1∞ ) will be any point on



the line α x 1 − of the yE1 -plane.

x k



E 1 μx β xy − 1+ ax = μeμτ1 +δ (eμτ1 −1 )x in the first quadrant

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

117

Fig. 6. Variation of the quantities: (i) x∞ , (ii) y∞ , (iii) E1∞ , (iv) E2∞ of the Bionomic equilibrium with the variation of age selections of both prey and predator (τ1 − τ2 ). Parameters are: c1 = 2.5, c2 = 3.5, p1 = 2.5, p2 = 3, E1 = 0.5 and E2 = 0.4. Other parameters are as in the text. p μx

Case D. If c1 < μeμτ1 +δ1(eμτ1 −1 )x and c2 < p2 e−γ τ2 y, then the revenues for both the species being positive so that the whole fishery will be in operation. In this case x∞ =

c1 μeμτ1 p1 μ−δ (eμτ1 −1 )c1

and

y∞ = p2 eγ τ2 (green region in Fig. 5). Substituting the values of x∞ 2 and y∞ in Eqs. (11) and (12) we get c

E1∞ =

1

μ

   x β y∞  (μeμτ1 + δ (eμτ1 − 1 )x∞ ) α 1 − ∞ −

and

E2∞ = eγ τ2

k

 cβ x ∞ 1 + ax∞

1 + ax∞



−γ .

Now, we observe the dynamical feature of net economic rent or net revenue R at any time at a constant harvesting effort of both population (E1 = 0.4, E2 = 0.5), if population abundance be varied (0 < x < 2, 0 < y < 3), we observe the variation of economic rent (Fig. 5). When 0 < x < x∞ and 0 < y < y∞ , then the costs are greater than the revenue for the prey and predator both, then harvesting of prey and predator is not practicable. Also, the prey and predator harvesting is stopped (Case A: red region of Fig. 5), when 0 < x < x∞ and y > y∞ , the cost is greater than the revenue for the prey, then harvesting of prey is not practicable. The prey harvesting is stopped and then only predator harvesting remains operational (Case B: black region of Fig. 5), when x > x∞ and 0 < y < y∞ , the cost is greater than the revenue for the predator, then harvesting of predator is not practicable. Also, the predator harvesting is stopped and then only prey harvesting remains operational (Case C: blue region of Fig. 5) and finally when x > x∞ and y > y∞ , then both prey and predator harvesting are operational (Case D: green region of Fig. 5). It is the most important case (Case D) because economically it is more significant than Case

B and C [31,17,18]. Now, we observe the variation of the quantities of bionomic equilibrium with the variation of age selection strategies (τ 1 , τ 2 ) [15,23,48] of prey and predator respectively (Fig. 6(i): x∞ , Fig. 6(ii): y∞ , Fig. 6(iii): E1∞ , Fig. 6(iv): E2∞ ). From the figure it is clear that x∞ and E1∞ gradually increases with increasing τ 1 , whereas they do not have any variation with the variation of τ 2 . But y∞ and E2∞ can show just opposite characteristics with the same variation of τ 1 and τ 2 . It is clear that the opposite variational modes of (x∞ − y∞ ) and (E1∞ − E2∞ ) clarify the stable bionomic equilibrium of the system. 5. Optimal harvesting policy Optimization is a technique which maximizes the profit without any harmful effect on the system. Optimal harvesting policy gives us a strategy which helps us in both maximizing the profit by minimum level of effort as well as eradicating the risk of extinction of the harvested population [27,14]. It also depends on age-selection strategy [16,19,13] of the harvested population because if we harvest the population at an early age then there is an extinction risk of the population whereas if we harvest the population at much older age then it is not economically profitable, so this policy ensures us the sensitive age class which neither has extinction risk nor has profit dampening effect. In this section, our objective is to maximize the objective functional form of the harvesting model, with the instantaneous annual rate of discount ξ is as follows:

J ( E1 , E2 ) =



∞ 0

e −ξ t



p 1 μx

μeμτ1 + δ (eμτ1 − 1 )x



− c1 E1



+( p2 e−γ τ2 y − c2 )E2 dt

(14)

118

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

subject to the state constraints (6) and the control constraints

0 ≤ Ei (t ) ≤ Ei max ;

i = 1, 2

(15)

Hamiltonian for the problem is given by



Using equilibrium conditions, Eq. (23) becomes





p 1 μx − c1 E1 + ( p2 e−γ τ2 y − c2 )E2 μeμτ1 + δ (eμτ1 − 1 )x    x E 1 μx β xy +λ1 α x 1 − − − k 1 + ax μeμτ1 + δ (eμτ1 − 1 )x  cβ xy +λ2 − γ y − E2 e−γ τ2 y , (16) 1 + ax

H = e −ξ t

where λ1 and λ2 are adjoint variables. We have

  ∂H p 1 μx λ1 μx = e −ξ t − c1 − μτ μτ ∂ E1 μe 1 + δ ( e 1 − 1 ) x μeμτ1 + δ (eμτ1 − 1 )x = σ1 (t ), (17) ∂H = e−ξ t ( p2 e−γ τ2 y − c2 ) − λ2 e−γ τ2 y ∂ E2 = σ2 (t ).

(18)

The optimal control Ei (t ) (i = 1, 2 ) must clearly satisfy the conditions



Ei (t ) =

Eimax 0

if if

σi (t ) > 0 σi (t ) < 0.

Since σ i (t) causes Ei (t ) (i = 1, 2 ) to switch between level 0 and Ei max , σi (t ) (i = 1, 2 ) are called switching function. Depending on the sign of the switching function σ i (t), the optimal control Ei (t) is a bang-bang switching from one extreme point to other. When σi (t ) = 0 (i = 1, 2 ), the Hamiltonian function H becomes independent of the control variable Ei (t ) (i = 1, 2 ) and the optimal control can not be determined by the above procedure. It is then called a singular control Ei ∗ (t), 0 < Ei ∗ (t) < Ei max . Hence the optimal harvesting policy is



Ei (t ) =

Eimax 0 Ei ∗

if if if

λ1 = e−ξ t

(i = 1, 2 ). From Eqs. (17) and

 c1 (μeμτ1 + δ (eμτ1 − 1 )x ) p1 − , μx

 λ2 = e−ξ t p2 −

c2 e−γ τ2 y



.

(20)

By Pontryagain’ s maximum principle [40], the adjoint equations are

dλ1 =− dt

∂H , ∂x

dλ2 =− dt

∂H . ∂y

(21)

(22)

From Eqs. (16) and (21) we have



dλ1 p1 μ2 eμτ1 = −e−ξ t E1 μτ 1 dt (μe + δ (eμτ1 − 1 )x )2 −λ1 −λ2





 2α x E1 μ2 eμτ1 βy α− − − 2 μτ μτ 2 k (1 + ax ) ( μe 1 + δ ( e 1 − 1 ) x )  cβ y (1 + ax )2

.

(23)



dλ1 p1 μ2 eμτ1 = −e−ξ t E1 dt (μeμτ1 + δ (eμτ1 − 1 )x )2   α x c1 (μeμτ1 + δ (eμτ1 − 1 )x ) β axy −e−ξ t p1 − − + μx k (1 + ax )2 E1 μδ (eμτ1 − 1 )x + μτ μτ 2 1 1 ( μe + δ ( e − 1 ) x )   c β y c2 −e−ξ t p2 − −γ τ . (25) e 2 y (1 + ax )2 Similarly we get





dλ2 λ1 β x cβ x = −e−ξ t E2 p2 e−γ τ2 + − λ2 − γ − E2 e−γ τ2 dt 1 + ax 1 + ax    c1 (μeμτ1 + δ (eμτ1 − 1 )x ) = −e−ξ t E2 p2 e−γ τ2 − p1 − μx  β x  × . (26) 1 + ax Integrating Eqs. (25) and (26) we get

λ1 =

λ2 =

 



p1 μ2 eμτ1 μτ 1 ξ (μe + δ (eμτ1 − 1 )x )2   α x c1 (μeμτ1 + δ (eμτ1 − 1 )x ) β axy + p1 − − + μx k (1 + ax )2    cβ y  c2 E1 μδ (eμτ1 − 1 )x + + p 2 − −γ τ , μτ μτ 2 e 2y ( μe 1 + δ ( e 1 − 1 ) x ) (1 + ax )2 (27) 1

e −ξ t E 1

1

e−ξ t E2 p2 e−γ τ2 − p1 −

ξ ×

(19)



Now putting the values of λ1 and λ2 from Eqs. (19) and (20) respectively in Eq. (24), we get

σi (t ) > 0 σi (t ) < 0 σi (t ) = 0

for i = 1, 2. For singular control σi (t ) = 0 (18) we get,





dλ1 p1 μ2 eμτ1 = −e−ξ t E1 dt (μeμτ1 + δ (eμτ1 − 1 )x )2  αx β axy E1 μδ (eμτ1 − 1 )x −λ1 − + + k (1 + ax )2 (μeμτ1 + δ (eμτ1 − 1 )x )2  cβ y −λ2 . (24) (1 + ax )2





 β x  1 + ax

c1 (μeμτ1 + δ (eμτ1 − 1 )x ) μx

.

 (28)

The monetary values assigned to the cost which are difficult to calculate or currently unknowable, referred as shadow price. The origin of these costs is typically due to an unwillingness to recalculate a system to account for marginal production or an externalization of costs. So, we neglect the constant of integration in order to shadow prices λi eξ t (i = 1, 2 ) of the two species are bounded. From Eqs. (19) and (27) we get,



e −ξ t p 1 − =

c1 (μeμτ1 + δ (eμτ1 − 1 )x ) μx

 





p1 μ2 eμτ1 ξ (μeμτ1 + δ (eμτ1 − 1 )x )2   α x c1 (μeμτ1 + δ (eμτ1 − 1 )x ) β axy + p1 − − + μx k (1 + ax )2    cβ y  c2 E1 μδ (eμτ1 − 1 )x + + p − , 2 e−γ τ2 y (μeμτ1 + δ (eμτ1 − 1 )x )2 (1 + ax )2 1

e −ξ t E 1

 μ( p μeμτ1 + xθ δ (eμτ1 − 1 )) 1 1 (μeμτ1 + xδ (eμτ1 − 1 )))2  (α x + kξ )(1 + ax )2 − aβ xy   cβ y  = θ1 − θ2 2 k(1 + ax ) (1 + ax )2 ⇒ E1

(29)

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

119

Fig. 7. Optimal harvesting effort of (i) prey and (ii) predator respectively with the variation of x and y with constant age selection τ1 = 0.5 and τ2 = 0.1 along with annual discount rate ξ = 1.5. Parameters are as in the previous figure.

where θ1 = p1 −

c1 (μeμτ1 +δ (eμτ1 −1 )x ) μx

and θ2 = p2 −

c2 . e −γ τ 2 y

Simi-

larly, from Eqs. (20) and (28) we get

e

−ξ t

=



1

ξ

p2 −





c2 e−γ τ2 y



e−ξ t E2 p2 e−γ τ2 − p1 −

c1 (μeμτ1 + δ (eμτ1 −1 )x ) μx

 β x  1 + ax

,

⇒ E2 p2 e−γ τ2     β x  c2 c1 (μeμτ1 + δ (eμτ1 −1 )x ) = ξ p 2 − −γ τ + p1 − . e 2y μx 1 + ax (30) Eqs. (29) and (30) give the optimal harvesting efforts as

E1Opt =

 μeμτ1 + xδ (eμτ1 − 1 ) 2 1 + ax

 θ (α x + kξ )(1 + ax )2 − β y(ax + ckθ )  1 2 × kμ( p1 μeμτ1 + xθ1 δ (eμτ1 − 1 )) and

E2Opt =

1 p2 e−γ τ2



+ p1 −

  ξ p2 −

c2 e−γ τ2 y



 c1 (μeμτ1 + δ (eμτ1 − 1 )x ) μx

(31)

β x  1 + ax

.

(32)

Hence solving steady state equations together with (31) and (32) we obtain the optimal solution (xξ , yξ ) and optimal harvesting efforts E1ξ and E2ξ . Now we are going through an important observation regarding the optimal harvesting policy of the system (6). At a constant age selection τ1 = 0.4, τ2 = 0.5 of prey and predator with annual discount rate ξ = 1.5, we obtain an opposite change of optimal harOpt Opt vesting policy of prey (E1 ) and predator (E2 ) with the variation of population abundance 0 < x < 6, 1 < y < 2 (Fig. 7). Whereas if we treat the same system with a constant discount rate ξ = 1.5 and apply this set up on a constant number of prey (4) and predator (3), we observe that optimal harvesting of both population increases with increasing age selection strategy of prey (0 < τ 1 < 6), but there is almost no variation of optimal harvesting of both population with the variation of predator age selectivity (1 < τ 2 < 2) (Fig. 8). Therefore, it is observed that optimal harvesting policy is size-specific as well as density dependent [4].

6. Conclusions and discussions Species interactions are not static but they are state-dependent (i.e. nonlinear) and change as ecosystem factors shift to new states [11]: e.g. fish populations show sensitivity to oceanographic conditions that increases when population decline [1]. The condition of maximum profit arrives before the MSY where the value of the landed species is dependent on the number of individuals landed. The situation of maximum profit is equivalent as the MSY when the value of the landed species is independent of harvest. The economic effects on communities dependent on fisheries are considered by social objectives [9,26,24]. In this paper, we analyzed the interaction between prey and predator fish population with age-selection harvesting. We carried out numerical simulations to clarify the effect of prey and predator age selection. First, we observed that the model system remains stable for either τ 1 ∈ (0, 1.2) or τ 2 ∈ (0, 2.91). As we increase the prey age-selection, we should decrease predator age-selection or its vice-versa to maintain the stability of the model system. If we select both the age-selection harvesting more than the above defined range simultaneously, both the population start to oscillates around its stable equilibrium point. For a fixed predator age-selection harvesting, we found that the MSY increases gradually with the increase of prey harvesting effort till the predator harvesting reaches upto 0.413. It is also observed that the bionomic equilibrium is directly affected by age-selections. (x∞ , E1 ∞ ) depends on prey age-selection whereas (y∞ , E2 ∞ ) decreases on predator age-selection. It is interesting to note that the optimal harvesting of both the population increases with the increase of prey age-selection, i.e., optimal harvesting increases independent of predator age-selection. Most emphasize on that the model formulation is new and results i.e. MSY, bionomic equilibrium and optimal harvesting policy are highly depending upon the age-selection of harvested population like the stability of biological equilibrium and these results point out that under particular selective harvesting policy, we can prevent extinction risk of population and biodiversity. Therefore, our economic purposes will be solved. Appendix A Case 1: Instantaneous harvesting In case of instantaneous harvesting (i.e., τ1 = 0 = τ2 ), the Eq. (8) becomes

120

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

Fig. 8. Optimal harvesting effort of prey (i) and predator (ii) respectively with the variation of τ 1 and τ 2 with constant population density x = 4 and y = 3 along with annual discount rate ξ = 1.5. Parameters are as in the previous figure.

λ2 + (A + C1 + C2 )λ + (B + D1 + D2 + E ) = 0.

All roots of the Eq. (A1) have negative real parts if and only if

(H1 ) A + C1 + C2 > 0 and B + D1 + D2 + E > 0. Therefore, the equilibrium point E∗ (x∗ , y∗ ) is locally asymptotically β y∗ ( γ +E ) stable when (H1 ) holds. Note that B + D1 + D2 + E = (1+ax∗ )22 is al-



aβ y∗

ways positive and A + C1 + C2 = x∗ αk − (1+ax∗ )2 if E2 > caβ



ak(α −E1 )−α ak(α −E1 )+α





− γ and E1 ≤ α 1 −



will be positive

cβ +aγ ak(cβ −aγ )



. The results

are summarized in the following theorem. E∗ (x∗ ,

y∗ )

Theorem 2.1. The interior equilibrium of the system (6) exists and becomes locally asymptotically stable in absence of delays if the following conditions hold:



γ

aγ (ak+1 ) ak−1

Let

(A1)

(i ) cβ > max aγ + k , , (ii ) ak > 1 , (α−E1 )−α (iii ) caβ ak − γ < E2 < ac+βγ − γ and ak(α −E1 )+α k

γ +E2 (iv ) E1 < α 1 − k(cβ −a (γ +E2 )) .

If Theorem 2.1 and (H3 ) hold then the Eq. (A4) has a unique positive root ω02 . Substituting ω02 into (A3), we have

τ2 n =

1

ω0



cos−1

 (D2 + E )(ω02 − B − D1 ) − (A + C1 )C2 ω02 2 nπ + , ω0 C22 ω02 + (D2 + E )2

n = 0, 1, 2, . . . .. Define the function θ (τ 2 ) ∈ [0, 2π ) such that cosθ (τ 2 ) is given by the right-hand sides of the last equation. Then the τ 2 we seek, at which stability switches occur, are the solutions of

S n ( τ2 ) = τ2 − τ2 n

(A5)

(H4 ) C22 −(A + C1 )2 + 2(B + D1 )>0, (B + D1 )2 −(D2 + E )2 >0 and [C22 − (A + C1 )2 + 2(B + D1 )]2 > 4[(B + D1 )2 − (D2 + E )2 ].

If prey is harvested instantaneously and predator is harvested selectively then τ1 = 0 and τ 2 = 0. In this case, the characteristic Eq. (8) becomes

(A2)

Let iω(ω > 0) be a root of the Eq. (A2). Then we have,

(D2 + E ) cos ωτ2 + C2 ω sin ωτ2 = ω2 − (B + D1 ), C2 ω cos ωτ2 − (D2 + E ) sin ωτ2 = −(A + C1 )ω.

( B + D1 ) 2 − ( D2 + E ) 2 < 0 .

Let,

Case 2: Selective harvesting of predator only

λ2 + (A + C1 )λ + B + D1 + (C2 λ + D2 + E )e−λτ2 = 0.

(H3 )

(A3)

2 If Theorem 2.1 and (H4 ) hold then (A4) has two positive roots ω+ 2. and ω− Substituting ω± into (A3), we get

τ = ± 2k

ω±



cos

−1

 (D2 + E )(ω±2 − B − D1 ) − (A + C1 )C2 ω±2 2kπ + , 2 + ( D + E )2 ω± C22 ω± 2

k = 0, 1, 2, . . . .. If λ(τ 2 ) be the root of (A2) satisfying Reλ(τ2n ) = 0 (respectively, Reλ(τ2± ) = 0) and Imλ(τ2n ) = ω0 (respectively, Imλ(τ2± ) = ω± ), we k k get



This leads to

1

ω4 − [C22 − (A + C1 )2 + 2(B + D1 )]ω2 + (B + D1 )2 − (D2 + E )2 = 0.

d Re(λ ) d τ2

(A4)



= τ2 =τ20 ,ω=ω0

It follows that the Eq. (A4) has no positive roots if the following conditions are satisfied:

(H2 ) (A + C1 )2 −C22 −2(B + D1 )>0 and (B + D1 )2 −(D2 + E )2 >0. Hence, all roots of the Eq. (A4) will have negative real parts when τ 2 ∈ [0, ∞) if conditions of the Theorem 2.1 and (H2 ) are satisfied.

ω 4 + ( D2 + E ) 2 − ( B + D1 ) 2 ω2 (C22 ω2 + (D2 + E )2 ) >

Similarly, we can show that



d Re(λ ) d τ2



> 0, τ2 =τ2+ ,ω=ω+ k

ω2

C22 ω2 + (D2 + E )2



d Re(λ ) d τ2

> 0, (by H3 ).

 < 0. τ2 =τ2− ,ω=ω− k

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

From Corollary (2.4) in Ruan and Wei (2003), we have the following conclusions. Theorem 2.2. Assume τ1 = 0, τ2 = 0 and conditions of the Theorem 2.1 are satisfied. Then the following conclusions hold: (i) If H2 holds then the equilibrium E∗ (x∗ , y∗ ) is asymptotically stable for all τ 2 ≥ 0. (ii) If H3 holds then the equilibrium E∗ (x∗ , y∗ ) is conditionally stable. It is locally asymptotically stable for τ2 < τ20 and unstable for τ2 > τ20 . Furthermore, the system (6) undergoes a Hopf bifurcation at E∗ (x∗ , y∗ ) when τ2 = τ20 , where

τ2 0 =

1

ω0



cos−1

 (D2 + E )(ω02 − B − D1 ) − (A + C1 )C2 ω02 . C22 ω02 + (D2 + E )2 0

m−1

0

1

, τ2+ ) and unstable when τ2 ∈ [τ2+ , τ2− ) ∪ (τ2+ , τ2− ) ∪ m

. . . .. ∪ (τ2+

0

1

1

) ∪ (τ2+m , ∞ ). Furthermore, the system undergoes a Hopf bifurcation at E∗ (x∗ , y∗ ) when τ2 = τ2+ , k = m−1

, τ2−

0

m−1

k

0, 1, 2, . . . ..

In the following, we assume that



(H6 )

d (Reλ ) d τ1



= 0. λ=iω0

Therefore, by the general Hopf bifurcation theorem, we have the following result on the stability and bifurcation of the system (6). Theorem 2.3. For the system (6), suppose parameters satisfy conditions of the Theorem 2.1, H3 and H5 and τ2 ∈ [0, τ20 ). Then the equilibrium E∗ (x∗ , y∗ ) is asymptotically stable when τ1 ∈ (0, τ10 ), unstable when τ1 > τ10 and a Hopf bifurcation occurs when τ1 = τ10 . Case 4: Selective harvesting of prey only

(iii) If H4 holds then there is a positive integer m such that the equilibrium is stable when τ2 ∈ [0, τ2+ ) ∪ (τ2− , τ2+ ) ∪ . . . .. ∪

(τ2−

121

Case 3: Selective harvesting of both the prey and predator If both the prey and predator are harvested selectively then τ 1 = 0 and τ 2 = 0. We consider Eq. (8) with τ 2 in its stable interval (τ2 ∈ [0, τ20 )) and regard τ 1 as a parameter (τ 1 ∈ (0, ∞)). Without loss of generality, we consider that the system parameters satisfy Theorem 2.1 and (H3 ). Let iω(ω > 0) be a root of Eq. (8) and we obtain

ω4 + Aω2 + B2 + D22 − D21 − E 2 + 2Bsin ωτ2 + 2Ccos ωτ2 = 0, (A6) where

 = A2 + C 2 − 2B − C 2 , A 2 1

Next, we assume that predator is harvested instantaneously but prey is harvested selectively so that τ 1 = 0 but τ2 = 0. The proof is similar to the Case 2 and we only summarize the results below. Theorem 2.4. Assume τ1 = 0, τ2 = 0 and Theorem 2.1 is satisfied. Then the following conclusions hold: (i) If (H7 ) (A + C2 )2 − C12 − 2(B + D2 ) > 0, (B + D2 )2 − (D1 + E )2 > 0 hold, then the equilibrium E∗ (x∗ , y∗ ) is locally asymptotically stable for all τ 1 ≥ 0. (ii) If (H8 ) (B + D2 )2 − (D1 + E )2 < 0 holds, then the equilibrium E∗ (x∗ , y∗ ) is locally asymptotically stable for τ1 < τ˜10 and unstable for τ1 > τ˜10 . Furthermore, the system (6) undergoes a Hopf bifurcation at E∗ (x∗ , y∗ ) when τ1 = τ˜10 , where τ˜10 = min{τ˜1 p | p = 0, 1, 2, . . .} and

1

τ˜1 p =

ω0

cos−1

  (D1 + E )(ω20 − B − D2 ) − (A + C2 )C1 ω20 2 pπ + , ω0 C12 ω20 + (D1 + E )2

p = 0, 1, 2, . . . .. (iii) If (H9 ) C12 − (A + C2 )2 + 2(B + D2 ) > 0, (B + D2 )2 − (D1 + E )2 > 0 and [C12 − (A + C2 )2 + 2(B + D2 )]2 > 4[(B + D2 )2 − (D1 + E )2 ] hold, then there is a positive integer l such that the equilibrium is stable when τ1 ∈ [0, τ˜1+ ) ∪ (τ˜1− , τ˜1+ ) ∪ . . . .. ∪ (τ˜1− , τ˜1+ ), and unsta0

0

1

l−1

l

ble when τ1 ∈ [τ˜1+ , τ˜1− ) ∪ (τ˜1+ , τ˜1− ) ∪ . . . .. ∪ (τ˜1+ , τ˜1− ) ∪ (τ˜1+ , ∞ ). 0

0

1

1

l−1

l−1

l

Furthermore, the system undergoes a Hopf bifurcation at (x∗ , y∗ ) when τ1 = τ˜1+ , p = 0, 1, 2, . . . ..

 = ωC1 E − ω3C2 − ωAD2 + ωBC2 and B

p

 = −D1 E − ω2 D2 + BD2 + ω2 AC2 . C

Case 5: Selective harvesting of both the prey and predator

We define,

Finally, we assume that both the prey and predator are harvested selectively then τ 1 = 0 and τ 2 = 0. We consider Eq. (8) with τ 1 in its stable interval and regard τ 2 as a parameter. The proof is similar to the Case 3 and we only summarize the results below.

F (ω ) = ω

4

ω2 + B2 + D2 − D2 − E 2 + 2B sin ωτ2 + 2C cos ωτ2 +A 2 1

and assume that

(H5 )

( B + D2 ) 2 − ( D1 + E ) 2 < 0 .

Then it is easy to check that F(0) < 0 and F (∞ ) = ∞. Thus, (A6) has finite positive roots ω1 , ω2 , . . . .., ωk . For every fixed ωi , i = j 1, 2, . . . ., k, there exists a sequence {τ1 | j = 1, 2, . . . .}, where i

1  M 2 jπ 1 τ = cos−1 + , i = 1, 2, . . . ., k; j = 1, 2, . . . . ωi N1 ωi

Theorem 2.5. For the system (6), suppose parameters satisfy conditions of the Theorem 2.1, H3 and H5 and τ1 ∈ [0, τ˜10 ). Then the equilibrium E∗ (x∗ , y∗ ) is locally asymptotically stable when τ2 ∈ (0, τ 20 ), unstable when τ2 > τ 20 and a Hopf bifurcation occurs when τ2 =

τ 20 , where τ 20 = min{τ t2s |s = 1, 2, . . . , q; t = 0, 1, 2, . . .},

j 1i

τ t2s =

where,

M1 = P1 S1 + P2 T1 + (Q1 S1 + R1 T1 )cosωi τ2 + (R1 S1 − Q1 T1 )sinωi τ2 , N1 = S12 + T12 , P1 = −ωi2 + B, P2 = Aωi , Q1 = D2 , R1 = C2 ωi , S1 = −(E cosωi τ2 + D1 ), T1 = E sinωi τ2 − C1 ωi , i = 1, 2, . . . ., k Let τ10 = min{τ1 | i = 1, 2, . . . ., k; j = 1, 2, . . . .}. Define the function j

1  M 2t π 2 cos−1 + , s = 1, 2, . . . ., q; t = 1, 2, . . . . ωs N2 ωs

where,

M2 = P1 S2 + P2 T2 + (Q2 S2 + R2 T2 )cosωs τ1 + (R2 S2 − Q2 T2 )sinωs τ1 , N2 = S22 + T22 , P1 = −ωs2 + B, P2 = Aωs , Q2 = D1 , R2 = C1 ωs , S2 = −(E cosωs τ1 + D2 ), T2 = E sinωs τ1 − C2 ωs , s = 1, 2, . . . ., q.

i

θ (τ 1 ) ∈ [0, 2π ) such that cosθ (τ 1 ) is given by the right-hand sides of the last equation. Then the τ 1 we seek, at which stability switches occur, are the solutions of

Tn (τ1 ) = τ1 − τ1n

(A7)

References [1] Anderson CNK, Hsieh CH, Sandin SA, Hewitt R, Hollowed A, Beddington J, May RM, Sugihara G, et al. Why fishing magnifies fluctuations in fish abundance. Nature 2008;452:835–9.

122

D. Jana et al. / Chaos, Solitons and Fractals 93 (2016) 111–122

[2] Arino J, Wang L, Wolkowicz GSK. An alternative formulation for a delayed logistic equation. J Theo Biol 2006;241:109–19. [3] Banerjee TK, Hazra US, Banerjee R. Studies on the incidence behavior and morphometry of hilsa fish, Tenualosa ilisha in relation to environmental attributes in the selected portion of the upper stretch of hooghly estuary. Spring 2013;1(1):226–47. [4] Botsford LW. Optimal fishery policy for size-specific, density-dependent population models. J Math Biol 1981;12(3):265–93. [5] Brauer F. Stability of some population models with delay. Math Biosci 1977;33:345–58. [6] Brauer F, Soudack AC. Coexistence properties of some predator-prey systems under constant rate harvesting and stocking. J Math Biol 1981;12:101–14. [7] Calder WA. Size, Function, and Life History. Cambridge, Massachusetts: Harvard University Press; 1984. [8] Chaudhuri KS, Saha Roy S. On the combined harvesting of a prey-predator system. J Biol Sys 1996;4(3):373–89. [9] Clark CW. Mathematical Bioeconomics: The optimal Management of Renewable Resources. New York: John Wiley & Sons; 1976. [10] Conover DO, Munch SB. Sustaining fisheries yields over evolutionary time scales. Science 2002;297:94–6. [11] Deyle ER, May RM, Munch SB, Sugihara G. Tracking and forecasting ecosystem interactions in real time. Proc R Soc B 2016;283:2015–258. [12] De DK, Datta NC. Age, growth, length-weight relationship and relative condition in hilsa, Tenualosa ilisha (hamilton) from the hooghly estuarine system. Indian J Fish 1990;37(3):199–209. [13] Diekert FK, Hjermann DØ, Nævdal E, Stenseth NC. Spare the young fish: Optimal harvesting policies for north-east arctic cod. Environ Resource Econ 2010;47:455–75. [14] Dubey B, Chandra P, Sinha P. A resource dependent fishery model with optimal harvesting policy. J Biol Sys 2002;10(1):1–13. [15] Fenberg PB, Roy K. Ecological and evolutionary consequences of size-selective harvesting: how much do we know? Molecular Ecol 2008;17:209–20. [16] Frederick SW, Peterman RM. Choosing fisheries harvest policies: when does uncertainty matter? Can J Fish Aquat Sci 1995;52:291–306. [17] Gardner C., Larkin S., Seijo JC.. Systems to maximize economic benefits in lobster fisheries. 2013. http://www.researchgate.net/publication/255787730. [18] Harris JM., Codur AM. Economics of fisheries. 2014. The Encyclopedia of Earth. [19] Hritonenko N, Yatsenko Y. Optimization of harvesting age in an integral age-dependent model of population dynamics. Math Biosci 2005;195:154–67. [20] Iversen ES. Living Marine Resources: Their Utilization and Management. New York, NY: Chapman and Hall; 1996. [21] Jana D, Pathak R, Agarwal M. On the stability and hopf bifurcation of a prey– generalist predator system with independent age-selective harvesting. Chaos Solit & Fract 2016;83:252–73. [22] Jhingran VG, Natarajan AV. Derivation of average lengths of different age groups in fishes. Fish Res Bd Canada 1969;26(11):3037–76. [23] Jørgensen C, Ernande B, Fiksen Ø. Size-selective fishing gear and life history evolution in the northeast arctic cod. Evolutionary Appl 2009;2(3):356–70. [24] Kar TK. Selective harvesting in a prey-predator fishery with time delay. Math Comp Model 2003;38:449–58. [25] Kar TK, Chaudhuri KS. On non-selective harvesting of two competing fish species in the presence of toxicity. Ecol Model 2003;161:125–37. [26] Kar TK, Matsuda H. A bioeconomic model of a single-species fishery with a marine reserve. J Environ Manag 2008;86:171–80. [27] Kim C. Optimal management of multi-species north sea fishery resources. Weltwirtschaftliches Archiv 1983;119(1):138–51. [28] Krˇivan V, Jana D. Effects of animal dispersal on harvesting with protected areas. J Theo Biol 2015;364:131–8. [29] Kuang Y, Freedman HI. Uniqueness of limit cycles in gause-type models of predator-prey systems. Math Biosci 1988;88:67–84. [30] Law R, Grey DR. Evolution of yields from populations with age-specific cropping. Evol Ecol 1989;3:343–59. [31] Lian C, Singh R, Weninger Q. Fleet restructuring, rent generation, and the design of individual fishing quota programs: Empirical evidence from the pacific coast groundfish fishery. Marine Resour Econ 2009;24(4):329–59.

[32] Martin A, Ruan S. Predator-prey models with delay and pray harvesting. J Math Biol 2001;43:247–67. [33] Mesterton-Gibbons M. A technique for finding optimal two-species harvesting policies. Nat Resource Model 1996;92:235–44. [34] Netessine S, Shumsky R. Introduction to the theory and practice of yield management. INFORMS Trans Education 2002;3(1). [35] Neubert MG, Herrera GE. Triple benefits from spatial resource management. Theo Ecol 2008;1:5–12. [36] Peters RH. The Ecological Implications of Body Size. Cambridge, UK: Cambridge University Press; 1983. [37] Pikitch EK, Santora C, Babcock EA, Bakun A, Bonfil R, Conover DO, Dayton P, Doukakis P, Fluharty D, Heneman B, Houde ED, Link J, Livingston PA, Mangel M, McAllister MK, Pope J, Sainsbury KJ. Ecosystem-based fishery management. Science 2004;305:346–7. [38] Pillay SR, Rao VK. Observations on the biology and fishery of the hilsa, Hflsa ilisha (hamilton) of river godavari. Proc Indo-Pacific Fish Coun 1958;10(2):37–61. [39] Pillay SR, Rao VK. Observations on the biology and fishery of the hilsa, Hflsa ilisha (hamilton) of river godavari. Proc Indo-Pacific Fish Coun 1962;10(2):37–61. [40] Pontryagin LS, Boltyonsku VG, Gamkrelidre RV, Mishchenko EF. The Math. Theo Optimal Process. New York: Wiley; 1962. [41] Pope JG, Rice JC, Daan N, Jennings S, Gislason H. Modelling an exploited marine fish community with parameters results from a simple size-based model. ICES J Marine Sci 2006;63:1029–44. [42] Quddus MMA, Shimizu M, Nose Y. Comparison of age and growth of two types of hilsa ilisha in bangladesh waters. Bull Jap Soc Sci Fish 1984;50(1):51–7. [43] Ramakrishnaiah M. Biology of Hilsa ilisha (hamilton) from the chilka lake with an account of its racial status. Indian J Fish 1972;19(1 & 2):35–53. [44] Rajyalakshmi T. The population characteristics of the godavari hilsa over the years 1963-67. Indian J Fish 1973;20(1):78–94. [45] Rothschild BJ, Keiley EF, Jiao Y. Failure to eliminate overfishing and attain optimum yield in the new england groundfish fishery. ICES J Marine Sci 2014;71(2):226–33. [46] Royce WF. Introduction to the Practice of Fishery Science. San Diego, CA: Academic Press; 1996. [47] Shepherd JG. A versatile new stock-recruitment relationship of fisheries and construction of sustainable yield curves. J du Conseil International pour l’Exploration de la Mer 1982;40:67–75. [48] Skonhoft A. The maximum economic yield management of an age-structured salmon population. laksedemografi bookchapter 2012. [49] Srinivasu PDN, Ismail S, Naidu CR. Global dynamics and controllability of a harvested prey-predator system. J Biol Sys 2001;9(1):67–79. [50] Svedang H, Hornborg S. Selective fishing induces density-dependent growth. Nature Commun 2014. doi:10.1038/ncomms5152. [51] Takeuchi Y. Global Dynamical Properties of Lotka-Volterra Systems. Singapore: World Scientific Publishing Company; 1996. [52] Walsh MR, Munch SB, Chiba S, Conover DO. Maladaptative changes in multiple traits caused by fishing: impediments to population recovery. Ecol Letter 2006;9:142–8. [53] Webster CR, Lorimer CG. Single-tree versus group selection in hemlock-hardwood forests: are smaller openings less productive? Can J Forest Res 2002;32(4):591–604. [54] Xia J, Liu Z, Yuan R, Ruan S. The effects of harvesting and the delay on predator-prey systems with holling type II functional response. SIAM J Appl Math 2009;70(4):1178–200. [55] Yender R, Michel J, Lord C. Managing Seafood Safety after an Oil Spill. National Oceanic and Atmospheric Administration, NOAA National Ocean Service, Office of Response and Restoration; 2002. [56] Zhou S, Smith ADM, Puntb AE, Richardson AJ, Gibbs M, Fulton EA, Pascoe S, Bulman C, Bayliss P, Sainsbury K. Ecosystem-based fisheries management requires a change to the selective fishing philosophy. Proc Nat Acad Sci USA 2010;107(21):9485–9.