# Metapopulation dynamics on the brink of extinction

###### Abstract

We analyse metapopulation dynamics in terms of an individual-based, stochastic model of a finite metapopulation. We suggest a new approach, using the number of patches in the population as a large parameter. This approach does not require that the number of individuals per patch is large, neither is it necessary to assume a time-scale separation between local population dynamics and migration. Our approach makes it possible to accurately describe the dynamics of metapopulations consisting of many small patches. We focus on metapopulations on the brink of extinction. We estimate the time to extinction and describe the most likely path to extinction. We find that the logarithm of the time to extinction is proportional to the product of two vectors, a vector characterising the distribution of patch population sizes in the quasi-steady state, and a vector – related to Fisher’s reproduction vector – that quantifies the sensitivity of the quasi-steady state distribution to demographic fluctuations. We compare our analytical results to stochastic simulations of the model, and discuss the range of validity of the analytical expressions. By identifying fast and slow degrees of freedom in the metapopulation dynamics, we show that the dynamics of large metapopulations close to extinction is approximately described by a deterministic equation originally proposed by Levins (1969). We were able to compute the rates in Levins’ equation in terms of the parameters of our stochastic, individual-based model. It turns out, however, that the interpretation of the dynamical variable depends strongly on the intrinsic growth rate and carrying capacity of the patches. Only when the growth rate and the carrying capacity are large does the slow variable correspond to the number of patches, as envisaged by Levins. Last but not least, we discuss how our findings relate to other, widely used metapopulation models.

###### keywords:

Metapopulations, stochastic dynamics, extinction, migration^{†}

^{†}journal: Theoretical Population Biology

## 1 Introduction

The habitats of animal populations are often geographically divided into many small patches, either because of human interference or because natural habitats are patchy. Understanding the dynamics of such populations is a problem of great theoretical and practical interest. Isolated small patches are often extinction prone, for example because of inbreeding in combination with demographic and environmental stochasticity (Hanski, 1999). When the patches are connected into a network by migration, local populations may still be prone to extinction, but the whole population may persist because empty patches are re-colonised by migrants from surrounding occupied patches. Moreover, when the migration rate is sufficiently large, local populations can be stabilised by an inflow of immigrants (the rescue effect). Given a model of such a population, the main concern is usually to find out how the different parameter values in the model affect the growth and persistence of the metapopulation.

In Levins’ model, the metapopulation dynamics is simplified by treating each patch as either occupied or empty (Levins, 1969). The rates of patch colonisation and extinction are expressed as functions of the total fraction of occupied patches in the population. In its simplest form Levins’ model describes the time change of the fraction of occupied patches:

(1) |

Here is the rate of successful colonisation of an empty patch, and is the rate at which occupied patches turn empty (the rate of patch extinction). How these rates are related to the life-history parameters determining the stochastic, individual-based metapopulation dynamics is not explicitly known. In spatially explicit models the parameters are chosen to be functions of for example the number of occupied neighbouring patches (Hui and Li, 2004; Roy et al., 2008). Eq. (1) is generally motivated by assuming a separation of time scales between colonisation and extinction of patches on the one hand, and the population dynamics within patches on the other hand (Levins, 1969; Hanski and Gyllenberg, 1993; Lande et al., 1998; Etienne, 2002). Many of the mechanisms that can be observed in natural populations (e.g. Allee and rescue effects) can be represented by suitable modifications of Eq. (1), and the qualitative behaviour of such models is well understood (Hanski and Ovaskainen, 2000; Etienne, 2000; Harding and McNamara, 2002; Zhou and Wang, 2004). See also Zhou et al. (2004); Taylor and Hastings (2005), and Martcheva and Bolker (2007).

Models that describe the metapopulation dynamics in terms of the fraction of occupied patches have no explicit connection to the population dynamics within each patch. In other words, it is not apparent how individual births, deaths, and migration events translate into colonisation and extinction rates (Etienne, 2002).

An alternative line of analysis focuses on the population dynamics within a single patch (Drechsler and Wissel, 1997; Lopez and Pfister, 2001; Newman et al., 2004). The advantage of this approach is that birth and death processes, emigration, and immigration can be modelled explicitly in terms of the number of individuals in the patch. This makes the model much more immediate in terms of biological processes (such as density dependence and Allee effects), but the rest of the metapopulation is reduced to a background source of immigrants (assumed to be stationary) and is not explicitly modelled.

Several authors have attempted to bridge the gap between detailed local population dynamics and the dynamics at the overall population level. Keeling (2002) estimated rates of colonisation and extinction events in Levins’ equation from individual-based simulations. This method rests on the assumption that Eq. (1) is an appropriate description of the metapopulation dynamics.

