## Abstract

Immigration can rescue local populations from extinction, helping to stabilize a metapopulation. Local population dynamics is important for determining the strength of this rescue effect, but the mechanistic link between local demographic parameters and the rescue effect at the metapopulation level has received very little attention by modellers. We develop an analytical framework that allows us to describe the emergence of the rescue effect from interacting local stochastic dynamics. We show this framework to be applicable to a wide range of spatial scales, providing a powerful and convenient alternative to individual-based models for making predictions concerning the fate of metapopulations. We show that the rescue effect plays an important role in minimizing the increase in local extinction probability associated with high demographic stochasticity, but its role is more limited in the case of high local environmental stochasticity of recruitment or survival. While most models postulate the rescue effect, our framework provides an explicit mechanistic link between local dynamics and the emergence of the rescue effect, and more generally the stability of the whole metapopulation.

## 1. Introduction

The increasing rate of habitat destruction and fragmentation poses a major threat to the persistence of many species [1,2]. Connectivity among populations has long been recognized to play a key role in preventing extinction [3,4]: immigration can lead to the recolonization of previously extinct patches, thus promoting the long-term persistence of the network of populations (i.e. the metapopulation) [3]. Also, in the absence of population turnover, immigration can reduce the probability of local extinction by boosting numbers of inhabitants under adverse conditions, thus buffering local populations against demographic and environmental stochasticity. This reduction in local extinction probability was termed the ‘rescue effect’ by Brown & Kodric-Brown [5], and ecologists have accumulated a wealth of empirical evidence of its potential importance (e.g. [1,2,6,7]).

While it is intuitively obvious that the rescue effect can play a pivotal role in determining the persistence of local populations as well as of the entire metapopulation, quantifying the link between local dynamics and the strength of the rescue effect is a challenging endeavour. Levins’s classical metapopulation model [8] considers patches as either occupied or empty, thus ignoring local population dynamics within patches. In an attempt to integrate the rescue effect in Levins’s model, several authors have assumed that the extinction probability of a given patch is a function of the number of occupied patches [9,10]. However, the shape of this function cannot be derived from first principles and has to be fitted to the data for each species of interest. Thus, the lessons learnt from this approach are specific to the dataset in question.

In order to investigate the factors that affect the strength of the rescue effect, it is necessary to explicitly model the local stochastic population dynamics. Agent-based models can be used for this purpose (e.g. [11]), but it is often difficult to understand the mechanisms that determine the global metapopulation dynamics, and to quantify the influence of the different factors characterizing the local population demography.

An analytical treatment of the metapopulation is much more desirable and informative, but it is not trivial. Lande *et al*. [12] were the first to develop such an analytical model based on a diffusion approximation for the local population dynamics. Their approach characterizes the metapopulation dynamics as a sequence of isolated, local extinction and recolonization events, and relies on the assumption of very low migration rates. This assumption is appropriate when a large fraction of suitable patches is not occupied near demographic equilibrium. This scenario is at odds with the observation that the majority of suitable patches tend to remain occupied in many metapopulations (e.g. [13–17]). Keeling [18] took a different approach, relying on a moment-closure technique, to investigate the persistence of stochastic metapopulations. This approach relies on knowing the distribution of population sizes, a problematic restriction if we are interested in investigating the dynamics of a metapopulation close to extinction, for which the shape of the distribution is likely to differ from the corresponding shape near demographic equilibrium.

In this paper, we develop an analytical framework, based on the identification of the slow and fast modes underlying the metapopulation dynamics. Our approach allows us to investigate the link between local stochastic dynamics and the strength of the rescue effect. We compare our results with spatially explicit simulations to explore the spatial scales at which our analytical results are relevant. While for this paper we focus on a model that describes the local dynamics in terms of a Ricker map [19], the approach can be applied to other models for population dynamics, formulated in terms of discrete non-overlapping generations or in terms of continuously occurring birth and death processes.

## 2. Metapopulation model

