## Ft tf 211

denotes the fraction of the population Y that will (in the steady state) eventually be consumed by the predator Z. The complementary parameter

ay Y

denotes the fraction of the population that will eventually die because of natural mortality. In terms of these scale parameters Eq. (2.9) can be written as y = ay[gy(x,y) - py9z(y,z) - fjym(y)]. (2.13)

In this example we have managed to find interpretable scale parameters by introducing one parameter that denotes the scale of the total turnover, ay, and subsequently measuring the relative contributions to this turnover. Even in much more complicated models this procedure generally succeeds to reveal easily interpretable scale parameters.9 In some cases it can be useful to introduce multiple levels of grouping. Suppose for instance that the equation of motion contained multiple loss terms that arise from the predation by different predators. In this case we could use one scale parameter ay to denote the total turnover, then another scale parameter to denote the relative contribution of the sum of all predation terms to the total turnover and finally a third parameter Yy,i to denote the relative contribution of the predation by a certain predator i to .

The second difficulty, that can arise in the construction of a generalized model, is that the individual terms in the model can be conceptually more complicated. Let us illustrate this situation by the well studied example of predation on multiple prey populations.12 In comparison to a single prey population, this situation is for two reasons more complicated. First, we know that some relations between the losses of the prey and the gain of the predator exist. While these relations should be reflected in the model, the losses of either prey are no longer directly proportional to the total gain of the predator. Second, some derivatives can arise which do not have a direct intuitive interpretation. For example it is not always intuitively clear how the loss rate of one prey population responds to a variation in the population density of the other. Both of these problems can be solved by including some additional mechanistic reasoning, which enters the model in the form of auxiliary variables and equations.

Let us denote the two prey populations by X and Y and the predator population by Z. We use the function G(X, Y, Z) to describe the gain of the predator by predation and the functions Lx (X, Y, Z) and Ly (X, Y, Z) to describe the predative losses of the prey populations. In addition we introduce the auxiliary variable P which denotes the total amount of prey that is perceived by the predator Z. Let us assume that P can be written as a sum

