Dynamical complexities in a tri-trophic hybrid food chain model with Holling type II and Crowley – Martin functional responses

Abstract. We study how predator behavior influences community dynamic s of predatorprey systems. It turns out that predator behavior plays a dom inant role in community dynamics. The hybrid model studied in this paper reveals tha t period-doubling and period-doubling reversals can generate short-term rec urrent chaos (STRC), which mimics chaotic dynamics observed in natural populations. S TRC manifests itself when deterministic changes in a system parameter interrupt chao ti behavior at unpredictable intervals. Numerical results reinforce an earlier suggest ion that period-doubling reversals could control chaotic dynamics in ecological models. In eco logical terms, the prey and intermediate predator populations may go to extinction in t he event of a catastrophe. The top predator is always a survivor. In contrast to this, th is is not the case when the constituent populations are interacting through Holling t ype II functional response. Even this top predator can go to extinction in the event of such cat astrophes.


Introduction
Understanding ecosystem's dynamics is one of the most challenging tasks.The biotic part of the ecosystem comprises living entities known as species.These species are spread over space.The key to understanding of ecosystem dynamics is population systems.The populations are groups of individuals of different species spread over a given geographical area.When population densities are high enough, one may regard the interaction of food chain species as "well-mixed"; at such densities, laws of chemical kinetics govern population dynamics.In what follows, we discuss how such population dynamic models are designed.
The population dynamic models are assembled by combining functional and numerical responses of predators in a suitably chosen scheme.Two such schemes (formulations of predator-prey interaction) are in vogue; (i) Volterra scheme (ii) Leslie scheme.The former assumes that the predator population dies out exponentially in the absence of its lone food; latter is based on the premise that the per capita growth of the predator population is limited by the per capita availability of prey.We present a model, which is designed using Volterra scheme.The model system employs Holling type II functional response for the first predator.The top predator's functional response is modeled by Crowley-Martin (CM) function.
The first predator-dependent model was given by Hassell and Varley [1], who proposed that the attack rate should decrease with increasing predator density.They proposed the functional response in the form Parameter m can be interpreted as an interference coefficient.This view is also plausible biologically [2].When m = 0 or Y = 1, the Hassell-Varley functional response reduces to Holling type II functional response.Holling type II response function assumes that the predator spends some time searching for prey and some time for processing each captured prey item.The instantaneous per capita feeding rate of the predator is given by where N is the prey density.This response function assumes that there exists no interference among individuals of predators and that depletion of prey causes competition among predators for food [3].Beddington-DeAngelis (BD) functional response assumes that individual predator allocates time not only to searching for prey, but also to engaging in encounters with other predators.The instantaneous per capita feeding rate is given by , where P is the predator's density and c is a positive constant describing the magnitude of interference among predators.The underlying assumption of BD model is that handling and interfering are exclusive activities.The CM model allows for interference among individuals of predators engaged in handling or searching at a given instant of time.The CM model adds an additional term in the denominator , (P − 1) is replaced by P when predator's interference is modeled as continuous variable or other mechanisms of predator dependence such as prey behavior that depends upon predator's density.The basic difference between these two functional responses is that the BD model predicts that the effects of predator interference become negligible at high prey densities.While the CM model maintains that the interference effects on feeding rate remain important at all densities.
Predators behave differently in different situations they are foraging.In spacelimited situations where predators can move freely but can only encounter prey in a region of limited size, predator interference effects are produced by spatial mechanisms [2].Common examples of space-limited situations arise when predators could only encounter prey at the edge of a prey refuge or if the prey were only exposed to predation in gaps in protective cover.It is expected that natural food chains could be better represented by models, which incorporate both Holling type II and CM functional responses.The latter assumes that the per capita feeding rate of the predators is a function of predator interference effects.In this paper, we explore dynamical complexities (attractors and their basins, bifurcations, etc.) of such a hybrid food chain model.Upadhyay and Naji [4] studied the local and global stability of this hybrid model system and also established the persistence criteria.
Hybrid ecological models do not have pristine history.Upadhyay and Rai [5] proposed and studied a hybrid food chain model to understand why deterministic chaos is rarely observed in natural populations.These models were designed by combining two formulations of predator-prey dynamics: Volterra and Leslie-Gower schemes.The former describes the population dynamics of a specialist predator and the latter that of a generalist predator.The authors have been successful in their attempt to develop a theory of ecological chaos based on these new classes of models [6,7].The present study is based on numerical computations.Stability analyses can be performed for simple systems of differential equations describing trophic structures.It is impossible to determine analytically the nature of unstable dynamics (regular vs chaotic oscillations).Recently, Gross et al. [8] proposed a method to investigate the potential for chaotic dynamics in general food chain models of variable length.It would be beneficial if this approach, based on bifurcations of higher co-dimension as indicators of chaos, could be extended to food web architectures.

