Limit theorems for discrete-time metapopulation models

We describe a class of one-dimensional chain binomial models of use in studying metapopulations (population networks). Limit theorems are established for time-inhomogeneous Markov chains that share the salient features of these models. We prove a law of large numbers, which can be used to identify an approximating deterministic trajectory, and a central limit theorem, which establishes that the scaled fluctuations about this trajectory have an approximating autoregressive structure.


Metapopulations
A metapopulation is a population confined to a network of geographically separated habitat patches that may suffer extinction locally and be recolonized through dispersal of individuals from other patches.The term was coined by Levins [41], but the idea goes back much further, to MacArthur and Wilson [42], Andrewartha and Birch [4] and Wright [67], and has been refined more recently by Hanski and others (see for example [26] and [27]).Levins [40] was the first to provide a succinct mathematical description of a metapopulation, proposing that the number n t of occupied patches at time t in a group of N patches should follow the law of motion with c being the colonization rate and e being the local extinction rate.This is Verhulst's model [63] for population growth and Levins used Pearl's rationale [51,52,53] to derive it.Furthermore, Levins was able to divine an explicit solution to (1) in the case where both c and e are time dependent, and he derived a diffusion approximation for n t (surprisingly, the time-inhomogeneous version of Levins' model has had no traction in the ecology literature, despite the obvious implications for biological control in varying environments).
The natural stochastic analogue of ( 1) is a birth-death process on S = {0, 1, . . ., N } with birth rates λ n = λn(N − n), where λ = c/N , and death rates µ n = en (in the usual notation); indeed (1) can be viewed as its large-N mean field limit [54].Known as the SIS (susceptible-infectious-susceptible) model in the epidemiology literature [64], it is an example of Feller's stochastic logistic model [23], which includes variants of the SIS model which are also of use in ecology, for example, models that account for colonisation from an external source [2,57].A significant drawback of the SIS model is that patches are assumed to be identical [50], although it does account for patch proximity (λ is inversely proportional to N ).A more immediate drawback is that seasonal variation is not taken into account.Seasonal variation could certainly be accommodated using a time-inhomogeneous version, allowing c and e to be time dependent, but a simpler approach is to suppose that the number of occupied patches follows a Markov chain in discrete time.Akçakaya and Ginzburg [1] and Day and Possingham [22] developed a Markov chain model that assumes colonisation and extinction events occur in distinct successive phases; one might envisage an annual cycle, with local populations being susceptible to extinction during winter, while new populations establish in empty patches during the spring.It is assumed that a census takes place either at the end of successive colonisation phases (EC model) or at the end of successive extinction phases (CE model) and the state of the Markov chain is the state of the population at a census time.This approach has become predominant in the applied metapopulation literature, because it provides a vehicle for parameter estimation [46] and permits control mechanisms to be investigated using simple optimisation tools such as dynamic programming [60,62].Indeed, discrete-time Markov chain models predominate in the ecology literature (even in cases where they are not faithful to population dynamics), perhaps due in part to a widespread misconception that a discrete time model is needed if populations are observed (and controlled) at discrete time points [59,58].Numerical methods and simulation are generally used to analyse discrete time metapopulation models, typically the EC case only [29,34,68], and until recently there have been few analytical studies [18,45].
Our purpose here is to present limit theorems that can be brought to bear in the study of chain-binomial metapopulation models (and indeed any discretetime Markovian model that shares their salient features).We present a law of large numbers, which is used to identify an approximating (discrete-time) deterministic trajectory, and a central limit theorem, which establishes that the scaled fluctuations about this trajectory have an approximating autoregressive structure.Limit theorems of this kind are standard fare in the context of continuous-time Markovian models (see Kurtz [35,36,37,38,39] and Barbour [11,12,13,14] and, more recently, Darling and Norris [20]); an approximating continuous-time deterministic trajectory is identified, and the fluctuations about this trajectory are approximated by a Gaussian diffusion.These methods have been exploited in the analysis of continuous-time metapopulation mod-els [2,5,6,7,15,56,57,54].Diffusion approximations can also be realised for discrete-time Markovian models [49], but these approximations do not respect the discrete-time structure.
We begin by describing a class of one-dimensional discrete-time metapopulation models, called homogeneous SPOMs (stochastic patch occupancy models) in the ecology literature-homogeneous because patches are assumed to be identical.Section 3 contains the basic limit theorems.Our major source is the work of Klebaner and Nerman [33] (see also Klebaner [32]), who studied a generalisation of the Galton-Watson process where the offspring distribution was allowed to depend on the current population size measured as a proportion of some threshold.We explain how their results can be extended to accommodate timeinhomogeneous discrete-time Markov chains with density dependent transition probabilities.These results are applied to our metapopulation models in Sections 4 and 5.We conclude by exploring the relationship between discrete-time metapopulation models and their counterparts in continuous time.

Discrete-time metapopulation models
Extinction and colonisation are assumed to occur in alternating Markovian phases.In the extinction phase occupied patches are assumed to go extinct independently, each with the same probability e (0 < e < 1), while in the colonisation phase empty patches are colonised independently, each with a probability c(x) that depends on the proportion x of occupied patches at the start of this phase.For example, we might have c(x) = cx, where c (0 < c ≤ 1) is the (hypothetical) probability that a single unoccupied patch would be colonised by the fully occupied network, or c(x) = c 0 + cx (c ≥ 0, c 0 > 0, c 0 + c ≤ 1) if we wish to account for an external source of colonisation.Or, we might have c(x) = 1 − exp(−βx) (β > 0), which effectively assumes (see [34,29]) that colonising individuals propagate from each occupied patch at rate β.However, we will assume only that c is continuous, increasing and concave, with c(0) ≥ 0 and c(x) ≤ 1.
A census is assumed to take place either at the end of successive colonisation phases (the EC model) or at the end of successive extinction phases (the CE model); we consider both scenarios.If N is the total number of patches and n t the number occupied at census t ∈ {0, 1, . . .}, then (n t , t ≥ 0) is a Markov chain taking values in S N = {0, 1, . . ., N } with the following chain-binomial structure: (We adopt the convention that Ber(p), Bin(n, p), Poi(µ) and N(µ, σ 2 ) denote random variables with the corresponding Bernoulli, Binomial, Poisson and Gaussian distribution.)Under the conditions we have imposed, S N is irreducible unless c(0) = 0, in which case there is a single absorbing state 0, corresponding to total extinction of the population, with the remaining states forming an irreducible transient class E N = {1, 2, . . ., N } from which 0 is accessible.
Models with a state-independent colonisation probability c(x) = c 0 (> 0) have appeared in different guise in the epidemiology literature (see Section 4.4 of [19]).They were studied by us in detail in [18].We found that n t has the same distribution as the sum of two independent binomial random variables, Bin(n 0 , p t ) and Bin(N − n 0 , q t ), whose success probabilities (p t ) and (q t ) could be exhibited explicitly.Indeed n . is equivalent (in law) to the urn model n t+1 = Bin(n t , p 1 ) + Bin(N − n t , q 1 ), so that p 1 and q 1 can be interpreted as 'effective' survival and colonisation probabilities.This special property of equivalent independent phases carries over to some extent in the general case: for the CE model (only) it is easy to show that n t+1 has the same distribution as the sum of two independent binomial random variables, Bin(n t , 1 − e) and Bin(N −n t , (1−e)c(n t /N )), so that (1−e)c(x) is the effective colonisation probability when the proportion of occupied patches is x.In [18] we also examined the proportion X N t = n t /N of occupied patches and considered what happens for fixed t as N gets large, assuming X N 0 → x 0 .We proved a law of large numbers that established existence of a limiting deterministic trajectory x .with initial value x 0 , which could be exhibited explicitly.The corresponding central limit law for the scaled fluctuations x t ) about this trajectory was also proved, assuming Z N 0 → z 0 , the limiting Gaussian distribution having a variance that could be exhibited explicitly, and we mooted that the process Z N .might converge (in the sense of finite-dimensional distributions) to a Gaussian Markov chain Z .with initial value z 0 .
These results can be extended to accommodate general discrete-time metapopulation models with state-dependent colonization probabilities.We will see that, under mild conditions, X N t P → x t for all t ≥ 0, where x . is determined by x t+1 = f (x t ) (t ≥ 0) with f specific to the model.It is also possible to identify conditions which ensure that if Z N 0 D → z 0 , then Z N .converges weakly (in the usual product topology) to the Gaussian Markov chain Z .defined by , where v is specific to the model.
We will also examine the following infinite-patch models, with n t ∈ S = {0, 1, . . .}: where m(n) (> 0) is the expected number of patches colonised during any one colonisation phase when the number occupied at the start of that phase is n.Even though there is now no ceiling on the number of occupied patches, the dependence of m on n t allows for regulation of the colonisation process.The case m(n) = mn (where the constant m > 0 is the expected number of colonisations by any one occupied patch) is a natural analogue of the N -patch models described above, for if c(0) = 0 and c has a continuous second derivative near 0, then, for fixed n, Bin(N − n, c(n/N ))

