Analysis of a predator-prey model for exploited fish populations with schooling behavior

Applied Mathematics and Computation 317 (2018) 35–48

Contents lists available at ScienceDirect

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

Analysis of a predator-prey model for exploited fish populations with schooling behavior Debasis Manna a, Alakes Maiti b, G.P. Samanta c,∗ a

Department of Mathematics, Surendranath Evening College, Kolkata 700009, India Department of Mathematics, Vidyasagar Evening College, Kolkata 700006, India c Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India b

a r t i c l e

i n f o

Keywords: Predator-prey model Fishing effort Schooling Stability Hopf bifurcation

a b s t r a c t In this paper, a predator-prey model for exploited fish populations is considered, where the prey and the predator both show schooling behavior. Due to this coordinated behavior, predator-prey interaction occurs only at the outer edge of the schools formed by the populations. Positivity and boundedness of the model are discussed. Analysis of the equilibria is presented. A criterion for Hopf bifurcation is obtained. The optimal harvest policy is also discussed using Pontryagin’s maximum principle, where the effort is used as the control parameter. Numerical simulations are carried out to validate our analytical findings. Implications of our analytical and numerical findings are discussed critically. © 2017 Elsevier Inc. All rights reserved.

1. Introduction It is a fact that many fish species live in groups. They enjoy many benefits from living in groups, such as higher success in finding mates, reduction in the risk of predation, increase in the foraging success, protection from bad weather, etc. Such collective behaviors are of two types, one is shoaling and the other is schooling. When members of a group swim independently in such a way that they stay connected, then this behavior is known as shoaling. If all members are swim in the same direction in a coordinated fashion with same speed then such a behavior is called schooling. Schooling fishes are usually of the same species and of the same size. Again shoalers and schoolers are of two types. Obligate shoalers or schoolers exhibit shoaling or schooling behavior all the time. Facultative shoalers or schoolers show collective behavior for finding mates or some other reasons [25]. Schools that are traveling can form thin lines, or squares or ovals or amoeboid shapes. Usually interactions of different species take several forms, depending on whether the influences are beneficial or detrimental to the species involved. Among these interactions, predator-prey relationship is considered to be an extremely important one. It is true that the preys always try to develop the methods of evasion to avoid being eaten. However, it is certainly not true that a predator-prey relationship is always harmful for the preys, it might be beneficial to both. Further, such a relationship often plays an important role to keep ecological balance in nature. Mathematical modeling of predator-prey interaction was started in the 1920s. Interestingly, the first predator-prey model in the history of theoretical ecology was developed independently by Lotka (a US physical chemist) and Volterra (an Italian mathematician) [17,29]. Subsequently, ∗

Corresponding author. E-mail addresses: [email protected] (D. Manna), [email protected] (A. Maiti), [email protected], [email protected]. ac.in, [email protected] (G.P. Samanta). http://dx.doi.org/10.1016/j.amc.2017.08.052 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

36

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

this model has been used as a machine to introduce numerous mathematical and practical concepts in theoretical ecology. Many refinements of the Lotka–Volterra model have also been made to overcome the shortcomings of the model and to get better insights of predator-prey interactions. We notice an important assumption in such modeling. In these models, it is assumed that the individuals of both predator and prey populations live independently so that any predator can interact with any prey. Consequently, the interaction term is proportional to the product of the density of both populations. This concept is used for many decades by many authors. Many fish species, which are in predator-prey relationship, have economic importance also. Therefore, developing their management strategies is an important area of research. The literature abounds with evidences of many such predator-prey systems. For example, Arctic cod (Gadus morhua) and its commercially important prey species: capelin (Mallotus villosus) or herring (Clupea harengus) [14]. Another example is tuna (Thunnus spp.), which feeds mainly on anchovy (e.g. Engraulis mordax) and sardine (e.g. Sardinops sagax) [5,22]. Most of the above mentioned fishes (both prey and predator) exhibit schooling behavior. For example, adult cods are usually found in dense schools [6]. Schooling behavior of juvenile tuna is well known [24]. On the other hand, small fishes like capelin, herring and sardine are among very spectacular schooling fishes [8,13,23]. Hence consideration of schooling behavior in predator-prey systems of fishes is of utmost importance. Although a number of predator-prey models for fish populations have been cultured in theoretical ecology (see [10,11,18,27] and references there in), but the effect of schooling behavior has remained unuttered. However, there should be no denying that, in case of many fish populations, ignoring schooling behavior is simply ignoring reality. From the above viewpoint, we consider a predator-prey model for fish populations, where both prey and predator live in schools. A number of interesting results on stability are established. Some results on extinction of the populations are obtained. A criterion for Hopf bifurcation is established. Optimal harvest policy (to be adopted by the fishery management) is discussed and dynamic optimization of the harvest policy is carried out using Pontryagin’s Maximum Principle. The paper is structured as follows. In Section 2, we present the mathematical model with basic considerations. Positivity and boundedness of solutions of the model are established in Section 3. Section 4 contains the detailed analysis of the equilibria, their stability analysis and a criterion for Hopf bifurcation. Bionomic equilibrium points, optimal harvest policy at equilibrium level, and dynamic optimization of the optimal harvest policy are presented in Section 5. To illustrate our analytical findings, computer simulations of variety of solutions of the system are performed; and the results are presented in Section 6. Section 7 contains the general discussion of the paper and biological significance of our analytical findings. 2. The mathematical model At time t, let x(t) denotes the density of the prey fish, and y(t) denotes the density of the predator fish. In the following, we explain the basic considerations that motivate us to introduce our basic mathematical model. (i) So far as the growth of the prey (in the absence of the predator) is concerned, many modelers have suggested the logistic growth function to be a logically reliable choice. The function is introduced in 1838 by the Belgian mathematician Verhulst [28]. The logistic growth equation is given by





