Stochastic logistic models reproduce experimental time series of microbial communities
Abstract
We analyze properties of experimental microbial time series, from plankton and the human microbiome, and investigate whether stochastic generalized LotkaVolterra models could reproduce those properties. We show that this is the case when the noise term is large and a linear function of the species abundance, while the strength of the selfinteractions varies over multiple orders of magnitude. We stress the fact that all the observed stochastic properties can be obtained from a logistic model, that is, without interactions, even the niche character of the experimental time series. Linear noise is associated with growth rate stochasticity, which is related to changes in the environment. This suggests that fluctuations in the sparsely sampled experimental time series may be caused by extrinsic sources.
Introduction
Microbial communities are found everywhere on earth, from oceans and soils to gastrointestinal tracts of animals, and play a key role in shaping ecological systems. Because of their importance for our health, humanassociated microbial communities have recently received a lot of attention. According to the latest estimates, for each human cell in our body, we count one microbe (Sender et al., 2016). Dysbiosis in the gut microbiome is associated with many diseases from obesity, chronic inflammatory diseases, some types of cancer to autism spectrum disorder (Gilbert et al., 2016). It is therefore crucial to recognize what a healthy composition is, and if unbalanced, be able to shift the composition to a healthy state. This asks for an understanding of the ecological processes shaping the community and dynamical modeling.
The dynamics of complex ecosystems can be studied by considering the number of individuals of each species, referred to as abundances, at subsequent time points. There are several ways to characterize experimental time series properties. Models typically focus on one specific aspect such as the stability of the community (May, 1972; Coyte et al., 2015; Levine et al., 2017; Grilli et al., 2017; Gavina et al., 2018; Gibbs et al., 2018), the neutrality (Fisher and Mehta, 2014; Washburne et al., 2016), or mechanisms leading to longtailed rank abundance distributions (Solé et al., 2002; Brown et al., 2002; McGill et al., 2007; Matthews and Whittaker, 2015). Different types of dynamical models have been proposed. A first distinction can be made between neutral and nonneutral models. Neutral models assume that species are ecologically equivalent and that all variation between species is caused by randomness. In such models, no competitive or other interactions are included. A second distinction is made between populationlevel and individualbased models. Generalized LotkaVolterra (gLV) models describe the system at the population level and assume that the interactions between species dictate the community’s time evolution. Both deterministic and stochastic implementations exist for gLV models. Stochastic models include a noise term. There are multiple origins of the noise: intrinsic noise captures the fluctuations due to small numbers, extrinsic noise models external factors such as changing immigration rates of species or changing growth rates mediated by a varying flux of nutrients. Individualbased or agentbased models include selforganized instability models (Solé et al., 2002) and the controversial neutral model of Hubbell, 2001; Rosindell et al., 2011. A classification scheme that assesses the relative importance of different ecological processes from time series has been proposed in Faust et al., 2018. The scheme is based on a test for temporal structure in the time series via an analysis of the noise color and neutrality. Applied to the time series of human stool microbiota, it tells us that stochastic gLV or selforganized instability models are more realistic. Here, we will however only focus on stochastic gLV models. The reason for this is twofold. First, one can encompass the whole spectrum of ecosystems from neutral to niche with gLV models (Fisher and Mehta, 2014). Second, we aim at describing dense ecosystems and even though an individualbased model might be more accurate, in the large number limit it will be captured by a Langevin approximation, that is, by the stochastic gLV model.
Our goal is to compare time series generated by stochastic gLV models with experimental time series of microbial communities. We aim at capturing all observed properties mentioned above—the rank abundance profile, the noise color, and the niche character—as well as the statistical properties of the differences between abundances at successive time points with one model. As is shown in Properties of experimental time series, the abundance distribution is heavytailed, which means that few species are highly abundant and many species have low abundances. Despite the large differences in abundances, the ratios of abundances at successive time points and the noise color are independent of these abundances and although the fluctuations are large, the results of the neutrality tests indicate that the experimental time series are in the niche regime. To sum up, we seek growth rates, interaction matrices, immigration rates, and an implementation of the noise in stochastic gLV models to obtain the experimental characteristics.
We simulated time series using gLV equations. The interaction matrices are random as was introduced by May, 1972. The growth rates are determined by the choice of the steadystate, which is set to either equal abundances for all species or abundances according to the rank abundance profiles found for experimental data. For the noise, we consider different implementations corresponding to different sources of intrinsic and extrinsic noise.
Our analysis constrains the type of stochastic gLV models able to grasp the properties of experimental time series. First, we show that there is a correlation between the noise color and the product of the mean abundance and the selfinteraction of a species. The noise color profile for such models will, therefore, depend on the steadystate. This implies that imposing equal selfinteraction strengths for all species, what can be done to ensure stability (Fisher and Mehta, 2014; Gibson et al., 2016), is incompatible with the properties of experimental time series. Second, from the differences between abundances at successive time points, we conclude that a model with mostly extrinsic (linear) noise agrees best with the experimental time series. Third, neutrality tests often result in the niche regime for time series generated by noninteracting species with noise. We, therefore, conclude that all stochastic properties of experimental time series are captured by a logistic model with large linear noise. However, interactions are not incompatible with those properties. This suggests using stochastic logistic models as null models to test for interactions. Our results go along the lines of the ones obtained by Grilli, 2019 which state that the stochastic logistic model can be interpreted as an effective model capturing the statistics of individual species fluctuations.
All codes are available online (see Additional files: Code).
Results
Properties of experimental time series
We study time series from different microbial systems: the human gut microbiome (David et al., 2014), marine plankton (MartinPlatero et al., 2018), and diverse body sites (hand palm, tongue, fecal) (Caporaso et al., 2011; Figure 1A). A study of the different characteristics for a selection of these data is represented in Figure 1. The complete study of all time series can be found in Supplementary file 1: Analysis of experimental data. We propose a detailed description of the properties of the experimental time series. They fall essentially into two categories. The stability and rank abundance are tightly connected to the deterministic part of the equations while the differences between abundances at successive time points and noise color explain the stochastic behavior. The neutrality is more subtle and depends on the complete system.
Box 1.
Definitions of the studied characteristics We study multiple characteristics of the dynamics of microbial communities.
We here define these characteristics. The labels (AF) denote the different figures of Figure 1 and Figure 4.
A time series represents the time evolution of the abundances of different species of the community.
The rank abundance distribution describes the commonness and rarity of all species. It can be represented by a rank abundance plot, in which the abundances of the species are given as a function of the rank of the species, where the rank is determined by sorting the species from high to low abundance. These curves can generally be fitted with power law, lognormal, or logarithmic series functions (Limpert et al., 2001; McGill et al., 2007; Brown et al., 2002).
The noise color describes the distribution of the frequencies of the fluctuations of a time series of a species. It is defined by the slope of a linear fit through the power spectral density. White, pink, brown and black noise correspond to slopes around 0,–1, −2 and −3 respectively. The more negative the slope is—this corresponds to darker noise—the more structure there is in the time series (Faust et al., 2018).
We study the mean absolute difference between abundances at successive time points $\u27e8\mid x(t+\delta t)x(t)\mid \u27e9$ as a function of the mean abundance $\u27e8x(t)\u27e9$. These values represent the jumps of the abundances from one time point to the next.
We measure the ratios of the abundances at two successive time points $x(t+\delta t)/x(t)$. The advantage of this method is that it captures the direction of a jump between two time points: for ratios higher than one the jump is positive, for ratios lower than one the jump is negative. The distribution of these ratios fits a lognormal curve with a mean at one as the fluctuations occur around steadystate and the width of the distribution tells how large the fluctuations of a time series are. The goodness of the fit is defined by the pvalue of the KolmogorovSmirnov test. Higher pvalues denote a better fit. We use the width as a characteristic and compare the widths of different species. Examples of the fitted lognormal curve can be found in Supplementary file 1: Supporting results.
The KullbackLeibler divergence measures how different the multivariate distribution of species abundances is from a distribution constructed under the assumption of ecological neutrality. The idea of the neutral covariance test is to compare the time series with a WrightFisher process. A WrightFisher process is a continuous approximation of Hubbell’s neutral model for a large and finite community. In particular, it tests the invariance with respect to grouping. More about the validity of these neutrality measures can be found in the Supplementary file 1: Supporting results.
The time series show fluctuations over time
The experimental time series show large fluctuations over time. We can ask the question whether the origin of this variation is biological or technical, and assume that most of the variation can be contributed to biological processes. This hypothesis is supported by the results of Silverman et al., 2018 for microbial communities of an artificial gut. Here, the biological variation becomes five to six times more important than the technical variation for the sampling interval of a day. Also, Grilli, 2019 shows the time correlation of experimental time series which is nonzero. In the case where the variation is mostly due technical errors, we expect to see no correlation. Because no experimental errorbars are available for most of the data and because we assume most variation has a biological origin, we did not consider the errors on the species abundances.
The abundance distribution is heavytailed
The first aspect of community modeling that has been widely studied during the last years is the stability of the steadystates. Large random networks tend to be unstable (May, 1972). This problem is often solved by considering only weak interactions, sparse interaction matrices (May, 2001) or by introducing higherorder interactions (Grilli et al., 2017; Gavina et al., 2018; Sidhom and Galla, 2019). Although the stability of gLV models decreases with an increasing number of participating species, the stability only depends on the interaction matrix and not on the abundances (Gibbs et al., 2018). The abundance distribution of the experimental data is heavytailed. This means that there are few common and many rare species. The distribution of the steadystate values can also be represented by a rank abundance curve (see Box 1B). Although the abundances show large fluctuations over time, the rank abundance remains stable (Figure 1B).
The differences between abundances at successive time points are large and linear with respect to the species abundance
Time series can be described by the differences between abundances at successive time points. We propose to focus on two specific representations of the information contained in those differences. First, we consider the mean absolute difference between abundances at successive time points $\u27e8\mid x(t+\delta t)x(t)\mid \u27e9$ as a function of the mean abundance $\u27e8x(t)\u27e9$ (see Box 1D). For the experimental data, the relation between these variables is a monomial—this means that it is linear on the loglog scale (Figure 1D). The fact that the slope of this line is almost one hints at a linear nature of the noise.
Second, we examine the distribution of the ratios of the abundances at two successive time points $x(t+\delta t)/x(t)$ (see Box 1E). The width of this distribution tells how large the fluctuations are. To measure this width, we fit the distribution with a lognormal curve for which the mean is fixed to be one as the fluctuations occur around steadystate. For most of the species of experimental data (except for the stool data), the fit of the distribution to a lognormal curve is good (Figure 1E). Furthermore, we notice that the distribution is wide—in the order of 1—and that the width does not depend on the mean abundance of the species (Figure 1E).
The noise color is independent of the mean abundance of the species
The noise of a time series can be studied by considering the distribution of the frequencies of the fluctuations. This distribution can be defined by its slope, which is interpreted as the noise color (see Box 1C). We notice that there is no correlation between the noise color and the mean abundance of the species for experimental time series (Figure 1C).
Experimental time series are in the niche regime
In neutral theory, it is assumed that all species or individuals are functionally equivalent. It is challenging to test whether a given time series was generated by neutral or niche dynamics. We use two definitions of neutrality measures: the KullbackLeibler divergence as used in Fisher and Mehta, 2014 and the neutral covariance test as proposed by Washburne et al., 2016 (see Box 1F). Both neutrality measures indicate that most experimental time series are in the niche regime (Figure 1F).
Reproducing properties of experimental time series from stochastic generalized LotkaVolterra models
We find that the aforementioned characteristics of experimental time series can be reproduced by stochastic logistic equations. We first explain how to choose the growth rate to obtain the heavytailed experimental abundance distribution. Next, we discuss how the noise color determines the selfinteraction of a species given its abundance and how the implementation of the noise determines the slope of the mean absolute increment $\u27e8\mid x(t+\delta t)x(t)\mid \u27e9$ and the mean abundance $\u27e8x(t)\u27e9$ (such as in Figure 1D). In the end, by using the appropriate choice for the selfinteractions, growth rates, and noise implementation, we conclude that a stochastic logistic model can reproduce all the stochastic properties, including the niche regime for the neutrality tests although the model does not include any interactions.
The rank abundance distribution can be imposed by fixing the growth rate
Random matrix models do typically not give rise to heavytailed abundance distributions. Neither is it known which properties of the interaction matrix and growth rates are required to obtain a realistic rank abundance distribution. We can however enforce the desired rank abundance artificially by solving the steadystate of the gLV equations. Given the steadystate abundance vector ${\overrightarrow{x}}^{*}$ and interaction matrix ω, we impose the growth rate $\widehat{g}=\omega {\overrightarrow{x}}^{*}$. One model that results in heavytailed distributions is the selforganized instability model proposed by Solé et al., 2002.
For logistic models, the growth rate is equal to the product of the selfinteraction and mean abundance. The noise color and the width of the distribution of ratios $x(t+\delta t)/x(t)$ depend on this product. To obtain given characteristics—a predefined noise color and width of the distribution of ratios $x(t+\delta t)/x(t)$—the choice of the growth rate will dictate the choice of the remaining free parameters, the sampling time step δt and the noise strength σ.
The noise color is determined by the mean abundance and the selfinteraction of the species
To study the noise color, we first consider a model where the species are not interacting. The noise color is independent of the implementation of the noise but depends on the product of the mean abundance and the selfinteraction of the species (Figure 2A). For noninteracting species, the growth rate equals the product of the selfinteraction and the steadystate abundance. Because we consider fluctuations around steadystate, the mean and the steadystate abundance are nearly equal and the xaxis of Figure 2A; Figure 2B; Figure 2C; can be interpreted as the growth rate. Also, the strength of the noise does not change its color (Figure 2C). A parameter that is important for the noise color is the sampling rate: the higher the sampling frequency the darker the noise becomes (Figure 2B). This is in agreement with the results of Faust et al., 2018. Darker noise corresponds to more structure in the time series. The more frequent the abundances are sampled the more details are visible and the underlying interactions become more visible. We conclude that the noise color is only dependent on the mean abundance, the selfinteractions, and the sampling rate. Figures of the dependence on the mean abundance and selfinteraction separately can be found in Supplementary file 1: Supporting results.
For interacting species, increasing the strength of the interactions makes the color of the noise darker in the high mean abundance range (Figure 2D; Figure 2E). Importantly, for interacting species with a lognormal rank abundance, the correlation between the noise color and mean abundance is preserved (Figure 2E). The data can be fit to obtain a bijective function between the product of the mean abundance and the selfinteraction, and the noise color. Assuming this model is correct, we can obtain an estimate for the selfinteraction coefficients given the mean abundance and noise color by fixing the sampling rate and the interaction strength. The uncertainty on the estimates is larger where the fitted curve is more flat (slopes of the power spectral density around −1.7 and 0), but many experimental values of the stool microbiome data lie in the pink region where the selfinteraction can be estimated for this model.
The implementation of the noise determines the correlation between the mean absolute increment $\u27e8\mid x(t+\delta t)x(t)\mid \u27e9$ and the mean abundance $\u27e8x(t)\u27e9$
Next, we study the differences between abundances at successive time points (see Figure 1D). From the results of the noise color, we can estimate the selfinteraction for the dynamics of the experimental data. We use the rank abundance and the selfinteraction inferred from noise color of the microbiome data of the human stool to perform simulations and calculate the characteristics of the distribution of differences between abundances at successive time points. We here assume that there are no interactions. More results for dynamics with interactions are in Supplementary file 1: Supporting results. We first study the correlation between the mean absolute difference between abundances at successive time points $\u27e8\mid x(t+\delta t)x(t)\mid \u27e9$ and the mean abundance $\u27e8x(t)\u27e9$. For linear multiplicative noise, the slope of the curve of the logarithm of the mean absolute difference between abundances at successive time points ${\mathrm{log}}_{10}\left(\u27e8\mid x(t+\delta t)x(t)\mid \u27e9\right)$ as a function of the logarithm of the mean abundance ${\mathrm{log}}_{10}\left(\u27e8x(t)\u27e9\right)$ is one. For multiplicative noise that scales with the square root of the abundance, the slope is around 0.66 and for additive noise, the slope is zero. By combining both linear noise and noise that scales with the square root of the abundance, slopes with values between 0.6 and 1 can be obtained (Figure 3A). The slopes of experimental data range between 0.84 and 0.99, we therefore conclude that linear noise is a relatively good approximation to perform stochastic modeling of microbial communities.
The strength of the noise determines the width of the distribution of ratios $x(t+\delta t)/x(t)$
Next, we examine the distribution of the ratios of abundances at successive time points (see Box 1E). As expected, for significant noise, this distribution can be approximated by a lognormal curve and the width of the distribution becomes larger for increasing noise strength (Figure 3B). In order to have widths that are of the same order of magnitude as the ones of the experimental data, the noise must be sufficiently strong. Another way of increasing the width is through interactions, this effect is only moderate. These results are presented in Supplementary file 1: Supporting results.
Stochastic logistic models capture the properties of experimental time series
By using all previous results and imposing the steadystate of experimental data, we find that it is possible to generate time series with identical characteristics to the ones seen in the experimental time series (Figure 4). Furthermore, these time series can be generated without introducing any interaction between the different species, but their neutrality measures can still be in the niche regime (Figure 4F). Out of 100 simulations, 62 had a pvalue smaller than 0.05 for the neutral covariance test which means they are in the niche regime. The colors of the noise fix the selfinteraction values (Figure 4C), next the rank abundance distribution is imposed by calculating the growth vector $\widehat{g}$ (Figure 4B). The slope of the curve of the mean absolute difference between abundances at successive time points as a function of the mean abundance is one by using linear multiplicative noise (Figure 4D) and the width of the fluctuations is tuned by choosing a large noise size σ (Figure 4E). In most experimental time series, only the fractional abundances of species can be measured per time point and not the absolute ones. Because the total abundance of all species remains nearly constant in time series generated by a stochastic logistic equation, our results still hold for time series with fractional abundances (see Supporting results). Similar results can be obtained for models with interactions (see Supporting results), but we want to stress that interactions are not needed to reproduce the properties of experimental time series.
Discussion
Recent research has focused on different aspects of experimental time series of microbial dynamics, in particular the rank abundance distribution, the noise color, the stability, and neutrality. Within the framework of stochastic generalized LotkaVolterra models, we studied the influence of growth rates, interactions between species, and the different sources of stochasticity on the observed characteristics of the noise and on neutrality. Our observations are:
Even when we consider the case without interactions between species, the result of the neutrality test on the time series is often niche. We should, therefore, be careful in the interpretation of the results of neutrality tests.
For a given sampling step δt, the noise color depends on the product of the selfinteraction and the mean abundance, which for noninteracting species reduces to a dependence on the growth rate. Assuming the model can be used for microbial communities, the selfinteraction coefficients can be estimated given the mean abundance, noise color, and sampling rate. Low sampling rates result in larger errors (Figure 2B). For sparsely sampled experimental data, the standard deviation of the selfinteraction inferred using the noise color will be larger. For the experimental time series (plankton, gut, and human microbiome) the selfinteraction strengths range over several orders of magnitude. The convention of equalling all selfinteractions to −1 used in several studies (Fisher and Mehta, 2014; Gibson et al., 2016), cannot be adopted for stochastic models of communities with a heavytailed abundance distribution.
The exponent of the mean absolute differences between abundances at successive time points with respect to the mean abundances is slightly smaller than one for experimental time series. Linear multiplicative noise results in a value of one, square root noise results in lower values (0.6). A mix of linear and square root noise can result in slopes with intermediate values.
A large multiplicative linear noise is in agreement with both the distribution of the ratios of abundances at successive time points and the relation between the differences between abundances at successive time points and mean values.
To conclude, characteristics of experimental time series, from plankton to gut microbiota, can be reproduced by stochastic logistic models with a dominant linear noise. We expect, however, that for higher sampling rates, modeling the interactions between microbes would be necessary to explain the properties of the time series. For gut microbial time series, the system is sampled only once a day and therefore dominated by the noise in the growth terms corresponding to a linear noise.
Predictive models for the dynamics of microbial communities will certainly require a more indepth description of the system. Nutrients and spatial distribution of microbes should play a role to dictate the evolution of the community, as well as the interaction with the environment. Synthetic microbial communities are currently being developed and will hopefully provide a more comprehensive view on the complexity of microbial communities (Vrancken et al., 2019).
Materials and methods
Modeling generalized LotkaVolterra equations
Request a detailed protocolIn a microbial community different species interact because they compete for the same resources. Moreover, they produce byproducts that can affect the growth of other species. Depending on the nature of the byproducts, harmful, beneficial, or even essential, the interaction strength will be either negative or positive. To describe the dynamics of interacting species, one can use the generalized LotkaVolterra equations:
where x_{i}, λ_{i} and g_{i} are the abundance, the immigration rate, and the growth rate of species i respectively, and ${\omega}_{ij}$ is the interaction coefficient that represents the effect of species $j$ on species i. The diagonal elements of the interaction matrix ${\omega}_{ii}$, the socalled selfinteractions, are negative to ensure stable steadystates. The offdiagonal elements of the interaction matrix ${\omega}_{ij}$ are drawn from a normal distribution with standard deviation α (${\omega}_{ij}\sim \mathcal{N}(0,{\alpha}^{2})$). The gLV equations only consider pairwise effects and no saturation terms, or other higherorder terms. Due to this drawback, these models sometimes fail to predict microbial dynamics (Momeni et al., 2017; Levine et al., 2017). However, they are among the most simple models for interacting species and therefore widely studied and used. Noninteracting species can be described by the logistic model, which is a special case of the gLV model obtained by setting all offdiagonal elements of the interaction matrix to zero.
Implementations of the noise
Request a detailed protocolThere exist two principal types of noise: intrinsic and extrinsic noise. Extrinsic noise arises due to external sources that can alter the values of the different variables: the immigration rate and growth rate fluctuate in time through colonization of species or a changing flux of nutrients. These processes give rise to additive and linear multiplicative noise respectively. The remaining parameters, inter and intraspecific interactions can also, change depending on the environment. The formulation of this noise is more subtle (used in Zhu and Yin, 2009). Intrinsic noise is due to the discrete nature of individual microbial cells. Thermal fluctuations at the molecular level determine the fitness of the individual cells. Therefore, cell growth, cell division, and cell death can be considered as stochastic Poisson processes. For large numbers of microbes, these fluctuations will be averaged out.
We first consider the extrinsic noise. If the time series is calculated by ${x}_{i}(t+dt)={x}_{i}(t)+d{x}_{i}(t)$, the implementation of the linear multiplicative noise is as follows,
where dW is an infinitesimal element of a Brownian motion defined by a variance of dt ($dW\sim \sqrt{dt}\mathcal{N}(0,1)$). Changes in immigration rates of microbial species can be modeled with additive noise,
with $d{W}_{\text{const}}\sim \sqrt{dt}\mathcal{N}(0,1)$. Our main motivation is to model the gut microbiome in the colon. Here, we ignore the immigration of species for two reasons. First, the number of microbes in the colon is orders of magnitude larger than the number of microbes in the other parts of the gut (Marteau et al., 2001; Gorbach et al., 1967)—therefore, the flux of incoming microbes in the colon is small. Second, we only consider systems around steadystate, for which we assume immigration does not play an important role. For perturbed systems, which are far from equilibrium, immigration rates cannot be ignored. Ignoring immigration may be too restrictive for some microbial systems such as the skin microbiome or plankton.
To derive the form of intrinsic noise in generalized LotkaVolterra equations, we can consider every species abundance making a random walk in one dimension. The average displacement is zero and the variance of displacement is the sum of the rate of growth (jumping to the right) and the rate of death (jumping to the left). For the generalized LotkaVolterra equations, this results in a noise term
with ω the interaction matrix and where functions f and h each decouple the growth and death terms. In the generalized LotkaVolterra model no difference is made between negative interactions as a result of slowing down the growth rate or increasing the death rate, only the resulting net rates are used. This distinction must however be made to implement the intrinsic noise for gLV. In our analysis, we use the simpler logistic models where the resulting variance of the noise is proportional to the square root of the abundance $\sqrt{x}$. One must be careful not to use this noise for values that are smaller than one because this derivation relies on Poisson statistics which is defined for integer numbers.
We implement the intrinsic noise by a term that scales with the square root of the species abundance (Walczak et al., 2012; Fisher and Mehta, 2014),
with $d{W}_{\text{sqrt}}$ again an infinitesimal element of a Brownian motion defined by a variance of dt ($d{W}_{\text{sqrt}}\sim \sqrt{dt}\mathcal{N}(0,1)$). The size of this noise ${\sigma}_{i,\text{sqrt}}$ is determined by the cell division (g^{+}) and death rates (g^{}) separately, which are in our model combined to one growth vector ($g={g}^{+}{g}^{}$, ${\sigma}_{i,\text{sqrt}}=\sqrt{{g}^{+}+{g}^{}}$), for large division and death rates the intrinsic noise will be larger.
To sum up, we focus on linear multiplicative noise because: (a) extrinsic noise is dominant as microbial communities contain a very large number of individuals and (b) we ignore the immigration of individuals in our analysis.
We verified that our analysis is robust with respect to the multiple possibilities for the discretization of these models. We also compare our populationlevel approach with individualbased modeling approaches. Details can be found in the Supplementary file 1: Supporting results.
Neutrality measures
Request a detailed protocolThere is no consensus on the definition of neutrality. In general, ecosystems are considered neutral if the dominating cause of fluctuations is random birth and death processes and not fitness advantages of species.
Different neutrality measures focus on different aspects of neutrality. The KullbackLeibler divergence verifies whether all species are equal (equal abundances and equal covariances). The neutrality covariance test studies the grouping invariance of species in time series.
Given two distributions P and Q, the KullbackLeibler divergence is defined as
where ${E}_{P}$ is the expectation value using the probabilities of distribution P. The density function of a multivariate Gaussian distribution is
where μ and K are the mean and covariance matrix of the distribution respectively. The KullbackLeibler divergence for two multivariate Gaussian distributions in ${\mathbb{R}}^{n}$ is (Duchi, 2007)
For every time series, we can calculate the mean μ and covariance matrix K, and define values ${\mu}_{N}$ and ${K}_{N}$ for a corresponding neutral time series in which all species are equal (Fisher and Mehta, 2014). The distance to neutrality ${D}_{KL}(P{P}_{N})$ can thus be calculated by computing the probability distribution of the original time series P and the associated neutral distribution ${P}_{N}$ with mean values ${\mu}_{N}={S}^{1}{\sum}_{i=1}^{S}{\mu}_{i}$ and ${K}_{P,ii}={S}^{1}{\sum}_{i=1}^{S}{K}_{ii}$ and ${K}_{P,ij}={S}^{1}{(S1)}^{1}{\sum}_{i=1}^{S}{\sum}_{j=1,i\ne j}^{S}{K}_{ij}$ with S the number of species.
The neutral covariance test was designed by Washburne et al., 2016. We used a python translation of the code developed by this author.
Noise color
Request a detailed protocolThe color of the noise in a time series is determined by the slope of the power spectral density in a loglog scale. This slope can be determined by a linear fit through the spectrum. A different technique to estimate this slope has been put forward by Faust et al., 2018. There, it is argued that the power spectral density does not have a constant slope and that, therefore, a nonlinear curve must be fitted. They choose a spline fit and consider the minimal value of its derivative to be the value of the noise color. Because the minimal value of the slope of the fit is taken, the noise color tends to be darker when using this technique. For our time series, however, we see that the spline fit only deviates from the linear fit for low frequencies (Figure 5). We ignore the low frequencies for fitting because of the windowing effect. Therefore, we opt for a linear fit after omitting the values for low frequencies (one order of magnitude of the lowest frequencies).
The correspondence between the colors and slopes is here:
Slope  Color 

