Introduction

Ecological and epidemiological cycles have long been the subject of interest and controversy (Soper 1929; Elton and Nicholson 1942; Levy and Wood 1992; Ricker 1997; Trenberth 1997; Kendall et al. 1998; Grover et al. 2000; Krebs et al. 2001; Myers and Cory 2013). The underlying mechanisms causing such cycles are often hard to identify, in particular when the only available data is in form of a recorded time series. Commonly, suspected mechanisms are noise-sustained oscillations (Royama 1992; Kaitala et al. 1996; Nisbet and Gurney 2004), limit cycles arising from intrinsic ecological interactions (May 1972; Kendall et al. 1999; Gilg et al. 2003), and cyclic environmental forcing (Elton 1924; London and Yorke 1973; Sinclair et al. 1993; Hunter and Price 1998; García-Comas et al. 2011).

Noise-sustained oscillations can in principle be distinguished from intrinsic limit cycles and forced oscillations based on the system’s probability distribution (Pineda-Krch et al. 2007), which describes how much time system trajectories spend in different areas of phase space. Discerning between limit cycles and forced oscillations is much harder, partly because the system’s probability distribution is similar in both cases, namely concentrated along the cycle. Conventional approaches include autoregression analysis (Royama 1992), fitting generic, or system-specific mechanistic models to data (Turchin and Taylor 1992; Kendall et al. 1999), as well as cross-correlation and cross-spectral analysis of time series with hypothesized environmental driving forces (Sinclair et al. 1993; Stenseth et al. 2002). The first two methods implicitly assume an intrinsic mechanistic origin of the observed cycles, and it is not always clear at which point this assumption should be dropped in favor of an explanation involving periodic forcing. The third method requires a concrete suspected driving force, to which the focal time series is to be compared.

Investigation of ecological cycles often involves their power spectrum (PS) or autocovariance (ACV), both standard concepts in time series analysis (Platt and Denman 1975; Royama 1992). The general perception is that forced cycles are characterized by a non-decaying autocovariance, i.e., a long-term temporal coherence, in contrast to intrinsic limit cycles (Berryman 2002; Nisbet and Gurney 2004). But this criterion implicitly assumes that any potential driving force would be temporally coherent or phase-remembering in Nisbet and Gurney's (2004) terminology. This is the case for the seasonal cycle, but is certainly not the general rule. For example, climatic cycles such as the El Niño-Southern Oscillation (Burgers 1999; Stenseth et al. 2002), the solar cycle (Sinclair et al. 1993; Klvana et al. 2004), and population cycles possibly causing fluctuations of coexisting focal species (Bulmer 1974) are typically phase forgetting, that is, exhibit low temporal coherence. Ecological variables driven by such cycles will clearly also be phase forgetting. Therefore, classifying ecological cycles as limit cycles solely based on a decaying autocovariance can lead to wrong conclusions.

In this article, we argue that a time series-based distinction between intrinsic limit cycles (LC) and forced oscillations (FO) is in principle possible when the focal variable is subject to random perturbations. The reason is that LCs and FOs respond qualitatively differently to noise, which shapes their power spectrum and autocovariance in characteristic ways. To the best of our knowledge, these qualitative but subtle differences have not been pointed out before in the context of distinguishing between LCs and FOs. Based on our observations, we propose a novel approach for classifying cyclic ecological variables as either part of a limit cycle or simply forced by one. We assume that other possibilities like noise-sustained oscillations and quasi-periodic chaos have been ruled out, and that the time series shows a single dominant spectral peak. We propose a set of unitless summary variables, V PS and V ACV, that are to be extracted from the estimated PS and ACV, and based on which a distinction between LC and FO becomes possible. We demonstrate their predictive power using numerical simulations of simple generic models for limit cycles and forced oscillations.

Theoretical background

Perturbations should reveal periodic forcing

In the following, we identify a periodic force or environmental parameter with the temporal variation it causes to the position of a stable node of the forced system. For example, the hypothetical effects of sockeye salmon population cycles (Ricker 1997) on zooplankton communities shall be seen as a periodic change in zooplankton equilibrium densities due to varying predation pressure.

