## Abstract

Populations occasionally experience abrupt changes, such as local extinctions, strong declines in abundance or transitions from stable dynamics to strongly irregular fluctuations. Although most of these changes have important ecological and at times economic implications, they remain notoriously difficult to detect in advance. Here, we study changes in the stability of populations under stress across a variety of transitions. Using a Ricker-type model, we simulate shifts from stable point equilibrium dynamics to cyclic and irregular boom–bust oscillations as well as abrupt shifts between alternative attractors. Our aim is to infer the loss of population stability before such shifts based on changes in nonlinearity of population dynamics. We measure nonlinearity by comparing forecast performance between linear and nonlinear models fitted on reconstructed attractors directly from observed time series. We compare nonlinearity to other suggested leading indicators of instability (variance and autocorrelation). We find that nonlinearity and variance increase in a similar way prior to the shifts. By contrast, autocorrelation is strongly affected by oscillations. Finally, we test these theoretical patterns in datasets of fisheries populations. Our results suggest that elevated nonlinearity could be used as an additional indicator to infer changes in the dynamics of populations under stress.

## 1. Introduction

Traditionally, ecosystem management has been based on mechanistic models that attempt to predict ecosystem dynamics. Unfortunately, however, these models are usually based on limited mechanistic knowledge and have been notoriously poor at prediction [1,2]. This shortcoming is evidenced by the fact that ever more elaborate models do not necessarily improve out-of-sample prediction of complex ecological dynamics [3,4]. For example, in fisheries, increasing the complexity of a mechanistic model by including variables that are believed to be influential, such as temperature, can actually reduce predictability [1]. Thus, there is an urgent need for improved methods when it comes to forecasting ecological responses [4–7] especially under current levels of anthropogenic stress [8].

No doubt, this is not a trivial task as the uncertainty in ecological responses to anthropogenic stresses poses a major challenge for ecological forecasting [4,7]. One main reason for that is that stress may induce radical shifts in ecological dynamics. Such shifts are mathematically related to the crossing of bifurcation points: thresholds where the stability properties and, thus, the behaviour of a system, changes fundamentally [9]. Although ecologists have a long history in studying bifurcations in models [10], it is difficult to empirically predict if and when such bifurcation points will be crossed. More importantly, crossing particular bifurcations can have serious ecological consequences when they involve catastrophic shifts between alternative states [11]. Examples of such catastrophic shifts (also termed critical transitions [12]) are the overgrowth of coral reefs by macroalgae [13,14], unexpected pest outbreaks like boom-and-bust cycles of spruce budworm beetles [15] or the abrupt shifts in the composition of pelagic marine communities [16,17]. Similar bifurcation points may also lie behind less dramatic transitions, like the occurrence of algal blooms in lakes due to eutrophication [18], population extinction in deteriorating environments [19] or the elimination of infectious diseases [20].

Coping with such fundamental sources of unpredictability requires alternative modelling strategies that go beyond traditional approaches [6,21]. One set of such strategies focuses on early warning signs that attempt to detect abrupt transitions based on the statistical properties of the observed time series without requiring any mechanistic model [12]. Such changes are caused by critical slowing down (CSD), a phenomenon by which stable systems close to local bifurcation points respond slowly to disturbances [22,23]. Variance and autocorrelation are the most studied statistical indicators that arise from CSD. Rising autocorrelation has been shown to be an indicator of the increasing risk of extinction in stable laboratory experiments with yeast [24] and zooplankton populations [19], whereas increasing variance can mark the shift of lake dominance from piscivorous to planktivorous fish [25,26].

