- Email: [email protected]

Metapopulation Models for Extinction Threshold in Spatially Correlated Landscapes Otso Ovaskainennw, Kazunori Satoz, Jordi Bascompte y8 and Ilkka Hanskin n

Metapopulation Research Group, Department of Ecology and Systematics, PO Box 17, Arkadiankatu 7, FIN-00014 University of Helsinki, Finland, zDepartment of Systems Engineering, Faculty of Engineering, Shizuoka University 5-1, Johoku 3-chome, Hamamatsu 432-8561 Japan and yNational Center for Ecological Analysis and Synthesis, University of California at Santa Barbara, 735 State St. Suite 300, Santa Barbara, CA 93101, U.S.A (Received on 29 January 2001, Accepted in revised form on 22 November 2001)

Simple analytical models assuming homogeneous space have been used to examine the eﬀects of habitat loss and fragmentation on metapopulation size. The models predict an extinction threshold, a critical amount of suitable habitat below which the metapopulation goes deterministically extinct. The consequences of non-random loss of habitat for species with localized dispersal have been studied mainly numerically. In this paper, we present two analytical approaches to the study of habitat loss and its metapopulation dynamic consequences incorporating spatial correlation in both metapopulation dynamics as well as in the pattern of habitat destruction. One approach is based on a measure called metapopulation capacity, given by the dominant eigenvalue of a ‘‘landscape’’ matrix, which encapsulates the eﬀects of landscape structure on population extinctions and colonizations. The other approach is based on pair approximation. These models allow us to examine analytically the eﬀects of spatial structure in habitat loss on the equilibrium metapopulation size and the threshold condition for persistence. In contrast to the pair approximation based approaches, the metapopulation capacity based approach allows us to consider species with long as well as short dispersal range and landscapes with spatial correlation at diﬀerent scales. The two methods make dissimilar assumptions, but the broad conclusions concerning the consequences of spatial correlation in the landscape structure are the same. Our results show that increasing correlation in the spatial arrangement of the remaining habitat increases patch occupancy, that this increase is more evident for species with short-range than long-range dispersal, and that to be most beneﬁcial for metapopulation size, the range of spatial correlation in landscape structure should be at least a few times greater than the dispersal range of the species. r 2002 Elsevier Science Ltd.

Introduction Habitat loss is widely considered to be the principal threat to the long-term survival of wAuthor to whom correspondence should be addressed. E-mail: [email protected]ﬁ Current address: Department of Zoology, University of Cambridge, Downing Street, Cambridge CB2 3EJ, U.K. ! Biologica ! * 8Current address: Estacion de Donana, CSIC. Apdo 1056, E-41080 Sevilla, Spain. 0022-5193/02/050095 + 14 $35.00/0

species both locally and worldwide (Barbault & Sastrapradja, 1995). As with other complex ecological processes, direct observations (e.g. Andre! n 1994, 1996) and experiments (e.g. Debinski & Holt, 2000) are badly needed but they are hardly suﬃcient to decipher all the long-term and large-scale biological consequences of habitat loss and fragmentation, for which purpose robust theory is a necessity. r 2002 Elsevier Science Ltd.

96

O. OVASKAINEN ET AL.

The current population ecological theory on habitat loss and fragmentation is still rudimentary, being largely based on Lande’s (1987) extension of Levins’s (1969) familiar model of classic metapopulation dynamics. The latter model assumes an inﬁnitely large network of identical and equally connected habitat patches, with the fraction of occupied patches being the only dynamic variable. Lande added the assumption that only a fraction h of the patches is suitable for occupancy, though migrating individuals still continue to arrive at all patches, including the unsuitable ones. Habitat loss is thus represented by 1 h in Lande’s model, and the rate of change in the fraction of occupied patches p (out of all patches) is given by dp=dt ¼ cpðh pÞ ep;

ð1Þ

where e and c are, respectively, the extinctionand colonization-rate parameters. Increasing the fraction of unsuitable patches reduces the rate of establishment of new populations. Eventually, when habitat availability is reduced below a threshold value, the metapopulation goes extinct even if a fraction of the patches is still suitable for occupancy. Modelling studies based on eqn (1) include May (1991), Nee & May (1992) and Tilman et al. (1994). An apparently useful prediction of this simple model is that, at equilibrium, the fraction of empty but suitable patches is constant, E n ¼ h pn ¼ e=c ¼ d; for all values of habitat loss allowing deterministic persistence. Therefore, the threshold level of habitat availability needed for long-term persistence (h ¼ d) can be estimated simply from the knowledge of the equilibrium fraction of empty habitat ( pn ) in landscapes with h > d (Lande, 1987; Nee, 1994; Lawton et al., 1994). Hanski et al. (1996) have coined the term Levins rule for this model prediction. Some attempts to test the Levins rule have been made by Doncaster et al. (1996), Andre! n (1997) and Hanski (1999). The above-described model represents a very simple caricature of habitat loss and its population ecological consequences. To start with, the Levins model lacks any structure in terms of local population sizes and patch qualities. Analysing such structured models, Gyllenberg

& Hanski (1997) showed that whether dE n =dh equals zero (the Levins rule) depends on the relative values of existing and added/removed patches to support metapopulation growth. In the Levins model, all patches are equal, hence the Levins rule applies. Including the rescue eﬀect into the Levins model (Hanski et al., 1996; Gyllenberg & Hanski, 1997) leads to dE n =dho0; because now local extinction rate decreases with increasing p and hence the value of a patch depends on the state of the entire metapopulation. Another structure missing from the Levins model is the spatial population structureFthe Levins model assumes that all patches and populations are equally connected. This is a reasonable assumption only for species with large dispersal range in relation to landscape size or for landscapes with little spatial heterogeneity in habitat patch characteristics (area, density, etc.). Furthermore, to apply the Levins rule implies that habitat loss occurs via random loss of entire patches, retaining spatial homogeneity in landscape structure. In reality, of course, landscapes tend to be heterogeneous and habitat loss occurs non-randomly, for instance being greater in some parts of the landscape than in some other parts, and habitat may be lost via reduction in the areas of the existing patches rather than via complete loss of patches. Assuming realistically that dispersal distances are limited, we might expect that it makes a diﬀerence in which manner habitat is actually destroyed (McCarthy et al., 1997). This expectation is conﬁrmed by simulation studies (Dytham, 1995; Bascompte & Sole! , 1996; Tilman et al., 1997; With & King, 1999; Hill & Caswell, 1999) and an analytical study comparing three diﬀerent landscape patterns (Neuhauser, 1998). Generally, increasingly localized dispersal makes colonization more diﬃcult and shifts the extinction threshold to lower values of habitat loss (Bascompte & Sole! , 1996), while spatially aggregated habitat loss may allow the species to persist in the untouched part of the landscape and shifts the extinction threshold to higher values of habitat loss (With & King, 1999; Hill & Caswell, 1999). In this paper, we present two diﬀerent analytical approaches to study the

