Authors: Ashish B. George, James O’Dwyer
Categories: Biological Sciences, Physical Sciences, fluctuations, microbial ecology, science of cities, stochastic population dynamics, urban dynamics, 414, 426
Source: Proceedings of the National Academy of Sciences of the United States of America
Universal abundance fluctuations across microbial communities, tropical forests, and urban populations
Authors: Ashish B. George, James O’Dwyer
Abstract
Large fluctuations in complex populations, such as microbial communities and urban populations can trigger catastrophic ecological and economic collapse. Understanding these fluctuations is crucial for mitigating these rare events. However, models of these systems are often too detailed to fully parameterize and analyze, and not all these details affect the statistical behavior of the population. Here, we analyze time-series data from three disparate complex populations—gut microbial species, employees in US cities, and forest tree species—to demonstrate that emergent behavior in all these populations is described by a single simplified model, resembling the universality observed in many physical systems. Our analysis predicts the distributions of population abundances and fluctuations to improve risk estimation and forecasting of rare events.
The dynamics of complex populations is studied in fields ranging from microbiology to economics. These studies have culminated in theoretical and computational models, with various levels of system-specific detail, that have made progress toward explaining system behavior and developing conceptual understanding. This includes the consumer–resource and Lotka–Volterra models of microbial communities; economic, econometric, and physics-based models of urban dynamics; and ecological models of forests based in niche and neutral theory (1–6). Quantitative analysis, inspired by these models, has also unearthed statistical patterns in the data, hinting at simple emergent behavior of each population (7–10).
What is lacking, however, is an investigation of emergent behavior using a common framework across different populations and an analysis of the temporal fluctuations in these populations. Classic models in statistical physics, such as diffusion and the Ornstein–Uhlenbeck process, have successfully described the behavior of diverse physical systems (11–13). Their success demonstrates that some emergent properties are determined by only a few key underlying details of the dynamics. Efforts at applying this philosophy in other fields have sometimes succeeded (4, 14–16) but are often hindered by the lack of data and the incorporation of many system-specific details that make models difficult to analyze. In this paper, we undertake a statistical physics-inspired investigation of three different complex microbial communities, tropical forests, and human cities, through a common framework. Our analysis will aim to uncover the key underlying similarities and differences between the populations; this will not only deepen our understanding of each system but also facilitate the interconnection of tools and techniques between research fields.
All three complex populations we consider fluctuate in time. Large fluctuations in these populations are associated with catastrophic events such as disease, economic or ecological collapse (17–20); hence, understanding these fluctuations is crucial for risk assessment, quantitative biological methods, and forecasting (20–24). Yet, many models of these populations study the equilibrium and steady-state properties such as the average population abundance (1, 2, 19, 25–28). This restriction is in part because 1) analyzing dynamical properties of complex models is harder than analyzing their steady-state behavior and 2) temporal data required to fit and validate complex models have been lacking. We will tackle these two challenges by 1) adopting a minimal model which makes analysis of fluctuations feasible and 2) using data from three separate systems to fit and validate the model. Generating a reliable null model for population fluctuations that enables improved risk estimation is therefore a major goal of the study.
Fitting models using time series analysis has a precedent in each of the fields we consider here. For example, multiple studies have analyzed the possibility that simple models (making a range of assumptions) can reproduce empirical patterns in time-series data of gut microbiome communities (8, 29, 30), tropical forests (31–34), and urban populations (4, 14) (SI Appendix, section S4). Despite this history of previous analyses, we do not know of any comparison made across these data types using the same model, hence putting all three datasets on the same footing and amenable to direct comparison. We will thus go beyond earlier individual studies by analyzing both time-series and temporal snapshots of data simultaneously, using the same model applied across all three different systems, and introducing ecologically-meaningful categorizations of each set of population data into distinct variant types motivated by domain-specific knowledge.
While these three systems span spatial scales from a single body to an entire country, temporal scales from days to years, and are studied in separate research fields, we harness the similarities in structure of the data to identify the same emergent features in all three systems. We compare the observed features to predictions from a model of stochastic population dynamics with just three population-specific parameters. This simplicity allows us to make analytical predictions for emergent features that we can fit and validate using the limited data. Remarkably, the model is able to capture most of the observed variation, despite its simplicity.
Model and Data
To develop a unified understanding of complex population dynamics across spatio-temporal scales, we analyze time-series population data from three disparate microbial communities in the human gut, trees in a tropical forest, and employment in US cities (Fig. 1A). Traditionally studied by different fields, we interrogate the data using a shared framework by taking advantage of the similarities in the structure of the data. Each dataset contains the relative abundance of population categories within each community at a sequence of time points (Fig. 1B).
More specifically, the gut microbiome dataset records the relative abundances of microbial species in the human gut sampled at daily intervals for almost a year (35). A microbial species or operational taxonomic unit (OTU) was defined based on genomic similarity (35). The data from the Barro Colorado Island (BCI) forest record the number of tree species within a 50-hectare plot on the island, sampled at 5-y intervals for two decades (36). We group trees based on clustering along trait axes into four clusters. There is a long history of grouping species by their maximum height (37), based on the idea that species with access to similar levels, variability, or horizontal uniformity of light are more likely to compete strongly with each other (38). Precisely which species to assign to which height cluster has since been put on a more quantitative footing (39), with four distinct height clusters (shrubs, understory treelets, midstory trees, and canopy trees) identified on BCI. The employment dataset (https://www.bls.gov/cew/) records the number of people employed in different economic sectors, classified according to the North American Industry Classification System (NAICS), in 383 US cities [Metropolitan Statistical Areas (MSAs)]. The data are sampled monthly for 17 y.
We analyze the relative abundances in all three systems. Relative employment in a sector (employment in the sector in a given city, divided by total employment in the city) removes the effect of large variation in total population sizes of cities (5). Similarly, we analyze the relative abundance of a tree species within a height cluster. Using relative abundances in both these systems allows us to treat data from all three systems on the same footing, normalizes out any temporal changes in the total population sizes stemming from overall population growth or decline (SI Appendix, Fig. S10), and eliminates the effects of the large differences in city sizes. SI Appendix, section S1 gives a mathematical definition of relative abundance in each system and SI Appendix, section S12 and Figs. S10–S12 provide our analysis of absolute abundance data. Although working in terms of relative abundances introduces a constraint (that all relative abundances sum to equal one), for diverse communities, this constraint has a relatively insignificant effect. We might also expect departures from our model for very high relative abundances, given that this constraint will tend to change the fluctuation properties for such variants as they approach the size of the entire system. But in practice, this is a small effect; it is rare in these data for any relative abundance to even approach 0.5. For the forest data, we show that our results are robust to the choice of relative vs absolute abundances by analyzing absolute abundances in SI Appendix, section S12 and Figs. S11 and S12.
In contrast to the static predictions from most deterministic models, the data in Fig. 1B show that population abundances continue to fluctuate in time. Furthermore, the strength of these fluctuations differs between the three systems with the largest fluctuations observed in the microbiome. To quantify the strength of fluctuations, we measured the logarithmic fold-change in abundance over a time interval Δt, defined [1]λi=logxi(t+Δt)xi(t),
where xi(t) denotes the (relative) abundance of a species/sector i at time t. The empirical distribution of logarithmic fold-change or the LFD is shown for a single species/sector from each system in Fig. 1C. The comparison with the fit of the logarithmic fold-change by the normal distribution (equivalent to fitting fold-change by the lognormal distribution) illustrates that fluctuations cannot be understood as an outcome of environmental noise without any additional structure or mechanism (SI Appendix). Further, the fit by the normal distribution indicates that fluctuations much larger than expected from a normal distribution may occur in some of these complex populations. These large fluctuations have a major impact on risk estimates and time-series forecasting, and so understanding their distribution is a key step towards advancing these goals.
We now develop a simple model that is capable of predicting statistical features of all three datasets, including the LFD. To achieve our goal of describing all three datasets, we keep the model both as generally applicable and simply formulated as possible. The model assumes that the abundance of a species/sector i, xi fluctuates around an equilibrium value, xi∗, determined by the metabolic, ecological, or economic niche the species/sector occupies. We do not explore the system-specific mechanisms (resource competition, metabolic/economic interactions, etc.) that determine the particular equilibrium value (1, 2, 27, 40–42). Fluctuations then occur due to the stochastic processes governing population growth and decline. Deviations from the equilibrium value result in a linear restoring force, described by (xi∗−xi)/τi, where τi is the timescale of return to equilibrium, where we neglect additional contributions from species interactions. Based on these assumptions, we call the model the stochastic linear-response model (SLRM) of population dynamics.
Assuming that population growth and decline occur in proportion to the abundance, we can write down the stochastic differential equation of the SLRM governing population [2]dxidt=xi∗τi−xiτi+2σixiη(t),
where σi captures the strength of population fluctuations, and η(t) is delta-correlated Gaussian noise or white noise. The model resembles the classic Ornstein–Uhlenbeck process, used to describe many stochastic quantities in physics and finance (11, 13, 43), with one crucial the scaling of the stochastic fluctuations with the square root of xi. This scaling reflects our assumption on the stochastic processes underlying population growth and decline; specifically, the square root scaling arises when growth and decline is a consequence of many independent, random events and the number of such events is proportional to xi, and is commonly used in many stochastic population models (44). We note that this model has been proposed and analytical solutions were obtained previously in the context of singular diffusion processes by Feller (45), bond interest rates by Cox and coauthors (46), and in forest ecology (as a birth–death–immigration model) by Azaele et al. (31). Here, we go beyond these studies by applying a single model across three different systems. The SLRM has not to our knowledge been applied individually to employment data or microbiome data, and our application to forest ecology uses populations divided into niche-based categories.
In ecological terms, the SLRM resembles a scenario where species occupy well-separated niches, with equilibrium abundance xi∗ for species i. This idea of niche-separation determining model parameters is reflected in our assumptions. Specifically, in the microbiome, each species is described by its own SLRM parameters; in urban employment, each sector is described by its own SLRM parameters (independent of city); in forests, each trait-based cluster is described by its own SLRM parameters (independent of species index within the cluster). This assumption is reflected in our choice of showing the empirical LFD for a single sector across all cities, a single microbial species in the gut, and all tree species within a trait cluster in Fig. 1C. While interspecies interactions may occur in each system, incorporating interactions drastically increases the model complexity hindering parameter inference. Further, species interactions may not be necessary to describe population fluctuations in these systems, as we will show.
The minimal nature of model allows us to make analytical predictions for two key quantities that characterize its long-term behavior. First, at long times (t>τi>0), the distribution of population abundances will converge to a steady-state distribution. The steady-state distribution arises from the balance of stochastic fluctuations that kick the population from its deterministic equilibrium value and the restoring force toward the equilibrium. The analytical form of the distribution is (see refs. 31, 45, and 46 and SI Appendix, section S4)[3]Pss,i(x)=(σiτi)−xi∗/σiτiΓ(xi∗/σiτi)xxi∗σiτi−1e−xσiτi,
which is a gamma distribution defined by two σiτi and the equilibrium to noise ratio (ENR), xi∗/σiτi, and Γ() is the gamma function. σiτi quantifies the effective strength of the i.e., the fluctuations over the equilibration time period. The ENR, xi∗/σiτi, quantifies the relative strength of deterministic restoring force to stochastic fluctuations.
Second, we can analytically derive how abundances at steady-state fluctuate in time. Specifically, we quantify the temporal fluctuations by measuring the logarithmic fold-change of population abundances between two time points Δt time units apart, denoted as λi=logxi(t+Δt)xi(t). The distribution of logarithmic fold-changes λi or LFD in steady-state was derived in earlier studies (see refs. 46 and 31 and SI Appendix, section S4)[4]Pfluc,i(λ,Δt)=Γxi∗σiτi+122πΓxi∗σiτieΔtτi−1xi∗σiτieΔt2τicoshλ21eΔtτi(coshλ2)2−1xi∗σiτi+12.
Despite its imposing appearance, the equation depends only on two the ENR, xi∗/σiτi, and the equilibration timescale, τi, which scales the time interval Δt. This distribution quantifies the likelihood of large abundance fluctuations in the model. In particular, the tails of the distribution at large λ decay exponentially at a rate proportional to the ENR, i.e., ∼e−xi∗σiτi|λ|, which ensures the moments of the distribution are well defined. Note that this exponential decay for the logarithmic fold-change at large λ corresponds to a power-law decay for the fold-change, ρ=xi(t+Δt)xi(t), at large ρ with exponent −xi∗σiτi−1.
Results
Stochastic Linear Response Model (SLRM) Reproduces Empirical Abundance Fluctuations across Populations.
In line with our goal of characterizing abundance fluctuations in complex populations, we plot the empirical LFD in the three systems in Fig. 2. Fig. 2A–C plots the LFD of urban employment, aggregated across all US cities, for three different sectors. Fig. 2D–F plots the LFD of three different species in the human gut microbiome. Fig. 2 G and H plots the LFD of all the tree species belonging to two height cluster in the forest. The fold-change in abundance was calculated for time intervals of 1 y, 1 d, and 5 y for the city, microbiome, and forest data, respectively.
All LFDs were roughly centered around λi=0 (average λi was 0.001, 0.004, and 0.05 for city, microbiome, and forest) but exhibited a large variation relative to the mean (coefficient of variation ≫1; see SI Appendix, Fig. S1). This relative dominance of fluctuations over systematic trends for relative abundances not only makes the study of fluctuations easier but also underscores the need to understand fluctuations better. The SLRM, which neglects any systematic trends in the mean, is therefore a good candidate model of population fluctuations.
The red lines in Fig. 2 show the fit of the observed LFDs with the SLRM prediction, Pfluc (Eq. 4). A separate fit is performed for each employment sector, microbial species, and forest niche, allowing the inference of corresponding population parameters. Despite having only two free parameters, the SLRM is able to fit the observed fluctuations across these diverse systems. We compared the SLRM fit to fits by two other candidate distributions using the Akaike Information Criterion (47), as an additional test of the fit. The two candidate distributions (normal distribution and Laplace distribution) were chosen based on their success at fitting empirical fluctuations and alternative models (SI Appendix and ref. 7, 44, and 48). The SLRM prediction outperformed the other candidate distributions in the majority of cases (SI Appendix, Fig. S2 and Datasets S1–S3). We note that our forest analysis demonstrates the importance of the distinct niches, given that inferred timescales range from 400 to 1,700 y, and are significantly different between shrubs and other categories (SI Appendix, section S11 and Fig. S13).
We then further quantified the fits, independent of comparisons with other models, by testing whether the observed data were likely to have been generated by the model. Specifically, for each fit and corresponding parameter estimates, we compared the likelihood of the observed data to the likelihood of 100 random samples of the same size as the data generated from the fitted distribution. This kind of “exact” statistical significance test follows a precedent from testing neutral ecological models and population genetics models (49, 50). If the likelihood of the observed data was smaller than the 95% of the samples, we concluded that the data were unlikely to have been generated by the fitted distribution and rejected the fit. The percentile score in Fig. 1I quantifies the percentage of random samples with a larger likelihood than the observed data. All fits had a percentile score <95%, meaning that SLRM passes our goodness of fit test in all cases. Thus, Pfluc serves as a reasonable null expectation for the distribution of population fluctuations.
Comparing the three systems in Fig. 2, we notice that the scale of fold-change on the x-axis appears larger for the microbiome (Fig. 2D–F) than the two macroscopic systems (Fig. 2A–C, G, and H). These differences in the data are reflected in the parameter values of the fit, and we explore this parameter variation in the next section. Importantly, despite these differences in scale and shape of the data, Pfluc provides a good fit to data from all three systems.
Universal Fluctuations Appear Different in Micro- and Macroscopic Systems due to Different Generation Times.
Fitting Pfluc to the empirical LFD provides us with maximum likelihood estimates for the equilibration timescale τi and the ENR xi∗/σiτi. We compare the dimensionless versions of these inferred parameters across the three systems. Specifically, we compare the ratio of the equilibration timescale to time interval of observation τ/tobs and the ENR across the three systems.
Fig. 3A plots the ENR and τ/tobs across the three systems. While the ENR spans a similar range across all three datasets, the micro- and macroscopic systems differ in the inferred values of τ/tobs. For microbial populations, τ/tobs<1, which means that the duration over which we observe the system, about a year, is longer than its equilibration timescale which is around 10 d. In contrast, cities and forests have τ/tobs>1, implying that the duration of observation, two decades, is much shorter than the equilibration timescale. This difference in τ/tobs is responsible for the difference in shapes of the fitted Pfluc in Fig. 2. Thus, the duration of observation is an important difference between microscopic (microbiome) and macroscopic systems (cities and forests).
An alternative way to compare these populations is to use the ratio of the equilibration timescale to the generation time, τ/tgen (Fig. 3B). The generation times in the three systems were estimated as 4.2 h, 10 y, and 55 y for microbes, employment, and forests respectively (SI Appendix, section S2). Remarkably, when viewed in terms of generation time, all three systems occupy a similar region of parameter space. Hence, fluctuations in the three systems are described by the same distribution over a similar parameter range when time is measured in generations, highlighting the similarities in the emergent behavior across the population types. We note that a small subset of microbes in Fig. 3B appear as outliers, with τ/tgen<1. For these microbes, the empirical LFD is better-fit by a normal distribution than the model prediction Pfluc (SI Appendix, Fig. S3). For the remaining data, the equilibration time scale is on the order of 10 generations for all three systems. Therefore, when viewed in terms of generation time rather than physical time, emergent fluctuations in the three systems are highly similar.
In addition to examining how the inferred timescale, τ, varied across the three systems, we also examined how τ varied within each system. Analyzing employment data, we found that the most abundant sectors in cities such as healthcare and retail trade had the longest timescales while the least abundant sectors in cities such as agriculture and mining had the shortest timescales. Quantitatively, we found that the median relative abundance of a sector across cities was correlated with the inferred timescale (Pearson’s r =0.86, Spearman r =0.66, P<0.05). We found a similar relationship for microbes, with the most abundant species (belonging to Bacteroides) having the longest timescales (Pearson’s r =0.93, Spearman r =0.68, P<0.05) (SI Appendix, Fig. S5). For forest clusters, we found that shrubs had a significantly shorter timescale than the other three, taller height clusters (SI Appendix, section S11 and Fig. S13). We will return to interpret these correlations shortly.
SLRM Reproduces Empirical Distribution of Species Abundances.
From Fig. 3, we see that the observation duration for the microbiome is longer than its equilibration time scale (τ/tobs<1). Hence, the temporal trajectory of abundances of each microbial species should trace out the corresponding steady-state distribution, which is a gamma distribution described by two parameters Eq. 3. We plot the distribution of relative abundance of microbial species and fit it with the two-parameter gamma distribution (black line), via maximum likelihood in Fig. 4A–C.
For city and forest data, the observation timescale is longer than the equilibration timescale (τ/tobs>1), and so the temporal trajectory of abundances will not converge to the steady-state distribution. Instead, we plot the cross-sectional distribution of abundances, i.e., the relative abundance of a sector in a city aggregated across all cities and the abundances of tree species within a height cluster, at a random time point (Fig. 4D–H). The fit with the gamma distribution is the black line.
To test whether the observed data could have been generated by the model, we repeat the procedure used in Fig. 2. Specifically, we compared the likelihood of the observed data to the likelihood of randomly generated samples from the fitted distribution. All datasets had a likelihood comparable to a random sample, as evidenced by the percentile scores (from 1,000 random samples of the same size) shown in Fig. 4I. Two additional goodness of fit tests were also performed and the majority of species passed both tests (SI Appendix, section S9 and Figs. S6 and S7). Together, these observations indicate that the abundances in these systems are well described by the gamma distribution. This conclusion is further supported by previous research that used a gamma distribution to successfully fit cross-sectional microbial abundances across microbiomes (8).
The steady-state distribution (SSD) and LFD predicted by the model (Eqs. 3 and 4) share a parameter, the ENR (xi∗/σiτi). We compare the maximum likelihood estimates of the ENR obtained by the separate fits of the empirical abundance distribution and LFD in Fig. 5. The discrepancy in inferred ENR values could arise from limitations of the data or the model.
To check whether the model is able to simultaneously fit both the LFD and abundance distribution, we performed a modified version of our goodness of fit test we used previously. First, using the ENR inferred from fitting the abundance, we compare the likelihood of fitting the LFD by Pfluc to the likelihood of 100 random samples from the fitted distribution of the same size. Then, we fit the other way around, i.e., using the ENR inferred from fitting the LFD, we compare the likelihood fitting the abundance by Pss to the likelihood of 1,000 random samples from the fitted distribution of the same size. If the likelihood of the observed data was within 95% of the likelihood of the samples for either of these comparisons, then we conclude that the inferred parameters were consistent, i.e., the ENR obtained from fitting one distribution is able to provide a reasonable fit of the other distribution. The data points deemed consistent are depicted by filled markers in Fig. 5. If the data likelihood was smaller than 95% in both cases, we term the inferred parameters as inconsistent. The data points deemed inconsistent are shown as unfilled markers in Fig. 5.
The ENR estimates for the majority of the data were consistent (Fig. 5). A minority of the microbes (28/85) were rejected as being inconsistent. This could partly be due to temporal correlations in the data used to fit the abundance distribution, which we neglected. These correlations exist over timescales of τi and only vanish when Δt≫τi. For cities and forests, we used cross-sectional data, which do not suffer from this drawback. All employment data, including the apparently large outliers in Fig. 5, were consistent. The consistency of these large outliers was because although the abundance distribution of some sectors was well fit by large ENR values, substantially smaller ENR values also provided a reasonable fit and could not be rejected using the few hundred observations. Overall, the majority of the observed species/sectors in the three systems followed the expected relationship between the predicted distributions for abundances and fluctuations.
Variation of Model Parameters Explains Taylor’s Law.
In addition to providing a simple 2-parameter null model for fluctuations and abundances in complex populations, the SLRM can also help understand other empirical patterns in the data. Prior research in microbiome data has revealed an approximate power-law scaling (with exponent 2) of the variance vs. mean of the abundance, called Taylor’s law (7, 8) (Fig. 6), which may arise due to various mechanisms (51, 52). Examination of inferred model parameters provides an alternative explanation for this empirical observation.
To investigate why Taylor’s law arises, we compare the variation of the inferred model parameters across microbial species, using the coefficient of variation (Fig. 6A). The ENR (xi∗/σiτi) has a substantially lower coefficient of variation than the other parameters, and so can be considered to be approximately constant. For the gamma distribution, the mean abundance of a species i is xi∗, and the ratio of variance to mean squared is the ENR, xi∗/σiτi. In microbiome data, the ENR is constant, and so the ratio of variance to mean squared remains fixed while the mean abundance varies, which leads to the observation of an approximate power-law scaling of the variance with mean. Hence, the approximate constancy of the ENR provides an alternative explanation for the observed power-law known as Taylor’s law in microbiome data. The approximate constancy of ENR implies that the distribution of fluctuations (LFD) of different microbial species will be similar when time intervals are measured in terms of τi. Furthermore, the approximate constancy may explain one of our earlier findings—that relative abundance correlates with inferred timescale (SI Appendix, Fig. S5). Mathematically, since the ENR does not vary significantly, relative abundance in our model must be proportional to σiτi. Future work may uncover the mechanisms behind why ENR is approximately constant across these systems.
Comparing SLRM and a Model with Environmental Noise.
The SLRM incorporates “square-root” fluctuations, referred to as demographic noise, and commonly used in many population dynamics models (44). Demographic noise captures the fluctuations arising from accumulation of small, independent random growth and death events, and mathematically, it can be identified by the square-root scaling of the noise term with population size (Eq. 2). An alternative form of noise used in population dynamics models is environmental noise (44). Environmental noise aims to capture fluctuations arising from random fluctuations of the overall growth and death rates of the population and has been used in the analysis of both microbiome data (8) and local forest communities (32–34). The latter involve a range of choices of model specification, including the way competition is imposed in the local community and how dispersal is modeled from regional pool to local patches. This complexity tends to yield models without analytical solutions for LFD and SSD. On the other hand, a relatively simple implementation of environmental stochasticity is the stochastic logistic model (SLM), which has been applied to recapitulate the observed abundance distributions in another of our three data microbiome communities (8, 29). All of these models are characterized by the same linear scaling of the noise term in the population size (Eq. 5), and so comparing our model with an environmental noise model provides an initial test of whether environmental stochasticity will inevitably tend to provide a better description of fluctuations than demographic noise alone. Therefore, in this section, we compare the SLRM, which incorporates demographic noise, with the SLM.
The SLM is also a three-parameter model defined by the following equation for the relative species abundance xi:[5]dxidt=xiτi′1−xiKi+Siτi′xiηi(t),
with parameters Ki, which describes the carrying capacity of the population, Si, which captures the strength of fluctuations, and τi′, which sets the timescale of growth and equilibration. ηi is delta-correlated Gaussian noise or white noise.
The steady-state abundance distribution (SSD) of the SLM is also a gamma distribution (8), like the SLRM:[6]Pss,i(SLM)(x)=x2Si−1−2Γ(2Si−1−1)2KiSi2Si−1−1e−2xKiSi.
It is parameterized by combinations of the two parameters Ki and Si. However, unlike the SLRM, there is no analytical prediction for the LFD of the SLM.
Therefore, to facilitate a direct comparison between the SLRM and SLM, we adopted the following first, we fixed two of the three parameters in both models by fitting the gamma-distributed SSD predicted by the models to the observed abundance distribution. Then, we simulated the SLM for a range of values of remaining parameter, τ′, to obtain a series of predicted LFDs from SLM simulations. We obtained the LFD for the same range of τ values of the SLRM through the analytical predictions. Finally, we compared the disagreement between the two sets of predicted LFDs and the empirical LFDs in the three systems by computing the Jensen–Shannon distance (JSD) between them (see SI Appendix, section S3 for further details).
In Fig. 7 A and B, we illustrate this procedure applied to data on the employment in the management sector in US cities. The distributions predicted by the two models as the timescale parameters (τ′, τ) are varied are shown by the colored lines alongside the observed data (black circles). The disagreement between the model prediction and the observed data is measured using the JSD and shown in the insets. The model with a lower JSD better explains the observed data. We repeated this analysis for all microbial species, employment sectors, and forest niches, and found that the SLRM provided a better fit in the majority of the cases (Fig. 7C). Specifically, the SLRM had a lower JSD in 17 out of the 18 employment sectors, 72 out of 85 microbial species, and 2 out of the 4 forest clusters. Interestingly, it is the forest data where a body of evidence exists to demonstrate the importance of environmental noise in explaining other aspects of population fluctuations (34, 53), where the SLRM and SLM were most similar in their performance. But outside of these cases, the SLRM provides the better description of the data.
Discussion
Understanding the behavior of diverse complex populations is challenging, but critical for progress in many fields. Here, we analyze all three of microbiome data, forest data, and employment data on an equal footing. This unified approach to studying the three populations complements detailed investigations of each specific population. In the microbiome context, where many studies focus on species interactions (54, 55), our work demonstrates how statistical features of the data can arise from stochastic fluctuation, with minimal interactions. For forest data, our work finds some independent support for prior analysis that groups trees into distinct height niches (39). And for employment data, where many studies have focused on detailed econometric models and analyses (6, 40, 56), our work provides a unique perspective from the point of view of a simple dynamical model. Moreover, by showing that the emergent distributions of population abundances and fluctuations in all three populations can be reproduced by a single model, our analysis highlights emergent properties that are independent of system-specific details. When using generations to measure time, all three systems occupy a similar range of parameters, suggesting that fluctuations in cities and forests over decades and centuries would closely resemble microbiome fluctuations. In broader terms, the fact that our analysis provides a good description across a variety of systems may point to the existence of universal properties in the fluctuations of complex systems.
In quantitative terms, our analysis provides a simple two-parameter functional form for the distribution of fold-changes in the populations and relates this to the abundance distribution. The predicted distribution of fluctuations Pfluc (Eq. 4) is able to fit fluctuations in all three systems despite having only two parameters, and moreover, we find that when measured in terms of generation time, parameter fits for all three systems collapse into a narrow window of fitted values. Further, Pfluc arises from a well-defined model and suggests a plausible mechanism; it is not simply chosen from the vast library of statistical distributions historically examined. The observation that we find all three data types are described by similar parameter values suggests important, deeper connections that may yet be uncovered in future work.
The predicted distribution of fluctuations Pfluc can serve as an important null model of the fluctuations in complex populations, where understanding the likelihood of large fluctuations is crucial. For employment fluctuations, large fluctuations impact urban planning and economic stability; for forests, large fluctuations impact ecological management strategies; for microbiome, large fluctuations can cause dysbiosis, which affects the health of the host (17, 19, 20). A two-parameter null model for fluctuations in complex populations is useful in practical, data-limited settings; it can estimate the risk of large fluctuations more accurately (SI Appendix, Fig. S4) and improve quantitative methods, such as those utilizing Bayesian inference from time-series data to classify ecosystem states (57) and priors for priors for decision-making and modeling (58–60). In SI Appendix, we show how Pfluc fits empirical data from the three systems when fluctuations are measured over different time intervals (SI Appendix, section S6 and Figs. S8 and S9).
The SLRM can be understood as the linearization of a more complex nonlinear model around its equilibrium when species interactions are neglected. This is shown in our SI Appendix, section S5, where we also demonstrate an example where a model with interspecies interactions actually reduces exactly to the SLRM. In general though, this simplification drastically reduces the number of parameters, making parameter inference from available data feasible. Deviations from the model predictions, could indicate the presence of species interactions, which are often modeled by Lotka–Volterra, consumer–resource, and other models of higher order species interactions (1, 2, 54), autocorrelated noise, or other mechanisms. Such models could potentially be parameterized by using specialized methods with additional data (22, 61). In SI Appendix, we discuss how the SLRM can be extended into a stochastic model incorporating species interactions in a linear regime. Analysis of this extended model could pave the way for novel inference methods that account for the stochastic fluctuations in observational data.
While we compared our model predictions with a range of classic distributions, we also note that environmental noise (44) has been proposed as an explanation for fluctuations in abundance across different complex systems, including multiple forest datasets (32–34) and microbiome data (8). To capture this mechanism and compare our model to its predictions, we tested the performance of the SLRM to the SLM (8, 29), a three-parameter model that combines nonlinear logistic growth with environmental stochasticity. The SLM lacks an analytical solution for the LFD, but through numerical simulations, we compared the fits of the SLRM and SLM to the empirical data. We found that the SLRM outperformed the model with environmental noise in the majority of our data (Fig. 7). While there are multiple other types of environmental noise model, for example, those that have provided a good description of local forest community fluctuations (33, 34), this comparison demonstrates that environmental stochasticity does not necessarily provide a better description of fluctuations in complex populations. More general models of environmental stochasticity tend to lack simultaneous analytical solutions for the LFDs and SSD, making numerical comparison more challenging. However, it is certainly possible that, just as with species interactions, more general kinds of noise should form part of the basis for extending our model, and future analysis will likely shed light on this question.
Framing the SLRM as a useful base model for further research, we note that it can be easily augmented with additional mechanisms, including environmental fluctuations and species interactions, which could be tested with additional data (SI Appendix). Other modifications could help understand evolving populations. For the timescales examined, we assumed that the equilibrium value x∗ remains constant. Over longer timescales, however, the equilibrium value could change due to biological evolution, climate change, or socio-technological revolution. Investigating the model when x∗ changes in time could help understand emergent dynamics in complex populations over evolutionary timescales and presents an interesting direction for future research.
To butcher two well-worn phrases, all models are wrong, some are useful, and some are unreasonably effective. We believe the SLRM falls into the latter two categories and that its surprising effectiveness across such diverse datasets points to something universal about the way complex populations fluctuate. The SLRM also provides valuable two parameter null models for the distributions of abundances and fluctuations in complex populations, which are of particularly utility in data-limited scenarios for forecasting and risk analysis. The unified analysis of the three population types highlights both similarities and differences between the systems and paves the way for a fruitful exchange of tools, techniques, and interpretations between these very different fields.
Materials and Methods
Data Processing.
City data.
Public domain employment data from 2003 to 2019 were obtained from Quarterly Census of Employment and Wages (https://www.bls.gov/cew/). The data list employment classified into industrial sectors by the NAICS at the county level in the United States (SI Appendix, Table S1). We aggregated data at the county level to 383 MSAs, which we call cities. MSAs are independent statistical units defined by the Census Bureau that encompass a central city and surrounding connected areas. We used the list of counties in each MSA in 2017 from US Census Bureau, County Business Patterns program, to maintain a consistent definition of MSAs across all time points. Although employment data are recorded at monthly intervals, we used a Δt of 1 y for calculating the empirical LFD since many industries, such as agriculture and accommodation, display seasonal employment trends within a year. The duration of observation, tobs, was 17 y. Relative sectoral employment was calculated by dividing sectoral employment in a city with the total employment in the city. We plot and fit data with only non-zero abundance. See SI Appendix, section S1 for further details.
Microbiome data.
Microbiome data from ref. 35 was obtained and processed as in ref. 7. We consider only the gut microbiome data for individual M3 since it was substantially longer than other time series. The data were collected in time intervals of one day with some sampling gaps for 336 time-points; hence, Δt=1 d for microbiome data. Due to sampling gaps, we used an approximate duration of observation, tobs, of 300 d. To process the data, first, read counts at each time point were normalized to obtain relative abundance. Only prevalent (present in more than half of the time points) or abundant species (average relative abundance >10−3) were used. A total of 85 species met these criteria. We plot and fit data with non-zero relative abundance. Each species has its own set of SLRM parameters.
Forest data.
The BCI forest data was obtained from the Center for Tropical Forest Science website (https://forestgeo.si.edu) (36). Abundance data collected at 5-y intervals, in years 1990, 1995, 2000, and 2005, were used. Only trees that were alive and had a diameter at breast height >10 cm were counted. The trees were grouped into four height shrubs, understory treelets, midstory trees, and canopy trees based on ref. 39. There were 87, 75, 60, and 63 species in the four clusters. Relative abundance of a species in a cluster was calculated as the absolute abundance of the trees of the species at that time point divided by the absolute abundance of all tree trees within the cluster. We assume that all species within a height cluster are highly similar in trait values and can be described by the same set of SLRM parameters. Each height cluster is fit separately, like the different employment sectors. Time interval for calculating LFD, Δt, was 5 y. The duration of observation, tobs, was 20 y.
Fitting and Sampling Procedures.
The data were fit and parameters were estimated by maximum likelihood estimation from the Scipy package in Python. In addition to the probability distribution Pfluc, the corresponding cumulative distribution (SI Appendix, section S4) was defined to make sampling from the distribution more efficient.
In the city data, we had substantially more data points in each LFD (∼30,000 points) than for the empirical abundance distribution (∼300 points). This made performing the goodness of fit test computationally harder for the LFD than the abundance distribution, since the test required sampling from the distribution of the same size as the data. Hence, for each fit of the abundance distribution, we obtained 1,000 samples to compare with the data likelihood, whereas for each fit of the LFD, we obtained only 100 samples.