Metapopulation dynamics: effects of habitat quality and landscape architecture
Atte Moilanen
INTRODUCTION
Ongoing fragmentation and destruction of natural habitats are of great concern to conservationists and ecologists throughout the world. A small local population isolated from other conspecific populations is prone to local extinction, but the species may have a chance of survival in a network of habitat patches connected by dispersal. Ecologists have asked about the most appropriate and effective approaches to the study of populations in fragmented landscapes (Durrett and Levin 1994, Tilman and Kareiva 1997). The currently most popular approach is based on the metapopulation concept (Levins 1969) and on the study of metapopulation dynamics (for reviews see Gilpin and Hanski 1991, Hanski 1994b, 1997, Hastings and Harrison 1994, Hanski and Gilpin 1997).
The term metapopulation has been used in a variety of meanings (Hanski and Simberloff 1997). In our application a metapopulation is an assemblage of local populations inhabiting spatially distinct habitat patches. A conceptually important assumption is that all local populations have a significant risk of extinction. This assumption leads to a focus on extinction-colonization dynamics, in other words, the metapopulation is assumed to persist in a stochastic equilibrium between local extinctions and colonizations of currently empty but suitable habitat patches. Migration is assumed to be distance dependent (Hanski 1994a), and hence the spatial configuration of the landscape affects metapopulation dynamics unlike in the well-known Levins model (Levins 1969), in which only the fraction of presently occupied habitat patches has an effect on colonization. Also, unlike the Levins model, we allow for differences in the areas of habitat patches. In contrast to another commonly studied metapopulation scenario, the mainland-island metapopulation structure (MacArthur and Wilson 1967, Harrison 1991, 1994), the structure assumed here has no single-habitat patch large enough to guarantee the long-term survival of the respective local population and thereby the metapopulation as a whole.
Recently, we have studied metapopulation processes with a stochastic patch occupancy model called the incidence function model (Hanski 1994a, 1997, Hanski et al. 1995b, Moilanen and Hanski 1995, Moilanen et al. 1998). This model is practical in the sense that it has only a few parameters, and realistic in the sense that the model includes the biologically most significant processes: the probability of local extinction depends on population size via patch size, migration and hence colonization of empty patches is distance dependent, the spatial structure of the landscape is explicitly accounted for, and the so-called rescue effect (Brown and Kodric-Brown 1977) is included in the model. Furthermore, the basic incidence function model uses only two landscape variables, patch area and isolation, which are both easy to measure in practice. The incidence function model has been successfully applied to a large data set, covering the entire geographical range of the butterfly Melitaea cinxia in Finland, with [greater than] 1600 habitat patches within an area of 50 x 70 [km.sup.2] (Hanski et al. 1995a, 1996b). The model has also been applied to other species (Hanski 1992, Cook and Hanski 1995, Wahlberg et al. 1996, Moilanen et al. 1998), and it has been used in the study of several particular questions of practical and theoretical importance (Moilanen and Hanski 1995, Hanski et al. 1996a, b).
While attempting to predict the observed pattern of occurrence of the butterfly M. cinxia in its fragmented environment, some regional differences between model predictions and observations were found (Hanski et al. 1996b). There are several possible reasons for such discrepancies. (1) If the landscape structure has been recently altered, the metapopulation may not have had time to reach a new equilibrium (Lamberton et al. 1992, Hanski et al. 1996a). This phenomenon was termed nonequilibrium metapopulation dynamics by Harrison (1991) and Hanski et al. (1996a), and it may lead to what was dubbed “debt of extinction” by Tilman et al. (1994), i.e., additional extinctions that will occur in the future because of past habitat destruction. (2) The level of temporal and spatial variation in the occurrence of the butterfly in the habitat patches may depend on the interaction with a specialist braconid parasitoid, Cotesia melitaearum (Lei and Hanski 1997). This interaction may ultimately necessitate a metapopulation approach explicitly accounting for host-parasitoid metapopulation dynamics (Taylor 1991, Harrison and Taylor 1997, Nee et al. 1997). (3) The theory of structured metapopulation dynamics (Hanski 1985, Gyllenberg and Hanski 1992, Hanski and Gyllenberg 1993, Hastings and Harrison 1994, Gyllenberg et al. 1997) predicts that metapopulations with strong rescue effect may exhibit two alternative stable equilibria, one with most of the available habitat occupied and the other one corresponding to metapopulation extinction. In the presence of environmental stochasticity a metapopulation with multiple equilibria may suddenly crash to extinction from a state with practically all habitat occupied. The M. cinxia data set has revealed convincing evidence for alternative equilibria, which present general problems for the parameterization of metapopulation models (Hanski et al. 1995b). (4) There may be spatial variation in habitat quality and in the landscape structure affecting extinctions and colonizations, which is the focus of the present study.
The starting point of this study, the basic incidence function model, takes implicitly into account some environmental factors and assumes a certain landscape structure. A piece of landscape is assigned as a “habitat patch” based on the ecological requirements of the focal species, whereas the rest of the landscape is assumed to be uniform nonhabitat. It follows that if habitat quality and landscape structure are homogenous across the patch network, the assumptions of the basic incidence function model are satisfied, and the model should perform well. However, if there are regional differences in habitat quality or landscape structure, it would be useful to find out whether these differences are substantial enough to affect the dynamics of the metapopulation as a whole. In other words, a fundamental question is to what extent metapopulation dynamics can be explained by the two basic variables in metapopulation dynamics, habitat patch area and isolation, and to what extent other factors possibly improve the predictive power of metapopulation models.
Spatial variation in the environment can conceivably affect metapopulation dynamics in many different ways (Murphy et al. 1990). At the level of a single patch, habitat quality may affect both the probability of local extinction and the probability of the colonization of an empty patch. Landscape structure between the patches may affect the migration rate in general (Aberg et al. 1995), and it may affect the connectivity of particular pairs of patches (Wiens 1997). One potentially significant mechanism affecting population turnover is interaction between habitat quality and climatic perturbations. For instance, a meadow that is generally too dry for caterpillars may be just optimal during a wet year (Murphy et al. 1990, Weiss et al. 1993, Hanski et al. 1996b). Discovering and especially parameterizing such interactions for the purpose of modeling will be difficult in practice. Other environmental factors have simpler effects that can be modeled more easily. For instance, our previous studies have suggested that the probability of occurrence of M. cinxia is smaller in a grazed than ungrazed meadow (Hanski et al. 1996b). This effect is probably due to grazing animals accidentally eating and trampling larvae, which can be modeled simply as an increased risk of extinction of populations in grazed meadows.
This study consists of two parts. We first describe a method of adding environmental factors affecting habitat patch quality into the basic incidence function model. A set of such environmental factors in the M. cinxia data set is examined, and the significance of these factors is investigated both statistically and by comparing the results with independent biological information about these factors. Secondly, we attempt to include information about the landscape structure as revealed by satellite images in the incidence function model. The idea here is to use knowledge of the type of the habitat between the patches to modify effective interpatch distances, which in the model translates to modified patch connectivity.
METHODS
The incidence function model
The incidence function model (Hanski 1994a) is based on a first-order, linear Markov chain model for changes in the state of one patch, which may be either occupied or empty. Changes in the state of patch i are determined by extinction and colonization probabilities [E.sub.i] and [C.sub.i], which are calculated independently for every patch i in each time unit. The model uses patch areas [A.sub.i] and patch isolations [S.sub.i] as variables that determine [E.sub.i] and [C.sub.i]. Specifically, the area of a patch is assumed to determine the extinction probability, [Mathematical Expression Omitted], where e and x are two parameters. The colonization probability of patch i is defined by a sigmoid function [Mathematical Expression Omitted], where the isolation (actually connectivity, but we use isolation to confer to common practice) [S.sub.i] of patch i is defined by
[Mathematical Expression Omitted] (1)
where [p.sub.j] equals 1 for occupied and 0 for empty patches, [d.sub.ij] is the distance between patches i and j, and [Alpha] is a constant setting the distance-dependent migration rate. [A.sub.j’]s are included in Eq. 1 because the sizes of the source populations affect the absolute migration rate to the focal patch i, and the exponent b is used to scale patch areas to expected population sizes. In this study we use b = 0.5 based on Hanski et al. (1996b). The model thus has five parameters, x, y, e, b, and [Alpha].
The model parameters x, y, and e can be estimated from a “snapshot” patch occupancy pattern by regressing [p.sub.i] against the steady-state expression for the probability of presence in patch i, the incidence [J.sub.i] (Hanski 19948):
[Mathematical Expression Omitted]. (2a)
With the rescue effect (a population may be rescued from extinction by immigration enhancing local population sizes) added to the model, the equation for [J.sub.i] becomes (Hanski 1994a)
[Mathematical Expression Omitted] (2b)
where e[prime] = e[y.sup.2]. Some additional information is now needed to separately estimate the values of e and y in Eq. 2b. For instance, one may estimate the minimum patch area [A.sub.0], for which [Mathematical Expression Omitted], or use information about turnover events between two or more time steps. Note also that in model simulations that include the rescue effect, [E.sub.i] is replaced with (1 – [C.sub.i])[E.sub.i]. For more details about the model and parameter estimation see Hanski (1994a, 1997) and Moilanen et al. (1998).
Adding environmental factors to the incidence function model
In the context of the incidence function model, a natural way to think about the impact of additional environmental factors is to assume that they affect either the extinction or the colonization probability in the model. A factor affecting patch quality, such as the amount of food or the number of nesting places, will affect the expected size of the local population in a patch. In the model such an effect will translate to a modified extinction risk, and a practical way to achieve this is to use patch quality to modify effective habitat patch area. The advantage of this approach is that it will simultaneously take into account the effect of patch quality on the contribution that the patch makes to migration. A factor specifically affecting emigration from or immigration to a patch can be included via modified patch isolation. In the following model, we therefore replace [A.sub.i] by a modified area [U.sub.i][A.sub.i], and the contribution of patch j to the isolation of patch i, [S.sub.ij], is replaced by a modified isolation [V.sub.i][S.sub.ij][V.sub.j]. The modifiers for area and isolation, [U.sub.i] and [V.sub.i], are functions of the environmental factors measured in patch i.
Let us denote by [G.sub.r], r = 1, 2, . . ., R and [H.sub.s], s = 1, 2, . . ., S the functions that map R factors affecting patch area and S factors affecting isolation to the corresponding modifiers for patch area and isolation. Similarly let [g.sub.ri] and [h.sub.si] be the values of the factors measured in patch i. We used third-degree polynomials (e.g.,f(x) = a[x.sup.3] + b[x.sup.2] + cx + d) to model [G.sub.r] and [H.sub.s]. In the simplest case, where only one environmental factor affects patch area and/or isolation (R, S = 1), the functions for the multipliers, [U.sub.i] and [V.sub.i], become
[Mathematical Expression Omitted] (3)
and
[Mathematical Expression Omitted]. (4)
We return to the case with R or S [greater than] 1 below. Eqs. 1 and 2b augmented by environmental factors now become
[Mathematical Expression Omitted] (5)
and
[Mathematical Expression Omitted]. (6)
Note that here the factors affecting isolation have been measured in the patches, which assumes knowledge of habitat quality in the habitat patches and their immediate vicinity only. Thus we modify [S.sub.i] itself, not the term exp(-[Alpha][d.sub.ij]) in Eq. 5, which would require knowledge of the intervening habitat composition between the patches (this will be done in the second variant of the basic model, below). In Eq. 5 a factor affecting colonization will affect similarly both emigration and immigration, which is justified in many cases. For instance, forest surrounding a patch is likely to decrease both immigration and emigration in much the same way. If a factor affects only patch detection during migration, and hence immigration to a patch, Eq. 5 must be trivially modified.
Third-degree polynomials were used to model the functions [G.sub.r] and [H.sub.s] for all the environmental factors used in this study. Third-degree polynomials, unlike for example the logistic function, can accommodate most plausible functional forms, including linear, asymptotic, and humped curves. One could, of course, use tailored functional forms for each environmental factor, but it is convenient from a practical point of view to use the same general function for all factors. Here we are mostly interested in the shapes of the fitted functions, not the parameters of the functions [G.sub.r] and [H.sub.s] themselves.
Parameters x and e[prime] and the coefficients of the functions [G.sub.r] and [H.sub.s] were estimated by maximum likelihood regression of the empirically observed patch occupancies [p.sub.i] against the predicted incidences [J.sub.i] (Eq. 6):
[Mathematical Expression Omitted]. (7)
Eq. 7 was minimized with simulated annealing using 4000 function evaluations for the basic model and an additional 7000 evaluations for each additional environmental factor (for details see Moilanen 1995). Note that parameter [Alpha] can be estimated by setting it a free parameter and using Eq. 7.
We evaluated each environmental factor as follows. (Note that below P denotes statistical significance, p the proportion of occupied patches, and [p.sub.i] the observed incidence of patch i. These symbols have not been changed because they are commonly used in the literature.) Let us denote by [[Theta].sub.b] and [[Theta].sub.a] the error in the fit of the basic model and the model including the additional environmental factor, respectively. We define the significance of a factor as the probability that a model fit error as low as [[Theta].sub.a] is not a chance effect given the distribution of the environmental factor among the patches. Let [[Theta].sub.rk] be the model fit error when the empirically observed values of the environmental factor have been randomized among the patches. Randomization was repeated 100 times to obtain values [[Theta].sub.rk], = 0, 1, . . ., 100. The approximate significance (P value) of a factor is the number of times the original fit error, [[Theta].sub.a], is higher than that obtained from the randomized distribution for the environmental factor, that is, P = #{[[Theta].sub.a] [greater than] [[Theta].sub.rk], k = 0, 1, . . ., 100}/100. The decrease in the error following the addition of the environmental factor, [[Theta].sub.b] – [[Theta].sub.a], was used to measure the power of the factor. Note that if an environmental factor with no significance at all is added to the model, the fitting error of the model will still decrease, as the model has become more complex (more parameters). The effect of adding unnecessary complexity to the model can be seen by examining [[Theta].sub.r], which we define as the average error in the fit of the model after the environmental factor has been randomized among the patches, that is, [[Theta].sub.r] is the average of [[Theta].sub.rk] values. One should also recognize that even if the model would predict the expected long-term patch occupancy probabilities perfectly, the model fit error [Theta] does not go to zero. For example, if a patch is occupied half of the time then it necessarily makes a contribution of -2 x 0.5 x ln 0.5 [approximately equal to] 0.69 to [Theta]. The expected minimum error for the model fit, [[Theta].sub.min], can be calculated from Eq. 7 by assuming perfect prediction ability, that is [p.sub.i] = [J.sub.i]. By using the parameter values obtained for the basic model ([Alpha] = 0.332, x = 0.359, y = 9.10, e = 0.176 with error [[Theta].sub.b] = 793.3) we obtained [[Theta].sub.min] = 599.0. The power of the factor is now defined as the proportion of the explainable error that the factor actually explains, [r.sub.e] = ([[Theta].sub.b] – [[Theta].sub.a])/([[Theta].sub.b] – [[Theta].sub.min]). Noting that we expect the model fit error to drop from [[Theta].sub.b] to [[Theta].sub.r] even if the factor is completely insignificant, we can calculate the power of the significant improvement in the fit of the model as [r.sub.s] = ([[Theta].sub.r] – [[Theta].sub.a])/([[Theta].sub.r] – [[Theta].sub.min]).
Finally, the parameter estimates obtained for the augmented model were used to simulate the Markov chain model to check the dynamic behavior of the model. Model fit error during iteration, [[Theta].sub.it], gives a measure of how well the observed and predicted patch occupancy patterns match during model iteration. [[Theta].sub.it] was calculated as the average error in 10 replicate simulations. Each simulation produced an error score that was calculated with Eq. 7, using J values calculated as the fraction of time that each patch was occupied during the simulation run. The advantage of calculating J values from simulation is that unlike Eq. 6, model simulation takes into account any variation in the [C.sub.i] values, resulting from extinction-colonization stochasticity. Each simulation lasted for 600 time units, with the first 100 time units omitted to allow the metapopulation to settle close to the stochastic equilibrium. For a summary of the different errors used in calculations see Table 1.
When multiple environmental factors are added to the model, there is an additional problem in parameter estimation, because the factors occur in product terms in Eqs. 8 and 9 (below) and therefore cannot be estimated independently. Suppose, for example, that two factors affect patch quality (area). If the function for the first factor is multiplied by 0.5 and the function for the latter by 2, the model still gives the same result as with the original values. Thus with multiple factors the shapes of the functions G and H are significant, whereas the absolute values are not. This situation may lead to some problems with the convergence of the numerical optimization algorithm. One way to resolve the problem is to force the multiplier equal to 1 for a patch with the median value for the factor. This fixes the functions and unambiguous parameter values can be found. However, since the true multiplier for a median-valued patch is not necessarily I for all the factors, one should add parameters ([[Lambda].sub.1], [[Lambda].sub.2],) setting the absolute levels of the products. With these modifications, equations for including multiple environmental factors become
[Mathematical Expression Omitted] (8)
and
[Mathematical Expression Omitted]. (9)
[TABULAR DATA FOR TABLE 1 OMITTED]
Combining GIS data with the incidence function model
The basic incidence function model assumes a simple landscape structure in the form of distinct suitable habitat patches with uniform quality surrounded by entirely unsuitable habitat. In the second variant of the basic model we used standard GIS data to allow for heterogeneous landscape surrounding the patches. It is natural to assume that different habitat types have differently permeabilities to a migrating individuals (Wiens 1997). For example, the Glanville fritillary butterfly does not fly in dense forests, and as the landscape in Aland is a mosaic of forests and open areas, the habitat composition between pairs of patches is likely to affect migration between the patches. Here we assume that the habitat composition between pairs of patches will affect effective interpatch distances in the calculation of isolation. In other words, “difficult” habitat between patches will increase the effective interpatch distance, whereas “easy” habitat will decrease it. The modified equation for isolation is
[Mathematical Expression Omitted] (10)
where [[Alpha].sub.h] is the migration coefficient for habitat type h. The distance between patches i and j crossing habitat type h, [d.sub.hij], was determined by scoring different habitat types in 25-m Landsat data along straight lines connecting the two patches. Note that in this case P, [[Theta].sub.r], and [r.sub.s] cannot be calculated, because the landscape cannot be “randomized” like the values of the environmental factors characterizing the patches can themselves (see Adding environmental factors to . . . above).
The Landsat data available for Aland islands consisted of 44 habitat classes. Many of these classes were represented by an insignificantly small amount of habitat and including them into Eq. 10 would achieve little else than complicate the model. Therefore we first reduced the number of habitat types by classifying them into groups that consisted of habitat types that are ecologically similar and had similar effects in the following analysis. We first scored the proportion of each habitat type within a 400-m circle around the patch. Then, using logistic regression, the observed patch occupancies ([p.sub.i]) were explained by the incidences predicted by the basic model ([J.sub.i]) and by the habitat surrounding the patch. The habitat types that did not have a significant effect in this analysis were subsequently grouped together to represent a “standard” background habitat.
Empirical data
The models were parameterized and tested using a three-year data set (1993-1995) on the Glanville fritillary butterfly, Melitaea cinxia. The data were collected from a large metapopulation on the main Aland island, southwest Finland, with 1264 habitat patches, small dry meadows with the larval host plants Plantago lanceolata and Veronica spicata ([ILLUSTRATION FOR FIGURE 1 OMITTED]; Hanski et al. 1995a, 1996b). The data set includes patch areas, spatial locations, and presence/absence data on the butterfly in the three years. The [p.sub.i] values used in parameter estimation are the average occupancies for the three years. For the small number of patches discovered since the first large survey in 1993, data for only one or two years were used in the calculation of [p.sub.i] values.
In the first part of this study, environmental factors were used to modify habitat patch quality (area). Data for this study was gathered in 1993, when a number of habitat variables apart from area were scored for each habitat patch (Hanski et al. 1996b). In a preliminary analysis, nine factors were identified as potentially important (Table 2), and these factors were recorded again for all the patches during the survey in 1995, Grazing by cattle and sheep was scored as a binary variable; the remaining variables have values ranging from 0 to 100% (Table 2). During numerical optimization, the values of all variables were scaled to the interval [0, 1].
In the second part of this study, the quality of the interpatch habitat was taken into account in the calculation of isolation. The data used for this purpose were based on Landsat TM 25-m satellite images. We used a standard interpretation of these data that was available from the Finnish Geodetic Institute. This interpretation is mainly intended for forestry use and the habitat types it differentiates include water, agricultural fields, rocky outcrops, various kinds of swamps, clear-cut areas and pine, spruce, and deciduous forests of different densities. The data are further described in Vuorela (1996).
RESULTS
Table 3 reports the results for models with one extra environmental factor added to the basic incidence function model. Four factors, grazing, percentage of patch perimeter bordering cultivated fields, percentage of patch perimeter bordering and semi-open habitats, and the amount of nectar-producing plants were found to have statistically significant effects. Of these four factors only grazing ([r.sub.s] = 0.10) and patch perimeter bordering cultivated field ([r.sub.s] = 0.15) have a power that is large enough to merit attention.
[TABULAR DATA FOR TABLE 2 OMITTED]
Fig. 2 shows the fitted functions for the four factors that were found to have a significant effect. Grazing by cattle or sheep greatly decreased effective patch area and hence occupancy. Patch boundary with cultivated fields increased isolation, while patch boundary with semi-open habitat decreased isolation. The results for flowering plants were more complicated. In a preliminary study, when examining the effect of flowering plants on effective patch area and isolation, conflicting results were obtained, with a large improvement in model fit but poor simulation results, indicating a poor match between the observed and predicted patterns of patch occupancy. The problem turned out to be the assumption that flowering plants had a similar effect on immigration and emigration (Eq. 5). Sensible results were obtained when assuming that abundant flowering plants increase immigration but decrease emigration; for simplicity we assumed here that the effect on immigration is the inverse of that on emigration.
Fig. 3 shows the functions for grazing, flowering plants, and patch perimeter with cultivated fields when estimated simultaneously. This model had [[Theta].sub.b] – [[Theta].sub.a] = 42.0 and power = 0.22. These results are consistent with those in Fig. 2, where the functions for the factors were estimated separately. The fourth individually significant factor (Table 3), patch perimeter with semi-open habitat, was not included into this model because it was negatively correlated with patch perimeter bordering cultivated fields.
The results for the incidence function model augmented with GIS data were disappointing. We initially used logistic regression to help in the grouping of the numerous habitat types present in the original interpretation of the satellite data. Table 4 shows the best logistic regression model for patch occupancy ([p.sub.i]) explained by the incidence ([J.sub.i]) as predicted by the basic model and by the habitat composition surrounding the [TABULAR DATA FOR TABLE 3 OMITTED] patch. Several significant effects of habitat types were detected (Table 4), but many common habitat types, such as pine forests of different densities, were found to have no effect. When the significant habitat types were added one at a time to the incidence function model (Eq. 10) only very small improvements to the fit of the basic model were obtained (Table 5). The two habitat types with the largest effects, cultivated fields and open rocky outcrops, had powers of only 7.4% and 2.6%, respectively.
Discussion
We have here described and applied methods of adding environmental factors (habitat quality) and landscape structure to a spatially realistic metapopulation model, the incidence function model (Hanski 1994a). This sort of approach has been called for to bridge the gap between metapopulation ecology, which is patch oriented, and landscape ecology, which is focused on the landscape structure (Wiens 1997). As it turned out, when tested with data from the well-studied large metapopulation of the Glanville fritillary butterfly (Hanski et al. 1994, 1995a, 1996b), the additional information on patch quality and the structure of the surrounding landscape did not greatly improve our ability to quantitatively [TABULAR DATA FOR TABLE 4 OMITTED] predict the dynamics of the butterfly in its highly fragmented environment. Below, we discuss the effects of habitat quality and landscape structure in turn, as well as some general issues in the implementation of the present type of approach in practice.
Habitat quality
The present results on the effects of habitat quality on metapopulation dynamics of the Glanville fritillary largely agree with the results of a previous analysis employing logistic regression (Hanski et al. 1996b). This is encouraging, since the two analyses are based on entirely different models: a static regression analysis of patch occupancy and (in the present study) a dynamic metapopulation model. The only substantial difference in the results was the effect of flowering plants on migration detected in the present study but not in the regression analysis. We found that for moderate values (0-0.6) the density of flowering plants decreased emigration but increased immigration. Our confidence in this result is increased by Kuussaari et al.’s (1996) experimental study, which demonstrated the same effects. The actual mechanism is also known: abundant flowers tend to keep the butterflies where they are, and to attract migrating butterflies (Kuussaari et al. 1996). [TABULAR DATA FOR TABLE 5 OMITTED] High density of flowering plants ([greater than]0.6) seems to discourage colonization in the few patches where this many flowers were observed. This result is possibly due to the fact that patches with this high density of flowering plants often have high vegetation with relatively little of the larval host plants, which makes the patch less suitable for breeding and hence decreases the probability of successful colonization. It is encouraging that the dynamic modeling revealed these qualitative effects.
The effects of patch boundary with cultivated fields and semi-open habitats were opposite, which is at first surprising. The negative effect of open fields is most likely due to elevated emigration rate from patches with much open boundary (Kuussaari et al. 1996), which is apparently not compensated for by increased immigration rate. Patch occupancy is thus less likely if the patch has a highly permeable boundary, as suggested by much of the landscape ecological literature (reviewed in Wiens 1997). In contrast, the percentage of patch perimeter bordering semi-open habitats had a positive effect on between-patch migration (decreased isolation). A static regression analysis detected the same effect (Hanski et al. 1996b). We do not know for sure why this effect is observed, but suspect that it is related to various semi-open habitats being structurally similar to the meadows themselves, and thereby possibly functioning as movement corridors and facilitating migration.
Grazing by cattle and sheep decreases patch quality, because of increased larval morality due to foraging and trampling by grazing animals (M. Kuussaari, personal communication). In our analysis, grazing greatly decreased effective patch area and thereby substantially increased the probability of local extinction. Nonetheless, it is worth emphasizing that grazing has a potentially positive long-term effect in keeping the vegetation relatively open and thus favoring the growth of the larval host plant P. lanceolata, a relatively poor competitor in the study area (Hanski et al. 1996b). Similar effects of grazing have been observed in other studies on butterflies (Murphy et al. 1990, Thomas and Jones 1993, Thomas and Hanski 1997); grazing is generally beneficial for species inhabiting early-successional habitats.
The environmental factors that we have investigated here with the exceptionally large data set on the Glanville fritillary turned out to have relatively little effect on the overall performance of the incidence function model. The most significant factor had a power of only 0.15. This result implies that in this case no single environmental factor had a great effect compared with the primary effects of patch area and isolation. In a sense, this is good news, because a key justification for the basic incidence function model is simplicity and ease of parameterization while maintaining realism. Each additional environmental factor would make the gathering of data so much more difficult and expensive and would also increase the size of the data set needed for reliable parameter estimation. However, we add two important caveats to these conclusions.
The results could have been very different for some other species and other landscapes; there are no theoretical reasons to expect that habitat quality or landscape structure have no major effects on metapopulation dynamics. The Glanville fritillary in Aland exemplifies a metapopulation living in a highly fragmented landscape, in a network of generally very small patches, with a median size of only 300 [m.sup.2] but a large range from 12 [m.sup.2] to 6 ha (Hanski et al. 1995a). One could expect that the effects of patch area and isolation are especially strong in this sort of a metapopulation (Andren 1994). The habitat patches were originally identified on the basis of the presence of the larval host plants, but they also happen to be rather uniform in quality. In other species and landscapes, habitat quality may vary more among patches and hence have a greater effect on the dynamics.
Our second caveat is about the frequency of different types of habitat patches in the network. When considering the potential effect of an environmental factor one should observe both the shape of its estimated function and the empirical distribution of the values for the factor. For grazing a power of 0.11 was obtained with only 15% of the patches being affected by grazing [ILLUSTRATION FOR FIGURE 2 OMITTED]. On the other hand, the effective area of a grazed patch is only 6% of the effective area of an ungrazed patch of the same size [ILLUSTRATION FOR FIGURE 2 OMITTED], which indicates a strongly decreased probability of occupancy for grazed patches. Thus grazing would have a very strong effect on metapopulation dynamics in this species if grazed patches were more numerous. In this sense the effect of grazing is stronger than the effect of patch boundary with cultivated fields even if the latter explains a larger proportion of the error in the model fit in this particular metapopulation (Table 3). As things stand, however, grazing is common enough to affect the persistence of only limited subregions of the M. cinxia metapopulation.
Landscape structure
It is both surprising and disappointing that the inclusion of GIS data into the basic incidence function model failed to improve model fit to any significant degree. There are several possible reasons for this result. Landsat data have a grid size (25 m) that is large in comparison with the spatial scale that is important for the butterfly. The median patch size in Aland is only 300 [m.sup.2] (Hanski et al. 1995a), which is less than half of the area of one grid cell in Landsat data. Second, the habitat classification used for the GIS data may not be optimal for this study. The abundance and distribution of flowering plants are evidently important for the migration of the butterfly, but there is practically no way to extract this kind of information from the current satellite data. These problems could possibly be alleviated by using a tailored habitat classification, aerial photographs, and field work to construct a highly detailed map of the study area. But this will be very expensive and time consuming for large areas. We are in the process of attempting a detailed mapping of an area of 25 [km.sup.2], but to extend this to the entire Aland island (1500 [km.sup.2] of land) is simply not economically possible.
Even if problems of data quality could be solved, the largest problem of all would still be left: one has to make assumptions about the movement behavior of individuals if one wishes to use detailed data on landscape structure in these kinds of models. In this study we scored the habitat composition along a straight line between two patches and used this information to modify effective interpatch distances. We also tried elliptical and rectangular buffer zones instead of straight lines with a similar lack of positive results. What one would really want is a measure of the difficulty of migration between any two patches. Unfortunately, this measure will depend not only on the proportions of different habitat types between the patches but also on the exact spatial configuration of the landscape mosaic. To fully utilize such information one will have to quantify species-specific characteristics such as transition probabilities between pairs of habitat types, movement speeds, and mortalities in different habitats and possibly various edge effects. While being hard to model in a simple and realistic manner such effects are likely to be extremely difficult to quantify and demonstrate in empirical studies.
In a recent paper, Gustafson and Gardner (1996) studied the effect of habitat matrix heterogeneity on dispersal and habitat patch colonization with an individual-based model that used self-avoiding random walkers. In their study 89% of the variability in dispersal success was explained by differences in the sizes and isolations of the habitat patches, though habitat structure could significantly affect emigration success from a particular patch. Gustafson and Gardner (1996) also found that dispersal corridors were often diffuse and difficult to identify from the structure of the landscape. Our results with the Glanville fritillary data set are in agreement with the simulation results of Gustafson and Gardner (1996). Taken together, these results underline the principal role of habitat patch area and isolation in the modeling of spatial dynamics in highly fragmented landscapes.
The one encouraging aspect in the GIS study was that the factor with the largest effect was the amount of cultivated fields (Table 4), which had a similar negative effect on migration in the first part of this study, where the quality of the surrounding habitat was scored only in the immediate vicinity of the patch (Table 3). In the GIS study this effect was weaker (power = 0.074) than in the first analysis (power = 0.15). One possible reason for this difference is that when habitat surrounding the patch is scored by visual inspection from the patch, the inaccuracies of satellite data are avoided. Also, Gustafson and Gardner (1996) observed that movement corridors between habitat patches often become diffuse and hard to identify. This may imply that habitat in the immediate vicinity of a patch is more important in determining migration rates than habitat further away from the patches, because habitat immediately surrounding the patch must be crossed during migration, whereas habitat farther away is less likely to be part of the realized migration route of a particular individual. Hence, using information about the habitat immediately surrounding the patch to modify the effective isolation of a patch may be a useful first step in using landscape information in a metapopulation model. The advantages of this approach are that it is easy to implement, it only requires data that are relatively easy to obtain, and it avoids the general difficulties in the modeling of individual migration. Conversely, any approach aiming to use information on the entire landscape mosaic in a quantitative metapopulation model will necessarily face the formidable obstacles described above.
Parameter estimation
Even though our study was not primarily focused on parameter estimation of metapopulation models, this subject deserves some attention since there are common pitfalls that one should be aware of. When metapopulation dynamics are modeled as a stochastic extinction-colonization process, it follows that the proportion of occupied patches, p, fluctuates up and down from year to year as shown in Fig. 4. Such fluctuations in the absence of any environmental noise can be quite large and they pose an obvious problem for parameter estimation. One has practically no way of knowing from a snapshot occupancy pattern whether a low, average, or high proportion of patches was occupied with respect to the long-term average. For this reason, and if data are available for one or a few years only, parameter estimation from few (especially small) patch networks will be unreliable (Hanski et al. 1996b). In this paper we used a 3-yr data set from a very large patch network, which includes several smaller semiindependent patch networks. One can reasonably assume that this data set is large enough to average out network-specific fluctuations. On the other hand, little can be done to account for possible large-scale and long-term spatially correlated fluctuations that can result from, e.g., climatic perturbations.
Another problem is that there is no sure way of knowing whether all of the semi-independent networks used in parameter estimation were even close to an equilibrium (Hanski et al. 1995b). If some of the networks included in parameter estimation are badly off balance, parameter estimates may be thrown off mark. There are ways to allow for bad data in parameter estimation, one of which is to calculate [Theta] values from only a part of the data, leaving out for example the 5% of the patches with worst individual log-likelihood errors. This was tested with the current data set and very similar parameter values were obtained when omitting 5, 10, or even 25% of the data in parameter estimation. Another method for graphically examining the reliability of model fit is based on repeated parameter estimations using bootstrap replicate data sets. An example of this method is shown in Fig. 5.
Our main conclusion is that we failed in this study to produce compelling evidence for the need to add environmental factors into the incidence function model, even though statistically significant results were obtained for the effects of some factors and hence more was learned about the biology of the species. A probable reason for the rather small improvements in model fit subsequent to the addition of extra environmental factors is that the habitat patches have been delineated carefully using extensive biological information and thus a lot of information on “environmental factors” has already been implicitly taken into account in the area variable of the basic model. Even so, there were indications, such as the strong effect of grazing, which suggest that including environmental factors into the basic incidence function model could well be necessary in other metapopulation settings. The inclusion of detailed landscape structure into the basic model produced very little, not fulfilling the hopes that have been placed on such an approach (Wiens 1997). This failure is due to reasons such as inaccuracies in satellite data, generic problems in the modeling of migration routes and the fact that the landscape across Aland is fairly homogeneous. Finally, we emphasize the constructive facet of these results: the study and modeling of many species at the metapopulation level may be profitably based on the effects of just the two variables in the heart of the classical metapopulation concept, the areas and the isolations of the habitat patches.
ACKNOWLEDGMENTS
We thank Mikko Kuussaari, Juha Poyry, and two anonymous referees for useful comments on the manuscript and Juha Poyry for providing the Aland data set in a conveniently recorded electronic form. All students who participated in the field surveys in 1993-1995 are thanked for their contributions.
LITERATURE CITED
Aberg, J., G. Jansson, J. E. Swenson. and P. Angelstam. 1995. The effect of matrix on the occurrence of hazel grouse (Bonasa bonasa) in isolated habitat fragments. Oecologia 103:265-269.
Andren, H. 1994. Effects of habitat fragmentation on birds and mammals in landscapes with different proportions of suitable habitat: a review. Oikos 71:355-366.
Brown, J. H., and A. Kodric-Brown. 1977. Turnover rates in insular biogeography: effect of immigration on extinction. Ecology 58:445-449.
Cook, R. R., and I. Hanski. 1995. On expected lifetimes of small and large species of birds on islands. American Naturalist 145:307-315.
Durrett, R., and S. A. Levin. 1994. Stochastic spatial models: a user’s guide to ecological applications Philosophical Transactions of the Royal Society, London B 343:329-350.
Efron, B., and R. J. Tibshirani. 1993. An introduction to the bootstrap. Chapman and Hall, New York, New York, USA.
Gilpin, M., and I. Hanski, editors. 1991. Metapopulation dynamics: empirical and theoretical investigations. Academic Press, London, UK.
Gustafson, E. J., and R. H. Gardner. 1996. The effect of landscape heterogeneity on the probability of patch colonization. Ecology 77:94-107.
Gyllenberg, M., and I. Hanski. 1992. Single-species metapopulation dynamics: a structured model. Theoretical population Biology 42:35-61.
Gyllenberg, M., I. Hanski, and A. Hastings, 1997. Structured metapopulation models. Pages 93-122 in I. Hanski and M. Gilpin, editors. Metapopulation biology: ecology, genetics and evolution Academic Press, London, UK.
Hanski, I. 1985. Single-species spatial dynamics may contribute to long-term rarity and commonness Ecology 66: 335-343.
—–. 1992. Inferences from ecological incidence functions. American Naturalist 139:657-662.
—–. 1994a. A practical model of metapopulation dynamics. Journal of Animal Ecology 63:151-162.
—–. 1994b. Spatial scale, patchiness and movement on land. Philosophical Transactions of the Royal Society of London, Series B 343:19-25.
—–. 1997. Metapopulation dynamics: from concepts and observations to predictive models Pages 69-92 in I. Hanski and M. Gilpin, editors Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
Hanski, I., and M. Gilpin, editors. 1997. Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
Hanski, I., and M. Gyllenberg. 1993. Two general metapopulation models and the core-satellite hypothesis. American Naturalist 142:17-41.
Hanski, I., M. Kuussaari, and M. Nieminen. 1994. Metapopulation structure and migration in the butterfly Melitaea cinxia. Ecology 75:747-762.
Hanski, I., A. Moilanen, and M. Gyllenberg. 1996a. Minimum viable metapopulation size. American Naturalist 147: 527-541.
Hanski, I., A. Moilanen, T. Pakkala, and M. Kuussaari. 1996b. Metapopulation persistence of an endangered butterfly: a test of the quantitative incidence function model. Conservation Biology 10:578-590.
Hanski, I., T. Pakkala, M. Kuussaari, and G. Lei. 1995a. Metapopulation persistence of an endangered butterfly in a fragmented landscape. Oikos 72:21-28.
Hanski, I., J. Poyry, T. Pakkala, and M. Kuussaari. 1995b. Multiple equilibria in metapopulation dynamics. Nature 337:618-621.
Hanski, I., and D. Simberloff. 1997. The metapopulation approach, its history, conceptual domain and application to conservation. Pages 5-26 in I. Hanski and M. Gilpin, editors. Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
Harrison, S. 1991. Local extinction in metapopulation context: an empirical evaluation. Pages 73-88 in M. Gilpin and I. Hanski, editors. Metapopulation dynamics: empirical and theoretical investigations. Academic Press, London, UK.
—–. 1994. Metapopulations and conservation. Pages 111-130 in P. J. Edwards, R. M. May, and N. R. Webb, editors. Large scale ecology and conservation biology. Blackwell, Oxford, UK.
Harrison, S., and A. Taylor. 1997. Empirical evidence for metapopulation dynamics. Pages 27-42 in I. Hanski and M. Gilpin, editors. Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
Hastings, A., and S. Harrison. 1994. Metapopulation dynamics and genetics. Annual Review of Ecology and Systematics 25:167-188.
Kuussaari, M., M. Nieminen, and I. Hanski. 1996. An experimental study of migration in the Glanville fritillary butterfly Melitaea cinxia. Journal of Animal Ecology 65:791-801.
Lamberson, R. H., K. S. McKelvey, R. J. Noon, and C. Voss. 1992. A dynamic analysis of northern spotted owl viability in a fragmented forest landscape. Conservation Biology 6:505-512.
Lei, G., and I. Hanski. 1997. Metapopulation structure of Cotesia melitaearum, a parasitoid of the butterfly Melitaea cinxia. Oikos 78:91-100.
Levins, R. 1969. Some demographic and genetic consequences of environmental heterogeneity for biological control. Bulletin of the Entomological Society America 15: 237-240.
MacArthur, R. H., and E. O. Wilson. 1967. The theory of island biogeography. Princeton University Press, Princeton, New York, USA.
Moilanen, A. 1995. Parameterization of a metapopulation model: an empirical comparison of several different genetic algorithms, simulated annealing and TABU search. Pages 551-556 in Proceedings of the Second IEEE International Conference on Evolutionary Computation. IEEE Press, Piscataway, New Jersey, USA.
Moilanen, A., and I. Hanski. 1995. Habitat destruction and competitive coexistence in a spatially realistic metapopulation model. Journal of Animal Ecology 64:141-144.
Moilanen A., A. Smith, and I. Hanski. 1998. Long-term dynamics in a metapopulation of the American pika. American Naturalist, in press.
Murphy, D. D., K. E. Freas, and S. B. Weiss. 1990. An environment-metapopulation approach to population viability analysis of a threatened invertebrate. Conservation Biology 4:41-51.
Nee, S., R. M. May, and M. P. Hassell. 1997. Two-species metapopulation models. Pages 123-148 in I. Hanski and M. Gilpin, editors. Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
Taylor, A. 1991. Studying metapopulation effects in predator-prey systems. Pages 305-323 in M. E. Gilpin and I. Hanski, editors. Metapopulation dynamics. Academic Press, London, UK.
Thomas, C. D., and I. Hanski. 1997. Butterfly metapopulations. Pages 359-386 in I. Hanski and M. Gilpin, editors. Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
Thomas, C. D., and T. M. Jones. 1993. Partial recovery of a skipper butterfly (Hesperia comma) from population refuges: lessons for conservation in fragmented landscapes. Journal of Animal Ecology 62:472-481.
Tilman, D., and P. Kareiva, editors. 1997. Spatial ecology. Princeton University Press, Princeton, New Jersey, USA.
Tilman, D., R. M. May, C. L. Lehman, and M. A. Nowak. 1994. Habitat destruction and the extinction debt. Nature 371:65-66.
Vuorela, A. 1996. Satellite image based land cover and forest classification of Finland. In R. Kuittinen, editor. Finnish-Russian remote sensing research seminar. Reports of the Finnish Geodetic Institute X:96.
Wahlberg, N., A. Moilanen, and I. Hanski. 1996. Predicting the occurrence of species in fragmented landscapes. Science 273:1536-1538.
Weiss, S. B., D. D. Murphy, P. R. Ehrlich, and C. F. Metzler. 1993. Adult emergence phenology in checkerspot butterflies: the effects of macroclimate, topclimate, and population history. Oecologia 96:261-270.
Wiens, J. 1997. Metapopulation dynamics and landscape ecology. Pages 43-62 in I. Hanski and M. Gilpin, editors. Metapopulation biology: ecology, genetics and evolution. Academic Press, London, UK.
COPYRIGHT 1998 Ecological Society of America
COPYRIGHT 2000 Gale Group
