Global dynamics and bifurcation in a stage structured prey–predator fishery model with harvesting

Applied Mathematics and Computation 218 (2012) 9271–9290

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Global dynamics and bifurcation in a stage structured prey–predator fishery model with harvesting Kunal Chakraborty a,⇑, Soovoojeet Jana b, T.K. Kar b a Information Services and Ocean Sciences Group, Indian National Centre for Ocean Information Services, Hyderabad ‘‘Ocean Valley’’, Pragathi Nagar (BO), Nizampet (SO), Hyderabad 500 090, India b Department of Mathematics, Bengal Engineering and Science University, Shibpur, Howrah 711 103, India

a r t i c l e

i n f o

Keywords: Prey–predator fishery Stage structure Harvesting Hopf bifurcation Optimal control Global stability

a b s t r a c t This paper describes a prey–predator model with stage structure for predator and selective harvesting effort on predator population. The Holling type II functional response function is taken into consideration. All the equilibria of the proposed system are determined and the behavior of the system is investigated near them. Local stability of the system is analyzed. Geometric approach is used to derive the sufficient conditions for global stability of the system. The occurrence of Hopf bifurcation of the model system in the neighborhood of the co-existing equilibrium point is shown through considering maximal relative increase of predation as bifurcation parameter. Fishing effort used to harvest predator population is considered as a control to develop a dynamic framework to investigate the optimal utilization of the resource, sustainability properties of the stock and the resource rent earned from the resource. Pontryagin’s maximum principle is used to characterize the optimal control. The optimal system is derived and then solved numerically using an iterative method with Runge–Kutta fourth-order scheme. Simulation results show that the optimal control scheme can achieve sustainable ecosystem.  2012 Elsevier Inc. All rights reserved.

1. Introduction The study of the dynamics of prey–predator systems is one of the dominant subjects in mathematical ecology due to its universal existence and importance. Thus predator–prey models have been in the focus of ecological science from the early days of this discipline. It has been turned out very soon that predator–prey systems can show different dynamical behaviors such as steady states, oscillations, bifurcations and chaos depending on the model parameters. There has been an increasing interest during last few decades to study predator–prey system in population dynamics and have studied by both mathematicians and ecologists [1–4]. On the other hand, the fundamental objective of most fishery management is to ensure conservation of fishery resource into the future and to provide a sustainable flow of benefits to human society. In order to control the stock, catches and fishing effort of the fishery, it is possible to derive some useful management procedure that can be imposed to the fishery to get protection from overexploitation. Therefore, it requires some strong scientific research in the field of fisheries biology and economics. Economic aspects of a resource can be easily examined with the help of a suitable model using the empirical data of the resource [5,6]. The dynamical analysis of the predator– prey model plays an important role in mathematical biology. Though many biologists believe that if unique positive equilibrium point of a predator–prey system is locally asymptotically stable, then it is globally asymptotically stable, ⇑ Corresponding author. E-mail addresses: [email protected] (K. Chakraborty), [email protected] (S. Jana), [email protected] (T.K. Kar). 0096-3003/$ - see front matter  2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.03.005

9272

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

however it is not always true. It is proved that a unique positive locally asymptotically stable equilibrium point has at least one limit cycle surrounding the equilibrium point under suitable condition. Thus many mathematicians try to use some well known methods to find the conditions for global stability for the equilibrium point of predator–prey system [7– 10]. The optimal management of renewable resources as fishery, which has a direct relationship to sustainable development, has been studied extensively by Clark [11]. The subject of harvesting in predator–prey systems is considered as a multidisciplinary area of research which includes economists, ecologists and natural resource managers. Most research has focused attention on optimal exploitation, guided entirely by profits from harvesting. Clark [12], Hannesson [13], Ragozin and Brown [14], and Ströbele and Wacker [15] derive golden rules for optimal steady state harvesting in a multi-species context. In addition, Ragozin and Brown [14] study the approach path towards the optimal steady state. Kar and Matsuda [16] investigated a prey–predator model with Holling type of predation and harvesting of mature predator species. They proved that harvesting effort can be used as control to break the cyclic behavior of the system and drive the system to a required steady state. Kar and Matsuda [17] studied a harvested prey–predator system with Holling type III response function. They analyzed the effects of harvesting efforts on the prey–predator system. Saito et al. [18] discussed the effects of cannibalism on a basic stage-structured population model for a single species. Kar and Chaudhuri [19] investigated a dynamic reaction model in the case of a prey–predator type fishery system, where only the prey species is subjected to harvesting taking taxation as a control instrument. In this paper, they critically compare the uniqueness of their developed model with other researchers like Chaudhuri and Johnson [20], Ganguly and Chaudhuri [21], Anderson and Lee [22], Krishna et al. [23] and Pradhan and Chaudhuri [24]. Jang and Diamond [25] derived a discrete size-structured fishery model to study the population-level effects of the timing of density dependence. Ecological systems are characterized by the interaction between species and their natural environment. An important type interaction which affects population dynamics of all species is predation. Thus predator–prey models have been in focus of ecological science since the early days of this discipline [26]. Great attention has been paid to the dynamics properties of the predator –prey models which have significant biological background. Many excellent and interesting results have been obtained in Xu and Ma [27], Gao et al. [28]; Kuang and Takeuchi [29]; Song and Guo [30]. Chakraborty et al. [31] described a prey–predator model with stage structure for prey. They considered the harvesting of adult prey and predator populations and discussed the dynamical behavior of the system. The trophic interactions of food webs are significantly dependent on the capacity of predators to find, kill and consume prey population. It plays an important role in shaping the entire ecosystem. The success rate of an individual predator depends on several factors which include the existence of the populations with stage structure. The mutual interactions of the species are required to be incorporated in the predation model to understand the dynamics of the species involved. The most important component in such models is the density of prey that determines the functional response. Here our model is significantly differs from the other models as we have considered the ability of a mature predator individual to adjust its feeding rate to changes in both prey density as well as in immature predator density. Moreover, it is necessary to include the mutual interference between searching predators during its interaction with prey population together with immature predator population. In general, functional response depends on the spatial scale on which it is measured. Thus, the intrinsic growth rate of predators depends on the density of both prey population and immature predator population. Thus, the study of population dynamics with harvesting is a very interesting subject of mathematical bioeconomics. Therefore, our objective is to provide some results which are theoretically beneficial to maintaining the sustainable development of the prey–predator system as well as keeping the economic interest of harvesting at an ideal level. In the present paper, we describe the dynamics of a stage structured prey–predator fishery system. The asymptotic behavior of the proposed system is analyzed and the sufficient conditions are derived for the global stability of the system at the interior equilibrium point using geometric approach. Attempt has made to interpret the obtained results in terms of ecology to the extent possible. The control parameter is characterized and the optimal control problem is formulated. The control problem is solved using an iterative method with Runge–Kutta fourth-order scheme. The numerical results provide more realistic features of the model system. 2. The model and its qualitative properties We consider a prey–predator model with stage structure for predator and it is assumed that only mature predator population is harvested. Let us assume x, y and z are respectively the size of prey population, immature predator and mature predator population at time t. It is assumed that the populations are growing in closed homogeneous environment. At any time (t > 0), the birth rate of prey population is assumed to follow logistic growth and the birth rate of the immature predator population is assumed to be proportional to the density of existing mature predator population with proportionality parameter a. The mature predator population catches the prey population and the immature predator population at Holling type-II functional response which are respectively mxz/(a + x) and nyz/(b + y) where m and n are the respective maximal predator per capita consumption rate, i.e. the maximum number of prey and immature predator that can be eaten by a predator in each time unit and a, b are the half capturing saturation constant, i.e. the number of prey and immature predator necessary to achieve one-half of the maximum rate m and n [32]. Keeping these aspects in view, the dynamics of the system may be governed by the following system of differential equations:

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

 dx x  mxz ¼ rx 1   ; dt K aþx dy nyz ¼ az   dy; dt bþy   dz cz  hðtÞ; ¼ sz 1  dt xþy

