Emil Mallmin<fnref n="1" />*, Arne Traulsen<fnref n="1" />, and Silvia De Monte<fnref n="1" />,<fnref n="2" />

Revised version, February 2024

**Abstract:** The composition of ecological communities varies not only between different locations but also in time. Understanding the fundamental processes that drive species towards rarity or abundance is crucial to assessing ecosystem resilience and adaptation to changing environmental conditions. In plankton communities in particular, large temporal fluctuations in species abundances have been associated with chaotic dynamics. On the other hand, microbial diversity is overwhelmingly sustained by a 'rare biosphere' of species with very low abundances. We consider here the possibility that interactions within a species-rich community can relate both phenomena. We use a Lotka-Volterra model with weak immigration and strong, disordered, and mostly competitive interactions between hundreds of species to bridge single-species temporal fluctuations and abundance distribution patterns. We highlight a generic chaotic regime where a few species at a time achieve dominance, but are continuously overturned by the invasion of formerly rare species. We derive a focal-species model that captures the intermittent boom-and-bust dynamics that every species undergoes. Although species cannot be treated as effectively uncorrelated in their abundances, the community's effect on a focal species can nonetheless be described by a time-correlated noise characterized by a few effective parameters that can be estimated from time series. The model predicts a non-unitary exponent of the power-law abundance decay, which varies weakly with ecological parameters, consistent with observation in marine protist communities. The chaotic turnover regime is thus poised to capture relevant ecological features of species-rich microbial communities.

**Keywords:** ecological chaos | rare biosphere | plankton ecology | disordered Lotka-Volterra model

# 1 Introduction

The dynamic nature of ecological communities of species has long been recognized [1]. Fluctuations in species' abundances can have fundamentally different origins depending on the spatial and temporal scales considered, and the particular community of interest [1, 2]. For example, if environmental filtering shapes community composition, abundance fluctuations may reflect changing external conditions. Another possibility is that stochastic effects of demography, dispersal, and diversification dominate over the ecological differences between species in driving turnover [3]. As communities are increasingly monitored in the wake of

climate change and biodiversity decline, there is growing opportunity and need to understand why abundances fluctuate, and how fluctuations relate to patterns of biodiversity and biogeography.

An alternative to environmental and stochastic effects as the main determinants of abundance fluctuations is the hypothesis that they reflect intrinsically chaotic dynamics arising from the complexity of ecological interactions. Mathematically, chaos is the phenomenon whereby a deterministic, nonlinear dynamical system (e.g. describing the populations of interacting species) generates bounded but apeironic trajectories that depend sensitively on initial conditions [4]. While chaos can be readily identified in simple mathematical models, its presence in empirical time series is challenging to ascertain, and the relevance of chaos for natural communities has been controversial [5, 6]. However, recent methodological advances and systematic assessment of a large ecological time series database using validated, non-parametric methods showed that ecological chaos is generally not rare [7, 8], and is particularly prevalent in planktonic communities, where it was found in ~80% of time series.

The biodiversity of microbial communities such as plankton is overwhelmingly sustained by the 'rare biosphere' revealed by recent methods of high-throughput genomic sequencing [9–11]—an extreme instance of a near-universal observation that ecological communities harbour a few highly abundant, dominant species, and a much larger number of low-abundance, rare species. In plankton profit communities sampled from multiple distant locations in the world oceans, the number of rare species increases as a power-law as lower abundances are considered; a pattern that is quantitatively uniform across samples despite their strong compositional differences [12]. In addition to spatial variations, strong temporal turnover has been observed for plankton, where species abundances can change dramatically on a short time scale, even when abiotic conditions do not vary substantially [13–15]. A role for intrinsic ecological dynamics in driving such complex oscillations is supported by mesocosm experiments, where sustained abundance fluctuations have been observed even under stable external conditions, both for plankton [16, 17] and other microbes [18, 19].

The conditions enabling ecological chaos can be investigated with mathematical models. Traditionally, models of population dynamics have considered only a handful of taxa. There, chaos tends to occur only within a narrow parameter range [20, 21]. In contrast, high-dimensional dynamical systems (involving dozens or hundreds of interacting degrees of freedom) seem to display chaos more generically [22]. Robust fluctuating states (variably iden-

arXiv:2306.11031v2 [q-bio.PE] 21 Nov 2024

<footnote n="1">¹Max Planck Institute for Evolutionary Biology, Plön, Germany</footnote>

<footnote n="2">²Institut de Biologie de l'ENS (IBENS), Département de Biologie, Ecole Normale Supérieure,
CNRS, INSERM, Université PSL, 75005 Paris, France</footnote>

<footnote>*mallmin@evolbio.mpg.de</footnote>
tified as chaos) were found in models of species-rich communities with competitive [23–25], predator-prey [26, 27], or consumer-resource [28] interactions. Some studies reported that chaotic regimes tended to sustain a higher diversity than equilibria due to the availability of more (spatial) temporal niches [20, 24, 27, 29]. Nonetheless, stabilizing mechanisms are required to prevent large abundance fluctuations from causing diversity-limiting extinctions. A metacommunity structure (a network of patches connected through dispersal) offers one plausible solution. Under conditions where patches' abundance dynamics do not synchronize, local extinction can be compensated by migration from another patch, and the fluctuations persist on time scales much longer than the local dynamics [24–26, 30].

Here, we are instead interested in characterizing the within-patch chaotic dynamics in order to relate two complementary perspectives: that of the fluctuating abundance time series of individual species, and that of local community-level statistics—such as the instantaneous distribution of abundances across species and the overall strength of interactions. We consider a general model for large communities where strong ecological interactions encompass, as in microbial ones, vigorous competition between the composing species, but also facilitation. We simplify spatially structured models by considering a local, well-mixed community with a constant, small immigration. Like a 'seed bank' [31] or the effect of metacommunity dispersal, this prevents irreversible loss of species. Following a well-established 'disordered' approach to complex communities [32–35], we consider pairwise interactions drawn from a random distribution, but focus on the little studied regime where only a handful of species can dominate the community at any time, while most other species are rare, consistent with the rare biosphere pattern.

We show that in this general setting a broad range of model parameters allows species to alternate chaotically between rarity and abundance on a characteristic timescale such that the community composition moves through a succession of low-diversity states. The distribution of abundances attained by any given species over a long time series largely overlaps with the distribution of abundances found in the whole community at any given time, which is a power-law across many orders of magnitude in abundance values. This correspondence suggests an equivalence among different species despite their clear ecological differences and short-term competitive exclusion dynamics. We derive a stochastic focal-species model that captures, in a statistical sense, the dynamical features common to all species, and also identify the origin of species-specific deviations in the propensity to dominate the community.

## 2 Model

We describe a community of S species by their time-dependent absolute abundances $x_i(t)$, with $i = 1, 2, ..., S$ the index of a species. Microbial communities have been described by deterministic equations where changes in abundance relate to competition within species and pairwise species interactions [19, 36, 37]. According to the Lotka-

Volterra equations [38], the abundance of any species in isolation grows logistically: if initially the species is rare, its abundance grows exponentially at a maximum rate $r$, doubling every $(\ln 2)/r$ time units. Eventually, it saturates to a carrying capacity $K$ set by resources, predation, and abiotic conditions, assumed constant and not modelled explicitly. For simplicity, we set $r$ and $K$ to unity for all species, but discuss heterogeneity in these parameters in Supplementary Note S3. The interaction coefficients $\alpha_{ij}$ (real numbers) quantify the effect of species $j$ on the growth rate of species $i$; by convention, detrimental when $\alpha_{ij} > 0$, and facilitative when $\alpha_{ij} < 0$. We include a small rate of immigration $\lambda \ll 1$ into the community, constant and equal for each species, to set a lowest level of rarity and prevent extinctions. Abundances thus change in time as

$$ \dot{x}_i(t) = x_i(t) \left( 1 - x_i(t) - \sum_{j=1(\ne i)}^{S} \alpha_{ij} x_j(t) \right) + \lambda. \quad (1) $$

In species-rich communities, the number of potential interactions—$S \times S$—is very large, and their values hard to estimate in natural settings. A classic approach is therefore to model the set of interaction coefficients as a realization of a random interaction matrix $A$ [32–35]. When $S$ is large, patterns of ecological interest are expected to depend on the summary statistics of $A$ rather than its particular realization. We consider for simplicity Gaussian statistics $A_{ij} \sim N(\mu, \sigma^2)$ ($i \ne j$). A correlation $\gamma$ between diagonally opposed elements can be introduced, biasing interactions toward predator-prey ($\gamma = -1$) or symmetric competition ($\gamma = 1$); here, we focus on independent interaction coefficients ($\gamma = 0$) and discuss other cases in Supplementary Figure S5.

The interaction coefficients for distinct species $i, j$ can be represented in terms of the mean $\mu$ and standard deviation $\sigma$ of the interaction matrix, as

$$ \alpha_{ij} = \mu + \sigma z_{ij}, \quad (2) $$

where the $z_{ij}$ are realizations of random variables with zero mean and unit variance. We note that, by convention, we have separated the self-interaction term from the intra-specific interaction terms in Eq. (1). The diagonal element $\alpha_{ii}$ therefore does not appear in the sum, and is not defined.

Equation (1) with randomly sampled interactions defines the disordered Lotka-Volterra (dLV) model. By tuning the ecological parameters $S, \mu, \sigma, \lambda$, it exhibits a number of distinct dynamical behaviours which have been thoroughly explored in the weak-interaction regime, where the interaction between any particular pair of species is negligible, but a species' net competition term from all other species is comparable to its (unitary) self-interaction. If species are near their carrying capacities, the net competition is approximately

$$ \sum_{j(\ne i)} \alpha_{ij} = S\mu + \sqrt{S\sigma} z_i, \quad (3) $$

where the net interaction bias

$$ z_i := \frac{1}{\sqrt{S}} \sum_{j(\ne i)} z_{ij} \quad (4) $$
is a realization of a random variable ~ N(0, 1). A finite net competition in the limit of a large species pool requires

$$ \mu = \frac{\tilde{\mu}}{S}, \quad \sigma^2 = \frac{\tilde{\sigma}^2}{S}, \qquad (5) $$

where $\tilde{\mu}$, $\tilde{\sigma}$ do not grow with S. Under this scaling, methods from statistical physics (dynamical mean-field theory [34, 39–41], random matrix theory [32, 42], and replica theory [43, 44]) allow exact analytical results in the limit of $S \to \infty$, although in practice $S \sim 100$ is sufficient for good agreement between theory and simulations. Sharp boundaries were shown to separate a region where species coexist at a unique equilibrium and one with multiple attractors, including chaotic steady-states [34, 39–41].

Since we are here interested in the scenario of large differences in species abundance (rare biosphere pattern) and rapid turnover dynamics, we instead consider the strong-interaction regime where the statistics of the interaction matrix do not scale with species richness S according to Eq. (5). For $S\mu \gg 1$, the overall competitive pressure makes it impossible for all species to simultaneously attain abundances close to their carrying capacities. Abundant species tend to exclude one another, resulting in instability and complex community dynamics. Arguably, strong interactions are more plausible than weak ones for microbial communities, where metabolic cross-feeding, toxin release, phagotrophy, and competition over limited nutrients lead species to depend substantially on one another’s presence [45, 46].

## 3 Results

In the strong-interaction regime, numerical simulations of the disordered Lotka-Volterra model show that the community can display several different classes of dynamics, from equilibrium coexistence of a small subset of species, to different kinds of oscillations, including chaos. In sections 3.1–3.4, we focus on reference value of the interaction statistics ($\mu = 0.5, \sigma = 0.3$) representative of chaotic dynamics, and describe its salient features. In Section 3.5 and 3.6, we describe how the dynamics depends qualitatively on the statistical parameters $\mu$ and $\sigma$. Unless otherwise stated, simulations use $S = 500$ and $\lambda = 10^{-8}$. Further details on the numerical implementation are presented in Appendix A.

### 3.1 A chaotic turnover of rare and abundant species

For a broad range of parameters in the strong-interaction regime, the community undergoes a chaotic turnover of dominant species. As illustrated by the time series of stacked abundances in Figure 1A, the overwhelming share of the total abundance at any given time is due to just a few species. Which species are abundant and which are rare changes on a characteristic timescale, $\tau_{dom} \approx 30$ time units, comparable to the time it would take an isolated species to attain an abundance on the order of its carrying capacity starting from the lowest abundance set by immigration. While the total abundance fluctuates moderately around a well-defined time average, individual species follow a ‘boom-bust’ dynamics. If this simulation represented a natural microbial community, only the most abundant species—that

<figure><img src="image_1.png" /><figcaption>**Figure 1: Turnover of the dominant component.** A The stacked abundances of all species under steady-state conditions: there is a turnover of species such that only the dominant component is visible at any given time (each species has a distinct random colour). B Bray-Curtis index of community composition similarity between the dominant component of the community at time t, and the composition if it were isolated from the rare species and allowed to reach equilibrium: the community appears to approach the composition of few-species equilibria before being destabilized by invasion from the pool of rare species.</figcaption></figure>

