Detecting alternative attractors in ecosystem dynamics
Detecting alternative attractors in ecosystem dynamicsWe use empirical dynamical modeling, a set of equation-free tools for analyzing non-linear time series (for a review and assumptions see25,26, respectively), to test if the temporal dynamics of alternative dynamical regimes are qualitatively different. Empirical dynamic modeling builds fundamentally on Takens embedding theorem, which shows that attractors of multi-dimensional dynamical systems can be reconstructed using higher order lags of its embedded time series27. However, if a dynamical system has gone through a bifurcation, or switched to an alternative basin of attraction, attractors are qualitative dissimilar in the two regimes. Theoretically, this infers that it should be possible to reconstruct the attractor of one regime using information from the same regime, but not from the other regime. In practice, this implies that if a model (attractor reconstruction) based on one dynamical regime is used to predict the dynamics of variables from the same dynamical regime predictions should be accurate (i.e. low prediction errors), whereas if an attractor reconstruction based on one dynamical regime is used to predict the dynamics of variables of another attractor predictions should be less accurate (i.e. high prediction errors). We make use of this idea by specifically testing if prediction errors of across and within regime predictions are different. As explained below this idea can be used for both univariate and multivariate time series data.Univariate approachUnivariate attractor reconstructions can be found using the simplex algorithm28,29. First, for a given dynamical regime, a time series can be split into a library of vectors, and each vector is described by$${underline{y}}_{A}(t)= < {Y}_{A}(t),{Y}_{A}(t-1),{Y}_{A}(t-2),ldots ,{Y}_{A}(t-(E-1)) > ,$$
(1)
where ({Y}_{A}(t)) is an observation of variable Y at time t in dynamical regime A and E is the reconstructed attractors embedding dimension. Using the simplex projection algorithm, a one-step ahead forecast is produced as follows:$${hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{B}=mathop{sum}limits_{m=1ldots E+1}{w}_{m}{Y}_{B}({t}_{m}+1),$$
(2)
where tm is a time index of an observation in dynamical regime B, E is the embedding dimension of regime B, and wm is an exponential weighting described by:$${w}_{m}={u}_{m}/mathop{sum}limits_{n=1,ldots ,E+1}{u}_{n},$$
(3)
where n and m belongs to the set of the E+1 nearest neighbors of vector ({underline{y}}_{A}(t)) in the set of vectors ({{underline{y}}_{B}({t}_{m})}), ({u}_{m}=exp {-d[{underline{y}}_{A}(t),{underline{y}}_{B}({t}_{m})]/d[{underline{y}}_{A}(t),{underline{y}}_{B}({t}_{1})]}), and (d[{underline{y}}_{A}(t),{underline{y}}_{B}({t}_{1})],)is the Euclidean distance between the prediction vector ({underline{y}}_{A}(t)) and its nearest neighbor ({underline{y}}_{B}({t}_{1})) in the set ({{underline{y}}_{B}({t}_{m})}).The only parameter that is estimated using the simplex algorithm is the embedding dimension E. This parameter is estimated by optimizing the correlation between observations (({Y}_{A}(t+1))) and predictions (({hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{A})) using a leave-one-out cross validation approach (See Supplementary Discussion). The embedding dimension E and its corresponding set of E-dimensional vectors (Eq. 1) constitutes the reconstructed attractor, MA, of a given dynamical regime A. This reconstructed attractor (MA) is then used to predict data for both the same dynamical regime (({hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{A})), and the contrasting dynamical regime ({hat{Y}}_{B}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{A}). Likewise, the reconstructed attractor MB can be used to predict time series dynamics from both dynamical regimes; that is, ({hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{B}) and ({hat{Y}}_{B}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{B}), respectively.Multivariate approachA multivariate time series describes a number of simultaneously evolving variables. For example, a bivariate time series can be described by variables X and Y. For such time series, Sugihara et al.30 developed an approach for testing if two variables (time series) are dynamically coupled. Their methodology builds on the fact that a reconstructed attractor should map 1:1 to the original attractor on which the reconstruction is based. This infers that two attractor reconstructions (based on two different variables) should also map 1:1 to each other30. Practically, this means that if two variables are dynamically coupled one-time series should be predictable based on an attractor reconstruction of another variable. However, if a dynamical system has gone through a bifurcation, or potentially switched to an alternative basin of attraction, a new set of rules will govern the dynamics of the system. Hence, a new attractor should have emerged. Now, since this new attractor is most likely governed by a new set of rules it should be difficult to predict the dynamics of this new alternative attractor based on information from the former attractor. Thus, if one variable in one dynamical regime is used to predict another variable in another dynamical regime, predictions should be biased. Yet, if one variable from one dynamical regime is used to predict another variable from the same regime predictions should be more accurate.The simplex algorithm can be used to make predictions of a variable Y using a time series of another variable X30. Predictions are produced as follows:$${hat{Y}}_{{{{{{boldsymbol{A}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{B}=mathop{sum}limits_{m=1ldots E+1}{w}_{m}{Y}_{B}({t}_{m}),$$
(4)
where tm is the time series index of a vector of variable X of dynamical regime B, wm is an exponential weighting based on variable X:$${w}_{m}={u}_{m}/mathop{sum}limits_{n=1,ldots ,E+1}{u}_{n},$$
(5)
where n and m belongs to the set of the E+1 nearest neighbors of ({underline{x}}_{A}(t)) in ({{underline{x}}_{B}({t}_{m})}), ({u}_{m}=exp {-d[{underline{x}}_{A}(t),{underline{x}}_{B}({t}_{m})]/d[{underline{x}}_{A}(t),{underline{x}}_{B}({t}_{1})]}), and (d[{underline{x}}_{A}(t),{underline{x}}_{B}({t}_{1})],)is the Euclidean distance between the prediction vector(,{underline{x}}_{A}(t)) and its nearest neighbor ({underline{x}}_{B}({t}_{1})) in dynamical regime (B).The reconstructed attractors, MA and MB, for each variable and regime are found using the univariate simplex algorithm described above28,29,30. Similar to the univariate case, the reconstructed attractor (MA) is used to predict data from the same dynamical regime (({hat{Y}}_{{{{{{boldsymbol{A}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{A})), and to predict time series of a contrasting dynamical regime (({hat{Y}}_{{{{{{boldsymbol{A}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{B})). Yet, it is important to stress that MA here reflects an attractor reconstruction based on a variable that is not being predicted (that is, variable X is used to predict variable Y). This prediction approach thus infers that predictions are made on data that was not used to fit the model (X predicts Y and vice versa). Thus, neither across nor within regime predictions are made on data used to fit a model.Test statisticWe used mean absolute prediction errors to test for difference between across and within regime predictions. Alternative metrics, such as mean sum of square errors, can also be used. However, since our approach gives skewed prediction errors we used mean absolute prediction errors to reduce the impact of extreme values. Further, since the absolute prediction errors are non-normally distributed we used a permutation test. The null hypothesis that is tested reads:$$H0:{{{{{rm{MAP{E}}}}}}}_{A} < {{{{{rm{MAP{E}}}}}}}_{w},$$
(6)
where MAPEA is the mean absolute prediction error for across regime predictions (that is, ({{{{{rm{MAP{E}}}}}}}_{A}=frac{1}{n}mathop{sum}limits_{t=1:n}{{{{{rm{abs}}}}}}({hat{Y}}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{B}-{Y}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t))), and ({{{{{rm{MAP{E}}}}}}}_{w}) is the mean absolute prediction error for within regime predictions (that is, ({{{{{rm{MAP{E}}}}}}}_{w}=frac{1}{n}mathop{sum}limits_{t=1:n}{{{{{rm{abs}}}}}}({hat{Y}}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{A}-{Y}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t))). A test is consider significant if observed difference in across and within regime mean prediction errors is larger than the 95th percentile of 1000 permuted data sets.Food-chain modelWe used a food-chain model parameterized as in McCann and Yodzis31 to simulate food-chain dynamics:$$frac{{{{{{rm{d}}}}}}R}{{{{{{rm{d}}}}}}t}=Rleft(1-frac{R}{K}right)-frac{{x}_{c}{y}_{c}CR}{R+{R}_{0}}$$
(7)
$$frac{{{{{{rm{d}}}}}}C}{{{{{{rm{d}}}}}}t}={x}_{c}Cleft(-1+frac{{y}_{C}R}{R+{R}_{0}}right)-frac{{x}_{P}{y}_{P}PC}{C+{C}_{0}}$$$$,frac{{{{{{rm{d}}}}}}P}{{{{{{rm{d}}}}}}t}={x}_{P}Pleft(-1+frac{{y}_{P}C}{C+{C}_{0}}right),$$where R is the resource density, C consumer density, and P predator density. All parameters, except half-saturation constants R0 (here set to 0.16129) and C0 (here set to 0.5), and resource carrying capacity K, are derived from bioenergetics and body size allometry30 (xc = 0.4, yc = 2.009, yp = 2.876, R0, r = 1, xp = 0.08).This model can display a rich set of dynamics depending on parameter values31. Here we alter resource carrying capacity K in order to simulate the dynamics (using the deSolve package32 in R) of qualitatively different attractors (See Supplementary Fig. 1; K = 0.78, equilibrium; K = 0.85; two-point limit cycle; K = 0.92, four-point limit cycle; K = 0.997, chaotic dynamics). Every fifth time step of the simulated dynamics, corresponding to a sampling frequency of ≈10 samples per cycle for the 2-point limit cycle, was sampled. Observation noise was thereafter added to the deterministic dynamics produced by the model:$${N}_{l}(t)={N}_{l}^{prime}(t)+rho * e(t);e(t) sim N(0,{sigma }_{N^{prime_{l}}}),$$
(8)
where (N_{l}^{prime}(t)) is the abundance of species l (P, C or R) simulated by the food-chain model at time point t, (rho) is the level of observation noise and ({sigma }_{N_{l}^{prime}}) is the standard deviation of the deterministic dynamics of species l produced by the food chain model.In order to investigate how time series length and observation noise affects the probability of detecting alternative attractors we derived probability landscapes. These were derived by testing the null-hypothesis (H0:(|{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{i}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}| > |{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{j}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|); See Test statistic above) across 100 replicates for each combination of time series length and level of observation noise, (rho). Time-series length was varied from 10 to 100 in steps of 10, and observation noise, (rho), was varied from 0.01 to 0.3 in steps of 0.01, in total yielding 300 combinations of observation noise and time series length, for each combination of dynamical regimes i and j. Predator dynamics was used to predict consumer and resource dynamics using the multivariate approach described above (results for the cases where consumer or resource dynamics are used to predict the other species´ dynamics are presented in Supplementary Figs. 2, 3). All time series were standardized ((mu =0;sd=1)) prior testing for dynamical difference.Experimental data setThe experimental data set was given by Fussman et al.7. This data set contains 14 time series of a predator Brachionus calyciflorus and its prey Chlorella vulgaris derived from chemostat experiments. Time series for different dilution rates were produced by keeping the dilution rate fixed in different chemostats (Supplementary Figs. 3–11). Brachionus calyciflorus and Chlorella vulgaris time series were used to predict Chlorella vulgaris and Brachionus calyciflorus time series, respectively, using the multivariate approach described above. We tested for qualitative difference in the temporal dynamics across all time series, which were standardized ((mu =0;sd=1)) prior testing.Alternative stable state modelWe used a stochastic version of a well-known alternative stable state model4,33 to produce alternative stochastic dynamical regimes. The model is described by:$${{{{{rm{d}}}}}}x=left(xleft(1-frac{x}{{{{{{rm{K}}}}}}}right)+frac{c{x}^{2}}{1-{x}^{2}}right){{{{{rm{d}}}}}}t+sigma {{{{{rm{d}}}}}}w,$$
(9)
where K is the carrying capacity (here set to 11), c is a harvest rate, and σ (here set to 0.01) is the magnitude of noise which is described by a Wiener process (dw).The model was simulated for fixed harvest rates (c) assuming that the system state resides in either of its two basins of attraction. The initial value for the simulation was set to the equilibrium of the noise-free model skeleton for fixed harvest rates c, and σ is set low in order to avoid stochastic flips, so-called flickering, between alternative basins of attraction. Dynamics was integrated (Δt = 0.01) using the matlab-package SDE-Tools34.In order to investigate how time-series length and harvest rate, c, affects the probability of detecting alternative attractors in stochastic regimes we derived probability landscapes.These were derived by testing the null-hypothesis H0:(|{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{i}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}| > |{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{j}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|) (permutation test p = 0.05) across 100 simulated data sets for each combination of time series length and harvest rate, c. Time-series length was varied between 50 and 150 in steps of 10, and c was varied between 1.83 and 2.73 in steps of 0.05, in total yielding 209 combinations of time series length and harvest rate. Each time series was standardized ((mu =0;sd=1)) prior testing for difference in temporal dynamics of contrasting regimes.Natural time-series dataIn a previous study on early warning signals of impending regime shifts, Gsell et al.18 used breakpoint analysis to identify two potential alternative dynamical regimes. We here test if these two-time series segments constitute alternative dynamical attractors. Prior analysis, we imputed a few missing observations (n = 24) using a kalman smoother35. The two time series segments, i.e. pre- and post-breakpoint time series, were standardized ((mu =0;sd=1)) prior testing for dynamical difference.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
