How Levins’ dynamics emerges from a Ricker metapopulation model on the brink of extinction

F. Elías-Wolff, A. Eriksson, A. Manica B. Mehlig,

1 Department of Physics, University of Gothenburg, SE-41296 Gothenburg, Sweden

2 Department of Zoology, University of Cambridge, Cambridge, CB2 3EJ, UK

3 The Linnaeus Centre for Marine Evolutionary Biology, University of Gothenburg, SE-405 30 Gothenburg, Sweden

E-mail: Bernhard.Mehlig@physics.gu.se

## Abstract

Understanding the dynamics of metapopulations close to extinction is of vital importance for management. Levins-like models, in which local patches are treated as either occupied or empty, have been used extensively for this purpose, but they ignore the important role of local population dynamics. In this paper, we consider a stochastic metapopulation model where local populations follow a Ricker dynamics, and use this framework to investigate the behaviour of the metapopulation at the brink of extinction. As long as dispersal rates are not too large, the system is shown to have a time evolution consistent with Levins’ dynamics. We derive analytical expressions for the colonisation and extinction rates ( and ) in Levins-type models in terms of reproduction, survival, and dispersal parameters of the local populations, providing an avenue to parameterising Levins-like models from the type of information on local demography that is available for a number of species. To facilitate applying our results, we provide a numerical implementation for computing and .

Keywords: colonisation, extinction, Levins’ equation, Ricker model

## Introduction

Most species live in systems of populations connected by dispersal, also known as metapopulations. While small local populations can go extinct frequently (due to stochastic events, competition, predation and other factors), the whole metapopulation persists due to dispersal allowing for the recolonisation of empty patches and the rescuing of patches close to local extinction.

Levins [1] provided the first mathematical representation of metapopulations. His approach takes the simplifying assumption of treating individual patches as either occupied or empty, thus ignoring the local patch dynamics. The total fraction of occupied patches changes as a function of time, and is described by Levins’ equation:

(1) |

Here is the rate (per patch) of successful colonisation of empty patches, and is the corresponding extinction rate. Levins’ equation has formed the basis of most metapopulation modeling, and the basic model has been extended to account for a large number of additional factors, such as allowing for different patch sizes (so-called core-satellite models), Allee effects, and multi-species interactions.

Despite the extensive literature on the behaviour of Levins-like models (summarised e.g. by Etienne [2]), a fundamental question has received limited attention: under which conditions, that is to say for which type of local dynamics do real metapopulations exhibit a Levins-like behaviour? Keeling [3] investigated numerically a metapopulation dynamics based on local coupled maps, in order to determine whether Levins’ model is a good approximation. He found that for certain parameter values, the results of numerical simulations of the metapopulation fitted well to the form of Levins’ equation (similar results were obtained in [4]). Fronhofer et al. [5] have performed numerical simulations of a spatially structured population model, to answer the question under which conditions it behaves qualitatively as a metapopulation, Eq. (1), but make no attempt at determining the rates and . Fronhofer et al. [5] base their conclusions on a number of criteria (e.g. levels of patch occupancy and turnover) and conclude that on the brink of extinction the dynamics of their model is goverend by a Levins-like balance between patch colonisations and extinctions.

Even for conditions under which metapopulations can be shown to follow a Levins-like behaviour [6], it remains a challenge to parameterise the corresponding colonisation and extinction rates, as it requires long time series of censuses of a number of patches. While it is straightforward to generate such data from simulations, they are hard to obtain for empirical systems; indeed, this type of comprehensive data is rare in the ecological literature. See Ref. [6] for a study of an endangered butterfly population shown to behave Levins-like.

A possible solution to this problem is to develop an analytical formulation of the local stochastic dynamics within a metapopulation framework, thus providing a direct link between local demographic parameters and the colonisation and extinction rates of patches at the metapopulation level. Using such a framework, Eriksson et al. [7] estimated metapopulation extinction times for the case of local dynamics with time-continuous birth-death dynamics. In a related approach, Lande et al. [8] have proposed an elegant stochastic generalisation of Levins’ approach for continuous population dynamics, in order to quantify how the local patch dynamics affects the properties of the quasi-steady state of the metapopulation.

However many population models commonly used in the literature are based on discrete time evolution in terms of non-overlapping generations of individuals. In such models, the dynamics is usually represented by a collection of coupled maps [9, 10, 11]. Examples are the logistic map, the Ricker map [12], and the Hassell map [13].

