Automated Techniques for the Qualitative Analysis of Ecological Models: Continuous Models

Lynn van Coller¹

• Abstract
• Introduction
• Ratio-Dependent Model
• Resilience
• Conclusion
• Acknowledgments
• Responses
• Literature Cited

The mathematics required for a detailed analysis of the behavior of a model can be formidable. In this paper, I demonstrate how various computer packages can aid qualitative analyses by implementing techniques from dynamical systems theory. Because computer software is used to obtain the results, the techniques can be used by nonmathematicians as well as mathematicians. In-depth analyses of complicated models that were previously very difficult to study can now be done. Because the paper is intended as an introduction to applying the techniques to ecological models, I have included an appendix describing some of the ideas and terminology. A second appendix shows how the techniques can be applied to a fairly simple predator-prey model and establishes the reliability of the computer software. The main body of the paper discusses a ratio-dependent model. The new techniques highlight some limitations of isocline analyses in this three-dimensional setting and show that the model is structurally unstable. Another appendix describes a larger model of a sheep-pasture-hyrax-lynx system. Dynamical systems techniques are compared with a traditional sensitivity analysis and are found to give more information. As a result, an incomplete relationship in the model is highlighted. I also discuss the resilience of these models to both parameter and population perturbations.

KEY WORDS: bifurcation; computer analysis; dynamical systems; ecological models; qualitative analysis; ratio-dependent model; XPPAUT.

Caution: If your browser does not handle graphics, the mathematical models and all Greek characters will be unreadable in this online version. Please email the author to obtain a complete hard copy.

A fundamental aim of ecological models is to help us better understand ecological systems. If this aim is to be fulfilled, then we need to know the range of possible behavior that a model can exhibit, and which relationships or mechanisms give rise to this behavior. Although numerous mathematical techniques exist for analyzing systems of model equations, in practice, the application of these techniques to anything other than the very simplest models is often impossible, or at least sufficiently difficult that it is seldom done. The development of dynamical systems computer software has changed this situation. It is now possible for both mathematicians and nonmathematicians to study the behavior of complicated models and to obtain useful, practical information. This is illustrated in the remainder of the paper.

Most parameter values in nature are not known with certainty (Hadamard 1952, Frank 1978). It is also well known that different parameter values in a model can give rise to very different dynamics, such as stable equilibrium behavior, oscillatory behavior, or a decline to extinction. However, finding parameter combinations that give rise to these different dynamics is often a matter of trial and error and involves numerous simulations. The dynamical systems software provides an easier and more systematic way to do this.

In this paper, I substantiate this claim by using the software XPPAUT (Ermentrout 1994), which contains an interface to AUTO86 (Doedel 1981), to analyze three continuous ecological models of varying complexity. Appendix 5 contains a list of other packages that have capabilities similar to XPPAUT or AUTO86. In the main body of this paper, I consider only one of these models in detail. This model is a recent example from the controversial area of ratio-dependent modeling. It has three state variables that represent a plant, a herbivore, and a predator. In addition to highlighting some of the limitations of an isocline analysis in this three-dimensional setting, the analysis shows that the model is structurally unstable. Even small perturbations to the ratio-dependent terms alter the dynamics.

Appendix 4 discusses the application of the software to a system dynamics model that has 10 state variables and a large number of parameters. Analytical work done by hand and isocline analyses are of little use in such situations. Traditionally, computers have been used to obtain numerical solutions corresponding to a fixed parameter set and to implement sensitivity analyses. In Appendix 4 , I show how dynamical systems techniques can be used to increase our understanding of the relationships between different components in the model. In particular, bifurcation diagrams give more information than sensitivity analyses. These diagrams also highlight an incomplete relationship in the model and lead to an improvement in the formulation of the equations. Thus, the packages can help us to formulate more plausible models.

Because this paper is intended as an introduction to the application of dynamical systems techniques to ecological models, I include a number of appendices. Appendix 2 introduces some of the ideas and terminology and also provides a glossary of the terms used in the paper and appendices. Those readers unfamiliar with dynamical systems ideas may benefit from a quick reading of this appendix before continuing. A more comprehensive exposition, also specifically written for ecologists, can be found in Yodzis (1989). Another excellent introduction to nonlinear dynamics can be found in Strogatz (1994). Appendix 3 shows how the dynamical systems techniques can be applied to a fairly simple predator-prey model using XPPAUT. I show how I reproduced and improved on Bazykin's analytical results (Bazykin 1974) for the model. This appendix establishes some confidence in the reliability of the computer software.

I also discuss the topic of model resilience, a concept introduced by Holling (1973), to parameter and population perturbations.

Gutierrez et al. (1994) developed and analyzed a tritrophic model of a plant-herbivore-predator system. Currently, there is widespread debate over the validity of such ratio-dependent models (Appendix 1 summarizes the arguments). A numerical study of how such a model behaves under extreme conditions would test its validity. To do this, I introduce a small modification to the ratio-dependent model of Gutierrez et al. (1994) and study this modified model in conjunction with the original one. My analysis shows that the original model is structurally unstable, because a small perturbation to the ratio-dependent terms substantially alters the dynamics.

I also claim that the technique used by Gutierrez et al. (1994) to study the model, namely isocline analysis (for technique and examples, see Edelstein-Keshet 1988), is unsuitable in this three-dimensional setting. This is discussed in detail in section 2.2.4. Although isocline analyses have been employed in many settings and with considerable success (Fitzhugh 1961, Rosenzweig and MacArthur 1963, Gilpin 1974, Hethcote 1976, Ludwig et al. 1978, Fairen and Velarde 1979, Murray 1981), in more complicated, higher dimensional models for which the categorization of variables as slow vs. fast is not possible, their application is limited. (If the state variables in a model vary on different time scales, it is often possible to approximate the system by a two-dimensional model representing either the slow or the fast dynamics. An isocline analysis can then be done using the reduced system.) An isocline analysis allows two variables, at most, to vary simultaneously. This means that, for the current model, one variable is held fixed, resulting in a partly static representation of the dynamics. However, the dynamical systems techniques allow all three state variables to vary simultaneously, thus permitting a more accurate analysis.

2.1 Model Equations

The basic model equations and a description of the parameters can be found in Appendix 1 . There are three equations describing the dynamics of a plant, a herbivore, and a predator. In Appendix 1 , I describe how I nondimensionalized the equations and introduced a small modification to the ratio-dependent terms. The resulting equations are as follows:

In the next section, I compare results from the original model, (a_i = 0) with those from the modified equations (a_i = 0.001). I chose this latter value for the a_i's because it only affects the isocline configurations (see Appendix 1 ) at low values of the state variables, which is where the difficulties with the ratio-dependent model are encountered. I also investigate a few other values.

2.2 Model analysis

2.2.1 One-parameter studies

2.2.2 Analysis of the predator parameters

Bifurcation diagrams showing the effects of varying for both the original and the modified model are shown in Fig. 1 , with M₁ plotted on the y-axis. These diagrams were obtained using AUTO through XPPAUT. In both Fig. 1 diagrams, the position of the equilibrium point does not vary with , supporting the previous observation that does not affect the isocline configuration. However, does affect the stability of the system. In both cases, the equilibrium point is unstable for very low values of (low predator growth rate) and the stable attractor is a limit cycle for these values. The original model gives an example of hard loss of stability (see Appendix 2) so that, for certain values of , there are two stable attractors, a sink and a stable limit cycle (as in Fig. 11 ). The initial values of the state variables determine which final state is reached. Also, perturbations to the system may cause a rapid change from one stable attractor to the other if the disturbance is sufficiently large. For a very small range of values near the Hopf bifurcation point HB in Fig. 1 a, there are two stable limit cycles. The range of values is so small, however, that it is not of much biological significance.

FIGURE 1:     One-parameter bifurcation diagrams showing the effects of varying for (a) the original model and (b) the modified model. The state variable is plotted on the y-axes.

---

FIGURE 11:     Isoclines in the (a) and (b) planes for , , , , , , , , and . Both the stable equilibrium point and the limit cycle are shown.

---

