Viability analysis for endangered metapopulations: a diffusion approximation approach

E. E. Holmes1
National Marine Fisheries Service
Northwest Fisheries Science Center
2725 Montlake Blvd E
Seattle, WA 98112

B. Semmens
Zoology Department
University of Washinton
Seattle, WA 98195

1 Corresponding author
e-mail: eli.holmes@noaa.gov
phn: 206-860-3369
fax: 206-860-3467

23.1.  Introduction
Population viability analysis (PVA) assesses the rate of population decline and the risks of extinction or quasi-extinction over a defined time horizon for a population of concern (Gilpin and Soule 1986; Boyce 1992; Morris and Doak 2002).  Although the techniques employed to conduct PVA are varied, they typically involve building quantitative models that are parameterized by demographic and environmental data.  PVA was first used in the early 1980s (Shaffer 1981), and in the past decade it has gained broad acceptance in the conservation community as a useful tool for assessing and managing ‘at risk’ species (Beissinger 2002; Morris and Doak 2002; Reed et al. 2002).  This is particularly true for demogaphic PVAs, due in large part to the advancements in monte carlo techniques and desktop computers (Beissinger 2002).  The International Union for the Conservation of Nature (IUCN)’s Red List Criteria, probably the most widely applied set of decision rules for determining the status of at risk species, is partially defined by metrics that require some form of PVA (IUCN 1994).  For instance under one of the Red List criteria, a taxon may be classified as endangered if a “reduction of at least 50%, projected or suspected to be met within the next ten years or three generations” is predicted.
Although many PVAs are focused on single populations in single sites, there are often needs for spatially explicit PVAs: many populations of conservation concern are distributed across multiple sites and additionally, the primary anthropogenic threats facing at risk species are habitat destruction and alteration, which are fundamentally spatial processes (Wilcove et al. 1998).  Several software packages have been written for spatially explicit PVA, including RAMAS® Metapop (Akçakaya 1997) and RAMAS® GIS (Boyce 1996), ALEX (Possingham and Davies 1995), and VORTEX (Lacy 1993).  These models incorporate a diversity of demographic and spatial attributes such as distance-dependent migration, allee effects, social population structure, habitat quality and spatial arrangement, and genetic variability.  The development of flexible sophisticated PVA software packages such as these has made the construction and simulation of spatially explicit PVA models feasible for those who are not highly skilled programmers and has greatly increased the number of managers and scientists capable of using spatially realistic PVA models.  
As the use of PVA has grown in conservation science, so have concerns that PVAs are often overextended given limited data sets (Reed et al. 2002).  Beissinger and Westpahl (1998) suggested that PVA should be limited to assessing relative risks over short time frames using the simplest models that can reasonably be justified.  For single species with spatially simple structure, data needs can often be met when Beissinger and Westpahl’s call for model moderation and simplicity are heeded.  When one is faced with species with more complex spatial structure, a much larger amount of data is needed to parameterize the dynamics of individual local populations, the levels and patterns of dispersal, and the spatial pattern of temporal correlations among local populations (e.g. Ralls et al. 2002).  Unfortunately, collection of the data needed to parameterize a spatial model is rare for species of conservation concern, at least in the United States (Morris et al. 2002), and there is a disconnect between the parameter requirements for spatially explicit PVA models, and the willingness and/or ability of management agencies to collect the types of data needed to appropriately apply such tools.  Since it is usually impossible to retroactively fulfill data requirements for a spatial PVA and there will always be cases where collection of spatial data is infeasible, managers require PVA tools that that can help guide conservation of metapopulations in the absence of large amounts of spatial data.

Diffusion approximation for metapopulations
One approach to the problem of limited population data is to find a diffusion approximation that correctly models the long-run statistical properties of a complex population process.  This approach has been used successfully for single population models (Karlin and Taylor 1981; Lande and Orzack 1988; Lande 1993; Dennis et al. 1991; Hill et al. 2002; see also Morris and Doak 2002, Ch. 3, and Lande et al. 2003) and reduces the problem of parameterizing a large model with many parameters to the much simpler task of parameterizing a two-parameter diffusion model.  One of the main practical implications of the diffusion approximation approach is that it is not necessary to know the multitude of parameters describing the local dynamics, dispersal levels, spatial patterns of dispersal, and spatial synchrony between local populations in order to make basic predictions about the statistical distribution of the long-term metapopulation or local population trajectories.  The relevant two parameters for the diffusion approximation can be estimated from a simple time series of counts from the population process.  
In this chapter the diffusion approximation approach is used to model the long-run behavior of spatially structured populations.  Our focus is on stochastic metapopulations characterized by structured population size, density-independent local dynamics, and, in keeping with the assumption of density independence, a metapopulation that is declining as a whole.  Local populations are assumed to have patch-specific structured local dynamics and dispersal rates, with spatial structure among local populations both in terms of their local dynamics and dispersal patterns.  Description of the long-run statistical behavior of the metapopulation trajectories using a diffusion approximation allows the estimation of PVA risk metrics such as the long-term rate of metapopulation decline and the probability of reaching different threshold declines over different time horizons (i.e. probabilities of extinction or quasi-extinction).  These methods for estimating metapopulation PVA metrics are illustrated using data from two chinook salmon metapopulations in the U.S. Pacific Northwest.  

23.2.  A stochastic metapopulation model
Our focus is on declining metapopulations, and thus what has been termed non-equilibrium metapopulations.  We model a collection of local populations connected by dispersal where local populations have density-independent local dynamics, which may be “sources” or “sinks”, but the metapopulation as a whole is declining.  Dispersal levels could very low, resulting in basically independent local populations, or extremely high, resulting in essentially one population.  From a practical standpoint, this approach is most appropriate when dispersal is not insignificant (e.g. above 2% per year localized dispersal or 0.1% global dispersal), otherwise parameterization of the model requires inordinately long time series.  Data from this type of metapopulation would be characterized by fluctuating local population trajectories, but actual extinctions would be unusual until the metapopulation has very few individuals.  Our model assumes no density-dependence nor carrying capacities within the individual local populations.  Such a model is only appropriate in cases where the population is declining and all local populations are well below their carrying capacities.  Our example using data on chinook salmon illustrates a situation that is likely to be well modeled as this type of metapopulation.  
In the following section, a rather parameter intensive mathematical description of a stochastic, declining metapopulation is given.  However, the reader should keep in mind that this model will not be parameterized.  Rather the asymptotic behavior of this model’s trajectories will be derived and that information will be used to develop a diffusion approximation of the process.  Time series data will then used to parameterize the diffusion approximation.

The model
Consider an individual local population i with stochastic yearly growth and stochastic dispersal to and from other local populations.  Such a local population’s numbers in year t, Ni(t), could be described as follows:

    [1]

where zi(t) is the stochastic growth rate of local population i in year t and is a random variable with some unspecified statistical distribution with mean mi and variance .  The mi term will be referred to as the local population’s intrinsic growth rate; it will not be observed since the local population is subject to immigration and emigration.  Some fraction of individuals, di(t), leaves local population i at year t and disperses to other local populations, and dispersal into local population i occurs from other local populations.  The fraction of dispersers from local population j that go to local population i in year t is aji(t) and can vary depending on the destination, i, thus allowing for spatially-structured dispersal.  The dispersal parameters, di(t), and aji(t), are assumed to be temporally random variables from some unspecified statistical distribution.

The model in matrix form
The model for the entire metapopulation can be written using a random transition matrix, A(t), which encapsulates both the dispersal and local growth: 

                                                 [2]

where

                            [3]

The ‘(t)’ on the d’s, a’s, and z’s have been left off to remove clutter.  There may be any level or spatial pattern of temporal correlation between the intrinsic local growth rates, zi’s, dispersal rates, di’s, and dispersal patterns, aji’s.
In the matrix model, each row represents 1 unit of habitat.  Local populations with multiple units of habitat appear as multiple rows with very high dispersal between the units of habitat in that local population.  The habitat units within a local population could vary in quality (i.e. habitat within a local population need not be uniform), and different local populations certainly differ in the number of habitat units they contain.  The di’s and aji’s are assumed to be drawn from some distribution that can be different for each local population or local population pair.  Although the di’s, aji, and zi’s are temporally random variables, they are assumed to be stationary, i.e. that there is no overall change in the mean values over time.  For the purposes of this paper, it will be assumed that the di’s, aji’s, and zi’s are all strictly postitive, which this means that all local populations are connected to each other to some (although possibly very low) degree and that mean yearly geometric growth rates, exp(mi)’s, while possibly very small are not zero.  These assumptions imply that the A(t) describe an ergodic set of matrices (Caswell 2001, page 375).  The assumption of strict positivity is not strictly necessary.  It is possible for A(t) to describe an ergodic set if some elements of A are zero; it depends on the pattern of zeros within A (cf. Caswell 2001 for a discussion of the conditions under which matrices are ergodic).
The model is very general, allowing some sites to be dispersal sources and others to be dispersal targets, allowing any spatial pattern of dispersal or spatially correlated local growth rates, allowing any pattern of temporal correlation amongst local growth rates and allowing any combination or pattern of habitat sizes of local sites.

