# «Generalized modeling of ecological population dynamics Justin D. Yeakel · Dirk Stiefs · Mark Novak · Thilo Gross Received: date / Accepted: date ...»

Journal of Theoretical Ecology manuscript No.

(will be inserted by the editor)

Generalized modeling of ecological population dynamics

Justin D. Yeakel · Dirk Stiefs · Mark

Novak · Thilo Gross

Received: date / Accepted: date

**Abstract**

Over the past years several authors have used the approach of generalized

modeling to study the dynamics of food chains and food webs. Generalized models come

close to the eﬃciency of random matrix models, while being as directly interpretable

as conventional diﬀerential-equation-based models. Here we present a pedagogical introduction to the approach of generalized modeling. This introduction places more emphasis on the underlying concepts of generalized modeling than previous publications. Moreover, we propose a shortcut that can signiﬁcantly accelerate the formulation Justin D. Yeakel Department of Ecology and Evolutionary Biology 1156 High Street University of California Santa Cruz, CA 95064 USA E-mail: jdyeakel@gmail.com Dirk Stiefs Max-Planck Institute for the Physics of Complex Systems N¨thnitzer Str. 38 o 01187 Dresden Germany E-mail: stiefs@pks.mpg.de Mark Novak Long Marine Lab 100 Shaﬀer Road Santa Cruz, CA 95060 USA E-mail: mnovak1@ucsc.edu Thilo Gross Max-Planck Institute for the Physics of Complex Systems N¨thnitzer Str. 38 o 01187 Dresden Germany E-mail: thilo.gross@physics.org of generalized models and introduce an iterative procedure that can be used to reﬁne existing generalized models by integrating new biological insights.

Keywords Omnivory · Generalized Modeling · Bifurcation · food chain · food web 1 Introduction Ecological systems are fascinating because of their complexity. Not only do ecological communities harbor a multitude of diﬀerent species, but even the interaction of just two individuals can be amazingly complex. For understanding ecological dynamics this complexity poses a considerable challenge. In conventional mathematical models, the dynamics of a system of interacting species are described by a speciﬁc set of ordinary diﬀerential equations (ODEs). Because these equations are formulated on the level of the population, all complexities arising in the interaction of individuals must be cast into speciﬁc functional forms. Indeed, several important works in theoretical ecology present derivations of functional forms that include certain types of individual-level eﬀects [16, 30, 3, 8, 6]. Although these allow for a much more realistic representation than, say, simple mass-action models, they cannot come close to capturing all the complexities existing in the real system. Even if detailed knowledge of the interactions among individuals were available and could be turned into mathematical expressions, these would arguably be too complex to be conducive to a mathematical analysis.

In this light the functional forms that are commonly used in models can be seen as a compromise, reﬂecting the aim of biological realism, the need to keep equations simple, and often the lack of detailed information.

Because of the many unknowns that exist in ecology, it is desirable to obtain results that are independent of the speciﬁc functional forms used in the model. This has been achieved by a number of studies that employed general models, in which at least some functional forms were not speciﬁed [7, 21, 4, 23, 20, 24, 44]. These works considered not speciﬁc models, but rather classes of models comprising simple, commonly used, functions, as well as the whole range of more complex alternatives.

That ecological systems can be analyzed without restricting the interactions between populations to speciﬁc functional forms is in itself not surprising–in every mathematical analysis the objects that are analyzed can be treated as unknown. The results of the analysis will then depend on certain properties of the unknown objects. In a general ecological model we thus obtain results that link dynamical properties of the model, e.g. the presence of predator-prey oscillations to properties of the (unknown) functions describing certain processes, e.g. the slope of the functional response evaluated at a certain point. Accordingly, the analysis of general models reveals the decisive properties of the functional forms that have a distinctive impact on the dynamics.

Whether such results are ecologically meaningful depends crucially on our ability to attach an ecological interpretation to the decisive properties that are identiﬁed.

