Chaotic Dynamics in a Three Species Aquatic Population Model with Holling Type II Functional Response

Abstract. A three-trophic model for marine community is proposed and i nvestigated by means of numerical bifurcation analysis. The proposed mo el based on a modified version of the Leslie-Gower scheme, incorporates mutual in terference in all the three populations and generalizes several other known models in t he ecological literature. We investigate the dynamical behavior of the model system by co nsidering the Holling type II functional response of toxin liberation process. Bifurcat ion diagram and two-dimensional parameter scan suggest that chaotic dynamics is robust to va riations in toxin production by phytoplankton. Our study suggests that toxic substances rel ased by TPP population may act as bio-control by changing the state of chaos to order . The mutual interference also induces chaos and acts as both stabilizing and destabil izing factors.


Introduction
Ecological systems have all the necessary characteristics (nonlinearity, high-dimensions, etc.) to support chaotic dynamics [1].Chaotic dynamics and limit cycles are common in tri-trophic food chain model and are of common interest to both the theoretical and experimental population biologists/ecologists.To assess the ecological implications of chaotic dynamics in different natural system, it is important to explore changes in the dynamics when structural assumptions of the system are varied.One such approach to the study of the dynamics of marine ecological community is its food web and the coupling of interacting species with each other [2].Upadhyay and Rai [3] provided new examples of a chaotic population system in a simple tri-trophic food chain with Holling type II functional response.Aziz-Alaoui [4] revisited the Upadhyay and Rai model and found that the chaotic dynamics is observed via sequences of period-doubling bifurcation of limit cycles which however suddenly break down and reverse giving rise to a sequence of period-halving bifurcation leading to limit cycles.Upadhyay and Chattopadhyay [5] modified the model of Upadhyay and Rai [3], by introducing an extra mortality term in middle predator and interpreted the system for aquatic environment consisting of TPP-Zooplankton-Molluscs food chain model.They observed that an increase in the strength of toxic substance released by toxin producing phytoplankton population reduces the propensity of chaotic dynamics and changes the state of chaos to limit cycle and finally settles down to stable focus.Ruxton [6] also showed that the system of linked populations has a stabilizing effect on tri-trophic food chain model.Further study [7] reveals that the rate of toxin production by TPP plays an important role for controlling oscillations in the plankton system.
Many studies investigated the effect of mutual interference on the population dynamics.DeAngelis et al. [8] studied the dynamical properties of a continuous-time autonomous model system incorporating their interference model.This model was studied by Hwang [9] to establish that the periodic orbits, if they exist, are unique.The models considered for interference have different mathematical expressions and different conceptual foundations [10].From theoretical studies and empirical evidences, a consensus has been reached to conclude that interference has a stabilizing influence on population dynamics [11], although Hassell and May [12] pointed out that there was an upper limit on the interference constant beyond which the dynamics becomes unstable.Predator-prey models incorporating mutual interference were first proposed in Hassell [13] and Rogers and Hassell [14].A model incorporating density-dependent death rates was considered by Levin [15].Freedman and Rao [16] considered the Gause-type model incorporating mutual interference among predators and a density-dependent predator death rate.Erbe and Freedman [17] applied it to the simple food chain, modeled by Lotka -Volterra dynamics.Here, I have used the concept for modelling it with modified version of Leslie-Gower scheme in a simple food-chain model modelling marine ecosystems.Motivated by the above studies, we show that the chaotic behaviour as described by Upadhyay and Rai [3,18] could be controlled by an auto-control mechanism.
In this paper, we propose a generalized model of aquatic ecological system by introducing mutual interference in all the three populations, an extra mortality term in zooplankton population and also taking into account the toxin liberation process of TPP population.This model generalizes several other known models in the literature like Upadhyay and Rai model [3,18] and Hastings and Powell model [19].One of the main objectives of this study is to examine the roles of mutual interference parameters and the parameter θ, the rate of toxin release by TPP population on the dynamics of the model system.Different types of toxin release function f 1 (x 1 ), which represents the toxin liberation process of TPP population is considered.The results reported in this paper are only for Holling type II functional response.
This paper is organized as follows: in Section 2, we present the details of the model system.The methodology used is presented in Section 3, to help us in selecting the parameter values to perform simulation experiments.Numerical results are summarized in Section 4 and conclusions are presented in Section 5.