9273

ð2:1Þ

where r is the intrinsic growth rate of prey population, K is the environmental carrying capacity, d is the death rate of immature predator population, s is the intrinsic growth rate of mature predator population and c is the equilibrium ratio between both the prey and predator population. Cod fishery is a good example of the above model. An ecological interpretation of cod cannibalism is that the cod has adopted with the cyclic recruitment pattern of capelin by eating its own offspring when the capelin stock is depleted and is in a state of rebuilding. Harvesting has a strong impact on the dynamic evaluation of a population subjected to it. First of all, depending on the nature of the applied harvesting strategy, the long run stationary density of population may be significantly smaller than the long run stationary density of a population in the absence of harvesting. Therefore, while a population can in the absence of harvesting be free of extinction risk, harvesting can lead to the incorporation of a positive extinction probability and therefore, to potential extinction in finite time. Secondly, if a population is subjected to a positive extinction rate then harvesting can drive the population density to a dangerously low level at which extinction becomes sure no matter how the harvester affects the population afterwards. We take the harvest rate h(t) in the form

hðtÞ ¼ qEz;

ð2:2Þ

where q is the catchability co-efficient, E is the fishing effort used to harvest the predator population. Thus our system becomes

 dx x  mxz ¼ rx 1   ; dt K aþx dy nyz ¼ az   dy; dt bþy   dz cz ¼ sz 1   qEz; dt xþy

ð2:3Þ

with initial conditions x(0) P 0, y(0) P 0, z(0) P 0. 3. Positivity and boundedness of the solution In this section, we try to find the conditions for which the solutions of the system (2.3) are positive as well as bounded. Theorem 3.1. If y(t) is always nonnegative then all possible solutions of the system (2.3) are positive.

Proof. From the first equation of the system (2.3) we can write,

dx ¼ x

   x mz dt; r 1  K aþx

which implies,

dx x

n  mz o ¼ /ðx; zÞdt where /ðx; zÞ ¼ r 1  Kx  aþx .

Now integrating above differential equation in the region [0, t] we have,

xðtÞ ¼ xð0Þe

R

/ðx;zÞdt

> 0;

8t:

Again, from the third equation of the system (2.3) we can write,

dz ¼ z

    cz s 1  qE dt; xþy

which implies,

dz z

n   o cz  qE . ¼ wðx; y; zÞdt where wðx; y; zÞ ¼ s 1  xþy

Now integrating above differential equation in the region [0, t] we get,

R zðtÞ ¼ zð0Þe wðx;y;zÞdt > 0 8t: Hence, in consequence to our assumption y(t) P 0 "t, it may be concluded that all the solutions of the system (2.3) are always positive. h

9274

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

In the rest part of our analysis we assume that y(t) is always nonnegative so that the solutions of the system (2.3) are always positive. In the next theorem, we try to find some sufficient condition for which the solutions of the system (2.3) are bounded. Theorem 3.2. If a < dc then the solutions of the system (2.3) are bounded above.

Proof. From the first equation of the system (2.3) we may conclude that,

xðtÞ 6 K;

8t:

Again, from the third equation of the system (2.3) we can write,

czðtÞ 6 ðxðtÞ þ yðtÞÞ; 8t or; zðtÞ 6

K

c

þ

yðtÞ

c

; 8t:

Using above equation in the second equation of the system (2.3) we obtain

  dy K y a  dy ¼ ny þ K 1 where n ¼ d  ; 6a þ c dt c c or;

K1 ¼

aK c

dy þ ny 6 K 1 : dt

Now if a < dc then n > 0 and then integrating above equation we get,

yðtÞ 6

nyð0Þ  K 1 nt e þ K1 n

i.e. y(t) 6 K1 as t ! 1. Using above relation we may write that,

zðtÞ 6

1

c

ðK þ K 1 Þ as t ! 1:

Hence, all the solutions of the system (2.1) are bounded above if a < dc. h

Note 1. If the ratio of the proportionality parameter of the immature predator population to the density of existing mature predator population and death rate of immature predator population is less than the equilibrium ratio between both of prey and predator populations then the solutions of the proposed system are bounded above. This is the ecological interpretation of the above theorem. Note 2. The above condition is the sufficient condition for boundedness of the system (2.3) but not necessary. Again, if all the solutions of the system (2.3) are positive then we say that the above condition is the sufficient condition for the solutions of the system (2.3) are bounded (i.e. both bounded below and above). Definition 3.1. The solution of the system (2.3) will said to be Invariance if the solution always remain in the same region from which it start. Theorem 3.3. The axes of x and y, the planes xy and yz and the first octant are invariant for the solution of the system (2.3) provided y(t) is always nonnegative. Proof. If we start the solution of the system (2.3) from the x-axis i.e. if the initial conditions are x(0) > 0, y(0) = z(0) = 0, then the solution of the system (2.3) always remain on the x-axis i.e. we may say x-axis is invariant for the solution of the system (2.3). Similar result will arise for the axis y. Again if the solution starts from the xy plane then initially z(0) = 0. Using the result of the theorem (3.1) we have z(t) = 0, for all t and so in this case the solution of the system (2.3) always remain on the x y plane. Similarly we show that the yz plane is also invariant for the solution of the system (2.3). Lastly if the solution starts from first octant, then x(0), y(0), z(0) P 0 and the result of the theorem (3.1) implies that the solution of our system always remains in the positive octant. Hence the theorem. h

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

9275

4. Analysis of the model at its interior equilibrium To analyze the system (2.3) at its interior equilibria we first try to find all possible non negative equilibria of the system though we, particularly, are interested about the analysis around the interior equilibrium point A⁄(x⁄, y⁄,z⁄) to the system (2.3) for its natural importance. Clearly, the system has three possible non negative equilibria namely, (i) The boundary equilibrium A0(K, 0, 0). (ii) The prey free equilibrium A1(0, y1, z1) where y1 = b(a1  d)/{(n  ba)a1 + bd}, z1 = a1y1 and a1 = (s  qE)/sc. (iii) The interior equilibrium A⁄(x⁄, y⁄, z⁄) where (x⁄, y⁄, z⁄) are the positive root of system x_ ¼ y_ ¼ z_ ¼ 0. At the interior equilibrium point, we must have,

Aðx þ yÞ ¼

 dyðb þ yÞ r x s  qE ¼ z; where A ¼ : ¼ ða þ xÞ 1  K sc ab þ ða  nÞy m

The above expression can be further written as

y ¼ A1 x2 þ A2 x þ A3 r where, A1 ¼ AmK ; A2 ¼ ðK  aÞA1  1; A3 ¼ aKA1 and

B1 y2 þ B2 y þ ðB1 þ dÞxy þ ðB2 þ bdÞx ¼ 0 where, B1 ¼ Aða  nÞ  d ¼ Bb2  An. Finally, it is possible to derive the following equation,