Model system
We consider the following system as a model representing a tritrophic food chain.The model employs both Holling type II and CM type of functional responses.It is described by following system of differential equation, where X(T ) is the population density of the lowest trophic level species (prey) at time T .Y (T ) is the population density of the middle trophic level species (intermediate predator) at time T and Z(T ) is the population density of highest trophic level species (top predator) at time T .The intermediate predator Y feeds on the prey X according to the Holling type II functional response, however the top predator (Z) feeding rate on Y varies according to the CM type functional response [4].
The prey X grows with intrinsic per capita growth rate a 1 and carrying capacity K in the absence of predation.D measures the extent to which the environment provides protection to prey X. w is the maximum value which per capita reduction rate of X can attain, w 1 /w represents the conversion rate of an eaten prey.The constants w 2 , w 3 , b, and d are the saturating CM type functional response parameters, in which b measures the magnitude of interference among predator.Further, a 2 is the death rate of the intermediate predator and a 3 is the death rate of the top predator.All the model parameters are assuming only positive values.Obviously, model system (1) is a three species simple food chain involving a hybrid type of prey-dependent and a predatordependent functional response.
The model system is rendered dimensionless using the following variables and parameters: The model equations in dimensionless form are: Clearly, the non-dimensional system (3) has eight parameters in all.Obviously the right hand sides of system (3) are continuous smooth functions on R 3 + = {(x, y, z) ∈ R 3 : x ≥ 0, y ≥ 0, z ≥ 0}.Indeed, they are Lipschizian on R 3 + and then the solution of the system (3) exists and is unique.

Methodology
Computer simulations were performed on MATLAB 7.0.Model system (3a)-(3c) is integrated numerically using the fourth-order Runge-Kutta method with time step size 0.001.We explore dynamical complexities of the proposed model system in two-dimensional parameter spaces.Bifurcation diagrams are computed by treating crucial parameters of the model system as bifurcation parameters.Non-linear dynamical systems exhibit dynamical complexities, which are either of deterministic or of stochastic origin.2D scans help us to examine the former and the latter which are investigated by computing basin boundary diagrams.We also generate phase portraits and time series of the model system for certain parameter sets.
We now mention briefly how basin boundary calculations are performed.First, we define the basin of infinity.Let SD denote the diameter of the computer screen.It may be possible that the point at infinity is an attractor.Since we cannot examine rigorously whether the trajectory of a point goes to infinity, we conclude that a trajectory diverges or is diverging if it leaves the computer screen area, that is, it goes to left or to right of the screen by more than one SD width of the screen, or goes above or below the screen area by more than one SD screen height.The basin of infinity is the set of initial points whose trajectories are diverging.Maryland Chaos group has done pioneering work in this area and have developed a tool to calculate basin boundary structures.We have used the research version of the software which accompanied the book entitled "Dynamics: Numerical Explorations" authored by Nusse and Yorke [9].We have used the BAS (basins and attractors structure) method for all the computations.This method divides the basin into the following two groups: (i) The basin of attraction A whose points will be plotted, (ii) The basin B whose points will not be plotted.A generalized attractor is the union of finitely many attractors, and a generalized basin is the basin of a generalized attractor.The BAS routine does not plot the bowl lying outside.The strategy is to test each grid point which is the centre of the grid box.In the event that the centre of a grid box is in basin A, then the same is plotted (colored).In the default case, basin A is the set of points whose trajectories are diverging, while basin B is empty.Therefore, BAS routine will plot a grid box if the trajectory of its centre is diverging.The important aspects of the basin boundary calculations are to specify the basins A, B and to find the radius RA, where RA stands for radius of attraction for storage vectors which help to specify the basins A and B. The value of RA will be different for different dynamical systems.It must be set appropriately in order to avoid any misleading basin picture.

