Abstract
Background
Research on the ecological consequences of global climate change has elicited a growing interest in the use of time series analysis to investigate population dynamics in a changing climate. Here, we compare linear and nonlinear models describing the contribution of climate to the density fluctuations of the population of wolves on Isle Royale, Michigan from 1959 to 1999.
Results
The nonlinear self excitatory threshold autoregressive (SETAR) model revealed that, due to differences in the strength and nature of density dependence, relatively small and large populations may be differentially affected by future changes in climate. Both linear and nonlinear models predict a decrease in the population of wolves with predicted changes in climate.
Conclusions
Because specific predictions differed between linear and nonlinear models, our study highlights the importance of using nonlinear methods that allow the detection of nonlinearity in the strength and nature of density dependence. Failure to adopt a nonlinear approach to modelling population response to climate change, either exclusively or in addition to linear approaches, may compromise efforts to quantify ecological consequences of future warming.
Background
Several recent studies have suggested that species from many taxa and environments are already responding to regional climate change with, for example, poleward and elevational range shifts, changes in abundance, changes in body size, and changes in the timing of life history events [13]. Because these changes may influence important stabilizing mechanisms in communities and ecosystems, and, thus, the goods and services they provide to our society, and because the rates of such changes may continue to increase [3], they have moved to the forefront of political and environmental concern [1]. Nonetheless, there is still considerable uncertainty regarding the magnitude of future changes in climate on both local and regional scales and the impact they will have on populations, communities, and ecosystems. Clearly, studying the consequences of climate change for the dynamics and stability of populations will improve our ability to predict these effects and respond accordingly.
One growing body of research investigating this issue has used time series analysis of longterm datasets to quantify the contribution of largescale climatic fluctuation to population dynamics (reviewed in [4] and [5]). The approach commonly used in such studies is loglinear autoregressive modeling. The dynamical properties of this approach make it a logical, tractable, and appropriate candidate for examining and quantitatively predicting the influence of climate change on mean population density [6].
While appropriate in some situations, linear models may obscure potentially important properties of the population, which may ultimately influence the accuracy of predictions drawn from them. For example, linear models cannot detect nonlinear density dependence and the differential influence of abiotic factors above and below certain density thresholds (see, e.g., [7]). While several studies have indicated the potential for nonlinear densitydependence in natural populations, the consequences of incorrectly assuming linear densitydependence for the accuracy of model predictions is not well understood.
Here, we present a comparative analysis of linear and nonlinear models describing the contribution of climate to the interannual density fluctuations of the population of wolves (Canis lupus) on Isle Royale, Michigan, USA, from 1959 to 1999. Our analysis builds on previous studies indicating that the dynamics of the focal wolf population may be influenced by largescale climatic variation [8,9], and that the population may exhibit phasedependence [1012] that could generate nonlinearity in the dynamics [13]. Thus, there are indications that climate change may have important implications for the dynamics of wolves on Isle Royale, and, in addition, there are convincing arguments suggesting that nonlinear modeling may be appropriate for investigating this issue. We proceed by first developing a nonlinear selfexcitatory threshold autoregressive (SETAR) model for the focal population. We then develop a linear autoregressive model using similar methods and compare the stability properties, statistical fit, and predictions obtained by both models. Our analysis indicates that linear and nonlinear models may, in fact, lead to significantly different predictions of population response to climate change.
Results
The fit of the models
The most parsimonious nonlinear model identified here includes the number of packs in the previous year and both temperature and snow in the current year as covariates in the lowdensity regime; and pack size in the previous year, the number of packs in the previous year, and snow of the current year as covariates in the highdensity regime (Table 1). Piecewise nonlinear modeling indicated that the optimal density threshold for the SETAR model was 3.4, which on a raw scale is a density of approximately 30 wolves (Table 1). The nonlinear SETAR model is discontinuous at this threshold. In the lower regime, only the parameters for the delayed density and temperature covariates are significant at α = 0.1 Table 1). In contrast, only the parameters for delayed density and delayed number of packs covariates are significant in the upper regime (Table 1).
Table 1. The most parsimonious SETAR models of the population dynamics of wolves on Isle Royale, Michigan, 1959–99. Covariates included winter snow accumulation (SW_{t}), Northern Hemisphere temperature anomalies (T_{t}), number of packs in the previous year (PA_{t1}), and mean pack size in the previous year (PS_{t1}). X_{t }is log_{e}transformed density, and a_{i,j }are statistical parameters (i = 1 and 2 corresponds to lower and upper regimes, respectively, j = 0, 1, 2, 3 correspond to the constant, lag1 density coefficients, number of packs coefficients, and climatic coefficients, respectively). n indicates the number of data points in each regime. Parameter estimates were obtained by the method of conditional least squares. X* is the equilibrium point on the logscale.
The most parsimonious linear model includes the number of packs in the previous year and temperature of the current year as covariates and the parameters for these covariate are significant at α = 0.1 (Table 1). Given only the first observed value of the density time series and the observed values of snow, temperature, the number of packs, and pack size, both models simulate dynamics reasonably well (Fig. 2). However, AIC scores indicate a better fit for the nonlinear model (Table 1).
Stability properties
The stability of the dynamics within the each regime of the nonlinear model and for the linear model were determined by the values of the coefficients of density dependence (a_{1,1}, a_{2,1}, and a_{1 }in Table 1). Ninetyfive percent confidence intervals for the coefficient of direct density dependence in the lower regime (0.05 ± 0.28) and the upper regime (2.63 ± 0.13) do not overlap suggesting that the strength of density dependence differs between the two regimes. The equilibrium point in the lower regime is stable, while that for the upper regime is not (Table 1). The equilibrium point in the linear model is stable (Table 1). The location of the equilibrium points in both regimes of the nonlinear model and in the linear model, however, depends on climate (Table 1). If the average values of pack size (8.24), number of packs (3.7), temperature (0.10), and snow (190.4) from the observed time series are used, then the equilibrium points of the lower and upper regime of the nonlinear model and of the linear model are approximately 11, 36, and 20 wolves on the raw scale, respectively (Table 1).
Discussion
Nonlinear SETAR model
Previous work has revealed that density dependence in the wolf population on Isle Royale is phasedependent [12], which may contribute to nonlinearity in the dynamics of this population. Importantly, our results also indicate that there are differences in the nature of density dependence between low and highdensity regimes (Table 1, Fig. 1). Analysis of the model suggests that populations below the density threshold experience negative density dependence and exhibit monotonic damping, while populations above the threshold display positive density dependence.
Figure 1. Population growth rate (R_{t }= X_{t } X_{t1}) of wolves versus lntransformed density in the previous year (X_{t1}). Solid lines indicate nonlinear spline functions, and dashed lines indicate 95% confidence bands, estimated with generalized additive models. These plots suggest a threshold at ~3.25–3.50 wolves on a log_{e }scale. The linear relationship below this threshold suggests weak density dependence while the relationship above this threshold suggests that density dependence is stronger at higher densities.
The social behavior of wolves provides a convincing biological explanation for the presence of a density dependence threshold and positive density dependence above the threshold. Because they live on an island and, thus, have a finite amount of space, growth in the population must be accompanied by an increase in the number of wolves per pack and/or the number of packs on the island. A closer examination of the data demonstrates that at low densities, population growth on Isle Royale is accompanied by an increase in average pack size (Fig. 3). The number of packs on the island remains relatively constant (3 to 5 packs). When the total population size reaches 30 wolves, however, the number of packs begins to increase (Fig. 3). The average pack size remains relatively constant at 10–15 wolves per pack. These data suggest that there may be an upper limit to the size of a pack, possibly due to difficulties obtaining enough to food to feed all pack members or due to an increase in aggressive behaviors. This upper limit in pack size may generate a density threshold. In addition, this shift in behavior may lead to positive density dependence at high densities. An increase in the number of packs corresponds to an increase in the number of breeding pairs on the island. As density continues to increase, more packs are formed, leading to higher population growth. However, and importantly, this positive feedback loop is inherently unstable. As the number of packs on the island increases, the amount of territory per pack decreases and the frequency of aggressive encounters among packs must increase. Together, these processes should lead to increased mortality. In fact, the observed dynamics of the population support this theory. The first time the population reached 30 wolves in 1965, it declined over the next four years, and the two times it has crossed this threshold, between 1978 and 1980, it has crashed precipitously [12].
Figure 2. Dynamics predicted by nonlinear SETAR model and linear model. Predicted dynamics were calculated from the full models in Table 1 using one seed density estimates from 1960 and the observed climate and pack time series (one step ahead simulations). (a). Dynamics predicted by SETAR model with winter snow accumulation, SW_{t}, as a covariate. (b). Dynamics predicted by linear model with the number of packs, PA_{ti}, and the Northern Hemisphere winter temperature anomaly, T_{t}, as covariates. (c). Predicted versus observed values for the nonlinear model. (d). Predicted versus observed values for the linear model.
Figure 3. (a) Number of packs (PA_{t}) versus Number of wolves (N_{t}). (b) Average pack size (PS_{t}) versus Number of wolves (N_{t}). Threshold for nonlinear model is ~30 wolves on the raw scale.
Population response to climate change
Differences in the nature of density dependence above and below the threshold suggest that the consequences of climate change for the dynamics of Isle Royale wolves differ depending on if the population is relatively small or large. At low densities dynamics are largely determined by within pack dynamics, specifically their ability to capture prey and feed all pack members. Because changes in climate influence the ability of wolves to capture prey, climate change should have significant effects when the population is small. Because dynamics maybe determined more by interpack dynamics and less by the ability of individual packs to capture prey at high densities, climate and climate change should be less important in determining population dynamics when the population is large. Climate may influence the dynamics of large populations, but the signature of climatic effects may be diluted by the effects of interpack aggression. These predictions are supported by the fact that though the most parsimonious nonlinear model contains climate as covariates in both regimes, the only significant coefficient of a climatic variable occurs in the lowdensity regime (Table 1). Thus, our results suggest that the dynamics of this population, when its density is low, are stable, determined mainly by direct density dependence, and affected by changes in climate. At higher densities, the dynamics are intrinsically unstable and the population is less influenced by climate.
In addition to predicting when climate and climate change is important, the nonlinear SETAR model suggests that predicted changes in climate will cause the population of wolves on Isle Royale to decline, possibly to extinction. While wolf density is negatively related to both temperature and snow in the low density regime, the coefficient for snow is not significant (Table 1). This suggests that below the density threshold, changes in snow depth will have a lesser effect on the population than changes in temperature. Hence, predicted increases in temperature may lead to a decrease in the population of wolves on Isle Royale, possibly contributing to extinction risk. At densities above the threshold, the population is unlikely to be systematically affected directly by changes in climate. However, because dynamics are unstable at high densities, large populations are likely to crash, falling below the threshold where predicted changes in climate will cause the population to decline.
While increases in mean temperature may cause the population to decline, this warming trend may also be accompanied by increases in variability that could have additional effects on the population. Extreme variations in climate may cause the population of wolves on Isle Royale to fluctuate above and below the density threshold. This may cause extreme fluctuations in the size of the population greatly increasing the chance of extinction.
Comparison to linear model
Qualitatively, the linear model obtained here predicts similar responses to climate change in comparison to the nonlinear model: increasing temperatures will lead to declines in population density. However, the linear model indicates negative density dependence for the whole range of observed densities and suggests that the number of packs is important in determining the dynamics of the population. Importantly, the absence of a density threshold and positive density dependence in this model suggests that the population should experience relatively small fluctuations in density in constant and variable environments. Negative density dependence should regulate the population at all densities and buffer the population from changes in climate. Thus, in contrast to the nonlinear SETAR model, the linear model predicts similar effects of climate and climate change on relatively small and large populations, predicts that the population is relatively stable, and suggests that population fluctuations will be relatively small.
In terms of statistical fit to the data, the nonlinear SETAR model performs better than the linear model, indicated by the fact that the AIC of the linear model is almost double that of the nonlinear model (Table 1). Because these are phenomenological models, this criterion does not necessarily mean that the nonlinear model more accurately portrays the biology of the system or that it is better for the purposes of predicting response to climate change. In fact, we urge caution in interpretations of both models. For the nonlinear model in particular, there are only eight data points, but relatively many parameters in the upper regime. This decreases certainty surrounding parameter estimates and thus, our conclusions regarding stability and response to climate change for large populations. In addition, despite differences in AIC scores, the linear model still provides a reasonably good fit to the data. It may, therefore, give just as reliable predictions as the nonlinear model. The main contribution of our analysis is its demonstration of how linear versus nonlinear time series modeling approaches may influence conclusions regarding population response to climate change.
The importance of nonlinear timeseries modelling
Using time series data for examining population dynamics has both disadvantages and advantages, regardless of whether linear or nonlinear approaches are applied. Our ability to model the dynamics of a population and the reliability of the predictions that can be drawn from it depend on the data available and the length of the time series. As mentioned above, trophic interactions were not directly included in either model but could be very important in determining not only the dynamics of the population of wolves on Isle Royale, but also its response to climate change. In addition, though this data set represents one of the longest time series of a threetrophic level system, it spans only 40 years and teasing out the relative influences of intrinsic versus extrinsic processes on the dynamics of longlived species may require longer time series [14]. In addition, it is possible that the factors most important in determining dynamics vary through time. For example, it is possible that the population was exposed to canine parvovirus in the 1980s [11]. The dynamics of the population during this time period may differ in important ways that cannot be captured in a simple model of the full time series lacking this information.
Despite the potential difficulties, the analysis of time series data offers several advantages. First, these methods provide one of the only means for examining population dynamics when experimentation and detailed demographic analyses are impractical or impossible, as they are on Isle Royale. Second, these methods allow analyses of population dynamics when only limited data are available. Finally, time series analysis is useful when examining processes that occur over large spatiotemporal scales such as climate change.
Conclusions
While time series methods may allow us to predict the effects of climate change on population dynamics, the reliability of those predictions and our mechanistic understanding of population dynamics is influenced by whether linear or nonlinear approaches are taken. Linear models may be more appropriate for determining response of the whole system to changes in climate, however, linear models assume that the strength of density dependence and, thus, the susceptibility of the population to climate change, are independent of density (sensu [6]). The nonlinear model presented here, however, suggests that the nature and strength of density dependence may be nonlinear (i.e. differs between high and lowdensity regimes) (Table 1) and, because of this, the manner in which the population responds to climate change may differ depending on if the population is at relatively high or relatively low densities. By specifically examining the possibility of nonlinearity, we were able to suggest mechanisms through which changes in wolf behavior may influence population dynamics and their response to changes in climate. Thus, we suggest that because density dependence determines, in part, population response to climate change, it is important to investigate methods that allow the detection of nonlinearity in the nature and strength of this force and other potentially important processes.
Methods
Developing the model
Wolves on Isle Royale are the top predator in a straightchain, three trophiclevel system including moose (Alces alces) and balsam fir (Abies balsamea) [16]. Therefore, the dynamics of the wolves in this system are most broadly influenced by trophic interactions with moose and intrinsic density dependence. Because wolves prey primarily on young and old moose, their population dynamics may depend more on the density of these age groups rather than the total density of moose. Unfortunately, there is very little information regarding the age structure of the moose population. In addition, estimates of moose density are based upon cohort reconstruction and, therefore, may be less accurate than those of wolf density, which are based on direct counts. For these reasons, in addition to previous studies suggesting that a univariate equation in delay coordinates may indirectly incorporate trophic interactions [17], we begin with a simple singlespecies linear autogressive model:
X_{t }= α_{0 }+ α_{1}X_{t1 }+ α_{2}X_{t2 }+ .... + α_{n}X_{tn }+ ε_{t } (1)
where X_{i }is log_{e}transformed wolf density at time i and ε_{i }is white noise. Although this simplification compromises detail in terms of the ecology behind wolf dynamics, it may facilitate further manipulations while, at some level, incorporating trophic interactions between wolves and moose (see, e.g. [13,18]).
Partial autocorrelation function (PACF) plots and recent first and secondorder linear and piecewise nonlinear modeling of the wolf time series indicate that a firstorder model provides a better fit in terms of statistical parsimony than secondorder or higher models [9,12]; A. Ellis, unpublished results) giving:
X_{t }= α_{0 }+ α_{1}X_{t1 }+ ε_{t } (2)
To examine the structure of density dependence we set X_{t } X_{t1 }= R_{t}, simplify coefficients, and express equation (2) in terms of the population growth rate of wolves:
R_{t }= β_{0 }+ β_{1}X_{t1 }+ ε_{t } (3)
The scatter plot of R_{t }versus X_{t1 }fit with a spline function (Fig. 1) indicates that the relationship between population growth rate in the current year and density in the previous year may be nonlinear. While this inference is supported by the reduction in residual error achieved by applying a generalized additive model [19] with a nonlinear spline function (deviance = 1.44) compared to a linear regression (deviance = 2.19) between R_{t }and lagged density, more data points are needed to verify the robustness of this claim. For the purposes of this study, this plot suggests that the strength of density dependence may differ between relatively small and relatively large populations around some threshold (in X_{t1}) and that this nonlinearity may appropriately be expressed with simple linear piecewise models.
A firstorder selfexcitatory threshold autoregressive (SETAR) model with a low and highdensity regime and a threshold in X_{t1 }has the general form:
X_{t }= a_{1,0 }+ a_{1,1}X_{t1 }+ ε_{1,t } X_{t1 }< θ (4)
X_{t }= a_{2,0 }+ a_{2,1}X_{t1 }+ ε_{2,t } X_{t1 }≥ θ
where, X_{t }and X_{t1 }are log_{e}transformed densities of the current and previous years, respectively, a_{i,j }are statistical parameters (i = 1 and 2 corresponds to lower and upper density regimes, respectively, j = 0 and 1 correspond to the constant and lag1 coefficients, respectively), ε_{i,t }is white noise (i = 1 and 2 corresponds to lower and upper density regimes, respectively), and θ is the threshold set in terms of density [7,20].
Covariates
Clearly, wolf dynamics are not determined exclusively by the direct influence of density dependence and the indirect influence of moose density. Wolf dynamics may also be influenced by social structure [9,12,21], local climate [22], and largescale climate [9]. There is, for example, typically one breeding pair per pack [23], so that the birth rate of the population depends on the number of packs in the population. A closely related aspect of social structure, pack size, may also influence dynamics because pack size may affect ability to capture prey [9,12]. Climate, specifically snow depth and the North Atlantic Oscillation (NAO) [24], may also influence wolf dynamics through its influence on prey vulnerability and kill rate [8,22]. Because snow depth itself is influenced by winter temperature (R^{2 }= 0.14 for the period of 1959 to 1999), and because the NAO may represent the vehicle through which largescale warming occurs [25], it is reasonable to suggest that winter temperature may also influence wolf dynamics on Isle Royale though, to our knowledge, this influence has not been examined.
Hence, we examined all possible nonlinear SETAR models of the form specified in equation (4) with covariates of mostrecent lags including pack size of the previous year, number of packs of the previous year, cumulative winter snow depth of the current year (Ontario Climate Center), the Northern Hemisphere winter temperature anomaly of the current year, and the North Atlantic Oscillation of the current year [26]. Because estimates of cumulative snowfall and temperature on Isle Royale are not available, we used snowfall measurements from nearby Thunder Bay, Ontario and the Northern Hemisphere winter temperature anomaly.
Threshold selection
Before obtaining parameter estimates and comparing models, we determined the optimal threshold for each possible model. The optimal threshold for a particular model based upon the skeleton (Eq. 4) is the lagone density term (X_{t1}) from the observed unique values of X_{t1 }that minimizes the negative loglikelihood value for the entire model. This was determined by analyzing each model with θ in every unique value of X_{t1}, except for the lowest and highest values to ensure minimum degrees of freedom to perform a linear regression. The threshold that minimized the negative loglikelihood was the optimal threshold for that particular model.
Model selection
All possible SETAR models with the form of equation (4) and the covariates mentioned above were analyzed to find the most statistically parsimonious model. For comparison, all possible linear models of the form of equation (2) and the covariates mentioned above were also analyzed and the most parsimonious model selected. We examined models with the separate effects of snow, temperature, or the NAO and models that included two or three climate variables to investigate the interaction between them. The most parsimonious model was that which minimized the corrected Akaike's Information Criterion (AIC) score [27]. Selected models were then tested for multicollinearity. Pairs of covariates resulting in variance inflation factors greater than 10 were not simultaneously included in the model, excluding the covariate with lower explanatory power when necessary.
Stability
Stability theory predicts that the greater the strength of density dependence, the more likely a population is to remain within a range of densities despite environmental change [28]. We inferred the stability of the system on the basis of the value of the coefficient of direct density dependence (a_{1,1}, a_{2,1 }in Eq. (4) and α_{1 }in Eq. (2)) in the regime where the equilibrium point lies [18,29,30].
Authors' contributions
AE designed the study, conducted the analyses, and did most of the writing. EP conceived of the study, supervised its design, and contributed to the writing. Both authors read and approved of the final manuscript.
Acknowledgements
We graciously thank Ottar Bjørnstand for his insight and contributions to this study. His help was invaluable. We also wish to thank Mads C. Forchhammer, Joe Kiesecker, Katriona Shea, James Lettieri, and Jeff Haeni for their helpful comments. This research was supported by a graduate fellowship to A.E. and a grant from NSF to E.P.
References

