Chaos is not rare in natural ecosystems

Chaotic dynamics are thought to be rare in natural populations but this may be due to methodological and data limitations, rather than the inherent stability of ecosystems. Following extensive simulation testing, we applied multiple chaos detection methods to a global database of 172 population time series and found evidence for chaos in >30%. In contrast, fitting traditional one-dimensional models identified <10% as chaotic. Chaos was most prevalent among plankton and insects and least among birds and mammals. Lyapunov exponents declined with generation time and scaled as the −1/6 power of body mass among chaotic populations. These results demonstrate that chaos is not rare in natural populations, indicating that there may be intrinsic limits to ecological forecasting and cautioning against the use of steady-state approaches to conservation and management. Previous meta-analyses have found that few ecological time series are actually chaotic; here, the authors use recently developed chaos detection methods and a larger number of time series to show that 30% of studied time series are characterized as chaotic.

C haos was introduced to ecology nearly 50 years ago 1,2 to provide an explanation for widespread fluctuations in abundance of natural populations. The defining characteristics of chaos are bounded, deterministic, aperiodic dynamics that depend sensitively on initial conditions. If common, chaos would offer the promise of short-term predictability while setting hard limits on long-term forecasting 3 . It would also mean that the 'stable ecosystem' paradigm-the theoretical justification for linear statistical models of ecological dynamics 4 and steady-state management policies 5would need rethinking. Chaos has been observed in many ecological models [6][7][8][9] , demonstrated in laboratory experiments (for example, insects 10 , microbes 11 and plankton 12 ) and detected in a handful of well-studied field systems [13][14][15][16] . However, most meta-analyses assessing the prevalence of chaos in natural field populations have found chaos to be absent or rare [17][18][19][20][21] . The most recent global meta-analysis concluded that only 1 out of 634 ecological time series was chaotic 18 .
The apparent rarity of chaos in free-living natural populations is a mystery for several reasons. Ecosystems involve tens to thousands of species and large complex systems are prone to chaos 9,22,23 . Nonlinear dynamics, a necessary condition for chaos, are also common in ecological time series 24 and many abiotic drivers of population dynamics are themselves chaotic 25 . In light of this, we hypothesize that the dearth of evidence for ecological chaos reflects methodological and data limitations, rather than genuine rarity. Importantly, many meta-analyses in ecology have tested for chaos by fitting one-dimensional, parametric population models to time series [17][18][19][20] -models in which the current state depends only on the last state and this dependency is constrained to a particular functional form. We know from theory and focused empirical studies that overcompensation is not the only mechanism that can generate chaos and that chaos often arises through ecological interactions 16,26 . Although seminal models of ecological chaos were one-dimensional 1 , using one-dimensional models to classify natural populations treats ecological complexity (for example, species interactions) as noise, thereby hindering chaos detection 27,28 . Non-parametric, multidimensional methods for chaos detection that make minimal assumptions about the dynamics 27,29,30 are more mathematically robust and potentially more accurate, particularly in cases where the underlying dynamics are complex and not well understood. In contrast to the one-dimensional parametric studies, the last global meta-analysis to use flexible, higher-dimensional methods, published in 1995, found evidence for chaos in 11% of 27 field time series (excluding four series on measles cases), noting that this was probably an underestimate 30 . The question of chaos prevalence using comparable methods has not been revisited and, in the interim, many more time series of sufficient length have become available-a critical factor for detecting chaos 31 . In addition, new chaos detection tools are also now available 32-35 but how these methods, developed outside ecology, will perform on ecological time series is unknown.
Here, we revisit whether chaos is, in fact, rare in ecological systems using a suite of flexible, higher-dimensional approaches. The definitive and most widely used index of chaos is the Lyapunov exponent (LE), which measures the average rate of divergence between nearby points in phase space 36 (Supplementary Note 2). Positive LE values are indicative of chaotic dynamics. We selected two methods of estimating LEs (direct 37 and Jacobian 31 ) and four additional chaos detection algorithms (recurrence quantification analysis 32 , permutation entropy 33 , horizontal visibility graphs 34 and the chaos decision tree 35 ).
We tested the six methods on data simulated with a variety of chaotic, periodic and stochastic models to benchmark misclassification rates under ecologically relevant time series lengths and levels of noise (Supplementary Notes 2-4). We tested the generality of this classification accuracy using two additional suites of simulation models. Three methods had error rates >0.5 and so were not pursued further (Table 1 and Extended Data Fig. 1). This included the direct LE method, which was unable to differentiate divergence due to chaos from divergence due to noise 38 . We applied the remaining three methods (Jacobian, recurrence quantification and permutation entropy) to time series from the Global Population Dynamics Database (GPDD) 39 . The GPDD aggregates 4,471 time series from 1,891 taxa. Previous analyses of the GPDD have concluded that many of these time series are too noisy to permit accurate modelling 40,41 . Therefore, we restricted our attention to the subset of the GPDD where chaos could be detected if present; that is, relatively long time series of good data quality without any major gaps (Methods). Applying these criteria produced a dataset of 172 time series representing 138 different taxa from 57 distinct locations with between 30 and 197 observations. To confirm the prevalence Chaos is not rare in natural ecosystems Tanya L. Rogers 1 ✉ , Bethany J. Johnson 2 and Stephan B. Munch 1,2,3 ✉ Chaotic dynamics are thought to be rare in natural populations but this may be due to methodological and data limitations, rather than the inherent stability of ecosystems. Following extensive simulation testing, we applied multiple chaos detection methods to a global database of 172 population time series and found evidence for chaos in >30%. In contrast, fitting traditional one-dimensional models identified <10% as chaotic. Chaos was most prevalent among plankton and insects and least among birds and mammals. Lyapunov exponents declined with generation time and scaled as the −1/6 power of body mass among chaotic populations. These results demonstrate that chaos is not rare in natural populations, indicating that there may be intrinsic limits to ecological forecasting and cautioning against the use of steady-state approaches to conservation and management.
of chaos among plankton in the GPDD, which were sampled largely from one location, we also analysed an independent dataset of 34 zooplankton time series from three lakes with between 138 and 639 observations.
We then explored how Jacobian LE values varied among taxa and depended on intrinsic timescale (generation time), body size (mass), time series length (generations sampled) and embedding dimension (E) which we define here as the number of lags needed to reconstruct the dynamics (Supplementary Note 1). The Jacobian method estimates LEs from a local linear model fit to lags of the time series 42,43 .

