Qualitative behavior of three species food chain around inner equilibrium point : spectral analysis ∗

Abstract. This work deals with analytical investigation of local qualitative temporal behavior around inner equilibrium point of a model for three species food chain, studied earlier by Hastings and Powel and others. As an initial step towards the spectral analysis of the model, the governing equations have been split into linear and nonlinear parts around arbitrary equilibrium point. The explicit parameter dependence of eigenvalues of Jacobi matrix associated to the linear part have been derived. Analyzing these expressions in conjunction with some pedagogical analysis, a lot of predictions on stable, unstable or chaotic change of species have been highlighted. Agreement of predictions of this work with available numerical or semi-analytical studies suggest the utility of analytical results derived here for further investigation/analysis of the model as desired by earlier works.


Introduction
One aspect to study the dynamics of ecological community is its food web and the interaction among constituent species.Field and laboratory studies are hard to design their evolution.As a consequence, dynamical behavior of food web are studied with the help of mathematical models governed mostly by nonlinear differential or difference equations.Theoretical studies of food webs in ecology begins with the study of ditrophic food chains after the pioneering work of Lotka and Volterra [1,2].However, an ecological system with only two interacting species can account for only a small number of phenomena that are commonly exhibited in nature.In order to explain more and more new observations, such model thus extended subsequently by introducing additional number of variables as well as varieties of nonlinear interaction terms coupled with some parameters which control the evolution of the system.Following this spirit, Maity et al. [3], Upadhyay and his coworkers [4,5] recently studied the stability and bifurcation behavior of models describing evolution of three species food chain consist of Holling type response functions.From these study it appears that the presence of increasing number of variables and the additional parameters raise considerable problems for analysis of the model, both for the experimenter and theoretician.Nonetheless, such models need to be analyzed in detail since those models are found to be efficient to explain more exotic process of food web.At the begining of nineties, Hastings and Powell (HP) [6] produced an example of a chaotic system in a three species food chain modelled with type II functional responses and found to be most widely accepted [7][8][9][10][11][12][13] for defining the movement of aquatic ecosystem from its equilibrium state to chaotic states.
On the other hand Chattopadhyay and Sarkar [14] modified HP model [6] by introducing extra mortality term in one of the three species and showed that the chaotic behavior less likely occurs in a real food chain dynamics.In their works the solutions of the nonlinear equations governing the time development of three species are determined numerically by choosing some initial values.The choice of the initial value settles-down to a steady state solution depicting a situation which is at the limit cycle oscillation of the system around an equilibrium point.
In the recent past, Lonngren, Bai and Ahmet [15] showed that the system described by the HP model can be treated as the "master" in the synchronization with another "food chain" system, commonly known as "slave", that is described with a similar set of equations but containing a different set of model parameters.In their investigation they used both the systems have an identical form with one of the numerical parameters that appears in the model assuming a different value and noted that, the results would critically depend upon the values that were chosen.Their final observation was that, the dynamical behavior of the model needs to be investigated further in order to gain a complete understanding and control of the chaotic behavior.
The present work is the first step in the direction of complete understanding of the temporal behavior of the model as desired by Longrenn et al. [15] and to carry out the investigations which are yet to be explored as pointed out by Klebanoff and Hastings in their work [11], viz., desireness of the extension of their partial analysis around the equilibrium point on the plane face of the positive quadrant to inner equilibrium points and for wide range in parameter space which may support to describe other food webs, although the algebraic manipulation become more complex.
Starting from the HP model in terms of dimensionless variables and parameters, local qualitative behavior of the system around the inner equilibrium point has been studied rigorously and compared with results obtained by numerical methods in Section 2. The theoretical prediction of Sil'nikov chaos in certain domain in parameter space and its comparison with the result obtained by numerical methods are discussed in Section 3. The last section deals with the conclusion of this work.
2 Behaviors of the solution of the system around inner equilibrium point (x 0 > 0, y 0 > 0, z 0 > 0) The basic mathematical model to be studied here is a system of three first order coupled nonlinear ordinary differential equations with non-polynomial vector field describing the temporal changes of three species of food chain.In the dimensionless units, the HP model [6] is given by where the definition of dimensionless variables x, y, z; parameters even the equilibrium points of the system of Eq. ( 1) are discussed rigorously by Klebanoff and Hastings [11].
For the spectral analysis of the system around an arbitrary equilibrium point (x 0 , y 0 , z 0 ), we translate the origin of the variables (x, y, z) to (x 0 , y 0 , z 0 ) through the transformation and get the equations for new variables ξ, η, ζ as   ξ η ζ with and Here, Eqs. ( 3)-( 5) describe the behavior of the dynamical system, Eqs.(1), around the equilibrium point (x 0 , y 0 , z 0 ) irrespective of any restriction on parameter space R 6 .
The inner equilibrium point of the dynamical system described by Eq. ( 1) is In order to study the behavior of the system around the inner equilibrium point, E * , we first note that this point will be ecologically meaningful if the parameters satisfy the conditions The Jacobi matrix at E * is given by To find the eigenvalues of J * we recast the characteristic equation for it in the form where the coefficients A 1 , A 2 , A 3 are given by Transforming λ → ρ = λ + A 1 , one gets the characteristic equation in the form of standard cubic as The explicit expressions for coefficients H and G in terms of parameters are given by, and The expression is known as discriminant of the cubic equation, plays the decisive role to determine the nature of roots of the Eqs.( 9) or (11).Following the standard procedure for solving cubic equation one can obtain the expression for roots of ( 9) in terms of parameters through their presence in H, G, ∆ as where 2 and ω 2 = ω * are the cube roots of unity.