x dx = rx 1 − , dt k

(2.1)

where r is the intrinsic per capita growth rate and k is the carrying capacity of the environment. The logic behind this is very simple. As the resources (e.g., space, food, essential nutrients) are limited, every population grows into a saturated phase from which it cannot grow further; the ecological habitat of the population can carry just so much of it and no more. This indicates that the per capita growth rate is a decreasing function of the size of the population, and reaches zero as the population achieves a size k (in the saturated phase). Further, any population reaching a size that is above this value will experience a negative growth rate. The term −rx2 /k may also be regarded as the loss due to intraspecific competition. (ii) Usually, it is assumed that the individuals of both predator and prey populations live independently so that any predator can interact with any prey. Consequently, the interaction term is proportional to the product of the density of both populations. This concept is used for many decades by many authors. If predators form a school and they attack a prey population which lives independently; then the predators who occupy the edge of the school would get maximum benefit. The idea of such community behaviors of predators had been given by Cosner et al. [12]. According to this idea, the interaction term should be proportional to the product of density of prey and the square root of the predator density. Similarly, when preys form a school and predators live independently, then the most harmed preys during predator hunting are those staying on the boundary of the school; and so the interaction term should be modified in an analogous manner. Unfortunately, such an innovative idea has not been used by the researchers for about a decade. The work of Chattopadhyay et al. [9] may be regarded as a strong recognition of this concept. Then came the works of Ajraldi et al. [1] and Braza [7], which have given such modeling a new dimension. Their works have stimulated developments in the modeling of group behaviors among various populations. The idea is very interesting. Let P be the density of a population that gathers in groups, and suppose that group occupies an area A. The number of individuals staying at outermost positions in the group is proportional to the length of the perimeter of the patch where√the group is located. Clearly, √ its length is proportional to A. Since P is distributed over a two-dimensional domain, P would therefore count the

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

37

individuals at the edge of the group. This is the main idea of Ajraldi et al. [1] and Braza [7]. These works have motivated research on the modeling of group behavior, and a number of works have been done on this area in last three or four years [2–4,19–21,30]. Here we have assumed that both prey and predator exhibit schooling behavior. Accordingly, when attack of a schooling predator (of density y) on a schooling prey (of density x) is to be modeled, the interaction term will be proportional to the product of the square-roots of the densities of the two populations, because the interaction occurs √ √ only at the outer edges of the schools. In other words, the interaction term will be b x y, where b is the consumption rate of the predator. We also assume that the predator depends only on this specific prey. (iii) We assume that both the prey and the predator are subject to a combined (non-selective) harvesting effort E. Then, based on the catch-per-unit-effort hypothesis [11], the catch rate functions for the prey and predator are q1 Ex and q2 Ey, respectively, where q1 and q2 are the catchability coefficients of the prey and the predator species, respectively. On the basis of the above considerations, we introduce the following predator-prey model within the framework of the following set of nonlinear ordinary differential equations:





√ √ x dx = rx 1 − − b x y − q1 Ex, x(0 ) > 0, dt k √ √ dy = c x y − dy − q2 Ey, y(0 ) > 0. dt

(2.2)

The parameter r is the intrinsic growth rate of the prey, k is the carring capacity of the environment for the prey. d represents the death rate of the predator and E denotes the combined harvesting effort. Here b is the consumption rate and c is the conversion rate for the predator. The parameters q1 , q2 are catchability coefficients of the two species. It is an obvious assumption that all the parameters are non-negative. 3. Positivity and boundedness Positivity and boundedness of a model guarantee that the model is biologically well behaved. For positivity of the system (2.2), we have the following theorem. Theorem 3.1. All solutions of the system (2.2) that start in R2+ remain positive forever. The proof is simple and therefore it is omitted. The following theorem ensures the boundedness of the system (2.2). Theorem 3.2. All solutions of the system (2.2) that starts in R2+ are uniformly bounded. Proof. Let (x(t), y(t)) be any solution with positive initial conditions (x(0), y(0)). Since



x dx ≤ rx 1 − dt k



So, by a standard comparison theorem, we have

lim sup x(t ) ≤ γ , t→∞

where γ = max{x(0 ), k}. Let

W = cx(t ) + by(t ). Then

W˙ = cx˙ + by˙





x − cq1 Ex − bdy − bq2 Ey k ≤ c (r + 1 )x − cx − bdy

= crx 1 −

≤ c ( r + 1 )x − τ W ≤ c ( r + 1 )γ − τ W where

τ = min{1, d}. That is

W˙ + τ W ≤ c (r + 1 )γ . So,

0 ≤ W (x(t ), y(t )) ≤

c ( r + 1 )γ

τ

+ e−τ t W (x(0 ), y(0 )).

38

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

When t → ∞, then

0≤W ≤

c ( r + 1 )γ

τ

.

Thus, all solutions of (2.2) enter into the region



D = (x, y ) : 0 ≤ W ≤ Hence the theorem.

c ( r + 1 )γ

+  , for any

τ

 >0 .



4. Stability and bifurcation analysis 4.1. The boundary equilibria The system (2.2) always has the equilibrium point (0, 0). The equilibrium ( xˆ, 0 ) exists if E < r/q1 . If this condition is ¯ ¯ q E rq −dq rq −dq satisfied, then xˆ = k(1 − 1r ). Let E¯ = 12 q2 q 1 + 2q l q , E = 12 q2 q 1 − 2q l q , l¯ = (rq2 + dq1 )2 − 4bcq1 q2 . 1 2

1 2

1 2

1 2

The system (2.2) is not linearizable about (0, 0) and (xˆ, 0 ). So local stability of them can not be studied. However, we have studied some global behaviors of (0, 0) and (xˆ, 0 ). Theorem 4.1. If E >