In the following we develop an analytical formulation of a metapopulation where the local population dynamics is governed by a stochastic version of the Ricker map. We mention that the spatial structure of metapopulations with local Ricker dynamics was investigated in Refs. [14, 15, 16].

In this paper we investigate under which conditions a metapopulation model based on Ricker’s map exhibits a Levins-like behaviour. We then estimate patch colonisation and extinction rates (the two key parameters in Levins’ equation) in terms of reproduction, survival and dispersal parameters of the local populations.

The remainder of this paper is organised as follows. We first describe the metapopulation model. Then we derive a deterministic law for the evolution of this population, in the limit of many patches. We then proceed to show how this dynamics simplifies when the population is close to extinction, and how a discrete version of Levins’ equation (1) is obtained. Finally, in a concluding section, we discuss our results and put them in context. Mathematical details are deferred to appendices.

## Stochastic metapopulation model

Our model consists of a population distributed among patches evolving in discrete, non-overlapping generations. It is closely related to the model studied numerically by Keeling [3]. During each generation, three processes take place, in this order: reproduction, survival to maturity, and dispersal. The first two phases occur within patches, and are described by a stochastic Ricker map. The Ricker map as it is commonly used corresponds to the deterministic dynamics:

(2) |

Here denotes the number of adults in a local population in generation , the parameter characterises the average offspring size, that is, the fecundity rate. The parameter determines the probability of density-dependent survival (such as the incidence rate of cannibalism by adults).

Stochasticity is added in the following manner. Consider a patch with individuals in generation . We take the number of individuals after reproduction to be Poisson–distributed with mean . During the second phase, individuals survive to maturity independently with probability . Thus the number of adults in the next generation is a random variable with mean equal to the right-hand side of Eq. (2). During the third phase, the dispersal phase, individuals independently move away from the patch with probability . Finally, the migrants from all patches are gathered in a common pool from where the target patch of each migrant is independently and randomly chosen, and thus the migrants are uniformly randomly distributed among all patches. Fig. 1 illustrates this migration process. The pseudocode for the direct numerical simulations of the model described here is given in algorithm 1. Note that because we start from a Poisson distributed number of offspring and use binomial sampling in the survival and dispersal phases, it follows that the patch population sizes are Poisson distributed in each phase [17] (and that the population size is statistically independent between patches, conditional on the patch population sizes in the previous generation ).

## Dynamics in the limit of infinitely many patches

In this section we analyse the metapopulation dynamics in the limit of infinitely many patches. Since all patches are equivalent (the local population dynamics is given by the parameters and ), it is convenient to describe the state of the metapopulation by recording the frequencies of patches with a given number of individuals. In other words, the state of the metapopulation is described by the vector , where denotes the number of empty patches, the number of patches containing one individual, etc. In the limit of we expect the frequencies to approach -independent limits, and change deterministically from generation to :

(3) |

Here is the probability that a patch inhabited by individuals in generation has population size in generation . Since the fractions must sum to unity, it is convenient characterise the state of the metapopulation using the vector with components . Note that depends on . This dependence is determined by the migration scheme. In appendix A we derive the explicit form of . For the following discussion the details do not matter.

The long-term behaviour of Eq. (3) is decided by the values of the parameters , , and by the value of the emigration rate . For sufficiently large emigration rates, the metapopulation converges to a quasi-steady state, described by the steady-state solution of the dynamics (3). In a finite metapopulation, the state cannot persist ad infinitum since finite populations must eventually become extinct unless they continue to grow [18]. It is nevertheless of interest to study the limit because stable steady states in this limit correspond to very long-lived, quasi-steady states in finite but large populations.

As the emigration rate decreases (at fixed values of and ), the infinitely large metapopulation approaches the brink of extinction: there is a critical value where the infinite metapopulation becomes unstable and ceases to persist. For , the steady state is stable, thus extinction is certain.

How does depend on the life-history parameters and ? In order to answer this question we must investigate the linear stability of by linearising the dynamical equation (3). This procedure is briefly described in appendix B. The result of the analysis is shown in Fig. 2a for . We see that for small values of . The line separating stable from unstable steady states is shown as a solid line in the --plane in Fig. 2a. Above this line the non-trivial steady states are linearly stable (the infinite metapopulation persists ad infinitum), below this line they are unstable.