General structure: Density dependence
Let (n N t , t ≥ 0) be a family of Markov chains indexed by N , each taking values in a subset set S N of Z + .Suppose that the family is density dependent in that there are sequences of non-negative functions (f t ) and (v t ) such that Then, the 'density process' (X N t , t ≥ 0), defined by Our first result is a law of large numbers that establishes convergence of the density process to a deterministic trajectory.
Theorem 1. Suppose that, for all t ≥ 0, f t (x) and v t (x) are continuous in x and such that f t (X N t ) and v t (X N t ) are a.s.uniformly bounded.Then, if Proof.We will use mathematical induction.Suppose X ) is a.s.uniformly bounded, and so f t (X N t ) r → f t (x t ) (convergence in r-th mean) for all r ≥ 1, which entails, in particular, that But, for all ǫ > 0, Pr(|X , and the proof is complete. We anticipate applying Theorem 1 in cases where X N t itself is bounded, for example, when X N t is a proportion (as it is in our N -patch models).To accommodated cases where X N t is unbounded (as it can be in our infinite-patch models) we relax uniform boundedness in favour of a Lipschitz condition, but at the expense of requiring a more stringent initial condition, that X N 0 converges to x 0 in mean square.
Theorem 2. Suppose that, for all t ≥ 0, f t (x) and v t (x) are Lipschitz continuous in x.
Proof.We will again use mathematical induction.Suppose , for some positive constant κ t .On taking expectations, we see that ).We deduce that Var X N t+1 → 0, for, as noted earlier, Var and so Having established convergence in probability to a limiting deterministic trajectory x . ,we next consider the 'fluctuations process' (Z N t ) defined by → z 0 , we aim to identify conditions under which Z N .converges weakly to a Gaussian Markov chain Z . .Additional structure is needed.We will assume that where r N t = N r t (n N t /N ) and g N t = N g t (n N t /N ) with r t (x) and g t (x) being continuous in x, and ξ N jt (j = 1, . . ., r N t ) are iid having a distribution that depends only on t and on n N t /N , and which has bounded third moment.In particular, we assume that there are functions m t (x) and Of course all of our ingredients must be such that r N t and g N t are positive integers and, then, that n N t+1 ∈ S N .Notice that E(n N t+1 |n N t ) = N f t (n N t /N ), where f t (x) = g t (x) + r t (x)m t (x), and Var (n N t+1 |n N t ) = N v t (n N t /N ), where v t (x) = r t (x)σ 2 t (x).This setup is similar to that of Klebaner and Nerman [33] (see also Klebaner [32]) who studied a generalisation of the Galton-Watson process where the offspring distribution was allowed to depend on the current population size measured as a proportion of some threshold N .Their model was time homogeneous and had g t ≡ 0 and r t (x) = x.None the less, many of the results in Section 3 of their paper carry over to the present case with only minor changes.We content ourselves with the following central limit law.Theorem 3. Suppose that, for all t ≥ 0, f t (x) is twice continuously differentiable in x with bounded second derivative and that X N t P → x t , where x .satisfies . converges weakly to the Gaussian Markov chain Z .defined by where (E t ) are independent with E t ∼ N(0, v t (x t )).
Proof.First observe that we may rewrite (2) as ([ • ] denotes integer part), noting that, for fixed x, η N t (x) is independent of X N t .Next, for fixed t and x, η N t (x) D → N(0, v t (x)) as N → ∞.To see this, apply the Feller-Lindeberg Theorem (see for example Theorems 27.2 and 27.3 of Billingsley [16]) to the zero-mean triangular array (W N j , j = 1, . . ., r N ), where W N j = ξ N jt − m t (x) and r N = [N r t (x)], noting that E(W 2 N j ) = σ 2 t (x) and the Lyapounov condition, that → 0 (for some δ > 0), is satisfied here with δ = 1 because b t (x), the third centred moment of ξ N jt , is bounded in x.
We have assumed that X N t P → x t for all t ≥ 1, where x . is determined by we also have, by Taylor's theorem, that for some θ N t between X N t and x t , and so, from (4), Since f ′′ t (x) is bounded in x, we may thus write where P → 0 as N → ∞.To establish weak convergence of Z N .to Z .it is sufficient to establish conver- gence of the finite-dimensional distributions.To this end, consider the characteristic function Since, for fixed t and x, η N t (x) this iteration defines the characteristic function of (Z 0 , . . ., Z t ), where Z . is the proposed limiting Gaussian process.This completes the proof.
Remarks.(i) An alternative approach to proving Theorems 1 and 3 would be to adapt the results of Karr [30] to a time inhomogeneous setting.He considered a sequence of time-homogeneous Markov chains (X N t , t ≥ 0) on a general state space (E, E) N with transition kernels (K x in E, then the corresponding probability measures (P N πN ) on (E, E) N also converge: P N πN ⇒ P π .Karr's main result (Theorem 1 of [30]) remains true, with obvious modifications, for a time inhomogeneous Markov chain.Applying this to our density process (X N t ), and then to the two-dimensional process (X N t , Z N t ), would establish the required convergence, and the form of the limiting processes would become apparent once the limiting kernels were evaluated.
(ii) The chain binomial models described in Section 2 are time homogeneous, yet we need the present level of generality to accommodate them because it is not always possible to establish density dependence, even on occasions when we anticipate the kind of asymptotic (large N ) behaviour exhibited in Theorems 1 and 3.For our models (or, more generally, for Markov chains with two alternating phases), it is natural to construct a time-inhomogeneous Markov chain (n N t , t ≥ 0) by setting n N 2t = n t and n N 2t+1 = ñt .Then, density dependence in the phases is often enough to establish the required asymptotic behaviour.Of course, this programme extends to models with more than two density-dependent phases.For example, our methods can easily accommodate the (three-phase) extinction-reproduction-settlement model of Klok and De Roos [34].
(iii) Klebaner and Nerman (Theorem 3 via Theorem 6 of [33]) anticipated a special case of Theorems 1 and 3: when the time homogeneous version of ( 4) is in force with r t (x) = x.
(iv) The mean and covariance function of Z .are easy to evaluate by iterat- ing (3): and where and (Here and henceforth empty products are to be interpreted as being equal to 1.) Furthermore, it is clear that, for any t ≥ 1, , and so these formulae can be used to approximate the mean and covariance function of n N . .Indeed, the joint distribution of n N t1 , . . ., n N tn , where t 1 , . . ., t n is any finite set of times, can be approximated by an n-dimensional Gaussian distribution with . is time homogeneous, and hence the approximating deter- ministic model has the form x t+1 = f (x t ), the full range of long-term behaviour, including chaos, is possible.For example, if f has a stable fixed point x * and X N 0 P → x * , then, assuming f and v are continuous, and f (X N t ) and v(X N t ) are a.s.uniformly bounded, we will have ) and, if f is twice continuously differentiable with bounded second derivative, then, assuming Z N 0 D → z 0 , the limit process Z .will be an au- toregressive (AR-1) process with the representation and |a| < 1, then there will be a sequence of times (tN) such that More generally, if f admits a stable limit cycle, x * 0 , x * 1 , . . ., x * d−1 , and The distribution of Y 0 , and both the coefficient matrix A and the covariance matrix Σ, are determined using ( 6)-( 9) (obtained by iterating (3)) with x nd+j = x * j (n ≥ 0, j = 0, . . ., d − 1).This was done explicitly by Klebaner and Nerman (Theorem 4 of [33]) for the population dependent branching process, and observed to be true more generally (Theorem 6 of [33]).Indeed, their Theorems 4 and 5 hold word for word, but with obvious adjustments in the definitions of f and v, in the present more general context of a time-homogeneous density dependent family.

