Dynamics of the interaction of plankton and planktivorous fish with delay

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.


ABOUT THE AUTHORS
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.

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.

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 (1991) and Pal and Chatterjee (2012). 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 (2000) 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 (2008). Realistic patterns of patchiness in plankton-fish dynamics have been studied by Upadhyay, Kumari, and Rai (2009). 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 (2003), 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. (2009).
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 (1992), 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, 1977;Freedman & Rao, 1983). Chemostat-type models incorporating discrete delays have been investigated by Freedman, So, and Waltman (1989). Also, delayed models in population biology have been discussed in Gopalsamy (1984), Kuang (1993), Khare, Misra, Singh, andDhar (2011) andLiu, Beretta, andBreda (2010). A discrete time delay term to the model of Beretta, Bischi, and Solimano (1991) was introduced in Ruan (1995). 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 (2008). Rehim and Imran (2012) 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, 2012). In particular, Gazi and Bandyopadhyay (2006) 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 (2008), Maity, Patra, and Samanta (2007) and Yongzhen, Min, and Changguo (2011).
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 nondelay models. Direction of Hopf bifurcation is discussed. Numerical simulations under a set of parameter values have been performed to support our analytical results.

The mathematical model
Let P(t) be the concentration of the phytoplankton at time t with carrying capacity K and constant intrinsic growth rate r. Let Z(t) be the concentration zooplankton biomass and F(t) be the concentration of planktivorous fish biomass present at any instant of time t.
Let 1 be the maximal zooplankton ingestion rate and 2 maximal zooplankton conversion rate for the growth of zooplankton, respectively ( 2 ≤ 1 ). Let 1 be the maximal planktivorous fish ingestion rate and 2 ( 2 ≤ 1 ) be the maximal planktivorous fish conversion rate due to grazing of herbivorous zooplankton. Further, 1 , 2 , 3 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 K 1 and K 2 as half saturation constant.
The basic equations with all of the parameters are: The system (Equation 1) can be written in compact form as Ẋ = J(X). Its Jacobian is

Positive invariance
The system (Equation 1) is not homogeneous. However, it is easy to check that whenever choosing X(0) ∈ R + 3 with X i = 0, for i = 1, 2, 3, then J i (X) | X i = 0 ≥ 0. This ensures that the solution remains within the positive orthant, ensuring the biological well posedness of the system.

Boundedness of the system
Theorem 1 All the solutions of (Equation 1) are ultimately bounded.
Proof We define a function For the large value of t, we get dw dt ≤ rP − D 0 w, where D 0 = min{ 1 , 2 , 3 }.

Planktivorous fish free equilibrium
At E 1 , the population levels are This equilibrium is feasible if it is satisfying At E 1 , the Jacobian (Equation 2) factorizes to give one explicit eigenvalue 2 Z 1 (K 2 + Z 1 ) −1 − 3 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 E 1 is then ensured by

The coexistence equilibrium
The coexistence equilibrium E * = (P * , Z * , F * ) cannot be found explicitly since Z * = 3 K 2 Thus, the condition for the existence of the interior equilibrium point E * (P * , Z * , F * ) is given by, P * > 0, Z * > 0, F * > 0. For feasibility, we certainly need 2 > 3 and 0 < r < min
Therefore, according the Routh-Hurwitz criteria, all roots of the above cubic equation have negative real parts satisfying A 1 > 0, A 3 > 0, and A 1 A 2 − A 3 > 0. Then, the system becomes locally asymptotically stable around E * . Thus, depending upon system parameters, the system may exhibit stable or unstable behavior in this case.
The analytical results are summarized in Table 1.

