Dynamical Complexity in Some Ecological Models : Effect of Toxin Production by Phytoplankton

Abstract. We investigate dynamical complexities in two types of chaotic tri-trophic aquatic food-chain model systems representing a real situation in the mar ine environment. Phytoplankton produce chemical substances known as toxins to reduce grazing pressure by zooplankton [1]. The role of toxin producing phytoplankton (TPP) on the chaotic behavior in these food chain systems is investigated. Holling type I, II, and III functional response forms are considered to study the interference between phy toplankton and zooplankton populations in the presence of toxic chemical. Our study show s that chaotic dynamics is robust to changes in the rates of toxin release as well a s the toxin release functions. The present study also reveals that the rate of toxin p r duction by toxin producing phytoplankton plays an important role in controlling oscillation s in the plankton system. The different mortality functions of zooplankton due to to xin producing phytoplankton have significant influence in controlling oscillations, coexiste nce, survival or extinction of the zooplankton population. Further studies are needed to a scert in if this defence mechanism suppresses chaotic dynamics in model aquatic sys tem .


Introduction
It is a challenge to understand the dynamical complexities of an ecological system.Field and laboratory studies are hard to design and implement.One of the main reasons for this is the fact that there has not been a theory which can guide the experiments and one finds hardly any experiment to corroborate the predictions of an ecological model.In recent years, the success of efforts by an interdisciplinary team [2] has made us believe that complex dynamics in ecological data could be the result of simple rules.The tribolium project indeed corroborates May's hypothesis [3].This encourages us to repose more faith in models based on a set of established ecological principles.The resulting models are non-linear and deterministic and can be explored using concepts and notions from non-linear dynamics.
Mathematical models have been designed and studied since the pioneering work of Sir Robert May [3,4].In population ecology, the practice has been to design either the difference or differential equation models.The difference equation models describe the evolution of biological populations with non-overlapping generations.On the other hand, the differential equation models correspond to populations with overlapping generations.The success of these models depends on the underlying ecological principles.The identification of general ecological principles in itself is a great challenge.Of late, there have been some thoughts on this issue [5,6].Two ecological principles that form the skeleton of the model systems that we study in this paper are: 1.A specialist predator population decays exponentially fast in the absence of its lone prey.
2. The generalist predator switches to an alternative food option as and when it faces difficulty to find its favorite prey.The per capita growth of a generalist predator is limited by dependence on its favorite preys and severity of this limitation is inversely proportional to per capita availability of preys at any instant of time.
In the present paper, we study and compare the dynamical complexity of two nonlinear deterministic prey-predator models of aquatic ecosystems.The first model has both kinds of predators: specialist as well as generalist.The second one has only one kind of predator; that is the specialist predator.In our earlier work [7], similar type of model systems are considered for terrestrial ecosystems and suggested that the biology of the top predator would be a crucial factor for the determination of food chain dynamics.In the proposed work, we modify the model system of Upadhyay and Rai [7,8] and Hastings and Powell [9] for aquatic environment by introducing the toxin liberation process of TPP population for which the mortality of zooplankton increases.It should be noted that we assume these populations to be submerged in a homogeneous environment, therefore, diffusive processes present in marine environments are ignored [10].To observe the role of TPP, Holling type I, II and III functional response forms are considered for the description of consumption of prey by its predator.We investigate the dynamical complexity in these model systems with the help of bifurcation study.In order to defend themselves against grazing by zooplankton, phytoplankton release toxins.These toxins weaken the rate at which zooplankton graze [1].Recently, it has been shown that the toxin production by phytoplankton suppresses chaotic dynamics [11,12].These authors added an extramortality term to the rate equation for the middle predator in the Hastings and Powell model [9] and Upadhyay and Rai model [8] to incorporate the effect of toxin release by phytoplankton.We study the dynamical complexities of such model systems in detail.One of the main objectives of our study is to investigate how dynamical complexities of a model given by Upadhyay and Rai [8] and Hasting and Powell [9] changes their basic character in response to different types of toxin release functions.
It is observed that toxin producing phytoplankton (TPP) population do not release toxic chemical always, release only in the presence of dense zooplankton population around it.This phenomenon has been included in the interaction of Holling type I and type II functional form respectively.But as liberation of toxin reduces the growth of zooplankton, causes substantial mortality of zooplankton and in this period TPP population is not easily accessible, hence a more common and intuitively obvious choice is of Holling type II or type III functional form to describe the grazing phenomena in the presence of toxic substances [13,14].The Holling type II and III predation form is an obvious choice to represent the hunting behavior of predator [15,16].A realistic description of zooplankton grazing with functional responses to phytoplankton abundance was introduced by Ivlev [17] with a modification by Mayzaud and Poulet [18].Holling type response term are also in use [19,20].Edwards and Brindley [21] observed that the choice of functional form and mortality of zooplankton has a major influence in the dynamics of excitable nature of blooms.Now we propose two model of TPP-Zooplankton -fish/molluscs interaction.
The plan of the paper is as follows: in the next section, we present details of these two model systems.The third section describes stability analysis and Hopf-bifurcation for the given model systems.Numerical results are presented in Section 4. Section 5 concludes this paper.

