Effect of salinity and fish predation on zooplankton dynamics in Hooghly-Matla estuarine system, India

1.1 MB·

by  Ujjwal Roy

Effect of salinity and fish predation on zooplankton
dynamics in Hooghly-Matla estuarine system, India
Ujjwal Roy(1), S. Sarwardi(2)∗, N. C. Majee(1), Santanu Ray(3)
(1) Department of Mathematics, Jadavpur University
6
1 Kolkata-700 032, West Bengal, India
0
2 email: [email protected]
n (2) Department of Mathematics, Aliah University, IIA/27, New Town
a Kolkata - 700 156, West Bengal, India
J
email: [email protected]
4
(3)Ecological Modelling Laboratory, Department of Zoology, Visva-

]
S Bharati University, Santiniketan-731 235, West Bengal, India.
D
email: [email protected]
.
h
t
a
m
Abstract
[

The Hooghly-Matla estuarine complex is the unique estuarine system of the world. Nutrient
1
v from the litterfall enrich the adjacent estuary through tidal influence which in turn regulate the
3
phytoplankton, zooplankton and fish population dynamics. Environmental factors regulate the
6
biotic components of the system, among which salinity plays a leading role in the regulation of
2
1 phytoplankton, zooplankton and fish dynamics of the estuary. In this article, a PZF model is
0 considered with Holling type-II response function. The present model considers salinity based
.
1 equations onplanktondynamics ofthe estuary. The interiorequilibriumis consideredasthe most
0
important equilibrium state of this model. The model equations are solved both analytically and
6
1 numericallyusing the realdata baseofHooghly-Matlaestuarinesystem. The essentialmathemat-
: ical features of the present model have been analyzed thorough local and global stability and the
v
i bifurcationsarisinginsomeselectedsituations. Acombinationofsetofvaluesofthesalinityofthe
X
estuaryareidentifiedthathelpedtodeterminethesustenanceoffishpopulationinthesystem. The
r
a ranges of salinity under which the system undergoes Hopf bifurcation are determined. Numerical
illustrations are performedin orderto validate the applicability of the model under consideration.
Mathematics Subject Classification: 92D25, 92D30, 92D40
Keywords: PZF model; Equilibria; Local stability; Global stability; Hopf bifurcation; Numerical
simulation
∗Authorto whom all correspondence should be addressed
1 Introduction
The delta of Hooghly-Matla estuarine system is networked by seven major river along with their
tributaries and creeks. This deltaic system harbours luxuriant mangroves and constitutes Sundarban
mangrove ecosystem. Themangroves arethemajorresources ofdetritusandnutrientstotheadjacent
estuary that make up favourable habitat for the growth of shell fish and fin fish (cf. Mandal et al.,
(2009)). This estuary act as route and refuge areas for a variety of migratory fish species. The
estuary supports 53 species of pelagic fish belonging to 27 families and 124 species under 49 families
of demersal fish (cf. Hussain and Acharya, (1994)). Fishery in this estuarine water contributes a part
in the economy of the state of West Bengal, India. The organic matter that passes from litterfall to
the adjacent estuary supports both the grazing and detritus food chains of this lotic system.
Among the chemical components studied so far, the salinity plays a crucial role in the abundance (cf.
Bhunia, (1979)) and dynamics of zooplankton of the estuary (cf. Ghosh, (2001), Ketchum, (1951)).
Because, this community lies in the middle of grazing food chain; phytoplankton-zooplankton-fish
(PZF). In addition, perturbation to this trophic level may trigger imbalance in the food chain which
in turn affect the phytoplankton and fish community. The dynamics of salinity is season dependent.
During monsoon and early post monsoon, huge fresh water enters in the estuary from the upstream
resultingintheloweringofsalinity. Inthepre-monsoon,freshwaterrunofffromtheupstreambecomes
very less and due to tidal influence of the adjacent Bay of Bengal, the salinity increases. Throughout
the year, a gradient of salinity is observed between upstream and downstream area of the estuary (cf.
Mandal et al., (2009)).
The abundance of different species of zooplankton varies according to the salinity of the estuary
throughout the year, as a result the grazing rate also changes with seasons. In PZF system, the
grazing rate of zooplankton is one of the most sensitive parameter in Hooghly-Matla estuarine system
as the dynamics of zooplankton also depends on the lower trophic level of phytoplankton as well
as fish population (cf. Mandal et al., (2009)). Moreover, fish predation exhibits top-down effect on
zooplankton community. Therefore, in the PZF system, two important parameters that shape up
the zooplankton dynamics of the estuary are salinity dependent grazing rate of zooplankton and fish
predation rate (cf. Dube et al., (2010)).
Therehave been only a few models based on the effects of salinity on plankton dynamics. Few studies
have been done in the relationship between salinity levels to the types of species that can occur in
algal blooms (cf. Griffin et al., (2001), Marcarelli et al., (2006), Quinlan and Phlips, (2007)). Many
PZF model are constructed on estuarine system (cf. Ray et al., (2001); Cottingham, (2004); Dube
and Jayaraman, (2008); Dube et al., (2010)) and top-down effect (cf. Morozov et al., (2005); Irigoien
et al., (2005); Calbet and Saiz, (2007)) but noneof the models have considered the salinity dependent
grazing rate of zooplankton which plays an important rule in the zooplankton dynamics of estuary.
The present account deals with a PZF model, where salinity dependent grazing rate of zooplankton
is taken into consideration along with the variation of upstream and downstream salinity.
2
1.1 Mathematical model formulation
Let P, Z, F denote the populations of phytoplankton, zooplankton and fish respectively. In this
present PZF model, light and temperature dependant photosynthesis rate of phytoplankton and
salinity induced grazing rate of zooplankton have been incorporated into the model proposed by
Mandal et al., (2011). The modified model under consideration is as follows:
dP = m P 1 P g PZ
dt 1 − kP − sP+kZ
 ddZt = agsPP+(cid:16)ZkZ −gf(cid:17)ZZ+FkF −m2Z (1)
 dF = g ZF m F.
dt fZ+kF − 3


The salinity induced grazing of zooplankton is g = δg , where δ = su , the dilution factor and
s Z su−sd
m , s , s , g , k , k are net growth rate of phytoplankton, downstream salinity, upstream salinity,
1 d u Z Z P
grazing rate of zooplankton, half saturation constant of zooplankton on phytoplankton grazing by
zooplankton and half saturation constant grazing by phytoplankton respectively. Since estuary is a
transition zone of river and sea, so there is always fluctuation of salinity throughoutthe year, which is
dueto dilution by upstream river water and/or mixing by downstream sea tidal water. δ is calculated
by following the equation of Ketchum, (1951). Besides grazing the abundance of zooplankton is also
dependent on loss due to E = ZE , respiration R = Zr , fish predation F = Zr , and
zoo zo zoo zo p rf
mortality M = ZM , R is governed by r .
zoo Z zoo zo
The net mortality rate of zooplankton is m = (E +r +r +M ).TheabundanceofF isgoverned
2 zo zo fp Z
by many processes in the estuary. Fish predation F on Z follows Michaelis-Menten kinetics (Holling
p
type II) which enriches the fish pool of the estuarine system. F depends on k and g . F population
p Z f
is reduced by mortality rate of fish M , respiration Fr , harvest by fishing FH and excretion of
f f cf
fish FE . FH are controlled by R and H respectively. Respiration rate of Fish and H Harvest
f f f f f
rate of Carnivorous Fish respectively.
The net mortality rate of fish is m = (E +M +r +H ),wherer ,r and M arefishpredation
3 f f f f fp zo rz
rate, zooplankton respiratory rate and mortality rate respectively.
1.2 Existence and positive invariance
Letting X (P,Z,F)T, f : R3 R3,f =(f ,f ,f )T, the system (1) can be rewritten as X˙ = f(X).
1 2 3
≡ →
Here f C∞(R) for i = 1,2,3, where f = m P(1 P ) g PZ , f = ag PZ m Z ZF ,
i ∈ 1 1 − kP − sP+kZ 2 sP+kZ − 2 − Z+kF
f = g ZF m F. Since the vector function f is a smooth function of the variables P, Z and F in
3 fZ+kF − 3
the positive octant Ω0 = (P,Z,F);P > 0,Z > 0,F > 0 , the local existence and uniqueness of the
{ }
system hold.
1.3 Boundedness of the system
Boundedness of a system guarantees its biological validity. The following theorem establishes the
uniform boundedness of the system (1).
Theorem 1. All the solutions of the system (1) which start in R3 are uniformly bounded.
+
3
Proof. Let (P(t),Z(t),F(t)) be any solution of the system with positive initial conditions. From the
real ecological field study one can consider max m , m < 1.
2 3
{ }
Now let us we define the function X = aP +Z +F. The time derivative of X gives
dX dP dZ dF
= a + + ,
dt dt dx dt
P PZ PZ ZF ZF
= am P(1 ) ag +ag g m Z +g m F
1 s s f 2 f 3
− k − P +k P +k − Z +k − Z +k −
P Z Z F F
P
= am P(1 ) m Z m F. (2)
1 2 3
− k − −
P
Now for each v > 0, we have
dX P
+vX = am P(1 ) m Z m F +v(aP +Z +F)
1 2 3
dt − k − −
P
P
= am P(1 )+avP +(v m )Z +(v m )F
1 2 3
− k − −
P
P
am (1 )+avP, if v min m , m ,
1 2 3
≤ − k ≤ { }
P
a(m +v)2k
1 P
. (3)
≤ 4
Using the theory of differential inequality (cf. Birkoff and Rota, (1982); Sarwardi et al. (2011)), one
can easily obtain
a(m +v)2k
1 P
limsupX(t) = ρ. (4)
t→+∞ ≤ 4
Therefore, all the solutions of the system (1) enter into the compact region = (P,Z,F) R3 :
B ∈ +
aP +Z +F a(m1+v)2kP , which completes the proof. n
≤ 4
o
1.4 Equilibria and their feasibility
To determine the equilibrium points of the system of(1) we put P˙ = Z˙ = F˙ = 0. The equilibria of
the system(1) are (i) the null equilibrium E = (0,0,0);
0
(ii) the axial equilibrium E = (k ,0,0);
1 P
(iii) the boundary equilibrium E = ( m2kZ , (am1kZ)(akPgs−kPm2−m2kZ),0); and
2 (ags−m2) (ags−m2)2kP
(iv) the interior equilibria Ei = (P∗,Z∗,F∗); i = 1,2,
∗ i i
where P∗=−m1(gf−m3)(kP−kZ)±√(m1(gf−m3)(kP−kZ))2−4m1(m3−gf)(kPkZm1(gf−m3)−m3kfgskP),
i 2m1(m3−gf)
Z∗= m3kF and F∗ = m1kFa(P∗−kP) + kF(ags−m2).
i (gf−m3) i m3kP (gf−m3)
1.5 Existence of planer and interior equilibria and their stability
Equilibria of the model (1) can be obtained by solving P˙ = Z˙ = F˙ = 0. Though the system (1) has
several non-negative steady state but we have mentioned the existence of interior equilibrium point
4
only for its biological importance.
∗ ∗ ∗
The interior equilibrium point E∗ = (P ,Z ,F ) of the system (1) exists if (gf m2) > 0 . When
−
∗
these conditions are satisfied, we have P as the positive root of the following quadratic equation:
A P∗2 A P∗ A = 0, (5)
0 1 2
− −
where
A = m (g m ), A = (g m m m )(k k ), A =k k m (g m ) m k g k.
0 1 f 3 1 f 1 3 1 P Z 2 P Z 1 f 3 3 f s
− − − − − −
Since g > m , A is positive. Therefore, one positive root of (5) can be found as
f 3 0
1
P∗ = [A + A2+4A A ], (6)
2A 1 1 0 2
0
q
which exists if A > 0 and A > 0. Therefore, the equation (5) has a unique positive solution if
2 0
∗
g > m . Once we get the unique positive solution of P from equation (5), it is easy to obtain the
f 3
other components of the interior equilibrium E∗.
Feasibility: It is clear that the axial equilibrium is feasible. The interior equilibrium is feasible if
the condition (i) (g m ) > 0 and (ii) (k k ) < 0 hold good.
f 3 P z
− −
1.6 Local stability of equilibrium
In order to find out about the stability of the equilibrium points we need to linearize the system (1).
Sincethe axial and boundaryequilibria are less important on this system (1). That’s why we don’tgo
for detailed analysis of the system around the equilibria other than interior equilibrium. The stability
of interior equilibrium, we find the Jacobian matrix. For this purpose we define the system (1) as
follows:
dX
= F(X), (7)
dt
m P(1 P ) g PZ P
1 − kP − sP+kZ
where f(X)= ag PZ g ZF m Z , and Y = Z .
 sP+kZ − fZ+kF − 2   
g ZF m F F
 fZ+kF − 3   
   
The variational matrix of the system at any arbitrary point ( P˜,Z˜,F˜ ) is
2g Z˜ g Z˜k g P˜
s s Z s
m 0
P˜ +k − 1− (P˜ +k )2 −(P˜ +k )
 Z Z Z 
ag k Z˜ g F˜ g F˜k g Z˜
J˜= ∂f = s Z f f F f .
∂X |(P˜,Z˜,F˜)  (P˜ +kZ)2 (Z˜+kF) − (Z˜+kF)2 −(Z˜+kF) 
 gfF˜kF 
 0 0 
 (Z˜+k )2 
 F 
 
Therefore, the Jacobian matrix for the system (1) at E∗ is given by
R R m K 0
1 2 1 1
− − −
J∗ =  R2 K2 R3 m3 ,
− −
0 R 0
3
 
 
5
where R = 2gsZ∗ , R = gsZ∗kZ , R = gfF∗kF , K = gsP∗ , K = gfF∗ .
1 (P∗+kZ) 2 (P∗+kZ)2 3 (Z∗+kF)2 1 (P∗+kZ) 2 (Z∗+kF)
The characteristic equation corresponding to the Jacobian matrix J∗ is
∆(λ) = λ3+D λ2+D λ+D = 0. (8)
1 2 3
For thestability of thesolutions of the system (1)all of theroots of thecharacteristic equation for the
Jacobian matrix at E∗ ie. J∗ should have negative real parts. This can be achieved without actually
solving for all the roots of the characteristic equation, but by applying the Routh-Hurwitz condition
for negative real parts of the characteristic roots. According to the Routh-Hurwitz conditions, the
solutions of (1) will be asymptotically stable if the conditions D > 0, D D > D , D > 0 are
1 1 2 3 3
satisfied.
2 Global Stability
Now we examine the global stability issue. Here we describe the general method described by Li and
Muldowney, (1996) to show an n-dimensional autonomous dynamical system f : D Rn,D Rn,
→ ⊂
an open and simply connected set and f C1(D), where the dynamical system is given by
∈
dx
= F(x) (9)
dt
is globally stable undercertain parametric conditions (cf. Jin and Haque, (2005); Haque et al., (2009)
and Buonomo et al., (2008))
Now we assume the following incidents:
∗
(A1) The autonomous dynamical system has a unique interior equilibrium point x in D.
(A2) The domain D is simply connected.
(A3) There is a compact absorbing set Ω D.
⊂
Now we give a definition due to Haque et al., (2005)
∗
Definition 1. Theuniqueequilibriumpointx ofthedynamicalsystem(9)isglobally 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=(J ) be the variational matrix of the system (9) and J[2] be second additive compound
ij 3x3
matrix with order nC nC . In particular for n=3 we can write
2 2
×
V +V V V
11 22 23 13
J[2] = ∂f[2] = V V +V −V Let M(x) C1(D), be an nC nC matrix
∂X  32 11 33 12  ∈ 2 × 2
V V V +V
31 21 22 33
 − 
valued function. Moreover, we consider B as a matrix such that B=M M−1+MJ[2]M−1 where the
f
matrix M is represented by
f
∂M t
ij
(M (x)) = .f(x)= M .f(x). (10)
ij f ij
∂x ∇
(cid:18) (cid:19)
6
Again due to Jr. Martin et al., (1974), we consider the Lozinskiˇi measure Γ of B with respect to a
vector norm . in Rn . N = nC , we have
2
||
I +hB 1
Γ(B)= lim | |− .
h→0+ h
Li and Muldowney (1996) showed that the system (9) will be globally stable if the conditions (A1),
(A2) and (A3) together with
1 t
limsupsup Γ(B(x(s,x )))ds < 0 (11)
0
t→∞ t Z0
hold simultaneously.
The above stated condition not only assures that there are no orbits (i.e., homoclinic orbits, hetero-
clinic cycles and periodic orbits) which gives rise to a simple closed rectifiable curve in D, invariant
for the system (9) butalso it is a robustBendixson criterion for (10) in addition with (11) is sufficient
for local stability thesystem (9)around E∗. Now weuse theabove discussionto show that our system
(1) is globally stable around its interior equilibrium. We make the following transformation of the
variables x P−1, y Z, z F, for which the autonomous system (1) is transformed to the
→ → →
following one:
dX
= f(X), (12)
dt
m1(1 1 ) gsy x2 x
x − xkP − 1+xkZ −
where f(X) = agsy gfyz m y and X = y . The variational matrix of the
 1+xkZ − y+kF − 2   
gfyz m z z
 y+kF − 3   
system (12) canbe written as   
xg y(2+xk ) g x2
s Z s
m 0
(1+xk )2 − 1 (1+xk )
 Z Z 
ag yk ag ag y g yz g y
V = ∂f = s Z s s + f f .
∂X  −(1+xk )2 (1+xk ) − (1+xk ) (y+k )2 −(y+k ) 
 Z Z Z F F 
 gfZkF 
 0 0 
 (y+k )2 
F
 
 
If V[2] be the second additive compound matrix of V then due to Buonomo et al. (2008) we can write
g y
f
a 0
11
−(y+k )
F
 g Zk g x2 
f F s
V[2] = (y+k )2 a22 (1+xk ) , where
 F Z 
 ag yk 
s Z
 0 a 
 −(1+xk )2 33 
 Z 
 
xg y(2+xk ) g yz
s Z f
a = m + ,
11 (1+xk )2 − 1 (y+k )2
Z F
xg y(2+xk )
s Z
a = m ,
22 (1+xk )2 − 1
Z
g yz
f
a = .
33 (y+k )2
F
We consider P(X) in C1(D) in a way that P = diag 1, y,yx2 then we have P−1 = diag 1, z, z ,
z z y yx2
(cid:0) (cid:1) (cid:0) (cid:1)
7
P P−1 = diag 0, y˙ z˙, 2x˙ + y˙ z˙ and
f { y − z x y − z}
a gfy 0
11 −(y+kF)
PV[2]P−1 = gfkFz a gsx2 . After some algebraic calculation, we get
 (y+kF)2 22 (1+xkZ) 
0 agsykZ a
 −(1+xkZ)2 33 
 
B B
B = P P−1+PV[2]P−1 = 11 12 ,
f
B B
" 21 22 #
where
B = a = xgsy(2+xkZ) m + gfyz ,
11 11 (1+xkZ)2 − 1 (y+kF)2
T a + y˙ z˙ gsx2
B = gfy 0 , B = gfzkF 0 , B = 22 y − z (1+xkZ) .
12 h −y+kF i 21 h (y+kF)2 i 22 " −(1a+gsxykkZZ)2 a33+ 2xx˙ + yy˙ − zz˙ #
Now letusdefinethefollowing vector norminR3 as (u,v,w) =max u, v + w ), where (u,v,w) is
| | | | | | | | | |
the vector norm in R3 and it is denoted by Γ, the Lokˇinski measure with respect to this norm. There-
(cid:0)
fore,Γ(B) sup l ,l ,wherel = Γ (B )+ B fori=1,2andi=j,here B and B arethematrix
1 2 i 1 ii ij 12 21
≤ { } | | 6 | | | |
norms with respect to the L1 vector norm and Γ is the Lokˇinski measure with respect to that norm.
1
Therefore, we can easily obtain the following terms: Γ (B ) = xgsy(2+xkZ) m + gfyz , B =
1 11 (1+xkZ)2 − 1 (y+kF)2 | 12|
max gfy , 0 ,
{y+kF }
B = max gfzkF , 0 ,
| 21| {(y+kF)2 }
Γ (B ) = y˙ z˙+max gsy(2x+x2kZ−akZ), gsx2 + gfyz + 2x˙) .Usingthesystemequation (12),
1 22 y−z (1+xkZ)2 (1+xkZ) (y+kF)2 x
we have l = xgsy(2+xkZn) m + gfyz + gfy , l = y˙ z˙ 2mo(1 1 )+ gsyx + gfykF .
1 (1+xkZ)2 − 1 (y+kF)2 y+kF 2 y − z − 1 − xkP 1+xkZ (y+kF)2
Now from the expression of l and l one can obtain
1 2
Γ(B)= max xgsy(2+xkZ) m + gfyz + gfy ,y˙ z˙ 2m (1 1 )+ gsyx + gfykF .
(1+xkZ)2 − 1 (y+kF)2 y+kF y − z − 1 − xkP 1+xkZ (y+kF)2
(cid:26) (cid:27)
= y˙ z˙ 2m (1 1 )+ gsyx + gfykF , if 1< xk < 2.
y − z − 1 − xkP 1+xkZ (y+kF)2 P
Thus,
Γ(B) y˙ 2gs (gs + gf )ρ (m +m )
≤ y− ρ+kF − kZ kF − 1 3
(cid:20) (cid:21)
y˙ µ, where µ = 2gs (gs + gf )ρ (m +m )> 0.
≤ y − ρ+kF − kZ kF − 1 3
Taking average value of Γ(B) within the time scale 0 and t one can have
1 t 1 y(t)
Γ(B)ds log µ,which gives
t ≤ t y(0) −
Z0
1 t
lim sup sup Γ(B(x(s,x )))ds < µ < 0, (13)
0
t→∞ t −
Z0
whenever 1 > 1 > 1.
xkP 2
Now we are in a position to state the following theorem:
8
Theorem 2.1 The model system (1) is globally asymptotically stable around its interior equilibrium
if µ > 0, i.e 2gs (gs + gf )ρ (m +m ) > 0.
ρ+kF − kZ kF − 1 3
2.1 Numerical simulation
Forthepurposeofmakingqualitativeanalysis,numericalsimulationshavebeencarriedoutbymaking
use of MATHLAB-R2010a and Maple-18. In Table 1, we summarize the used parameters with their
admissible values and their biological interpretations. All these results have also been verified by
numerical simulations,of which we reportsome in the figures. In order to compare our modelwith the
corresponding model. we have run simulations using the standard MATLAB differential equations
integrator for the Runge-Kutta method, i.e, MATLAB routine ODE45. The numerical experiments
performed on the system (1) using the experimental data taken from different sources to confirm
our theoretical findings. Specially we focus here on the following values of the system parameters
as shown in the table-I. The problem described by the system (1) is well posed. The P, Z, F axes
are invariant under the flow of the governing system. In the present system the total environmental
population is bounded above (cf. Subsection 1.3). Therefore, any solution starting in the interior
of the first octant never leaves it. This mathematical fact is consistent with the biologically well
behaved system and is common in the modern research work on mathematical biology. We took a
set of parameter values: m = 0.6, m = 0.0698, m = 0.324, a = 0.8, g = 0.75, g = 0.6894, s =
1 2 3 s f d
12.30, s = 8.23, k = 12, k = 38, k = 10.1. For this set of parameter values the system possessed
u P Z F
a unique equilibrium point E∗(1.809,8.964,3.112). We have established the sufficient conditions for
global stability of the coexistence equilibrium. The coexistence coexistence equilibrium point E∗ has
been found through numerical simulations whose global asymptotical stability has been depicted in
Figures 1(a) and Figure 2(a). It is observed from the Figure 1:(b)-(c) and Figure 2(b)-(c) that the
fish population gradually decreases as the upstream salinity increases from s = 8.33 to s = 8.51
u u
while the downstream salinity is keeping fixed at s = 12.30. For higher upstream salinity (s = 8.51)
d u
Figure 1(c) in the estuary is detrimental to the present zooplankton population. Zooplankton in
this particular region are stenohaline thus the fish population cannot persist in such high salinity
and the population gradually decreased to zero. Therefore, the salinity has an important role on
regulating fish population in the present estuarine system. Representative numerical simulations are
shown in Figures: 3–5. In Figure 3(a), the solutions plots describing the Hopf bifurcation and Figure
3(b) exhibits the corresponding phase space diagram. The present dynamical system confirms the
existence of chaotic nature which has been presented in Figure 4:(a)-(b). Figure 5:(a)-(f) shows
different phase space diagram for different values of upstream salinity within the range 4.50–7.50,
while downstream salinity has been kept fixed at s = 10.30 and the remaining parameter values are
d
same as taken in Figure 1. It is observed the catastrophic change in the behaviour of the solution
plots as the parameter s varies from 4.50 to 8.0 (cf. Figure 5:(a)-(f)). More specifically, Figure 5(a)
u
shows periodic orbit with period 4; Figure 5(b) shows the chaotic orbit; Figure 5(c) shows periodic
orbit with period 6; Figure 5(d) shows the period doubling bifurcation; Figure 5(e) demonstrates a
limit cycle and finally Figure 5(f) ensures the stability of the system around an interior equilibrium
point (2.622,8.878,10.12) of the system (1).
9
Table 1: The set of parameter values and their biological terminologies
Descriptions Fixed Values Units References
Parame-
ters
Growth rate of phyto- m 0.6 day−1 Mandal et al., (2012)
1
plankton
Half saturation constant k 12 dimensionless Mandal et al., (2012)
P
for nutrient uptake by
phytoplankton
Grazing rate of zoo- g 0.75 day−1 Ray and Straskraba,
z
plankton (2001)
Optimal light intensity E 0.04 day−1 Ray et al., (2001)
zo
Respiration rate of zoo- r 0.0153 day−1 Ray et al., (2001)
zo
plankton
Fish predation r 0.2 day−1 Ray et al., (2001)
fp
Mortality rate of zoo- M 0.0145 day−1 Ray et al., (2001)
z
plankton
Excretion of fish E 0.049 day−1 Mandal et al., (2012)
f
Mortality rate of fish M 0.021 day−1 Mandal et al., (2012)
f
Respiration rate of Fish r 0.0125 day−1 Mandal et al., (2012)
f
Harvest rate of Carnivo- H 0.1090 day−1 Mandal et al., (2012)
f
rous Fish
Growth rate of Carnivo- g 0.6894 day−1 Mandal et al., (2012)
f
rous Fish
Half saturation constant k 10 dimensionless Mandal et al., (2012)
F
for zooplankton grazing
by Carnivorous Fish
Half saturation constant k 38 dimensionless Mandal et al., (2012)
Z
for phytoplankton graz-
ing by zooplankton
10