## Abstract

Critical population phase transitions, in which a persistent population becomes extinction-prone owing to environmental changes, are fundamentally important in ecology, and their determination is a key factor in successful ecosystem management. To persist, a species requires a suitable environment in a sufficiently large spatial region. However, even if this condition is met, the species does not necessarily persist, owing to stochastic fluctuations. Here, we develop a model that allows simultaneous investigation of extinction due to either stochastic or deterministic reasons. We find that even classic birth–death processes in spatially extended ecosystems exhibit phase transitions between extinction-prone and persistent populations. Sometimes these are first-order transitions, which means that environmental changes may result in irreversible population collapse. Moreover, we find that higher migration rates not only lead to higher robustness to stochastic fluctuations, but also result in lower sustainability in heterogeneous environments by preventing efficient selection for suitable habitats. This demonstrates that intermediate migration rates are optimal for survival. At low migration rates, the dynamics are reduced to metapopulation dynamics, whereas at high migration rates, the dynamics are reduced to a multi-type branching process. We focus on species persistence, but our results suggest a unique method for finding phase transitions in spatially extended stochastic systems in general.

## 1. Introduction

The stability and survival of populations depend on climate, resources, predation, dispersal and other factors. Environmental changes may not only lead to gradual changes in population abundance but may also lead to sudden shifts and even to sudden extinctions [1]. Identification of such phase transitions between persistent and extinction-prone populations, promoted by environmental changes, is fundamentally important in ecology for several reasons. First, extinction is a fundamental phenomenon that shapes ecosystems [2]. Second, estimation of extinction risks is a key factor in evaluation of the effects of environmental changes [3]. Third, this estimation may suggest conservation management activities, such as building corridors between fragmented habitats to enhance migration [4].

Populations in nature commonly are subdivided. Each individual is spatially restricted and interacts locally with its neighbours for long periods of time, but eventually the individual migrates and then interacts with its new neighbours [5]. Migration may occur continuously over time, such as in infection diseases where bacteria or viruses migrate from host to host, or it may occur during a dispersal stage, such as with pollen of plants, with the larvae of benthic marine life, with the adolescents of many vertebrates and during the post-teneral migratory phase of adult insects [6,7]. The spatial structure of a population plays a major role in determining its persistence. For example, space may promote longer survival periods as migrants may re-establish locally extinct populations [5] (rescue effect). However, a major question is how to distinguish between persistent and extinction-prone populations. Populations are inherently discrete and finite, and therefore, sooner or later, all populations are prone to extinction due to their stochastic nature. The *mean time to extinction* (MTE) changes smoothly with environmental parameters [8], and therefore it is hard to recognize a specific extinction threshold, namely to put the finger on a certain environmental change that triggers extinction.

A common approach to address this question and to identify and analyse phase transitions is metapopulation modelling, in which infinitely many spatial sites are considered, each of which has two possible states, either populated or empty, and the dynamical variables are the fractions of occupied sites [9–11]. Because these models describe spatially infinite systems, they may exhibit phase transitions. Specifically, when the mortality rate, at which occupied sites are emptied, is sufficiently low, the steady-state fraction of occupied sites is greater than zero, and the population is persistent. Nevertheless, this fraction gradually decreases as mortality increases, until it approaches zero where mortality rate is higher than the rate of colonization of empty sites by migrants [9]. This demonstrates a second-order phase transition.

The metapopulation approach has lead to major advances in the study of population dynamics, but at the same time, it still suffers some disadvantages and may lead to misleading conclusions in several cases, because the assumption that each site has only two possible states is sometimes over-simplifying and within-site dynamics may be important. First, local fixation dynamics in which a site's population with initially a single individual either becomes established or becomes extinct, is not necessarily much faster than extinction of already established site's populations, and therefore, migration during local fixation is sometime important. Second, the effect of migration rate on extinction is not explicit in the metapopulation dynamics, and therefore negative consequences of migration such as dilution of the source population, and also positive impacts such as ongoing immigration are commonly ignored. Third, metapopulation models do not take into account more realistic dynamics such as a local Allee effect, in which locally diluted populations decrease in size and a population can grow only when the local density is higher than a certain density.

