# Population regulation in snowshoe hare and Canadian lynx: Asymmetric food web configurations between hare and lynx

## Abstract

The snowshoe hare and the Canadian lynx in the boreal forests of North America show 9- to 11-year density cycles. These are generally assumed to be linked to each other because lynx are specialist predators on hares. Based on time series data for hare and lynx, we show that the dominant dimensional structure of the hare series appears to be three whereas that of the lynx is two. The three-dimensional structure of the hare time series is hypothesized to be due to a three-trophic level model in which the hare may be seen as simultaneously regulated from below and above. The plant species in the hare diet appear compensatory to one another, and the predator species may, likewise, be seen as an internally compensatory guild. The lynx time series are, in contrast, consistent with a model of donor control in which their populations are regulated from below by prey availability. Thus our analysis suggests that the classic view of a symmetric hare–lynx interaction is too simplistic. Specifically, we argue that the classic food chain structure is inappropriate: the hare is influenced by many predators other than the lynx, and the lynx is primarily influenced by the snowshoe hare.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

The cyclic changes in abundance of the snowshoe hare (

*Lepus americanus*Erxleben, 1777) and the Canadian lynx (*Lynx canadensis*Kerr, 1792) are well known (1–4). These 9- to 11-year fluctuations are commonly discussed in ecology texts (e.g., refs. 5 and 6) as examples of coupled predator–prey cycles (7–13).Even though the biodiversity of the boreal forest is low (14), it is still a too complex ecosystem to be modeled intelligibly (15). By focusing on those species directly connected with the hare, a smaller and more tightly interlinked food web emerges (Fig. 1

*A*). This is even more so when focusing on the species directly connected to the lynx (Fig. 2*A*).Figure 1

Figure 2

To estimate the number of key interactions determining the dynamics of the snowshoe hare and the lynx, we have analyzed time series data on these species (Figs. 1

*B*and 2*B*) and interpreted the results on the basis of recent ecological data [primarily from the Kluane Ecosystem Project (4, 16)]. The main statistical result of our investigation is that the embedding dimension for the lynx series is roughly two, whereas that of the hare series is closer to three. We discuss these statistical results from the point of view of ecological interactions: in spite of the apparent complexity of the underlying food web, the dynamics of the hare and the lynx seem to be rather simple (low dimensional). On the basis of the documented pattern, we propose two ecological models—one for the lynx and one for the snowshoe hare.## The Data

The time series data derive from the compilation of fur records on hares (17) and lynx (2) carried out by Charles Elton, Helen Chitty, and others of the

*Canadian Snowshoe Rabbit Enquiry*(ref. 18; see also ref. 19). It was partly through this inquiry that ecologists were convinced that the vertebrate cycles of the boreal forest were not just an artifact of the trapping or smoothing of random numbers (20).The snowshoe hare data derive from the main drainage of Hudson Bay, whereas the lynx data are from 10 different regions across boreal Canada (the two hare series correspond regionwise to the combined James Bay and Lakes lynx series; L13 in Table 1). Although benchmark data on the abundance of cyclic species, they have some drawbacks. (

*i*) The hare data presented by MacLulich (17) as one time series are really two different sets of data (1844–1904 represents fur records, whereas 1905–1935 derives from questionnaires; Fig. 1*B*); they thus ought to be analyzed as two series. (*ii*) The data on the different regions of the lynx are of variable lengths (extending between 1821 and 1939; most end in 1934) and contain occasional missing values (for most series, data are lacking for 1892–1896 and 1914). (*iii*) The regions to which the lynx series refer vary somewhat in demarcation over time (2). Among the lynx series, the one from northwestern Canada (the MacKenzie River district adjacent to the area studied in the Kluane Ecosystem Project (61°N, 138°W)) is the one most commonly quoted (L3 in Table 1).Table 1

No. | Time series | Years | General model: X = _{t}F(X_{t−1}, X_{t−2}, ⋯ , X_{t−d}) + ɛ_{t} | Additivity (P value) | The General Additive Model: X = _{t}f_{1}(X_{t−1}) + f_{2}(X_{t−2}) + ⋯ + f(_{d}X_{t−d}) + ∂_{t} | |||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Dimension | Linearity (P value) | Dimension | Slope at equilibrium of the GAM (linear coefficients ± 2 SEM) | |||||||||

N-W* | LL reg* | H-T | Tsay | Linear* | GAM* | lag1 (X_{t−1}) | lag2 (X_{t−2}) | lag3 (X_{t−1}) | ||||

H1 | Snowshoe hare | 1844–1904 | 1 (0.67;0.74) | 1 (0.70;0.79) | √ (0.18) | √ (0.07) | √ (0.11) | 4 (0.69;0.76) | 4 (0.59;0.66) | 0.82 (0.62 ± 0.26) | 0.13 (0.1 ± 0.30) | −0.42 (−0.23 ± 0.26) |