Model systems
Consider a situation where TPP population (prey) of size x is predated by individuals of specialist predator zooplankton population y.This zooplankton population, in turn, serves as a favorite food for the generalist predator molluscs population of size z.This interaction is represented by the following system of a simple prey -specialist predatorgeneralist predator interaction [12] where a 1 , a 2 , b 1 , w, w 1 , w 2 , w 3 , D, D 1 , D 2 , c and θ are positive constants.The detailed description of the model system is given in the paper by Upadhyay et al. [12].Since the generalist predator z in (1c) are assumed to be sexually reproducing species, their growth has two phases: a linear phase and a quadratic phase [22].So in this case, the last equation (1c) is modified to where D 3 represents the residual loss in z population due to severe scarcity of its favorite food y.The typical situation represented by this model is presented in Fig. 1.Equations (1a), (1b) and (1d) represent model system (1).We choose to study the nonlinear phase (model system (1)) as the linear phase does not support the chaotic behaviour at all.The sexually reproducing population are covered by this phase when they are under Allee effect [23].Hastings and Powell [9] produced a new example of chaotic population system in a simple tri-trophic food chain model with Holling type II functional responses.Chattopadhyay and Sarkar [11] modified the Hastings and Powell model by introducing an extra mortality term in zooplankton population and studied the reduction of propensity of chaotic dynamics as described by Hastings and Powell.Both Hastings and Powell [9] and Chattopadhyay and Sarkar [11] have used half saturation constant of zooplankton as key parameter in their model to study the system dynamics from order to chaos.Mandal et al. [24] studied the modified Hastings and Powell model by considering different body sizes of zooplankton and accordingly the growth rate and half saturation constant have changed.This interaction is represented by the following system of a simple preyspecialist predator -specialist predator interaction.The basic mathematical model can now be represented by a set of three ordinary differential equations describing the rate of change of TPP, zooplankton and fish population over time [11,24] given by where w 2 measures the maximum value attainable by the per-capita functional response of the specialist predator z which feeds only on y.The parameter c is the decay rate of the predator z in absence of its prey y and w 3 is a measure of its assimilation efficiency.In case of θ = 0, the model has been studied by many researchers [5,7,9].These equations (2a)-(2c) describe model system (2).The real world example for this model is presented in Fig. 2.
To characterize interface between phytoplankton and zooplankton populations in the presence of toxic chemical, Holling type I, II, and III functional responses for f (x) have been considered to study the behavior of the system [13,16,25].