Here, we overcome these difficulties using a novel approach. We extend the metapopulation model by allowing any number of individuals in each site. We assume a spatially infinite system of locally interacting individuals, while still allowing detailed stochastic local dynamics and spatially heterogeneous environments. We derive our model equations using a master equation approach, which takes into account the dynamics within each site explicitly. We first introduce our method for a fully connected, homogeneous environment, and later we extend the theory to allow mixed connectivity and heterogeneous environments. This approach allows simultaneous and comprehensive investigation of both positive and negative effects of migration, as well as of other dynamical rates, on the dynamics. We show that higher migration rates promote spatial synchronization that prevents sudden extinctions, but at the same time, very large migration rates prevent selection of habitats of higher quality, thereby promoting extinction. Therefore, intermediate migration rates are optimal for persistence. We also show that the spatial dynamics of populations with several metastable local configurations (e.g. with an Allee effect) may exhibit first-order phase transitions. Thus, environmental changes may promote a mass reduction in the population—a process that will not be reversed even after the environment gets back to its vital state.

## 2. Intermediate migration rates promote persistence

We consider a population that is subdivided into habitable sites. Each site is occupied by a finite number of individuals, and therefore, the dynamics within each site are stochastic. Here, in our first example, we assume identical sites and consider the following classical birth–death dynamics that comprise the following four processes (figure 1). For every site *i*, (i) a resident individual may give birth to another individual, (ii) a resident individual may die, and (iii) any interaction between two individuals may lead to a death due to competition. In addition, (iv) for every two sites *i* and *j*, an individual may migrate from *i* to *j*. These processes are summarized in the following rate equations: for all *i* and *j*,
2.1(*i*, *j* = 1, 2, … , *M*), where *A _{i}* is an individual located at site

*i*,

*β*is the birth rate,

*δ*the natural death rate,

*γ*is the strength of competition,

*λ*is the migration rate,

*c*is the probability that a migrant lands on a viable site and

*M*is the number of sites. Note that, in this example, all sites are identical and have the same rates.

First, consider no migration (*λ* = 0), so the system encompasses *M* independent sites. In a given site *i*, the probability for birth is given by *βn*_{i}, where *n _{i}* is the number of individuals in the site and the probability for mortality owing to either death or competition is given by

*δn*

_{i}+

*γn*

_{i}(

*n*− 1). The mortality probability in site

_{i}*i*becomes equal to the birth probability when

*n*approaches the metastable population density, For and sufficiently short timescales, the

_{i}*mean-field*approximation is valid, and

*n*evolves approximately according to [12] which implies that

_{i}*n*approaches

_{i}*N*

_{0}. However, owing to its discrete nature, the population ultimately becomes extinct (

*n*= 0) with probability one, in which case it will never recover (absorbing state; the MTE of a given site in (2.1) increases exponentially with

_{i}*N*

_{0}if [8]). Therefore, the mean-field equation is invalid where long timescales are considered. Instead, the single-site dynamics without migration is fully described by the following master equation, which is exact and describes the dynamics of

*P*(

_{n}*t*)—the probability that a given site is occupied by exactly

*n*individuals at time

*t*[12]: The population ultimately becomes extinct as this equation comprises a single steady state that is a global attractor: and for

*n*≥ 1.

Next, assume that the sites are coupled by migrants (*λ* > 0). The master equation becomes increasingly complex as the number of coupled sites, *M*, increases because an appropriate master equation, in this case, describes the dynamics of —the probability that site 1 is occupied by exactly *n*_{1} individuals, site 2 by exactly *n*_{2} individuals, etc. Still, for any finite *M*, the MTE changes smoothly with any of the rates (no phase transitions), and the population still ultimately becomes extinct as the only absorbing state corresponds to extinction in all sites (see [13]).

Here, we extend the master equation approach to include infinitely many sites, which may allow phase transitions. As in the metapopulation approach, we consider only relative abundance and not the entire configuration. Specifically, we examine the dynamics of *S _{n}*(

*t*)—the fraction of sites occupied by exactly

*n*individuals at time

*t*. The time evolution of

*S*is derived from the process (2.1) with infinitely many identical sites: 2.2where is the average population per site at time

_{n}*t*.