where Cx and Cy are general positive functions that describe the contribution of the populations X and Y depending on the respective population sizes. While it is in many cases reasonable to assume that these functions are linear, they can be nonlinear if, for instance, the predators can improve the success rate of attacks with practice.12 We can normalize auxiliary equations, like Eq. (2.14), by applying the same normalization procedure that we have used for the differential equations. In the notation introduced above the normalized auxiliary equation reads as p(a

We identify the constant factor p = Cx*/P* as a scale parameter which denotes the relative contribution of population X to the total amount of available prey, while the complementary variable p = 1 — p = C*/P* denotes the fraction contributed by population Y. This allows us to write the normalized amount of available prey as

Let us now investigate how the losses of population X relate to the gain of Z (the losses of Y are completely analogous and hence will not be treated separately). Since population X contributes a fraction Cx/P to the available amount of prey, it can be assumed that it contributes the same fraction to the captured amount of prey. From this we deduce the form of the corresponding loss rate as

The normalization of this equation yields

Thus, by introducing the auxiliary variable P we have managed to determine the relation between the predative losses of the prey populations and the gain of the predator. However, the main advantage lies in the fact that the derivatives of the auxiliary variables with respect to the normalized state variables have a much more direct interpretations. In the Jacobian all terms relating to predation can be expressed as exponent parameters by the following derivatives iP9{p>z)

Here, the exponent parameter gp describes the nonlinearity of the predation rate with respect to prey density, while gz describes its dependence

on the predator density. These two parameters are completely analogous to the parameters gx and gy that have been introduced in the previous section to describe the predator-prey system with a single prey population. The ability to describe structurally different systems with directly comparable parameters is one of the advantages of generalized modeling. The two new parameters cx and cy describe the nonlinearity of the contributions of the two prey populations to the total amount of prey. For example the case of passive prey switching corresponds to cx x = cy y = 1 while active prey switching can lead to larger values.

This example of predation on multiple prey populations illustrates that additional constraints can be taken into account in generalized models by including auxiliary equations. The introduction of such auxiliary equations is often useful since it makes room for additional theoretical reasoning, which can greatly enhance the interpretability of a given model without introducing too many new assumptions.

### 2.5. A generalized spatial model

The investigation of generalized models proposed here is not limited to models that are formulated in the language of ordinary differential equations, but can be extended for example also to systems of partial differential equations (PDEs). In ecology PDEs are frequently used to describe ecological populations in physical space. The underlying assumption in these models is that, at a certain scale, the evolution of population densities is captured by a diffusion equation. It is well known that in reaction-diffusion systems instabilities with respect to spatially inhomogeneous perturbations with a certain wavenumber k can exist.13 The corresponding qualitative transition in the phase portrait of the system is known as Turing bifurcation and wave instabilities. Beyond a Turing bifurcation spatially inhomo-geneous patterns form spontaneously from an initially homogenous state. This transition has been extensively studied in conventional models.14-18 More recently it has been discovered as the driving force of pattern formation in certain ecological systems.19-21

To illustrate these ideas, we consider a system of partial differential equations (PDEs) that was recently studied in Ref. 22, in which the simple predator-prey system Eq. (2.1) is extended to a spatial system. Thus we describe the dynamics in one point of physical space by the equations

X = S(X) - G(X, Y)+ DxAX, Y = r/G(X, Y) - M (Y) + Dy AY, [ ' '

where A denotes the Laplace operator, and Dx and Dy are diffusion, or dispersal , constants.

At first glance it seems that our analysis of the generalized model is complicated by diffusion. The diffusion term is neither a pure gain nor a pure loss term, but a mix of both. In particular, in a homogenous equilibrium it vanishes. This means that the normalization procedure described above cannot be applied to the diffusion term. However, recall that the main purpose of the normalization was to map unknown rates of the processes in the steady state to a known position. Since we know that the diffusion term vanishes in a homogeneous state, we can consider this case without normalizing the diffusion term. Moreover, the vanishing diffusion term does not interfere with the normalization of the other states in the model. We therefore obtain the normalized equations

y = ay (g(x, y) - m(y)) + Dy Ay, where Dx and Dy now act as scale parameters describing the diffusion.22 In order to investigate the stability of this system one considers the Jacobian with respect to perturbations with a wavenumber k which is given by22

Note, that for homogenous perturbations (k = 0) this Jacobian is identical to the one of the well-mixed system given in Eq. (2.6).

2.6. Local stability in small and intermediate models

In the previous sections the formulation and normalization of generalized models has been discussed. In the following we will be concerned with some ways in which information can be extracted from the resulting models. The Jacobian matrices computed from generalized models are in general simple in the sense that they do not contain complicated terms that usually arise in conventional models from the computation of steady states. In systems of small (dimension N < 4) or intermediate (N < 10) size, it is therefore often possible to compute the local bifurcations analytically.

In systems of ODEs local bifurcations of steady states occur if the variation of a parameter causes the real part of one or more eigenvalues of the Jacobian to change sign.1 Eigenvalues generally either cross the imaginary axis as a pair of two complex conjugate eigenvalues, or pass through the origin of the complex plane as a single real eigenvalue. The first case corresponds to a Hopf bifurcation which, at least transiently, gives rise to oscillations as the stability of the steady state is lost. The latter case corresponds to bifurcations of saddle-node type (e.g. fold, transcritical or pitchfork bifurcations) in which the number and/or stability of steady states changes. It is interesting to note that the direct computation of both of these types of bifurcations is in general less difficult than the computation of the eigenvalues themselves or the computation of steady states in a conventional model. The computation of eigenvalues involves the factorization of a polynomial of degree N which analytically is in general only possible for N < 4. By contrast, a test function that describes the local bifurcation points can always be constructed. The determinant of the Jacobian is a convenient test function that vanishes in (and in general only in) bifurcation points of saddle-node type. By applying slightly more involved techniques analogous test functions for the computation of Hopf bifurcations can be constructed.8'23'24 While these techniques can in principle be applied in systems of any size, the resulting expressions become too long to handle analytically in large systems (N > 10).

In small and intermediate systems the analytical computation of local bifurcations of steady states is a very efficient tool for the investigation of generalized models. For instance in the predator-prey model proposed in Sec. 2.3 we find bifurcation points of saddle-node type at

my and Hopf bifurcation points at ay (my - 9y) f «x(my - gy) io oa\

In order to find the Turing bifurcation points one formulates a condition for the existence of a positive eigenvalue and then considers the wavenumber for which this condition is first satisfied. This calculation is shown in detail in Ref. 22. As a result we find that the Turing bifurcation points are located

Fig. 2.1. Three-parameter bifurcation diagrams of generalized predator-prey systems. Left: the ODE model from Eq. (2.1), right: the spatial PDE model from Eq. (2.20). The bifurcation surfaces are shown in dependence on the prey sensitivity gx, the timescale separation r = ay /ax and the exponent of closure my. In both diagrams the steady state under consideration is stable in the top-most volume of parameter space. If gx is decreased destabilization occurs in a Hopf bifurcation (red surface) or in a bifurcation of saddle-node type (blue surface) or-in the spatial model-in a Turing bifurcation (green surface). On the lines, on which two surfaces meet, codimension-2 bifurcation points are located. Other parameters: sx = 0.5, gy = 1 and d = 30.

Fig. 2.1. Three-parameter bifurcation diagrams of generalized predator-prey systems. Left: the ODE model from Eq. (2.1), right: the spatial PDE model from Eq. (2.20). The bifurcation surfaces are shown in dependence on the prey sensitivity gx, the timescale separation r = ay /ax and the exponent of closure my. In both diagrams the steady state under consideration is stable in the top-most volume of parameter space. If gx is decreased destabilization occurs in a Hopf bifurcation (red surface) or in a bifurcation of saddle-node type (blue surface) or-in the spatial model-in a Turing bifurcation (green surface). On the lines, on which two surfaces meet, codimension-2 bifurcation points are located. Other parameters: sx = 0.5, gy = 1 and d = 30.

= 3 Vdy - \ mY + for sx + r(gy - my) < < ~{gy - my), at

These results are visualized in a three-parameter bifurcation diagram shown in Fig. 2.1, where we have assumed intermediate nutrient availability sx = 0.5 and the absence of intraspecific competition between predators gy = 1. In the three dimensional parameter space the bifurcation points of Hopf and saddle-node type form surfaces, which divide the parameter space into regions of qualitatively different long-term dynamics. The normalized steady state is stable in the topmost volume of the parameter space. As the prey sensitivity is lowered the steady state loses its stability as a Hopf bifurcation point (red surface), a bifurcation point of saddle-node type (blue surface) or the Turing bifurcation (green surface) is encountered.

In small and intermediate systems one can obtain a good impression of the full local bifurcation structure of the system by considering several of such three-parameter bifurcation diagrams with different axes. If analytical expressions for the bifurcation surfaces are available, then these diagrams can be generated without much effort. By visual inspection of the bifurcation diagrams one can usually tell the way in which the individual parameters effect the dynamics of the system. Once such an intuition is gained it can be verified mathematically. For instance, in the case of our general predator-prey system the sensitivity of the predator gx has a strong stabilizing effect. By increasing the value of gx one can always stabilize, but never destabilize a steady state. Furthermore, since my > gy in almost all systems, Eq. (2.24) shows that the critical value of gx at which the Hopf bifurcation occurs decreases as r = ay/ax is increased. This result is counter-intuitive since it implies that oscillations are less likely if the timescale separation indextimescale!separation between predator and prey is small.

In Ref. 22 a similar way of reasoning was used to identify the conditions under which the spontaneous formation of spatial and spatio-temporal patterns in predator-prey systems is likely. In particular it was shown that high nutrient supply, low competition for nutrients among prey, high abundance of prey and predators, strong intraspecific competition in the predator population and density dependent predator mortality promote the spontaneous pattern formation. Since all of these are typically found in enriched systems, these results indicate that anthropogenic eutrophication could lead to the formation of spatial or spatio-temporal patterns in natural predator-prey systems. A similar conclusion was reached in Ref. 25 based on the investigation of a conventional model.

Another interesting effect connected to eutrophication is the so-called paradox of enrichment. This paradox revolves around the observation that many ecological systems can be destabilized by increasing the supply of nutrients or prey.26'27 While this was initially felt to be counter-intuitive, the effect is now well understood. From a modern perspective the true paradox lies in the fact that many ecological systems observed in experiments are stabilized by an increase of nutrients or prey while almost all models predict a destabilization.28-30 In the past several solutions to this paradox have been pointed out,30-33 among them is the formation of spatio-temporal patterns mentioned above.25

Our work on generalized models suggests a different solution: The purely destabilizing effect of enrichment that is observed in many models may be an artifact, that is produced because of the specific functional forms that are usually employed in modeling.5 Let us focus on the functional response G(X, Y) and the corresponding parameter gx. We have already discovered that high values of gx have a stabilizing effect on the