Dynamics of populations in fish farm: Analysis of stability and direction of Hopf-bifurcating periodic oscillation

Applied Mathematical Modelling 36 (2012) 2118–2127

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Dynamics of populations in fish farm: Analysis of stability and direction of Hopf-bifurcating periodic oscillation Nurul Huda Gazi Department of Mathematics, Aliah University, DN-41, Sector V, Salt Lake, Kolkata 700 091, India

a r t i c l e

i n f o

Article history: Received 14 October 2009 Received in revised form 21 July 2011 Accepted 9 August 2011 Available online 19 August 2011 Keywords: Fish farm Nonlinear differential equation model Stability Time delay Hopf-bifurcation

a b s t r a c t The paper deals with the dynamical behavior of fish and mussel population in a fish farm where external food is supplied. The ecosystem of the fish farm is represented by a set of nonlinear differential equations involving nutrient (food), fish and mussel. We have listed some results already obtained. We have analyzed for the direction of Hopf-bifurcation, stability of the Hopf-bifurcating periodic orbits, and the period of the periodic orbits by using Poincare’ normal form and center manifold theory. We have performed numerical simulation to support the analytical results. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction The industrial wastes, pesticides and fertilizer round-off and fecal matters are the pollutants which can deteriorate water body affecting the wildlife and ecosystem surrounding a water body. With the believe that water would dilute any substance, the lakes and any water body have become destinations of dumping grounds of pollutants. The water pollution is a concern for fish farm. In many fish farms the fish depends on natural food and also on the commercial feeds. The external food is sometimes a threat to the fish farm. The heavy organic load in the fish farm causes depletion of dissolved oxygen in the water column. The water current speed is sometimes decreased for the uncontrolled proliferation of milk fish cages which causes the depletion of the dissolved oxygen in the water column [1]. The harmful algal bloom in the open coastal water affects severely on the aquaculture. The shellfish, seaweeds and oyster are used as bio-filter. These species usually used as bio-filter which can consume excess feed not dissolved into water. These, in turn, can enhance the water quality. The artificial external food can cause rapid growth of the fish population. This, again, may have negative effect when it is overfed leading to a low oxygen condition. The shellfish plays important role in controlling the water quality. The characteristic of the shellfish is that they are filter feeders, feeding on bacteria and nutrients. The Barnacles, mussels and other shellfish are used in an innovative way so that they behave like a bio-filtering system to clean up fish farms [2]. Mussels can help in pollution problem in two ways: mussels can digest particles of feed waste and fish excreta and also consume phytoplankton which also thrive on inorganic nutrients such as phosphorous and nitrogen. The present work is an extension of the earlier work by Gazi et al. [3]. We have considered the model comprising fish and mussel. The external food is supplied constantly. The incorporation of gestational delay reveals oscillation of the system from stable state of the equilibria [3]. In many works the delay has been incorporated leading to delay differential equation. Delay

E-mail address: [email protected] 0307-904X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.08.009

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

2119

differential equation mathematical models of population biology have been analyzed in Cushing [4], Gopalsamy [5], Kuang [6], May [7] and Beretta [8], and the references cited therein. In this work we have mainly studied the stability of the periodic oscillation of different species. The Hopf-type bifurcating oscillation is identified and we have studied the direction and period of the oscillation using center manifold theory and Poincare’ normal form [9]. We now consider the following model system representing a lake ecosystem with two components consisting food (x) and fish (y):

dx ¼ K  ða þ byÞx; dt dy ¼ ½L þ My  Nxðt  sÞy: dt

ð1:1Þ

Here K, a, b, represent the external food, outflow or sedimentation rate of nutrient, nutrient uptake rate of fish population respectively. L, M and N are respectively the death rate, the intra-specific competition, and the proportion of nutrient contributing to the biomass of fish. s is the delay due to gestation of the fish for consumption of food. For higher level of external food (K), the fish population can not assimilate it, the excess of the food makes the fish farm eutrophic. It has great impact on the ecosystem functioning. At this stage, an another species feeding on the by-product in the fish farm could balance the ecosystem where both the species co-exist. The inclusion of mussel species is traditional way to control eutrophication in any water body. The growth equation of the mussel population (z) may be taken as

dz ¼ ða1  b1 xÞz; dt

ð1:2Þ

where a1 and b1 are the death rate and the conversion efficiency of the mussel population. Thus the external feed not assimilated into the fish biomass could reach mussel in the form of particulate organic matter and mussel can easily consume it. So, we consider the following set of differential equations representing the fish farm with three components: nutrient, fish and mussel

dx ¼ K  ða þ by þ czÞx; dt dy ¼ ½L þ My  Nxðt  sÞy; dt dz ¼ ða1  b1 xÞz; dt

ð1:3Þ