We emphasize that this equation is exact and rigorously derived from the process (2.1) where *M* → *∞*. The terms in the first three lines are equivalent to the respective terms in the single-site master equation, and result from within-site birth, death and competition: the probability that a process, say birth, is occurring within a given site between *t* and *t* + d*t* is given by *βn*d*t*, and this is also the fraction of sites occupied by *n* individuals that become occupied by *n* + 1 individuals during that time (each site increases by one with probability *βn*d*t*, and there are infinitely many independent events, because there are infinitely many sites occupied by exactly *n* individuals). The term in the fourth line is somewhat trickier and results from migration. The number of individuals in a site may either decrease if a migrant leaves, or increase if a migrant arrives. The key observation is that these two events are independent in the infinite-sites limit. Therefore, the rate at which migrants leave is given by *λn*, and the rate at which migrants arrive is given by *λcν*, because migrants may arrive from all other sites with equal probabilities.

Note that, despite the similarity in the form of the equations for *P _{n}* and for

*S*, they have different meanings:

_{n}*P*is the

_{n}*probability*that the site is occupied by

*n*individuals (in a single-site system), whereas

*S*is the

_{n}*actual fraction*of sites occupied by

*n*individuals (in a multiple-sites system). In addition, note that if

*M*is large but finite, then for short timescales, the dynamics of

*S*follow equation (2.2) plus fluctuations that are proportional to (

_{n}*MS*)

_{n}

^{−}^{1/2}. By contrast, for long timescales, the population ultimately becomes extinct for any parameter set, but the MTE increases exponentially with

*M*if the population is persistent (see [13]).

The global extinction state (*S*_{0} = 1 and *S _{n}* = 0 for

*n*≥ 1) is a steady-state solution of equation (2.2) for any parameter set, but there may be other steady-state solutions as well. To examine the dynamics, we solved equation (2.2) numerically as a system of

*n*

_{max}coupled differential equations, where

*S*is simulated for all 0 ≤

_{n}*n*≤

*n*

_{max}. Here,

*n*

_{max}is a truncation threshold, which is sufficiently large such that

*S*is effectively zero where

_{n}*n*>

*n*

_{max}, and consequently our results are indistinguishable from those of the actual infinite system in equation (2.2). Our results show that for sufficiently small migration rates,

*λ*<

*λ*

_{c}

_{1}, the population ultimately becomes extinct as the only steady-state solution of equation (2.2) is the global extinction state. However, as

*λ*surpasses the threshold

*λ*

_{c}

_{1}, a sustainable population steady state emerges (figure 2

*a*). The dynamics of any initial population approach this sustainable steady state and the global extinction state become dynamically unreachable. For

*λ*just above the threshold, most sites at the sustainable steady state are empty, and a few sites are occupied by nearly

*N*

_{0}individuals (figure 2

*a*, insets). The fraction of occupied sites gradually increases with

*λ*, which demonstrates a second-order phase transition at

*λ*=

*λ*

_{c}

_{1}in the limit

*M → ∞*, owing to the synchronization of sites [9,14,15]: the migration rate is large enough to ‘feed’ empty sites fast enough to ensure that they become repopulated (rescue effect), and in turn, each populated site has nearly

*N*

_{0}individuals (figure 2

*a*). Nevertheless, if the environment is heterogeneous (

*c*< 1), then as

*λ*further increases above some point, the average population per site at steady state decreases, because many migrants do not land on a viable site. Finally, as

*λ*surpasses a second threshold,

*λ*

_{c}

_{2}, the population becomes extinction-prone again, with global extinction again the only steady state (figure 2

*a*).

To calculate the first threshold, *λ*_{c}_{1}, note that when , local fixation duration is much shorter than the MTE of an established site's population. Moreover, when , migration has only a weak influence during local fixation. In addition, the metastable population density, *N*_{0}, is barely changed by migrants. Therefore, the dynamics are reduced to metapopulation dynamics, and we may assume that each site is either empty or populated by *N*_{0} individuals. The extinction rate equals a single site's (MTE)^{−1}, which was derived by Assaf & Meerson [16] using the Wentzel–Kramers–Brillouin (WKB) approach, and colonization rate equals *λN*_{0}*c* times the fixation probability of a single individual. Because the population can stably persist if and only if the average number of established migrants that are being sent from a single site before its population becomes extinct is greater than one, it follows that
(see electronic supplementary material). This implies that *λ*_{c}_{1} decreases exponentially with *N*_{0}.