N -patch models
Here we explain how the results of Section 3 can be used to study our N -patch metapopulation models.We first identify an approximating deterministic model and describe the structure of the approximating Gaussian process, and then give a detailed account of long-term behaviour.

Limit theorems for the proportion of occupied patches
Let X N t = n t /N be the proportion of occupied patches at census t.We first prove a law of large numbers for the process X N . ,thus establishing the existence of an approximating deterministic trajectory, and a central limit law for the fluctuations about this trajectory.For the N -patch CE model, we may evaluate E(n t+1 |n t ) and Var (n t+1 |n t ) explicitly, because E(n t+1 |ñ t ) and Var (n t+1 |ñ t ) are both linear functions of ñt : = e E(n implying that We may apply Theorem 1 directly since both f and v are continuous, and f (X N t ) and v(X N t ) are bounded because 0 ≤ X N t ≤ 1 and c(x) ≤ 1.For the N -patch EC model we cannot evaluate the conditional mean and variance explicitly, and it will be clear that the model is not density dependent unless we impose further restrictions on c.For the N -patch EC model (and indeed, alternatively, for the CE model) we apply Theorem 1 to the time-inhomogeneous Markov chain (n N t , t ≥ 0) obtained by setting n N 2t = n t and n N 2t+1 = ñt .Then, for the EC model, All of these functions are continuous, and f t (X N t ) and v t (X N t ) are uniformly bounded because 0 ≤ X N t ≤ 1 and c(x) ≤ 1.The limiting trajectory satisfies (in particular) x 2(t+1) = f (x 2t ), where f = f 1 • f 0 .Thus we have the following simple result.Theorem 4. For the N -patch metapopulation models with parameters e and c(x), let X N t be the proportion of occupied patches at census t.
To obtain the corresponding central limit law for , first observe that our N -patch models can be represented as )), (CE model) where the (Ber j (p)) are collections of iid Bernoulli random variables with success probability p.We may thus apply Theorem 3 to the time-inhomogeneous Markov chain (n N t , t ≥ 0) obtained (as above) by setting n N 2t = n t and n N 2t+1 = ñt , because this chain will have the form (2) with the (± ξ N jt ) being appropriate sequences of iid Bernoulli random variables.For both models, and Thus, for the CE model the roles of f 0 and f 1 , and v 0 and v 1 , are reversed.If c is twice continuously differentiable, then (remembering that c is increasing and concave with c(0) ≥ 0 and c(x) ≤ 1) in both cases f t (x) will be twice continuously differentiable in x with bounded second derivative, v t (x) will be continuous in x, and b t (x) will be bounded in x.
We have already seen that our deterministic trajectory satisfies x 2(t+1) = f (x 2t ), where f = f 1 • f 0 for the EC model and f = f 0 • f 1 for the CE model, with f 0 and f 1 as given in (10).Similarly, it is clear that our limiting Gaussian Markov chain Z .must satisfy 1 for the CE model, with v 0 and v 1 as given in (11), noting that Thus we arrive at the following result.Theorem 5.For the N -patch metapopulation models with parameters e and c(x), suppose that c is twice continuously differentiable.Let , where X N t is the proportion of occupied patches at census t and where x . is determined by x t+1 = f (x t ) (t ≥ 0) with f given as in Theorem 4.Then, if . converges weakly to the Gaussian Markov chain Z .defined by Since our limiting Gaussian Markov chain Z .satisfies (3), now with f given as in Theorem 4, an immediate consequence of Theorem 5 is that where (even though the constructed Markov chain n N . is time-inhomogeneous).Also, the joint distribution of numbers of occupied patches, observed at census times t 1 , . . ., t n , can be approximated by an n-dimensional Gaussian distribution with means N x ti + √ N µ ti and covariances N c ti, tj , where c t, s