Our model considers *N* patches connected by global dispersal, and we specifically focus on metapopulations consisting of a large number of patches. Each patch can be inhabited by a population that evolves in discrete, non-overlapping generations. Each generation has two phases, the first corresponding to local recruitment and the second to dispersal to other populations. For biological realism, we use the model by Melbourne & Hastings [20] to represent local recruitment and density-dependent survival. This model is a stochastic version of the classical deterministic Ricker model [19] that incorporates the effect of both demographic and environmental stochasticity within populations. In our model, each individual has a Poisson-distributed number of progeny, with mean given by the recruitment rate *R*_{i} for individual *i*. The individual recruitment rate *R*_{i} is itself random (independent between generations), giving rise to demographic and environmental recruitment stochasticity. The environmental recruitment stochasticity determines a mean (environmental) recruitment rate, *R*_{E}, for the local population (*R*_{E} is gamma-distributed with mean *R* and shape parameter *k*_{E}). Given the environmental state of the local population, individual recruitment rates *R*_{i} are independent gamma-distributed stochastic variables with mean *R*_{E} and shape parameter *k*_{D} (giving rise to demographic recruitment stochasticity). Finally, each progeny survives (independently) to adulthood with density-dependent probability exp(-*αn*), where *n* is the number of adults in the population and *α* regulates the strength of the density dependence. As shown by Melbourne & Hastings [20], the expected number of progeny surviving to adulthood is the same as in the classical Ricker dynamics (with recruitment *R* and density dependence *α*), but the variance of the distribution is strongly affected by the amount of demographic and environmental recruitment stochasticity (defined by the shape parameters *k*_{E} and *k*_{D}). Putting it all together (following [20]), the probability of having, after the first phase, *j* individuals in a population with *i* adults (*P _{ij}*) can be written as
2.1a

The original formulation by Melbourne & Hastings [20] detailed above describes a system where environmental stochasticity acts on recruitment. We also consider an alternative model in which demographic stochasticity of recruitment is the same as in equation (2.1*a*), but where environmental stochasticity affects survival. Survival stochasticity is modelled by dividing *α* by a gamma-distributed variable with unit mean and shape parameter *k*_{A}. Thus, for this model we have
2.1b

In the following, the derivations do not rely on the specific form of *P _{ij}* in equations (2.1

*a*) and (2.1

*b*), and we will therefore use ‘equation (2.1)’ when either equation (2.1

*a*) or (2.1

*b*) can be used interchangeably.

In the second phase, the dispersal phase, each surviving progeny disperses with probability *m*, and otherwise stays in the local population. For simplicity, we assume that a dispersing individual is equally likely to move to any other patch; if the destination patch is currently unoccupied, the individual colonizes it and starts a new local population. We later investigate the effect of the range of dispersal using spatially explicit simulations of the dynamics.

As a result of demographic and environmental stochasticity, individual populations undergo substantial fluctuations in size (figure 1). When the number of populations is large, by contrast, we expect that fluctuations in local populations have little effect upon the global distribution of population sizes in the metapopulation. In this limit of many patches, we formulate the metapopulation dynamics in terms of the fraction *f _{i}* of populations with

*i*individuals. The fraction

*f*after local recruitment and density dependence (denoted as ) can be written as 2.2where

_{i}*P*is given by equation (2.1). After dispersing individuals have left their natal populations, but before they have arrived to their new populations, the fraction is given by 2.3

_{ij}Finally, as there is no mortality during dispersal (mortality only occurs within the natal patch in our model), the mean number of individuals arriving to each patch is equal to the mean number of individuals dispersing out of each patch. Because of the large number of patches, the distribution of individuals arriving to each patch is approximately Poisson-distributed, allowing us to write the dispersal rate (i.e. the expected number of individuals arriving to a patch), *I*, as
2.4

We obtain the population size distribution for the next generation as 2.5

Equations (2.1)–(2.5) can be written as a single equation for the dynamics in matrix form 2.6

where ** f** is the vector of population size frequencies (

*f*) and

_{i}**A**is a matrix that depends on the dispersal rate

*I*(

*t*).

Note that equations (2.1)–(2.4) describe processes that occur independently in each patch; only in the dispersal stage (equation (2.5)) have we taken advantage of the assumption that the number of patches is large and treated *f _{i}*(

*t*) as deterministic rather than stochastic. As we will show later using simulations, this approximation is well justified. We also point out that, because

*I*(

*t*) depends upon

*f*(

_{i}*t*) in equation (2.4), the left-hand side of equation (2.6), Δ

**(**

*f**t*), depends nonlinearly on the distribution of population sizes in generation

*t*. This nonlinearity is important in determining how we can analyse the behaviour of the system in §3.

### (a) Spatial scale relevant for our spatially implicit model