Higgins (2009) has investigated the extinction risk of metapopulations subject to strong dispersal, that is, in a limit where migration is faster than the local population dynamics within patches. As pointed out above, the converse is commonly assumed in writing Eq. (1). Within his model, Higgins (2009) addressed the question of how the degree of fragmentation of the metapopulation affects its risk of extinction.

Chesson (1981, 1984), Hanski and Gyllenberg (1993), Casagrandi and Gatto (1999, 2002), Nachman (2000), Metz and Gyllenberg (2001) and others have analysed the equilibrium states and possible persistence of metapopulations in models consisting of infinitely many patches with local dynamics coupled via a migration pool (reviewed in Massol et al., 2009). The assumption that there are infinitely many patches is crucial in these studies, as it makes it possible to pose the question under which circumstances the metapopulation persists ad infinitum, that is, for which choice of parameters the metapopulation dynamics reaches a non-trivial stable equilibrium. For finite metapopulations, by contrast, there is no stable equilibrium corresponding to persistence (for a review of stochastic extinctions in biological populations, see e.g. Ovaskainen and Meerson, 2010). As is well known, finite populations must eventually become extinct (unless they continue to grow). This fact is referred to as the ‘merciless dichotomy of population dynamics’ by Jagers (1992). Moreover, it remains unclear how the equations employed in these studies are related to Eq. (1).

In this paper we characterise the dynamics of metapopulations with a large, but finite, number of patches. We derive, from first principles, how the stochastic metapopulation dynamics determines the distribution of individuals over patches. The only critical assumption in our derivation is that the number of patches is large enough. Especially, it is not necessary to assume that the typical number of individuals per patch is large. Neither is it necessary to assume a time-scale separation between local and migration dynamics. We show that our results represent the typical transient of the metapopulation towards a quasi-steady state (if such a state exists).

On the brink of extinction, that is, close to the bifurcation point where an infinite metapopulation ceases to persist, we use a systematic expansion of an exact master equation in powers of (where is the number of patches) to find the most likely path to extinction, as well as the leading contribution to the time to extinction. In this case (close to the bifurcation) the metapopulation dynamics simplifies considerably: it can be approximated by a simple, one-dimensional dynamics. This fact is a consequence of a general principle (Guckenheimer and Holmes, 1983) stating that a dynamical system close to a bifurcation exhibits a ‘slow mode’: a particular linear combination of the dynamical variables is found to relax slowly, and the remaining degrees of freedom relax much more quickly and may be assumed to be in local equilibria. In other words, when the dynamical system is in a perturbed state, the slow mode evolves towards the equilibrium state on a longer time-scale than the fast variables do. This renders the dynamics effectively one-dimensional. We find the slow mode of the metapopulation dynamics and show how it depends on the properties of the local dynamics (given by the local growth rate and the local carrying capacity). The slow mode determines the stochastic dynamics of finite metapopulations as well as the deterministic dynamics of metapopulations consisting of infinitely many patches. In the latter case we find that the slow mode obeys Eq. (1). We derive how the parameters and depend on the parameters determining the life history of the local populations. We show, however, that the variable is in general not given by the fraction of occupied patches as envisaged by Levins. But it turns out that approaches the fraction of occupied patches in the limit of large carrying capacities and large local growth rates (this is the limit of time-scale separation mentioned above).

In other words, we have derived Levins’ model, Eq. (1), from a stochastic, individual-based model of a finite metapopulation. We find that Eq. (1) is still valid (close to the bifurcation) even when there is no time-scale separation between the local and the migration dynamics. This is the consequence of the existence of a slow mode.

Lande et al. (1998) have suggested an elegant stochastic generalisation of Levins’ model, in an attempt to compute how the local patch dynamics affects the properties of the quasi-steady state of the metapopulation (and the average time to extinction of this population). Their main idea is to connect the extinction and colonisation rates in Eq. (1) to local processes. The extinction rate is calculated as the inverse expected time to extinction of a single patch at the carrying capacity (which is determined self-consistently), and colonisation is defined as the rate of a single migrant arriving to an empty patch, seeding a population that grows to the carrying capacity. This scheme is persuasive but not rigorous. The predictions of Lande et al. (1998) have, to our knowledge, never been tested by comparisons to results of simulations of stochastic, individual-based metapopulation models. Therefore it is important that our approach allows us to compute the rates of extinction and colonisation from first principles. In the limit of large patch carrying capacities and close to the bifurcation we obtain expressions (exact in the limit we consider) that are very similar, but not identical, to the relations proposed by Lande et al. (1998).

In summary, we characterise the stochastic dynamics of metapopulations on the brink of extinction. Using an expansion of the exact master equation describing the stochastic dynamics with the number of patches as a large parameter, we identify a slow mode in the metapopulation dynamics, regardless of whether there is a time-scale separation between local and global dynamics or not. We show under which circumstances widely used metapopulation models provide accurate descriptions of metapopulation dynamics.