METAPOPULATIONS IN CORRELATED LANDSCAPES

metapopulation dynamic consequences of spatially structured habitat loss in a patch network represented by a regular lattice of inﬁnite size. Our goal is to present an improved general understanding of the dependence of the extinction threshold on the spatial structure of the landscape and on the spatial pattern of habitat loss. The ﬁrst approach is based on the recent work by Ovaskainen and Hanski (Hanski & Ovaskainen, 2000; Ovaskainen & Hanski, 2001), who modelled the rates of change in the probabilities of patch occupancy in a ﬁnite patch network. This approach combines the idea of classic metapopulation dynamics (as in the Levins model) with assumptions about how the habitat patch areas and connectivities inﬂuence the rates of local extinction and colonization (as in the incidence function model; Hanski, 1994). It turns out that the landscape eﬀects on population turnover can be captured in a single number, the dominant eigenvalue of an appropriate ‘‘landscape’’ matrix (Adler & Nuernberger, 1994), which number we have termed the metapopulation capacity of a fragmented landscape (Hanski & Ovaskainen, 2000). With metapopulation capacity, we are able to demonstrate how the species dispersal range and the range of spatial correlation in landscape structure aﬀect metapopulation persistence. The metapopulation capacity based approach can also be applied to irregular networks of habitat patches (see Discussion and Hanski & Ovaskainen, 2000). Our second approach is based on pair approximation, a technique originally developed in the physical sciences to correct mean-ﬁeld models [such as eqn (1)] when the explicit spatial structure should not be ignored (Dickman, 1986; Katori & Konno, 1991). After an early version of the present form (Matsuda, 1987; Matsuda et al., 1987), pair approximation was introduced in ecology and epidemiology as a way of exploring spatial structure in the contact process and in similar spatially explicit models (Matsuda et al., 1992; Harada & Iwasa, 1994; Sato et al., 1994; Satulovsky & Tome! , 1994; Tainaka, 1994a, b; Harada et al., 1995). Subsequently, the applications of pair approximation and other moment closure techniques in ecology have

97

quickly increased (Sato & Konno, 1995; Kubo et al., 1996; Levin & Durret, 1996; Nakamaru et al., 1997, 1998; Takenaka et al., 1997; Iwasa et al., 1998; Levin & Pacala, 1998; Pacala & Levin, 1998; Ives et al., 1998; Bolker & Pacala, 1999; Dieckmann et al., 2000; Hiebeler, 2000). The basic idea is to include in the diﬀerential equations for the mean occupancy second-order terms deﬁning the densities of pairs of points (occupied–empty, occupied–occupied, and so on). The dynamics of pairs are described by additional diﬀerential equations, which contain yet higher-order terms, etc. The usual strategy is to close the equations by appropriately neglecting correlations beyond pairs. Pair approximation has recently been used to analyse spatial structure of habitat loss in metapopulation dynamics (Hiebeler, 2000), though for a discrete-time model in contrast with the continuous time model in the present paper. Ives et al. (1998) have previously modelled the correlation among suitable sites in their application of pair approximation to the study of how the pattern of land-cover aﬀects the dynamics of herbs in forests. The Metapopulation Capacity based Approach The metapopulation capacity lM is a measure that describes the capacity of a highly fragmented landscape to support a viable metapopulation (Hanski & Ovaskainen, 2000). As a ﬁrst approximation, lM measures the amount of suitable habitat, like h in Lande’s version of the Levins model [eqn (1)], but additionally lM takes into account the actual spatial conﬁguration of the fragmented landscape. Metapopulation capacity was initially developed for patch-based models (Ovaskainen & Hanski, 2001), but it can be applied equally well for lattice models as is done here. We model metapopulation dynamics with a spatially explicit version of the Levins model. The full stochastic model is a Markov process, but we will restrict the analysis to the deterministic version, deﬁned as the ﬁrst-order approximation of the master equations. In the deterministic model, the rate of change in pi ; the probability that a suitable site i is occupied, is given by (Hanski &

98

O. OVASKAINEN ET AL.

Gyllenberg, 1997)

behaves as

dpi ¼ Ci ðpÞð1 pi Þ Ei pi ; dt

ð2Þ

where p is a vector with components pi ; Ci ðpÞ is the colonization rate of site i (when empty) and Ei is the extinction rate of site i (when occupied). We will assume that the colonization rate Ci ðpÞ is given as Ci ðpÞ ¼ cSi ðpÞ; where c is a speciesspeciﬁc parameter and Si ðpÞ ¼ Sja i sij pj gives the metapopulation-dynamic connectivity of site i when the system is in state p (Ovaskainen & Hanski, 2001), the component sij giving the contribution that patch j makes (when occupied) to the connectivity of patch i: The exact functional form of Si ðpÞ depends on further model assumptions and will be speciﬁed below. As all sites are assumed to be identical, we will assume that the local extinction rate of site i is given as Ei ¼ e; where e is a species-speciﬁc parameter. Under these assumptions, the threshold condition for persistence is (Adler & Nuernberger, 1994; Ovaskainen & Hanski, 2001) lM > d;

ð3Þ

where the metapopulation capacity lM is given as the leading eigenvalue of an appropriate ‘‘landscape’’ matrix M describing the dependence of the colonization and extinction rates on the spatial conﬁguration of habitat patches, and d is the ratio between the species-speciﬁc extinction and colonization rate parameters, d ¼ e=c (Hanski & Ovaskainen, 2000). With the above model assumptions, the elements of matrix M are simply given as mij ¼ sij : Under the structural model assumptions made here, persistence in the sense of a stable non-trivial equilibrium state is equivalent with invasibility in the sense of instability of the trivial equilibrium state (Ovaskainen & Hanski, 2001). The condition set by eqn (3) follows easily by stability analysis. In Lande’s version of the Levins model [eqn (1)], the equilibrium fraction of occupied patches is given by pn ¼ 1 d=h: This equation has a counterpart in the spatially realistic Levins model, where a weighted average of the equilibrium values pni in eqn (2), denoted by pnl ;