Using random theory to understand the model’s statistical behavior
Together, Eqns. 2 and 3 describe a quite generic model of a declining metapopulation with density-independent local dynamics.  From a viability analysis perspective, one might ask the question: “Can one predict the viability of the total metapopulation?”  In more precise terms, this is asking, what are the statistical properties of the metapopulation trajectories of this type of connected collection of local populations (of the form in Eqns. 2 and 3)?  Clearly, the matrix A(t) has a large number of parameters that would be difficult, if not impossible, to estimate for any given metapopulation of conservation concern.  However, using random theory, it can be shown that the long-term dynamics can be described by only two parameters and that it is unnecessary to know the multitude of other parameters for the purpose of projecting long-run dynamics.
To use this random theory, we first need to recognize that this stochastic metapopulation model falls into the class of random processes that involve products of ergodic random matrices, in this case products of A(t), which can be seen by using Eqn. 2 to project the vector of local population sizes forward:

                                                       [4]

where N(t) is the column vector of Ni’s at time t, in Eqn. 2.  Products of random ergodic matrices have a well-established theoretical foundation and have certain well-studied asymptotic statistical properties.  A brief review of two of the key results from this theory is provided in Box 23.1 and a simulated example is shown to illustrate these results.  As described in Box 23.1, the theory demonstrates that this stochastic, density-independent, metapopulation will have an asymptotic growth rate and that the metapopulation, , the individual Ni(t)’s and sets of Ni(t)’s representing the units of habitat comprising a semi-independent local population will be lognormally distributed with the same parameters:

                                    [5]

Figure 23.1 shows an example of this behavior.  A metapopulation is simulated (described in Box 23.1) and over time, the metapopulation declines at a constant rate and all Ni(t)’s have the same long-term fate.  When viewed over short time frames, t small in the figure, the local sites show different growth rates with some declining more or less than the long-term rate, but over the long-term their rate of decline is the same.
 The model studied here approximates the local dynamics by a simple exponential growth (or decline) model.  However, it has been shown that the results from random theory (presented in Box 23.1) also apply to a more complicated metapopulation model where local dynamics are described by stochastic age-structured Leslie matrices (Sanz and Bravo de la Parra 1998, pgs. 286-294).  Essentially, this occurs because even when the local dynamics are described by a local matrix model, the system can still be described by products of random matrices.

23.3.  Diffusion approximation
The asymptotic distribution of log M(t) in Eqn. 5 has the same properties as the distribution of a diffusion process with drift; it is normal and the mean and variance of the distribution of log M(t) increase linearly with time, t.  This observation in the context of age-structured matrix population models (Lande and Orzack 1988; Dennis et al. 1991) led to the use of a diffusion approximation to enable parameterization using simple time series and to enable calculation of extinction probabilities.  Diffusion approximation methods for single population populations are an important and established method for approximating stochastic trajectories (Lande and Orzack 1988; Dennis et al. 1991; chapter 3 in Morris and Doak 2002; chapter 5 in Lande et al. 2003).  The models for single populations are mathematically analogous to the models used here for metapopulations with a stochastic process involving products of random matrices.  However, in single population models, the matrix represents a life history matrix rather than a growth and dispersal matrix, and the N(t) vector (in Eqn. 2) represents different age or stage classes while in the metapopulation matrix, it represents different local sites and populations.  
	A diffusion approximation with drift is a stochastic process with the following properties (cf. Karlin and Taylor 1981, chap. 7):

                                    [6]

For any non-overlapping pair of time periods, t1 < t2 and t3 < t4, X(t2)-X(t1) and X(t4)-X(t3) are independent random variables.  X(t+t) is a random variable with distribution .  Correspondingly the probability density function for X(t+t) given log X(t) is

             [7]

Behavior of metapopulation trajectories versus diffusion trajectories
The diffusion approximation is based on the behavior of log M(t) as t goes to infinity, however in PVA settings the time frame of interest is substantially less than infinity and typically in the range of 25 to 100 years.  How well does the diffusion approximation do over these finite time periods?  To explore this, a collection of 50 local populations were simulated that were connected by global dispersal ranging from 0.1 to 5% per year and that had correlated local dynamics, zi(t), drawn from a Normal(mean = -0.05, variance = 0.09) and a temporal covariance of (0.2)(0.09) between the zi(t)’s of local populations in any given year.  
If the log metapopulation trajectories behave like a diffusion process, then if we repeatedly generate a large sample of replicate metapopulation trajectories, the mean and variance of from those trajectories, should be a constants over the time period of interest.  Additionally,  should be normally distributed.  To examine whether the metapopulation trajectories had these properties, the simulations were started from a distribution of local population sizes selected from the equilibrium set of local population distributions and then run forward for 200 years.  This was repeated (using the same initial distribution of local populations) 1000 times to estimate the distribution of .  This process was repeated for four randomly chosen initial distributions of local population sizes.  The mean and variance of  is denoted as , respectively, in the figure and discussion below.
Figure 23.2 illustrates the results.  For dispersal levels 1% or higher, the trajectories behaved like a diffusion process with  roughly constant and the distributions approximately normal according to a Kolmogorov-Smirnov test at P = 0.05.  For low dispersal, 0.1%, the trajectories did not behave like a diffusion process for t less than 200 at least.  Neither  nor  were constant except for t > 150 and the normality assumption was generally violated except again at large t.  This means that when dispersal is very low, diffusion approximations for this metapopulation would more approximate than for metapopulations with higher dispersal.
	Figure 23.2 illustrates results from one particular model.  Repeating this process for a number of different models indicated some general behaviors.  The higher the dispersal levels, the more trajectories behaved like a diffusion process.  Global dispersal levels of at least 2 to 5% were generally high enough to result in diffusion-like behavior within a short time frame.  Note that localized dispersal has the effect of lowering the effective dispersal rates.  The higher the amount of temporal covariance between local populations in terms of their yearly growth rates, the more the trajectories behaved like a diffusion process.  The simulations were done with the local population sizes within the equilibrium set of local population distributions – indeed the theory is predicated on the local populations being near equilibrium.  For metapopulations with 2 to 5% dispersal, the local populations equilibrated fairly quickly starting from all local populations with equal numbers.  However, at very low dispersal, equilibration took thousands of time steps.  This suggests that the assumption of equilibrium should be viewed cautiously for metapopulations that have very low dispersal rates between local populations.

23.4.  Estimating the parameters
Maximum likelihood estimates of  and  can be calculated using the diffusion approximation for log M(t).  Denote the observed time series as M = M(0), M(1), M(2),…, M(n).  If we approximate log M(t) as a diffusion process, the likelihood function  is given by the product of the probability function distributions for the transitions from log M(t+1) to log M(t), which is Eqn. 7 with t = 1, over t = 0, 1, 2, …, n.  Thus the log likelihood function is:

                                [8]

The maximum likelihood estimates are obtained by solving for the  and  which maximize Eqn. 8:

	                                         [9]

Note that the unbiased estimator for  would use (n-1) rather than n.  The  and  are analogous to the estimates of mean and variance from n samples from a normal distribution, and confidence intervals on mm and  are analogous:

	                                     [10]

where ta,q is the critical value of a t distribution at P = a and q degrees of freedom and  is the critical value of a chi-square distribution at P = a and q degrees of freedom.  See Dennis et al. (1991), page 120, for a more in-depth discussion of maximum likelihood estimates for diffusion processes.  Following Dennis et al.’s monograph, parameter estimation based on the diffusion approximation has been widely used for the analysis of single population trajectories.  For a discussion of parameter estimation that is not based on the diffusion approximation, the reader is referred to Heyde and Cohen (1985).  
The maximum likelihood estimates assume that the metapopulation has reached a stochastic equilibrium and thus that the diffusion approximation is reasonable.  When exploring these methods using simulations, it’s important to allow the system to equilibrate, after starting the simulation with something peculiar like all local populations at the same size.  Equilibruim can be monitored by waiting for the variance of (log(N(t)i)-log[mean(Ni(t))]) to stabilize.  In simulations done for this paper, the distribution stabilized relatively quickly when dispersal was non-zero.  If dispersal is zero, however, the distribution never stabilizes and the variance of (log(Ni(t))-log[mean(Ni(t))]) continually increases.  For an actual metapopulation, for which one wants to conduct a PVA, it is also critical to test the appropriateness of the diffusion approximation for one’s time series data.  Dennis et al. (1991), page 121-122, and Morris and Doak (2002), pages 73-79, review how to do this which is based on diagnostic procedures for evaluating the appropriateness of linear models.