The remainder of this paper is organised as follows. In section 2 we describe the individual-based stochastic metapopulation model investigated here. We summarise how our numerical experiments were performed and briefly describe how to represent the metapopulation dynamics in terms of a master equation, and how to expand this equation in powers of . Our results are described in section 3. We first discuss the limit , demonstrate under which circumstances metapopulations persist in this limit, and analyse the metapopulation dynamics. Second, we turn to stochastic fluctuations of finite metapopulations, and summarise our results on fluctuations in the quasi-steady state, the most likely path to extinction as well as the average time to extinction. Section 4 contains our conclusions. Appendices A, B, and C summarise details of our calculations.

## 2 Methods

In this section we define the stochastic, individual-based metapopulation model that is analysed in this paper. We describe how our numerical experiments are performed and derive a master equation, Eq. (2.3), that is the starting point for our mathematical analysis of metapopulation dynamics. In general it is not possible to solve this equation in closed form. We demonstrate how an approximate solution can be obtained by expanding the master equation using the number of patches as a large parameter. In the limit of , the dynamics reduces to a deterministic model. The equilibrium properties of this model were analysed by Casagrandi and Gatto (1999, 2002), and by Nachman (2000).

### 2.1 Stochastic, individual-based metapopulation model

The model consists of a population distributed amongst patches, as illustrated in Fig. 1. In each patch, the local population dynamics is a birth-death process with birth rates , and death rates . Here denotes the population size in a given patch (and ). To simplify the discussion we assume that the rates are the same for all patches, but more general cases can be treated within the approach described in this paper. In the following we illustrate our results for a particular choice of birth and death rates:

(2) | ||||||

(3) |

The parameter is the birth rate per individual, is the density-independent per capita mortality, and determines the carrying capacity of a single patch. For simplicity, we take hereafter (this corresponds to measuring time in units of the expected life-time of individuals in the absence of density dependence). The parameters occurring in Eqs. (2,3) are listed in Table 1, which summarises the notation used in this article.

In addition to the local population dynamics, the number of individuals in each patch can change because some individuals emigrate from their patch to other patches, or because immigrants arrive from other patches.

Following Hanski and Gyllenberg (1993), the migration process is modelled as follows: individuals emigrate from a patch at rate , where is the population size in the patch in question (). If the individuals migrate independently and with constant rates, is proportional to :

(4) |

More complex migration patterns can be incorporated. If for example individuals moved to avoid overcrowding, the emigration rate would be density dependent.

The emigrants enter a common dispersal pool, containing the migrants from all patches that have not yet reached their target patch. Each migrant stays an exponentially distributed time in the pool, with expected value , before reaching the target habitat, which is chosen with equal probability among all patches. This process is illustrated in Fig. 1. Migration may fail if the individuals die before reaching the new habitat; this is modelled by the rate of dying during dispersal. In practice, the probability of successful migration depends on the background mortality of the individuals, on the time the migrating individual spends in the dispersal pool, and on additional perils individuals may be exposed to during dispersal (e.g. increased risk of predation due to lack of cover, etc.). In summary, if there are migrants, individuals leave the dispersal pool per unit of time, and the rate of immigration to a given patch is .

### 2.2 Numerical experiments

In the direct numerical simulations, the population evolves in the following manner. First, at any given time, the local rates, Eqs. (2-4), sum to a rate of the next locally generated event, , for patch . This rate is simply the sum of the local rates, Eqs. (2-4). If patch contains individuals then we have:

(5) |

The sum of these rates over all patches, , for , yields the rate for the next event occurring in the population. Thus the simulation proceeds by generating an exponentially distributed random number with expected value . Second, at each time step, a patch is chosen with probability . Third, the type of event is chosen randomly: birth with probability , death with probability or emigration with . Fourth, the numerical experiments described below were performed in the limit and for . This means that emigrating individuals are immediately assigned to a randomly chosen patch (possibly the one they come from).

### 2.3 Master equation

A natural and commonly adopted approach to describe stochastic population dynamics is to derive a master equation (van Kampen, 1981) for the change in time of the probability of observing individuals in the first patch, individuals in the second patch, and of observing individuals in the migrant pool. Recently this approach was adopted by Meerson and Sasorov (2011) to describe the dynamics of local birth-death processes coupled by nearest-neighbour interactions (diffusion).

In the following we pursue a different approach. Since the local population dynamics is assumed to be the same within all patches, it is sufficient to count the number of patches with a given number of individuals rather than keeping track of the number of individuals in each patch. Let denote the number of patches with inhabitants at a given time. The state of the population is described by the variables and . In the master equation, the time derivative of the probability of observing the system in a given state at time , is the rate of arriving to a given state from other states, minus the rate of transitions to other states. Thus, we find the master equation for the probability by considering all possible transitions between the states of the population, contributing to the change of .