we call the dominant component of the community—would be detectable by morphological inspection or shallow sequencing.

We wish to characterize the dominant component, and understand how it relates to the pool of rarer species. In order to quantify the notion of dominance, we define the effective size of the community as Simpson’s (reciprocal) diversity index [47],

$$ S_{\text{eff}}(t) := \frac{1}{\sum_{i} p_{i}^{2}(t)}, \qquad (6) $$

where $p_i = x_i / \sum_j x_j$ denote relative abundances. $S_{eff}$ approaches its lowest possible value of 1 when a single species is responsible for most of the total abundance, and its maximum S when all species have similar abundances. Its integer approximation provides the richness, i.e. number of distinct species, of the dominant component.

The effective size $S_{eff}$ of the community in our reference simulation fluctuates around an average of 9 dominant species, which make up 90% of the total abundance. The relative abundance threshold for a species to be in the dominant component fluctuates around 3%, which is comparable to the arbitrary 1%-threshold used in empirical studies [48]. In Supplementary Figure S4 we show that the number of dominant species grows slowly (but super-logarithmically) with S, up to about 15 for $S = 10^4$. Thus, strong interactions limit the size of the dominant component, and the vast majority of species are rare at any point in time.

The turnover of dominant species is not periodic; indeed, even over a large time-window, where every species is found on multiple occasions to be part of the dominant component, its composition never closely repeats (Supplementary Figure S3). This aperiodicity suggests the presence of chaotic dynamics. We give numerical evidence for sensitive dependence on initial condition and positive maxi-
mal Lyapunov exponent in Supplementary Figure S1 and S2. The turnover dynamics has the character of moving, chaotically, between different quasi-equilibria corresponding to different compositions of the dominant community (c.f. 'chaotic itinerancy' [49]). To reveal this pattern, we measure a 'closeness-to-equilibrium', defined as the similarity in composition between the observed dominant component at a given time, and the equilibrium that this dominant component would converge to if it were isolated from the rare component and allowed to equilibrate. As a similarity metric we use the classical Bray-Curtis index (Appendix B), which has also been used to measure variations in community composition in plankton time series [13]. In Figure 1B, we see that the similarity at times slowly approaches 100%, followed by faster drops, towards about 50%, indicating the subversion of a coherent dominant community by previously rare invaders.

The fact that the community composition is not observed to closely repeat is arguably due to the vast number of possible quasi-equilibria that the chaotic dynamics can explore. In the weak-interaction regime, a number of unstable equilibria exponential in $S$ has been confirmed [50, 51]. It is therefore conceivable that the number of quasi-equilibria in our case is also exponentially large. The LV equations for $\lambda = 0$ admit up to one coexistence fixed point (not necessarily stable) for every chosen subset of species [38]. Hence, we expect on the order of $\sim S^{S_{eff}}$ quasi-equilibria, which for $S = 500$ and $S_{eff} \approx 9$ evaluates to $10^{24}$. If the dynamics explores the astronomical diversity of such equilibria on trajectories which depend sensitively on the initial conditions, the dominant component may look as if having been assembled 'by chance' at different points in time.

The composition of the dominant community is not entirely arbitrary, though. While the abundance time series of most pairs of species have negligible correlations, every species tends to have a few other species with a moderate degree of correlation. In particular, if $(\alpha_{ij} + \alpha_{ji})/2$ is significantly smaller than the expectation $\mu$, and hence species i and j are close to a commensal or mutualistic relationship, these species tend to 'boom' one after the other (Supplementary Figure S6).

## 3.2 Species' abundance fluctuations follow a power-law

In a common representation of empirical observations, where relative abundances are ranked in descending order (a *rank-abundance plot* [52]), microbial communities display an overwhelming majority of low-abundance species [10]. Our simulated community reproduces this feature; Figure 2A. The exact shape of the plot changes in time, as does the rank of any particular species, but the overall statistical structure of the community is highly conserved. An alternative way to display the same data is to bin abundances and count the frequency of species occurring within each bin, producing a *species abundance distribution* (SAD) [52]. The histogram in Figure 2B illustrates the 'snapshot' SAD for the rank-abundance plot in Figure 2A of abundances sampled at a single time point. Whenever observations are available for multiple time points, it is also possible to plot, for a given species, the histogram of its abundance in

<figure><img src="image_2.png" /><figcaption>**Figure 2: Statistical features of abundance variations across species and in time.** A Snapshot rank-abundance plot for the relative abundances in the reference simulation: most species have orders of magnitude smaller abundances than the top ranks. Different lines represent observations at well-separated time points. B Species abundance distribution (SAD, blue histogram) corresponding to the blue rank-abundance plot; overlaid, abundance fluctuation distribution (AFD), averaged over all species (black line) with $\pm$ one standard deviation across species shaded in grey: the snapshot SAD appears to be a subsampling of the average AFD, indicating an equivalence, but de-synchronization, of species in their abundance fluctuations. The one bar missing from the SAD is the effect of finite species richness, as high-abundance bins only ever contain a couple of species for $S = 500$. The vertical dashed line indicates the immigration level which determines a lower limit to abundances.</figcaption></figure>

