Metapopulations in multifractal landscapes: on the role of spatial aggregation

Javier G.P Gamarra


The use of fractals in ecology is currently pervasive over many areas. However, very few studies have linked fractal properties of landscapes to generating ecological mechanisms and dynamics. In this study I show that lacunarity (a measure of the landscape texture) is a well suited ecologically scaled landscape index that can be explicitly incorporated in metapopulation models such as the classical Levins equation. I show that the average lacunarity of an aggregated landscape is linearly correlated to the habitat that a species with local spatial processed information may perceive. Lacunarity is a computationally feasible index to measure, and is related to the metapopulation capacity of landscapes. A general approach to multifractal landscapes has been conceived, and some analytical results for self-similar landscapes are outlined, including the specific effect of landscape heterogeneity, decoupled from that of contagion by dispersal. Spatially explicit simulations show agreement with the semi-implicit method presented.


1. Introduction

The quest for the relationship between patterns and processes that for so long has puzzled ecologists has a cornerstone in the field of spatial ecology. In particular, it is the aim of many landscape ecologists to discern the dynamic functionality of landscapes, but little theory has been produced to date (With & King 1999a; McIntyre & Wiens 2000; Vos et al. 2001; Halley et al. 2004). Life history and dispersal attributes of species have been formerly coupled with neutral landscape models to infer plausible mechanisms building up realistic patterns of species distributions (With & King 1999b). It is the interplay of space and ecological traits what seems to affect growth and extinction of populations (Hanski 1998). Since the seminal work by Levins (1969) much attention has been given to the value of space for the persistence of populations. In general, the different approaches that focus on space as a mechanism and not just a pattern are implicit (Levins 1969; Lande 1987; Hess 1996), lattice-based (Bascompte & Solé 1996; Hill & Caswell 1999; With & King 1999a), analytical (Ives et al. 1998; Hiebeler 2000; Bascompte 2001) and patch-based (Hanski & Ovaskainen 2000), the last ones coining the discipline of metapopulation theory. Since the 1990s most of the studies have agreed on the importance of the quantity of habitat and its spatial arrangement on population persistence (Cantrell & Cosner 1991; Ritchie 1997).

Whatever the approach followed, all of them need to include some specific features defining the spatial correlation of sites in the experimental or analytical dynamic approximation of the ecological processes under study. Analytical approaches are based on pair-approximation (Ives et al. 1998; Hiebeler 2000; Ovaskainen et al. 2002) or nearest-neighbour correlation techniques (Bascompte 2001). Spatially explicit simulations run on correlated lattices (Hill & Caswell 1999; With & King 1999a; King & With 2002). Patch-based (incidence function) metapopulation models rely on patch network structures where population dynamics takes place at patches of different sizes, and then, have to deal with specific information from every patch (Hanski 1998).

The inability of some of these models to scale up or down within the landscape is nonetheless one of their main drawbacks. Such processes as foraging and dispersal ability are clearly species-specific and may be related to body-size or trophic position (Ritchie 1998; Haskell et al. 2002). Foraging or dispersion responses of species depend on the way they perceive resources on the landscapes (Milne 1991; Ritchie 1997; McIntyre & Wiens 2000). While incidence function techniques (e.g. Hanski & Ovaskainen 2000) have been successful in dealing with species-specific dispersal attributes by assigning different dispersal kernels to the species at play, the use of matrices in the definition of patch attributes usually requires a large computational task.