0  white 
1  pink 
2  brown 
3  black 
References

1
The fractal nature of nature: power laws, ecological complexity and biodiversityPhilosophical Transactions of the Royal Society of London. Series B: Biological Sciences 357:619–626.https://doi.org/10.1098/rstb.2001.0993

2
Moving pictures of the human microbiomeGenome Biology 12:R50.https://doi.org/10.1186/gb2011125r50
 3
 4

5
Logistic models, version 37398e9Github.
 6
 7
 8
 9
 10

11
On the origins and control of community types in the human microbiomePLOS Computational Biology 12:e1004688.https://doi.org/10.1371/journal.pcbi.1004688
 12

13
Studies of intestinal microfloraGastroenterology 53:874–880.https://doi.org/10.1016/S00165085(19)341241
 14
 15

16
The Unified Neutral Theory of Biodiversity and Biogeography. No. 32 in Monographs in Population BiologyPrinceton University Press.
 17
 18

19
Comparative study of bacterial groups within the human cecal and fecal MicrobiotaApplied and Environmental Microbiology 67:4939–4942.https://doi.org/10.1128/AEM.67.10.49394942.2001
 20

21
REVIEW: On the species abundance distribution in applied ecology and biodiversity managementJournal of Applied Ecology 52:443–454.https://doi.org/10.1111/13652664.12380
 22

