- Email: [email protected]
Chaotic population dynamics and biology of the top-predator
Chaos, Solitons and Fractals 21 (2004) 1195–1204 www.elsevier.com/locate/chaos
Chaotic population dynamics and biology of the top-predator Vikas Rai a, Ranjit Kumar Upadhyay a
b,*
Department of Applied Mathematics, Delhi College of Engineering, Bawana Road, Delhi 110 042, India b Department of Applied Mathematics, Indian School of Mines, Dhanbad, Jharkhand 826 004, India Accepted 8 December 2003
Abstract We study how the dynamics of a food chain depends on the biology of the top-predator. We consider two model food chains with specialist and generalist top-predators. Both types of food chains display same type of chaotic behavior, short-term recurrent chaos; but the generating mechanisms are drastically different. Food chains with specialist top-predators are dictated by exogenous stochastic factors. On the contrary, the dynamics of those with the generalist top-predator is governed by deterministic changes in system parameters. The study also suggests that robust chaos would be a rarity. Ó 2004 Elsevier Ltd. All rights reserved.
1. Introduction Ecology is a parent discipline of chaos. Since the seminal work of Sir Robert May [1] model ecosystems have been designed and studied [2,3]. These studies conveyed that chaos might be a dominant mode of ecosystem evolution. As such data requirements to capture chaos in the wild are stringent. One needs a string of precise measurements of species abundance for a length of time, which cannot be fixed a priori. On contrary to this, the truth of the real world is that ecological time-series are short and noisy. This severely limits the possibility of chaos getting detected even if it is there. A couple of years ago, it occurred to us that the origin of the crisis might be intrinsic to organization of ecological interactions. Since then, we have made sustained efforts to understand why it is so difficult to observe chaos in nature. Earlier authors had suggested [4,5] that the narrowness of the boundary between chaos and noise is responsible for the crisis. All natural systems are noisy. Chaotic dynamics in appearance is noise-like. The noisy character of these systems are caused either by intrinsic factors or extrinsic forces. When a system is subjected to external influences (deterministic or stochastic) complex dynamical phenomena (chaos and noisy chaos) emerge. These phenomena are often indistinguishable. This leads to failure of attempts to observe chaos. The work [6,7] of the present authors has opened up a new possibility that the reason for chaos getting poor evidential support might lie elsewhere. Moreover, it also suggested that the biology of the top-predator has an important role to play as far as the dynamical complexity of a food chain is concerned. In the present paper, we take up this issue in detail. The paper is organized as follows. Section 2 presents the model systems that have been chosen for detailed dynamical investigations. The methodology of the present study is presented in Section 3. We discuss the main findings in Section 4 and record the messages of the study in Section 5.
*
Corresponding author. Tel.: +91-326-2200817; fax: +91-326-2210028. E-mail address: [email protected] (R.K. Upadhyay).
0960-0779/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2003.12.065
1196
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
2. Model systems Among four basic interactions (predator–prey, competition, interference and mutualism), predator–prey interaction is the most common. There are two types of predators: (i) specialists and (ii) generalists. A specialist predator is the one, which dies out exponentially fast when its favorite food is absent or is in short supply. The latter predators switch over to alternative food options when its most preferred food is in short supply. The Rosenzweig–MacArthur (RM) model [8] is the one, which describes the dynamics of a specialist predator and its prey: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ
ð1aÞ
dY w1 YX ¼ a2 Y þ ; dt ðX þ D1 Þ
ð1bÞ
where X is a prey for specialist predator Y with Holling type II functional response. a1 is the rate of self-reproduction for the prey. The parameter a2 measures how fast the predator Y will die when there is no prey to capture, kill and eat. b1 measures the intensity of competition among individuals of species X for space, food etc. and D measures the efficiency of the prey in evading a predator’s attack. It depends on the protection afforded by the environment to the prey. D1 has similar meaning as that of D. A model given by Holling and Tanner (HT) [18] describes the dynamics of a generalist predator and its prey: dZ Z w3 UZ ; ð2aÞ ¼ AZ 1 dt K ðZ þ D3 Þ dU w4 U 2 ¼ cU ; dt Z
ð2bÞ
where Z is the most favorite food for the generalist predator U . In this model, prey and predator both grow logistically. A and K are respectively the rate of self-reproduction and carrying capacity for the prey Z. c is the growth rate of the generalist predator due to sexual reproduction. The last term in Eq. (2b) describes how loss in species U depends on per capita availability of its prey (Z). The other parameters have their usual meaning. Eqs. (2a) and (2b) together represent this model. In what follows, we attempt to arrive at a graphical representation, which gives us an idea about the parameter regimes displaying distinct dynamical possibilities. A predator–prey system qualifies as a Kolmogorov (K) system if it can be cast into the following form: dX ¼ XF ðX ; Y Þ; dt dY ¼ YGðX ; Y Þ; dt where F and G are continuous and analytic functions in the domain X P 0, Y P 0. It can be easily seen that both the models qualify as a K-system. Kolmogorov analysis of RM model yields following conditions: w1 > a2 ;
ð3aÞ
a1 D1 a2 : > b1 ðw1 a2 Þ
ð3bÞ
The local stability analysis puts the additional constraint D1 a2 2b1 þ b1 D a1 > 0: w1 a2
ð4Þ
One obtains a stable limit cycle in the phase space, if inequality (4) is violated and those from Kolmogorov analysis [9,10] are honored while making the choices of the parameter values. The two conditions can be combined to give a straight line with slope ð2D1 =ðw1 a2 ÞÞ. This line is located at a distance D from the origin of the coordinate system. a1 =b1 and a2 are the ordinate and abscissa. The straight line divides the space between two regions: the region above presents parameter values, which correspond to stable limit cycle solutions. The region below the dividing line and the straight line (a1 =b1 ¼ D) is populated by stable equilibrium solutions (cf. Fig. 1a).
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
1197
Fig. 1. (a) Two stability regions are separated by a straight line with slope 2D1 =ðw1 a2 Þ and meets the abscissa at D [15]. (b) Graphical representation of criterion (5). The figure is drawn for A=c ¼ w3 =w4 . Similar stability boundaries exist for other relationships. This figure has been adapted from May [9].
Local stability analysis of HT model [9] provides a criterion c 2ða RÞ > ; A 1þaþbþR
ð5Þ
where a¼
w3 c ; w4 A
b¼
D3 K
and
R ¼ ½ð1 a bÞ2 þ 4b1=2 :
ð6Þ
The criterion when graphically represented, provides an idea about what parameter values will lead to a limit cycle. A graphical representation of the criterion for the case Ac ¼ ww34 is given in Fig. 1b. The unstable region signifies stable limit cycle solutions, which emanate from stable ones through super critical Hopf bifurcations. The assumed relationship is biologically realistic and was chosen for the shake of simplicity. Similar stability boundaries exist for other relationships. It should be noted that a Kolmogorov analysis of this system does not put any constraint on the parameter values of the system (Eqs. (2a) and (2b)). When parameters are chosen in such a way that both RM and HT systems display regular persistent periodic behavior, these systems act as oscillators. Let us call RM oscillator as LCO (1) and HT oscillator as LCO (2). When suitably coupled, these systems will force each to generate chaos. An ecologically sound coupling of these systems gives us the following model: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ dY w1 YX w2 YU ; ¼ a2 Y þ dt ðX þ D1 Þ ðY þ D2 Þ dZ Z w3 UZ ; ¼ AZ 1 dt K ðZ þ D3 Þ dU w4 U 2 ; ¼ cU ðY þ ZÞ dt
ð7aÞ ð7bÞ ð7cÞ ð7dÞ
where parameter a1 is the rate of self-reproduction of species X , b1 measures the intensity of competition among individuals of species X for space, food, etc. w=ðX þ DÞ is the per capita rate of removal of species X by Y . Y is a specialist predator, i.e., X is the only food for it. Therefore, Y dies out exponentially in the absence of X . w is the maximum value that function wX =ðX þ DÞ can attain. a2 is the rate at which Y dies out exponentially in the absence of its prey X . w1 X =ðX þ D1 Þ denotes gain in the specialist predator population due to proportionate loss in its prey. A and K are respectively the rate of self-reproduction and carrying capacity for prey Z. The last term in Eq. (7b) represents the
1198
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
functional response of the predator U , which is a generalist predator. It switches its prey whenever its favorite food option Z is in short supply. The first term in the evolution equation for U describes the growth of U due to selfreproduction. The last term in Eq. (7d) describes how loss in species U depends on per capita availability of its preys (Z and Y ). This gives us a model which was studied by us in Upadhyay et al. [11] with the functional response of the vertebrate predator. w2 serves as the coupling parameter. Let us imagine a situation that arises when one of the prey species (Z) leaves out the patch for a better habitat. The resulting system is described by the following equations: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ
ð8aÞ
dY w1 YX w2 YZ ; ¼ a2 Y þ dt ðX þ D1 Þ ðY þ D2 Þ
ð8bÞ
dZ w3 Z 2 ¼ cZ : dt ðY Þ
ð8cÞ
It should be noted that generalist predator is denoted by Z, instead of U . This is for the shake of the continuity of the notations for different species in the model. The generalist predator Z (in Eq. (8c)) is a sexually reproducing species. For most sexual species over most population densities, reproduction is determined primarily by female numbers. Because of the female-biased demography of the sexual species the population growth rate is linear in population size, which is well below the carrying capacity. However, the growth of a sexually reproducing population is proportional to the square of the number of individuals present, in situations when the population is at very low densities. This is known as the Allee effect [12]. In such species, the decline of the population to low densities occurs due to a variety of causes [13]. Most important of these is inbreeding depression, which results from physical or chemical modification of the environment of organisms by social interaction or by density-dependent mating success. For example, 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. The Allee effect variant of this model can be written by replacing Z by Z 2 in the first term of Eq. (8c) and adding an additional constant (D3 ) in the denominator of the last term of the same equation. In this case, the last equation of the model system is modified to: dZ w3 Z 2 ¼ cZ 2 : dt ðY þ D3 Þ
ð9Þ
Some insect top-predators very often switch to alternative prey in situations when their favorite food is in short supply. Eqs. (8a)–(8c) define the linear phase of model 1. The non-linear phase is described by Eqs. (8a), (8b) and (9). These equations describe model 1. We chose to study the non-linear phase as the linear phase does not support the chaotic behavior at all. The sexually reproducing populations are covered by this phase when they are under Allee effect. A dynamical representation of these model systems are given in Fig. 2. Limit cycle oscillator of the first kind (LCO (1)) is obtained when parameters of Rosenzweig–McArthur model (RM; first subsystem) [8] are set at limit cycle values derived from an application of pseudo–prey method [6]. Likewise, the
prey
Specialist predator
pseudo-prey
Generalist predator
prey
Specialist predator
pseudo-prey
Generalist predator
(a) L. C. O.(1)
L.C.O.(2)
(b) Fig. 2. (a) A representation of the biological system. (b) Two coupled limit cycle oscillators.
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204 L. C. O.(1)
1199
L.C.O.(1)
Fig. 3. Dynamical representation of model 2.
other subsystem (obtained by cutting the link between the specialist and generalist predators in the second diagram of Fig. 2a when set at appropriate parameter values serves as a limit cycle oscillator of the second kind. The predator–prey interaction in RM model is formulated in Volterra scheme, which assumes exponential decay of a predator population in the absence of its sole prey. In contrast to this, the second subsystem is a valid representation of predator–prey interaction when the predator has other food options in addition to its most favorite prey. When we couple two LCO (1)’s, we get a model proposed by Rosenzweig and MacArthur [8] and studied by Hastings and Powell [19] and by Rinaldi et al. [3]. The model is described by the following equations: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ
ð10aÞ
dY w1 YX w2 YZ ; ¼ a2 Y þ dt ðX þ D1 Þ ðY þ D2 Þ
ð10bÞ
dZ w3 YZ ¼ cZ þ : dt ðY þ D3 Þ
ð10cÞ
We call this model 2. The same can be dynamically represented as done in Fig. 3. 3. The methodology The two models dynamically differ in one essential way: model 2 represents two coupled oscillators of type 1 whereas model 1 results from a coupling of LCO (1) and LCO (2). Model parameters are selected in accordance with a method given in Upadhyay et al. [6]. For 2D scans, two parameters are varied on both sides of their limit cycle values. An important part of the methodology of selection of parameter values for simulation experiments is the following theorem: Two coupled K-systems in the oscillatory mode would yield either cyclic (stable limit cycles and quasi-periodicity) or chaotic solutions depending on the strength of coupling between the two. If key parameters of a model system are known, one can study the dynamics by varying these system parameters on both sides of their basic value for which both subsystems yield a stable limit cycle in the phase space. Rest of the parameters are kept constant. For example, a1 and c are two key parameters for model 1. The ranges for variation of these parameters in simulation experiments are decided on the available biological information [16].
4. Results and discussion As far as observability of chaos in nature is concerned, it can be divided into two distinct categories: (i) robust chaos, and (ii) short-term recurrent chaos. We give precise definitions for both of them and present results, which suggest that robust chaos should be rarely found in nature. 4.1. Robust chaos If chaotic dynamics exists in a region in 2D parameter space spanned by two crucial parameters of a model system and also if the basin of attraction for the associated attractor is large; one understands that the chaotic dynamics displayed by the system is robust. If this is the case, chaos would be a dominant dynamical mode of the corresponding real system. 4.2. Short-term recurrent chaos There exist two mechanisms which generate this kind of chaotic behavior. If chaotic dynamics is displayed in a region in 2D parameter space of a model system, it can still support short-term recurrent chaos (strc) provided the
1200
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
0.046 0.044 0.042 0.040 0.038 c 0.036 0.034 0.032 0.030 0.028 1.92
1.93
1.94
1.95
1.96
1.97
1.98
a1 Fig. 4. Distribution of points, where chaos was observed in model 1 with invertebrate top-predator. a1 was varied in the range ½0:5; 2:5. The parameter c was varied in the range ½0:005; 0:05. a2 was set at 1.
basins of attraction of different attractors are intermixed. Alternatively, it can be generated when chaos exists at discrete and isolated points in 2D parameter space. As chaos (exponential sensitivity of system’s trajectories on ‘‘initial conditions’’) takes time (inversely proportional to the max. Lyapunov exponent of the system) to develop, system gets needed amount of time to lock itself on to a non-chaotic attractor before it settles to a fully developed chaotic state. In the former case the agency which generates strc is exogenous stochastic influences. It is purely deterministic changes in the crucial parameters of the system, which cause this behavior in the latter. In the present paper, we explore chaotic dynamics of two model systems, which differ from each other in one essential way; i.e., the top-predator in model 1 is a generalist predator whereas that of model 2 is a specialist one. In these model systems, the top-predator is assumed to be an invertebrate. Fig. 4 present results of simulation experiments on model 1. It evident that chaotic dynamics exist at discrete and isolated points in these systems. The secular changes in model parameters cause frequent interruptions in chaotic dynamics, at which time, the system evolves on non-chaotic attractors (periodic and fixed point). Chaotic dynamics in model 1 was observed only when a2 was fixed at 1. For other values of this parameter, no chaos was found in the model. This suggests that there does not exist a window in parameter space for chaotic dynamics. Chaotic dynamics does not exist for other values of a2 . Fig. 4 suggests that there does not exist a window in parameter space for chaotic dynamics. We present basin boundary calculations for chaotic attractors supported by model 1 at bottom-left corner ða1 ¼ 1:92; c ¼ 0:03Þ, top-right corner ða1 ¼ 1:975; c ¼ 0:045Þ and middle ða1 ¼ 1:95; c ¼ 0:035Þ of chaotic region in Figs. 5–7, respectively. All the three figures are the XY view of the basin boundary structure (100 6 X 6 100, 100 6 Y 6 100) of chaotic attractor (shown in yellow colour 1). 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 [20] for all the basin boundary calculations. The basin boundary calculations are presented as colour diagrams (for detail description see [11,21]). It is clear from Figs. 5–7 that the chaotic attractor is the dominant one. No other attractor with appreciable basin size exist. The encroachment into the basin of chaotic attractor by basin of attractor at infinity (shown in green colour here) can be observed. The basin boundary diagram also shows that there is no intermixing of basins of different attractors. This establishes that strc in this model is generated by deterministic changes in the crucial parameters of the system. Fig. 8 present results of 2D parameter scan on model 2. It can be seen from this figure that chaos exists in a region (non-zero area) in a1 w3 parameter space. But the chaotic region is intruded by limit cycle attractors. We show results of basin boundary calculations at the top-left corner ða1 ¼ 1:75; w3 ¼ 3:75Þ, top-right corner ða1 ¼ 3:0; w3 ¼ 3:75Þ and in the middle ða1 ¼ 2:5; w3 ¼ 3:0Þ of the chaotic region in Figs. 9–11, respectively. All these figures reveal that the basins of different attractors are intermixed. In such a case, the dynamics of the system is decided by external stochastic
1
For interpretation of color in Figs. 5–11, the reader is referred to the web version of this article.
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
1201
Fig. 5. Model 1. Basin boundary calculation for bottom-left corner point (a1 ¼ 1:92; c ¼ 0:03) in Fig. 4 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:55, D2 ¼ 10, w3 ¼ 1, D3 ¼ 20.
Fig. 6. Model 1. Basin boundary calculation for top-right corner point (a1 ¼ 1:975; c ¼ 0:045) in Fig. 4 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:55, D2 ¼ 10, w3 ¼ 1, D3 ¼ 20.
influences and the result is strc; a form of chaotic behavior, wherein it is interrupted by other kind of dynamics at irregular and unpredictable intervals. Therefore, we see that short-term recurrent chaos can be caused either by deterministic changes in the system parameters (model 1) or by exogenous stochastic influences (model 2). In the former, chaotic dynamics is interrupted by smooth changes in system’s parameters whereas stochastic influences affect such interruptions in the latter one.
5. Conclusions Thus, it is clear that robust chaos is less likely to be found in nature. Instead, what nature seems to abound in has recently been called strc [14,15]. The strc can occur through two mechanisms. The first mechanism involves deterministic changes in the system parameters. The other invokes exclusive role of abrupt changes in initial conditions. The only systems, wherein robust chaos has been observed so far, are diffusively coupled predator–prey systems [17]. Diffusion couples stable limit cycle oscillators on a spatial gradient. These coupled oscillators force each other at ‘‘incommensurate’’ frequencies to generate chaos. The chaotic dynamics exists in a region of non-zero measure in 2D
1202
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
Fig. 7. Model 1. Basin boundary calculation for middle point (a1 ¼ 1:95; c ¼ 0:035) in Fig. 4 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:55, D2 ¼ 10, w3 ¼ 1, D3 ¼ 20.
4.5
SF
3.75
LC Chaos
w3
3 2.25 1.5 0.75 0 0
0.5
1
1.5
2
2.5
3
3.5
a1
Fig. 8. Model 2. 2D scan diagram in (a1 ; w3 ) parameter space. The parametric values for limit cycle solution are (using pseudo-prey method) chosen as a1 ¼ 2, b1 ¼ 0:05, w ¼ 1, D ¼ 10 ¼ D1 ¼ D2 , a2 ¼ 1, w1 ¼ 2, w2 ¼ 1:5c ¼ 0:7, w3 ¼ 1, D3 ¼ 20. The parameter space for the simulation experiment is chosen as 0:5 6 a1 6 3, 0:1 6 w3 6 4.
Fig. 9. Model 2. Basin boundary calculation for top-left corner point (a1 ¼ 1:75; w3 ¼ 3:75) in Fig. 8 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:5, D2 ¼ 10, c ¼ 0:7, D3 ¼ 20.
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
1203
Fig. 10. Model 2. Basin boundary calculation for top-right corner point (a1 ¼ 3:0; w3 ¼ 3:75) in Fig. 8 for the chaotic attractor. The values of other parameters are as in Fig. 9.
Fig. 11. Model 2. Basin boundary calculation for middle point of chaotic region (a1 ¼ 2:25; w3 ¼ 3:0) in Fig. 8 for the chaotic attractor. The values of other parameters are as in Fig. 9.
parameter space and there is no other competing attractor in the initial condition space. Thus, there is no intermixing of the attractor basins. In our earlier work [6], we had suggested that the biology of the top-predator would be a crucial factor for the determination of food chain dynamics. The results of the present study indicate that chains ending at specialist predators are dictated by exogenous stochastic influences. On the contrary, those with generalist top-predators are governed by deterministic changes in system parameters. We opine that the natural systems with first kind of food chains would present difficult challenges as far as program of quantification of their dynamical complexity is concerned. The other kind of systems seems to allow such a program to be implemented smoothly. This conjecture is to be tested in the laboratory and in the field.
Acknowledgements The research described in this paper is supported by Department of Science and Technology, New Delhi, grant under SERC Fast Track scheme for young scientists 2001–2002, to the second author (RKU).
1204
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
References [1] May RM. Simple mathematical models with very complicated dynamics. Nature 1977;261:459–67. [2] Gilpin ME. Spiral chaos in a predator–prey model. Am Nat 1979;113:306–8. [3] Rinaldi S, Casagrandi R, Gragnani A. Reduced order models for the prediction of the time of occurrence of extreme episodes. Chaos, Solitons & Fractals 2001;12:313–20. [4] Ellner S, Turchin P. Chaos in a noisy world: new methods and evidence from time-series analysis. Am Nat 1995;145:343–75. [5] (a) Turchin P, Hanski I. An empirically based model for latitudinal gradient in vole population dynamics. Am Nat 1997;149(5):842–74; (b) Smith RH, Daniels S, Simkiss K, Bell ED, Ellner SP, Forest MB. Blowflies as a case study in non-linear population dynamics. In: Perry JN, Smith RH, Woiwod IP, Morse D, editors. Chaos in real data. Dordrecht, Netherlands: Kluwer Academic Publishers; 2000. p. 137–72. [6] Upadhyay RK, Rai V. Why chaos is rarely observed in natural populations? Chaos, Solitons & Fractals 1997;8:1933–9. [7] Upadhyay RK, Iyengar SRK, Rai V. Chaos: an ecological reality? Int J Bifurcat Chaos 1998;8:1325–33. [8] Rosenzweig ML, MacArthur RH. Graphical representation and stability conditions of predator–prey interactions. Am Nat 1963;97:205–23. [9] May RM. Stability and complexity in model ecosystems. Princeton, New Jersey: Princeton University Press; 1974. [10] Yodzis P. Introduction to theoretical ecology. New York, USA: Harper and Row; 1989. [11] Upadhyay RK, Rai V, Iyengar SRK. How do ecosystems respond to external perturbations? Chaos, Solitons & Fractals 2000;11:1963–82. [12] Allee WC, Emerson AE, Park O, Schmidt T. Principles of animal ecology. London: Saunders; 1949. [13] Lande R. Genetics and demography in biological conservation. Science 1988;241:1455–60. [14] Turchin P, Ellner SP. Living on the edge of chaos: population dynamics of Fennoscandian voles. Ecology 2000;81:3099–116. [15] Rai V. Chaos in natural populations: edge or wedge? Ecological Complexity, 2004, in press. [16] Jorgensen SE, editor. Handbook of environmental data and ecological parameters. Pergamon Press; 1979. p. 142–216. [17] Rai V, Jayaraman G. Is Diffusion induced chaos robust? Curr Sci 2003;84(7):925–9. [18] Pielou EC. Mathematical ecology. New York: John Wiley and Sons; 1977. [19] Hastings A, Powell T. Chaos in three species food chain. Ecology 1991;72:896–903. [20] Nusse HE, Yorke JA. Dynamics: Numerical exploration. New York: Springer-Verlag; 1994. [21] Upadhyay RK. Multiple attractors and crisis route to chaos in a model food-chain. Chaos, Solitons & Fractals 2003;16:737–47.
Chaotic population dynamics and biology of the top-predator Vikas Rai a, Ranjit Kumar Upadhyay a
b,*
Department of Applied Mathematics, Delhi College of Engineering, Bawana Road, Delhi 110 042, India b Department of Applied Mathematics, Indian School of Mines, Dhanbad, Jharkhand 826 004, India Accepted 8 December 2003
Abstract We study how the dynamics of a food chain depends on the biology of the top-predator. We consider two model food chains with specialist and generalist top-predators. Both types of food chains display same type of chaotic behavior, short-term recurrent chaos; but the generating mechanisms are drastically different. Food chains with specialist top-predators are dictated by exogenous stochastic factors. On the contrary, the dynamics of those with the generalist top-predator is governed by deterministic changes in system parameters. The study also suggests that robust chaos would be a rarity. Ó 2004 Elsevier Ltd. All rights reserved.
1. Introduction Ecology is a parent discipline of chaos. Since the seminal work of Sir Robert May [1] model ecosystems have been designed and studied [2,3]. These studies conveyed that chaos might be a dominant mode of ecosystem evolution. As such data requirements to capture chaos in the wild are stringent. One needs a string of precise measurements of species abundance for a length of time, which cannot be fixed a priori. On contrary to this, the truth of the real world is that ecological time-series are short and noisy. This severely limits the possibility of chaos getting detected even if it is there. A couple of years ago, it occurred to us that the origin of the crisis might be intrinsic to organization of ecological interactions. Since then, we have made sustained efforts to understand why it is so difficult to observe chaos in nature. Earlier authors had suggested [4,5] that the narrowness of the boundary between chaos and noise is responsible for the crisis. All natural systems are noisy. Chaotic dynamics in appearance is noise-like. The noisy character of these systems are caused either by intrinsic factors or extrinsic forces. When a system is subjected to external influences (deterministic or stochastic) complex dynamical phenomena (chaos and noisy chaos) emerge. These phenomena are often indistinguishable. This leads to failure of attempts to observe chaos. The work [6,7] of the present authors has opened up a new possibility that the reason for chaos getting poor evidential support might lie elsewhere. Moreover, it also suggested that the biology of the top-predator has an important role to play as far as the dynamical complexity of a food chain is concerned. In the present paper, we take up this issue in detail. The paper is organized as follows. Section 2 presents the model systems that have been chosen for detailed dynamical investigations. The methodology of the present study is presented in Section 3. We discuss the main findings in Section 4 and record the messages of the study in Section 5.
*
Corresponding author. Tel.: +91-326-2200817; fax: +91-326-2210028. E-mail address: [email protected] (R.K. Upadhyay).
0960-0779/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2003.12.065
1196
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
2. Model systems Among four basic interactions (predator–prey, competition, interference and mutualism), predator–prey interaction is the most common. There are two types of predators: (i) specialists and (ii) generalists. A specialist predator is the one, which dies out exponentially fast when its favorite food is absent or is in short supply. The latter predators switch over to alternative food options when its most preferred food is in short supply. The Rosenzweig–MacArthur (RM) model [8] is the one, which describes the dynamics of a specialist predator and its prey: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ
ð1aÞ
dY w1 YX ¼ a2 Y þ ; dt ðX þ D1 Þ
ð1bÞ
where X is a prey for specialist predator Y with Holling type II functional response. a1 is the rate of self-reproduction for the prey. The parameter a2 measures how fast the predator Y will die when there is no prey to capture, kill and eat. b1 measures the intensity of competition among individuals of species X for space, food etc. and D measures the efficiency of the prey in evading a predator’s attack. It depends on the protection afforded by the environment to the prey. D1 has similar meaning as that of D. A model given by Holling and Tanner (HT) [18] describes the dynamics of a generalist predator and its prey: dZ Z w3 UZ ; ð2aÞ ¼ AZ 1 dt K ðZ þ D3 Þ dU w4 U 2 ¼ cU ; dt Z
ð2bÞ
where Z is the most favorite food for the generalist predator U . In this model, prey and predator both grow logistically. A and K are respectively the rate of self-reproduction and carrying capacity for the prey Z. c is the growth rate of the generalist predator due to sexual reproduction. The last term in Eq. (2b) describes how loss in species U depends on per capita availability of its prey (Z). The other parameters have their usual meaning. Eqs. (2a) and (2b) together represent this model. In what follows, we attempt to arrive at a graphical representation, which gives us an idea about the parameter regimes displaying distinct dynamical possibilities. A predator–prey system qualifies as a Kolmogorov (K) system if it can be cast into the following form: dX ¼ XF ðX ; Y Þ; dt dY ¼ YGðX ; Y Þ; dt where F and G are continuous and analytic functions in the domain X P 0, Y P 0. It can be easily seen that both the models qualify as a K-system. Kolmogorov analysis of RM model yields following conditions: w1 > a2 ;
ð3aÞ
a1 D1 a2 : > b1 ðw1 a2 Þ
ð3bÞ
The local stability analysis puts the additional constraint D1 a2 2b1 þ b1 D a1 > 0: w1 a2
ð4Þ
One obtains a stable limit cycle in the phase space, if inequality (4) is violated and those from Kolmogorov analysis [9,10] are honored while making the choices of the parameter values. The two conditions can be combined to give a straight line with slope ð2D1 =ðw1 a2 ÞÞ. This line is located at a distance D from the origin of the coordinate system. a1 =b1 and a2 are the ordinate and abscissa. The straight line divides the space between two regions: the region above presents parameter values, which correspond to stable limit cycle solutions. The region below the dividing line and the straight line (a1 =b1 ¼ D) is populated by stable equilibrium solutions (cf. Fig. 1a).
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
1197
Fig. 1. (a) Two stability regions are separated by a straight line with slope 2D1 =ðw1 a2 Þ and meets the abscissa at D [15]. (b) Graphical representation of criterion (5). The figure is drawn for A=c ¼ w3 =w4 . Similar stability boundaries exist for other relationships. This figure has been adapted from May [9].
Local stability analysis of HT model [9] provides a criterion c 2ða RÞ > ; A 1þaþbþR
ð5Þ
where a¼
w3 c ; w4 A
b¼
D3 K
and
R ¼ ½ð1 a bÞ2 þ 4b1=2 :
ð6Þ
The criterion when graphically represented, provides an idea about what parameter values will lead to a limit cycle. A graphical representation of the criterion for the case Ac ¼ ww34 is given in Fig. 1b. The unstable region signifies stable limit cycle solutions, which emanate from stable ones through super critical Hopf bifurcations. The assumed relationship is biologically realistic and was chosen for the shake of simplicity. Similar stability boundaries exist for other relationships. It should be noted that a Kolmogorov analysis of this system does not put any constraint on the parameter values of the system (Eqs. (2a) and (2b)). When parameters are chosen in such a way that both RM and HT systems display regular persistent periodic behavior, these systems act as oscillators. Let us call RM oscillator as LCO (1) and HT oscillator as LCO (2). When suitably coupled, these systems will force each to generate chaos. An ecologically sound coupling of these systems gives us the following model: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ dY w1 YX w2 YU ; ¼ a2 Y þ dt ðX þ D1 Þ ðY þ D2 Þ dZ Z w3 UZ ; ¼ AZ 1 dt K ðZ þ D3 Þ dU w4 U 2 ; ¼ cU ðY þ ZÞ dt
ð7aÞ ð7bÞ ð7cÞ ð7dÞ
where parameter a1 is the rate of self-reproduction of species X , b1 measures the intensity of competition among individuals of species X for space, food, etc. w=ðX þ DÞ is the per capita rate of removal of species X by Y . Y is a specialist predator, i.e., X is the only food for it. Therefore, Y dies out exponentially in the absence of X . w is the maximum value that function wX =ðX þ DÞ can attain. a2 is the rate at which Y dies out exponentially in the absence of its prey X . w1 X =ðX þ D1 Þ denotes gain in the specialist predator population due to proportionate loss in its prey. A and K are respectively the rate of self-reproduction and carrying capacity for prey Z. The last term in Eq. (7b) represents the
1198
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
functional response of the predator U , which is a generalist predator. It switches its prey whenever its favorite food option Z is in short supply. The first term in the evolution equation for U describes the growth of U due to selfreproduction. The last term in Eq. (7d) describes how loss in species U depends on per capita availability of its preys (Z and Y ). This gives us a model which was studied by us in Upadhyay et al. [11] with the functional response of the vertebrate predator. w2 serves as the coupling parameter. Let us imagine a situation that arises when one of the prey species (Z) leaves out the patch for a better habitat. The resulting system is described by the following equations: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ
ð8aÞ
dY w1 YX w2 YZ ; ¼ a2 Y þ dt ðX þ D1 Þ ðY þ D2 Þ
ð8bÞ
dZ w3 Z 2 ¼ cZ : dt ðY Þ
ð8cÞ
It should be noted that generalist predator is denoted by Z, instead of U . This is for the shake of the continuity of the notations for different species in the model. The generalist predator Z (in Eq. (8c)) is a sexually reproducing species. For most sexual species over most population densities, reproduction is determined primarily by female numbers. Because of the female-biased demography of the sexual species the population growth rate is linear in population size, which is well below the carrying capacity. However, the growth of a sexually reproducing population is proportional to the square of the number of individuals present, in situations when the population is at very low densities. This is known as the Allee effect [12]. In such species, the decline of the population to low densities occurs due to a variety of causes [13]. Most important of these is inbreeding depression, which results from physical or chemical modification of the environment of organisms by social interaction or by density-dependent mating success. For example, 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. The Allee effect variant of this model can be written by replacing Z by Z 2 in the first term of Eq. (8c) and adding an additional constant (D3 ) in the denominator of the last term of the same equation. In this case, the last equation of the model system is modified to: dZ w3 Z 2 ¼ cZ 2 : dt ðY þ D3 Þ
ð9Þ
Some insect top-predators very often switch to alternative prey in situations when their favorite food is in short supply. Eqs. (8a)–(8c) define the linear phase of model 1. The non-linear phase is described by Eqs. (8a), (8b) and (9). These equations describe model 1. We chose to study the non-linear phase as the linear phase does not support the chaotic behavior at all. The sexually reproducing populations are covered by this phase when they are under Allee effect. A dynamical representation of these model systems are given in Fig. 2. Limit cycle oscillator of the first kind (LCO (1)) is obtained when parameters of Rosenzweig–McArthur model (RM; first subsystem) [8] are set at limit cycle values derived from an application of pseudo–prey method [6]. Likewise, the
prey
Specialist predator
pseudo-prey
Generalist predator
prey
Specialist predator
pseudo-prey
Generalist predator
(a) L. C. O.(1)
L.C.O.(2)
(b) Fig. 2. (a) A representation of the biological system. (b) Two coupled limit cycle oscillators.
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204 L. C. O.(1)
1199
L.C.O.(1)
Fig. 3. Dynamical representation of model 2.
other subsystem (obtained by cutting the link between the specialist and generalist predators in the second diagram of Fig. 2a when set at appropriate parameter values serves as a limit cycle oscillator of the second kind. The predator–prey interaction in RM model is formulated in Volterra scheme, which assumes exponential decay of a predator population in the absence of its sole prey. In contrast to this, the second subsystem is a valid representation of predator–prey interaction when the predator has other food options in addition to its most favorite prey. When we couple two LCO (1)’s, we get a model proposed by Rosenzweig and MacArthur [8] and studied by Hastings and Powell [19] and by Rinaldi et al. [3]. The model is described by the following equations: dX wYX ¼ a1 X b1 X 2 ; dt ðX þ DÞ
ð10aÞ
dY w1 YX w2 YZ ; ¼ a2 Y þ dt ðX þ D1 Þ ðY þ D2 Þ
ð10bÞ
dZ w3 YZ ¼ cZ þ : dt ðY þ D3 Þ
ð10cÞ
We call this model 2. The same can be dynamically represented as done in Fig. 3. 3. The methodology The two models dynamically differ in one essential way: model 2 represents two coupled oscillators of type 1 whereas model 1 results from a coupling of LCO (1) and LCO (2). Model parameters are selected in accordance with a method given in Upadhyay et al. [6]. For 2D scans, two parameters are varied on both sides of their limit cycle values. An important part of the methodology of selection of parameter values for simulation experiments is the following theorem: Two coupled K-systems in the oscillatory mode would yield either cyclic (stable limit cycles and quasi-periodicity) or chaotic solutions depending on the strength of coupling between the two. If key parameters of a model system are known, one can study the dynamics by varying these system parameters on both sides of their basic value for which both subsystems yield a stable limit cycle in the phase space. Rest of the parameters are kept constant. For example, a1 and c are two key parameters for model 1. The ranges for variation of these parameters in simulation experiments are decided on the available biological information [16].
4. Results and discussion As far as observability of chaos in nature is concerned, it can be divided into two distinct categories: (i) robust chaos, and (ii) short-term recurrent chaos. We give precise definitions for both of them and present results, which suggest that robust chaos should be rarely found in nature. 4.1. Robust chaos If chaotic dynamics exists in a region in 2D parameter space spanned by two crucial parameters of a model system and also if the basin of attraction for the associated attractor is large; one understands that the chaotic dynamics displayed by the system is robust. If this is the case, chaos would be a dominant dynamical mode of the corresponding real system. 4.2. Short-term recurrent chaos There exist two mechanisms which generate this kind of chaotic behavior. If chaotic dynamics is displayed in a region in 2D parameter space of a model system, it can still support short-term recurrent chaos (strc) provided the
1200
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
0.046 0.044 0.042 0.040 0.038 c 0.036 0.034 0.032 0.030 0.028 1.92
1.93
1.94
1.95
1.96
1.97
1.98
a1 Fig. 4. Distribution of points, where chaos was observed in model 1 with invertebrate top-predator. a1 was varied in the range ½0:5; 2:5. The parameter c was varied in the range ½0:005; 0:05. a2 was set at 1.
basins of attraction of different attractors are intermixed. Alternatively, it can be generated when chaos exists at discrete and isolated points in 2D parameter space. As chaos (exponential sensitivity of system’s trajectories on ‘‘initial conditions’’) takes time (inversely proportional to the max. Lyapunov exponent of the system) to develop, system gets needed amount of time to lock itself on to a non-chaotic attractor before it settles to a fully developed chaotic state. In the former case the agency which generates strc is exogenous stochastic influences. It is purely deterministic changes in the crucial parameters of the system, which cause this behavior in the latter. In the present paper, we explore chaotic dynamics of two model systems, which differ from each other in one essential way; i.e., the top-predator in model 1 is a generalist predator whereas that of model 2 is a specialist one. In these model systems, the top-predator is assumed to be an invertebrate. Fig. 4 present results of simulation experiments on model 1. It evident that chaotic dynamics exist at discrete and isolated points in these systems. The secular changes in model parameters cause frequent interruptions in chaotic dynamics, at which time, the system evolves on non-chaotic attractors (periodic and fixed point). Chaotic dynamics in model 1 was observed only when a2 was fixed at 1. For other values of this parameter, no chaos was found in the model. This suggests that there does not exist a window in parameter space for chaotic dynamics. Chaotic dynamics does not exist for other values of a2 . Fig. 4 suggests that there does not exist a window in parameter space for chaotic dynamics. We present basin boundary calculations for chaotic attractors supported by model 1 at bottom-left corner ða1 ¼ 1:92; c ¼ 0:03Þ, top-right corner ða1 ¼ 1:975; c ¼ 0:045Þ and middle ða1 ¼ 1:95; c ¼ 0:035Þ of chaotic region in Figs. 5–7, respectively. All the three figures are the XY view of the basin boundary structure (100 6 X 6 100, 100 6 Y 6 100) of chaotic attractor (shown in yellow colour 1). 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 [20] for all the basin boundary calculations. The basin boundary calculations are presented as colour diagrams (for detail description see [11,21]). It is clear from Figs. 5–7 that the chaotic attractor is the dominant one. No other attractor with appreciable basin size exist. The encroachment into the basin of chaotic attractor by basin of attractor at infinity (shown in green colour here) can be observed. The basin boundary diagram also shows that there is no intermixing of basins of different attractors. This establishes that strc in this model is generated by deterministic changes in the crucial parameters of the system. Fig. 8 present results of 2D parameter scan on model 2. It can be seen from this figure that chaos exists in a region (non-zero area) in a1 w3 parameter space. But the chaotic region is intruded by limit cycle attractors. We show results of basin boundary calculations at the top-left corner ða1 ¼ 1:75; w3 ¼ 3:75Þ, top-right corner ða1 ¼ 3:0; w3 ¼ 3:75Þ and in the middle ða1 ¼ 2:5; w3 ¼ 3:0Þ of the chaotic region in Figs. 9–11, respectively. All these figures reveal that the basins of different attractors are intermixed. In such a case, the dynamics of the system is decided by external stochastic
1
For interpretation of color in Figs. 5–11, the reader is referred to the web version of this article.
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
1201
Fig. 5. Model 1. Basin boundary calculation for bottom-left corner point (a1 ¼ 1:92; c ¼ 0:03) in Fig. 4 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:55, D2 ¼ 10, w3 ¼ 1, D3 ¼ 20.
Fig. 6. Model 1. Basin boundary calculation for top-right corner point (a1 ¼ 1:975; c ¼ 0:045) in Fig. 4 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:55, D2 ¼ 10, w3 ¼ 1, D3 ¼ 20.
influences and the result is strc; a form of chaotic behavior, wherein it is interrupted by other kind of dynamics at irregular and unpredictable intervals. Therefore, we see that short-term recurrent chaos can be caused either by deterministic changes in the system parameters (model 1) or by exogenous stochastic influences (model 2). In the former, chaotic dynamics is interrupted by smooth changes in system’s parameters whereas stochastic influences affect such interruptions in the latter one.
5. Conclusions Thus, it is clear that robust chaos is less likely to be found in nature. Instead, what nature seems to abound in has recently been called strc [14,15]. The strc can occur through two mechanisms. The first mechanism involves deterministic changes in the system parameters. The other invokes exclusive role of abrupt changes in initial conditions. The only systems, wherein robust chaos has been observed so far, are diffusively coupled predator–prey systems [17]. Diffusion couples stable limit cycle oscillators on a spatial gradient. These coupled oscillators force each other at ‘‘incommensurate’’ frequencies to generate chaos. The chaotic dynamics exists in a region of non-zero measure in 2D
1202
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
Fig. 7. Model 1. Basin boundary calculation for middle point (a1 ¼ 1:95; c ¼ 0:035) in Fig. 4 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:55, D2 ¼ 10, w3 ¼ 1, D3 ¼ 20.
4.5
SF
3.75
LC Chaos
w3
3 2.25 1.5 0.75 0 0
0.5
1
1.5
2
2.5
3
3.5
a1
Fig. 8. Model 2. 2D scan diagram in (a1 ; w3 ) parameter space. The parametric values for limit cycle solution are (using pseudo-prey method) chosen as a1 ¼ 2, b1 ¼ 0:05, w ¼ 1, D ¼ 10 ¼ D1 ¼ D2 , a2 ¼ 1, w1 ¼ 2, w2 ¼ 1:5c ¼ 0:7, w3 ¼ 1, D3 ¼ 20. The parameter space for the simulation experiment is chosen as 0:5 6 a1 6 3, 0:1 6 w3 6 4.
Fig. 9. Model 2. Basin boundary calculation for top-left corner point (a1 ¼ 1:75; w3 ¼ 3:75) in Fig. 8 for the chaotic attractor. The values of other parameters are a2 ¼ 1:0, w ¼ 1, D ¼ 10, b1 ¼ 0:05, w1 ¼ 2, D1 ¼ 10, w2 ¼ 1:5, D2 ¼ 10, c ¼ 0:7, D3 ¼ 20.
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
1203
Fig. 10. Model 2. Basin boundary calculation for top-right corner point (a1 ¼ 3:0; w3 ¼ 3:75) in Fig. 8 for the chaotic attractor. The values of other parameters are as in Fig. 9.
Fig. 11. Model 2. Basin boundary calculation for middle point of chaotic region (a1 ¼ 2:25; w3 ¼ 3:0) in Fig. 8 for the chaotic attractor. The values of other parameters are as in Fig. 9.
parameter space and there is no other competing attractor in the initial condition space. Thus, there is no intermixing of the attractor basins. In our earlier work [6], we had suggested that the biology of the top-predator would be a crucial factor for the determination of food chain dynamics. The results of the present study indicate that chains ending at specialist predators are dictated by exogenous stochastic influences. On the contrary, those with generalist top-predators are governed by deterministic changes in system parameters. We opine that the natural systems with first kind of food chains would present difficult challenges as far as program of quantification of their dynamical complexity is concerned. The other kind of systems seems to allow such a program to be implemented smoothly. This conjecture is to be tested in the laboratory and in the field.
Acknowledgements The research described in this paper is supported by Department of Science and Technology, New Delhi, grant under SERC Fast Track scheme for young scientists 2001–2002, to the second author (RKU).
1204
V. Rai, R.K. Upadhyay / Chaos, Solitons and Fractals 21 (2004) 1195–1204
References [1] May RM. Simple mathematical models with very complicated dynamics. Nature 1977;261:459–67. [2] Gilpin ME. Spiral chaos in a predator–prey model. Am Nat 1979;113:306–8. [3] Rinaldi S, Casagrandi R, Gragnani A. Reduced order models for the prediction of the time of occurrence of extreme episodes. Chaos, Solitons & Fractals 2001;12:313–20. [4] Ellner S, Turchin P. Chaos in a noisy world: new methods and evidence from time-series analysis. Am Nat 1995;145:343–75. [5] (a) Turchin P, Hanski I. An empirically based model for latitudinal gradient in vole population dynamics. Am Nat 1997;149(5):842–74; (b) Smith RH, Daniels S, Simkiss K, Bell ED, Ellner SP, Forest MB. Blowflies as a case study in non-linear population dynamics. In: Perry JN, Smith RH, Woiwod IP, Morse D, editors. Chaos in real data. Dordrecht, Netherlands: Kluwer Academic Publishers; 2000. p. 137–72. [6] Upadhyay RK, Rai V. Why chaos is rarely observed in natural populations? Chaos, Solitons & Fractals 1997;8:1933–9. [7] Upadhyay RK, Iyengar SRK, Rai V. Chaos: an ecological reality? Int J Bifurcat Chaos 1998;8:1325–33. [8] Rosenzweig ML, MacArthur RH. Graphical representation and stability conditions of predator–prey interactions. Am Nat 1963;97:205–23. [9] May RM. Stability and complexity in model ecosystems. Princeton, New Jersey: Princeton University Press; 1974. [10] Yodzis P. Introduction to theoretical ecology. New York, USA: Harper and Row; 1989. [11] Upadhyay RK, Rai V, Iyengar SRK. How do ecosystems respond to external perturbations? Chaos, Solitons & Fractals 2000;11:1963–82. [12] Allee WC, Emerson AE, Park O, Schmidt T. Principles of animal ecology. London: Saunders; 1949. [13] Lande R. Genetics and demography in biological conservation. Science 1988;241:1455–60. [14] Turchin P, Ellner SP. Living on the edge of chaos: population dynamics of Fennoscandian voles. Ecology 2000;81:3099–116. [15] Rai V. Chaos in natural populations: edge or wedge? Ecological Complexity, 2004, in press. [16] Jorgensen SE, editor. Handbook of environmental data and ecological parameters. Pergamon Press; 1979. p. 142–216. [17] Rai V, Jayaraman G. Is Diffusion induced chaos robust? Curr Sci 2003;84(7):925–9. [18] Pielou EC. Mathematical ecology. New York: John Wiley and Sons; 1977. [19] Hastings A, Powell T. Chaos in three species food chain. Ecology 1991;72:896–903. [20] Nusse HE, Yorke JA. Dynamics: Numerical exploration. New York: Springer-Verlag; 1994. [21] Upadhyay RK. Multiple attractors and crisis route to chaos in a model food-chain. Chaos, Solitons & Fractals 2003;16:737–47.