time. As time gets large (practically, we considered 100'000 time units after the transient), the histogram converges to a smooth distribution, that we call the *abundance fluctuation distribution* (AFD) [53]. Its average shape across all species is also displayed in Figure 2B.

Several conclusions can be drawn by comparing SADs and AFDs. First, a snapshot SAD appears to be a subsampling of the average AFD. Therefore, SADs maintain the same statistical structure despite the continuous displacement of single species from one bin to another. Second, *ev ery* species fluctuates in time between extreme rarity ($x \approx \lambda = 10^{-8}$) and high abundance ($x \gtrsim 10^{-1}$). This variation is comparable to that observed, at any given time, between the most abundant and the rarest species. Third, species are largely equivalent with respect to the spectrum of fluctuations in time, as indicated by the small variation in AFDs across species. We evaluate the regularities and differences of single-species dynamics more thoroughly in Section 3.4.

The most striking feature of these distributions, however, is the power-law $x^{-\nu}$ traced for intermediate abundances. This range is bounded at low abundances by the immigration rate and at high abundances by the single-species carrying capacity. The power-law exponent is $\nu \approx 1.18$ for the reference simulation, but it varies in general with the ecological parameters, as we discuss further in the following sections.

The regularity of the abundance distributions across species suggests that it may be possible to describe the dynamics of a 'typical' species in a compact way—this is the
<figure><img src="image_3.png" /><figcaption>**Figure 3: Comparison of the stochastic focal-species model and the chaotic dLV model.** A Time series of one arbitrary species in the disordered Lotka-Volterra (dLV) model (blue), and one realization of the stochastic focal-species model (Eq. (7)) with parameters as in Eq. (9): the time series are statistically similar. B Comparison of the average abundance fluctuation distribution (AFD) from Figure 2 (black), and the AFD of the focal-species model (pink): excellent agreement is found for the power-law section. The 'unified coloured noise approximation' solution for the focal-species model's AFD (dashed, pink line) predicts the correct overall shape of the distribution, but not a quantitatively accurate value for the power-law exponent.</figcaption></figure>

goal of the next section.

## 3.3 A stochastic focal-species model reproduces boom-bust dynamics

Fluctuating abundance time series are often fitted by one-dimensional stochastic models [7]; for example, stochastic logistic growth has been found to capture the statistics of fluctuations in a variety of datasets on microbial abundances [53, 54]. The noise term encapsulates variations in a species' growth rate whose origin may not be known explicitly. In our virtual Lotka-Volterra community, once the interaction matrix and initial abundances have been fixed, there is no uncertainty; nonetheless, the chaotic, high-dimensional dynamics results in species' growth rates fluctuating in a seemingly random fashion. We are therefore led to formulate a model for a single, focal species, for which explicit interactions are replaced by stochastic noise. Because we have found species to be statistically similar, its parameters do not depend on any particular species, but reflect the effective dynamics of any species in the community.

Following dynamical mean-field-like arguments and approximations informed by our simulations (Appendix E), we derive the focal-species model

$$ \dot{x}(t) = x(t)(g(t) - x(t)) + \lambda, \quad (7a) $$

$$ g(t) = -k + u \eta(t), \quad (7b) $$

where $g(t)$ is a stochastic growth rate with mean $-k$, and fluctuations of magnitude $u$ and correlation time $\tau$. The process $\eta(t)$ is a coloured Gaussian noise with zero mean

and an autocorrelation that decays exponentially;

$$ \langle \eta \rangle = 0, \quad \langle \eta(t) \eta(t') \rangle = e^{-|t-t'|/\tau}, \quad (8) $$

where brackets denote averages over noise realizations. The connection between the ecological parameters $S, \mu, \sigma, \lambda$ and the resulting dynamics of the disordered Lotka-Volterra model in the chaotic phase is then broken down into two steps: how the effective parameters $k, u, \tau$ relate to the ecological parameters; and how the behaviour of the focal-species model depends on the effective parameters.

For the first step we find

$$ k = \mu \bar{X} - 1 \quad \text{and} \quad u = \sigma \frac{\bar{X}}{\sqrt{S_{\text{eff}}}}, \quad (9) $$

where $X$ is the total community abundance of the original dynamics Eq. (1), the effective community size $S_{\text{eff}}$ is as in Eq. (6), and an overline denotes a long-time average. Equation (9) relates the focal species' growth rate to the time-averaged net competition ($\approx \mu \bar{X}$) from all other species. We find in simulations of Eq. (1) in the chaotic phase that competition is strong enough to make $k > 0$. The second relation captures the variation in the net competition that a species experiences because of turnover of the dominant community component. Due to sampling statistics, this variation is larger when the dominant component tends to have fewer species; hence the dependence on $(S_{\text{eff}})^{-1/2}$. The third effective parameter, the timescale $\tau$, controls how long the focal species stays dominant, once a fluctuation has brought it to high abundance. This timescale is essentially equal to the turnover timescale $\tau_{\text{dom}}$ of the dominant component (defined more precisely by autocorrelation functions in Appendix E). In the weak-interaction regime, where any pair of species can be treated as effectively independent at all times, self-consistency relations such as $S\langle x \rangle = \bar{X}$ allow to implicitly express the focal-species model in terms of the ecological parameters. For strong interactions, however, the disproportionate effect of the few dominant species on the whole community invalidates this approach; we therefore relate the effective parameters to the community-level observables $\bar{X}, S_{\text{eff}}, \tau_{\text{dom}}$ which are obtained from simulation of Eq. (1) at given values of the ecological parameters.

For the second step, we would like to solve Eqs. (7) for general values of the effective parameters. While this is intractable due to the finite correlation time of the noise, the equations can be simulated and treated by approximate analytical techniques. In Figure 3A we compare the time series of an arbitrary species in the dLV model with a simulation of the focal-species model. By eye, the time series appear statistically similar. The typical abundance of a species can be estimated by replacing the fluctuating growth rate in Eq. (7) with its typical value (i.e. $\eta = 0$), yielding the equilibrium $\lambda/k$ if $k > 0$, as indeed confirmed by the simulation. Thus the typical abundance value is on the order of the immigration threshold. Figure 3B shows that the average AFD of the dLV agrees remarkably well with the stationary distribution of the focal-species model, in particular for the power-law section. Using the unified coloured noise approximation
<figure><img src="image_4.png" /><figcaption>**Figure 4:** Species differences in dominance. A Example of a long abundance time series for the three species who are ranked first, median, and last, with respect to the 'dominance bias' (fraction of time spent in the dominant component relative to the species median). Some species 'boom' more often than others. B The scaling of median fraction of time spent in the dominant component against reciprocal species pool size: increasing S results in a proportional decrease in median dominance time. C Distribution of dominance biases against relative dominance rank for a range of S: there appears to be convergence towards a non-constant limiting distribution, implying that net species differences are not due to small-S effects. Note that, by definition, the dominance bias is 1 for the middle rank, indicated by the dashed line separating positively from negatively biased species. D Scatter of dominance bias against the net interaction bias, zᵢ Eq. (4): lower net competition correlates with higher dominance bias. Species in the tails of the zᵢ distribution are also less 'typical', with typicality quantified by the index θᵢ, Eq. (16), representing the similarity of a species AFD to the species-averaged AFD. Panel A and D are both for S = 500.</figcaption></figure>

[55] (Appendix F), one predicts that the stationary distribution, for $\lambda \ll x \ll 1$, takes the power-law form $x^{-\nu}$, where the exponent

$$v = 1 + \frac{k}{u^2 \tau} \qquad (10)$$

is strictly larger than one—the value predicted for weak interactions [41] and for neutral models [56]. Even if Eq. (10) is not quantitatively precise (Figure 3B), this formula suggests a scaling with the effective parameters that we will discuss later on.

## 3.4 Species with lower net competition are more often dominant

The similarity of all species' abundance fluctuation distributions in Figure 2 is reflected in the focal-species model's dependence on collective properties like the total abundance. However, the logarithmic scale downplays the variance between species' AFDs, particularly at higher abundances. Indeed, while all abundances fluctuate over orders of magnitude, some species are observed to be more often dominant (or rare). Such differences are reminiscent of the distinction between 'frequent' and 'occasional' species observed in empirical time series [57, 58].

In order to assess the nature of species differences in simulations of chaotic dLV, we rank species by the fraction of time spent as part of the dominant component. Observing the community dynamics on a very long timescale of tens of thousands of generations (400 times longer than in Figure 1), the first-ranked species appears to boom much more often than the last (Figure 4A). The frequency of a species is chiefly determined by the number of booms rather than their duration, which is comparable for all species. The median dominance time decreases with the total species richness

(Figure 4B): a doubling of S leads to each species halving its dominance time fraction. As the community gets crowded—while its effective size hardly increases, as remarked in Section 3.1—all species become temporally more constrained in their capacity to boom. Yet some significant fraction of species is biased towards booming much more often or rarely than the median, regardless of community richness. We quantify this trend by plotting in Figure 4C the dominance bias—the dominance time fraction normalized by the median across all species—against the relative rank (i.e., rank divided by S). For high richness ($S \sim 10^3$), the distribution of bias converges towards a characteristic, nonlinearly decreasing shape, where the most frequent species occur more than four times as often as the median, and the last-ranked species almost zero.

The persistence of inter-species differences with large S may seem to contradict the central limit theorem, as species' sets of interaction coefficients converge towards statistics that are identical for every species. In the chaotic regime, however, even the smallest differences in growth rates get amplified during a boom. As we show in Appendix D, if Eq. (1) is rewritten in terms of the proportions $p_i$, the relative advantage of species i is quantified by a selection coefficient whose time average scales as $-S^{-1/2}\sigma z_i$. Correspondingly, the relative, time-averaged growth rate is proportional to the net interaction bias $z_i$ (defined in Eq. (4)), resulting in species with larger $z_i$ to have positive dominance bias (Figure 4D). Outliers of the scatter plot, i.e. species that have particularly high or low dominance ranks, are also the species whose AFD is furthest from the average AFD of the community, as quantified by the typicality index $\theta_i \in [0, 1]$, defined in Appendix B.

In conclusion, the relative species-to-species variation in
<figure><img src="image_5.png" /><figcaption>**Figure 5: Dynamical phases of the disordered Lotka-Volterra model as a function of the interaction mean and standard deviation.** Probability of persistent chaos in long-time simulations: for each $\mu$ and $\sigma$ (with 0.01 increment), 30 simulations were made, each with a different random initial condition $x_i \sim U(\lambda, 2/S)$ and realization of the interaction matrix. Parameters yielding divergence every time are marked with grey. The boundary separating the chaotic phase from the rest of the multiple-attractor phase (in which cycles and multi-stable fixed point are common in addition to chaos) is not sharp, unless probed adiabatically in the way explained in Supplementary Figure S8. The unique fixed-point phase has been studied analytically in the weak-interaction regime ($\mu \sim 1/S$). When inter-specific competition is in general stronger than intra-specific competition, a single species (identity depending on initial condition) dominates, in line with the classical competitive exclusion principle [59].</figcaption></figure>

the total interaction strength drives the long-term differences in the dynamics of single species in the community. While the focal-species model emphasizes the similarity of species, species differences can also be taken into account by employing species-specific effective parameters. In particular, replacing $k$ with a distribution of $k_i$'s would create a dominance bias, and is in fact motivated upon closer examination of our focal-species model derivation (Figure 7D in Appendix E).

## 3.5 Interaction statistics control different dynamical phases

Hitherto, we have focussed on reference values of the interaction statistics $\mu$ and $\sigma$ that produce chaotic turnover of species abundances. We now broaden our investigation to determine the extent of validity of our previous analysis when the interaction statistics are varied. For every pair of $(\mu, \sigma)$ values, we run 30 independent simulations, each with a different sampling of the interaction matrix and uniformly sampled abundance initial condition. After a transient has elapsed, we classify the trajectory as belonging to one of four different classes: equilibrium, cycle, chaos, or divergence. Figure 5 displays the probability of observing chaos, demonstrating that it does not require fine tuning of parameters, but rather occurs across a broad parameter range.

The parameter region where chaos is prevalent, the ‘chaotic phase’, borders on regions of qualitatively different community dynamics. For small variation in inter-

action strengths (below the line connecting $(0, \sqrt{2/S})$ to $(1, 0)$), the community has a unique, global equilibrium that is fully characterized for weak interactions (c.f. Fig. 2. of [34]). The transition from equilibrium to chaos has been investigated with dynamical mean-field theory [41]. For low interaction variance, but with mean exceeding the unitary strength of intra-specific competition, a single species comes to dominate, as expected by the competitive exclusion principle [59]. Adiabatic simulations, implemented by continuously rescaling a single realization of the interaction matrix (details in Supplementary Figure S8), reveal that lines radiating from the point $(\mu, \sigma) = (1, 0)$ separate sectors where stable fixed points have different numbers of coexisting species. Traversing these sectors anti-clockwise, $S_{eff}$ increases by near-integer steps from one (full exclusion) up to about 8. From thence, a sudden transition to chaos occurs at the dashed line in Figure 5. We note, however, that the parameter region between chaos and competitive exclusion contains attractors of different types: cycles and chaos, coexisting with multiple fixed points, resulting in hysteresis (Supplementary Figure S8B). This ‘multiple attractor phase’ [34, 41] is a complicated and mostly uncharted territory whose detailed exploration goes beyond the scope of this study. Finally, for large variation in interactions, some abundances diverge due to the positive feedback loop induced by strongly mutualistic interactions, and the model is biologically unsound.

Across the phase diagram, community-level observables such as the average total abundance $\bar{X}$ and effective community size $\bar{S}_{eff}$ vary considerably (Supplementary Figure S9). The weak-interaction regime (whether in the equilibrium or chaotic phase) allows for high diversity, so $\bar{X}$ and $\bar{S}_{eff}$ are of order $S$; strong interactions, on the other hand, imply low diversity, with $\bar{S}_{eff}$ and $\bar{X}$ of order unity. An explicit expression for how these community-level observables depend on the ecological parameters $(S, \mu, \sigma, \lambda)$ is intractable (although implicit formulas exist in the weak-interaction regime [34]). Nonetheless, an approximate formula that we derive in Appendix C allows to relate community-level observables to one another and to $\mu$ and $\sigma$:

$$ \bar{X} \approx \left[ \mu + \frac{1-\mu}{\bar{S}_{eff}} - \sigma \bar{\rho} \right]^{-1}, \quad (11) $$

in which we introduce the collective correlation

$$ \bar{\rho} := - \sum_{ij} z_{ij} \overline{p_i p_j}, \quad (12) $$

involving the time-averaged product of relative abundances weighted by the their normalized interaction coefficient Eq. (4). By construction, the collective correlation is close to zero when all species abundances are uncorrelated over long times, as would follow from weak interactions. On the contrary, it is positive when pairs of species with interactions less competitive than the average tend to co-occur, and/or those with more competitive interactions tend to exclude one another.

Eq. (11) is particularly useful in understanding the role of correlations in the chaotic phase. As we observed in Section 3.3, the effective parameter $k = \mu \bar{X} - 1$ is positive in the
chaotic phase, implying that the growth rate of a species is typically negative, and abundances are therefore typically on the order of the small immigration rate rather than near carrying capacity. The existence of these two ‘poles’ of abundance values is key to boom-bust dynamics. By combining $k > 0$ with Eq. (11), we estimate a minimum, critical value of the collective correlation required for boom-bust dynamics:

$$ \rho_c = \frac{1 - \mu}{\sigma} \frac{1}{\bar{S}_{\text{eff}}}. \qquad (13) $$

Numerical simulations demonstrate that $\bar{\rho} \gtrsim \rho_c$ in the chaotic phase, where the critical value is approached at the boundary with the unique-equilibrium phase (Supplementary Figure S11). With this result in hand, Eq. (11) and Eq. (13) establish that $\bar{X} \gtrsim 1/\mu$ in the chaotic phase. For strong interactions, total abundances are predicted to be of order one, and for weak interactions $\bar{X} \approx S/\bar{\mu}$ (recall Eq. (5)), which recovers the observed scalings of these observables. As one moves deeper into the chaotic phase, the collective correlation increases continuously, as the effective community size drops, suggesting a seamless transition from a weak-interaction, chaotic regime amenable to exact treatment [41, 60], to the strongly correlated regime that we have analyzed by simulations and the approximate focal-species model.

## 3.6 Self-organization between community-level observables constrains abundance power-law variation

In Section 3.3 we established a focal-species model depending on the effective parameters $k, u$, and $\tau$, that were related to the ecological parameters $S, \mu, \sigma, \lambda$ indirectly via community-level observables $\bar{X}, \bar{S}_{\text{eff}}, \tau_{\text{dom}}$. Furthermore, in the previous section we studied how the latter vary in the chaotic phase. Putting these results together, we here examine the corresponding variation of the effective parameters and of the focal-species model’s predictions.

Because the trio $k, u, \tau$ ultimately derives from only two independent variables, $\mu, \sigma$ (considering fixed $S, \lambda$), they must be dependent. Figure 6A demonstrates that, across the chaotic phase, an approximate linear relationship holds between $k$ and $u$, as well as between $u$ and $\tau$. Because $k$ and $u$ are related to the mean and the variance of abundances via Eq. (9), their proportionality is reminiscent of the empirical Taylor’s law which posits a power-law relation between abundance mean and variance as they vary across samples [61]. The slope of the relationship of $u$ to $k$ is close to one (and varying little with $S$ and $\lambda$; Supplementary Figure S10), which implies with Eq. (9) that

$$ \bar{X} \approx \left[ \mu - \frac{\sigma}{\sqrt{\bar{S}_{\text{e}}}} \right]^{-1}. \qquad (14) $$

Comparison to Eq. (11) then yields that $\rho - \rho_c \approx \bar{S}_{\text{eff}}^{-1/2}$. This empirical relationship thus supports the aforementioned convergence—in the limit where $\bar{S}_{\text{eff}}$ is large, as for weak interactions—of the collective correlation to its critical value.

We find in Figure 6 that the slope $\nu_{\text{foc}}$ of the power-law trend obtained from simulation of the focal-species model finds good agreement with the value $\nu$ from the full dLV model. There is a narrow overall variation of the exponent; a consequence of the interdependency of the effective parameters. As can be intuited by the approximate expression Eq. (10) for the focal-species model, the exponent is strictly larger than 1, a value it approaches if the turnover time scale diverges, as indeed it does on the boundary to the unique equilibrium phase. The exponent increases as interactions become more competitive, up to about 1.4 at $(\mu, \sigma) = (1, 0)$. However, the exponent also depends on $S$ and $\lambda$, showing a constant slope against $\log S$ or $-1/\log \lambda$ (Supplementary Figure S7).

## 4 Discussion

Following growing empirical evidence for the presence of ecological chaos in natural and synthetic communities [7, 19], and increasing interest in the role of the rare biosphere [10, 62], we have sought a connection between the two through a minimal model of community dynamics: the disordered Lotka-Volterra (dLV) model with strong interactions and weak immigration. Our analysis of this model by extensive simulations, and through the derivation of an effective focal-species model, showed that: first, persistent chaos arises generically, and can drive fast and extensive turnover of rare and abundant species; second, a statistical equivalence between species emerges such that a single focal species’ fluctuation statistics predict the largely invariant power-law abundance distributions; third, deviations from this equivalence are associated with species differences in frequency of occurrence. In the following we discuss the generality of these results and their interpretation in the context of plankton ecology.

The chaotic turnover of rare and abundant species occurs because every subset of species that could stably coexist at high abundances is invadable by some rare species. This phenomenon should be robust to generalizations of the model as long as the dominant component remains exposed to a sufficient diversity of potential invaders, and the niche space that underlies species interactions contains enough trade-offs that no species can be a superior competitor across many biotic contexts. Our simplifying assumptions such as uniform growth rates and carrying capacities, and uncorrelated interactions can be relaxed (see our limited explorations in Supplementary Figure S5). Additional sources of modest noise should not cancel the deterministic contributions to fluctuations; indeed, the dynamical phases we have indicated are qualitatively similar to those arising in an individual-based version of the dLV model accounting for demographic stochasticity [23]. On the other hand, if the connectivity of the interaction network were reduced, lowering the exposure to competitors, one might expect a loss of persistent chaos at some critical connectance value [63]. Highly structured and hierarchical interactions would also undermine autonomous turnover on ecological timescales.

On a more technical note, the type of chaos we observe is likely ‘chaotic itinerancy’ [49, 64]. Lotka-Volterra systems
<figure><img src="image_6.png" /><figcaption>**Figure 6:** Relations between effective parameters in the chaotic phase. A Colour legend of the chaotic phase (boundaries from Figure 5). Each pair of $(\mu, \sigma)$ has been mapped to a distinct colour. B Co-dependence of the effective parameters $u, k, \tau$: the amplitude $u$ of growth-rate fluctuations approximately equals the absolute value $k$ of the negative growth rate (only weakly depending on $\lambda$ and $S$; Supplementary Figure S10); $u$ is roughly proportional to the inverse turnover time, but the slope of the relationship depends on $\lambda$ and $S$. C The exponent $\nu$ of the power-law section of the AFD for the chaotic dLV model plotted against the analogue $\nu_{\text{ foc}}$ obtained for the focal-species model: generally good agreement is found, with more outliers for parameters close to phase boundaries. A few outliers lie beyond the plotted range. Exponents have been estimated by fitting a power-law in the interval [100$\lambda$, 0.01] of the abundance distribution.</figcaption></figure>

without immigration admit heteroclinic networks [65–67]; equilibria with stable and unstable directions (i.e. saddle points) connected by orbits. Without immigration, such saddles are found on the system boundary, corresponding to some subset of species being extinct—in our case, these are the low-diversity equilibria reflected in the dominant component. The chaotic attractors appear when the saddles are ‘pushed off’ the boundary by the immigration term. Consistent with chaotic itinerancy in the dLV, characteristics of heteroclinic orbits—dynamical slowdown and ‘aging’—appear in the limit of vanishing immigration [60, 68].

While the assumption of disordered interactions may appear *ad hoc*, predictions for the onset of instability by the dLV model qualitatively match experiments in synthetic bacterial communities [19]. In a plankton context, we take the dLV to be a minimal yet relevant phenomenological representation of the relationships between species (or ‘operational taxonomic units’ from sequencing) of marine protists of a similar size class: the protistan interactome is largely uncharted [69], the ubiquity of mixoplankton blurs consumer–resource distinctions [70], and the effects of a diversity of zooplankton and viruses can manifest as apparent competition between species.

For rare plankton protists, the empirical snapshot SADs show a clear power-law trend, with an exponent around 1.6, varying little between different locations in the world oceans, despite large composition differences across samples [12]. The unified neutral theory of biodiversity, based on the interchangeability of individuals regardless of species identity, predicts a power-law tail of the SAD with exponent one [3, 56]. To approach the empirical value, previous studies augmented neutral theory with nonlinear growth rates [12] or chaotic mixing [71] to find an exponent dependent on the model parameters. However, for large census sizes such as that of plankton communities, neutral theories predict astronomically large turnover timescales [72, 73], inconsistent with observation. As we have shown, the dLV exhibits fast turnover when interactions are strong and sufficiently varied. For this model, $\nu \to 1$ as immigration tends to zero (Supplementary Figure S7; also shown in the weak-

interaction limit [28, 60]), but, if interactions are not weak, $\nu$ is substantially larger than one for small but finite values of immigration. The approximate solution to the focal-species model, Eq. (10), shows that the positive deviation from $\nu = 1$ depends on three inter-related effective parameters: the mean, amplitude, and timescale of fluctuations in each species’ net competition. As these fluctuations drive the turnover pattern, boom-bust dynamics comes to be associated with a larger-than-one exponent. The relatively weak variation of $\nu$ across the space of ecological parameters moreover suggests a reason for the limited geographical heterogeneity of the empirical value of the exponent.

A role for chaos in the plankton has long been advocated for [29, 74]. Proposed mechanisms include coupling of population dynamics to chaotically fluctuating environmental variables [75, 76], nonlinearity of low-dimensional zooplankton–phytoplankton dynamics [77], resource competition between phytoplankon species [20], the effect of marine viruses on populations of cyanobacteria strains [26]. These possibilities are not mutually exclusive, but relevant at different scales, from coarser to finer levels of taxonomic resolution. Adding a degree of structure to our species-level model to represent multiple functional groups would offer a way to investigate the connection between fluctuations at different scales. Empirical findings to replicate are the weakening signal for chaos as taxa are aggregated at higher orders [8], and more dynamical regularity and predictability in succession patterns at the level of functional groups [78]. In fact, even our unstructured model captures the feature that fluctuations are less severe at the aggregated level (e.g. total biomass, the envelope in Figure 1).

Besides explicit incorporation of structured interactions, an extension of our model with particular biological relevance would be to allow interactions to evolve, notably as they are reshaped by the appearance of novel species—a different scenario than our immigration term captures. Persistent turnover can then manifest on long timescales even if the ecological dynamics is—contrary to our case—at equilibrium. Such turnover has been shown in numerical models of evolving food webs structured by body size [79, 80] and
when adaptive dynamics occurs in high-dimensional trait spaces [81]. An open question is what evolutionary process may produce interactions that underpin chaotic turnover on ecological timescales. The observation that evolution sustains higher diversity under boom-bust ecology than under equilibrium ecology [82], together with the propensity of diversity to cause instability, suggests a possible role for evolutionary feedbacks.

Our approximate derivation of an explicit focal-species model demonstrates how ecological chaos comes to resemble noise. Parallel work to ours shows that an exact but implicitly defined effective model can be derived in the combined limit of weak interactions and infinitesimal immigration, where compositional turnover is slow [60]. In our model, the effective parameters could be used in fitting observational time series. Formally, Eq. (7) is similar to heuristic stochastic single-species logistic growth models that predict empirical distributions of microbial abundances [53, 54]. A notable difference lies in the negative mean growth rate we find, which together with noise-correlation and immigration yields fluctuations over many orders of magnitude, from rare to abundant. An insight from our model is that a species may be rare for an exceedingly long time, without rarity being a permanent character. On the other hand, species differences in the propensity to become abundant could reflect small differences in effective parameters that depend on a multitude of factors, which—like interaction rates—might not be individually measurable with precision. Together, these findings suggest that the abundances of particular species may not be easily explained by their traits, should fluctuations be determined by community complexity rather than a more direct coupling to environmental variables. In closing, a comparison of time series data to focal-species models could provide a complement to non-parametric methods [6, 7] in establishing the plausibility of ecological chaos as a driver of abundance fluctuations.

# Appendices

## A Numerical implementation

For Lotka-Volterra simulations we used a fixed time-step Euler scheme with $\Delta t = 0.01$, applied to the logarithm of abundances. This guarantees the positivity of all abundances at all times, regardless of immigration rate. To automatically classify the long-time behaviour of trajectories as fixed-points, cycles, or chaos, we used a heuristic method of counting abundance vector recurrences, validated against visual inspection of trajectories and calculated maximal Lyapunov exponent for a subset of trajectories. Further details are given in Supplementary Note S2.

All code was deposited on Zenodo: https://doi.org/10.5281/zenodo.10646601

## B Similarity metrics

The Bray-Curtis similarity index [83] is defined as

$$BC(\mathbf{x}, \mathbf{y}) = \sum_{i} w_{i} \frac{\min(x_{i}, y_{i})}{\text{mean}(x_{i}, y_{i})}, \quad (15)$$

where $w_i$ is the relative abundance of species i with respect to the joined abundances $\mathbf{x} + \mathbf{y}$. By definition, $BC(\mathbf{x}, \mathbf{y}) = 1$ iff $\mathbf{x} = \mathbf{y}$, and $BC \approx 0$ when, for each i, either $x_i \gg y_i$ or $y_i \gg x_i$; this makes it suitable for communities where abundances span orders of magnitude.

For the similarity graph Figure 1B, we have plotted $BC(\mathbf{x}^{\text{dom}}(t), \mathbf{y}^{*}(t))$, where $\mathbf{x}^{\text{dom}}(t)$ is the restriction of $\mathbf{x}(t)$ in the reference simulation to only the dominant species at time $t$, and $\mathbf{y}^{*}(t)$ is the fixed point reached from $\mathbf{x}^{\text{dom}}(t)$ as initial condition, with $\lambda = 0$.

To compare the similarity the AFD of species i, $P_i(\mathbf{x})$, to the species-averaged AFD $P = \sum_i P_i/S$, we define the index

$$\theta_i := 1 - \sup_{\mathbf{x}} |F(\mathbf{x}) - F_i(\mathbf{x})|, \quad (16)$$

where $F_i$ and $F$ are the cumulative distribution functions of $P_i$ and $P$, respectively; i.e., the index $\theta_i$ is based on the Kolmogorov-Smirnov distance [47] of the AFDs.

## C Derivation of time-averaged total abundance

Direct summation of Eq. (1) over i (assuming $r_i = 1$), and then division on both sides by $X(t) = \sum_i x_i(t)$, yields

$$\frac{d}{dt} \ln X(t) = 1 - X(t) \left[ \mu + (1 - \mu)/S_{\text{eff}} - \sigma \rho(t) \right] + \frac{S\lambda}{X(t)} \quad (17)$$

with $S_{\text{eff}}$ as in Eq. (6) and $\rho(t)$ as Eq. (12) but without the time average. Taking the long-time average of Eq. (17), the left-hand side becomes $\lim_{T \to \infty} (\ln X(T) - \ln X(0))/T$, which evaluates to zero on the assumption that no species diverges in abundance. The right-hand side contains terms such as $\overline{X}/S_{\text{eff}}$ and $\overline{X}\rho$. If the relative fluctuations in $X$, $S_{\text{eff}}$, $\rho$ are small (see Supplementary Figure S9), or these functions are at most weakly correlated to one another, then we obtain, approximately,

$$0 = 1 - \mu \overline{X} - (1 - \mu) \overline{X} / \overline{S}_{\text{eff}} + \sigma \overline{\rho} + O(S\lambda). \quad (18)$$

We neglect the small immigration term; solving for $\overline{X}$ then finds Eq. (11). The relative error in $\overline{X}$ between Eq. (11) for simulated values of the community-level observables in the right-hand side, and the simulated value of $\overline{X}$, is typically less than a few percent (see Supplementary Figure S14).

## D Selective advantage

The dynamics of the relative abundance $p_i = x_i/X$ is found by summing and differentiating Eq. (1) as

$$\dot{p}_i = X p_i \left\{ \sum_j p_j^2 - p_i - \sum_j p_j \left( \alpha_{ij} - \sum_k \alpha_{jk} p_k \right) \right\} + \frac{\lambda S}{X} \left( \frac{1}{S} - p_i \right). \quad (19)$$

Using Eqs. (3), (6), and (12) in defining

$$\pi(t) := 1/S_{\text{eff}}(t) - \sigma\rho(t), \quad s_i(t) := -\sigma \sum_j z_{ij} p_j, \quad (20)$$
we can write Eq. (19) as

$$\dot{p}_i = X p_i \{\pi + s_i - p_i\} + O(\lambda). \quad (21)$$

The term $s_i$ is responsible for the bias of species i against the reference proportion $\pi$. As a heuristic means of calculating the time-averaged bias, we suppose the $p_i$'s can be treated independently of the $z_{ij}$ and be replaces by $\bar{p}_i \approx 1/S$; then we obtain $\bar{s}_i \approx -\sigma z_i / \sqrt{S}$. On this basis, we expect $z_i$ to be indicative of a species' dominance bias.

## E Derivation of the stochastic focal-species model from dynamical mean-field arguments

We write Eq. (1) as

$$\dot{x}_i = x_i(g_i - x_i) + \lambda, \quad g_i = 1 - \sum_{j(\neq i)} \alpha_{ij} x_j. \quad (22)$$

If we suppose that the abundances $\{x_j(t)\}$ (or, rather, their statistical properties) are independent of the particular realization $[\alpha_{ij}]$ of the interaction matrix, then, for a given realization of $\{x_j(t)\}$,

$$g_i(t) \sim \mathcal{N} \left( 1 - \mu \sum_{j(\neq i)} x_j(t), \sigma^2 \sum_{j(\neq i)} x_j^2(t) \right), \quad (23)$$

based on the properties of sums of Gaussian variables. The time-varying mean and variance of $g_i$ means that, averaged over time, $g_i$ does not necessarily follow a Gaussian distribution. We introduce

$$a(t) := 1 - \mu \sum_{i} x_i(t), \quad b(t) := \sigma \sqrt{\sum_{i} x_i^2(t)}, \quad (24)$$

which are found to exhibit significant relative fluctuations, with skewed distributions (Figure 7A and B). However, once we shift and scale $g_i(t)$ into the “effective noise”

$$\eta_i(t) := \frac{g_i(t) - a(t)}{b(t)}, \quad (25)$$

we recover (closely) a $\mathcal{N}(0, 1)$ distribution, for both the set $\{\eta_i(t)\}_{1,...,S}$ at any given time $t$, and for the stationary distribution of $\eta_i(t)$, at least for typical species (Figure 7C and D). The empirical distribution of the $g_i$ across all species and times is closely approximated by the stationary distribution $\mathcal{N}(\bar{a}, \bar{b})$ (Figure 7E). Therefore, we suppose that, despite their fluctuations, we can replace $a(t)$ and $b(t)$ with their time-averages and model $g_i$ as a stochastic process $g(t) = \bar{a} + \bar{b}\eta(t)$, where $\eta(t)$ is a process with stationary distribution $\mathcal{N}(0, 1)$. The parameter correspondence in Eq. (9) follows by $k = -\bar{a}$, $u = \bar{b} \approx \bar{X}/\sqrt{\bar{S}_{\text{eff}}}$, and $\tau = \tau_{\eta}$, the correlation time of $\eta$.

Note that, up to neglecting a diagonal term of the sum, the effective noise can be written

$$\eta_i(t) = - \sum_{j(\neq i)} z_{ij} q_j(t), \quad (26)$$

with $z_{ij} \sim \mathcal{N}(0, 1)$, and $q(t) = x(t)/\|x(t)\|_2$. Given the chaotic turnover pattern, the latter is expected to perform

<figure><img src="image_7.png" /><figcaption>**Figure 7: Statistical properties of the effective noise.** A, B Time series and distribution of $a_{rel} = a/\bar{a} - 1$, etc. C, D Histograms of $\eta_i(t)$ across all species and time (grey), over just species for one random time (green), over all time for the first/mid/last-ranked species with respect to average abundance (blue/pink/yellow), with $\mathcal{N}(0, 1)$ (black, dashed) for reference. E The empirical distribution of $g$ in Eq. (22) over all species and times, compared to the distribution $\mathcal{N}(\bar{a}, \bar{b})$ assumed for $g$ in the focal-species model. F Autocorrelation functions: for every species (grey), first/mid/last-rank species (blue/pink/yellow)), and the average over all species (black). The left inset compares the ACFs of $x$ (green), $\eta$ (black), and the exponential fit to the latter (red); the right inset shows the distribution of the $\tau$ parameter in exponential fits to each species ACF.</figcaption></figure>