, the availability (or nutritional value) of the herbivore to the predator, also directly affects the predator. Bifurcation diagrams for the original and modified models, respectively, are shown in Fig. 2 and Fig. 3 . From these figures, it can be seen that as increases, there is a general increase in the M₁ equilibrium value, or limit cycle maximum. The larger , the greater the availability of the herbivore to the predator and the easier it is for the predator to control the herbivore. Obviously, the lower the herbivore population, the higher the plant equilibrium. As the M₁ equilibrium value approaches the M₁ carrying capacity in the original model, a Hopf bifurcation (see Appendix 2 ) occurs at = 0.09 (see Fig. 2 ). The periodic orbit associated with this Hopf bifurcation undergoes a number of period-doubling bifurcations (see Appendix 2 ), which leads to more complicated cycling behavior. An example of the temporal dynamics when is shown in Fig. 4 . The minima of the M₂ and M₃ cycles, in particular, get very small (on the order of 10^-15, according to XPPAUT's data window). From a practical viewpoint, these populations would be considered extinct due to statistical variation, in which case the plant population would increase to its carrying capacity for these values of . This is exactly the case for the modified model (see Fig. 3 ). Thus, the upper Hopf bifurcation may be an artifact of ratio-dependent models. This will be discussed in more detail later.

FIGURE 2:     One-parameter bifurcation diagrams showing the effects of varying on (a) , (b) and (c) for the original model.

---

FIGURE 3:     One-parameter bifurcation diagrams showing the effects of varying on (a) , (b) and (c) for the modified model.

---

FIGURE 4:     Time plots showing the behavior of (a) , (b) and (c) when in the original model. Two complete cycles are shown in the figures.

---

has a significant effect on the equilibrium values of the state variables. This is in agreement with Gutierrez et al. (1994), but the way in which they arrive at this conclusion is not entirely correct. Gutierrez et al. (1994) state that a less efficient parasitoid has a wider C-shaped M₂-isocline. It is true that if is decreased, the M₂-isocline widens. However, this is provided that M₃ is constant. If the system is integrated and M₃ is allowed to vary until a new equilibrium is reached, and the isoclines are plotted with this new equilibrium M₃ value, then the final M₂-isocline may, in fact, have a narrower C shape than before.

This analysis has shown that the properties of the predator affect the stability of the system as well as the equilibrium magnitudes of the herbivore (directly) and the plant (indirectly). The extent of these effects depends on the properties of both the plant and the herbivore. In the next section, a parameter affecting the herbivore is examined in more detail.

2.2.3 Analysis of a herbivore parameter

Varying , the assimilation or conversion efficiency of the herbivore when plants are abundant, gives the bifurcation diagrams in Fig. 5 . Comparing Fig. 5 b with Fig. 3 a reveals that they are almost mirror images. If the mean of the cycle maxima and minima in Fig. 5 a and Fig. 2 a is taken, then these two figures are also almost mirror images. That is, decreasing has a very similar effect to increasing . Both parameters can be thought of as affecting the resistance of the herbivore to the predator. Decreasing causes a decline in the condition of the herbivore, because it cannot convert food as effectively. As a result, the detrimental effect of the predator on the herbivore is greater. Increasing , the availability (nutritional value) of the herbivore to the predator, achieves the same result but more directly.

FIGURE 5:     One-parameter bifurcation diagrams showing the effects of varying for (a) the original model and (b) the modified model. The state variable is plotted on the y-axes.

---

FIGURE 6:     Two-parameter bifurcation diagrams showing the effects of varying both and on the positions of the Hopf bifurcations in

Figure 5. (a) and (b) correspond to the original and modified models respectively. The meaning of "less complex cycles" is made clear in Figure 8.

FIGURE 7:     Two-parameter bifurcation diagram of the Hopf bifurcation continuations in -space for the modified model with .

---

FIGURE 8:     Time plots corresponding to the points marked with *'s in figures 6(a), 6(b) and 7. (a) These plots were obtained using the original model with (i) , and (ii) , . (b) These plots were obtained using the modified model with and , and (ii) , . (c) These plots were obtained using the modified model with and (i) , and (ii) , .

---

2.2.4 The role played by the isocline configurations

FIGURE 9:     Examples of isocline configurations showing the different possibilities for the position of the tritrophic equilibrium.

---

In general, it is the proximity of the tritrophic equilibrium point to the peaks of the M₁ and M₂ isoclines in the M₁M₂ and M₂M₃ planes, respectively, that is important for determining the robustness of model behavior with respect to parameter perturbations. If the equilibrium point is close to one of these peaks, then a small parameter perturbation may change the qualitative structure of the isoclines and, hence, the dynamics. However, to obtain this information, the exact equilibrium isocline configuration for a given set of parameter values needs to be known.

By inference, these criticisms have all noted that if the exact positions of the isoclines and the tritrophic equilibrium were known in both phase planes, then we could obtain a fair amount of information from them. Using XPPAUT, this is possible. In particular, the effects of introducing nonzero values for the a_i's can be studied. Figure 10 shows the results obtained using the reference parameter set for model 1 (Eq. A1.1) of Appendix 1 with a_i = 0, a_i = 0.001, and a_i = 0.005 (i = 1, 2, 3).

FIGURE 10:     Isocline configurations for (a) , (b) and (c) together with the reference parameter set.

---

2.3 Discussion

The preceding analysis has shown how XPPAUT and AUTO can be used to study the effects of various parameters on the behavior of the ratio-dependent model of Gutierrez et al. (1994). Regions of stable equilibria, sustained oscillations, and multiple equilibria were located in this study. The limitations of an isocline analysis in this three-dimensional setting were also discussed, thus supporting the need for new techniques such as the dynamical systems ones.

The ease with which models may be studied using XPPAUT allowed the simultaneous study of the original model and a modified one. I concluded that the original model is structurally unstable, because a small perturbation to the ratio-dependent terms has a significant effect on the dynamics. This supports the argument that ratio-dependent models exhibit pathological behavior and are not valid near the axes (see the background information for ratio-dependent models in Appendix 1 ). Thus, they cannot be used to study extinction or situations in which one of the state variables attains low values. However, the model by Gutierrez et al. (1994) that has been analyzed is a biological control model whose aim is to suggest what kind of predator can keep herbivore numbers low.

This paper aims to show the usefulness of dynamical systems computer packages and to indicate the ease with which they may be used to gain insights into the behavior of a model and, thus, concentrates on these points. However, the arguments should not be taken out of context. My claim is not that the packages provide a comprehensive way in which to study ecological models. Instead, they should be seen as aids to be used in conjunction with other methods of analysis. Also, these packages do have their limitations. They can only be used to study systems of ordinary differential equations where all the functions are continuous (although simple step functions can often be approximated by continuous functions, as is done in Appendix 4 for the sheep-hyrax-lynx model). Other types of models require different methods of analysis. The packages also require an initial set of parameter values to be chosen. Thus, other techniques are required for theoretical studies that aim to find the exact parameter relationships giving rise to different types of behavior. However, numerical studies can be useful for indicating the relationships that do exist, and thus can provide direction for theoretical studies.

The concept of resilience of ecological models was introduced by Holling (1973), and is described by May (1981) as the magnitude of population perturbations a system will tolerate before collapsing into some qualitatively different dynamical regime. A slightly broader interpretation (Beddington et al. 1976), identifies resilience with the magnitude of the perturbation in any specific quantity that some property of a system can undergo before suffering qualitative changes.

Wollkind et al. (1988) use the latter definition to examine the resilience of their model system for predator-prey mite interactions to changes in food quality of prey for conversion into predator births. This resilience depends on the range of parameter values that give rise to a certain type of qualitative behavior. The wider this range, the greater the systemØs resilience to perturbations of the relevant parameter, provided the parameter is not too near the range endpoints.

Collings and Wollkind (1990) discuss the resilience of various stable phenomena in their model, which exhibits metastability. In this case, the resilience is in terms of the magnitude of population (rather than parameter) perturbations that the system can withstand. If a stable phenomenon has a small domain of attraction (see Appendix 2), then small changes to the populations can lead to dramatic changes in system behavior. Such a system is not resilient.

Both of these studies used bifurcation diagrams to obtain results. Thus, bifurcation diagrams can characterize the resilience of a model system to both parameter and population perturbations. In the sheep-hyrax-lynx model in Appendix 4, the stable equilibrium behavior predominates for wide ranges of the parameter values that were considered. Because there was only a single equilibrium point, the system is resilient to both population perturbations and perturbations in the sheep female culling normal (S_FCN), provided S_FCN is within the desirable range of values and the initial population values are reasonable.

The ratio-dependent model is less resilient, in that parameter perturbations can push the system from stable equilibrium behavior into oscillatory behavior and vice versa. However, Fig. 6 shows that the regions of stable and oscillatory behavior in space are fairly large. Hence, for parameter combinations that are well within these regions, the system is resilient to small perturbations in and .

This paper has shown how techniques from dynamical systems theory can help us to understand the behavior of ecological models. Previously, trial and error had to be used to locate parameter values that gave rise to different types of behavior, but we now have a more systematic approach. In addition, bifurcation diagrams provide a concise way of summarizing results.

The results obtained were compared with those from traditional approaches. AUTO improved on some approximations that Bazykin (1974) made in his analytical study of a predator-prey model. In the sheep-hyrax-lynx model (Swart and Hearne 1989), the bifurcation diagrams gave more information than a traditional sensitivity analysis. This highlighted a relationship that had been omitted from the model. Thus, the dynamical systems techniques can be helpful in constructing better models.

For the ratio-dependent model (Gutierrez et al. 1994), XPPAUT allowed two variations of the model to be studied simultaneously. The resulting conclusion was that the original model is structurally unstable, due to its ratio-dependent terms. The dynamical systems techniques also highlighted the limitations of an isocline analysis in a three-dimensional setting.

In addition, the bifurcation diagrams can be used to characterize the resilience of models to both parameter and population perturbations.

The most important point of this paper is that dynamical systems software was used to arrive at these conclusions. A complete understanding of the underlying mathematical techniques is not required to obtain useful information about the behavior of a model. The software also allows complicated models to be analyzed in greater depth than was previously possible. It is also possible to study discrete models, the topic of a second paper.

None of the models in this paper has a seasonal component. Seasonality plays an important role in the dynamics of many natural systems; thus, studies of models that include such components would be of interest. A few such studies using dynamical systems techniques have been done (e.g., Rinaldi et al. 1993, Gragnani and Rinaldi 1995), but the bifurcation structure is more complicated than that discussed in this paper and requires a continuation package that can handle higher codimension bifurcations, such as LOCBIF (Khibnik et al. 1993; see Appendix 5). Because this paper is intended to introduce the application of dynamical systems techniques to ecological models, I have not included one of these more complex examples.

Responses to this article are invited. If accepted for publication, your response will be hyperlinked to the article. To submit a comment, follow this link. To read comments already accepted, follow this link.

I would like to thank my supervisor, Don Ludwig, for all his help and guidance throughout the research, and for suggesting avenues of analysis that would be of most interest to ecologists. I would also like to thank A. Gutierrez, John Hearne, and Johan Swart for the use of their models and for responding to all my questions, and Eusebius Doedel and Bard Ermentrout for making the various packages available and for responding to numerous questions. I am also grateful to Colin Clark, Simon Levin, and Wayne Nagata for reading various versions of the manuscript, and to the reviewers, who gave many helpful suggestions. I woould like to acknowledge Emma Smith, Oppenheimer, and I.W. Killiam Trusts for their financial support, which made possible my studies at the University of British Columbia.

Abrams, P.A. 1994. The fallacies of "ratio-dependent'' predation. Ecology 75(6): 1842-1850.

Akçakaya, H.R., R. Arditi, and L.R. Ginzburg. 1995. Ratio-dependent predation: an abstraction that works. Ecology 76(3): 995-1004.

Allgower, E.L., and K. Georg. 1992. Continuation and path following. Acta Numerica 1-64.

Arditi, R., and H. Saïah. 1992. Empirical evidence of the role of heterogeneity in ratio-dependent consumption. Ecology 73(5):1544-1551.

Arnold, V.I. 1983. Geometrical methods in the theory of ordinary differential equations. Springer-Verlag, New York, New York, USA.

Aronson, D.G., M.A. Chory, G.R. Hall, and R.P. McGehee. 1982. Bifurcations from an invariant circle for two-parameter families of maps of the plane: a computer-assisted study. Communications in Mathematical Physics 83: 303-354.

Back, A., J. Guckenheimer, M. Myers, F. Wicklin, and P. Worfolk. 1992. dstool: computer-assisted exploration of dynamical systems. Notices of the American Mathematical Society 39: 303-309.

Bazykin, A.D. 1974. Volterra's system and the Michaelis-Menten equation. Pages 103-142 in V.A. Ratner, editor. Problems in mathematical genetics. U.S.S.R. Academy of Science, Novosibirsk, U.S.S.R.

Beddington, J.R., C.A. Free, and J.H. Lawton. 1975. Dynamic complexity in predator-prey models framed in difference equations. Nature 255: 58-60.

Beddington, J.R., C.A. Free, and J.H. Lawton. 1976. Concepts of stability and resilience in predator-prey models. Journal of Animal Ecology 45: 791-816.

Berlinski, D. 1976. On systems analysis: an essay concerning the limitations of some mathematical methods in the social, political, and biological sciences. Massachusetts Institute of Technology Press, Cambridge, Massachusetts, USA.

Berryman, A.A. 1992. The origins and evolution of predator-prey theory. Ecology 73(5): 1530-1535.

Brylinski, M. 1972. Steady-state sensitivity analysis of energy flow in a marine ecosystem Pages xxx-xxx in B.C. Patten, editor. Systems analysis and simulation in ecology, II. Academic Press, New York, New York, USA.

Collings, J.B., and D.J. Wollkind. 1990. A global analysis of a temperature-dependent model system for a mite predator-prey interaction. SIAM Journal of Applied Mathematics 50(5): 1348-1372.

Collings, J.B., D.J. Wollkind, and M.E. Moody. 1990. Outbreaks and oscillations in a temperature-dependent model for a mite predator-prey interaction. Theoretical Population Biology 38: 159-191.

Deuflhard, P., B. Fiedler, and P. Kunkel. 1987. Efficient numerical pathfollowing beyond critical points. SIAM Journal of Numerical Analysis 24: 912-927.

Doedel, E.J. 1981. AUTO: a program for the automatic bifurcation analysis of autonomous systems. Congressus Numerantium 30: 265-284.

Dowden, P.B., H.A. Jaynes, and V.M. Carolin. 1953. The role of birds in a spruce budworm outbreak in Maine. Journal of Economic Entomology 46: 307-312.

Edelstein-Keshet, L. 1988. Mathematical models in biology. Random House, New York, New York, USA.

Ermentrout, G.B. 1995. XPPAUT 1.63: the differential equations tool.

Fairen, V., and M.G. Velarde. 1979. Time-periodic oscillations in a model for respiratory process of a bacterial culture. Journal of Mathematical Biology 8: 147-157.

Fitzhugh, R. 1961. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal 1: 445-466.

Forrester, J.W.1961. Industrial dynamics. Massachusetts Institute of Technology Press, Cambridge, Massachusetts, USA.

Frank, P.M. 1978. Introduction to system sensitivity theory. Academic Press, New York, New York, USA.

Freedman, H.I., and R.M. Mathsen. 1993. Persistence in predator-prey systems with ratio-dependent predator influence. Bulletin of Mathematical Biology 55(4): 817-827.

Gerald, C.F., and P.O. Wheatley. 1989. Applied numerical analysis. Addison-Wesley, Reading, Massachusetts, USA.

Gilpin, M.E. 1974. A model of the predator-prey relationship. Theoretical Population Biology 5: 333-344.

Ginzburg, L.R., and H.R. Akçakaya. 1992. Consequences of ratio-dependent predation for steady-state properties of ecosystems. Ecology 73(5): 1536-1543.

Gleeson, S.K. 1994. Density dependence is better than ratio dependence. Ecology 75(6): 1834-1835.

Gragnani, A., and S. Rinaldi. 1995. A universal bifurcation diagram for seasonally perturbed predator-prey models. Bulletin of Mathematical Biology 57(5): 701-712.

Guckenheimer, J., and P. Holmes. 1983. Nonlinear oscillations, dynamical systems, and bifurcation of vector fields. Springer-Verlag, New York, New York, USA.

Gutierrez, A.P. 1992. Physiological basis of ratio-dependent predator-prey theory: the metabolic pool model as a paradigm. Ecology 73(5): 1552-1563.

Gutierrez, A.P., N.J. Mills, S.J. Schreiber, and C.K. Ellis. 1994. A physiologically based tritrophic perspective on bottom up-top down regulation of populations. Ecology 75: 2227-2242.

Hadamard, J. 1952. Lectures on Cauchy's problem in linear partial differential equations. Dover, New York, New York, USA.

Hale, J.K., and H. Kogak. 1989. Differential and difference equations through computer experiments. Springer-Verlag, New York, New York, USA.

Hanski, I. 1991. The functional response of predators: worries about scale. Trends in Ecology and Evolution 6(5): 141-142.

Hearne, J.W. 1987. An approach to resolving the parameter sensitivity problem in system dynamics methodology. Applied Mathematical Modelling 11: 315-318.

Hethcote, H.W. 1976. Qualitative analyses of communicable disease models. Mathematical Biosciences 28: 335-356.

Holling, C.S. 1965. The functional response of predators to prey density and its role in mimicry and population regulation. Memoirs of the Entomological Society of Canada 45: 1-60.

______ . 1973. Resilience and stability of ecological systems. Annual Review of Ecology and Systematics 4: 1-23.

Hultquist, P.F. 1988. Numerical methods for engineers and computer scientists. Benjamin Cummings, Menlo Park, California, USA.

Jeffers, J.N.R. 1978. An introduction to systems analysis with ecological applications. University Park Press, Baltimore, Maryland, USA.

Kaas-Petersen, C. 1989. PATH: user's guide. University of Leeds, Leeds, UK.

Khibnik, A., Y. Kuznetsov, V. Levitin, and E. Nikolaev. 1993. Continuation techniques and interactive software for bifurcation analysis of ordinary differential equations and iterated maps. Physica D 62: 360-371.

Kowal, N.E. 1971. A rationale for modeling dynamic ecological systems. Pages xxx-xxx in B.C. Patten, editor. Systems analysis and simulation in ecology, I. Academic Press, New York, New York, USA.

Ludwig, D., D.D. Jones, and C.S. Holling. 1978. Qualitative analysis of insect outbreak systems: the spruce budworm and forest. Journal of Animal Ecology 47: 315-332.

Lundberg, P., and J.M. Fryxell. 1995. Expected population density vs. productivity in ratio-dependent and prey-dependent models. American Naturalist 146(1): 153-161.

MacDonald, N. 1978. Time lags in biological models. Springer-Verlag, Berlin, Germany.

______ . 1989. Biological delay systems. Cambridge University Press, Cambridge, Massachusetts, USA.

May, R.M. 1973. Time delay vs. stability in population models with two and three trophic levels. Ecology 54: 315-325.

______ . 1981. Models for two interacting populations. Pages 78-104 in R.M. May, editor. Theoretical ecology: principles and applications. Blackwell, Oxford, UK.

Maynard Smith, J. 1974. Models in ecology. Cambridge University Press, Cambridge, UK.

McCarthy, M.A., L.R. Ginzburg, and H.R. Akçakaya. 1995. Predator interference across trophic chains. Ecology 76(4): 1310-1319.

Murray, J.D. 1981. A pre-pattern formation mechanism for animal coat markings. Journal of Theoretical Biology 88: 161-199.

Nusse, H.E., and J.A. Yorke. 1994. Dynamics: numerical explorations. Applied Mathematical Sciences 101. Springer-Verlag, New York, New York, USA.

Patten, B.C. 1971. A primer for ecological modeling and simulation with analog and digital computers. Pages xxx-xxx in B.C. Patten, editor. Systems analysis and simulation in ecology, I. Academic Press, New York, New York, USA.

Rheinboldt, W.C., and J.V. Burkardt. 1983a . A locally parameterized continuation process. ACM Transactions on Mathematical Software 9: 215-235.

Rheinboldt, W.C., and J.V. Burkardt. 1983b. Algorithm 596: a program for a locally parametrized continuation process. ACM Transactions on Mathematical Software 9: 236-241.

Rinaldi, S., S. Muratori, and Y. Kuznetsov. 1993. Multiple attractors, catastrophes, and chaos in seasonally perturbed predator-prey communities. Bulletin of Mathematical Biology 55(1): 15-35.

Rosenzweig, M.L., and R.H. MacArthur. 1963. Graphic representation and stability conditions of predator-prey interactions. American Naturalist 97: 209-223.

Royama, T. 1984. Population dynamics of the spruce budworm Choristoneura fumiferana. Ecological Monographs 54(4): 429-462.

Sanders, C.J., R.W. Stark, E.J. Mullins, and J. Murphy, editors. 1985. Recent advances in spruce budworm research. Proceedings of CANUSA Spruce Budworm Research Symposium, 16-20 September 1984, Bangor, Maine, Canadian Forestry Service, Ottawa, Ontario, Canada.

Sarnelle, O. 1994. Inferring process from pattern: trophic-level abundances and imbedded interactions. Ecology 75(6): 1835-1841.

Seydel, R. 1988. From equilibrium to chaos: practical bifurcation and stability analysis. Elsevier, New York, New York, USA.

______ . 1991. BIFPACK: a program package for continuation, bifurcation, and stability analysis. Version 2.3+. University of Ulm, Ulm, Germany.

Strogatz, S.H. 1994. Nonlinear dynamics and chaos. Addison-Wesley, Reading, Massachusetts, USA.

Swart, J. 1987. Sensitivity of a hyrax-lynx mathematical model to parameter uncertainty. South African Journal of Science 83: 545-547.

Swart, J., and J.W. Hearne. 1989. A mathematical model to analyze predation and competition problems in a sheep-farming region. System Dynamics Review 5: 35-50.

Tomovic, R. 1963. Sensitivity analysis of dynamic systems. Translated by D. Tornquist. McGraw-Hill, London, UK.

Waterloo Maple Software. 1990-1992. The maple roots report. Waterloo Maple Software, Waterloo, Canada.

Watt, K.E.F. 1966. The nature of systems analysis. Pages xxx-xxx in K.E.F. Watt, editor. Systems analysis in ecology. Academic Press, New York, New York, USA.

Wiggins, S. 1990. Introduction to applied nonlinear dynamical systems and chaos. Springer-Verlag, New York, New York, USA.

Williams, T., and C. Kelley. 1986. GNUPLOT: an interactive plotting program: user manual. http://www.cs.dartmouth.edu/gnuplot/gnuplot.html

Wollkind, D.J., J.B. Collings, and J.A. Logan. 1988. Metastability in a temperature-dependent model system for predator-prey mite outbreak interactions on fruit trees. Bulletin of Mathematical Biology 50(4): 379-409.

Yodzis, P. 1989. Introduction to theoretical ecology. Harper and Row, New York, New York, USA.

However, Akçakaya et al. (1995) state that these problems near the axes are only pathological in a mathematical sense; in biological terms, both species would increase initially, and then predators would consume all the prey and both species would become extinct. Because prey-dependent models cannot predict this outcome, they are pathological in a biological sense. In the next section, however, I will show that ratio-dependent models do not necessarily predict this outcome either. I will describe a ratio-dependent model that predicts oscillations of large amplitude in these "pathological" regions (see Fig. 2 ). Although these large-amplitude oscillations may be interpreted as signaling extinction from a practical viewpoint, this is not the prediction of the ratio-dependent model.

Akçakaya et al. (1995) state that:

A realistic model of prey-predator interactions should be able to predict the whole range of dynamics observed in such systems in nature. A ratio-dependent model can have stable equilibria, limit cycles, and the extinction of both species as a result of over exploitation.

In an attempt to clarify some of the arguments in this debate, I introduce a small modification to the ratio-dependent model of Gutierrez et al. (1994) and study this modified model in conjunction with the original one. The analysis to follow shows that the original model is structurally unstable, because a small perturbation to the ratio-dependent terms substantially alters the dynamics.

A1.2 Model equations

The model equations are functionally homogeneous (i.e., all three equations contain the same basic terms), because the authors argue that the same generalized functional and numerical responses must describe the search, acquisition, and conversion of all organisms as they seek to satisfy their metabolic requirements. Details of the model formulation can be found in Gutierrez et al. (1994). The final equations are:

  ,

where M₁ is the biomass of plants, M₂ is the biomass of herbivores, and M₃ is the biomass of predators. The parameters represent assimilation rates corrected for the efficiency of biomass conversion; represents the proportion of the resource that is available to its consumers (that is, its apparency); D_i is the per unit demand of the consumer for resources; r_i is the base respiration rate corrected for the efficiency of biomass conversion; and b_i is the degree of self-limitation. M₀ represents a biomass equivalent of the light energy incident in the growing space of the plants, and and (which both lie between 0 and 1) scale the potential photosynthesis rate to the realized rate.

The respiration term, r_iM_i^{1 + b_i}, requires further explanation. Respiration usually increases with population density (Gutierrez et al. 1994) and, thus, should be an increasing function of M_i. However, introducing such a functional dependence increases the model's complexity. Because the effect is usually small, the authors chose the simpler formulation r_iM_i^{1 + b_i}, where b_i has a value between 0.02 and 0.05. The disadvantage of this choice is that r_i must have rather unusual units (grams per b_i per day, where grams are the units of M_i) that depend on b_i. Thus, the term r_iM_i^{1 + b_i} has the same units as (namely, grams per day, where grams are the units of M_i). This dependence of the units on b_i is not satisfactory from a mathematical viewpoint, but since b_i is small, I chose to ignore this initially. In the next section, I discuss a small alteration to the model that takes care of the difficulty.

Gutierrez et al. (1994) use parameter values corresponding to a cassava (Manihotesculenta Crantz)-mealybug (Phenacoccuss manihoti Mat.-Ferr.)-parasitoid system in Africa. They claim that their analysis demonstrates that the parasitoid Epidinocarsis lopezi (De Santis) can control the mealybug (except on poor soils), whereas E. diversicornis (Howard) and native natural enemies cannot. However, the parameter values for the cassava-mealybug-parasitoid system given to me by the authors (only values for are reported in the paper) did not yield the isocline configurations or the behavior that they described. I then scaled the equations to get a better idea of the relative magnitudes of the parameters. The procedure involved is discussed in the next section.

A1.3 Nondimensionalization

The technique of nondimensionalizing, or scaling, is commonly used to simplify a system of equations because it has many other advantages. It illuminates which parameters are most important in determining the dynamics of the model (Edelstein-Keshet 1988) and gives insight into the relative magnitudes of the parameters required to produce biologically reasonable behavior. Also, the state variables are scaled so that they all have the same order of magnitude, say between 0 and 1. This is important when solving the equations numerically, as very different magnitudes can lead to computer round-off errors (Gerald and Wheatley 1989). In the cassava-mealybug-parasitoid system, the biomass of cassava is much larger than that of the mealybug and the parasitoid (cassava mass averages ~ 2 kg, whereas the mealybug and parasitoids have average masses of ~ 2 mg). Scaling, at least, is thus vital. (For an introduction to nondimensionalizing systems of ordinary differential equations, see Edelstein-Keshet 1988.)

Natural scalings for the state variables are given by their carrying capacities when resources are abundant. Gutierrez et al. (1994) calculate these to be

Replacing M_i by (i=1,2,3) , and t by (here K_i has the dimensions of M_i and the dimensions of t; hence, (i =1,2,3) and t* are dimensionless variables) and choosing the dimensionless combinations of parameters

,

where K₀ = M₀ gives the nondimensionalized Eqs. A1.2, where I have replaced by M_i and t* by t for convenience:

(I replaced the product by , because all three parameters have the same effect on the dynamics and, thus, do not need to be considered separately.)

The choice of dimensionless parameters is not unique. Other combinations would have led to slightly different final equations, but the choices here lend themselves to biological interpretation. For example, can be thought of as the potential per unit biomass growth rate (Gutierrez et al. 1994), or as the conversion efficiency of the consumer in converting the resource into biomass. can be thought of as the availability of the resource to the consumer, or perhaps the nutritional value of the resource. is made up of a ratio of quantities. The numerator can be thought of as the maximum amount of resource available to the consumer, and the denominator as the maximum demand of consumers for resource. More simply, gives a measure of the ratio of supply to demand. Results in Gutierrez et al. (1994) are based on the relationship between and , since is just a scaling factor. Hence, results using Eqs. A1.2 are comparable with those in Gutierrez et al. (1994).

I mentioned earlier that the parameters r_i in the original model have units that depend on b_i. This can be prevented by replacing the terms r_iM_i^{1 + b_i} with terms of the form , where T_i has the same units as M_i. T_i can be thought of as a threshold value above which self-limitation becomes noticeable. With this modification, the units of r_i are per day, and r_i can be interpreted as a respiration rate, as was originally intended. The new carrying capacities are given by

It can be shown that setting the K_i 's equal to these new values and scaling the equations as above results in system (A1.2) once again. Thus, the problems with the respiration term can be ignored in the remainder of the paper.

Having scaled the equations, I needed to choose values for the parameters. An advantage of scaling is that there are now 11 parameters instead of the original 18. For comparison with the results of Gutierrez et al. (1994), I wanted to find values that resulted in isocline configurations similar to those in their paper. XPPAUT calculates and displays isoclines in two dimensions, and parameter values can be altered interactively. Since the competition effect is very small, but difficult to quantify, I followed A. Gutierrez (personal communication) and chose B_i = 0.02 (i = 1,2,3). Values of , , , , , , , gave the isocline configurations shown in Fig. 11 . These isoclines have similar shapes to those in Gutierrez et al. (1994). A noticeable difference is that the M₃ -isocline intersects the M₂ -isocline to the left of the M₂ -peak. In fact, using the techniques in Gutierrez et al. (1994), it would have been impossible to conclude that the tritrophic equilibrium is stable for the isocline configuration shown in Fig. 11, because of the position of the intersection point in the plane. For this parameter set (which I shall call the reference set), there is also a stable limit cycle. The initial values of M₁,M₂ and M₃ determine whether the system approaches the stable equilibrium or the limit cycle.

Biological considerations suggest that these values for the are rather high and that the value for is rather low. However, in the absence of better information and since this parameter set has a nontrivial, stable equilibrium point, it is a convenient starting point.

Before beginning an analysis of the model, I will introduce a small modification to the equations. Because the problems associated with ratio-dependent models described in section A1.1 involve low population densities, it would be interesting to know what effects a small modification to the model, which prevents the denominators of the ratios from getting too close to zero, would have on the dynamics. Abrams (1994) states that modifications to ratio-dependent models cannot be made biologically realistic because the original models have no clear mechanistic derivation. However, some form of modification that prevents the ratios from tending to infinity may be useful for revealing any spurious behavior near the axes that may result from the ratio dependence.

The ratios in model (A1.2) have the form

.

The difficulties are experienced when M_i approaches zero. Adding a constant in the denominator, that is, replacing the ratio by

would alleviate the problem. Although this addition may appear difficult to justify biologically, Gutierrez (1992) used an exponent of this form in his functional response term. That model is physiologically based, as is the present one.

In model A.1.2 ], the ratio-dependent terms have the form

  ,