McCarthy JJ, Canziani OF, Leary NA, Dokken DJ, White KS, eds: Climate change: Impacts, adaptation, and vulnerability. Contribution of Working Group II to the Third Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge: Cambridge University Press; 2001.

Hughes L: Biological consequences of global warming: is the signal already apparent?
TREE 2000, 15:5661. PubMed Abstract  Publisher Full Text

Walther GR, Post E, Convey P, Menzel A, Parmesan C, Beebee TJC, Fromentin JM, Guldberd OH, Barlein F: Ecological responses to recent climate change.
Nature 2002, 416:389395. PubMed Abstract  Publisher Full Text

Post E, Forchhammer MC, Stenseth NC: Population ecology and the North Atlantic Oscillation (NAO).

Ottersen G, Planque B, Belgrano A, Post E, Reid PC, Stenseth NC: Ecological effects of the North Atlantic Oscillation.
Oecologia 2001, 128:114. Publisher Full Text

Ives AR: Predicting the response of populations to environmental change.

Grenfell BT, Wilson K, Finkenstädt B, Coulson TN, Murray S, Albon SD, Pemberton JM, CluttonBrock TH, Crawley MJ: Noise and determinism in synchronized sheep dynamics.
Nature 1998, 394:674677. Publisher Full Text

Post E, Peterson RO, Stenseth NC, McLaren BE: Ecosystem consequences of wolf behavioural response to climate.
Nature 1999, 401:905907. Publisher Full Text

