Plankton nutrient interaction model with effect of toxin in presence of modified traditional Holling Type II functional response

A three-dimensional plankton-nutrient interaction model is proposed and analysed which mediated by a toxin-determined functional response. The new functional response is a modification of the traditional Holling Type II functional response by explicitly including a reduction in the consumption of phytoplankton by the zooplankton due to chemical defenses. Our analysis leads to different thresholds in terms of model parameters to find out different steady-states behaviour. It is found that constant nutrient input and dilution rate of nutrient influence the plankton ecosystem model and maintain stability around the coexistent equilibrium. Our observations indicate that if the constant nutrient input crosses a certain critical value, the system enters into Hopf bifurcation. In addition, we have studied the direction of Hopf bifurcation by applying the normal form method. The maximal amount of toxin of phytoplankton species plays a crucial role to change the steady-state behaviour. Computer simulations illustrate the results.


Introduction
In all marine ecosystems, the interactions between nutrient and plankton are important. These interactions are partially understood due to their complexities. Several scientists have discussed the role of competition in phytoplankton population for the occurrence and control of plankton bloom and the effect of toxin on nutrient-phytoplanktonzooplankton (NPZ) model (Bairagi, Pal, Chatterjee, and Chattopadhyay, 2008;Chakraborty, Bhattacharya, Feudel, and Chattopadhyay, 2012;Chatterjee, Pal, and Chatterjee, 2011;Mukhopadhyay and Bhattacharyya, 2005;Pal, Chatterjee, and Chattopadhyay, 2007). The adverse effect on plankton ecosystem model and downward movement of the equilibrium level of zooplankton biomass due to presence of toxin have been examined by Khare, Misra, and Dhar (2010). As stated by Banerjee and Venturino (2011), the toxic phytoplankton don't lead the zooplankton population to become extinct in plankton marine ecosystem. Jang, Baglama, and Rick (2006) investigated that when phytoplankton and zooplankton can uptake toxin separately, the growth rates of plankton population decrease. But Chattopadhyay, Sarkar, and Pal (2004) showed that there may not be release of toxic chemicals by all the plankton species. Ruan (1993Ruan ( , 2001 has pointed out the importance of nutrient for growth of plankton population in phytoplankton-zooplankton interaction model. The effects of multi nutrient and harmful algal bloom on NPZ model have been investigated by some researchers with variable * Corresponding author. Email: samaresp@yahoo.co.in outcomes (Mitra, 2006;Mitra and Flynn, 2006). The role of delay on plankton ecosystem in the presence of a toxin producing phytoplankton has been examined by Khare, Misra, Singh, and Dhar (2011). This study has also shown the effect of delay in the digestion of nutrient by phytoplankton biomass. In a phytoplankton-zooplankton-fish model with fish affecting the migrations of zooplankton, in turn affects the migrations of motile phytoplankton, as studied by Samanta, Chowdhury, and Chattopadhyay (2013).
Intake of forage is raised by the nutrient (nitrogen, phosphorous, carbohydrates, etc.) as herbivores require them and often are limiting in food supply. The herbivores either stops eating poisonous food or die, as by satiating the detoxification system of herbivores, herbivory is reduced by toxins (Li, Feng, Swihart, Bryant, and Huntly, 2006). Consumption of toxic food stops before the occurrence of lethal poisoning and hence poisoning is not so common. The zooplankton functional response, the linker of zooplankton to phytoplankton population, is a key component in phytoplankton-zooplankton interaction model. The instantaneous change in the rate of prey intake by herbivore in response to changing prey biomass is referred as the said functional response. Holling Type II response, featured by a monotonic increase in intake and asymptotic towards a maximum rate of biomass consumption, is a commonly used functional response. The traditional Holling type II functional response assuming the growth rate of herbivore as a monotonically increasing function of plant density is one of the most commonly used descriptions of the plant consumption by herbivores (Liu, Feng, Zhu, and L DeAngelis, 2008).
However, this will not be appropriate if the chemical defense of phytoplankton is considered, where the negative effect of phytoplankton toxin on zooplankton can lead to a decrease in the growth rate when the phytoplankton density is high. The Holling Type II functional responses is characterized by a monotonic increase in phytoplankton consumption as phytoplankton biomass increases asymptotically towards a maximum consumption rate, at which the zooplankton is satiated. To examine the effect of toxicity on feeding rates, we constructed a simple toxin-determined functional response model (TDFRM) with phytoplankton-zooplankton species.
The toxin-determined functional response is a modification of the traditional Holling Type II response to include the negative effect of toxin on zooplankton growth, which can overwhelm the positive effect of biomass ingestion at sufficiently high phytoplankton toxicant concentrations. The model adds one parameter to the Holling Type 2 response; G 1 , which describes the negative effect that phytoplankton toxins have on the zooplankton population's time-dependent per capita intake of phytoplankton biomass. This addition allows one to dynamically model how toxin-determined zooplankton affects when it consumes toxic phytoplankton in plankton dynamics. Thus, the model is a potential tool for resource managers working in ecosystems containing toxic phytoplankton.
As observed by Swihart, L DeAngelis, Feng, and Bryant (2009), three analytical modifications should be considered when incorporating the effects of toxins on plant-herbivore dynamics. Toxin-mediated functional responses should (a) explicitly account for the negative effects of phytoplankton toxins on zooplankton growth; (b) permit herbivores to regulate intake of toxins; and (c) allow for intake of multiple phytoplankton that are detoxified with independent pathways. Here an open system is considered with three interacting components consisting of phytoplankton (P), zooplankton(Z) and dissolved limiting nutrient (N ). In this paper a nutrient-phytoplankton-zooplankton model with the effect of toxin producing phytoplankton is described. The stability of equilibrium point is described. We have also derived the conditions for instability of the system around the interior equilibrium and direction of Hopf bifurcation. Numerical simulations under a set of parameter values have been performed to support our analytical result.