pnl ¼ 1 d=lM :

ð4Þ

The type of the weighted average is based on the dynamical values of the individual sites, and may be derived as follows. Let x denote the right leading eigenvector of matrix M; and let r denote * deﬁned the left leading eigenvector of matrix M; n by m* ij ¼ mij ð1 pi Þ=d: Scaling the two vectors as rT x ¼ 1; and deﬁning the weights wi by wi ¼ ri xi ; it may be shown that pnl deﬁned by pnl ¼ Si wi pni satisﬁes eqn (4). To illustrate the nature of this weighted average, we ﬁrst note that x is proportional to the equilibrium state at the extinction threshold (Ovaskainen & Hanski, * is a density2001). Second, the matrix M dependent version of matrix M; r representing the ‘‘reproductive value’’ of patch i for a species described by the parameter d: Taking into account the exact relationship between pnl and lM given by eqn (4), we conclude that the weighting w combines information from the threshold condition and the equilibrium state in an appropriate way. As illustrated by Fig. 1, well-connected sites attain the highest weights, whereas the weights assigned to isolated patches can be negligible, especially if the species has a short-range dispersal kernel. In this paper, we consider a landscape consisting of inﬁnitely many local sites (‘‘habitat patches’’) arranged as a regular two-dimensional lattice. We will analyse how the metapopulation capacity of this inﬁnite-size lattice changes with habitat loss. We will assume that spatial correlation in the landscape structure decays exponentially with parameter b: The dispersal kernel of the focal species is also exponential with parameter a: In the analysis we need to know the distribution of suitable sites around a site that is known to be suitable, leading to conditional probabilities of sites being suitable. Let Hk denotes the proposition ‘‘Site k is suitable habitat’’. We will utilize the notation hk ¼ 1 and hk ¼ 0 corresponding to the cases Hk ¼ True and Hk ¼ False: It is neither a trivial nor a well-deﬁned task to construct a landscape with ‘‘long-range spatial correlation’’. One approach commonly used by physicists is the Fourier ﬁltering method (FFM)

METAPOPULATIONS IN CORRELATED LANDSCAPES

99

Fig. 1. An illustration of the weighting w utilized in eqn (4). Panel (a) depicts a 30 30 lattice with periodic boundary conditions, habitat being indicated by black sites. In panels (b) and (c), the sizes of the dots are proportional to the weights associated with the individual sites in the lattice. In panel (b), the dispersal of the species is restricted to nearest neighbours (1=a ¼ 0), whereas in panel (c), the dispersal kernel has a long tail (a ¼ 1). The species parameter d was set to d ¼ 0:3; and the lattice is characterized by the parameters h ¼ 0:4; r ¼ 0:1 and b ¼ 1:

and its extensions, which are based on ﬁltering the Fourier components of an uncorrelated sequence of random numbers with a suitable correlation ﬁlter (Makse et al., 1996). However, in the present case this approach is unnecessarily complicated, as the goal is mainly to demonstrate whether and how longrange correlation in the landscape structure matters. Furthermore, in the analysis we will need only the one- and two-point conditional probabilities, higher-point conditional

probabilities having no eﬀect on our approximation. In order to keep the analysis as simple as possible, we characterize the landscape directly through these two conditional probabilities as PðHi jHk Þ ¼ h þ r 2bð1dki Þ ;

ð5Þ

PðHi -Hj jHk Þ ¼ PðHi jHk ÞPðHj jHk ÞPðHi jHj Þ=h; ð6Þ

100

O. OVASKAINEN ET AL.

where dki denotes the distance between sites k and i; measured as dki ¼ jxk xi j þ jyk yi j: Note that the parameters h and r must be chosen so that probabilities (5) and (6) are always between 0 and 1. Equation (5) is based on the assumption that the unconditional probability for site i being suitable is h; 1 h measuring thus the amount of habitat loss. If site k at distance dki ¼ 1 from site i is known to be suitable, the conditional probability for site i being suitable is h þ r: When dki changes between 1 and N; the conditional probability for site i being suitable decays exponentially from h þ r to h: The amount of correlation is thus measured by r; whereas the range of correlation is measured by 1=b: If r > 0; there is a positive correlation in the landscape structure, whereas ro0 corresponds to a negative correlation. With real landscapes, the case r > 0 is often relevant, whereas the case ro0 is mainly of theoretical interest. Small 1=b corresponds to short-range correlation, and large 1=b corresponds to long-range correlation. Equation (6) is an extension of eqn (5), and it is consistent in the sense that landscape structure is invariant with respect to space. Some landscape patterns illustrating eqns (5) and (6) are shown in Fig. 2. The algorithm with which these patterns have been generated is described in Appendix A. We emphasize that while this is not the most general way to incorporate spatial correlation into landscape structure, this is an analytically tractable approach that is reasonable for the purpose of demonstrating the consequences of habitat loss in spatially correlated landscapes. It should also be noted that the metapopulation capacity based approach is not dependent on the exponential assumption in eqn (5), and it could be as readily applied to any other structural assumptions. We assume that the dispersal behaviour of the species is described by an exponentially decaying dispersal kernel. Dispersal is deﬁned in terms of Si ðpÞ; the metapopulation dynamic connectivity of site i; given that the system is in state p; Si ðpÞ ¼ CðaÞSjAAðiÞ 2að1dij Þ pj ;

ð7Þ

where AðiÞ is the set of all suitable sites other than site i: The dispersal kernel is scaled by

CðaÞ ¼ 22ð1þaÞ ð2a 1Þ2 to give Si ð1Þ ¼ 1 in the case of a fully occupied lattice with all sites being suitable habitat. Here, 1 denotes a vector with all components equal to 1. As the lattice and thus the ‘‘landscape’’ matrix are inﬁnitely sized, the leading eigenvalue lM is not analytically tractable. Before proceeding to analyse habitat loss with metapopulation capacity, we have to derive an approximation for lM : To do this, we deﬁne Ri ¼ Si ð1Þ; which corresponds to all patches being occupied in eqn (7). Ovaskainen & Hanski (2001) showed that Rl ¼ EðRÞ þ VarðRÞ=EðRÞ is a good approximation for the metapopulation capacity lM ; and we rely on this approximation in this paper. The approximation is based on the ﬁrst terms given by the ‘‘power method’’, which is an iterative scheme for computing the leading eigenvalue (Golub and van Loan, 1990). It is possible to evaluate analytically the inﬁnite sums needed to compute EðRÞ; leading to EðRÞ ¼ h þ