To calculate the second threshold, *λ*_{c}_{2}, note that if , then sudden extinctions do not play a significant role, because colonization of empty sites is rapid. Then, the average population per site at steady state, *ν**, decreases with *λ*, because migrants may land on non-viable sites when *c* < 1. The average population size per site becomes *ν** ≈ (*β* − *δ* − *λ*(1 − *c*))/*γ* (figure 2*a*). The second threshold appears where *ν** = 0, namely where linear growth equals zero:

Above this threshold, the average growth rate of small populations is negative even without sudden extinctions. Note that if migrants always land on viable sites (*c* = 1), then *λ*_{c}_{2} → *∞*, which implies that the steady-state population size increases with migration rate and asymptotically approaches its maximum as *λ* → *∞*. This implies that the mechanism for the reduced sustainability at higher migration rates is the inability to select for and to stay in the suitable habitats when the environment is heterogeneous.

## 3. Population collapse: first-order extinction transition

In §2, we demonstrated a second-order phase transition, in which an extinction-prone population gradually increases its persistence as *λ* surpasses *λ*_{c}_{1}. However, first-order phase transitions, in which sudden population collapses following environmental changes, may characterize birth–death processes as well. Consider the following processes: for all *i* and *j*,
3.1(*i*, *j* = 1, 2, … , *M*). Without migration, the corresponding mean-field approximation is

For some choices of the parameters, the mean-field potential has two minima, denoted by *N*_{1} and *N*_{2}, and the number of individuals at each site stays close to one of these metastable population abundance, *N*_{1} or *N*_{2}, for some time, before it switches to the other owing to some stochastic event. For *α* = 0, the lower abundance corresponds to extinction (*N*_{1} = 0, which becomes an absorbing state of the system), and the dynamics become a version of a birth–death process with the Allee effect, in which populations require to exceed a certain density to have a positive mean growth rate within a given site.

Next, we allow migration and examine the dynamics with infinitely many sites: 3.2

For sufficiently small *λ*, the population is asynchronous as some sites are occupied by nearly *N*_{1} individuals and the rest by nearly *N*_{2} individuals (figure 2*b*). However, as *λ* surpasses *λ*_{c}_{1}, the population exhibits a first-order phase transition and becomes synchronous at either *N*_{1} or *N*_{2}, depending on the initial configuration: the extinction (low abundance) state exists in parallel with the persistent (high-abundance) state for *λ* > *λ*_{c}_{1}, and the fate of the population depends on the initial conditions (figure 2*b*). The second transition at *λ*_{c}_{2} is still exhibited, and is also a first-order one. This demonstrates that a change in *λ*, which may be induced by some environmental change, may result in an abrupt, irreversible population collapse. More generally, the dynamics of *S _{n}* may entail more than two bifurcations if the mean-field potential has more minima. We emphasize that the more complicated site dynamics in this model preclude the reduction to metapopulation model. Moreover, the second threshold is nonlinear and cannot be calculated by the procedures used in the previous section.

## 4. Heterogeneous environments

In the previous sections, we considered fully connected, identical sites. However, sites in nature may vary significantly in their quality, and moreover, the transition between sites is not uniform due to demography and to preference of sites with higher quality. Here, we generalize our model and allow a variety of site types, where a site of type *k* appears with frequency *q _{k}*, and has its unique characteristic rates (hereafter denoted by a subscript

*k*). The transition probabilities between different types are given by the connectivity matrix

*D*, namely the probability that a migrant from a

*k*-site arrives to an

*l*-site is

*λ*d

*, where*

_{kl}*Σ*

*d*

_{l}*= 1 for all*

_{kl}*k*and

*l*. If migration is uniform, then

*d*=

_{lk}*q*, but we allow the more general case. However, we still assume global connectivity in the sense that for all

_{k}*k*and

*l*, each site of type

*k*is connected to many sites of types

*l*, such that in the limit of many sites (

*M*→

*∞*), the probability that a migrant arrives to a given type depends only on the average site occupancies of each type and on

*D*(see the electronic supplementary material). Note that similar assumption is required in metapopulation models that allow several site types [10,11]. In analogy to previous sections, we formulated an equation for the dynamics of

*S*

_{k,n}, the fraction of sites of type

*k*that occupy exactly

*n*individuals (equation (S1) in electronic supplementary material). This equation has the same structure as equations (2.2) and (3.2), except that a subscript

*k*is added and

*c*→

*ν**Σ*

*(*

_{l}*q*/