Panels b to d in Fig. 2 show a comparison between the dynamics obtained from Eq. (3), solid lines, and the results of direct stochastic simulations of the metapopulation model (symbols) for (panel b) and (panels c and d) patches. The dashed lines correspond to Levins’ dynamics (7), discussed in the following section.

## Levins’ equation

Points on the brink of extinction in Fig. 2a (solid line) correspond to parameter values where the metapopulation dynamics (3) suffers a qualitative change as mentioned in the preceding section. Slightly above this line, the metapopulation dynamics simplifies considerably, it becomes essentially one-dimensional [19, 20, 7]. In this section we describe the form of this one-dimensional equation, show that it corresponds to a time-discrete version of Levins’ equation (1), and determine the parameters and in this equation in terms of local demographic parameters.

Our results are valid for metapopulations close to (and above) the critical line in Fig. 2a. We use the parameter

(4) |

to characterise how close the population is to the brink of extinction (). For small values of the metapopulation dynamics (3) exhibits a so-called ‘slow mode’: a linear combination of the that changes slowly in time. The remaining ‘fast variables’ (different linear combinations of the ) rapidly relax to steady states that depend on the instantaneous value of the slow mode. We denote the slow mode by , it is given by a linear combination of the :

(5) |

In appendix C we show how to compute the coefficients . The details do not matter, but it is important to note the in the limit of large values of and small values of . In other words, when the local dynamics is much faster than the dispersal dynamics, the slow mode corresponds to the fraction of occupied patches

(6) |

See Fig. 3(b). In this limit we find the following equation for the change of from one generation to the next:

(7) |

This is a discrete version of Levins’ equation (1). The colonisation and extinction rates and in Eq. (7) are derived in appendix C. They depend on and as well as on :

(8) |

The coefficients and depend on the values of and . This dependence is illustrated in Fig. 3c. Eq. (8) show how the rates and in Levins’ equation (7) depend upon the life history of the local populations (determined by and in the case considered here, and upon the migration rate).

In panels b to d of Fig. 2 we compare Levins’ dynamics (described by Eqs. (7) and (8)), to results of direct simulations of the metapopulation model. As pointed out above, and as discussed in appendix C, we derived these equations assuming that the metapopulation is on the brink of extinction (solid line in panel a), and that the local dynamics is fast compared to migration. Both assumptions are satisfied in panel d, and we observe excellent agreement between Eqs. (7,8) and the results of direct stochastic simulations, except for a rapid initial transient. In panel b, the agreement is not as good. This is due to the fact that the parameter is too large, implying that the local and global time scales are not as clearly separated as in panel d. Panel c corresponds to parameter values much further from the critical line. In this case, Eqs. (7,8) fail to describe the metapopulation dynamics.

The dependence of the colonisation and extinction rates on the emigration rate is shown in Fig. 4a. Solid lines correspond to the analytical result (8) and (C.8,C.9) in appendix C. The results of direct stochastic simulations (symbols) were obtained as follows. First, for each set of parameters we determined the dynamics of by direct stochastic simulations. Second, assuming a clear time-scale separation between the local patch dynamics and migration, the resulting was projected according to . Third, the resulting values of were averaged over 100 independent simulations. Fourth, from this average, we computed as a function of , and found the coefficients of Eq. (7) by a least-squares polynomial fit to the parabola on the right-hand side of Eq. (7). Fits are shown in panel b.

As the migration rate approaches its critical value, Fig. 4a shows that the colonisation rate approaches the same value as the extinction rate in agreement with Eq. (8) for . For larger values of , the extinction rate is approximately independent of . This is a consequence of the fact that in the limit of time-scale separation (see Fig. 3c). The colonisation rate is seen to scale as , as predicted by Eq. (8). As the emigration rate increases further, we observe deviations between the theoretical prediction Eq. (8) and the numerical results. This is expected: the approximations leading to Eqs. (7,8) are valid on the brink of extinction and fail further away from the critical migration line.

## Discussion