Parameter bias
The estimators are unbiased maximum likelihood estimators for the diffusion approximation, X(t).  It is important to understand whether and how these estimates are biased when working with short time series of metapopulation trajectories, M(t), as opposed to an actual diffusion process.  In particular,  is certain to be biased to some degree since it relies on the diffusion approximation holding for t = 1 in log M(t+t)/M(t), regardless of the length of the time series used for estimation.  This is not the case for, , which is also an unbiased predictor for M(t) given a long time series (Heyde and Cohen 1985).
To numerically explore parameter bias from short time series, simulations were used to look at the difference between  and  from a 20-year time series versus their true values mm and .  An example metapopulation of 50 local sites was simulated with global dispersal and correlated local growth rates, zi(t), drawn yearly from a normal distribution with mean = mi, variance = f and covariance of 0.2*f between any two local growth rates  Two versions of the simulation were run; one to model uniform site quality (spatially uniform mi = -0.05) and one to model highly variable site quality (spatially variable mi’s).  To explore biases over a range of different dispersal and variability levels, the models were run with dispersal between 0.1 and 5% per year and local variability, f, between 0.1 and 0.5.  These parameters translated to metapopulation-level rates, mm, in the range of 0.01 to –0.05 and metapopulation-level variability, , in the range of 0.001 to 0.08.  For each dispersal and local variability pair, 1000 replicate metapopulation trajectories were simulated, each with an initial distribution of local population sizes selected randomly from the equilibrium set.
The mean difference between  and  over the dispersal and local variability parameter space was very low, < 0.0015, for both the uniform and variable mi simulations.  Overall the lack of bias in  supports metrics that rely primarily on this parameter, such as the metapopulation l (next section).  For most of the parameter space explored, , representing a 0 to 20% under- or overestimation of .  Larger biases, , representing a >20% under- or overestimation, were seen for some parameter combinations.  The impact of this bias depended on where  was used.  For instance, the effect on the estimated confidence intervals on  (Eqn. 10) was minimal with the width of the interval changing by a median 0.002.  The effect on estimated passage probabilities was higher although not dramatic.  For example, the estimated probability that the metapopulation will be ≤10% of current levels at the end of 50 years was decreased by 0 to 0.04 (on a scale from 0 to 1) for the uniform mi simulation and increased by 0 to 0.04 for the variable mi simulation.  The estimated probability that the metapopulation will pass below 10% of current levels at any point during the next 50 years was changed by 0 to 0.09.  Overall, the effect of  bias was low in these simulations, but this will depend on the particular metapopulation and will need to be investigated for individual cases of interest.
In practical applications, one must contend with other factors which can lead to parameter bias, but which are outside the scope of this paper.  In particular observation error, non-equilibrium local population distributions and temporal autocorrelation can lead to parameter bias.  Such problems are being studied in the context single population PVA.  Much of this work is likely to be relevant for metapopulation PVA.  See Morris and Doak (2003), chapters 3-5, for a review and discussion of current work in this area.  

23.5.  Metapopulation viability metrics
One of the most basic viability metrics is the long-term geometric rate of decline (or growth) of a population, termed generally l in the PVA literature.  If l is less than 1.0, the population ultimately declines to extinction and 100*(1-l) is roughly the average yearly percent decline.  The metapopulation l is exp(mm) and its estimate is then

.                                                                  [11]

This definition of l follows Caswell’s use of the symbol ls as the long-term average stochastic growth rate: ls = (N(t)/N(0))1/t as  (Caswell 2001, Eqn. 14.44).  This is the long-run geometric growth rate that would be observed in almost every trajectory.  Defined this way, if l < 1, the population goes extinct with certainty, eventually.  This differs from Dennis et al.’s use of the symbol l where l is used for  and the long-term average geometric growth rate is instead denoted by a = exp(m).  The maximum likelihood estimate of l is a biased estimator; since  is normally distributed, the median value of  is  but the mean value is not.  Dennis et al. (1991, Eqn. 79) give an unbiased estimator (mean() = l) based on Shimizu and Iwase (1981), although Dennis et al. found negligible differences between the biased and unbiased estimators in their examples.
From the asymptotic distribution of log M(t), Eqn. 5, the probability that the metapopulation is below a threshold b at the end of y years can be calculated:

                                    [12]

Although this uses the asymptotic distribution, this is mitigated by the fact that is used for the distribution at the end of y years but not at any time before that.  The estimate of  replaces  and  by their estimates  and .  Like the estimate of l, the median estimate of  is equal to the true value, but not the mean.
Some metapopulations can have a low long-term risk of being below a threshold due to a l near 1.0, but high short-term risks of hitting that threshold due to high variability.  Such quasi-extinction or extinction probabilities are commonly used and very important PVA metrics.  The diffusion approximation for log M(t) can be used to estimate these probabilities for the metapopulation.  The probability of that the diffusion process, X(t), experiences a decline below a threshold log b at some time T less than y years is calculated by integrating over the probability density function for first passage times for a diffusion process with drift (Karlin and Taylor 1981, Theorem 5.3).  Lande and Orzack (1988), page 7419, goes through the calculation, which leads to:

                       [13]

F( ) is the cumulative distribution function for a standard Normal(mean = 0, variance =1).  The estimate of  for the metapopulation uses  and  with X(0) = log M(0).  The estimated probability of extinction (to 1 individual) is calculated using Eqn. 13 and setting b equal to 1.  The reader is cautioned that estimates of extinction are problematic and estimates of quasi-extinction (e.g. some threshold greater than 1 individual) are more robust (cf. Morris and Doak 2002, pages 43-44).  Also Eqn. 13 uses the diffusion approximation over short time scales since it calculates the probability of hitting a threshold at any time, including short times, before y years.  This makes Eqn. 13 more approximate than other metrics.
Other viability metrics based on the diffusion approximation such as the mean time to extinction and median time to extinction are discussed in Lande and Orzack (1988) and Dennis et al. (1991).  Karlin and Taylor (1981) 

Risk metric uncertainty
The 100(1-a)% confidence intervals are often used as characterizations of uncertainty.  These can be calculated for the risk metrics using the estimated distributions of  and .  The confidence intervals for  are

.                                [14]

where ta,q is the critical value of a t distribution at P = a and q degrees of freedom.  The corresponding significance level, a, for a hypothesis test, such as “Is l < b” is the a such that:

.                                                     [15]

Confidence intervals on  and  can be calculated by parametric bootstrapping from the estimated distributions of  and : Normal() and Gamma(shape = (n-1)/2, scale = ).  A large number of (,) pairs are randomly generated by sampling from these distributions and the risk metric Y is calculated (Eqn. 13 or 12) for each pair.  The range of Y over the (,) bootstrapped pairs, for which both parameters are within their respective 100(1-a)% confidence intervals, defines the 100(1-a)% confidence interval for Y.  This and other methods for calculating confidence intervals for diffusion approximation risk metrics are discussed in Dennis et al. (1991), page 127-128.
An alternate way to present the level of uncertainty is to estimate the data support for different values of a risk metric.  There are both frequentist and Bayesian approaches for this (see Wade 2001 for a review geared towards conservation applications).  Holmes (2003) presents a Bayesian approach, which uses posterior probability distributions to illustrate data support.  That method is adapted here for estimating the level of data support for the metapopulation risk metrics.  Let Y be a risk metric.  The probability that Y is greater than some threshold j given the data is 

              [16]

where  is the likelihood function given ,  is the likelihood function given ,  are the priors on mm and , and the normalizing constant is 

            [17]

The posterior distribution of Y is  over all j.  Examples of this calculation for l and the probability of being below thresholds at the end of 25 years are shown in the salmon examples.  Holmes (2003) supplies Splus code for these calculations.

23.6.  A simulated example
In this example, a collection of 49 local populations in a 7 x 7 grid was simulated with neighborhood dispersal.  Local populations were specified with variable mean local growth rates; thus, some mi’s were much larger than others.  The local growth rates in any given year were slightly correlated between sites.  Thus all sites were more likely than random to have good and bad years together.  Dispersal was variable between 5 and 10% from year to year and was mainly to the four nearest neighbors (or two and three for corner or edge sites).  In specific terms, A(t) was specified with zi(t)’s drawn from a normal distribution with mean = mi and variance of 0.0625.  The mi were different for each local population and were chosen randomly between -0.22 and -0.01.  Each year, new zi(t)’s were selected from the normal distribution for that local population.  The zi(t)’s were correlated among the local populations such that the covariance of zi(t) and zj(t) was (0.1)(0.0625).  The di(t) varied from year to year.  Each year and separately for each local population, di(t) was selected from a uniform random distribution between 0.05 and 0.1; thus the dispersal varied from year to year and between local populations in any given year.  Most of this dispersal, 80%, was to nearest neighbors.  Thus for nearest neighbors, aji = 0.80 dj(t)/nn, where nn is the number of nearest neighbors, and for non-nearest neighbors, aji = 0.2 dj(t)/nnn; where nnn is the number of non-nearest neighbors.
	The simulation was started from a set of local population sizes randomly drawn from the stochastic equilibrium, and the starting sizes were drawn anew from this distribution for each replicate of the simulation.  For each replicate, a 25-year time series was generated, and from this time series, mm and  were estimated using Eqn. 9.  From the estimates, the probability that the metapopulation would be below different thresholds (50% or 75% of starting levels) at the end of 25 years was predicted and compared to the actual probabilities obtained by repeatedly (1000 times) running the simulation for 25 years starting from the point where the initial 25-year time series stopped.  This simulation was replicated 500 times to generate the distribution of estimated probabilities of 50 and 75% decline in 25 years versus the true probability.  Also, from each 25-year simulation, the metapopulation l was estimated and compared to the actual value calculated by running a 10000-year simulation.  For each estimated risk metric, confidence intervals were estimated via the methods in section 23.5.