H2 | Snowshoe hare | 1905–1935 | 3 (0.39) | 3 (0.33) | √ (0.39) | √ (0.52) | √ (0.13) | 3 (0.30) | 3 (0.22) | 1.02 (0.95 ± 0.33) | 0.14 (0.01 ± 0.49) | −0.56 (−0.49 ± 0.33) |

L1 | West | 1825–1856 | 2 (0.38) | 3 (0.36;0.38) | √ (0.09) | √ (0.31) | √ (0.77) | 4 (0.35;0.39) | 4 (0.34;0.44) | 1.00 (1.12 ± 0.33) | −0.93 (−0.55 ± 0.33) | |

L2† | West | 1897–1934 | 2 (0.28) | 2 (0.25) | √ (0.87) | √ (0.16) | √ (0.08) | 2 (0.24) | 4 (0.23;0.25) | 1.38 (1.34 ± 0.26) | −073 (−0.67 ± 0.26) | |

L3‡ | MacKenzie River | 1821–1934 | 2 (0.19) | 4 (0.15;0.16) | (0.01) | (0.00) | √ (0.40) | 4 (0.17;0.18) | 4 (0.14;0.15) | 1.38 (1.38 ± 0.13) | −0.95 (−0.75 ± 0.13) | |

L4 | Athabasca Basin | 1821–1891 | 2 (0.31) | 2 (0.24) | (0.05) | (0.00) | √ (0.88) | 3 (0.35;0.41) | 3 (0.26;0.26) | 1.16 (1.03 ± 0.22) | −0.72 (−0.41 ± 0.22) | |

L5† | Athabasca Basin | 1897–1934 | 4 (0.47;0.50) | 2 (0.20) | √ (0.44) | (0.00) | √ (0.20) | 2 (0.22) | 2 (0.18) | 1.23 (1.33 ± 0.26) | −0.67 (−0.74 ± 0.26) | |

L6 | West Central | 1821–1891 | 2 (0.15) | 2 (0.09) | √ (0.76) | √ (0.08) | √ (0.92) | 2 (0.10) | 2 (0.09) | 1.43 (1.46 ± 0.14) | −1.03 (−0.82 ± 0.14) | |

L7† | West Central | 1897–1934 | 4 (0.38;0.48) | 4 (0.23;0.30) | √ (0.32) | (0.01) | √ (0.11) | 2 (0.30) | 2 (0.32) | 0.92 (1.09 ± 0.34) | −0.45 (−0.45 ± 0.34) | |

L8 | Upper Saskatch. | 1821–1891 | 2 (0.21) | 3 (0.17;0.20) | √ (0.94) | √ (0.27) | √ (0.37) | 3 (0.19;0.21) | 3 (0.20;0.21) | 1.30 (1.32 ± 0.18) | −0.83 (−0.68 ± 0.18) | |

L9 | Winnipeg Basin | 1821–1891 | 3 (0.19;0/21) | 4 (0.14;0.16) | √ (0.86) | √ (0.98) | √ (0.31) | 4 (0.14;0.15) | 4 (0.15;0.16) | 1.42 (1.40 ± 0.15) | −0.91 (−0.78 ± 0.15) | |

L10 | North Central | 1821–1891 | 2 (0.35) | 2 (0.35) | (0.01) | (0.00) | √ (0.85) | 4 (0.39;0.44) | 3 (0.28;0.29) | 1.15 (1.02 ± 0.21) | −0.71 (−0.52 ± 0.21) | |

L11† | James Bay | 1895–1939 | 2 (0.20) | 2 (0.11) | √ (0.74) | √ (0.86) | √ (0.11) | 2 (0.12) | 3 (0.11;0.12) | 1.37 (1.47 ± 0.20) | −0.82 (−0.84 ± 0.20) | |

L12† | Lakes | 1897–1939 | 2 (0.26) | 2 (0.20) | √ (0.41) | √ (0.21) | (0.03) | 2 (0.18) | 2 (0.20) | 1.26 (1.34 ± 0.26) | −0.79 (−0.56 ± 0.26) | |

L13† | James B.+Lakes | 1897–1939 | 2 (0.20) | 2 (0.13) | √ (0.76) | √ (0.49) | √ (0.65) | 2 (0.13) | 2 (0.14) | 1.50 (1.46 ± 0.20) | −0.85 (−0.76 ± 0.20) | |

L14† | Gulf | 1897–1939 | 2 (0.55) | 2 (0.39) | √ (0.80) | √ (0.60) | √ (0.20) | 2 (0.37) | 2 (0.40) | 1.02 (1.05 ± 0.29) | −0.73 (−0.47 ± 0.29) |