In this paper we have analysed the dynamics of a stochastic metapopulation model based on coupled Ricker maps, and have shown how Levins’ dynamics emerges on the brink of extinction. Our derivation of a discrete version of Levins’ model is valid in the limit of many patches, and when the local dynamics is faster than the global migration dynamics. Levins’ model emerges on the brink of extinction because the deterministic dynamics (valid in the limit of many patches) exhibits a unique slow mode. Our finding that Levins’ equation is applicable at the brink of extinction, based on an analytical study of a stochastic metapopulation model, is in line with results obtained by other authors who have attempted to bridge the gap between the results of direct numerical simulations of individual-based stochastic metapopulation models and Levins’ model. Keeling [3] did so by fitting results from numerical simulations to the rates and in Levins’ model. Fronhofer et al. [5], also through numerical simulations, have explored under which circumstances individual-based models qualitatively behave as metapopulations. In this paper, we provide a convenient analytical framework that allows us to define formally under which circumstances Levins’ model provides an adequate approximation of the full metapopulation dynamics.

Another important outcome of our analytical treatment of stochastic metapopulations is our derivation of Levins’ key parameters and in terms of life history parameters from the Ricker dynamics that describes the local populations, see Eq. (8) in the main text and Eqs. (C.8), (C.9) in appendix C. A numerical implementation of our method for computing and is available in the supplementary material.

While Levins’ model has been used extensively by ecologists to explore the behaviour of metapopulations, it is very difficult to fit the model to ecological data. Ideally, one would need long time series tracking individual populations and recording extinction and colonisation. However, such datasets are rare, and unlikely to be available for species of management interest because of their risk of global extinction. Some of those rare examples are presented in Refs. [6, 21], where the authors have additionally fitted the ecological data in question to a metapopulation model in order to use it predictively with some success. In this paper, we have provided an alternative solution to this problem, as demographic parameters are often easily available from short term studies of individual populations, especially for larger mammals in which population turnover is likely to be low [22]. In addition, since the parameters in Levins’ equation depend upon the local population dynamics, as well as on the dispersal between patches, parameters obtained by fitting to a given data set do not necessarily transfer to other contexts. In contrast, our mechanistic connection between local population dynamics and the coefficients of Levins’ equation allows using a fitted model to investigate which aspects of the local population dynamics has the strongest effect on extinction and colonisation rates. For example, is the population most sensitive to adult or juvenile mortality, or is it mainly affected by fecundity?

We have focused on the Ricker map, which was originally devised for fish [12], but is also applicable to a variety of taxa, such as mammals [23], birds [24, 25], and insects [26, 27]. Furthermore, the methods used in this paper are not restricted to the Ricker map, but can be used to obtain estimates of and based on any other type of local stochastic population dynamics, as long as dispersal is still random and uniform (but the number of dispersing individuals can have any dependence on the local population size). The ability to fit Levins’ model to metapopulations at the brink of extinction should provide a potentially powerful tool for applied ecologists, who can now, with minimal information about the species of interest, leverage several decades of detailed analysis of the dynamics of Levins’ model.

## Acknowledgments

Financial support by Vetenskapsrådet, by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine, and by CeMEB are gratefully acknowledged. FEW wishes to thank the Centro Internacional de Ciencias in Cuernavaca for its hospitality during an extended stay.

## Appendix A

In this appendix we derive an expression for in Eq. (3). We calculate as follows. First consider reproduction and survival to adulthood. In our model the number of individuals after reproduction is Poisson distributed with mean . Survival is a binomial process with survival probability . Therefore the number of survivors is Poisson distributed with expected value . After dispersing individuals have left the patch (but before any dispersers arrive), the number of individuals is Poisson-distributed with expected value . Lastly, since dispersing individuals choose their destination independently, the number of individuals dispersing into a given patch is Poisson-distributed with expected value

(A.1) |

Bringing it all together, the probability in Eq. (3) of finding individuals in the patch in the next generation starting from individuals in the current generation is

(A.2) |

We took advantage of the fact that the number of resident individuals and the number of dispersing individuals arriving at the patch are two independent Poisson-distributed numbers (the sum of which is also Poisson distributed). Note that depends upon through the number of immigrants, . This non-linearity is necessary for a quasi-stable steady state to exist. Note also that the components are not independent, because the number of patches is equal to . We incorporate this constraint explicitly by substituting for . This results in the following dynamics for

(A.3) |

## Appendix B

In order to determine how depends upon the life history parameters and , we investigate the linear stability of the steady state by determining the eigenvalues of the stability matrix

(B.1) |

