Health and Disease Imprinted in the Time Variability of the Human Microbiome.

The animal microbiota (including the human microbiota) plays an important role in keeping the physiological status of the host healthy. Research seeks greater insight into whether changes in the composition and function of the microbiota are associated with disease. We analyzed published 16S rRNA and shotgun metagenomic sequencing (SMS) data pertaining to the gut microbiotas of 99 subjects monitored over time. Temporal fluctuations in the microbial composition revealed significant differences due to factors such as dietary changes, antibiotic intake, age, and disease. This article shows that a fluctuation scaling law can describe the temporal changes in the gut microbiota. This law estimates the temporal variability of the microbial population and quantitatively characterizes the path toward disease via a noise-induced phase transition. Estimation of the systemic parameters may be of clinical utility in follow-up studies and have more general applications in fields where it is important to know whether a given community is stable or not. IMPORTANCE The human microbiota correlates closely with the health status of its host. This article analyzes the microbial composition of several subjects under different conditions over time spans that ranged from days to months. Using the Langevin equation as the basis of our mathematical framework to evaluate microbial temporal stability, we proved that stable microbiotas can be distinguished from unstable microbiotas. This initial step will help us to determine how temporal microbiota stability is related to a subject's health status and to develop a more comprehensive framework that will provide greater insight into this complex system.

INTRODUCTION

The quest to understand the factors that influence human health and cause disease has always been one of the major driving forces of biological research. With growing evidence of the new “holobiont” and “hologenome” concepts (1, 2), research focuses not only on human physiology but also on the associated microbial population, although these concepts are still under debate (3). Research has revealed that the human microbiome is intimately linked to our physiology through the metabolism of bile acids (4), choline (5), and key metabolites such as short-chain fatty acids (6, 7), which are also involved in immune system maturation (8, 9). The human microbiota is plausibly related to diseases such as type 2 diabetes (10), cardiovascular disease (11), irritable bowel syndrome (IBS) (12), and Crohn’s disease (13) and some afflictions like obesity (14, 15) and malnutrition (16), as well as many other diseases (17). Current studies have revealed that gut microbes also influence brain function and behavior and are related to neurological disorders like Alzheimer’s disease through the gut-brain axis (18, 19). Recently, chronic fatigue syndrome (CFS), a subtle but devastating condition often cited as a psychosomatic disease, has been associated with a reduced diversity and altered composition of the gut microbiome (20).

This research area has progressed greatly thanks to high-throughput methods for microbial 16S rRNA gene sequencing and SMS (shotgun metagenomic sequencing), which reveal the composition of archaeal, bacterial, fungal, and viral communities located in and on the human body. Modern high-throughput sequencing and bioinformatics tools provide a powerful means of understanding how the human microbiome contributes to health and its potential as a target for therapeutic interventions (21). Research is under way to establish normal host-gut microbe interactions and understand how microbiota compositional changes can cause certain diseases (22–24).

Biology has recently acquired new technological and conceptual tools to investigate, model, and understand living organisms at a systems level, thanks to progress in quantitative techniques, large-scale measurement methods, and joint experimental and computational approaches. In particular, systems biology strives to reveal the general laws governing the complex behavior of microbial communities (25–27), including a proposal for universal dynamics (28). Microbiota can be approached in the light of ecological theory, which includes general principles like Taylor’s law (29, 30) relating the spatial or temporal variability of the population with its mean. This law, also known as the fluctuation scale law, is ubiquitous in the natural world and can be found in several systems such as random walks (31), stock markets (32, 33), tree (34) and animal (30, 35, 36) populations, gene expression (37), and the human genome (38). Taylor’s law has been applied to microbiotas spatially by Zhang et al. (39), with results showing that this population tends to be an aggregated one rather than having a random distribution. Despite its ubiquity, this law has only been tested in experimental settings (40, 41) and has never been applied in follow-up studies on microbiota, despite major efforts to infer the community structure from a dynamic point of view (42–44).

This paper presents the hallmarks of health status (healthy or diseased) in the macroscopic properties of the microbiota by studying its temporal variability. We analyzed over 40,000 time series of taxa from the gut microbiomes of 99 subjects obtained from publicly available high-throughput sequencing data related to different conditions, i.e., diseases, diets, trips, obesity, antibiotic therapy, and healthiness. On finding that all of the cases followed Taylor’s law, we used this empirical fact to model how the relative abundances of taxa evolved over time by using the Langevin equation, similarly to the approach applied by Blumm et al. (45). We used this mathematical framework to explore the temporal stability of the microbiotas under different conditions to understand how this is related to the health status of the subjects.

RESULTS

