To: cestvr@ces.iisc.ac.in http://www.consecol.org/Journal/vol1/iss1/art5/ --------------7C9372155F5F Content-Type: text/html; charset=iso-8859-1 Content-Transfer-Encoding: quoted-printable Content-Base: "http://www.consecol.org/Journal/vol1/i ss1/art5/" Conservation Ecology: Automated Techniques for the Qualitative Analysis of Ecological Models - Continuous Models Copyright © 1997 by the Ecological Society of America

van Coller, L. 1997. Automated techniques for the qualitative analysis of ecological models: continuous models. Conservation Ecology [online]1(1): 5. Available from the Internet. URL: http://www.consecol.org/vol1/iss1/art5
A version of this article in which text, figures, tables, and appendices = are separate files may be found by following this link.
Research

Automated Techniques for the Qualitative Analysis of Ecological Model= s: Continuous Models

Lynn van Coller1
Department of Mathematics,University of British Columbia

ABSTRACT

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.

1. INTRODUCTION

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.

2. RATIO-DEPENDENT MODEL

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 =2E 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, (ai =3D 0) with those from the modified equations (ai =3D 0.001). I chose this latter value for the ai'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

Gutierrez et al. (1994) did a partial qualitative analysis of the original model, but used different techniques (namely, isocline analysis), so it will be informative to compare some of the results. I chose dimensionless parameters with this in mind (see Appendix 1 ). Because one of the main conclusions in Gutierrez et al. (1994) concerns the relative efficacy of two parasitoids in controlling the cassava mealybug population, I begin by studying two parameters that affect the third trophic level, namely and = . I then discuss = some parameters affecting the herbivore. Only a few examples are included here. Results, together with possible biological interpretations, are described in the next section.

2.2.2 Analysis of the predator parameters

The parameter = can be thought o= f as the potential predator biomass growth rate when prey are abundant, or as the predator's conversion efficiency in the presence of abundant prey. Note that the value of = does not affect = the isoclines. The M3 zero isocline is given by = , and the solutio= n of this equation is independent of = (see equations (1). Hence, an isocline analysis similar to that done in Gutierrez et al. (1994) would not give any insight into how this parameter influences the behavior of the system.

Bifurcation diagrams showing the effects of varying = for both the ori= ginal and the modified model are shown in = Fig. 1 , with M1 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 s= table 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) pla= nes for , , , , , , , , and . Both the stable equilibrium point and the limit cycle are shown.


Analysis of the temporal dynamics of the system (using XPPAUT) for different values of = shows that large= r values of = decrease the tim= e taken to reach equilibrium. Increasing = , the potential g= rowth rate of the predator, has a stabilizing influence on the system. This seems biologically plausible, because higher values of = suggest that the= predator is better adapted to controlling its prey. Interestingly, this trait does not affect any of the equilibrium biomasses. The modified model has a much smaller range of parameter values over which cycles occur, and the amplitudes of these cycles are smaller than for the original model (see Fig. 1 ). Thus, even though the ai' s have small values, they appear to have a stabilizing influence on the dynamics.

, 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 =2E From these figures, it can be seen that as = increases, there= is a general increase in the M1 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 M1 equilibrium value approaches the M1 carrying capacity in the original model, a Hopf bifurcation (see Appendix 2 ) occurs at =3D 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 =2E The minima of the M2 and M3 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) <= IMG ALIGN=3DBOTTOM SRC=3D"_2448_tex2html_wrap1378.gif">, (b) and (c) for the original model.


FIGURE 3:     One-parameter bifurcation diagrams showing the effects of varying on (a) <= IMG ALIGN=3DBOTTOM SRC=3D"_2448_tex2html_wrap1378.gif">, (b) and (c) for the modified model.


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