something like a random walk on the $S$-sphere, with a de-correlation time corresponding to the turnover of dominant species. This timescale is inherited by the effective noise. More precisely, we compare autocorrelation functions (ACF). The ACF of a function $f$ is defined as

$$\phi_f(t_{\text{lag}}) := \text{mean}_t[\delta f(t) \cdot \delta f(t + t_{\text{lag}})]/\text{var} f, \quad (27)$$

with $\delta f = f - \bar{f}$. By definition $\phi_f(0) = 1$. For each species' effective noise we compute numerically $\phi_{\eta_i}(t_{\text{lag}})$, as shown in Figure 7F. Due to the small number of 'booms' per species, even over a large simulation time, ACFs are slightly irregular. In order to make estimations more accurate, we consider the averaged ACF $\phi_{\eta} := S^{-1} \sum_i \phi_{\eta_i}$. The decay of correlation is well-approximated by the ex-
ponential $\exp(-t_{lag}/\tau_{\eta})$, where the parameter $\tau_{\eta}$ (fitted by least squares) represents the noise correlation timescale for a 'typical' species.

The approximately $N(0, 1)$ distribution and exponential autocorrelation function of the effective noise $\eta$ suggest that it can be modelled as an Ornstein-Uhlenbeck process, the only Markov process with these two properties;