Microbiome temporal variability was analyzed to extract the global properties of the system. As fluctuations in total counts are plagued by systematic errors, we worked on the temporal variability in the relative abundance of each taxon. Our first finding was, without exception, that changes in the relative abundances of taxa followed a ubiquitous pattern, known as the fluctuation scaling law (46) or Taylor’s power law (30). In other words, the microbiota of all of the taxa detected followed , a power law dependence between mean relative abundance x and dispersion σ. The law seems to be ubiquitous, even spanning 6 orders of magnitude in the observed relative abundances. As shown in Fig. 1, where V corresponds to the y intercept and β is the slope of the fit, the most abundant species were less volatile in relative terms than the less abundant ones. The fit to the power law was always robust (R2 > 0.88) and did not depend on microbiome condition. The power law (or scaling) index β and the variability V (hereafter Taylor’s parameters) appear to be correlated with community stability. Accordingly, we assume that Taylor’s parameters behave as proxies for stability. On the one hand, β is a scaling index that provides information about the statistical properties of the ecosystem. If it is 1/2, the system behaves like a Poisson distribution. If β is 1, the system behaves like an exponential distribution. Generally speaking, metagenomes undergo time course variations with β between these two universal classes. In our study, the fact that β was less than 1 indicates that the most abundant taxa in the microbial community were less susceptible to perturbations than the other taxa. On the other hand, the variability V is a direct estimator of the amplitude of fluctuations over time. V represents the maximum variability attainable by a hypothetically dominant genus (with relative abundance close to 1). It is an important parameter that characterizes the type of system. If V is small, the ranking is stable. For example, this would be the case for the number of diagnoses of a particular disease recorded in Medicare during a month (47). If V is large, as occurs for metagenomic samples, the ranking might be unstable, as it would be for the number of hourly page views of articles in Wikipedia (45, 46). Interestingly, the Taylor parameters were related to the health status of the host, which is the main finding of this study.

FIG 1 

x-weighted power law fits of the standard deviations (SD) versus the mean values for each bacterial genus monitored over time. The fit is shown for samples from a healthy subject (top) and from a subject with a diagnosis of IBS (bottom) studied in our lab (12). Taylor’s power law seems to be ubiquitous, spanning 6 orders of magnitude. V corresponds to the y intercept, and β corresponds to the slope of the line. The error bars (mean axis) are the SEM.

Taylor’s parameters describing the temporal variability of the gut microbiome in our sampled subjects are shown in Tables S1 to S4 in the supplemental material. Our results are indicative of ubiquitous behavior. First, the variability (corresponding to the maximum amplitude of fluctuations) was large, which suggests resiliency of the microbiota. Second, the scaling index was always smaller than 1, which means that the more abundant taxa were less volatile than the less abundant ones. In addition, Taylor’s parameters for the microbiome of healthy subjects in different studies (12, 48–53) were compatible with estimated errors. This enabled us to define an area in the Taylor parameter space that we called the healthy zone.

Taylor’s parameters of subjects with either animal-based (A) or plant-based (P) diets (52). Previous to the diet, the population sampled was described by V̄ = 0.09 ± 0.05 and β̄ = 0.77 ± 0.04. Download TABLE S1, PDF file, 0.03 MB.

Taylor’s parameters for subjects taking antibiotics (Ab) in the antibiotic study (49), persons with a diagnosis of IBS in the IBS study (12) and for special intervals concerning gut microbiota in the host lifestyle study (HLS) (53). Prior to antibiotic intake, the population sampled in the antibiotic study (49) was described by V̄ = 0.12 ± 0.05 and β̄ = 0.75 ± 0.04. Healthy subjects sampled in the IBS study (12) were characterized by V̄ = 0.135 ± 0.010 and β̄ = 0.692 ± 0.024. The healthy and quotidian periods in the host lifestyle study (53) are characterized by V̄ = 0.25 ± 0.09 and β̄ = 0.777 ± 0.025. Download TABLE S2, PDF file, 0.03 MB.

Taylor’s parameters for the healthy-subject (DH) and kwashiorkor (DK) parts of the discordant-twin studies (51). The population of healthy twins is characterized by V̄ = 0.25 ± 0.10 and β̄ = 0.863 ± 0.028. Download TABLE S3, PDF file, 0.03 MB.

Taylor’s parameters for subjects with various degrees of overweight and obesity (50). Healthy people in this study, who were not obese, are characterized by V̄ = 0.19 ± 0.06 and β̄ = 0.806 ± 0.034. Download TABLE S4, PDF file, 0.03 MB.