23
Stability and Complexity in Model Ecosystems. 1st Princeton Landmarks in Biology Ed Ed. Princeton Landmarks in BiologyPrinceton University Press.
 24
 25

26
The unified neutral theory of biodiversity and biogeography at age tenTrends in Ecology & Evolution 26:340–348.https://doi.org/10.1016/j.tree.2011.03.024
 27
 28
 29

30
Selforganized instability in complex ecosystemsPhilosophical Transactions of the Royal Society of London. Series B: Biological Sciences 357:667–681.https://doi.org/10.1098/rstb.2001.0992

31
Synthetic ecology of the human gut microbiotaNature Reviews Microbiology 17:754–763.https://doi.org/10.1038/s4157901902648

32
Analytic Methods for Modeling Stochastic Regulatory NetworksIn: X Liu, M. D Betterton, editors. Computational Modeling of Signaling Networks, 880. Totowa, NJ: Humana Press. pp. 273–322.
 33

34
On competitive Lotka–Volterra model in random environmentsJournal of Mathematical Analysis and Applications 357:154–170.https://doi.org/10.1016/j.jmaa.2009.03.066
Decision letter

Sandeep KrishnaReviewing Editor; National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India

Aleksandra M WalczakSenior Editor; École Normale Supérieure, France
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This study shows that a very simple noninteracting model, with an appropriate choice of noise, can explain various summary statistics of microbial community time series. This serves as an important null hypothesis, that can be used to rebut overinterpretations of such time series data when they are used to infer underlying interactions.
Decision letter after peer review:
Thank you for submitting your article "Stochastic logistic models reproduce experimental time series of microbial communities" for consideration by eLife. Your article has been reviewed by Aleksandra Walczak as the Senior Editor, a Reviewing Editor, and two reviewers. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID19 (https://elifesciences.org/articles/57162). Specifically, when editors judge that a submitted work as a whole belongs in eLife but that some conclusions require a modest amount of additional new data, as they do with your paper, we are asking that the manuscript be revised to either limit claims to those supported by data in hand, or to explicitly state that the relevant conclusions require additional supporting data.
Our expectation is that the authors will eventually carry out the additional experiments and report on how they affect the relevant conclusions either in a preprint on bioRxiv or medRxiv, or if appropriate, as a Research Advance in eLife, either of which would be linked to the original paper.
Summary:
Both reviewers found points of interest in your manuscript, and appreciated the usefulness of a simple model that is consistent with the data. However, both have also raised some important issues which need to be addressed, which range from clarity of the presentation, to questions about the methods, and robustness of the results. Please address each of the issues raised below by the reviewers in your revised manuscript.
In addition, it is clear from the reviews and our reading that you have not ruled out other more complex models (e.g. models with interactions). Therefore, we would like you to remove or rewrite statements that make overly strong claims about the real systems, and be more clear which claims are about your model and which about the real system. For instance, the last sentence of the Abstract – "This suggests that fluctuations in the sparsely sampled experimental time series are caused by extrinsic sources." – is too strong a claim to make about the real system. Another example from the Abstract is "We show that the noise term should be large and that it is a linear function of the species abundance, while the strength of the selfinteractions should vary.…", which should be modified to make it clearer that these are requirements for the model to fit the data.
Reviewer #1:
Deschmeemaeker et al., use a stochastic generalized LotkaVolterra model to demonstrate that the properties of human microbes can be explained using a logistic modeling approach. Moreover, their results suggest that the niche state of community dynamics can be explained by the intraspecies fluctuations!, this outcome I find particularly interesting. The manuscript highlights a novel approach to studying complex microbial communities and the findings are generally interesting. However, I would suggest some changes. Mostly, these concerns can be considered minor, but answering them is crucial for more transparency and clarity on the work presented in this manuscript.
1) Palm and tongue, as well as gut microbiomes, change over time. I would have imagined that these changed cause changes in the rank abundance profile of a community over time. The data, therefore, is interesting (because it is derived from humanassociated microbial communities), though not surprising.
2) Readers will benefit if the authors can include information from the methods sections. Such as definitions of noise, noise color, abundance and time series can be included in the main text. This change will help the message from this manuscript to be more accessible to a broader audience.
3) Since the frequency of sampling has an effect on noise color. Can authors elaborate a bit more ( subsection “The selforganized instability model can be reproduced by the stochastic generalized LotkaVolterra model”) on how the limitation of having experimental time series from only a handful of time points might influence their inference?
4) Statistics: Figure 1(AE), Figure 4(A,B,D,E,F) error not mentioned.
5) Figure 1A, legends for this figure is not clear to me. These details are part of the main text however it is important that legends of the figure are described in greater detail. For example, what is the meaning of different colors in the time series? Are these randomly selected species? or species with specific ranks are selected for this analysis?
Reviewer #2:
This study aims to find the simplest model which can explain various summary statistics of microbial community time series. The key result is that various statistics of interest can be obtained using a noninteracting logistic model with appropriate choice of noise. This serves as a form of null hypothesis, and could be used to rebut overinterpretations of such time series.
Overall, I am sympathetic to the goals of the study. I have the following suggestions:
Specific comments:
1) For nonspecialists it would be good to start with a summary of the neutral vs niche regime, as this point is rather central to the discussion. (If the time series seem to be operating in the niche regime, then it is not surprising per se that a noninteracting model can capture some summary statistics of the time series.)
2) This comes across mainly a study of a model, rather than a study of datasets which produces specific insights into microbial communities of different types. This is partly the point of the paper, that the real time series are somehow not sufficient to reveal the true interactions. Of course, the interactions in any real system are neither zero, nor equal, nor randomly sampled, but can be discerned with sufficiently sensitive experiments. So one expects that, after full study of this type, the noninteracting model proposed by the authors will not be supported. The authors do note that interactions cannot be ruled out in their approach. But could some text be included to say what the reader should take away from the study: is it meant to be a null model, or predictive?
3) The results are obtained for a particular choice of growth rates to match given abundances. This particular choice is indeed shown to replicate the observed time series summary statistics. However, it does mean that the other observations made by the authors (specific relations of noise / abundances etc.) could depend on this growth rate choice. How robust are the various observations to changing the particular mechanism of setting abundances?
4) I had difficulty distinguishing between terms "dt", "δt", "δ", "delay", and "sampling time" which the authors seem to use interchangeably through the text and supplement. For example, in Figure 2B should it be "Sampling dt" or "Sampling δt"? In Figure S4 they use the term "δ = 0" which I think should be "δt = 0"? The term "time delay" is not defined elsewhere. More broadly, in a correct implementation of a stochastic simulation, the result should be asymptotically independent of integration timestep. The authors provide detailed comparison of different discretizations, this is appreciated (Figure S5). However, in the supplement I could not see how the "dt" in the SDE is related to the integration time for numerical solution of the Langevin equation "δt", unless they were the same term. I could not at all follow Figure 2B. Is the "sampling time" different from the numerical integration timestep? Sampling is correctly done by first integrating with a small timestep and showing results independent of timestep, and then sampling at the desired sampling intervals. If the sampling is related to the bandwidth limitation, that should be analyzed separately from the numerical integration.
https://doi.org/10.7554/eLife.55650.sa1Author response
Summary:
Both reviewers found points of interest in your manuscript, and appreciated the usefulness of a simple model that is consistent with the data. However, both have also raised some important issues which need to be addressed, which range from clarity of the presentation, to questions about the methods, and robustness of the results. Please address each of the issues raised below by the reviewers in your revised manuscript.
In addition, it is clear from the reviews and our reading that you have not ruled out other more complex models (e.g. models with interactions). Therefore, we would like you to remove or rewrite statements that make overly strong claims about the real systems, and be more clear which claims are about your model and which about the real system. For instance, the last sentence of the Abstract – "This suggests that fluctuations in the sparsely sampled experimental time series are caused by extrinsic sources." – is too strong a claim to make about the real system. Another example from the Abstract is "We show that the noise term should be large and that it is a linear function of the species abundance, while the strength of the selfinteractions should vary.…", which should be modified to make it clearer that these are requirements for the model to fit the data.
We changed the overly strong claims of the Abstract.
Reviewer #1:
Deschmeemaeker et al., use a stochastic generalized LotkaVolterra model to demonstrate that the properties of human microbes can be explained using a logistic modeling approach. Moreover, their results suggest that the niche state of community dynamics can be explained by the intraspecies fluctuations!, this outcome I find particularly interesting. The manuscript highlights a novel approach to studying complex microbial communities and the findings are generally interesting. However, I would suggest some changes. Mostly, these concerns can be considered minor, but answering them is crucial for more transparency and clarity on the work presented in this manuscript.
1) Palm and tongue, as well as gut microbiomes, change over time. I would have imagined that these changed cause changes in the rank abundance profile of a community over time. The data, therefore, is interesting (because it is derived from humanassociated microbial communities), though not surprising.
Indeed, although the species abundances fluctuate over time, the rank abundance profile remains stable over time, as can be seen in Figure 1B. A comment to acknowledge this fact was added to the text (Results section).
2) Readers will benefit if the authors can include information from the methods sections. Such as definitions of noise, noise color, abundance and time series can be included in the main text. This change will help the message from this manuscript to be more accessible to a broader audience.
We added sentences introducing the concepts ’abundance’ and noise in the Introduction. The technical aspects of noise can be found in subsection “Implementations of the noise”. We also created a feature box with the definitions of all the concepts, such that the reader can find this information easily. The features are labeled with the letters of Figure 1 and Figure 4, of which we changed the order to remain consistent.
3) Since the frequency of sampling has an effect on noise color. Can authors elaborate a bit more (subsection “The selforganized instability model can be reproduced by the stochastic generalized LotkaVolterra model”) on how the limitation of having experimental time series from only a handful of time points might influence their inference?
Low sampling rates result in larger errors (Figure 2B). For sparsely sampled experimental data, the standard deviation of the selfinteraction inferred using the noise color will be larger (Discussion section).
4) Statistics: Figure 1(AE), Figure 4(A,B,D,E,F) error not mentioned.
We did not represent experimental error bars on Figure 1(AB) as this data is not available. We expect the errors due to measurements or sampling to be smaller than the fluctuations due to biological processes. This hypothesis is supported by the results of [Silverman et al., 2018] for microbial communities of an artificial gut. Here, the biological variation becomes five to six times more important than the technical variation for the sampling interval of a day (Figure 3B of [Silverman et al., 2018]). Also, [Grilli, 2019] shows the time correlation of experimental time series which is nonzero while if technical variation is the main source of variation, no correlation should be noticeable. We added a paragraph in the manuscript about these errorbars (Results section). For Figure 1C, we added the error of the slope of the power spectral density (noise color). For Figure 1D, we added the error on the exponents. For Figure 1E, the pvalue corresponds to the KolmogorovSmirnov test. Concerning Figure 4, similar errors have been added. We did not include errors for Figure 4(AB) as those figures correspond to one simulation.
5) Figure 1A, legends for this figure is not clear to me. These details are part of the main text however it is important that legends of the figure are described in greater detail. For example, what is the meaning of different colors in the time series? Are these randomly selected species? or species with specific ranks are selected for this analysis?
In the first submission, we selected species with evenly spread ranks. In the figure of this submission, we used another selection of species with ranks 1, 5, 10, 50 and 100. We added a legend denoting the rank number of all species.
Reviewer #2:
This study aims to find the simplest model which can explain various summary statistics of microbial community time series. The key result is that various statistics of interest can be obtained using a noninteracting logistic model with appropriate choice of noise. This serves as a form of null hypothesis, and could be used to rebut overinterpretations of such time series.
Overall, I am sympathetic to the goals of the study. I have the following suggestions:
Specific comments:
1) For nonspecialists it would be good to start with a summary of the neutral vs niche regime, as this point is rather central to the discussion. (If the time series seem to be operating in the niche regime, then it is not surprising per se that a noninteracting model can capture some summary statistics of the time series.)
We added extra sentences about the concept of niche and neutrality (Introduction). A negative neutrality test is sometimes used to state that there must be interactions between species which are subsequently inferred, as for instance in [Faust et al., 2018]. We say that a noninteracting model can be classified in the niche regime by a neutrality test.
2) This comes across mainly a study of a model, rather than a study of datasets which produces specific insights into microbial communities of different types. This is partly the point of the paper, that the real time series are somehow not sufficient to reveal the true interactions. Of course, the interactions in any real system are neither zero, nor equal, nor randomly sampled, but can be discerned with sufficiently sensitive experiments. So one expects that, after full study of this type, the noninteracting model proposed by the authors will not be supported. The authors do note that interactions cannot be ruled out in their approach. But could some text be included to say what the reader should take away from the study: is it meant to be a null model, or predictive?
The logistic model with linear noise can be considered as a null model (Introduction). We do not claim that there are no interactions between species. We show that logistic models are able to reproduce most of the characteristics of experimental time series, what implies that no interactions are needed in the model to obtain the statistical properties.
3) The results are obtained for a particular choice of growth rates to match given abundances. This particular choice is indeed shown to replicate the observed time series summary statistics. However, it does mean that the other observations made by the authors (specific relations of noise / abundances etc.) could depend on this growth rate choice. How robust are the various observations to changing the particular mechanism of setting abundances?
For logistic models, the growth rate is equal to the product of the selfinteraction and mean abundance. The noise color and the width of the distribution of ratios x(t + δt)/x(t) depend on this product. In order to obtain given characteristics—a predefined noise color and width of the distribution of ratios x(t+δt)/x(t)—the choice of the growth rate will therefore dictate the choice of the remaining free parameters, the sampling time step δt and the noise strength σ. This clarification was added to the main paper (Results section).
4) I had difficulty distinguishing between terms "dt", "δt", "δ", "delay", and "sampling time" which the authors seem to use interchangeably through the text and supplement. For example, in Figure 2B should it be "Sampling dt" or "Sampling δt"? In Figure S4 they use the term "δ = 0" which I think should be "δt = 0"? The term "time delay" is not defined elsewhere. More broadly, in a correct implementation of a stochastic simulation, the result should be asymptotically independent of integration timestep. The authors provide detailed comparison of different discretizations, this is appreciated (Figure S5). However, in the supplement I could not see how the "dt" in the SDE is related to the integration time for numerical solution of the Langevin equation "δt", unless they were the same term. I could not at all follow Figure 2B. Is the "sampling time" different from the numerical integration timestep? Sampling is correctly done by first integrating with a small timestep and showing results independent of timestep, and then sampling at the desired sampling intervals. If the sampling is related to the bandwidth limitation, that should be analyzed separately from the numerical integration.
We are sorry for the confusion. We use “dt” for the integration time step and “δt” for the sampling step, i.e. the time step in the time series which is a multiple of the integration time step. Figure 2B was therefore wrongly labeled and we changed this in the new version of the manuscript. We integrated with small time steps dt, which we previously called the fundamental time steps but are now labeled as integration time steps for clarity. We compared our results for different values of this integration time step dt (Figure 5E of the supplement) and concluded that our chosen time step dt was small enough such that there was no dependence on this step size. The terminology of “delay” and “δ” was introduced in the context of the autocorrelation function which is defined as a function of time delay.
https://doi.org/10.7554/eLife.55650.sa2Article and author information
Author details
Funding
Vrije Universiteit Brussel (SRP31)
 Lana Descheemaeker
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Karoline Faust for interesting discussions, Pankaj Mehta for clarifications about his work on the nicheneutral transition and Stefan Vet and Lendert Gelens for careful reading of the manuscript.
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewing Editor
 Sandeep Krishna, National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India
