Metapopulation Models for Extinction Threshold in Spatially Correlated Landscapes

Metapopulation Models for Extinction Threshold in Spatially Correlated Landscapes

J. theor. Biol. (2002) 215, 95–108 doi:10.1006/jtbi.2001.2502, available online at on Metapopulation Models for Extinction...

358KB Sizes 0 Downloads 8 Views

J. theor. Biol. (2002) 215, 95–108 doi:10.1006/jtbi.2001.2502, available online at on

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 effects 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 effects of landscape structure on population extinctions and colonizations. The other approach is based on pair approximation. These models allow us to examine analytically the effects 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 different 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 beneficial 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]fi 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 sufficient 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.



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 infinitely 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;


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 effect 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 difference in which manner habitat is actually destroyed (McCarthy et al., 1997). This expectation is confirmed by simulation studies (Dytham, 1995; Bascompte & Sole! , 1996; Tilman et al., 1997; With & King, 1999; Hill & Caswell, 1999) and an analytical study comparing three different landscape patterns (Neuhauser, 1998). Generally, increasingly localized dispersal makes colonization more difficult 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 different analytical approaches to study the


metapopulation dynamic consequences of spatially structured habitat loss in a patch network represented by a regular lattice of infinite 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 first 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 finite 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 influence the rates of local extinction and colonization (as in the incidence function model; Hanski, 1994). It turns out that the landscape effects 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 affect 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-field 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


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 differential equations for the mean occupancy second-order terms defining the densities of pairs of points (occupied–empty, occupied–occupied, and so on). The dynamics of pairs are described by additional differential 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 affects 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 first 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 configuration 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, defined as the first-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 &



Gyllenberg, 1997)

behaves as

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


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 speciesspecific 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 specified 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-specific parameter. Under these assumptions, the threshold condition for persistence is (Adler & Nuernberger, 1994; Ovaskainen & Hanski, 2001) lM > d;


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 configuration of habitat patches, and d is the ratio between the species-specific 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 :


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 * defined the left leading eigenvector of matrix M; n by m* ij ¼ mij ð1  pi Þ=d: Scaling the two vectors as rT x ¼ 1; and defining the weights wi by wi ¼ ri xi ; it may be shown that pnl defined by pnl ¼ Si wi pni satisfies eqn (4). To illustrate the nature of this weighted average, we first 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 infinitely many local sites (‘‘habitat patches’’) arranged as a regular two-dimensional lattice. We will analyse how the metapopulation capacity of this infinite-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-defined task to construct a landscape with ‘‘long-range spatial correlation’’. One approach commonly used by physicists is the Fourier filtering method (FFM)



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 filtering the Fourier components of an uncorrelated sequence of random numbers with a suitable correlation filter (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 effect 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 Þ ;


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



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 defined 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 ;


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 infinitely 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 define 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 first 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 infinite 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 define 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



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 infinitely 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



sites, where i; jAf0; 1; 2g: These variables, called doublet densities, represent the probability that two randomly chosen nearest neighbour sites are in the configuration 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 definition qi=j ¼ rij =rj : In a system with global mixing, rij ¼ ri rj ; and hence qi=j ¼ ri : Similar to eqn (5), we measure effective correlation among nearby sites as r þ 2r12 þ r22  h: r ¼ 11 h


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


dr1 ¼ r2  r2 q1=2 =d; dt


(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


dr00 ¼ 0; dt


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


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


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


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 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  3 16ð1  hÞð1  q0=0 Þd=h þ ½3  4ð1  hÞð1  q0=0 Þ=h 2 : n

dr0 ¼ 0: dt


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


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-field version of


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


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 effect of changing landscape structure on the fraction of occupied patches out of suitable ones, defined 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)


show the effect 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 fixed amount of habitat loss, increasing correlation in the spatial arrangement of the remaining habitat increases metapopulation size.

Fig. 3. The effect 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 figure 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:



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 difference between the two models completely. An additional factor is that the metapopulation capacity is nonetheless a mean-field 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 effect 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 infinitely large. In contrast, real landscapes are finite, and both the details of the exact lattice configuration and stochastic population fluctuations are likely to affect 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 difference between the deterministic mean-field model and the full stochastic model. The stochastic model predicts

a lower level of habitat occupancy due to the fact that the mean-field 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 finding 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 infinite size of the lattice. If the size of the lattice would be finite, the case h-0 together with r fixed 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 infinity, 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 infinity we have h-0 with r fixed 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 infinitely sized lattice must be interpreted cautiously, especially at the limits of the landscape parameters. Panels (c) and (d) illustrate the effect 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


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) fixed, 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 beneficial 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 effect 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-field 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.


However, there are situations where this simplification 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 effect 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 configuration (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 sufficient to clearly show the effect 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 effect 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 defining and using the effective correlation parameter r:



The metapopulation capacity based approach relies on a different 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 effect 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 affected 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 significance of the spatial scale of habitat loss is influenced by the spatial scale of dispersal of the focal species, as illustrated in Fig. 3(a) and (b). By assuming long-range dispersal, we effectively make the landscape more homogeneous for the species and reduce the significance 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 effect on the persistence of a species in a fragmented landscape. As habitat loss in real landscapes tends to be non-random, this effect 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 financial 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). Effects of habitat fragmentation on birds and mammals in landscapes with different 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 different 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-field 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 effect 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.


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 Scientific 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 simplification and scaling of spatially distributed processes. In: Spatial Ecology. The Role of Space in Population Dynamics and Interspecific Interactions


(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 Scientific 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 Interspecific 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.



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 first 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 infinite 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 find that Pðkj0Þ ¼ h þ ð1  hÞ ðr=ð1  hÞÞk þ Oðrkþ1 Þ: Thus, the correlation structure is slightly different from that specified by eqn (5).