Dimension refers to the optimal order using CV, adopting the Nadaraya–Watson kernel (N-W), the locally linear model (LL-reg), the additive model (GAM), or the linear autoregression model (Linear). The corresponding CV value are given in parentheses. The CV value for dimension two (and three) for all lynx (and hare) series are given as the second number in the parentheses whenever not optimal (not boldface). Two methods were used when testing for nonlinearity [H-T (21); Tsay (22)]. For both methods,

*P*values are given in parentheses. Additivity gives*P*values from a Lagrange multiplier test for additivity (23). The estimated partial derivatives of the various*f*functions (see text) at equilibrium of the skeleton model is given under lag1, lag2, and lag3. The corresponding linear autoregressive parameters (±2 SEM) are given in parentheses. No detrending has been carried out in spite of possible trends in some of the time series. Check marks indicate tests for which the null hypothesis was_{i}*not*rejected. Boldface indicates order estimates corresponding to the structure and complexity hypothesized in the main text of the paper.^{*}

CV value in parenthesis. First number gives the CV for the optimal dimension; second number (if present) gives the CV for dimension three for hare, and dimension two for lynx. The degrees of freedom for each lag is restricted to be 0 (the lag is not present), 1 (linear), 2 (nonlinear), or 4 (strongly nonlinear).

^{†}

Series has been interpolated for the missing observation in year 1914.

^{‡}

The L3 series in the linear case is selected to have dimension 11 based on the AIC (24) criterion when there is no upper limit for the possible dimensions to be selected. This is exceptional compared with the other lynx series.

## The Structure of the Time Series

### Statistical Modeling.

Data were log-transformed to stabilize the variance (25). This transformation is biologically suitable due to the multiplicative nature of the population dynamics process involving birth and death processes (26). The log-transformed series {

*X*} were scaled to have zero mean and variance equal to one._{t}In our search for parsimonious models for the time series under study, we start with a general model and simplify it as far as statistically permissible. Throughout we rely on recent developments in time series analysis (refs. 21, 23, 27–32; see references therein for details).

We follow the long standing conjecture (24, 33–37) that the transformed series may adequately be modeled in delay coordinates as a general autoregressive model of order [or dimension;

*sensu*Royama (38)]*d*; that is,*X*=_{t}*F*(_{d}*X*_{t−1},*X*_{t−2}*, . . . , X*_{t−d}) + ɛ_{t}. We assume that {ɛ_{t}} is time- and state-independent white noise (30). As suggested by Cheng and Tong (27), the function*F*determining the laws of the population dynamics can be estimated for a given choice of_{d}*d*with minimal assumptions, using nonparametric regression. Provided*F*_{d}is continuous, Cheng and Tong (28, 29) showed that the sample size requirement is exponential in*d*for estimation of the functional form of*F*_{d}but is at most quadratic in*d*for the estimation of the dimension. To find the appropriate dimension for the series (and their underlying process), we employ cross-validation (39–42). That is, for each*d*considered (*d*= 1, . . . , 4), we estimate*F*based on all but one data point, which is subsequently predicted. This is repeated for each data point. The sum of squared out-of-sample prediction errors [the cross-validation (CV) value; a rough estimate of the percentage noise is used as a measure of the model suitability. The_{d}*d*corresponding to the autoregressive model minimizing the CV value provide a consistent estimate of*d*, under minimal assumptions (25–28). To estimate the functions*F*, we employ two different nonparametric regression methods: (_{d}*i*) the Nadaraya–Watson kernel regression (21, 27, 28) and (*ii*) a locally linear regression (31, 32, 43). We employ a Gaussian product kernel for both. The consensus arising (Table 1) is that the snowshoe hare may parsimoniously be described by a model of dimension three (note that H2 appears cleaner in terms of predictability, based on its CV value, than H1), whereas the series for the lynx is more likely to be governed by a process of dimension two.§ The advantage of employing flexible models for*F*is that conclusions reached are not biased by any parametric prejudices. The disadvantage is that such general models have a high number of degrees of freedom (ref. 44, chapter 3) causing loss in precision (high statistical variance). These models may also be victims of the curse of dimensionality (45). Indeed, the statistical uncertainty associated with any estimate may be paramount when the number of data points is relatively low. The high degree of consistency across the different data sets (Table 1) is, however, encouraging._{d}To increase precision we investigate the possibilities of imposing some restrictive assumptions. One such simplifying assumption involve assuming that there is no interaction, on the logarithmic scale, between years (44, 46):On the basis of the test for additivity by Chen

\[ \begin{equation*}X_{t}=F_{d}({\cdot})=f_{{\mathrm{l}}}(X_{t-1})+{\cdots}+f_{d}{\vert}X_{t-d})+{\varepsilon}_{t}.\end{equation*}\]

[1]