with initial densities x(h) = n(h) where nðhÞ 2 Cð½s; 0; Rþ Þ and n(h) P 0, n(0) P 0, y(0) > 0 and z(0) > 0. Here c is the maximum nutrient uptake rate of the mussel and x(t), y(t) and z(t) are the densities of nutrient, fish and mussel populations biomass at time ‘t’ and N < b, b1 < c. We want to study the delay differential equation (1.3) for its dynamical behavior. We organize the paper as follows: In Section 2, we have listed the results obtained in the earlier work of Gazi et al. [3]. In Section 3, we have analyzed the direction and stability of the Hopf bifurcation. The periodicity of the periodic oscillation is calculated. In the last conclusion section we make critical analysis of the analytical results. 2. Steady states and stability results For the above model system (1.3), following are the equilibrium states: 1. E1 ðx1 ; 0; 0Þ; x1 ¼ Ka .

    qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   M bL 2 bN 2. E2 ðx2 ; y2 ; 0Þ; x2 ¼ 2bN  a  bL þ 4K þ a  ; y2 ¼ M1 ½Nx2  L. M M M   3. E3 ðx3 ; 0; z3 Þ; x3 ¼ ba11 ; z3 ¼ 1c xK3  a .   4. E ðx ; y ; z Þ; x ¼ ba11 ; y ¼ M1 ðNx  LÞ; z ¼ 1c xK  a  by . Let

K1 ¼ a

L ; N

K2 ¼ a

a1 ; b1



 bL K 2 K3 ¼ K2 1 þ 1 : aM K 1

ð2:1Þ

The feasibility of the equilibria are discussed in Section 3 of Gazi et al. [3]. Lemma 2.1. All solutions of the system (1.3) that initiate in R3þ are eventually uniformly bounded within the region B, where B ¼ fðx; y; zÞ 2 R3þ : 0 6 x þ y þ z 6 K=minfa; L; a1 g þ , for any  > 0}. Proof. We assume that the right-hand sides of system of Eq. (3.1) are smooth functions of (x(t), y(t), z(t)) of tðt 2 Rþ Þ. Let x(t), y(t) and z(t) be any solution at time ‘t’ with initial condition x(0) > 0, y(0) > 0, z(0) > 0.

2120

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

From the first equation of (1.3) we have dx < K  ax. Thus, according to the standard comparison theorem x 6 K for t P 0. dt Let us now construct a time dependent function W(t) = x(t) + y(t) + z(t) of ‘t’. The time derivative of W(t) along the solution of model system (1.3) is

dW ¼ K  ða þ by þ czÞx  ðL þ My  NxÞy  ða1  b1 xÞz: dt Since N < b, b1 < c, the above expression reduces to

dW < K  ðax þ Ly þ a1 zÞ < K  PW; dt where P = min{a, L, a1}. Thus

dW þ PW 6 K: dt Applying a theorem in differential inequalities [10] we obtain

0 6 Wðx; y; zÞ 6

K þ Wðxð0Þ; yð0Þ; zð0ÞÞePt ; P

ð2:2Þ

and for t ! þ1; 0 6 Wðx; y; zÞ 6 KP . Thus the model system is dissipative. The asymptotic bound of the system is K/P. Therefore, all solutions of system (3.1) initiated at (x(0), y(0), z(0)) enter into the compact region B ¼ fðx; y; zÞ 2 R3þ : 0 < x þ y þ z 6 K=P þ , for any  > 0}. Thus all solutions of system (1.3) is uniformly bounded initiated at (x(0), y(0), z(0)). This completes the proof. h The model system (1.3) with s = 0 is locally asymptotically stable if the external food supply K attains different level at different equilibria: at E1: K < K1,K < K2; at E2: K < K3; at E3: K2 < K1; and at E⁄: Always [3]. According to Gazi et al. [3] the model system (1.3) with s = 0 is globally asymptotically stable around E2 and E⁄ (Section 3.2 in [3]). Theorem 2.1 (According to Gazi et al. [3]). The model system (1.3) is locally asymptotically stable about E⁄ if the external food supply K satisfies the following: (i)K1 < K2, (ii)K> K3 and4 h   i (iii) xK þ My > 4 2 xK þ My b1 cMx y z þ ðbNx y Þ2 . 2.1. Oscillatory behavior The model system (1.3) has equilibrium state E⁄ about which we want to study the oscillatory behavior. The study had been carried out in Gazi et al. [3]. Here we need some results of the previous work for further study of the system (1.3). Here the Jacobian matrix about E⁄ of the linearized system of (1.3) is

0

1 ða þ by þ cz Þ bx cx B C UðE ðx ; y ; z ÞÞ ¼ @ Ny ðL þ 2My  Nx Þ 0 A: b1 z 0 ða1  b1 x Þ

ð2:3Þ