where k is a parameter or combination of parameters. When the exponent is small, we have

.

In order to preserve this property when using the modified term (Eq. A1.3 ), I replaced Eq. A1.4 by

,

where a_i is a small constant, say 0.001. The resulting equations are:

  .

If a_i = 0 (i = 1, 2, 3) , then the above model is equivalent to system (A1.2).

It is important to note that bifurcation diagrams summarize the behavior associated with a range of parameter values. They do not represent the dynamics corresponding to a continually varying parameter (Wiggins 1990). To read a bifurcation diagram, fix the parameter at a particular value and mentally draw a vertical line at that value. Each crossing of this line with a curve in the diagram corresponds to an equilibrium point or a periodic orbit (limit cycle). The local stability properties of a particular phenomenon are given by the type of curve: solid, dotted, or open or closed circles. For example, the phase portraits in Fig. 14b (i) and (ii) were obtained by mentally drawing vertical lines at the parameter values µ = µ₁ and µ = µ₂, respectively, in Fig. 14a.

FIGURE 14:     (a) A bifurcation diagram of a transcritical bifurcation, (b)(i) a phase portrait corresponding to and (b)(ii) a phase portrait corresponding to . The stabilities of the two equilibrium points interchange for parameter values on either side of the bifurcation point.

---

FIGURE 19:     Two-parameter bifurcation diagram of -parameter space for , , and . The regions (i), (ii) and (iii) correspond to those in

