Complex Dynamics of a Three Species Food-chain Model with Holling Type Iv Functional Response

In this paper, dynamical complexities of a three species food chain model with Holling type IV predator response is investigated analytically as well as numerically. The local and global stability analysis is carried out. The persistence criterion of the food chain model is obtained. Numerical bifurcation analysis reveals the chaotic behavior in a narrow region of the bifurcation parameter space for biologically realistic parameter values of the model system. Transition to chaotic behavior is established via period-doubling bifurcation and some sequences of distinctive period-halving bifurcation leading to limit cycles are observed.


Introduction
Ecological systems (natural systems, prey-predator communities etc.) are dynamic, complex and non-linear in nature.In non-linear systems, the effect is not proportional to the cause.This non-linearity in the presence of dissipation (frictional forces in physics and non-trophic interactions in ecology) gives rise to different attractor types as representations of asymptotic states of these dynamical systems.The attractors which are frequently observed in the model ecological systems are: stable foci, stable limit cycles and chaos [1].These are obtained in the state space (or phase space) of an ecological system when an intrinsic attribute (control parameter) of the system is varied; e.g., the quantity and quality of food supplied to a chemostat.If the first two are most frequently observed in a model ecological system, the system's dynamics is classified to be regular.On the other hand, if chaotic attractors are detected in a dense region in two-dimensional parameter space, then the system is termed to be chaotic.Chaotic states are characterized by exponential sensitivity of system's dynamics with respect to initial conditions which can not be fixed or determined without committing a finite amount of error either in the laboratory or in the field.In a recent work [2], we studied how different attractor types c Vilnius University, 2011 appear and disappear in a suit of five model ecological systems when control parameters are varied in two-dimensional parameter spaces.The sensitivity of the dynamics of model ecological systems to variations in control (key) parameters are gauged in this way.But, this is just one aspect of the dynamical behavior of a model ecological system.In the natural world, dynamical behavior exhibited by various ecological systems is very complicated.For instance, the number of wild animals and plants are continuously variable quantities and these variations are irregular in period and also in amplitude [3].Prince et al. [4] remarked that the study of community behavior with the help of mathematical models must be based on at least three trophic levels and hence more focus should be made to study the complex behavior exhibited by deterministic models consists of three and more trophic levels.Hastings and Powell [5] examined the complex nonlinear behavior of three-species continuous-time ecological model.His study initiated a concerted effort among theoretical ecologists who sought to analyze the subtle dynamics of these nonlinear multi-trophic models [6][7][8][9][10].
Predator-prey communities are embedded in complex food webs.These webs are woven by trophic (predation) as well as non-trophic interactions (e.g., competition, mutualism and interference).The latter forms weak links with food chains in the web and is important as far as the stability of food webs are concerned.In the food chain the predator which appears at the bottom of the chain is known as top-predator.The response of food webs to external perturbations (e.g., forest fire, epidemics, invasion by alien species, etc.) is dependent on the structure of interconnected food chains.In an earlier work [11], we have shown that the response of a food chain to exogenous environmental perturbations is dependent on whether the top-predator is a specialist or generalist.
Various mathematical techniques, like bifurcation analysis, extensive numerical simulations have been used to detect chaotic dynamics in multi-dimensional deterministic models of ecological systems.There is no unique mathematical tool to detect the parametric domain for which the model under consideration will exhibit chaotic oscillations and also what type of non-linear coupling terms are required to be present within the mathematical models which have potential to exhibit chaotic dynamics.
In the history of population ecology, both mathematicians and ecologists have a great interest in the Holling type predator-prey models [12][13][14]; including Holling type I-III, originally due to Holling [15] and Holling type IV, suggested by Andrews [16].The Holling type IV functional response is of the form As the parameter i (which is a measure of the predator's immunity from or tolerance of the prey) decreases, the predator's foraging efficiency decreases.In the limit of large i, it reduces to a type II functional response [17].The parameters w and a can be interpreted as the maximum per capita predation rate and the half saturation constant in the absence of any inhibitory effect.This response function describes a situation in which the predator's per capita rate of predation decreases at sufficiently high prey densities.Using the Floquet theory of impulsive equations and small perturbation, Baek [18] obtained conditions for the stability analysis of a food chain system with Holling type IV functional response.
Lian and Xu [19] investigated the affects of time delay on the dynamics of the food chain systems.Existence of local Hopf bifurcations is analyzed by choosing the delay as a control parameter.Shen [20] obtained the sufficient conditions for the permanence and global attractivity of the food chain system with Holling type IV functional response are obtained using comparison theorem.Van Gemerden [21] fitted a type IV functional response to the uptake of hydrogen sulphide by purple sulphure bacteria.Collings [22] used the response function in a mite predator-prey interaction model and called it a Holling type IV function.
Recently Upadhyay et al. [23] studied the model system under the influence of environmental driving forces.These driving forces are unable to drive the system from a regime of deterministic chaos towards a stochastically stable situation.In the present paper, we investigate the existence of complex dynamical behavior and persistence in a three trophic levels food chain involving Holling type IV functional response.In the next section, we present details of the model system.The existence of solution and the dissipativeness are presented in Section 3. The stability analysis of the model system is carried out and the persistence conditions are established in Section 4. Section 5 enumerates results of our simulation experiments.Discussions and conclusions are presented in the last section.

The food chain model system
The three species food chain models are consisting of one prey (x) and two predators (y and z).The middle predator y feeds on the prey x and the top predator z preys upon y, both according to the Holling type IV functional response.The three species food chain model can be represented mathematically by the following system of differential equations: where a 1 , b 1 , a 2 , w, w 1 , w 2 , w 3 , c, i, i 1 and a are positive constants.a 1 is the intrinsic growth rate of the prey population x, a 2 be the intrinsic death rate of predator population y in the absence of the only food x, the parameter c is the decay rate of the top predator z in the absence of its prey y and the ratio w 3 /w 2 is measure of the assimilation efficiency.w is the maximum value which per capita reduction rate of x can attain, w 1 has similar meaning to w. w 2 is the maximum value which per capita reduction rate of y can attain, w 3 has similar meaning to w 2 .b 1 is the strength of intra-specific competition among the prey species x.The parameter a can be interpreted as the half-saturation constant in the absence of any inhibitory effect.The parameters i and i 1 , in turn, are a direct measure of the middle and top predator's immunity from, or tolerance of the prey respectively.
Clearly the model system (1) has eleven parameters in all and it can be written as the following form: 3 Existence and dissipativeness Obviously the interaction functions g i (i = 1, 2, 3) of the model system (2) are continuous and have continuous partial derivatives on R 3 + = {(x, y, z) ∈ R 3 : x ≥ 0, y ≥ 0, z ≥ 0}.Therefore the solution of the model system (2) with non-negative initial condition exists and is unique, as the solution of system (2) initiating in the non-negative octant is bounded.Further more, the system is said to be dissipative if all population initiating in R 3  + are uniformly limited by their environment [24].Accordingly, the following theorem establishes the dissipativeness of model system (2).
Proof.From Eq. (2a) of the model system (2) we have x(0) − b 1 be the constant of integration.Hence for the large value of time we get Let m(t) = w1 w x + y + w2 w3 z, hence it is easy to verify that dm dt where p = min(a 2 , c).So by comparison lemma we obtain, for all t ≥ T ≥ 0, Thus m(t) = w1 w x + y + w2 w3 z ≤ w1a1(a1+p) wb1p , and then all species are uniformly bounded for any initial value inR 3  + .Note that, for a biologically realistic model system (2) has to be dissipative (i.e., all population are uniformly limited in time by their environments).Therefore, according to the above theorem we assume that there exists (η 1 , η 2 , η 3 ) > 0 such that Ω(x 0 , y 0 , z 0 where Ω(x 0 , y 0 , z 0 ) is the omega limit set of the orbit initiating at (x 0 , y 0 , z 0 ).Thus, the model system (2) is dissipative.

Stability analysis of 3D system and persistence
In this section, the existence and local stability analysis of the non-negative equilibrium point of model system (2) are investigated.Also the persistence condition of the system (2) is established.
There are at most four non-negative equilibrium point of model system (2).The existence and the stability condition for them are as follows: 1.The trivial equilibrium point E 0 = (0, 0, 0) always exists.
3. In the absence of top predator the middle predator can survive on its prey.Hence the equilibrium point E 2 = ( x, y, 0) exists in the interior of positive quadrant of xy-plane, where x and y are given as follows: 4. Due to the extinction scenario of middle predator there is no equilibrium point in the xz-plane.Moreover, neither y nor z can survive in the absence of prey species x, hence there is no equilibrium point in the yz-plane.
5. The positive equilibrium point E 3 = (x * , y * , z * ) exists in the interior of the first octant if and only if there is a positive solution to the following algebraic nonlinear system: Nonlinear Anal.Model.Control, 2011, Vol.16,No. 3, Straight forward computation show that x * is a positive root of the cubic equation and this equation can be written as where Therefore the positive equilibrium point E 3 = (x * , y * , z * ) exists under the following conditions: Now, in order to study the behavior of solution near the equilibrium points, we need to compute the variational matrix of the model system (2).Assume that V (x, y, z) denotes the variational matrix of the system (2) at the point(x, y, z).Then where Accordingly, the linear stability analysis about the equilibrium points E i , i = 0, 1, 2, 3, gives the following results: 1.The equilibrium point E 0 = (0, 0, 0) is a saddle point, indeed near E 0 the prey population grows while both the predator population decline.
2. The equilibrium point and then E 1 is a saddle point with stable manifold in the xz-plane and with unstable manifold in y-direction.Moreover, near E 1 the prey population remains in the neighborhood of a 1 /b 1 , while the middle predator increase when the top predator decline.
3. The root of the characteristic equation p 3 (λ) = 0 of the above variational matrix ( 8) about E 2 = ( x, y, 0) satisfy the following: The equilibrium point E 2 = ( x, y, 0) is stable or unstable in the positive direction orthogonal to the xy-plane, i.e., z-direction depending on whether λ 3 = w3 y y 2 /i1+ y+a − c is negative or positive, respectively.
4. The variational matrix of the model system (2) at the positive equilibrium E 3 = (x * , y * , z * ) is given as follows: where α = y * 2 /i 1 + y * + a, and β = x * 2 /i + x * + a.So the characteristic equation of the variational matrix ( 9) is written as: where According to the Routh-Hurwitz criterion, where Therefore, depending on the above analysis the following theorem can be proved easily.
Theorem 2. Suppose that the positive equilibrium point E 3 = (x * , y * , z * ) exists in the interior of the positive octant, and then E 3 is locally asymptotically stable provided that conditions (10)-( 13) hold.
In the following theorem, the global stability analysis of E 3 = (x * , y * , z * ) is investigated.
Theorem 3. Assume that the positive equilibrium point E 3 = (x * , y * , z * ) is locally asymptotically stable.Then it is a globally asymptotically stable in the interior of positive octant (i.e., Int R 3 + ) provided that where η 1 = a1 b1 , η 2 = w1a1 wb1 (1 + a1 4a2 ), η 3 = w1a1 wb1 (1 + a1 p ) and p = min(a 2 , c).Proof.The proof is based on a Lyapunov direct method.Consider the following positive definite function Obviously, V is a continuous function on Int R 3 + , C j (j = 1, 2, 3) are the positive constants to be determined.Now in order to investigate the global dynamics of the nonnegative equilibrium point E 3 of the model system (2) the derivative of V with respect to time along the solution of system ( 2) is computed as Since, Substituting (16a)-(16c) in ( 15) we obtain The above equations can be write as sum of the quadratics where , , , , Sufficient conditions for dV dt to be negative definite are that the following inequalities hold: Since a 13 = 0, condition (20) is automatically satisfied.It can be seen that under condition (14a), a 11 > 0. Under condition (14b), a 22 > 0; and under condition (14c), 1 /i+η1+a) and C 3 = C2w2y * aw3z * , then it can be checked that conditions (18) and ( 19) are automatically satisfied.
In the following, we shall find the conditions for the persistence of the food chain system (2).Theorem 4. Assume that there are no non-trivial periodic solutions in the xy-plane, and then the necessary condition for the persistence of system (2) is and the sufficient condition for the persistence of the system is Proof.Since the boundedness is proved and λ 3 is the eigenvalue which gives the stability in the positive direction orthogonal to the xy-plane (i.e., z-direction).Hence if there are no non-trivial periodic solution in the xy-plane and condition (21a) is violated (i.e., λ 3 < 0).Then E 2 is a stable equilibrium point, which means that there is orbit converges to it in the positive cone.Therefore, condition (21a) is the necessary condition for the persistence.Now for the sufficient condition of persistence of the model system (2), we shall apply the abstract theorem of Freedman and Waltman [25].According to the growth functions g 1 , g 2 and g 3 of system (2) the following hypothesizes are satisfied: (A2) The prey species x grows to carrying capacity in the absence of predator that is g 1 (0, 0, 0) = a 1 > 0 and g 1 ( a1 b1 , 0, 0) = 0. While, due to the intraspecific competition within prey species we have ∂g1 ∂x (x, 0, 0) = −b 1 < 0. (A3) There is no equilibrium point in xz-plane and yz-plane.
(A4) The intermediate predator can survive on its prey in the absence of top predator, that is there exist an equilibrium point E 2 = ( x, y, 0) in the xy-plane such that g 1 ( x, y, 0) = 0 = g 2 ( x, y, 0).However the top predator cannot survive on the prey x in xy-plane.
Theorem 5. Suppose that condition (21b) holds and there are a finite number of limit cycles in the xy-plane.Then for each limit cycle (θ(t), v(t)) in the xy-plane the persistence condition for system (2) takes the form Here T is the time period of the limit cycle.Proof.Assume that there exists a limit cycle in the xy-plane, then the variational matrix about the limit cycle, x(t) = θ(t), y(t) = υ(t), z(t) = 0, take the form where all the partial derivative and g j (j = 1, 2, 3) in ( 23) are computed at the limit cycle (θ(t), υ(t), 0).Consider now a solution of the model system (2) with positive initial condition (t, α 1 , α 2 , α 3 ) sufficiently close to the limit cycle.
From the variational matrix ( 23), ∂z ∂α3 is a solution of the system dz dt = g 3 (θ(t), υ(t), 0)z with z(0) = 1.Thus, Hence by using Taylor expansion, we have Then z increase or decrease according to T 0 g 3 (θ(t), υ(t), 0) dt is positive or negative, respectively.Since E 3 and those limit cycles are the only possible limit in the xy-plane of trajectories with positive initial condition, hence the trajectories go away from the xy-plane if conditions (21b) and ( 22) hold.

Numerical results
The global dynamical behavior of the non-linear model system (2), in the positive octant, is investigated numerically.A numerical integration for model system (2) is carried out for various choices of biologically feasible parameter values and for different sets of initial conditions.In all the cases being considered here the data is chosen such that the persistence conditions are satisfied.Therefore, as the solution for system is bounded we expect that model system (2) have a rich dynamics including limit cycle, quasi-periodic or even chaotic dynamic.Model system (2) is solved numerical using the predictor-corrector method with sixth order Runge-Kutta method [26].
The chaotic attractor and its corresponding time series of the model system (2) are decided on the following data set (see Fig. 1):             Phase portrait in the xy, yz and xz-planes of a strange attractor are presented in Figs.1(a)-1(c).Model system (2a)-(2c) display sensitivity to the initial condition for the small change in population densities are presented in Figs.1(d)-1(f).The prey and predator densities plots are redrawn for a small change of size 0.0001 in the initial condition, which intersects the initial plot significantly.The temporal variations in prey and predator densities show highly irregular oscillations in the model system.This behaviour is observed when the time history of the population density is plotted with small changes in the initial condition.The two time trajectories intersect each other confirming the dynamics of the system to be chaotic for that particular set of parameter values.Therefore, a detailed study of the system dynamics must include the SIC property of the system leading to confirmation of chaos.
The data set given in Eq. ( 24), we only change the value of the strength of intraspecific competition among prey species, b 1 from 1.0 to 1.5, we obtain the limit cycle attractor of period 2 and its corresponding time series are presented in Fig. 2. The limit cycle attractor of period 1 and its corresponding time series of model system (2) are presented in Fig. 3 for b 1 = 1.75.The values of other parameters are same as given in Eq. ( 24).
If we set the strength of intra-specific competition among prey species at the value b 1 = 2.2 and the value of other parameters are same as given in Eq. ( 24), we obtain the stable focus for the model system (2) (see Fig. 4).
The simulation experiments were done to determine the regions in the parameter spaces, which support different dynamical behaviors for the model system (2).The computed results are given in Table 1, which are the main results of the complex dynamical behavior of this model system.From Table 1, it is found that chaos was observed in the ranges 1.4 ≤ a 1 ≤ 2.7, 0.7 ≤ b 1 ≤ 1.4, 1.9 ≤ w 1 ≤ 2.3, 0.08 ≤ a 2 ≤ 0.089, 0.093 ≤ a 2 ≤ 0. Extinction of top-predator (z) is also observed for most of the parameter values and in some cases the extinction of the middle predator y is also observed.For fixed i (i = 0.3) and i 1 ∈ [0.2, 4.0], chaos is observed but if we fixed i 1 (i 1 = 0.3) then chaos is observed in the range i ∈ [0.3, 0.9] and stable limit cycle is observed in the range        Under the bifurcation analysis of the model system (2), very rich and complex dynamics are observed, presenting various sequences of period-doubling bifurcation leading to chaotic dynamics or sequences of period-halving bifurcation leading to limit cycles.For bifurcation diagram of model system (2) presented in column (a) of Fig. 5, the value of the intrinsic growth rate of prey population parameter a 1 varies in the range [0.62, 2.0] and the values of other parameter are as follows: The blow-up bifurcation diagrams are presented in column (b) of Fig. 5, show that the model system possesses rich variety of complex dynamic behavior for bifurcation parameter a 1 in the range [0.63, 0.75].
The bifurcation diagrams are generated for the successive maxima of the prey population x in the range [0.001, 2.0], the intermediate predator population y and the top predator population z in the ranges [0.001, 2.5] and [0.001, 0.3] respectively as a function of intrinsic growth rate of prey population a 1 in the range 0.62 ≤ a 1 ≤ 2.0.The blowup bifurcation diagrams show that the model system possesses rich variety of dynamical behavior if the bifurcation parameter a 1 varies in the range 0.63 ≤ a 1 ≤ 0.75 and successive maxima of the prey and predator populations in the ranges [0.01, 0.8], [0.01, 0.6] and [0.005, 0.12] respectively.
The bifurcation diagrams are also generated for the successive maxima of the prey population x in the range [0.001, 1.8] as a function of strength of intra-specific competition among prey population b 1 in the range 0.64 ≤ b 1 ≤ 2.2.The blow-up bifurcation diagram is generated in the range x max ∈ [0.001, 0.7] and the bifurcation parameter b 1 in the range 1.4 ≤ b 1 ≤ 2.2.Similarly for the successive maxima of the predator populations y and z and the control parameter b 1 , the bifurcation diagrams and the blow -up bifurcation diagrams are presented in Fig. 6.The bifurcation diagrams are drawn for the parameter values as given in Eq. ( 24) with a 1 = 1, b 1 = 1, and presented by Figs. 5 and 6.These Figures (Figs. 5 and 6) show clearly the evidence of the route to chaos through the cascade of period-doubling and period-halving bifurcation respectively.The b 1 -bifurcation parameter is varied in the range [0.64, 2.2].A period-halving cascade is observed.After the accumulation point, the behavior settles down onto a chaotic attractor which is structured on a skeleton of periodic orbit [27].The period-doubling phenomena leading to chaos is a well known feature of a range of nonlinear systems of biological populations.For a 1 bifurcation parameter, a period-doubling cascade is observed which is given in Fig. 5.The increase in size of a chaotic attractor as the system parameter is varied is considered to be the hallmark of the crisis (sudden destruction of a chaotic attractor) route to chaotic dynamics [28,29].The crisis occurs precisely at the point (a 1 = 0.687, b 1 = 1.622)where the unstable period 3 orbit created at the original saddle -node bifurcation intersects with the narrow chaotic region.It shows the evidence of the route to chaos through the cascade of period halving and other important change in the chaotic set include interior crisis in which a chaotic attractor undergoes a sudden increase in the size [27] along with the appearance or sudden enlargement of a fractal structure in the basin boundary.

Discussions and conclusions
In this paper, we have investigated the dynamical complexities of a three-species food chain model with Holling type IV predator response.The boundedness of the trajectories, existence of an attracting set, as well as existence of equilibrium points is established.The local and global stability for non-negative equilibrium point has been analyzed.Condition for the persistence of the model system is obtained.Numerical simulation of model system (2) shows that, the food chain system has rich dynamics including periodic and chaotic dynamics.Extinction of the top or middle predator population is observed for the extreme values of the model parameters except for the death rate of the top-predator z for which extinction is observed in the intermediate ranges [0.34, 0.9].
From the bifurcation analysis, we observed that the asymptotic behavior of the model system is extremely sensitive to the value of intrinsic growth rate of prey population and the strength of intra-specific competition among the prey species parameters a 1 and b 1 in the ranges (0.63, 0.75) and (1.4,2.2) respectively.From Figs. 5 and 6, we observed the period-doubling route to chaos for the intrinsic growth of prey population parameter (a 1 ) and period-halving route to limit cycle for the strength of intraspecific competition parameter (b 1 ) among the prey species.Thus it is observed that even small variation in parameters a 1 and b 1 may cause a shift from limit cycles to chaos and vice-versa respectively.This study gives support to the view that multi-species model systems are able to generate unpredictable and complex behavior and points to the difficulties in understanding the observed dynamical behavior of even simple ecological systems in the absence of a reasonable a priori model of their growth and dynamics.

Fig. 1 .
Fig. 1.Phase plane diagram for the model system (2) for the parameter values given in Eq. (24) with initial condition [0.9302, 0.1795, 0.0132]: (a) phase portrait in the xyplane of a strange attractor; (b) phase portrait in the yz-plane of a strange attractor; (c) phase portrait in the xz-plane of a strange attractor; (d) sensitive dependence on initial condition (SIC) for the small change in prey population density x; (e) SIC for the small change in middle predator population density y; (f) SIC for the small change in top-predator population density z.
13 and at the discrete point a 2 = 0.091 and other parameters are fixed www.mii.lt/NA at the limit cycle value i.e., b 1 = 1.75.For this set of values no chaos was observed in the whole range of parameter c ∈ [0.01, 3.0] and w 3 ∈ [0.01, 4.0].But if we set the strength of intra-specific competition among prey species at the value b 1 = 1.0, we observe chaos for these parameter values in the ranges c ∈ [0.14, 0.3] and w 3 ∈ [1.4,4.0].

Fig. 5 .Fig. 6 .
Fig. 5. Period-doubling bifurcation of the model system (2) for the parameter values given in Eq. (24): (a) successive maxima of the population densities verses the control parameter a1 for the range 0.62 ≤ a1 ≤ 2.0; (b) blow-up bifurcation diagram of (a) for the control parameter a1 for the range 0.63 ≤ a1 ≤ 0.75.