Alternative stable states, nonlinear behavior, and predictability of microbiome dynamics

Microbiome dynamics are both crucial indicators and drivers of human health, agricultural output, and industrial bio-applications. However, predicting microbiome dynamics is notoriously difficult because communities often show abrupt structural changes, such as “dysbiosis” in human microbiomes. We here integrate theoretical and empirical bases for anticipating drastic shifts of microbial communities. We monitored 48 experimental microbiomes for 110 days and observed that various community-level events, including collapse and gradual compositional changes, occurred according to a defined set of environmental conditions. We then confirmed that the abrupt community changes observed through the time-series could be described as shifts between “alternative stable states” or dynamics around complex attractors. Furthermore, collapses of microbiome structure were successfully anticipated by means of the diagnostic threshold defined with the energy landscape analysis of statistical physics or that of a stability index of nonlinear mechanics. These results indicate that abrupt microbiome events in complex microbial communities can be forecasted by extending classic ecological concepts to the scale of species-rich microbial systems.

changes (see "abruptness" index in Fig. 1b) were observed in Water/Medium-A, 120 Soil/Medium-A, and Water/Medium-B treatments (abruptness > 0.5). Within these treatments, 121 taxonomic compositions and timing of abrupt shifts in community structure varied among 122 replicate communities (Extended Data Fig. 3). Large shifts of community compositions 123 through time were observed as well in Soil/Medium-B treatment, albeit the community shifts 124 were more gradual (maximum abruptness through time-series, 0.36 ~ 0.57; Extended Data 125 Fig. 3). In contrast, Medium-C condition yielded relatively steady microbiome dynamics with 126 continuously low taxonomic diversity (e.g., dominance of Aeromonas in Water/Medium-C 127 treatment), although shifts of dominant taxa were observed latter in the experiment in some 128 replicate communities (Extended Data Fig. 3). In all the six treatments, the a-diversity 129 (Shannon diversity) of ASVs displayed fluctuations, but not monotonic decrease, through 130 time (Extended Data Fig. 1e). 131 132 Framework 1: energy landscape analysis 133 By compiling the microbiome time-series data, we examined the distributions of stable states 134 within the multidimensional space of community structure based on an energy landscape 135 analysis 24 . Because no variation in environmental conditions was introduced through the 136 time-series in our experiment, a fixed "energy landscape" of community states was assumed 137 for each of the six treatments. On this assumption, shifts between alternative stable states are 138 attributed to perturbations to variables (i.e., population density of microbial ASVs) but not to 139 "regime shifts 34-36 ", which, by definition, requires energy landscape reorganization caused by 140 changes in environmental parameters (i.e., temperature). 141 In each experimental treatment, multiple stable states were estimated to exist (Fig. 2;142 Extended Data Fig. 4), indicating that the observed abrupt changes in community 143 compositions could be described as shifts between alternative stable states. Therefore, in this 144 approach of statistical physics 24-26 , community dynamics are divided into phases of 145 fluctuations around local equilibrium points and those of shifts into adjacent equilibria. In 146 other words, the presence of multiple equilibrium points (Extended Data Fig. 4) We next analyzed the time-series data based on the framework of empirical dynamic 152 modeling. We first focused on the population dynamics (increase/decrease) of the microbial 153 ASVs constituting the microbial communities. In ecology, population dynamics data have 154 often been analyzed with methods assuming linear dynamics (i.e., without considering "state 155 dependency 37 "). Meanwhile, a series of empirical dynamics modeling approaches applicable 156 to nonlinear time-series processes, such as simplex projection 20 and sequential locally 157 weighted global linear maps 19 (S-map), have been increasingly adopted to capture key 158 properties lost with linear dynamic assumptions (Fig. 3a). We found that ca. 85 % of the 159 microbial populations in our experiments exhibited nonlinear behavior (i.e., nonlinearity 160 parameter q > 0; Fig. 3b). This result suggests the predominance of nonlinear dynamics over 161 linear dynamics in microbial populations 32 , in line with populations of other organismal 162 groups such as fish 38 and plankton 21 . 163 We then reconstructed the attractors of nonlinear dynamics based on Takens' embedding 164 theorem 39 (Fig. 3a). To examine the performance of the attractor reconstruction, we conducted 165 forecasting of the population dynamics of respective microbial ASVs by means of simplex 166 projection and S-map (Fig. 3c). The population density (16S rRNA copy concentration) of an 167 ASV in a target replicate community at time point t + p (p represents time steps in 168 forecasting) was forecasted based on the ASV's population density at time point t and time-169 series data of other replicate communities (hereafter, reference replicate communities; see 170 Methods for details; Fig. 3a). For many microbial ASVs, predicted population densities was 171 positively correlated with observed ones Extended Data Fig. 5). As expected, 172 correlation between predicted and observed population size increased with increasing number 173 of reference replicate communities, suggesting dependence of forecasting skill on the size of 174 state-space reference databases (Extended Data Fig. 6). 175 By assembling the forecasting results of respective ASVs at the community level, we 176 further conducted forecasting of microbiome compositions ( Fig. 4a; Extended Data Fig. 7). 177 The forecasting precision of community-level dynamics varied depending on culture media 178 and the dissimilarity (b-diversity) of community structure between target and reference 179 replicates (Fig. 4b). Despite the utility of the forecasting platform, we observed high 180 prediction error immediately after the peak of abrupt community changes ( Based on the frameworks of the energy landscape analysis and empirical dynamic modeling, 192 we explored ways for anticipating abrupt events in community dynamics. In the former 193 framework roach, the reconstructed energy landscapes were used to estimate "energy gap" 194 and "stable-state entropy" indices, which represented stability/instability of community 195 states 24 (Fig. 5a). In the latter framework, the inferred Jacobian matrices of the multi-species 196 time-series dynamics (see Methods) were used to calculate "local Lyapunov stability 40 " and 197 "local structural stability 41 ". We examined how these indices could help us forecast large 198 community-compositional shifts such as those observed in Medium-A and Medium-B 199 treatments (Fig. 1b). 200 Among the signal indices examined, energy gap or stable-state entropy of community 201 states ( Fig. 5a) was significantly correlated with the degree of future community changes in 202 Medium-A and Medium-B treatments (FDR < 0.05 for all treatments; Fig. 5b; Extended Data 203 Fig. 9). In the seven-day-ahead forecasting of abrupt community-compositional changes 204 (abruptness > 0.5), for example, stable-state entropy showed relatively high diagnostic 205 performance on the two-dimensional surface of detection rate (sensitivity) and false detection 206 rate (1 -specificity) as represented by receiver operating characteristic (ROC) curve 42 . 207 Specifically, although the small number of points with abruptness greater than 0.5 (Extended 208 Data Fig. 10) precluded the application of the ROC analysis in Soil/Medium-B treatment, 209 diagnostic performance as evaluated by area under the ROC curve (AUC) ranged from 0.726 210 to 0.957 in other Medium-A and Medium-B treatments (Fig. 6a). 211 Local Lyapunov or structural stability was correlated with the degree of community 212 changes as well, but the correlations were less consistent among experimental treatments than 213 energy gap and stable-state entropy (Fig. 5b;Extended Data Fig. 9). Meanwhile, local 214 structural stability exhibited exceptionally high diagnostic performance in Water/Medium-A 215 stability can be used as signs of future microbiome collapse, although further technical 217 improvement in the state space reconstruction of species-rich communities (e.g., multi-view 218 distance regularized S-map 43 ) may be needed to gain consistent forecasting performance 219 across various types of microbiomes. 220 By further utilizing the frameworks of the energy landscape analysis and empirical 221 dynamic modeling, we next examined the availability of diagnostic thresholds for anticipating 222 community collapse. For this aim, we first focused on stable-state entropy because its 223 absolute values in the unit of well-known entropy index (Fig. 5a) were comparable across 224 diverse types of biological communities. Based on the ROC curve representing all the stable-225 state entropy data of Medium-A and Medium-B treatments, the balance between detection and 226 false-detection rates were optimized with the Youden index 42 . With a relatively high AUC 227 score (0.848), the threshold stable-state entropy was set as 1.343 (Fig. 6b). In the same way, 228 we calculated the threshold value for local Lyapunov stability because this index originally 229 had a tipping value (= 1) for diagnosing community-level stability/instability 40 . Indeed, the 230 estimated threshold of local Lyapunov stability on the ROC curve was 0.9802, close to the 231 theoretically expected value (Fig. 6b). 232 233

234
By compiling datasets of experimental microbiome dynamics under various environmental 235 (medium) conditions, we here tested whether two lines of ecological concepts could allow us 236 to anticipate drastic compositional changes in microbial communities. Despite decades-long 237 discussion on alternative stable/transient states of community structure [15][16][17]35,36 , the 238 application of the concept to empirical data of species-rich communities has been made 239 feasible only recently with the computationally intensive approach of statistical physics 240 (energy landscape analyses 24 ). On the other hand, the concept of dynamics around complex 241 forms of attractors has been applicable with the emerging framework of nonlinear 242 mechanics 27,40,41 (e.g., empirical dynamic modeling), microbiome time-series data satisfying 243 the requirements of the analytical frameworks remained scarce 32 . Thus, this study, which was 244 designed to apply both frameworks, provided a novel opportunity for fuel feedback between 245 empirical studies of species-rich communities and theoretical studies based on 246 classic/emerging ecological concepts. 247 Our analysis showed drastic events in microbiome dynamics, such as those observed in 248 dysbiosis of human-gut microbiomes 13,14 , could be forecasted, at least to some extent, by 249 framing microbiome time-series data as shifts between alternative stable states or dynamics 250 around complex attractors. In the forecasting of abrupt community changes observed in our 251 experimental microbiomes, the former concept (model) seemingly outperformed the latter 252 . This result is of particular interest, because concepts or models more efficiently 253 capturing dynamics of empirical data are expected to provide more plausible planforms in not 254 only prediction but also control of biological community processes. Nonetheless, given the 255 ongoing methodological improvements of nonlinear mechanics frameworks for describing 256 empirical time-series data 43 , further empirical studies comparing the two concepts are 257 necessary. 258 A key next step for forecasting and controlling microbial (and non-microbial) 259 community dynamics is to examine whether common diagnostic thresholds could be used to 260 anticipate collapse of community structure. This study provided the first empirical example 261 that the tipping value theoretically defined in noncolinear mechanics 40 (local Lyapunov 262 stability = 1) could be actually used as a threshold for alerting microbiome collapse. 263 Likewise, although estimates of diagnostic thresholds can vary depending on the definition of 264 community collapse (e.g., abruptness > 0.5 in this study), stable-state entropy scores greater 265 than 1.3 can be used to anticipate undesirable community events (dysbiosis) across medical, 266 agricultural, and industrial applications.  Fig. 1c-d). 486 The five standard DNAs were designed to be amplified with a primer set of 515f 50 and 487 806rB 51 but not to be aligned to the V4 region of any existing prokaryote 16S rRNA. Note 488 that the number of 16S rRNA copies per genome generally varies among prokaryotic taxa 52 489 and hence 16S rRNA copy concentration is not directly the optimal proxy of cell or biomass points t -4 to t and those from t + p to t + p + 4 (i.e., 5-day time-windows) were calculated 556 before evaluating degree of community changes for time point t and time step p in each 557 replicate community. Dissimilarity of community compositions between the time windows 558 before (from t -4 to t) and after (t + p to t + p + 4) each target time point t with a given time 559 lag p was calculated based on Bray-Curtis b-diversity as a measure of abrupt (sudden and 560 substantial) community changes (hereafter, "abruptness" of community-compositional 561 changes). A high value of this index indicates that abrupt community-compositional changes 562 occurred around time point t, while a low value represents a (quasi-)stable mode of 563 community dynamics. We also evaluated temporal changes of community compositions using 564 nonmetric multidimensional scaling (NMDS) with the R package vegan. 565 566