Figure 23.3 shows the distribution of l estimates and the estimated probabilities of 50 and 75% decline versus the true values.  As expected, the median estimate of l was equal to the true value ( is an unbiased estimator of mm).  The median estimate of l was 0.97 compared to the true value of 0.97.  The median estimates of 50 and 75% decline were 0.63 and 0.14 compared to the true values of 0.62 and 0.12, respectively.  Although the median estimates were very close to the true values, the estimates were variable.  The estimates of l ranged between 0.9 and 1.0.  The estimates of declines to thresholds were also variable.  The variability depended on the threshold and the time frame.  In this example, there was low variability around the estimate of 50% decline in 25 years, but high variability in the estimate of 75% decline.  The true values for each of the metrics are shown by the solid lines in the middle of the distributions.  
The variability of the estimates is due to the stochastic nature of the process and is not a fault of the estimation methods per se; by chance short trajectories will appear to have underlying parameters that are different than the true underlying parameters and this leads to variability in the estimated viability metrics.  When estimates are inherently variable, it is critical that the confidence intervals for the estimates be correctly estimated.  Figure 23.3 (right panels) confirms that the estimated confidence intervals properly characterize the uncertainty for the estimate risk metrics: e.g. 100(1-a)% of the time the 100(1-a)% confidence intervals contain the true values.

23.7.  Salmon as metapopulations
Salmonid populations (Oncorhynchus spp.) show strong spatial structuring and they have often been referred to as metapopulations (Reiman and McIntyre 1995; Policansky and Magnuson 1998; Cooper and Mangel 1999; Hill et al. 2002).  Spawning and rearing habitats of different salmon stocks occur on discrete and physically separated river or stream sections.  Salmon have a well-known and strong tendency to return to their natal streams with a low (1 to 20%) dispersal to other stocks (Fulton and Pearson 1981; Mathews and Waples 1991; Quinn 1993).  Within the U.S. Pacific Northwest, collections of anadromous salmon stocks have been divided into “Evolutionary Significant Units” or ESUs (Waples 1991) which represent substantially reproductively isolated conspecific groups that can be distinguished based on their coherence on a genetic level and known dispersal between the stocks.  Salmon within a stock spawn on individual streams or river sections and the majority of offspring return to spawn in their natal stream or river.  Straying of returning adults to non-natal streams is spatially structured and occurs more frequently within subbasins.  The stocks within an ESU have some level of synchrony due to exposure to common migratory corridors between the ocean and natal stream and also due to exposure to similar large-scale ocean dynamics (Pearcy 1992; Ware 1995; Mantua et al 1997).  However stocks also show a great deal of asynchrony due to exposure to their independent spawning and juvenile rearing habitats and variability in migration timing between stocks (e.g. PSTRT 2001).  Throughout the Pacific Northwest, most salmonid populations show regional decline with the majority of individual stocks showing steady declines with densities well below historical levels (Rieman and Dunham 2000; McClure et al. 2003)

23.8.  Snake River Spring/Summer Chinook ESU
The Snake River spring/summer chinook ESU (Fig. 23.4) includes all spring and summer chinook spawning within the subbasins of the Tucannon River, Grande Ronde River and the South, Middle and East Fork Salmon Rivers, which flow into the Snake River below the Hells Canyon Dam (Mathews and Waples 1991).  Juvenile fish rear in the mountain streams and then migrate down the Snake and Columbia Rivers to the ocean.  After maturing in the ocean, adult fish return to spawn at variable ages between 3 to 5 years (mean = 4.5 years).  Tagging experiments in Columbia River basin (which the Snake River basin is a part of) have found that the proportion of individuals that disperse and spawn away from their natal sites is on the order of 1-3% for wild-born individuals (Quinn 1993).  
The Snake River spring/summer chinook ESU was listed as threatened under the U.S. Endangered Species Act in 1992.  The stocks within this large and complex basin, like salmon stocks throughout the Pacific Northwest, are negatively impacted by a variety of factors (Wissman et al. 1994) and many have experienced substantial declines (Myers et al. 1998; McClure et al. 2003).  There is habitat degradation in many areas related to forestry, grazing, mining and irrigation practices, resulting in lack of pools, high temperatures, low flows, poor overwintering conditions, and high sediment loads in many areas.  At the same time, a substantial portion of the ESU is protected as part of federally designated wilderness (Mathews and Waples 1991).  The official ESU designation does not include salmon in the Clearwater basin since chinook in this subbasin originate from hatchery fish that were stocked in the subbasin after the original natural fish were extirpated in the 1940s.  However from a metapopulation dynamic perspective, the current stocks in the Clearwater River basin interact with the stocks within other subbasins.  Thus, in this analysis all stocks in the entire Snake River basin were analyzed together.  
A total metapopulation-level time series was available for this ESU from counts of the total number of wild-born spawners returning through the Ice Harbor dam at the downstream end of the ESU (Fig. 23.4).  Returning spawners can be either wild-born or hatchery-born as hatcheries have been operating in the basin since the early 1970s.  McClure et al. (2003) discuss the effects of hatchery production on viability analyses.  Effectively by focusing on the wild-born spawner time series and not incorporating a correction for hatchery production, the in-stream viability metrics assume that hatchery-born fish all return to the hatchery and do not spawn in-stream (which would produce wild-born offspring).  As discussed by McClure et al., this means our viability metrics are optimistic upper bounds, since some unknown fraction of hatchery fish do stray to the wild spawning grounds and potentially reproduce.  
In addition to the metapopulation-level dam count, time series of redds-per-mile (rpm), which are indices of the density of gravel egg nests made by spawning females, were available for the majority of stocks within the Snake River basin.  Redds-per-mile are an index of the redds (and consequently returning spawners) trend within a stock, but the total redds are unknown since the total spawning habitat is not surveyed.  The majority of the rpm and dam data are available in the digital appendices of McClure et al. (2003).

Parameter estimation
Our Ice Harbor dam time series starts in 1962 and ends in 1999.  The wild-born component of the dam count is denoted M(0), M(1), M(2), … M(37) where M(0) is the 1962 count and M(37) is the 1999 count.  The maximum likelihood estimates presented in Eqn. 9 assume that the data do not contain sampling error or other non-process error, however salmon data typically have high levels of sampling error and boom-bust cycles that confound estimation of mm and especially  (Holmes 2001).  An alternate approach uses data filtering and examination of the rate at which variance increases within the time series to improve parameter estimation and separate out sampling error variance from the time series (Holmes 2001; Holmes and Fagan 2002; cf also Morris and Doak 2002 pgs. 171-179).  These methods have been extensively cross-validated with salmon data (Holmes and Fagan 2002; Fagan et al. 2003) and are used here to estimate parameters.  First, the data are transformed using a running-sum:

                                            [18]

The estimates of mm and  are then

                [19]

The estimate of  uses the property that the variance of the underlying stochastic process should increase linearly with time: .  The confidence intervals for  and  using  are slightly different than Eqn. 10 (Holmes and Fagan 2002):

	.                         [20]

where L is the number of counts summed together for the running sum and 
df ≈ 0.333 + 0.212 (n+1) – 0.387 L = 6.84 (L = 4 and n = 38 here).  The estimated 95% confidence intervals on mm and  are (-0.133, 0.020) and (0.017, 0.111), respectively.
	 
Metapopulation viability metrics
The estimate of l for the Snake River spring/summer chinook ESU is .  To the extent that long-term trends continue, the expected population size in 25 years is 21% of current levels (=).  The point estimate of the probability of that the ESU drops below 10% of current levels at any time over the next 25 years is 

                           [21]

The corresponding estimate of 90% decline over the next 50 years is 0.74.  The probability of extinction was not estimated since this requires an estimate of the total population size.  The number of returning spawners is not the total population size since non-mature fish remain in the ocean.  However, if the true l of the metapopulation is less than 1.0, the population will eventually go extinct.
	The posterior probability density functions (Eqn. 16) for the estimated metrics are shown in Fig. 23.5.  The posterior probability distributions give an indication of the degree to which the data support different risk levels.  The distribution for l shows considerable data support for a l < 1, indicating a declining metapopulation.  There is also strong data support for a high risk of 90% decline over the next 50 years, however the estimate of 90% decline over 25 years is very uncertain.  The mean value is 0.23, but the probability distribution is very broad over the 0 to 1 range.  This illustrates that uncertainty in estimates of probabilities of quasi-extinction can vary widely depending on the time frame over which one is interested.

23.9.  Puget Sound Chinook ESU
The Puget Sound ESU is a subset of the major chinook salmon group in Washington’s northern coastal basins and Puget Sound.  The ESU (Fig. 23.6) includes all spring-, summer-, and fall-runs in the Puget Sound region from the North Fork Nooksack River to the Elwha River on the Olympic Peninsula (Myers et al. 1998).  The Elwha and Dungeness coastal basins of the Strait of Juan de Fuca, Hood Canal, and the Puget Sound area north to the northern Nooksack River Basin and the U.S. Canadian Border are all a part of the Puget Sound ESU.  Basin-to-basin dispersal rates have been observed at between 0.1 and 6% based on recoveries of tagged juveniles returning as adults (PSTRT 2001).  Fish in this ESU typically mature at ages 3 and 4 and are coastally-oriented during the ocean phase of their life history.  The Puget Sound ESU does not include Canadian or coastal Washington populations.  The Puget Sound ESU was listed as threatened under the Endangered Species Act in March of 1999.  Trends in abundance throughout the ESU are predominantly downward with several populations exhibiting severe short-term declines.  Degraded spawning and rearing habitats as well as access restrictions to spawning grounds and migration routes have all likely contributed to population declines.  Salmon in this ESU do not migrate through a hydropower system as the Columbia River ESUs do.
Our data for this ESU consist of yearly estimates of the total returning spawners (wild- plus hatchery-born) to the 44 separate river and creek systems feeding into the Puget Sound (Fig. 23.6).  These time series were compiled by the National Marine Fisheries Service (Seattle, WA) based on a variety of data: redd counts, carcass counts, in-stream harvest records, weir counts and hatchery return counts.  An independent metapopulation-level count was not available; unlike spawners returning to the Columbia River basin, spawners here do not pass through a hydropower system where they can be enumerated.  Instead a 1979-1997 index of the metapopulation was constructed by added together the 29 time series for the local populations with data over the 1979-1997 period.  As for the Snake River analyses, our viability metrics implicitly assumes that hatchery fish have not been reproducing and will be optimistic if some hatchery fish do not return to the hatchery and instead successfully spawn in the wild.