To test the importance of our assumption of a common dispersal pool, we compared the predictions of our spatially implicit model with spatially explicit simulations of the stochastic model on a two-dimensional grid of local populations connected by dispersal, varying the range of dispersal to explore the spatial scale at which the analytical framework is relevant (see the electronic supplementary material, appendix A for details). We quantified the behaviour of the metapopulation by estimating steady-state mean population sizes for a range of dispersal probabilities. The results of our analytical framework were in excellent agreement with the simulations (figure 2), except for dispersal restricted to the nearest neighbour population in the lattice (yellow triangles in figure 2). Even at this most local form of dispersal, the analytical framework was able to capture qualitatively the relationship between our observable (mean population size) and the probability of dispersal *m*.

## 3. Dynamical analysis

### (a) Identifying the slow mode underlying the metapopulation dynamics

As mentioned above, individual populations exhibit substantial population size fluctuations (figure 1). But the distribution of population sizes across the whole metapopulation (graphs at top of figure 1) changes relatively slowly through time, because demographic events in individual populations have very little effect as long as the number of populations is large, indicating a separation of time scale between local and global dynamics, and the existence of a slow mode. We also find that when the metapopulation recovers from different population sizes, it retraces the same trajectory to the global stationary state (except for brief transients, lasting a few generations after the perturbation; see the electronic supplementary material, figure S1). This separation of time scales is a common property of large interconnected nonlinear systems, and there are well-established methods for analysing the fast and slow components of such systems (e.g. [21]).

In general, fast and slow components of dynamical systems are found by analysing the eigenvectors and eigenvalues of the linearized dynamics. Eigenvalues of small magnitude correspond to slow components. In our case, the structure of the dynamics is determined by the dispersal rate *I* (equation (2.6)). Therefore, this quantity is a natural observable for parameterizing the slow mode. Dispersal between patches renders the metapopulation dynamics nonlinear (the *I*-dependence of **A** makes equation (2.6) nonlinear). Therefore, the eigenvectors and eigenvalues of the linearized dynamics (determined by the structure of the Jacobian) depend on the state of the metapopulation (e.g. the degree of occupancy, its proximity to extinction, etc.). To accommodate this important fact, we have to solve locally with respect to *I* for the fast and slow components; and determine how they interact (see the electronic supplementary material, appendix B).

Numerically, the slow mode corresponding to a given rate of dispersal *I* (i.e. the local solution) can be found by spectral decomposition of the linearized dynamics (see the electronic supplementary material, appendix B for details). The distribution of population sizes must (i) be normalized such that it sums to unity and (ii) have a rate of dispersal that matches the given value of *I* (two linear constraints). To uniquely determine the slow mode, we use the adiabatic approximation that the population size distribution is unchanged along the directions of the fast degrees of freedom (as determined from the spectral decomposition).

Thus, we can reconstruct the full distribution of population sizes as a function of *I*, and thereby also estimate Δ*I*, giving a closed dynamics for *I*. Figure 3*a* shows Δ*I* as a function of *I* from simulations (thin line) and from theory (thick line). In this example, several higher eigenmodes each give small but significant contributions to Δ*I* and are tied (‘slaved’) to the first eigenmode. By taking these higher eigenmodes into account, we obtained a very close match between the simulations and the theory. We note that the theory is exact at *I* = 0 (trivially) and in the stable steady state (see the electronic supplementary material, appendix B). In figure 3*b*, we show how the distribution of local population sizes evolves as the metapopulation recovers from initially very small local population sizes. The theory (lines) matches the distributions obtained from simulations (symbols) very well.

Thus, the adiabatic solution to the system provides us with the ability to predict the detailed state of the metapopulation based on the effect of dispersal on the local dynamics. From the predicted states, it is possible to estimate rate of local extinctions (in other words, the strength of the rescue effect). This approach has the important advantage of being derived from the full stochastic description of the system (in other words, the rescue effect emerges mechanistically from the interaction of the local stochastic populations).

### (b) Understanding the slow mode in biological terms

