Nonlinear time series analysis unravels underlying mechanisms of interspecific synchrony among foliage‐feeding forest Lepidoptera species

Funding information Operational Programme Research, Development and Education, Grant/Award Number: CZ.02.1.01/0.0/0.0/16_019/ 0000803; JSPS KAKENHI, Grant/Award Numbers: 15H04613, 18K14797, 16K18625, 16H04846 Abstract Interspecific synchrony, that is, synchrony in population dynamics among sympatric populations of different species can arise via several possible mechanisms, including common environmental effects, direct interactions between species, and shared trophic interactions, so that distinguishing the relative importance of these causes can be challenging. In this study, to overcome this difficulty, we combine traditional correlation analysis with a novel framework of nonlinear time series analysis, empirical dynamic modeling (EDM). The EDM is an analytical framework to identify causal relationships and measure changing interaction strength from time series. We apply this approach to time series of sympatric foliage-feeding forest Lepidoptera species in the Slovak Republic and yearly mean temperature, precipitation and North Atlantic Oscillation Index. These Lepidoptera species include both free-feeding and leaf-roller larval life histories: the former are hypothesized to be more strongly affected by similar exogenous environments, while the latter are isolated from such pressures. Correlation analysis showed that interspecific synchrony is generally strongest between species within same feeding guild. In addition, the convergent cross mapping analysis detected causal effects of meteorological factors on most of the free-feeding species while such effects were not observed in the leaf-rolling species. However, there were fewer causal relationships among species. The multivariate S-map analysis showed that meteorological factors tend to affect similar free-feeding species that are synchronous with each other. These results indicate that shared meteorological factors are key drivers of interspecific synchrony among members of the free-feeding guild, but do not play the same role in synchronizing species within the leaf-roller guild.

Here we explore the causes of interspecific synchrony using a system studied by  and Raimondo, Turcáni, et al. (2004) consisting of sympatric foliage-feeding forest Lepidoptera species in the Slovak Republic. This insect community can be categorized into two groups based on larval feeding biology, that is, freefeeding species and leaf-rollers Raimondo, Turcáni, et al., 2004). These two species guilds are likely to be differently influenced by environmental conditions and/or predators that use search images within same feeding group. Applying cross-correlation analysis to lepidopteran abundance time series and using a model analysis of prey dynamics with shared predators, they found that (a) significant synchrony occurs more frequently among species exhibiting similar feeding strategies, and (b) pairs of prey species yielding similar search images to predators are easily synchronized compared to species with different images (Raimondo, Turcáni, et al., 2004). Based on these results they concluded that generalist predators are likely causal drivers of observed interspecific synchrony among forest Lepidoptera (Raimondo, Turcáni, et al., 2004). However, their conclusion dependent upon inference from the correlation, which is not necessarily a good indicator of causality (Sugihara et al., 2012), and inference from the behavior of a population model, whose validity in the real world is unknown. Thus, we still cannot rule out the possibility of the other factors, such as shared environmental factors and direct interactions between lepidopteran species.
A recently developed approach to nonlinear time series analysis, empirical dynamic modeling (EDM), is a promising method for overcoming the inference problems described above. The approach of EDM is based on theories of nonlinear dynamical systems (e.g., Deyle & Sugihara, 2011;Takens, 1981), and enables testing for causal relationships between observed time series. In particular, we focused on the two methods of EDM, convergent cross mapping (CCM) and the multivariate S-map procedure. The CCM method is used to determine the direction of causal effect (Sugihara et al., 2012;Ye et al., 2015), and the S-map analysis is used to measure the strength of the causal interaction (Deyle, May, Munch, & Sugihara, 2016). The EDM is applicable to a wide variety of problems involving inference of causal relationships from ecological time series, and its application in ecology is rapidly increasing (e.g., Kawatsu & Kishi, 2018;Sugihara et al., 2012;Ushio et al., 2018). For example, Kawatsu and Kishi (2018) leveraged the CCM method and the multivariate S-map analysis to identify critical interactions in the competition dynamics between two bean beetle species, where two different interspecific interactions, that is, resource competition between larvae and reproductive interference between adults, can co-occur. Similarly, application of EDM approach to lepidopteran time series and candidate causal series (e.g., weather conditions and predators) could be expected to provide more definitive identification of the drivers of interspecific synchrony. For example, since freefeeding larvae are likely exposed to similar environmental conditions while leaf-rolling larvae are likely more isolated from exogenous environments (Raimondo, Turcáni, et al., 2004), we hypothesize that the free-feeding species may be more strongly synchronized by shared meteorological conditions (i.e., the direction of causality and the sign of climate effects) than would leaf-rollers.
In this study, we reanalyze the dynamics of the foliagefeeding forest Lepidoptera community studied by Raimondo, Turcáni, et al. (2004) applying EDM analysis. The analysis consists of the following four steps. First, to confirm the degree of interspecific synchrony, we quantify crosscorrelation between each species. Second, to test the importance of the weather conditions on interspecific synchrony, we performed CCM analysis between the lepidopteran time series and meteorological time series. Third, to explore the possibility of direct interactions, we performed the CCM analysis between the lepidopteran time series. Finally, to infer the sign of the detected causal effects of climatic factors on lepidopteran species, we applied multivariate S-map analysis. The results shed light on the underlying mechanisms causing synchrony among populations of the foliage-feeding forest lepidopteran community.