In ecological processes, when scaling issues are involved, techniques concerning fractal approaches are often useful. By their very nature, fractals involve in many cases invariant properties across many orders of magnitude (Kunin 1998; Ritchie 1998; Ritchie & Olff 1999; Keitt 2000; Haskell et al. 2002; Halley et al. 2004). This invariance has often proved to be related to the existence of single mechanisms acting over different scales, i.e. patch selection is clearly scale-dependent, and resources are more or less available depending on the local knowledge a species has on its environment (Ritchie 1997). If fractal assumptions hold over a wide range of scales, one can define the density of the local landscape via simple analytical relationships involving a window size for perception, the fractal dimension of the landscape and some normalization constant. This normalization constant, the so-called lacunarity, is an expression of the structural contagion of the landscape (Olff & Ritchie 2002) and has been previously stated to correlate with extinction thresholds for metapopulations (With & King 1999b), although no explicit incorporation of this parameter has been done so far. Many criticisms have been made against the main assumption involving monofractal relationships, stating that instead, real landscapes are quasi-fractals, where some typical scaling holds over a finite range (Lennon et al. 2002) or multifractals (Halley et al. 2004), where the momentums of fractal dimensions are not the same and some scales dominate others.

Hanski & Ovaskainen (2000) defined a neat relationship between population persistence, extinction thresholds and life-history traits on one side (such as colonization, extinction and dispersal capacity parameters in metapopulation models) and habitat quantity and contagion on the other. In their incidence function approach, they coined the term metapopulation capacity as the leading eigenvalue in a matrix of spatial connectivities between patches. Ovaskainen et al. (2002) further proved the same relationship in lattice-based correlated landscapes. In this study I show that this term is simply related to the lacunarity of the landscape. Furthermore, lacunarity can be explicitly included into the equation for the equilibrium population and the extinction threshold condition. This is not only applicable to strongly fragmented metapopulations, but to any lattice-based multifractal landscape at use.

2. The model

Following former approaches (Lande 1987; Noon & McKelvey 1998; Hanski 1999; Hill & Caswell 1999; Hanski & Ovaskainen 2000), the basic foundation of this study stands on the modified Levins model of patch occupancy in fragmented landscapesEmbedded Image(2.1)where c and e are colonization and extinction parameters, h represents the fraction of the landscape composed of suitable patches (occupied+empty), and p is the fraction of suitable habitat h that is occupied. This mean-field model assumes global mixing and, after solving for stability, gives the solution p*=1−δh−1 for the equilibrium fraction of occupied patches, where δ=ec−1 defines a species-specific life-history ratio (sensu Hanski & Ovaskainen 2000).

The colonization rate C=ch(1−p) depends on the quantity of habitat h that a species may detect to disperse into empty sites from occupied sites. There is also a critical fraction of habitat under which populations become extinct, which is given by hc=δ, the so-called extinction threshold (Lande 1987; Dytham 1995; Bascompte & Solé 1996).

However, as earlier mentioned, non-random spatial configurations of habitat may have a huge effect on the dynamic behaviour of this model (Hill & Caswell 1999; Bascompte 2001; Ovaskainen & Hanski 2001). Thus, one can relax Levins' assumption by assuming: (i) local information and (ii) non-random habitat distribution. Let us first consider a landscape consisting of a lattice where species present local resource perception abilities. For simplification, instead of using an exponential decreasing function (Hanski 1999), the species will be able to disperse only to those patches inside some square window with an area a×a surrounding the source patch (a cell in the lattice), where a is the linear length of any of its sides expressed as a fraction of the linear dimension of the whole landscape. a is a species-specific parameter defining the spatial scale of patch interaction. In a purely random landscape in discrete time, the distribution of habitat cells within the neighbourhood a×a is Poisson-determined (Hill & Caswell 1999).

On the other hand, in non-random, correlated landscapes, fractal theory is a straightforward approach to the study of metapopulation dynamics (Mandelbrot 1983; Plotnick et al. 1993; Ritchie 1997; Haskell et al. 2002). In such landscapes, the habitat that a species may detect at local scales follows well known scaling relationships.