While the separation in fast and slow modes detailed in §3*a* is accurate and mathematically tractable, the formulation of slow mode reconstructed above remains difficult to interpret in biological terms. At the expense of losing some accuracy, we can express the dynamics in terms of the demographic response of a single patch to a constant influx of dispersing individuals. Rather than a full spectral analysis of the dynamics, we express the slow mode as a linear combination of eigenvectors of the matrix **A** and keep only the terms corresponding to the two eigenvalues closest to zero (one eigenvalue is zero; see the electronic supplementary material, appendix C for details). This is equivalent to assuming a separation of time scales in the local population dynamics [22]. The coefficients of the two eigenvectors are uniquely determined by the constraints that the distribution of patch sizes must sum to unity and that the rate of dispersal *I* must match a given value. Putting it all together, we find a closed equation for the change in dispersal rate, *I*, from one generation to the next:
3.1where *D*(*I*) is the average number of individuals dispersing from a population with a stable stream of *I* individuals (on average, Poisson-distributed) arriving to the population each generation, and *λ*_{2}(*I*) is the smallest non-zero eigenvalue of **A** in equation (2.6). The first factor in equation (3.1), *D*(*I*) − *I*, expresses the imbalance between dispersal from the population and immigration into the population at the local steady state. The second factor in equation (3.1), –*λ*_{2}(*I*), reflects the elasticity of Δ*I* with respect to this imbalance (i.e. how quickly the demographics of the local populations respond to the imbalance). We note that because single-population dynamics with immigration and dispersal is a Markov process with a unique steady state, the elasticity is always positive. Therefore, the dispersal rate (and hence the mean population size in the metapopulation) increases when the imbalance (the first factor) is positive and decreases when it is negative (reflecting whether patches are sources or sinks for the dynamics, on average). The magnitude of the change in dispersal rate is affected by both factors. Numerically, equation (3.1) gives solutions that are qualitatively similar to the slow-mode theory and the simulations (see electronic supplementary material, figure S2). The main difference between the simplified and full reconstruction is the magnitude of Δ*I*: especially, the sign of Δ*I* is still given by the first term in equation (3.1), the imbalance of dispersal and immigration in the local population model (because both theories are exact when Δ*I* = 0).

It is instructive to compare our approximate one-dimensional dynamics to the classical equations in Levins’s model. In the simplest version, Levins’s equation predicts a simple quadratic relationship between a variable (traditionally the fraction of occupied patches) and its change from one generation to the next
3.2Here, *p** is the fraction of occupied patches at the global steady state and *c* is the colonization probability of empty patches. Eriksson *et al.* [23] showed that under certain circumstances (specifically, when the metapopulation is on the brink of extinction) a Levins-like dynamics may emerge because the distribution of population sizes in occupied patches has a rigid shape, so that only the fraction of occupied patches *p* changes. In this case, the immigration rate *I* and the mean population size are both proportional to *p* (because all three quantities are linear combinations of the underlying distribution of population sizes), and therefore they too have Levins-like dynamics (but with different coefficients in equation (3.2)).

In the current model, by contrast, the shape of the distribution of population sizes changes markedly during recovery to the global steady state (figure 1). As expected, the fraction of occupied patches shows very strong deviations from Levins’s equation caused by the rescue effect (see electronic supplementary material, figure S3A, where the expectation from the Levins’s equation would be a constant line given by *c*). Surprisingly, however, we find that *I* and mean population size () both follow rather closely the predictions from Levins’s equation (see electronic supplementary material, figure S3B,C). In order to understand this apparent contradiction, we need to consider the dispersal–immigration imbalance (*D*(*I*) − *I*) and its elasticity (given by the factor −*λ*_{2}(*I*)) in equation (3.1). It turns out that the changes in elasticity for different values of *I* (which are not accounted for by Levins’s equation) are counterbalanced by corresponding changes in the dispersal–immigration imbalance (see electronic supplementary material, figure S3D). Repeating the analysis for the model with environmental stochasticity on survival gave qualitatively similar results (see electronic supplementary material, figure S4).

## 4. The rescue effect in action

We explored the behaviour of the metapopulation after a perturbation in the form of a reduction in the population size of all patches (red lines in figure 4*a*,*b*). Because environmental survival stochasticity gave qualitatively similar results, we focus on environmental recruitment stochasticity in figure 4. The fraction of occupied patches (figure 4*a*) displayed a non-monotonic transient during recovery, with an initial loss of a number of patches due to stochastic extinctions of the diminished populations followed by a progressive recolonization from the surviving patches. The average population size of occupied patches, on the other hand, recovered monotonically towards the equilibrium value (figure 4*b*). The behaviour was different when the system was perturbed by increasing population sizes in all patches. In this case, both the fraction of occupied patches and mean population size returned monotonically to their equilibrium values (blue lines in figure 4*a*,*b*).