To jointly visualize and compare the results of subjects from the above-mentioned studies (12, 48–53), their Taylor parameters were standardized, with standardization meaning that the mean value was subtracted from each parameter and the result was divided by the standard deviation of the group of healthy subjects for every study independently. Because of the different systematics in each study, we defined a healthy region for each of them, standardized to a mean of 0 and a variance of 1, then we computed the mean and variance of the “unhealthy” with this standardization (for details of the procedure, see the standardization subsection in Materials and Methods). Therefore, different studies were isolated so that the subjects in any study did not affect the results for the “unhealthy” subjects in the other studies. We think this statistical approach was safer, as we avoided combining data with very different systematic errors. The healthy zone and the standardized Taylor parameters for the temporal variability of the gut microbiome in subjects whose gut microbiota was compromised (i.e., they had IBS, kwashiorkor, an altered diet, antibiotic intake, or a Salmonella infection or had gone on a trip abroad) are shown in Fig. 2. The variability in children with kwashiorkor was smaller than that of their healthy twins. A meat/fish-based diet significantly increased variability compared to a plant-based diet. All other cases presented increased variability, which was particularly severe and statistically significant at an over 95% confidence level (CL), for grade III obese patients on a diet, subjects taking antibiotics, the subject who had a Salmonella infection, the subject who had traveled abroad, and the patients with a diagnosis of IBS. One global property that emerges from these comprehensive data is that the Taylor parameters characterized the statistical behavior of microbiome changes. Furthermore, we verified that our conclusions were robust to systematic errors resulting from taxonomic assignment (see taxon level selection in Materials and Methods).

FIG 2 

Taylor’s law parameter space. All of the data studied in this work are compiled here. The large light-blue circle corresponds to the 68% CL region of healthy subjects in the Taylor parameter space, while the dashed circumference around it delimits the 98% CL region. Points with errors place the gut microbiome in the Taylor parameter space for each subject whose microbiota was compromised. It should be noted that the parameters were standardized (standard deviation units) to the healthy group for every single study independently for purposes of demonstration and comparison.

Taylor’s power law has been explained in terms of various effects, though without a general consensus. It has its origin in mathematical convergence, which is similar to the central limit theorem, and thus, virtually any statistical model designed to produce a Taylor law converges to a Tweedie distribution (54), providing a mechanistic explanation based on the statistical theory of errors (55–57). To reveal the generic mechanisms that drive different scenarios in the β-V space, we modeled the system by assuming that taxon relative abundance followed a Langevin equation with, on the one hand, a deterministic term that captured the fitness of each taxon and, on the other hand, a randomness term associated with Gaussian random noise (45). Both terms were modeled by power laws, with coefficients that can be interpreted as taxon fitness F and variability V (see model in Materials and Methods). F captures the time scale that the system needs to reach equilibrium (the size of V may or may not allow equilibrium to be reached). F has dimensions of 1/time and roughly corresponds to the half-life of the system when it is decaying to a stable state. In fact, it is exactly the half-life if β is 1 and V is negligible. In this model, when V is sufficiently low, abundances are stable in time. Differences in V can induce a noise-induced phase transition in the relative abundances of taxa. The Fokker-Planck equation governs the temporal changes in the likelihood that a given taxon has an abundance x, given its fitness. The results of solving this equation show that stability is best captured by a phase space determined by the fitness (F) and the amplitude of fluctuations (V) (Fig. 3).

FIG 3 

Microbiota states can be placed in the phase space F-V. The light-blue-shaded region corresponds to the stable phase, while the gray-shaded region is the unstable phase (the phase transition line is calculated for α = β = 0.75). We placed healthy subjects (green) and subjects whose gut microbiota is threatened (antibiotics, IBS) in the phase space fitness-variability. The gut microbiotas of healthy subjects over a long time span show a quasiperiodic variability (the central period is 10 days). We show that taking antibiotics (AB1 and AB2 correspond to the first and second treatments, respectively) induces a phase transition in the gut microbiota that has an impact on future changes. We also show a patient with a diagnosis of IBS transitioning from the unstable to the stable phase.

The model predicted two phases for the gut microbiome, a stable phase with large variability that enabled some changes in the relative abundances of taxa and an unstable phase with even larger variability, above the phase transition, where the order of abundant taxa varies significantly over time. The phase transition is continuous (of second order), as is the crossing of the boundary. The state variable is the composition. Any disturbance modifies the composition of the microbiota, with different compositions encoding different F and V values. We have shown that effective perturbations significantly change V and lead the microbiota to a transition from the ordered phase to the noise-induced one. Our model can be solved analytically, which allows for a simple understanding of the different regimes and, in particular, to calculate the formula of the transition region. The order parameter is the composition x that maximizes the probability distribution; 0 < x < 1 defines the ordered phase, while x > 1 defines the disordered phase. If V is sufficiently smaller than F, likelihood peaks in the physical region (relative compositions larger than 0 and smaller than 1); i.e., there is a best-composition solution of the differential equation, which is the ordered solution. Conversely, if V is sufficiently larger than F, the likelihood peaks outside the physical region; i.e., the best-composition solution of the differential equation is at the boundaries (either 0 or 1) and all physical solutions have comparable likelihoods, i.e., the noise-induced phase. The microbiomes of healthy subjects were found to be in the stable phase, while the microbiomes of several other subjects were in the unstable phase. In particular, subjects taking antibiotics and patient P2 (with a diagnosis of IBS) had the most severe symptoms. In this phase diagram (Fig. 3), each microbiota state is represented by a point at its measured V and inferred F. The model predicted high average fitness for all taxa; i.e., taxa were narrowly distributed in F. The fitness parameter was chosen with different values as a demonstration. Fitness was larger for the healthiest subjects and smaller for the patients with a diagnosis of IBS.