_{l}*q*)

_{k}*d*, where the last term is the average flow into a site of type

_{lk}*ν*_{l}*k*as

*ν*

_{l}is the average per-site population at sites of type

*l*. Note that

*c*is omitted, because the original model is reconstructed when the fraction of habitable sites is

*c*and the fraction of lethal sites is 1 −

*c*, namely

*q*

_{1}=

*c*,

*q*

_{2}= 1 −

*c*and

*δ*

_{2}→

*∞*.

Equation (S1) in the electronic supplementary material can be solved numerically, but several analytic results can still be derived despite its complexity. For sufficiently low migration rates, sites from all types barely send migrants before fixation and non-viable sites (sinks) barely send migrants before their population becomes extinct. Therefore, for a single-well mean-field potential (e.g. equation (2.1)), the dynamics can be reduced to metapopulation dynamics if the metastable state of viable sites with positive linear growth (sources) is sufficiently large. In this case, *λ*_{c}_{1} can be derived from the combination of *D*, of the average extinction rates of all sources, and of the average colonization rates of all sources (see §2), using the formalism given in [11]. These rates can be calculated using the same technique given in §2 and in the electronic supplementary material, where for each type, the metastable population size of the type substitutes *N*_{0}.

To calculate *λ*_{c}_{2} and examine the deterministic stability of the extinction state, we consider the opposite limit where sudden extinctions do not play a significant role, namely sites are synchronized (). Then, linear growth determines the survival of small populations, and the transition can be calculated using the theory of multi-type branching processes [17,18], so the growth of small populations is given by the largest eigenvalue of the linear growth matrix A. For example, for the dynamics given by equation (2.1), *A _{ii}* =

*β*

_{i}−

*δ*

_{i}−

*λ*(1 −

*d*) and

_{ii}*A*=

_{i≠j}*λd*

_{ij}. Therefore,

*λ*

_{c}

_{2}zeroes the largest eigenvalue of

*A*. Note, however, that

*λ*

_{c}

_{2}does not necessarily exist: in the limit

*λ*→

*∞*, the distribution among site types is fully determined by the connectivity matrix

*D*, and therefore, the second extinction transition exists if and only if where is the

*k*th component of the positive eigenvector of

*D*−

*I*(provided that

*D*−

*I*is a regular matrix; note that for uniform migration). Namely,

*λ*

_{c}

_{2}exists if and only if the average growth rate of a small population is negative where different types are weighted by . Otherwise,

*ν*decreases and approaches some positive value asymptotically as

*λ*→

*∞*.

## 5. Extreme migration and global synchronization

Next, we derive the asymptotic behaviour at *λ* → *∞* (*λ* much larger than the other rates) for the general case. In this ‘hydrodynamic limit’ [19], the dynamics are reduced to those of a single, infinite site that occupies the entire population and exhibits well-mixed dynamics. Each individual has the same probability to interact with any given other individual, and therefore, when environment is homogeneous, the population density follows the mean-field equation, whereas if the environment is heterogeneous, then the density follows some effective equation, which has the same form as the mean-field equation, but its rates are given by a weighted average of the actual rates. In the following, we derive this equation from our model.

Because migration is much faster than local processes, separation of timescales is possible. First, we find the steady-state solution of the dynamics with migration being the only process, namely for all *k*
where . This implies a Poisson distribution of site occupancies for each type, while keeping the overall average number of individuals per site, *ν*, constant. Namely, at all *t* and for all *n* and *k*,
and the only unknown variable is *ν*(*t*), which continuously evolves owing to the intra-site (slow) birth–death processes, according to the effective potential *U*_{eff}:

Next, to find *U*_{eff}, we substitute *S _{k}*

_{,n}into the corresponding equation: for homogeneous environments, we substitute