Linear stability analysis
The requirement of the equilibrium point is to be node or saddle is that all the roots of the characteristic equation, Eq. ( 9) are to be nonzero real and distinct and, the point is focus or center if Eq. ( 9) has a pair of complex roots.From the pedagogical analysis of Eq. ( 9) one can thus conclude: Proposition 1.The equilibrium point E * of Eq. ( 1) is to be node or saddle (focus or center), i.e., hyperbolic if (i) the discriminant ∆ of ( 14) is strictly less (greater) than zero, and, The proof of this proposition can be outlined by mentioning the relation between roots and the discriminant as: if λ 1 , λ 2 , λ 3 are three roots of Eq. ( 9), then From ( 16), it is clear that if G 2 + 4H 3 < 0, all the roots are real, distinct and the system is known as irreducible case of Cardan's solution.In this case the solutions of Eq. ( 9) can be found as In addition, condition (ii) ensures non-vanishing of all the roots of Eq. ( 9).Consequently, the ecologically meaningful equilibrium point E * will be node or saddle for the parameters lying in the domain (7).As a corollary one can further conclude that Proposition 2. The system of equations in (1) is locally conjugate, i.e, have the same qualitative behavior) to the linear equation those points in the parameter space where the eigenvalues are resonant.
The proof of this proposition follows from Poincaré's theorem.We have kept the analysis of the resonant case beyond this work.
Propositions 1 and 2 mentioned above deal with situation when the eigenvalues of the Jacobi matrix J are distinct and nonzero.In this case the qualitative behavior of the solution of the system of Eqs.(1) changes smoothly for infinitesimal variation of parameters.But this feature does not maintained whenever real part of one or more eigenvalues of J are zero.The smooth variation of parameters around such critical values of parameters for which eigenvalues of J moves around both sides of imaginary axis, produces entirely different qualitative behavior of the solution of the system, commonly known as bifurcation.