We derive the characteristic equation about E⁄(x⁄, y⁄, z⁄) as under

Dðk; sÞ  k3 þ L1 k2 þ L2 k þ L3 keks þ L4 ¼ 0;

ð2:4Þ

where L1 = K/x⁄ + My⁄, L2 = KM y⁄/x⁄ + b1c x⁄z⁄, L3 = bNx⁄y⁄ and L4 = Mx⁄y⁄z⁄. Let s – 0, then from (2.4), if k = a + ix is a root, then for periodic solution the characteristic Eq. (2.4) has a purely imaginary root. If s be the bifurcation parameter and s = s0 be the value of s at which the stable solution becomes unstable, then a(s0) = 0 and x(s0) = x0 – 0 [9]. Putting k = ix0 in the Eq. (2.4), we get

  Dðix0 ; s0 Þ ¼ L1 x20 þ L3 x0 sin x0 s0 þ L4 þ i x30 þ L2 x0 þ L3 x0 cos x0 s0 ¼ 0:

ð2:5Þ

From the above equation we get two equations

 L1 x20 þ L3 x0 sin x0 s0 þ L4 ¼ 0;  x30 þ L2 x0 þ L3 cos x0 s0 ¼ 0:

ð2:6Þ

From these two equations we arrive at









x60 þ L21  2L2 x40 þ L22  2L1 L4  L23 x20 þ L24 ¼ 0:

ð2:7Þ

2121

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

Thus one of the values of x0 is given by (2.7). The corresponding delay parameter value s0j is given (from (2.6)) by

s0j ¼

1

x0

arctan



 L1 x20  L4 jp þ ; x30  L2 x0 x0

j ¼ 0; 1; 2; . . . :

ð2:8Þ

So, j = 0 in (2.8) gives

s0 ¼

1

x0

arctan



 L1 x20  L4 ; x30  L2 x0

ð2:9Þ

which is the smallest time delay for which the stable solution becomes unstable. Thus the characteristic Eq. (2.4) has a pair of purely imaginary roots for x = x0 and s = s0. So, we are to verify transversality condition ddas jðix0 ;s0 Þ – 0. Here ddas jðix0 ;s0 Þ – 0 (For details see [3]). ddas jðix0 ;s0 Þ is given by

   i da x2 h ¼ 20 3x40 þ 2 L21  2L2 x20 þ L22  2L1 L4  L23 ;  ds ðix0 ;s0 Þ D

ð2:10Þ

where D2 ¼ D21 þ D22 . D1 and D2 are given as

  D1 ¼ 3x20 þ L2 cos x0 s0  2L1 x0 sin x0 s0 þ L3 ;   D2 ¼ 3x20 þ L2 sin x0 s0 þ 2L1 x0 cos x0 s0  L3 x0 s0 :

ð2:11Þ

Theorem 2.2 (According to [3]). (i) If K1 < K2, K > K3 and (K/x⁄ + My⁄)4 > 8(K/x⁄ + My⁄) M b1c x⁄y⁄z⁄ + 4(bN x⁄y⁄)2, then the equilibrium E⁄ is asymptotically stable for all s P 0. If the above is satisfied then E⁄ is locally asymptotically stable for s 2 [0, s0). (ii) If K1 < K2, K > K3 and (KMy⁄/x⁄ + b1c x⁄z⁄)2 < 2(K/x⁄ + M y⁄) M b1c  x⁄y⁄z⁄ + (bN x⁄y⁄)2, then model system (1.3) is unstable for s P s0. Hopf-bifurcation occurs at s = s0. Thus the stable solution bifurcates from E⁄ into periodic oscillations as s passes through s0. In the next section we will analyze the stability of the periodic oscillation. 3. Stability and direction of periodic solution In this section we want to study the stability, direction and the period of the Hopf-bifurcating periodic solution obtained in Gazi et al. [3]. We have used the Poincare’ normal form and the center manifold theory to determine the above-mentioned qualitative behavior about E⁄ [9]. Same method has been used in [11–13]. We have obtained k = ± ix0 at the critical parameter value s = s0. We use this value of x and s in this section.  i ¼ ui ðssÞ; s ¼ s0 þ l, then dropping the bars (we have replaced s by t) We take t ¼ ss; u1 ¼ x  x4 ; u2 ¼ y  y4 ; u3 ¼ z  z4 ; u we get,



K u_ 1 ðtÞ ¼ ðs0 þ lÞ  u1  bx4 u2  cx4 u3  bu1 u2  cu1 u3 ; x4   u_ 2 ðtÞ ¼ ðs0 þ lÞ Ny4 u1 ðt  sÞ  My4 u2 þ Nu1 ðt  sÞu2  Mu22 ;

ð3:1Þ