For example, consider the effect of local births on the probability of finding the system in a given, ‘focal’ state at time . A birth event in a patch with individuals corresponds to the transition and . Thus, has the contribution

(6) |

from births in the focal state. This contribution is negative since any such birth in the focal state moves the system away from this state. In order to obtain the positive contributions to , we calculate the rate of arriving to the focal state by a birth in a different state:

(7) |

A compact way of writing these contributions is in terms of the raising and lowering operators , defined by their effects on a function (van Kampen, 1981):

(8) | ||||

In our example, the total contribution to from birth events can thus be written as

(9) |

Adding up the contributions due to death, emigration, and immigration we find:

(10) |

To simplify the discussion we assume that migration is instantaneous, this corresponds to taking the limit and . In this limit, the immigration rate to a patch is given by:

(11) |

and the corresponding master equation takes the form:

(12) |

The last two terms on the right-hand side of Eq. (12) describe instantaneous migration where . The terms involving Kronecker -symbols (see Table 1) on the right-hand side of Eq. (12) arise from enforcing that, for the migration of an individual, emigration must precede immigration. These terms are of higher order in .

The number of patches is given by . The following discussion is simplified by making this constraint explicit in the master equation (12). This is achieved by considering , the number of empty patches, to be be a function of :

(13) |

Using Eq. (13) we find the following master equation:

(14) |

Here the components of the vector denote the number of patches with individuals, and . Eq. (2.3) describes the stochastic population dynamics of the metapopulation model considered here.

In the next section we describe the approximate method of solving Eq. (2.3) adopted in the following. It corresponds to a systematic expansion of the master equation, using the number of patches in the population as a large parameter. This approach does not require that the number of individuals per patch is large, neither is it necessary to assume a time-scale separation between local population dynamics and migration. Our approach makes it possible to accurately describe the dynamics of metapopulations consisting of many small patches. In the limit of infinitely many patches, to lowest order in the expansion, a deterministic metapopulation dynamics is obtained. Stochastic fluctuations in metapopulations with a large but finite number of patches are described by the leading order of the expansion.

### 2.4 Expansion of the master equation

When is large we expect that the probability of observing a given distribution of changes only little when patches change in population size. It is important to emphasise that no assumption is made concerning the size of the changes to the population size of any single habitat. It is merely assumed that there are sufficiently many patches that they form a statistical ensemble. In this case it is perfectly possible to have big jumps in the population sizes of individual patches in the model without violating the assumption that the distribution of is smooth, and to allow the local population dynamics to depend sensitively on the local patch population size when there are few individuals in the patch — modelling for example the effect of abundance on mating success (Sæther et al., 2004; Melbourne and Hastings, 2008).

It is convenient to express the the components of the vector in terms of the scaled frequencies . The probability of observing is related to the probability by

(15) |

In contrast to the distribution of , the distribution of is expected to become approximately independent of for a large number of patches. A change of a single count leads only to a change of in . Since the probability is an approximately continuous function of in the limit of large values of , we seek an approximate solution of Eq. (2.3) by expanding this equation in powers of .

In the limit of we expect that stochastic fluctuations are negligible, so that the metapopulation dynamics becomes deterministic. It takes the form of a kinetic equation for :

(16) |

The right-hand side of this equation, , depends upon the details of the model analysed. For the model described in Section 2.1, the form of is given in Eq. (18). In the limit of , the question whether or not the metapopulation may persist is answered by finding the steady states of the system (16), given by . Persistence corresponds to the existence of a stable steady state with positive components for some values of . If, by contrast, the only stable steady state is , then the metapopulation will definitely become extinct. The time to extinction is determined by Eq. (16), and by the initial conditions.

For large (but finite) values of the metapopulation fluctuates around the stable steady states of Eq. (16). These fluctuations can be described by expanding the master equation (2.3) to leading order in . The stochastic fluctuations are expected to be small when is large, but they are crucial to the metapopulation dynamics as mentioned in the Introduction. When the number of patches is finite, the metapopulation must eventually become extinct (the state is the only absorbing state of the master equation). When is large, the time to extinction from a stable steady state of the deterministic dynamics is expected to be large. By analogy with standard large-deviation analysis, we expect that the time to extinction increases exponentially with increasing , giving rise to a long-lived quasi-steady state. Below we analyse its properties, and estimate the time to extinction of the metapopulation. Our results are consistent with the above expectation, and we show that the time to extinction depends sensitively on the parameters of the model (, and ).

#### 2.4.1 Deterministic dynamics in the limit of

To simplify the notation, we drop the tilde in Eq. (15) so that is the probability of observing, at time , a fraction of patches with one individual, a fraction of patches with two individuals, and so forth. In deriving the expansion of the master equation (2.3) we use the approach described by van Kampen (1981). Assuming that is a smooth function of , the action of the raising and lowering operators, Eq. (8), can be written as

(17) |