The methods that we propose for classifying cyclic ecological variables as either part of intrinsic cycles or forced oscillations are based on the fundamentally different effects that random perturbations have on cycles in these two cases. Perturbations of forced oscillations are reversed by the forced system, which tends to stabilize on the equilibrium currently defined by the external driving force (Fig. 1a). Thus, perturbations of the focal ecological variable do not feed back into the driving force, and in particular do not influence the cycle’s phase. In contrast, perturbations of limit cycles can cause random phase shifts along the cycle that are not reversed by the system’s dynamics (Fig. 1b). The resulting stochastic phase drift leads to an increased spectral bandwidth and a decay of the system’s autocovariance at larger time lags. These effects, known as jitter in electronic signal theory (Demir et al. 2000), result in the temporal decoherence of the cycle and a deviation from true periodicity. It is on these grounds that temporal decoherence has been perceived as a signature of limit cycles, as opposed to driven oscillations (Berryman 2002; Nisbet and Gurney 2004). However, as already pointed out, forced oscillations can also be decoherent when driven by decoherent cyclic forces. Hence, temporal decoherence by itself is generally not a sufficiently informative property.

Fig. 1
figure 1

On the fundamentally different effects of perturbations (first half of thick arrow) on forced oscillations (a) and intrinsic limit cycles (b). In the former case, perturbations tend to be reversed by the intrinsic dynamics of the forced system. In the latter case, perturbations have long-term effects on the cycle’s phase

Conceptually, deciding whether the focal ecological variable is part of a limit cycle or not translates to deciding whether the cycle’s phase is causally affected by changes of the variable or not. Provided that the latter is subject to random perturbations, this question is equivalent to determining the extent to which these perturbations contribute to the cycle’s decoherence. Stable systems driven by a cyclic force will typically exhibit power spectra and autocovariances that relate to the ones of the driving force by means of their so called frequency transfer function and response function, respectively (Schetzen 2003; Nisbet and Gurney 2004). However, as discussed above, additional random perturbations are expected to modify their power spectra and autocovariances in characteristic ways that differ from perturbations that cause decoherence in the underlying cycle. It is these signatures that we attempt to quantify in the following section.

Considered models

As a limit cycle model, we consider the normal form of a supercritical Hopf bifurcation (Guckenheimer and Holmes 1985, Section 3.4), perturbed by additive Gaussian white noise. The Hopf bifurcation describes the emergence of a limit cycle around an unstable focus. It appears in a wide range of ecological and epidemiological models (Rosenzweig and MacArthur 1963; Greenhalgh 1997; Fussmann et al. 2000; Pujo-Menjouet and Mackey 2004), which is why we chose it as a generic model structure. We concentrate on a single model variable, denoted by L, because we wish to draw an analogy to ecological time series, which are often only available for a small subset of system variables, such as the population density of a single species. More precisely, we let

$$\begin{array}{@{}rcl@{}} \frac{d}{dt}\left(\begin{array}{l} L\\ y\end{array}\right) &=& \left(\begin{array}{ll}\lambda_{\mathrm{l}}/2 & -\alpha_{o}\\ \alpha_{o} & \lambda_{\mathrm{l}}/2\end{array}\right)\cdot\left(\begin{array}{l} L\\ y\end{array}\right)\notag\\ && + \frac{(L^2+y^2)}{A^2} \left(\begin{array}{l}-\lambda_{\mathrm{l}} L/2 - (\alpha-\alpha_{o}) y\\-\lambda_{\mathrm{l}} y/2 + (\alpha-\alpha_{o}) L\end{array}\right)\notag\\ && + \sigma_{\mathrm{l}} \left(\begin{array}{l}\varepsilon_1\\\varepsilon_2\end{array}\right), \end{array} $$
(1)