The mathematical model
Consider a situation where a prey population x 1 is predated by individuals of a population x 2 .The population x 2 , in turn serves as a favourite food for individuals of a population x 3 .This interaction is represented by the following system (prey -specialist predatorgeneralist predator interaction) of ordinary differential equations: where m i > 0 for i = 1, 2, 3, a 1 , a 2 , b 1 , θ, w 0 , w 1 , w 2 , w 3 , c and D 0 , D 1 , D 2 , D 3 , D 4 are the positive constants and f i ∈ C 0 (R + ) for i = 1, 2. The parameters m i for i = 1, 2, 3 are mutual interference parameters that model the intraspecific competition among predators when hunting for prey [13,16,17,20,21].
In this model, TPP population (prey) of size x 1 serves as the only food for the specialist predator (zooplankton) population of size x 2 .This zooplankton population, in turn, serves as a favorite food for the generalist predator (molluscs) population of size x 3 .The equations for rate of change of population size for prey and specialist predator have been written following the Volterra scheme that is, predator population dies out exponentially in the absence of its lone prey.The interaction between this predator x 2 and the generalist predator x 3 is modeled by the modified version of the Leslie-Gower scheme, where the loss in a predator population is proportional to the reciprocal of per capita availability of its most favorite food.a 1 is the intrinsic growth rate of the prey population x 1 , a 2 is the intrinsic death rate of the predator population x 2 in the absence of the only food x 1 , c measures the rate of self-reproduction of generalist predator x 3 , w 0 , w 1 , w 2 , w 3 are the maximum values which per capita growth rate can attain.b 1 measures the strength of intra-specific competition among the individuals of the prey species x 1 .D 0 and D 1 quantify the extent to which environment provides protection to the prey x 1 and may be thought of as a refuge or a measure of the effectiveness of the prey in evading a predator's attack.D 2 is the value of x 2 at which per capita removal rate of x 2 becomes w 2 /2.For m 1 = 1 the coefficient w 0 /(x 1 + D 0 ), of the third term on the right hand side of (1a) is obtained by considering the probable effect of the density of the prey's population on predators attack rate.If this coefficient is multiplied by x 1 (the prey population at any instant of time), it gives the attack rate on the prey per predator.Denote p(x 1 ) = w 0 x 1 /(x 1 + D 0 ), when x 1 → ∞, p(x 1 ) → w 0 which is the maximum value that it can reach.Some aquatic organisms condition their medium by releasing substances that stimulate growth of species, which have similar genetic make-up.Sparse populations rarely provide sufficient opportunities for social interaction necessary for reproduction.
Here, f 1 (x 1 ) represents the toxin liberation process of TPP population for which the mortality of zooplankton increases and as a result, the grazing pressure of zooplankton on TPP population decreases.The parameter θ is the rate of toxin release by TPP population.w 3 measures the limitation on the growth of the generalist predator x 3 by its dependence on per capita availability of its most favorite prey x 2 represented by the function where D 3 represents the residual loss in x 3 population due to severe scarcity of its favourite food x 2 .Equations (1a)-(1c) describe the proposed model system.
A model system could be more realistic from ecological point of view and interesting from mathematical point of view if one considers different predator's functional response and compares the dynamic effects of these functional responses.Since functional response encapsulates attributes of both prey and predator biology, so handling time, search efficiency, encounter rate, prey escape ability, etc. should alter, in general, the functional responses [22].Therefore, predator's functional response may be different when a particular predator preys different prey having different escape ability and if a particular prey is predated by different predators having different hunting ability.The structure of prey habitat is also responsible to alter the functional response.Thus, a predator which follows type II functional response in homogeneous habitat may follow type III in a heterogeneous medium.Anderson [23] experimentally observed in a kelp bass-kelp parch predatorprey interaction for none and medium amounts of habitat structure, the type II functional response had a better fit than linear models.However, for the highest amount of habitat structure a type III functional response had a better fit.In reality, the raptorial behaviour of copepods is highly complex and exhibits a hunting behaviour [24] and hence type II or type III is an appropriate choice.To characterize interface between phytoplankton and zooplankton populations in the presence of toxic chemical, Holling type II functional responses for f 1 (x 1 ) is considered to study the dynamical behaviour of the model system.
A separate investigation is required when the parameters m i are sub-linear (0 < m i < 1).In this case, we make the following assumptions: Assumption 1.There exist functions h j continuous on R 3 + , where + .As in Erbe et al. [17], we consider the following change of variables for (1a)-(1c) (2) The system (1a)-(1c) transforms as The change of variables given in (2) transforms the sublinear system (3c) into (3b) in which no sublinearities are present.Biologically, this amounts to requiring that the mutual interferences are not too strong.
The above discussions may be summarized as follows.

