H. Resit Akcakaya

Linking GIS with Models of Ecological Risk Assessment for Endangered Species

A model that links GIS to models for viability analysis and risk assessment is applied to endangered species, including the Spotted Owl in the northwestern US, the Red-cockaded Woodpecker in Louisiana, and the California Gnatcatcher in Orange County, California. The model integrates landscape data on habitat requirements with demographic data to analyze risks of extinction, evaluate management options, and assess human impact on wildlife populations. Other applications of the model involve design of nature reserves, wildlife management, and population viability analysis. The model analyzes habitat data exported from a GIS, and identifies the patches of habitat that can support a population. The structure of these patches, including their locations, sizes and distances from each other, define the spatial structure of the metapopulation. The spatial structure is combined with demographic data and other information on the ecology of the species to complete a metapopulation model, which incorporates age or stage structure and density dependence for each population, spatial correlation and dispersal among populations, environmental and demographic stochasticity and catastrophes. The model performs a risk analysis, and runs multiple simulations, automatically changing parameters to analyze the sensitivity of risks to input data.

INTRODUCTION

One of the main tasks of conservation biologists is to evaluate the viability of endangered and threatened species under different natural conditions, and under alternative options for wildlife management, reserve design, and habitat protection plans. These evaluations usually ask questions about the predicted future abundance, risk of extinction, or chance of recovery of the species; and are addressed by population viability analysis (PVA), which is a systematic examination of interacting factors that place a population or species at risk of extinction (Gilpin and Soule 1986, Shaffer 1990). The factors that a PVA examines may be both natural and anthropogenic in origin, and their analysis often involves mathematical or computer models that predict the future changes in the abundance and distribution of the species in question, given information about its ecology and demography (Burgman et al. 1993).

Habitat loss and fragmentation are among the most common threats facing endangered species, making GIS-based evaluations an essential component of population viability analyses. Often habitat loss and fragmentation, combined with the natural heterogeneity of landscapes, forces species to exist in multiple populations inhabiting relatively isolated habitat patches. Such a collection of populations of the same species is called a metapopulation. The existence of multiple populations usually introduce complexities that make it impossible to evaluate the viability of the species based on PVAs performed on separate populations, and necessitate a metapopulation modeling approach.

This paper describes a computer program for building metapopulation models and performing GIS-based PVAs, and discusses its application to cases involving endangered species.

THE MODEL

The PVA program RAMAS/GIS (Akcakaya 1994a) is designed to link GIS-generated landscape data with a detailed metapopulation model for extinction risk assessment, viability analysis, reserve design and wildlife management. Descriptions of the first version of the program can be found in Akcakaya (1994b) and Akcakaya et al. (1995a); it was also reviewed by Kingston (1995). This section describes the second version of the program, which is being tested with applications to various endangered species (see the section on "Applications" below).

The program operates in four steps. First, landscape data are analyzed and the patch structure of the habitat is exported to a metapopulation model. Second, temporal changes in habitat characteristics are modeled. Third, a metapopulation model is built by combining spatial and demographic information. Fourth, simulations are run to estimate risks of extinction or decline, and to predict the abundance and distribution of individuals in the metapopulation. The essential aspects of these four model components are summarized below.

Landscape data

The function of the Landscape data component is to analyze GIS data to determine the spatial structure of the metapopulation as well as several population-specific demographic parameters that may depend on habitat quality. This component works in the 6 steps described below.

1. Imports landscape data from a GIS

RAMAS/GIS operates under MS-DOS, and can import landscape data in the form of raster maps exported in ASCII format from ARC/INFO, GRASS, IDRISI, and any GIS that supports bare grid ASCII format; it can also import maps in binary format from IDRISI. The requirements are: (i) maps must be in raster format, with the cells arranged in a square (not hexagonal) format; (ii) map values must be numerical; (iii) all maps must describe the landscape with the same precision (i.e., number of cells in both north-south and east-west directions must be the same); (iv) coverage of all maps must be identical (the corners of the rectangular area described by all maps must be the same); (v) maximum map size depends on the available memory, but cannot exceed 5000 rows and 5000 columns in version 2.0.

2. Creates a habitat map