Post E, Forchhammer MC: Pervasive influence of largescale climate in the dynamics of a terrestrial vertebrate community.
BMC Ecology 2001, 1:5. PubMed Abstract  BioMed Central Full Text

Peterson RO, Thomas NJ, Thurber JM, Vucetich JA, Waite TA: Population limitation and the wolves of Isle Royale.

Peterson RO, Page RE: The rise and fall of Isle Royale wolves, 1975–1986.

Post E, Stenseth NC, Peterson RO, Vucetich JA, Elllis AM: Phase dependence and population cycles in a largemammal predatorprey system.

Stenseth NC: Population cycles in voles and lemmings: density dependence and phase dependence in a stochastic world.

Peterson RO, Page RE, Dodge KM: Wolves, moose, and the allometry of population cycles.

Ginzburg LR, Ferson S, Akçakaya HR: Reconstructibility of density dependence and the conservative assessment of extinction risks.

McLaren BE, Peterson RO: Wolves, moose, and tree rings on Isle Royale.

Forchhammer MC, Asferg T: Invading parasites cause a structural shift in red fox dynamics.
Proceedings of the Royal Society of London, Series B 2000, 267:779786. Publisher Full Text

Royama T: Analytical population dynamics. London: Chapman & Hall; 1992.

Hastie TJ, Tibshirani RJ: Generalized additive models. London: Chapman & Hall/CRC; 1999.