Rank stability of the taxa.

The rank dynamics and stability plots in Fig. 4 and 5 show the variations in rank over time for the most dominant taxa and their calculated rank stability index (RSI; as discussed in Materials and Methods) for the gut microbiome taxa of a healthy subject, namely, subject A in the host lifestyle study (53). Figure 4 covers the period when the subject was traveling abroad, and Fig. 5 covers the subsequent period. The taxa are listed according to their accumulated frequency over the time series, with the y axis being the overall dominant axis for each sample set. Generally speaking, we observed that the most dominant taxa had the highest rank stability.

FIG 4 

Rank variation over time of the 50 most dominant elements (taxa) and their calculated RSI, RV, and DV, as detailed in rank stability and variability in Materials and Methods, for a special period (days 72 to 122, traveling abroad) belonging to subject A in the host lifestyle study (53).

FIG 5 

Rank variation over time of the 50 most dominant elements (taxa) and their calculated RSI, RV, and DV, as detailed in rank stability and variability in Materials and Methods, for an ordinary period (days 123 to 256, after the trip) belonging to subject A in the host lifestyle study (53).

For the trip abroad in Fig. 4, beyond the differences in the dominance of particular taxa, we observed that the most dominant ones were also the most rank stable ones. Moreover, the medium-ranked taxa were quite rank unstable, mostly because of transient (often one or two consecutive samples), albeit dramatic, declines in their relative abundance, which usually occurred more than twice during their time series.

Nevertheless, in the case of the next period (Fig. 5), the one subsequent to the trip, some taxa showed higher stability than other more dominant taxa, forming “rank stability islands” for medium-ranked taxa, and displaying a moderately stable index (RSI roughly over 70%). In particular, this was the case for the genera Actinomyces, Leuconostoc, Lachnobacterium, Eggerthella, Clostridium, and Collinsella. For these genera, both the overall rank and the RSI were clearly lower during the trip (RSI under 70%). Actinomyces and Lachnobacterium are not shown in Fig. 4 because they sank to positions 56 and 77, respectively. In contrast, Leuconostoc was the least sensitive to the lifestyle change. Interestingly Lachnobacterium showed anticorrelation over time compared to the vast majority of the taxa classified in this study.

We also found those “rank stability islands” for medium-ranked taxa in the other periods belonging to subject A in the host lifestyle study (53) (see Fig. S1 and S2 for the corresponding rank plots). See Table S5 for details of the rank and RSI of the above-mentioned taxa over the different periods considered.

Rank variation over time of the 50 most dominant elements (taxa) and their calculated RSI, RV, and DV, as detailed in rank stability and variability in Materials and Methods, for an ordinary period (days 0 to 70, before the trip) belonging to subject A in the host lifestyle study (53). Download FIG S1, EPS file, 1.1 MB.

Rank variation over time of the 50 most dominant elements (taxa) and their calculated RSI, RV, and DV, as detailed in rank stability and variability in Materials and Methods, for an ordinary period (days 257 to 364, further after the trip) belonging to subject A in the host lifestyle study (53). Download FIG S2, EPS file, 1.2 MB.

Rank and RSI, as discussed in Materials and Methods, over different periods for the taxa listed as rank stability islands regarding the gut microbiome of subject A in the host lifestyle study (53). Download TABLE S5, PDF file, 0.04 MB.

Time dependence of model parameters.

Finally, we studied the time dependence of variability V and the power law index β (see model in Materials and Methods) by using a sliding-window approach. The total number of time points was divided into subsets of five points, where the following subset was defined by adding the next sampling time and eliminating the earliest one. Both parameters were calculated for each subset against the average time lapse. Figure 6 shows variability V as a function of time for the two subjects in the study of Caporaso et al. (48) corresponding to the gut microbiota of a male (upper plot) and a female (lower plot). Both samples showed changes in V with quasiperiodic behavior peaking at about 10 days. Variability grew more for the gut microbiota of the male and shared a minimal value of around 0.1 with the gut microbiota of the female.

FIG 6 