*et al.*(23), this restrictive assumption appears permissible for 15 of 16 time series (Table 1); the rejection level of the deviating series (L12) is not far from the nominal 0.05 level. Additivity may, of course, be an artifact of low power due to the small sample size. However, the conclusion is biologically reasonable (see below) and the sample size is, after, all not particularly low. Estimating the optimal dimension by using CV of the additive model (Table 1) yields conclusions consistent with the analyses based on the more general models. Note, specifically that the additive models are usually as good as the more complicated models for prediction (see CV values in Table 1).Since we have demonstrated approximate additivity, a linear model (on the logarithmic scale) represents the next level of simplification:Two tests of nonlinearity were employed: (

\[ \begin{equation*}X_{t}={\alpha}_{{\mathrm{l}}}X_{t-1}+{\alpha}_{2}X_{t-2}+{\alpha}_{3}X_{t-3}+{\cdots}+{\alpha}_{d}X_{t-d}+{\varepsilon}_{t}.\end{equation*}\]

[2]

*i*) The nonparametric test by Hjellvik and Tjøstheim (21) and (*ii*) the Tukey one-degree-of-freedom test by Tsay (22). The null hypothesis of linearity is rejected in 5 of the 14 lynx series and is close to a nominal 0.05 level for one of the hare series (H1). This (31% at a nominal 0.05 level and 50% at a 0.1 level) is more than expected under linearity (assuming the different series to be independent). However, it is less than one may expect under extreme nonlinearity (see ref. 47). The estimates of*d*based on cross-validation of the linear models are generally higher than that based on the nonlinear models, as expected for nonlinear systems (48). Thus, it appears that the hare series may be represented an additive process on a logarithmic scale of order around three, while the order for the lynx series is around two. The lynx is certainly a nonlinear process (as previously concluded for series L3; see, for example, refs. 24 and 49). The same seems to be the case for the hare. Note, however, that the nonlinear function appear monotonic and not highly curved (Fig. 3).Figure 3

Despite significant nonlinearities, the monotoniety of the

*f*_{i}functions justifies the linear autoregressive model (Eq. 2) as a useful average of the state dependent (functional) autoregressive coefficients of the dynamics, because these may be interpreted as coefficients of statistical density dependence (13, 47, 51).Both the log-linear and the log-additive model indicate negative direct density dependence for the hare (α

_{1}-1; see refs. 13 and 51) ranging from −0.38 to −0.05 with a mean of −0.23. There is essentially no dependency on the second lag (α_{2}ranging from 0.06 to 0.11 with a mean of 0.09), but a strong negative dependency on the third lag is observed (α_{3}ranging from −0.52 to −0.24 with a mean of −0.38) (Table 1 and Fig. 3*A*). The three-dimensional structure of the hare system is consistent with earlier theoretical arguments (17, 52, 53) and with experimental results (4). For the lynx, the direct density dependency is positive (α_{1}-1 ranging from 0.01 to 0.48 with a mean of 0.26); the lagged dependency is strongly negative (α_{2}ranging from −0.84 to −0.27 with a mean of −0.63) (Table 1 and Fig. 3*B*). The two-dimensional structure for the lynx (L3) has been deduced in several earlier studies including Moran’s (ref. 33; but see §).All additive skeleton models (with ɛ

_{t}≡ 0 for all values of*t*) shown in Table 1 give rise to dampened oscillations (as judged by the Schur–Cohn stability criterion for linear difference equations; refs. 54 and 55). Thus, we do not observe limit cycles as reported by Tong (36) [based on the SETAR(2;2,2) model for the L3 lynx series]. The dampened oscillations will, however, be sustained in the presence of environmental stochasticity. The difference between limit cycles and weakly dampened cycles are, therefore, not so conspicuous in stochastic systems (see, e.g., ref. 47).### Ecological Interpretations.

Recent studies have shown that the hares feed on a variety of food plants and are eaten by an array of predators (Fig. 1

*A*) (56–59). The lynx is a specialist on hares but may utilize other prey species (Fig. 2*A*) (60, 61). When prey are scarce some predators may act as top predators on the lynx (62).The present statistical results are consistent with the three-level trophic model for the hare (4, 56, 58). Other models involving spatial structure or age–structure are also possible. Due to recent field experiments and observations (e.g., ref. 4), we emphasis the trophic hypothesis as the more plausible scenario. In contrast to the more complex regulation of the hare, the lynx dynamics are thought to be food-driven (60, 63–65). Thus, the classic view of the lynx–hare cycle as a simple and symmetric predator–prey interaction seems to be an oversimplification (refs. 13 and 64; see also below).