While the CSD indicators are model-free, they mostly posit stable dynamics and transitions between stationary attractors. Parallel work has focused on forecasting non-equilibrium and chaotic population dynamics that arise mechanistically from unstable attractors based on an approach called empirical dynamical modelling (EDM) [1,27]. Similar to the statistical approach taken by CSD indicators, this approach is also equation-free. EDM is based on reconstructing an attractor manifold directly from time series (https://www.youtube.com/watch?v=rs3gYeZeJcw). This reconstruction (essentially, a multi-dimensional scatterplot of the points in the time series) can be done with multiple time series of interest (each axis as a time series) or with a single time series (each axis as a lag of the time series) [28]. In either case, how well the attractor is reconstructed can be verified by its ability to forecast future states. In fact, EDM has been shown to outperform equation-based approaches at forecasting recruitment in sockeye salmon populations [1], dynamics of Pacific sardines [29] and the fate of experimental flour beetle populations [30]. With its capacity to forecast the future state of a system, EDM could potentially be a useful alternative approach for detecting the proximity of a system to nearby bifurcations.

Theoretically, when a dynamical system is approaching a bifurcation (i.e. where attractors and/or dynamics fundamentally change), stochastic events may more easily push a system across attractor boundaries or into a region of state space where dynamics may be also affected by a different attractor. This implies that, close to a bifurcation, the realized dynamics become increasingly state-dependent. State-dependence means that the future evolution of a system is determined by the dynamics around its current state. EDM can evaluate state-dependence (or what we call nonlinearity) by comparing forecast performance obtained when using global versus local information to model the system [31,32]. Local information refers to the local neighbourhood of points on the attractor closest to the system's current state, while global information refers to all points on the attractor. If local information gives a better forecast of system state compared with global information, the behaviour of the system is deemed state-dependent or nonlinear. In principle, this concept can be generalized to a number of bifurcations: across stable equilibria, from stable states to cyclic attractors and to chaos. In that way, while CSD indicators rely mostly on changes in stability between stable attractors [33], EDM does not, so that in an unknown system, it is likely that the union of these two approaches may be more informative for anticipating critical transitions in ecosystems under stress.

Here, we study whether nonlinearity, a measure of state-dependence quantified by EDM, can provide warning of imminent shifts in the dynamics of populations under stress. We generate time series using a simple stochastic Ricker-type population model. We estimate nonlinearity and we compare it to CSD indicators (variance and autocorrelation) across different dynamical transitions. As an illustration, we test all indicators in stressed and non-stressed populations from two long-term fisheries datasets. Our aim is to explore the capacity of nonlinearity to be an indicator of loss of stability and to show how the combination of these model-free approaches can be helpful for detecting dynamical changes in populations under stress.

## 2. Material and methods

### 2.1. Simulated data

We used a discrete Ricker-type model that describes the logistic growth of population *N* with density-dependence and an extra loss term (e.g. due to exploitation). We chose this model not only because it has been generically used for describing population dynamics for a variety of organisms (fisheries, insects and birds, e.g. [34–36]), but also because it can exhibit different types of bifurcations that are frequently found in most ecological models. The model reads
2.1where *N _{t}* is population biomass,

*r*is the intrinsic growth rate,

_{t}*b*defines the strength of density-dependence (=

*r*/

_{t}*K*, where

*K*is the carrying capacity set by the environment (=10)), and exploitation follows a sigmoid functional response (

*p*= 2), with half-saturation

*h*(=0.75), and maximum harvesting rate

*F*. We assumed process error in the model to represent environmental stochasticity with a Gaussian term

*ɛ*of zero mean and

_{t}*σ*(=0.25) standard deviation. We also considered demographic stochasticity in the growth rate

_{E}*r*by using exponential filtering at each time step (), where

_{t}*r*

_{0}is the mean and

*σ*(=0.1) the standard deviation of the Gaussian noise term

_{r}*ɛ*

_{r}_{,t}[37].

We considered two scenarios where changes in external conditions can trigger changes in the dynamics of a population. In the first scenario, we hypothesized that anthropogenic stress may lead to changes in population demographic traits. For example, size-selective harvesting of large individuals in a population may cause small size individuals to mature at an earlier age [38,39]. Such earlier age-at-maturation can be associated with an overall increase in the intrinsic growth rates of a population and it has been suggested as an explanation for affecting the dynamics of overharvested fish populations [40]. Changes in abiotic conditions (e.g. climate warming), pesticide resistance or decreasing competition from resident populations can increase intrinsic growth rates of pests or alien species leading to outbreaks [41,42]. On the other hand, toxicity, habitat degradation or food depletion may lead to a decrease in the intrinsic growth rate of a population favouring the chances of stochastic extinction [43–45]. We mimic all these effects by gradually changing growth rate *r*_{0} (=[0.01, 3]), while setting the overall harvesting rate *F* to zero (figure 1*a*). We call this scenario the ‘demographic change scenario’. Changing growth rate in the Ricker model leads to a well-known series of transitions from stable equilibrium to cycles through a flip bifurcation, and eventually transitions to chaos through a series of period-doubling bifurcations [9,46]. At zero growth rate, the population goes to extinction through a transcritical bifurcation.

In the second scenario, we hypothesized that a population runs the risk of collapse due to direct exploitation pressure [47]. We labelled this the ‘exploitation scenario’ and we simulated it by progressively increasing harvesting rate *F* (=[0, 3]), starting from a stable non-fluctuating population (*r*_{0} = 0.75) (figure 1*b*). Intensifying harvesting leads to the abrupt collapse of the population due to the crossing of a fold bifurcation that forces the population to shift from an underexploited to an overexploited state of low biomass.

In both scenarios, we gradually increased the bifurcation parameters, *r*_{0} and *F*, respectively, in 100 equidistant steps. At each step, we burned-in the models for a period of 100 time steps to discard transients, and we simulated another 100 points to use as time series for analysis. For each step of the bifurcation parameter, we produced 1000 replicate time series that were used to estimate nonlinearity and CSD indicators. We also compared the stationary distributions of the indicators to a transient simulation scheme, where growth rate *r*_{0} (=[0.01, 3]) and harvesting rate *F* (=[0, 3]) continuously increased over 200 time steps. In such a transient scheme, we estimated all indicators within a moving window equal to half the size of the time series (i.e. 100 points).

### 2.2. Nonlinearity derived from empirical dynamical modelling

We use nonlinearity (i.e. quantification of state-dependence) as a potential new indicator for anticipating shifts in ecological dynamics. To determine whether a time series reflects changes in state-dependence, we compared the out-of-sample forecast skill of a linear model (i.e. relying on global information for forecasts) versus an equivalent nonlinear model (i.e. relying on local information). This involves state-space reconstruction (also known as EDM) using lagged coordinate embeddings with a two-step procedure as follows. First, we used simplex projection [27] to determine the embedding dimension (*E*) of the system. *E* represents the number of independent variables (axes) needed to reconstruct the system state-space and is operationalized as the number of lagged coordinates used to reconstruct the system attractor. A set of trial values for *E* is used to test different attractor shapes. The best *E* for a given time series is selected according to forecast skill, calculated as the correlation coefficient between actual and model forecast values from the time series. The best *E* also defines the number of points for the neighbourhood around the current system state for forecasting.

Second, using this embedding dimension, S-map [31] is applied to compare linear versus nonlinear forecasting models, by tuning a nonlinear weighting parameter *θ*. *θ* is a parameter that progressively weighs the neighbourhood around the current system state. Forecasts in S-map are made using all points in the attractor but with different weights depending on their proximity to the current system state. At *θ* = 0, the weighting surface is flat and all points in the attractor are equally used for forecasting [32]. At *θ* > 0, the weighting surface gives exponentially greater weight to the local neighbourhood around the current system state. As *θ* increases, points far on the trajectory from the current state receive ever-diminishing weight until they do not contribute any information to the forecast algorithm. If the forecast skill of the nonlinear model (*θ* > 0) outperforms that of the linear model (*θ* = 0), the observed dynamics in the system are classified as nonlinear (or state-dependent). As in simplex projection, forecast skill is evaluated based on the correlation between S-map predicted out-of-sample values and the actual observed values in the time series. A full description of this established methodology with related algorithms can be found in [1,32,48,49].

We applied the above two-step procedure after first-differencing and standardizing both simulated and empirical time series [32]. As EDM requires a time series of at least 30 observations [50], we used simulated time series of 100 time steps and selected empirical records that contained at least 30 points. We produced simplex projections using lagged coordinates of one time step (*τ* = 1) for a series of different embedding dimensions *E*. We used a range of *E*: 1 through 10 for the empirical data and 1 through 3 for the simulated data—we used a smaller range of *E* for the simulated data as we knew the dimensionality of the Ricker model attractor is three according to Whitney's theorem that *n* *≤* *E* *≤* *2n* + 1, where *n* is the dimensionality of the system (*n* = 1 for our model). We chose the best embedding dimension *E* for each time series by estimating the Pearson correlation (*ρ*) between observed and forecast values. The best *E* was then used for fixing the embedding space in the S-map. In the S-map, we estimated forecasting skill as above for a range of *θ* values [0, 0.001, 0.005, 0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 1, 2, 5]. Lastly, we quantified the improvement in forecasting skill of the nonlinear over the linear model as the difference in Pearson correlation: Δ*ρ* = max(*ρ _{θ}* −

*ρ*

_{θ}_{=0}): the maximum difference between the correlation