22b ð2a 1Þ2 r: ð2aþb 1Þ2

Unfortunately, the expression for VarðRÞ ¼ EðR2 Þ EðRÞ2 is too complicated to be solved analytically and we have to rely on numerical approximation. To do this, we deﬁne E k ðR2 Þ ¼ CðaÞ2 Si;jAQk 21d0i 21d0j Eðhi hj Þ; where Qk is the square fk; y; kg fk; y; kg with the origin (i ¼ 0) excluded, and Eðhi hj Þ ¼ PðHi jH0 Þ for i ¼ j and Eðhi hj Þ ¼ P ðHi -Hj jH0 Þ for iaj: A straightforward way to evaluate EðR2 Þ ¼ limk-N E k ðR2 Þ would be to compute E k ðR2 Þ for k ¼ 1; 2; y until convergence to desired accuracy had been reached. Unfortunately, especially with small a (or b), in which case the dispersal kernel (or correlation in the landscape structure) has a long range, convergence is achieved only when k is very large. Noting that the computation of E k ðR2 Þ in the two-dimensional lattice involves four successive sums from k to k; it is clear that the method is computationally not feasible. To speed up convergence, we assume that the sum converges as 2minða; bÞk ; and write E kþ1 ðR2 Þ ¼ E k ðR2 Þ þ cðkÞ2minða; bÞk ; where cðkÞ is assumed to

METAPOPULATIONS IN CORRELATED LANDSCAPES

101

Fig. 2. Spatial pattern of habitat loss. The fraction of suitable habitat (non-destroyed) is h ¼ 0:5 in all four panels. Habitat is indicated by black sites: (a) random pattern (r ¼ 0), (b),(c), positive correlation (r ¼ 0:2) and (d), negative correlation (r ¼ 0:1). Panels (b) and (d) correspond to landscapes used in pair approximation (1=b ¼ 1), whereas in panel (c) the range of spatial correlation is longer (1=b ¼ 2).

approach a constant c: Numerical experiments demonstrate that our assumption is correct, and with moderate k thus EðR2 Þ may be extrapolated P minða; bÞl 2 : as EðR2 ÞEE k ðR2 Þ þ c N l¼k The Pair Approximation based Approach In this approach we model metapopulation dynamics on two-dimensional inﬁnitely sized regular lattice using a basic contact process with colonization restricted to one of the four nearest

neighbour (NN) sites, corresponding to the case 1=a-0 in the metapopulation capacity based approach. Each lattice site is in one of the following states: destroyed (0), empty but suitable (1), or occupied (2). We denote the fractions of sites that are destroyed, empty but suitable and occupied by r0 ; r1 and r2 ; respectively, so that h ¼ 1 r0 ¼ r1 þ r2 represents the fraction of suitable habitat. These variables are called singlet densities. Let us further denote by rij the fraction of (ij) pairs of

102

O. OVASKAINEN ET AL.

sites, where i; jAf0; 1; 2g: These variables, called doublet densities, represent the probability that two randomly chosen nearest neighbour sites are in the conﬁguration ij: Finally, qi=j is the conditional probability of a site being in state i given that one of its NN sites is in state j: In other words, qi=j is a measure of NN correlation, which is equivalent to the mean crowding index (Lloyd, 1967; Iwao, 1968; Harada & Iwasa, 1994). By deﬁnition qi=j ¼ rij =rj : In a system with global mixing, rij ¼ ri rj ; and hence qi=j ¼ ri : Similar to eqn (5), we measure eﬀective correlation among nearby sites as r þ 2r12 þ r22 h: r ¼ 11 h

ð8Þ

Expression (8) represents the increase (or decrease) in correlation with respect to what we would expect in a random distribution of suitable and destroyed sites in a landscape (Ives et al., 1998; Bascompte, 2001). As shown in Appendix A, this does not exactly correspond to eqn (5) used with the metapopulation capacity based approach. In the comparison between the two approaches, we have approximated the landscape for the pair approximation by choosing b ¼ 1 in eqn (5). We normalize the extinction rate e by setting it equal to 1, in which case the normalized colonization rate is given by c ¼ e=d ¼ 1=d: The time derivatives for singlet densities are (Matsuda et al., 1992) dr2 ¼ r2 q1=2 =d r2 ; dt

ð9Þ

dr1 ¼ r2 r2 q1=2 =d; dt

ð10Þ

(permanently destroyed). For doublet densities we can write: dr20 1 ¼ 1 q2=10 r10 =d r20 ; ð12Þ z dt dr10 1 ¼ r20 1 q2=10 r10 =d; z dt

ð13Þ

dr00 ¼ 0; dt

ð14Þ

dr11 1 ¼ 2r12 2 1 q2=11 r11 =d; z dt

ð15Þ

dr22 1 1 ¼ 2 þ 1 q2=12 r12 =d 2r22 ; z z dt

ð16Þ

dr12 1 ¼ r22 þ 1 q2=11 r11 =d z dt 1 1 þ 1 q2=12 r12 =d r12 ; z z

ð17Þ

where z is the number of nearest neighbours (z ¼ 4 in the present case). Note that the equations for doublets contain higher-order terms. Pair approximation or doublet decoupling approximation consists of neglecting such higher-order terms by assuming that qi=jk Eqi=j : Once the previous system of equations is closed, we have only three independent variables and we can solve them at the steady state. An analytical result for the equilibrium patch occupancy (out of suitable patches) pn ¼ rn2 =h can be explicitly obtained as follows (with z ¼ 4):

1 p ¼ 15 12ð1 hÞð1 q0=0 Þ=h 8d 6 2d qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3 16ð1 hÞð1 q0=0 Þd=h þ ½3 4ð1 hÞð1 q0=0 Þ=h 2 : n

dr0 ¼ 0: dt

ð11Þ

Equation (11) is based on the assumption that the fraction of destroyed sites is ﬁxed

ð18Þ

This analytical result has previously been shown to be in good correspondence with spatially explicit simulations [Harada et al., 1995, see Fig. 3(a)]. The mean-ﬁeld version of

METAPOPULATIONS IN CORRELATED LANDSCAPES

this model, corresponding to the Levins model, can be written as dr2 ¼ r2 r1 =d r2 : dt