| Data
The lepidopteran data was collected by the Slovak Forest Research Institute, Zvolen, the Slovak Republic. Insect counts were made from 1955 to 1964 and from 1966 to 1982 in early May (usually between May 10 and 15) at 20 sites throughout Slovakia. The size of forest stands used for monitoring varied from several hectares to several hundred hectares. Lepidopteran insects were collected from each site by selecting lower canopy branches from 20 oak, Quercus, trees and beating them onto a net to dislodge insects. For details of sampling, refer to Raimondo, Turcáni, et al. (2004). As in Raimondo, Turcáni, et al. (2004), we analyzed data for the 10 most abundant species, representing four families, for the analysis. Of these species, six were primarily free feeders as larvae and four were leaf-rolling larvae (all species names are cited in Table S1). The lepidopteran time series were ln(x + 0.5) transformed (Yamamura, 1999) for each species for all analyses. All data were missing from 1965 so values were filled using a spline curve ("smooth spline" function in R; Figure 1).
We also used three meteorological time series, monthly average precipitation, mean monthly minimum temperature recorded in Bratislava, the Slovak Republic (http://www. ecad.eu/download/ensembles/download.php), and the North Atlantic Oscillation Index (NAOI, http://www.cpc/noaa.gov/ data/teledoc/nao.shtml) from 1955 to 1982 for the causality analysis. The NAOI is a climatic index which measures the difference of atmospheric pressure at sea level in the North Atlantic Ocean and has a well-known correspondence to broad-scale weather conditions in regions near the Atlantic Ocean (Hurrell, 1995).

| Spatial covariance analysis
We hypothesized that interspecific synchrony among forest Lepidoptera species would be different between the freefeeding and the leaf-rolling species. Specifically, we anticipated that interspecific synchrony is stronger between species pairs within the same feeding group than those between different groups. To test this prediction, we evaluated the correlation among all possible pairs of time series of the 10 species. Interspecific synchrony was quantified by calculating the spatial cross-correlation function (SCCF, Bjørnstad, Ims, & Lambin, 1999) for all possible species pairs; the average correlation and its confidence interval were obtained from 1,000 bootstrap samples. The significance of interspecific F I G U R E 1 Log-transformed abundance dynamics of 10 lepidopteran species. Different colors indicated the 20 sites throughout the Slovak Republic from 1955 to 1982. The panels (a)-(j) correspond to the dynamics of Agriopis aurantiaria, Alsophila aescularia, Erannis defoliaria, Operophtera brumata, Orthosia cruda, Lymantria dispar, Aleimma loefilingia, Archips xylosteana, Eudemis profundana and Totrix virida, respectively. The five species of the upper row and L. dispar are categorized into the free-feeder, and the other species are categorized into the leafroller group [Color figure can be viewed at wileyonlinelibrary.com] synchrony was judged based on the calculated SCCF values with a Bonferroni correction. The correlation analysis was performed by using "ncf" package for R version 3.5.1 (R Core Team, 2018).

| EDM analysis
The EDM is based on Takens' theorem (Takens, 1981) and also on the related nonlinear dynamical theories (e.g., Deyle & Sugihara, 2011). A multivariate system is geometrically represented as an attractor (or a manifold), a set of states and their trajectories in a multidimensional space where each coordinate represents to each system variable. Takens' theorem proves the original attractor (e.g., the community dynamics of interacting Lepidoptera species) can be reconstructed with time-delayed coordinates of an observed time series even if all system variables cannot be observed (Takens, 1981). The reconstruction of the original attractor with a time-delay embedding is called state space reconstruction (SSR), and the CCM leverages the SSR theories to test causal interaction. We here explain the algorithm by using a simple example. If a meteorological factor X affects the dynamics of a Lepidoptera species Y then the time series of Y contains dynamical information of the causality from X. Therefore, the attractor reconstructed from the Y's time series alone preserves similar mathematical properties to that reconstructed from the X's time series, such as neighbor points and their trajectories. These dynamical similarities improve with an increase in the number of points on the SSR attractor of Y so that a causality test can be accomplished by checking the prediction skill which is evaluated by similarity with the observed X and predicted X with Y's attractor (cross-mapping; Sugihara et al., 2012). In contrast, due to the lack of causal effects of Lepidoptera Y on meteorological factor X, there are no information of Y in X's time series and the CCM analysis would fail in prediction.
Although EDM is ordinarily applied to a single long time series (>30 time points), a lepidopteran time series from a site contains at most 27 time points so that it might be insufficient to extract dynamical information even with EDM methods. To avoid this problem, we combined time-series data from different sites to yield single longer time series for each lepidopteran species, prior to the analysis. This procedure is reasonable because time series of a same species is expected to contain similar dynamical information (Clark et al., 2015;Hsieh, Anderson, & Sugihara, 2008;Kawatsu & Kishi, 2018). Then, these combined time-series data were standardized with those means and standard errors. Note that an SSR was performed in order that a reconstructed attractor did not cross gaps between different time series. All EDM analyses were tested with leave-one-out (i.e., jackknife) cross-validation.
We conducted the CCM analysis between the Lepidoptera time series and the meteorological time series, then applied the CCM analysis for all possible pairs of the Lepidoptera species. First, we determined the optimal embedding dimension E i for each Lepidoptera species, based on the prediction skill which is estimated by simplex projection (Sugihara & May, 1990) of each time series. We searched the dimension up to E i = 15 to ensure a sufficient data sizes of the reconstructed attractor and calculated the cross-map skill ρ (Pearson's correlation coefficient) with the minimum library length and the maximum total data size (to equalize the number of points in the reconstructed attractor of each Lepidoptera species, we used E i + 1 as the minimum length). In addition, we generated null time-series data of a putative causal factor by using twin-surrogate method (Thiel, Romano, Kurths, Rolfs, & Kliegl, 2006), which enables the generation of time-series data that preserve the nonlinearity of the focused time series but destroy the dynamic coupling between variables. We judged the causality as significant if the cross-map skill at the maximum data is higher than both the minimum and the surrogate using Fisher's Δρ test (Chang, Ushio, & Hsieh, 2017).
To measure the effect of causal drivers on lepidopteran time series, we applied multivariate S-maps to the Lepidoptera-weather dataset. The S-map (sequential locally weighted global linear map) procedure sequentially calculates the best linear model by leveraging information from all points in the reconstructed manifold (Sugihara, 1994). Specifically, an S-map predicts the future system state at a target time point by giving greater weight to neighboring points on the attractor with a nonlinearity parameter θ (Sugihara, 1994). Applying the S-map procedure to a multivariate embedding with a causal variable led to a good approximation of the effects of causal interaction (Deyle et al., 2016). We created multivariate embedding of the lepidopteran time series as follows. First, based on optimal E i determined above, a univariate embedding was reconstructed with one lepidopteran time series. Then, the last coordinate of the univariate embedding was replaced by one of the causal weather time series. The effects of the causality were sequentially obtained as the E i length coefficient vector. We performed this procedure for the subset of pairs where the causality of weather conditions was detected by CCM analysis. Furthermore, we performed the multivariate S-map analysis to the relationship between lepidopteran species with causal interaction detected by the above CCM analysis. EDM analyses described above were performed using the R version 3.5.1 (R Core Team, 2018) and the EDM package (rEDM version 0.7.2; doi: 10.5281/ zenodo.1081784) for R. We attached the R code and the time series data set that reproduce our EDM results (Supplementary materials).

| Correlation analysis
Because we could not find any spatial trends in our spatial covariance analyses (results not shown), the regional average cross-correlations between every pair of the species are summarized in Figure 2. These results are qualitatively similar to those of Raimondo, Turcáni, et al. (2004) who analyzed the same data. The shaded pairs in Figure 2a indicate the species pairs with significant interspecific synchrony (significant regional average cross-correlation) and the darker shading indicates greater correlation values. As found by Raimondo, Turcáni, et al. (2004), interspecific synchrony was generally greater for pairs of species within the same feeding group (i.e., within the free-feeders and within the leaf-rollers). A hierarchical clustering based on the similarity of SCCF values in all scales supports this inference (R function "hclust" with method "average"). The lepidopteran species belonging to the free-feeder are grouped into relatively similar clusters compared to those of the leaf-rollers (Figure 2b).

| EDM analysis
Before the CCM and S-map analyses, we identified the embedding dimension and the nonlinearity for each Lepidoptera time series, by use of simplex projection and univariate S-map, respectively (the embedding dimension E i and the nonlinearity parameter θ i ranges from 1 to 15 and from 0.0 to 3.5, respectively; for more detail, see Table S1). Using the estimated E i , we performed the CCM analysis of the Lepidoptera time series between the meteorological time series (i.e., testing causality from meteorological factors to each lepidopteran species). The results demonstrated that the cross-map skill at the maximum data are significantly higher than both that at the minimum data and that of the surrogate one (p ≤ .05) for many cross-mapping pairs in the freefeeders (Figure 3a-f; 10 out of 18 pairs). The species Lymantria dispar was exceptional in that no convergence in cross-map skill was detected (Figure 3f). In contrast to the free-feeding species, there were no significant convergence in meteorological times series cross-mappings the leafrollers (0/12 pairs). That is, a significant causality of the meteorological time series was detected more frequently in the free-feeders than in the leaf-rollers (Fisher's exact test, p < .001). In addition, the cross-map skills at the maximum data were higher in the free-feeders than in the leaf-rollers (ANOVA, p < .001). In contrast, our additional analysis of meteorological time series cross-mapping Lepidopteran time series shows no causal effects of Lepidoptera on climatic factors ( Figure S1). The results of the CCM analysis between the Lepidoptera time series are summarized in Figure 4. The CCM analysis detected 10 significant causalities (p ≤ .05) from the 90 possible species pairs in total. More specifically, the significant causalities are distributed in 4/30 pairs in the cross-mapping between free-feeders, 3/24 pairs in the free-feeders crossmapping leaf-rollers, 2/24 pairs in the leaf-rollers crossmapping free-feeders and 1/12 pairs in the cross-map between leaf-rollers ( Figure 4). These differences in the cross-mapping pairs are not statistically significant (Fisher's exact test).
The results of the multivariate S-map analysis are summarized in Figure 5. The results demonstrated that the estimated effects of meteorological time series on the freefeeding group are different among the Lepidoptera species (Figure 4a-d). However, for a pair of Lepidoptera species with strong synchrony, that is, Erannis defoliaria, and Operophtera brumata Figure 2b), the signs (i.e., negative or positive) and effect sizes of the climate effects are similar for in all climatic factors. For Agriopsis aurantiaria, the effect of NAOI has a similar sign and effect size, but the other two factors have different signs for those of E. defoliaria and O. brumata. The results of multivariate S-maps for the species pairs with causal relationship are summarized in F I G U R E 3 Summary of the convergent cross mapping (CCM) analysis for each lepidopteran species cross-mapping the meteorological time series. The labels of X-axis, Temp, Prec, NAOI, represent the three meteorological factors, yearly mean temperature, precipitation and North Atlantic Oscillation Index, respectively. The red, blue and gray points and regions indicate the mean value and the 95% CI of the cross-map skill with the maximum library size, the minimum library size and the surrogate data, respectively. The panels (a)-(j) correspond to the CCM results of Agriopis aurantiaria, Alsophila aescularia, Erannis defoliaria, Operophtera brumata, Orthosia cruda, Lymantria dispar, Aleimma loefilingia, Archips xylosteana, Eudemis profundana and Totrix virida, respectively. The five species of the upper row and L. dispar are categorized into the free-feeder, and the other species are categorized into the leaf-roller group [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 4 Summary of the convergent cross mapping (CCM) analysis between lepidopteran time series. Predictive skills (Pearson's correlate coefficient, ρ) of species of Y-axis cross-mapping species of X-axis are plotted for all the pair of the species. Cross-map skills increase as the plotted color becomes darker. The dashed vertical-and horizontal-line indicate the border of the different feeding groups (species upper/left the lines are the free-feeders). The species IDs are same as Figure 2 Figure 5e. As shown in the figure, the estimated interaction effects tend to be larger in the species pairs categorized into the same feeding group (i.e., "Free-Free" or "Leaf-Leaf") than in the pairs categorized into different feeding groups ("Free-Leaf" or "Leaf-Free").

| DISCUSSION
Interspecific synchrony can arise via several possible mechanisms including common environmental effects, direct species interactions and shared trophic interactions, so that discriminating the relative importance of these causes is challenging. In this study, we tackled this issue by analyzing a time-series dataset of meteorological factors and the dynamics of several lepidopteran species using EDM, a nonlinear time series analysis method. The results indicate clear evidence that exogenous climatic conditions can be causal factors of interspecific synchrony among Lepidoptera species. Raimondo, Turcáni, et al. (2004) also employed the Lepidopteran data and found that there was interspecific synchrony within the same feeding groups. However, due to the lack of the data, they could not evaluate the relationship between climatic conditions and species dynamics in the framework of newly developed causality analyses but instead based their inference on correlation (not causal) analysis. We believe that our study adds useful new insight to studies of synchronous phenomena in ecology.
Specifically, the interspecific synchrony occurred more frequently in species pairs belonging to the same feeding group than between species from different groups (Figure 2). F I G U R E 5 Summary of multivariate S-map analyses. The top four panels are the results of S-map where multivariate embeddings of freefeeding species with causal meteorological factors are analyzed (a-d are those of Agriopis aurantiaria, Erannis defoliaria, Operophtera brumata and Orthosia cruda, respectively). The labels Temp, Prec and NAOI represent yearly mean temperature, precipitation, and North Atlantic Oscillation Index, respectively. The boxplots represent the distribution of S-map coefficients corresponding to the effects of each meteorological factor. (e) The results of lepidopteran species with causal relationships. The X-axis represents the type of species pairs analyzed by the multivariate S-maps (e.g., "Free-Free" indicates the estimated effects of a leaf-roller species on a free-feeding species) In addition, the CCM analysis detected causal effects of meteorological factors on several species categorized into the free-feeding group, while such effects were not observed in the leaf-rolling species (Figure 3). The multivariate S-map analysis demonstrated that the meteorological factors exert similar effects on a pair of free-feeding species with strong synchrony (E. defoliaria and O. brumata), and such agreements of the demographic effect would yield strong synchrony between the dynamics of different species (Figure 5b,c). Although we found that temperature and precipitation significantly affect the dynamics of A. aurantiaria, E. defoliaria and O. brumata, the multivariate S-map analysis showed that the effects of the two climatic factors on A. aurantiaria seem to be opposite to those of the other species (Figure 5a). The reason why A. aurantiaria showed the dynamics similar to the other free-feeding species and it requires further studies that leverage specific modeling or explore proximate mechanisms of climatic factors. In summary, these results suggest that meteorological factors are key drivers of interspecific synchrony of the free-feeding guild, but not for the leaf-rolling guild.
Combining correlation analysis between time series of different lepidopteran species and the exploration using a simple population model, Raimondo, Turcáni, et al. (2004) hypothesized that shared generalist predators could be causing the interspecific synchrony observed among foliagefeeding forest Lepidoptera species. Our EDM results do not necessarily rule out this hypothesis but rather reinforce its validity at least for the leaf-roller group. The foliage-feeding lepidopteran species exhibit relatively similar behavior and morphology within the same feeding group (Raimondo, Turcáni, et al., 2004). In particular, larvae of the free-feeding group are generally larger and are more exposed to exogenous environments, while leaf-rolling larvae are smaller and spend the majority of their time within their leaf shelter where they feed (Stehr, 1998). This biological information suggests that the free-feeding group may be affected by both weather conditions and shared predators while leaf-rolling species are not as strongly affected by environmental conditions. Therefore, we predicted that meteorological factors have stronger causal influences on the dynamics of freefeeders than the leaf-rollers and that interspecific synchrony is more prominent in the free-feeding group than the leaf-rollers: both the SCCF and CCM results are consistent with this prediction (Figures 2 and 3).
The EDM analysis also suggested the possible influence of biological factors on synchrony. If any pair of lepidopteran species share generalist predators, the dynamics of the two species may be linked by their shared nonlinear impacts, for example, apparent competition (Bonsall & Hassell, 1997;Holt, 1977) and switching predation with shared search image (Oaten & Murdoch, 1975). In fact, we found causal relationships in the five species pairs belonging to the same feeding group (Figure 4), and the estimated interaction effects of those causal relationships nonlinearly fluctuate ( Figure 5). Interestingly, of those causal relationships, the pair of leaf-rolling species, Eudemis profundana and Archips xylosteana, shows synchronized dynamics (Figure 2a) but this was apparently not caused by shared meteorological factor impacts (Figure 4). The CCM analysis also identified causal interactions between pairs of species that belong to different feeding group (e.g., E. profundana was affected by A. aurantiaria and E. defoliaria, Figure 4). These results may imply that the existence of generalist predators is involved in the interspecific synchrony. However, these results of course do not rule out the possibility that other synchronizing biological factors, such as direct competition between lepidopteran species and the effect of shared host plants. Application of the causality test to host and predator time series would be useful for clarifying these effects.
Our results demonstrate how CCM analysis can be used to identify unique features of individual species. Specifically, spatial cross-correlation analysis indicated that L. dispar was generally not synchronous with other Lepidoptera species (Figure 2). Further, causality test did not detect the causality from meteorological factors to L. dispar dynamics (Figure 3f). In fact, in spite of being a free-feeding insect, L. dispar is unique in that larvae are covered with dense setation (skin-hair like structure), possibly protecting them from many predators (Raimondo, Turcáni, et al., 2004). In addition, L. dispar is well known as an outbreaking species: populations periodically reach highdensity levels with outbreaks common in European, Asian and North American countries (Johnson, Liebhold, Bjørnstad, & McManus, 2005). However, the processes driving L. dispar outbreaks are not entirely understood (Liebhold, Elkinton, Williams, & Muzika, 2000;Turcáni, Novotny, & Zúbrik, 2001). Application of the EDM analysis to both lepidopteran and predator time series may resolve this problem. The CCM results here reveal that meteorological factors do not affect the dynamics of leaf-rolling species and L. dispar (Figure 3f), and thus suggest that predators or parasitoids are more likely mechanisms driving the population dynamics in these lepidopteran species. Thus, if L. dispar are free from strong predation pressures compared to the other leaf-roller species, this lack of strong predator/ parasitoid effects may contribute to the outbreak behavior of L. dispar dynamics. To test this possibility, further studies are needed that apply the EDM analysis to Lepidoptera system with time-series data of natural enemy species.
We note here that, although the causality tests were done with the rigorous criteria (i.e., cross-map skill of the maximum library size is compared to the minimum library's and the surrogate data's one), our EDM analyses suffer from the use of multiple comparison and some conclusions may be based upon false discovery of causalities. Application of simple correction methods, such as the Bonferroni correction, would result in excessively conservative alpha-levels because of the two causality criteria; ultimately, to avoid the problem of multiple comparison, it may be necessary to develop new criteria for significant causality. However, this problem would not qualitatively alter the main finding of this study: causal effects of meteorological conditions occur more frequently in the free-feeding group than in the leafroller group (Figure 3). Thus, the EDM approach is a useful way for the study of interspecific synchrony. The EDM may also be powerful in exploring the other ecological and practical issues. Originally, the EDM framework was developed to forecast nonlinear time series (Sugihara, 1994;Sugihara & May, 1990). In addition, recent studies reported that the method successfully identified critical interactions behind ecological phenomena from the observed time series (Deyle et al., 2016;Kawatsu & Kishi, 2018). Ushio et al. (2018) also applied the EDM framework to time-series data of aquatic fishes to measure the changeable community stability. Thus, like the Lepidoptera insects in our study, applying the EDM framework to the accumulated time-series dataset in ecology will pave the new way for understanding hidden ecological relationships, measuring changing community stability and forecasting future ecological system states.