Stability analysis and Hopf bifurcation
The model system (1) has the following equilibrium points: • The trivial equilibrium point E 0 = (0, 0, 0) always exists.
• The equilibrium point E 1 = (a 1 /b 1 , 0, 0) exists on the boundary of the first octant.
• E 2 = (x, ȳ, 0) is the planer equilibrium point on x − y plane, where x is the positive root of the equation Clearly E 2 exists provided the following condition satisfied • The nontrivial equilibrium E 3 = (x * , y * , z * ) exists if and only if there is a positive solution to the following set of equations: Straight forward computations show that x * is a positive root of the quadratic equation Therefore, E 3 = (x * , y * , z * ) is positive equilibrium point under the following conditions: Now, in order to investigate the local behavior of model system (1) around each of the equilibrium points, the variational matrix V of the point (x, y, z) is computed as Let V i , i = 0, 1, 2, 3 denotes the variational matrix at E i , i = 0, 1, 2, 3 respectively.Hence From the above variational matrix, it is observed that there is an unstable manifold along x-direction and a stable manifold along y-direction.Therefore, the equilibrium point E 0 is a saddle point.The variational matrix for E 1 is , From the variational matrix V 1 , it is found that the equilibrium point E 1 is locally asymptotically stable provided w1a1 a1+b1D1 < a 2 + θf ( a1 b1 ).The variational matrix about another equilibrium point E 2 is .
Therefore, the linearized systems about E 0 , E 1 and E 2 have zero eigen values.Therefore, these are non-hyperbolic points and hence the dynamical behavior near them can be stable, periodic, or even chaotic.
However, for the positive point E 3 = (x * , y * , z * ) the variational matrix is .
According to Routh-Hurwitz criterion, E 3 = (x * , y * , z * ) is locally asymptotically stable provided the following conditions are satisfied A 1 > 0, A 3 > 0, and A 1 A 2 > A 3 ; where A i , i = 1, 2, 3 are the coefficients of the characteristic equation of V 3 = ⌊a ij ⌋, i, j = 1, 2, 3: Straight forward computations show that, A 1 > 0 and A 3 > 0 if and only if the following condition is satisfied: Also since, Therefore, the necessary condition for Further more, by substituting the values of a ij , i, j = 1, 2, 3, we get Clearly, M 1 > 0 provided that condition ( 6) is satisfied and M 2 > 0 provided that condition ( 5) is satisfied.Hence, if conditions ( 5) and ( 6) hold, then the necessary and sufficient condition for Finally, due to the above analysis, the local stability analysis of the positive equilibrium point E 3 = (x * , y * , z * ) can be summarized as the following theorem.
Now, in order to investigate the Hopf bifurcation of model system (1), we will follow the Liu approach [26].According to Liu approach, the simple Hopf bifurcation at µ = µ * can occur provided that

2.
dΨ(µ) dµ µ=µ * = 0. Now, let c, the growth rate of the generalist predator, be the bifurcation parameter.Therefore, if conditions ( 5) and ( 6) hold together with the following condition Then, obviously Further it is easy to verify that Accordingly, the following theorem establishes the Hopf bifurcation conditions: Theorem 2. Under the conditions (5), (6), and (8), there is a simple Hopf bifurcation of the positive equilibrium point E 3 = (x * , y * , z * ) of model system (1) at some critical value of the parameter c given by (8).