The master equation (2.3) is expanded as follows. We replace by in Eq. (2.3), insert Eq. (17), and expand in powers of . Keeping only the lowest order in , we arrive at an equation for that corresponds to deterministic dynamics of the form (16). We find that the components of are given by:

(18) | ||||

Here is the rate of immigration into a given patch, corresponding to Eq. (11). Since depends upon , the deterministic dynamics (16,18) is nonlinear. We note that Eqs. (16,18) correspond to the metapopulation model suggested by Casagrandi and Gatto (1999) and Nachman (2000). Here we have derived it by a systematic expansion of the exact master equation in powers of where is the number of patches. Our derivation emphasises the fact that Eqs. (16,18) approximate the metapopulation dynamics by a deterministic equation. This approximation improves as becomes larger, and becomes exact as . Arrigoni (2003) has shown that this limit holds under quite general assumptions for the underlying stochastic model, namely that the time-evolution of the probability measure over the states of the model converges to the deterministic time evolution as . Casagrandi and Gatto (1999), Nachman (2000), and others have studied how the stability of the steady states of the deterministic dynamics depends upon the parameters of the model. However, Eqs. (16,18) cannot be used to determine how stochastic population dynamics affects metapopulation persistence. In order to take the stochastic fluctuations into account, it is necessary to consider the next order in .

#### 2.4.2 Quasi-steady state distribution at finite but large values of

Consider a stable steady state of the metapopulation in the limit of infinitely many patches, that is, a stable steady state of the deterministic dynamics (16).

Metapopulations consisting of a finite number of patches exhibit random fluctuations caused by the random sequence of birth-, death-, and migration events. When the number of patches is large, these fluctuations are expected to be small, and the components of the vector are expected to fluctuate closely around those of the stable steady state . If the fluctuations around are small, the metapopulation may persist for a very long time (but extinction of this finite metapopulation is certain, as explained above.). This situation is commonly referred to as a ‘quasi-stable’ steady state. We analyse its properties by expanding the master equation to leading order in , using a standard method that is referred to as ‘WKB analysis’ (Wilkinson et al., 2007), as the ‘eikonal approximation’ (Dykman et al., 1994), or as the ‘large-deviation principle’ in the mathematical literature (Freidlin and Wentzell, 1984).

We briefly outline this method in the remainder of this subsection. For a comprehensive description, the reader is referred to Altland and Simons (2010), Elgart and Kamenev (2004), and Dykman et al. (1994). The method has been successfully used to describe fluctuations in finite biological populations, describing for example the spreading of epidemics (Dykman et al., 2008) or the risk of extinction of biological populations (see Ovaskainen and Meerson (2010) for a review).

In the quasi-steady state we expect and seek a solution of the master equation of the form

(19) |

The function is commonly referred to as the ‘action’. It not only depends upon , but also on the parameters of the problem (, , and in our case). This dependence is not made explicit in Eq. (19). In the limit of a large number of patches, this function determines the form of the probability distribution in the quasi-steady state, and its sensitive dependence upon the parameters , , and .

As mentioned above, is expected to be strongly peaked at in the limit of large values of . In other words, the quasi-steady state distribution concentrates on the deterministic fixed point as the number of patches tends to infinity. It is therefore convenient to define the action function such that . When the number of patches is large, we expect the distribution to be Gaussian in the vicinity of the steady state. Correspondingly, is expected to be approximated by a quadratic function of . By contrast, non-Gaussian tails of this distribution, corresponding to large deviations of from , reflect the particular properties of the extinction dynamics of the metapopulation.

The form of the function is determined by inserting the ansatz (19) into the master equation (2.3), making use of Eq. (17), and expanding in . One finds a first-order partial differential equation for :

(20) |

In our case, for the master equation (2.3), we obtain:

(21) | ||||

Here , for , and .

The solution of Eqs. (20,21) is found by recognising that Eq. (20) is a so-called ‘Hamilton-Jacobi equation’ (Freidlin and Wentzell, 1984). As a consequence, the action can be determined by solving the set of equations

(22) |

These equations have the form of Hamilton’s equations in classical mechanics with configuration-space variables . The variables

(23) |

are therefore referred to as ‘momenta’. An important property of Eq. (22) is that the function remains constant along the ‘trajectories’ , that is, along the solutions of Eq. (22).

How is the form of determined by these solutions? The quasi-steady state solution of Eq. (2.3) corresponds to solutions of Eq. (22) satisfying

(24) |

To every such solution corresponds an action obtained by integrating the momentum along the path :

(25) |

Here denotes the transpose of the vector .

The boundary conditions (24) are motivated as follows. Eq. (20) enforces . Further, observe that the stable steady state of the deterministic dynamics (16) corresponds to a steady state of Eq. (22). This can be seen by expanding the function to second order in around

(26) |