Figure 16.

---

FIGURE 16:     Bazykin's two-parameter bifurcation diagram of -parameter space.

A2.1.2 Bifurcation point

A bifurcation point is a point in parameter space at which the qualitative behavior of the system changes. A stable equilibrium may become unstable at this point, or there may be a change from a stable equilibrium to oscillatory behavior. Examples can be found in sections A2.1.6, A2.1.8, A2.1.11, and A2.1.19.

A2.1.3 Continuation branch

A solution branch, or continuation branch, is a curve of equilibrium points (or limit cycles or bifurcation points) that indicates how the position and properties of the equilibrium point (or limit cycle or bifurcation point) change as a parameter (or parameters) is altered. Together, a number of these branches make up a bifurcation diagram.

A2.1.4 Domain of attraction

Suppose the system in which we are interested has a stable equilibrium point (see section A2.1.9). Then the collection of all initial state variable values from which the system tends toward this equilibrium, as time progresses, is the domain (or basin) of attraction of the equilibrium point. The equilibrium point is called an "attractor." Any stable phenomenon, such as a stable limit cycle, also has a domain of attraction and is referred to as an attractor.

A2.1.5 Hard loss of stability

In the case of a Hopf bifurcation, this occurs when there is a rapid change from stable equilibrium behavior to stable limit cycles of large amplitude. This results from a subcritical Hopf bifurcation in which the unstable orbit reaches a limit point and turns around to stabilize. An example is shown in Fig. 12. When the parameter µ is increased beyond the Hopf bifurcation at µ^*, the system rapidly changes to limit cycles of large amplitude instead of starting off with small limit cycles that grow in size as µ increases, as shown in Fig.15 . Also, as µ is decreased,there is a jump from large-amplitude cycles to a zero-amplitude equilibrium point, but this takes place at µ₁, which is less than µ^*. For µ₁ < µ < µ^*, there are two stable attractors, an equilibrium point and a limit cycle. Examples of hard loss of stability arise in the analysis of the ratio-dependent model.