The state is linearly stable provided all eigenvalues are less than unity in modulus. The bifurcation in the dynamics, that is, the point where the extinction steady state changes from being unstable to being stable, occurs when the leading eigenvalue (that is, the largest in absolute value) approaches unity: as . This criterion allows us to compute the dependence of upon . We do this numerically in the following manner. At the bifurcation we have and the matrix , the elements of which is given in appendix D. We then numerically compute the eigenvalues of for a suitable range of , thus constructing a function . Finally, from this function we interpolate the value of such that .

## Appendix C

In this section we show how Levins’ dynamics Eqs. (7,8) emerges from the deterministic dynamics (A.3). In order to show how the dynamics in Eq. (A.3) simplifies for small positive values of , we follow Refs. [19, 20, 7] and expand the right-hand side of Eq. (A.3) in powers of , noting that the components of are of order . To simplify the notation we define and drop the generation dependence :

(C.1) |

Here are the elements of the stability matrix evaluated at . The coefficients , , , and are given in appendix D.

In order to identify the slow mode, we diagonalise the linearisation of the deterministic dynamics (A.3) at :

(C.2) |

We take the left and right eigenvectors to be bi-orthogonal and ordered as

(C.3) |

The condition (C.3) does not determine the normalisation of the eigenvectors. Fig. 3(a) shows that the components approach a -independent constant. We take this constant to be unity, thus fixing the normalisation of the leading left eigenvector, and thus by orthonormality fixing also the normalisation of . The limiting value is approached more rapidly for larger values of and smaller values of .

At we have that . For small values of we find that is real and slightly larger than unity, , and for . Thus, for small values of , perturbations away from grow slowly in the direction of , and decay rapidly in the directions , . Thus, the slow mode is given by

(C.4) |

This equation implies that the slow mode corresponds to the fraction of occupied patches if . Figure 3a shows the components of for different values of . We adopted a normalisation convention leading to as becomes large and tends to zero. This ensures that can be interpreted as the fraction of occupied patches in this limit as shown in Fig. 3c. The fast variables are given by

(C.5) |

The fast variables quickly approach local equilibria that depend on the instantaneous value of the slow variable, . In order to derive an equation for the dynamics of we start from (Appendix C). Inserting (C.4) and (C.5) into (Appendix C) we find an expansion of in powers of . We expect that is of order . Let us first consider this expansion to second order in :

(C.6) |

The fast variables quickly approach local equilibria given by . This yields

(C.7) |

On the right-hand side of Eq. (Appendix C), terms involving the fast variables are neglected. This is consistent since Eq. (C.7) implies that are of order close to the critical line. The coefficients and are given by

(C.8) | ||||

(C.9) |

Figure 3c shows how and depend on and .

We now make a further simplification, assuming time-scale separation between the fast local processes and slow emigration. In keeping with this approximation we neglect terms of order and higher. To orders and we find from Eqs. (Appendix C), (C.4), and the results of appendix D that

(C.10) |

Eq. (C.10) has the form of Levins’ equation (for discrete time evolution):

(C.11) |

The colonisation and extinction rates are given by:

(C.12) |

We note that is positive while is negative. Moreover, as becomes large and , the coefficient is approximately equal to . This implies that the coefficient is approximately independent of , while scales as , see Fig. 4a in the main text. A numerical implementation of our method for computing and is available as a MATLAB script in the supplementary material. The fixed points of Eq. (C.11) are the extinction point and the steady-state . For and , the extinction point is stable and the steady state is unstable. At , both fixed points meet and exchange their stability for .

At a period-doubling bifurcation occurs [28], and as increases still further, chaos ensues. See [14] for the study of chaotic transients. However, this behaviour is outside the range of applicability of Eq. (C.11), that is, too far from the critical line in Fig. 2a.

We conclude by summarising the three assumptions that were necessary to obtain Eqs. (C.11) and (C.12). First, it was assumed that the number of patches is large. Second, it was assumed that the metapopulation is on the brink of extinction (allowing us to assume that the parameter is small). Third, it was assumed that the local dynamics is much faster than the global dynamics (allowing us to assume that is small).

## Appendix D

In this section we give the coefficients occurring in the expansion (Appendix C): the matrix elements of and , as well as the coefficients , and . All coefficients are evaluated at , where the steady state is . The stability matrix has elements

for , | (D.1) |

with . The elements of the matrix are given by

for . | (D.2) |

The coefficients are given by

for . | (D.3) |

The coefficients are given by

for . | (D.4) |

Finally, the coefficients are given by