r q1

then

lim x(t ) = 0.

t→∞

Proof. We have from the first equation of (2.2) that





√ x dx = rx 1 − − b xy − q1 Ex, dt k ≤ − ( q1 E − r )x < 0. This implies that limt→∞ x(t ) = 0.



Theorem 4.2. The equilibrium point (0, 0) is globally stable if E¯ < E < r/q1 . Proof. When E¯ < E < r/q1 , then the system (2.2) has only two equilibria, namely, (0, 0) and (xˆ, 0 ). Also





x˙ = x (r − q1 E ) − b < bx



c − d + q2 E



c d+q2 E

>

r−q1 E . b

So,



r y − x , x k



rx y − . x bk

(4.1)

Let x0 , y0 denote values of x(t), y(t) at t = 0 respectively. Now, we consider the following three cases. √ √ c x (i) If (x0 , y0 ) lies in the region where (d+ < y, then from (4.1), we have x˙ < 0. Also y˙ < 0. That is, both x(t), y(t) decrease q2 E ) with time. So after certain time, the solution curve intersects the line y = ( d+cq E )2 x. Then y(t) is constant but x˙ < 0. So 2 √ √ c x the trajectory comes back to the region (d+ < y, and finally it moves towards the origin. q2 E ) √ √ r−q E √ (ii) If (x0 , y0 ) lies in the region where b1 x < y < d+cq E x, then x˙ < 0, y˙ > 0. That is, x(t) decreases and y(t) increases 2

with time. So after some time, the solution curve intersects the line y = ( d+cq E )2 x, whence y(t) is constant but x˙ < 0. 2 √ √ c x Therefore, the trajectory again returns back to the region (d+ < y , and finally it approaches the origin. q2 E ) √ r−q E √ (iii) If (x0 , y0 ) lies in the region where y < b1 x, then y˙ > 0. That is, y(t) increases with time. So after certain time, the solution curve crosses the line y = ( √ √ c x the region (d+ < y. q E)

r−q1 E 2 ) x. b

Then x(t) is decreases and y(t) increases with time, and finally it moves to

2

Thus, the trajectories starting from any point in R2 will eventually approach (0, 0). Hence (0, 0) is globally stable.



Theorem 4.3. The equilibrium point (xˆ, 0 ) of the system (2.2) is always unstable. √ √ Proof. The equilibrium point (xˆ, 0 ) lies in region D2 = {(x, y ) : c x − (d + q2 E ) y ≥ 0} and y˙ ≥ 0 for all (x, y) ∈ D2 ; the equal√ √ ity occurs when (x, y) lies on the line c x − (d + q2 E ) y = 0 or on x axis. Now consider an arbitrary  neighborhood of (xˆ, 0 ) (say D ) given by D (xˆ, 0 ) = {(x, y ) ∈ R2+ , (x − xˆ)2 + y2 <  2 ,  > 0}. If we choose the initial position (x0 , y0 ) = (xˆ,  /2 ) ∈ D , then y˙ at (x0 , y0 ) is positive. So y does not approach to zero. Arbitrariness of  implies that in every neighbourhood of (xˆ, 0 ), we can find a point (x0 , y0 ) such that any trajectory starting at (x0 , y0 ) will not approach (xˆ, 0 ). In other words, (xˆ, 0 ) is an unstable critical point. 

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

39

To determine the behavior near the origin for min{0, E } < E < E¯ , we consider the following two cases. √ Case 1: x, y are very small (close to zero) and y < < x < 1. In this case, it is clear that, xy << x so that near the origin,





r dx ≈ x ( r − q1 E ) − x , dt k √ dy = c xy − dy − q2 Ey. dt

(4.2)

Since E < E¯ < r/q1 , so for any x0 < k(r − q1 E )/r, x˙ > 0 at (x0 , y0 ). Where x0 , y0 denote initial value of x(t) and y(t) at time 2 k(r−q1 E ) t = 0. After certain time, solution curve intersect the line x = . Then x˙ = 0. Finally y → (d+cq xE )2 when t → ∞. r 2 √ Case 2: x, y are very small (close to zero), x < < y < 1 and E < E¯ . In this case, it is clear that, 1 − kx ≈ 1 and xy << y so that near the origin,

√ dx ≈ rx − b xy − q1 Ex, dt dy ≈ −dy − q2 Ey. dt

(4.3)

Solving Eq. (4.3), we get

√ x=

√ √ √  (r−q E )t (d+q2 E )t b y0 b y0 1 e− 2 + x0 − e 2 , r − q1 E + d + q2 E r − q1 E + d + q2 E

y = y0 e−(d+q2 E )t .

(4.4) (4.5)

√ √ When (x0 , y0 ) lies in the region D1 = {(x, y ) : (r − q1 E + d + q2 E ) x − b y > 0}, then as t → ∞, y → 0 and x(t) → ∞. That is y(t) decreases and x(t) is increases with time. If we choose (x0 , y0 ) from the region y > x, then after certain time, the solution curve crosses the line y = x and finally enters the region where y < x. Then from case-1 we conclude that the solution curve will not ultimately return to the point (0, 0). √ √ When (x0 , y0 ) lies in the region where (r − q1 E + d + q2 E ) x − b y < 0, then as t → ∞, y → 0. That is y(t), x(t) both decrease with time and x(t) vanishes for some positive value of y, after which the predator y declines to zero (as dy/dt < 0). So the trajectory ultimately approaches the point (0, 0). Such observations are similar to those obtained by Braza [7]. 4.2. Analysis of the interior equilibrium Theorem 4.4. The interior equilibrium (x∗2 , y∗2 ) exists if any one of the following conditions hold: (i) bc − rd < 0 and 0 ≤ E ≤ E¯ , (ii) bc − rd > 0, (rq2 + dq1 )2 > 4bcq1 q2 , rq2 − dq1 > 0 and E ≤ E ≤ E¯ . If either (i) or (ii) is satisfied, then x∗2 , y∗2 are given by √ c x∗ q E kbc x∗2 = (k(1 − 1r ) − r (d+ ), y∗2 = ( d+q 2E )2 . q E) 2