V of the two subjects in the study of Caporaso et al. (48) as a function of time. Samples of the gut microbiomes of a male (top) and a female (bottom) are shown.

Figure 7 shows time course changes in V for patient P2 in the IBS study (13) (upper plot) and patient D in the antibiotic study (49) (lower plot). The variability of the gut microbiota of P2 decreased from over 0.3 to below 0.2, showing a slow tendency to increase the order of the system. Antibiotic intake led to a quick increase in variability that lasted a few days until V recovered a low level. The second antibiotic treatment showed some memory traits (smaller increase in variability) with slower recovery.

FIG 7 

V as a function of time for patient P2 in the IBS study (12) (top) and patient D in the antibiotic study (47) (bottom). The blue vertical lines in the lower plot show the periods of antibiotic treatment.

DISCUSSION

One of the highlights of this study is that it shows, independently of its condition, that the microbiota follows Taylor’s law. We have seen that, in each case, the value of the scaling index is always less than unity (using the standard deviation as the dispersion measurement), which provides us with information about the community structure. This means that, in relative terms, the most abundant elements in the population are less volatile to perturbations than the less abundant ones. The explanation for this universal pattern is not clear, although some hypotheses have been tested in other studies, such as the presence of negative interactions in the population (58), and a demonstration that this may depend on reproductive correlation (59). Nevertheless, none of these explanations is sufficient when it comes to the microbiota, as the term reproduction is diffuse and the interactions between its components are not only based on competition (60–62). Moreover, even such a negative interaction may not effectively yield values of less than unity when referring to a bacterial species (41). Nonetheless, the values obtained in all cases were very similar to each other, which may suggest that the community structure is preserved throughout the different scenarios studied here.

The second parameter provides information about noise and can be directly linked to the variability or fluctuation amplitude of the population over time. It is a direct estimator of the stability of the system under study. As we have shown above, the healthy subset in each study has less variability than the nonhealthy subset when adult subjects are considered. Interestingly, the variability parameter was higher in the healthy subset in the study of discordant twins suffering from kwashiorkor disease (51). In this respect, research has shown that the infant microbiota needs to develop toward a definite, adult state (63). This implies that temporal variability is greater in children than in healthy adults, who should be temporally stable. Thus, our results could indicate that this variability is necessary to reach that adult state. Furthermore, as we wanted to see how this variability shifted over time, we calculated the changes in this parameter for the samples that had enough sampling coverage. As shown in Fig. 6, the variability of the microbiota fluctuated over time. Interestingly Fig. 7 shows how this parameter reflected the two antibiotic intakes by one of the patients in the study by Dethlefsen and Relman (49), particularly the apparent resilience of the microbiota due to the reduced increase in variability during the second antibiotic intake.

The primary hypothesis of this work is that, in adults, having a healthy microbiota means that the microbial population is stable over time. This stability means that the microbiota does not shift and become susceptible to external or internal perturbations causing dysbiosis. To use the valuable information provided by the empirical law of Taylor’s work, here we have proposed the use of Langevin’s equation to model how stability ranking changed over time. While the system noise component can be directly measured as its variability, the other main term needs to be inferred from the model. This term, which we have named “fitness,” enables the system to remain stable when confronted with potential perturbations. In ecological terms, this could represent the nature of the interactions present among bacteria, between bacteria and other minority populations such as fungi or archaea, and between bacteria and the viral component of the microbiota and interactions between the host and the whole microbiota. As this is a first step in modeling the temporal stability of the microbiota and given its complex nature, we calculated fitness by using the fluctuation dissipation theorem as a first approximation (64). Thus, further work is required to model the fitness of the microbiota to provide a more accurate model with greater predictive power.

By solving Langevin’s differential equation, we obtained a phase diagram in which each microbiota sample can be placed by its fitness and variability into one of two phases, according to the stability ranking of the system. As shown by the phase space in Fig. 3, three different conditions can occur. The first is a healthy microbiota with some fluctuations, as shown by one of the subjects in the study of Caporaso et al. (48). As this case would have good fitness, its temporal variability would not place the microbiota in the unstable phase of the diagram. Second, we have a subject from the study of Dethlefsen and Relman (49) whose microbiota was perturbed twice by antibiotic intake, undergoing sufficient change to lose its stability and hence be placed in the unstable part. In this location, it is more sensitive to potential perturbations such as, for example, opportunistic infections. In the third and last condition, the subject was already in the unstable phase because of a health issue, i.e., IBS. This can be observed in one of the patients in the study of Durbán et al. (12). In addition, it was shown that this subject’s health status improved during the experiment, implying that his/her microbiota also recovered stability. Interestingly, in the study of David et al. (53), the subject who had a Salmonella infection during the experiment underwent a significant shift in variability with eventual recovery from the perturbed state (see Fig. S3).