ð19Þ

The equilibrium fraction of occupied patches (out of suitable patches) is pn ¼ rn2 =h ¼ 1 d=h; consistent with eqn (1). Results Figure 3 illustrates the eﬀect of changing landscape structure on the fraction of occupied patches out of suitable ones, deﬁned by pn ¼ rn2 =h in the pair approximation based approach, and by pn ¼ 1 d=lM in the metapopulation capacity based approach. Panels (a) and (b)

103

show the eﬀect of habitat loss and fragmentation for a species with dispersal restricted to nearest neighbour sites (1=a-0). Here, we also assume that the range of spatial correlation in the landscape structure corresponds to pair approximation, which is approximated by 1=b ¼ 1 in the metapopulation capacity based approach. Panel (a) depicts the result for pair approximation, whereas panel (b) shows the result for metapopulation capacity based approach. From these results we may draw the following conclusions. First, with increasing habitat loss the size of the metapopulation at equilibrium declines. Second, with a ﬁxed amount of habitat loss, increasing correlation in the spatial arrangement of the remaining habitat increases metapopulation size.

Fig. 3. The eﬀect of landscape structure on the fraction of occupied patches. Panels (a) and (b) correspond to the assumptions of the pair approximation both with respect to the dispersal kernel of the species (1=a ¼ 0) and with respect to the spatial correlation in the landscape structure (approximated by 1=b ¼ 1 in the metapopulation capacity based approach). Panels (b) (1=a ¼ 1=b ¼ 1) and (c) (1=a ¼ 1; h ¼ 0:6) depict the metapopulation capacity result for a species with long-range dispersal kernel. The lines in panel (a) are the predictions of the pair approximation (pn ¼ rn2 =h), whereas the lines in panels (b)–(d) are based on the metapopulation capacity based approach (pnl E1 d=Rl ). In all four panels, the lines correspond (from top to down) to the parameter values r ¼ 0:2; r ¼ 0:1; r ¼ 0 and r ¼ 0:1: In order to assess the validity of the approximations, a set of 30 30 lattices with periodic boundary conditions was generated (one such lattice is illustrated in Fig. 1). The circles depict deterministic equilibrium states and the crosses depict stochastic simulation results for this set of lattices. The initial state for the simulation was chosen to correspond to the deterministic equilibrium state, and the simulation was run (after an initial phase of 500 events) until 5000 events (extinctions or colonizations) had occurred. The value 0 in the ﬁgure indicates that the simulation went extinct. Both the circles and the crosses are depicted for the values r ¼ 0 and 0.1. In all four panels, the value of the species parameter d was set to d ¼ 0:3:

104

O. OVASKAINEN ET AL.

Pair approximation predicts a smaller metapopulation size than the metapopulation capacity based approach [Fig. 3(a) and (b)]. This is largely due to the fact that in pair approximation, the result is shown for the fraction of occupied patches ( pn ), whereas for metapopulation capacity, the result is shown for a weighted fraction of occupied patches ( pnl ). As the weights are largest for the best connected sites, it is evident that the weighted average of the occupancy probabilities is higher than just the simple fraction of occupied patches. However, this does not explain the diﬀerence between the two models completely. An additional factor is that the metapopulation capacity is nonetheless a mean-ﬁeld approximation in the sense that it ignores the part of spatial correlation in patch occupancy that is generated by population dynamics, though it takes explicitly into account the spatial structure of the landscape. In contrast, pair approximation includes the eﬀect of population dynamics on spatial correlation. To see that this factor lowers the predicted fraction of occupied patches for the pair approximation, we note that the aggregation of occupied sites lowers the probability that an empty but suitable site has occupied sites as its neighbours, thus reducing the frequency of successful colonization events. Nonetheless, the qualitative predictions made by the two models are similar [Fig. 3(a) and (b)], suggesting that the results are not likely to be heavily modeldependent. Both approximations are based on a deterministic model, which assumes that the lattice is inﬁnitely large. In contrast, real landscapes are ﬁnite, and both the details of the exact lattice conﬁguration and stochastic population ﬂuctuations are likely to aﬀect metapopulation size. This is illustrated by the circles and crosses in Fig. 3, which shows the deterministic equilibrium state of eqn (2) and results from stochastic simulations. As the relation pnl ¼ 1 d=lM is exact, the discrepancy between the lines and the circles in Fig. 3(b) –(d) represents solely the error in the approximation Rl ElM for the sample landscapes. The discrepancy between the circles and the crosses represents the diﬀerence between the deterministic mean-ﬁeld model and the full stochastic model. The stochastic model predicts

a lower level of habitat occupancy due to the fact that the mean-ﬁeld model ignores the part of spatial correlation generated by population dynamics (see above). Note that we have to restrict the range of habitat loss (1 h) considered in some of the cases depicted in Fig. 3. For example, the combination 1 ho0:2 and r ¼ 0:2 is not feasible, as eqn (5) would give a conditional probability greater than 1 for ﬁnding a suitable site next to a known suitable site. The seemingly strange behaviour in Fig. 3(b) of the positively correlated cases with nearly all habitat lost is explained by the inﬁnite size of the lattice. If the size of the lattice would be ﬁnite, the case h-0 together with r ﬁxed to a positive value would not be possible, as h-0 indicates that almost all habitat is lost and thus there cannot be an area with a cluster of habitat patches. However, if the size of the lattice is allowed to approach inﬁnity, one may conceive such a situation. Consider a relatively large cluster of habitats arranged as a square. Let this cluster be located in the middle of a lattice, and let all the other sites in the lattice be destroyed. Now letting the size of the lattice approach inﬁnity we have h-0 with r ﬁxed to a positive value. The counterintuitively good performance of the species at the limit h-0 [Fig. 3(a)] is explained by this example. There are no important ecological conclusions to be drawn from this example, but it demonstrates that results derived from an inﬁnitely sized lattice must be interpreted cautiously, especially at the limits of the landscape parameters. Panels (c) and (d) illustrate the eﬀect of habitat loss and fragmentation for a species with long-range dispersal kernel. Below, we call this species as species B, and the species with shortrange dispersal [panel (b)] as species A. As pair approximation considers only nearest neighbours, it cannot be used here, and the results are given only for the metapopulation capacity based approach with 1=a ¼ 1: In panel (c), we keep the correlation in the landscape structure at the same level as in panel (b). There are two conclusions to be drawn. First, uncorrelated habitat loss (r ¼ 0) is even more detrimental for species B than it was for species A. Second, increasing correlation in landscape structure (r > 0) helps species B less than it helps species