Hopf bifurcation at coexistence
Theorem 2 When the carrying capacity K crosses a critical value, say K * , 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 bifurca- The transversality condition, the third (iii), is obtained observing that the eigenvalues of the characteristic equation, of the form i = u i + iv i , must satisfy the transversality condition We will now verify the Hopf bifurcation condition (iii), putting = u + iv in the characteristic equation, we get On separating the real and imaginary parts and eliminating v, we get It is clear from the above that u( ≠ 0, providing the third condition (iii).
This ensures that the above system has a Hopf bifurcation around the positive interior equilibrium E * .
Next, we perform detailed analysis of the bifurcation solutions to study the nature of Hopf bifurcation. We start by rewriting system (Equation 1) Now, let the normalized eigenvectors of M and M T corresponding to the eigenvalues i and −i be q and p, respectively, where . Let us define the following multi-linear functions for the three With these, F(X) will be of the form F(X) = 1 2 B(X, X) + 1 2 C(X, X, X).
Solving the corresponding linear system, we get where Again, where E = (2i − n 11 )(2i − n 12 ) − n 12 n 21 . Therefore, Using the invariant expression for the first lyapunov coefficient l 1 (0), we get the first lyapunov coefficients as Mukhopadhyay and Bhattacharyya (2011).
Now, if l 1 (0) < 0, the Hopf bifurcation will be of supercritical nature and there will be emergence of stable limit cycles.

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.
The system (Equation 8) has to be analyzed with the following initial conditions,

Qualitative analysis of the model with time delay
Remark The characteristic equation for the variational matrix V 0 about the plankton free steady states E 0 and E 01 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 E 0 and E 01 .
To find conditions for the other local asymptotic stability of system (Equation 8), we use the following theorem from Gopalsamy (1992).
Theorem 3 If D E ( , ) = 0 denotes the characteristic equation, then a set of necessary and sufficient conditions for the equilibrium to be asymptotically stable for all ≥ 0 are (i) The real parts of all the roots of D E ( , 0) = 0 are negative.
(ii) For all real and any ≥ 0, The characteristic equation for the variational matrix V 1 about E 1 takes the following form .
Here, D E 1 ( , 0) = 0 has roots with negative real parts, provided system (Equation 8) is locally asymptotically stable about equilibrium E 1 (for the condition, see Section 3.3.3).
and separating the real and imaginary parts, we get −l 2 + n = −f 2 cos + h cos + g cos , − 3 + m = f 2 sin − h sin + g cos .
Squaring and adding the above two equations, we get Sufficient conditions for the non-existence of a real number satisfying D E 1 (i , ) = 0 can be writ-

Estimation of the length of delay to preserve stability at E 1
In this section, we assume that in the absence of delays, E 1 is locally asymptotically stable. Let u 1 (t), v 1 (t), and w 1 (t) be the respective linearized variables of this model. The system (Equation 8) becomes Let û 1 (L), v 1 (L), and ŵ 1 (L) be the Laplace transform of u 1 (t), v 1 (t), and w 1 (t), respectively. Taking the Laplace transformation of system (Equation 11), we have The inverse Laplace transform of û 1 (L), v 1 (L), and ŵ 1 (L) will have terms which exponentially increase with time, if û 1 (L), v 1 (L), and ŵ 1 (L) have pole with positive real parts. Since E 1 needs to be locally asymptotically stable, it is necessary and sufficient that all poles of û 1 (L), v 1 (L), and ŵ 1 (L) have negative real parts. We shall employ the Nyquist criterion which states that if L is the arc length of a curve encircling the right half-plane, the curve of û 1 (L), v 1 (L), and ŵ 1 (L) will encircle the origin a number of times equal to the difference between the number of poles and the number of zeroes of û 1 (L), v 1 (L), and ŵ 1 (L) in the right half-plane.
ReG(iv 2 ) = 0, ImG(iv 2 ) > 0, To get an estimation on the length of delay, we utilize the following conditions, Therefore, E 1 will be stable if the above inequality holds at v = v 2 , where v 2 is the first positive root of Equation 14. We shall now estimate an upper bound v ++ of v 2 , which would be independent of . Then, we estimate so that (Equation 15) holds for all 0 ≤ v ≤ v ++ , and hence in particular at v = v 2 .
, then from (Equation 16) Then, the Nyquist criterion holds for 0 ≤ ≤ 1 , where 1 = 1 gives an estimate for the length of delay for which stability is preserved.
Theorem 4 If there exist in 0 ≤ ≤ 1 such that A 2 + B + C > 0, then 1 is the maximum value (length of delay) of for which E 1 is asymptotically stable.

Qualitative analysis of the model at E * with gestation time delay ( ≠ 0)
The characteristic equation for the variational matrix V * about E * takes the following form: −lv 2 + n ≤ fv 2 + h + gv.
(17) Let D E * (ī, ) = 0 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 D E * = 0 can be written as ̄6 +Q 1̄4 +Q 2̄2 +Q 3 > 0, which can be transformed to Therefore, a sufficient condition for E * to be stable is (a) Q 1 > 0 and (b) Q 3 >Q 2 2 4Q 1 .

Estimation of the length of delay to preserve stability at E *
In this section, we assume that in the absence of delays, E * is locally asymptotically stable. Let u(t), v(t), and w(t) be the respective linearized variables of this model. The system (Equation 8) becomes Let ū(L), v(L), and w(L) be the Laplace transform of u(t), v(t), and w(t), respectively. Taking the Laplace transformation of system (Equation 19), we havē Then, E * is locally asymptotically stable if ReF(iv 0 ) = 0 and ImF(iv 0 ) > 0 i.e.
where v 0 be the smallest positive root of Re(iv 0 ) = 0.
To get an estimation on the length of delay, we use the following conditions, Therefore, E * will be stable if above inequality holds at v = v 0 , where v 0 is the first positive root of Equation 20. We shall now estimate an upper bound v + of v 0 , which would be independent of . Then, we estimate so that (Equation 21) holds for all 0 ≤ v ≤ v + , and hence in particular at v = v 0 .
Thus, the unique positive solution of (f Then, the Nyquist criterion holds for 0 ≤ ≤ + , where + = 1 2A (−B + √ B 2 + 4A C ) 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 0 ≤ ≤ + such that A 2 + B + C > 0, then + is the maximum value (length of delay) of for which E * is asymptotically stable.

Bifurcation analysis in the presence of discrete delay
Let us consider ≠ 0 and assume 1 = + i in Equation 18. Then, separating the real and imaginary parts, we get a system of transcendental equations as Let us consider 1 , and hence and are functions of . We want to know the change of stability of E * which will occur at the values of =̂ for = 0 and ≠ 0.
(2) If H(z)is decreasing (increasing) at all its positives roots, stability (instability) is preserved.
We note the following facts: (i) For the existence of ̂, H(z) must have at least one positive real root.
(ii) Since H(z) is cubic in z , (iii) If H(z) has a unique positive root, then it must increase at that point to satisfy (ii).
(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 B 3 < 0, then H(z) has only one root.
(vi) If B 1 > 0, B 3 > 0, and B 2 < 0, then (i) will be satisfied. Now, if B 1 > 0, B 3 > 0, and B 2 < 0, then the minimum of H(z) will exist at z min = Thus, for Equation 30 to hold, it is sufficient that 27B 3 + 2B 3 . Now, we can state the following theorem.

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 2, Chatterjee & Pal, 2011) for which the existence condition of the coexistence equilibrium point E * is satisfied and the coexistence equilibrium point E * = (0.40, 0.10, 0.54) is locally asymptotically stable in the form of a stable focus with eigenvalues −0.8003, −0.0129 ± i0.2064 (cf. Figure 1).

Effect of r
We have observed that the system shows oscillatory behavior around the positive interior equilibrium E * for increasing the value of r = 1.6 with eigenvalues −1.2305, 0.0013 ± i0.0.2057 (cf. Figure 2).
Decreasing the value of r from 1.2 to 0.3 with the same set of parametric values in Table 2, the system shifts to planktivorous fish free equilibrium E 1 = (0.0250, 0.0766, 0) in the form of a stable focus with eigenvalues −0.0173, −0.0026 ± i0.1188 (cf. Figure 3). Figure 4(a-c) depicts the different steadystate behaviors of phytoplankton, zooplankton, and planktivorous fish in the system (Equation 1) for the parameter r. Here, we see a Hopf bifurcation point at r 0 c =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.

Effects of K
For K = 0.8, leaving all other parameters unaltered, the system exhibits oscillation around the positive interior equilibrium E * with eigenvalues −0.8599, 0.0069 ± i0.2370 (cf. Figure 5). Figure 6(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 K * = 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.

Effects of 3
Keeping the other parameters fixed and decreasing the value of 3 from 0.08 to 0.02, the system exhibits oscillatory behaviors around the positive interior equilibrium E * with eigenvalues −1.0833, 0.0002 ± i0.1084 (cf. Figure 7). But the system shifts to planktivorous fish free equilibrium  Table 2 (cf. Figure 8).

Combined effect of r and 3
It has already been shown that the system depicts oscillatory behavior around the positive interior equilibrium E * for r = 1.6. If the mortality rate of planktivorous fish 3 in increased from 0.08 to 0.12, the  Figure 9). But if the mortality rate of planktivorous fish 3 is increased from 0.12 to 0.3, the system shifts to planktivorous fish free equilibrium E 1 = (0.025, 0.4625, 0) with eigenvalues −0.0284, −0.0104 ± i0.2918 (cf. Figure 9).

Combined effect of K and 3
It has already been shown that the system depicts oscillatory behavior around the positive interior equilibrium E * for K = 0.8 (cf. Figure 5). But the system shifts to asymptotically stable behavior at E * = (0.5960, 0.1590, 0.7282) in the form of a stable focus with eigenvalues −0.6461, −0.0065 ± i0.2876 for increasing the value of 3 from 0.08 to 0.12 (cf. Figure 10).

Effects of
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  Figure 11) for same set of parametric values as in Figure 1. We see that decreasing the value of r from 1.2 to 1, the system shifts to stable behavior around the positive interior  Figure 12). 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 13 and cf. Figure 14) for same set of parametric values as in Figure 11.

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 15). We have also plotted another single parameter bifurcation diagram for r (cf. Figure 16). Now, we have plotted a two parameters bifurcation diagrams, K-r, K − 3 and r − 3 (cf. Figures 17-19). Here, we see a generalized Hopf (GH) point, where the first Lyapunov coefficient (l 1 ) 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  20).

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 1). 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 E * 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 E * will remain asymptotically stable, where (d) High value of mortality rate of planktivorous fish prevents instability behavior in presence of gestation delay.
Throughout the article of this model (analytically and numerically), an attempt has been made to search for a suitable way to control the system and maintain a stable coexistence between all the species. It has been observed that to control the fracturing population and to maintain stability around the coexistence equilibrium, we have to balance the rate of carrying capacity of phytoplankton, mortality rate of planktivorous fish population, and constant intrinsic growth rate. Theorem 5). Moreover, we analyzed conditions for bifurcation of the positive interior equilibrium. We observed that a positive interior equilibrium remains stable if B 1 > 0 and B 3 > 0 for all > 0 (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.