Taylor’s law parameter space for intervals concerning gut microbiota in the host lifestyle study (53). We observe that subject B, who suffered a Salmonella infection during the experiment, had a relevant shift in the parameters from “_before” to “_infection” and a final recovery from the perturbed state to “_afterinfec,” which lies in the parameter area compatible with the healthy and stable intervals (see Table S2). Subject A also had a shift in variability from “_before” to “_abroad” and back to “_returned,” also in the proximity zone of healthy and stable periods. Download FIG S3, EPS file, 0.4 MB.

Specifically, in the host lifestyle study (53), the presence of “rank stability islands” among medium-ranked taxa is an interesting feature revealed by the analysis of rank stability in different time periods in subject A. Interestingly, this stability was compromised when the period was not an ordinary one, suggesting that those taxa were sensitive to lifestyle changes. Among the genera identified as “rank stability islands,” Lachnobacterium and Clostridium were catalogued as genera predictive of dysbiosis in the work of Larsen and Dai (65), which analyzed the same data set (53). Furthermore, research has recently confirmed a clear relationship between Actinomyces and conventional adenoma (66), one of the two main precursors of colorectal cancer. Finally, Eggerthella is an opportunistic pathogen often associated with serious gastrointestinal pathology (67).

One might question the role of these taxa as key players in the phase transition of the microbiota and wonder whether they are more susceptible to perturbations than the most abundant taxa. The types of interactions that could sustain this particular behavior are unclear, as these nonabundant taxa are usually excluded from dynamic studies to obtain a community matrix. Further experiments and data analysis are needed to clarify whether “rank stability islands” are a widespread feature of microbiotas and whether they appear at lower taxonomic levels too.

Notwithstanding the above, we should be aware that the hypothesis is too simplistic to directly apply to reality. Indeed, the situation is more complex than the idea that healthy people can be distinguished from nonhealthy people in solely compositional terms, as highlighted by Moya and Ferrer in their recent review (17). There are several feasible scenarios in which we can consider the microbiota to be stable, irrespective of its compositional shifts over time. For example, it may depend on the ability of the microbiota to recover its initial composition (resilience) or its ability to recover its original function despite its composition (functional redundancy). What we have shown in this work could be explained as the transition of a stable microbiota into a state of dysbiosis.

This is a first step toward understanding microbiota stability, although the model presents some limitations and thus further research is required. From a biological perspective, many questions arise from this work. We have observed the same pattern in Taylor’s parameters under all of the conditions studied, but a pertinent question is whether this is really a universal feature in the hugely diverse microbial niches. Furthermore, another relevant question relates to mechanisms involved in the maintenance of population structure. Undoubtedly, the nature of the interactions between community components has a great bearing on this issue, and this is related to community fitness, as mentioned above. How we should address community fitness remains unclear, but studies like the one by Tikhonov (68) could point us in the right direction and help us to unravel the complexity of the microbiota and its relationship to host health.

MATERIALS AND METHODS

Model.

We modeled microbial abundances over time along the lines of Blumm et al. (45). The dynamics of taxon relative abundances was described by the Langevin equation, where F captured the fitness of the taxon i, V corresponded to the noise amplitude, and ξ(t) was a Gaussian random noise with a mean of zero, 〈ξ(t)〉 = 0, and variance uncorrelated in time, 〈ξ(t)ξ(t′)〉 = δ(t′ − t). The function ϕ(t) ensured normalization at all times, ∑x(t) = 1, and corresponded to . The time course changes in the probability that a taxon i has a relative abundance x(t), P(x,t), was determined by the Fokker-Planck equation:

The microbiota evolved toward a steady state with a time-independent probability depending on the values of α, β, F, and V. For α < 1 (otherwise, systems are always unstable), the steady-state probability was localized in a region around a preferred value or broadly distributed over a wide range, depending on whether fitness F dominated or was overwhelmed by the noise amplitude V. The steady-state solution of the Fokker-Planck equation was given by: if 2β ≠ 1 + α and if 2β = 1 + α, where and C and C are integrals solved numerically for the parameters of interest. The ordered phase occurred when the solution had a maximum in the physical interval (0 < x < 1). For larger V values, the transition to a disordered phase happened when the maximum shifted to the unphysical region, x < 0, which set the phase transition region V(α, β, F). The phase transition region was calculated analytically in particular cases, if β = α ≠ 1 and if 2β = 1 + α, where the first case simplified to F = 3V2 if β = 0.75 and the fitness of this taxon dominated in ϕ0. In many physical systems (Brownian motion is the classic example [69]), the two terms of the Langevin equation are related. The fluctuation dissipation theorem states a general relationship between the response to an external disturbance and the internal fluctuations of the system (64). The theorem can be used as the basic formula to derive fitness from the analysis of fluctuations in the microbiota, assuming that it is in equilibrium (the ordered phase).