Publication history
 Received: January 31, 2020
 Accepted: July 20, 2020
 Accepted Manuscript published: July 20, 2020 (version 1)
 Version of Record published: August 6, 2020 (version 2)
Copyright
© 2020, Descheemaeker and de Buyl
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 599
 Page views

 96
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Computational and Systems Biology
Sepsis is not a monolithic disease, but a loose collection of symptoms with diverse outcomes. Thus, stratification and subtyping of sepsis patients is of great importance. We examine the temporal evolution of patient state using our previouslypublished method for computing risk of transition from sepsis into septic shock. Risk trajectories diverge into four clusters following early prediction of septic shock, stratifying by outcome: the highestrisk and lowestrisk groups have a 76.5% and 10.4% prevalence of septic shock, and 43% and 18% mortality, respectively. These clusters differ also in treatments received and median time to shock onset. Analyses reveal the existence of a rapid (30–60 min) transition in risk at the time of threshold crossing. We hypothesize that this transition occurs as a result of the failure of compensatory biological systems to cope with infection, resulting in a bifurcation of low to high risk. Such a collapse, we believe, represents the true onset of septic shock. Thus, this rapid elevation in risk represents a potential new datadriven definition of septic shock.

 Computational and Systems Biology
Predicting gene expression from DNA sequence remains a major goal in the field of gene regulation. A challenge to this goal is the connectivity of the network, whose role in altering gene expression remains unclear. Here, we study a common autoregulatory network motif, the negative singleinput module, to explore the regulatory properties inherited from the motif. Using stochastic simulations and a synthetic biology approach in E. coli, we find that the TF gene and its target genes have inherent asymmetry in regulation, even when their promoters are identical; the TF gene being more repressed than its targets. The magnitude of asymmetry depends on network features such as network size and TFbinding affinities. Intriguingly, asymmetry disappears when the growth rate is too fast or too slow and is most significant for typical growth conditions. These results highlight the importance of accounting for network architecture in quantitative models of gene expression.