Information in different maps are then combined to make a habitat suitability (HS) map, with a user-defined habitat suitability function. This function must be estimated outside the program. Most methods of estimating the HS function involve statistical procedures, using species occurrence or abundance at each location as the dependent variable and the habitat characteristics as the set of independent variables. The statistical procedures most commonly used are stepwise multiple regression methods (including logistic regression) and stepwise discriminant function analysis. The function is then entered in RAMAS/GIS, which creates a habitat suitability map by calculating the value of the function for each cell of the raster map.

3. Finds habitat patches

The program employs a patch-recognition algorithm to find clusters or groups of nearby cells that have HS values higher than or equal to a threshold habitat value, and labels them as patches. For this procedure, the program uses two parameters that determine how the species perceives (or reacts to) the patchiness of the habitat. These are threshold HS value and neighborhood distance. Threshold HS is the minimum HS value (as defined by the HS function) below which the habitat is not suitable for reproduction and/or survival (although individuals may disperse or migrate through habitat that has a lower HS than this threshold). Neighborhood distance is used to identify nearby cells that belong to the same patch. Suitable cells (as defined by the threshold parameter) that are separated by a distance less than or equal to the neighborhood distance are regarded to be in the same patch. For an animal species, the neighborhood distance parameter may represent the foraging distance.

4. Calculates demographic parameters

Next, the program calculates five demographic parameters for each patch. These are carrying capacity, initial abundance, maximum growth rate, relative fecundity and relative survival. The meaning of these parameters will be described below (see the section "Metapopulation model"). The program allows these five parameters to be calculated as any arbitrary function of patch characteristics, such as total habitat suitability in the patch, average habitat suitability, area as the number of cells in the patch, length of the perimeter (edge) of the patch, average values of input maps for the patch. In addition to these landscape variables, the functions for the five demographic parameters can also include any of the standard mathematical functions (such as min, max, log, exp, etc.)

5. Determines the spatial structure

In addition to the population-specific demographic parameters described above, the program calculates several parameters related to the spatial structure of the metapopulation. These include the location of patches, the distances among patches and migration (dispersal) rates, and spatial correlations based on these distances. These two metapopulation-level parameters, dispersal and correlation, will be discussed below. The rate of dispersal between two population can depend on one of three distance measures: edge to edge (minimum distance between the two patches), center to edge (distance from the center of source patch to the closest edge of the target patch), or center to center (distance between the centers of the two patches).

6. Exports the results to a metapopulation model

Finally, the program exports the patch coordinates, demographic parameters (carrying capacities, growth rates, etc.), the dispersal (migration) rates and the correlation structure to the metapopulation model.

Habitat dynamics

The second component of the program is designed to be used in modeling temporal changes in habitat. It allows the calculation of a time-series of carrying capacities for each population. It reads data files saved from the Landscape data subprogram, and outputs a set of data files for the Metapopulation model subprogram. One of these data files contains the main metapopulation model, and others contain the temporal changes in carrying capacities. The Metapopulation model subprogram inputs these files of carrying capacities and vital rates, and uses them to model habitat dynamics.

Metapopulation model

The main component of the program is where the spatial and demographic parameters calculated by the Landscape data component are combined with other ecological and demographic information about the species to develop a metapopulation model. The model may incorporate the factors and parameters discussed below.

Age structure or Stage structure within populations is modeled by a matrix model (Caswell 1989) that incorporates age- or stage-specific vital rates (survival rates and fecundities). Each population in the model can have a different stage matrix, and a different initial age or stage structure (initial number of individuals in each age or stage). The initial structures can be specified as the stable age or stage distributions. The population-specific stage matrix can be specified to change through time, by reading two files (one each for fecundities and survival rates) that contain the temporal change in the relative values of the vital rates.

Density dependence in population dynamics is modeled by modifying the mean values of survival rates and fecundities as a function of the population size (N). Density-dependent population growth may involve a simple ceiling model, logistic-like functions that describe contest- or scramble-type intraspecific competition (including Ricker and Beverton-Holt functions), Allee effects (i.e., density dependence at low population sizes), or Allee effects combined with density dependence at high population sizes. All density dependence functions are parameterized with the same set of parameters that include maximal growth rate (Rmax) and carrying capacity or ceiling (K), random variation in K, and temporal trend in K . Each population can have a different set of parameters.

Habitat change, e.g., habitat loss (as a result of human impact) or habitat increase (as a result of vegetation growth) can be modeled by specifying how the carrying capacity of each population changes through time. This can be done either with a constant rate, or as a time-series of carrying capacities saved in a disk file.