FIGURE 12:     Hard loss of stability (adapted from Seydel (1988), p.74). This phenomenon gives rise to sudden changes between stable equilibrium behavior and limit cycles of large amplitude.

A2.1.6 Hopf bifurcation

A Hopf bifurcation (also known as a Poincaré-Andronov-Hopf bifurcation; Arnold 1983, Wiggins 1990) is a bifurcation point at which an equilibrium point alters stability and a limit cycle (period orbit) is initiated. The example in Fig. 15 has a stable limit cycle surrounding an unstable equilibrium point. It is also possible to have an unstable limit cycle encircling a locally stable equilibrium point. Unstable periodic orbits are indicated by open rather than solid circles.

FIGURE 15:     Diagram showing the dynamics corresponding to a range of parameter values surrounding a Hopf bifurcation at .

A2.1.7 Limit cycle

(See periodic orbit, section A2.12)

A2.1.8 Limit point

A limit point or saddle-node bifurcation occurs when there are two equilibrium points on one side of the bifurcation point, but none on the other side. Figure 13a shows an example of a bifurcation diagram of a limit point (LP). For µ < µ^*, there are no equilibrium points at which both populations are nonzero. µ^* is thus the limiting value of µ for which equilibrium points exist; hence, the name limit point. A possible phase portrait in two dimensions for µ > µ^* is shown in Fig. 13b for the particular value µ = µ₁. A is a locally stable equilibrium point and B is a saddle point. The initial values of x₁ and x₂ determine the subsequent behavior of the system. If the initial point is in the domain of attraction of A (to the right of point B in Fig. 13b), then the system will approach A. However, if the initial point lies on the other side of B, it will be repelled away from B in the opposite direction to A. Notice how the size of the domain of attraction of A, in terms of the state variable x₁, decreases as µ decreases towards µ^* (see Fig. 13a).