where y is a second variable coupled to L that is needed for the complete description of a limit cycle. The first row in Eq. 1 corresponds to the linear dynamics near the unstable focus, while the second row captures the nonlinearities giving rise to the limit cycle. Here, α>0, A>0, and λ l < 0 are the limit cycle’s angular frequency, amplitude, and negative Lyapunov exponent (or resilience), respectively, while α o >0 is the system’s angular frequency in the proximity of the focus. ε 1 and ε 2 are independent standard Gaussian white noise processes, and σ l>0 is the noise amplitude. σ l quantifies the random perturbations of the cycle, such as weather events or demographic stochasticity in predator–prey cycles.

The power spectrum (PS) and autocovariance (ACV) of L can be approximated analytically when fluctuations around the limit cycle are weak (\(\sigma _{\mathrm {l}}^2\ll A^2\lambda _{\mathrm {l}}\)) (Louca 2013, work in review). At any angular frequency ω, the PS is given by

$$\begin{array}{@{}rcl@{}} PS[L](\omega) &\approx& A^{2} F\left[\left(\frac{\sigma_{\mathrm{l}}}{A}\right)^{2},\alpha,\omega\right]\notag\\ &&+\frac{\sigma_{\mathrm{l}}^2}{2\lambda_{\mathrm{l}}}F\left[{\left(\frac{\sigma_{\mathrm{l}}}{A}\right)^{2} + 2\lambda_{\mathrm{l}},\alpha,\omega}\right], \end{array} $$
(2)

where we defined

$$F[s,\alpha,\omega]:=\frac{2s[{4(\alpha^2+\omega^2) + s^2}]}{[{4(\alpha-\omega)^2 + s^2}][{4(\alpha+\omega)^2 + s^2}]}. $$

The autocovariance of L at time lag τ is approximately

$$\begin{array}{@{}rcl@{}} \text{ACV}[L](\tau) &\approx& \frac{1}{2}{\left[A^{2} + \frac{\sigma_{\mathrm{l}}^2}{2\lambda_{\mathrm{l}}}e^{-\lambda_{\mathrm{l}} |{\tau}|}\right]}\cos(\alpha \tau)\notag\\ &&\times \exp\left[{-\frac{|{\tau}|\sigma_{\mathrm{l}}^{2}}{2A^{2}}}\right]. \end{array} $$
(3)

The exponential decay expressed in the second row of Eq. 3 quantifies the temporal decoherence of the cycle due to a Brownian motion-like phase drift along the deterministic trajectory. Figures 2a and b show the power spectrum and autocovariance estimated from a typical times series generated by such a noisy limit cycle, along with the expected values in Eqs. 2 and 3, respectively.

Fig. 2
figure 2

a and b Power spectrum and autocovariance, respectively, estimated from time series generated using model (Eq. 1) for a noisy limit cycle of frequency α. c and d Power spectrum and autocovariance, respectively, estimated from time series using model (Eq. 4) for a stable linear system forced by a noisy limit cycle. Black dashed curves represent theoretical formulas (Eqs. 2 and 5) for the power spectra and the autocovariances, respectively. Parameter values given in Online Resource 1, Table 2

A forced ecological variable will respond approximately linearly to the variation of its equilibrium induced by the driving force, if this variation is sufficiently weak. We consider the simple case of a one dimensional linear system with resilience λ s whose equilibrium is, after appropriate rescaling, identified with L:

$$ \frac{dX}{dt} = -\lambda_{\mathrm{s}}(X - L(t)) + \sigma_{\mathrm{s}} \varepsilon. $$
(4)

Hence, X is driven by the cyclic variable L(t), which is part of a perturbed limit cycle (Eq. 1) with variable perturbation strength. The amplitude σ s of the additive Gaussian white noise in Eq. 4 quantifies the strength of non-decohering perturbations, such as the effects of short-term weather events on seasonally driven phytoplankton populations. In the deterministic limit (σ l=σ s=0), both the limit cycle in Eq. 1 and the driven system in Eq. 4 generate similar signals. In the stochastic case, it can be shown that the process given in Eq. 4 has the power spectrum

$$ \text{PS}[X](\omega) = \frac{\lambda_{\mathrm{s}}^2}{\omega^2+\lambda_{\textbf{s}}^2}\cdot \text{PS}[L](\omega) + \frac{\sigma_{\mathrm{s}}^2}{\omega^2+\lambda_{\mathrm{s}}^2} $$
(5)