Thus, f(a), the average local fraction of habitat that a species locally perceives, can be explicitly incorporated into the average colonization rate within equation (2.1). Now, it reads as C=c(1−p)f(a), whereEmbedded Image(2.2)with Φ being a normalization factor defining the inverse of the number of boxes of size a occupied at least in one cell of minimum resolution size ϵ, and k an index for every box of size a. For landscapes with large amounts of habitat available and low contagion, Φa−2, where the exponent is defined by the embedding dimension (in this case, a two-dimensional landscape). That is, Φ is in that case the total number of boxes of size a filling the landscape, therefore f(a)≃h. In this step, we relax the assumption of a zero-variance distribution under the implicit model, and modify the quantity of local available habitat to incorporate into the general metapopulation dynamic equation (2.1). It can be shown (see Appendix A) that f(a) is related to the lacunarity of the landscape, which has been usually considered as a constant normalization factor defining the texture or spatial aggregation of the landscape in a typical scaling relationship. Lacunarity, in fact, is an index of landscape texture accounting for second moments in the probability mass distribution M(n, x) (number of occupied cells per box of size a×a in the set; Plotnick et al. 1993)Embedded Image(2.3)where Embedded Image and Varn(x) are the mean and variance of the number of occupied sites per box.

Thus, it incorporates explicit space into a simple index. When multifractal landscapes are at play, lacunarity, as a measure of heterogeneity, shows scale-dependence (Plotnick et al. 1993). In such a case, the approach shows that f(a) can be explicitly defined through the computation of lacunarity on the real landscapes. Thus, the final equation results in a semi-implicit approach. Now the solutionEmbedded Image(2.4)is straightforward. Likewise, the extinction threshold becomes fc(a)=δ. Substituting equation (A 4) into equation (2.4) the averaged extinction threshold (for p*=0) now readsEmbedded Image(2.5)It is interesting to note the astonishing resemblance between equation (2.4) and that given by Hanski & Ovaskainen (2000) for the fraction of population at equilibrium, where f(a) is replaced by the metapopulation capacity λM. This parameter is the leading eigenvalue of the connectivity matrix in patch-based (Hanski & Ovaskainen 2000) and lattice-based landscapes (Ovaskainen et al. 2002). In the latter, they found thatEmbedded Image(2.6)where R stands for the local connectivity a species may perceive, and is a function of some dispersal parameter α defining the rate at which dispersion of propagules decays with distance. Thus, R∈{0,1} is clearly related to f(a).

Observe that the bracketed term is a good approximation to equation (2.3). When dealing with largely decorrelated landscapes, landscapes with a large amount of suitable habitat, or species with high dispersal capabilities, E(R)≃h. From this approximation, together with equation (A 4), it follows the simple form λMf(a). That is, metapopulation capacity and landscape lacunarity are formally equivalent and directly applicable to the dynamic metapopulation equations.

A special case holds when self-similar, so-called monofractal or homogeneous landscapes are concerned. Those landscapes, although never accomplished in the real world (Halley et al. 2004), are a good basis in the quest for formally based ecological scaling laws (Ritchie 1997; Ritchie & Olff 1999; Haskell et al. 2002). Because of the specific relationship by which the box-counting dimension holds constant over all scales (i.e. the well known scale-invariant features), there is a tight relationship with strong covariance between the box-counting dimension and the habitat available (Olff & Ritchie 2002; Halley et al. 2004). The very nature of this relationship allows one to discern formal relationships for the populations at equilibrium or the extinction thresholds (see Appendix B). By taking equation (B 2) it is easy to get the solutions for f(a) (see figure 1b)Embedded Image(2.7)which can be further applied to obtain equilibrium populations (figure 1a)Embedded Image(2.8)Finally, the average extinction threshold isEmbedded Image(2.9)

Figure 1

Response to habitat availability of lacunarity-dependent variables in monofractal, self-similar landscapes: (a) equilibrium occupancy, (b) habitat perceived. The responses are shown for different window sizes 2i−1256−1, where i represents the number of the line in top-left/bottom-right direction. Then the top-left line is assigned for a=ϵ and the bottom-right for a=1.

3. Results