2

Remark 1. It is easy to notice that E¯ < qr when (x∗2 , y∗2 ) exists. Therefore, if max{0, E } < E < E¯ , then both (xˆ, 0 ), (x∗2 , y∗2 ) 1 exist. On the other hand, if E¯ < E < qr then only (xˆ, 0 ) exists. In all other cases, both the non-trivial equilibrium points do 1 not exist. These ranges of effort are very useful from practical point of view. Theorem 4.5. If (x∗2 , y∗2 ) exist with u < E < v, then (x∗2 , y∗2 ) is locally asymptotically stable, where v = (r+d )q2 −q1 d−l ,l= q2 (2q1 −q2 )



(r+d )q2 −q1 d+l , u= q2 (2q1 −q2 )

(rq2 + dq1 )2 − 3bcq2 (2q1 − q2 ).

Proof. The variational matrix V2 at E2∗ is given by

⎛ 

b y∗2

r ∗ ⎜ 2 x ∗ − k x 2 ⎜ 2  V2 = ⎜ ⎝ c y∗2



2 x∗2

 ⎞

− −

b x∗2

 ⎟

2 y∗2 ⎟  ⎟. c x∗2 ⎠



2 y∗2

The characteristic equation of V2 is

λ

2

cr ∗ − λP + x 2k 2



x∗2 = 0, y∗2

where

P=

1 (E − u )(E − v ). 6 b( d + q 2 E )

(4.6)

40

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

Since, u < E < v, we have P < 0. Thus, sum of the roots of the Eq. (4.6) is negative and product of the roots is positive. Hence the theorem.  4.3. Hopf bifurcation Theorem 4.6. If (x∗2 , y∗2 ) exists with q2 − 2q1 > 0; d2 + 2rd − 3bc < 0 then a Hopf bifurcation occurs at E = E ∗ , where E∗ is the unique positive real root of q2 (q2 − 2q1 )E 2 + (2d (q2 − q1 ) + 2rq2 )E + d2 + 2rd − 3bc = 0. Proof. Since q2 − 2q1 > 0; d2 + 2rd − 3bc < 0, the equation

q2 (q2 − 2q1 )E 2 + (2d (q2 − q1 ) + 2rq2 )E + d2 + 2rd − 3bc = 0 has only one positive real root E∗ (say). Then the characteristic Eq. (4.6) has a pair of purely imaginary roots. Again



m dP  = − {q2 (q2 − 2q1 )E + (d (q2 − q1 ) + rq2 )} = 0.  dE E =E ∗ c Thus all conditions of the Hopf bifurcation theorem are satisfied and hence a simple Hopf bifurcation occurs at E = E ∗ .



5. Bionomic equilibrium points, optimal harvest policy and dynamic optimization 5.1. Bionomic equilibrium points Let cˆ be the cost per unit effort and p1 , p2 are the price per unit biomass of prey and predator, respectively. Then the net economic revenue to the society is given by

π (x, y, E, t ) = ( p1 q1 x + p2 q2 y − cˆ)E

(5.1)

The term bionomic equilibrium is an amalgamation of the concepts of biological equilibrium as well as economic equilibrium. In case of bionomic equilibrium, x˙ = y˙ = π = 0. Here the bionomic equilibrium points are

   c2     cˆ  bc bc r (0, 0, 0 ), (k, 0, 0 ), k 1 − , 2k 1− ,0 , , 0, 1− rd

d

rd

p1 q1

q1

cˆ p1 q1 k



.

(5.2)

Also, there exists a non-zero bionomic equilibrium (x∞ , y∞ , E∞ ), which satisfy the following equations



x r 1− k





− q1 E = b

y , x



c

x = d + q2 E, y

p1 q1 x + p2 q2 y = cˆ.

(5.3)

One can easily obtain

x∞ =

cˆ(d + q2 E∞ )2 , y∞ = p1 q1 ( d + q2 E∞ )2 + p2 q2 c 2



c d + q2 E∞

and E∞ is the positive real root of the equation



p1 q21 q32 E 4 + (dq1 − rq2 ) p1 q1 q22 + 2dp1 q21 q22 +



2

x∞ ,

(5.4)



r cˆ 3 3 q E k 2

+ (bc − rd ) p1 q1 q22 + q1 q2 ( p1 q1 d2 + p2 q2 c2 ) + 2(dq1 − rq2 ) p1 q1 q2 d +



+ (bc − rd )2d p1 q1 q2 + ( p1 q1 d2 + p2 q2 c2 )(dq1 − rq2 ) + +

r cˆ 3 d + (bc − rd )( p1 q1 d2 + p2 q2 c2 ) = 0. k





r cˆ 3dq22 E 2 k

r cˆ 2 3d q2 E k (5.5)

5.2. Optimal harvest policy Here we discuss the optimal harvest policy to be adopted by the fishery management to maximize the present value J of a continuous time-stream of revenues from the fishery using effort as a control variable when the system is at equilibrium level. Then



J=

0

∞

L(x, E, t )dt ,

where

L(x, E, t ) = e−δt ( p1 q1 x + p2 q2 y − cˆ)E.

(5.6)

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

41