Metapopulation viability metrics
Parameters were estimated as for the Snake River.  The parameter estimates are , and .  The estimate of l for the Puget Sound Chinook ESU is .  The point estimate of the probability that the ESU drops below 10% of current levels at any time over the next 25 years is 0.000 and over the next 50 years is 0.001.  
	The posterior probability distributions (Fig. 23.5, right panels) illustrate the high uncertainty, given the data, as to whether this ESU is declining, stable, or increasing.  The most that can be said from the data is that there is low data support for a severely declining (l < 0.9) or increasing (l > 1.1) metapopulation.  Interestingly, the low support for small l’s translates into high data support for a low risk of 90% decline in the short-term (over 25 years).  Over the longer term, however, the uncertainty as to whether the metapopulation is declining or increasing gives rise to a U-shaped distribution meaning that the data give the most support to a probability of 0 or 1, reflecting that l could be either less than or greater than 1.0.  This example illustrates that while the data may be equivocal on some questions of conservation concern, such as “Is l < 1?”, the data may still give information on other questions, such as “Are the short-term risk of severe decline high?”.

23.10.  Using the stochastic metapopulation model to investigate effects of management
Determining how to distribute effort in order to recover an at-risk species is a routine, and challenging, task of conservation managers.  For salmon, management actions tend as a generality to affect an entire ESU or multiple ESUs or to affect individual stocks.  Management actions such as harvest reductions or increases to survival during migration (between spawning areas and the ocean) or improvements to estuarine environments are examples of actions that will tend to improve conditions for all stocks within an ESU or multiple ESUs.  Habitat improvements or protections that affect spawning areas and management of in-stream water levels are examples of actions that tend to affect individual stocks.  Without knowing the local stock dynamics or dispersal rates, one can still give certain types of guidance about how much effort is required for recovery of a declining metapopulation and about how effort should be distributed across all local populations.  

Metapopulation-level actions
When management actions affect all local populations roughly equally, it can be estimated how change would change the metapopulation l.  Mathematically, this means that all mi’s increase by some dm.  An absolute dm change in all mi’s is equivalent to multiplying all elements in A(t) by a constant = exp(dm).  The mean of the distribution of log M(t)/M(0) becomes .  Thus .  The change, dm, can be translated into currency that is more meaningful from a management standpoint by using the relationship l = R01/T, between l, the net reproductive rate, R0, and the mean generation time, T (Caswell 2001, Eqn. 5.73).  This is illustrated here for harvest and hydropower effects on salmon in the Snake River spring/summer chinook ESU (cf. McClure et al. 2003).

Harvest:  In the Pacific Northwest, harvest rates for salmon are generally expressed in terms of the fraction of spawners that did not return to the spawning grounds but that would have without harvest, e.g. a harvest rate of 0.8 indicates that the actual number of returning spawners is 20% of what it would have been if there had been no harvest.  Harvest rates are expressed in this way so that harvest that occurs in the stream versus in the ocean can be compared via a common currency.  We can write the net reproductive rate using fecundity and age-specific survival (cf. Caswell 2001, Eqn. 5.61)

      [22]

where h is the harvest rate, si is the survival from age i-1 to i, Fi is the fraction of spawners that return at age i, and f is the mean offspring per spawner.  Using Eqn. 22, the change in the l from a change in h alone is:

                                               [23]

Hydropower:  Juvenile salmon from the Snake River basin must migrate through the mainstem of the Snake River, enter the Columbia River and descend down the Columbia River on their journey to the ocean.  This migration, and the return migration of spawning adults, involves passage through four large hydropower dams on the Columbia River and four Snake River hydropower dams.  Improving the survival of both juvenile and adult fish migrating through the Columbia and Snake River hydropower systems has been the focus of much effort and is one of the human impacts that has been relatively well-quantified.  
Following a strategy similar to that used for harvest, the effect of changes in survival through the hydropower system on the rate of decline at the ESU-level can be estimated.  Denoting by cd and cu the proportional increase in down- and upstream passage survival due to improvement in the hydropower system, the improved net reproductive rate is 

.          [24]

Thus, for the assessing the impacts of increased survival through the hydropower system:

                                               [25]
Estimates of the impacts of harvest and hydropower changes to the Snake River ESU:  The mean ocean and in-river 1980-1999 harvest rate for the Snake River spring/summer chinook ESU was h = 0.08 (McClure et al. 2003).  By setting hnew = 0, we can examine the effect of successful selective harvest management that would substantially eliminate harvest impacts on salmon in this ESU.  Using Eqn. 23 and a mean generation time of 4.5 years, the estimated increase in l with hnew = 0 is roughly 2%.  NMFS (2000) has recently required that agencies operating the Federal Columbia River Power System implement a variety of activities, including increased spill, improved passage facilities, and increased barging of salmon around the dams as a means of improving survival through the system.  The estimated improvement in passage survival from the improvements proposed by NMFS are on the order of 5-6% (i.e. cdcu = 1.05-1.06) for the Snake River spring/summer chinook (McClure et al. 2003).  This translates into a 1% improvement in l for this particular ESU, using Eqn. 25.  Thus if the combined effects of substantially reduced harvest and the proposed passage improvements are additive, then roughly a 3% increase in l is estimated for these actions.  If the true l is less than 0.97, a 3% increase would not be sufficient to achieve l > 1.  Figure 23.5 indicates that the data cannot rule out that the l in this is ESU is greater than 0.97, but the data certainly give more support to a lower l.  This suggests that other recovery actions such as improvements at the stock-level will also be necessary.

Local population-level actions
The effects of changes to individual units of habitat are harder to quantify than the effects of metapopulaion-level changes.  The change in l achieved by a change at the level of a specific unit of habitat depends on the level of dispersal, the spatial pattern of dispersal, whether that habitat is connected to source or sink habitat, the level and pattern of synchrony between sites, etc.  In other words, it depends on the type of detailed information that has traditionally been difficult to obtain for metapopulations of conservation concern.  Interestingly although it is difficult to determine how much change in l can be achieved, it appears possible to estimate where the largest dl from a given dm change (per unit of habitat) in the local growth rate is achieved, even though the size of the resultant dl cannot be determined.  
Recall that each row of A represents a unit of habitat, and that a local population is comprised of some set of units of habitat with high connectivity.  When the intrinsic growth rate, mj, in a unit of habitat j is changed by dm ,  to exp(mj +dm), all the aij elements of column j in matrix A(t) are multiplied by exp(dm).  The goal is to calculate the total change in l from this dm change to all elements in column j by summing over rows i:

                                        [26]

The term  is the elasticity of l.  Caswell (2001), section 14.4, presents the calculation for the elasticity of l for products of stochastic matrices:

                                        [27]

where R(t) is the relationship between the right eigenvector and A(t), .  Thus, the dl from a dm change in a unit of habitat j can be solved for by summing Eqn. 27 over i.

                                [28]

This can be translated into matrix form:  

                                [29]

where “°” denotes the scalar product.  Using the relationship between the left eigenvector and A(t), ,

                                [30]

Thus dl from a change in a unit of habitat j is a weighted temporal average of the reproductive value of local population j times it’s density:

                                                 [31]

where c(t) is a constant that depends on t but not j.  A local population a is comprised of units of habitat in the set a = {a1, a2, a3, …, am}, where {a1, a2, a3, …, am} denotes which rows of A correspond to the units of habitat in local population a.  The dl per dm per unit habitat for a particular local population a is , where m is the number of units of habitat in local population a.  In words, this means the change in l is proportional to the product of the ‘average’ density of individuals in a particular local population times the ‘average’ reproductive value of its units of habitat.  
Although reproductive values are unknown, there are many cases where the product vjwj is a positive function of vj as long as dispersal is not too unidirectional (meaning, dispersal from A to B but not B to A).  This can be shown analytically in three extreme cases: a) 100% uniform and equal dispersal, b) all mj’s equal, or c) dispersal extremely low.  In cases a) and b), the reproductive values are all equal and vjwj = (a constant) x vj.  In case c),  and .  However, this positive relationship was also found in simulations with variable local growth rates, neighborhood dispersal, and dispersal sources and targets.  An obvious exception to this positive relationship is if dispersal is unidirectional, for example a linear chain of local populations with dispersal via a steady directional wind or ocean current.  However, as the simulations below illustrate, the general relationship can still hold even when dispersal is strongly, although not strictly, directional.