Environmental stochasticity is modeled by random fluctuations in vital rates and in carrying capacities. The random fluctuations can be normal- or lognormal-distributed, and can be correlated among populations. They are assumed to be perfectly correlated among age classes or stages within each population.

Demographic stochasticity is modeled by sampling the number of survivors from a binomial and the number of offspring from a Poisson distribution (Akcakaya 1991).

Catastrophes are rare events that can affect abundances (a proportion of all individuals die), vital rates (survival rates and fecundities are reduced after a catastrophe), or carrying capacities (which are reduced after a catastrophe). The spatial extent of catastrophes may be local or regional. The impact of catastrophes can be population-specific (some populations may be more prone to, or more affected by catastrophes than other populations), or they may be stage-specific (some stages, or even certain vital rates may be more affected by catastrophes than others).

Density dependence in migration is modeled by making the total rate of emigration or dispersal a function of the size of the population. The rate of dispersal/emigration can be specified either to increase or decrease as N increases. Migration rates can also be stage-specific.

Geographic configuration is specified by the coordinates of each population. The distance of populations from each other and their relative positions are utilized in modeling the effects of the two spatial factors discussed below.

Dispersal (migration) describes the movement of individuals among populations. In RAMAS/GIS, migration is modeled by specifying the proportion of individuals that move from each population to each other at every time step. These rates are input in the form of a migration matrix. In most cases, the rate of dispersal may be a function of the distance between source and target populations. RAMAS/GIS allows users to specify a function that describes the dependence of dispersal rates on distance. The matrix can filled according this to function, and can be edited to account for habitat corridors (by increasing the rate between specific pairs of populations) and for obstacles or geographic barriers to migration (by decreasing the rate). The migration rates may also be specified to be dependent on population size to allow for density-dependent migration (see above), or on the age or stage of the individuals to allow for age- or stage-specific dispersal tendencies.

Correlations among populations describe the similarity of environmental patterns experienced by each population. This factor is important in the "rescue effect" in metapopulations: when fluctuations are spread over a number of separate populations, the overall risk faced by the metapopulation is reduced. If the fluctuations in the environment are at least partially independent (uncorrelated), it will be less likely that all populations go extinct at the same time than if the fluctuations are dependent (i.e., synchronous or correlated). In uncorrelated environments, extinct populations will have a chance to be recolonized. Thus, correlation of environmental fluctuations has important effects on metapopulation persistence and viability (Gilpin 1988; Akcakaya and Ginzburg 1991; Burgman et al. 1993).

RAMAS/GIS models correlations by sampling the vital rates of each population from a normal or lognormal distribution which is correlated with the vital rates of other populations according to a correlation matrix specified by the user. This matrix can be specified as a function of the distance between populations (as closer populations are more likely to experience similar environmental patterns).

Risk assessment

In this component, metapopulation models built in RAMAS are used to simulate the dynamics and to predict the future of the metapopulation. RAMAS/GIS summarizes the results of a simulation with 16 types of output, some of which are superimposed. Most of these results are related to risk assessment and report risk analytical measures such as risk of extinction and time to extinction (Akcakaya 1992). Types of output include

* Risk of extinction (or decline to any level) any time during, or at the end of, the simulated time period;
* Probability of exceeding a range of population sizes any time during, or at the end of, the simulated time period;
* Distribution of times to extinction (or decline to a specific level), and cumulative time to extinction or decline;
* Distribution of times that the metapopulation will exceed a specified threshold, and cumulative time to increase;
* Abundance (of the metapopulation and each of its populations) through time;
* Metapopulation occupancy (number of extant populations through time);
* Local occupancy (for each population, the number of time steps that it was extant);
* Local extinction duration (for each population, the maximum number of consecutive time steps during which a population remained unoccupied during a typical replication);
* Final stage structure (for each population, histogram of the number of individuals in each stage at the end of the simulation);
* Metapopulation structure (histogram of the number of individuals in each population at a specific time step).

The model can be run several times, to analyze the sensitivity of results to input parameters by varying them automatically, to compare management options, or to assess anthropogenic impact by comparing outputs from simulations with parameters for impacted and non-impacted situations. The sensitivity analysis facility also allows the user to superimpose graphs from different simulations to make the comparisons easier.