Here x = [x, y]T is the vector of state variables which evolve according to the state equations (2.2) written as x˙ = f(x, E ), x(0 ) = x0 and E(t) is the control variable. We notice that δ is the instantaneous annual discount rate. We consider that the present value of capital flow over time depends on discount rate δ . The discount rate would therefore determine the stock level, maximizing the present value of the flow of the resource rent over time. This reference point is sometimes referred to as the optimal economic yield biomass. Our intention is to maximize (5.6) subject to the state equations (2.2) and to the control constraint

0 ≤ E (t ) ≤ Emax

(5.7)

using Pontryagin’s Maximum Principle [26]. 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 allows us to give the existence of the optimal control. The Hamiltonian function is formulated by

H (x, E, t ) = L(x, E, t ) + λ (t )f(x, E ), T

where λ(t) is the vector of costate or adjoint variables. √ √ For our problem, f(x, E ) = [rx(1 − kx ) − b xy − q1 Ex, c xy − dy − q2 Ey]T and we take λ1 (t), λ2 (t) as the adjoint variables so that the adjoint variable vector becomes λ(t ) = [λ1 (t ), λ2 (t )]T . Therefore, the Hamiltonian is given by

 

H (x, y, E, t ) = e−δt ( p1 q1 x + p2 q2 y − cˆ)E + λ1 (t ) rx 1 −

x k







√ √ − b xy − q1 Ex + λ2 (t ) c xy − dy − q2 Ey



(5.8)

For H to be maximum on the control set 0 ≤ E(t) ≤ Emax , we must have

∂H = 0. ∂E That is,

e−δt F4 (x, y ) − λ1 (t )q1 x − λ2 (t )q2 y = 0,

(5.9)

where

F4 (x, y ) = p1 q1 x + p2 q2 y − cˆ. Pontryagin’s minimum principle states that the optimal state trajectory, optimal control, and corresponding adjoint variable vector must satisfy the following: T

−λ˙ (t ) = Hx , which are known as the adjoint equations. Obviously, in our problem, the adjoint equations become

   by  ∂H x c y = e−δt p1 q1 E + λ1 (t ) r 1 − 2 − − q1 E + λ2 (t ) , ∂x k 2 x 2 x  c x  dλ2 ∂H b x − = = e−δt p2 q2 E − λ1 (t ) + λ2 (t ) − d − q2 E . dt ∂y 2 y 2 y −

dλ1 = dt

(5.10)

(5.11)

For non-zero optimal equilibrium solutions, x˙ = 0 = y˙ (that is, x, y are independent of t) and we have



x r 1− k



c





−b

y − q1 E = 0, x

(5.12)

x − ( d + q2 E ) = 0. y

(5.13)

From above, it is clear that E also independent of t. Using this, we get



dλ1 rx b − = e−δt p1 q1 E + λ1 (t ) − + dt k 2 −

dλ2 b = e−δt p2 q2 E − λ1 (t ) dt 2



  y x

x c − λ2 (t ) y 2

c + λ2 (t ) 2





y x

x y

(5.14)

(5.15)

Differentiating (5.9) and substituting value of λ˙1 , λ˙2 we have

  r b√ c√ λ1 eδt − q1 x2 + yx(q1 − q2 ) + λ2 eδt xy(q1 − q2 ) = δ F4 − q21 xp1 E − p2 q22 Ey. k

2

2

(5.16)

42

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

Now (5.9) and (5.16) give

  √ δ q2 y − 2c xy(q1 − q2 ) − q2 yE ( p1 q21 x + p2 q22 y ) λ1 eδt = q q √ ,   2 1 2 xy(by + cx ) − q1 q2 r xky − 2b yx (q2 y )2 − 2c yx (q1 x )2 2   √ F4 − δ q1 x + 2b xy(q1 − q2 ) − qk1 r x2 + q1 xE ( p1 q21 x + p2 q22 y ) λ2 eδt = q q √ .   2 1 2 xy(by + cx ) − q1 q2 r xky − 2b yx (q2 y )2 − 2c yx (q1 x )2 2 F4

(5.17)

(5.18)

Now eliminating E from (5.12) and (5.13), we get



x r 1− k





−b

 x

y q1 = c x q2

y



−d .

(5.19)

This gives the optimal trajectory in the steady state described by the optimal populations x = xδ , y = yδ . Substituting λ1 and λ2 in (5.15), we get optimal equilibrium level of effort given by



Eδ =



p2 q22 q1 xy −



rx k

+

b 2

y x

y

x

−

√ crx x 2ky

√ − δ − q2 p2 q21 2c x xy − p1 q21 x q1 x

δ+

F4 q1 x − δ 2 + δ −



rx k

+

b 2

x

−

c 2

y





 c 2

x y

√ − 2b q2 xy

.

(5.20)

As δ → ∞, the Eq. (5.20) leads to the obvious result that F4 = 0. That is,

p1 q1 x∞ + p2 q2 y∞ = cˆ ⇒ π (x∞ , y∞ , E ) = 0, which implies that an infinite discount rate leads to complete dissipation of economic revenue. This conclusion is also drawn by Clark [11]. 5.3. The dynamic optimization Now our intention is to maximize (5.6) subject to the state equations (2.2) and to the control constraint

0 ≤ E (t ) ≤ Emax using Pontryagin’s Maximum Principle when the control variable E(t) is a dynamic (i.e. time-dependent) variable. We also determine optimal control E∗ (t) and the corresponding optimal response x∗ (t), y∗ (t) which maximize (5.6) subject to the conditions stated above. The Hamiltonian function for the dynamic optimization problem is

 

H (x, y, E, t ) = e−δt ( p1 q1 x + p2 q2 y − cˆ)E + λ1 (t ) rx 1 −

x k







√ √ − b xy − q1 Ex + λ2 (t ) c xy − dy − q2 Ey



(5.21)