Tong H: Nonlinear time series analysis. Oxford: Clarendon Press; 1990.

Vucetich JA, Peterson RO, Waite TA: Effects of social structure and prey dynamics on extinction risk in grey wolves.
Conservation Biology 1997, 11:957965. Publisher Full Text

Fuller TK: Effect of snow depth on wolf activity and prey selection in north central Minnesota.

Mech LD: The wolf: the ecology and behavior of an endangered species. Minneapolis: University of Minnesota Press; 1970.

Hurrell JW: Decadal trends in the North Atlantic Oscillation: regional temperatures and precipitation.

Visbeck MH, Hurrell JW, Polvani L, Cullen HM: The North Atlantic Oscillation: past, present and future.
Proceedings of the National Academy of Sciences, USA 2001, 98:1287612877. Publisher Full Text

University of East Anglia, Norwich, NR 2002.
http://www.cru.uea.ac.uk/ftpdata/tavenh.dat webcite, http://www.cru.uea.ac.uk/~timo/projpages/naoupdate.htm webcite

Hurvich CM, Tsai CL.: Regression and time series model selection in small samples.

Ives AR, Gilchrist G: Climate change and ecological interactions. In Biotic interactions ang global change. Edited by Karieva PM, Kingsolver JG, Huey RP. Boston: Sinauer Associates; 1993:120146.

May RM: Stability in randomly fluctuating versus deterministic environments.
The American Naturalist 1973, 107:621650. Publisher Full Text

May RM, Conway GR, Hassell MP, Southwood TRE: Time delays, densitydependence and singlespecies oscillations.