Here is given by Eq. (18), and the symmetric matrix has elements . The elements are given in appendix A. Eq. (26) shows that the dynamics for corresponds to the deterministic dynamics given in Eq. (16). The stability of the steady state is determined by the eigenvalues of the matrix

(27) |

The matrix has elements evaluated at . It is the stability matrix of the stable steady state of the deterministic dynamics (16). Its elements can be obtained from Eq. (A.2) in appendix A. Similarly, is evaluated at . Assuming that the steady state is stable, the eigenvalues of must have negative real parts. In our case it turns out that the eigenvalues are in fact negative. We write . The eigenvalues of occur in pairs , and . The steady state of Eq. (22) is thus a saddle. In other words, stochastic fluctuations allow metapopulations consisting of a finite number of patches to escape from the steady state (that is stable in the limit ) to extinction (). In general there are (infinitely) many such paths satisfying the boundary conditions (24). Freidlin and Wentzell (1984) formulated a variational principle for the most likely escape path: in the limit of large values of , the metapopulation goes extinct predominantly along this path. The quasi-steady state distribution reflects this property: configurations along this path are assumed with higher probability. According to the principle described by Freidlin and Wentzell (1984), the most likely escape path is the one with extremal action, Eq. (25).

The picture summarised above is schematically depicted in Fig. 2, for the case where the metapopulation persists in the limit of . Fig. 2a shows the - plane. As explained above, the point corresponds to the stable steady state of the deterministic dynamics, where the infinitely large metapopulation persists. The point corresponds to extinction, it is unstable. Solving Eq. (22) and inserting yields the deterministic dynamics (16), connecting these two fixed points. In the limit of , the metapopulation dynamics is constrained to the -axis and must approach , as the arrow on the -axis in Fig. 2a indicates. In other words, in the situation depicted in Fig. 2a extinction never occurs in the limit of .

In finite metapopulations, the situation is entirely different. Fig. 2a illustrates that there is a path from to , reaching at the so-called ‘fluctuational extinction point’ . Along this path, the momenta assume non-zero values. These variables characterise the sensitivity of the finite metapopulation to stochastic fluctuations ( in the deterministic limit as mentioned above). The form of (Fig. 2b) is determined by evaluating Eq. (25) along this path, satisfying conditions (24). In the limit of large values of , the average time to extinction of the metapopulation scales as (Dykman et al., 1994)

(28) |

The coefficient may depend on , as well as ,and . In the argument of the exponent, only the leading -dependence is explicit in Eq. (28). The argument of the exponential determines the sensitive dependence of the time to extinction. In one-dimensional problems with a single component (the case discussed in the review by Ovaskainen and Meerson (2010)), Eq. (25) shows that is given by the shaded area in Fig. 2a. In the case of our metapopulation model, by contrast, the vector has infinitely many components. It is therefore not possible, in general, to find the most likely path from to the fluctuational extinction point explicitly. In practice one may truncate the dynamics by only considering a finite number of variables (up to a maximal value of ). This is expected to be a good approximation when is taken to be much larger than the carrying capacity , since for . In our subsequent analysis of the problem we make use of the fact that the dynamics simplifies considerably in the vicinity of a bifurcation of the deterministic dynamics (Guckenheimer and Holmes, 1983; Dykman et al., 1994).

As mentioned above, near the steady state , the distribution (19) is expected to be Gaussian. This corresponds to an action quadratic in :

(29) |

The matrix , which is the covariance matrix of the distribution (19) multiplied by , can be obtained from the linearised dynamics, Eq. (22). Using we have

(30) |

According to Eq. (24), the dynamics must obey . It follows from Eq. (26) that the linearised dynamics satisfies this constraint provided

(31) |

This equation determines the covariance matrix in Eq. (29) in terms of the matrices and .

## 3 Results and discussion

In this section we summarise our results for the population dynamics of the metapopulation model described in Section 2.1, using the expansion of the master equation outlined in Section 2.4. This section is divided into two parts. We first discuss the limit of infinitely many patches. Second, we analyse the most likely path to extinction in finite metapopulations. We also summarise our results for the average time to extinction of the metapopulation, and demonstrate that it depends sensitively on upon the parameters of the model (, and ).

### 3.1 Infinitely many patches

It was shown in Section 2.4 that in the limit of infinitely many patches, the metapopulation dynamics is described by the deterministic equation (16), corresponding to the model proposed by Casagrandi and Gatto (1999) and Nachman (2000). In this subsection we briefly summarise our results on the persistence and the relaxation behaviour of the metapopulation model introduced in section 2.1.

#### 3.1.1 Metapopulation persistence in the limit of infinitely many patches ()

For sufficiently large emigration rates , the deterministic dynamics (16,18) has two viable steady states. The state is unstable, and there is a second steady state given by

(32) |