where λ1 (t), λ2 (t) are the adjoint variables. Adjoint equations are

   by  ∂H x c y −δ t = e p1 q1 E + λ1 (t ) r 1 − 2 − − q1 E + λ2 (t ) , ∂x k 2 x 2 x  c x  dλ2 ∂H b x − = = e−δt p2 q2 E − λ1 (t ) + λ2 (t ) − d − q2 E . dt ∂y 2 y 2 y dλ1 − = dt

(5.22)

(5.23)

For H to be maximum on the control set 0 ≤ E(t) ≤ Emax , we must have

∂H = 0. ∂E That is,

HE = e−δt F4 (x, y ) − λ1 (t )q1 x − λ2 (t )q2 y = 0,

(5.24)

where

F4 (x, y ) = p1 q1 x + p2 q2 y − cˆ. From Eq. (5.24), we get another relation

dHE = Me−δt + (N − δ q1 x )λ1 (t ) + dt

c√ 2

xy(q1 − q2 ) − δ q2 y

 λ2 (t ) = 0,

(5.25)

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

where

 

x k

M = p1 q1 rx 1 −



43



√ √ rx2 b√ − b xy + p2 q2 (c xy − dy ), N = − q1 + xy(q1 − q2 ). k 2

Solving Eqs. (5.24) and (5.25) we get

λ1 (t ) = λ2 (t ) =

Be−δt

(q1 − q2 )(q1 cx − bq2 y ) Ae−δt

(q1 − q2 )(q1 cx − bq2 y )

where



x B = p1 q1 q2 xy r 1 − k



 −b

√ xy 2

+ kr q1 q2 x2 y

√ xy 2

+ kr q1 q2 x2 y

c y − x 2

  x y

(5.26)

(5.27)

+ p2 q22

 cy√xy 2



− dy2 +

√ c xy { p1 q21 x + p2 q1 q2 y − (q1 − q2 )cˆ} − δ q2 yF4 , 2 (5.28)

A=

p1 q21 x2





x −r 1−2 k



  y

b + 2



b rx + p2 q2 q1 xy (d + ) − k 2



x √  rx b  y  b xy + { p1 q1 x + p2 q2 y − cˆ} − cˆq1 x − + δ q1 F4 x. 2 k 2 x

y −c x

  x y (5.29)

Again differentiating (5.25), we get

d2 HE = U (x, y )E (t ) + V (x, y ) = 0, dt 2





U (x, y ) = q21 p1 x r − 1 +

 rx2

+λ1 (t ) q21

k

+

3x k



(5.30)





√ + δ + q22 p2 yd + q1 q2 xy( p1 b − p2 c ) e−δt







c√ b√ xy(q1 − q2 )2 + λ2 (t ) − xy(q1 − q2 )2 + δ p2 q22 ye−δt 4 4





√ 5rx p1 q1 r2 x 2 2 (k − x )(k − 2x ) + by − b xy 3r − 2 k k

V (x, y ) =







√ − bcx + bd xy e−δt







p2 q2 c √ rx2 q1 3b + xy(rk − rx − 3dk ) − bcy + c2 x + 2yd2 e−δt + λ1 (t ) 2 k k 2



√ br + xy(q1 − q2 ) 4

√



−1+

 rc

3x k



+ (b − c )

c 4



(5.31)

x d ( 2c − b ) + y 4





y −r x





cd rx2 (q1 + q2 ) + − δ 2 (q1 xλ1 (t ) 4 4 2k   −2r √ +q2 yλ2 (t )) + δ x2 q1 λ1 (t ) + xy(q1 − q2 )(bλ1 (t ) + cλ2 (t )) k +λ2 (t )

xy(q1 − q2 )

−b+

(5.32)

For optimality of singular control, it is necessary to satisfy the generalized Legendre condition. That is,

∂  d2 HE  ≥ 0. ∂ E dt 2

(5.33)

The condition (5.33) holds if the following conditions are satisfied



k δ 1− 3 r



< x,

(5.34)

p1 c ≥ , p2 b

(5.35)

λ1 c > . λ2 b

(5.36)

44

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

Fig. 1. Nullclines and equilibria of the system (2.2) when r = 6, k = 44, c = 2, d = 0.1, q1 = 0.3, q2 = 0.3, b = 3, E = 6.

Fig. 2. Here r = 6, k = 44, c = 2, d = 0.1, q1 = 0.3, q2 = 0.3, b = 3. (a) Stability behavior of (0, 0) for different choices of x(0) and y(0) when E¯ < E = 18 < r/q1 . (b) Stability behavior of (0, 0) for different choice of x(0) and y(0) when E = 21 > r/q1 .

The Eq. (5.34), shows that, the bionomic optimal (x∗ , y∗ ) is established when prey population is greater than one third of the carrying capacity. It is also noted from the Eq. (5.36) that, if biomass conversion rate of prey population to predator population is greater than 1 then the shadow price of the prey exceeds that of the predator at the optimal level of fishery. In modern terminology, λ1 (t), λ2 (t) denote the shadow price of the capital stock x, y at time t. This term reflects the fact that the value of the capital asset is calculated from its future productivity instead of taking it as its direct sale value. From the Eq. (5.35) we also conclude that the existence of an optimal singular control is ensured if the actual price of each prey fish species be greater than or equal to the conversion factor (of prey to predator) times actual price of the predator species. Eq. (5.30)gives the optimal singular control E∗ (t) as

E ∗ (t ) = −

V ( x∗ , y∗ ) . U ( x∗ , y∗ )

Solving (2.2), we get singular optimal trajectory (where (x∗ , y∗ ) lies) after substituting E = E ∗ .

(5.37)

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

45