I tested the numerical predictive power of the lacunarity approach for the metapopulation model by building two-dimensional correlated landscapes with the use of the midpoint displacement algorithm (MDA) for three-dimensional rugged landscapes (Saupe 1988; Hill & Caswell 1999; Halley et al. 2004) and then using slices cutting off those landscapes horizontally in order to fit the desired quantity of habitat available in the landscape. The method recursively draws imperfect, rugged surfaces based on fractional Brownian motion controlled by two input parameters: the Hurst exponent H=2−D, where D is the box-counting dimension, and the quantity of suitable habitat h. I simulated 10 replicates and computed the lacunarity and f(a) for each of the replicates. Average values are shown in figure 2 for different window linear sizes. Several points can be observed from the figure. (i) The MDA is unable to produce monofractal landscapes, since there is a constraint motivated by the need of both input parameters D and h, and by finite size effects, so the imposed scaling relationship from monofractals must be relaxed. That means that D only holds for some small part of the range a∈{ϵ, 1}, and that this scaling is highly dependent on the h chosen. Then, the dominant fractal dimension gradually changes to avoid hitting geometrical constraints earlier or later than expected (a=1), where equation (A 4) becomes f(a)=h. This is the largest window size, where all approaches collapse. Further discussion on this topic can be found in Allain & Cloitre (1991), Lennon et al. (2002) and Halley et al. (2004). (ii) The quantity of available habitat perceived for colonization decreases as the window size increases. This decrease is more pronounced at very local scales in random landscapes, and at the larger ones in the more aggregated landscapes.

Figure 2

Fraction of habitat perceived by species in multifractal landscapes for different aggregation scenarios, with ϵ=256−1. Recall that H=2−D. Decay of the fraction f(a) (and the lacunarity) as a function of window size a for two different h scenarios and different levels of aggregation.

Lacunarity values computed from the generated landscapes were used to determine equilibrium populations under equation (2.4). For every block of 10 replicates, confidence intervals (α<0.05) were also computed. Since the MDA stands on stochastic settings, different initial seeds may produce very different landscapes under the same D and h, which produces very different landscapes and lacunarity values. Average values of the trivial equilibrium populations p* resulting from the mean-field lacunarity-based (MFLB) model, and their respective confidence intervals can be seen in figure 3. Furthermore, I run a spatially explicit simulation (SES) version of the metapopulation model on those same landscapes in order to test the accuracy of the MFLB model. Ten different runs of the spatially explicit simulation model also accompany the former semi-analytic approach in figure 3.

Figure 3

Effect of habitat availability on the fraction of occupied patches for different aggregation scenarios and window sizes. Thin dashed lines depict the results of the original implicit, mean-field model approach. Black lines and grey lines stand for the average and confidence intervals (α<0.05) in the MFLB model described in this study. Circles are the results for 10 different replicates for every value of h in the spatially explicit simulation approach. MFLB and SES model results are based on stochastic landscapes generated with the MDA. See more details in the text.

I computed landscapes with ϵ=256−1. The metapopulation model had parameter values c=0.2 and e=0.125. SES simulations were stopped when the absolute value of the slope in the least-square regression line between the successive last 20 arbitrary time steps was less than 10−4. Lattice boundary conditions were wraparound or periodic (Hiebeler 2000). In order to account for the specific effect of landscape heterogeneity (and hence lacunarity), decoupled from that of contagion by local dispersal, I used the infinite-dispersal approximation, where habitat perception is local, but propagules randomly disperse across the landscape after every time step of integration (Durrett & Levin 1994; Hiebeler 2000). This shuffling of propagules ensure avoidance of spatial correlations in cell occupancy generated by the population dynamics (Hiebeler 1997). Since the metapopulation model with explicitly incorporated lacunarity is still a mean-field model, it does not account for dispersal-generated correlations, but only for those associated to the local landscape heterogeneous structure. Therefore, the model is better suited for host–parasitoid systems where propagules may well be localized in hosts that are aggregated in small-scale areas, and then disperse through larger areas (J. M. Montoya, personal communication). Thus, simulations must ensure that spatial correlations are purely structural or landscape-generated (Ives et al. 1998; Hiebeler 2000; Ovaskainen et al. 2002).