Bifurcation behavior
A bifurcation is said to be a steady state bifurcation (SSB) [20] whenever a eigenvalue of the Jacobi matrix J * at the equilibrium point is real and vanishes and the real part of all other eigenvalues are non zero.Proposition 3. The system of equations, Eq. (1) undergo steady state bifurcation around the equilibrium point E * provided the parameters contained in the domain The proof is straightforward since the first inequality constraints the roots are to be real and distinct and, the second condition suggests that one root must be zero.Fig. 1(a) demonstrates the variation of a 1 and b 1 satisfying the conditions A 3 = 0 and ∆ < 0 for the fixed values of a 2 , b 2 , d 1 , d 2 at 0.1, 2.0, 0.4, 0.01, respectively.
The Takens-Bagdanov bifurcation [20] is a bifurcation of an equilibrium point of a system of autonomus ordinary differential equations at which the Jacobi matrix has a zero eigenvalue of multiplicity two.Following the similar reasoning one can easily state that: Proposition 4. The system of equations in (1) undergo Takens-Bogdanov bifurcation (TBB) around the bifurcation point E * provided the critical values of the parameters belongs to the subset This statement can be proved by using the relation (16) between roots and the discriminant or putting the values of A 2 and A 3 into Eq.( 9).The occurrence of Hopf bifurcation (there is a pure imaginary pair of eigenvalues) of Eq. ( 1) can be identified from the following proposition: The proof of this statement follows from the fact that the first condition, viz., the inequality confirms the existence of pair of complex roots.However, the real part of that pair may not necessarily vanish.The non-vanishing of the coefficients A 1 , A 2 , A 3 as well as the discriminant ensure the existence of a non-zero real root.Finally, the equation A 3 = 9A 1 A 2 ensures the vanishing of the real part of pair of complex eigenvalues which can be justified by comparing the coefficients of like powers of λ in the product (λ − δ)(λ − iω)(λ + iω) with the L.H.S. of Eq. ( 9).Fig. 1(c) exhibits partially region D P AHB of critical values in a 1 − b 1 space for parameters a 2 , b 2 , d 1 , d 2 fixed at (0.1, 2.0, 0.4, 0.01) and variable parameters a 1 , b 1 within the limit 0 < a 1 < 8, 0 < b 1 < 6.From this figure it is evident that whenever 0 < a 1 < 5.75, curve of bifurcation points in parameter space consists of two branches whereas for 5.75 < a 1 < 8, it consists of four branches.So with the help of equation A 3 = 9A 1 A 2 one can easily extract the critical values of parameters associated to the bifurcation points in the codimension two parameter space around E * .The occurrence of Hopf bifurcation (there is a pure imaginary pair of eigenvalues) of Eq. ( 1) can be identified from the following proposition: The proof of this statement follows from the fact that the first condition, viz., the inequality confirms the existence of pair of complex roots.However, the real part of that pair may not necessarily vanish.The non-vanishing of the coefficients A 1 , A 2 , A 3 as well as the discriminant ensure the existence of a non-zero real root.Finally, the equation A 3 = 9A 1 A 2 ensures the vanishing of the real part of pair of complex eigenvalues which can be justified by comparing the coefficients of like powers of λ in the product (λ − δ)(λ − iω)(λ + iω) with the L.H.S. of Eq. ( 9).Fig. 1(c) exhibits partially region D P AHB of critical values in a 1 − b 1 space for parameters a 2 , b 2 , d 1 , d 2 fixed at (0.1, 2.0, 0.4, 0.01) and variable parameters a 1 , b 1 within the limit 0 < a 1 < 8, 0 < b 1 < 6.From this figure it is evident that whenever 0 < a 1 < 5.75, curve of bifurcation points in parameter space consists of two branches whereas for 5.75 < a 1 < 8, it consists of four branches.So with the help of equation A 3 = 9A 1 A 2 one can easily extract the critical values of parameters associated to the bifurcation points in the codimension two parameter space around E * .To illustrate the utility of analytical results derived here representative of some of our theoretical predictions of node, saddle, focus, center etc., of HP model are summarized in Table 1 in some domain in the parameter space.For comparison of our result with the existing results obtained by numerical techniques we fix the parameters (a 2 , b 2 , d 1 , d 2 ) at (1.0, 2.0, 0.4, 0.01), change b 1 discretely by 0.5, 2.2,and 3.0 and vary a 1 continuously in appropriate range depending upon the choice of b 1 .From Table 1 it appears that for b 1 = 0.5, the domain of a 1 is (0.62229 ≤ a 1 ≤ 7.777) and six bifurcation points within this domain in the parameter space.Out of them two are SSB points and the other four are PAHB points.Depending upon the sign of the real part of the eigenvalues as obtained from (15), the nature of the equilibrium points have been predicted as: it is stable nodes in (0.6223 ≤ a 1 ≤ 0.8897) and saddle within (0.8897 ≤ a 1 ≤ 1.5650) and so on.
Till now we have investigated the qualitative behavior of deterministic nature of temporal evolution around E * of species described by HP model.Observing the appearance of their chaotic change in the earlier semi-analytical studies [11][12][13] around the equilibrium point on the plane face of the positive quadrant of the state space, we will continue our spectral analysis for investigating the chaotic motion around the inner equilibrium point E * .

