Microbial communities perform essential ecosystem services and carry out important functions in their hosts. For instance, the healthy human gut microbiome protects its host from pathogens, expands the host’s digestive capacities and contributes to human immune system maturation [1]. Thanks to recent advances in sequencing technology, we now have access to data sets that capture the dynamics of entire microbial communities at a high phylogenetic resolution over long periods of time. Such densely sampled long-term time series data have been collected for instance for skin and gut [2, 3], but also for lakes and oceans [4, 5]. These studies have illustrated that proportions of microbial taxa do not always remain constant, but can vary over time, sometimes considerably so, for instance in cases of seasonality [5] and succession [6].

Fluctuations in microbial community composition have been linked to a variety of inter-dependent factors, including ecological interactions between community members, environmental conditions, immigration from adjacent ecosystems, the history of the community, and the evolution of community members [7, 8]. The importance of these factors varies across ecosystems; hence, a better characterization of their contribution will enable selecting suitable community models and improve our understanding of microbial community functions.

A popular strategy to learn more about microbial community structure and dynamics is to fit community models to sequencing data of microbial communities and to analyze the parameterized models. Two frequently selected community models are Hubbell’s neutral model [9, 10] and the generalized Lotka-Volterra (gLV) model [11, 12], where the former assumes that species are ecologically equivalent and community dynamics is governed by local extinction and immigration from a metacommunity, whereas the latter describes the change of species abundances as a function of growth rates and species interactions. The parameters (and therefore the species interactions) of the deterministic gLV model and its stochastic, time discrete version, the Ricker model, have been inferred directly from time series data by several authors (e.g., [13–18]). The continuous form of the Hubbell model developed by Sloan [19] served as a null model against which over- or under-representation of microbial taxa was tested [20, 21]. The self-organized instability (SOI) model [22] proposed by Solé and colleagues combines aspects of the gLV and neutral model, namely interactions with stochastic immigration and extinction. Alternative ways to integrate interactions and immigration have also been suggested [23–26].

Although these models emphasize different aspects of community dynamics, they can be seen as realizations of an encompassing community model framework spanning a gradient from stochastic to deterministic and from unstructured to temporally structured in time (Fig. 1). Temporal structure is defined here in the sense of a Markovian structure, which means that the state of the system at a certain time point depends on its state at preceding time points. Structured models, which include all models considered here except the DM, encode rules that describe how the community composition changes from one time step to the next. In contrast, unstructured models, which are stochastic, successively sample independent microbial community composition from a probability distribution, such as the Dirichlet-multinomial distribution [27, 28].

Temporally structured stochastic models, such as the neutral and the SOI model, describe microbial community dynamics at the level of individuals. With increasing number of individuals, their behavior approaches that of deterministic models describing community dynamics at the level of populations, such as the gLV and its discrete variant, the Ricker model [14]. These deterministic models in turn can move towards the stochastic end of the gradient by integrating intrinsic noise or random environmental fluctuations of increasing strength.

In this work, we aim to provide guidelines that help distinguish the different generating processes directly from the time series, thereby complementing and guiding standard model fitting procedures. More precisely, our goal is to determine whether the inference of a microbial network from a time series data set is meaningful.

Our proposed classification scheme consists of three steps: (i) test for temporal structure, (ii) test for neutrality, (iii) fit an interaction model (Fig. 2). This classification scheme is complementary to the one proposed by Gibbons and colleagues, which does not test for the presence of temporal structure or neutrality [29]. In the following, we will apply this scheme to simulated and real-world microbial time series data.

In order to demonstrate the efficacy of our classification scheme, we generated time series with the Dirichlet-multinomial (DM) distribution, Hubbell’s neutral model, the SOI model, the Ricker model with varying levels of noise, and the generalized Lotka-Volterra model. In addition, we included time series from two individuals (A and B) from a long-term study of the human gut microbiota [3]. In total, we generated 60 community time series with various parameter settings, each including 100 species followed over 3000 time points. Since such long time series are not yet available in practice, we also repeated every test for the first 100 time points of each time series. An overview of the parameter settings, time series properties, and test results is given in Additional file 1: Table S1. Time series generation is summarized in “Methods”; model details are provided in Additional file 2.