Density and l sensitivity
Three different types of metapopulation models were used to look at the relationship between average densities in units of habitat versus the dl from a small increase in the local growth rate in each unit of habitat.  In each model, dispersal was non-uniform among the local populations so that some sites were dispersal sources (more dispersal out than in) and others dispersal targets (more dispersal in than out). In the first model, local growth rates were equal among all sites and dispersers were globally distributed among all sites.  In the second model, local growth rates were variable so that some sites had much higher local growth rates than others and dispersal was mainly to nearest neighbors.  In the third model, local growth rates were again variable and dispersal mainly to the two south and east neighbors, however a small proportion of dispersers were distributed globally.  Thus the three examples illustrate global, local and directional dispersal.
A hundred randomly generated matrix models in each of these three categories were made and dlj calculated via Eqn. 30.  Figure 23.7 shows the relationship between the average density of a local population and the dl from increasing the local growth rate in that unit of habitat.  The x-axis ranks the dlj, thus “1” indicates the local population with the highest dlj in any simulation and “49”, the lowest.  The y-axis shows the corresponding mean density rank of that local population; “1” indicates the population had the highest density amongst the 49 sites and “49”, the lowest.  The results from the 100 randomly generated models are summarized by showing a box plot, which shows the median and range of all density ranks for the sites with a given dlj/dm rank.  Thus, the box plot at the x-axis position “1” shows the range of density ranks for the units of habitat with the highest dl/dm in each model.  The model results show a strong positive relationship between the relative density rank within a unit of habitat and which unit of habitat produced the largest increases in the metapopulation l for a given dm.  The two to three units of habitat with the highest average densities were consistently the units that produced the largest dl for a given dm.  This suggests that plotting the distribution of the relative densities within local populations in a metapopulation could give a rapid indication of the sensitivity of the metapopulation to changes to individual local populations.
One application of this would be to estimate where negative impacts would lead to the greatest decrease in l , thus suggesting where protection in most critical.  It would also suggest where local improvements would be most effective for a given increase in the local growth rate.  However in actual management situations where improvements are being sited, one is generally trying to maximize the “bang per buck”, .  The cost, d$, is actual monetary cost or some combination of monetary, logistical and political costs and dmj/d$ is the cost of a unit improvement to a unit of habitat j.  Thus dlj/dm is one part of the equation, and the other part, the cost of a unit improvement in different habitats, would have to come from a specific analysis of the costs and estimated effects of management actions on different local populations.

Example using the Snake River ESU:  The overall level of salmon dispersal between and among stocks within this ESU is known to be fairly low and spatially localized (Mathews and Waples 1991; Quinn 1993).  In addition, there is high variability in the habitat quality between stocks, with some stocks relatively pristine and protected within wilderness areas while others are exposed to high and multiple impacts (such as stream degradation and disturbance, pollution, in-stream harvest, and irrigation impacts).  Figure 23.8 (top) shows the distribution of average normed redds-per-mile per mile for 50 Snake River spring/summer chinook stocks.  For each year between 1980-1995, the redds-per-mile count for each stock was divided by the maximum count among the 50 stocks in that year.  The average over the 16 years was then used as an estimate of the average normed redds-per-mile.  The long-tailed distribution is the expectation from theory given low dispersal and high variability in stock habitat quality.  
Estimation of the average normed redds-per-mile was repeated using a variety of different time periods.  Regardless of the time period or number of years used for averaging, six stocks consistently appeared among the top five stocks with the highest density of redds: Johnson Cr., Poverty Cr. and Secesh R. in the South Fork of the Salmon R., the Lostine R. in the Grande Ronde subbasin, Marsh Cr. in the Middle Fork of the Salmon R., and the Imnaha R.  Perhaps not surprisingly, all of these are in relatively low impacted regions of the ESU.  At a subbasin level, the overall highest redd density was in the South Fork of the Salmon R. where summer-run chinook primarily occur (Fig. 23.8, bottom).  The other regions are primarily spring-run chinook.  The South Fork of the Salmon R. is relatively pristine and few hatchery fish have been released into this subbasin, and the stocks presumably have experienced relatively low interbreeding with hatchery-reared stocks.  In addition, the later run-timing may somehow be associated with less straying, lower harvest or lower hydropower impacts.
This analysis predicts that the l of the Snake R. spring/summer chinook ESU would be most sensitive to changes to the summer-run stocks in the South Fork of the Salmon River and to the spring-run stocks, the Lostine R., Imnaha River, and Marsh Creek and should be preferentially protected from impacts.  This can be counterintuitive in some situations.  For example, imagine making choices about where to allow a limited catch-and-release sport fishery.  Sites with the highest density would seem to be the prime candidates whereas the analysis of dl/dm indicates just the opposite.  In terms of determining where to direct improvements, the dl/dm suggests that these pristine sites are where a given dm would produce the greatest metapopulation l, however the regions where dl/dm is the highest are not necessarily the regions where m is most easily improved.  Indeed a given unit of improvement may be more difficult in pristine sites.  Choosing where to direct stock improvements requires consideration of the cost and difficulty of a given dm for different stocks in combination with the estimate of the sensitivity of l to local changes.

23.11.  PVA in practice
The purpose of this chapter is to present a theoretical framework for metapopulation PVA using time series data and diffusion approximations.  These methods are then illustrated using data from two salmon metapopulations.  The salmon analyses are intended as an example of how to calculate the diffusion parameters and metrics.  An actual PVA must grapple with other important issues that are outside the scope of this chapter, but which anyone contemplating an actual PVA must be aware.  Morris and Doak (2002, chapter 12) give a review of the criticisms and caveats surrounding the use of PVA and outline general recommendations and cautions when conducting a PVA.  In the context of diffusion approximation methods in particular, Holmes (2003) outlines an approach using matrix models to conduct sensitivity analyses in order to chose among different parameterization methods and metrics for a specific PVA application.  
One of the issues that is especially pertinent for our chapter is the issue of variability in estimated risk metrics.  A number of recent PVA cross-validations using actual data on a large number of different populations have shown that careful PVA analyses give unbiased risk estimates (Brook et al. 2000; Holmes and Fagan 2002; Fagan et al. 2003).  Although this is very encouraging, a difficult issue is the high inherent variability associated with estimated probabilities (such as probability of extinction) – even though they may be unbiased (Ludwig 1996, 1999; Fieberg and Ellner 2000; Holmes 2001; Ellner et al. 2002).  How to properly use risk metrics that have high variability is currently being debated within the field with arguments ranging from “don’t use them” (Ludwig 1996, 1999; Fieberg & Ellner 2000), to “use to estimate risks within collections of populations” (Fagan et al. 2001; Holmes and Fagan 2002), to “use where data are extensive and high quality” (Coulson et al. 2001), to “PVA metrics based on data, even if variable, are better than the alternatives” (Brook et al. 2002).  An encouraging aspect of diffusion approximation methods is that cross-validations using real time series data have indicated that the uncertainty in the estimated metrics appears to be properly characterized (Holmes and Fagan 2002).  Nonetheless, how to use and present metrics with high variability, albeit well characterized, is not an easy question to answer.  Presentation of 100(1-a)% is an oft used approach, but experience in the forum of salmon recovery planning in the Pacific Northwest has shown that it is easy to misinterpret confidence intervals.  For example, it is easy to interpret 95% confidence intervals for l that overlap 1.0 as an indication that the data are equivocal as to whether the population is declining or increasing; whereas, there may be considerable data support for a declining population.  Graphical presentations of data support for different risk levels have been more compelling and informative, although translating levels of data support into numbers that policymakers can use to take uncertainty into account in policy decisions has been challenging.

23.12.  Discussion
This chapter has focused on the calculation of metapopulation PVA metrics, however there are other more general PVA insights from an examination of stochastic metapopulations and of this specific class of declining density-independent metapopulations.  First, by definition the trajectory of a stochastic metapopulation is subject to random processes and thus the metapopulation trajectory observed in any one snippet of time is unlikely to capture the long-term dynamics.  The shorter the time frame, the farther the observed trend is likely to be from the long-term trend.  Thus the trends in any two adjacent time periods are unlikely to be identical, and the difference indicates not necessarily a change in the underlying rate of decline, can be due simply to chance.  The variability of observed rates of decline can be estimated from the level of the variability driving the long-term dynamics, and thus statistical tests performed to determine the likelihood that an apparent change in trend occurred due to the stochastic nature of the process rather than an underlying change in conditions.
Second, the local populations within a metapopulation are linked and experience the same long-term growth rates, regardless of the underlying difference in local population conditions (i.e. whether they are ‘sources’ or ‘sinks’).  However, measured over a short time period, there will be differences in the observed local population trends due to chance and local conditions.  This means that over a given time period, local populations will appear to be declining at different rates, but this is not an indication the long-term trends and not necessarily related local conditions being better or worse than other areas.  That the long-term trends of the individual local populations are the same as the metapopulation has a direct impact on PVA for local populations within a metapopulation.  The rate of decline observed among the different local populations will differ as will the apparent level of variability in the local time series.  Thus if an individual viability analysis is done using parameters estimated from local population time series alone, it will appear that there is tremendous variability among the local populations risk levels when in fact their long-term risks are similar.  When looking at the long-term risks, use of metapopulation-level parameters leads to better estimates of the long-term local population risks.  Short-term risks, however, are still strongly influenced by local conditions.  Clearly estimates of both short-term and long-term risks are needed to capture the whole viability picture for a metapopulation.  Although local populations within the type of metapopulations modeled here will be eventually repopulated by dispersal if they undergo extreme declines, the resulting loss of genetic diversity leads to a gradual erosion of the genetic health of the metapopulation.  Indeed this has happened for salmon species throughout the Pacific Northwest.  
Recovery planning for endangered and threatened species typically requires determining where to put the most effort.  Rarely is it the case that maximal effort can be applied everywhere.  Using the stochastic metapopulation model, a sensitivity analysis was used to look for local characteristics that predict where local changes would produce the biggest change in the metapopulation growth rate.  Interestingly, local density (not absolute numbers) was a strong predictor of the where a unit change in local growth rates led to the largest metapopulation growth rate.  This relationship was observed even in simulations with dispersal sources and targets and strongly directional dispersal, although it will break down when dispersal is strictly unidirectional.  Determining which local populations are best suited for restoration efforts also requires assessing the feasibility, cost, and acceptance of restoration efforts.  Indeed when it comes to actually implementing recovery actions, optimizing the efficiency of effort in terms affecting recovery requires solving a complex function of biological, economic and political information.  However, understanding the population dynamics of the species of concern and gaining insight regarding how the demography of the species will respond to alternative management actions are fundamental and primary components of this conservation equation.

