

Home  Archives  About  Login  Submissions  Notify  Contact  Search  
Copyright © 1997 by The Resilience Alliance^{*} 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/ To print the entire article (including any tables, figures, or appendices), use this link to the unified version of the article. To print separate parts of the article, click in the frame containing that part (text, figure, table, or appendix) before choosing File, Print. Research Automated Techniques for the Qualitative Analysis of Ecological Models: Continuous Models Lynn van Coller^{1} Department of Mathematics,University of British Columbia
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. Indepth 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 predatorprey model and establishes the reliability of the computer software. The main body of the paper discusses a ratiodependent model. The new techniques highlight some limitations of isocline analyses in this threedimensional setting and show that the model is structurally unstable. Another appendix describes a larger model of a sheeppasturehyraxlynx 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; ratiodependent 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 ratiodependent 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 threedimensional setting, the analysis shows that the model is structurally unstable. Even small perturbations to the ratiodependent terms alter the dynamics. 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 predatorprey 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. 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. 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 plantherbivorepredator system. Currently, there is widespread debate over the validity of such ratiodependent 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 ratiodependent 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 ratiodependent 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 EdelsteinKeshet 1988), is unsuitable in this threedimensional 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 twodimensional 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 ratiodependent 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 ratiodependent model are encountered. I also investigate a few other values.
2.2 Model analysis
2.2.1 Oneparameter studiesGutierrez 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 parametersThe parameter can be thought of 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 M_{3} zero isocline is given by , and the solution 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 original and the modified model are shown in Fig. 1, with M_{1} plotted on the yaxis. 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. 1a, there are two stable limit cycles. The range of values is so small, however, that it is not of much biological significance. Analysis of the temporal dynamics of the system (using XPPAUT) for different values of shows that larger values of decrease the time taken to reach equilibrium. Increasing , the potential growth 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 a_{i}' 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. From these figures, it can be seen that as increases, there is a general increase in the M_{1} 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_{1} equilibrium value approaches the M_{1} 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 perioddoubling 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_{2} and M_{3} 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 ratiodependent models. This will be discussed in more detail later. 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 become difficult to evaluate numerically as M_{i} 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 do 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 Cshaped M_{2}isocline. It is true that if is decreased, the M_{2}isocline widens. However, this is provided that M_{3} is constant. If the system is integrated and M_{3} is allowed to vary until a new equilibrium is reached, and the isoclines are plotted with this new equilibrium M_{3} value, then the final M_{2}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 parameterVarying , the assimilation or conversion efficiency of the herbivore when plants are abundant, gives the bifurcation diagrams in Fig. 5. Comparing Fig. 5b with Fig. 3a reveals that they are almost mirror images. If the mean of the cycle maxima and minima in Fig. 5a and Fig. 2a 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. Twoparameter diagrams in (, )space can be generated by continuing the Hopf bifurcations in Fig. 5 in two parameters. The results are shown in Fig. 6. These diagrams show clearly that decreasing or increasing has a similar effect, 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 oneparameter bifurcation diagrams give us the additional information that, for certain parameter ranges, limit cycles are possible. 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. 6a; this problem does not occur in Fig. 6b. A closer investigation reveals that the equilibrium values for the state variables are close to zero in the upper left triangle of the twoparameter space, and this results in numerical problems when using the original model. Figure 6b 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 M_{1} 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 twoparameter space are less complex. These claims are made clearer in Fig. 7 and Fig. 8. Figure 7 shows the (, )space for the modified model with a_{i} =0.002 (i = 1,2,3). The region of stable equilibria is even larger than in Fig. 6b, resulting in smaller regions of cycles. The presence of the a_{i} '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 lefthand section of the upper region of cycles in Figs. 6a, b, and Fig. 7. The lefthand portion of this region is where the values of M_{2} and M_{3} are low and, hence, where the nonzero a_{i} 's have most effect. Clearly, nonzero a_{i} '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 (a_{i} =0, i = 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 configurationsFrom 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 M_{1}M_{2}(M_{2}M_{3}) phase plane depends on the value of M_{3}(M_{1}) 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 M_{1}M_{2} plane corresponds to a tritrophic equilibrium. Figure 9 shows three possibilities. Two of these (Fig. 9b and c) appear in Gutierrez et al. (1994), who assumed, incorrectly, that the equilibrium point in Fig. 9b was unstable.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 unstable 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 M_{1} and M_{2} isoclines in the M_{1}M_{2} and M_{2}M_{3} 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). A comparison of Fig. 10a and b demonstrates that introducing the a_{i}'s prevents the M_{1} and M_{2} 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 a_{i}'s from 0.001 to 0.005 reduces the humped shape of the M_{1} and M_{2} isoclines. The result is an even more robust model. In fact, when a_{i} = 0.005 (i = 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 ratiodependent 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 ratiodependent 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 threedimensional 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 ratiodependent terms has a significant effect on the dynamics. This supports the argument that ratiodependent models exhibit pathological behavior and are not valid near the axes (see the background information for ratiodependent 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 sheephyraxlynx 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 predatorprey 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 sheephyraxlynx 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 ratiodependent 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 predatorprey model. In the sheephyraxlynx 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 ratiodependent 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 ratiodependent terms. The dynamical systems techniques also highlighted the limitations of an isocline analysis in a threedimensional 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.
Address of Correspondent: Lynn van Coller 28 Woodside Avenue Cowies Hill, 3610, South Africa phone: 2731728013 fax: 2731728013 lvcoller@nbsmail.nbs.co.za ^{*}The copyright to this article passed from the Ecological Society of America to the Resilience Alliance on 1 January 2000.


Home  Archives  About  Login  Submissions  Notify  Contact  Search  