METAPOPULATIONS IN CORRELATED LANDSCAPES

A. These results may seem counterintuitive, as it should be easier for the more dispersive species B to locate suitable habitat in increasingly fragmented landscapes. However, one should keep in mind that in panel (b) the scale of habitat fragmentation is rather short [as shown in Fig. 2(b)]. Thus, in the fragmented landscape, the largest remaining habitat blocks as perceived by species A are larger than those perceived by species B. As both species do equally well in an undisturbed environment (h ¼ 1), it is not surprising that species B is more sensitive to the type of habitat loss considered here. Panel (d) is again for species B, but instead of varying the amount of suitable habitat as in panel (c), we now vary the range of correlation in the landscape structure. These results further elucidate the reason for the relatively poor performance of species B in panel (c). While keeping both the amount of suitable habitat (h) and the level of correlation (r > 0) ﬁxed, increasing the range of landscape correlation makes a substantial contribution to the persistence of species B. As a rule of thumb, we conclude that it is beneﬁcial for the species if the range of spatial correlation in landscape structure is at least a few times greater than the dispersal range of the species. Discussion We have used mathematical models to study the eﬀect of spatially correlated habitat loss on the equilibrium metapopulation size, which has previously been investigated with numerical simulations only, with the exception of Hiebeler (2000). In this paper, we used two approaches, one based on the metapopulation capacity measure and the other one based on pair approximation. Pair approximation provides one simple way of introducing spatial correlation into analytical models and improving predictions based on mean-ﬁeld models that assume homogeneous space. As with all approximations, a number of assumptions are made. In particular, we neglect higher-order correlations than those among pairs. This is reasonable in the context of basic contact process models, in which interactions are restricted to the nearest neighbour sites only.

105

However, there are situations where this simpliﬁcation may break down, for instance near the extinction threshold, which is the equivalent in physics of a critical point, where a phase transition takes place. In general, deviations between pair approximations and simulation results are greatest near the extinction threshold (Matsuda et al., 1992; Harada & Iwasa, 1994; Sato et al., 1994; Harada et al., 1995; Levin & Durret, 1996). Here, we have used an ‘‘ordinary’’ pair approximation to investigate the qualitative eﬀect of spatial correlation on the population dynamic consequences of habitat destruction. We note that the quantitative approximations could be improved e.g. by taking into account the higher-order terms of the spatial conﬁguration (Dickman, 1986; Sato et al., 1994; van Baalen, 2000). Other approaches may avoid this problem by estimating r from spatially explicit simulation and substituting the value thus obtained into a mathematical expression (Bascompte, 2001). However, the present model is tractable analytically, and its level of accuracy is suﬃcient to clearly show the eﬀect of spatial correlation in the pattern of habitat loss on equilibrium metapopulation size. The same type of conclusions that we drew from results based on pair approximation [Fig. 3(a)] were drawn by Hiebeler (2000) for a discrete-time model, with one exception. He concluded that ‘‘pn depends mainly on the habitat-clustering parameter q00 ; and only very slightly on the amount of habitat available, p0 ’’ (p0 ¼ h and 1 q00 ¼ ð1 hÞq0=0 =h in our notation), which seems to be in contradiction with Fig. 3(a), where the fraction of occupied patches strongly depends on the amount of habitat loss. However, the discrepancy is not in the models but in the interpretation of the results. Hiebeler (2000) studied the eﬀect of habitat loss and fragmentation with respect to parameters h and q00 ; whereas we have done it in terms of h and r: Hiebeler’s (2000) counterintuitive conclusion that the fraction of occupied patches would depend only very slightly on the amount of habitat available is explained by the high correlation between h and q00 ; a problem which was here solved by deﬁning and using the eﬀective correlation parameter r:

106

O. OVASKAINEN ET AL.

The metapopulation capacity based approach relies on a diﬀerent kind of approximation than pair approximation. The great advantage of the metapopulation capacity is that it can readily be applied to landscapes with arbitrarily long (or short) correlation range, and it allows one to make structural assumptions other than the exponential decay used in eqn (5). The metapopulation capacity was originally developed for irregular networks of habitat patches, in which both the areas and connectivities of habitat patches vary from patch to patch. Many real metapopulations live in such networks (Hanski & Gilpin, 1997; Hanski, 1999). Hanski & Ovaskainen (2000) have studied the eﬀect of habitat loss in the context of a spatially realistic metapopulation model and an irregular patch network. In agreement with our results, Hanski & Ovaskainen (2000) found that if habitat is lost in large blocks rather than in a random pattern, the metapopulation capacity lM of the patch network decreases less than in proportion to the area lost, as metapopulation dynamics in the remaining habitat are aﬀected relatively little. In the present case, persistence is enhanced by positive correlation in landscape structure, which implies that the remaining habitat occurs in large blocks. Hanski & Ovaskainen (2000) also showed that a highly detrimental form of habitat loss in the model that they studied consists of each patch losing a constant fraction of its area, in agreement with our results for the case of negatively correlated (ro0) landscape pattern. The population dynamic signiﬁcance of the spatial scale of habitat loss is inﬂuenced by the spatial scale of dispersal of the focal species, as illustrated in Fig. 3(a) and (b). By assuming long-range dispersal, we eﬀectively make the landscape more homogeneous for the species and reduce the signiﬁcance of short-range landscape correlation. The opposite happens when the dispersal range of the species is short but the spatial correlation in landscape structure is long. The overall conclusion from our analysis is that the spatial pattern in habitat loss can have a major eﬀect on the persistence of a species in a fragmented landscape. As habitat loss in real landscapes tends to be non-random, this eﬀect should not be ignored.

OO and IH thank the Academy of Finland for funding (Grant 50165 and the Finnish Centre of Excellence Programme (2000–2005), Grant 44887). KS and JB are grateful to the Spatial Ecology Program of the Division of Population Biology, University of Helsinki, for ﬁnancial support and intellectual stimulus during a short visit. JB was a Postdoctoral Associate at the National Center for Ecological Analysis and Synthesis, funded by NSF (Grant DEB-94-21535), the University of California at Santa Barbara, and the State of California.