In the present paper we speciﬁcally consider the approach called generalized modeling. This approach constitutes a procedure by which the local dynamics in models can be analyzed in such a way that the results are almost always interpretable in the context of the application. Generalized modeling was originally developed for studying food chains [11, 9, 10] and was only later proposed as a general approach to nonlinear dynamical systems [12]. Subsequently, generalized modeling was used in systems biology, where it is sometimes called structural-kinetic modeling [34, 36, 45, 28] and is covered in recent reviews [33, 40, 17, 35, 29, 32]. In ecology generalized models have been employed in several recent studies [38, 14, 13, 42, 2], for instance for explorating the effects of food-quality on producer-grazer systems [38] and for identifying stabilizing factors in large food webs [14]. The latter work demonstrated that the approach of generalized modeling can be applied to large systems comprising 50 diﬀerent species and billions of food web topologies.

In the present paper we present a pedagogical introduction to generalized modeling and explain the underlying idea on a deeper level than previous publications. Furthermore, we propose some new techniques that considerably facilitate the formulation and analysis of generalized models. The approach is explained using a series of ecological examples of increasing complexity, including a simple model of omnivory that has so far not been analyzed by generalized modeling.

We start out in Sec. 2 with a brief introduction to fundamental concepts of dynamical systems theory. In Sec. 3, we introduce generalized modeling by considering the example of a single population. In contrast to previous generalized analyses of this system we use a shortcut that accelerates the formulation of generalized models. This shortcut is also used in Sec. 4, where we apply generalized modeling to a predator-prey system. Our ﬁnal example, shown in Sec. 5, is a simple omnivory scenario involving three species. This example already contains all of the diﬃculties that are also encountered in larger food webs.

**2 Local analysis of dynamical systems**

Generalized modeling builds on the tools of nonlinear dynamics and dynamical systems theory. Speciﬁcally, information is typically extracted from generalized models by a local bifurcation analysis. Mathematically speaking, a bifurcation is a qualitative transition in the long-term dynamics of the system, such as the transition from stationary (equilibrium) to oscillatory (cyclic) long-term dynamics. The corresponding critical parameter value at which the transition occurs is called the bifurcation point.

In this section we review the basic procedure for locating bifurcation points in systems of coupled ODEs. This analysis is central to the exploration of both generalized and conventional models and is also covered in many excellent text books, for instance [19, 15].

In the following we consider systems of N coupled equations

where |∗ indicates that the derivative is evaluated in the steady state.

Because the Jacobian is a real matrix, its eigenvalues are either real or form complex conjugate eigenvalue pairs. A given steady state is stable if all eigenvalues of the corresponding Jacobian J have negative real parts. When the function f (x) is changed continuously, for instance by a gradual change of parameters on which f (x) depends, the eigenvalues of the corresponding Jacobian change continuously as well.

Local bifurcations occur when a change in parameters causes one or more eigenvalues to cross the imaginary axis of the complex plane. In general, this happens in either of two scenarios: In the ﬁrst scenario, a real eigenvalue crosses the imaginary axis, causing a saddle-node bifurcation. In this bifurcation two steady states collide and annihilate each other. If the system was residing in one of the steady states before the transition, the variables typically change rapidly while the system approaches some other attractor. In ecology crossing a saddle-node bifurcation backwards can, for instance, mark the onset of an Allee eﬀect. In this case one of the two steady states emerging from the bifurcation is a stable equilibrium, whereas the other is an unstable saddle, which marks the tipping point between long-term persistence and extinction.

In the second scenario, a complex conjugate pair of eigenvalues crosses the imaginary axis, causing a Hopf bifurcation. In this bifurcation the steady state becomes unstable and either a stable limit cycle emerges (supercritical Hopf) or an unstable limit cycle vanishes (subcritical Hopf). The supercritical Hopf bifurcation marks a smooth transition from stationary to oscillatory dynamics. A famous example of this bifurcation in biology is found in the Rosenzweig-MacArthur model [31], where enrichment leads to destabilization of a steady state in a supercritical Hopf bifurcation. By contrast, the subcritical Hopf bifurcation is a catastrophic bifurcation after which the system departs rapidly from the neighborhood of the steady state.