These low minima also lead to numerical difficulties, resulting from the way in which the model is formulated, particularly the dependence of many of the terms on the ratio = . These ratios be= come difficult to evaluate numerically as Mi approaches zero, causing the ratio to tend to infinity. XPPAUT fails to calculate zero isoclines, whereas AUTO often enters an infinite loop if such a situation arises and may crash. Setting a fairly low total number of steps for a continuation sometimes allows AUTO to break out of the loop and signal nonconvergence. Manually stopping a continuation when one of the state variables gets very close to zero also prevents the package from crashing. This explains why the limit cycles in Fig. 2 are only calculated up to = . Such problems d= o not occur when using the modified model.

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 M2-isocline. It is true that if = is decreased, th= e M2-isocline widens. However, this is provided that M3 is constant. If the system is integrated and M3 is allowed to vary until a new equilibrium is reached, and the isoclines are plotted with this new equilibrium M3 value, then the final M2-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 =2E 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 simil= ar 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 availabilit= y (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.


Two-parameter diagrams in (, )-space can be ge= nerated by continuing the Hopf bifurcations in Fig. 5 in two parameters. The results are shown in = Fig. 6 =2E These diagrams show clearly that decreasing = or increasing = has a similar ef= fect, since the Hopf bifurcation curves occur roughly along the diagonal that has one end at the origin. Part of this result could have been predicted from Gutierrez et al. (1994), because they note that it is the ratio of and = that determines = the nature of the herbivore isocline. This inverse relationship is thus expected. However, the one-parameter bifurcation diagrams give us the additional information that, for certain parameter ranges, limit cycles are possible.

= FIGURE 6:     Two-parameter bifurcation diagrams showing the effects of = varying both and= on the position= s of the Hopf bifurcations in Figure 5. (a) and (b) c= orrespond to the original and modified models respectively. The meaning of "less complex cycles" is made clear in Figure 8.


Some important observations can be made by comparing the Fig. 6 diagrams. First, AUTO could not calculate beyond the point denoted by MX in Fig. 6 a; this problem does not occur in Fig. 6 b. A closer investigation reveals that the equilibrium values for the state variables are close to zero in the upper left triangle of the two-parameter space, and this results in numerical problems when using the original model. Figure 6 b gives a more complete picture of the dynamics. There are two regions that correspond to stable cycles, whereas stable tritrophic equilibria occur in the other region. Low values of give rise to high equilibrium values of M1 and high values of give rise to low equilibrium values, an ecologically important distinction. A comparison of the two Fig. 6 diagrams shows that the dynamics for the modified model, although very similar to those for the original model, are more stable in general. That is, the region corresponding to tritrophic equilibria is larger and the cycles in the upper half of the two-parameter space are less complex. These claims are made clearer in = Fig. 7 and Fig. 8 =2E

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) , .


Figure 7 shows the (, )-space for the m= odified model with ai =3D0.002 (i =3D 1,2,3). The region of stable equilibria is even larger than in F= ig. 6 b, resulting in smaller regions of cycles. The presence of the ai 's seems to have a stabilizing effect on the dynamics of the system. Figure 8 shows time plots corresponding to points marked with *'s in the left-hand section of the upper region of cycles in Figs. 6 a, b, and Fig. 7 =2E The left-hand portion of this region is where the values of M2 and M3 are low and, hence, where the nonzero ai 's have most effect. Clearly, nonzero ai 's reduce the complexity of the cycles (even for very small values), and increasing their values also reduces the cycle amplitude. Additionally, the cycles for the original model (ai =3D0, i =3D 1,2,3) undergo long periods of extremely low values, which is unrealistic from an ecological viewpoint.

2.2.4 The role played by the isocline configurations =46rom their isocline analysis, Gutierrez et al. (1994) concluded that the parasitoid Epidinocarsis lopezi could control the cassava mealybug, whereas E. diversicornis could not. However, with three state variables all having similar time scales, these deductions are not as straightforward as they may seem. First, the equilibrium isocline configuration in the M1M2(M2M<= sub>3) phase plane depends on the value of M3(M1) as well as the parameter values. Thus, noting how an isocline changes as a parameter is varied does not give a complete picture. Second, it is impossible to tell from the qualitative structure of the isoclines which intersection point in the M1M2 plane corresponds to a tritrophic equilibrium. = Figure 9 shows three possibilities. Two of these (Fig. 9 b and c) appear in Gutierrez et al. (1994), who assumed, incorrectly, that the equilibrium point in Fig. 9 b was unstable.

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


Even if the exact position of the equilibrium point is known, it is impossible to tell, from the qualitative structure of the isoclines, whether this point is stable or unstable and whether or not limit cycles occur. For example, although altering the parameter = has no effect on= the isoclines, low values of = give rise to uns= table fixed points and stable limit cycles, and high values give rise to a stable equilibrium (see Fig. 1 ). Thus, numerical computation is needed to determine the exact location as well as the local stability of an equilibrium point for the current model.