and autocovariance

$$\begin{array}{@{}rcl@{}} \text{ACV}[X](\tau) &=& \frac{\lambda_{\mathrm{s}}}{2}\int_{\mathbb{R}} \text{ACV}[L](u)\cdot e^{-\lambda_{\mathrm{s}}|{u-\tau}|}\ du\notag\\ &&+ \frac{\sigma_{\mathrm{s}}^{2}}{2\lambda_{\mathrm{s}}}e^{-\lambda_{\mathrm{s}}|{\tau}|} \end{array} $$
(6)

(Online Resource 1, Section 1). The modulation of the force’s spectrum in Eq. 5 is due to the finite resilience λ s of the forced system, which dampens high-frequency variations in L. Similarly, the convolution of ACV[L] in Eq. 6 expresses the memory of the forced system due to its finite resilience. For promptly responding systems (λ s), the first terms in Eqs. 5 and 6 become identical to PS[L] and ACV [L], respectively. In the general case, the integral in Eq. 6 becomes a complicated function of the model parameters (given in Online Resource 1, Section 1.3).

We are particularly interested in the right-most terms in Eqs. 5 and 6, originating in the non-decohering perturbations. These terms represent promising signatures of forced oscillations. Figures 2c and d show the power spectrum and autocovariance estimated from a typical time series generated by the model in Eq. 4, together with their theoretical values. The increased power at low frequencies seen in Fig. 2c, as well as the upward shift of the autocovariance at small time lags in Fig. 2d, are characteristic of forced oscillations, subject to non-decohering perturbations (Nisbet and Gurney 2004).

Classifier variables

The formulas for the expected power spectrum Eq. 5 and autocovariance Eq. 6 allow, at least in principle, an estimation of A,α, σ l,λ l, σ s, and λ s from a cyclic time series, for example, through ordinary least squares fitting to the estimated power spectrum and autocovariance. Ideally, the right-most additive terms in Eqs. 5 and 6 should be estimated close to zero if the focal variable is part of an intrinsic limit cycle, and far from zero if the cyclic pattern results from cyclic forcing. In view of these expectations, we propose two summary variables of cyclic ecological time series:

$$V_{\text{PS}}=\frac{\sigma_{\mathrm{s}}^2}{\lambda_{\mathrm{s}}^2PS(0)},\quad V_{\text{ACV}}=\frac{\sigma_{\mathrm{s}}^2}{2\lambda_{\mathrm{s}}\text{ACV}(0)}. $$

V PS, estimated for example by fitting Eq. 5, relates the power due to non-decohering perturbations to the total process power at zero frequency. V ACV, estimated by fitting Eq. 6, relates the variance originating in non-decohering perturbations to the total variance. Both V PS and V ACV are expected to be distributed at higher values (≲1) in the case of forced oscillations (FO), and at lower values (≈0) for limit cycles (LC). Figure 3a, b show the estimated power spectrum and autocovariance of an El Niño-Southern Oscillation index (Trenberth 1997), along with fitted versions of Eqs. 5 and 6. The same is shown in Fig. 3 c, d for a summer temperature profile of Vancouver, clearly forced by the diurnal cycle. We chose these two textbook examples because their true nature is unambiguous. As seen in the figures, the two time series give values for V PS and V ACV compatible with our predictions.

Fig. 3
figure 3

a Power spectrum and b autocovariance of the El Niño N3.4 index between January 1871 and December 2007. The low values for V PS and V ACV (given in the figure legends) are in accordance with the oscillatory character of sea surface temperature, as part of the El Niño Southern Oscillation. c Power spectrum and d autocovariance of Vancouver temperature during July and August, 2013. The high values for V PS and V ACV reflect the forced character of the diurnal temperature oscillation. The dashed curves in all plots show the fitted formulas for the power spectrum in Eq. 5 and autocovariance in Eq. 6, from which V PS and V ACV were obtained. a, b were estimated from monthly data provided by NCAR CAS National Centre for Atmospheric Research (2007). c, d were estimated from hourly data provided by Environment Canada (2013)