Numerical simulations
In order to understand the dynamics of the model system (3), we turn to numerical simulations.It is observed that the model system (3a)-(3c) has a chaotic solution for the following set of parameter values: The parameter values are chosen on the basis of our previous studies [4].Time series corresponding to the chaotic attractor after transients were died out were recorded for a typical set of parameter values and are presented in Fig. 1.The chaotic attractor for the model system (3a)-(3c) is presented in Fig. 1(a).The time series representations for chaotic dynamics in the model system are presented in Figs.1(b)-1(d) which shows that while populations x and y reach extinction-sized densities, the population density of the top-predator, z does not.Extensive numerical simulations are carried out for various values of parameters and for different sets of initial conditions.Two different control parameters are discussed, the death rate of intermediate predator scaled to per-capita reproduction rate of its prey, w 5 and scaled death rate of top predator w 11 .These values are selected after a thorough study of the asymptotic dynamics of Kolmogorov subsystem, in particular the system parameters are chosen in such a way that subsystems perform their journies on stable limit cycles.The bifurcation diagram of system (3) for the successive maxima of the top predator population z as a function of w 11 in the range 0 < w 11 < 0.1 with step size 0.01, is plotted in Fig. 2. According to this bifurcation diagram the solution of system (3) has different type of attracting sets in the range 0 < w 11 < 0.055 including periodic and chaotic ones.However, as w 11 increases further 0.055 < w 11 < 0.1 a stable limit cycle is observed.Moreover, for w 11 > 0.1, the top predator population z declines and reaches to extinction.It is evident from Fig. 2 that PD and PDR bifurcations take place when parameter w 11 is decreased.Deterministic chaos shows up in (0.01, 0.02).The magnified bifurcation diagrams of Fig. 2 are blown up in Fig. 3