Long-term behaviour
Next we look at stationarity/quasi stationarity, and begin by examining the long-term (t → ∞) behaviour of our deterministic models.First notice that, in an obvious notation, f CE ((1 − e)x) = (1 − e)f EC (x), and so the fixed points of f CE and f EC are related by x * CE = (1 − e)x * EC (again in an obvious notation).Furthermore, x * CE and x * EC will have the same stability properties, because , where r(x) = ρ x/(1 − x) and ρ = e/(1 − e), the function r having slope ρ at x = 0 and increasing strictly from 0 to ∞ as x increases from 0 to 1. But, recall that c is strictly increasing from c(0) ≥ 0 and concave with c(x) ≤ 1.Therefore, we always have precisely one stable fixed point and x t approaches this point monotonically.We have the following three cases.
(i) Stationarity: c(0) > 0. There is a unique fixed point x * in [0, 1] and this satisfies 0 < x * < 1. Moreover it is stable, because clearly c ′ (x * CE ) < r ′ (x * CE ), and a simple calculation shows that this entails (iii) Quasi stationarity: c(0) = 0 and c ′ (0) > ρ.There are two fixed points in [0, 1], 0 and x * ∈ (0, 1); 0 is unstable because now f ′ CE (0) = (1 − e)(1 + c ′ (0)) > 1 and x * is stable by the same argument as used in (i).Thus, in Case (i) we expect the unique stationary distribution of the process n N . to be centred near N x * .In Cases (ii) and (iii) there is an absorbing state 0 which is reached in a finite time.However, in Case (ii) we expect the process to be absorbed quickly (even for N quite large), while in Case (iii) a quasi-equilibrium will be reached; we expect the unique quasi-stationary distribution of n N . (being the limiting conditional state probabilities, conditional on non-extinction) to be centred near N x * .
Notice that x * CE < x * EC .Indeed, since f CE (x) < f EC (x), the deterministic trajectory x .will be uniformly smaller for the CE model than for the EC model.This is to be expected, for our models differ only in when the census is taken; numbers observed following an extinction phase would likely be smaller than numbers observed following a colonisation phase.These remarks are supported by illustrations in Figure 1.Simulations are depicted for both the EC and CE models with c(x) = cx (see Example 1 below), as well as the corresponding deterministic trajectories and quantities relating to the limiting Gaussian processes.The quasi-stationary distribution, . was eval- uated as the normalized left eigenvector of the transition matrix restricted to E N corresponding to its Perron-Frobenius eigenvalue (see Darroch and Seneta [21]), and this was compared with the approximating Gaussian pdf with mean N x * and variance N V * , where V * = v(x * )/(1 − f ′ (x * ) 2 ) (see Corollary 1 below).We note that, whilst the quasi-stationary distributions for our two models cannot be exhibited explicitly, they are always related by π EC = π CE C (in an obvious notation) with common Perron-Frobenius eigenvalue, where C denotes the matrix C (the colonisation transition matrix) restricted to E N .
Our next result is obtained from Theorems 4 and 5 on setting x 0 = x * .It established that in the stationary and quasi-stationary cases, where there is a positive stable deterministic equilibrium x * , the fluctuations . about x * can be approximated by an AR-1 process whose parameters can be exhibited explicitly.
Corollary 1.For the N -patch metapopulation models with parameters e and c(x), let X N t be the proportion of occupied patches at census t and let Suppose that, in addition to c being twice continuously differentiable, c(0) > 0 or c(0) = 0 and c ′ (0) > e/(1 − e), and let x * be the stable fixed point of f (f given as in Theorem 4).Then, . converges weakly to the AR-1 process Z .defined by Z t+1 = f ′ (x * )Z t + E t (Z 0 = z 0 ), with iid errors , where v is given as in Theorem 5.
One consequence of the corollary is that, for all t ≥ 1, then there will be a sequence of times (tN) such that Z N tN D → N(0, V * ), where V * = v(x * )/(1 − a 2 ).Also, we expect that if the process has reached equilibrium/quasi equilibrium then the joint distribution of the numbers of occupied patches, observed at census times t 1 , . . ., t n , can be approximated by an n-dimensional Gaussian distribution with means N x ti + √ N µ ti and covariances N c ti, tj , where c t, s := Cov (Z t , Z s ) = V t a |s−t| .It would be of interest to determine how closely, for how long, and over what ranges, X N .is faithfully ap- proximated.To this end, we might look at the time τN = inf{t ≥ 1 : . from an interval containing x * , where eN → ∞.Based on Theorem 1 of Barbour [12] (who considered this problem for the continuous time analogue-a limiting Ornstein-Uhlenbeck process), we conjecture that if eN does not grow too quickly, X N .is asymptotically equally likely to leave to the right of the interval as to the left, and, conditional on (say) leaving to the right, the exit time is asymptotically geometrically distributed.
We have already noted the simple relationship between the deterministic equilibria of our two models, x * CE = (1 − e)x * EC , and that the decay rates are the same: . The stationary variances of the approximating AR-1 processes are also related.First, because c(x * ) = r(x * ), it is easy to prove that And, since it can also be shown that and therefore (1

Examples
We now illustrate these results by looking at particular instances of c(x).
Example 1. Suppose that c(x) = cx (0 < c ≤ 1).In this case we may write f (x) = x(1 + r(1 − x/x * )), where r = c(1 − e) − e for both models and x * is the appropriate equilibrium: x * EC = r/(c(1 − e) 2 ) or x * CE = r/(c(1 − e)) (both being strictly positive, and then stable, if and only if c > e/(1 − e)).Thus, our limiting deterministic trajectory follows the discrete logistic model (see for example Section 3.2 of Renshaw [55]), with r being the 'natural growth rate' and x * being the 'carrying capacity' (expressed as a proportion of the ceiling N ).Of course the logistic model is well known to exhibit a wide range of dynamic behaviour, but we emphasise that here 0 < 1 + r = (1 − e)(1 + c) < 2. Our limiting Gaussian Markov chain has error variance (CE model) In the quasi-equilibrium case (r > 0), the limiting AR-1 process is defined by Z t+1 = aZ t + E t , where a = 1 − r (0 < a < 1), with E t ∼ N(0, v * ), where, from ( 12) and ( 13 that this decreases with a, so the faster the decay in the mean, the smaller the stationary variance. Example 2. Suppose that c(x) = c 0 , where 0 < c 0 ≤ 1.This case was studied in detail by us in [18].As mentioned earlier, the metapopulation behaves as if, at every census, each occupied patch remains occupied with probability p, and, independently, each unoccupied patch is colonised with probability q, where Notice that for the EC model the 'effective' extinction probability is e(1 − c 0 ).This accords with Hanski's [25] interpretation of the 'rescue effect', a term coined by Brown and Kodric-Brown [17] to describe the deleterious effect of colonisation on extinction when colonisation is frequent; Hanski argued that the extinction probability should be (1 − c 0 )e, where e is the extinction probability in the absence of migration.
In [18] we proved that, for all t ≥ 1, n N t has the same distribution as the sum of two independent random variables, Bin(n N 0 , p t ) and Bin(N − n N 0 , q t ), with success probabilities q t = q * (1 − a t ) and p t = q t + a t (t ≥ 0), where a = p − q = (1 − e)(1 − c 0 ) (the same for both EC and CE) and q * = q/(1 − a).The proportion X N t of occupied patches at time t has mean and variance given by So, on the one hand, as t → ∞, EX N t converges to q * at geometric rate a (note that 0 < a < 1) and N Var (X N t ) → q * (1 − q * ) (indeed, n N . has a Bin(N, q * ) stationary distribution).On the other, letting N → ∞ with t fixed, EX N t → x t (x 0 ) and N Var ( t is the sum of two independent binomial random variables, it is clear that will converge in distribution to Gaussian random variables if their initial values converge.Theorem 5 and Corollary 1 provide more detailed information.Since f (x) = px + q(1 − x), and hence f ′ (x) = a, and v . converges weakly to an AR-1 process Y .with EY t = a t y 0 and Cov (Y t , Y s ) = q * (1 − q * )a |s−t| (1 − a 2t ); the error variance here is v(q * ) = q * (1 Finally, we remark that if n N 0 follows the stationary Bin(N, q * ) law, our representation n t+1 = Bin(n t , p 1 ) + Bin(N − n t , q 1 ) is termed binomial autoregressive [43,44,65,66].Our results establish a connection with standard autoregressive processes.The error variance can be evaluated, but omitted for brevity's sake.The limiting AR-1 process has stationary variance v * /(1 − a 2 ), where v * is given by ( 12) or (13).In the quasi-stationary case the deterministic equilibria cannot be exhibited explicitly, but can be evaluated numerically by iterating the map The limiting AR-1 process has stationary variance v * /(1 − a 2 ), where v * is evaluated using (12) or (13).A simple calculation reveals that .

Infinite-patch models
As previously, let n t be the number of occupied patches at time t, but suppose now that (n t , t ≥ 0) is a Markov chain taking values in S = {0, 1, . . .} that evolves as follows: where m(n) ≥ 0. So, as before, extinction and colonisation occur in alternating phases, and occupied patches go extinct independently with probability e (0 < e < 1), but now the number of colonisations follows a Poisson law and the expected number of colonisations is a function of the number of patches presently occupied.
Before embarking on the general case, let us examine the important special case where the expected number of colonisations is a linear function of the number of patches presently occupied.

The infinite-patch branching model
Suppose that m(n) = mn, where m > 0. The parameter m can be interpreted as the expected number of colonisations by any one occupied patch.As we noted earlier, this is the natural analogue of our N -patch models, for recall that, if c(0) = 0 and c has a continuous second derivative near 0, then Bin(N − n, c(n/N )) D → Poi(mn) as N → ∞, where m = c ′ (0).But, this infinite-patch scheme has a simplifying feature, namely branching, which makes both models much simpler to analyse.Notice that if there are n occupied patches at the beginning of any given phase, then the number occupied at the end of that phase has the same distribution as the sum of n independent copies of either B = Ber(1 − e) (extinction phase) or P + := 1 + Poi(m) (colonisation phase).And, since the phases are conditionally independent, the net effect is that n t+1 will have the same distribution as the sum of n t independent copies of Y , where Y is either B independent copies of P + (the EC model) or P + independent copies of B (the CE model).We therefore make the following simple observation.
Proposition 1.The process (n t , t ≥ 0) is a Galton-Watson process whose offspring distribution has pgf G(z) given by Thus, we think of the census times as marking the 'generations' of our branching process, the 'particles' being the occupied patches, and the 'offspring' being the occupied patches that they notionally replace in the succeeding generation.For the EC model, it is as if each occupied patch becomes empty with probability e, or otherwise colonises Poi(m) patches, while, for the CE model, it is as if each occupied patch survives with probability 1 − e or becomes empty with probability e, but, whatever happens, it colonises Poi(m(1 − e)) patches.
We may now invoke the encylopaedic theory of branching processes [8,9,10,28] to prove results for this important special case of the model; it is just a matter of which questions are of interest.For example, it is easy to prove that offspring distribution has mean µ = (1 + m)(1 − e) (the same for both models), and so E(n t |n 0 ) = n 0 µ t (t ≥ 1).Also, our branching process is subcritical , critical or supercritical according as m is less than, equal to or greater than the critical value ρ = e/(1 − e).This accords immaculately with our earlier criteria for evanescence (c ′ (0) ≤ ρ) versus quasi stationarity (c ′ (0) > ρ) of our N -patch models with c(0) = 0. We also have the following simple result concerning the probability that the metapopulation becomes extinct (totally extinct), starting with n 0 patches occupied.
Corollary 2. For the infinite-patch branching model, total extinction occurs with probability 1 if and only if m ≤ ρ; otherwise total extinction occurs with probability η n0 , where η is the unique fixed point of G on the interval (0, 1), with G given as in Proposition 1.
The extinction probability η cannot be exhibited explicitly, but can of course be obtained numerically by iterating the map G.
Finally, the variance of the offspring distribution is , depending on which model, and so A simple extension of the branching model is obtained by setting m(n) = m 0 + mn, where now m ≥ 0, and the new parameter m 0 (> 0) is to be interpreted as the expected number of colonisations from an external source.It can be derived from our N -patch models if they are modified so that colonisation probability c(n/N ) is replaced by m 0 /N + c(n/N ) (imagine that an external colonization potential m 0 is apportioned equally among all N patches), for then Bin The effect is to introduce an additional (independent) Poi(m 0 ) number of colonisations in each colonisation phase.It is easy to see that the resulting process (n t , t ≥ 0) will be the Galton-Watson process identified in Proposition 1, but modified so that there are Poi(d) immigrant particles in each generation, where d = m 0 for the EC model and d = (1 − e)m 0 for the CE model.Again we can invoke general theory.
For example, on applying Theorem VI.7.2 of [10], we learn that n .has a proper limiting distribution if and only if m < ρ.

The infinite-patch model with regulated colonisation
Returning now to the general case, where the expected number of colonisations depends arbitrarily on the number of occupied patches, let us consider what First observe that our infinite-patch models have the equivalent representation where the (Poi j ( • )) are collections of iid Poisson random variables with mean µ(n t /N ).Thus, mirroring our argument leading to Theorem 5 for our N -patch models, we may apply Theorem 3 to the time-inhomogeneous Markov chain (n N t , t ≥ 0) obtained by setting n N 2t = n t and n N 2t+1 = ñt .This chain also has the form (2), but with the (± ξ N jt ) now being appropriate sequences of iid Poisson random variables.For both models, and , leading to the same expressions for f t , v t and b t , but with the roles f 0 and f 1 , v 0 and v 1 , and b 0 and b 1 , reversed.Since µ is twice continuously differentiable with bounded second derivative, in both cases f t (x) will be twice continuously differentiable in x with bounded second derivative, v t (x) will be continuous in x, and b t (x) will be bounded in x.
We have just seen that the limiting deterministic trajectory satisfies x 2(t+1) = f (x 2t ), where f = f 1 • f 0 for the EC model and f = f 0 • f 1 for the CE model, with f 0 and f 1 as given in (14).Similarly, it is clear that our limiting Gaussian Markov chain Z .should take the form 1 for the CE model, with v 0 and v 1 as given in (15), noting that f ′ 0 (x) = 1 − e and f ′ 1 (x) = 1 + µ ′ (x).Thus we arrive at the following result.Theorem 7.For the infinite-patch metapopulation models with parameters e and µ(x), suppose that µ is twice continuously differentiable with bounded second derivative.Let X N t be the proportion of occupied patches at census t, and suppose that X N 0 2 → x 0 , so that X N t P → x t for all t ≥ 0, where x . is determined by . converges weakly to the Gaussian Markov chain Z .defined by Z t+1 = f ′ (x t )Z t + E t (Z 0 = z 0 ), with (E t ) independent and Equilibrium behaviour is richer and more interesting than for the earlier Npatch models, because now the limiting deterministic models can exhibit the full range of long-term behaviour, and we cannot be as precise in classifying this behaviour as we were earlier.First notice that, in the notation adopted earlier, f CE ((1 − e)x) = (1 − e)f EC (x), and so fixed points of f CE and f EC are related by x , where ρ = e/(1 − e).So, if µ(0) = 0 then 0 is a fixed point; it is stable if µ ′ (0) < 1 and unstable if µ ′ (0) > 1 (if µ ′ (0) = 1 its stability is determined by higher derivatives of µ near x = 0).However, even when µ(0) = 0, there might be other (conceivably many) fixed points; our conditions on µ do not preclude this.Certainly if there is a unique positive fixed point x * , it will be stable if µ ′ (x * ) < 1 and unstable if µ ′ (x * ) > 1 (again we need to consider higher derivatives when µ ′ (x * ) = 1).Finally, notice that the d-th iterates of our maps are also related by f Ricker growth dynamics.Suppose that µ(x) = x exp(r(1 − x)), where r > 0, so that the colonisation potential of the occupied patches is greatest when their number is close to N/r; the parameter r can be interpreted as the growth rate and N the carrying capacity of the metapopulation in the absence of extinction., implying that f ′′ CE (0) = −2(1 − e)re r .Therefore, if r ≤ r 0 , 0 is the unique non-negative fixed point, and it is stable.If r > r 0 , then x * CE is an additional positive fixed point; it is stable because f ′ CE (x * CE ) = 1 − e(r − r 0 ) < 1 (and 0 is unstable).However, if r is sufficiently large, we get limiting cycles with period doubling towards chaos, as illustrated in Figure 2.
Figure 3 illustrates some of the range of behaviour exhibited by the stochastic CE model with Ricker growth dynamics.Four cases are depicted.Evanescence (a): for r between 0 and r 0 (≃ 0.8473), corresponding to 0 being the unique stable fixed point, the process dies out quickly.Quasi stationarity (b): for r between r 0 and r 1 ≃ 0.7, corresponding to x * CE (≃ 1 − 0.8473/r) being the unique stable fixed point, the process exhibits (quasi-) equilibrium behaviour; r 1 is the point of first period-doubling (we have not been able to determine r 1 analytically).Oscillation (c) and (d): for r bigger than r 1 the process 'tracks' the limit cycles of the deterministic model (period 2 and 4, respectively).To be sure, the theory predicts the kind of behaviour shown in Figure 3(c), but it is no less remarkable when one witnesses it; one could not mistake it for, say, bistability.→ x * 0 .Corollary 3 follows from Theorem 7 on setting x 0 = x * , so that then x t = x * for all t ≥ 0, and evaluating the stationary error variance v(x * ) in both cases.Corollary 4 follows from Theorem 7 on setting x 0 = x * 0 , so that then x .tracks the limit cycle, that is, x nd+j = x * j (n ≥ 0, j = 0, . . ., d − 1).The representation of Z .as a d-variate AR-1 process Y . ,and in particular the form of the coefficient matrix A and the error covariance matrix Σ d , follow by iterating Z t+1 = f ′ (x t )Z t + E t (Z 0 = z 0 ), with (E t ) independent N(0, v(x t )) random variables: using expressions (6) to (9) with we obtain a representation of Z nd+j (j = 1, . . ., d) in terms of Z nd (n ≥ 0), as well as the stationary covariance matrix V .Notice that, because Z . is Markovian, only Z (n+1)d−1 contributes to the drift in Y n .

Continuous-time analogues
When extinction and colonisation events happen in random order, rather than in alternating phases, it is natural to take (n t , t ≥ 0), where n t is the number occupied patches at time t, to be a Markov chain in continuous time.But, what is the most appropriate model?If the probability of colonisation in a small time interval were independent of the number of occupied patches, then the SIS model would seem to be the most appropriate N -patch model, and the Levins model (1) could be used, in much the same way as above, to draw conclusions about its long-term behaviour.Evanescence and quasi stationarity could be distinguished by examining the stability of the equilibrium points, 0 and n * = N (1 − e/c); if c ≤ e, 0 would be stable and the population would have genuine evanescent character, while if c > e, n * would be strictly positive and stable, and the population would persist (in this latter case the quasi-stationary distribution would be centred near n * [47,48]).
If, as envisaged here, the probability of colonisation were to depend on the proportion of patches currently occupied, the natural N -patch continuous-time model would be a birth-death process on S = {0, 1, . . ., N } with birth rates λ n = c(n/N )(N − n) and death rates µ n = en, where c(x) is as above (assumed to be continuous, increasing and concave, with c(0) ≥ 0 and c(x) ≤ 1).To see this, suppose that an occupied patch becomes empty in a time interval of length h with probability eh + o(h), or remains occupied with probability 1 − eh + o(h), while if there are n occupied patches at time t, then any given unoccupied patch becomes occupied in the interval (t, t + h] with probability c(n/N )h + o(h), or remains unoccupied with probability 1 − c(n/N )h + o(h).Suppose also that the chance of two or more events of either kind happening in (t, t + h] is o(h).Then, a transition from n to n + 1 in time h is effected by having exactly one colonization and no extinctions (there are other ways, but all have probability o(h)), and the chance of this happening in time h is Similarly, and as well as Pr(n t+h = m|n So, analogous to Examples 1, 2 and 3 above, if c(x) = cx we get the stochastic SIS model, if c(x) = c 0 we get the continuous-time Ehrenfest model (see Section 1.4 of [31]), while if c(x) = c 0 + cx we get the mainland-island model of Alonso and McKane [2] (see also [57]), all being instances of Feller's stochastic logistic model [23].Whatever the form of c(x), we can obtain continuous-time analogues of Theorems 4 and 5 and Corollary 1, because our birth-death model is density dependent in the sense of Kurtz [35,36].It follows immediately from Theorem 3.1 of [35] that the proportion X N t = n t /N of occupied patches at time t converges (uniformly in probability over finite time intervals) to a deterministic trajectory (x t , t ≥ 0) satisfying the law of motion dx dt = F (x) (t ≥ 0), where F (x) = c(x)(1 − x) − ex (0 ≤ x ≤ 1), (16) assuming of course that X N 0 → x 0 (our conditions on c imply that F is Lipschitz continuous on [0, 1]).Furthermore, if we let Z N t = √ N X N t − x t , then, assuming Z N 0 → z 0 , Theorem 3.5 of [36] can be used to show that process (Z N t , t ≥ 0) converges weakly in D[0, t] (the space of right-continuous left-limits functions on [0, t]) to a Gaussian diffusion (Z t , t ≥ 0) with initial value Z 0 = z 0 and with mean E Z t = M t z 0 , where M t = exp( t 0 B u du) and B t := F ′ (x t ), and variance V t = M 2 t t 0 M −2 u G(x u ) du, where G(x) = F (x) + 2ex.In the important special case where x 0 is taken to be an equilibrium point x * of ( 16), usually a stable equilibrium, the approximating diffusion is an Ornstein-Uhlenbeck process, and more precise results are available [61].For example, if F (x * ) = 0 and B := F ′ (x * ) < 0, then E Z t = e −at z 0 , where a = −B, and V t = V * (1 − e −2at ), where V * = G(x * )/(−2B) = ex * /a (the stationary variance).
There is a one to one correspondence between the equilibria of ( 16) and the fixed points  ).Thus, stationarity happens when c(0) > 0, evanescence when c(0) = 0 and c ′ (0) ≤ e, and quasi stationarity when c(0) = 0 and c ′ (0) > e.One might have expected that, in stationary and quasi-stationary regimes, the stable equilibrium would lie between x * EC and x * EC , but this turns out not to be the case.
To illustrate this, and compare other equilibrium characteristics, suppose that c(x) = cx.We will assume that ρ < c ≤ 1, so that also c > e.Then,   )) and V * SIS = e/c (in an obvious notation).The natural continuous-time analogue of our infinite-patch models is the birth-death process on S = {0, 1, . . .} with birth rates λ n = m(n) and death rates µ n = en.When m(n) = mn, where m > 0, we get the simple immigrationdeath process (see Section 6.1 of [24]), which can be interpreted as a Markov branching process with binary splitting.Even with the inclusion of an immigration term, m(n) = m 0 + mn, the model is completely tractable (explicit results are catalogued in Section 3.2 of [3]).If, more generally, m(n) = N µ(n/N ), for some threshold N , where µ is continuous with bounded first derivative (implying that F is Lipschitz continuous), Theorem 3.1 of [35] guarantees that if the number X N t = n t /N of occupied patches at time t, measured relative to the threshold, converges (uniformly in probability over finite time intervals) to a deterministic trajectory (x t , t ≥ 0) satisfying dx dt = F (x) (t ≥ 0), where F (x) = µ(x) − ex (x ≥ 0), (17) whenever X N 0 → x 0 .Furthermore, Theorem 3.5 of [36] guarantees that (Z N t , t ≥ 0), where Z N t = √ N X N t − x t , converges weakly in D[0, t] to a Gaussian diffusion (Z t , t ≥ 0) with initial value Z 0 = z 0 , assuming Z N 0 → z 0 .Its mean is E Z t = M t z 0 , where M t = exp( t 0 B u du) and B t := F ′ (x t ) = µ ′ (x t ) − e, and its variance is u G(x u ) du, where G(x) = F (x) + 2ex.When x 0 is taken to be a stable equilibrium point x * of (17) Z . is an Ornstein-Uhlenbeck process with E Z t = e −at z 0 , where a = e − µ ′ (x * ), and V t = V * (1 − e −2at ), where V * = ex * /a is the stationary variance.It will be clear that the equilibria of ( 17), and their classification, will be the same as that laid out in Section 5.2 for the discrete-time models, but with ρ replaced by e.

Fig 1 .
Fig 1. N -patch metapopulation model with N = 30, e = 0.2 and c(x) = cx.Simulation (solid circles) of the EC Model with (a) c = 0.2 (evanescence) and (b) c = 0.6 (quasi stationarity); deterministic trajectories are shown (solid), together with ±2 standard deviations of the Gaussian approximation (dotted).(c) Simulation of the EC (solid circles) and CE (open circles) models with c = 0.6; both deterministic trajectories shown (solid).(d) Quasi-stationary distribution (bars) of the EC Model with c = 0.6 and the stationary Gaussian pdf (dotted).

Fig 3 .
Fig 3. Simulation (open circles) of the infinite-patch CE model with Ricker growth dynamics, together with the corresponding limiting deterministic trajectories (small solid circles).Here e = 0.7 and N = 200, and, (a) r = 0.84, (b) r = 1 (c) r = 4 and (d) r = 5.In (a), (b) and (c), the dotted lines indicate ±2 standard deviations of the Gaussian approximation (in (c) every second point is joined to indicate variation about each of the two limit cycle values).

Fig 4 .
Fig 4. N -patch metapopulation model with N = 30, e = 0.2 and c(x) = cx.(a) x * CE (solid), x * SIS (dotted) and x * EC (dashed) versus c.(b) pdf of the stationary distribution of the limiting Gaussian process for the CE model (solid), the SIS model (dotted) and the EC model (dashed) with c = 0.6, together with the corresponding quasi-stationary distributions (open circles, diamonds and solid circles).

Figure 4 (
a) shows the equilibria plotted against c for fixed e (e = 0.2); note that here ρ = 0.25 and ρ + e = 0.45.