$$\dot{\eta}(t) = -\frac{1}{\tau}\eta(t) + \sqrt{\frac{2}{\tau}}\xi(t), \quad (28)$$

where $\xi(t)$ is a Gaussian white noise; $\langle \xi \rangle = 0$, $\langle \xi(t)\xi(t') \rangle = \delta(t - t')$. The timescale referred to as $\tau_{dom}$ in the main text can be defined as $\tau_x$, the decay time of the exponential fit to the ACF of the abundance vector. For a vector-valued function, Eq. (27) gives

$$\phi_x = \frac{1}{S} \sum_i w_i \phi_{x_i}, \quad w_i = \frac{\mathrm{Var}[x_i]}{\frac{1}{S} \sum_j \mathrm{Var}[x_j]} \quad (29)$$

$\phi_x$ and $\phi_{x_i}$ match very well (inset of Figure 7F) for the reference simulation; as do the associated timescales $\tau_x$ and $\tau_{\eta}$ for all $(\mu, \sigma)$ in the chaotic phase (Supplementary Figure S13). This observation motivates identifying $\tau_{\eta}$ of the focal species model with the turnover timescale $\tau_{dom}$. Thus, the focal species model and its parameters have been fully specified.

We point out a critical difference between our approach and dynamical mean-field theory applied in the weak-interaction regime [34, 41]. Under strong interactions, $a(t)$ and $b(t)$ are determined by a small number ($S_{eff}$) of dominant species that, during the time of co-dominance, have strong effects on each other. Therefore, they can not be determined 'self-consistently' from the focal-species equations by assuming that every species can described simul-tanesouly as an independent realization of it. For example, the self-consistency relation for $k$ in Eq. (7b) is $k = \mu S \langle x \rangle - 1$. In our reference simulation $k = 0.26$, whereas $\mu S \langle x \rangle - 1 = -0.31$ even has the wrong sign. This discrepancy is due to the neglected inter-species correlations needed for the collective correlation $\bar{\rho}$, Eq. (12), to exceed the critical value Eq. (13) associated with $k > 0$ and boom-bust dynamics.

## F Steady-state solution of the focal-species model under the unified coloured noise approximation

The unified coloured noise approximation [55] assumes overdamped dynamics to replace a process $\dot{x} = F(x) + G(x)\eta$, driven by Gaussian correlated noise $\eta$ of correlation-time $\tau$, with a process driven by white noise. The approximation is exact in the limits $\tau \to 0$ or $\tau \to \infty$. The stationary distribution of the corresponding white-noise process is