Using our ability to predict the rate of local extinctions from the distribution of patch sizes (as detailed in the electronic supplementary material, appendix B), we explored the probability of extinction for any given patch following perturbations. This allows us to unravel the mechanisms that force the system to return to its steady state. In the absence of a rescue effect, there should be no relationship between local extinction and the fraction of occupied patches. However, we observed a clear inverse relationship between the probability of patch extinction and the fraction of occupied patches (figure 4*c*). We compared this with the relationship expected from the slow mode (thick shaded line in figure 4*c*) and found that the dynamics very quickly converges to this line, despite starting off the slow mode. We also found a steep inverse relationship between the probability of patch extinction and mean population size (figure 4*d*), and also for this relationship a quick convergence to the slow mode. These relationships reflect both the rescue effect and the increasing resilience to demographic stochasticity of larger populations. Thus, in our model, the rescue effect emerges mechanistically from the interaction of the stochastic dynamics of individual populations.

Compared with models in which the rescue effect is represented by an arbitrary function linking extinction rates with the proportion of occupied patches, the mechanistic nature of our framework allows us to investigate the relative roles of demographic and environmental stochasticity in affecting the strength of the rescue effect. We considered the effect of each type of stochasticity separately by setting the strength of the other kinds to zero (figure 5). As we were interested in the rescue effect, we conditioned our simulations to have the same mean within-patch occupancy (10 individuals, thus standardizing the metapopulation state). While increasing either type of stochasticity had a negative effect on patch persistence, environmental stochasticity had a much stronger impact on the rescue effect than demographic stochasticity (figure 5). This difference arises because environmental stochasticity tends to affect all individuals within the same patch simultaneously, thus giving little room for the rescue effect to come into play. On the other hand, demographic stochasticity will only affect part of the population, thus leaving small populations that can be rescued. Among the two types of environmental stochasticity, recruitment stochasticity had a larger negative effect than survival stochasticity. Interestingly, very low levels of environmental survival stochasticity led to a slight decrease in extinction probability (a phenomenon also observed by Chesson [24]).

## 5. Discussion

Our analytical approach provides a powerful way to explore the link between local and metapopulation dynamics. While the distribution of local population sizes in a metapopulation undergoes complex changes while the population approaches the global steady state, it was possible to find a close approximation to the full population dynamics using the instantaneous dispersal rate as the dynamic variable. Simulations showed that our analytical approach is perfectly adequate in describing metapopulations with localized dispersal, as long as the range of dispersal extends beyond the immediate neighbours. This robustness to the exact scale of dispersal makes our approach widely applicable, providing a convenient alternative to individual-based simulations when investigating metapopulations with explicit local dynamics. An important property highlighted by our framework is the emergence of the rescue effect from the interaction of local dynamics mediated via dispersal.

Our framework makes it possible to quantify the rescue effect, as determined by the local dynamics, in models where the strength of the effect has to be defined *a priori*. For example, when investigating extinction debt using the framework developed by Ovaskainen & Hanski [25], the rescue effect is represented by introducing nonlinearities in the colonization and extinction functions. The shape and scale of such nonlinear functions could be quantified only if several long time series collected under similar circumstances were available, a luxury that is usually not available. Within our framework, the same feat can be achieved simply by taking into consideration the local population dynamics, something that is often studied and quantified by ecologists. A further advantage is that our approach also allows us to investigate how changing the local dynamics might affect the whole metapopulation.

Local stochasticity, and in particular its environmental component, was shown to play an important role in determining the metapopulation dynamics, especially during the recovery from perturbation. Local environmental stochasticity not only had a direct effect on the local population size, but also had the potential to negate the rescue effect from the rest of the metapopulation. These results are in line with much of the literature on this topic [12,20,26–28] and stem from the fact that local environmental stochasticity affects all individuals within the same population simultaneously, and thus has a much stronger effect on the persistence time of the population compared with demographic stochasticity. Global environmentally driven fluctuations (affecting all populations simultaneously), which were not modelled in this paper, have the potential of also affecting recolonization rates because of synchrony within the metapopulation [26,28]. It will be interesting in the future to explore the link between local demography, global environmental stochasticity and metapopulation persistence time.

The key take-home message of our analysis is that local population dynamics can play a very important role in determining metapopulation behaviour, directly as well as via the rescue effect. We have provided a convenient analytical framework that allows this role to be investigated, and found our approach to be surprisingly robust to the spatial details to dispersal (except for very localized movement to adjacent patches). Interestingly, the dynamics of the metapopulation can be reconciled with a Levins-like dynamics with appropriately modified rates. This provides an opportunity to establish which local dynamics are compatible with previous models based on Levins’s equation and its generalizations.

## Funding statement

B.M. gratefully acknowledges financial support by Vetenskapsrådet, by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine and by the Centre for Evolutionary Marine Biology at Gothenburg University.

- Received November 28, 2013.
- Accepted January 14, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.