We argue that V PS and V ACV are promising classifier variables that can be used in standard binary classification schemes for discerning limit cycles from forced oscillations (McLachlan 2004; Fradkin and Muchnik 2006). In practice, the distributions of V PS and V ACV for the two cases, LC and FO, are expected to overlap as a result of finite time series and estimation errors, therefore compromising the accuracy of any classifier based on them.

Numerical evaluation of V PS and V ACV

Methods

Using Monte Carlo simulations (Kroese et al. 2011), we evaluated the suitability of the summary variables V PS and V ACV for discerning between time series generated by the limit cycle model in Eq. 1 and the forced system in Eq. 4. Specifically, we examined the behavior of a simple threshold-based classification scheme based on any one of V∈{V PS,V ACV}, defined as follows: Depending on the value of V, estimated from some given time series, we

  • classify the time series as a LC if VT l (where T l is some fixed threshold between 0 and 1),

  • classify the time series as a FO if V>T u (where T u is some fixed threshold between T l and 1), and

  • reject the case as undetermined if T l<VT u.

Similar univariate classification schemes are widespread in diagnostic decision making (Dudoit et al. 2002; Pepe 2004).

The false FO detection rate (FDR) of such a classifier is defined as the fraction of LCs misclassified as FOs. Its true FO detection rate (TDR) is the fraction of FOs correctly classified as FOs. Its accuracy is the fraction of non-rejected cases that are correctly classified as either LC or FO. Varying the thresholds T l and T u changes the FDR, TDR, the rejection rate, and the classification accuracy. For example, a zero rejection rate is only achieved if T u=T l. In general, the FDR, TDR, rejection rate, and classification accuracy depend on the assumed probability of encountering a LC or FO. For our investigations, we took a non-biased approach and assumed both LC and FO to be equally likely. In principle, if the distributions of V PS and V ACV are known for LCs and FOs, the classification accuracy can be maximized by choice of T l and T u, whose difference should be sufficiently small for the rejection rate not to exceed a certain acceptable value.

Even if the rejection rate is zero (T l=T u), a tradeoff exists between reducing the FDR and increasing the TDR. The relationship between the FDR and TDR, as the classification threshold is varied, can be visualized using the so called receiver operating characteristic (ROC) (Brown and Davis 2006). The ROC curve is a widespread tool for the evaluation of threshold-based classifiers in signal detection theory and diagnostic decision making (Zweig and Campbell 1993). The area under the ROC curve (AUC), which ranges between 0.5 and 1.0, is often taken as a simple measure of the discriminative power of the classification variable, independently of any particular threshold. An AUC =0.5 corresponds to no discriminative power at all and an AUC >0.8 is generally considered excellent (Hosmer 2013).

We sampled a large part of the entire 7-dimensional parameter space for the LC model and the FO model, generating an equal number of LC and FO time series of fixed size and duration. The considered parameter ranges are summarized in Online Resource 1, Table 1. We discarded cases that passed the D’Agostino-Pearson K 2 normality test (D’Agostino et al. 1990) with a significance above 0.05. This test evaluates the deviation of the phase-space distribution of trajectories from the normal distribution. A significant deviation from normality has been proposed as a criterion for ruling out noise-sustained oscillations as the cause of cycles (Pineda-Krch et al. 2007). We only considered cases whose estimated spectral peak had a statistical significance P<0.05 (Scargle 1982; Horne and Baliunas 1986). For each run, we calculated V PS and V ACV by least squares fitting of Eqs. 5 and 6. We estimated the probability distributions of V PS and V ACV based on 1,000 runs. Using either V PS or V ACV, we calculated the maximum achievable accuracy of the classification scheme described above, by suitable choice of the thresholds T l and T u and for various rejection rates. We also constructed the ROC curves for the non-rejecting classifier (T l=T u) using either V PS or V ACV. Technical details are given in Online Resource 1, Section 2.

Results