In addition to the generic local bifurcation scenarios, discussed above, degenerate bifurcations can be observed if certain symmetries exist in the system. In many ecological models one such symmetry is related to the unconditional existence of a steady state at zero population densities. If a change of parameters causes another steady state to meet this extinct state, then the system generally undergoes a transcritical bifurcation in which the steady states cross and exchange their stability. The transcritical bifurcation is a degenerate form of the saddle-node bifurcation and is, like the saddle-node bifurcation, characterized by the existence of a zero eigenvalue of the Jacobian.

3 Density Dependent Growth of a Single Species In the following we demonstrate how the approach of generalized modeling can be used to ﬁnd local bifurcations in general ecological models. We start with the simplest example: the growth of a single population X. A generalized model describing this type

Let us now discuss the interpretation of the other two parameters sx and dx. For this purpose, note that these parameters are deﬁned as logarithmic derivatives of the original functions. Such parameters are also called elasticities, because they provide a nonlinear measure for the sensitivity of the function to variations in the argument.

For any power-law aX p the corresponding elasticity is p. For instance all constant functions have an elasticity of 0, all linear functions an elasticity of 1, and all quadratic functions have an elasticity of 2. This also extends to decreasing functions such as a/X for which the corresponding elasticity is -1. For more complex functions the value of the elasticity can depend on the location of the steady state. However, even in this case the interpretation of the elasticity is intuitive. For instance the Holling Type-II functional response is linear for low prey density and saturates for high prey density [16]. The corresponding elasticity is approximately 1 in the linear regime, but asymptotically decreases to 0 as the predation rate approaches saturation. A similar comparison for the Holling type-III function is shown in Fig. 1.

Elasticities are used in several scientiﬁc disciplines because they are directly interpretable and can be easily estimated from data [5]. In particular we emphasize that elasticities are deﬁned in the state that is observed in the system under consideration, and thus do not require reference to unnatural situations, such as half-maximum values or rates at saturation that often cannot be observed directly. We note that in previous publications the elasticities have sometimes been called exponent parameters and have been obtained by a normalization procedure. In comparison to this previous procedure the application of Eq. (6), proposed here, provides a signiﬁcant shortcut.

We now return to the discussion of the example system. So far we have managed to express the Jacobian determining the stability of all steady states by the three parameters α, sx, and dx. Because this simple example contains only one variable, the Jacobian is a 1-by-1 matrix. Therefore, the Jacobian has only one eigenvalue which is directly λ = α (sx − dx ). (12) The steady state under consideration is stable if λ 0, or equivalently

In words: In every system of the form of Eq. (4) a given steady state is stable whenever the elasticity of the mortality in the steady state exceeds the elasticity of reproduction.

A change in stability occurs when the elasticities of gain and loss become equal,

If this occurs the eigenvalue of the Jacobian vanishes and the system undergoes a saddle-node bifurcation.

For gaining a deeper understanding of how the generalized analysis relates to conventional models it is useful to consider a speciﬁc example. One model that immediately comes to mind is logistic growth, which is characterized by linear reproduction and quadratic mortality. However, based on our discussion above, it is immediately apparent that linear reproduction must correspond to sx = 1 and quadratic mortality to dx = 2. Without further analysis we can therefore say that steady states found for a single population under logistic growth must always be stable regardless of the other parameters.

Fig. 1 A. In the specifc example system a reporoduction rate, S(X) of the form of a Holling type-III functional response, aX 2 /(k2 +X 2 ), is assumed. This function starts out quadratically at low values of the population density X, but saturates as X increases. B. The corresponding elasticity, sx, is close to two near the quadratic regime χ = k/X ∗ ≈ 0, but approaches zero as saturation sets in.

A more interesting example is obtained when one assumes a reproduction rate following a Holling type-III kinetic and linear mortality,