Fig. 3. Here r = 6, k = 44, c = 2, d = 0.1, q1 = 0.3, q2 = 0.3, b = 1.3, E = 2.3 < E¯ . Phase portrait of the system (2.2) for different initial choices, showing that (38.94, 0) is an unstable equilibrium point.

Fig. 4. Here r = 6, k = 44, c = 2, d = 0.1, q1 = 0.3, q2 = 0.3, b = 3, E = 6. (a) Phase portrait of the system (2.2) showing that the interior equilibrium (7.6421, 8.4671) is locally asymptotically stable when x(0 ) = 7 and y(0 ) = 8. (b) Stable behavior of x(t), y(t) with time.

Remark 2. If δ = 0, then

dH (x, y, E, λ1 , λ2 ) = dt

∂H ∂H ∂H ∂H ˙ ∂H ˙ ∂H ˙ + x˙ + y˙ + λ + λ + E ∂t ∂x ∂y ∂λ1 1 ∂λ2 2 ∂ E ∂H ∂H ˙ = − λ˙ 1 x˙ − λ˙ 2 y˙ + x˙ λ˙ 1 + y˙ λ˙ 2 + E ∂t ∂E ∂H ∂H ˙ = + E ∂t ∂E

(5.38)

dH vanishes (as Hamiltonian has no explicit dependence on time, the second partial derivative vanishes for an dt optimal control or E˙ vanishes for a boundary solution). This implies that the value of Hamiltonian is constant over time along the optimal trajectory. In this case, the optimal control and its responses must satisfy the following condition [15,16]:

Clearly,

46

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

Fig. 5. Here all the parameter values are same as in Fig. 4 except E = 5. (a) Phase portrait of the system (2.2) showing a periodic orbit growing out of the interior equilibrium when x(0 ) = 6.8 and y(0 ) = 8.7. (b) Oscillations of x(t), y(t) with time.

Fig. 6. Here r = 6, k = 10, c = 2, d = 0.18, q1 = 0.3, q2 = 0.3, b = 3, p1 = 1, p2 = 1, cˆ = 2. Bionomic equilibrium points of the systems (2.2) and (5.1).

H (x∗ (t ), y∗ (t ), E ∗ (t ), λ1 (t ), λ2 (t )) = 0.

(5.39)

The equation of singular optimal trajectory is given by

 

B1 rx 1 −

x k





√ √ − b xy + A1 (c xy − dy ) = 0,

(5.40)

where B1 and A1 denote the values of B, A, respectively, at δ = 0. Optimal singular control is obtained from the Eq. (5.37) after substituting δ = 0. 6. Numerical simulation In this section we present computer simulation of the system (2.2) to illustrate our analytical findings. We choose the parameters of the system (2.2) as r = 6, k = 44, c = 2, d = 0.1, q1 = 0.3, q2 = 0.3, b = 3, E = 15 < E¯ , then all the three equilibrium points of the system exist (see the nullclines presented in Fig. 1).

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

47