Noise types differentiate between unstructured and structured models

The first step of the classification scheme distinguishes between the presence of temporal structure, i.e., rules that govern the dynamics, and its absence. To test for the presence of temporal structure, we exploit the fact that it will introduce a dependency between time points that is absent otherwise. Different measures exist that can detect dependency on previous community states, such as autocorrelation, the Hurst exponent, and noise types. We focus on the latter here, but also tested the two other measures with similar results (see Additional file 3: Figure S1). The noise type of a species is obtained by decomposing its abundance fluctuations (or relative abundance fluctuations in most cases) into spectral densities at specific frequencies using a Fourier transform. When spectral densities tend to be high at low frequencies (i.e., long intervals) and vice versa, they signal a strong dependency on previous time points, which is visible as a negative slope of densities versus frequencies in the periodogram. Depending on the value of this slope in log-log scale, one can distinguish black (below − 2), brown (around − 2), pink (around − 1), and white noise (no negative slope). A shift in the noise color from black to pink indicates a reduced dependency on previous community states. In brief, the dependency on previous time points is strongest for black noise and weakest for pink noise, while it is absent for white noise.

We computed the noise type for each species and reported the percentage of black, brown, pink, and white species in each community time series (Fig. 3). While structured community models (gLV, Ricker, Hubbell, and SOI) generate time series with no or only a small percentage of white noise species, the DM distribution gives rise to mostly white noise species, as expected for unstructured data. Deterministic models such as gLV and Ricker with low intrinsic noise are dominated by black noise species, whereas Hubbell time series are mostly composed of brown noise species and SOI time series of pink noise species. However, we found that longer intervals and shorter time series increase the number of white noise species in the stool data and SOI time series, that higher connectance increases the proportion of black taxa at the cost of pink ones in Ricker but not SOI time series (Additional file 4: Figure S2 and Additional file 5: Figure S3), that high mortality rates favors pink instead of brown noise in the Hubbell time series (Additional file 6: Figure S4), and that increasing intrinsic noise increases the percentage of pink and brown noise in Ricker time series. These delicate shifts in noise-type profiles illustrate the effect of confounding factors such as sampling interval and noise. While confounders complicate recognizing an underlying model on the basis of its noise-type profile alone, they do not affect the distinction between temporally structured and unstructured community models.

This noise-type-based test for structure is also robust to compositionality and Poisson and multinomial noise (Additional file 7: Figure S5 and Additional file 8: Figure S6) and increases in accuracy with increasing number of time points (Additional file 9: Figure S7 and Additional file 10: Figure S8). While it is more effective for the first 100 time points in the transient part of the time series, it also distinguishes temporally structured from unstructured data in the last 100 time points (Additional file 7: Figure S5).

Neutrality test distinguishes neutral from non-neutral dynamics

To distinguish neutral from non-neutral structured models, we applied a test previously developed to recognize neutrality in time series data [30]. Neutrality, in this test, is defined as a general feature in which demographic rates and performance are independent of species identities. The idea of the neutrality test is to test for group invariance, which means that time series properties are not affected by summing species into groups that represent higher-level taxa such as genera or families. In particular, the group invariance of inter-species quadratic covariation is tested; thus, not only group invariance of covariance is tested, but also its correspondence to a type of covariance common to all neutral processes. When validating this test on our community time series, the DM and Hubbell time series are classified as neutral (p values > 0.05) as expected, whereas the other simulated time series are correctly classified as non-neutral (Fig. 4a, b). The time interval between samples is an important variable when defining neutrality, since for stool data, time series sampled at larger intervals are classified as neutral. Thus, neutrality is defined relative to the time-scale of an investigation. Furthermore, the addition of Poisson noise introduces false positives (time series falsely classified as non-neutral; Additional file 11: Figure S9).

Microbial network inference through deterministic model fitting