Formulation of the mathematical model
Let N (t) be the concentration of the nutrient at time t. In addition, P(t) and Z(t) be the concentration of phytoplankton and zooplankton population, respectively, at time t. N 0 be the constant input of nutrient concentration, D is the dilution rate of nutrient (Fredrickson and Stephanopoulos, 1981). The constant 1/D has the physical dimension of a time and represents the average time that nutrient and waste products spend in the system (Smith, 1981). Let α 1 and α 2 be the nutrient uptake rate for the phytoplankton population and conversion rate of nutrient for the growth of phytoplankton population, respectively (α 2 ≤ α 1 ). e 1 is the rate of encounter per unit of phytoplankton species. The parameter e 1 is the resource encounter rate, which depends on the movement velocity of the zooplankton and its radius of detection of food items. The parameter σ 1 (0 < σ 1 ≤ 1) is the fraction of food biomass encountered that the zooplankton ingests, while h is the handling time per unit biomass, which incorporates the time required for the digestive tract to handle the item.
The second part of the toxin determined functional response, that is, 1 − f (N )/4G 1 , accounts for the negative effect of toxin, where G 1 = M 1 /T 1 . The parameter M 1 is a measure of the maximum amount of toxicant per unit time that the zooplankton can tolerate, and T 1 is the amount of toxin per unit phytoplankton biomass. G 1 decreases when the concentration of toxin tolerated by the zooplankton decreases (high toxicity) or the concentration of toxin in the phytoplankton biomass increases. Therefore, the smaller the value of G 1 , the larger the effect the toxin has on intake.
If we assume that M 1 is the same for all individuals within a zooplankton population, then the variation in G 1 will result from the different values of T 1 of the phytoplankton species encountered. The factor 4 is just a multiplier that simplifies the form of the peak value of C(P) as a function of P. Let D 1 be the mortality rate of the phytoplankton population and D 2 be the mortality rate of the zooplankton population. Let β 1 be the maximal zooplankton conversion rate. We have considered Holling type II functional response to describe the grazing phenomena with K 1 as half saturation constant. Now following the above assumption, we consider the mathematical model where We take the initial conditions for the system (1) as

Positive invariance
Due to the constant supply of nutrients N 0 , the system (1) is not homogeneous. By considering together with X (0) = X 0 ∈ R + 3 . It is easy to verify that whenever choosing 3 for all t > 0. This ensures that the solution remains within the positive orthant, ensuring the biological wellposedness of the system.

Boundedness of the system
THEOREM 1 All the solutions of (1) are bounded.
The proof is obvious.

Equilibria
The system (1) possesses the following equilibria: we calculate third equation of system (1) by putting zero in right hand side we get For Case 2 and Case 3, we get unique positive real root and Thus the conditions for the existence of the unique interior equilibrium point E * (N * , P * , Z * ) are given by, N * > 0, P * > 0, Z * > 0 with and Case 4: If A > 0, B < 0, the quadratic equation has two positive real roots. Let Then and