FIGURE 13:     (a) Bifurcation diagram showing a limit point (LP) and (b) a phase portrait corresponding to . For there are no equilibrium points at which both populations are nonzero.

A2.1.9 Local stability

Suppose the system in which we are interested is disturbed slightly from its equilibrium point. For example, a week of warmer weather may cause an insect's growth rate to increase slightly. If, after the disturbance is removed, the system returns to its original equilibrium, then the equilibrium point is said to be locally stable and is called an "attractor." Otherwise, it is said to be unstable and is a "repeller." Locally stable equilibrium points are called sinks and locally unstable ones are called saddle points or sources.

A2.1.10 Parameter

A parameter is a quantity, e.g., fecundity rate or predation rate, that is used in describing the dynamics of a state variable. Although a state variable evolves with time, a parameter is kept constant as time progresses. In this paper, parameter values are varied across ranges of values to see how they affect the qualitative behavior of the state variables. For example, increasing the fecundity rate of a population that is at equilibrium may cause the population to start oscillating.

A2.1.11 Period-doubling bifurcation

A period-doubling bifurcation occurs when a limit cycle or periodic orbit undergoes a bifurcation and there is an exchange of stability to orbits having double the period. See Seydel (1988) for a detailed explanation.

A2.1.12 Periodic orbit

The terminology periodic orbit is used to describe a state variable (e.g., population density) that is oscillating in a regular, repetitive manner. For continuous models having dimension of at least two, such an orbit is also called a limit cycle, provided that the orbit is isolated (that is, trajectories starting at points near the orbit spiral either toward or away from the cycle; Strogatz 1994).

A2.1.13 Phase portrait

Suppose our system is continuous and has two state variables, say a prey (x₁) and a predator (x₂). We can represent the behavior over time of both populations in a single diagram called a phase portrait. Examples are shown in Figs. 13b and 14b . An illustration of the relationship between time plots and phase portraits can be found in Holling (1973: 3).

A2.1.14 Saddle point

A saddle point is an equilibrium point that attracts in certain directions and repels in others. Such an equilibrium point has unstable and stable manifolds. Initial points lying on these manifolds are repelled from, or attracted toward, the equilibrium point, respectively. Other initial points first may be attracted and then repelled (see Edelstein-Keshet 1988).

A2.1.15 Sink

A sink is a locally stable equilibrium point and may be either a stable node or a spiral attractor. In the case of a spiral attractor, the state variables oscillate with decreasing amplitude as they approach the equilibrium point.

A2.1.16 Soft loss of stability

This occurs at a Hopf bifurcation when there is a continuous change from stable equilibrium behavior to limit cycles of small amplitude. The amplitude of these cycles increases gradually for parameter values further from the Hopf bifurcation. Figure 15 gives an example of soft loss of stability. See also section A2.1.5.

A2.1.17 Source

A source is an equilibrium point that is locally unstable. Any disturbance to the system will cause the state variables to move away from this point. In the case of a spiral repeller, this repulsion may be oscillatory.

A2.1.18 State variable

Suppose we are interested in a system consisting of plants, herbivores, and predators. Then the"state" of the system can be described by the relative biomasses or densities of these populations. The variables that are used in a mathematical model of the system to represent these biomasses or densities are called state variables.

A2.1.19 Transcritical bifurcation

At a transcritical bifurcation point, two equilibrium points coincide and exchange stabilities. An example is shown in Fig. 14. There are two equilibrium points, A and B, at each value of the parameter µ. A is stable for µ < µ^* and a saddle point for µ > µ^*. The situation is reversed for B. Figures 14b(i) and (ii) show possible phase portraits in two dimensions for µ = µ₁ and µ = µ₂. Note that the bifurcation diagram in Fig. 14a only indicates the positions of the equilibrium points in terms of one of the state variables, x₁.

A2.2 Introduction to dynamical systems concepts

In this paper, I am more interested in obtaining qualitative than quantitative information about model behavior. That is, I am more interested in finding out how varying a parameter value affects the modes of behavior of a model, rather than finding the exact parameter values at which these changes occur. Berlinski (1976) notes that qualitative insights are at a much greater depth than partially quantitative results. Because of the uncertainty associated with parameter values in nature (Hadamard 1952, Frank 1978, Swart 1987), and because models are necessarily simplifications of reality, the exact quantitative predictions of a model have little significance. It is the type or mode of behavior that is most important.

Qualitative analyses are concerned mainly with long-term rather than transient behavior. Transient dynamics vary with the initial values of variables and the time period over which solutions are calculated, whereas qualitative studies deal with the eventual behavior of the system once the initial transients have died away.

There are many different modes of behavior that characterize the long-term dynamics of a model. Suppose that the modelØs state variables represent interacting populations. In the long term, any given population may become extinct, reach a stable equilibrium, or oscillate in time. Another mathematical possibility is that a population may grow indefinitely, although this is not biologically plausible except for short time periods. As one or more parameters are varied, it is possible to change from one mode of behavior to another.

Consider the following two-dimensional model due to Bazykin (1974):

Here, x represents prey density and y predator density. The term ax describes the exponential growth of the prey population in the absence of predators, and takes into account intraspecific competition among prey. In the absence of the predator, the prey population experiences logistic growth. The term - cy describes the exponential decline in the predator population in the absence of prey. The Holling type-II functions (Holling 1965) and describe the interaction between the two populations.

If intraspecific competition in the prey population is sufficiently intense (i.e., is sufficiently large), then the prey population reaches a stable equilibrium. However, the equilibrium value is too low to support the predator population and this forces the predator into extinction.

If is decreased, thus reducing competition among prey, then the predator-prey interaction becomes relatively more important and both populations reach stable equilibria.

If is decreased further, the dynamics become more complicated. Instead of reaching a stable equilibrium, the two populations oscillate in time in a sustained manner. This is known as stable limit cycle behavior. For any positive initial values, the predator and prey populations oscillate, with the amplitudes of the oscillations increasing or decreasing until the regular oscillatory pattern corresponding to the limit cycle is reached.

The exact values of at which the previously mentioned changes in qualitative behavior occur are known as bifurcation points. In particular, the transition between equilibrium behavior and limit cycle behavior is known as a Hopf bifurcation (at least in this case; limit cycles can also be generated in other ways).

Figure 15 is an example of a one-parameter bifurcation diagram. These diagrams summarize the phase portraits corresponding to different (fixed) values of a parameter, showing how the qualitative behavior of a system of equations changes as a parameter is varied across a range of values. The parameter in question is plotted on the x -axis and one of the state variables (in this case, x) is plotted on the y -axis. Solid lines represent stable equilibria, dotted lines represent unstable equilibria, and solid circles indicate the maxima and minima of stable limit cycles. In order to interpret a bifurcation diagram, imagine the phase portraits corresponding to different values of the parameter. This is done by choosing a particular value of the parameter and mentally drawing a vertical line at that value. The points at which this line crosses the curves in the bifurcation diagram will indicate the points (in terms of the quantity on the y -axis) at which equilibria or limit cycles occur, as well as the stability of these phenomena. This information can be used to construct a possible phase portrait.

From Fig. 15, we can see that for x reaches a stable equilibrium from any initial value, as indicated by the arrows. But for the equilibrium is unstable and the system is attracted instead toward the stable limit cycle. The maxima and minima of the limit cycle (with respect to x) are indicated by the solid circles. Notice that the amplitude of the limit cycle oscillations increases as decreases.

Whereas phase portraits and time plots represent the dynamics of the system for a fixed set of parameter values, bifurcation diagrams summarize the dynamics for a whole range of parameter values. They show where stable equilibria or limit cycles can be expected, and can also indicate whether multiple stable states are possible at a given parameter value, or whether extinction can be expected. Thus, bifurcation diagrams provide more concise and complete information than other graphical methods (Collings et al. 1990). Each parameter can be varied in turn to see what effect it has on the dynamics. Hence, instead of using trial and error to locate the parameter combinations giving rise to oscillations or other behavior, we now have a more systematic method.

A drawback of the dynamical systems approach is that the mathematics required to obtain a bifurcation diagram can be quite formidable, especially when large models are involved. Fortunately, in recent years, many computer packages have become available to simplify the task. I used XPPAUT (Ermentrout 1994) and DSTOOL (Back et al. 1992) to obtain the results in this paper. XPPAUT has a graphical interface for the continuation and bifurcation package AUTO86 (Doedel 1981). I shall refer to AUTO rather than XPPAUT when making use of this interface. In Appendix 3, I demonstrate how these packages can be used and illustrate the reliability of AUTO's routines by reproducing and, in fact, improving upon the results that Bazykin obtained by hand for his model (Eq. A2.1).