(a) and 3(b).
A typical long-term attracting set of system (3) is drawn in Figs.4(a)-4(d).The blown up bifurcation diagram Fig. 3, and the attracting sets given in Fig. 4 show clearly the route to chaos through the cascade of periodic-doubling bifurcations.In fact for the range of 0 < w 11 < 0.055, the solution of system (3) has rich dynamics generated by perioddoubling and period-halving bifurcations.For the same parameters given by equation ( 4), if we fix w 5 = 0.25, the strange chaotic attractor given in Fig. 4(d), revealed via a perioddoubling cascade the diagram of which is given by Fig. 3(a).This phenomenon of perioddoubling can suddenly break down and reverse, giving rise to period-halving bifurcations (cf.Fig. 3(b)) leading to stable limit cycles.It is observed that for w 11 = 0.01, 0.022, 0.025 and 0.035 the solution of system (3) approaches to chaotic attractor, period-2 attractor, period-4 attractor and chaotic attractor respectively.
Two bifurcation diagrams of system (3) for the successive maxima of the intermediate predator population y as a function of scaled death rate of intermediate predator w 5 in the range 0.15 < w 5 < 0.5, keeping other parameters as given in equation ( 4) with w 11 = 0.03, are drawn in Fig. 5(a).The bifurcation diagram given in Fig. 5(a) shows that the solution of system (3) is very sensitive to the death rate of intermediate predator in the range 0.18 < w 5 < 0.44.This figure shows clearly the presence of the cascade of periodic doubling bifurcations leading to chaos.
Similar bifurcation diagrams are drawn in Fig. 5(b), as those in Fig. 5(a), for the parameter values given in equation ( 4) with 0.25 < w 5 < 0.5 and w 11 = 0.06.Although the food chain system has rich dynamics including chaos, an increase in the death rate w 11 from 0.03 to 0.06, causes the reduction in the zones of chaos and extension in the periodic windows in the range 0.28 < w 5 < 0.44.Finally, from Figs. 5(a) and 5(b), it is easy to check that for w 5 = 0.375 with the rest of parameters given by equation ( 4), the solution of system (3) approach to chaotic attractor at w 11 = 0.03 while it approaches to a periodic attractor for w 11 = 0.06.Fig. 5(a) shows that PD bifurcations appear as parameter w 5 is increased.The system displays chaotic dynamics when this parameter takes a value more than 0.3.The system does not support either regular or chaotic dynamics beyond 0.35 and below 0.4 (cf.Fig. 5(b)).Fig. 6(a) shows that the system dynamics shuttles between regular and chaotic motion.This shuttling is caused by changes in either of the two parameters (w 5 and w 11 ).The size of the chaotic attractor increases with the decrease in parameters w 5 and w 11 .The basin boundary of this chaotic attractor with that of the attractor shown in blue is complex.The other attractors have basins riddled with basins of several co-existing attractors.The complexity of the basin boundary structure decreases as one move to the bottom-left point in the 2D parameter scan (Fig. 6(a)).Fig. 6(b) suggests that a change in parameter w 5 affects a transition in the dynamical behavior of the system.Basin boundaries shown in Fig. 8 suggest that the system dynamics is sensitive to changes in initial conditions.This implies that stochastic external influences dictate the dynamical behavior of this hybrid system.Since chaotic dynamics exists in a narrow strip (cf.Fig. 6(c)), w 5 and w 6 are two sensitive parameters of this model system.The size of the chaotic attractor decreases when values of both the parameters slide down.Number of co-existing attractors is diminished for the top-middle point (Fig. 9(b)).The dynamical behavior of this system under the influence of exogenous factors may be unpredictable as intertwined basin boundaries are common.
We also present 2D scan diagrams in various parameter spaces (see Figs. 6(a)-6(c)) and the basin boundary calculations for chaotic attractors of the model system (3) with respect to various parameter spaces (see Figs. 7-9).In most of the figures we have presented the yz-view and in some figures the xz-view of the basin boundary structure of chaotic attractor (shown in yellow color).The basin boundary calculations are performed using the basins and attractors structure (BAS) routine developed by Maryland Chaos group.We have used the dynamics software package of Nusse and Yorke [9] for all the basin boundary calculations.It is clear from these figures that basin boundaries of the chaotic attractor are fractal which show the dynamical complexities of the hybrid model system (3).It is also seen that basin of attraction of different attractors are intermixed.The encroachment into the basin of chaotic attractor by basin of attractor at infinity (shown in green color) can be observed in Figs. 7, 8 and 9.It appears between the first attractor (shown in green color) and its basin (shown in sky blue color).The interesting feature in the model system ( 3) is that the riddled basin with fractal boundary lies in the basin of repeller which has many rectangular and square holes created by chaotic attractor.This complicated basin boundary structure suggests that the system dynamics may have loss of even qualitative predictability in the case of external disturbances.Distribution of points in the parameter space (c.f.Fig. 6(a)) suggests that for model system (3) displays STRC.It is characterized by chaotic bursts repeated at unpredictable intervals.The reason for existence of chaos at discrete isolated parameter values is that chaotic dynamics is abruptly terminated by period-doubling reversals.Figs.7(a)-7(c) show that attractors with different geometries coexist at the same set of parameter values.This is known as multi-stability in non-linear dynamics.Basin boundaries are riddled.This suggests that the system loses qualitative predictability if it interacts with external stochastic influences.Chaotic dynamics is confined to a wedge-shaped region in the parameter space (cf.Fig. 6(c)).This is certainly not a manifestation of STRC.The robust chaos is localized to a narrow region of the parameter space.The system switches to other kinds of attractors (mainly periodic) once the parameter values move outside the narrow region.Figs.9(a) and 9(c) present more complex basin boundary structures.The complexity of the basin boundaries is reduced at the parameter value, which lies in the middle of the wedge shaped region.External stochastic perturbations would have reduced effect on the system's dynamics in this case.