23.13.  Acknowledgements
The graduate students in the Winter 2000 seminar on Spatial Ecology and Conservation Biology at the University of Washington were the initial impetus of this work.  Special thanks are given to Amanda Stanley who did many of the simulations that eventually led to the realization of strong patterns within the dynamics of stochastic metapopulations.  Jon Hoekstra (NMFS) provided us with data on the Puget Sound chinook ESU.  Discussions with Michelle McClure, John McElhany, and Mary Ruckleshaus (leads of the NMFS Technical Recovery Teams for Pacific NW salmon) have been instrumental at various stages of this work.  Thanks to Michael McCarthy whose detailed review and critiques led to many improvements in the presentation of the theory.

Literature Cited
Akçakaya, H. R.  1997.  RAMAS ® Metapop: viability analysis for stage-structured metapopulations. Version 2. Applied Biomathematics, Setauket, New York, USA.
Beissinger, S.R.  2002.  Population viability analysis: past, present, future. In S.R. Beissinger and D.R. McCullough. Editors. Population viability analysis. University of Chicago Press, Chicago.
Beissinger, S. R., and M. I. Westphal.  1998.  On the use of demographic models of population viability in endangered species management. Journal of Wildlife Management 62:821-841
Boyce, M. S.  1992.  Population viability analysis. Annual Review of Ecology and Systematics 23:481-506.
Boyce, M.  1996.  Review of RAMAS/GIS. Quarterly Review of Biology 71:167-168.
Brook, B.W., Burgman, M.A., Akcakaya, H.R., O'Grady, J.J., and Frankham, R.  2002.  Critiques of PVA ask the wrong questions: Throwing the heuristic baby out with the numerical bath water.  Conservation Biology 16: 262-263.
Brook, B. W., J. J. O'Grady, A. P. Chapman, M. A. Burgman, H. R. Akçakaya, and R. Frankham.  2000.  Predictive accuracy of population viability analysis in conservation biology. Nature. 404:385-387.
Caswell, H.  2001.  Matrix population models: construction, analysis, and interpretation.  Sinauer Press, Sunderland, MA, USA.
Cooper, A. B. and M. Mangel.  1999.  The dangers of ignoring metapopulation structure for the conservation of salmonids.  Fisheries Bulletin 97:213–226.
Coulson, J.C., G. M. Mace, E. Hudson, and H. Possingham.  2001.  The use and abuse of population viability analysis. Trends in Ecology and Evolution. 16: 219-221.
Dennis, B., P.L. Munholland and J.M. Scott.  1991.  Estimation of growth and extinction parameters for endangered species. Ecological Monographs 64:205-224.
Ellner, S. P., J. Fieberg, D. Ludwig, and C. Wilcox.  2002.  Precision of population viability analysis.  Conservation Biology 16:258-261
Fagan, W. F., E. E. Holmes, J. J. Rango, A. Folarin, J. A. Sorensen, J. E. Lippe, and N. E. McIntyre.  2003.  Cross-validation of quasi-extinction risks from real time series: an examination of diffusion approximation methods.  Submitted; 2nd revision in review.
Fagan, W.F., E. Meir, J. Prendergast, A. Folarin, and P. M. Kareiva.  2001.  Characterizing vulnerability to extinction for 758 species. Ecology Letters 4:132-138.
Fieberg, J., and S.P. Ellner.  2000.  When is it meaningful to estimate an extinction probability?  Ecology 81: 2040-2047.
Fulton, L. A., and R. E. Pearson.  1981.  Transplantation and homing experiments on salmon, Oncorhynchus spp., and steelhead trout, Salmo gairdneri, in the Columbia River system: Fish of the 1939-44 broods. U.S. Dep. Commer., NOAA Tech. Memo. NMFS F/NWC-12, 97p.
Furstenberg, H. and Kesten, H.  1960.  Products of random matrices.  Annals of Mathematical Statistics 31:457-469.
Gilpin, M.E. and M.E. Soule.  1986.  Minimum viable populations of species extinction. Pp. 19-34 In M.E. Soule (ed.). Conservation biology: the science of scarcity and diversity. Sinauer Associates, Inc., Sunderland, Massachusetts.
Heyde, C. C. and J. E. Cohen.  1985.  Confidence intervals for demographic projections based on products of random matrices.  Theoretical Population Biology 27:120-153.
Hill, M. F., A. Hastings, and L. W. Botsford.  2002.  The effects of small dispersal rates on extinction times in structured metapopulation models.  American Naturalist 160:389-402.
Holmes, E. E.  2001.  Estimating risks in declining populations with poor data.  Proceedings of the National Academy of Science 98:5072-5077.
Holmes, E. E. and W. F. Fagan. 2002.  Validating population viability analysis for corrupted data sets.  Ecology 83:2379-2386.
Holmes, E. E.  2003.  Beyond theory to application and evaluation: diffusion approximations for population viability analysis.  Manuscript conditionally accepted, Ecological Applications.
IUCN.  1994.  IUCN Red List Categories, Prepared by the IUCN Species Survival Commission As approved by the 40th Meeting of the IUCN Council, Gland, Switzerland.
Karlin, S. and H. M. Taylor.  1981.  A second course in stochastic processes.  Academic Press, New York, NY, USA.
Lacy, R.C.  1993.  VORTEX: A computer simulation model for Population Viability Analysis. Wildlife Research 20:45-65.
Lande, R.  1993.  Risks of population extinction from demographic and environmental stochasticity and random catastrophes. American Naturalist 142:911–927.
Lande, R., S. Engen, and B. Saether.  2003.  Stochastic population models in ecology and conservation: an introduction.  Oxford University Press, Oxford, UK.
Lande, R. and S. H. Orzack.  1988.  Extinction dynamics of age-structured populations in a fluctuating environment.  Proceedings of the National Academy of Science USA  85:7418-7421.
Ludwig, D.  1996.  Uncertainty and the assessment of extinction probabilities. Ecological Applications 6:1067-1076.
Ludwig, D.  1999.  Is it meaningful to estimate a probability of extinction?  Ecology 80:298-310.
Mantua, N. J., S. R. Hare, Y. Zhang, J. M. Wallace, and R. C. Francis.  1997.  A Pacific interdecadal climate oscillation with impacts on salmon production. Bulletin of the American Meteorological Society 78:1069-1079.
Mathews, G. M. and R. S. Waples.  1991.  Status review of Snake River spring and summer chinook salmon. U.S. Dept. Commer., NOAA Tech. Memo. NMFS F/NWC-200.
McClure, M., E. Holmes, B. Sanderson, and C. Jordan.  2003.  A large-scale, multi-species risk assessment:  anadromous salmonids in the Columbia River Basin.  In press in Ecological Applications.
Morris, W.F., P. L. Bloch, B. R. Hudgens, L. C. Moyle, and J. R. Stinchcombe.  2002.  Population viability analysis in endangered species recovery plans: past use and future improvements.  Ecological Applications 12:708-712.
Morris, W. F., and D. F. Doak.  2002.  Quantitative conservation biology: theory and practice of population viability analysis. Sinauer Associates, Inc, Sunderland MA.
Myers, J.M., R.G. Kope, G.J. Bryant, D. Teel, L.J. Lierheimer, T.C. Wainwright, W.S. Grand, F.W. Waknitz, K. Neely, S.T. Lindley, and R.S. Waples.  1998.  Status review of chinook salmon from Washington, Idaho, Oregon, and California. U.S. Dept. Commer., NOAA Tech. Memo. NMFS-NWFSC-35, 443pp.
NMFS (National Marine Fisheries Service).  2000.  Operation of the federal Columbia River power system including the juvenile fish transportation program and the Bureau of Reclamation's 31 projects, including the entire Columbia Basin Project. U.S. Department of Commerce, NOAA, National Marine Fisheries Service, Northwest Region, NW Region, Seattle.
Pearcy, W.G.  1992.  Ocean Ecology of North Pacific Salmon. University of Washington Press, Seattle, WA. 179p.
Policansky, D. and J. J. Magnuson.  1998.  Genetics, metapopulations, and ecosystem management of fisheries.  Ecological Applications 8:S119-S123.
Possingham, H. P., and I. Davies.  1995.  ALEX: a model for the viability of spatially structured populations. Biological Conservation 73:143-150.
Puget Sound Technical Recovery Team (PSTRT).  2001.  Independent populations of chinook salmon in Puget Sound.  Public review draft.  April, 2001.  Northwest Fisheries Science Center, NOAA, Seattle, WA, USA.
Quinn, T. P. 1993.  A review of homing and straying of wild and hatchery-produced salmon.  Fisheries Research 18:29-44.
Ralls, K., S. R. Beissinger, and J. F. Cochrane.  2002.  Guidelines for using population viability analysis in endangered-species management. Pages 521-550 in Population Viability Analysis (S. R. Beissinger and D. R. McCullough, Eds.). University of Chicago Press, Chicago, IL.
Reed, J. M., Mills, L. S.; Dunning Jr, J. B., Menges, E. S., McKelvey, K. S., Frye, R., Beissinger, S. R., Anstett, M. C., Miller, P.  2002.  Emerging issues in population viability analysis. Conservation Biology 16:7-19.
Rieman, B. E. and J. B. Dunham.  2000.  Metapopulations and salmonids: a synthesis of life history patterns and empirical observations.  Ecology of Freshwater Fish 9:51-64.
Rieman, B. E. and J. D. McIntyre.  1995.  Occurrence of bull trout in naturally fragmented habitat patches of varied size.  Transactions of the American fisheries Society 124:285-296.
Sanz, L. and R. Bravo de la Parra.  1998.  Variables aggregation in time varying discrete systems.  Acta Biotheoretica 46:273-297.
Shimizu, K. and K. Iwase.  1981.  Uniformly minimum variance unbiased estimation in lognormal and related distributions.  Communications in Statistics, Part A – Theory and Methods 10:1127-1147.
Shaffer, M.L.  1981.  Minimum population sizes for species conservation. BioScience 31:131-134.
Tuljapurkar, S. D.  1990.  Population dynamics in variable environments.  Lecture Notes in Biomathematics. Springer-Verlag, New York, NY, USA.
Tuljapurkar, S. and Orzack, S. H.  1980.  Population dynamics in variable environments. I. Long-run growth rates and extinction.  Theoretical Population Biology 18: 314-342.
Wade, P. R.  2001.  The conservation of exploited species in an uncertain world: novel methods and the failure of traditional techniques.  In Conservation of exploited species, (J. D. Reynolds and G. Mace, eds.), Cambridge University Press.
Waples, R. S.  1991.  Definition of "species" under the Endangered Species Act: Application to Pacific salmon. U.S. Department of Commerce, NOAA Tech. Memo. NMFS F/NWC-194, 29 pp.
Ware, D. M.  1995.  A century and a half of change in the climate of the NE Pacific. Fisheries Oceanography 4:267-277.
Wilcove, D. S., D. Rothstein, J. Dubow, A. Phillips and E. Losos.  1998.  Quantifying threats to imperiled species in the United States. Bioscience 48:607-615.
Wissmar, R. C., J. E. Smith, B. A. McIntosh, H. W. Li, G. H. Reeves and J. R. Sedell.  1994.  A history of use and disturbance in riverine basins of eastern Oregon and Washington (Early 1800s - 1900s). Northwest Science 68:1-35.