(see Appendix 2) are summarized in the two-parameter bifurcation diagram in Fig. 16. This figure shows three regions in -parameter space, each corresponding to a different form of qualitative behavior. In region (i), there is a spiral attractor (say A) at which both x and y are nonzero, as well as a saddle point (say B) on the x -axis to the right of A. In region (ii), A is an unstable source surrounded by a stable limit cycle and B is still a saddle point on the x -axis. In region (iii), B is now the attractor and A is a saddle point, but A has moved into the fourth quadrant.

Since A and B exchange stability in crossing from regions (i) and (ii) to region (iii), we expect a transcritical bifurcation to occur as the line with negative slope in Fig. 16 is crossed from regions (i) and (ii) into region (iii). In crossing from region (i) to region (ii), point A changes from stable to unstable and a limit cycle is initiated. Thus, we expect a curve of Hopf bifurcations to divide these regions.

I used XPPAUT to reproduce Bazykin's results. (A comprehensive, interactive tutorial on how to use XPPAUT is available on the World Wide Web at the address

ftp://mthbard.math.pitt.edu/pub/bardware/xpptut/start.html.

The tutorial can also be obtained via anonymous ftp from ftp.math.pitt.edu in the directory pub/bardware). The required driving program containing the model equations and parameter values is shown at the end of this appendix A3.1. As can be seen, it is very straightforward. In order to use AUTO to create a bifurcation diagram, we need a starting point that must be an equilibrium point. By choosing values of 0.3 for and 0.1 for , we can either determine such a point analytically (as Bazykin did), or we can use XPPAUT to perform the task numerically.

I then used AUTO to vary . AUTO locates a transcritical bifurcation (see Appendix 2) at and a Hopf bifurcation at , together with stable limit cycles (see Fig. 17). Since is fixed in Fig. 17, this one-parameter diagram describes the dynamics along a horizontal line at, say, in Fig. 16.

FIGURE 17:     One-parameter bifurcation diagram obtained by varying in Bazykin's model with , , , and . The labels A and B mark the continuation branches corresponding to the equilibrium points, HB stands for Hopf bifurcation and BP for bifurcation point (transcritical in this case).

---

These are obtained by setting the right-hand sides in the predator-prey equations (see Eqs. A2.1) equal to zero and solving for x and y. As expected, does not appear in the x -coordinate for A, but does appear in that for B.

In the range , point A is unstable and B is a saddle point (see Appendix 2). There is also a stable limit cycle surrounding A. Hence, these values of correspond to region (ii) in Fig. 16. For , A is a stable equilibrium point and B is again a saddle point. This configuration corresponds to region (i) in Fig. 16. For , B is now stable and A is a saddle point, but the numerical output from AUTO shows that the y -coordinate for A is negative for these values of . This corresponds to region (iii) in Fig. 16.

AUTO can be used to continue the Hopf bifurcation at in as well as , resulting in a two-parameter bifurcation diagram. This is shown in Fig. 18.

FIGURE 18:     Two-parameter continuation of the Hopf bifurcation shown in Figure 17.

---

Note that the curve of Hopf bifurcations is very different from Bazykin's straight line in Fig. 16. I will return to this point shortly. A second observation is that, for the given values of a, b, c, and d , a Hopf bifurcation (and hence limit cycle behavior) is only possible if and , which is a fairly small region of parameter space and, hence, probably not of much biological significance. (AUTO slows down considerably as increases toward 0.5 and tends to 0 and never actually reaches this point, although the curve does get very close if AUTO is left to run for a sufficiently long time period. It can be verified analytically that the curve does pass through (0,0.5). However, complex behavioral changes occur at this point, which is why AUTO has computational difficulties.)

It is not possible to continue transcritical bifurcations in two parameters using AUTO, as they are structurally unstable unless the system possesses intrinsic symmetry. However, if is fixed at a number of different values and one-parameter bifurcation diagrams similar to Fig. 17 are created by varying in each case, then the values corresponding to transcritical bifurcations can be recorded. An approximation to the two-parameter curve can then be drawn through these points. Figure 19 shows the resulting curve, together with the Hopf bifurcation continuation.

Bazykin's model is simple enough that the curves in Fig. 19 can be determined analytically, although the algebraic manipulations are rather tedious. It can be shown that transcritical bifurcations occur along the straight line

,

and Hopf bifurcations occur along the curve