Criterion for the extinction of phytoplankton
THEOREM 2 Let the inequality Theorem (2) demonstrates that if the ratio of maximal nutrient conversion rate for the growth of phytoplankton and the loss rate of the phytoplankton is less than unity then the phytoplankton population will eliminate from the system. The proof is obvious and hence omitted.

Eigenvalue analysis
In this section, we investigate the local stability behaviour of the system (1) around each of the equilibria by computing the jacobian matrix. LetĒ = (N ,P,Z) be any arbitrary equilibrium. Then the jacobian matrix forĒ is given bȳ Proof By computing the jacobian matrix V 0 for the equilibrium E 0 of the system (1) we get three eigen- . Clearly E 0 is asymptotically stable if and only if R 0 < 1. When R 0 > 1, plankton free steady state is unstable and a feasible zooplankton free steadystate E 1 exists.
LEMMA 2 There exists a feasible zooplankton free steady-state E 1 of the system (1) which is unstable if Proof Further the eigenvalues of the jacobian matrix V 1 around the equilibrium E 1 of the system (1) are λ 1 , λ 2 which are the roots of the equation Clearly λ 1 and λ 2 have negative real parts because all the coefficients of above quadratic equation is positive. Therefore λ 3 > 0 if condition (9) is satisfied and hence E 1 is unstable.
We summarized the analytical results in Table 1. Stability analysis of the positive interior equilibrium of the system (1).
The jacobian matrix of system (1) around the positive equilibrium E * = (N * , P * , Z * ) is The characteristic equation is where Then the system becomes locally asymptotically stable around E * . Thus depending upon system parameters, the system may exhibit stable or unstable behaviour in this case.

Remark 1
The system could have a Hopf-bifurcation at the coexistence equilibrium if the following two conditions are satisfied,

Hopf bifurcation at coexistence
In fact, the bifurcation occurs when the characteristic Equation (10) has to have two purely imaginary roots ±η 2 i, where i denotes the imaginary unit, in addition to the real one e. It follows that it can be rewritten as (y 2 + η 2 2 )(y − e) = 0. Expanding, we find y 3 − ey 2 + η 2 2 y − eη 2 2 = 0. Comparing the coefficients with those of Equation (10) the first condition (11) follows. The transversality condition, the second (11), is obtained observing that the eigenvalues of the characteristic equation, of the form λ i = u i + iη i must satisfy the transversality condition du i /dN 0 | N 0 =N 0 c = 0. Substituting λ i into the characteristic equation, separating the real and imaginary parts and eliminating η, using Differentiating this cubic with respect to N 0 , we find providing the second condition (11). This ensures that the system (1) has a Hopf-bifurcation around the positive interior equilibrium E * . An analytic verification of Equation (11) appears in these conditions a very difficult task, as the coefficients depend on the population levels at equilibrium E * . However, these remarks help in guiding the numerical simulations, that in fact reveal the existence of sustained population oscillations.
THEOREM 3 The parameter μ 22 determines the direction of the Hopf-bifurcation. If μ 22 > 0 (< 0) then the Hopfbifurcation is supercritical (subcritical) and the bifurcating periodic solutions exists for N 0 > N 0 c . The stability and the period of the bifurcating periodic solutions are, respectively, determined by the parameters β 2 and τ 2 defined in the proof. The solutions are orbitally stable (unstable) if β 2 < 0 (> 0) and the period increases (decreases) if τ 2 > 0 (< 0).
Proof We just outline the sketch of the proof, because the computational details are very complicated. The method is based on the normal form theory (Hassard, Kazarinoff, and Wan, 1981). In the conditions of the Hopf bifurcation given in the above Remark 1, let us denote by a bar the system parameters; the eigenvector corresponding to the eigenvalue σ = iη 2 is The eigenvector corresponding to the eigenvalue e is Using the transformation where as a shorthand we have introduced the following quantities: The equilibrium point of the new system (12) is now the origin. At it, the Jacobian of (12) simplifies since several entries vanish: Further, the following auxiliary quantities can be explicitly calculated in terms of the system parameters, but we omit the explicit formulae in view of their excessive length.

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, reference Feng, Qiu, Liu, and L DeAngelis, 2011) for which the existence criterion at equilibrium of coexistence E * = (2.086, 0.3929, 2.314) is satisfied. In this case coexistence equilibrium point E * is locally asymptotically stable in the form of a stable focus (cf. Figure 1).