APPLICATIONS

The second version of the program is currently being tested with applications to four endangered birds (California gnatcatcher, cactus wren, red-cockaded woodpecker, and northern spotted owl). These applications are briefly described below.

California gnatcatcher and cactus wren

The program is applied to California gnatcatcher Polioptila californica and cactus wren Campylorhynchus brunneicapillus metapopulations in central and coastal regions of Orange County, California. The California gnatcatcher (a federally-listed threatened species) and cactus wren are two of the three "target species" in the State of California's Natural Community Conservation Planning (NCCP) program. The target species are considered as surrogates for the conservation of coastal sage scrub, a plant community that has decreased substantially from its historic coverage as a result of urban and agricultural development, and has become fragmented.

The application started with a statistical analyses of observation locations, and GIS maps for elevation, slope, and vegetation (exported from ARC/Info) to characterize the habitat of the two birds. The habitat maps were then analyzed to determine the spatial structure of the metapopulation model. These data were combined with demographic parameters estimated from field work conducted by Atwood et al. (1995) to produce spatially-explicit, stage-structured, stochastic metapopulation models for the two species (Akcakaya and Atwood, in prep.). Simulations and sensitivity analyses with these models pointed out to the relative importance of temporal variation in vital rates (including the frequency and the magnitude of catastrophic declines in vital rates as was observed in the winter of 1994-95).

Red-cockaded woodpecker

The goal of this application of the program was to evaluate the impact of timber prescriptions (forest management practices) on the viability of a red-cockaded woodpecker Picoides borealis metapopulation in Louisiana. The habitat of the red-cockaded woodpecker was characterized based on the stand type (i.e, the dominant tree species), stand condition, and basal area of pine species. These variables were input into the program in the form of GIS maps exported from GRASS (Akcakaya et al. 1995b). Currently the habitat characterization is being extended to include other habitat variables such as distance to streams, and age of dominant trees in the stand. The habitat descriptions are used to evaluate the impact of forest management practices, such as a planned cut, on the carrying capacities of habitat patches and other parameters of the red-cockaded woodpecker populations inhabiting these patches. The model will then be run with different management options to compare their impacts in terms of the predicted risk of extinction, or rate of recovery of the metapopulation.

Northern spotted owl

This application focused on factors affecting the viability of the Northern Spotted Owl Strix occidentalis caurina throughout its range in the U.S. A second goal in this analysis was to incorporate two sources of variability in determining the threat the species faces. The study incorporated natural variation (resulting from temporal fluctuations in environmental factors) in the form of randomly distributed vital rates (survivals and fecundities). In addition, demographic stochasticity was modeled to describe chance variations in reproduction, survival and dispersal. These types of natural variation (environmental and demographic) were used to express the model results in probabilistic terms such as the viability of the species (for example in terms of the chance of survival or risk of extinction). Uncertainties that result from a lack of knowledge were incorporated in the form of parameter ranges, and were used to estimate upper and lower bounds on the estimated viability of the species.

Based on the habitat maps provided by the Forest Service, the program found 42 habitat patches. The size distribution of the patches was very skewed, with the 4 largest patches making up about 94% of the total area of all patches, and the seven largest making up about 96%. Because of the large differences in sizes of neighboring populations, the model results (risk of decline) were not very sensitive to the rate of inter-patch dispersal of juvenile spotted owls. The model predicted a large difference between upper and lower bounds on the viability of the northern spotted owl, based on the best-case and worst-case scenarios, which were parameter combinations that resulted in best and worst chance for survival. According to sensitivity analyses, the viability of the metapopulation was most sensitive to the set of vital rates used (the dependence of fecundities and survival rates on habitat), and also sensitive to the degree of spatial correlation among vital rates of the populations, and to the carrying capacities of the populations. In addition, metapopulation occupancy was sensitive to dispersal and Allee effects (Akcakaya, unpublished).

Other applications

The methodology used in the models, and the relative effects of factors discussed above on metapopulation dynamics have been discussed in several previous publications (Akcakaya 1991, 1992; Akcakaya and Ginzburg 1991; Burgman et al. 1993). Previous versions and precursors of the program have been applied to a wide variety of organisms. RAMAS/GIS (version 1.0) was used in assessing the effectiveness of translocating individuals as a management option for the endangered helmeted honeyeater (Lichenostomus melanops cassidix) in Australia (Akcakaya et al. 1995a).