*ρ*at each

_{θ}*θ*to the correlation

*ρ*

_{θ}_{=0}found for

*θ*= 0. In other words, Δ

*ρ*defines how much the best model outperforms the linear model [40]. If the linear model is the best model, Δ

*ρ*is equal to zero. Δ

*ρ*served as the indicator of nonlinearity.

### 2.3. Critical slowing down indicators

CSD is defined as the decrease in recovery rate upon small perturbations in the vicinity of local bifurcation points [9]. It is a generic property of dynamical systems that undergo transitions between different attractors when a stress parameter crosses a threshold. In mathematical terms, CSD is associated with a diminishing dominant eigenvalue *λ*, where *λ* defines the rate of exponential decay of a perturbation close to equilibrium (Δ*x* *=* e* ^{−λt}*) [22]. The consequence of this slow decay is that both variance and autocorrelation of the recorded ecosystem dynamics will increase close to a transition point [12]. We estimated variance as coefficient of variation (

*CV*= standard deviation/mean), and autocorrelation at lag-1 (

*AR1*) as the Pearson correlation between lagged time series at one time step [51].

### 2.4. An empirical example: comparing stability across fish populations

We tested the indicators in empirical fisheries datasets collected from scientific surveys on the Northeast Shelf (NES) in the northwest Atlantic and in the southern California Current Ecosystem (CCE) in the eastern Pacific. NES data were collected through the Northeast Fisheries Science Center. These data are relative biomass estimates generated from an annual fall bottom trawl survey and include 29 stocks of demersal fishes sampled from 1963 to 2008. Of these, 20 stocks were exploited (subject to fishing pressure). CCE data were collected through the California Cooperative Oceanic Fisheries Investigations (CalCOFI). These data are relative biomass estimates generated from regular ichthyoplankton tows and include 29 coastal-neritic fish species [52] sampled from 1951 to 2007. Among the 29 species, 16 were exploited and 13 were unexploited. We omitted missing values in estimating the indicators. The NES data are available from the Northeast Fisheries Science Center, National Marine Fisheries Service (NMFS), National Oceanic and Atmospheric Administration (NOAA), USA (http://www.nefsc.noaa.gov/nefsc/saw/). The CCE fish data are available from the Southwest Fisheries Science Center, NMFS, NOAA, USA (http://coastwatch.pfeg.noaa.gov/erddap/search/index.html?page=1&itemsPerPage=1000&searchFor=calcofi).

We compared indicators from the fisheries data to distributions of indicators of exploited and unexploited populations from our simulated scenarios. In the demographic change scenario, we marked populations with *r* > 1.48 as exploited, and populations with *r* < 1.48 as unexploited, while in the exploitation scenario, exploited populations were harvested with rate *F* between 0.57 and 1.18, and unexploited populations with rate *F* below 0.57. We then randomly selected 2000 exploited and unexploited population time series.

All simulated data and the estimation of the nonlinearity and CSD indicators were produced using Matlab v. 2015a (Mathworks). For open source options, nonlinearity indices can be estimated using the rEDM package (https://cran.r-project.org/web/packages/rEDM/index.html), and CSD indicators can be computed using the R package earlywarnings (https://cran.r-project.org/web/packages/earlywarnings/index.html).

## 3. Results

In our demographic change scenario, increases in growth rate led to a sequence of well-understood transitions: a shift from stable point to excitable point equilibrium dynamics, followed by a transition to regular oscillations, and finally a shift to chaos. In addition, there was also a transition to extinction when growth rate decreased to 0. Obviously, stochasticity makes the boundaries between these attractors difficult to identify, as exemplified by the simulated time series in figure 1*a*. We found that CSD indicators (*CV*, *AR1*) and nonlinearity (Δ*ρ*) could partly differentiate the transitions between these dynamical regimes (figure 1*c*). Around the boundary to spiralling dynamics, *CV* and Δ*ρ* were the lowest, while *AR1* changed sign from positive to negative values. These patterns reflected the fact that close to this boundary disturbances die off fast but in an oscillating manner, determined by the decreasing local eigenvalue that shifts sign in the deterministic model (see the electronic supplementary material, figure S1). Before the onset of cycles, all indicators increased, reflecting the progressive loss of stability in the deterministic model. In particular, *AR1* reached its highest (negative) value, signalling the birth of the deterministic cyclic attractor. As the instability increased further towards the transition to chaos, *CV* and Δ*ρ* continued to strongly increase, but *AR1* decreased weakly implying that correlated patterns were disrupted due to the chaotic dynamics. Lastly, in the opposite direction, the approach to extinction was marked by a clear increase in all indicators that reflected a considerable loss of stability*.*

In the exploitation scenario, a gradual increase in harvesting rate *F* caused populations to experience a slight decrease until collapsing to an overexploited state (figure 1*b*). Owing to stochasticity, the collapse occurred earlier than the actual transition of the deterministic model (figure 1*b*). Approaching the transition, the stability of the point equilibrium gradually diminished (electronic supplementary material, figure S1*d*). Consequently, we found that both CSD indicators (*CV*, *AR1*) increased (figure 1*d*). Strikingly, nonlinearity (Δ*ρ*) increased as well before the collapse (figure 1*d*). Once the transition was crossed, the pattern reversed: *CV*, *AR1* and Δ*ρ* all dropped, marking the progressive gain in stability and the decrease in nonlinearity of the overexploited state.

Nonlinearity Δ*ρ* is estimated after the empirical reconstruction of the underlying dynamical attractor (best embedding dimension *E*) and the fitting of a model with varying levels of state-dependence (best weighting factor *θ*). In figure 2, we explicitly monitored both these quantities in our simulations. In the demographic change scenario, we found that the fraction of nonlinear models with best embedding dimension *E* = 2 peaked around the onset of the chaotic regime, while the proportion of linear models with the best forecasts (*θ* = 0) dropped to 0 almost at the shift to cyclic dynamics (figure 2*a*,*c*). In the exploitation scenario (figure 2*b*,*d*), the patterns were noisy, although we noticed an increasing trend and a peak in the most nonlinear model (*θ* = 5) at the shift to overexploitation.

We also compared results from the stationary distributions to a situation where external conditions change continuously with time. We simulated this more realistic scenario by linearly increasing growth and harvesting rates over 200 time steps, and we measured CSD and nonlinearity indicators within a moving window along the time series (figure 3). We found similar patterns to those obtained in the case of stationary distributions (figure 1). However, we found strong uncertainties in the patterns of nonlinearity especially for the demographic change scenario (figure 3*g*,*h*).

Ideally, for populations that have a history of monitoring, there exist biomass estimates at different levels of stress (environmental conditions, exploitation) that help compare stability among populations using the above presented indicators. Unfortunately, in our empirical fisheries records, we have no explicit values of demographic rates (such as *r* in our model), nor of harvesting rate *F* as in the simulations. Instead, we can only discriminate populations based on whether they were commercially fished (i.e. exploited) or not (i.e. unexploited or bycatch). Using this discrimination as a proxy for exploitation stress, we tested for differences in *CV*, *AR1* and Δ*ρ* between exploited versus unexploited populations and we compared them to the patterns found from our scenarios (figure 4*a*–*c*). We found that mean Δ*ρ* and *AR1* were higher for exploited populations than unexploited populations in both datasets (figure 4*e*,*f*), whereas mean *CV* was higher in exploited populations in the CCE but not in the NES dataset (figure 4*d*). However, only Δ*ρ* and *CV* were significantly higher in the CCE exploited populations (Mann–Whitney *U*-test). Taking both datasets together (electronic supplementary material, table S1), we found negative correlations between *CV* and *AR1* for both exploited and unexploited populations, stronger positive correlation between *CV* and Δ*ρ* for unexploited than exploited populations, and also stronger positive correlation between *AR1* and Δ*ρ* for unexploited than exploited populations. However, similar to our simulated time series, only the correlations between *CV* and Δ*ρ* were significant.

## 4. Discussion

Driven by chaotic dynamics, stochasticity or alternative attractors, our ability to accurately describe nonlinear systems and to efficiently manage them remains limited. Here, we explored whether EDM can complement CSD indicators for detecting transitions in the dynamics of ecological systems. Both approaches belong to a class of model-free methods for describing ecological dynamics. EDM has never been explored as a means of measuring ecosystem stability or detecting the risk of dynamical transitions in the generic sense as CSD indicators do. Our results show that elevated nonlinearity Δ*ρ* based on EDM can be potentially used as an indicator for the proximity to transitions between alternative states as well as to transitions from stable dynamics to irregular oscillations.

Nonlinearity (or state-dependence), measured as Δ*ρ*, basically captures changes in the deterministic stability of the Ricker model in terms of degree of state-dependence (figure 1 and electronic supplementary material, figure S1). In the demographic change scenario, the gradual increase of Δ*ρ* from a stable equilibrium to oscillations and finally chaos is expected as populations become increasingly unstable and unpredictable. This implies that future states are better predicted by using local information: that is, points that are closer together in the reconstructed state space of dimension 2 (figure 2*a*). On the other hand, the rise of Δ*ρ* before the shift to the alternative state is less intuitive in the exploitation scenario. Owing to critical slowing, there is a build-up of memory in the produced time series. This implies that the best reconstructed attractor is of dimension 1 (figure 2*b*), which basically means that future states are better predicted by using local information that corresponds to points closer in time.

One may observe a large difference in the behaviour of the indicators between the two scenarios. In the exploitation scenario, all three indicators peak at the transition to the alternative state (figure 1*d*), while in the demographic change scenario, indicators change in a smooth way (figure 1*c*). This observation raises challenges in the operationalization of these indicators. Ideally, an indicator should produce a distinct mark of changing dynamics. This is true for the exploitation scenario, because the change in the dynamics is related to a true discontinuity in the stability of the system (electronic supplementary material, figure S1*d*). In this case, all three indicators serve as independent signals of the approaching discontinuity. In the demographic change scenario, however, changes in dynamics are continuous (electronic supplementary material, figure S1*c*). Therefore, it is expected that these indicators (as any other indicator) would also change in a continuous way. For this reason, interpretations of changes in dynamics could only be based on parallel monitoring and comparison among multiple indicators [51].

For example, we find that elevated nonlinearity and rising variability are strongly correlated in both scenarios (Pearson correlation: demographic change scenario 0.86, exploitation scenario 0.88). This observation implies that similar changes in nonlinearity (Δ*ρ*) and variability (*CV*) may be a good proxy for detecting shifts in population dynamics that are characterized by non-fixed point attractors and strong stochasticity. This is because rising variance and nonlinearity can be understood as fingerprints of increasing state-dependence [31] that is not affected by the nature of the underlying deterministic attractor. Population dynamics are typically the result of a mix of transients across stable and unstable equilibria affected by environmental and demographic stochasticity. CSD indicators can capture changes in these dynamics but only when it comes to stable equilibria in the presence of weak stochasticity [12]. Identifying transitions across chaotic attractors or, more generally, in systems with nonlinear dynamics may be difficult with CSD indicators [53,54]. By contrast, EDM-derived nonlinearity may be broader in its application as it can capture changes in dynamics beyond stable attractors typical of the dynamics encountered in natural systems.

Our findings suggest that an ecosystem manager could estimate both CSD and EDM metrics in order to rank populations according to their stability [55]. For instance, Krkosek & Drake [56] looked at patterns of *CV* and *AR1* for Pacific salmon populations and found that they were higher for pink salmon stocks that had a population growth rate close to zero. In this work, the authors assumed that salmon populations would suffer a transcritical bifurcation due to growth rates approaching zero. Trends in *CV* and *AR1* in our analyses are in line with these findings (figure 1). Moreover, we also find that nonlinearity increases at decreasing growth rates approaching to zero. Our results imply that regardless of the type of transition (be it a transcritical, a fold or a flip bifurcation), a simultaneous increase in *CV* and nonlinearity could signal proximity to a transition in ecological dynamics. When tested in the empirical data, *CV* and nonlinearity also show a consistent positive relationship (electronic supplementary material, table S1) that follows our theoretical expectations. However, we failed to find an unequivocal signature in the indicators that discriminate the exploited from the unexploited group. For example, we found elevated nonlinearity and *AR1* in the exploited populations for both NES and CCE datasets, but higher variability only for the CCE dataset (figure 4). Although the datasets we analysed are sampled under the same environmental conditions, differences in life-history traits between species can modify population responses to the same stress [52]. If information on such differences is available, correcting for them across populations may improve the interpretation of the indicators.

The indicator patterns we identified in our Ricker-type model will most probably also hold for other discrete or continuous population models with similar bifurcations. Nonetheless, a number of factors can confound the detection capacity of the indicators. For example, we found strong fluctuations in nonlinearity estimates, especially in the exploitation scenario (electronic supplementary material, figure S2f). Changes in variance and autocorrelation can also be unreliable in the presence of short time series [51], high levels of stochasticity [57], fast changing stress drivers [54,58] or due to portfolio effects [56] and life-history strategies [59]. Further research is needed to find if similar constraints hold for the elevated nonlinearity indicator we propose here.

In the quest for understanding and anticipating ecosystem responses to stress, testing novel and alternative approaches is of high priority. The indicators we examined here contribute to equation-free, data-driven approaches that aim at quantifying differences in the stability of populations under increasing environmental stress. Translating such differences to a risk assessment scheme might be a useful tool for improving ecosystem management in the face of global environmental change.

## Authors' contributions

V.D., S.G., C.H.H. and G.S. conceived the study. V.D. and C.H.H. performed analyses. S.G. and C.H.H. contributed data. V.D. wrote the paper with contributions from all co-authors.

## Competing interests

We decalre we have no competing interests.

## Funding

V.D. acknowledges support by a Marie Curie EU IntraEuropean fellowship (2013–2015) and funding from the Center for Adaptation to a Changing Environment (ACE) at ETH Zürich. C.H.H. is supported by the National Center for Theoretical Sciences, Foundation for the Advancement of Outstanding Scholarship, and Ministry of Science and Technology of Taiwan.

## Acknowledgements

We would like to thank Florian Grziwotz and Arndt Telschow for valuable comments on the manuscript, as well as the four anonymous reviewers and the editor who all helped improve the quality of this work. We also thank the California Cooperative Oceanic Fisheries Investigations, Southwest Fisheries Science Center and Northeast Fisheries Science Center for the use of the data. Hui Liu assembled the NES data used in the study.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3691876.

- Received October 21, 2016.
- Accepted February 3, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.