Methods of investigation
The model system presented above is a multi-parameter system.Model parameters are selected in accordance with a method given in upadhyay et al. [3,18].A few hundred parameter combinations (choosing two at a time) are possible.This is simply not feasible for any one to scan the system in all the parameter spaces.Application of non-linear dynamics is in unison with the knowledge of biology of the system and enables one to choose parameter combinations for simulation experiments.The most crucial part of the present methodology is the following conjecture: Two coupled Kolmogorov systems in oscillatory mode would yield either cyclic (stable limit cycles and quasi-periodic) or chaotic solutions depending on the strength of coupling between the two.
In the present case, the set of parameter values for which the system admits a limit cycle solution is found to be There is one more important aspect of these simulation experiments i.e., choosing the step size for the variation of a system parameter from a parameter combination within the chosen range.It depends on the nature of the parameter concerned: whether it is a slow varying or fast varying one.
The most useful way to study such a dynamical system is to monitor the amplitude (maxima) of the subsequent oscillations as the control parameter of the system is varied.A small change in parameter values may lead to a bifurcation: an abrupt, qualitative change in the dynamics.There are number of ways to detect chaotic dynamics in dynamical systems.We have used in our study the phase space representation, bifurcation diagram and two dimensional scan.

Numerical results
Model system is integrated numerically using six-order Runge-Kutta method along with predictor corrector method.It is observed that the model system (1a)-(1c) has a chaotic solution at the following set of parameter values (see Fig. 1).The parameter values are selected on the basis of previous studies [4,7] and correspond to quantitative measures of attributes of the TPP-Zooplankton-Molluscs food chain.To confirm the existence of chaos, the dynamics of the model system is studied by constructing bifurcation diagram.For Holling type II functional response form for toxin liberation process, we have plotted the successive maxima of top predator x 3 as a function of the parameter θ (rate of toxin substances release by TPP population) keeping other parameters fixed as given in equation ( 5) for model system (1a)-(1c).The Fig. 2 represents the bifurcation diagrams of model system (1a)-(1c) with Holling type II functional response.This figure shows the transition from chaos to order through sequences of period-halving bifurcation.From this bifurcation diagram, it is observed that an increase in the value of toxic substances released by TPP population has a stabilizing effect.The blow-up bifurcation diagram (see Fig. 2(b)) shows that the model system possesses rich variety of dynamical behaviour for bifurcation parameter θ in the range [0, 0.06].A period -doubling cascade is observed.After the accumulation point, the behaviour settles down onto a chaotic attractor.When θ, the bifurcation parameter is decreased, new periodic orbits are created.The chaotic attractor emanating from the main one is destroyed by a boundary crisis with the unstable periodic orbit created by the saddle-node bifurcation.A saddle-node bifurcation is merge and disappearance of two steady states one of them is saddle and other is node.Two co-existing period -doubling cascades are then observed.Dynamical behavior of model system (1a)-(1c) depending on the results of bifurcation diagrams given in Fig. 2 is presented in Table 1.From this result, we observe stable focus, different order limit cycles and strange chaotic attractor in different ranges of θ, the rate of toxic substance released by TPP.Also, we conclude that for the model system, the increase in the value of toxic substances released by TPP has a stabilizing effect.These observations indicate that to maintain the order of an ecosystem functioning, Holling type II functional form for toxin liberation process is more appropriate.Table 1.Dynamical behavior of model system (1a)-(1c) depending on the results of bifurcation diagrams given in Fig. 2. Pi -limit cycle of period i for (i = 2, 4, 5, 6), SF -stable focus, LC -limit cycle, LP -long period, SCA -strange chaotic attractor Results of model (1a)-(1c) for Holling type II functional response We have also investigated the role of mutual interference parameters on the dynamics of trophic system in detail.The values of mutual interference parameters were chosen on the basis of the values reported in Katz [26].We have observed stable focus, limit cycles and chaotic dynamics phenomena in the model system by changing the mutual interference parameters m i , i = 1, 2, 3 and the rate of toxin release by TPP population θ, in the fixed range.We have also reported the function error or argument domain error, the region in the parameter space, where no dynamics is observed.In this domain, the values of mutual interference parameters are not conducive for simulation experiment i.e., in real situation, no species can attain these values of mutual interference.Our approach is to fix m 1 and m 2 then vary m 3 in the interval [1,3] and θ in the interval [0, 1) and then observe the exchange of states (stability -limit cycle -period doubling -chaos) in the model system for three different cases of m i (m i >, =, < 1).
It is found that in most of the cases, x 2 becomes extinct and (x 1 , x 3 ) rests on stable focus for higher values of θ.For lower values of θ, all the three populations rest on stable focus and limit cycle attractor in the phase plane.It is also observed that for m 1 = m 2 = 1.25, 1.5, 1.75 and for whole range of the parameter space (m 3 , θ) (i.e., 1 ≤ m 3 ≤ 3, 0 ≤ θ ≤ 1), the model system (1a)-(1c) predicts no dynamics.The simulation results show function error or argument domain error.