where is the Gamma function, is the rate of immigration into a patch in the steady state, and . Eq. (32) is most easily understood by recognising that the deterministic dynamics (16,18) takes the form of a master equation for the probabilities that a patch is occupied by individuals. Eq. (32) is equivalent to Eqs. (5a,b) in (Nachman, 2000) and to Eqs. (5,6) in (Casagrandi and Gatto, 2002). In Eq. (32), the factor is a normalisation factor (equal to the frequency of empty patches in the steady state) determined by the requirement that . Note that the factors in the product in Eq. (32) depend upon the rate . This gives rise to a self-consistency condition for the rate of immigration in the steady state:

(33) |

The steady state (32) is stable provided the eigenvalues of (this matrix is introduced in the previous section and its elements are given in appendix A) have negative real parts. This is the case when the steady-state immigration rate is larger than zero. The steady-state immigration rate is expected to decrease when the emigration rate decreases. Consider for instance decreasing , defined in Eq. (4), while keeping all other parameters constant. As approaches a critical value, , we observe that . At the same time all components of tend to zero, and . Below this critical point, that is for , the stable steady state ceases to exist (and turns stable). In the limit of infinitely many patches (), the persistence condition

(34) |

ensures that a stable steady state exists, with for some non-zero values of . This persistence criterion for infinitely large metapopulations is equivalent to the criterion suggested by Chesson (1984) and re-derived by Casagrandi and Gatto (2002), namely that the metapopulation persists provided that the expected number of emigrants from a patch with initially one individual and into which immigration is excluded is greater than unity. In a slightly different form this principle is quoted by Hanski (1998), namely that the expected number of successful colonisations out of a given patch during its lifetime in an otherwise empty metapopulation should be larger than unity. We emphasise that this criterion (or any of the equivalent criteria, such as (34)) cannot be used to determine how stochastic fluctuations in finite metapopulations affect their persistence.

The critical value is obtained by analysing the steady-state condition (33) for close to . We write

(35) |

and expand the steady-state immigration rate in powers of (note that vanishes at ):

(36) |

The constant is given in appendix B. Expanding Eq. (33), we find to lowest order in a condition for :

(37) |

Fig. 3 shows how the solution of Eq. (37) depends upon the carrying capacity for the model introduced in section 2.1. The critical migration rate is shown as a solid line in the left panel of Fig. 3, in the following referred to as the ‘critical line’. Above this line, a metapopulation consisting of an infinite number of patches persists. Expanding the condition (37) for large carrying capacities we find that (see appendix C):

(38) |

Fig. 3 also shows as a function of for six stable steady states corresponding to different values of and . Results similar to those shown in the six panels on the right-hand side of Fig. 3 are given in Fig. 2 in (Nachman, 2000), and in Fig. 2 a-c in (Casagrandi and Gatto, 2002).

One observes how the number of empty patches tends towards unity as the critical line is approached. To leading order in we have:

(39) |

The constant is given in appendix B. Similarly, the expected number of individuals per patch approaches zero:

(40) |

as . The constant is given in appendix B. Fig. 4 shows the average number of individuals per patch, and as a function of . Shown are numerical solutions of the steady-state condition (33), solid lines, of Eqs. (39) and (40), dashed lines, and of direct numerical simulations, symbols, of the metapopulation model described in section 2.2. The numerical experiments were performed for and patches and averaged over an ensemble of different realisations. In principle, in a finite metapopulation no stable steady states exist with for some . But the larger the number of patches, the larger the average time to extinction is expected to be. , it turns out, is sufficiently large for the parameter values chosen that extinction did not occur during the simulations. Consequently we observe good agreement between the direct numerical simulations and the numerical solution of the deterministic steady-state condition. When the number of patches is small, by contrast, this is no longer the case. In order to determine the quasi-steady state in such cases, the simulations must be run for a time large enough for the initial transient to die out, but shorter than the expected time to extinction. Stochastic realisations leading to extinction during the simulation time must be discarded.

#### 3.1.2 Relaxation dynamics in the limit of and Levins’ model

We continue to discuss metapopulations with infinitely many patches. In this section we analyse how the metapopulation relaxes to the stable steady state . This question is important for two reasons. First, when the relaxation time is much smaller than the expected time to extinction (that is, when the number of patches is large enough), Eq. (16) describes the relaxation dynamics well. Second, and more importantly, the deterministic dynamics exhibits a ‘slow mode’ in the vicinity of the critical line (in Fig. 3): for small values of , it turns out, there is a particular linear combination of the variables that relaxes slowly, with a rate proportional to . This slow mode dominates the relaxation dynamics which is thus essentially one-dimensional for small values of . This fact is well known in the theory of dynamical systems (Guckenheimer and Holmes, 1983), see also Dykman et al. (1994). For our model, this fact has important implications. Essentially, the deterministic dynamics close to the critical line in Fig. 3 reduces to Levins’ model, as we demonstrate in this section.