Because several food species are eaten by the hare and the hare itself is eaten by many predators (Fig. 1

*B*), the deduced three-dimensional structure of the time series lead us to hypothesize that the food species act as a group from the hare’s point of view and that the predators act as another group; if not, we might expect higher dimensionality. There is evidence from the Kluane Ecosystem Project that the predator community is strongly affected by intraguild predation, particularly as hares decline in abundance (62). Additional circumstantial evidence for this hypothesized compensation exists: Until about 1980, lynx were extensively trapped in northern Canada, but this predator removal has had no reported gross effect on the hare cycle (C.J.K., unpublished observations). During the 1986–1994 hare cycle, lynx were a major predator at Kluane. During this same period, northern goshawk numbers at Kluane were below normal (66). Finally, on Anticosti Island in the St. Lawrence River, there are no lynx present, yet the hare cycle persists (3). Both red fox and great-horned owls are particularly abundant on Anticosti Island and appear to compensate for the absence of lynx (2, 3). Keith (ref. 3, p. 115) documents cases of the extermination of lynx in southern Canada with continued hare cycles.The vegetation appears to segregate into palatable and nonpalatable species for the hare (56, 57). Among the palatable species (primarily

*Betula glandulosa*,*Salix glauca*, and*Picea glauca*), the hare has a mixed diet (66). Preferred foods are typically reduced during the peak of the hare cycle, whereas less preferred (but still palatable) species remain common (59). The impact of hares on their food plants is transient (58) and the recovery of the preferred food plants is substantial even before the decline is complete. At Kluane essentially no effect on the herbaceous vegetation of excluding hares from a 4-ha plot for 8 years has been found (R. Turkington, personal communication).Snowshoe hares constitute a major part of the lynx’s diet. Lynx may themselves be eaten by several top predators (Fig. 1

*A*; refs. 60 and 65 and M. O’Donoghue, personal communication). Studies in North America suggest, however, the hare to be the one major factor influencing the population dynamics of the lynx (60, 63, 67–69). The statistical models for the time series support the view of the dynamics of the lynx as a hare–lynx interaction.Through our analyses we have simplified the bewildering complexity portrayed in Figs. 1

*A*and 2*A*: If we are to study the dynamics of the lynx, we need to focus on the hare—in addition to studying the lynx. If we are to study the dynamics of the snowshoe hare, we need to focus on the guild of predators as a unit and the vegetation as a unit—in addition to the hare. The analyses have further documented detailed patterns of statistical density-dependence (Figs. 1*B*and 2*B*and Table 1). A viable hypothesis for the dynamics should be able to account for these patterns as well as for the dimensionality. In the following section, we suggest two simple mathematical models for trophic interactions within the boreal ecosystem. We investigate which constraints are required on the ecological interactions to generate the observed patterns of statistical density dependence. The detailed dynamic behavior of these models are deemed outside the scope of this report.## Two Ecological Models

### A Vegetation–Hare Predator Model.

Let where ɛ

*H*be the abundance of hares at time_{t}*t*,*V*_{t}be the abundance of the vegetation at time*t*, and*P*_{t}be the abundance of predators at time*t*; notice that*P*_{t}does not consist of lynx only but the combination of a variety of predators preying upon the hare. The functions*F*_{h},*F*_{v}, and*F*_{p}describe*per capita*ecological interactions such that\[ \begin{equation*}V_{t+1}=V_{t}F_{{\mathrm{v}}}(V_{{\mathrm{p}}},\hspace{.167em}H_{{\mathrm{p}}},\hspace{.167em}{\varepsilon}_{{\mathrm{v}}}),\end{equation*}\]

[3a]

\[ \begin{equation*}H_{t+1}=H_{t}F_{{\mathrm{h}}}(V_{{\mathrm{p}}},\hspace{.167em}H_{{\mathrm{p}}},\hspace{.167em}P_{{\mathrm{p}}},\hspace{.167em}{\varepsilon}_{{\mathrm{h}}}),\end{equation*}\]

[3b]

\[ \begin{equation*}P_{t+1}=P_{t}F_{{\mathrm{p}}}(H_{{\mathrm{p}}},\hspace{.167em}P_{{\mathrm{p}}},\hspace{.167em}{\varepsilon}_{{\mathrm{p}}}),\end{equation*}\]

[3c]

_{v}, ɛ_{h}, and ɛ_{p}represent stochastic influences. There is no general agreement about the impact of hares on edible vegetation (see also above): Wolff (56) and Keith (57) have measured stronger impacts on vegetation at the peak of the cycle than measured in the Kluane Ecosystem Project (R. Turkington, unpublished results). However, rapid vegetation recovery after the cyclic peak have been observed in all studies (56–58). The effect of the hare on the vegetation seems by and large negligible, thus allowing us to assume ∂*F*_{v}/∂*h*≈ 0 (*h*= ln*H*). The total predator community is strongly affected by the hare cycle, and the abundance of all the major hare predators—lynx, coyote, great-horned owl, and northern goshawk—follow the hare cycle closely (58, 70), suggesting ∂*F*_{p}/∂*h*> 0 and ∂*F*_{h}/∂*p*< 0 (where*p*= ln*P*).Approximating the where α