The first two tests classify the stool data as resulting from a temporally structured, non-neutral community. Thus, gut microbial community dynamics is to a large part shaped by interactions between community members, which is in agreement with our knowledge of the numerous cross-feeding and competition interactions taking place between gut organisms [31]. A number of algorithms have been developed to parameterize the gLV or Ricker model directly from time series data [14, 18, 32], thereby inferring ecological interactions. Here, we employ LIMITS [14] to test this parameterization step. We applied it to the top 60 most abundant species in our test time series and measured inference accuracy as the mean correlation between the known and the inferred interaction matrix. The LIMITS algorithm was found to perform well on time series data with known interaction matrices (Fig. 4c). Interestingly, the accuracy of LIMITS is also relatively high for SOI time series, although LIMITS assumes a different community model. These results were reproduced with lower accuracy for short time series and selected time series with Poisson noise (Additional file 11: Figure S9 and Additional file 12: Figure S10).

The goodness of fit to the Ricker model was computed as the mean correlation between the community time series predicted with the parameterized Ricker in a step-wise manner and the test time series (Fig. 4d). Not surprisingly, the goodness of fit is inversely correlated to the strength of the intrinsic noise in Ricker (Spearman’s rho − 0.76, p value < 0.0001). We also applied LIMITS to neutral time series. Although the Hubbell model does not assume direct ecological interactions between species, it does introduce competition between all species. Notably, the goodness of fit of Hubbell time series to the Ricker model was high, but the absolute interaction strengths were significantly smaller than those inferred for all other time series, even when randomly sub-sampling from the latter to account for different sample numbers (Wilcoxon rank sum test p value < 0.0001).

After having tested its performance, we applied LIMITS to explore the interaction networks of the two stool time series. While gLV/Ricker models have been parameterized on stool time series previously [14, 17], the LIMITS algorithm has to our knowledge not yet been applied to the stool time series collected by David and colleagues [3] (see Additional file 13: Figure S11 for a time series plot). We noticed that repeated rarefactions alter the noise-type classification of a few taxa (Additional file 14: Figure S12) and therefore inferred interaction networks from the 30 top-abundant OTUs consistently classified as pink, brown, or black across several rarefactions. The resulting interaction networks have a large percentage of negative links (70% for individual A, 63% for individual B), in agreement with the theoretical expectation that a high percentage of negative interactions is needed to stabilize the community [33]. The same Faecalibacterium OTU (OTU_165924) forms negative hubs in both networks (with nine and five negative links, respectively; Fig. 5). In both networks, Firmicutes and Bacteroidetes OTUs form significantly more links within their phylum than expected at random, with the majority of within-phylum links being negative, while inter-phylum link numbers are not higher than expected at random. Thus, for top-abundant representatives of Bacteroidetes and Firmicutes, intra-phylum competition may be more intense than competition with members of other phyla. Our finding that Faecalibacterium forms hub nodes is in agreement with the results of an alternative network inference method (Granger causality) applied to the same data [29]. However, not all links represent ecological interactions, since networks inferred from real-world data likely contain more false positive links due to underlying environmental variables.

Our results show that several criteria can be combined to discriminate processes underlying the dynamics of microbial communities from time series.

The tendency of ecological time series to display pink, brown, or black noise was noticed previously [34] and linked to the environment [35]. However, since none of our tested models takes environmental factors into account, our results indicate that such noise types can result from the intrinsic community dynamics alone. Self-organized critical systems, to which the SOI model is closely related, in particular are known for their pink (1/f) noise [36], whereas the neutral model was stated previously to generate brown noise [37]. The noise types measure dependency between time points, which can be referred to as a memory effect. The strongest memory effect is found in the gLV and noise-free Ricker time series with their black noise. In general, the gLV time series quickly reach a steady state and thus represent the extreme case of perfect memory. However, for two of the noise-free Ricker time series, all taxa are classified as black while displaying periodicities. Adding a noise term in the Ricker model weakens the memory effect, thereby shifting the noise type from black via brown to pink.