Fig. 5 shows the relaxation of the average number of individuals and of the frequency of empty patches to their values in the stable steady state , for three different initial conditions. Shown are numerical solutions of Eqs. (16,18), solid lines, as well as results of direct numerical simulations, symbols. The numerical experiments were performed for , and patches, and averaged over realisations. Even for patches, and for the parameter values chosen in Fig. 5, the average time to extinction is large enough so that extinction did not occur in the numerical simulations. The agreement between the direct numerical simulations for and the solution of Eqs. (16,18) is good. For and we can observe deviations to the solution of Eqs. (16,18). Here the effect of the finite number of patches becomes apparent.

In general Eqs. (16,18) must be solved numerically. We now show that the deterministic dynamics simplifies considerably close to the critical line in Fig. 3, when is small. At we have that , and consequently . For small positive values of , the components of are of order , and we expand Eq. (16) in powers of and (Dykman et al., 1994):

(41) |

Here are the elements of the matrix evaluated at the critical line (). They are obtained from (A.2):

(42) | ||||

The slow mode exists because this matrix has one zero eigenvalue, . All other eigenvalues are negative. At finite but small values of the corresponding matrix (A.2) evaluated at the steady state has one small eigenvalue, . To identify the slow mode we diagonalise given by (42). Since this matrix is not symmetric, its left and right eigenvectors differ:

(43) |

Here denotes a left eigenvector with components . denotes a right eigenvector with components . We take the eigenvectors to be bi-orthonormal, . As before, the eigenvalues are ordered as . Multiplying Eq. (41) from the left with and making use of Eq. (16) yields the equations of motion of (that is, the equations governing the time-evolution of ):

(44) | ||||

At (that is, for ) we have that . For small values of we find for . While , , approach constants as , tends to zero in this limit. This implies

(45) |

In the following, is therefore referred to as ‘slow mode’, whereas for are termed fast variables. As explained by Dykman et al. (1994), the fast variables rapidly approach quasi-steady states

(46) |

that depend on the instantaneous value of the slow variable, . Terms including fast variables ( for ) are not kept on the right-hand side of Eq. (46) because they are of higher order in . This is a consequence of the fact that close to the steady state, is of order . Eq. (46) is correct to order . The coefficients and are given by

(47) | ||||

for . The elements and are given in appendix A. The equation of motion for the slow mode is to third order in :

(48) | ||||

In the remainder of this subsection we neglect the -terms and higher-order terms in Eq. (48). The -terms are discussed in subsection 3.2. To second order in , the equation of motion for is:

(49) |

According to Eq. (47), the coefficients and are given in terms of the components of the left eigenvector and the components of the right eigenvector of . For the right eigenvector we find, from Eqs. (42) and (43),

(50) |

the first component being . For the left eigenvector, we obtain the following recursion (with the boundary condition for ):

(51) |

This recursion is solved by:

(52) | |||

We show in appendix C that the elements of approach the following limiting form as :

(53) |

We choose such that approaches unity for as . This corresponds to the choice , resulting in

(54) |

in the limit of . In this limit, and for large values of , we see that approaches the vector . The convention leading to Eq. (54) also fixes which must be chosen so that . Explicit expressions for and , obtained from Eqs. (47), (50), and (52) are given in appendix B.

Eq. (49) has two steady states and . Note that is positive since . These two steady states correspond to the steady states and of Eq. (16). Comparison with Eq. (32) shows that is proportional to to lowest order in (that is, to lowest order in ). In fact we have

(55) |

Fig. 6 illustrates how the deterministic dynamics relaxes to the stable steady state . Shown are the expected number of individuals per patch, and the fraction of empty patches as functions of time, determined from numerical solutions of Eqs. (16,18). Curves for three different initial conditions are shown (solid lines). The corresponding solutions of Eq. (49) are shown as dashed lines. Initially the variable is not much slower than the -variables for , and the approximate one-dimensional dynamics (49) is not a good approximation. But the solution of Eqs. (16,18) rapidly relaxes to a form where becomes a slow mode, and Eq. (49) accurately describes the slow approach to the steady state.

Comparing Eqs. (49) and (1) shows that the slow mode obeys Levins’ equation. The results of this subsection allow us to compute and in terms of , , and . We find to order

(56) | ||||

The contributions to of order can be computed explicitly using the formulae derived above. For the sake of brevity we do not specify these contributions here. In Sec. 3.2.3 we give an explicit formula valid in the limit of large .

We see that that Levins’ model describes the dynamics of the metapopulation regardless of whether the patches are strongly coupled or not, provided the metapopulation is sufficiently close to criticality (that is, close to the critical line in Fig. 3). But we emphasise that in general the slow variable is not the fraction of occupied patches as envisaged by Levins. The variable is given by . Fig. 7 shows how the components