## Abstract

Biologists use mathematical functions to model, understand and predict nature. For most biological processes, however, the exact analytical form is not known. This is also true for one of the most basic life processes: the uptake of food or resources. We show that the use of several nearly indistinguishable functions, which can serve as phenomenological descriptors of resource uptake, may lead to alarmingly different dynamical behaviour in a simple community model. More specifically, we demonstrate that the degree of resource enrichment needed to destabilize the community dynamics depends critically on the mathematical nature of the uptake function.

## 1. Introduction

Simple community models predict that increasing the availability of resources will destabilize community dynamics from equilibria to oscillatory dynamics, a phenomenon termed the ‘paradox of enrichment’ (Rosenzweig 1971; Gilpin 1972; May 1972; Myerscough *et al*. 1996). Attempts to establish this effect in experiments or in the wild have met with only partial success (McCauley & Murdoch 1987, 1990; Persson *et al*. 2001), indicating that real communities respond to enrichment in a more complicated way than simple models suggest. Environmental conditions and properties of the community (web-like structure, shift to inedible prey, inducible defences) have been offered as explanations for why communities might fail to destabilize as a consequence of enrichment (Persson *et al*. 2001; Vos *et al*. 2004). Here, we offer a different explanation that is related to the properties of the mathematical model which was first used to put forth the ‘paradox’.