Figures
Figure 23.1.  Illustration of two of the main results from random theory:  (A) All local population go towards the same long-term rate of population growth (or decline) as t gets large.  (B) The distribution of log M(t) is a normal distribution with mean given by the long-term rate of growth (or decline) multiplied by t and the variance given by t multiplied by the rate that variance increases in an individual trajectory, i.e.  for t not overly small.  Here the variance was estimated from one time series using t = 10 and this used to predict the distribution at t = 100;

Figure 23.2.  Illustration of the performance of a diffusion approximation for modeling the behavior of a metapopulation with 50 local populations and uniform 0.1%, 1%, or 5% yearly dispersal.  The diffusion approximation performs well for a given time frame when  and  are constants over that time frame and when log M(t)/M(0) is normal.

Figure 23.3.  Estimated viability metrics and their estimated confidence intervals versus the true values for a 49-site metapopulation in a 7 x 7 grid with 5-10% dispersal to the closest neighboring sites.  The left panels compared the true metrics to the distribution of estimated metrics from 500 simulations starting from the same initial conditions.  The right panels show the performance of the estimated confidence intervals by looking at the fraction of estimated 100(1-a)% confidence intervals that contain the true values.

Figure 23.4.  Map of the Snake River spring/summer chinook ESU.  The ESU includes the stocks from the Snake River and its tributaries between the Ice Harbor and Hells Canyon dams.  The Hells Canyon hydropower dam has no passage facilities and blocks the migration of salmon into their historical habitat in the upper Snake River basin.

Figure 23.5.  Estimated posterior probability distributions for l, and the probability of 90% decline in 25 and 50 years.  The posterior probability distributions, which were calculated using uniform priors on  and , indicate the relative levels of data support for different risk metric values.  The left and right panels show the distributions for the Snake River Spring/Summer chinook and Puget Sound ESUs, respectively.

Figure 23.6.  Map of the Puget Sound chinook ESU.

Figure 23.7.  Relationship between the influence of a given habitat unit on the metapopulation l versus the average density in that habitat unit.  One hundred 7x7 metapopulations with spatially variable dispersal rates (some sites dispersal sinks and others targets) were randomly generated in each of three classes: 1) spatially uniform growth rates and global dispersal, 2) spatially variable growth rates with neighborhood dispersal and 3) spatially variable growth rates with directional neighborhood dispersal to the S and E two neighbors only.  The x-axis shows the rank in terms of dl/dm and the y-axis shows a box plot of the distribution of density ranks for sites with a given dl/dm rank across all 100 models in each class.  Thus the box plot at x = 1 shows the distribution of ranks for the sites with the highest dl/dm in each model.  The line in each box shows the median density rank for the sites with a given dl/dm rank, the box encloses 50% of the ranks, and the whiskers show the range from all 100 randomly generated models.

Figure 23.8.  Distribution of densities of redds in the Snake River spring/summer chinook ESU at a stock and subbasin level.  The average normed redds densities (upper plot) are shown for the 50 stocks with 1980-1995 data (the years were chosen to maximum the number of stocks with data).  For each stock the normed redd density was averaged over the 16 years to get an estimate of the normed average density.  In the lower plot, relative average densities over all stocks within different basins are shown (with the number of stocks in each basin shown above the bars).  The basin designations are: GR = Grande Ronde, I = Imnaha, SFS = South Fork Salmon, MFS = Middle Fork Salmon, US = Upper Salmon, C = Clearwater.  The redds due to hatchery fish released into stocks were removed before doing these analyses; otherwise the density will be artificially high simply due to hatchery fish releases.  This correction could not be done for the Upper Salmon or Clearwater regions since the fraction of spawners that are hatchery strays were unknown; however, the hatchery releases are very high in these basins and thus the corrected relative densities would be much lower than shown.

Box 23.1  Key results from random theory
Here are reviewed two of the fundamental results from the theory of products of random matrices and interpreted in the context of our metapopulation model.  The reader is referred to Caswell (2001) section 14.3.3 and Tuljapurkar (1990) section 2.1 for other reviews interpreted in the context of demographic, single population models.

The metapopulation and local populations decline at the same rate
One of the basic results from Furstenberg and Kesten’s “Products of random matrices” (1960), is that the product of ergodic random matrices asymptotically goes to an equilibrium.  Say that Xt is an ergodic random ‘k x k’ matrix and that Y (also a k x k matrix) denotes the product of n of the X matrices: X1, X2, X3,… Xn.  Then Furstenberg and Kersten’s results say that, Y goes an equilibrium state such that:

= a constant which is the same for all a               [1]

We can use this result to show that the long-run exponential growth rate of the metapopulation and the local populations will be the same.  

Thus from Eqn. 1,

== a constant = mm

The distribution of local population and metapopulation sizes is lognormally distributed
One of the most powerful results, for our purposes at least, concerns the statistical distribution of the metapopulation and local trajectories.  This tells us what distribution of sizes we would see if we ran our model over and over again, and allows us to make population viability analyses for metapopulations since we have a prediction about the likelihood of different metapopulation futures.  Random theory (Furstenberg and Kersten 1960, Tuljapurkar and Orzack 1980) shows that any sum of the Ni(t)’s, such as the total metapopulation (all i’s), a single local population (one i), or any other subset, goes to the same distribution:

                                           [2]

where the sum of local populations is denoted in matrix terms as  and c is a column vector with 0’s and 1’s to show which local populations to sum together.

Example
These results are simple to see with simulations.  Here is shown an example of a linear chain of ten local populations connected via 2% yearly dispersal to their nearest neighbors and 0.2% to non-nearest neighbors.  The local dynamics were  where zi is a normally distributed random variable, Normal(mi,).  The local growth rates, mi’s, for local populations 1 to 10 were respectively: 0.97, 1.00, 0.96, 0.83, 0.88, 1.00, 1.00, 0.89, 0.99, 0.81.  Figure 1A shows that the long-run growth rate of the local population and metapopulations is equal to the same constant.  Figure 1B shows that the distribution of metapopulation size after 100 years is .  The expected distribution was estimated using the maximum likelihood (ML) estimates for mm and  (Eqn. 9) from a single 1000-year time series of metapopulation counts.  The ML estimate for  relies on an assumption of normality for t = 1, although strictly speaking normality only holds for t large.  However, it does quite well as can be seen in Figure 1B.