Discussion and conclusion
In an attempt to understand difficulties in detecting chaos in natural populations, present authors have proposed a theory [6,7].These authors discovered that deterministic chaos manifests itself as short-term recurrent chaos.The chaotic behavior is interrupted by non-chaotic dynamics at unpredictable intervals and we thought that a special kind of bifurcation process was responsible for STRC [10].Period-doubling (PD) and perioddoubling reversing (PDR) bifurcations generate and terminate chaotic dynamics in this system (cf.Figs.5(a) and 5(b)).Fig. 5(a) shows that sustained chaos exists in the range w 5 ∈ (0.31, 0.39).Short-term recurrent chaos exists in the range (0.225, 0.255).Fig. 5 In [7], we have studied the model system in which both the intermediate and top predators have Holling type II functional response.The chaotic dynamics is sensitive to changes in both the parameters: per capita death rate of the intermediate predator and the per capita intrinsic growth and death rate of the top predator.In case of the present model system in which the top predator has a CM functional response, the chaotic dynamics is distributed over larger regions of the parameter spaces (Figs.6(a)-6(c)).In these figures the involved system parameters have the following meanings: w 5 is the ratio of the specific death rate of the intermediate predator to the per capita rate of the self-reproduction of the prey.w 4 is the ratio of the half-saturation constant to the carrying capacity of the environment for the prey.In ecological terms, chaotic dynamics for the present system means that prey and intermediate predator populations may go to extinction in the event of a catastrophe.The top predator z is always a survivor (cf.Fig. 1).In contrast to this, this is not the case when the constituent populations are interacting through Holling type II functional response.Even this predator can go to extinction in the event of such catastrophes.
The present study suggests that deterministic changes in system parameters can cause transitions between ordered and chaotic system states (Figs.6(a)-6(c)).Ordered states are represented by stable foci and stable limit cycle attractors.The ratio of the per capita death rate of the intermediate predator to per capita growth rate of its prey is a crucial parameter of this model system.Another sensitive parameter is the ratio of per capita death rate of the top-predator to per capita growth rate of the prey.Another source of dynamical complexity is abrupt changes in initial conditions of the system, which might be brought in by ecological catastrophes (e.g., forest fires, flood, drought, etc.).
Non-linear dynamics and disturbance ecology have unraveled a suite of behaviors including multiple equilibria, basins of attraction and transient dynamical behavior.Analyzing a data set accumulated over a period of twelve years subsequent to catastrophic disturbance of a rain forest located on the western seaboard of Nicaragua, Vandermeer et al. [12] concluded the existence of multiple basins of attraction and non-equilibrium community structure in this forest system.The post-succession dynamics of a forest can be studied quantitatively.One measure which is particularly useful in disturbance ecology is resilience.Resilience of a system consists in the magnitude of a perturbation the system can receive without changing its structure and function.

Fig. 2 .
Fig. 2. Bifurcation diagram of the system (3), the successive maxima of z as a function of w11 is plotted in the range 0 < w11 < 0.1 for the parameters given in equation (4).
[11]stablishes that PD and PDR bifurcations are responsible for STRC.PD bifurcations initiate chaos and period reversing bifurcations terminate it in ecological time.It is surprising that PD and PDR are capable of generating this behavior for realistic parameter values.The results reinforce an earlier suggestion that period-doubling reversals could control chaotic dynamics in ecological models.Stone[11]has studied dynamics of logistic map in the presence of immigration and found that period-doubling reversals are responsible for suppression of chaotic dynamics.Fig.7suggests that chaotic dynamics is sensitive to value of the ratio of per capita death rate of the intermediate predator to the value of the per capita growth rate of the prey.Figs.7(b)and7(c)presentlesscomplex basin boundary structure.Figs.7(a)-7(c)giveanidea of complexity of basin boundaries of different attractors.Basin boundaries are riddled and interwoven.This gives rise to unpredictability in system's dynamics.Basin Boundaries presented in Figs.7-9confirm that a dominant source of unpredictability in this model system is external influences.