Numerical results
Our primary interest is to explore the occurrence of chaotic dynamics in two model systems which differ from each other in one essential way i.e., the top predator in model system (1) is a generalist predator whereas that of model system (2) is a specialist predator.We also try to observe the role of toxin producing phytoplankton on the chaotic dynamics in such ecosystems.The role of TPP for controlling blooms or for decreasing grazing pressure is now well known but the functional forms for releasing toxic substances are not known [16].For this reason, we have considered Holling types I, II, and III functional forms to describe the liberation of toxin production process, motivated by the literatures available in this field [13,14,16,25,27].For the following form of the functional responses: (Holling type II) (Holling type III).
Model system (1) and model system (2) are integrated numerically using six-order Runge-Kutta method along with predictor corrector method.It is observed that model system (1) has a chaotic solution at the following set of parameter values (see Fig. 3).These parameter values are selected on the basis of paper by Letellier and Aziz-Alaoui [28].However, Hastings and Powell [9] and Rai and Upadhyay [7] observed that, for the following set of parameter values, model system (2) exhibits a chaotic dynamics (see Fig. 4).Our approach is to observe the exchanging of the states (chaos, period doubling, limit cycle, stable point) for different values of θ and also for different types of functional response forms.In real life situation, it has been observed that increasing the strength of toxic substances has a stabilizing effect.So, in this paper we like to see whether this is true or not in our considered model systems and then discuss the effect of different types of functional forms on the liberation of toxin production process.
There exists empirical evidence that the level of nutrients in an ecological system is coupled with toxin release by phytoplankton [29].Thus, it is important to examine the dynamical complexity of the model system.The common tool used for the study of dynamical complexity of two considered model systems is bifurcation diagrams with θ as the bifurcation parameter.For different types of functional response forms, we have plotted the successive maxima of top predator z as a function of the parameter θ (rate of toxin substances released by TPP population) keeping other parameters fixed as given in equation (9) for model system (1) and equation (10) for model system (2).
The figures, Fig. 5, Fig. 6, and Fig. 7, are representing the bifurcation diagrams of model system (1) with functional response of Holling type I, II, and III respectively.However, for model system (2), the bifurcation diagrams are also drawn for functional response given by Holling type I, II, and III, and are presented in Fig. 8, Fig. 9, and Fig. 10 respectively.All these figures show clearly the transition from chaos to order through sequence of period halving bifurcation.Therefore, for both the models and for different forms of toxic substance liberation process (i.e.f (x)), it is observed that, increase of value of toxic substances released by TPP population (i.e.θ) has a stabilizing effect.The blow-up bifurcation diagrams show that the model system possesses rich variety of dynamical behaviour for bifurcation parameter θ in the ranges [0, 0.0029] for Holling type I (see Fig.  Further, Table 1 and Table 2 describe the dynamical behavior of model system (1) and model system (2) respectively according to the above mentioned bifurcation diagrams.From Table 1, it is found that chaos was observed in the ranges Table 1.Dynamical behavior (DB) of model system (1) depending on the results of bifurcation diagrams given in Fig. 5, Fig. 6 and Fig. 7. Pi -limit cycle of period i (i = 2, 3, 4, 5, 6), SF -stable focus, LC -limit cycle, LP -long period, SCA -strange chaotic attractor, EX -extinction Results of model ( 1) for Results of model ( 1) for Results of model ( 1) for Holling type I: Holling type II: Holling type III: From Tables 1 and 2, it has been also observed that the gradual increase of toxin values of phytoplankton in the model systems which turn the system dynamics from chaos to doubling state to different order limit cycles and the system finally settles down to an equilibrium state.If we increase more toxin to the model systems then the dynamics shows the extinction of zooplankton.It will be a subject of further study to examine if chaotic dynamics is robust against changes in other system parameters in model ecological systems which include the effect of toxin production by phytoplankton.Dynamical complexities of these model systems as manifested in bifurcation diagrams (Figs.5-10) possess many subtleties.Because of the coexistence of non-chaotic attractors, i.e., two, three, four or six-periodic attractors, in ecological interaction, the presence of chaotic attractor cannot be judged solely on the basis of the sensitivity dependence of initial conditions.One can clearly see that chaotic dynamics is robust to changes against the rates in toxin production by phytoplankton as it exists in large range of θ values.Table 2. Dynamical behavior of model system (2) depending on the results of bifurcation diagrams given in Fig. 8, Fig. 9 and Fig. 10.Pi -limit cycle of period i (i = 2, 3, 4, 5, 6), SF -stable focus, LC -limit cycle, LP -long period, SCA -strange chaotic attractor, EX -extinction Results of model ( 2) for Results of model ( 2) for Results of model ( 2) for Holling type I: Holling type II: Holling type III:

Conclusions
According to the above discussion for the dynamical behavior of model system (1) and model system (2) (see the bifurcation diagrams and Tables 1 and 2), it is clear that increasing the strength of toxic chemical released by TPP population reduce the prevalence of chaos, therefore, our observation is that the toxin released by TPP population acts as a bio-control by changing the state of chaos to order.Further, From Table 1 and Table 2, extinction has been observed for higher values of θ in Holling type II and III functional responses.However, in case of Holling type I, extinction was observed very earlier at θ = 0.017 for model system (1) and θ = 0.016 for model system (2).Therefore, the rate of toxic substance released by TPP is to be high for type I functional form than those of type II and type III functional form.These observations indicate us that to maintain the order of an ecosystem functioning, type II or type III functional form for toxin liberation process is more appropriate.We have also observed that for different ranges of θ, the rate of toxin release by TPP, the systems locally stable as well as oscillate around the equilibria related to planktonic blooming.The phenomenon of dynamic stabilization of a locally unstable equilibrium point can also be observed for certain rate of toxin production.Different Holling type mortality function of zooplankton due to TPP has significant influence in controlling oscillations, coexistence, survival or extinction of zooplankton.The rate of toxin production by TPP plays an important role for controlling oscillation in plankton systems.
The consequences of toxin production by phytoplankton are of great interest both for intrinsic scientific merit and also because of its possible detrimental effect on fisheries.In our earlier work [7], we had suggested that the biology of the top predator would be a crucial factor for the determination of dynamical complexity in food chain models.The dynamical complexity depends on the nature (e.g., inter-specific competition vs predation) and degree (e.g., the number of interacting components) of complexity of trophic interactions.On the basis of present study, we opine that the natural systems with first kind of food chains (i.e., chain ending at specialist predator as in model system (1)) would present difficult challenges as far as program of quantification of their dynamical complexity is concerned.The other kind of systems (the chain ending at generalist predator as in model system (2)) seems to allow such program to be implemented smoothly.This conjecture is to be tested in the laboratory and in the field.