![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
Abstract
This paper is devoted to the study of a plankton–fish ecosystem model. The model represents the interaction between phytoplankton, zooplankton, and fish with Holling II functional response consisting of carrying capacity and constant intrinsic growth rate of phytoplankton. It is observed that if the carrying capacity of phytoplankton population crosses a certain critical value, the system enters into Hopf bifurcation. We have introduced discrete time delay due to gestation in the functional response term involved with the growth equation of planktivorous fish. We have studied the effect of time delay on the stability behavior. In addition, we have obtained an estimate for the length of time delay to preserve the stability of the model system. Existence of Hopf bifurcating small amplitude periodic solutions is derived by considering time delay as a bifurcation parameter. It is observed that constant intrinsic growth rate of phytoplankton and mortality rate of planktivorous fish play an important role in changing one steady state to another steady state and oscillatory behavior of the system. Computer simulations illustrate the results.
Keywords:
Public Interest Statement
Economically important fish species have long been regarded in isolation from each other and their habitat. In order to comprehensively assess the impacts of fisheries, the entire habitat must be considered. Only then will a sustainable and economic fishery system be possible. Predation by fish determines the abundance of herbivorous zooplankton which in turn regulates the level of phytoplankton. Also, changes in the abundance of planktivorous fish do affect both the phytoplankton and zooplankton. The relationship between phytoplankton, zooplankton, and fish culture is of paramount importance in determining the water quality, on one hand, and the natural productivity and the fish production, on the other. It is observed that constant intrinsic growth rate of phytoplankton and mortality rate of planktivorous fish play an important role in changing one steady state to another steady state and oscillatory behavior of the system.
1. Introduction
It is well known that plankton plays an integral role in marine ecosystem. Phytoplankton and zooplankton are the main two types of plankton. Many models have already been built to simulate zooplankton–phytoplankton interactions. In marine ecosystem, much attention has been paid to the effect of fish and plankton biomass. Top-down and bottom-up effects on plankton–fish dynamics have been discussed by Scheffer (Citation1991) and Pal and Chatterjee (Citation2012). The authors have discussed that at high fish densities, zooplankton is controlled by fish predation and algal biomass is light or nutrient limited, whereas at low fish densities, zooplankton is food limited and phytoplankton density is controlled by zooplankton grazing. Scheffer, Rinaldi, and Kuznetsov (Citation2000) have discussed about isocline analysis and simulations to show how an increase in fish predation may affect plankton dynamics in the model and to explain which type of bifurcations may arise. Role of gestation delay in a plankton–fish model under stochastic fluctuations has been discussed by Mukhopadhyay and Bhattacharyya (Citation2008). Realistic patterns of patchiness in plankton–fish dynamics have been studied by Upadhyay, Kumari, and Rai (Citation2009). The authors also studied both one-dimensional and two-dimensional reaction–diffusion model of nutrient–phytoplankton–zooplankton–fish interaction. As studied by Biktashev, Brindley, and Horwood (Citation2003), the food source for fish larvae depends on the zooplankton dynamics, which in turn is coupled to larval populations through the zooplankton mortality term. A conceptual mathematical model of fish and zooplankton populations inhabiting the lakes Naroch and Myastro has been developed and studied by Medvinsky et al. (Citation2009).
In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay can cause destabilization of equilibria and can induce various oscillations and periodic solutions. As observed by Bischi (Citation1992), the effects of the time delay involved nutrient recycling on resilience, that is, the rate at which a system returns to a stable steady state following a perturbation. The researchers have suggested that characterizing a system by oscillatory behavior, an increase in the distributed time delay, can have a stabilizing effect. But it has been found that the introduction of time delays is a destabilizing process, in the sense that increasing the time delay could cause a stable equilibrium to become unstable and/or cause the populations to fluctuate (Cushing, Citation1977; Freedman & Rao, Citation1983). Chemostat-type models incorporating discrete delays have been investigated by Freedman, So, and Waltman (Citation1989). Also, delayed models in population biology have been discussed in Gopalsamy (Citation1984), Kuang (Citation1993), Khare, Misra, Singh, and Dhar (Citation2011) and Liu, Beretta, and Breda (Citation2010). A discrete time delay term to the model of Beretta, Bischi, and Solimano (Citation1991) was introduced in Ruan (Citation1995). This discrete time delay term may be considered as delay due to gestation. The author also allowed the washout rates for nutrient and plankton to be different. Recently, effect of delay on nutrient cycling in phytoplankton–zooplankton interaction model has been studied by Das and Ray (Citation2008). Rehim and Imran (Citation2012) have induced a discrete time delay to both the consume response function and distribution of toxic substance terms in phytoplankton–zooplankton model. They have found out a range of gestation delays which initially impart stability, then induce instability and ultimately lead to periodic behavior. Due to gestation of prey, a discrete time delay toxin producing phytoplankton–zooplankton system shows rich dynamic behavior including chaos and limit cycles (Singh & Gakkhar, Citation2012). In particular, Gazi and Bandyopadhyay (Citation2006) have introduced discrete time delay due to recycling of dead organic matters and gestation of nutrients to the growth equations of various trophic levels. They have studied the effect of time delay on the stability behavior and obtained an estimate for the length of time delay to preserve the stability of the model system. Effects of time delay on a harvested predator–prey model have discussed in Gazi and Bandyopadhyay (Citation2008), Maity, Patra, and Samanta (Citation2007) and Yongzhen, Min, and Changguo (Citation2011).
Here, an open system is considered with three interacting components consisting of phytoplankton (P), zooplankton (Z), and planktivorous fish (F). In this paper, a plankton–fish interaction model is described. The stability of equilibrium point is described. We have derived the conditions for instability of the system around the interior equilibrium and Hopf bifurcation in presence of delay and non-delay models. Direction of Hopf bifurcation is discussed. Numerical simulations under a set of parameter values have been performed to support our analytical results.
2. The mathematical model
Let be the concentration of the phytoplankton at time t with carrying capacity K and constant intrinsic growth rate r. Let
be the concentration zooplankton biomass and
be the concentration of planktivorous fish biomass present at any instant of time t.
Let be the maximal zooplankton ingestion rate and
maximal zooplankton conversion rate for the growth of zooplankton, respectively (
). Let
be the maximal planktivorous fish ingestion rate and
be the maximal planktivorous fish conversion rate due to grazing of herbivorous zooplankton. Further,
,
,
denote the mortality rates of the phytoplankton, zooplankton, and planktivorous fish biomass, respectively. We choose Holling type II functional form to describe the grazing phenomena with
and
as half saturation constant.
The basic equations with all of the parameters are:(1)
(1)
Set and
, with
. The system (Equation 1) can be written in compact form as
. Its Jacobian is
(2)
(2)
3. Some preliminary results
3.1. Positive invariance
The system (Equation 1) is not homogeneous. However, it is easy to check that whenever choosing with
, for
, then
. This ensures that the solution remains within the positive orthant, ensuring the biological well posedness of the system.
3.2. Boundedness of the system
Theorem 1
All the solutions of (Equation 1) are ultimately bounded.
Proof
We define a function
We have which gives
as
, where
is constant.
For the large value of t, we get , where
.
Let be the solution of
, satisfying
.
Then, as
. By comparison, it follows that lim
, proving the theorem.
3.3. Equilibria
The system (Equation 1) possesses the following four equilibria: the plankton free equilibrium , the zooplankton free equilibrium
, the planktivorous fish free equilibrium
, and possibly the coexistence of the three populations
.
3.3.1. Plankton free equilibrium
is always feasible; the Jacobian (Equation 2) evaluated at this equilibrium has the eigenvalues
,
and
. Clearly,
is always unstable.
3.3.2. Zooplankton free equilibrium
At , the Jacobian (Equation 2) factorizes to give three explicit eigenvalues
,
and
, where
is defined below. Thus,
is locally asymptotically stable if and only if
(3)
(3)
3.3.3. Planktivorous fish free equilibrium
At , the population levels are
This equilibrium is feasible if it is satisfying(4)
(4)
At , the Jacobian (Equation 2) factorizes to give one explicit eigenvalue
and the quadratic
. The Routh–Hurwitz conditions for the latter are easily seen to hold, but only when (Equation 4) is satisfied. Thus, stability of
is then ensured by
(5)
(5)
3.3.4. The coexistence equilibrium
The coexistence equilibrium cannot be found explicitly since
,
and
, where
,
,
.
Thus, the condition for the existence of the interior equilibrium point is given by,
. For feasibility, we certainly need
Stability analysis of the positive interior equilibrium of the system (Equation 1)
The variational matrix of system (Equation 1) around the positive equilibrium is
where ,
,
,
,
.
The characteristic equation is ,
where ,
,
.
Case1: If
, it implies that
. Thus,
is unstable.
Case2: If
, it implies that
.
Also, if
where
,
and
is negative. So,
if
. Therefore, according the Routh–Hurwitz criteria, all roots of the above cubic equation have negative real parts satisfying
, and
. Then, the system becomes locally asymptotically stable around
. Thus, depending upon system parameters, the system may exhibit stable or unstable behavior in this case.
Table 1. The table representing thresholds and stability of steady states
The analytical results are summarized in Table .
3.4. Hopf bifurcation at coexistence
Theorem 2
When the carrying capacity K crosses a critical value, say , the system (Equation 1) enters into a Hopf bifurcation around the positive equilibrium, which induces oscillations of the populations.
Proof
The necessary and sufficient conditions for the existence of the Hopf bifurcation for , if it exists, are (i)
,
(ii)
and (iii)
. The transversality condition, the third (iii), is obtained observing that the eigenvalues of the characteristic equation, of the form
, must satisfy the transversality condition
.
We will now verify the Hopf bifurcation condition (iii); putting in the characteristic equation, we get
(6)
(6)
On separating the real and imaginary parts and eliminating , we get
(7)
(7)
It is clear from the above that if
. Further, at
,
is the only root since the discriminate
if
.
Further, differentiating (Equation 7) with respect to K, we have , where
.
At , we have
, so that above equation becomes
, providing the third condition (iii).
This ensures that the above system has a Hopf bifurcation around the positive interior equilibrium .
Next, we perform detailed analysis of the bifurcation solutions to study the nature of Hopf bifurcation. We start by rewriting system (Equation 1) in the form where
:
,
,
and
.
Here, ,
where ,
,
.
Now, let the normalized eigenvectors of M and corresponding to the eigenvalues
and
be q and p, respectively, where
,
. Let us define the following multi-linear functions for the three vectors
,
, and
.
where ,
With these, will be of the form
.
Solving the corresponding linear system, we get
where
Again,
where .
Therefore,
where ;
;
.
Using the invariant expression for the first lyapunov coefficient , we get the first lyapunov coefficients as Mukhopadhyay and Bhattacharyya (Citation2011).
.
Now, if , the Hopf bifurcation will be of supercritical nature and there will be emergence of stable limit cycles.
4. The mathematical model with time delay
We have already discussed the usefulness of time delay in realistic modeling of ecosystems. Here, we consider the following modification of the model (Equation 1) incorporating discrete time delay in it. In this paper, we have introduced a discrete time delay term to the above model; this term may be considered as delay due to gestation. Here, is the discrete time delay.
(8)
(8)
The system (Equation 8) has to be analyzed with the following initial conditions,(9)
(9)
4.1. Qualitative analysis of the model with time delay
Remark 1
The characteristic equation for the variational matrix about the plankton free steady states
and
remains the same as obtained for the non-delayed system (Equation 8). Thus, in our model system, the delay has no effect on the stability nature of the system about
and
.
To find conditions for the other local asymptotic stability of system (Equation 8), we use the following theorem from Gopalsamy (Citation1992).
Theorem 3
If denotes the characteristic equation, then a set of necessary and sufficient conditions for the equilibrium to be asymptotically stable for all
are
(i) | The real parts of all the roots of | ||||
(ii) | For all real |
The characteristic equation for the variational matrix about
takes the following form
(10)
(10)
where ,
,
,
,
,
.
Here, has roots with negative real parts, provided system (Equation 8) is locally asymptotically stable about equilibrium
(for the condition, see Section 3.3.3).
For ,
.
Now, for ,
.
Let =0 and separating the real and imaginary parts, we get
,
.
Squaring and adding the above two equations, we get
, where
,
,
.
Sufficient conditions for the non-existence of a real number satisfying
can be written as
,
which can be transformed to .
Therefore, a sufficient condition for to be stable is (i)
and (ii)
.
4.2. Estimation of the length of delay to preserve stability at ![](//:0)
![](//:0)
In this section, we assume that in the absence of delays, is locally asymptotically stable. Let
be the respective linearized variables of this model. The system (Equation 8) becomes
(11)
(11)
where ,
,
,
,
,
.
Let ,
, and
be the Laplace transform of
,
, and
, respectively. Taking the Laplace transformation of system (Equation 11), we have
where .
The inverse Laplace transform of ,
, and
will have terms which exponentially increase with time, if
,
, and
have pole with positive real parts. Since
needs to be locally asymptotically stable, it is necessary and sufficient that all poles of
,
, and
have negative real parts. We shall employ the Nyquist criterion which states that if
is the arc length of a curve encircling the right half-plane, the curve of
,
, and
will encircle the origin a number of times equal to the difference between the number of poles and the number of zeroes of
,
, and
in the right half-plane.
Let from (Equation 10). Then,
is locally asymptotically stable if
(12)
(12)
(13)
(13)
where is the smallest positive root of Equation 12 (Freedman, Erbe, & Rao, Citation1986).
Now, (Equation 12) and (Equation 13) become.
To get an estimation on the length of delay, we utilize the following conditions,(14)
(14)
(15)
(15)
Therefore, will be stable if the above inequality holds at
, where
is the first positive root of Equation 14. We shall now estimate an upper bound
of
, which would be independent of
. Then, we estimate
so that (Equation 15) holds for all
, and hence in particular at
.
Maximizing , subject to
,
, we obtain
(16)
(16)
Thus, the unique positive solution of is denoted by
. Hence, if
, then from (Equation 16) we have
which ensures that
is independent of
. Then, we estimate
so that (Equation 15) holds for all values
; now, rearranging (Equation 15) by
and
, we have
(17)
(17)
Thus, (Equation 15) will be satisfied if , where
,
,
.
Then, the Nyquist criterion holds for , where
and
gives an estimate for the length of delay for which stability is preserved.
Theorem 4
If there exist in
such that
, then
is the maximum value (length of delay) of
for which
is asymptotically stable.
4.3. Qualitative analysis of the model at ![](//:0)
with gestation time delay (![](//:0)
)
The characteristic equation for the variational matrix about
takes the following form:
(18)
(18)
Now, .
Let and separating the real and imaginary parts, we get
Squaring and adding the above two equations, we get
where
The sufficient conditions for the non-existence of a real number satisfying
can be written as
, which can be transformed to
Therefore, a sufficient condition for to be stable is (a)
and (b)
.
4.4. Estimation of the length of delay to preserve stability at ![](//:0)
![](//:0)
In this section, we assume that in the absence of delays, is locally asymptotically stable. Let
be the respective linearized variables of this model. The system (Equation 8) becomes
(19)
(19)
where ,
,
,
,
,
,
,
.
Let ,
, and
be the Laplace transform of
,
, and
, respectively. Taking the Laplace transformation of system (Equation 19), we have
where ,
.
Let .
Then, is locally asymptotically stable if
and
i.e.
where be the smallest positive root of
.
To get an estimation on the length of delay, we use the following conditions,(20)
(20)
(21)
(21)
Therefore, will be stable if above inequality holds at
, where
is the first positive root of Equation 20. We shall now estimate an upper bound
of
, which would be independent of
. Then, we estimate
so that (Equation 21) holds for all
, and hence in particular at
.
Maximizing , subject to
,
, we obtain
(22)
(22)
Thus, the unique positive solution of is denoted by
. Hence, if
, then from above equation we have
which indicates that
is independent of
. Then, we estimate
so that (Equation 21) holds for all values
; now, rearranging (Equation 21) by
and
, we have
(23)
(23)
Thus, Equation 21 will be satisfied if
Then, the Nyquist criterion holds for , where
and
gives an estimate for the length of delay for which stability is preserved. Thus, we are now in a position to state the following theorem.
Theorem 5
If there exist in
such that
, then
is the maximum value (length of delay) of
for which
is asymptotically stable.
5. Bifurcation analysis in the presence of discrete delay
Let us consider and assume
in Equation 18. Then, separating the real and imaginary parts, we get a system of transcendental equations as
(24)
(24)
(25)
(25)
Let us consider , and hence
and
are functions of
. We want to know the change of stability of
which will occur at the values of
for
= 0 and
.
Then, the above two equations become(26)
(26)
Now, eliminating , we have
(27)
(27)
From (Equation 26), we get . The smallest
is given by
and we take it as
. In order to establish the Hopf bifurcation at
, we need to show that
at
. We differentiate Equations 24 and 25 with respect to
and setting
,
and
, we get
(28)
(28)
where
Solving (Equation 28),
we get ,
where has the same sign as that of
.
Substituting the values of and using (Equation 26), we get
Let , where
,
,
which is left side of (Equation 27) with
. Then,
and we note that
(29)
(29)
Hence, we can describe the criteria for the preservation of stability (instability) geometrically as follows:
(1) | If the polynomial | ||||
(2) | If |
(i) | For the existence of | ||||
(ii) | Since | ||||
(iii) | If | ||||
(iv) | If H(z) has two or three distinct positive real roots, then it must decrease at one root and increase at other; hence, (ii) is not satisfied. | ||||
(v) | If | ||||
(vi) | If |
and (i) will be satisfied if , i.e.
(30)
(30)
or .
Since (since
,
, and
), and
, hence
.
Thus, for Equation 30 to hold, it is sufficient that .
Table 2. A set of parametric values
Now, we can state the following theorem.
Theorem 6
If ,
, and relation (Equation 30) hold, then the stable positive equilibrium
remains stable for all
.
6. Numerical simulations
In this section, we focus our attention on the occurrence and termination of the fluctuating plankton population. We begin with a parameter set (see Table , Chatterjee & Pal, Citation2011) for which the existence condition of the coexistence equilibrium point is satisfied and the coexistence equilibrium point
is locally asymptotically stable in the form of a stable focus with eigenvalues
(cf. Figure ).
6.1. Effect of ![](//:0)
![](//:0)
We have observed that the system shows oscillatory behavior around the positive interior equilibrium for increasing the value of
with eigenvalues
(cf. Figure ). Decreasing the value of r from 1.2 to 0.3 with the same set of parametric values in Table , the system shifts to planktivorous fish free equilibrium
in the form of a stable focus with eigenvalues
(cf. Figure ). Figure (a–c) depicts the different steady-state behaviors of phytoplankton, zooplankton, and planktivorous fish in the system (Equation 1) for the parameter r. Here, we see a Hopf bifurcation point at
=1.583292 (denoted by a red star (H)) with first Lyapunov coefficient being –0.2999282. Clearly, the first Lyapunov coefficient is negative. This means that a stable limit cycle bifurcates from the equilibrium, when it loses stability.
6.2. Effects of ![](//:0)
![](//:0)
For K = 0.8, leaving all other parameters unaltered, the system exhibits oscillation around the positive interior equilibrium with eigenvalues
(cf. Figure ). Figure (a–c) depicts the different steady-state behaviors of phytoplankton, zooplankton, and planktivorous fish in the system (Equation 1) for the parameter K. Here, we see a Hopf bifurcation point at
=0.711977 (denoted by a red star (H)) with first Lyapunov coefficient being –0.3651997. Clearly, the first Lyapunov coefficient is negative. This means that a stable limit cycle bifurcates from the equilibrium, when it loses stability.
Figure 3. The figure depicts stable behavior around the planktivorous fish free equilibrium point of system (Equation 1) for decreasing r, from 1.2 to 0.3 with other parametric values as given in Table .
![Figure 3. The figure depicts stable behavior around the planktivorous fish free equilibrium point E1 of system (Equation 1) for decreasing r, from 1.2 to 0.3 with other parametric values as given in Table 2.](/cms/asset/2044eb99-a537-465c-b26f-b6ea8e3cb629/oama_a_1074337_f0003_oc.gif)
Figure 4. (a) The figure depicts different steady-state behaviors of phytoplankton for the effect of r. (b) The figure depicts different steady-state behaviors of zooplankton for the effect of r with other parametric values as given in Table . (c) The figure depicts different steady-state behaviors of planktivorous fish for the effect of r with other parametric values as given in Table .
![Figure 4. (a) The figure depicts different steady-state behaviors of phytoplankton for the effect of r. (b) The figure depicts different steady-state behaviors of zooplankton for the effect of r with other parametric values as given in Table 2. (c) The figure depicts different steady-state behaviors of planktivorous fish for the effect of r with other parametric values as given in Table 2.](/cms/asset/5ede528e-9402-42e4-b5ce-571f09486363/oama_a_1074337_f0004_oc.gif)
Figure 5. The figure depicts oscillatory behavior around the positive interior equilibrium point of system (Equation 1) for increasing K, from 0.5 to 0.8 with other parametric values as given in Table .
![Figure 5. The figure depicts oscillatory behavior around the positive interior equilibrium point E∗ of system (Equation 1) for increasing K, from 0.5 to 0.8 with other parametric values as given in Table 2.](/cms/asset/b79059bf-0082-4c0b-8fb6-99ceb0c6b027/oama_a_1074337_f0005_oc.gif)
Figure 6. (a) The figure depicts different steady-state behaviors of phytoplankton for the effect of K. (b) The figure depicts different steady-state behaviors of zooplankton for the effect of K with other parametric values as given in Table . (c) The figure depicts different steady-state behaviors of planktivorous fish for the effect of K with other parametric values as given in Table .
![Figure 6. (a) The figure depicts different steady-state behaviors of phytoplankton for the effect of K. (b) The figure depicts different steady-state behaviors of zooplankton for the effect of K with other parametric values as given in Table 2. (c) The figure depicts different steady-state behaviors of planktivorous fish for the effect of K with other parametric values as given in Table 2.](/cms/asset/fd6f84dd-0233-4d13-8bb2-d6a628d0d1f1/oama_a_1074337_f0006_oc.gif)
Figure 7. The figure depicts oscillatory behavior around the positive interior equilibrium point of system (Equation 1) for decreasing
, from 0.08 to 0.02 with same set of parametric values as use in Table .
![Figure 7. The figure depicts oscillatory behavior around the positive interior equilibrium point E∗ of system (Equation 1) for decreasing μ3, from 0.08 to 0.02 with same set of parametric values as use in Table 2.](/cms/asset/601251d7-d518-4d09-83a0-74478ec5ea89/oama_a_1074337_f0007_oc.gif)
Figure 8. The figure depicts stable behavior around the planktivorous fish free equilibrium point of system (Equation 1) for increasing
, from 0.08 to 0.3 with same set of parametric values as use in Table .
![Figure 8. The figure depicts stable behavior around the planktivorous fish free equilibrium point E1 of system (Equation 1) for increasing μ3, from 0.08 to 0.3 with same set of parametric values as use in Table 2.](/cms/asset/86243cd3-7de8-4799-9e0e-d75e099ca834/oama_a_1074337_f0008_oc.gif)
Figure 9. The figure depicts asymptotically stable behavior at positive interior equilibrium point of system (Equation 1) for r = 1.6 and
with other parametric values as given in Table ( Blue solid line). The trajectory (Black dotted line) shows stable behavior around the planktivorous fish free equilibrium point
of system (Equation 1) for r = 1.6 and
with other parametric values as given in Table .
![Figure 9. The figure depicts asymptotically stable behavior at positive interior equilibrium point E∗ of system (Equation 1) for r = 1.6 and μ3=0.12 with other parametric values as given in Table 2 ( Blue solid line). The trajectory (Black dotted line) shows stable behavior around the planktivorous fish free equilibrium point E1 of system (Equation 1) for r = 1.6 and μ3=.3 with other parametric values as given in Table 2.](/cms/asset/f47a78a2-bfbd-42cb-9f97-7778b2b221b1/oama_a_1074337_f0009_oc.gif)
Figure 10. The figure depicts stable behavior around the positive interior equilibrium point of system (Equation 1) for
and
with other parametric values as given in Table .
![Figure 10. The figure depicts stable behavior around the positive interior equilibrium point E∗ of system (Equation 1) for K=0.8 and μ3=0.12 with other parametric values as given in Table 2.](/cms/asset/431b8eae-7f61-477c-995f-ca84a6343cb8/oama_a_1074337_f0010_oc.gif)
6.3. Effects of ![](//:0)
![](//:0)
Keeping the other parameters fixed and decreasing the value of from 0.08 to 0.02, the system exhibits oscillatory behaviors around the positive interior equilibrium
with eigenvalues
(cf. Figure ). But the system shifts to planktivorous fish free equilibrium
with eigenvalues
, when
increase from 0.08 to 0.3 with the same set of parametric values as in Table (cf. Figure ).
6.4. Combined effect of ![](//:0)
and ![](//:0)
![](//:0)
It has already been shown that the system depicts oscillatory behavior around the positive interior equilibrium for
. If the mortality rate of planktivorous fish
in increased from 0.08 to 0.12, the system becomes asymptotically stable at positive interior equilibrium
in the form of a stable focus with eigenvalues
(cf. Figure ). But if the mortality rate of planktivorous fish
is increased from 0.12 to 0.3, the system shifts to planktivorous fish free equilibrium
with eigenvalues
(cf. Figure ).
6.5. Combined effect of K and ![](//:0)
![](//:0)
It has already been shown that the system depicts oscillatory behavior around the positive interior equilibrium for
(cf. Figure ). But the system shifts to asymptotically stable behavior at
in the form of a stable focus with eigenvalues
for increasing the value of
from 0.08 to 0.12 (cf. Figure ).
6.6. Effects of ![](//:0)
![](//:0)
Now, we are to introduce gestation delay, , in system (Equation 1). Keeping the other parameters fixed and increasing the value of the delay parameter
from 0 to 1, we observe that the solution of system (Equation 8) changes from stable behavior to oscillatory behavior around the positive interior equilibrium
(cf. Figure ) for same set of parametric values as in Figure . We see that decreasing the value of r from 1.2 to 1, the system shifts to stable behavior around the positive interior equilibrium (cf. Figure ). Similar cases happen when decreasing the carrying capacity and increasing mortality rate of planktivorous fish from 0.5 to 0.35 and 0.08 to 0.11, respectively, (cf. Figure and cf. Figure ) for same set of parametric values as in Figure .
6.7. Hopf bifurcation
For a clear understanding of the dynamical changes due to change in constant nutrient input, K, a bifurcation diagram is plotted with this parameter as the bifurcation parameter with other three species (cf. Figure ). We have also plotted another single parameter bifurcation diagram for r (cf. Figure ).
Now, we have plotted a two parameters bifurcation diagrams, K–r, and
(cf. Figures –). Here, we see a generalized Hopf (GH) point, where the first Lyapunov coefficient (
) vanishes but the second lyapounov coefficient is non-zero. The bifurcation point separates branches of sub and supercritical Andronov–Hopf bifurcations in the parameter plain. The Bogdanov–Takens (BT) point is the common point for the limit point curve and the curve corresponding to equilibria. Actually, at BT point, the system has an equilibrium with a double zero eigenvalues. For nearby parameter values, the system has two equilibria (a saddle and a non-saddle) which collide and disappear via a saddle-node bifurcation. The non-saddle equilibrium undergoes an Andronov–Hopf bifurcation, generating a limit cycle. Finally, a bifurcation diagram is plotted with
as the bifurcation parameter with other three species (cf. Figure ).
7. Discussion
A phytoplankton–zooplankton–fish interaction model is considered with nutrient uptake functions. Firstly, the model is studied analytically and the threshold conditions for the existence and stability of various steady states are worked out (see Table ). It is observed that for constant intrinsic growth rate r, there is a chance for the population to fluctuate. But further increase in the value of constant intrinsic growth rate r may lead to extinction of planktivorous fish population.
It is observed that the system shifts to oscillatory behavior in presence of high value of carrying capacity. We have also observed that there is a chance of extension of planktivorous fish population for high value of mortality rate of planktivorous fish. Our study indicates that low value of mortality rate of planktivorous fish may lead to oscillatory behavior around the positive interior equilibrium.
Next, we have established that the system remains locally asymptotically stable and all solutions approach toward in presence of discrete delay whenever its magnitude lies below some threshold condition. Further, we have also observed the maximum value (length of delay) of
(i.e.
) for which a locally asymptotically stable interior equilibrium
will remain asymptotically stable, where
(see Theorem 5). Moreover, we analyzed conditions for bifurcation of the positive interior equilibrium. We observed that a positive interior equilibrium remains stable if
and
for all
(see Theorem 6). Further, numerical simulations demonstrate the following conclusions:
(a) | The system exhibits dynamic instability due to higher gestation delay. | ||||
(b) | Low levels of constant intrinsic growth rate r induce stability around the positive interior equilibrium in presence of gestation delay. | ||||
(c) | Low levels of carrying capacity of phytoplankton K help the system (Equation 8) to approach stability at | ||||
(d) | High value of mortality rate of planktivorous fish prevents instability behavior in presence of gestation delay. |
Acknowledgements
The authors would like to thank the referees very much for their valuable comments and suggestions.
Additional information
Funding
Notes on contributors
Samares Pal
Our research group is headed by Samares Pal and includes Anal Chatterjee as a post-doctoral researcher. Research in the group addresses a wide range of questions broadly concerning the dynamics of plankton ecosystems under environmental conditions. The main emphasis of our work includes mathematical modeling of marine ecosystems affected by fish induced, overfishing of zooplankton. The research reported in this paper will help in studying the interspecies between phytoplankton and zooplankton in presence of fish and the subsequent changes the stability behavior.
Notes
Corresponding author: Samares Pal, Department of Mathematics, University of Kalyani, Kalyani 741235, India. E-mail: [email protected].
References
- Beretta, E., Bischi, G. I., & Solimano, F. (1991). Stability in chemostat equations with delayed nutrient recycling. Journal of Mathematical Biology, 28, 99–111. doi:10.1007/BF00171521
- Biktashev, V. N., Brindley, J., & Horwood, J. W. (2003). Phytoplankton blooms and fish recruitment rate. Journal of plankton Research, 25, 21–33. doi:10.1093/plankt/25.1.21
- Bischi, G. I. (1992). Effects of time lags on transient characteristics of a nutrient cycling model. Mathematical Biosciences, 109, 151–175. doi:10.1016/0025-5564(92)90043-V
- Chatterjee, A., & Pal, S. (2011). Effect of dilution rate on the predictability of a realistic ecosystem model with instantaneous nutrient recycling. Journal of Biological Systems, 19, 629–650. doi:10.1142/S021833901100410X
- Cushing, J. M. (1977). lntegrodifferential equations and delay models in population dynamics (Lecture notes in biomathematics, Vol. 20). Berlin: Springer-Verlag. doi:10.1007/978-3-642-93073-7
- Das, K., & Ray, S. (2008). Effect of delay on nutrient cycling in phytoplankton--zooplankton interactions in estuarine system. Ecological Modelling, 215, 69–76. doi:10.1016/j.ecolmodel.2008.02.019
- Freedman, H. I., Erbe, L. H., & Rao, V. S. H. (1986). Three species food chain models with mutual interference and time delays. Mathematical Biosciences, 80, 57–80. doi:10.1016/0025-5564(86)90067-2
- Freedman, H. I., & Rao, V. S. H. (1983). The trade-off between mutual interference and time lags in predator-prey systems. Bulletin of Mathematical Biology, 45, 991–1004. doi:10.1007/BF02458826
- Freedman, H. I., So, J., & Waltman, P. (1989). Coexistence in a model of competition in the chemostat incorporating discrete delay. SIAM Journal on Mathematical Analysis, 49, 859–870. doi:10.1137/0149050
- Gazi, N. H., & Bandyopadyay, M. (2006). Effect of time delay on a detritus-based ecosystem. International Journal of Mathematics and Mathematical Sciences, Article ID: 25619, 1–28. doi:10.1155/IJMMS/2006/25619
- Gazi, N. H., & Bandyopadhyay, M. (2008). Effect of time delay on a harvested predator-prey model. Journal of Applied Mathematics and Computing, 26, 263–280. doi:10.1007/s12190-007-0015-2
- Gopalsamy, K. (1984). Delayed responses and stability in two-species systems. The Journal of the Australian Mathematical Society Series B Applied Mathematics, 25, 473–500. doi:10.1017/S0334270000004227
- Gopalsamy, K. (1992). Stability and oscillations in delay differential equations of population dynamics (Mathematics and its applications, Vol. 74). Dordrecht: Kluwer. doi:10.1007/978-94-015-7920-9
- Khare, S., Misra, O. P., Singh, C., & Dhar, J. (2011). Role of delay on planktonic ecosystem in the presence of a toxic producing phytoplankton. International Journal of Differential Equations, Article ID: 603183, doi:10.1155/2011/603183
- Kuang, Y. (1993). Delay differential equations with applications in popular dynamics. New York, NY: Academic Press.
- Liu, S., Beretta, E., & Breda, D. (2010). Predator-prey model of Beddington--DeAngelis type with maturation and gestation delays. Nonlinear Analysis: Real World Applications, 11, 4072–4091. doi:10.1016/j.nonrwa.2010.03.013
- Maity, A., Patra, B., & Samanta, G. P. (2007). Persistence and stability in a ratio-dependent predator-prey system with delay and harvesting. Natural Resource Modeling, 20, 575–600. doi:10.1111/j.1939-7445.2007.tb00221.x
- Medvinsky, A. B., Rusakov, A. V., Bobyrev, A. E., Burmensky, V. A., Kriksunov, A. E., Nurieva, N. I., & Gonik, M. M. (2009). A conceptual mathematical model of aquatic communities in lakes Naroch and Myastro (Belarus). Biophysics, 54, 90–93. doi:10.1134/S0006350909010151
- Mukhopadhyay, B., & Bhattacharyya, R. (2008). Role of gestation delay in a plankton--fish model under stochastic fluctuations. Mathematical Biosciences, 215, 26–34. doi:10.1016/j.mbs.2008.05.007
- Mukhopadhyay, B., & Bhattacharyya, R. (2011). A stage-structured food chain model with stage dependent predation: Existence of codimension one and codimension two bifurcations. Nonlinear Analysis: Real World Applications, 12, 3056–3072. doi:10.1016/j.nonrwa.2011.05.007
- Pal, S., & Chatterjee, A. (2012). Role of constant nutrient input and mortality rate of planktivorous fish in plankton community ecosystem with instantaneous nutrient recycling. Canadian Applied Mathematics Quarterly, 20, 179–207.
- Rehim, M., & Imran, M. (2012). Dynamical analysis of a delay model of phytoplankton--zooplankton interaction. Applied Mathematical Modelling, 36, 638–647. doi:10.1016/j.apm.2011.07.018
- Ruan, S. (1995). The effect of delays on stability and persistence in plankton models. Nonlinear Analysis, Theory, Methods and Applications, 24, 575–585. doi:10.1016/0362-546X(95)93092-I
- Scheffer, M. (1991). Fish and nutrients interplay determines algel biomass: A minimal model. Oikos, 62, 271–282. Retrieved from http://www.jstor.org/stable/3545491
- Scheffer, M., Rinaldi, S., & Kuznetsov, A. Y. (2000). Effects of fish on plankton dynamics: A theoretical analysis. Canadian Journal of Fisheries and Aquatic Sciences, 57, 1208–1219.
- Singh, A., & Gakkhar, S. (2012). Analysis of delayed toxin producing phytoplankton--zooplankton system. International Journal of Modeling and Optimization, 2, 677–680. doi:10.7763/IJMO.2012.V2.208
- Upadhyay, R. K., Kumari, N., & Rai, V. (2009). Wave of chaos in a diffusive system: Generating realistic patterns of patchiness in planktonfish dynamics. Chaos, Solitons and Fractals, 40, 262–276. doi:10.1016/j.chaos.2007.07.078
- Yongzhen, P., Min, G., & Changguo, L. (2011). A delay digestion process with application in a three-species ecosystem. Communications in Nonlinear Science and Numerical Simulation, 16, 4365–4378. doi:10.1016/j.cnsns.2011.03.018