u_ 3 ðtÞ ¼ ðs0 þ lÞðb1 z4 u1 þ b1 u1 u3 Þ: The system (3.1) is transformed into a functional differential equation in C = C([ 1, 0], R3) as

u_ ¼ Ll ðut Þ þ Fðl; ut Þ;

ð3:2Þ

where u(t) = (u1(t), u2(t), u3(t))T 2 R3, and Ll: C ? R, F : R  C ? R are given by

0 B Ll ð/Þ ¼ ðs0 þ lÞ@

K=x4

bx4

0

My4

b 1 z4

0

cx4

10

/1 ð0Þ

1

0

0

CB B C 0 A@ /2 ð0Þ A þ ðs0 þ lÞ@ Ny4 /3 ð0Þ 0 0

0 0

10

/1 ð1Þ

1

CB C 0 0 A@ /2 ð1Þ A; 0 0

ð3:3Þ

/3 ð1Þ

and

0 B Fðl; /Þ ¼ @

b/1 ð0Þ/2 ð0Þ  c/1 ð0Þ/3 ð0Þ N/1 ð1Þ/2 ð0Þ  M/22 ð0Þ

1 C A;

ð3:4Þ

b1 /1 ð0Þ/3 ð0Þ where /(h) = (/1(h),/2(h), /3(h))T 2 C. Here Ll is a one parameter family of bounded linear operators in C[ 1, 0] ? R3.

2122

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions g(h, l) for 2 h 2 ½1; 0 ! R3 , such that

Ll / ¼

Z

0

/ðhÞdgðh; lÞ for / 2 C:

ð3:5Þ

1

For this representation we can choose g(h, l) as follows:

0

 xK4

gðh; lÞ ¼ ðs0 þ lÞB @ 0

bx4

cx4

0

0

B C 0 AdðhÞ  ðs0 þ lÞ@ Ny4 0 0

My4

b 1 z4

1

0

0 0

1

C 0 0 Adðh þ 1Þ;

ð3:6Þ

0 0

where d is the Dirac delta function. For / 2 C1[ 1, 0], we define