Standardization.

To properly show all of the studies under common axes, we decided to standardize the Taylor parameters by using the group of healthy subjects for every single study independently. By this approach, all of the studies can be visualized in a shared plot with units of Taylor-parameter standard deviation on their axes.

For a Taylor parameter, e.g., V, the estimate of the mean (V̂) of the healthy subpopulation, composed of h subjects, is as , since ω are normalized weights calculated as where σ is an estimate of the uncertainty in V obtained together with V from the x-weighted power law fit for healthy subjects.

Likewise, the estimation of the standard deviation for the healthy population is where , which finally yields

Selection and methods.

The bacterial and archaeal taxonomic assignments were obtained by analyzing 16S rRNA sequences, which were clustered into operational taxonomic units (OTUs) with 97% sequence identity by QIIME (70). SMS data (51) were analyzed and assigned at the strain level by the Livermore metagenomic analysis toolkit (LMAT) (71), according to their default quality threshold. The genus, with the best balance between error assignment and the number of taxa, was chosen as our reference taxonomic level. We verified that our conclusions were not significantly affected by selecting the family or species as the reference taxonomic level (see Fig. S4).

Overview of the comparison of different approaches based on adjacent taxonomic levels using plots in the Taylor parameter space. The former row of subfigures is for 16S rRNA, where the levels are family (blue circles) versus genus (purple triangles), whereas the latter row of subfigures is for SMS, where the levels are genus (blue circles) versus species (purple triangles). The left column shows the raw results, and the right column plots the standardized results (see standardization in Materials and Methods). Download FIG S4, EPS file, 0.4 MB.

Sample selection.

We chose studies about relevant pathologies containing metagenomic sequencing time data series of bacterial populations from humans in various healthy and nonhealthy states. Only subjects for whom three or more time points of data were available in databases were selected. The study by Caporaso et al. (48) was selected because it featured two healthy subjects measured over a very long time span, with almost daily sampling. The study of Faith et al. (50) was selected because of the body mass index (BMI) differences among the subjects. Moreover, some of them were on diets, which could be treated as system perturbations. Only those subjects whose BMI was normal or who were overweight were considered healthy. The study by Smith et al. (51) was selected for both the age of the patients and the rare disease. Regarding kwashiorkor, we considered only the discordant twins and deemed subjects unaffected by kwashiorkor as being healthy in each pair of patients. The study of David et al. (52) was selected for its differential diets. The healthy period was considered to be the initial samples of each subject before starting the diet, while the remaining time points were considered perturbations. The study of Dethlefsen and Relman (49) was selected because of the interesting treatment of two intakes of the same antibiotic by three different subjects. The healthy period was considered to correspond only to those times before any antibiotic treatment, whereas the periods during and after antibiotic intake were considered perturbations. The work of David et al. (53) was selected because of the comprehensive longitudinal data that it provides plus its complete metadata and the interesting events experienced by both subjects (an infection and a trip abroad). The healthy period was taken from time points before and after each event. Finally, we also considered a study carried out by Durbán et al. (12) in which the healthy subjects were considered those who did not suffer from IBS, while the patients who had this disease were taken as perturbations.

The metadata of each study are provided in Tables S1 to S4. The studies all used 16S rRNA gene sequencing, except for the study of the discordant kwashiorkor twins (51), in which both SMS and 16S rRNA data were used. In the latter case, we chose to work with SMS data to show that our method is valid, regardless of the source of taxonomic information. Each of the data sets was treated as follows.

16S rRNA sequence processing.

Reads from the selected studies were first quality filtered with the FastX toolkit (72), allowing only those reads scoring over 25 for quality in 75% of the complete sequence. 16S rRNA reads with 97% nucleotide sequence identity were then clustered into OTUs with the QIIME software package (70) (version 1.8). We used an open reference OTU picking workflow in all cases. The clustering method used was UCLUST, and the OTUs were matched against the Silva database (73) (version 111, July 2012) and were taxonomically assigned by a UCLUST-based consensus taxonomy assigner. The parameters used in this step were a similarity of 0.97, a prefilter percent identity of 0.6, 20 maximum accepts, and 500 maximum rejects.

Metagenomic sequence processing.