Our simulations verified our predictions on the separate probability distributions of V PS (or V ACV) for LCs and FOs (Fig. 4a–c). The classification accuracy itself strongly depends on the quality (resolution and length) of the time series. For high-quality data (200 points across 20 cycle periods), the classification accuracy can reach almost 100 % for any of the two variables, given a sufficiently high rejection rate and an optimal choice of the classification thresholds T l,T u. For example, at a zero rejection rate, the maximum achievable accuracy exceeds 80 % for both V PS and V ACV. Accordingly, the ROC curves for the non-rejecting classifier exhibit a high AUC (>0.85 for both V PS and V ACV, Fig. 5a), underlining the strong discriminative power of V PS and V ACV. At a 30 % rejection rate, the maximum accuracy even exceeds 90 %.

Fig. 4
figure 4

a Scatterplot of the summary variables (V PS,V ACV) obtained from Monte Carlo simulations of the limit cycle model (Eq. 1) (squares) and the forced oscillation model (Eq. 4) (crosses) as described in the main text. b Corresponding estimated probability distributions of V PS for limit cycles (continuous line) and forced oscillations (dashed line). c Corresponding estimated probability distributions of V ACV for limit cycles (continuous line) and forced oscillations (dashed line). ROC curves calculated using the distributions shown in (b) and (c) are shown in Fig. 5a. Time series consisted of 200 points spanning across 20 cycle periods

Fig. 5
figure 5

Receiver operating characteristics for the non-rejecting (i.e., single-threshold) LC/FO classifiers using V PS and V ACV, for high-quality (a) and low-quality (b) time series. The area under the curve (AUC) is given in the figure legends and is a measure of the discriminative power of the considered variable. Red circles mark points of maximum classifier accuracy, with corresponding classification thresholds given adjacently. Time series consisted of 200 points and spanned 20 cycle periods in (a) and consisted of 25 points spanning 5 cycles in (b)

For low-quality data (25 points across 5 cycle periods), the accuracy can drop by as much as 20 % and the AUC of the non-rejecting classifier can fall below 0.7 (Fig. 5b). The optimal classification thresholds were also found to slightly depend on the quality of the time series (Fig. 5), as well as on the considered range of model parameters. This is not surprising, since the probability distributions of V PS and V ACV depend not only on the distributions of model parameter values, but also on the accuracy of their estimation and therefore the quality of the time series.

It should be noted that power spectra estimated from time series are expected to show a flat background, should the sampling frequency be much lower than the system’s resilience. This becomes clear in expression Eq. 5 for the power spectrum. In that case, fitted parameter values (in particular σ l, σ s) exhibit a high variance, thus compromising the predictive power of V PS. Similar issues arise for V ACV when sampling durations are shorter than the typical correlation times of the limit cycle or of the forced system. Both of these effects have been observed in our Monte Carlo simulations.

Discussion

We have investigated the possibility of discerning intrinsic limit cycles from environmentally forced oscillations based on a single time series. Our proposed approach is superior to conventional procedures when neither an underlying limit cycle nor a periodic forcing can be excluded, and a concrete suspicion for a driving force as well as a mechanistic understanding of the studied system are both lacking. In our analysis, we did not require that cyclic environmental forces are temporally coherent, as is often (and often wrongfully) implicitly assumed. Hence, our reasoning also holds for cases where the periodic forcing itself results from a perturbed limit cycle exhibiting low temporal coherence.

We emphasize that care should be taken with the interpretation of signals classified as limit cycles. By construction, classifiers based on the summary variables V PS and V ACV can not differentiate between ecological variables undergoing forced oscillations but not subject to detectable perturbations, and dynamic variables bidirectionally coupled to a limit cycle. For example, the concentration of snowshoe hare droppings in northern Canada might be classified as part of a limit cycle, should a suitable time series be available (Krebs et al. 2001). However, this would be a result of the strong effects that hare populations have on dropping concentrations (i.e., droppings are a proxy for hares), and not because droppings themselves are part of a (putative) limit cycle.