REFERENCES Adler, F. R. & Nuernberger, B. (1994). Persistence in patchy irregular landscapes. Theor. Popul. Biol. 45, 41–75. Andre! n, H. (1994). Eﬀects of habitat fragmentation on birds and mammals in landscapes with diﬀerent proportions of suitable habitat: a review. Oikos 71, 355–366. Andr e! n, H. (1996). Population responses to habitat fragmentation: statistical power and the random sample hypothesis. Oikos 76, 235–242. Andr e! n, H. (1997). Habitat fragmentation and changes in biodiversity. Ecol. Bull. 46, 171–181. van Baalen, M. (2000). Pair approximation for diﬀerent spatial geometries. In: Geometry of Ecological Interactions. Simplifying Spatial Complexities (Dieckmann, U., Law, R. & Metz, J. A. J., eds.), pp. 359–387. Cambridge: Cambridge University Press. Barbault, R. & Sastrapradja, S. D. (1995). Generation, maintenance and loss of biodiversity. In: Global Biodiversity Assesment (Heywood, W. H., ed), pp. 193–274. Cambridge: Cambridge University Press. Bascompte, J. (2001). Aggregate statistical measures and metapopulation dynamics. J. theor. Biol. 209, 373–379, doi:10.1006/jtbi.2001.2275. Bascompte, J. & Sol e! , R. V. (1996). Habitat fragmentation and extinction thresholds in spatially explicit models. J. Anim. Ecol. 65, 465–473. Bolker, B. M. & Pacala, S. W. (1999). Spatial moment equations for plant competition: understanding spatial strategies and the advantages of short dispersal. Am. Nat. 153, 575–602. Debinski, D. M. & Holt, R. D. (2000). A survey and overview of habitat fragmentation experiments. Conserv. Biol. 14, 342–355. Dickman, R. (1986). Kinetic phase transitions in a surfacereaction model: mean-ﬁeld theory. Phys. Rev. A 34, 4246–4250. Dieckmann, U., Law, R. & Metz, J. A. J. (2000). The Geometry of Ecological Interaction: Simplifying Spatial Complexity. Cambridge: Cambridge University Press. Doncaster, C. P., Micol, T. & Jensen, S. P. (1996). Determining minimum habitat requirements in theory and practice. Oikos 75, 335–339. Dytham, C. (1995). The eﬀect of habitat destruction pattern on species persistence: a cellular model. Oikos 74, 340–344. Golub G. N. & van Loan, C. E. 1990. Matrix Computations. Baltimore and London: The Johns Hopkins University Press.

METAPOPULATIONS IN CORRELATED LANDSCAPES

Gyllenberg, M. & Hanski, I. (1997). Habitat deterioration, habitat destruction and metapopulation persistence in a heterogeneous landscape. Theor. Popul. Biol. 52, 198–215, doi:10.1006/tpbi.1997.1333. Hanski, I. (1994). A practical model of metapopulation dynamics. J. Anim. Ecol. 63, 151–162. Hanski, I. (1999). Metapopulation Ecology. New York: Oxford University Press. Hanski, I. & Gilpin, M. E. (1997). Metapopulation Biology: Ecology, Genetics, and Evolution. San Diego: Academic Press. Hanski, I. & Gyllenberg, M. (1997). Uniting two general patterns in the distribution of species. Science 275, 397–400. Hanski, I. & Ovaskainen, O. (2000). The metapopulation capacity of a fragmented landscape. Nature 404, 755–758. Hanski, I., Moilanen, A. & Gyllenberg, M. (1996). Minimum viable population size. Am. Nat. 147, 527–541. Harada, Y. & Iwasa, Y. (1994). Lattice population dynamics for plants with dispersing seeds and vegetative propagation. Res. Popul. Ecol. 36, 237–249. Harada, Y., Ezoe, H., Iwasa, Y., Matsuda, H. & Sato, K. (1995). Population persistence and spatially limited social interaction. Theor. Popul. Biol. 48, 65–91, doi:10.1006/ tpbi.1995.1022. Hiebeler, D. (2000), Populations on fragmented landscapes with spatially structured heterogeneities: landscape generation and local dispersal. Ecology 81, 1629–1641. Hill, M. F. & Caswell, H. (1999). Habitat fragmentation and extinction thresholds on fractal landscapes. Ecol. Lett. 2, 21–127. Ives, A. R., Turner, M. G. & Pearson, S. M. (1998). Local explanations of landscape patterns: can analytical approaches approximate simulation models of spatial processes? Ecosystems 1, 35–51. Iwao, S. (1968). A new regression method for analyzing the aggregation pattern of animal populations. Res. Popul. Ecol. 10, 1–20. Iwasa, Y., Nakamaru, M. & Levin, S. A. (1998). Allelopathy of bacteria in a lattice population: competition between colicin-sensitive and colicin-producing strains. Evol. Ecol. 12, 785–802. Katori, M. & Konno, N. (1991). Upper bounds for survival probability of the contact process. J. Stat. Phys. 63, 115–130. Kubo, T., Iwasa, Y. & Furumoto, N. (1996). Forest spatial dynamics with gap expansion: total gap area and gap size distribution. J. theor. Biol. 180, 229–246, doi:10.1006/jtbi.1996.0099. Lande, R. (1987). Extinction thresholds in demographic models of territorial populations. Am. Nat. 130, 624–635. Lawton, J. H., Nee, S., Letcher, A. J. & Harvey, P. H. (1994). Animal distributions: patterns and process. In: Large-scale Ecology and Conservation Biology (Edwards, P. J., May, R. M. & Webb, N. R., eds), pp. 41–58. Oxford: Blackwell Scientiﬁc Press. Levin, S. A. & Durrett, R. (1996). From individuals to epidemics. Philos. Trans. R. Soc. London B 351, 1615–1621. Levin, S. A. & Pacala, S. W. (1998). Theories of simpliﬁcation and scaling of spatially distributed processes. In: Spatial Ecology. The Role of Space in Population Dynamics and Interspeciﬁc Interactions

107