$$P^*(x) \propto \exp \left\{ \int^x \nu(x') dx' \right\}, \quad (30)$$

with

$$\nu = (\tau^{-1/2} H_{\tau}) \frac{F}{G^2} + \left( \ln \frac{H_{\tau}}{G} \right)', \quad (31)$$

and $H_{\tau}$ a function of $F,G$, and $\tau$. For Eq. (7)

$$F(x) = -x(k + x) + \lambda, \quad G(x) = ux, \quad (32)$$

$$H_{\tau}(x) = \tau^{-1/2} + \tau^{1/2}(x + \lambda x^{-1}). \quad (33)$$

With these functions, the integral in Eq. (30) can be performed exactly, yielding

$$P^*(x) = \frac{1}{\mathcal{N}} e^{-[q_+(x)+q_-(\lambda/x)]} x^{-\nu} \left( \tau^{-1} + x + \frac{\lambda}{x} \right), \quad (34)$$

where $\nu$ is given by Eq. (10), and

$$q_{\pm}(y) := \frac{(y + (\tau^{-1} \pm k))^2}{2u^2}. \quad (35)$$

## Acknowledgements

We thank Matthieu Baron and Giulio Biroli for sharing preliminary work on the strong-interaction regime, and Jules Fraboul, Giulio Biroli, and Matthieu Barbier for valuable discussions in the course of this project and feedback on the manuscript.

## References

[1] D. L. DeAngelis and J. C. Waterhouse, Equilibrium and nonequilibrium concepts in ecological models, *Ecological Monographs* **57**, 1 (1987).

[2] M. Vellend, *The theory of ecological communities* (Princeton University Press, 2016).

[3] S. P. Hubbell, *The unified neutral theory of biogeography* (Princeton University Press, 2001).

[4] H. G. Schuster, *Deterministic Chaos*, 3rd (VCH, Weinheim, 1995).

[5] A. A. Berryman and J. A. Millstein, Are ecological systems chaotic – and if not, why not?, *Trends in Ecology & Evolution* **4**, 26 (1989).

[6] S. B. Munch, T. L. Rogers, B. J. Johnson, U. Bhat, and C.-H. Tsai, Rethinking the prevalence and pelevance of chaos in ecology, *Annual Review of Ecology, Evolution*, and *Systematics* **53**, 227 (2022).

[7] T. L. Rogers, B. J. Johnson, and S. B. Munch, Chaos is not rare in natural ecosystems, *Nature Ecology & Evolution* **6**, 1105 (2022).

[8] T. L. Rogers, S. B. Munch, S. ichiro S. Matsuzaki, and C. C. Symons, Intermittent instability is widespread in plankton communities, *Ecology Letters* **26**, 470 (2023).

[9] M. L. Sogin et al., Microbial diversity in the deep sea and the underexplored “rare biosphere”, *Proceedings of the National Academy of Sciences* **103**, 12115 (2006).

[10] M. D. J. Lynch and J. D. Neufeld, Ecology and exploration of the rare biosphere, *Nature Reviews Microbiology* **13**, 217 (2015).
[11] F. Pascoal, R. Costa, and C. Magalhães, The microbial rare biosphere: current concepts, methods and ecological principles, *FEMS Microbiology Ecology* **97**, 10.1093/femsec/fiaa227 (2021).

[12] E. Ser-Giacomi et al., Ubiquitous abundance distribution of non-dominant plankton across the global ocean, *Nature Ecology & Evolution* **2**, 1243 (2018).

[13] J. A. Fuhrman, J. A. Cram, and D. M. Needham, Marine microbial community dynamics and their ecological interpretation, *Nature Reviews Microbiology* **13**, 133 (2015).

[14] A. M. Martin-Platero et al., High resolution time series reveals cohesive but short-lived communities in coastal plankton, *Nature Communications* **9**, 266 (2018).

[15] J. A. Gilbert et al., Defining seasonal marine microbial community dynamics. *The ISME journal* **6**, 298 (2012).

[16] E. Benincà et al., Chaos in a long-term experiment with a plankton community, *Nature* **451**, 822 (2008).

[17] I. V. Telesh et al., Chaos theory discloses triggers and drivers of plankton dynamics in stable environment, *Sci. Rep.* **9**, 20351 (2019).

[18] L. Becks, F. M. Hilker, H. Malchow, K. Jürgens, and H. Arndt, Experimental demonstration of chaos in a microbial food web, *Nature* **435**, 1226 (2005).

[19] J. Hu, D. R. Amor, M. Barbier, G. Bunin, and J. Gore, Emergent phases of ecological diversity and dynamics mapped in microcosms, *Science* **378**, 85 (2022).

[20] J. Huisman and F. J. Weissing, Biodiversity of plankton by species oscillations and chaos, *Nature* **402**, 407 (1999).

[21] P. Schippers, A. M. Verschoor, M. Vos, and W. M. Mooij, Does "supersaturated coexistence" resolve the "paradox of the plankton"?, *Ecol. Lett.* **4**, 404 (2001).

[22] I. Ispolatov, V. Madhok, S. Allende, and M. Doebeli, Chaos in high-dimensional dissipative dynamical systems, *Sci. Rep.* **5**, 12506 (2015).

[23] D. A. Kessler and N. M. Shnerb, Generalized model of island biodiversity, *Phys. Rev. E* **91**, 042705 (2015).

[24] F. Roy, M. Barbier, G. Biroli, and G. Bunin, Complex interactions can create persistent fluctuations in high-diversity ecosystems, *PLOS Comput. Biol.* **16**, e1007827 (2020).

[25] J. D. O'Sullivan, J. C. D. Terry, and A. G. Rossberg, Intrinsic ecological dynamics drive biodiversity turnover in model metacommunities, *Nature Communications* **12**, 3627 (2021).

[26] M. T. Pearce, A. Agarwala, and D. S. Fisher, Stabilization of extensive fine-scale diversity by ecologically driven spatiotemporal chaos, *PNAS* **117**, 14572 (2020).

[27] P. Rodríguez-Sánchez, E. H. van Nes, and M. Scheffer, Neutral competition boosts cycles and chaos in simulated food webs, *Royal Society Open Science* **7**, 191532 (2020).

[28] I. Dalmedigos and G. Bunin, Dynamical persistence in high-diversity resource-consumer communities, *PLOS Computational Biology* **16**, e1008189 (2020).

[29] M. Scheffer, S. Rinaldi, J. Huisman, and F. J. Weissing, Why plankton communities have no equilibrium: solutions to the paradox, *Hydrobiologia* **491**, 9 (2003).

[30] J. Denk, S. Martis, and O. Hallatschek, Chaos may lurk under a cloak of neutrality, *PNAS* **117**, 16104 (2020).

[31] J. T. Lennon and S. E. Jones, Microbial seed banks: the ecological and evolutionary implications of dormancy, *Nature Reviews Microbiology* **9**, 119 (2011).

[32] R. May, Will a large complex system be stable?, *Nature* **238**, 413 (1972).

[33] S. Allesina and S. Tang, The stability-complexity relationship at age 40: a random matrix perspective, *Population Ecology* **57**, 63 (2015).

[34] G. Bunin, Ecological communities with Lotka-Volterra dynamics, *Phys. Rev. E* **95**, 042414 (2017).

[35] M. Barbier, J.-F. Arnoldi, G. Bunin, and M. Loreau, Generic assembly patterns in complex ecological communities, *PNAS* **115**, 2156 (2018).

[36] A. F. Ansari, Y. B. S. Reddy, J. Raut, and N. M. Dixit, An efficient and scalable top-down method for predicting structures of microbial communities, *Nature Computational Science* **1**, 619 (2021).

[37] A. Skwara et al., Statistically learning the functional landscape of microbial communities, *Nature Ecology & Evolution*, 10.1038/s41559-023-02197-4 (2023).

[38] J. Hofbauer and K. Sigmund, *Evolutionary games and population dynamics* (Cambridge University Press, 2002).

[39] M. Barbier and J.-F. Arnoldi, The cavity method for community ecology, *bioRxiv:10.1101/147728*, 10.1101/147728 (2017).

[40] T. Galla, Dynamically evolved community size and stability of random Lotka-Volterra ecosystems, *Europhys. Lett.* **123**, 48004 (2018).

[41] F. Roy, G. Biroli, G. Bunin, and C. Cammarota, Numerical implementation of dynamical mean field theory for disordered systems: application to the Lotka-Volterra model of ecosystems, *J. Phys. A: Math. Theor.* **52**, 484001 (2019).

[42] J. W. Baron, T. J. Jewell, C. Ryder, and T. Galla, Breakdown of random-matrix universality in persistent Lotka-Volterra communities, *Phys. Rev. Lett.* **130**, 137401 (2023).

[43] G. Biroli, G. Bunin, and C. Cammarota, Marginally stable equilibria in critical ecosystems, *New J. Phys.* **20**, 083051 (2018).
[44] A. Altieri, F. Roy, C. Cammarota, and G. Biroli, Properties of equilibria and glassy phases of the random Lotka-Volterra model with demographic noise, Phys. Rev. Lett. 126, 258301 (2021).

[45] O. S. Venturelli et al., Deciphering microbial interactions in synthetic human gut microbiome communities, Molecular Systems Biology 14, e8157 (2018).

[46] D. Machado et al., Polarization of microbial communities between competitive and cooperative metabolism, Nature Ecology & Evolution 5, 195 (2021).

[47] S. Xu, L. Böttcher, and T. Chou, Diversity in biology: definitions, quantification and models, Phys. Biol. 17, 031001 (2020).

[48] R. Logares et al., Patterns of Rare and Abundant Marine Microbial Eukaryotes, Current Biology 24, 813 (2014).

[49] K. Kaneko and I. Tsuda, Chaotic itinerancy, Chaos 13, 926 (2003).

[50] G. B. Arous, Y. V. Fyodorov, and B. A. Khoruzhenko, Counting equilibria of large complex systems by instability index, PNAS 118, e2023719118 (2021).

[51] V. Ros, F. Roy, G. Biroli, G. Bunin, and A. M. Turner, Generalized Lotka-Volterra equations with random, nonreciprocal interactions: the typical number of equilibria, Phys. Rev. Lett. 130, 257401 (2023).

[52] T. J. Matthews and R. J. Whittaker, Fitting and comparing competing models of the species abundance distribution: assessment and prospect, Frontiers of Biogeography 6, 67 (2014).

[53] J. Grilli, Macroecological laws describe variation and diversity in microbial communities, Nat. Commun. 11, 4743 (2020).

[54] L. Descheemaeker and S. de Buyl, Stochastic logistic models reproduce experimental time series of microbial communities, eLife 9, e55650 (2020).

[55] P. Jung and P. Hänggi, Dynamical systems: a unified colored-noise approximation, Phys. Rev. A 35, 4464 (1987).

[56] D. A. Kessler and N. M. Shnerb, Neutral selection, arXiv:1211.3609, 10.48550/arXiv:1211.3609 (2012).

[57] A. E. Magurran and P. A. Henderson, Explaining the excess of rare species in natural species abundance distributions, Nature 422, 714 (2003).

[58] W. Ulrich and M. Ollik, Frequent and occasional species and the shape of relative-abundance distributions, Diversity and Distributions 10, 263 (2004).

[59] G. Hardin, The competitive exclusion principle, Science 131, 1292 (1960).

[60] T. A. de Pirey and G. Bunin, Many-species ecological fluctuations as a jump process from the brink of extinction, arXiv:2306.13634, 10.48550/arXiv:2306.13634 (2023).

[61] Z. Eisler, I. Bartos, and J. Kertész, Fluctuation scaling in complex systems: taylor's law and beyond, Advances in Physics 57, 89 (2008).

[62] X. Jia, F. Dini-Andreote, and J. F. Salles, Community assembly processes of the microbial rare biosphere, Trends in Microbiology 26, 738 (2018).

[63] G. B. Stav Marcus Ari M. Turner, Local and extensive fluctuations in sparsely-interacting ecological communities, arXiv:2308.01828v1, 10.48550/arXiv:2308.01828 (2023).

[64] K. Hashimoto and T. Ikegami, Heteroclinic chaos, chaotic itinerancy and neutral attractors in symmetrical replicator equations with mutations, J. Phys. Soc. Japan 70, 349 (2001).

[65] J. Hofbauer, Heteroclinic cycles in ecological differential equations, Tatra Mountains Math. Publ. 4, 105 (1994).

[66] V. Afraimovich, I. Tristan, R. Huerta, and M. I. Rabinovich, Winnerless competition principle and prediction of the transient dynamics in a Lotka-Volterra model, Chaos 18, 043103 (2008).

[67] C. Bick and M. I. Rabinovich, On the occurrence of stable heteroclinic channels in Lotka-Volterra models, Dynamical Systems 25, 97 (2009).

[68] T. A. de Pirey and G. Bunin, Aging by near-extinctions in many-variable interacting populations, Physical Review Letters 130, 098401 (2023).

[69] M. F. M. Bjorbækmo, A. Evenstad, L. L. Røsæg, A. K. Krabberød, and R. Logares, The planktonic protist interactome: where do we stand after a century of research?, The ISME Journal 14, 544 (2019).

[70] N. C. Millette et al., Mixoplankton and mixotrophy: future research priorities, Journal of Plankton Research 45, 576 (2023).

[71] P. Villa Martín, A. Koldaeva, and S. Pigolotti, Coalescent dynamics of planktonic communities, Phys. Rev. E 106, 044408 (2022).

[72] D. Kessler, S. Suweis, M. Formentin, and N. M. Shnerb, Neutral dynamics with environmental noise: age-size statistics and species lifetimes, Physical Review E 92, 022722 (2015).

[73] M. J. Behrenfeld, R. O'Malley, E. Boss, L. Karp-Boss, and C. Mundt, Phytoplankton biodiversity and the inverted paradox, ISME Communications 1, 10.1038/s43705-021-00056-6 (2021).

[74] M. Scheffer, Should we expect strange attractors behind plankton dynamics – and if so, should we bother?, J. Plankton Res. 13, 1291 (1991).

[75] A. Bracco, A. Provenzale, and I. Scheuring, Mesoscale vortices and the paradox of the plankton, Proceedings of the Royal Society of London. Series B: Biological Sciences 267, 1795 (2000).

[76] F. d'Ovidio, S. De Monte, S. Alvain, Y. Dandonneau, and M. Lévy, Fluid dynamical niches of phytoplankton types, Proceedings of the National Academy of Sciences 107, 18366 (2010).
[77] A. B. Medvinsky et al., Chaos far from the edge of chaos: a recurrence quantification analysis of plankton time series, *Ecological Complexity* **23**, 61 (2015).

[78] C. M. Mutshinda, Z. V. Finkel, C. E. Widdicombe, and A. J. Irwin, Ecological equivalence of species within phytoplankton functional groups, *Functional Ecol-ogy* **30**, 1714 (2016).

[79] K. T. Allhoff, D. Ritterskamp, B. C. Rall, B. Drossel, and C. Guill, Evolutionary food web model based on body masses gives realistic networks with permanent species turnover, *Scientific Reports* **5**, 10.1038/ srep10955 (2015).

[80] M. Hamm and B. Drossel, The concerted emergence of well-known spatial and temporal ecological patterns in an evolutionary food web model in space, *Scientific Reports* **11**, 10.1038/s41598-021-84077-0 (2021).

[81] M. Doebeli and I. Ispolatov, Chaos and unpredictability in evolution, *Evolution* **68**, 1365 (2014).

[82] M. Doebeli, E. C. Jaque, and Y. Ispolatov, Boom-bust population dynamics increase diversity in evolving competitive communities, *Communications Biology* **4**, 502 (2021).

[83] J. R. Bray and J. T. Curtis, An ordination of the upland forest communities of southern wisconsin, *Ecological Monographs* **27**, 325 (1957).
# Supplementary material

## Chaotic turnover of rare and abundant species in a strongly interacting model community

Emil Mallmin<fnref n="1" />*, Arne Traulsen<fnref n="1" />, and Silvia De Monte<fnref n="1" />,<fnref n="2" />

Revised version, February 2024

## Contents

S1 Supplementary Figures 1

S2 Numerical implementation of model simulations 9

S3 Simplification of the disordered Lotka-Volterra model with mixing 10

S4 Invasion analysis of few-species equilibria 11

S5 Solution of intermittency model under unified coloured noise approximation 12

## List of Figures

<table><tr><td>S1</td><td>Sensitive dependence on model parameters</td><td>2</td></tr><tr><td>S2</td><td>Convergence to positive maximum Lyapunov exponent (MLE)</td><td>2</td></tr><tr><td>S3</td><td>Decay of community similarity with time.</td><td>3</td></tr><tr><td>S4</td><td>Scaling of effective community size with richness</td><td>3</td></tr><tr><td>S5</td><td>Robustness of turnover dynamics under model variations</td><td>4</td></tr><tr><td>S6</td><td>Pairwise correlations in species abundances</td><td>5</td></tr><tr><td>S7</td><td>Scaling of AFD power-law exponent with S and λ</td><td>5</td></tr><tr><td>S8</td><td>Phase diagram form adiabatic simulations</td><td>6</td></tr><tr><td>S9</td><td>Variation of community-level observables across the phase diagram</td><td>7</td></tr><tr><td>S10</td><td>Dependence of the effective parameters u, k on S, λ</td><td>7</td></tr><tr><td>S11</td><td>Collective correlation</td><td>8</td></tr><tr><td>S12</td><td>Power-law exponent ν in the chaotic phase</td><td>8</td></tr><tr><td>S13</td><td>Comparison of autocorrelation times</td><td>8</td></tr><tr><td>S14</td><td>Error in approximate formula for X&#772;</td><td>9</td></tr></table>

## S1 Supplementary Figures

As in the main text, unless otherwise stated, reference parameters in simulations are $S = 500$, $\mu = 0.5$, $\sigma = 0.3$, $\lambda = 10^{-8}$.

<footnote n="1">¹Max Planck Institute for Evolutionary Biology, Plön, Germany</footnote>

<footnote n="2">²Institut de Biologie de l'ENS (IBENS), Département de Biologie, Ecole Normale Supérieure, CNRS, INSERM, Université PSL, 75005 Paris, France</footnote>

<footnote>*mallmin@evolbio.mpg.de</footnote>
<figure><img src="image_8.png" /></figure>

**Figure S1: Sensitive dependence on model parameters.** A chaotic system exhibits sensitive dependence on initial conditions, and hence also on any model parameters or numerical implementation details that affect the dynamic variables. A Reference simulation, showing stacked abundances, similar to Main Text Figure 1A. B A change of integration scheme, with respect to the reference; C a perturbation of the interaction coefficients by $O(10^{-6})$; D a perturbation of the initial abundances by $O(10^{-8})$. E Each type of perturbation leads to completely different community composition compared to the reference (measured as Bray-Curtis similarity) after a few hundred time units.

<figure><img src="image_9.png" /></figure>

**Figure S2: Convergence to positive maximum Lyapunov exponent (MLE).** A The dominant finite-time Lyapunov exponent (FTLE) over a few integration time steps ($n = 2$) fluctuates along a trajectory, indicating the alternation of periods of phase-space expansion (boom) and contraction (bust). B The cumulative average of the FTLE converges towards a limit that is the maximal Lyapunov exponent. Its positive value (0.02) indicates that the trajectory is chaotic.
<figure><img src="image_10.png" /></figure>

**Figure S3: Decay of community similarity with time.** A The temporal similarity matrix $\mathcal{T}$ has elements given by the Bray-Curtis similarity between the abundance vectors at two time points, $\mathcal{T}(t, t') = BC(x(t), x(t'))$. Because only the diagonal elements are far from zero, and the similarity index is mostly determined by the overlap of dominant species, we conclude that the dominant component is not closely repeated (unless, perhaps, after an exceedingly long time). The aberration around $t \approx t' \approx 8000$ reflects a time when some dominant component persisted for an unusually long time. B For a few well-separated time points $t$ (one graph each), we show how $\mathcal{T}(t, t')$ decays over time $t'$ on a timescale of 200 time units (top panel), and how it fluctuates around a small value over a longer time scale of 5000 time units. Thus, community composition decorrelates quickly in time, with some residual low peaks in similarity reflecting that one or a few species will eventually reappear in a dominant community that is otherwise differently composed.