Effects of N 0
Increasing the value of constant nutrient input N 0 from 3 to 5.5, the equilibrium switches from stable to oscillatory behaviour around the E * for same set of parametric values in Table 2 (cf. Figure 2). But in presence of low value of constant nutrient input, N 0 = 0.04, system goes to zooplankton free equilibrium E 1 (cf. Figure 3).

Effects of D
For D = 0.85, leaving all other parameters unaltered, the system exhibits oscillation around the positive interior equilibrium E * (cf. Figure 4). But decreasing the value of D from 0.3 to 0.02, the system goes to zooplankton free equilibrium E 1 (cf. Figure 5).

Effects of G 1
We have observed that the system shows oscillatory behaviour around the positive interior equilibrium E * for low value of G 1 = 0.25 (cf. Figure 6). If G 1 is further decreasing from 0.25 to 0.15, the system becomes stable at zooplankton free equilibrium E 1 in the form of a stable focus (cf. Figure 7).

Effect of D 2
If the mortality rate of zooplankton D 2 is increased from 0.04 to 0.4, the system goes to zooplankton free equilibrium at E 1 (cf. Figure 8).   Table 2.  Table 2.  Table 2.  Table 2.

Hopf-bifurcation
System behaviour can be demonstrated more prominently if we plot the bifurcation diagram for constant  Table 2.  Table 2.  Table 2.  Table 2. nutrient input, N 0 , as a bifurcation parameter with other three species. Here we see a Hopf bifurcation at N 0 = 5.225469 with first Lyapunov coefficient −1.986688 e −003 (cf. Figure 9(a)-(c). Clearly the value of first Lyapunov coefficient is negative. This means that a stable limit cycle bifurcates from the equilibrium, when it looses stability. In Figure 9 We have plotted another diagram with dilution rate of nutrient, D, with other three species (cf. Figure 10(a)-(c)). Again we see that a Hopf bifurcation occurs at D = 0.807954 with first Lyapunov coefficient is −2.295238 e −003 . This means that a stable limit cycle bifurcation from the equilibrium, when it looses stability. In Figure 10 Next we have plotted a two parameters bifurcation diagram N 0 − D for showing a clear picture of stable and unstable zone in this system (cf. Figure 11(a)). Finally, we have plotted a two parameters bifurcation diagram, N 0 − G 1 . 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 common point for the limit point curve and the curve corresponding to equilibria. The system has an equilibrium with a double zero eigenvalues at BT point. 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 (cf. Figure 11(b)).  Table 2.  Table 2.  Table 2 of the system (2.1). (b) The two parameter bifurcation diagram for N 0 -G 1 with all parametric values as given in Table 2 of the system (2.1).

Discussion
A nutrient-plankton interaction model in presence of toxin producing phytoplankton is considered with Holling type II nutrient uptake function. We study the dynamics of the system (1) with a functional response mediated by phytoplankton toxicity. This toxin-determined functional response model (TDFRM) is a modification of the traditional consumer resource model with Holling Type II functional response. 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. Note that our recent model is different from other models (Feng et al., 2011;Li et al., 2006;Liu et al., 2008). Our discussion has focused on nutrient-plankton ecosystem apart from plant-herbivore system. Recent modelling efforts implicate nutrient input as potentially key drivers of change in plankton ecosystem. The simulation result confirmed the existence of a stable periodic solution when the stable interior equilibrium loses its stability. It is observed that there is a chance for the fluctuating population for high value of constant nutrient input. Thus, our analysis shows that if the constant nutrient input crosses a critical value, say N 0 c , the system (1) enters into Hopf bifurcation that induces stable limit cycle periodic solutions around the positive equilibrium. Then all the populations start to oscillate. But for low values of constant nutrient input it may lead to extinction of zooplankton population. Moreover, change in stability of interior equilibrium is demonstrated using several key parameters including the dilution rate of nutrient, mortality rate of zooplankton and phytoplankton toxicity.
Our simulations also helped us to indicate that, in presence of high dilution rate of nutrient there is chance for oscillatory behaviour of plankton population. But at low value of dilution rate of nutrient, zooplankton goes to extinction.
Also our observation indicates that there is a chance for oscillatory behaviour in presence of low value of toxicity level of phytoplankton species. Further, decreasing the value of toxicity level of phytoplankton species, the zooplankton goes to extinction. Similar case arises for high value of mortality rate of zooplankton population.
Thus our analysis suggests that to control the fluctuating population and to maintain stability around the coexistence of equilibrium, we have to control the constant nutrient input, dilution rate of nutrient and toxicity level of phytoplankton population.