(Tilman, D. & Kareiva, P., eds), pp. 271–295. Princeton Monographs in Population Biology, Princeton. Levins, R. (1969). Some demographic and genetic consequences of environmental heterogeneity for biological control. Bull. Entomol. Soc. Am. 15, 237–240. Lloyd, M. (1967). Mean crowding. J. Anim. Ecol. 36, 1–30. Makse, H. A., Havlin, S., Schwartz, M. & Stanley, E. (1996). Method for generating long-range correlations for large systems. Phys. Rev. E 53, 5445–5449. Matsuda, H. (1987). Condition for the evolution of altruism. In: Animal Societies: Theories and Facts. (Ito, Y., Brown, J. & Kikkawa, J., eds), pp. 67–80. Tokyo: Japan Scientiﬁc Society Press. Matsuda, H., Tamachi, N., Ogita, N. & Sasaki, A. (1987). A lattice model for population biology. In: Mathematical Topics in Biology (Teramoto, E. & Yamaguti, M., eds), pp. 154–161, Lecture Notes in Biomathematics, Vol. 71. New-York: Springer-Verlag. Matsuda, H., Ogita, N., Sasaki, A. & Sato, K. (1992). Statistical mechanics of populationFthe lattice Lotka–Volterra model. Prog. Theor. Phys. 88, 1035–1049. May, R. M. (1991). The role of ecological theory in planning reintroduction of endangered species. Symp. Zool. Soc. London 62, 145–163. McCarthy, M. A., Lindenmayer, D. B. & Dreschsler, M. (1997). Extinction debts and risks faced by abundant species. Conserv. Biol. 11, 221–226. Nakamaru, M., Matsuda, H. & Iwasa, Y. (1997). The evolution of cooperation in a lattice-structured population. J. theor. Biol. 184, 65–81, doi:10.1006/ jtbi.1996.0243. Nakamaru, M., Nogami, H. & Iwasa, Y. (1998). Scoredependent fertility model for the evolution of cooperation in a lattice. J. theor. Biol. 194, 101–124, doi:10.1006/ jtbi.1998.0750. Nee, S. (1994). How populations persist. Nature 367, 123–124. Nee, S. & May, R. M. (1992). Dynamics of metapopulations: habitat destruction and competitive coexistence. J. Anim. Ecol. 61, 37–40 Neuhauser, C. (1998). Habitat destruction and competitive coexistence in spatially explicit models with local interactions. J. theor. Biol. 193, 445–463, doi:10.1006/ jtbi.1998.0713. Ovaskainen, O. & Hanski, I. (2001). Spatially structured metapopulation models: metapopulation capacity and threshold conditions for persistence. Theor. Popul. Biol., in press. Pacala, S. W. & Levin, S. A. (1998). Biologically generated spatial pattern and the coexistence of competing species. In: Spatial Ecology. The Role of Space in Population Dynamics and Interspeciﬁc Interactions (Tilman, D. & Kareiva, P., eds), pp. 204–232. Princeton Monographs in Population Biology, Princeton. Sato, K. & Konno, N. (1995). Successional dynamical models on the 2-dimensional lattice space. J. Phys. Soc. Jpn. 64, 1866–1869. Sato, K., Matsuda, H. & Sasaki, A. (1994). Pathogen invasion and host extinction in lattice structured populations. J. Math. Biol. 32, 251–268. Satulovsky J. & Tome¤ , T. (1994). Stochastic lattice gas model for a predator–prey system. Phys. Rev. E 49, 5073–5079.

108

O. OVASKAINEN ET AL.

Tainaka, K. (1994a). Vortices and strings in a model ecosystem. Phys. Rev. E 50, 3401–3409. Tainaka, K. (1994b). Intrinsic uncertainty in ecological catastrophe. J. theor. Biol. 166, 91–99, doi:10.1006/ jtbi.1994.1007. Takenaka, Y., Matsuda, H. & Iwasa, Y. (1997). Competition and evolutionary stability of plants in a spatially structured habitats. Res. Popul. Ecol. 39, 67–75. Tilman, D., May, R. M., Lehman, C. L. & Nowak, M. A. (1994). Habitat destruction and the extinction debt. Nature 371, 65–66. Tilman, D., Lehman, C. L. & Yin, C. (1997). Habitat destruction, dispersal, and deterministic extinction in competitive communities. Am. Nat. 149, 407–435. With, K. A. & King, A. W. (1999). Extinction thresholds for species in fractal landscapes. Conserv. Biol. 13, 314–326.

Appendix A Generation of Correlated Landscapes

Here, we describe the algorithm used for generating the landscapes utilized in Figs. 1–3. The ﬁrst step is to assign suitable habitat to randomly located sites on an initially empty lattice until fraction h has been used. Next the conditional probabilities PðdÞ for a patch being suitable, conditioned on a patch at distance dðd ¼ 1; y; dmax ) being suitable, are computed. To do this, one has to average across all sites on the lattice, taking into account all sites at distance d from the focal site. The values PðdÞ; d ¼ 1; y; dmax are then compared with the conditional probabilities given by eqn (5), which indicates how much and in which direction landscape correlation at each spatial scale should be adjusted. The adjustment is done by randomly choosing two sites, one of which is

suitable and the other one destroyed, and changing their status if that would make the PðdÞ values, averaged across the spatial scales d; more similar to the target values given by eqn (5). The procedure is continued until the PðdÞ values are correct with desired accuracy. In principle, one should set dmax ¼ N to exactly mimic eqn (5). As there are 4d sites at distance d; the computational work needed is proportional to n2 dmax ; and dmax has thus in practice to be set to a relatively low value. In the computation used to generate Figs 1–3, we used dmax ¼ 3: To illustrate that the correlation structure assumed with pair approximation does not exactly correspond to the one described by eqn. (5), consider a one-dimensional lattice of inﬁnite size. Denoting by PðajbÞ the probability that a patch at location a is suitable, conditioned on that a patch at location b is known to be suitable, pair approximation assumes that Pð1j0Þ ¼ h þ r: There is no direct formula for Pð2j0Þ; as pair approximation is closed by neglecting higher-order correlations. However, we may evaluate Pð2j0Þ as Pð2j0Þ ¼ Pð2j1Þ Pð1j0Þ þ Pjð2j/1Þ Pjð/1j0Þ; where a/ means that the site at location a is destroyed. By elementary probability theory, Pð/1j0Þ ¼ 1 Pð1j0Þ and Pð2j/1Þ ¼ Pð/1j2Þ Pð2Þ=ð1 Pð1ÞÞ ¼ ð1 ðh þ rÞÞh= ð1 hÞ; where PðaÞ is the unconditional probability for location a being habitat. By induction, we ﬁnd that Pðkj0Þ ¼ h þ ð1 hÞ ðr=ð1 hÞÞk þ Oðrkþ1 Þ: Thus, the correlation structure is slightly diﬀerent from that speciﬁed by eqn (5).