for . | (D.5) |

## References

- 1. Levins R (1969) Some demographic and genetic consequences of environmental heterogeneity for biological control. Bull. Entomol. Soc. Am. 15: 237–240.
- 2. Etienne R (2002) A scrutiny of the Levins metapopulation model. Comments Theor. Biol. 7: 257–281.
- 3. Keeling MJ (2002) Using individual-based simulations to test the Levins metapopulation paradigm. J. Animal Ecol. 71: 270–279.
- 4. Wennberg B, Hårding K, Ripa J (2013). Unpublished.
- 5. Fronhofer E, Kubisch A, Hilker F, Hovestadt T, Poethke H (2012) Why are metapopulations so rare? Ecology 93: 1967–1978.
- 6. Hanski I, Pakkala T, Kuussaari M, Lei G (1995) Metapopulation persistence of an endangered butterfly in a fragmented landscape. Oikos 72: 21–28.
- 7. Eriksson A, Elías-Wolff F, Mehlig B (2013) Metapopulation dynamics on the brink of extinction. Theor. Pop. Biol. 83: 101–122.
- 8. Lande R, Engen S, Sæther BE (1998) Extinction times in finite metapopulation models with stochastic local dynamics. Oikos 83: 383–398.
- 9. Ylikarjula J, Alaja S, Laakso J, Tesar D (2000) Effects of patch number and dispersal patterns on population dynamics and synchrony. J. Theor. Biol. 207: 377–387.
- 10. Ranta E, Kaitala V, Lindström J, Lindén H (1995) Synchrony in population dynamics. Proc. Roy. Soc. London B 262: 113–118.
- 11. Wysham DB, Hastings A (2008) Sudden shifts in ecological systems: intermittency and transients in the coupled Ricker population model. Bull. Math. Biol. 70.
- 12. Ricker WE (1954) Stock and recruitment. J. Fish. Res. Board Can. 11: 559–623.
- 13. Hassell MP (1975) Density-dependence in single-species populations. J. Animal Ecol. 44: 283–295.
- 14. Hastings A, Higgins K (1994) Persistence of transients in spatially structured ecological models. Science 263: 1133–1136.
- 15. Ruxton GD, Doebeli M (1996) Spatial self-organization and persistence of transients in a metapopulation model. Proc. Roy. Soc. London B 263: 1153–1158.
- 16. Dey S, Dabholkar S, Joshi A (2006) The effect of migration on metapopulation stability is qualitatively unaffected by demographic and spatial heterogeneity. J. Theor. Biol. 238: 78–84.
- 17. Karlin S, Taylor AD (1998) An introduction to Stochastic Modeling. San Diego, California: Academic Press.
- 18. Jagers P (1992) Stabilities and instabilities in population dynamics. J. Appl. Probab. 29: 770–780.
- 19. Guckenheimer J, Holmes P (1983) Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Appl. Math. Sci. New York: Springer, second edition.
- 20. Dykman MI, Mori E, Ross J, Hunt PM (1994) Large fluctuations and optimal paths in chemical kinetics. J. Chem. Phys. 100: 5735–5750.
- 21. Hanski I (1994) A practical model of metapopulation dynamics. J. Animal Ecol. 63: 151–162.
- 22. Elmhagen B, Angerbjörn A (2012) The applicability of metapopulation theory to large mammals. OIKOS 94: 89–10.
- 23. Sherratt T, Lambin X, Petty S, MacKinnon J, Coles C, et al. (2000) Use of coupled oscillator models to understand synchrony and travelling waves in populations of the field vole microtus agrestis in Northern England. J. Appl. Ecol. 37: 148–158.
- 24. Holyoak M, Baillie S (1996) Factors influencing detection of density dependence in British birds. Oecologia 108: 47–53.
- 25. Wiklund C (2001) Food as a mechanism of density-dependent regulation of breeding numbers in the merlin, falco columbarius. Ecology 82: 860–867.
- 26. Barlow N, Beggs J, Barron M (2002) Dynamics of common wasps in New Zealand beech forests: a model with density dependence and weather. J. Animal Ecol. 71: 663–671.
- 27. Melbourne B, Hastings A (2008) Extinction risk depends strongly on factors contributing to stochasticity. Nature 454: 100–143.
- 28. Ott E (2002) Chaos in dynamical systems. Cambridge: Cambridge Univ. Press, second edition.