Energy landscape analysis 567
On the assumption that drastic changes in microbiome dynamics are described as shifts 568 between local equilibria (i.e., alternative stable states), we reconstructed the structure of the 569 "energy landscape 24,25 " in each experimental treatment (tutorials of energy landscape analyses 570 are available at https://community.wolfram.com/groups/-/m/t/2358581). Because external 571 environmental conditions (e.g., temperature) was kept constant in the experiment, a fixed 572 "energy landscape" of community states was assumed for each of the six experimental 573 treatments. Therefore, probabilities of community states p are given by 574 where = ( 4 (5) , 6 (5) , … , 8 (5) ) is a community state vector of k-th community state and S 577 is the total number of taxa (e.g., ASVs, species, or genera) examined. 9 (5) is a binary 578 variable that represents presence (1) or absence (0) of taxon i: i.e., there are a total of 2 8 579 community states. Then, the energy of the community state is given by 580 where ℎ A is the net effect of implicit abiotic factors, by which i-th taxon is more likely to 582 present (h i > 0) or not (h i < 0), and AC represents a co-occurrence pattern of i-th and j-th taxa. 583 Since the logarithm of the probability of a community state is inversely proportional to 584 ( ), a community state having lower H is more frequently observed. To consider 585 dynamics on an assembly graph defined as a network whose 2 8 nodes represent possible 586 community states and the edges represents transition path between them (two community 587 states are adjacent only if they have the opposite presence/absence status for just one species), 588 we assigned energy to nodes with the above equation, and so imposed the directionality in 589 state transitions. Then, we identified the stable state communities as the energy minima of the 590 weighted network (nodes having the lowest energy compared to all its neighbors), and 591 determined their basins of attraction based on a steepest descent procedure starting from each 592 node. The data of ASV-level compositions were used in the calculation of community state 593 energy using Mathematica v12.0.0.0. The "energy" estimates were then plotted against the 594 NMDS axes representing community states of the microbiome samples in each experimental 595 treatment. Spline smoothing of the landscape was performed with optimized penalty scores 596 using the mgcv v.1.8-40 package 61 of R. 597 598

Empirical dynamic modeling 599
In parallel with the energy landscape analysis assuming the presence of local equilibria, we 600 also analyzed the microbiome time-series data by assuming the presence of complex 601 attractors. In this aim, we applied the framework of "empirical dynamic modeling 19,20,29,40 ". In 602 general, biological community dynamics are driven by a number of variables (e.g., abundance 603 of respective species and abiotic environmental factors). In the analysis of a multi-variable 604 dynamic system in which only some of variables are observable, state space constituted by 605 time-lag axes of observable variables can represent the whole system as shown in Takens' 606 embedding theorem 39 . Thus, for each ASV in each experimental treatment, we conducted 607 Takens' embedding to reconstruct state space which consisted of time-delayed coordinates of 608 the ASV's absolute abundance (e.g., 16S rRNA copy concentration estimates). The optimal 609 number of embedding dimensions 29,39 (E) was obtained by finding E giving the smallest root-610 mean-square error (RMSE) in pre-run forecasting with simplex projection 20 or S-map 19 as 611 detailed below. Taking into account a previous study examining embedding dimensions 62 , 612 optimal E was explored within the range from 1 to 20. Prior to the embedding, all the 613 variables were z-standardized (i.e., zero-mean and unit-variance) across the time-series of 614 each ASV in each replicate community. 615 616

Population-level forecasting 617
Based on the state space reconstructed with Takens' embedding, simplex projection 20 was 618 applied for forecasting of ecological processes in our experimental microbiomes. For each 619 target replicate community, univariate embedding of each ASV was performed using the data 620 of the seven remaining replicate communities. Therefore, the reference databases for the 621 embedding did not include the information of the target replicate community (Fig. 2a), 622 providing platforms for evaluating forecasting skill. 623 In simplex projection, a coordinate within the reconstructed state space was explored at 624 a focal time point (t * ) within the population dynamics of a focal ASV in a target replicate 625 community (e.g., replicate community 8): the coordinate can be described as [x target_rep (t * ), 626 x target_rep (t * -1), x target_rep (t * -2)] when E = 3. For the focal coordinate, E + 1 neighboring 627 points are explored from the reference database consisting of the seven remaining replicate 628 communities (e.g., replicate communities 1-7; Fig. 2a). For each of the neighboring points, 629 the corresponding points at p-time-step forward (p-days ahead) are identified. The abundance 630 estimate of a focal ASV in the target replicate community at p-time-step forward [e.g., 631 target_rep (t * + p)] is then obtained as weighted average of the values of the highlighted p-time-632 step-forward points within the reference database (Fig. 2a). The weighting was decreased with 633 Euclidean distance between x target_rep (t * ) and its neighboring points within the reference 634 database. This forecasting of population dynamics was performed for each ASV in each target 635 replicate community at each time point. The number of time steps in the forecasting (i.e., ) 636 was set at 1 (one-day-ahead forecasting) and 7 (seven-day-ahead forecasting). Meanwhile, is an n ´ E dimensional matrix given by 661 The weighting function is defined as 663 where q is the parameter representing the nonlinearity of the data, while mean Euclidean 665 distance between reference database points and the target point in the target experimental 666 replicate is defined as follows: 667 where T ref denotes the set of A ¢ . In our analysis, the optimal value of q was explored among 669 eleven levels from 0 (linearity) and 8 (strong nonlinearity) for each ASV in each target 670 replicate based on the RMSE of forecasting (optimal q was selected among 0, 0.001, 0.01, 671 0.05, 0.1, 0.2, 0.5, 1, 2, 4, and 8). The local linear model was estimated for each time point 672 in the target replicate data. 673 We then performed direct comparison between linear and nonlinear approaches of 674 forecasting based on empirical dynamics modeling. Specifically, to assume linear dynamics in 675 S-map method, the nonlinearity parameter q was set 0 for all the ASVs. We then compared 676 forecasting results between linear (q = 0) and nonlinear (optimized q) approaches. For the 677 forecasting of ASVs in a target replicate community, the data of the remaining seven 678 communities (reference databases) were used as mentioned above. 679 For each ASV in each of the 48 experimental replicates, Spearman's correlation 680 coefficients between predicted and observed abundance (16S rRNA copy concentrations) 681 were calculated for each of the nonlinear/linear forecasting methods [simplex projection, S-682 map with optimized q, and S-map assuming linearity (q = 0)]. We also examined null model 683 assuming no change in community structure for a given time step. The time points (samples) 684 excluded in the data-quality filtering process (see Bioinformatics sub-section) were excluded 685 from the above evaluation of forecasting skill. 686 687

Reference database size and forecasting skill 688
To evaluate potential dependence of forecasting skill on the size of reference databases, we 689 performed a series of analyses with varying numbers of reference replicate communities. For 690 replicate community for a target replicate community, a fixed number (from 1 to 7) of other 691 replicate communities within each experimental treatment were retrieved as reference 692 databases: all combinations of reference communities were examined for each target replicate 693 community. For each microbial ASV in each target replicate community, forecasting of 694 population size was performed based on S-map with optimized q as detailed above. 695 Spearman's correlation between predicted and observed population size across the time-series 696 was then calculated for each ASV in each target replicate community. The correlation 697 coefficients were compared between different numbers of reference database communities 698 based on Welch's t-test in each experimental treatment. 699 700

Community-level forecasting 701
The above population-level results based on empirical dynamics modeling were then used for 702 forecasting community-level dynamics. For a focal time point (day) in a target experimental 703 replicate, the 16S rRNA copy concentration estimates (predicted abundance) of respective 704 ASVs were compiled, yielding predicted community structure (i.e., predicted relative 705 abundance of ASVs). The predicted and observed (actual) ASV compositions (relative 706 abundance) of respective target replicates were then plotted on a NMDS surface for each of 707 the six experimental treatments. In addition, we evaluated dependence of community-level 708 forecasting results on experimental conditions (source inocula and media), a-diversity 709 (Shannon's H¢) of ASVs, and mean b-diversity against other replicates in a multivariate 710 ANOVA model of predicated vs. observed community dissimilarity. 711 712 Anticipating abrupt community shifts 713 We then examined whether indices derived from the energy landscape analysis and/or 714 empirical dynamics modeling could be used to anticipate drastic changes in community 715 structure. 716 In the framework of energy landscape analysis, we calculated two types of indices 717 based on the estimated landscapes of microbiome dynamics (Fig. 3a). One is deviation of 718 current community-state energy from the possible lowest energy within the target basins 719 (hereafter, energy gap; Fig. 3a): this index represents how current community states are 720 inflated from local optima (i.e., "bottom" of basins). The other is "stable-state entropy 24 ", 721 which is calculated based on the random-walk-based simulation from current community 722 states to bottoms of any energy landscape basins (i.e., alternative stable states). A starting 723 community state is inferred to have high entropy if reached stable states are variable among 724 random-walk iterations: the stable-state entropy is defined as the Shannon's entropy of the 725 final destinations of the random walk 24 . Therefore, communities approaching abrupt structural 726 changes are expected to have high stable-state entropy because they are inferred to cross over 727 "ridges" on energy landscapes. For an analysis of a target replicate community, energy 728 landscapes were reconstructed based on the data of the remaining seven replicate 729 communities. 730 In the framework of empirical dynamics modeling (nonlinear mechanics), we calculated 731 "local Lyapunov stability 40 " (local dynamic stability) and "local structural stability 41 " based 732 on Jacobian matrices representing movements around reconstructed attractors 27 . Specifically, 733 based on convergent cross-mapping 22,32 and multivariate extension of S-map 63 , local 734 Lyapunov stability and structural stability were estimated, respectively, as the absolute value 735 of the dominant eigenvalue and trace (sum of diagonal elements) of the Jacobian matrices 736 representing the time-series processes 40 . For a target replicate community, the remaining 737 seven replicate communities were used for inferring Jacobian matrices. Note that a high score 738 of local Lyapunov/structural stability represents a potentially unstable community state. In 739 particular, local Lyapunov scores reflect whether trajectories at any particular time are 740 converging (local Lyapunov score < 1) or diverging (1 < local Lyapunov score) 40 . 741 For each of the above indices, linear regression of abruptness scores of community-742 compositional changes was performed for each replicate sample in each experimental 743 treatment (seven-day-ahead forecasting). The time points (samples) excluded in the data-744 quality filtering process (see Bioinformatics sub-section) were excluded from this evaluation 745 of signal indices. 746 We also examined the diagnostic performance of the signal indices based on the 747 receiver operating characteristic (ROC) analysis. In seven-day-ahead forecasting, detection 748 rates (sensitivity) and false detection rates (1 -specificity) of large community-compositional 749 changes (abruptness > 0.5) were plotted on a two-dimensional surface for each experimental 750 treatment, yielding area under the ROC curve (AUC) representing diagnostic performance 42 . 751 The optimal threshold value of each signal index for anticipating abrupt community-752 compositional changes (abruptness > 0.5) was then calculated with the Youden index 42 for 753 each experimental treatment. In addition, for stable-state entropy and local Lyapunov stability, 754 we calculated optimal threshold values after assembling all the data of Medium-A and 755 Medium-B treatments.       was forecasted with S-map (optimized q) based on reference databases (Fig. 2a). The 935 forecasting was performed for each number of reference databases defined on the horizontal 936 axis. Spearman's correlations between predicted and observed population size (Fig. 2c) were 937 calculated for each microbial ASV in each replicate community. An asterisk represents 938 significant differences in forecasting skill (forecasting performance) between different 939 numbers of reference databases in each experimental replicate: i.e., false discovery rate (FDR) 940 based on Welch's t-tests. Replicate community