(A description of the methods involved in calculating these curves can be found in Guckenheimer and Holmes (1983) or Wiggins (1990). I used MAPLE (Waterloo Maple Software 1990-1992) for some of the algebraic manipulations. I have only included these expressions here to show the accuracy of AUTO's results. The availability of packages such as AUTO precludes the need for the reader to perform such mathematical manipulations.) These are exactly the curves shown in Fig. 19. Hence, AUTO's results are more accurate than Bazykin's linear approximation (Fig. 16) for the Hopf bifurcation curve. The new regions (i), (ii), and (iii) are shown in Fig. 19. XPPAUT or DSTOOL can be used to obtain phase portraits and time plots corresponding to different parameter combinations in Fig. 19. Due to space considerations, I have not included any such diagrams here. Time plots are useful for indicating the speed with which the stable equilibrium or limit cycle is attained and the period of the limit cycle, if applicable. If a system takes a long time to approach an attractor, then the transient dynamics may be of greater practical importance than the long-term behavior.

This section has shown how to rederive Bazykin's work (1974) on the system of Eqs. A2. more accurately and without having to completely understand the mathematical techniques involved. It has also established some confidence in AUTO's results.

3.1 Computer listing for Bazykin model

# Model equations.

dx/dt=a*x-b*x*y/(1+alp*x)-eps*x*x

dy/dt=-c*y+d*x*y/(1+alp*x)

# Initial values.

x(0)=0.2

y(0)=0.2

# Parameters and nominal values.

param a=0.6,b=0.3,c=0.4,d=0.2,alp=0.3,eps=0.001

done

Note that AUTO encounters points beyond which it cannot calculate: S_FCN=0.54 is such a point. Using XPPAUT to solve the system of equations numerically for parameter values on either side of this limiting point, we can determine the causes of the difficulties. Output from a numerical integration is shown in XPPAUT's data window, and from this we see that the sheep population dies out for S_FCN >0.54. Beyond the limiting value, there is no equilibrium at which all three populations are present and, hence, AUTO cannot continue the equilibrium branch any further.

A second observation can be made by looking at the bifurcation diagram for pasture in Fig. 20. As the equilibrium number of sheep declines (as S_FCN increases and approaches the value corresponding to sheep extinction), the equilibrium value for pasture increases rapidly (and AUTO encounters numerical difficulties). Clearly, this is unrealistic and suggests that some modification should be made to the model equations to limit pasture growth. I return to this in section A4.5.4. Before modifying the equations to limit pasture growth, it is informative to study the roles of the existing density-dependent functions. This can be viewed as another form of sensitivity analysis, in which we investigate the system's response to whole functions instead of single parameters.

FIGURE 20:     One-parameter bifurcation diagrams obtained from varying (nominal value = 0.28/yr). In (a) revenue is plotted on the y-axis, in (b) the state variable for lynx females is used, and in (c) the state variable for pasture is used.

A4.5.3 The effects of density dependence

dP/dt = pasture production - pasture grazing = A.P_PN - T_SSU.G_N.G_M(P_A),

whereAis the area of the farming region under study; P_PN is the pasture production normal (average pasture growth rate); T_SSU is total small stock units (a representative value for the number of sheep); G_N is the grazing normal (average amount of pasture grazed per unit stock); and G_M is the grazing multiplier, which is a function of P_A, the pasture availability index. P_A is given by , where P is an average pasture density. The grazing multiplier is an increasing function of P_A, which reaches a maximum of 1.5 around P_A = 3 and has a value of 1 at P_A = 1.

To remove a density-dependent function from the model, the simplest approach is to replace it by a constant function. For example, we can set G_M equivalent to 1. This was done for each multiplier function in turn, and bifurcation diagrams obtained from the altered model were compared with those from the original model in each case.

Removing the grazing multiplier, G_M, had the greatest effect on model dynamics. Even after I altered a number of parameter values, the system did not reach equilibrium. The sheep population either increased indefinitely or decreased to extinction for each parameter set that was tried. This is not surprising, since the grazing multiplier affects pasture grazing and sheep fecundity as well as sheep juvenile deaths. Without G_M, the density dependence of pasture grazing and sheep dynamics on pasture availability ( G_M is a function of P_A) is removed. The results show that this density dependence is critical for regulating the sheep-pasture subsystem.

The effects of the other functions in the model were quantitative rather than qualitative. That is, removing them from the model tended to decrease the range of parameter values over which sheep, hyrax, and lynx coexist at equilibrium, but did not alter the qualitative dynamics. This does not imply that these other functions play an insignificant role in the model. They have a regulatory effect and reduce the impact of parameter changes on model behavior. This resilience to disturbances is very desirable (Holling 1973) and is expected of many natural systems. I observed earlier that the model lacks a feedback relationship that would limit pasture growth when sheep densities are very low. Also, pasture availability has a significant influence on pasture grazing and sheep fecundity and, hence, on the predictions of the model. In the next section, Eq. A4.1 is modified to include density-dependent growth.

A4.5.4 Adding density dependence to pasture growth

dP/dt = pasture production - pasture grazing = A.P_PN.P_M(P_A) - T_SSU.G_N.G_M(P_A),

where P_M is a function of pasture availability, P_A. For high values of P_A, we expect pasture growth to slow down and saturate. We also expect a decline in growth when P_A is very low. A function having this general form is the Ricker function (see Fig. 21).

FIGURE 21:     The pasture multiplier () as a function of pasture availability ().

---

To test the effects of this new function, I generated a number of bifurcation diagrams and compared them with those from the original model. Figure 22 shows those diagrams that correspond to Fig. 20.

In Fig. 22a, total revenue declines to zero as S_FCN increases. This is more mathematically satisfactory than Fig. 20a, where the curve stops abruptly at a positive revenue value, because it indicates clearly where a positive sheep population is no longer possible. Figure 22c shows the maximum pasture density that occurs when a positive sheep equilibrium (and hence, a positive value for revenue) is impossible. This density is lower than in Fig. 20c, but depends on the exact nature of the pasture multiplier.

FIGURE 22:     One-parameter bifurcation diagrams obtained from varying (nominal value = 0.28/yr) for the new model which includes a pasture limiting multiplier. In (a) revenue is plotted on the y-axis, in (b) the state variable for lynx females is used, and in (c) the state variable for pasture is used.

---

Another observation from Fig. 22 is that, instead of AUTO being unable to converge at very low S_FCN (sheep female culling normal) values, a Hopf bifurcation occurs, and the bifurcation diagrams show that no stable equilibrium at which all three populations coexist is possible for S_FCN < 0.059. For these values, there is insufficient pasture to support the high sheep population. Thus, introducing the pasture multiplier has improved the dynamics at low values of S_FCN as well as high values, and has solved the problem of revenue increasing rapidly as in Fig. 20a.

Using XPPAUT to integrate the system numerically gives insight into the dynamics corresponding to the different regions in the bifurcation diagrams. In particular, the temporal dynamics show that the pasture multiplier slows and limits pasture growth as desired. Comparing the dynamics of the original and modified models at S_FCN = 0.5 shows that limiting pasture growth has a stabilizing influence on the dynamics (see Fig. 23 ). The oscillatory approach to equilibrium by the original model (Fig. 23a) is replaced by a smooth approach in Fig. 23b.

FIGURE 23:     Time plots obtained using (a) the original model and (b) the modified model with .

---

The limiting value for pasture in Fig. 22c is still rather high, but I had no experimental data on which to base the form of the pasture multiplier, and did not think it worthwhile to fiddle with the function to obtain a more plausible value. Effects of introducing the multiplier have been adequately demonstrated already.

A closer look at the behavior of the model may suggest further modifications of the equations. I have shown just one example of how bifurcation diagrams can help in the process of model building. The next section describes how two-parameter studies can be used to obtain useful summaries of model behavior.

A4.5.5 A summary of the effects of culling both hyrax and lynx

Using the modified model, the bifurcation diagrams for the lynx culling normal L_CN have the form shown in Fig. 24 . For L_CN below the limit point, the sheep population dies out under too much predation by lynx. Using AUTO, the limit point can be continued in H_CN, the hyrax culling normal, as well as L_CN. That is, we can see how the position of this limit point varies as a function of both H_CN and L_CN. This gives the two-parameter bifurcation diagram in Fig. 25 . From this figure, it can be seen that sheep become extinct as a result of the combined effect of lynx predation and competition with hyrax for pasture, since both L_CN and H_CNare low in the region where sheep die out.

FIGURE 24:     One-parameter bifurcation diagram of revenue as a function of for the modified model. The change in stability (denoted by the change from a solid to a dotted line) occurs at a limit point.

Figure 24. (To determine the behavior corresponding to a particular region in this two-parameter diagram, choose values for and in this region and then use XPPAUT to integrate the system numerically.)

FIGURE 26:     Two-parameter bifurcation diagram of the and parameter space for the modified model.

---

Figure 26 shows that all three populations coexist at equilibrium for a large set of culling rates. The diagram would be of even greater use if the revenue value corresponding to each point in this two parameter space were known. This can be done by recording information given by XPPAUT and using some other graphics package to plot a three-dimensional surface.

Using the modified model developed in the previous section, I fixed the value of H_CN, chose L_CN = 0.15 (say), and used XPPAUT to find the equilibrium point numerically. I then used the AUTO interface to vary L_CN in both directions. Using the GRAB feature of XPPAUT to move along the branch of equilibrium points, I recorded the revenue values at regular intervals along the curve. I did this for a number of H_CN values and plotted the results using the public domain graphics package GNUPLOT (Williams and Kelley 1986). A surface plot and corresponding contour plot are shown in Fig. 27 . As can be seen from the figure, there is a large region of parameter space over which revenue does not vary much, indicating that the model is very robust to changes in the culling rates in this region. This is a desirable property when it comes to developing management strategies.

FIGURE 27:     (a) Surface plot and (b) contour plot of revenue as a function of the hyrax and lynx culling normals. The arrow in (b) indicates the direction of increasing revenue.

A4.6 Conclusions

The analysis of the previous sections has led to a number of insights into the sheep-pasture-hyrax-lynx system. Figure 22 shows that altering the number of ewes that are culled only affects the sheep-pasture subsystem. However, the effects on this subsystem are considerable. Culling too many ewes will obviously cause sheep numbers to decline. More important is the effect on revenue. Significant increases in revenue are possible if the farmer culls fewer ewes, because there is a greater return if these sheep are allowed to reproduce than if they are taken to market (provided that the sheep stock is not too large for the pasture to support it, i.e., provided S_FCN > 0.059.) However, there is a trade-off: culling fewer ewes results in a lower cash flow. Another trade-off results from the decrease in pasture availability that accompanies a larger sheep stock. This is already reflected in the model by the dependence of sheep fecundity on pasture availability. However, an additional quantity representing the quality of sheep may be useful, as this will affect the returns from wool and mutton sales and, hence, revenue. This presents another opportunity for improving the model.

In section A4.5.3, I showed that density dependence of sheep and pasture dynamics on pasture availability is critical for regulating the sheep-pasture subsystem. In fact, this encouraged me to modify pasture growth to include density dependence. This modification restricts both pasture and revenue values from increasing indefinitely (compare Figs. 20 and 22) and also stabilizes the temporal dynamics (see Fig. 23). Other modifications, such as including the effects of hyrax grazing in the pasture equation, can also be made. In the interests of space, I have not discussed such a modification, but the approach in preceding sections indicates how such a modification can be introduced and studied.

Finally, the effects of culling both hyrax and lynx were summarized using two-parameter diagrams. In particular, Fig. 27 shows that the model is robust to changes in culling rates, provided that these rates are sufficiently high.

A4.7 Mathematical details for the sheep-hyrax-lynx model

A4.7.1 Modelling delays in system dynamics models

To represent this in the model, I use averaged or smoothed versions of certain quantities to calculate growth rates. Three quantities are averaged in this model: hyrax density (H_D), prey abundance (A_P), and the grazing multiplier (G_M). Sheep fecundity and juvenile mortality are functions of the grazing multiplier average, _, instead of G_M. I use a first-order distributed delay (May 1973, MacDonald 1978) to calculate the averaged version, e.g.,

d /dt = ( ) /t_del,

where t_del is the average delay time. Examples of this type of delay equation are found in Forrester (1961). MacDonald (1978) provides a detailed explanation of the mathematics underlying this differential equation. This ordinary differential equation is added to the original system of model equations, thus increasing the dimension of the system to be solved. Fortunately, this increase in dimension does not pose a problem when solving the system numerically.

Further discussions of delays in biological systems and their effects on system behavior are found in May (1973) and MacDonald (1989). For the sheep-hyrax-lynx model, the delays did not affect stability, even for large average delay times.

A4.7.2 Rescaling model equations

dv_i/dt = f_i(v₁, ... , v_i, ... , v_m).

Let v_i = s_i , then

= v_i/s_i

and

d /dt = (1/s_i)(dv_i/dt)= (1/s_i)f_i(s₁ , ... ,s_i , ... , s_m ).

Dropping the bars for convenience, we get

dv_i/dt = (1/s_i)f_i(s₁v₁, ... ,s_iv_i, ... , s_mv_m),

as required. The state variables have been replaced by s_iv_i and the differential equation for v_i is divided through by s_i, as stated in section A4.4.

7. PITCON (Rheinboldt and Burkardt 1983a, 1983b).This is a Fortran subprogram for continuation of equilibrium points and for detecting limit points.