The metapopulation model component was used to build a model for a land snail Arianta arbustorum metapopulation in northeastern Switzerland. The model was used to investigate the effect of population subdivision on the persistence of a land snail metapopulation and to analyze the interaction between spatial factors, population subdivision, and catastrophes (Akcakaya and Baur, in press).

A precursor to the metapopulation model component (RAMAS/space; Akcakaya and Ferson 1990) has been used to model metapopulation dynamics the California spotted owl (LaHaye et al. 1994), a subspecies of the northern spotted owl. Gibbs (1993) used RAMAS/space to model dynamics of spotted salamanders, bullfrogs, snapping turtles, swamp sparrows and water shrews in a mosaic of wetlands.

REFERENCES

Akcakaya, H.R. (1991) A method for simulating demographic stochasticity. Ecological Modelling 54:133-136.

Akcakaya, H.R. (1992) Population viability analysis and risk assessment. Pages 148-157 in Wildlife 2001: Populations. D.R. McCullough and R.H. Barrett, eds. Elsevier Publishers, London.

Akcakaya, H.R. (1994a) RAMAS/GIS: Linking Landscape Data With Population Viability Analysis (version 1.0). Applied Biomathematics, Setauket, New York.

Akcakaya, H.R. (1994b) GIS enhances endangered species conservation efforts. GIS WORLD Vol.7, November 1994, pp. 36-40.

Akcakaya, H.R., and Baur, B. (in press). Effects of population subdivision and catastrophes on the persistence of a land snail metapopulation. Oecologia .

Akcakaya, H.R., and Ginzburg, L.R. (1991) Ecological risk analysis for single and multiple populations. Pages 73-87 in: Species Conservation: A Population-Biological Approach. A. Seitz and V. Loeschcke (eds.) Birkhaeuser Verlag, Basel.

Akcakaya, H.R., and Ferson, S. (1990) RAMAS/space User Manual: Spatially Structured Population Models for Conservation Biology. Applied Biomathematics, New York.

Akcakaya, H.R., McCarthy, M.A., and Pearce, J. (1995a) Linking landscape data with population viability analysis: management options for the helmeted honeyeater Biological Conservation 73:169-176.

Akcakaya, H.R., Jackson, J., and Ginzburg, L.R. (1995b) Modeling the Red-cockaded Woodpecker Populations at Fort Polk. Report to U.S. Army Corps of Engineers.

Atwood, J.L., Fugagli, M.R., Reynolds, C.H., Luttrell, J.C., and Tsai, S.H. (1995) Distribution, population dynamics, and dispersal of California Gnatcatchers on the Palos Verdes Peninsula, 1993-1995. Proceedings of CalGnat'95 Symposium (in press).

Burgman, M.A., Ferson, S., and Akcakaya, H.R. (1993) Risk Assessment in Conservation Biology. Chapman and Hall, Population and Community Biology Series.

Caswell, H. (1989) Matrix Population Models: Construction, Analysis, and Interpretation. Sinauer Associates, Sunderland, Massachusetts.

Gibbs, J.P. (1993) Importance of small wetlands for the persistence of local populations of wetland associated animals. Wetlands 13:25-31.

Gilpin, M.E. (1988) A comment on Quinn and Hastings: extinction in subdivided habitats. Conservation Biology 2:290-292.

Gilpin, M.E., and Soule, M.E. (1986) Minimum viable populations: processes of species extinction. Pages 19-34 in Conservation Biology: the science of scarcity and diversity, M.E. Soule, ed. Sinauer Associates, Sunderland, Massachusets.

Kingston, T. (1995) Valuable modeling tool: RAMAS/GIS. Conservation Biology 9:966-968.

LaHaye, W.S., Gutierrez, R.J., and Akcakaya, H.R. (1994) Spotted owl metapopulation dynamics in southern California. Journal of Animal Ecology 63:775-785.

Shaffer, M.L. (1990) Population viability analysis. Conservation Biology 4: 39-40.

H. Resit Akcakaya
Applied Biomathematics
100 North Country Road
Setauket, New York 11733, USA
Tel: 516-751-4350
Fax: 516-751-3435
email: risk@life.bio.sunysb.edu