It is current practice in microbial network inference to filter out rare taxa, where rareness of a taxon is defined arbitrarily by applying prevalence or abundance thresholds. Filtering taxa by noise type may constitute an interesting alternative. Taxon abundances that are not dependent on previous abundances belong to taxa that either do not interact or are too rare to reveal interactions and will thus introduce false interactions in the network model. While our work is a first step in this direction, a number of challenges have to be overcome: first, differentiating white from non-white taxa involves a threshold choice and is affected by confounding factors such as noise. For instance, we found that while overall noise-type percentages were conserved, a few (abundant) taxa in the stool data changed their noise type upon re-rarefaction. Second, it is conceivable that although white taxa are not affected by other taxa, they may influence other taxa, for instance through cross-feeding. Extensions of the gLV have been proposed previously to deal with external factors such as antibiotic treatment [14, 17] and environmental factors [4]. White taxa suspected to affect other taxa could be treated as external factors in the gLV parameterization. Tools such as MDSINE allow such a treatment directly, whereas we modified the LIMITS algorithm to support external factors. Cross-sectional microbial network inference on the time series directly or regression on the residuals after gLV model fitting may be able to pinpoint such white taxa.

In this context, we would like to emphasize that the absence of temporal structure in time series only has implications for analyses that assume its presence, such as the inference of interaction networks. Analyses that do not rely on the presence of temporal structure, such as ordination or the comparison of microbial composition between conditions, are not concerned.

In our evaluation of LIMITS, we found high correlation between inferred and known (simulated) interaction matrices. However, in real-world applications, one does not know the true interaction matrix for assessment. Instead, the observed time series is usually compared to the predicted time series [17]. Based on this criterion, LIMITS performed well on Hubbell time series, although their generating model assumes the absence of specific ecological interactions (the replacement step introduces a generic competition between all species). In general, the goodness of fit criterion is biased by the amount of memory present in the time series, since it is easier to predict a time series that is highly autocorrelated than one less so. This strongly highlights the need to test for non-neutrality before network inference. The fact that an agreement between observed and predicted time series can be misleading was also discussed in [38].

Interestingly, in our tests, LIMITS reached reasonable accuracies (quantified as the correlation between known and inferred interaction matrix) not only for Ricker but also for gLV and SOI time series. Thus, the inference is robust with respect to discrete or continuous implementation of a model (Ricker versus generalized Lotka-Volterra) and even to the resolution (populations in the Ricker and gLV model, individuals in the SOI model). In addition, LIMITS performance was not much affected by Poisson noise. However, as also pointed out by Cao et al. [39], network inference accuracy drops with an increasing sampling interval.

The first step in our proposed classification scheme, i.e., the computation of noise types, is robust to compositionality, the presence or absence of transient dynamics, the presence of Poisson and multinomial noise and to the tested sampling intervals and works for relatively short time series, yet it is currently unknown how sensitive it is to the impact of the environment. While the environment may influence community dynamics [40], we may reasonably assume that it will only in extreme cases introduce temporal structure when it is absent otherwise. Moreover, as we remove the linear trend from the time series prior to power spectrum estimation, our approach should be robust to slowly varying environmental conditions. However, this point requires further investigation. It is of note that all the models with temporal structure presented in this work can be modified to account for the influence of external drivers. Although a number of confounding factors (strength of intrinsic noise, external noise, sampling interval) can alter the noise type and hence prevent the identification of the underlying model from the noise-type profile, these processes do not introduce non-white taxa when these are absent. Thus, the presence of non-white taxa is a robust indicator of temporal structure.

The second step, which tests for neutrality, is more sensitive to sampling interval and external noise than the first. While a high percentage of brown noise taxa gives an additional and independent hint for neutral dynamics, this indicator can also be biased by the effects of intrinsic noise, a high death rate and the sampling interval. This highlights the importance of sufficiently dense sampling in longitudinal studies.

According to our tests, the stool data are temporally structured and non-neutral, the latter in agreement with previous results [19, 41]. The noise-type profile of the stool data is closest to the noisy Ricker and SOI time series, indicating the presence of stochasticity. This overall deterministic behavior with a stochastic component agrees well with our expectation for a microbial community whose members are known to engage in various cross-feeding and competitive interactions while being exposed to daily perturbations through the host (e.g., diet).

In conclusion, we have demonstrated that knowledge about the fundamental properties of microbial community dynamics is required for unbiased community study and model inference. We have proposed a classification scheme for microbial sequencing time series that first differentiates between unstructured and structured models, then tests for neutrality and finally determines the goodness of fit to a deterministic model. To make these tests accessible, we have implemented a new R package called seqtime. While several hurdles still have to be overcome, this classification scheme is a first step towards model identification from time series data.