<figure><img src="image_11.png" /></figure>

**Figure S4: Scaling of effective community size with richness.** The time-average of the effective community size, $\bar{S}_{\text{eff}}$ (Main Text Eq. (6)), increases slowly (but super-logarithmically) with the overall richness $S$. That is, even if we add thousands of new species to the community, the dominant component at given time would just have a species or two more than before.
<figure><img src="image_12.png" /></figure>

**Figure S5: Robustness of turnover dynamics under model variations.** We here illustrate that chaotic dynamics is observed even when we relax the simplifying assumptions we made on model parameters in the Main Text; however, we leave a systematic investigation of these generalized scenarios for future work. **A** Non-uniform growth rates: we sample $r_i \sim U(0, 1)$. **B** Non-uniform carrying capacities: $K_i \sim \text{LogNorm}(0, 0.1)$. **C** Sparse interactions: each interaction has a 0.1 chance to be non-zero. **D** Symmetric bias: $\gamma = 0.2$ correlation between diagonally opposed interaction coefficients. **E** Predator-prey bias: $\gamma = -0.3$ correlation between diagonally opposed interaction coefficients.
<figure><img src="image_13.png" /></figure>

**Figure S6: Pairwise correlations in species abundances.** While most of the $S(S-1)$ pairs of species do not have meaningful levels of correlation over long times (here 100'000 time units), every species has some other species with which its correlation is substantial and non-spurious. The vertical axis has the correlation coefficient with lag time $t_{lag}$ $\phi_{ij}(t_{lag})$, and the horizontal axis has the rescaled interaction coefficient $z_{ij} = (\alpha_{ii} - \mu)/\sigma$. The inset show that most zero-lag correlation coefficients are close to zero; all zero-lag correlations are scattered in grey in the main plot. Blue (darker) points shows the values of maximum correlations $\max_j \phi_{ij}(0)$ for every species $i$; in order to see if correlations are stringer if we optimize over the delay time, we show in light blue $\max_{j,t_{lag}<200} \phi_{ij}(t_{lag})$. It is seen that the maximal correlations are around 0.25 in size, and clearly associated with $z_{ij} < 0$, i.e. a less-than-averagely negative (even positive) effect of species $j$ on $i$. Similarly, the extremal anti-correlations (pink for zero time lag, and light pink optimizing over time lag) are associated with $z_{ij} > 0$, i.e. a particularly negative effect of $j$ on $i$.

<figure><img src="image_14.png" /></figure>

**Figure S7: Scaling of AFD power-law exponent with S and $\lambda$.** From simulations, we have extracted the slope of the power-law section of the abundance fluctuation distribution (AFD). **A** For varying $S$, we find empirically that the exponent depends linearly on the logarithm of species richness, with coefficients that depend on the system's other parameters: $\nu = \nu_0 + c \log S$, where $\nu_0 = \nu_0(\mu, \sigma, \lambda)$, $c = c(\mu, \sigma, \lambda)$. **B** For varying $\lambda$, the exponent appears to follow $\nu = 1 - d/\log \lambda$, where $d = d(\mu, \sigma, S)$. The values of $\lambda$ are 10 to the power of negative 8, 12, 16, 20, 24, 28, 32, and 128 in order to extrapolate towards zero immigration. The dashed line connects $\nu = 1$ with the value at $\lambda = 10^{-8}$.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

**Figure S8:** Phase diagram for adiabatic simulations. Adiabatic simulations allow to track, in a numerically efficient fashion, the attractors of the dynamics as model parameters are changed slowly and continuously. To make the interaction statistics $\mu$ and $\sigma$ continuous parameters of the model, we use as interaction matrix $\alpha_{ij}(\mu, \sigma) = \mu + \sigma\zeta_{ij}$ where $\zeta$ is a single, fixed realization of a standard Gaussian random matrix. **A** For each value of $\sigma$, we initialized separate simulation runs starting at $\mu = 1.4$, and let their abundances evolve until an attractor was found. For each run, we then changed $\mu$ by small increments $\delta\mu = -0.1$, allowing enough time between each change for the abundances to relax from their previous state. This relaxation would either result in a small perturbation of the previous attractor, or instigate a jump to a different attractor. If a state diverged, the initial abundances for the next value of $\mu$ were set as the most recent non-divergent attractor. Thus, each simulation traced a sequence of attractors from $\mu = 1.4 \to -0.1$, corresponding to a horizontal line in the phase diagram. The colour quality reflects the class of the attractor, and the colour gradation indicates the effective community size, revealing the following features: First, we find mostly fixed points in the multiple attractor region. This is because, once a fixed point is converged to, it is "hold on to" until it vanishes or changes stability. If, instead, every simulation at given $\mu, \sigma$ would start from newly sampled initial abundances and interaction matrix, we would find different attractors every time, and the diagram becomes more heterogeneous (compare Main Text Figure 5). Second, clear lines radiate from $(\mu, \sigma) = (1, 0)$ and delineate sectors characterized by the number of high-abundance species coexisting at a fixed-point. In section S4 we show that an invasion analysis predicts such sectors, but not the right scaling of the lines' slope with $\bar{S}_{\text{eff}}$. Third, the jump from fixed-point to chaotic attractors occurs along a sharply defined line. **B** Stacked abundances of the attractor found in an adiabatic sequence $\mu = 1.4 \to 0.6$ (top panel, right to left) and the reverse $0.6 \to 1.4$ (bottom panel, left to right) at $\sigma = 0.3$. One can see sudden jumps to new equilibria involving more (or less) species. In the upper panel, reading right to left, a three-species equilibrium is found at $\mu = 1.15$, which jumps to a 6-species equilibrium by the invasion of three more species at $\mu = 1.11$; another two species displace one of the previous at $\mu = 0.9$; and at $\mu = 0.72$ a sudden jump onto a chaotic attractor occurs. Reversing the adiabatic protocol, the transition from chaos to fixed point occurs only at $\mu = 0.81$, and the sequence of equilibria is not identical to the forward direction (hysteresis). A systematic investigation of the multiple attractor phase and the transition to the chaotic phase is left for future work.
<figure><img src="image_17.png" /></figure>

**Figure S9:** Variation of community-level observables across the phase diagram. For the community-level observables in Main Text Eq. (11) we show: **A–C** their time-averaged values; **D–F** their relative relative fluctuations. The data comes from the adiabatic simulation detailed in the caption to Figure S8. An arrow on the end of the colour bar implies the range has been capped for clarity.

<figure><img src="image_18.png" /></figure>

**Figure S10:** Dependence of the effective parameters *u*, *k* on *S*, *λ*. The empirical, approximate relationship *u* ∝ *k*, found across the range of *μ*, *σ* in the chaotic phase, has a proportionality constant that depends relatively weakly on *S* and *λ*.
<figure><img src="image_19.png" /></figure>

**Figure S11: Collective correlation.** Within the bounds to the chaotic phase indicated in Figure S8A, we have run a long simulation for each parameter point with random initial condition and interaction matrix realization. Statistics were recorded for persistently chaotic trajectories; non-chaotic trajectories were discarded, and the parameter point rerun to obtain a long chaotic trajectory, up to five times, else the point was omitted (chaos probability was shown in Main Text Figure 5). **A** The collective correlation. **B** The critical value of the collective correlation as defined by Main Text Eq. (13). **C** The ratio $\bar{\rho}/\rho_c$ tends towards 1 at the boundary to the equilibrium phase. Note that $\bar{\rho}$ changes continuously across this boundary (Figure S9C). The arrow at the upper end of the colour bar implies the range has been capped for clarity.

<figure><img src="image_20.png" /></figure>

**Figure S12: Power-law exponent $\nu$ in the chaotic phase.** **A** Variation of the AFD power-law exponent across the chaotic phase. Apart from outliers, we find an exponent larger than one. **B** To test the accuracy of the focal-species model in predicting the exponent, we measure the relative error in $\delta\nu - 1$ (since we expect $\nu > 1$) with respect to the value obtained from simulations of the disordered Lotka-Volterra model. Data from the simulations described in Figure S11. The arrow at the upper end of the colour bar implies the range has been capped for clarity.

<figure><img src="image_21.png" /></figure>

**Figure S13: Comparison of autocorrelation times.** We compare the autocorrelation time $\tau_x$ of the abundance vector $\mathbf{x}$ and the autocorrelation time $\tau_\eta$ of the effective noise $\eta$. These two parameters are obtained by the exponential fit $e^{-t/\tau}$ applied to the respective autocorrelation functions. Across the chaotic phase, these timescale are quantitatively close, for reasons explained in Main Text Appendix E.
<figure><img src="image_22.png" /><figcaption>**Figure S14: Error in approximate formula for $\bar{X}$.** Here is shown that the approximate Main Text Eq. (11) generally gives an accurate prediction (small relative error) of $\bar{X}$ compared to its simulated value (adiabatic data; Figure S8), if given the values of $\bar{S}_{eff}$ and $\bar{\rho}$ from the same simulation. Except for close the divergent phase, the error is within ±2%. Since the approximations involved in deriving the formula amount to neglecting fluctuations, it is expected to most accurate when fluctuations are small (compare Figure S9); in particular, at fixed-points of the dynamics it becomes exact up to an amount proportional to the negligibly small immigration rate. An arrow on the end of the colour bar implies the range has been capped for clarity.</figcaption></figure>

## S2 Numerical implementation of model simulations

All numerical procedures of this work were carried out in python using the ` numpy ` and ` scipy ` packages. To simulate the disordered Lotka-Volterra model, we have opted for a fixed time integration with a small time step of $\Delta t = 0.01$. By default, we have used the logarithmic integration scheme defined below, whose implementation we have validated against a 'logistic' scheme (Figure S1) and the standard ODE solver of `scipy` using RK45. For simulations of the stochastic differential equation of the focal-species model we have simulated the coloured noise (Ornstein-Uhlenbeck process) under a Euler-Mayurama scheme with $\Delta t = 0.01$, and the abundance dynamics under the logistic scheme with the same time step.

In phase diagram simulations, we have used an expedient numerical heuristic to classify the long-term dynamics of trajectories, which we have validated against visual inspection, and measurement of the maximal Lyapunov exponent for a sample set of trajectories. First, if trajectories diverged, they tended to do so early in the simulation. Otherwise, after a transient interval $[t_0, t_1]$ of fixed duration, we hypothesised that the following time interval $[t_1, t_2]$ would contain stationary dynamics, assuming $t_2-t_1$ to be longer than any periodicity of the dynamics, if present. We then counted how many times $n$ within this interval the abundance vector $x(t)$ crossed a threshold of $1-\varepsilon$ similarity to the final vector $x(t_2)$, where $\varepsilon$ is a small tolerance, and the similarity metric is $d(x, y) = 1 - ||x - y||/(x^2 + y^2)$. If $n = 0$, then abundances were constant in the interval and we assume a stable fixed point has been reached; if $n = 1$, then the final composition was one not seen before in the interval, which we classified as chaos; if $n > 1$ then periodicity or quasi-periodic. For trajectories classified as chaotic, a subsequent long time interval $[t_2, t_3]$ was simulated and used to gather relevant statistics; finally, the $n$-classification was applied to a final interval $[t_3, t_4]$ to ascertain that chaotic dynamics were not lost during the previous interval of measurements.

### Definition of integration schemes

We consider the numerical integration of

$$ \dot{x}_i(t) = r_i x_i(t) (g_i(t) - x_i(t)), \quad g_i(t) = 1 - \sum_{j(\neq i)} \alpha_{ij} x_j(t). \qquad (S1) $$
Because abundances may become exponentially small, yet eventually recover, it is potentially problematic if numerical error can cause zero or negative abundances. To avoid this issue, we can consider the exact identity

$$x_i(t + \Delta t) = x_i(t) \exp \left( \int_t^{t+\Delta t} (g_i(t') - x_i(t')) dt' \right). \qquad (S2)$$

Since the exponential is always positive, any numerical integration scheme $x_i(n\Delta t) \rightarrow x_i[n]$ based on approximating the integral will preserve positivity of abundances. For a $\Delta t$ much smaller than the turnover timescale of the dominant community, the integrand can be treated as approximately constant, yielding the scheme

$$x_i[n+1] = x_i[n] \exp \left( \tilde{r}_i - \sum_j \tilde{\alpha}_{ij} x_j[n] \right), \qquad (S3)$$

where, for compactness, we have defined $\tilde{r}_i = r_i \Delta t$, and $\tilde{\alpha}_{ij} = \tilde{r}_i \alpha_{ij}$ ($i \neq j$), $\tilde{\alpha}_{ii} = \tilde{r}_i$. This 'logarithmic' scheme is equivalent to applying a standard Euler scheme to the evolution of log-abundances, $y_i(t) = \ln x_i(t)$. Indeed, any integration scheme applied to log-abundances will preserve positivity.

Another approach is based on the formal solution

$$x_i(t + \Delta t) = x_i(t) \cdot \frac{G_i(t + \Delta t | t)}{1 + r_i x_i(t) \int_t^{t+\Delta t} G_i(t' | t) dt'} \qquad G_i(t' | t) = e^{r_i \int_{t'}^{t} g_i(s) ds} \qquad (S4)$$

The difference to Eq. (S2) is that the right-hand-side does not explicitly depend on $x_i$. Any approximation of the integral in $G_i$ will preserve positivity of abundances. Choosing $G_i(t'|t) \approx \exp(r_i g_i(t)(t' - t))$ for $t' - t$ small, then performing the integral in the denominator of Eq. (S4), we obtain

$$x_i(t + \Delta t) \approx \frac{g_i(t)x_i(t)}{g_i(t)e^{-r_i g_i(t)\Delta t} + x_i(t)(1 - e^{-r_i g_i(t)\Delta t})} \qquad (S5)$$

The resulting 'logistic' integration scheme is thus (after some rearrangements)

$$x_i[n+1] = g_i[n] \times \frac{x_i[n]}{x_i[n] + (g_i[n] - x_i[n])e^{-r_i g_i[n] \Delta t}}, \quad g_i[n] = 1 - \sum_{j(\neq i)} \alpha_{ij} x_j[n]. \qquad (S6)$$

This scheme was proposed by Jules Fraboul [Fraboul2023].

To either integration scheme we can add a term $+\lambda\Delta t$ for the immigration.

## S3 Simplification of the disordered Lotka-Volterra model with mixing

We consider a well-mixed volume $V$ containing $S$ species with instantaneous absolute abundance $N_i(t)$, and constant nominal carrying capacities $K_i$ and growth rates $R_i$. The growth dynamics follows the standard Lotka-Volterra form. We add the effect of mixing with an external environment containing the same set of species but at abundances $N_i^{\text{ext}}(t)$: a fraction $\Lambda$ per unit time of the volume $V$ is exchanged with an equal volume from the external environment that gets instantaneously mixed in with the volume $V$. In total, the dynamics of the abundances is

$$\dot{N}_i(t) = R_i N_i(t) \left( 1 - \frac{N_i(t) + \sum_{j(\neq i)} \beta_{ij} N_j(t)}{K_i} \right) + \Lambda(N_i^{\text{ext}}(t) - N_i(t)). \qquad (S7)$$

Note that if we let $\Lambda \gg \max_i R_i$, we will force $N_i(t) \approx N_i^{\text{ext}}(t)$. Instead, we consider the slow-mixing scenario $\Lambda \ll \min_i R_i$, since our purpose in adding mixing is mainly to prevent extinction of rare species. We introduce
rescaled parameters

$$r_i = R_i - \Lambda \approx R_i, \qquad (S8)$$

$$K'_i = K_i \left(1 - \frac{\Lambda}{r_i}\right) \approx K_i, \qquad (S9)$$

$$\alpha_{ij} = \frac{\beta_{ij}K'_{j}}{K'_{i}}, \qquad (S10)$$

$$\lambda_i(t) = \Lambda \frac{N_i^{\text{ext}}(t)}{K'_i}, \qquad (S11)$$

$$x_i(t) = \frac{N_i(t)}{K'_i}, \qquad (S12)$$

so that

$$\dot{x}_i(t) = r_i x_i(t) \left( 1 - x_i(t) + \sum_{j(\neq i)} \alpha_{ij} x_j(t) \right) + \lambda_i(t). \qquad (S13)$$

Since the mixing occurs slowly, we suppose that it is justified to replace $N_i^{\text{ext}}$ with its time-average $\bar{N}_i^{\text{ext}}$. A parsimonious distribution for these abundances is that they result from effectively independent species (in particular in the chaotic phase; see main text) constrained by a roughly constant total community biomass $N^*$ independent of the number of species:

$$\bar{N}_i^{\text{ext}} = N^* \frac{K_i}{\sum_j K_j}. \qquad (S14)$$

Introducing $\tilde{\lambda} = \Lambda N^*/K^*$ with $K^* = S^{-1} \sum_j K_j$ we have

$$\lambda_i(t) \equiv \lambda = \frac{\tilde{\lambda}}{S}. \qquad (S15)$$

Upon the simplifying assumptions $r_i \equiv 1$ and Gaussian distribution of $\alpha_{ij}$ (resulting from a combination of some distribution of $K_i$ and $\beta_{ij}$; see [1]), we obtain the model of the main article.

## S4 Invasion analysis of few-species equilibria

Let $\mathcal{R}(t), \mathcal{D}(t)$ be the sets of rare and dominant species. The dynamics of a rare species $i \in \mathcal{R}(t)$ that successfully invades is

$$\dot{x}_i \approx x_i \left( 1 - \sum_{j \in \mathcal{D}} \alpha_{ij} x_j \right) := \gamma_i x_i, \qquad (S16)$$

while while $\lambda \ll x_i \ll 1$. The growth rate $\gamma_i$ must be (mostly) positive in the time interval where invasion occurs. Thus, no rare species is expected to be able to invade the dominant component while

$$\gamma_{\max} := \max_{i \in \mathcal{R}} \gamma_i < 0. \qquad (S17)$$

Suppose then that the $|\mathcal{D}| \approx \bar{S}_{\text{eff}}$ dominant species are in a few-species equilibrium that would be stable but for a potential invasion. We approximate the abundances of dominant species as equal (we expect them at least to be of the same order) and recall that they make up the overwhelming share of the total abundance. Thus, $x_j \approx \bar{X}/\bar{S}_{\text{eff}}$
for $j \in \mathcal{D}$. We suppose further that $\{\alpha_{ij}\}_{j \in \mathcal{D}}$ for a random $i \in \mathcal{R}$ can be treated as independent. Then the sum $\sum_{j \in \mathcal{D}} \alpha_{ij} \sim \mathcal{N}(\bar{S}_{\text{eff}}\mu, \bar{S}_{\text{eff}}\sigma^2)$. It follows that, for random $i \in \mathcal{R}$,

$$ \gamma_i \sim \mathcal{N} \left( 1 - \mu \bar{X}, \frac{\sigma^2 \bar{X}^2}{\bar{S}_{\text{eff}}} \right). \qquad (S18) $$

Using extreme value theory [2], and the fact that $|\mathcal{R}| \approx S$,

$$ \gamma_{\max} \approx 1 - \bar{X} \left( \mu - \frac{\sigma}{\sqrt{\bar{S}_{\text{eff}}}} h(S) \right), \qquad (S19) $$

where $h(S)$ is a random variable<fnref n="1" /> that scales approximately as $\sqrt{\ln S}$. In particular, $h(500) \approx 3.04 \pm 0.45$ and $h(10'000) \approx 3.85 \pm 0.35$. Using Main Text Eq. (11),

$$ \gamma_{\max} \approx \frac{1 - \mu - (\bar{S}_{\text{eff}}\bar{\rho} + \sqrt{\bar{S}_{\text{eff}}h(S)})\sigma}{1 + (\bar{S}_{\text{eff}} - 1)\mu - \bar{S}_{\text{eff}}\sigma\bar{\rho}}. \qquad (S20) $$

Condition Eq. (S17) for non-invadability then amounts to

$$ \sigma < \frac{\mu - 1}{\bar{S}_{\text{eff}}\bar{\rho} + \sqrt{\bar{S}_{\text{eff}}h(S)}}. \qquad (S21) $$

This predict lines radiating from $(\mu, \sigma) = (1, 0)$. However, with $\bar{\rho} \sim (\bar{S}_{\text{eff}})^{-1/2}$, the slope of the lines are less steep for sectors with more species in the equilibrium, contrary to observation.

## S5 Solution of intermittency model under unified coloured noise approximation

The unified coloured noise approximation was put forth by Jung & Hänggi [3] (see also Fox [4, 5]) to solve SDEs of the form

$$ \dot{x}(t) = F(x(t)) + G(x(t))\eta(t), \qquad (S22) $$

with $\eta(t)$ a coloured Gaussian noise: $\langle\eta\rangle = 0$, $\langle\eta(t)\eta(t')\rangle = \exp(-|t-t'|/\tau)$. By differentiating Eq. (S22), rearranging terms, and rescaling time into $\hat{t} = \tau^{-1/2}t$ (with $\hat{x}(\hat{t}) = x(\tau^{1/2}\hat{t})$) one obtains the exact equation

$$ \frac{d^2 \hat{x}}{d\hat{t}^2} + \frac{G'}{G} \left( \frac{d\hat{x}}{d\hat{t}} \right)^2 + H_{\tau} \frac{d\hat{x}}{d\hat{t}} = F + \sqrt{2\tau^{1/4}} G \hat{\xi}, \qquad (S23) $$

where $\hat{\xi}$ is a standard Gaussian white noise and

$$ H_{\tau}(x) = \tau^{-1/2} - \tau^{1/2} \left( F'(x) - \frac{G'(x)}{G(x)} F(x) \right). \qquad (S24) $$

As either $\tau \to \infty$ or $\tau \to 0$, one finds $H_{\tau} \to \infty$, and then Eq. (S23) can be replaced by the overdamped limit

$$ \frac{d\hat{x}}{d\hat{t}} = \frac{F}{H_{\tau}} + \sqrt{2\tau^{1/4}} \frac{G}{H_{\tau}} \hat{\xi}. \qquad (S25) $$

<footnote n="1">¹The maximum M of N i.i.d. standard Gaussian RVs tends, as $N \to \infty$, to $M = a_N + \xi/a_N$, where $a_N = \sqrt{2 \ln N - \ln(4\pi \ln N)}$ and $\xi$ follows a standard Gumbel distribution, whose mean is the Euler-Mascheroni constant ($\approx 0.577$) and variance is $\pi^2/6$.</footnote>
It is hoped that this approximation is accurate for intermediate $\tau$ as well.

The overdamped equation Eq. (S25) can be solved for its steady state $P^*(x)$ by conventional techniques [6], e.g. under Stratonovich interpretation of the noise: the associated stationary probability current is then

$$J^* = \frac{F}{H_{\tau}} P^* - \tau^{1/2} \frac{G}{H_{\tau}} \frac{d}{dx} \frac{G}{H_{\tau}} P^*. \qquad (S26)$$

In one dimension, $J^*(x)$ must be constant. Since $x$ is a non-negative abundance in our case, we must impose a boundary at $x = 0$ through which probability cannot flow. Therefore $J^* \equiv 0$. The solution for $P^*$ is then

$$P^*(x) \propto \exp\left\{\int^x v(x') dx'\right\}, \qquad (S27)$$

with

$$v = (\tau^{-1/2} H_{\tau}) \frac{F}{G^2} + \left(\ln \frac{H_{\tau}}{G}\right)' \qquad (S28)$$

In the stochastic intermittency model, we have

$$F(x) = -x(k + x) + \lambda, \qquad (S29)$$

$$G(x) = ux, \qquad (S30)$$

$$H_{\tau}(x) = \tau^{-1/2} + \tau^{1/2}(x + \lambda x^{-1}). \qquad (S31)$$

With these functions, the integral in Eq. (S27) can be performed exactly, yielding the result of the main article:

$$P^*(x) = \frac{1}{\mathcal{N}} e^{-[q_+(x)+q_-(\lambda/x)]} x^{-\nu} \left(\tau^{-1} + x + \frac{\lambda}{x}\right), \qquad (S32)$$

$$q_{\pm}(y) = \frac{1}{2u^2} [y + (\tau^{-1} \pm k)]^2, \qquad (S33)$$

$$\nu = 1 + \frac{k}{\tau u^2}. \qquad (S34)$$

The normalization constant $\mathcal{N}$, however, we do not have in closed form; instead we calculate it numerically, after the change-of-variables $y = \ln x$, as

$$\mathcal{N} = \int_{-\infty}^{\infty} dy e^{y(1-\nu)-q_{+}(e^y)-q_{-}(\lambda e^{-y})} (\tau^{-1} + e^y + \lambda e^{-y}). \qquad (S35)$$

## References

[1] M. Barbier and J.-F. Arnoldi, The cavity method for community ecology, bioRxiv:10.1101/147728, 10.1101/147728 (2017).

[2] P. Vivo, Large deviations of the maximum of independent and identically distributed random variables, Eur. J. Phys. 36, 055037 (2015).

[3] P. Jung and P. Hänggi, Dynamical systems: a unified colored-noise approximation, Phys. Rev. A 35, 4464 (1987).

[4] R. F. Fox, Functional-calculus approach to stochastic differential equations, Phys. Rev. A 33, 467 (1986).

[5] R. F. Fox, Uniform convergence to an effective Fokker-Planck equations for weakly colored noise, Phys. Rev. A 34, 4525 (1986).

[6] C. W. Gardiner, Stochastic methods, A handbook for the natural and social sciences (Springer, 2009).