1ðxÞ ¼ B1 A21 x4  A1 ð2A2 B1 þ ðB1 þ dÞÞx3 þ ðB1 A22  2A1 A3 B1  A1 B2 þ A2 ðB1 þ dÞÞx2 þ ð2A2 A3 B1 þ A2 B2 þ A3 ðB1 þ dÞ þ B2 þ bdÞx þ B1 A23 þ B2 A3 ¼ 0: Now, the sufficient conditions for the existence of the interior equilibrium point of system (2.3) are as follows: (i) All  ofA2, B1, B2 are non negative and A2 < 0 and 1(K) > 0 provided A2 < 2A1K. (ii) 1 2A 1 At first, we analyze the stability criterion of the system (2.3) at the boundary equilibria then we examine the local and global stability of the system followed by bifurcation analysis at its interior equilibrium A⁄(x⁄, y⁄, z⁄). Theorem 4.1. The boundary equilibrium A0 (K, 0, 0), is locally asymptotically stable if biotechnical productivity of the mature predator species is less than the fishing effort used to harvest i.e., s < qE. Proof. The characteristic equation of the system (2.3) at A0(K, 0, 0), can be written as (k + r)(k + d)(k  s + qE) = 0. Therefore, the eigen values to the system (2.3) at A0(K, 0, 0), are given by

k ¼ r; d; ðqE  sÞ: Thus, if s < qE then the boundary equilibrium A0(K, 0, 0), is stable otherwise it is unstable. h Theorem 4.2. The sufficient conditions for the system (2.3) to be locally asymptotically stable, at its prey free equilibrium A1(0, y1, z1), are s < qE, ar < mz1 and a < ny1/(b + y1). Proof. The characteristic equation of the system (2.3) at A1(0,y1, z1), can be written as

 mz1  2 krþ ðk þ a2 k þ a3 Þ ¼ 0; a where,

a2 ¼ d þ

bnz1 ðb þ y1 Þ

2

 s þ qE þ

2scz1 y1

and

(

a3 ¼

dþ

bnz1 ðb þ y1 Þ2

!

qE  s þ

  ) 2scz1 csz2 ny1 þ 21 : a y1 y1 b þ y1

9276

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

Clearly, if ar < mz1, s < qE and a < ny1/(b + y1) then all the eigenvalues to the system (2.3), at its prey free equilibrium, are negative consequently the prey free equilibrium is locally asymptotically stable. Hence the theorem. h Now, the characteristic equation of the system (2.3) around its interior equilibrium point can be written as:

k3 þ c1 k2 þ c2 k þ c3 ¼ 0;

ð4:1Þ

where, 

c1 ¼

rx scz bnz mx z þ þd ; þ  K x þ y ðb þ y Þ2 ða þ x Þ2

ð4:2Þ

   ny 2 2 2   sc bþy   a z drx brnx z dscz rscx z mscx z bnscz þ þ þ þ þ þ c2 ¼ K Kðb þ y Þ2 x þ y Kðx þ y Þ ða þ x Þðx þ y Þ2 ðb þ y Þ2 ðx þ y Þ ðx þ y Þ2 





2

dmx z

bmnx z

ða þ x Þ

ða þ x Þ2 ðb þ y Þ

 2

 2

mscx z

2

ða þ x Þ2 ðx þ y Þ

ð4:3Þ

;

   ny  2 2 3 2 rsc bþy   a x z drscx z dmscx z ða  y Þ bmnscx z ða  y Þ bnrscx z þ þ þ þ c3 ¼ Kðx þ y Þ ða þ x Þ2 ðx þ y Þ2 ða þ x Þ2 ðb þ y Þ2 ðx þ y Þ2 Kðb þ y Þ2 ðx þ y Þ Kðx þ y Þ2    ny  3 msc bþy   a x z  : ða þ x Þ2 ðx þ y Þ2

ð4:4Þ

At the interior equilibrium point we have, scz⁄/(x⁄ + y⁄) = s  qE and 

ny =ðb þ y Þ  a ¼ dy =z : Let us take d1 = d + s  qE, then c1, c2 and c3 can be expressed as:

c1 ¼ a1  a2 m; where, a1 ¼

rx K

þ

bnz ðbþy Þ2

c2 ¼ a3  a4 m;

c3 ¼ a5  a6 m;

þ d1 ,

 

a2 ¼ a3 ¼

xz

ða þ x Þ2

;

    d1 rx bnz rx dðs  qEÞ2 y ; þ þ s  qE þ 2 K K scz ðb þ y Þ 

ðs  qEÞ2 x ; scða þ x Þ ða þ ða þ þ !  rðs  qEÞx dðx þ 2y Þ bnz ; þ a5 ¼ ðx þ y Þ K ðb þ y Þ2 !  ðs  qEÞ2 x bnz  a6 ¼ : dð2y  aÞ þ scða þ x Þ2 ðb þ y Þ2 a4 ¼

d1 x z

x Þ2



2

bnx z

x Þ2 ðb

y Þ2



Now using Routh–Hurwitz criterion we can state and prove the following theorem for the local stability of the system (2.3) around its interior equilibrium point. Theorem 4.3. The sufficient condition for the system (2.3) to be asymptotically locally stable around its interior equilibrium is min{a2/a1, a6/a5} > m > m⁄ where m⁄ is the largest root of the equation /(m) = a2a4m2 + (a6  a1a4  a2a3) m + (a1a3  a5) = 0  provided a 6 bnz  2 þ 2y . dðbþy Þ

Proof. If the interior equilibrium of the system (2.3) exists then its characteristic equation at the interior equilibrium point is given by Eq. (4.1). Now by using Routh–Hurwitz criterion we can say that all the eigen values of (2.3) at A⁄ have negative real part i.e. all the roots of (4.1) has negative real part if c1, c3 > 0 and c1c2 > c3.  Now for m < min{a2/a1, a6/a5} both of c1 and c2 are positive provided a 6 bnz  2 þ 2y . Also for m > m⁄ we have c1c2 > c3. dðbþy Þ

Hence, for min{a2/a1 ,a6/a5} > m > m⁄ the system is locally asymptotically stable around its interior equilibrium point. h Prey–predator models with constant parameters are often found to approach a steady state in which the species coexist in equilibrium. But if parameters used in the model are changed, other types of dynamical behavior may occur and the critical

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

9277

parameter values at which such transitions happen are called bifurcation points. The purpose of this study is to determine the stability behavior of the system in presence of different density-dependent factors of the prey–predator interactions. To study the transition of the system with respect to the small changes in the density dependent factors, we have considered, m, n and a, as bifurcation parameters and m⁄, n⁄ and a⁄ represent the critical value or the bifurcating value of the concerned bifurcation parameter. Theorem 4.4. The system (2.3) undergoes through a Hopf bifurcation at its interior equilibrium for m = m⁄ where m⁄ can be obtained from the equation /(m) = 0. Proof. For m = m⁄ we have c1c2  c3 = 0 therefore it is possible to obtain the eigen values of the forms k1 = c1 and pffiffiffiffiffi k2;3 ¼ i c2 . However, let us assume that k1 = l1(m) and k2,3 = l(m) ± ix(m). Now it is very easy to show that dl/dm is nonzero at the point m = m⁄. Hence, using the conditions given in Venkatsubramanian et al. [33], it is easy to show that our considered system (2.3) undergoes a Hopf bifurcation at its interior equilibrium for m = m⁄. h Note: It may also be noted by using theorem (4.3) that the system (2.3) is locally stable for m > m⁄ and unstable for m < m ⁄. Again let us rewrite c1, c2 and c3 from (4.2)–(4.4) in the following form:

c1 ¼ b1 þ b2 n; where, b1 ¼

b2 ¼

rx K

c2 ¼ b3 þ b4 n;



mx z ðaþx Þ2

bz



ðb þ y Þ2

c3 ¼ b5  b6 n;

þ d1 ,

;

! d1 rx d1 x z ðs  qEÞ2 x dðs  qEÞ2 y þ  ; m b3 ¼ K scz ða þ x Þ2 scða þ x Þ      2 bz rx bmx z b4 ¼ ; þ s  qE þ ðb þ y Þ2 K ða þ x Þ2 ðb þ y Þ2 rdðs  qEÞðx þ 2y Þx mdðs  qEÞ2 ð2y  aÞx ;  Kðx þ y Þ scða þ x Þ2 ! bðs  qEÞx z r mðs  qEÞ b6 ¼  : K scða þ x Þ2 ðb þ y Þ2 b5 ¼

Theorem 4.5. The system (2.3) undergoes a bifurcation at its interior equilibrium for n = n⁄ where n⁄ = b5/b6 provided n⁄ > max{0, b1/b2, b3/b4}. Proof. According to the statement of the theorem for n = n⁄ we have c3 = 0. Thus, the characteristic equation is reduced to k3 + c1k2 + c2k = 0. It is clear from Eq. (4.1) that both c1 and c2 are positive since n⁄ > max{0, b1/b2, b3/b4}. Consequently, the characteristic equation has a simple zero eigen value around its interior equilibrium for n = n⁄. Hence, the system (2.3) passes through a bifurcation at n = n⁄ around its interior equilibrium point. h Theorem 4.6. The system (2.3) also undergoes a bifurcation at its interior equilibrium for a = a⁄ where a⁄ can be obtained by solving c3 = 0 provided a⁄ is positive and c1, c2 are also positive. Proof. It is clear from the statement of the theorem that when a = a⁄ the characteristic equation of the system (2.3) is reduced to k3 + c1k2 + c2k = 0 since c3 = 0. Consequently, the system has a simple zero eigenvalue around its interior equilibrium when a = a⁄. Hence, the system (2.3) passes through a bifurcation at a = a⁄ provided a⁄ is positive. h Note. It is easy to show that the system is locally asymptotically stable for a > a⁄ and unstable for a < a⁄. 5. Global stability Here, we illustrate the general method described by Li and Muldowney [34] to show an n-dimensional autonomous dynamical system f : D ! Rn ; D  Rn ; an open and simply connected set and f 2 C1(D), where the dynamical system is given by

dx ¼ f ðxÞ dt

ð5:1Þ

9278

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

is globally stable under certain parametric conditions. For detailed calculations one can see Haque et al. [35], Bunomo et al. [36], Kar and Mondal [37]. Now, we assume the following conditions. (A1) The autonomous dynamical system (5.1) has a unique interior equilibrium point x⁄ in D. (A2) The domain D is simply connected. (A3) There is a compact absorbing set X  D. Definition 5.1. The unique equilibrium point x⁄ of the dynamical system (5.1) is globally asymptotically stable in the domain D if it is locally asymptotically stable and all the trajectories in D converges to its interior equilibrium point x⁄. Let J = (Jij)n be the variational matrix of the system (5.1) and Jj2j be the second additive compound matrix with order n C2  nC2. In particular, for n = 3 we can write

0 J j2j ¼

j2j

@f B ¼@ @x

1

V 11 þ V 22

V 23

V 13

V 32

V 11 þ V 33

V 12

V 31

V 21

V 22 þ V 33

C A:

Let M(x) in C1(D), Li and Muldowney [34], be the n C2  n C2 matrix valued function. Moreover, we also consider B be a matrix such that B = MfM1 + MJj2jM1where the matrix Mf is represented by

ðM ij ðxÞÞf ¼

 t @M ij  f ðxÞ ¼ rM ij  f ðxÞ: @x

Again, we consider the Lozinskii measure C of B (Martin [38]) with respect to a vector norm j.j in RN ; N ¼ n C 2 then we have

CðBÞ ¼ limþ h!0

jl þ hBj  1 : h

If, (A1), (A2) and (A3) hold then Li and Muldowney [34] show that

lim sup sup

1 t

Z

t

CðBðxðs; x0 ÞÞÞds < 0:

ð5:2Þ

0

The condition (5.2) ensures that there are no orbits (i.e. homoclinic orbits, heteroclinic cycles and periodic orbits) which give rise to a simple closed rectifiable curve in D, invariant for the system (5.1). It is also a robust Bendixson criterion. Now, we use the above discussion to show that our system (2.3) is globally stable around its interior equilibrium. The autonomous system (2.3) can be written in the following form,

dX ¼ f ðXÞ; dt

ð5:3Þ

 1 0 1 rx 1  Kx  mxz ; aþx x nyz B az  bþy  dy; C where f ðXÞ ¼ @  A and X ¼ @ y A.  cz z  qEz sz 1  xþy 0

Then, the variational matrix, V(x, y, z), of the system (2.3) can be written as

0 V¼

B @f B ¼B @X @

mxz þ ðaþxÞ  rx 2 K

0 2

scz  ðxþyÞ 2

mx aþx

0

1

C ny C bnz  ðbþyÞ a þ bþy 2  d C: A scz2 scz  ðxþyÞ2  xþy

If Vj2j be the second additive compound matrix of V then Vj2j, Bunomo et al. [36], can be expressed as

0 V j2j

B B ¼B @

mxz ðaþxÞ2

bnz  rx  ðbþyÞ 2  d K 2

scz  ðxþyÞ 2 scz2 ðxþyÞ2

ny a þ bþy mxz ðaþxÞ2

scz  rx  xþy K

0

mx  aþx

0 scz bnz  ðbþyÞ 2  d  xþy

We consider M(X) in C1(D) in such a way that M = diag {x/z, x/z, x/z}. Then, we have M1 = diag {z/x, z/x, z/x} and

Mf ¼

dM _  ðx=z2 Þz_ ; x=z _  ðx=z2 Þz_ g: _  ðx=z2 Þz_ ; x=z ¼ diagfx=z dX

Thus, it is easy to show that,

1 C C C: A

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

9279

_  z_ =z; x=x _  z_ =z; x=x _  z_ =zg and MV j2j M 1 ¼ V j2j : Mf M 1 ¼ diagfx=x We have,

B ¼ Mf M1 þ MV j2j M 1 ¼



B11

B12

B21

B22

 ;

where

mxz rx bnz x_ z_ B11 ¼  þ    d; x z ða þ xÞ2 K ðb þ yÞ2 B12 ¼ B21 ¼

 

ny a þ bþy

2

mx aþx

scz2 ðxþyÞ2

scz  ðxþyÞ 2

 ; t

and

0_ B22 ¼ @

x x

scz mxz rx  zz_ þ ðaþxÞ 2  K  xþy

0

0 x_ x

scz bnz  zz_  ðbþyÞ 2  d  xþy

1 A:

Now let us define the following vector norm in R3 as