The classical Rosenzweig–MacArthur (R–M) model (Rosenzweig & MacArthur 1963) is probably the simplest formulation of a trophic community able to produce realistic dynamic behaviour (Turchin 2003). The model describes the changes over time of two populations coupled by predation:The prey population, *x*, grows logistically at the rate *g*(*x*)=*rx*(1−*x*/*K*), where *r* is the growth rate of the prey and *K* is the carrying capacity. The predator, *y*, consumes the prey and grows according to the nonlinear uptake function (functional response), *f*(*x*), and has mortality, *m*. As is common practice (Rosenzweig 1971; Murdoch *et al*. 1998), enrichment is simulated in this model by increasing *K*, i.e. by allowing the prey to grow to higher densities in the absence of predators.

We investigated whether the specific mathematical formulation of the functional response affects the community dynamics that the R–M model predicts. The surprising result is that the degree of destabilization caused by enrichment is extremely sensitive to the mathematical nature of the uptake function, that is, even functional response curves that are undistinguishable for all practical purposes may produce qualitatively and quantitatively different dynamics.

## 2. Methods

Minimum requirements for realistic uptake functions, *f*(*x*), are that the function be zero at zero resource concentration, be monotonically increasing with resource density, and be saturating when resource density goes to infinity (Myerscough *et al*. 1996). To maximize similarity among different functional response curves we restrict ourselves to functions with negative curvature over the whole prey range (table 1; figure 1*a*). Ivlev's function (Ivlev 1961) and Holling's type II function (Holling 1959) are the most widely used functions that fulfil these requirements, but others, e.g. trigonometric functions, have been proposed (Jassby & Platt 1976). Given the error with which resource uptake by real organisms is measured, it is usually unjustified to identify a best-fitting model, and structurally different analytical forms may be used interchangeably. For our theoretical investigation, we chose a generic parameterization of Ivlev's functional response (*a*_{I}=1, *b*_{I}=2; but see electronic Appendix A) and used nonlinear least squares to maximize the phenomenological similarity with Holling's type II response and a response curve based on a trigonometric function (table 1; figure 1*a*).

Graphical analysis revolves around plotting the curves (isoclines) in the prey–predator phase plane that denote zero growth of the model predator and prey populations. For Ivlev's and Holling's functional responses, the R–M model produces vertical predator and hump-shaped prey isoclines, although the left portion of the hump may be hidden in the region of negative prey densities (figure 1*b*,*c*). Intersections of predator and prey isoclines mark equilibrium points. If the intersection occurs in the right decreasing portion of the prey isocline, the equilibrium is locally stable and population densities starting in the neighbourhood of the equilibrium will converge in this point. Intersections in the left increasing portion of the prey isocline mark unstable equilibria, which give rise to oscillatory, limit cycle dynamics (Gilpin 1972; May 1972). Using the trigonometric response function within the R–M model results, for certain parameter sets, in a prey isocline with multiple decreasing portions, which renders analysis more difficult. Besides graphical analysis, the strength of stability is evaluated by calculation of the eigenvalues of the community matrix at equilibrium (e.g. Edelstein-Keshet 1988). Further, we use numerical integration of the R–M model to confirm the results from the stability analysis and also to determine global stability of the equilibrium.

## 3. Results

It becomes immediately apparent that almost identical resource uptake curves (figure 1*a*) give rise to very differently shaped prey isoclines (figure 1*b*,*c*), which have drastic consequences for the dynamic stability of the system. First, we study the case where the carrying capacity is set to unity, *K*=1 (figure 1*b*). In this case, only Holling's function has the potential to destabilize the population dynamics because only its prey isocline has an increasing portion for positive prey densities. With the two other functional responses, the prey isocline is decreasing and predator–prey dynamics settle on stable equilibria in both cases. This is also obvious from the simulated time-series of predator and prey densities starting from initial values close to the equilibrium (Electronic Appendix A). Using Holling's functional response, the equilibrium is unstable and the trajectories settle into a stable limit cycle with large amplitudes. By contrast, the trigonometric function leads to a stable equilibrium and the trajectories quickly spiral into the steady state. Ivlev's function produces a dynamical behaviour that is intermediate between these extremes. The equilibrium is stable but only very weakly so (see also figure 2); initially the trajectories undergo weakly damped oscillations which, however, eventually settle into the equilibrium.

To conclude, even though the resource uptake curves are nearly identical, the resulting time-course of predator and prey densities is very different in the three cases. This also has consequences for the extinction risk. Using the trigonometric response curve, both predator and prey abundance are always far from zero and therefore the populations have a high expectation of persistence. By contrast, with Holling's function the oscillating densities pass through very small values, putting them at a high risk of extinction. Again, the Ivlev response curve leads to intermediate behaviour. Initially, while the time solutions are still oscillatory there is a moderate extinction risk. However, with increasing time, the amplitude of the oscillations, and therefore the extinction risk, are more and more reduced.

With enrichment (by raising the carrying capacity), all three functions become potentially destabilizing, but not to the same degree (table 1; figure 2). Enrichment, by a fourfold increase of the carrying capacity (*K*=4) and leaving all other parameters unchanged, leads to limit cycles if Holling's or Ivlev's responses are used (figure 1*c*). However, with the same amount of enrichment the equilibrium remains locally stable in a system that is based on the trigonometric function. Nevertheless, for an appropriate choice of initial values, limit cycle oscillations can also be observed in this system, because in this range of intermediate enrichment, the trigonometric function results in multi-stability with coexistence of stable equilibrium and oscillatory dynamics (table 1). The system needs to be enriched even more (*K*>10.12) before destabilization occurs globally (table 1). Thus, the enrichment level at which the equilibrium is destabilized varies by a factor of more than 20 (table 1; figure 2). Although enrichment eventually leads to destabilization in all models, the vastly differing conditions at which it occurs will be disconcerting to anyone using mathematical models as a predictive tool.

The three functional responses can be ranked according to their potential to destabilize the dynamics of the R–M model (Holling II>Ivlev>trigonometric function). We found this pattern to be very general and extremely robust over wide ranges of the parameters *a*_{I} and *b*_{I}, which determine the steepness and saturation level of the uptake function. The ranking remained largely constant even when we constrained our fits either to have the same slope at the origin or the same saturation level as Ivlev's function (Electronic Appendix A).

## 4. Discussion

The dynamic stability of the R–M model and other models has been shown to depend on the *type* of functional response used (Oaten & Murdoch 1975; Armstrong 1976; Williams & Martinez 2004). It is important to understand that our result is of a very different nature. We did not study the effect on system stability of different functional response curves that represent mechanistically motivated alterations of a basic function, e.g. Holling's type I, II and III responses (Holling 1959). On the contrary, our goal was to investigate the effect of response functions that are as similar as possible phenomenologically, but are derived from entirely different mathematical foundations. What we found is that three functions, which are identical for all practical purposes, give completely different outcomes in terms of model dynamics, a phenomenon called sensitivity to model structure (Wood & Thomas 1999).

Sensitivity to model structure has been described for the R–M (Myerscough *et al*. 1996) and other ecological models (Wood & Thomas 1999; Gross *et al*. 2004), but we offer a simple explanation for this striking phenomenon based on the structural similarity of the mathematical functions that occur in the R–M model. Here, logistic prey growth, *g*(*x*), and resource uptake, *f*(*x*), are structurally very similar at small prey numbers, *x*. This has consequences for the prey isocline, (Recall that the stability of the equilibrium depends on the slope of at equilibrium.) In the extreme case that the two functions are exactly identical, the isocline is a constant and has a slope of zero everywhere. Thus, linear stability is not well defined and the system is structurally unstable. But assume that in some range close to the equilibrium point the two functions differ slightly, except for a constant: *cf*(*x*)=*g*(*x*)+ϵ(*x*), where ϵ(*x*) is a small function. Now the isocline can be approximated as and in this range the sign of the slope of entirely depends on the difference, ϵ(*x*). Therefore arbitrary small deviations of resource uptake, *f*(*x*), from the prey growth rate, *g*(*x*), determine the stability of the equilibrium.

We emphasize that the whole argument relies on the fact that, in the relevant interval, prey growth, *g*(*x*), and resource uptake, *f*(*x*), are structurally very similar functions. In the R–M model this is always the case for small prey numbers because both *g*(*x*) and *f*(*x*) start from zero as negatively curved functions (this becomes apparent from a Taylor expansion of *f*(*x*) up to second order). Thus, whenever the equilibrium is at small prey levels, e.g. for small mortality, *m* (as in figure 1), the R–M model is sensitive to minor variations in the form of the functional response curve. By contrast, if the equilibrium is at large values of *x*, where prey growth and resource uptake are significantly different functions, e.g. at large or intermediate levels of mortality, *m*, the effect of sensitivity to the model structure is not observed.

The same mechanism lies at the heart of one of the major drawbacks of the original Lotka–Volterra model, which is a special case of the R–M model where prey growth and resource uptake are linear functions, i.e. *g*(*x*)=*ax* and *f*(*x*)=*bx*. In this model the expression *g*(*x*)/*f*(*x*) is constant per definition, which leads to structural instability in the whole parameter range. With the introduction of nonlinear logistic prey growth and saturating functional response, Rosenzweig and MacArthur tried to circumvent these problems. Indeed this works out for most parameter ranges. However, as shown above, the same difficulties of sensitivity to infinitesimal variations in the model structure are still inherent in the R–M model and are able to enter through the back door in, for example, cases of small mortality.

Sensitivity to model structure may be responsible for the failure to observe destabilization as a result of enrichment in real communities. This is not to say that mechanistic explanations should generally be ruled out. McCauley *et al*. (1999), for example, showed convincingly that competition between inedible and edible prey can reduce the effective carrying capacity of the prey and thereby stabilize the community dynamics. We caution, however, that such conclusions should not be drawn prematurely, based on a mismatch between a particular theoretical model and empirical results. Instead, we advise an evaluation of the robustness of model predictions for alternative mathematical formulations whenever the exact mechanistic nature of the resource uptake is not known. Currently, ecological modellers use Holling's type II function as a standard in consumer–resource models although the true uptake mechanism may vary across and within communities (Jeschke *et al*. 2002; Mols *et al*. 2004). Our analysis has demonstrated that a much higher degree of enrichment may be required for destabilization than such standard models predict. We conclude that, unless the exact mechanistic nature of the relationship between consumer and food is known, caution should be used when predicting predator–prey dynamics and the effect of enrichment without considering the effects of sensitivity to model structure.

## Acknowledgements

We thank A. Gonzalez, F. Guichard and two anonymous referees for constructive comments on an earlier draft of this paper.

## Footnotes

↵† Present address: Biology Department, McGill University, 1205 Avenue Docteur Penfield, Montreal, Quebec, H3A 1B1, Canada.

The supplementary Electronic Appendix is available at http://dx.doi.org/10.1098/rsbl.2004.0246 or via http://www.journals.royalsoc.ac.uk.

- Received July 30, 2004.
- Accepted August 26, 2004.

- © 2005 The Royal Society