*S*(Poisson distribution) into equation (2.2) or (3.2), which implies that

_{n}*U*

_{eff}equals the mean-field potential; for heterogeneous environments, we substitute

*S*

_{k}_{,n}into equation (S1) in the electronic supplementary material, which implies where

*U*

_{MF,k}is the mean-field potential for a site of type

*k*. Therefore, each local minimum of

*U*

_{eff}corresponds to a steady state of the system in the large migration limit.

## 6. Discussion

We have demonstrated that spatially extended populations are persistent within a certain range of intermediate migration rates. Very low migration rates prevent synchronization, thereby promoting extinction at each site with only limited ability to recover. Very large migration rates prevent the population from selecting the better habitats: individuals do not stay in sites with higher quality, and instead they travel along a wider range that also includes sites with lower quality. We emphasize that our method applies only for systems with high connectivity. Nevertheless, our direct Monte Carlo simulations demonstrate that this result also applies to lower dimensional systems with local connectivity (nearest neighbours; figure 3).

This result is consistent also with previous studies that demonstrated maximal MTE for intermediate migration rates [14,20–22]. The transition from extinction-prone populations at low migration rates to persistence populations at higher migration rates is evident in metapopulation models [9], owing to the rescue effect where migrants from yet populated sites ‘rescue’ and repopulate sites where the population locally becomes extinct. This transition was found also in several experiments [23,24], and also in models of disease propagation, which found that disease may outbreak only above a certain migration or infection rate [25]. This transition was studied also in systems with local connectivity, and found, using renormalization methods, to be analogous to the well-known percolation transition [15]. On the other side of the coin, reduced persistence at higher migration rates was demonstrated in several studies, where two distinct underlying mechanisms were suggested. First, some studies considered homogeneous environment and suggested that over-synchronization may be fatal if the dynamics are either oscillatory or chaotic [20,21,26,27], or cyclic with two or more species [28], because, then the ability of a homogeneous environment to host different modes of the dynamics may promote persistence, but is masked by migration. Second, some other studies suggested that not all migrants arrive to their destination, and consequently too high mutation rate is fatal [14,29,30]. This mechanism is mathematically very similar to our model when only two site types exist, habitable and lethal (§§2 and 3).

When a sufficiently large population is required to initiate a sustainable population at a given site (‘local’ Allee effect), we demonstrated a first-order transition, in which the fate of the population depends on the initial configuration. Therefore, the entire population may irreversibly collapse either if migration rate decreases below some threshold or if it increases above some other threshold (figure 2*b*). Similar behaviour was revealed in metapopulation models that allow several site types (heterogeneous environments) [10,31], and it is also consistent with the experiment of Molofsky & Perdy [24] that examined extinction probabilities of the annual plant *Cardamine pensylvanica* and found a sudden jump in the steady-state population abundance below a certain migration rate. This type of transition is particularly important in ecology. First, it is much harder to foresee it, because the reduction in the steady-state population is not gradual. Second, after extinction has occurred, it is much harder to restore the system back, because for that, active introduction of sufficiently many individuals is necessary. One cannot overcome the local Allee effect solely by locally introducing many individuals, because, as we demonstrated, this will not solve the problem on the global scale, in which a ‘global’ Allee effect has emerged. Therefore, sufficiently many sites must be restored to allow recovery.

We demonstrated a novel method for transforming master equations characterizing the dynamics of well-mixed stochastic processes into deterministic differential equations characterizing the dynamics of infinitely many such systems that are coupled by migration. These transformed equations entail phase transitions as parameters vary, in contrast to the original master equation. Using this method, we showed that our model may be reduced into two existing theories in two different regimes of the migration rate. For low migration rates, the simpler version of our model (no Allee effect, equation (2.1)) may be reduced to metapopulation models that assume only two site-states (empty and populated), and exhibits ‘stochastic extinction transition’ from a population that is prone to extinction by demographic stochasticity to a persistent population [10,11,32,33]. For high migration rates, the population is highly synchronized, sudden extinctions of a single site's population are less important, and systems with small populations are reduced to deterministic models of population dynamics on graphs, namely multi-type branching processes [17,18]. These may undergo ‘deterministic extinction transition’ when deterministic linear growth becomes negative.