In general, it is the proximity of the tritrophic equilibrium point to the peaks of the M1 and M2 isoclines in the M1M2 and M2M3 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 ai'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 ai =3D 0, ai =3D 0.001, and ai =3D= 0.005 (i =3D 1, 2, 3).

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


A comparison of Fig. 10 a and b demonstrates that introducing the ai's prevents the M1 and M2 isoclines from passing through the origin. Hence, the equilibrium values for the state variables do not approach zero as rapidly as for the original model, and the modified model is more robust to parameter variations in this region of low biomasses. Increasing the ai's from 0.001 to 0.005 reduces the humped shape of the M1 and M2 isoclines. The result is an even more robust model. In fact, when ai =3D 0.005 (i =3D 1, 2, 3), no regions of cycling behavior are encountered when the various parameters are altered. Because values of 0.005 are still small, this suggests that the model is structurally unstable and, hence, predictions from ratio-dependent models should be treated with caution.

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.

3. RESILIENCE

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=D8s 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 (SFCN), provided SFCN 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, F= ig. 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 = .=

4. CONCLUSION

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

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.

Acknowledgments:

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.

LITERATURE CITED

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.

APPENDIX 1

= DETAILS FOR THE RATIO-DEPENDENT MODEL

A1.1 Background

Ratio-dependent models assume that the functional response terms depend o= n ratios of the state variables rather than on absolute values or product= s of variables, as is the case for classical models. Although not a new i= dea, the concept of ratio dependence in predator-prey interactions has be= en approached with fresh interest in ecological theory in recent years (B= erryman 1992).

One advantages of such models is that they prevent the paradoxes of enric= hment and biological control predicted by classical models (Berryman 1992= ). Experimental observations of Arditi and Saïah (1992) suggest that= prey-dependent models are appropriate in homogeneous situations and rati= o-dependent models are appropriate in heterogeneous situations. Based on = their work, Ginzburg and Akçakaya (1992) and McCarthy et al. (1995= ) also conclude that natural systems are closer to ratio dependence than = to prey-dependence. Gutierrez (1992) develops a physiological basis for t= he theory.

Gleeson (1994), however, questions the assumptions of ratio-dependent mod= els, noting that direct density dependence, or self-regulation, in the to= p consumer is sufficient to preclude the paradox of enrichment from class= ical models. From his work on whether patterns among trophic levels are a= reliable way to distinguish between prey- and ratio dependence, Sarnelle= (1994) concludes that the ratio-dependent approach should be applied onl= y when predator and prey are the top two trophic levels in an ecosystem. = Abrams (1994) argues that patterns and experimental results used in suppo= rt of ratio-dependent predation are consistent with numerous other explan= ations, and that these other explanations do not suffer from pathological= behaviors and a lack of plausible mechanism, as do ratio-dependent model= s. Lundberg and Fryxell (1995) note the difficulty of distinguishing betw= een competing hypotheses without a proper mechanistic understanding of th= e processes involved.

In a recent paper, Akçakaya et al. (1995) respond to some of these= criticisms. The argument relevant to my present paper concerns their ref= utation of the pathological behavior of ratio-dependent models. Freedman = and Mathsen (1993) note that ratio-dependent models are invalid near the = axes (that is, where the state variables are close to zero), because the = ratios tend to infinity in these regions. As a result, even when prey (re= source) densities are very low, ratio-dependent models predict a positive= rate of predator (consumer) increase, provided that predator densities a= re low enough, because the number of prey available per predator increase= s to infinity as predator density declines to zero (Abrams 1994, Gleeson = 1994). In terms of isoclines, the problem stems from the fact that the pr= edator isocline passes through the origin in ratio-dependent models, whic= h means that, even at low prey densities, a sufficiently small predator p= opulation can increase. According to Hanski (1991), this is against intui= tion and many field observations. It also means that ratio-dependent mode= ls cannot be used to study the extinction of species.

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-dependen= t models cannot predict this outcome, they are pathological in a biologic= al 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 th= ese "pathological" regions (see Fig. 2 ). Although these large-amplitude oscillations may be interpreted as sign= aling extinction from a practical viewpoint, this is not the prediction o= f the ratio-dependent model.

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

A realistic mod= el 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.
However, = a few sentences later, they agree that ratio-dependent models are not val= id at very low densities (a precursor of extinction); earlier in the pape= r, they state that " it is at the extremes of low and high densities that= strict ratio dependence may not be valid."