jðu; v ; wÞj ¼ maxfjuj; jv j þ jwjg; where (u, v, w) is the vector in R3 and it is denoted by C, the Lozinskii measure with respect to this norm. Therefore, C(B) 6 sup{p1, p2} where pi = C1(Bii) + jBijj for i = 1, 2 and i – j, where jB12j, jB21j are matrix norms with respect to the L1 vector norm and C1 is the Lozinskii measure with respect to that norm. Consequently, we can obtain the following terms

z_ mxz rx bnz    d; z ða þ xÞ2 K ðb þ yÞ2     mx  ny  a; ; jB12 j ¼ max  bþy aþx x_ x

C1 ðB11 Þ ¼  þ

jB21 j ¼

scz2 ðx þ yÞ2

;

( ) scz mxz rx bnz x_ z_ þ max C1 ðB22 Þ ¼    ; d : x z xþy ða þ xÞ2 K ðb þ yÞ2

  cz  qE. Again, from the system (2.3), we can write zz_ ¼ s 1  xþy Using the above expression we have,

      ny  mx cz mxz rx bnz x_  þ qE þ  a; p1 ¼  s 1     d þ max  2 2 xþy K ðb þ yÞ bþy aþx x ða þ xÞ      _x scz mxz ny mx rx bnz þ  a;   þ max   d; ¼ þ qE  s þ x þ y ða þ xÞ2 bþy aþx K ðb þ yÞ2 x ( ) scz2 mxz rx bnz x_  s þ max  ; d p2 ¼ þ qE þ x ða þ xÞ2 K ðb þ yÞ2 ðx þ yÞ2 ( ) scz2 rx mxz bnz x_ ¼ þ qE þ  s  min ; þd :  K ða þ xÞ2 ðb þ yÞ2 x ðx þ yÞ2 Now, from the expressions of p1 and p2 we obtain

(     ny  mx rx rx bnz scz mxz x_  ; þ   CðBÞ 6 þ qE  s  min þ d   max a b þ y  a þ x ;K K ðb þ yÞ2 x þ y ða þ xÞ2 x ) mxz scz2 bnz scz2 :  þ ; þdþ 2 2 2 ða þ xÞ ðx þ yÞ ðb þ yÞ ðx þ yÞ2

9280

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

It is assumed that there exists a positive l1 2 R and t1 > 0 such that l1 = inf{x(t), y(t), z(t)} whenever t > t1. Also, we take





 nl  ml1 l2 ¼ max  1  a; b þ l1 a þ l1 and

l3



( ) r l1 bnl1 ml21 scl1 r l1 ml21 sc bnl1 sc ¼ min þd   l2 ; þ ; þdþ þ  : K K 4 ðb þ l1 Þ2 ða þ l1 Þ2 2l1 ða þ l1 Þ2 4 ðb þ l1 Þ2