The results depicted in figure 3 show that simulations and the mean-field lacunarity-based approach fit very well for less aggregated landscapes, larger window sizes and more habitat available. Exceptions hold with small window sizes and largely aggregated landscapes (H→1), where the mean-field lacunarity-based approach slightly overestimates the simulation results. This might be an artefact of demographic stochasticity in the simulations, since very small amounts of habitat are more sensible to stochastic deaths in small patches (Andrén 1994; Ives et al. 1998). Under those exceptions, there is almost no extinction threshold, but it could reappear again if Allee effects were at play (Noon & McKelvey 1998; Liebhold & Bascompte 2003). This effect has also been observed in metapopulation capacity and even pair-approximation based approaches (Ives et al. 1998; Hiebeler 2000; Ovaskainen & Hanski 2001). An important observation is that in random, independent landscapes, the three approaches (i.e. Levin's mean-field, lacunarity-based mean-field and stochastic simulation) match and have the same behaviour independently of the window size used. This is accomplished because in random landscapes the lacunarity of the set is the same for a wide range of window sizes. However, when aϵ, f(a) decreases very quickly (figure 2), one should expect no agreement in those random landscapes when dealing with species with very localized perception.

4. Discussion

Mean-field models are intended to replace all interactions between patches with an average or effective interaction. In this sense, I have stressed the difference between a purely implicit mean-field model and the lacunarity-based, mean-field approach, which is nothing but the computation of the average ‘effective’ interaction from local scales of perception, which are computed from local values building the specific fractal mass probability density functions at every window size. Several points may be considered.

(a) Lacunarity and allometric habitat perception

The approach I have developed through this study intends to show that the lacunarity of the landscape is not only a good landscape index related to extinction thresholds (Plotnick et al. 1993; With & King 1999b), but may also be explicitly incorporated into general dynamic metapopulation equations. It adequately incorporates the complexity of heterogeneous, multifractal landscapes into a single ESLI (ecologically scaled landscape index; sensu Vos et al. 2001) parameter which can then be used as a factor multiplying the quantity of habitat that species with different dispersal strategies may perceive. In this sense, fractal body-size scaling relationships may well be related to population dynamic approaches through very simple modelling exercises, since body-size has already been associated to the stride length or dispersal capacities of species (Ritchie 1997; Haskell et al. 2002).

The scale-dependence of lacunarity seems to further relax some of the assumptions taken in ecological models with body-size scaling relationships incorporated. Mostly, they tend to assume self-similar landscapes with constant lacunarity at all scales (Haskell et al. 2002). However, we have seen that, depending on the landscape, the decrease in perceived habitat may be very different, only constrained at very extreme geometrical configurations (Hiebeler 2000; Halley et al. 2004). This means that the tight use of constant fractal dimensions to elucidate general scaling laws may well need a revision.

(b) Lacunarity and metapopulation capacity

A straightforward advantage in the use of the lacunarity-based approach is its computational feasibility in opposition to the difficulties of computing leading eigenvalues in lattice-based or even patch area-based landscape matrices. I have shown that metapopulation capacity and the scale-dependent lacunarity index are formally equivalent under geometrical transformation of a continuous to a discrete habitat configuration. Thus, lacunarity can easily substitute the metapopulation capacity in metapopulation dynamic models. For incidence function models, it would be desirable to redraw the landscape in pixelized form to compute lacunarities and at the same time, avoid the use of individual patch areas, since, as I have shown, the normalized lacunarity in equation (2.3) is nothing but a good approximation to the average species-specific (ecologically scaled) connectivity of the landscape. This functional connectivity has already been stressed as an important component in predicting extinction thresholds (With & King 1999a), different from those found on models with global mixing (Adler & Nüernberger 1994).

Regarding pair approximation, the lacunarity-based approach has the same advantages as the metapopulation capacity. It can deal with a wide range of dispersal or structural distances, while pair-approximation techniques can only deal with nearest-neighbour correlations. On the other hand, the pitfalls of the lacunarity and the metapopulation capacity approaches are their inability to predict the dynamics of aggregation in populations with small dispersal abilities in very small, aggregated habitats. This is the main reason why the lacunarity-based approach can not be used in multitrophic levels of food web metapopulation models, where trophic correlations add up to behavioural effects (King & With 2002), such as foraging or dispersal, to depict patterns very different from the ones given by the structure of the landscape. I have generated simulations (not shown) in a lacunarity-based metapopulation predator–prey model (Holt 1997; Hanski 1999) and the results hardly fit in situations where local, dispersal-driven spatial correlations are of much importance (small local neighbourhoods, aggregated landscapes, or small quantities of habitat available).

On the other hand, the use of lacunarity may have strong ties to some other modelling approaches as the ones used in individual territory models, such as those using equations with logistic-like growth and multiple searching opportunities (Noon & McKelvey 1998; With & King 1999b; Gamarra & Solé 2002). Under this framework, there exists a trade-off between dispersal distance and colonization success or between search and encounter rates (Ritchie 1998; With 2001). This compromise could be mechanistically explained if reproductive constraints are at play in individual-based models under lacunarity-based habitat perception. A positive, monotonic, nonlinear relationship exists between both parameters, (i.e. lacunarity and number of searching trials), although the exact formulation still must be developed.

(c) Contagion by dispersal versus contagion by landscape structure

First, this study shows that there is a decoupling between the correlations induced by landscape heterogeneity and those induced by dispersal (those usually increase the value of p* even further). Second, contrary to the statements formulated by King & With (2002), I may conclude that spatial pattern is important in determining population success even in landscapes with h>0.4. The effect of the box size lowers the percolation threshold in most cases, and consequently, the extinction threshold (With & King 1999a) practically vanishes in very aggregated landscapes (Hill & Caswell 1999), as noted in figure 3. This pattern has been formerly found in multidimensional, scale-free networks dealing with epidemiological models (Pastor-Satorras & Vespignani 2001; Bascompte 2003). A recent study (Kallimanis et al. 2005) shows that, besides dispersal and landscape structure, disturbance structure is also important in determining extinction thresholds. A semi-implicit lacunarity-based model may be a good approach whenever correlated disturbances around occupied patches are included in metapopulation models. This would be achieved by including lacunarity in the otherwise Poisson-distributed extinction rate. Since both the extinction and colonization rates play opposite effects in the computation of the extinction threshold, a fair hypothesis would point to the ratio of disturbance–colonization lacunarities as the tuning parameter defining extinction thresholds. A larger ratio would increase Embedded Image, in agreement with Kallimanis et al. (2005).

Finally, Halley et al. (2004) explicitly stated that there is little theory relating lacunarity to generating processes (if any at all). Although more work is still needed, this work represents, at its best, a good starting point to understand lacunarity and other ecologically scaled landscape indices (Vos et al. 2001) in ecological dynamic frameworks.


Partial support was given by the Cornell University Biogeochemistry and Biocomplexity Initiative. David Alonso has contributed with some mathematical advice. I am deeply grateful for interesting discussions to Mark Ritchie, Ricard Solé and the people at the Complex Systems Lab at UPF in Barcelona. The Spatial Ecology Program of the Division of Population Biology, University of Helsinki, gave me some financial support for a short visit. Finally, The ‘core group’ of Spaniards at Cornell University gave me the necessary motivation to write the manuscript.


    • Received March 14, 2005.
    • Accepted April 17, 2005.


View Abstract