In our models, we assumed that random perturbations are temporally uncorrelated. This assumption of white noise is only valid if environmental perturbations have short correlation times compared to the characteristic time scales of the studied ecosystem. It has been pointed out by several authors that environmental noise can show significant temporal correlations, particularly so in marine ecosystems (Steele 1985; Halley 1996; Halley and Inchausti 2004; Vasseur and Yodzis 2004). The differences in the power spectrum and autocovariance of limit cycles and forced oscillations, due to colored noise, remain to be investigated (see (Teramae et al. 2009) for some first results).

Our analysis is based on the assumption that quasi-periodic chaos has been ruled out as a possible cause of decoherent cyclicity. However, chaos is likely to occur in natural populations and has been proposed as an explanation for population-dynamical patterns like fluctuating vole populations (Turchin and Ellner 2000) and insect outbreaks (Dwyer et al. 2004). Identifying deterministic chaos in populations and distinguishing it from noise poses significant challenges and has long been subject to investigations (Hastings et al. 1993; Ellner and Turchin 1995; Cencini et al. 2000). Chaos is commonly detected in time series based on estimated Lyapunov exponents (Turchin and Ellner 2000; Kendall 2001; Becks et al. 2005; Beninca et al. 2008). Such non-parametric tests should precede the analysis proposed in this article when the system is suspected of being chaotic. Moreover, a poor fit of the theoretical power spectrum Eq. 5 and autocovariance Eq. 6 should always serve as an indicator of alternative mechanisms.

Complications are also expected for systems exhibiting multiple oscillatory mechanisms. For example, a combination of induced plant defenses and natural predation seems to explain forest insect outbreaks of alternating magnitude (Elderd et al. 2013) and snowshoe hare cycles likely result from an interaction of predation and food supplies (Krebs et al. 2001). Furthermore, environmental seasonality can interact with nonlinear intrinsic processes and a complicated phase-space structure to produce irregular fluctuations, as demonstrated by models explaining seasonal measles outbreaks (Keeling et al. 2001) or the interannual variability of plankton community structure (Dakos et al. 2009). Prior estimation of the system’s dimension (Abbott et al. 2009) can help avoid erroneous simplistic assumptions on the dynamics. Moreover, in many cases, different mechanisms leave characteristic signatures on the time series, such as separate peaks or higher harmonics in the power spectrum (Bjørnstad and Grenfell 2001; Hammer 2007; Elderd et al. 2013). Prior subtraction of secondary peaks, similar to the subtraction of long-term trends (Yamamura et al. 2006), might enable the analysis of isolated oscillatory mechanisms. However, we do not know how robust our analysis is against such data transformations. Ideally, the theory introduced here should be extended to these more complicated scenarios.

The anticipated practical value of V PS and V ACV in any classification scheme presumes a certain universality in their distribution for different systems, at least for systems whose limit cycles emerge through a supercritical Hopf bifurcation. In fact, the formulas (Eqs. 5 and 6) used to define V PS and V ACV are, strictly speaking, only meaningful for systems sharing a similar structure to Eqs. 1 and 4. We thus anticipate a low accuracy for cycles involving variables and perturbations at different scales of magnitude. More general versions of V PS and V ACV will be required to characterize such cases.

Regardless of the discriminative power of V PS and V ACV, choosing optimal classification thresholds is a non-trivial task that falls into the general category of classifier training (Alpaydin 2004). Thresholds determined from simulations of generic models, such as the ones presented here (Fig. 5), should serve as a starting point. Better results might be obtained through supervised learning using pre-classified time series from similar systems (e.g., using data from the Global Population Dynamics Database (NERC Centre for Population Biology 2010)), as is common in machine learning for medical diagnosis. For single-threshold classifiers, this translates to estimating the probability distributions of V PS and V ACV for typical ecological time series, based on which optimal thresholds can be determined.

Conclusions

Our exploratory work shows that the signatures of noise on ecological cycles can be used to distinguish variables contributing to the cyclic behavior from variables merely driven by it. The differences become particularly apparent in the summary variables V PS and V ACV, which quantify the stochasticity not directly affecting the cycle’s phase. Moreover, the described approach has the potential to reveal the stability properties of the focal system as well as the amount of perturbations it is subject to. In the absence of experimental manipulations, naturally occurring perturbations and the transients they induce can therefore provide useful insights into a cycle’s dynamics.