Shotgun metagenomic sequences were analyzed with the LMAT software package (71) (version 1.2.4, with the February 2015 release of the LMAT Grand database). LMAT was run with a Bull shared-memory node belonging to the team’s high-performance computing cluster. It was equipped with 32 cores (64 threads available with Intel hyperthreading technology) since it had two Haswell-based Xeons (22 nm technology), E5-2698v3 2.3-GHz processors, sharing half a tebibyte of DRAM (dynamic random access memory). This node is also provided with a PCIe solid-state drive (SSD) card as NVRAM (nonvolatile random access memory) and the Micron P420m HHHL (half-height half-length) with 1.4 terabytes and 750,000 reading IOPS (input-output operations per second) of 4 kilobytes, achieving 3.3 GB/s. The computing node was supplied with a RAID-0 (striping) scratch disk area. We used the Grand database (74), released in February 2015, provided by the LMAT team, where Grand refers to a huge database that contains k-mers from all of the viral, prokaryotic, fungal, and protist genomes present in the NCBI database, the human reference genome (hg19), GenBank Human, and the 1000 Human Genomes Project (this represents about 31.75 billion k-mers occupying 457.62 GB) (74). Before any calculations were made, the entire database was loaded into the NVRAM. With this configuration, the observed LMAT sustained sequence classification rate was 20 kbp/s/core. Finally, it is worth mentioning that a complete set of Python scripts was developed as back and front ends of the LMAT pipeline to manage the added complexity of time series analysis (https://github.com/DLSteam/MAUS_scripts).

Taxon level robustness.

We selected the genus as the taxonomic level for the subsequent steps of our work. To ensure that there were no crucial differences between adjacent taxonomic levels that could still be of relevance after standardization (see the last subsection of Materials and Methods), we tested two different data sets. In the former, the antibiotic study (49) using 16S rRNA data, we tested the differences between the genus and family levels. The last data set tested was the kwashiorkor discordant-twin study (51) for both genus and species taxonomic levels. See Fig. S4 (overview) and S5 (detail) for plots comparing the studies (and so 16S rRNA and SMS data) and adjacent taxonomic levels.

Detail of the comparison of different approaches based on adjacent taxonomic levels using plots of x-weighted power law fits (see Materials and Methods). The former row of subfigures shows examples for 16S rRNA, whereas the latter row of subfigures plots examples for SMS. The left column shows results for the superior taxonomic level (family for 16S rRNA, genus for SMS), while the right column shows results for the inferior level (genus for 16S rRNA, species for SMS). Download FIG S5, EPS file, 1.9 MB.

x-weighted power law fit.

When fitting the power law of standard versus mean, we took into account that every mean has uncertainty and can be estimated for a sample size n by the SEM (standard error of the mean). Here, the uncertainties affected the independent variable, so the fit was not as trivial as a y-weighted fit, where the uncertainties affect the dependent variable. A standard approach to perform this fit is to (i) invert the variables before applying the weights, (ii) perform the weighted fit, and (iii) revert the inversion. This method is deterministic, but the approximate solution worsens with smaller coefficients of determination. To overcome this limitation, we developed a stochastic method with a bootstrapping-like strategy that avoided inversion and was applicable regardless of the coefficient of determination.

The basic idea of bootstrapping is that inference about a population from sample data can be modeled by resampling the sample data and performing inference about a sample from resampled data (75). To adapt this general idea to our problem, we resampled the x data array by using its error array. That is, for each replicate, a new x data array was computed on the basis of , where v is a Gaussian random variable with mean μ = 0 and standard deviation σ = SEM, as defined previously. For each replicate, a complete unweighted power law fit was performed, where to choose between fitting power laws (y = V·xβ) by linear regression on log-transformed data versus nonlinear regression, we mainly followed the general guidelines for the analysis of biological power laws (76). The parameters of the x-weighted fit were then estimated by averaging through all of the replicate fits performed, and their errors were estimated by computing the standard deviation for all of the fits. At the end of each step, the relative error was calculated by comparing the fit parameter estimation in the last step with the previous one. Finally, both the coefficient of determination of the fit and the coefficient of correlation between the fit parameters were estimated by averaging.

RSI and variability.

The RSI is shown as a percentage in a separate bar on the right of the rank matrix plot in Fig. 4 and 5 (see also Fig. S1 and S2). The RSI is strictly 1 for an element whose range never changes over time and is strictly 0 for an element whose rank oscillates between the extremes from time to time. So, the RSI is calculated, per element, as 1 minus the quotient of the number of true rank hops taken divided by the number of maximum possible rank hops, all to the p power: where D is the total number of rank hops taken by the element studied, N is the number of elements that have been ranked, and t is the number of time samples. The power index, p = 4, was arbitrarily chosen to increase the resolution in the stable region.

Finally, under the rank matrices in Fig. 4 and 5, there are plots relevant to the variability of the rank over time. On the one hand, the RV (rank variability) of a sampling point shows the absolute difference between every taxon’s rank and its accumulated abundance rank (the overall rank), averaged for all of the taxa shown. On the other hand, the DV (difference variability) of a sampling point shows the absolute difference between every taxon’s rank at that time and the value that it had at the previous sampling point, averaged for all of the taxa shown.