*F*functions in Eq.**3**by the first terms in a Taylor expansion in log-transformed abundances, we may, under minimal additional assumptions, write the log-linear model in delay coordinates (see Eq. 2):\[ \begin{equation*}h_{t}={\alpha}_{{\mathrm{l}}}h_{t-1}+{\alpha}_{2}h_{t-2}+{\alpha}_{3}h_{t-3}+{\varepsilon}_{{\mathrm{t}}},\end{equation*}\]

[4]

_{1}, α_{2}, and α_{3}are statistical parameters that may be written as functionals of the ecological system portrayed by Eq.**3**.¶ Interpreting the statistical results on the basis of Eqs.**3**and 4, the observed patterns of direct and delayed density dependence (α_{1}> 0, α_{2}≈ 0, and α_{3}< 0) may be shown to be highly plausible: given internal regulation through intraguild predation within the predator guild, (∂*F*_{p}/∂*p*≤ 0), the condition α_{3}< 0 is satisfied by assuming fairly weak self-regulation in both the vegetation and the hare population (∂*F*_{v}/∂*v*and ∂*F*_{h}/∂*h*are slightly negative, as documented in refs. 71–74). If ∂*F*_{p}/∂*p*is too negative, the relation α_{1}> 0 is violated: our statistical results, thus, introduces restrictions on the strength of the intraguild predator interaction—it must be strong, but not too strong. The relation α_{2}≈ 0 is fulfilled under a wide range of interaction strengths.### A Hare–Lynx Model.

Let Lynx occupy discrete territories during most of the hare cycle (63, 65, 67, 68) and only at the low phase (when hares are scarce) does the territorial system break down as many lynx go nomadic. A few lynx individuals survive the low phase by diversifying their diet with red squirrels and rodents. Hares represent the main food, however, even in the low phase (refs. 65, 67, and 68; M. O’Donoghue, personal communication; B. Slough, personal communication).

*H*′_{t}be the abundance of hares (with the possible inclusion of other herbivores preyed upon by the lynx;*H*′_{t}may be different from*H*used in the hare model) at time_{t}*t*and let*L*_{t}the abundance of lynx at time*t*. The functions*G*_{h}and*G*_{l}describe the ecological interactions such that\[ \begin{equation*}H^{\prime}_{t+1}=H^{\prime}_{t}G_{{\mathrm{h}}}(H^{\prime}_{{\mathrm{p}}}L_{{\mathrm{p}}}{\varepsilon}_{{\mathrm{h}}}),\end{equation*}\]

[5a]

\[ \begin{equation*}L_{t+1}=L_{t}G_{{\mathrm{l}}}(H^{\prime}_{{\mathrm{p}}}L_{{\mathrm{p}}}{\varepsilon}_{{\mathrm{l}}}).\end{equation*}\]

[5b]

Approximating the

*G*functions in Eq.**5**by the first terms in a Taylor expansion in log-transformed abundances, we may write:*l*= β_{t}_{1}*l*_{t−1}+ β_{2}*l*_{t−2}+ ɛ_{t}, where β_{1}and*β*_{2}are statistical parameters that may be written as functionals of the ecological system portrayed by Eq.**5**.**Due to the tight trophic interactions between the lynx and the hare (i.e., ∂

*G*_{h}/∂*l*< 0 and ∂*G*_{l}/∂*h*> 0), the empirical patterns of statistical density dependence seen in the parsimonious time series model for the lynx (Table 1; β_{1}> 0 and β_{2}< 0), are easily fulfilled in the predator–prey system if self-regulation within both the lynx and the hare populations are not too strong (the hare model requiring similar weak self-regulation in the hare).## Conclusion

Despite the complexity of the boreal food web, the realized dynamics of the snowshoe hare and the Canadian lynx are found to be of low dimension. We have furthermore found an asymmetry in the way the lynx and the hare are positioned within the ecosystem: The snowshoe hare appear to be regulated from below and above (by a variety of predators including the lynx). The lynx, in contrast, seems to be regulated only from below, and primarily by the hare. Thus, from the hare’s point of view, the food chain is a vegetation–hare–predator chain, whereas from the lynx point of view, the hare–lynx interaction dominates. The snowshoe hare certainly plays an interesting role in the ecological theatre of the Canadian boreal forest ecosystem.

## ABBREVIATION

- CV
- cross-validation

## Notes

**

The parameters in the autoregressive model for the lynx are given as follows: β

_{1}= 2 + ∂*G*_{h}/∂*h*+ ∂*G*_{l}/∂*l*and β_{2}= ∂*G*_{h}/∂*l*⋅∂*G*_{l}/∂*h*− (1 + ∂*G*_{l}/∂*l*)⋅(1 + ∂*G*_{h}/∂*h*).§

Cheng and Tong (27) found order three for the McKenzie River series (L3; Table 1). However, the fit of order two and order three are very similar in our analysis (Table 1) as well as in ref. 27. We are nevertheless convinced by the convergence on dimension two when all lynx series are seen together (Table 1).

¶