Prediction of chaos from spectral analysis
The characteristics of chaos and its presence in nature are much discussed in ecology [6,7,[16][17][18].The most important mathematical attribute of chaos is the absence of any stable equilibrium point or any stable limit cycle in system dynamics, for which the pattern never repeat themselves.Recent developments in dynamical system theory consider chaotic fluctuations of a dynamical system as highly desirable because fluctuations allow such a system to be easily controlled.We now intend to show that one can identify the domain in the parameter space for which the HP model may exhibits chaotic change in its state space for any choice of the parameter in that identified region.
The Sil'nikov criterion for the existence of chaotic motion (known as Sil'nikov chaos) [19,20] is: if the equilibrium point x e of the system is a saddle focus and the eigenvalues γ and α ± iβ satisfy the Sil'nikov inequalities: then the system may have Sil'nikov chaos in some neighbourhood of equilibrium point x e .Using expressions (15), the conditions in ( 19) can be restated as Proposition 6.The motion of the system of equations Eq. ( 1) is (Sil'nikov) chaotic in some neighbourhood of the inner equilibrium point E * if the parameters satisfy the conditions: and where H, G and ∆ are given in ( 12)- (14).
Proof.(i) The condition for the equilibrium point x e to saddle focus is that Eq. ( 9) must has a pair of complex root and a non zero real root.Consequently, relation (16) ensures that ∆ > 0.
(ii) Substitution of explicit expression for ω 2 ) and ω 2 (= ω * ) in expressions (15) and separation into real and imaginary parts gives the imaginary part complex eigenvalues around x e as (− G 2 + √ ∆ 2 ) . Therefore, the condition (ii) of Sil'nikov criterion yields the condition (ii) of Proposition 6.
(iv) With the help of (15a) and (21a), the absolute value of the real part of the complex roots of Eq. ( 9) can be written as Therefore, |γ| − |α| ≥ 0 implies the condition (iv) stated in Proposition 6.

Conclusion
In the process of study of the qualitative behavior of temporal evolution around the inner equilibrium point of a three species food chain described by HP model, we have first derived the dependence of eigenvalues of Jacobi matrix on the parameters involved in the system.From the judicial analysis of the conditions on parameters presented in

Conclusion
In the process of study of the qualitative behavior of temporal evolution around the inner equilibrium point of a three species food chain described by HP model, we have first derived the dependence of eigenvalues of Jacobi matrix on the parameters involved in the system.From the judicial analysis of the conditions on parameters presented in Propositions 1 to 5, the domain in the parameter space can be identified for desired temporal evolution around the equilibrium point of the system.In particular, Propositions 1 and 2 help to choose parameters for which the system of nonlinear Eqs. ( 1) is conjugate (topologically equivalent) to linear one around the equilibrium point.The set of parameters for which the system undergoes bifurcation behavior can be identified through the conditions presented in Proposition 3 to 5. Our predictions have been summarized in Table 1 and compared for a selective set of parameters for which the behavior of the system have been studied by using numerical techniques.Agreement of results for such parameter values suggest, condition on parameters derived here can be utilized for the prediction of the qualitative behavior of the temporal evolution of the species for other domain in parameter space.Furthermore, the domain in the parameter space for the Sil'nikov chaotic motion of the system can be identified from the conditions on parameters derived in Proposition 6.The reliability of our prediction has been tested by numerical solution of the system of Eqs.(1).Consequently, a knowledge provided here may be quite helpful in risk management of complex and highly interdependent system found in an environment.One may further initiate the systematic studies of the normal form reduction [22], limit cycle and other exotic behaviors like synchronization of chaos of HP model around inner equilibrium point from the results reported here.

Fig. 1 (
Fig. 1(b) demonstrates the variation of a 1 and b 1 satisfying the conditions A 3 = 0 and A 2 = 0 for the fixed values of a 2 , b 2 , d 1 , d 2 at 0.1, 2.0, 0.4, 0.01 respectively.From the acute observation of this figure it reveals that the Takens-Bogdanov bifurcation will occur for the parameter values a 1 = 8 and b 1 = 1.0004 in addition to the values of other parameters as stated above.Consequently, one may predict other values of a 1 and b 1 for given other set of parameter values a 2 , b 2 , d 1 , d 2 either by solving or plotting the criterion A 3 = 0 and A 2 = 0 simultaneously.

Proposition 5 .
The system of equations, Eq. (1) undergo Poincaré-Andronov-Hopf-(Hopf-steady state) bifurcation (PAHB) around the bifurcation point E * whenever the critical parameter values are contained in the domain D

Proposition 5 .
The system of equations, Eq. (1) undergo Poincaré-Andronov-Hopf-(Hopf-steady state) bifurcation (PAHB) around the bifurcation point E * whenever the critical parameter values are contained in the domain

Table 1 .
Nature of the bifurcation points and the nature of the orbit in between critical