We emphasize that our models' dynamics, equations (2.2), (3.2) and electronic supplementary material (S1), simulate exactly the corresponding dynamics of highly connected systems given by equations (2.1), (3.1) and (3.1) with mixed connectivity, respectively, in the limit of infinitely many sites. No further approximations are used. This may serve as a baseline for investigation of the dynamics, at least numerically, before further approximations are made. Casagrandi & Gatto [14] used a model somewhat similar to our homogeneous model in §2, but the authors used several approximations. Like us, they concluded that intermediate migration rates induce the optimal balance between synchronization and habitat selection, but at the same time they argued that extinctions mediated by demographic stochasticity are probable even when migration rate is large, whereas we showed that demographic stochasticity cannot lead to extinction for sufficiently large migration rates when the system is large (a difference that may result from one of the approximations used in [14]), which implies distinction between deterministic and stochastic extinctions. Other studies considered continuous intra-site population densities, together with some ad hoc assumptions about the effects of migrants [31,34]. In another study, Van den Broeck *et al.* [35] developed a similar model, but instead of the master equation, they used the Fokker–Planck approximation that encompasses continuous density in each site. The authors considered a fourth-order polynomial and demonstrated spontaneous breaking of symmetry as noise varied. In implementing this approach to population dynamics, however, one must be careful, because the Fokker–Planck approximation underestimates extinction risks (unlike the master equation that is exact) [8]. Moreover, the consistent properties at the large migration limit (§5) are not obtained by the Fokker–Planck approach, which may be misleading when analysing large migration rates.

Our method applies directly to highly connected populations, where each individual may migrate to any other site. This assumption is common in the ecology literature and is realized in a large variety of ecologically realistic scenarios. First, high connectivity may be promoted by dispersal mechanisms that are non-local [5–7]: pollens that are non-locally transported by insects, adult insects that travel distances to lay their eggs, etc. Second, it may emerge in infectious diseases for which propagation follows non-local transports as hosts (e.g. humans) travel long distances, or when the hosts under focus are well-interacting. Third, generally, high connectivity may be realized in systems where birth and death rates are much lower than the migration rate, such that an individual may explore the entire system before a birth or a death occurs (see §5). However, our method does not apply to systems with local connectivity, where migration is restricted to neighbouring sites. Methods for identification of phase transitions in systems with local connectivity are more complicated and include renormalization group methods [15], or repeated numerical calculations of the MTE for various system sizes, where for extinction-prone populations, the MTE typically increases linearly with the system's size, whereas for sustainable population, it increases exponentially with the system's size [36]. Note that the validity of the global connectivity assumption was exhaustively discussed in the physics literature, particularly in the context of the Becker–Döring equations, which are somewhat methodologically similar to our equations and describe the distribution of cluster sizes, assuming no spatial correlation between distinct clusters [37].

Our approach implies a unique measure for population persistence and performance in stochastic systems: the steady-state average population (local) density, in the limit where the system is spatially infinitely large. This measure has some advantages over measuring the average population density assuming deterministic dynamics and ignoring stochasticity [30]. While the two measures coincide either when steady-state population density is large or when migration rate is large, the measures differ significantly where both density and migration rate are intermediate or low and sudden extinctions are important. Particularly, the migration rate that optimizes population abundance in a deterministic system may approach zero when the chance that a mutant arrives to its destination, *c*, decreases below some value [30], whereas in a stochastic system, the migration rate that optimizes population abundance decreases with *c*, but never approaches zero. Another possible measure is the MTE [20,22]. Similarity with our measure, when sufficiently many sites are considered, the MTE is longer for migration rates that promote persistence in our model. Nevertheless, our measure is more general and applies in the absence of extinction as well (figure 2*b*). In summary, our measure of population persistence is relevant for both small and large migration rates, and thereby it bridges a wide range of ecologically relevant scenarios.

## Funding statement

This research was supported by NSF grant no. DEB 1009957 to A.H.

## Acknowledgements

We thank Nadav Shnerb, Baruch Meerson and four anonymous reviewers for helpful comments on the manuscript, and Sebastian Schreiber for helpful discussions.

- Received June 23, 2013.
- Accepted July 10, 2013.

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