Conclusions
In this article, we have attempted to find whether mutual interference and toxic substances released by TPP always stabilize the predator -prey dynamics in aquatic environment?Our simulation experiments support the conclusion that TPP stabilizes the predator -prey dynamics in aquatic environment.From the tables, it is observed that for different values of mutual interference parameters in different ranges, dynamics of the model system is also influenced by the functional form of toxin liberation process.For m i < 1, i = 1, 2, 3, no dynamics was observed in the range 0.25 ≤ m i ≤ 0.75, but if we take the values of any one of the interference parameters close to 1, the system dynamics converges to stable focus.In this case, the top predator becomes extinct as m 3 reaches 1.For m i > 1, most often the dynamics rests on stable limit cycle or stable focus.From Tables 2, it is found that for m 1 = m 2 = 1.05 (i.e., close to 1) and m 3 in the range [1,3], system dynamics settles down to limit cycle attractor.In this case, model system also supports chaotic dynamics only at a few discrete points.But for m 1 = m 2 = 2 and m 3 in the range [1,3], the system dynamics mostly settled on stable focus and middle predator becomes extinct.These results show that the interaction between predators is a stabilizing factor.Chaotic dynamics/situation may arise from an equilibrium state for different reasons in any ecosystem.But to overcome this chaotic situation sometimes system itself has some mechanism and self-adaptability.There are many ways by which system can be selfadjusted and one of such ways is toxin production by phytoplankton, which reduces the zooplankton grazing, helps the system to recover chaotic situation.In aquatic system of such condition it is reported in Mandal et al. [2], that toxins are produces by many phytoplankton and these toxins may turn the ecosystem into ordered state from chaos by reducing the grazing pressure of zooplankton.
From the tables and 2D scan diagram, it was also observed that the model system supports chaotic dynamics for m i ≥ 1.We also observe from bifurcation diagram that chaotic dynamics is robust to changes in changes against rates in toxin production by phytoplankton as it exists for large range of θ value.Period doubling bifurcations seem to be responsible for this kind of dynamical behaviour.
In real life situations, it has been observed that increasing the strength of toxic substances and mutual interference parameters has a greater stabilizing effect.Here, we like to see whether this is true or not in our considered model system.Our simulation results show that interference might actually strongly destabilize the dynamics as well leading to chaotic dynamic behaviour.Further studies are needed to ascertain if this defense mechanism suppresses chaotic dynamics in model aquatic systems.