( d/

1 6 h < 0;

; dh R0

AðlÞ/ ¼

/ðsÞdgðs; lÞ; h ¼ 0: 1

and

RðlÞ/ ¼



0;

1 6 h < 0;

Fðl; /Þ; h ¼ 0:

Thus the system (3.2) is equivalent to

u_ t ¼ AðlÞut þ RðlÞut ;

ð3:7Þ

where ut(h) = u(t + h) for 1 6 h 6 0. System (3.7) is more mathematically pleasing one as this equation involves a single unknown ut rather than both ut and u. Here the operator A depends on the parameter l. When l = 0, the Hopf bifurcation occurs. We define an adjoint operator A⁄ of A as follows: For w 2 C1([0, 1], (R3)⁄), ((R3)⁄ is the three dimensional vector space of row vectors).

( 

A wðsÞ ¼

 dwðsÞ ; 0 < s 6 1; ds R0 T wðtÞdg ðt; 0Þ; s ¼ 0: 1

For / 2 C([ 1, 0]), C2) and w 2 C([0, 1], C2), define a bilinear form

 < wðsÞ; /ðhÞ >¼ wð0Þ/ð0Þ 

Z

0

1

Z

h

  hÞdgðhÞ/ðnÞdn; v h wðn

ð3:8Þ

n¼0

where g(h) = g(h, 0). Then A(0) and A⁄(0) are adjoint operators. The domains of A(0) and A⁄(0) are C1[ 1, 0] and C1[0, 1] respectively. We know that ± ix0s0 are the eigenvalues of A(0). They are also eigenvalues of A⁄. We will compute the eigenvectors of A(0) and A⁄(0) corresponding to ix0s0 and ix0s0 respectively. The eigenvectors are given by

Að0ÞqðhÞ ¼ ix0 s0 qðhÞ; ix0 s0

4e 1 z4 then qðhÞ ¼ ð1; a; bÞT eix0 s0 h , where a ¼ Ny ; b ¼ bix . If q⁄ = D(1, a⁄, b⁄) be the eigenvector of A⁄(0) for  ix0s0, then ix0 þMy 0 4







A ð0Þq ðsÞ ¼ x0 s0 q ðsÞ; ⁄ 4 4 where a ¼ ix0bx and b ¼ icx þMy x0 . To find D we use (3.8), and since < q (s), q(h) > = 1 we get, 4

D¼

1 :   þ s0 a   þ bb   Ny4 eix0 s0 1 þ aa

In the remainder of the section we compute the local coordinates to describe the center manifold C0 at l = 0. Let ut be the solution of (3.7), when l = 0. We change the variable as follows

zðtÞ ¼< q ; ut >

ð3:9Þ

and also we define (on the center manifold C0),

ðhÞ ¼ ut ðhÞ  2RefzðtÞqðhÞg: Wðt; hÞ ¼ ut ðhÞ  zðtÞqðhÞ  zðtÞq On the center manifold C0,

Wðt; hÞ ¼ WðzðtÞ; zðtÞ; hÞ;

ð3:10Þ

2123

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

where

Wðz; z; hÞ ¼ W 20 ðhÞ

z2 z2 z3 þ W 11 ðhÞzz þ W 02 ðhÞ þ W 30 ðhÞ þ . . . 2 2 6

ð3:11Þ

 . where z and z are respectively local coordinates for the center manifold C0 in the direction of q⁄ and q For the real solution ut 2 C0 of (3.7), we have

 ð0ÞFð0; Wðz; z; 0Þ; 2RefzqðhÞgÞ ¼ ix0 s0 z þ q  F 0 ðz; zÞ: z_ ¼ ix0 s0 z þ q

ð3:12Þ

The above, we write in abbreviated form

z_ ¼ ix0 s0 z þ gðz; zÞ;

ð3:13Þ

where

 ð0ÞF 0 ðz; zÞ ¼ g 20 ðhÞ gðz; zÞ ¼ q

z2 z2 z2 z þ ... þ g 11 ðhÞzz þ g 02 ðhÞ þ g 21 ðhÞ 2 2 2

ð3:14Þ

and also from (3.10), we have

_ ¼ u_ t  z_ q  z_ q ¼ W



 ð0ÞF 0 qðhÞg; AW  2Refq h 2 ½1; 0Þ;  ð0ÞF 0 qðhÞg þ F 0 ; h ¼ 0: AW  2Refq

This we write in abbreviated form

_ ffi AW þ Hðz; z; hÞ; W

ð3:15Þ

where

Hðz; z; hÞ ¼ H20 ðhÞ

z2 z2 z2 z þ ... þ H11 ðhÞzz þ H02 ðhÞ þ H21 ðhÞ 2 2 2

ð3:16Þ

ðhÞ and qðhÞ ¼ ð1; a; bÞT eix0 s0 h , therefore, From (3.9) and (3.11), we have ut ¼ ðu1t ðhÞ; u2t ðhÞ; u3t ðhÞÞ ¼ Wðt; hÞ þ zqðhÞ þ zq

z2 z2 ð1Þ ð1Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ z þ z þ oðjðz; zÞj3 Þ; 2 2 z2 z2 ð2Þ ð2Þ ð2Þ z þ oðjðz; zÞj3 Þ; u2t ð0Þ ¼ W 20 ð0Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ az þ a 2 2 z2 z2 ð3Þ ð3Þ ð3Þ z þ oðjðz; zÞj3 Þ; u3t ð0Þ ¼ W 20 ð0Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ bz þ b 2 2 z2 z2 ð1Þ ð1Þ ð1Þ u1t ð1Þ ¼ W 20 ð1Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ zeix0 s0 þ zeix0 s0 þ oðjðz; zÞj3 Þ; 2 2 z2 z2 ð2Þ ð2Þ ð2Þ zeix0 s0 þ oðjðz; zÞj3 Þ; u2t ð1Þ ¼ W 20 ð1Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ azeix0 s0 þ a 2 2 z2 z2 ð3Þ ð3Þ ð3Þ zeix0 s0 þ oðjðz; zÞj3 Þ: u3t ð1Þ ¼ W 20 ð1Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ bzeix0 s0 þ b 2 2 ð1Þ

u1t ð0Þ ¼ W 20 ð0Þ

From the definition of F(l, ut), we have

0

B gðz; zÞ ¼ q ð0ÞF 0 ðz; zÞ ¼ s0 Dð1; a ; b Þ@ 





bu1t ð0Þu2t ð0Þ  cu1t ð0Þu3t ð0Þ Nu1t ð1Þu2t ð0Þ 

Mu22t ð0Þ

1 C A:

b1 u1t ð0Þu3t ð0Þ Using the above in this expression and then comparing the coefficients with (3.14), we get

 ; g   ðNaeix0 s0  M a2 Þ  b1 bb g 20 ¼ 2s0 D½ðba þ cbÞ  a 11  ; g  Þ  b1 Refbgb   ðNða  eix0 s0 þ aeix0 s0 Þ  2Maa ¼ s0 D½ðbRefag þ cRefbgÞ  a 02 b  ; g  a  þ cbÞ   ðNa  eix0 s0  Ma  2 Þ  b1 b ¼ 2s0 D½ðba 21 n o n o ð1Þ ð2Þ ð2Þ  ð1Þ ð0Þ þ 2bW ð1Þ ð0Þ þ 2W ð3Þ ð0Þ þ W ð3Þ ð0Þ  W ð1Þ ¼ s0 Db a ð0Þ þ 2 a W ð0Þ þ 2W ð0Þ þ W ð0Þ  s0 Dc bW 20 11 11 20 20 11 11 20 n o n o ð1Þ ð2Þ ð2Þ ð2Þ   N 2aW 11   M aW 11  W ð2Þ þ s0 D a ð1Þ þ 2W 11 ð0Þeix0 s0 þ 2W 20 ð0Þeix0 s0  2s0 Da ð0Þ þ a 20 ð0Þ h i  ð1Þ ð0Þ þ 2bW ð1Þ ð0Þ þ 2W ð3Þ ð0Þ þ W ð3Þ ð0Þ :  b1 bW þ s0 D b 20 11 11 20 To compute g21 we need to calculate W20(h) and W11(h). On the center manifold C0 near the origin

_ ¼ W z z_ þ W z z_ : W

ð3:17Þ

2124

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

_ By We use (3.11) for Wz and W z and equations (3.13) and (3.14), to replace z_ and z_ to obtain the second expression for W. comparison of this result with (3.15) we can derive equations for the coefficients Wij(h). These are

ð2ix0 s0  AÞW 20 ðhÞ ¼ H20 ðhÞ; AW 11 ðhÞ ¼ H11 ðhÞ:

ð3:18Þ

From (3.15), we know that for h 2 [ 1, 0),

ðhÞ ¼ gðz; zÞqðhÞ  gðz; zÞq ðhÞ:  ð0ÞF 0 qðhÞ  q ð0ÞF 0 q Hðz; z; hÞ ¼ q

ð3:19Þ

Comparing the coefficients with (3.16), we get

ðhÞ H20 ðhÞ ¼ g 20 qðhÞ  g02 q ðhÞ: and H11 ðhÞ ¼ g 11 qðhÞ  g11 q

ð3:20Þ

From (3.18) and (3.20), we get

_ 20 ðhÞ ¼ 2ix0 s0 W 20 ðhÞ þ g qðhÞ þ g02 q ðhÞ: W 20

ð3:21Þ

Therefore,

W 20 ðhÞ ¼

ig 20

x0 s0

qð0Þeix0 s0 h þ

ig02 ð0Þeix0 s0 h þ C 1 e2ix0 s0 h : q 3x0 s 0

ð3:22Þ

Also, from (3.18) and (3.20), we get

ðhÞ AW 11 ðhÞ ¼ g 11 qðhÞ þ g11 q or,

_ 11 ¼ g qð0Þeix0 s0 h þ g11 q ð0Þeix0 s0 h : W 11 Therefore,

W 11 ðhÞ ¼

ig 11

x0 s 0

qð0Þeix0 s0 h þ

ig11

x0 s 0

ð0Þeix0 s0 h þ C 2 : q

ð3:23Þ

Here C1 and C2 arbitrary constant vectors. From the definition of A and (3.18), we obtain

Z

0

W 20 ðhÞdgðhÞ ¼ 2ix0 s0 W 20 ðhÞ  H20 ð0Þ;

ð3:24Þ

W 11 ðhÞdgðhÞ ¼ H11 ð0Þ where gðhÞ ¼ gð0; hÞ:

ð3:25Þ

1

and

Z

0

1

From (3.15) and (3.16) for h = 0, we know that

ð0Þ þ F 0 ¼ gðz; zÞqð0Þ  gðz; zÞq  ð0ÞF 0 qð0Þ  q ð0ÞF 0 q ð0Þ þ F 0 Hðz; z; 0Þ ¼ q     2 2 2 z z2 z z ð0Þ þ F 0 : ¼  g 20 ð0Þ þ g 11 ð0Þzz þ g 02 ð0Þ qð0Þ  g20 ð0Þ þ g11 ð0Þzz þ g02 ð0Þ q 2 2 2 2 Comparing the coefficients of similar powers of z, we get

0

1 ba  cb B ð0Þ þ 2s0 @ Naeix0 s0  M C H20 ð0Þ ¼ g 20 qð0Þ  g02 q A; b1 b

ð3:26Þ

and

0

1 bRefag þ cRefbg B ð0Þ þ 2s0 @ Nðaeix0 s0 þ a C  eix0 s0 Þ  2Maa H11 ð0Þ ¼ g 11 qð0Þ  g11 q A: b1 Refbg Substituting (3.22) and (3.26) into (3.24), we get

Z 2ix0 s0 I 

0

1

0 1 ba  cb

B C e2ix0 s0 dgðhÞ C 1 ¼ 2s0 @ Naeix0 s0  M A b1 b

ð3:27Þ

2125

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

Or,

0

2ix0 s0 þ xK4

B @ Ny4 e2ix0 s0

bx4 My4 þ 2ix0 s0

b1 z4

10 ð1Þ 1 0 1 C1 ba  cb B C C B C ð2Þ C ix0 s0  M A: 0 AB @ C 1 A ¼ @ Nae ð3Þ b1 b 2ix0 s0 C1 cx4

0

ð3:28Þ

From (3.28) we get ð1Þ

1 ½2ix0 s0 ðba  cbÞðMy4 þ 2ix0 s0 Þ  2ix0 s0 bx4 ðN aeix0 s0  MÞ  b1 cbx4 ðMy4 þ 2ix0 s0 Þ; D1    1 ¼ 2ix0 s0 ðba  cbÞNy4 e2ix0 s0 þ ðNaeix0 s0  MÞ 4x20 s20 þ b1 cx4 z4  b1 cNbx4 y4 e2ix0 s0 ; D1 1 ¼ ½b1 z ðba  cbÞðMy4 þ 2ix0 s0 Þ  bb1 x z ðNaeix0 s0  MÞ þ b1 bf2ix0 s0 ðMy þ 2ix0 s0 Þ þ bNx y eix0 s0 g; D1

C1 ¼ ð2Þ

C1

ð3Þ

C1 where

  D1 ¼  4x20 s20  b1 cx z ðMy4 þ 2ix0 s0 Þ þ 2ibx x0 s0 e2ix0 s0 : Substituting (3.23) and (3.27) into (3.25), we get

0

K x4

bx4

B @ Ny4

My4

b1 z4

0

1 10 ð1Þ 1 0 C2 bRefag þ cRefbg C B CB ð2Þ C ix0 s0 C  eix0 s0 Þ  2M aa þa A: 0 AB @ C 2 A ¼ @ Nðae ð3Þ b1 Refbg 0 C2

cx4

ð1Þ

ð3:29Þ

ð2Þ

 =Mz and  eix0 s0 Þ  2M aa Therefore, C 2 ¼ Refbg=z ; C 2 ¼ ½Nðaeix0 s0 þ a ð3Þ   þ b1 RefbgðKMy =x  eix0 s0 Þ  2M aa C 2 ¼ ½b1 My z ðbRefag þ cRefbgÞ  bb1 x z ½Nðaeix0 s0 þ a

þ bNx y Þ=ðb1 cMx y z Þ:

b 5

6

4

Fish density

Nutrient density

a 7

5 4 3

2 1

2 1

3

0

50

100

150

0

200

0

50

Time (t)

100

150

200

Time (t)

d

c 10 8

Mussel

Mussel density

10

6

5

0 5

4

10 2

0

50

100

Time (t)

150

200

Fis

5

h

0

0

Nutrient

Fig. 1. Stable time series evolution of the species of model system (1.3) with parameter values K = 200, a = 3.0, b = 25, c = 9.0, L = 12.0, M = 4.0, N = 8.0, a1 = 4, b1 = 1.75, s = 1.2.

2126

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

From (3.22) and (3.23), W20(h) and W11(h) are known. Their expressions are given as follows:

ig02  þ C 1ð1Þ ; q 3x0 s 0 ig ig02 ð2Þ a þ C 1ð2Þ ; W 20 ð0Þ ¼ 20 a þ x0 s 0 3x0 s0 ig ig02  ð3Þ ð3Þ b þ C1 ; W 20 ð0Þ ¼ 20 b þ x0 s 0 3x0 s0 ig 11 ig11 ð1Þ ð1Þ þ þ C2 ; W 11 ð0Þ ¼ ig 20

ð1Þ

W 20 ð0Þ ¼

x0 s 0

x0 s 0

ð2Þ W 11 ð0Þ

¼

ð3Þ

W 11 ð0Þ ¼

ig 11

x0 s 0 ig 11

x0 s 0

þ

x0 s0

aþ bþ

ig11

x0 s0

a ð0Þ þ C 2ð2Þ ;

ig11  ð3Þ bð0Þ þ C 2 :

x0 s 0

Therefore, g21 is completely known. So we can find the first Lyapunov coefficient c1(0) given by

c1 ð0Þ ¼

i

jg j2 g 20 g 11  2jg 11 j  02 3 2

2x0 s0

!

þ

g 21 : 2

ð3:30Þ

Therefore, the coefficients b2, l2 and T2 are also known which are given as follows:

b2 ¼ 2Refc1 ð0Þg; Refc1 ð0Þg l2 ¼  ; Refk0 ðs0 Þg Imfc1 ð0Þg þ l2 Imfk0 ðs0 Þg ; T2 ¼ 

ð3:31Þ

x0 s0

b 5

6

4

5

Fish density

Nutrient density

a 7

4 3

2 1

2 1

3

0

50

100 Time (t)

150

0

200

0

50

100 Time (t)

c

150

200

d

10 10

6

Mussel

Mussel density

8

4

0 5

2 0

5

10 0

50

100 Time (t)

150

200

Fis

5

h

0

0

Nutrient

Fig. 2. Unstable oscillation of the species of model system (1.3) with parameter values K = 200, a = 3.0, b = 25, c = 9.0, L = 12.0, M = 4.0, N = 8.0, a1 = 4, b1 = 1.75, s > s0.

N.H. Gazi / Applied Mathematical Modelling 36 (2012) 2118–2127

2127

which determine the qualities of bifurcating periodic solution in the center manifold at the critical value s0. Here, b2 determines the stability of the bifurcating periodic solution: the bifurcating periodic solution is stable (unstable) if b2 < 0 (b2 > 0); l2 determines the direction of the Hopf bifurcation: if l2 > 0 (l2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solution exists for s > s0(s < s0). T2 determines the period of the bifurcating periodic solution: the period increases (decreases) if T2 > 0 (T2 < 0). 4. General discussion In this paper we have considered the deterministic model of a fish farm with three interacting components: food (nutrient), fish and mussel population. This is an extension of the earlier work by Gazi et al. [3]. The major, as well as, important results we had obtained for stability (local and global) of the differential equations and the delay differential equation model system in order to study the effect of gestational delay on various dynamic behavior with different level of external food input. Time delay have ability to drive a stable equilibrium point to an unstable one and it is also responsible for the oscillations of various trophic levels. We have stated different results obtained in Gazi et al. [3] in Section 2. In section three, we have carried out the main analysis: we have analyzed the stability of the periodic oscillation and the direction of Hopf bifurcation. Here b2 determines the stability of the periodic oscillation. We consider a set of some hypothetical parameter values as given in the figure caption. Then the positive equilibrium point of the system (1.3) is E⁄ = (2.2857, 1.5714, 5.0238). When s = 1.2, E⁄(2.2857, 1.5714, 5.0238) is asymptotically stable (See Fig. 1) and x0 = 2.7648, s0 = 1.3452, Re k0 (s0) = 582.7659, C1(0) = 11.0292  29.2669i. So, b2 < 0, l2 > 0. Thus E⁄ is stable when s < s0 (Fig. 1) and unstable when s > s0 (Fig. 2). Hopf bifurcation occurs and the stable solution becomes periodic oscillation when s > s0. Since b2 is negative, the periodic oscillation is stable and the Hopf bifurcation is supercritical. The direction of the bifurcation is s > s0. Acknowledgements The author is grateful to Prof. C.G. Chakrabarti, S. N. Bose Professor of Theoretical Physics, Department of Applied Mathematics, University of Calcutta, for his continuous help and guidance throughout the preparation of the paper. The research work is done under the UPE Project of Calcutta University. References [1] G.S. Jacinto, Fish cage farming in coastal waters-environmental, biological, social and governance issues, pp. 48–56, in: B.V.L. Querijero, C.R. Pagdilao and S.V. Ilagan (eds.). Guidelines in the Establishment of Fish Cages and Other Structures in Lakes and Coastal Waters. PCAMRD Book Series No. 36/ 2006. Los Baos, Laguna, Philippines. (2006) 174 p. [2] K.S.P. Shin, Shellfish used as Fish farm biofilter, Research Frontier, 10 (2005), Hong Kong, China. [3] N.H. Gazi, S.R. Khan, C.G. Chakrabarti, Mussel integration in fish farm: mathematical model and analysis, Nonlinear Anal.: Hybrid Syst. 3 (2009) 74–86. [4] J.M. Cushing, Integro-Differential Equations and Delay model in Population Dynamics, Springer-Verlag, Heidelberg, 1977. [5] K. Gopalsamy, Stability and Oscillation in Delay Differential Equations of Population Dynamics, Kluwer Academic Publisher, The Netherlands, 1992. [6] Y. Kuang, Delay differential equations with application in population dynamics, Academic Press, New York, 1993. [7] R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, New Jersey, 2001. [8] E. Beretta, Y. Kuang, Convergence results in a well-known delayed prey–predator system, J. Math. Anal.-Appl. 204 (1996). [9] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf-bifurcation, Cambridge University Press, Cambridge, 1981. [10] G. Birkhoff, G.C. Rota, Ordinary Differential Equation, Ginn. and Co., 1982. [11] X. Meng, J. Wei, Stability and bifurcation of mutual system with time delay, Chaos Solitons Fract. 21 (2004) 729–740. [12] C. Sun, Y. Lin, M. Han, Stability and Hopf bifurcation for an endemic disease model with delay, Chaos Solitons Fract. 30 (2006) 204–216. [13] Zhan-Ping Ma, Hai-Feng Huo, Chun-Ying Liu, Stability and Hopf bifurcation analysis on a predator–prey model with discrete and distributed delays, Nonlinear Anal.: Real World Appl. 10 (2009) 1160–1172.