In the general case, the parameters of Eq. 4 may be given as follows: α

_{1}= 3 + ∂*F*_{v}/∂*v*+ ∂*F*_{h}/∂*h*+ ∂*F*_{p}/∂*p*, α_{2}= −3 − 2⋅∂*F*_{v}/∂*v*− 2⋅∂*F*_{h}/∂*h*− ∂*F*_{v/}∂*v*∂*F*_{h}/∂*h*+ ∂*F*_{v}/∂*h*⋅∂*F*_{h}/∂*v*− 2⋅∂*F*_{p}/∂*p*− ∂*F*_{v}/∂*v*∂*F*_{p}/∂*p*− ∂*F*_{h}/∂*h*⋅∂*F*_{p}/∂*p*+ ∂*F*_{h}/∂*p*⋅∂*F*_{p}/∂*h*, and α_{3}= (1 + ∂*F*_{v}/∂*v*)⋅(1 + ∂*F*_{h}/∂*h*)⋅(1 + ∂*F*_{p}/∂*p*) − (1 + ∂*F*_{v}/∂*v*)⋅∂*F*_{h}/∂*p*⋅∂*F*_{p}/∂*h*− (1 + ∂*F*_{p}/∂*p*)⋅∂*F*_{v}/∂*h*⋅∂*F*_{h}/∂*v*.## Acknowledgments

Comments by Rudy Boonstra, Robert M. May, George Sugihara, Howell Tong, and two anonymous referees on earlier versions of the paper were very helpful. Grants to N.C.S. and O.N.B. from the Norwegian Research Council and the University of Oslo and to C.J.K. from the Natural Sciences and Engineering Research Council of Canada are acknowledged.

## References

1

C S Elton

*Voles, Mice and Lemmings*(Clarendon, Oxford, U.K., 1942).2

C S Elton, M Nicholson

*J Anim Ecol***11**, 215–244 (1942).3

L B Keith

*Wildlife’s Ten-Year Cycle*(Univ. of Wisconsin Press, Madison, 1963).4

C J Krebs, S Boutin, R Boonstra, A R E Sinclair, J N M Smith, M R T Dale, K Martin, R Turkington

*Science***269**, 1112–1115 (1995).5

C J Krebs

*Ecology: The Experimental Analysis of Distribution and Abundance*(Harper Collins, 4th Ed., New York, 1994).6

R L Smith

*Ecology and Field Biology*(Harper Collins, 5th Ed., New York, 1996).7

A J Lotka

*Elements of Physical Biology*(Williams and Wilkins, Baltimore, 1925).8

V Volterra

*Nature (London)***118**, 558–560 (1926).9

M L Rosenzweig, R H MacArthur

*Am Nat***97**, 209–223 (1963).10

R M May

*Science***177**, 900–902 (1972).11

M E Gilpin

*Am Nat***107**, 727–730 (1973).12

W M Schaffer

*Am Nat***124**, 798–820 (1984).13

T Royama

*Analytical Population Dynamics*(Chapman & Hall, London, 1992).14

R H MacArthur

*Geographic Ecology*(Harper & Row, New York, 1972).15

R M May

*Stability and Complexity in Model Ecosystems*(Princeton Univ. Press, Princeton, NJ, 1973).16

S A Boutin, C J Krebs, R Boonstra, M R T Dale, S J Hannon, et al.

*Oikos***74**, 69–80 (1995).17

D A MacLulich

*Fluctuations in the Number of the Varying Hare (Lepus americanus)*(Univ. of Toronto Press, Toronto, 1937).18

H Chitty

*J Anim Ecol***19**, 15–20 (1950).19

C H Smith

*Can Field-Nat***97**, 151–160 (1983).20

L C Cole

*J Wildl Manage***18**, 2–24 (1954).21

V Hjellvik, D Tjøstheim

*Biometrika***82**, 351–368 (1995).22

R S Tsay

*Biometrika***73**, 461–466 (1986).23

R Chen, J S Liu, R S Tsay

*Biometrika***82**, 369–383 (1995).24

H Tong

*Non-Linear Time Series: A Dynamical System Approach*(Clarendon, Oxford, U.K., 1990).25

A Sen, M Srivastava

*Regression Analysis: Theory, Methods and Applications*(Springer, New York, 1990).26

M Williamson

*The Analysis of Biological Populations*(Edward Arnold, London, 1972).27

B Cheng, H Tong

*J R Stat Soc B***54**, 427–449 (1992).28

B Cheng, H Tong

*Developments in Time Series Analysis*, ed T T Subba Rao (Chapman & Hall, London), pp. 183–206 (1993).29

B Cheng, H Tong

*Philos Trans R Soc London A***348**, 325–341 (1994).30

K S Chan, H Tong

*J R Stat Soc B***56**, 301–311 (1994).31

Q Yao, H Tong

*J R Stat Soc B***56**, 701–725 (1994).32