It is already mentioned that the system (2.2) is not linearizable about the equilibrium points (0, 0) and (xˆ, 0 ). However, we have tried to present some behaviors of (0, 0) and (xˆ, 0 ) numerically (see Figs. 2 and 3). For r = 6, k = 44, c = 2, d = 0.1, q1 = 0.3, q2 = 0.3, b = 3, E = 6, it can be shown that the co-existence equilibrium (7.6421, 8.4671) is locally asymptotically stable. The corresponding phase portrait is depicted in Fig. 4a. Clearly it is a stable spiral converging to the equilibrium. Fig. 4b shows that the prey and predator populations approach to their equilibrium values as time progresses. If we reduce the value of E to 5.4 and keep other parameter values fixed (as in Fig. 4), it is seen that (x∗2 , y∗2 ) loses its stability with the formation of a limit cycle. The corresponding phase portrait is presented in Fig. 5a. Sustained oscillations of the prey and predator are shown in Fig. 5b. Finally, in Fig. 6, we present the bionomic equilibrium points of the system. The magenta curve is given by x˙ = 0, p1 q1 x + p2 q2 y = cˆ. The black curve is given by y˙ = 0, p1 q1 x + p2 q2 y = cˆ. The cyan curve is given by p1 q1 x + p2 q2 y = cˆ, E = 0. The blue and red curves are given by x˙ = 0, E = 0 and y˙ = 0, E = 0, respectively. 7. Concluding remarks This paper deals with a predator-prey model for fish populations under harvesting. Here both the prey and predator show schooling behavior. It has been proved that the solutions of this system remains positive forever and they are uniformly bounded. These, in turn, imply that the system is biologically well behaved. We have established that the prey population dies off if harvesting effort is greater than the BTP of the prey population. Further, if the prey dies out, then it is impossible for the predator to survive. A number of interesting criteria for stability of the equilibrium points are obtained. We have made some observations about the behavior of the solutions near the origin, and their dependence on initial conditions. These observations also establish a difference from other models, where the origin is usually a saddle. A bifurcation analysis has been carried out treating the effort as the bifurcation parameter. Optimal harvest policy (equilibrium level) has been discussed using effort as a control parameter. It has been found that the total user cost of harvest per unit effort equals the discounted value of the future price. It has also been noted that an infinite discount rate yields zero economic rent (then the fishes remain unexploited). Dynamic optimization of the harvest policy is also carried out taking effort as a dynamic variable. Finally the optimal singular trajectory is derived using Pontryagin’s Maximum Principle. Acknowledgments The authors are grateful to the anonymous referees and Professor T. E. Simos (Editor-in-Chief) for their careful reading, constructive comments and helpful suggestions, which have helped to improve the presentation of this work significantly. References [1] V. Ajraldi, M. Pittavino, E. Venturino, Modelling herd behavior in population systems, Nonlinear Anal. RWA 12 (2011) 2319–2338. [2] M. Banerjee, B.W. Kooi, E. Venturino, An ecoepidemic model with prey herd behavior and predator feeding saturation response on both healthy and diseased prey, Math. Model. Nat. Phenom. 12 (2017) 133–161. [3] S.P. Bera, A. Maiti, G.P. Samanta, Modelling herd behavior of prey: analysis of a prey-predator model, World J. Model. Simul. 11 (2015) 3–14. [4] S.P. Bera, A. Maiti, G.P. Samanta, Stochastic analysis of a prey-predator model with herd behaviour of prey, Nonlinear Anal. Model. Control 2 (2016) 345–361. [5] A.M. Boustany, R. Matteson, M. Castleton, C. Farwell, B.A. Block, Movements of pacific bluefin tuna (thunnus orientalis) in the eastern north pacific revealed with archival tags, Prog. Oceanogr. 86 (2010) 94–104. [6] J. Brattey, Biological characteristics of atlantic cod (gadus morhua) from three inshore areas of northeastern newfoundland, NAFO Sci. Coun. Stud. 29 (1997) 31–42. [7] A.P. Braza, Predator-prey dynamics with square root functional responses, Nonlinear Anal. RWA 13 (2012) 1837–1843. [8] E.D. Brown, J.H. Churnside, R.L. Collins, T. Veenstra, J.J. Wilson, K. Abnett, Remote sensing of capelin and other biological features in the north pacific using lidar and video technology, ICES J. Mar. Sci. 59 (2002) 1120–1130. [9] J. Chattopadhyay, S. Chatterjee, E. Venturino, Patchy agglomeration as a transition from monospecies to recurrent plankton blooms, J. Theor. Biol. 253 (2008) 289–295. [10] C.W. Clark, Bioeconomic Modelling and Fisheries Management, John Wiley & Sons, New York, 1985. [11] C.W. Clark, Mathematical Bioeconomic: The Optimal Mmanagement of Renewable Resources, John Wiley & Sons, New York, 1990. [12] C. Cosner, D.L. DeAngelis, J.S. Ault, D.E. Olson, Effects of spatial grouping on the functional response of predators, Theor. Pop. Biol. 56 (1999) 65–75. [13] M. D’Elia, B. Patti, A. Sulli, G. Tranchida, A. Bonanno, G. Basilone, G. Giacalone, F. I., S. Genovese, C. Guisande, S. Mazzola, Distribution and spatial structure of pelagic fish schools in relation to the nature of the seabed in the sicily straits (central mediterranean), Mar. Ecol. 30 (2009) 151–160. [14] O. Flaaten, On the bioeconomics of predator and prey fishing, Fish. Res. 37 (1998) 179–191. [15] B.S. Goh, G. Leitmann, T.L. Vincent, Optimal control of a prey–predator system, Math. Biosci. 19 (1974) 263–286. [16] G. Leitmann, An Introduction to Optimal Control, Mc Graw-Hill, New York, 1966. [17] A. Lotka, Elements of Physical Biology, Williams and Wilkins, Baltimore, 1925. [18] A. Maiti, B. Patra, G.P. Samanta, Bionomic exploitation of a ratio-dependent predatorprey system, Int. J. Math. Ed. Sci. Tech. 39 (2008) 1061–1076. [19] A. Maiti, P. Sen, D. Manna, G.P. Samanta, A predator-prey system with herd behaviour and strong allee effect, Nonlinear Dyn. Syst. Theory 16 (2016) 86–101. [20] A. Maiti, P. Sen, G.P. Samanta, Deterministic and stochastic analysis of a prey-predator model with herd behaviour in both, Syst. Sci. Control Eng. 4 (2016) 259–269. [21] D. Manna, A. Maiti, G.P. Samanta, Analysis of a harvested predator-prey system with schooling behaviour, Int. J. Dyn. Control (2017), doi:10.1007/ s40435- 017- 0321- y. [22] A. Marzano, Harvesting the Sea: The Exploitation of Marine Resources in the Roman Mediterranean, Oxford University Press, Oxford, 2013. [23] A.O. Misund, J.C. Coetzee, P. Freon, M. Gardener, K. Olsen, I. Svellingen, I. Hampton, Schooling behaviour of sardine sardinops sagas in false bay, South Africa, S. Afr. J. Mar. Sci. 25 (2003) 185–193.

48

D. Manna et al. / Applied Mathematics and Computation 317 (2018) 35–48

[24] Y. Mitsunaga, C. Endo, R.P. Babaran, Schooling behavior of juvenile yellowfin tuna thunnus albacares around a fish aggregating device (FAD) in the philippines, Aquat. Living Resour. 26 (2013) 79–84. [25] B. Partridge, T. Pitcher, M. Cullen, J. Wilson, The three-dimensional structure of fish schools, Behav. Ecol. Sociobiol. 6 (1980) 277–288. [26] L.S. Pontryagin, V.S. Boltyanskii, R.V. Gamkrelidze, The Mathematical Theory of Optimal Processes, Wiley-Interscience, 1962. [27] G.P. Samanta, D. Manna, A. Maiti, Bioeconomic modeling of a three-species fishery with switching effect, J. Appl. Math. Comput. 12 (2003) 219–231. [28] P.F. Verhulst, Notice sur la loi que la population poursuit dans son accroissement, Corresp. Math. Phys. 10 (1838) 113–121. [29] V. Volterra, Variazioni e fluttuazioni del numero di individui in specie animali conviventi, Mem. Accd. Lincei. 2 (1926) 31–113. [30] C. Xu, S. Yuan, Stability and hopf bifurcation in a delayed predator-prey system with herd behavior, Abstr. Appl. Anal., 2014, 2014. 1–8. http://dx.doi. org/10.1155/2014/568943.