In an attempt to clarify some of the arguments in this debate, I introduc= e 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 uns= table, because a small perturbation to the ratio-dependent terms substant= ially alters the dynamics.

A1.2 Model equations

The model equations are functionally homogeneous (i.e., all three equatio= ns 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 the= ir metabolic requirements. Details of the model formulation can be found = in Gutierrez et al. (1994). The final equations are:

  ,

where M1 is the biomass of plants, M2 is the biomass of herbivores, and M3<= /B> is the biomass of predators. The parameters = represent assimil= ation rates corrected for the efficiency of biomass conversion; = represents the pr= oportion of the resource that is available to its consumers (that is, its= apparency); Di is the per unit demand of the co= nsumer for resources; ri is the base respiration= rate corrected for the efficiency of biomass conversion; and bi is the degree of self-limitation. M0 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, riMi1 + bi, requires further explanation. Respiration usual= ly increases with population density (Gutierrez et al. 1994) and, thus, s= hould be an increasing function of Mi. However, = introducing such a functional dependence increases the model's complexity= =2E Because the effect is usually small, the authors chose the simpler fo= rmulation riMi1 + bi, where bi has a value between 0.02 = and 0.05. The disadvantage of this choice is that ri= must have rather unusual units (grams per bi per day, where grams are the units of Mi) that= depend on bi. Thus, the term riMi1 + bi has the same un= its as = (namely, grams p= er day, where grams are the units of Mi). This d= ependence of the units on bi is not satisfactory= from a mathematical viewpoint, but since bi 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. diversic= ornis (Howard) and native natural enemies cannot. However, the para= meter values for the cassava-mealybug-parasitoid system given to me by th= e authors (only values for = are reported in t= he 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 r= elative 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 sim= plify a system of equations because it has many other advantages. It illu= minates which parameters are most important in determining the dynamics o= f the model (Edelstein-Keshet 1988) and gives insight into the relative m= agnitudes of the parameters required to produce biologically reasonable b= ehavior. Also, the state variables are scaled so that they all have the s= ame order of magnitude, say between 0 and 1. This is important when solvi= ng the equations numerically, as very different magnitudes can lead to co= mputer round-off errors (Gerald and Wheatley 1989). In the cassava-mealyb= ug-parasitoid system, the biomass of cassava is much larger than that of = the mealybug and the parasitoid (cassava mass averages ~ 2 kg, whereas th= e mealybug and parasitoids have average masses of ~ 2 mg). Scaling, at le= ast, is thus vital. (For an introduction to nondimensionalizing systems o= f ordinary differential equations, see Edelstein-Keshet 1988.)

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

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

,

where K0 =3D M0 gives t= he nondimensionalized Eqs. A1.2, where I have replaced = by <= B>Mi and t* by t for co= nvenience:

 

(I replaced the product = by = , because all thre= e 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 her= e lend themselves to biological interpretation. For example, = can be thought of= as the potential per unit biomass growth rate (Gutierrez et al. 1994), o= r 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 nutr= itional value of the resource. = is made up of a r= atio of quantities. The numerator can be thought of as the maximum amount= of resource available to the consumer, and the denominator as the maximu= m demand of consumers for resource. More simply, = gives a measure o= f the ratio of supply to demand. Results in Gutierrez et al. (1994) are b= ased on the relationship between = and = , since = is just a scaling= factor. Hence, results using Eqs. A1.2 are comparable with those in Guti= errez et al. (1994).

I mentioned earlier that the parameters ri in t= he original model have units that depend on bi. = This can be prevented by replacing the terms riMi1 + bi with terms of the form = , where Ti has the same units as Mi. <= I>Ti can be thought of as a threshold value above whic= h self-limitation becomes noticeable. With this modification, the units o= f ri are per day, and ri can be interpreted as a respiration rate, as was originally intended. T= he new carrying capacities are given by

It can be shown that setting the Ki 's equal t= o 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 ig= nored in the remainder of the paper.

Having scaled the equations, I needed to choose values for the parameters= =2E An advantage of scaling is that there are now 11 parameters instead o= f the original 18. For comparison with the results of Gutierrez et al. (1= 994), I wanted to find values that resulted in isocline configurations si= milar to those in their paper. XPPAUT calculates and displays isoclines i= n two dimensions, and parameter values can be altered interactively. Sinc= e the competition effect is very small, but difficult to quantify, I foll= owed A. Gutierrez (personal communication) and chose B= i =3D 0.02 (i =3D 1,2,3). Values of , = , = , , , , , gave the isocl= ine configurations shown in Fig. 11 =2E These isoclines have similar shapes to those in Gutierrez et al. (199= 4). A noticeable difference is that the M3 -isoc= line intersects the M2 -isocline to the left of = the M2 -peak. In fact, using the techniques in G= utierrez et al. (1994), it would have been impossible to conclude that th= e tritrophic equilibrium is stable for the isocline configuration shown i= n Fig. 11, because of the position of the inter= section point in the plane. For this parameter set (which I shall call th= e reference set), there is also a stable limit cycle. The initial values = of M1,M2 and M3 determine whether the system approaches the stable equi= librium or the limit cycle.

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

Before beginning an analysis of the model, I will introduce a small modif= ication to the equations. Because the problems associated with ratio-depe= ndent 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 clo= se to zero, would have on the dynamics. Abrams (1994) states that modific= ations to ratio-dependent models cannot be made biologically realistic be= cause the original models have no clear mechanistic derivation. However, = some form of modification that prevents the ratios from tending to infini= ty may be useful for revealing any spurious behavior near the axes that m= ay result from the ratio dependence.

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

.

The difficulties are experienced when Mi approac= hes zero. Adding a constant in the denominator, that is, replacing the ra= tio by =

  =

would alleviate the problem. Although this addition may appear difficult = to justify biologically, Gutierrez (1992) used an exponent of this form i= n his functional response term. That model is physiologically based, as i= s 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 t= he 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 ai is a small constant, say 0.001. The res= ulting equations are:

  .

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

APPENDIX 2

DYNAMICAL SYSTEMS BASICS

A2.1 Glossary of dynamical systems terminology

This appendix lists alphabetically those topics that are required in the = main body of the paper and in the appendices. Brief descriptions are give= n and diagrams are used to explain the concepts as simply and intuitively= as possible. A more comprehensive introduction to dynamical systems theo= ry can be found in most introductory dynamical systems texts (e.g., Gucke= nheimer and Holmes 1983, Seydel 1988). Yodzis (1989) explains the theory = using ecological examples, and Strogatz (1994) gives an informal introduc= tion to nonlinear dynamics, with an emphasis on geometric intuition and s= cientific applications.

A2.1.1 Bifurcation diagram

A one-parameter bifurcation diagram summarizes the qualitative behavior c= orresponding to different values of a parameter. A state variable, or som= e function of the state variables, is usually plotted on the y -ax= is and the parameter on the x -axis. The positions and local stabi= lities of equilibrium points, as well as periodic orbits, are indicated u= sing different line types. Solid curves are used to represent locally sta= ble equilibria and dotted curves are used for locally unstable equilibria= =2E Maxima and minima of periodic orbits are indicated using circles: sol= id for stable orbits and open for unstable orbits.

It is important to note that bifurcation diagrams summarize the behavior = associated with a range of parameter values. They do not represent the dy= namics corresponding to a continually varying parameter (Wiggins 1990). T= o 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 peri= odic orbit (limit cycle). The local stability properties of a particular = phenomenon are given by the type of curve: solid, dotted, or open or clos= ed circles. For example, the phase portraits in Fig= =2E 14b (i) and (ii) were obtained by mentally drawing vertical lines at the par= ameter values µ =3D µ1 and µ =3D µ2, 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.


A two-parameter bifurcation diagram shows how the positions of bifurcatio= n points change as two parameters are varied. For example, if a bifurcati= on point is encountered in a one-parameter bifurcation diagram, a second = parameter may be varied to see how it affects the position of the bifurca= tion point. An example is shown in Fig. 19.

= 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 qualitativ= e behavior of the system changes. A stable equilibrium may become unstabl= e 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, A= 2.1.11, and A2.1.19.

A2.1.3 Continuation branch

A solution branch, or continuation branch, is a curve of equilibrium poin= ts (or limit cycles or bifurcation points) that indicates how the positio= n and properties of the equilibrium point (or limit cycle or bifurcation = point) change as a parameter (or parameters) is altered. Together, a numb= er 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 po= int (see section A2.1.9). Then the collection of all initial state variab= le values from which the system tends toward this equilibrium, as time pr= ogresses, is the domain (or basin) of attraction of the equilibrium point= =2E The equilibrium point is called an "attractor." Any stable phenomenon= , such as a stable limit cycle, also has a domain of attraction and is re= ferred 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 chan= ge from stable equilibrium behavior to stable limit cycles of large ampli= tude. This results from a subcritical Hopf bifurcation in which the unsta= ble 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 =2E Also, as µ is decreased,there is a jump from large-amplitude cycles to a zero-amplit= ude equilibrium point, but this takes place at µ1, which is less than µ*. For µ1 < µ < µ*, 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-dep= endent 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 bifurca= tion; Arnold 1983, Wiggins 1990) is a bifurcation point at which an equil= ibrium point alters stability and a limit cycle (period orbit) is initiat= ed. The example in Fig. 15 has a stable limit c= ycle surrounding an unstable equilibrium point. It is also possible to ha= ve an unstable limit cycle encircling a locally stable equilibrium point.= Unstable periodic orbits are indicated by open rather than solid circles= =2E

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 equili= brium 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 µ =3D µ1. A is a locally stable equilibrium point and B is a saddle p= oint. The initial values of x1 and x2 determine the subsequent behavior of the system. If the in= itial point is in the domain of attraction of A (to the right of p= oint B in Fig. 13b), then the system wil= l approach A. However, if the initial point lies on the other side= of B, it will be repelled away from B in the opposite dire= ction to A. Notice how the size of the domain of attraction of = A, in terms of the state variable x1, decrea= ses 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 equilib= rium point is said to be locally stable and is called an "attractor." Oth= erwise, it is said to be unstable and is a "repeller." Locally stable equ= ilibrium points are called sinks and locally unstable ones are called sad= dle points or sources.

A2.1.10 Parameter

A parameter is a quantity, e.g., fecundity rate or predation rate, that i= s used in describing the dynamics of a state variable. Although a state v= ariable evolves with time, a parameter is kept constant as time progresse= s. 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 equilib= rium 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 varia= ble (e.g., population density) that is oscillating in a regular, repetiti= ve manner. For continuous models having dimension of at least two, such a= n orbit is also called a limit cycle, provided that the orbit is isolated= (that is, trajectories starting at points near the orbit spiral either t= oward 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 = (x1) and a predator (x2= ). We can represent the behavior over time of both populations in a singl= e diagram called a phase portrait. Examples are shown in Figs. 13b and 14b =2E An illustration of the relationship between time plots and phase port= raits can be found in Holling (1973: 3).

A2.1.14 Saddle point

A saddle point is an equilibrium point that attracts in certain direction= s 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 poi= nts 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 n= ode or a spiral attractor. In the case of a spiral attractor, the state v= ariables oscillate with decreasing amplitude as they approach the equilib= rium 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 ampli= tude of these cycles increases gradually for parameter values further fro= m 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 disturbanc= e to the system will cause the state variables to move away from this poi= nt. 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, a= nd predators. Then the"state" of the system can be described by the relat= ive biomasses or densities of these populations. The variables that are u= sed 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 14<= /a>b(i) and (ii) show possible phase portraits in two dimensions for µ =3D µ1 and µ =3D µ2. Note that the bifurcation diagram in Fig. 14a o= nly indicates the positions of the equilibrium points in terms of one of = the state variables, x1.

A2.2 Introduction to dynamical systems concepts=

In this paper, I am more interested in obtaining qualitative tha= n quantitative information about model behavior. That is, I am mo= re interested in finding out how varying a parameter value affects the mo= des of behavior of a model, rather than finding the exact parameter value= s at which these changes occur. Berlinski (1976) notes that qualitative i= nsights are at a much greater depth than partially quantitative results. = Because of the uncertainty associated with parameter values in nature (Ha= damard 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 mos= t important.

Qualitative analyses are concerned mainly with long-term rather than tran= sient behavior. Transient dynamics vary with the initial values of variab= les and the time period over which solutions are calculated, whereas qual= itative studies deal with the eventual behavior of the system once the in= itial transients have died away.

There are many different modes of behavior that characterize the long-ter= m dynamics of a model. Suppose that the model=D8s state variables represe= nt interacting populations. In the long term, any given population may be= come extinct, reach a stable equilibrium, or oscillate in time. Another m= athematical possibility is that a population may grow indefinitely, altho= ugh 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 predato= r 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 o= f the predator, the prey population experiences logistic growth. The term= - cy describes the exponential decline in the predator pop= ulation in the absence of prey. The Holling type-II functions (Holling 19= 65)