Results and discussion
Across three independent classification methods, at least one-third of the GPDD time series were classified as chaotic ( Table 1). The most conservative estimate (34%) was obtained with the Jacobian LE method, which in our simulations had the best performance on short time series, was the most robust to process noise and underestimated the frequency of chaos in the presence of substantial observation error (Supplementary Note 4). We focus most of our remaining analyses on the Jacobian LE estimates.
Noise and non-stationarity can affect the classification of time series and, on the basis of tests with low-dimensional parametric models, these were found to be present in the GPDD 40,44 . However, the time series we selected for chaos detection only partially overlap these previous studies. Hence, we needed to address the role that noise and non-stationarity play in our specific results. If noisy time series were being incorrectly classified as chaotic, we would expect a higher frequency of chaos among series with lower prediction accuracy. However, the fraction classified as chaotic by the Jacobian method did not vary with prediction R 2 (logistic regression, n = 172, X 2 d.f.=1 = 0.006, P = 0.9; Extended Data Fig. 2a) and series with high prediction error did not have higher LEs (Extended Data Fig. 2b). The frequency of chaos (34%) also did not change if only series with prediction R 2 > 0.25 were considered. So, although chaotic series were more variable than non-chaotic series (Fig. 1a), they were actually somewhat more predictable (Fig. 1b); hence observation error is not inflating the frequency of chaos.
If non-stationarity in the mean was driving the results, we should expect chaotic series to exhibit strong monotonic trends, exponential growth or nearly linear dynamics. Only six time series that had either strong monotonic trends and/or near-linear dynamics were misclassified as chaotic (Fig. 1c,d). Reclassifying these series (four birds, one mammal and one insect) as not chaotic reduced the frequency of chaos to 30%. The majority of chaotic series, however, were strongly nonlinear (Fig. 1c), did not display a strong monotonic trend ( Fig. 1d and Extended Data Fig. 2c,d) and had a median growth rate near 0. Although these metrics do not capture more subtle forms of non-stationarity, they suggest that non-stationarity and exponential growth are not responsible, by and large, for the observed frequency of chaos.
These observations in the GPDD time series are consistent with our simulations and previously published results 45 : the Jacobian method was less likely to find chaos as observation noise increased, was minimally affected by process noise, rarely classified long-term trends as chaotic and effectively discriminated between chaos and stochastic linear dynamics with seasonality (Supplementary Figs. [1][2][3][4][5]. Taken together, these analyses indicate that the frequency of ecological chaos is not an artefact. So, why is chaos more prevalent in our study than in previous meta-analyses? Whereas the methods used here make minimal assumptions about the dynamics, most earlier analyses classified series by fitting one-dimensional population models [17][18][19][20] . To evaluate the effect of these constraints on chaos detection, we first restricted the Jacobian method to E = 1, essentially fitting a one-dimensional non-parametric model. This reduced the apparent frequency of chaos in the GPDD from 34% to 9.9%, with reductions seen across all taxonomic groups (Fig. 2). Changes in classification were most common among populations in which the optimal E was high ( Supplementary Fig. 6), consistent with the hypothesis that reducing dimensionality inhibits chaos detection. When we fit a set of one-dimensional parametric models used in previous meta-analyses, this further reduced the apparent frequency of chaos to 6% or less (Supplementary Note 5 and Extended Data Fig. 3). Thus, the one-dimensional parametric assumption used in other meta-analyses probably explains the rarity of chaos detection in these studies.
Data limitation might also account for differences with other analyses, many of which used much shorter time series. To address this, we re-evaluated the prevalence of chaos in the 106 time series with at least 50 data points; 42% were chaotic. Restricting further to the 57 series with at least 70 data points increased the prevalence of chaos to 58%. Overlong sampling intervals might also bias our results because time series with sampling intervals larger than the Lyapunov horizon (timesteps on the order of LE −1 ) should appear effectively stochastic, producing false negatives. Of the 30 series with at least 70 data points and less than four generations per timestep (none of which were plankton), 40% were classified as chaotic. Hence, chaos seems to be more visible in longer series, but long sampling intervals do not appear to inflate the prevalence of chaos. Approaching this from the other direction, we found that 42% of the chaotic time series were no longer classified as chaotic when truncated to 30 data points (the minimum length used in our analysis). These results are consistent with our simulations results where data limitation increased the false negative rate much more than the false positive rate (Supplementary Note 4 and Supplementary Figs. 1-5). Thus, with longer time series we expect to see a greater fraction of populations classified as chaotic. These results probably explain why other meta-analyses using very short time series (<20 data points) found no evidence for chaos 21 .
Having allayed most reasonable qualms about statistical artefacts, we further explored the biological contexts in which chaos occurs. The frequency of chaos differed among taxonomic groups; phytoplankton had the greatest proportion of chaotic series (81%), followed by zooplankton (77%), insects (43%), bony fishes (29%), birds (18%) and mammals (16%) (Fig. 2). The prevalence of chaos decreased in species with longer generation times (logistic regression, n = 166, X 2 d.f.=1 = 26.7, P < 0.001; Fig. 3a) which tended to have lower LEs as well (Fig. 3b). E also tended to decrease with increasing generation time (Pearson r = −0.30) and was lowest among birds ( Supplementary Fig. 7). There are several possible explanations for this pattern. Long-lived species definitionally have lower average mortality rates. Hence, on a per unit time basis (but not per generation), we might expect long-lived species to have relatively weaker interactions with other species, leading to lower LE and E compared with short-lived taxa. Long-lived species may also be better insulated from chaotic environmental drivers 19,46 . Data limitations might also account for lower rates of chaos in long-lived species, as chaos detection depends on the time series length relative to the intrinsic timescale for the system. If having fewer generations sampled reduces the ability to detect chaos, we would expect data truncation to have a larger effect on species with longer generation times. There was a trend where species with longer generation times were more likely to be reclassified as non-chaotic when series were truncated to 30 data points; however, this result was not statistically significant (logistic regression, n = 49, X 2 d.f.=1 = 0.79, P = 0.38; Supplementary Fig. 8).
Recent evidence indicates that LEs scale with body mass in experimental demonstrations of chaos 47 . To determine how generally this applies to natural populations, we evaluated whether variation in LEs among chaotic species exhibited analogous mass (M) scaling. Combining LEs from our study and existing studies [10][11][12]16,[48][49][50][51]  suggesting that differences in mass account for most of the average differences in LE among broad taxonomic groups. To account for potential non-independence of LEs within locations, we repeated the regression of log 10 (LE) on log 10 (M) with general least squares and again using only location means. The resulting slopes differed by no more than 0.02. The consistency of LE scaling between laboratory and natural populations cannot readily be explained as a statistical artefact. Moreover, since the laboratory-derived LEs are not artefacts of observation noise or non-stationarity, consistent scaling with the field data provides additional evidence for chaos in natural populations.
The 172 GPDD series come from 57 distinct locations, raising the potential issue of non-independence (for example, 33 of 34 plankton and 15 of 17 fish time series each came from one location; Extended Data Fig. 4a). In a well-mixed system where all variables interact on similar timescales, there will be a unique maximum LE and the LEs estimated from each state variable should be similar, as has been observed in some laboratory studies 12 . However, in systems with weak coupling, modularization or separation of timescalesprocesses known to occur in natural ecosystems 52,53 -different LEs can be reconstructed from different variables, representing a range of dynamics corresponding to different subsystems. In fact, we hypothesize that this weak coupling and timescale separation contribute to the observed mass scaling of the LEs. The number of independent LEs at a given location will be somewhere between one and the number of species sampled. The exact number in our database is not resolvable, however, as only three locations had more than eight time series (Extended Data Fig. 4a) and tools to resolve these subsystems have yet to be developed and tested (an important next step). As a coarse first pass, examining the distribution of LEs suggests bimodality in two of the three locations with more than eight series (Extended Data Fig. 5); however, this distribution cannot be resolved for the other locations. Thus, we have presented our results at the time series level but recognize that further resolution is required. For the time being, we can conduct a conservative test of system-level chaos by assuming that each location represents one well-mixed system and asking whether the median LE for each location is chaotic (>50% of LEs are significantly >0 with 95% confidence). By this standard, 21% (12/57) of locations were chaotic and the prevalence of chaotic locations was 41% (5/12) in insects, 25% (4/16) in birds and 9% (2/22) in mammals (plankton and fish were not sufficiently well represented to compute this frequency). We can also use the error rates from our simulation study to compute the probability that a system is chaotic, accounting for the known asymmetry between false positive and false negative rates (with the important caveat that these error rates are dependent on the particular suite of simulations used). By this standard, 25% of locations had a >50% probability of being chaotic and 19% of locations had a >80% probability of being chaotic (Extended Data Fig. 4b). Both of these location-level results are consistent with our time-series level results and our finding that chaos is not rare.
At the time series level, our estimate of 34% is greater than that of the last meta-analysis to use higher-dimensional models and Jacobian LEs 30 (11%) and the difference may simply reflect an absence of plankton in their database. As most plankton in the GPDD are from a relatively open marine system, it is plausible that what appears here as chaos reflects advection of patchily distributed populations. To address this, and the fact that all zooplankton series came from a single location, we evaluated the frequency of chaos and mass scaling of LEs in additional time series of zooplankton from three lakes. The prevalence of chaos for the lake zooplankton time series was 47% and all lakes had a >60% location-level probability of chaos. Among the chaotic taxa, none of these new LE estimates was significantly different from the mass scaling derived from the GPDD (deviation of observed values from regression predictions, P > 0.05 for all; Fig. 4). It is exceedingly unlikely that Colour distinguishes broad taxonomic groups. Includes data from this study (gPDD and supplemental results from three lake systems) and positive LEs compiled by ref. 47 (Ag2020). the log-log scale is in keeping with previous studies 47 . Note that the lakes data (squares) were not used to fit the regression line. Supplementary Fig. 9 shows the same points with lower confidence intervals.
advection would result in LEs that scale with mass consistently across three datasets, although it may contribute to the frequency of chaos.

Conclusions
Single species models are routinely used to evaluate population status in applied fields such as fisheries 5 and conservation biology 54 . However, our results clearly show that scalar population models typically mischaracterize dynamics, treating complexity as noise and leading to the conclusion that chaos is rare [17][18][19]55 . As noted by Robert May, such models "do great violence to reality" 56 . More flexible methods (for example, refs. 30,43,45 ) are better able to characterize complex dynamics and integrating these into population status assessments is an important area for future research.
Reflecting on the frequency of chaos in natural populations, we note that birds and mammals, the least chaotic taxa, make up 58% of the time series we analysed but represent <1% of the species on earth 57 . Thus, chaos may be considerably more common than the one-third presented here. Diseases, genetic variants, species and statistical events are labelled 'rare' using thresholds ranging from 0.001% to 5%. By these standards, chaos in natural ecosystems is far from rare. This presents both challenges and opportunities for ecology as a predictive science; although short-term forecasting is feasible 58 , precise long-term prediction is likely to be impossible and management should avoid defining objectives in terms of equilibrium conditions. However, with increasing amounts of data and modern learning algorithms, new frontiers are open for characterizing the complex, non-equilibrium and high-dimensional dynamics of ecology which will advance both our understanding of natural variability and improve our ability to manage ecosystems.

Methods
Data. We obtained abundance time series data from the GPDD 39 accessed through the R package 'rgpdd' 59 . Our analyses required reasonably long and continuous time series and for organisms to be detected with sufficient frequency to reconstruct their dynamics. Consequently, we selected series with a reliability score of at least 2, at least 30 non-missing time points, at least 5 unique abundance values, <60% zeros and <22% missing time points (in our dataset, this resulted in time series having no more than 11 missing values). We used only field-collected survey data (we excluded laboratory and harvest data), excluded human diseases and excluded the shorter and lower-quality of six duplicate time series that passed our filtering.
Before analysis, all untransformed abundance time series were rescaled to unit variance by dividing by the standard deviation. This transformation was not strictly necessary but aided in visualization and diagnostics. To allow for log transformations and calculations of population growth rate, ln (xt/xt−τ ), all time series containing zeros were rescaled after adding a constant (1 if all values were integers, the minimum non-zero value if the series contained non-integers). Leaving the zeros intact and using only model forms that did not require log transformations produced similar results.
As a measure of organismal intrinsic timescale, generation time was obtained from published sources for all species in our dataset. We used the age at first reproduction as a proxy for generation time, unless direct estimates of generation time or doubling time were available. Wet body mass was obtained from published sources or, if unavailable, was estimated from volume, assuming that organisms have the same density as water. Generation time and mass data were not included for seven taxa that were not finely resolved enough taxonomically to obtain this information. Sources for generation time and mass were the following: birds, mammals, fish, insects from ref. 60 ; diatoms from refs. 61,62 ; insect masses not included in ref. 60 from refs. [63][64][65][66][67][68] ; copepods from refs. [69][70][71] ; and dinoflagellates from refs. 72,73 .
The plankton data in the GPDD are nearly all marine and thus from relatively open systems. Hence, it is possible that their dynamics reflect water movement in addition to population growth. However, these time series display seasonal peaks and troughs that persist for months, rather than the more ephemeral fluctuations expected from water mass movement and it seems reasonable to assume that these represent population dynamics over a large spatial area as opposed to fluid dynamics. Nevertheless, to assess the robustness of these plankton results, we performed a supplemental analysis on 34 monthly zooplankton time series data from three lake systems which are arguably more 'closed' than the marine environment. These systems were Lake Zurich (Wasserversorgung Zürich), Lake Geneva (SOERE OLA-IS, AnaEE-France, INRA of Thonon-les-Bains, CIPEL, 19 December 2019, developed by the Eco-Informatique ORE system of the INRA 74 ) and Oneida Lake 75 . We also required these series to have <60% zeros. Mass data for these species were obtained from refs. 63,[76][77][78][79] . If only dry mass was available, we assumed dry mass was 20% of wet mass for arthropods and 4% of wet mass for rotifers 76,80 .
Analysis. Our goal was to use a combination of modern and classical methods for detecting chaos to characterize ecological time series. However, most of these methods were developed in data-rich fields and tested on finely spaced time series with thousands to millions of observations. Therefore, we began by testing six chaos detection methods on simulated data from 37 stochastic, periodic and chaotic models with ecologically relevant time series lengths and levels of observation and process noise. These simulations included both a test set of models used for tuning and new sets of models used for validation and evaluation of generality. The specific classification methods we tested were the 'direct' method of estimating Lyapunov exponents (DLE) 37 , the Jacobian method of estimating Lyapunov exponents (JLE) 31,43 , recurrence quantification analysis (RQA) 32,81 , permutation entropy (PE) 33 , the horizontal visibility algorithm (HVA) 34,82 and the chaos decision tree (CDT) 35 . Note that the more traditional DLE and JLE have been tested previously 30,31,38,45,83 and we re-test them here for comparison with the newer methods. Supplementary Note 1 provides a brief background on time-delay embedding and a comparison of methods for selecting the embedding dimension and time delay which are used in many of the detection methods. Supplementary Note 2 provides the mathematical definition for LE and full details on our implementation of each detection method. Supplementary Note 3 provides details on the simulation models and Supplementary Note 4 summarizes results of the simulation testing.
Under the conditions of our simulations, DLE, HVA and CDT had either false positive or false negative rates >0.5 in both the test and validation datasets and so were not pursued further. We applied the remaining methods, which all had false positive rates <0.2, to the empirical dataset to estimate the frequency of chaos in natural populations.
The JLE method derives LEs from the Jacobian matrices of a local linear time-delay embedding model (Supplementary Note 2.2). We explored several methods for generating confidence intervals for the LE and selected the method that produced the best classification accuracy in the simulated data (Supplementary Table 2). The JLE method proved to be the most accurate index of chaos in the simulated test and validation datasets with the lowest false positive rate.
Note that LEs based on short time series, such as those used in this analysis, are more a reflection of the local, rather than global LE. While we may not be able to determine if a system, in the long run, is chaotic, these estimates characterize how the system behaved over this particular time period. Moreover, although short time series may not generate a precise estimate of the LE, our simulations indicate that, for the models tested, the sign of the LE can be accurately determined using relatively low sample sizes.
Since the LE also is the most widely used index of chaos and provides a quantitative, scale-invariant measure of divergence rate, we used the numeric values of the LEs to further explore the relationship between chaotic dynamics, intrinsic timescale, body size, time series length and embedding dimension in the GPDD dataset, after converting the LE from units of timestep −1 to units of month −1 . To calculate the mass scaling of the LE, we followed previous work 47, 84 and used ordinary least squares regression on log 10 transformed data for LE > 0.
To test the effect of time series length and sampling interval on the inferred LE (and subsequent classification), we examined the proportion of time series classified as chaotic if we were to restrict our analysis to only those series with at least 50 or 70 observations and to only those with at least 70 observations and at least four generations per timestep. We also truncated all time series that had been classified as chaotic to the last 30 observations and recomputed the LE. To account for the fact that some of the time series were sampled from the same locations and thus may (but not necessarily) represent series from the same dynamical system, we also computed the proportion of locations with a median LE that was chaotic (at least 50% of time series were chaotic with 95% confidence). To account for known asymmetry in the false positive and false negative rates, we also computed the probability that a location was chaotic using the error rates for the JLE method in our simulation study; however, these results should be interpreted with caution since we cannot know how representiative the error rates are of ecological reality. The corrected probability of chaos was calculated as (P − FPR)/(TPR − FPR), where P is observed proportion of time series classified as chaotic, FPR is false positive rate and TPR is true positive rate.
To test whether non-stationarity or long-term trends affected our results, we examined whether LEs were greater in series with stronger monotonic trends.
We assessed the degree of monotonic trend using the squared Spearman rank correlation between abundance and time. To test whether the restriction of dimensionality affects the inferred LE, we recomputed the LE with the embedding dimension set to 1. To test whether the restriction of model form, in addition to restricting dimensionality, affects the inferred LE, we fit a set of common one-dimensional population models with the form xt+1 = xtexp [f (xt, q)] to each time series and used the fitted model to estimate the LE (Supplementary Note 5).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The GPDD data are available on KNB with identifier https://doi.org/10.5063/ F1BZ63Z8. Zooplankton data were obtained for Oneida Lake from KNB (identifier kgordon.17.67), for Lake Zurich from Wasserversorgung Zürich and for Lake Geneva from the Observatory on LAkes (OLA-IS, AnaEE-France, INRA of Thonon-les-Bains, CIPEL; https://doi.org/10.4081/jlimnol.2020.1944). The simulated datasets and generating code are available in the code repository. The specific GPDD time series used and associated metadata (including compiled generation time and mass data) are available in the code repository.