 Research
 Open Access
 Published:
Alternative stable states, nonlinear behavior, and predictability of microbiome dynamics
Microbiome volume 11, Article number: 63 (2023)
Abstract
Background
Microbiome dynamics are both crucial indicators and potential drivers of human health, agricultural output, and industrial bioapplications. However, predicting microbiome dynamics is notoriously difficult because communities often show abrupt structural changes, such as “dysbiosis” in human microbiomes.
Methods
We integrated theoretical frameworks and empirical analyses with the aim of anticipating drastic shifts of microbial communities. We monitored 48 experimental microbiomes for 110 days and observed that various communitylevel events, including collapse and gradual compositional changes, occurred according to a defined set of environmental conditions. We analyzed the timeseries data based on statistical physics and nonlinear mechanics to describe the characteristics of the microbiome dynamics and to examine the predictability of major shifts in microbial community structure.
Results
We confirmed that the abrupt community changes observed through the timeseries 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.
Conclusions
The results indicate that abrupt microbiome events in complex microbial communities can be forecasted by extending classic ecological concepts to the scale of speciesrich microbial systems.
Background
Optimizing biological functions of speciesrich systems is a major challenge in both basic and applied sciences [1,2,3,4,5,6,7]. Managing the compositions of human gut microbiomes, for example, is essential for preventing diabetes [8, 9], infectious disease [10], and neuropsychiatric disorders [11]. Likewise, soil and plantassociated microbiomes drive nutrient cycling and pest/pathogen outbreaks in agroecosystems [5, 6], while highly controlled microbiomes facilitate stable and resourceefficient management in biofuel production [7] and water purification [12]. Nonetheless, it remains generally difficult to control microbial ecosystem functions because microbial communities often show drastic structural (compositional) changes [13, 14]. Thus, predicting such communityscale events remains an essential task for preventing unfavorable compositional changes and thereby keeping ecosystem functions in microbiome dynamics.
Drastic changes in biological community structure have been theoretically framed as transient dynamics towards a global equilibrium [15, 16], shifts between alternative equilibria [16, 17], or dynamics around complex forms of attractors [18,19,20]. Within a state space with a sole equilibrium point, drastic community compositional changes may be observed in the course of convergence to the global equilibrium [15]. In contrast, if multiple equilibria exist within a state space, abrupt community changes can be described as shifts between alternative stable states [17]. In other words, fluctuations in population densities of constituent species (variables) or changes in environments (parameters) can cause shifts of community states from a stable state to the other ones [16, 17]. Meanwhile, drastic community changes may be depicted as well in terms of dynamics around periodic/quasiperiodic attractors (i.e., limit cycle or torus) or dynamics around attractors with noninteger dimensions [18, 21,22,23] (i.e., chaos).
In analyzing empirical timeseries data of microbiome structure, these concepts of community dynamics are implemented with two lines of frameworks (Fig. 1a). One is the framework of energy landscape analyses in statistical physics [24,25,26], in which stability/instability of possible community states (compositions) are evaluated in the metric of “energy”. In energy landscape analyses, stable states within a state space are defined as community compositions whose energy values are lower than those of adjacent community compositions [24]. Thus, based on the reconstruction of energy landscapes, large community compositional changes are interpreted as transient dynamics towards an equilibrium or shifts between alternative equilibria (Fig. 1a). The other framework for describing abrupt community changes is based on nonlinear mechanics, which allows us to assume the presence of complex forms of attractors [19, 20, 22, 27]. The framework of empirical reconstruction of attractors (“empirical dynamic modeling [28, 29]”), in particular, provides a platform for interpreting community dynamics as deterministic processes around any forms of attractors (Fig. 1a). The two frameworks are potentially useful for framing microbial community processes. Nonetheless, it remains to be examined whether drastic changes in microbiome dynamics, such as dysbiosis in humanassociated microbiomes [14, 30, 31], could be predicted with either or both of the frameworks.
The major constraint preventing the comparison of the two frameworks is the lack of empirical datasets that simultaneously meet the basic requirements of energy landscape analyses and empirical dynamic modeling. In other words, although comparison of the two frameworks require information of populationsize dynamics of respective microbial species across tens or more of time points, majority of microbiome studies provide data of relative abundance for a small number of time points. Therefore, by developing a monitoring system of experimental microbiomes, we compile a series of microbiome timeseries data with substantial communitycompositional changes. By implementing an energy landscape analysis and empirical dynamic modeling, we examine whether the substantial community changes could be anticipated as transient dynamics towards global equilibria, shifts between stable states, or dynamics around complex attractors. Based on the results, we discuss how we can integrate empirical and theoretical studies for predicting and controlling speciesrich microbial systems.
Results
Experimental microbiome dynamics
To obtain timeseries datasets of microbiome dynamics, we constructed six types of microbiomes based on the combinations of two inoculum sources (soil and pond water microbiomes; hereafter, soil and water) and three media differing in chemical properties (oatmeal, oatmealpeptone, and peptone; hereafter, MediumA, B, and C, respectively), each with eight replicates (Additional file 1: Figure S1). We kept the experimental system at a constant temperature condition and sampled a fraction of each microbiome and added fresh media every 24 h for 110 days. For each of the six experimental treatment, 880 community samples were obtained (in total, 110 time points × 8 replicates × 6 treatments = 5280 community samples), providing rich information for exploring stable states of community structure by means of energy landscape analyses. In total, the dataset represented population dynamics of 264 prokaryote amplicon sequence variants (ASVs) belonging to 108 genera. Using quantitative amplicon sequencing [32] for estimating 16S ribosomal RNA gene (16S rRNA) copy concentrations of respective microbes in each microbiome, we determined the dynamics of both “relative” and “absolute (calibrated)” ASV abundance (Fig. 1b; Additional files 1, 2, 3 and 4: Figures S1–4). By estimating not only relative but also calibrated abundance, we were able to reconstruct respective ASVs’ population dynamics (increase/decrease), satisfying the requirements for applying empirical dynamic modeling [19, 20, 22].
The experimental microbiomes exhibited various types of dynamics depending on source inocula and culture media (Fig. 1b; Additional files 3 and 4; Figure S3–4). Specifically, sharp decline of taxonomic diversity [33] and abrupt (sudden and substantial) community structural changes (see “abruptness” index in Fig. 1b) were observed in Water/MediumA, Soil/MediumA, and Water/MediumB treatments (abruptness > 0.5). Within these treatments, taxonomic compositions and timing of abrupt shifts in community structure varied among replicate communities (Additional files 3 and 4: Figure S3–4). Large shifts of community compositions through time were observed as well in Soil/MediumB treatment, albeit the community shifts were more gradual (maximum abruptness through timeseries, 0.36~0.57; Additional files 3 and 4: Figure S3–4). In contrast, MediumC condition yielded relatively steady microbiome dynamics with continuously low taxonomic diversity (e.g., dominance of Aeromonas in Water/MediumC treatment), although shifts of dominant taxa were observed latter in the experiment in some replicate communities (Additional files 3 and 4: Figure S3–4). In all the six treatments, the αdiversity (Shannon diversity) of ASVs displayed fluctuations, but not monotonic decrease, through time (Additional file 1: Figure S1e).
Framework 1: energy landscape analysis
By compiling the microbiome timeseries data, we examined the distributions of stable states within the multidimensional space of community structure based on an energy landscape analysis [24]. Because no variation in environmental conditions was introduced through the timeseries in our experiment, a fixed “energy landscape” of community states was assumed for each of the six treatments. On this assumption, shifts between alternative stable states are attributed to perturbations to variables (i.e., population density of microbial ASVs) but not to “regime shifts [34,35,36]”, which, by definition, requires energy landscape reorganization caused by changes in environmental parameters (i.e., temperature).
Through the 110day dynamics of each experimental treatment, multiple stable states were inferred to exist (Fig. 2; Additional files 5 and 6: Figure S5–6). This result suggests that the observed abrupt changes in community compositions could be described as shifts between alternative stable states. Therefore, in this approach of statistical physics [24,25,26], community dynamics are divided into phases of fluctuations around local equilibrium points and those of shifts into adjacent equilibria. In other words, the presence of multiple equilibrium points (Additional files 5 and 6: Figure S5–6), by definition, means that the observed dynamics of the experimental microbiomes are not described as transient dynamics towards a sole equilibrium point.
Framework 2: empirical dynamic modeling
We next analyzed the timeseries data based on the framework of empirical dynamic modeling. We first focused on the population dynamics (increase/decrease) of the microbial ASVs constituting the microbial communities using the calibrated abundance data. In ecology, population dynamics data have often been analyzed with methods assuming linear dynamics (i.e., without considering “state dependency [37]”). Meanwhile, a series of empirical dynamic modeling approaches applicable to nonlinear timeseries processes, such as simplex projection [20] and sequential locally weighted global linear maps [19] (Smap), have been increasingly adopted to capture key properties lost with linear dynamic assumptions (Fig. 3a). We found that ca. 85% of the microbial populations in our experiments exhibited nonlinear behavior (i.e., nonlinearity parameter θ > 0; Fig. 3b). This result suggests the predominance of nonlinear dynamics over linear dynamics in microbial populations [32], in line with populations of other organismal groups such as fish [18] and plankton [21].
We then reconstructed the attractors of nonlinear dynamics based on Takens’ embedding theorem [38] (Fig. 3a; Additional file 7: Figure S7). To examine the performance of the attractor reconstruction, we conducted forecasting of the population dynamics of respective microbial ASVs by means of simplex projection and Smap (Fig. 3c). The population density (16S rRNA copy concentration) of an ASV in a target replicate community at time point t + p (p represents time steps in forecasting) was forecasted based on the ASV’s population density at time point t and timeseries data of other replicate communities (hereafter, reference replicate communities; see “Methods” section for details; Fig. 3a). For many microbial ASVs, predicted population densities were positively correlated with observed ones, although prediction accuracy decayed with time steps in forecasting (Fig. 3c, d; Additional file 8: Figure S8). As expected, correlation between predicted and observed population size increased with increasing number of reference replicate communities, suggesting dependence of forecasting skill on the size of statespace reference databases (Additional file 9: Figure S9).
By assembling the forecasting results of respective ASVs at the community level, we further conducted forecasting of microbiome compositions (Fig. 4a; Additional files 10 and 11: Figure S10–11). The forecasting precision of communitylevel dynamics varied depending on inoculum, culture media, αdiversity (Shannon’s H ′), and the dissimilarity (βdiversity) of community structure between target and reference replicates (Fig. 4b). Despite the utility of the forecasting platform, we observed high prediction error immediately after the peak of abrupt community changes (Fig. 4c; Additional file 12: Figure S12). Although the nonlinear method (Smap with optimized θ) captured the observed abrupt shifts of community compositions within a narrower time window than the linear method (Smap with θ = 0) (Fig. 4c), quantitative forecasting of abrupt community changes seemed inherently difficult.
Nonetheless, even if precise forecasting of community compositional dynamics remains challenging, prediction of the occurrence of abrupt community changes per se may be possible. Thus, we next examined whether potential of abrupt community changes could be evaluated through microbiome dynamics.
Anticipating abrupt community shifts
Based on the frameworks of the energy landscape analysis and empirical dynamic modeling, we explored ways for anticipating abrupt events in community dynamics. In the former framework, the reconstructed energy landscapes were used to estimate “energy gap” and “stablestate entropy” indices, which represented stability/instability of community states [24] (Fig. 5a). In the latter framework, the inferred Jacobian matrices of the multispecies timeseries dynamics (see “Methods” section) were used to calculate “local Lyapunov stability [39]” and “local structural stability [40]”. We examined how these indices could help us forecast large communitycompositional shifts such as those observed in MediumA and MediumB treatments (Fig. 1b).
Among the signal indices examined, energy gap or stablestate entropy of community states (Fig. 5a) was significantly correlated with the degree of future community changes in MediumA and MediumB treatments (FDR < 0.05 for all treatments; Fig. 5b; Additional files 13 and 14: Figure S13–14). In the 7dayahead forecasting of abrupt communitycompositional changes (abruptness > 0.5), for example, stablestate entropy showed relatively high diagnostic performance on the twodimensional surface of detection rate (sensitivity) and false detection rate (1–specificity) as represented by receiver operating characteristic (ROC) curve [41]. Specifically, although the small number of points with abruptness greater than 0.5 (Additional file 15: Figure S15) precluded the application of the ROC analysis in Soil/MediumB treatment, diagnostic performance as evaluated by area under the ROC curve (AUC) ranged from 0.726 to 0.957 in other MediumA and MediumB treatments (Fig. 6a).
Local Lyapunov or structural stability was correlated with the degree of community changes as well, but the correlations were less consistent among experimental treatments than energy gap and stablestate entropy (Fig. 5b; Additional file 13 and 14: Figure S13–14). Meanwhile, local structural stability exhibited exceptionally high diagnostic performance in Water/MediumA treatment (AUC = 0.788; Fig. 6a; Additional file 15: Figure S15). Thus, local Lyapunov or structural stability can be used as signs of future microbiome collapse, although further technical improvement in the state space reconstruction of speciesrich communities (e.g., multiview distance regularized Smap [42]) may be needed to gain consistent forecasting performance across various types of microbiomes.
By further utilizing the frameworks of the energy landscape analysis and empirical dynamic modeling, we next examined the availability of diagnostic thresholds for anticipating community collapse. For this aim, we first focused on stablestate entropy because its absolute values in the unit of wellknown entropy index (Fig. 5a) were comparable across diverse types of biological communities. Based on the ROC curve representing all the stablestate entropy data of MediumA and MediumB treatments, the balance between detection and falsedetection rates were optimized with the Youden index [41]. With a relatively high AUC score (0.848), the threshold stablestate entropy was set as 1.343 (Fig. 6b). In the same way, we calculated the threshold value for local Lyapunov stability because this index originally had a tipping value (= 1) for diagnosing communitylevel stability/instability [39]. Indeed, the estimated threshold of local Lyapunov stability on the ROC curve was 0.9802, close to the theoretically expected value (Fig. 6b).
Discussion
By compiling datasets of experimental microbiome dynamics under various environmental (medium) conditions, we here tested whether two lines of ecological concepts could allow us to anticipate drastic compositional changes in microbial communities. Despite decadeslong discussion on alternative stable/transient states of community structure [15,16,17, 35, 36], the application of the concept to empirical data of speciesrich communities has been made feasible only recently with the computationally intensive approach of statistical physics (energy landscape analyses [24]). On the other hand, the concept of dynamics around complex forms of attractors has been applicable with the emerging framework of nonlinear mechanics [27, 39, 40] (e.g., empirical dynamic modeling), microbiome timeseries data satisfying the requirements of the analytical frameworks remained scarce [32]. Thus, this study, which was designed to apply both frameworks, provided a novel opportunity for fueling feedback between empirical studies of speciesrich communities and theoretical studies based on classic/emerging ecological concepts.
Our analysis showed drastic events in microbiome dynamics, such as those observed in dysbiosis of humangut microbiomes [13, 14], could be forecasted, at least to some extent, by framing microbiome timeseries data as shifts between alternative stable states or dynamics around complex attractors. In the forecasting of abrupt community changes observed in our experimental microbiomes, the former concept (model) seemingly outperformed the latter (Figs. 5 and 6). This result is of particular interest, because concepts or models more efficiently capturing dynamics of empirical data are expected to provide more plausible planforms in not only prediction but also control of biological community processes. Nonetheless, given the ongoing methodological improvements of nonlinear mechanics frameworks for describing empirical timeseries data [42], further empirical studies comparing the two concepts are necessary.
A key next step for forecasting and controlling microbial (and nonmicrobial) community dynamics is to examine whether common diagnostic thresholds could be used to anticipate collapse of community structure. This study provided the first empirical example that the tipping value theoretically defined in noncolinear mechanics [39] (local Lyapunov stability = 1) could be actually used as a threshold for alerting abrupt changes in microbiome structure. Likewise, although estimates of diagnostic thresholds can vary depending on the definition of major microbiome shifts (e.g., abruptness > 0.5 in this study), stablestate entropy scores greater than 1.3 can be used to anticipate undesirable community events (dysbiosis) across medical, agricultural, and industrial applications.
Given that changes in environmental parameters were not incorporated into our experimental design, it remains another important challenge to reveal how distributions of stable states or forms of attractors are continually reshaped by changes in environmental parameters through community dynamics [17, 34, 35]. Experimental manipulation of “external” environmental parameters in microcosms, for example, will expand the target of research into microbiome systems potentially driven by regime shifts [34,35,36]. Likewise, environmental alternations caused by organisms per se [43,44,45] deserve further investigations as potential drivers of drastic community shifts. It is also important to explore potential difference in timeseries dynamics between real microbiomes and closed experimental systems. For example, while human gut microbiomes are considered to have “normal” or “standard” community states of healthy host individuals [3, 13, 46], variability of community dynamics in experimental microbiomes has just started to be explored [33, 43]. More comparative studies across diverse types of microbiomes are necessary for deepening our understanding of microbiome dynamics.
Conclusions
We showed that large shifts in microbial community structure is predictable based on emerging approaches of statistical physics and nonlinear mechanics, although further updates of mathematical and informatics platforms are necessary for increasing forecasting accuracy. We also found that abrupt microbiome changes could be interpreted as both shifts between alternative stable states and dynamics around complex attractors, while more empirical studies are required to discuss which model is more suitable for describing microbiome dynamics. Controlling biological functions at the ecosystem level is one of the major scientific challenges in the twentyfirst century [5, 47, 48]. Interdisciplinary approaches that further integrate microbiology, ecology, and mathematics are becoming indispensable for maximizing and stabilizing microbiomelevel functions, and for providing novel solutions to a broad range of humanity issues spanning from human health to sustainable industry and food production.
Methods
Continuousculture of microbiome
To set up experimental bacterial communities, we prepared two types of source inocula (soil and aquatic microbiomes) and three media (oatmeal, oatmealpeptone, and peptone): for each combination of source media and inocula (six experimental treatments), eight replicate communities were established (in total, two source microbiomes × three media × eight replicates = 48 experimental communities; Additional file 1: Figure S1a). We used natural microbial communities, rather than “synthetic” communities with predefined diversity, as source microbiomes of the experiment. One of the source microbiomes derives from the soil collected from the A layer (0–10 cm in depth) in the research forest of Center for Ecological Research, Kyoto University, Kyoto, Japan (34.972 °N; 135.958 °E). After sampling, 600 g of the soil was sieved with a 2mm mesh and then 5 g of the sieved soil was mixed in 30 mL autoclaved distilled water. The source microbiome was further diluted 10 times with autoclaved distilled water. The source aquatic microbiome was prepared by collecting 200 mL of water from a pond (“Shoubuike”) near Center for Ecological Research (34.974 °N, 135.966 °E). In the laboratory, 3 mL of the collected water was mixed with 27 mL of distilled water in a 50mL centrifuge tube. We then introduced the source soil or aquatic microbiomes into three types of media: oatmeal broth [0.5% (w/v) milled oatmeal (Nisshoku Oats; Nippon Food Manufacturer); MediumA], oatmealpeptone broth [0.5% (w/v) milled oatmeal + 0.5% (w/v) peptone (Bacto Peptone; BD; lot number: 7100982); MediumB], and peptone broth [0.5% (w/v) peptone; MediumC]. In our preliminary experiments, microbiomes cultured with MediumA (oatmeal) tended to show high species diversity, while those cultured with MediumC (peptone) were constituted by smaller number of bacterial species. Thus, we expected that diverse types of microbiome dynamics would be observed with the three medium conditions. Among the three media, MediumB had the highest concentrations of nonpurgeable organic carbon (NPOC) and total nitrogen (TN), while MediumA was the poorest both in NPOC and TN: MediumC had the intermediate properties (Additional file 1: Figure S1b).
In each well of a 2mL deep well plate, 200 μL of a diluted source microbiome and 800 μL of medium were installed. The deepwell plate was kept shaken at 1000 rpm using a microplate mixer NS4P (AS ONE Corporation, Osaka) at 23 °C for 5 days. After the 5day preincubation, 200 μL out of 1000μL culture medium was sampled from each of the 48 deep wells after mixing (pipetting) every 24 h for 110 days. In each sampling event, 200 μL of fresh medium was added to each well so that the total culture volume was kept constant. In total, 5280 samples (48 communities/day × 110 days) were collected. Note that on day 82, 200μL of fresh MediumB was accidentally added to samples of MediumA but not to those of MediumB. While the microbiomes under MediumA treatments experienced increase in total DNA copy concentrations late in the timeseries, relative abundance remained relatively constant from day 60 to 110 (Additional files 2, 3 and 4: Figure S2–4), suggesting limited impacts of the accidental addition of the medium on microbial community compositions.
To extract DNA from each sample, 25 μL of the collected aliquot was mixed with 50 μL lysis buffer (0.0025% SDS, 20 mM Tris (pH 8.0), 2.5 mM EDTA, and 0.4 M NaCl) and proteinase K (× 1/100). The mixed solution was incubated at 37 °C for 60 min followed by 95 °C for 10 min and then the solution was vortexed for 10 min to increase DNA yield.
Quantitative 16S rRNA sequencing
To reveal the increase/decrease of population size for each microbial ASV, a quantitative amplicon sequencing method [32, 49] was used based on Illumina sequencing platform. While most existing microbiome studies were designed to reveal the “relative” abundance of microbial ASVs or operational taxonomic units (OTUs), analyses of population dynamics, in principle, require the timeseries information of “absolute” abundance. In our quantitative amplicon sequencing, five standard DNA sequence variants with different concentrations of artificial 16S rRNA sequences (0.1, 0.05, 0.02, 0.01, and 0.005 nM) were added to PCR master mix solutions (Additional file 1: Figure S1a). The DNA copy concentration gradient of the standard DNA variants yielded calibration curves between Illumina sequencing read numbers and DNA copy numbers (concentrations) of the 16S rRNA region in PCR reactions, allowing estimation of original DNA concentrations of target samples [32, 49] (Additional file 1: Figure S1cd). The five standard DNAs were designed to be amplified with a primer set of 515f [50] and 806rB [51] but not to be aligned to the V4 region of any existing prokaryote 16S rRNA. Note that the number of 16S rRNA copies per genome generally varies among prokaryotic taxa [52] and hence 16S rRNA copy concentration is not directly the optimal proxy of cell or biomass concentration. However, in our study, estimates of 16S rRNA copy concentrations are used to monitor increase/decrease of abundance (i.e., population dynamics) within the timeseries of each microbial ASV. Even if the concentrations of PCR inhibitor molecules in DNA extracts vary among timeseries samples, potential bias caused by such inhibitors can be corrected based on the abovementioned method using internal standards (i.e., standard DNAs within PCR master solutions).
Prokaryote 16S rRNA region was PCRamplified with the forward primer 515f fused with 3–6mer Ns for improved Illumina sequencing quality and the forward Illumina sequencing primer (5′TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG[3–6mer Ns]–[515f]3′) and the reverse primer 806rB fused with 3–6mer Ns for improved Illumina sequencing quality [53] and the reverse sequencing primer (5′GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G [3–6mer Ns][806rB]3′) (0.2 μM each). The buffer and polymerase system of KOD One (Toyobo) was used with the temperature profile of 35 cycles at 98 °C for 10 s, 55 °C for 30 s, 68 °C for 50 s. To prevent generation of chimeric sequences, the ramp rate through the thermal cycles was set to 1 °C/s [54]. Illumina sequencing adaptors were then added to respective samples in the supplemental PCR using the forward fusion primers consisting of the P5 Illumina adaptor, 8mer indexes for sample identification [55] and a partial sequence of the sequencing primer (5′AAT GAT ACG GCG ACC ACC GAG ATC TAC AC[8mer index]TCG TCG GCA GCG TC3′) and the reverse fusion primers consisting of the P7 adaptor, 8mer indexes, and a partial sequence of the sequencing primer (5′ CAA GCA GAA GAC GGC ATA CGA GAT[8mer index]GTC TCG TGG GCT CGG3′). KOD One was used with a temperature profile: followed by 8 cycles at 98 °C for 10 s, 55 °C for 30 s, 68 °C for 50 s (ramp rate = 1 °C/s). The PCR amplicons of the samples were then pooled after a purification/equalization process with the AMPureXP Kit (Beckman Coulter). Primer dimers, which were shorter than 200 bp, were removed from the pooled library by supplemental purification with AMpureXP: the ratio of AMPureXP reagent to the pooled library was set to 0.6 (v/v) in this process. The sequencing libraries were processed in an Illumina MiSeq sequencer (271 forward (R1) and 231 reverse (R4) cycles; 10% PhiX spikein).
Bioinformatics
In total, 67,537,480 sequencing reads were obtained in the Illumina sequencing. The raw sequencing data were converted into FASTQ files using the program bcl2fastq 1.8.4 distributed by Illumina. The output FASTQ files were then demultiplexed with the program Claident v0.2. 2018.05.2 9[56]. The sequencing reads were subsequently processed with the program DADA2 [57] v.1.13.0 of R 3.6.0 to remove lowquality data. The molecular identification of the obtained ASVs was performed based on the naive Bayesian classifier method [58] with the SILVA v.132 database [59]. In total, 399 prokaryote (bacterial or archaeal) ASVs were detected. We obtained a sample × ASV matrix, in which a cell entry depicted the concentration of 16S rRNA copies of an ASV in a sample. In this process of estimating original DNA copy numbers (concentrations) of respective ASVs from sequencing read numbers in each sample, the samples in which Pearson’s coefficients of correlations between sequencing read numbers and standard DNA copy numbers (i.e., correlation coefficients representing calibration curves) were less than 0.7 (in total, 430 samples out of 5280 samples) were removed as those with unreliable estimates. Samples with less than 350 reads were discarded as well. Because missing values within timeseries data are not tolerated in some of the downstream analyses (e.g., empirical dynamic modeling), they were substituted by interpolated values, which were obtained as means of the time points immediately before and after focal missing time points. The ASVs that appeared 5 or more samples in any of the replicate communities were retained in the following analyses: 264 ASVs representing 108 genera remained in the dataset. From the sample × ASV matrix, we calculated αdiversity (Shannon’s H′) of the ASV compositions in each experimental replicate on each day. We also evaluated dissimilarity of community compositions in all pairs of sampling days in each replicate community using BrayCurtis metric of βdiversity as implemented in the vegan 2.5.5 package [60] of R. For each ASV in each replicate community, a parameter representing the nonlinearity of the population dynamics [18, 19] (θ) was estimated based on Smap analysis of calibrated abundance as detailed below in order to evaluate the overall nature of the timeseries data.
Community dynamics
We evaluated the degree of communitycompositional changes for time point t based on the BrayCurtis βdiversity through time. To remove effects of minor fluctuations and track only fundamental changes of community structure, average community compositions from time points t – 4 to t and those from t + p to t + p + 4 (i.e., 5day timewindows) were calculated before evaluating degree of community changes for time point t and time step p in each replicate community. Dissimilarity of community compositions between the time windows before (from t – 4 to t) and after (t + p to t + p + 4) each target time point t with a given time lag p was calculated based on BrayCurtis βdiversity as a measure of abrupt (sudden and substantial) community changes (hereafter, “abruptness” of communitycompositional changes). A high value of this index indicates that abrupt communitycompositional changes occurred around time point t, while a low value represents a (quasi)stable mode of community dynamics. We also evaluated temporal changes of community compositions using nonmetric multidimensional scaling (NMDS) with the R package vegan.
Energy landscape analysis
On the assumption that drastic changes in microbiome dynamics are described as shifts between local equilibria (i.e., alternative stable states), we reconstructed the structure of the “energy landscape [24, 25]“ in each experimental treatment (tutorials of energy landscape analyses are available at https://community.wolfram.com/groups//m/t/2358581). Because external environmental conditions (e.g., temperature) was kept constant in the experiment, a fixed “energy landscape” of community states was assumed for each of the six experimental treatments. Therefore, probabilities of community states \(\textrm{p}\left({\overrightarrow{\sigma}}^{\left(k\right)}\right)\) are given by
where \({\overrightarrow{\sigma}}^{\left(k\right)}=\left({\sigma}_1^{(k)},{\sigma}_2^{(k)},\dots, {\sigma}_S^{(k)}\right)\) is a community state vector of kth community state and S is the total number of taxa (e.g., ASVs, species, or genera) examined. \({\sigma}_{i}^{(k)}\) is a binary variable that represents presence (1) or absence (0) of taxon i: i.e., there are a total of 2^{S} community states. Then, the energy of the community state \({\overrightarrow{\boldsymbol{\sigma}}}^{\left(k\right)}\) is given by
where h_{i} is the net effect of implicit abiotic factors, by which ith taxon is more likely to present (h_{i} > 0) or not (h_{i} < 0), and J_{ij} represents a cooccurrence pattern of ith and jth taxa. Since the logarithm of the probability of a community state is inversely proportional to \(H\left({\overrightarrow{\sigma}}^{\left(k\right)}\right)\), a community state having lower H is more frequently observed. To consider dynamics on an assembly graph defined as a network whose 2^{S} nodes represent possible community states and the edges represents transition path between them (two community states are adjacent only if they have the opposite presence/absence status for just one species), we assigned energy to nodes with the above equation, and so imposed the directionality in state transitions. The parameters h_{i} and J_{ij} were explored based on a gradient descent algorithm [25, 61, 62]. For a model with parameters h^{∗} and J^{∗}, the expected probability of taxon i and that of cooccurrence are expressed, respectively, as
By iteratively adjusting \({\left\langle {\overrightarrow{\boldsymbol{\sigma}}}_i\right\rangle}^{\ast }\) and \({\left\langle {\overrightarrow{\boldsymbol{\sigma}}}_i{\overrightarrow{\boldsymbol{\sigma}}}_j\right\rangle}^{\ast }\) to the mean occurrence and cooccurrence calculated from the observational data [i.e., \(\left\langle {\overrightarrow{\boldsymbol{\sigma}}}_i\right\rangle\) and \(\left\langle {\overrightarrow{\boldsymbol{\sigma}}}_i{\overrightarrow{\boldsymbol{\sigma}}}_j\right\rangle\)], the parameters h_{i} and J_{ij} were estimated [24]. At each step, the parameters were updated as
where the learning rate α = 0.1 [24]. The maximum number of iterations was set to 50,000.
Based on the above equations, we identified the stable state communities as the energy minima of the weighted network (nodes having the lowest energy compared to all its neighbors), and determined their basins of attraction based on a steepest descent procedure starting from each node. Note that extinctions of species can generate “forbidden” directed paths between community states, potentially causing bias in the reconstruction of energy landscapes: treating extinctions in closed ecosystem dynamics is an important future challenge in energy landscape analyses. The data of ASVlevel compositions were used in the calculation of community state energy using Mathematica v12.0.0.0. Because ASVs appearing in majority of samples and those appearing in only a small number of samples are uninformative in the energy landscape analysis, the ASVs detected from 2–98% of samples were targeted. The data matrices of calibrated abundance were then converted into binary (1 or 0) format prior to energy landscape analyses. Thus, the value 0 in the input binary data does not necessarily mean complete absence of a taxon. By setting thresholds for binarizing the input data, energy landscape analyses can explore transitions between potential stable states. The “energy” estimates were then plotted against the NMDS axes representing community states of the microbiome samples in each experimental treatment. Spline smoothing of the landscape was performed with optimized penalty scores using the mgcv v.1.840 package [63] of R.
Empirical dynamic modeling
In parallel with the energy landscape analysis assuming the presence of local equilibria, we also analyzed the microbiome timeseries data by assuming the presence of complex attractors. In this aim, we applied the framework of “empirical dynamic modeling [19, 20, 29, 39]”. In general, biological community dynamics are driven by a number of variables (e.g., abundance of respective species and abiotic environmental factors). In the analysis of a multivariable dynamic system in which only some of variables are observable, state space constituted by timelag axes of observable variables can represent the whole system as shown in Takens’ embedding theorem [38]. Thus, for each ASV in each experimental treatment, we conducted Takens’ embedding to reconstruct state space which consisted of timedelayed coordinates of the ASV’s calibrated abundance (e.g., 16S rRNA copy concentration estimates). The optimal number of embedding dimensions [29, 38] (E) was obtained by finding E giving the smallest rootmeansquare error (RMSE) in prerun forecasting with simplex projection [20] or Smap [19] as detailed below. Taking into account a previous study examining embedding dimensions [64], optimal E was explored within the range from 1 to 20. Prior to the embedding, all the variables were zstandardized (i.e., zeromean and unitvariance) across the timeseries of each ASV in each replicate community.
Populationlevel forecasting
Based on the state space reconstructed with Takens’ embedding, simplex projection [20] was applied for forecasting of ecological processes in our experimental microbiomes. For each target replicate community, univariate embedding of each ASV was performed using the data of the seven remaining replicate communities. Therefore, the reference databases for the embedding did not include the information of the target replicate community (Fig. 3a), providing platforms for evaluating forecasting skill.
In simplex projection, a coordinate within the reconstructed state space was explored at a focal time point (t^{*}) within the population dynamics of a focal ASV in a target replicate community (e.g., replicate community 8): as time delay was set to 1 in our analysis, the coordinate can be described as [x_{target_rep}(t^{*}), x_{target_rep}(t^{*} – 1), x_{target_rep}(t^{*} – 2)] when E = 3. For the focal coordinate, E + 1 neighboring points are explored from the reference database consisting of the seven remaining replicate communities (e.g., replicate communities 1–7; Fig. 3a). For each of the neighboring points, the corresponding points at ptimestep forward (pdays ahead) are identified. The abundance estimate of a focal ASV in the target replicate community at ptimestep forward [e.g., \(\hat{x}\)_{target_rep}(t^{*} + p)] is then obtained as weighted average of the values of the highlighted ptimestepforward points within the reference database (Fig. 3a). The weighting was decreased with Euclidean distance between x_{target_rep}(t^{*}) and its neighboring points within the reference database. This forecasting of population dynamics was performed for each ASV in each target replicate community at each time point. The number of time steps in the forecasting (i.e., p) was set within the range from 1 (1dayahead forecasting) to 7 (7dayahead forecasting).
While simplex projection explores neighboring points around a target point, Smap [19] uses all the data points after weighting contributions of each point within a reference database using a parameter representing nonlinearity of the system. In Takens’ embedding, the state space of a target replicate for forecasting at time t is defined as
where E is embedding dimension. Values on the second and higher dimensions {z_{2, target _ rep}(t), …, z_{E, target _ rep}(t)} are represented by timedelayed coordinates of a focal ASV. Likewise, the state space of the remaining replicates (i.e., the reference database) is defined as
where t′ represents each of nonoverlapping time points within the reference database [e.g., {10001, 10002, …, 10110} and {20001, 20002, …, 20110} for reference replicate 1 and 2, respectively]. For a target time point t^{∗} within the timeseries data of a target replicate community, a local linear model C is produced to predict the future abundance of a focal ASV [i.e., z_{1, target _ rep}(t^{∗} + p)] from the statespace vector at time point z_{target _ rep}(t^{∗}) as follows:
This linear model is fit to the vectors in the reference databases. In the regression analysis, data points close to the target point z_{target _ rep}(t^{∗}) have greater weighting. The model C is then the singular value decomposition solution to the equation b = AC. In this equation, b is set as an ndimensional vector of the weighted future values of z_{1, ref}(t_{i}^{′} ) for each point (t_{i}^{′}) in the reference database (n is the number of points in the set of t_{i}^{′}): i.e.,
where ・ is Euclidean distance between two points in an E dimensional space. Meanwhile, A is an n × E dimensional matrix given by
The weighting function w is defined as
where θ is the parameter representing the nonlinearity of the data, while mean Euclidean distance between reference database points and the target point in the target experimental replicate is defined as follows:
where T_{ref} denotes the set of t_{i}^{′} . In our analysis, the optimal value of θ was explored among eleven levels from 0 (linearity) and 8 (strong nonlinearity) for each ASV in each target replicate based on the RMSE of forecasting (optimal θ was selected among 0, 0.001, 0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 4, and 8). The local linear model C was estimated for each time point in the target replicate data.
We then performed direct comparison between linear and nonlinear approaches of forecasting based on empirical dynamic modeling. Specifically, to assume linear dynamics in Smap method, the nonlinearity parameter θ was set 0 for all the ASVs. We then compared forecasting results between linear (θ = 0) and nonlinear (optimized θ) approaches. For the forecasting of ASVs in a target replicate community, the data of the remaining seven communities (reference databases) were used as mentioned above.
For each ASV in each of the 48 experimental replicates, R^{2} values between predicted and observed abundance (16S rRNA copy concentrations) were calculated for each of the nonlinear/linear forecasting methods [simplex projection, Smap with optimized θ, and Smap assuming linearity (θ = 0)]. We also examined null model assuming no change in community structure for a given time step. The time points (samples) excluded in the dataquality filtering process (see “Bioinformatics” section) were excluded from the above evaluation of forecasting skill.
Reference database size and forecasting skill
To evaluate potential dependence of forecasting skill on the size of reference databases, we performed a series of analyses with varying numbers of reference replicate communities. For replicate community for a target replicate community, a fixed number (from 1 to 7) of other replicate communities within each experimental treatment were retrieved as reference databases: all combinations of reference communities were examined for each target replicate community. For each microbial ASV in each target replicate community, forecasting of population size was performed based on Smap with optimized θ as detailed above. R^{2} values between predicted and observed population size across the timeseries was then calculated for each ASV in each target replicate community. The correlation coefficients were compared between different numbers of reference database communities based on Welch’s ttest in each experimental treatment.
Communitylevel forecasting
The above populationlevel results based on empirical dynamic modeling were then used for forecasting communitylevel dynamics. For a focal time point (day) in a target experimental replicate, the 16S rRNA copy concentration estimates (predicted abundance) of respective ASVs were compiled, yielding predicted community structure (i.e., predicted relative abundance of ASVs). The predicted and observed (actual) ASV compositions (relative abundance) of respective target replicates were then plotted on a NMDS surface for each of the six experimental treatments. In addition, we evaluated dependence of communitylevel forecasting results on experimental conditions (source inocula and media), αdiversity (Shannon’s H′) of ASVs, and mean βdiversity against other replicates in a multivariate ANOVA model of predicated vs. observed community dissimilarity.
Anticipating abrupt community shifts
We then examined whether indices derived from the energy landscape analysis and/or empirical dynamic modeling could be used to anticipate drastic changes in community structure.
In the framework of energy landscape analysis, we calculated two types of indices based on the estimated landscapes of microbiome dynamics (Fig. 3a). One is deviation of current communitystate energy from the possible lowest energy within the target basins (hereafter, energy gap; Fig. 3a): this index represents how current community states are inflated from local optima (i.e., “bottom” of basins). The other is “stablestate entropy [24]”, which is calculated based on the randomwalkbased simulation from current community states to bottoms of any energy landscape basins (i.e., alternative stable states). A starting community state is inferred to have high entropy if reached stable states are variable among randomwalk iterations: the stablestate entropy is defined as the Shannon’s entropy of the final destinations of the random walk [24]. Therefore, communities approaching abrupt structural changes are expected to have high stablestate entropy because they are inferred to cross over “ridges” on energy landscapes [24]. For an analysis of a target replicate community, energy landscapes were reconstructed based on the data of the remaining seven replicate communities.
In the framework of empirical dynamic modeling (nonlinear mechanics), we calculated “local Lyapunov stability [39]” (local dynamic stability) and “local structural stability [40]” based on Jacobian matrices representing movements around reconstructed attractors [27]. Specifically, based on convergent crossmapping [22, 32] and multivariate extension of Smap [65], local Lyapunov stability and structural stability were estimated, respectively, as the absolute value of the dominant eigenvalue and trace (sum of diagonal elements) of the Jacobian matrices representing the timeseries processes [39]. For a target replicate community, the remaining seven replicate communities were used for inferring Jacobian matrices. Note that a high score of local Lyapunov/structural stability represents a potentially unstable community state. In particular, local Lyapunov scores reflect whether trajectories at any particular time are converging (local Lyapunov score < 1) or diverging (1 < local Lyapunov score) [39].
For each of the above indices, linear regression of abruptness scores of communitycompositional changes was performed for each replicate sample in each experimental treatment (7dayahead forecasting). The time points (samples) excluded in the dataquality filtering process (see “Bioinformatics” section) were excluded from this evaluation of signal indices.
We also examined the diagnostic performance of the signal indices based on the receiver operating characteristic (ROC) analysis. In 7dayahead forecasting, detection rates (sensitivity) and false detection rates (1–specificity) of large communitycompositional changes (abruptness > 0.5) were plotted on a twodimensional surface for each experimental treatment, yielding area under the ROC curve (AUC) representing diagnostic performance [41]. The optimal threshold value of each signal index for anticipating abrupt communitycompositional changes (abruptness > 0.5) was then calculated with the Youden index [41] for each experimental treatment. In addition, for stablestate entropy and local Lyapunov stability, we calculated optimal threshold values after assembling all the data of MediumA and MediumB treatments.
Availability of data and materials
The 16S rRNA sequencing data are available from the DNA Data Bank of Japan (DDBJ) with the accession number DRA013352, DRA013353, DRA013354, DRA013355, DRA013356, DRA013368, and DRA013379. The microbial community data are deposited at the figshare repository (10.6084/m9.figshare.20653440) [to be released after the acceptance of the paper]. All the codes used to analyze the data are available at the figshare repository (10.6084/m9.figshare.20653440) [to be released after the acceptance of the paper].
Abbreviations
 ASV:

Amplicon sequence variants
 AUC:

Area under the curve
 DDBJ:

DNA Data Bank of Japan
 FDR:

False discovery rate
 NMDS:

Nonmetric multidimensional scaling
 ROC:

Receiver operating characteristic
 rRNA:

Ribosomal ribonucleic acid
References
Costello EK, Stagaman K, Dethlefsen L, Bohannan BJM, Relman DA. The application of ecological theory toward an understanding of the human microbiome. Science. 1979;2012(336):1255–62.
Cho I, Blaser MJ. The human microbiome: At the interface of health and disease. Nat Rev Genet. 2012;13:260–70.
Huttenhower C, Gevers D, Knight R, Abubucker S, Badger JH, Chinwalla AT, et al. Structure, function and diversity of the healthy human microbiome. Nature. 2012;486:207–14.
Wu GD, Chen J, Hoffmann C, Bittinger K, Chen YY, Keilbaugh SA, et al. Linking longterm dietary patterns with gut microbial enterotypes. Science. 1979;2011(334):105–8.
Toju H, Peay KGKG, Yamamichi M, Narisawa K, Hiruma K, Naito K, et al. Core microbiomes for sustainable agroecosystems. Nat Plants. 2018;4:247–57.
Busby PE, Soman C, Wagner MR, Friesen ML, Kremer J, Bennett A, et al. Research priorities for harnessing plant microbiomes in sustainable agriculture. PLoS Biol. 2017;15:e2001793.
Kazamia E, Aldridge DC, Smith AG. Synthetic ecology  A way forward for sustainable algal biofuel production? J Biotechnol. 2012;162:163–9.
Vatanen T, Franzosa EA, Schwager R, Tripathi S, Arthur TD, Vehik K, et al. The human gut microbiome in earlyonset type 1 diabetes from the TEDDY study. Nature. 2018;562:589–94.
Zhao L, Zhang F, Ding X, Wu G, Lam YY, Wang X, et al. Gut bacteria selectively promoted by dietary fibers alleviate type 2 diabetes. Science. 1979;2018(359):1151–6.
Kim YG, Sakamoto K, Seo SU, Pickard JM, Gillilland MG, Pudlo NA, et al. Neonatal acquisition of Clostridia species protects against colonization by bacterial pathogens. Science. 1979;2017(356):315–9.
Chu C, Murdock MH, Jing D, Won TH, Chung H, Kressel AM, et al. The microbiota regulate neuronal function and fear extinction learning. Nature. 2019;574:543–8.
Sato Y, Hori T, Koike H, Navarro RR, Ogata A, Habe H. Transcriptome analysis of activated sludge microbiomes reveals an unexpected role of minority nitrifiers in carbon metabolism. Commun Biol. 2019;2:179.
Carding S, Verbeke K, Vipond DT, Corfe BM, Owen LJ. Dysbiosis of the gut microbiota in disease. Microb Ecol Health Dis. 2015;26:26191.
Ravel J, Brotman RM, Gajer P, Ma B, Nandy M, Fadrosh DW, et al. Daily temporal dynamics of vaginal microbiota before, during and after episodes of bacterial vaginosis. Microbiome. 2013;1:29.
Hastings A, Abbott KC, Cuddington K, Francis T, Gellner G, Lai YC, et al. Transient phenomena in ecology. Science. 1979;2018:361.
Fukami T. Historical contingency in community assembly: integrating niches, species pools, and priority effects. Annu Rev Ecol Evol Syst. 2015;46:1–23.
Beisner BE, Haydon DT, Cuddington K. Alternative stable states in ecology. Front Ecol Environ. 2003;1:376–82.
Hsieh CH, Glaser SM, Lucas AJ, Sugihara G. Distinguishing random environmental fluctuations from ecological catastrophes for the North Pacific Ocean. Nature. 2005;435:336–40.
Sugihara G. Nonlinear forecasting for the classification of natural time series. Philos Transact R Soc London Ser A Phys Eng Sci. 1994;348:477–95.
Sugihara G, May RM. Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature. 1990;344:734–41.
Benincá E, Huisman J, Heerkloss R, Jöhnk KD, Branco P, van Nes EH, et al. Chaos in a longterm experiment with a plankton community. Nature. 2008;451:822–5.
Sugihara G, May R, Ye H, Hsieh C, Deyle E, Fogarty M, et al. Detecting causality in complex ecosystems. Science. 2012;338:496–500.
Strogatz SH. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. 2nd edition. New York: CRC Press; 2015.
Suzuki K, Nakaoka S, Fukuda S, Masuya H. Energy landscape analysis elucidates the multistability of ecological communities across environmental gradients. Ecol Monogr. 2021;91:1–21.
Watanabe T, Masuda N, Megumi F, Kanai R, Rees G. Energy landscape and dynamics of brain activity during human bistable perception. Nat Commun. 2014;5:4765.
Becker OM, Karplus M. The topology of multidimensional potential energy surfaces: Theory and application to peptide structure and kinetics. J Chem Phys. 1997;106:1495.
Deyle ER, May RM, Munch SB, Sugihara G. Tracking and forecasting ecosystem interactions in real time. Proc R Soc B Biol Sci. 2016;283:20152258.
Chang CW, Ushio M, Hsieh C, hao. Empirical dynamic modeling for beginners. Ecol Res. 2017;32:785–96.
Munch SB, Brias A, Sugihara G, Rogers TL. Frequently asked questions about nonlinear dynamics and empirical dynamic modelling. ICES J Mar Sci. 2019. https://doi.org/10.1093/icesjms/fsz209.
Lahti L, Salojärvi J, Salonen A, Scheffer M, de Vos WM. Tipping elements in the human intestinal ecosystem. Nat Commun. 2014;5:1–10.
David LA, Maurice CF, Carmody RN, Gootenberg DB, Button JE, Wolfe BE, et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature. 2014;505:559–63.
Ushio M. Interaction capacity as a potential driver of community diversity. Proc R Soc B Biol Sci. 2022;289:20212690.
Goldford JE, Lu N, Bajić D, Estrela S, Tikhonov M, SanchezGorostiaga A, et al. Emergent simplicity in microbial community assembly. Science. 1979;2018(361):469–74.
Scheffer M, Carpenter SR. Catastrophic regime shifts in ecosystems: Linking theory to observation. Trends Ecol Evol. 2003;18:648–56.
Scheffer M, Carpenter S, Foley JA, Folke C, Walker B. Catastrophic shifts in ecosystems. Nature. 2001;413:591–6.
May RM. Thresholds and breakpoints in ecosystems with a multiplicity of stable states. Nature. 1977;269:4717.
Ovaskainen O, Tikhonov G, Dunson D, Grøtan V, Engen S, Sæther BE, et al. How are species interactions structured in speciesrich communities? A new method for analysing timeseries data. Proc R Soc B Biol Sci. 2017;284:20170768.
Takens F. Detecting strange attractors in turbulence. In: Rand DA, Young LS, editors. Dynamical Systems and Turbulence: Springer; 1981. p. 366–81.
Ushio M, Hsieh CH, Masuda R, Deyle ER, Ye H, Chang CW, et al. Fluctuating interaction network and timevarying stability of a natural fish community. Nature. 2018;554:360–3.
Cenci S, Saavedra S. Nonparametric estimation of the structural stability of nonequilibrium community dynamics. Nat Ecol Evol. 2019;3:912–8.
Akobeng AK. Understanding diagnostic tests 3: Receiver operating characteristic curves. Acta Paediatrica Int J Paediatrics. 2007;96:644–7.
Chang C, Miki T, Ushio M, Ke P, Lu H, Shiah F, et al. Reconstructing large interaction networks from empirical time series data. Ecol Lett. 2021;24:2763–74.
Amor DR, Ratzke C, Gore J. Transient invaders can induce shifts between alternative stable states of microbial communities. Sci Adv. 2020;6:eaay8676.
Jones CG, Lawton JH, Shachak M. Organisms as ecosystem engineers. Oikos. 1994;69:373–86.
OdlingSmee FJ, Laland KN, Feldman MW. Niche construction: The neglected process in evolution; 2013.
Arumugam M, Raes J, Pelletier E, le Paslier D, Yamada T, Mende DR, et al. Enterotypes of the human gut microbiome. Nature. 2011;473:174–80.
Mee MT, Collins JJ, Church GM, Wang HH. Syntrophic exchange in synthetic microbial communities. Proc Natl Acad Sci U S A. 2014;111:E2149–56.
Vrancken G, Gregory AC, Huys GRB, Faust K, Raes J. Synthetic ecology of the human gut microbiota. Nat Rev Microbiol. 2019;17:754–63.
Ushio M, Murakami H, Masuda R, Sado T, Miya M, Sakurai S, et al. Quantitative monitoring of multispecies fish environmental DNA using highthroughput sequencing. Metabarcoding Metagenom. 2018;2:1–15.
Caporaso JG, Lauber CL, Walters WA, BergLyons D, Lozupone CA, Turnbaugh PJ, et al. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci U S A. 2011;108(SUPPL. 1):4516–22.
Apprill A, Mcnally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat Microb Ecol. 2015;75:129–37.
Klappenbach JA, Saxman PR, Cole JR, Schmidt TM. Rrndb: The ribosomal RNA operon copy number database. Nucleic Acids Res. 2001;29:181–4.
Lundberg DS, Yourstone S, Mieczkowski P, Jones CD, Dangl JL. Practical innovations for highthroughput amplicon sequencing. Nat Methods. 2013;10:999–1002.
Stevens JL, Jackson RL, Olson JB. Slowing PCR ramp speed reduces chimera formation from environmental samples. J Microbiol Methods. 2013;93:203–5.
Hamady M, Walker JJ, Harris JK, Gold NJ, Knight R. Errorcorrecting barcoded primers for pyrosequencing hundreds of samples in multiplex. Nat Methods. 2008;5:235–7.
Tanabe A. Claident v0.2.2018.05.29, a software distributed by author at http://www.fifthdimension.jp/. 2018.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: Highresolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
Wang Q, Garrity GM, Tiedje JM, Cole JR. Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol. 2007;73:5261–7.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: Improved data processing and webbased tools. Nucleic Acids Res. 2013;41:D590–6.
Oksanen J. The vegan package available at https://cran.rproject.org/web/packages/vegan/index.html. 2007.
Harris DJ. Inferring species interactions from cooccurrence data with Markov networks. Ecology. 2016;97:3308–14.
Watanabe T, Hirose S, Wada H, Imai Y, Machida T, Shirouzu I, et al. Energy landscapes of restingstate brain networks. Front Neuroinform. 2014;8:12.
Wood S. mgcv: Mixed GAM computation vehicle with automatic smoothness estimation available at https://cran.rproject.org/web/packages/mgcv/index.html. 2022.
Navarrete R. Embeddings and prediction of dynamical time series. Doctor thesis: The University of Michigan; 2018.
Cenci S, Sugihara G, Saavedra S. Regularized Smap for inference and forecasting with noisy ecological time series. Methods Ecol Evol. 2019;10:650–60.
Acknowledgements
We thank Sayaka Suzuki and Keisuke Koba for support in the experiment and Tadashi Fukami and anonymous reviewers for insightful comments on the manuscript. Computation time was provided by the SuperComputer System, Institute for Chemical Research, Kyoto University.
Funding
This work was financially supported by JST PRESTO (JPMJPR16Q6), Human Frontier Science Program (RGP0029/2019), JSPS GrantinAid for Scientific Research (20K20586), NEDO Moonshot Research and Development Program (JPNP18016), and JST FOREST (JPMJFR2048) to H.T., JSPS GrantinAid for Scientific Research (20K06820 and 20H03010) to K.S., and JSPS Fellowship to H.F. and A.C.
Author information
Authors and Affiliations
Contributions
H.T. designed the work with H.F.; H.F. performed experiments; H.F. analyzed the data with M.U., K.S., M.S.A., M.Y., and H.T.; H.F. and H.T. wrote the paper with all the authors. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1: Figure S1.
Experimental setting and microbiome data formats.
Additional file 2: Figure S2.
Dynamics of calibrated abundance without interpolation.
Additional file 3: Figure S3.
Dynamics of calibrated abundance.
Additional file 4: Figure S4.
Dynamics of relative abundance.
Additional file 5: Figure S5.
Distribution of stable states on the energy landscapes.
Additional file 6: Figure S6.
Transitions between alternative stable states on the energy landscapes.
Additional file 7: Figure S7.
Histogram of optimal embedding dimension (E) of ASVs.
Additional file 8: Figure S8.
Examples of populationlevel forecasting results.
Additional file 9: Figure S9.
Dependence of populationlevel forecasting results on reference database size.
Additional file 10: Figure S10.
Comparison of predicted and observed community structure (sevendayahead forecasting).
Additional file 11: Figure S11.
Distribution of prediction error in the communitylevel forecasting.
Additional file 12: Figure S12.
Comparison of nonlinear and linear forecasting approaches.
Additional file 13: Figure S13.
Candidates of signal indices for anticipating abrupt community changes.
Additional file 14: Figure S14.
Dependence of relationship between signal index values and observed communitycompositional changes on time steps in forecasting.
Additional file 15: Figure S15.
ROC analysis of diagnostic performance.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Fujita, H., Ushio, M., Suzuki, K. et al. Alternative stable states, nonlinear behavior, and predictability of microbiome dynamics. Microbiome 11, 63 (2023). https://doi.org/10.1186/s40168023014745
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s40168023014745
Keywords
 Alternative stable states
 Biodiversity
 Biological communities
 Chaos
 Community collapse
 Community stability
 Dysbiosis
 Empirical dynamic modeling
 Microbiome dynamics
 Nonlinear dynamics