J Fan, Q Yao, H Tong

*Biometrika***83**, 189–206 (1996).33

P A P Moran

*Aust J Zool***1**, 163–173 (1953).34

M J Campbell, A M Walker

*J R Stat Soc A***140**, 411–431 (1977).35

H Tong

*J R Stat Soc A***140**, 432–435, & 448–468. (1977).36

H Tong

*Threshold Models in Non-Linear Time Series Analysis, Lecture Notes in Statistics*(Springer, Heidelberg, No. 21. (1983).37

V Haggan, T Ozaki

*Biometrika***68**, 189–96 (1981).38

T Royama

*Ecol Monogr***47**, 1–35 (1977).39

C J Stone

*Ann Stat***5**, 595–645 (1977).40

G Sugihara, R M May

*Nature (London)***344**, 734–741 (1990).41

G Sugihara, B Grenfell, R M May

*Philos Trans R Soc London B***330**, 235–251 (1990).42

G Sugihara, W Allan, D Sobel, K D Allan

*Proc Natl Acad Sci USA***93**, 2608–2613 (1996).43

J Fan

*J Am Stat Assoc***87**, 998–1004 (1992).44

T Hastie, R Tibshirani

*Generalized Additive Models*(Chapman & Hall, London, 1990).45

D W Scott

*Multivariate Density Estimation: Theory, Practice, and Visualization*(Wiley, New York, 1992).46

R Chen, R S Tsay

*J Am Stat Assoc***88**, 298–308 (1993).47

N C Stenseth, O N Bjørnstad, T Saitoh

*Proc R Soc London B***263**, 1117–1126 (1996).48

D S Broomhead, R Jones

*Proc R Soc London A***423**, 103–121 (1989).49

K S Chan, H Tong

*J R Stat Soc B***52**, 469–476 (1990).50

*S-Plus for Windows Version 3.2 Supplement*(Statistical Sciences, Seattle, WA, 1993).

51

O N Bjørnstad, W Falck, N C Stenseth

*Proc R Soc London B***262**, 127–133 (1995).52

R M May

*Ecology***54**, 315–325 (1973).53

N C Stenseth

*Science***269**, 1061–1062 (1995).54

M Marden

*Geometry of Polynomials*(Am. Math. Soc., Providence, RI, 1966).55

R M May

*Stability and Complexity in Model Ecosystems*(Princeton Univ. Press, 2nd Ed., Princeton, NJ, 1975).56

J O Wolff

*Ecol Monogr***50**, 111–130 (1980).57

L B Keith

*Oikos***40**, 385–395 (1983).58

L B Keith

*Curr Mammal***2**, 119–195 (1990).59

A R E Sinclair, C J Krebs, J N M Smith, S Boutin

*J Anim Ecol***57**, 787–806 (1988).60

C H Nellis, S P Wetmore, L B Keith

*J Wildl Manage***36**, 320–329 (1972).61

L B Keith, A W Todd, C J Brand, R S Adamcik, D H Rusch

*Proc Int Congr Game Biol***13**, 151–175 (1977).62

M O’Donoghue, E Hofer, F I Doyle

*Nat Hist***104**, 6–9 (1995).63

R M R Ward, C J Krebs

*Can J Zool***63**, 2817–2824 (1985).64

H R Akcakaya

*Ecol Monogr***62**, 119–142 (1992).65

U Breitenmoser, B G Slough, C Greitenmoser-Wursten

*Oikos***66**, 551–554 (1993).66

R I Doyle, J N M Smith

*Stud Avian Biol***16**, 122–129 (1994).67

K G Poole

*Can J Zool***73**, 632–641 (1995a).68

K G Poole

*J Wildl Manage***58**, 608–618 (1995b).69

G M Koehler

*Can J Zool***68**, 845–851 (1990).70

D S Hik

*Wildlife Res***22**, 115–129 (1995).71

E John, R Turkington

*J Ecol***83**, 581–590 (1995).72

S A Boutin

*Oecologia***62**, 393–400 (1984a).73

S A Boutin

*J Anim Ecol***53**, 623–637 (1984b).74

S A Boutin, C J Krebs, A R E Sinclair, J N M Smith

*Can J Zool***64**, 606–610 (1986).## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

Copyright © 1997, The National Academy of Sciences of the USA.

#### Submission history

**Received**: May 30, 1996

**Accepted**: March 3, 1997

**Published online**: May 13, 1997

**Published in issue**: May 13, 1997

#### Keywords

#### Acknowledgments

Comments by Rudy Boonstra, Robert M. May, George Sugihara, Howell Tong, and two anonymous referees on earlier versions of the paper were very helpful. Grants to N.C.S. and O.N.B. from the Norwegian Research Council and the University of Oslo and to C.J.K. from the Natural Sciences and Engineering Research Council of Canada are acknowledged.

### Authors

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

#### Cited by

Loading...

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.