Then, from the above expression we can write

) ( r l1 bnl1 scl1 ml21 r l1 ml21 sc bnl1 sc x_ CðBÞ 6 þ qE  s  min þd   l2 ; þ ; þdþ þ  x K 2l1 ða þ l1 Þ2 K 4 ðb þ l1 Þ2 ða þ l1 Þ2 4 ðb þ l1 Þ2 x_ ¼ þ qE  s  l3 ; x x_ i:e: CðBÞ 6  ðs þ l3  qEÞ; x Z 1 t 1 xðtÞ CðBÞds 6 log  ðs þ l3  qEÞ; i:e: t 0 t xð0Þ Z t 1 CðBðxðs; x0 ÞÞÞds < ðs þ l3  qEÞ < 0: lim sup sup t!1 t 0 In consequence to the above analysis we have reached to state the following theorem. Theorem 5.1. The system (2.1) is globally asymptotically stable around its interior equilibrium if (s + l3) > qE where

l3 and

( ) r l1 bnl1 ml21 scl1 r l1 ml21 sc bnl1 sc ¼ min þd   l2 ; þ ; þdþ þ  K K 4 ðb þ l1 Þ2 ða þ l1 Þ2 2l1 ða þ l1 Þ2 4 ðb þ l1 Þ2

 n o  ml1 nl1 þ l2 ¼ max bþ l1  a; aþl1 with l1 2 R such that for t1 > 0 we havel1 = inf{x(t), y(t), z(t)} whenever t > t1.

6. Optimal control problem In commercial exploitation of renewable resources the fundamental problem, from the economic point of view, is to determine the optimal trade-off between present and future harvests. The emphasis of this section is on the profit-making aspect of fisheries. It is a thorough study of the optimal harvesting policy and the profit earned by harvesting, focusing on quadratic costs and conservation of fish population by constraining the latter to always stay above a critical threshold. The prime reason for using quadratic costs is that it allows us to derive an analytical expression for the optimal harvest; the resulting solution is different from the bang-bang solution which is usually obtained in the case of a linear cost function. It is assumed that price is a function which decreases with increasing biomass. Thus, to maximize the total discounted net revenues from the fishery, the optimal control problem can be formulated as

JðEÞ ¼

Z

tf

edt ½ðp  xqEzÞqEz  cEdt;

ð6:1Þ

t0

where p is the constant price per unit biomass, c is the constant cost of harvesting effort, x is the economic constants and d is the instantaneous annual discount rate. The problem (6.1), subject to population Eq. (2.3) and control constraints 0 6 E 6 Emax, can be solved by applying Pontryagin’s maximum principle. The convexity of the objective function with respect to E, the linearity of the differential equations in the control and the compactness of the range values of the state variables can be combined to give the existence of the optimal control. Suppose Ed is an optimal control with corresponding states xd, yd and zd. We take Ad (xd, yd, zd) as optimal equilibrium point. Here we intend to derive optimal control Ed such that

JðEd Þ ¼ maxfJðEÞ : E 2 Ug; where U is the control set defined by U = {E:[t0, tf] ? [0, Emax]jE is Lebesgue measurable}. Now the Hamiltonian of this optimal control problem is

  

 

x  mxz nyz cz  þ k2 az   dy þ k3 sz 1   qEz ; H ¼ ðp  xqEzÞqEz  cE þ k1 rx 1  K aþx bþy xþy

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

9281

where k1, k2, k3 are adjoint or costate variables. Here the transversality conditions give ki(tf) = 0, i = 1, 2, 3. Now, it is possible to find the characterization of the optimal control Ed. On the set {tj0 < Ed(t) < Emax}, we have

@H ¼ pqz  2xq2 z2 E  c  k3 qz: @E Thus at Ad(xd, yd, zd), E = Ed(t) and This implies that,

Ed ¼

@H @E

¼ pqzd  2xq2 z2d Ed  c  k3 qzd ¼ 0.

pqzd  c  k3 qzd : 2xq2 z2d

ð6:2Þ

Now, the adjoint equations at the point Ad(xd, yd, zd) are

" ! #  dk1 @H 2rxd amzd csz2d þ k3 ; ¼ dk1  k1 r  ¼ dk1   @x Ad dt K ða þ xd Þ2 ðxd þ yd Þ2 " # !  dk2 @H bnzd csz2d ; ¼ k2 d þ þ d  k3 ¼ dk2   2 @y Ad dt ðb þ yd Þ ðxd þ yd Þ2

       dk3 @H mxd nyd 2sczd þ k  pqEd  2xq2 E2d zd : ¼ k  k a  d  s þ qE þ ¼ dk3  1 2 3  @z Ad dt a þ xd b þ yd xd þ y d

ð6:3Þ ð6:4Þ

ð6:5Þ

Eqs. (6.3)–(6.5) are first order system of simultaneous differential equations and it is easy to get the analytical solution of the equations with the help of initial conditions ki(tf) = 0, i = 1, 2, 3. In this regard, it is to be noted that we have formulated the optimal control problem through considering fishing effort as control parameter and the optimal control problem will be numerically solved using a forward–backward sweep technique of 4th order Runge–Kutta method to pursue numerical simulations in the next section. Therefore, we summarize the above analysis by the following theorem: Theorem 6.1. There exists an optimal control Ed and corresponding solutions to the system (2.3) xd, yd and zd that maximizes J(E) over U. Furthermore, there exists adjoint functions k1, k2 and k3 satisfying Eqs. 6.3, 6.4, 6.5 with transversality conditions ki(tf) = 0, 3 qzd i = 1, 2, 3. Moreover, the optimal control is given by Ed ¼ pqz2d ck xq2 z2 . d

7. Numerical simulations For the purpose of simulation experiments we mainly use the software MATLAB 7.6 and Mathematica 5.2. As the problem is not a case study, the real world data are not available for this model. We, therefore, take here some hypothetical data with the sole purpose of illustrating the analytical results that we have established in the previous sections. Moreover, it may be noted that as the parameters of the model are not based on real world observations, the main features described by the simulations presented in this section should be considered from a qualitative, rather than a quantitative point of view. However, numerous scenarios covering the breath of the biological feasible parameter space were conducted and the results shown above display the gamut of dynamical results collected from all the scenarios tested.

P1 ¼ fr; K; a; a; n; b; d; s; c; q; Eg ¼ f1:1; 20; 4:8; 0:4; 6:9; 9:9; 1:3; 0:25; 0:2; 0:1; 0:76g: For this set of parametric value we see that the critical value of m = m⁄ is 0.9. According to our theoretical result the system is locally asymptotically stable for m > m⁄ and unstable for m < m⁄ around its interior equilibrium (1.64, 0.46, 6.44). In Figs. 1 and 3, we show the asymptotic local stability in the form of solutions curve and phase space diagram for m = 0.96 and m = 0.99 respectively. Also in Fig. 2 and in Fig. 4, we have shown the instability of the system for m = 0.86 and m = 0.8 respectively in the form of solutions curve and phase space diagram. Another interesting observation is the Hopf bifurcation to the system at its interior equilibrium for not only for the critical value of the parameter m but also for the critical value of almost every parameters present in the system. For example, we give respectively the stable and unstable solution curves for the system (2.3) by taking a as the bifurcation parameter. Now we choose a hypothetical set of parameters

P2 ¼ fr; K; a; m; n; b; d; s; c; q; Eg ¼ f1:1; 20; 4:8; 0:86; 6:9; 9:9; 1:3; 0:25; 0:2; 0:1; 0:76g: It is possible to find the critical value of a from Theorem 4.5. In Fig. 2, we have already drawn that the solutions curve to the system (2.3) for a = 0.4(a⁄) while a⁄ = 0.5 (approximately) with the help of parameter set P2. Clearly Fig. 5 shows that the system (2.3) is asymptotically stable around its interior equilibrium.

9282

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

14 prey immature predator mature predator

12

populations

10

8

6

4

2

0

0

200

400

600

800

1000

1200

1400

time Fig. 1. Solution curve of the system (2.3) for m(=0.96) > m⁄ with the parameter set P1.

14 prey immature predator mature predator

12

populations

10

8

6

4

2

0

0

200

400

600

800

1000

1200

1400

time Fig. 2. Solution curve of the system (2.3) for m(=0.86) < m⁄ with the parameter set P1.

We have obtained another interesting result around the prey free equilibrium point of the system though considering the following hypothetical set of parameters:

P3 ¼ fr; K; m; a; a; n; b; d; s; c; q; Eg ¼ f1:1; 20; 6:8; 3:8; 0:4; 6:9; 9:9; 0:3; 0:25; 0:2; 0:1; 0:76g: In Fig. 6, in the form of the solution curves, it is evident that the system is asymptotically locally stable around its prey free equilibrium point. The numerical simulation of optimal control [39] under various parameter sets can be done using a forward–backward sweep technique of 4th order Runge–Kutta method to solve the system state Eq. (2.3) and their corresponding adjoint

9283

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

14 12

z

10 8 6 4 2 3 15

2 10 1

5 0

y

0

x

Fig. 3. Phase space diagram of the system (2.3) for m(=0.99) > m⁄ with the parameter set P1.

20

z

15

10

5

0 3 15

2 10 1 y

5 0

0

x ⁄

Fig. 4. Phase space diagram of the system (2.3) for m(=0.8) > m with the parameter set P1.

Eqs. (6.3)–(6.5). In this method, initially we make a guess for optimal control and then solve the system of state Eq. (2.3) forward in time using the Runge–Kutta method with the initial conditions (x0, y0, z0). Then, using the state values, the adjoint Eqs. (6.3)–(6.5) are solved backward in time using the Runge–Kutta method with the transversality conditions. At this point, the optimal control is updated using the values for the state and adjoint variables. The updated control replaces the initial control and the process is repeated until the successive iterates of control values are sufficiently close. The convergence of such an iterative method is based on the work of Hackbush [40]. The intrinsic growth rate r of prey population has an important role on the dynamics of the system. It is clearly observed from Fig. 7 that, in presence of harvesting, the density of the populations is directly proportional to the intrinsic growth rate of prey population. It is also evident from the figure that the fishing effort used to harvest the populations is inversely proportional to the intrinsic growth rate of prey population. It may be noted that at the initial level fishing effort gradually

9284

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

10 prey immature predator mature predator

9 8

populations

7 6 5 4 3 2 1 0

0

200

400

600

800

1000

1200

1400

time Fig. 5. Solution curves of the system (2.3) for a = 0.6(>a⁄) with the parameter set P2.

1.8 1.6 prey immature predator mature predator

1.4

populations

1.2 1 0.8 0.6 0.4 0.2 0 -0.2

0

200

400

600

800

1000

1200

1400

time Fig. 6. Solution curves of the system (2.3) with the parameter set P3.

decreases with the increasing time but after a particular time of span it remains unchanged consequently the populations converge to a finite quantity. The interesting fact of the system is that the predator population gets increased with the increasing intrinsic growth rate of prey population and the existence of the mature predator population is ensured due to the increasing density of the prey as well as immature predator populations. The impacts of density dependent factors like predation are illustrated in Figs. 8 and 9. According to the formation of the model, the mature predator population consumes its members of its own type or kind and prey population. Thus, due to the effect of cannibalism, the density of mature predator population gets decreased for higher rate of predation of mature predator population from other populations. It is also observed that harvesting has great impact on the dynamics of the system

9285

15

14

Mature predator population

13

0

5

10 Time

15

20

0.465 0.46 0.455 0.45 0.445

0

10

20 30 Time

40

50

6.4 Fishing effort

Prey population

16

Immature predator population

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

6.2 6 5.8

0

20

40

0.27

0.26

0.25

60

0

50

Time

100 Time

150

200

Prey population

15

14.5

14

13.5

0

10

20

30

Immature predator population

Fig. 7. Variation of optimal prey and predator biomass and fishing effort with the increasing time. The solid line corresponds to r = 1, the dotted line to r = 1.1 and the dashed line to r = 1.3.

0.46

0.455

0.45

0.445

0

50 Time

100

6.4 0.28 6.2

Fishing effort

Mature predator population

Time

6 5.8 5.6

0.27 0.26 0.25

0

10

20 30 Time

40

50

0

50

100 Time

150

200

Fig. 8. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to m = 0.9, the dotted line to m = 0.96 and the dashed line to m = 1.1.

since the rate of harvesting exhibits solution effects with respect to the stock, it is interesting to observe that the fishing effort gets decreased with the increasing time and ultimately control the co-existence of the populations. Hence, the fully dynamic interaction between the predation and fishing effort plays an important role for the existence of the entire ecosystem. It is clear from Fig. 10 that the proportionality constant a of immature predator population has an important implication to the dynamics of the system. It is noted from the figure that for lower proportionality constant a, both the immature and

14.55

14.5

14.45 50 Time

100

0.5

0.45

0.4

0

5 Time

10

0

50 Time

100

0.28 6.1

Fishing effort

Mature predator population

0

Immature predator population

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

Prey population

9286

6.08 6.06 6.04 6.02 0

50 Time

0.275 0.27 0.265 0.26

100

14.55

14.5

14.45 50 Time

100

6.15

0.5

0.4

0.3 0

1

2 3 Time

4

5

0.28 Fishing effort

Mature predator population

0

Immature predator population

Prey population

Fig. 9. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to n = 6.0, the dotted line to n = 6.9 and the dashed line to n = 7.5.

6.1

6.05 0

50 Time

100

0.27

0.26 0

50 Time

100

Fig. 10. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to a = 0.3, the dotted line to a = 0.4 and the dashed line to a = 0.5.

mature predator populations gets decreased consequently the prey population gets increased. At the initial level fishing effort used to harvest decreases but after a definite span of time it is remain unaltered. The natural death rate of immature predator population plays an important role to the dynamics of the system. As the cannibalism of predator population is incorporated to the model therefore the density of mature predator population is dependent on the density of the immature predator population. Fig. 11 depicts that the density of the immature predator

9287

Prey population

14.55

14.5

14.45 50 Time

100

0.5

0.45

0.4 0

1

2 3 Time

4

5

6.1 0.275 Fishing effort

Mature predator population

0

Immature predator population

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

6.08 6.06

0.27 0.265

6.04 0

50

100

0.26

150

0

Time

50 Time

100

14.8 14.6 14.4 14.2 10

20 30 Time

40

50

0.46

0.455

0.45

0.445

0

50 Time

100

0

50 Time

100

0.28 6.2 Fishing effort

Mature predator population

0

Immature predator population

Prey population

Fig. 11. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to d = 1.0, the dotted line to d = 1.3 and the dashed line to d = 1.8.

6

0.27

0.26

5.8 0

20

40 Time

60

80

Fig. 12. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to s = 0.15, the dotted line to s = 0.25 and the dashed line to s = 0.4.

population decreases with the increasing death rate of immature predator population consequently the density of the mature predator population decreases. However, density of the prey population gets increased due to the decreasing density of the mature predator population. It is clear from Fig. 12 that the intrinsic growth rate of mature predator population has great influence on the dynamics of the system. It is noted from the figure that the density of the predator populations increases with the increasing intrinsic growth rate therefore prey population gets decreased. However, it is also clear from the figure that the fishing effort

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

Prey population

15 14.5 14

0

10

20 30 Time

40

50

6.5

Fishing effort

Mature predator population

13.5

Immature predator population

9288

6

0.46

0.44

0.42 0

5

10 Time

15

10

0.25

5.5 0

5 Time

0.2

20

0

50

100 Time

150

200

Prey population

14.6

14.5

14.4

14.3 50 Time

100

0.458

0.456

0.454

0.452

0

50 Time

100

6.2 6.15

Fishing effort

Mature predator population

0

Immature predator population

Fig. 13. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to c = 1.9, the dotted line to c = 2.2 and the dashed line to c = 2.5.

6.1 6.05 6 0

50 Time

100

0.4 0.3 0.2 0.1 0

0

1

2 3 Time

4

5

Fig. 14. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to q = 0.05, the dotted line to q = 0.1 and the dashed line to q = 0.2.

decreases with the increasing intrinsic growth rate and after a particular time of span it remains unchanged. Hence, the optimal control (fishing effort) can be suitably designed so that the populations may get protection from extinction. The equilibrium ratio between both the prey and predator populations is an important factor towards the cannibalism of predator population. The carrying capacity of the mature predator population is dependent on the density of the prey as well as immature predator populations therefore it is evident from Fig. 13 that the mature predator population should decrease with the increasing value of c and the consequent effects are clearly observed for the immature predator and prey populations. It is also obvious that the fishing effort should be increased with the increasing value of c.

9289

14.52 14.5

Mature predator population

14.48

0

50 Time

100

0.455

0.4545

0.454

0.4535

0

50 Time

100

0

50 Time

100

0.28 6.08

Fishing effort

Prey population

14.54

Immature predator population

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

6.06

0.27

0.26

6.04 0

50 Time

100

0.25

Fig. 15. Variation of optimal prey, predator and fishing effort with the increasing time. The solid line corresponds to d = 0.001, the dotted line to d = 0.01 and the dashed line to d = 0.03.

It is evident from Fig. 14 that the mature predator population decreases with the increasing rate of catchability coefficient consequently the density of the immature predator population also decreases but the density of the prey population gets increased. The economic parameters such as price per unit biomass of catch, fishing cost per unit effort and discount rate determine the stock level maximizing the present value of the flow of resource rent over time. Moreover, characterization of the system is extremely dependent on the sensitivity of the discount rate consequently the efficiency of the fishery in economic sense can only be achieved if the impacts of different discount rate is taken under consideration. It is obvious from Fig. 15 that the fishing effort gets increased with the increasing discount rates consequently the predator populations gets decreased and the prey population gets increased. 8. Concluding remarks This paper deals with a prey–predator type fishery model with stage structure for the predator and the harvesting of the mature predator species. The interactions among immature and mature species are based on the assumptions as follows: (i) the recruitment of the immature species depend on the size of the mature species (ii) two substocks interact via cannibalism and (iii) immature predator population which are the product of the mature predator also become mature. Though our paper is not based on a case study but North East Atlantic Cod fishery may be one good example for our model. Cod has adopted to the cycle recruitment patter of capelin by eating its own offspring when the capelin stock is depleted. To examine the dynamic behavior of the system we have first discussed the existence of possible steady states and their local as well as global stability. Global stability of the system is shown by using a suitable geometrical approach. It is observed that there are three possible equilibria exist, one as predator free, one as prey free and the most important one as the interior equilibrium. An important and one of the interesting question in mathematical ecology is how predator populations are able to survive without prey population as sometimes it is observed in real world situation. It is observed in our model and it is happened due to incorporation of cannibalism of predator population. Our result suggest that capture rate of mature predators (both on prey and immature predator) may lead to the stable system to become unstable through a Hopf bifurcation when these parameters passes their critical value. It is also observed that our system may exhibit cyclic behavior and due to this cyclic nature of the system some population show periodic fluctuation in abundance with periodic crashes. From the point of view of the bioeconomic context, optimal control problem is formulated and solved in both analytically and numerically. Sensitivity analysis is carried out with respect to some key parameters of the system. Though this work is not a case study but our analytical results may be helpful for the researchers in this field. It may also be pointed out that in this paper several important parameters such as ecological fluctuations, refuge, interaction with other species etc. are disregarded. Hence, further research is necessary to accomplish the needs in this field.

9290

K. Chakraborty et al. / Applied Mathematics and Computation 218 (2012) 9271–9290

Acknowledgements We are very grateful to the anonymous reviewers for their careful reading, constructive comments and helpful suggestions, which have helped us to improve the presentation of this work significantly. The first author also gratefully acknowledges Director, INCOIS for his cordial attitude and unconditional help. This is INCOIS contribution number 100. Research of Soovoojeet Jana is supported by University Grants Commission, Government of India (F. 11-2/2002 (SA-1) dated 19 August, 2011). References [1] C.C. Chen, C.Y. Hsui, Fishery policy when considering the future opportunity of harvesting, Mathematical Biosciences 207 (2007) 138–160. [2] T.K. Kar, K. Chakraborty, Effort dynamics in a prey–predator model with harvesting, International Journal of Information Systems Science 6 (3) (2010) 318–332. [3] Z. Sebestyén, Z. Varga, J. Garay, R. Cimmaruta, Dynamic model and simulation analysis of the genetic impact of population harvesting, Applied Mathematics and Computation 216 (2) (2010) 565–575. [4] D. Xiao, W. Li, W. Han, Dynamics in a ratio dependent predator prey model with predator harvesting, Mathematical Analysis and Applications 324 (1) (2006) 4–29. [5] M.D. Garza-Gill, M.M. Varela-Lafuente, J.C. Suris-Regueiro, European hake fishery bioeconomic management (Southern stock) applying effort tax, Fisheries Science 60 (2003) 199–206. [6] T.K. Kar, K. Chakraborty, A bioeconomic assessment of the Bangladesh shrimp fishery, World Journal of Modelling and Simulation 7 (1) (2011) 58–69. [7] H. Matsuda, K. Nishimori, A size-structured model for a stock-recovery program for an exploited endemic fisheries resource, Fisheries Science 1468 (2002) 1–14. [8] W. Wang, G. Mulone, F. Salemi, V. Salone, Permanence and stability of a stage-structured predator prey model, Journal of Mathematical Analysis and Applications 262 (2001) 499–528. [9] B. Leard, J. Rebaza, Analysis of predator–prey models with continuous threshold harvesting, Applied Mathematics and Computation 217 (12) (2011) 5265–5278. [10] K. Chakraborty, S. Das, T.K. Kar, Optimal control of effort of a stage structured prey–predator fishery model with harvesting, Nonlinear Analysis: Real World Applications 12 (6) (2011) 3452–3467. [11] C.W. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, Wiley Series, New York, 1990. [12] C.W. Clark, Mathematical models in the economics of renewable resources, SIAM Review 21 (1) (1979) 81–99. [13] R. Hannesson, Optimal harvesting of ecologically interdependent fish species, Journal of Environmental Economics and Management 10 (1983) 329– 345. [14] D.L. Ragozin, G.J. Brown, Harvest policies and nonmarket valuation in a predator prey system, Journal of Environmental Economics and Management 12 (1985) 55–168. [15] W.J. Ströbele, H. Wacker, The economics of harvesting predator–prey systems, Journal of Economics 61 (1995) 65–81. [16] T.K. Kar, H. Matsuda, Controllability of a harvested prey–predator system with time delay, Journal of Biological Systems 14 (2) (2006) 243–254. [17] T.K. Kar, H. Matsuda, Global dynamics and controllability of a harvested prey–predator system with Holling type III functional response, Nonlinear Analysis: Hybrid Systems 1 (2007) 59–67. [18] Y. Saito, C. Jung, Y. Lim, Effects of cannibalism on a basic stage structure, Applied Mathematics and Computation 217 (5) (2010) 2133–2141. [19] T.K. Kar, K.S. Chaudhuri, Regulation of a prey–predator fishery by taxation: a dynamics reaction model, Journal of Biological Systems 11 (2) (2003) 173–187. [20] K.S. Chaudhuri, T. Johnson, Bioeconomic dynamics of a fishery modeled as an S-system, Mathematical Biosciences 99 (1990) 231–249. [21] S. Ganguly, K.S. Chaudhuri, Regulations of a single species fishery by taxation, Ecological Modelling 82 (1995) 51–60. [22] L.G. Anderson, D.R. Lee, Optimal governing instrument, operation level and enforcement in natural resource regulation: the case of the fishery, American Journal of Agricultural Economics 68 (1986) 678–690. [23] S.V. Krishna, P.D.N. Srinivasu, B. Kaymakcalan, Conservation of an ecosystem through optimal taxation, Bulletin of Mathematical Biology 60 (1998) 569–584. [24] T. Pradhan, K.S. Chaudhuri, A dynamic reaction model of two species fishery with taxation as a control instrument: a capital theoretic analysis, Ecological Modelling 121 (1999) 1–16. [25] S.R.-J. Jang, S.L. Diamond, Population-level effects of density dependence in a size-structured fishery model, Applied Mathematics and Computation 139 (1) (2003) 133–155. [26] A.A. Berryman, The origins and evolutions of predator–prey theory, Ecology 73 (1992) 1530–1535. [27] R. Xu, Z. Ma, Stability and Hopf bifurcation in a ratio-dependent predator prey system with stage-structure, Chaos, Solitons and Fractals 38 (2008) 669– 684. [28] S.J. Gao, L.S. Chen, Z.D. Teng, Hopf bifurcation and global stability for a delayed predator–prey system with stage structure for predator, Applied Mathematics and Computation 202 (2008) 721–729. [29] Y. Kuang, Y. Takeuchi, Predator–prey dynamics in models of prey dispersal in two-patch environments, Mathematical Biosciences 120 (1994) 77–98. [30] X. Song, H. Guo, Global stability of a stage-structured predator–prey system, International Journal of Biomathematics 1 (3) (2008) 313–326. [31] K. Chakraborty, M. Chakraborty, T.K. Kar, Optimal control of harvest and bifurcation of a prey–predator model with stage structure, Applied Mathematics and Computation 217 (21) (2011) 8778–8792. [32] R. Arditi, L.R. Ginzburg, Coupling in predator–prey dynamics: ratio-dependence, Journal of Theoretical Biology 139 (1989) 311–326. [33] V. Venkatsubramanian, H. Schattler, J. Zaborszky, Local bifurcation and feasibility regions in differential–algebraic systems, IEEE Transactions on Automatic Control 40 (12) (1995) 1992–2013. [34] M.Y. Li, J.S. Muldowney, A geometric approach to global stability problems, SIAM Journal on Mathematical Analysis 27 (4) (1996) 1070–1083. [35] M. Haque, J. Zhen, E. Venturino, An ecoepidemiological predator prey model with standard disease incidence, Mathematical Methods in the Applied Sciences 32 (7) (2008) 875–898. [36] B. Bunomo, A. Onofrio, D. Lacitignola, Global stability of an SIR epidemic model with information dependent vaccination, Mathematical Biosciences 216 (1) (2008) 9–16. [37] T.K. Kar, P.K. Mondal, Global dynamics and bifurcation in a delayed SIR epidemic model, Nonlinear Analysis: Real World Applications 12 (2011) 2058– 2068. [38] R.H. Martin Jr., Logarithmic norms and projections applied to linear differential systems, Journal of Mathematical Analysis and Applications 45 (1974) 432–454. [39] J.T. Workman, S. Lenhart, Optimal Control Applied to Biological Models, Chapman and Hall/CRC, 2007. [40] W. Hackbush, A numerical method for solving parabolic equations with opposite orientations, Computing 20 (3) (1978) 229–240.