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

Significance A prominent feature of ecological communities is that a few species are abundant while most are rare. Using a standard community model, in which species interactions are assigned fixed random values, we show that chaos is a generic outcome if interactions are strong and immigration prevents extinction. Each species then alternates, in an effectively stochastic way, between long periods of rarity and shorter periods of high abundance; yet the overall distribution of species abundances remains conserved and qualitatively consistent with observations in marine plankton protists. Our model results contribute to a rekindled debate about the role of chaos in ecological communities.


Introduction
Scientists have long marvelled at the complexity of ecosystems, and wondered what ecological processes allowed them to remain diverse despite species' competition for shared resources and space.The coexistence of a large number of di erent taxa reaches its climax in microbial communities, where several thousands of Operational Taxonomic Units (referred to as 'species' in the following), can be detected when sequencing samples of soil or water [1].A general feature of microbial communities is that abundances vary widely among the species that co-occur at one time of observation.Generally, just a handful of abundant species make up most of the total community biomass [2][3][4], whereas the large majority of rare species are present in such low numbers as to only be detectable by su ciently deep genomic sequencing.
Multiple hypotheses have been put forward to explain the nature and origin of such a 'rare biosphere' [2,5].One possibility is that distinct taxa belong preferentially to the abundant sub-community or to the rare, either because of their traits, or because evolutions drives new, more adapted species to dominance [6,7].The notion that abundant and rare sub-communities are well-distinct is supported by the observation that their macroecological patterns appear to di er [8,9].Alternatively, diversity in community composition could be maintained by variations in the environment that allow for many spatio-temporal niches [10].The concept of a microbial seed bank encapsulates the idea that a small number of dominant taxa are maintained by environmental ltering, while most taxa remain dormant until they meet conditions favourable to their growth [11].In this view, the abundance or rarity of a given species is a function of the extrinsically-driven conditions at the time of sampling.A third hypothesis is that cycling of species between the rare and abundant component of the community is driven by intrinsic ecological uctuations that self-sustain even in the absence of environmental variation.Such oscillations have been found in controlled settings [12,13], and have been argued to be relevant also in natural communities [14][15][16].For instance, planktonic bacteria display a fast turnover of species within a season, even as abiotic conditions do not vary substantially [17][18][19].This supports the notion that species interactions may be central in determining which species are abundant, and when.
That ecological interactions can cause instability of species' abundances and drive chaotic uctuations is predicted by numerous theoretical models.Chaos can be found in models with just a handful of species [14,20,21], but may then require a ne-tuning of parameters [22].In contrast, instability and chaos appear to be generic features of high-dimensional systems, such as species-rich communities [23,24].The disordered Lotka-Volterra equations, for instance, representing intra-speci c and randomly assigned inter-speci c interactions have been broadly used to explore the collective ecological dynamics of species-rich communities [25][26][27][28][29][30][31][32].
Self-sustained oscillations that resemble the turnover of natural communities have been obtained by numerical simulations of species-rich meta-communities, where local abundances are coupled via dispersal [33][34][35][36].In this setting, species that are locally driven to extinction get replenished by migration of individuals from neighbouring patches, where they were not outcompeted.Among the key features of natural ecosystems that such models reproduce are large numbers of coexisting taxa, large uctuations in abundance, and skewed species abundance distributions with prominent power-law trends.However, given the complexity of the models, it seems inevitable that analytical descriptions are only possible under speci c assumptions, for instance on the scaling or (anti-)symmetry of interactions [25,32,33].
On the other hand, more phenomenological models, that assume from the outset the existence of intermittent population dynamics, also predict highly diverse communities where species alternate between rarity and abundance [18,37].Similarly, single-species stochastic equations account for many features of empirical abundance distributions [38,39].When species are also subject to evolution, booms and busts allow to maintain a larger number of species than expected by ecological dynamics alone [40].
Given the range of alternative theoretical descriptions, it is thus unclear what level of complexity-in terms of species richness, type and variation of interactions, spatial structure-is necessary to explain non-stationary patterns of highly diverse communities.
Here, we show that abundance patterns consistent with observations of plankton communities emerge generically in a disordered Lotka-Volterra model under simple assumptions: a single patch with well-mixed populations, small and constant immigration, and strong inter-species interactions.Under such conditions, we nd a broad range of model parameters where species alternate over time between rarity and abundance, in such a fashion that the community is at any time dominated by just a few abundant species.The turnover in composition of this dominant component occurs on a characteristic timescale, and resembles a succession of low-diversity equilibria.We compare the distribution of abundances observed in the community at a xed time with the frequency with which the abundance of any given species occurs in a long time series.These distributions all have the same shape, which across many orders of magnitude in abundance values is a power-law, suggesting an emergent equivalence among di erent taxa [6,33,[41][42][43].We therefore propose an approximate, e ective model for the ecological dynamics of a 'typical' focal species.Guided by this model, we characterize the region where the ecological dynamics is chaotic and point to scaling relations that may explain the weak geographical variation of plankton protist communities' abundance distributions [3].We also explore the reason why some species deviate from the typical abundance pattern.Species that have a smaller average interaction strength with all other species boom more frequently, highlighting the importance of relative interaction statistics in addition to absolute ones.

Model
We describe a community of species by their timedependent absolute abundances ( ), with = 1, 2, … , the index of a species.Deterministic equations that relate the changes in abundance to competition within species and interactions between pairs of species have been argued to be relevant descriptions of diverse microbial communities [13,44,45].According to the Lotka-Volterra equations [46], the abundance of any species in isolation grows logistically: if initially the species is rare, its abundance grows exponentially at a maximum rate ; later, it saturates to a carrying capacity set by resources, predators, and environmental conditions assumed constant and not modelled explicitly.For simplicity, we set and to unity for all species, but discuss heterogeneity in these parameters in Supplementary Note S3.The interaction coe cients (real numbers) quantify the e ect of species on the growth rate of species .We include a small rate of immigration ≪ 1 into the community; constant and equal for each species.This term prevents extinctions, and re ects immigration from a regional pool or the existence of as a 'seed bank' [11].Abundances thus change in time as In species-rich communities, the number of potential interactions-× -is very large, and their values hard to estimate in natural settings.Therefore, a classic approach [23,25,27,47] is to choose the set of interaction coe cients as a realization of a random interaction matrix , whose elements are Gaussian random variables; ∼ ( , 2 ) ( ≠ ).It is customary to allow a correlation between diagonally opposed elements, biasing interactions toward predator-prey ( = −1) or symmetric competition ( = 1); here, we focus on independent interaction coe cients ( = 0) and discuss other cases in Supplementary Figure S5.According to Eq. ( 1), positive denotes that species reduces the growth of species , as for instance when they compete for a common resource.A negative value indicates that species facilitates the growth of species .The interaction coe cients for distinct species , can be represented in terms of the mean and standard deviation of the interaction matrix, as where the 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 intraspeci c interaction terms in Eq. ( 1).The diagonal element therefore does not appear in the sum, and is not de ned.Equation ( 1) with randomly sampled interactions de nes the disordered Lotka-Volterra (dLV) model.By tuning the ecological parameters , , , , 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 where the net interaction bias is a realization of a random variable ∼ (0, 1).To achieve a nite net competition in the limit of a large species pool requires where ̃ , ̃ do not grow with .Under this scaling, methods from statistical physics (dynamical mean-eld theory [25,26,28,30], random matrix theory [23,48], and replica theory [29,31]) allow exact analytical results in the limit of → ∞, although in practice ∼ 100 is su cient 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 [25,26,28,30].
While vanishing interactions entail signi cant mathematical convenience, one can question how well the weak-interaction regime represents microbial communities.For instance, bacterial species engage in metabolic cross-feeding, toxin release, phagotrophy, and competition over limited nutrients, so that some species do depend substantially on another's presence [49,50].Moreover, some empirical species abundance distributions-notably those of plankton communities [3,4]-deviate qualitatively from those predicted for weak interactions [25,26].Finally, there is increasing evidence that both absolute and relative abundances of microbial species display large and frequent variation even on time scales of days, where environmental conditions are not expected to undergo dramatic variations [16][17][18].Weak interactions, instead, produce moderate uctuations, so that the total abundance can be assumed to be constant [30].For these reasons, as well as for completing the analysis of disordered Lotka-Volterra model, we consider here the strong-interaction regime where the statistics of the interaction matrix do not scale with species richness according to Eq. ( 5).For ≫ 1, the overall competitive pressure makes it impossible for all species to simultaneously attain abundances close to their carrying capacities, resulting in instability and complex community dynamics.

Results
In the strong-interaction regime, numerical simulations of the disordered Lotka-Volterra model show that the community can display several di erent 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 ( = 0.5, = 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 and .Unless otherwise stated, simulations use = 500 and = 10 −8 .Further details on the numerical implementation are presented in Appendix A.

A chaotic turnover of rare and abundant species
For a broad range of parameters in strong-interaction regime, the community undergoes a chaotic turnover of dominant species.As illustrated by the time series of stacked dances 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 , 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.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, dom ≈ 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 (− ln ≈ 18).While the total abundance uctuates moderately around a well-de ned time average, individual species follow a 'boom-bust' dynamics.If this simulation represented a natural community, only the most abundant species-that 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 de ne the e ective size of the community as Simpson's (reciprocal) diversity index [51], where = ∕ ∑ denote relative abundances.e approaches its lowest possible value of 1 when a single species is responsible for most of the total abundance, and its maximum when all species have similar abundances.Its integer approximation provides the richness, i.e. number of distinct species, of the dominant component.
The e ective size e of the community in our reference simulation uctuates 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 uctuates around 3%, which is comparable to the arbitrary 1%-threshold used in empirical studies [9].In Supplementary Figure S4 we show that the number of Chaotic turnover ... in a strongly interacting model community Mallmin, Traulsen, De Monte (2023) dominant species grows slowly (but super-logarithmically) with , up to about 15 for = 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 suggest the presence of chaotic dynamics.We give numerical evidence for sensitive dependence on initial condition and positive maximal Lyapunov exponent in Supplementary Figure S1 and S2.The turnover dynamics has the character of moving, chaotically, between di erent quasi-equilibria corresponding to di erent compositions of the dominant community (cf.'chaotic itinerancy' [52]).To reveal this pattern, we measure a 'closenessto-equilibrium', de ned 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 [17].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 a previously rare invader.
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 has been con rmed [53,54].It is therefore conceivable that the number of quasi-equilibria in our case is also exponentially large.The LV equations for = 0 admit up to one coexistence xed point (not necessarily stable) for every chosen subset of species [46].Hence, we expect on the order of ∼ e quasi-equilibria, which for e ≈ 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 di erent 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 ( + )∕2 is significantly smaller than the expectation , and hence species and are close to a commensal or mutualistic relationship, these species tend to 'boom' one after the other (Supplementary Figure S6).

Species' abundance uctuations follow a power-law
In a common representation of empirical observations, where relative abundances are ranked in descending order (a rank-abundance plot [55]), microbial communities display an overwhelming majority of low-abundance species [5].Our simulated community reproduces this feature; Fig- with ± 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 uctuations.The one bar missing from the SAD is the e ect of nite species richness, as high-abundance bins only ever contain a couple of species for = 500.The vertical dashed line indicates the immigration level which determines a lower limit to abundances.ure 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) [55].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 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 uctuation distribution (AFD) [38].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, every species uctuates in time between extreme rarity ( ≈ = 10 −8 ) and high abundance ( ≳ 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 uctuations in time, as indicated by the small variation in AFDs across species.We will evaluate the regularities and di erences of single-species dynamics more thoroughly in  7)) with parameters as in Eq. ( 9): the time series are statistically similar.B Comparison of the average abundance uctuation distribution (AFD) from Figure 2 (black), and the AFD of the focal-species model (pink): excellent agreement is found for the powerlaw section.The 'uni ed 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.
The most striking feature of these distributions, however, is the power-law − 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 ≈ 1.18 for the simulation analysed, 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 goal of the next section.

A stochastic focal-species model reproduces boom-bust dynamics
Fluctuating abundance time series are often tted by onedimensional stochastic models [15]; for example, stochastic logistic growth has been found to capture the statistics of uctuations in a variety of data sets on microbial abundances [38,39].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 xed, there is no uncertainty; nonetheless, the chaotic, high-dimensional dynamics results in species' growth rates uctuating 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 re-ect thee e ective dynamics of any species in the community.
Following dynamical mean-eld-like arguments and approximations informed by our simulations (Appendix E), we derive the focal-species model where ( ) is a stochastic growth rate with mean − , and uctuations of magnitude and correlation time .The process ( ) is a coloured Gaussian noise with zero mean and an autocorrelation that decays exponentially; where brackets denote averages over noise realizations.The connection between the ecological parameters , , , and the resulting dynamics of the disordered Lotka-Volterra model in the chaotic phase is then broken down into two steps: how the e ective parameters , , relate to the ecological parameters; and how the behaviour of the focalspecies model depends on the e ective parameters.
For the rst step we nd where is the total community abundance of the original dynamics Eq. ( 1), the e ective community size e is as in Eq. ( 6), and an overline denotes a long-time average.Equation (9) relates the focal species' growth rate to the timeaveraged net competition (≈ ) from all other species.We nd in simulations of Eq. ( 1) in the chaotic phase that competition is strong enough to make > 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 ( e ) −1∕2 .The third e ective parameter, the timescale , controls how long the focal species stays dominant, once a uctuation has brought it to high abundance.This timescale is essentially equal to the turnover timescale dom of the dominant component (de ned more precisely by autocorrelation functions in Appendix E).In the weak-interaction regime, where any pair of species can be treated as e ectively independent at all times, self-consistency relations such as ⟨ ⟩ = allow to implicitly express the focal-species model in terms of the ecological parameters.For strong interactions, however, the disproportionate e ect of the few dominant species on the whole community invalidates this approach; we therefore relate the e ective parameters to the community-level observables , e , 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 e ective parameters.While this is intractable due to the nite correlation time of the noise, the equations can be simulated and treated by approximate analytical techniques.In Figure 3A we compare the time series A Example of a long abundance time series for the three species who are ranked rst, 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 results in a proportional decrease in median dominance time.C Distribution of dominance biases against relative dominance rank for a range of : there appears to be convergence towards a non-constant limiting distribution, implying that net species di erences are not due to small-e ects.Note that, by de nition, 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 normalized sum of interaction coe cients, Eq. ( 4): lower net competition correlates with higher dominance bias.Species in the tails of the distribution are also less 'typical', with typicality quanti ed by the index , Eq. ( 16), representing the similarity of a species AFD to the species-averaged AFD.Panel A and D are both for = 500.
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 uctuating growth rate in in Eq. ( 7) with its typical value (i.e.= 0), yielding the equilibrium ∕ if > 0, as indeed con rmed 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 uni ed coloured noise approximation [56] (Appendix F), one predicts that the stationary distribution, for ≪ ≪ 1, takes the power-law form − , where the exponent is strictly larger than one-the value predicted for weak interactions [30] and for neutral models [57].Even if Eq. ( 10) is not quantitatively precise (Figure 3B), this formula suggests a scaling with the e ective parameters that we will discuss later on.

Species with lower net competition are more often dominant
The similarity of all species' abundance uctuation distributions in Figure 2 is re ected 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 uctuate over orders of magnitude, some species are observed to be more often dominant (or rare).Such di erences are reminiscent of the distinc-tion between 'frequent' and 'occasional' species observed in empirical time series [58,59].
In order to assess the nature of species di erences 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 (400 times longer than in Figure 1), the rst-ranked species appears to boom much more often than the last (Figure 4A).The frequency of a species is chie y 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 leads to each species halving its dominance time fraction.As the community gets crowded-while its e ective size hardly increases, as remarked in Section 3.1-all species become temporally more constrained in their capacity to boom.Yet some signi cant 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 speciesagainst the relative rank (i.e., rank divided by ).For high richness ( ∼ 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 di erences with large may seem to contradict the central limit theorem, as the vectors of the interaction coe cients converge towards statistics that are identical for every species.In the chaotic regime, however, even the smallest di erences in growth rates get ampli ed during a boom.As we show in Appendix D, if Eq. ( 1) is rewritten in terms of the proportions , the relative advantage of species is quanti ed by a selection  2∕ ) 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 multipleattractor phase (in which cycles and multi-stable xed point are common in addition to chaos) is not sharp, unless probed adiabatically in the way explained in Supplementary Figure S8.The unique xed-point phase has been studied analytically in the weak-interaction regime ( ∼ 1∕ ).When inter-speci c competition is in general stronger than intra-speci c competition, a single species (identity depending on initial condition) dominates, in line with the classical competitive exclusion principle [60].coe cient whose time average scales as − −1∕2 .Correspondingly, the relative, time-averaged growth rate is proportional to the net interaction bias (de ned in Eq. ( 4)), resulting in species with larger 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 quanti ed by the typicality index ∈ [0, 1], de ned in Appendix B.
In conclusion, the relative species-to-species variation in the total interaction strength drives the long-term di erences in the dynamics of single species in the community.While the focal-species model emphasizes the similarity of species, species di erences can also be taken into account by employing species-speci c e ective parameters.In particular, replacing with a distribution of '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).

Interaction statistics control di erent dynamical phases
Hitherto, we have focussed on reference values of the interaction statistics and 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 ( , ) values, we run 30 independent simulations, each with a di erent 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 di erent classes: equilibrium, cycle, chaos, or divergence.Figure 5 displays the probability of observing chaos, demonstrating that it does not require ne 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 di erent community dynamics.For small variation in interaction strengths (below the line connecting (−1∕ , √ 2∕ ) to (1, 0)), the community has a unique, global equilibrium that is fully characterized for weak interactions (cf.Fig. 2. of [25]).The transition from equilibrium to chaos has been investigated with dynamical mean-eld theory [30].For low interaction variance, but with mean exceeding the unitary strength of intra-speci c competition, a single species comes to dominate, as expected by the competitive exclusion principle [60].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 ( , ) = (1, 0) separate sectors where stable xed points have di erent numbers of coexisting species.Traversing these sectors anti-clockwise, e 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 di erent types: cycles and chaos, coexisting with multiple xed points, resulting in hysteresis (Supplementary Figure S8B).This 'multiple attractor phase' [25,30] 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 and e ective community size e vary considerably (Supplementary Figure S9).The weak-interaction regime (whether in the equilibrium or chaotic phase) allows for high diversity, so and e are of order ; strong interactions, on the other hand, imply low diversity, with e and of order unity.An explicit expression for how these community-level observables depend on the ecological parameters ( , , , ) is intractable (although implicit formulas exists in the weak-interaction regime [25]).Nonetheless, an approximate formula that we derive in Appendix C allows to relate community-level observables to one another and to and : in which we introduce the collective correlation involving the time-averaged product of relative abundances weighted by the their normalized interaction coe cient Chaotic turnover ... in a strongly interacting model community Mallmin, Traulsen, De Monte (2023) 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 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 e ective parameter = −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 > 0 with Eq. ( 11), we estimate a minimum, critical value of the collective correlation required for boom-bust dynamics: Numerical simulations demonstrate that ≳ 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 ≳ 1∕ in the chaotic phase.For strong interactions, total abundances are predicted to be of order one, and for weak interactions ≈ ∕ ̃ (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 e ective community size drops, suggesting a seamless transition from a weak-interaction, chaotic regime amenable to exact treatment [30], to the strongly correlated regime that we have analyzed by simulations and the approximate focalspecies model.

Self-organization between community-level observables constrains abundance power-law variation
In Section 3.3 we established a focal-species model depending on the e ective parameters , , and , that were related to the ecological parameters , , , indirectly via community-level observables , e , 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 e ective parameters and of the focal-species model's predictions.Because the trio , , ultimately hails from but two independent variables, , (considering xed , ), they must be dependent.Figure 6A demonstrates that, across the chaotic phase, an approximate linear relationship holds between and , as well as between and .Because and 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 to is close to one (and varying little with and ; Supplementary Figure S10), which implies with Eq. ( 9) that Comparison to Eq. ( 11) then yields that − ≈ e −1∕2 .This empirical relationship thus supports the aforementioned convergence-in the limit where e is large, as for weak interactions-of the collective correlation to its critical value.
We nd in Figure 6 that the slope foc of the power-law trend obtained from simulation of the focal-species model nds good agreement with the value from the full dLV model.There is a narrow overall variation of the exponent; a consequence of the interdependency of the e ective 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 ( , ) = (1, 0).However, the exponent also depends on and , showing a constant slope against log or −1∕ log (Supplementary Figure S7).

Discussion
We have sought a possible theoretical underpinning for macroecological patterns of dominance and rarity in species-rich communities.To this end, we studied a Lotka-Volterra model with strong interactions and weak immigration using numerical simulations and approximate analytical techniques.We characterized a parameter regime where species generically turn over between a small dominant component, and a large pool of temporarily rare species.In this process, each species' abundance undergoes a chaotic boom-bust dynamics, asynchronous with respect to most other species.The resulting distribution of abundances-of a single species over long times, or of the whole community at a single time-has a prominent power-law trend.
The phenomenology of the model-chaos, boom-bust dynamics, and a power-law shaped SAD with exponent larger than one-is consistent with observations of marine plankton communities [3,12,16,18,62].While the evidence for chaos in ecological time series has been generally ambiguous, a recent systematic assessment concludes that chaos is commonplace, especially for plankton [16].Experiments with closed plankton microcosms have revealed chaotic, high-amplitude uctuations sustained over many years [12,62].Abundance uctuations indicative of chaos were also seen in non-planktonic synthetic microbial communities [13,63], but of lower amplitude.That planktonic population sizes uctuate over many orders of magnitude in abundance is made especially poignant by algal blooms, which can become visible even from space.The timing of blooms, and the succession of functional groups within a season, are coupled to environmental factors such as nutrient concentrations.Yet, for the non-dominant taxa, the large di erences in a species' abundance between ocean samples show  5).Each pair of ( , ) has been mapped to a distinct colour.B The exponent of the power-law section of the AFD for the chaotic dLV model plotted against the analogue 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.C Co-dependence of the e ective parameters , , : the amplitude of growth-rate uctuations approximately equals the absolute value of the negative growth rate (only weakly depending on and ; Supplementary Figure S10); is roughly proportional to the inverse turnover time, but the slope of the relationship depends on and .little environmental signature [3].This suggests that the turnover of abundances might rather be driven by complex interactions or mixing dynamics.Empirical snapshot SADs of marine protists show a clear power-law trend for nondominant species within the same size class, with a largerthan-one exponent varying little between samples, as in our model.
The empirical value of the exponent of the SAD's powerlaw trend (around 1.6 [3]) is important because it rules out particular model assumptions.The chaotic phase in the Lotka-Volterra model produces a unitary exponent in the weak-interaction limit, and seemingly also for strong interactions if the immigration rate vanishes, → 0. Similarly, neutral theory predicts a power-law tail of the SAD with exponent one [57,64].To approach the empirical value, previous studies augmented neutral theory with nonlinear growth rates or chaotic mixing [3,65] to nd an exponent dependent on the model parameters.For dLV with strong interactions, we have shown that > 1 for all parameter combinations within the chaotic phase when interactions are strong.The approximate solution to the focalspecies model, Eq. ( 10), shows that the positive deviation from = 1 depends on three inter-related e ective parameters: the mean, amplitude, and timescale of uctuations in each species' net competition.As these uctuations drive the turnover pattern, boom-bust dynamics comes to be associated to a larger-than-one exponent.In our model, approaches the empirical value in the green region of Figure 6A.There, interactions are strong (smaller than but relatively close to self-interactions) and vary moderately from species to species, and the turnover time scale is large (but not yet divergent).We can speculate that similar features underpin plankton ecological dynamics, and may di erentiate these communities from other, more stable microbial assemblages.
The chaotic phase ends at the parameter point ( , ) = (1, 0) where all non-divergent dynamical phases meet (Figure 5).These phases appear to be qualitatively similar to those observed in an individual-based version of the dLV model, accounting for demographic stochasticity [66].
There, the phases meet at the "Hubble point" which recovers neutral theory; every pair of individuals compete equally, and species only come to di er through randomness in birth-death events.Despite being deterministic, our model shows some similarity to neutral theory (sometimes referred to as e ective neutrality [43,67,68]): species are largely equivalent in a statistical sense, and they uctuate (pseudo-)randomly and mostly independently.It is important, however, to highlight the di erences that allow distinguish between models.As argued before [33], the vast census sizes of planktonic species would imply enormous timescales of turnover if by demographic noise alone, the main driver of diversity in neutral models.In contrast, deterministic species interactions can produce rapid turnover.If the turnover time is long, or if there is no turnover, then species-speci c AFDs would di er substantially from one another, and not span the range of abundances observed across the whole community.It is therefore informative to measure AFDs in addition the more commonly measured SADs, whenever this is feasible.
Fluctuations, whether due to demographic or environmental noise, or complex interactions, can drive species extinctions.In our model, the immigration term represents some mechanism able to sustain abundances above the extinction threshold long enough for a species to rebound.Meta-community models of spatial patches connected by migration ows show how drastic levels of global extinctions can be avoided through a temporal storage effect [33][34][35][36]: if the patches' dynamics desynchronize, then a species that went extinct in one patch may eventually be reestablished there through immigration from another patch where it persisted.Within-patch turnover of composition was observed in the meta-community setting but attributed to the local interactions in a single patch [33,36].In particular, Ref. [33] focussed on disordered Lotka-Volterra dynamics where both within-and between-species interactions follow the same statistics.This required anti-symmetric, i.e. predator-prey-paired, interactions ( = −1) for sustained uctuations.We have instead followed a competitive paradigm (normalized self-interaction, and = 0), di-Chaotic turnover ... in a strongly interacting model community Mallmin, Traulsen, De Monte (2023) rectly extending the studied range of the -phase diagram of earlier work [25,30].Furthermore, our model assumptions have allowed the formulation of an explicit focal-species model, in terms of a stochastic logistic equation with coloured noise.It is similar to models that have been successfully tted to microbial time series [38,39], but a notable di erence lies in the negative mean growth rate we nd, which together with noisecorrelation and immigration yields uctuations over many orders of magnitude.Our focal-species model constitutes an ecologically motivated candidate for tting to ecological time-series with large uctuations.Its derivation also serves an important conceptual purpose, as it illustrates in an ecological context how complex dynamics can come to resemble a simple noise process.Indeed, it is notoriously di cult to distinguish random noise from chaos in empirical time series.The prevalence of ecological chaos may have been underestimated for this reason [15,16].
Chaos is a multifaceted phenomenon, and the question of how and why it arises in a given context is of great theoretical interest.It has long been recognized that Lotka-Volterra systems can admit heteroclinic networks [69][70][71]; saddle-points, i.e. equilibria with stable and unstable directions, connected by orbits.For LV equations without immigration, such saddles are found on the system boundary, corresponding to some subset of species having zero abundance (i.e.being extinct).One route to chaos is via the 'deformation' of a heteroclinic network, when, for instance, the introduction of an immigration term pushes saddles o the boundary.This can result in 'chaotic intinerancy', whereby trajectories traverse the vicinities of low-dimensional quasiattractors (formerly saddles on the boundary) via higherdimensional, chaotic phase-space regions [52,72].This picture ts well with our description of the chaotic turnover in relation to Figure 1.A further understanding of this mechanism may come from analytical investigations of disordered models related to, but more tractable than, Lotka-Volterra [73]; or focussing on the transition to chaos from the unique equilibrium regime in the weak-interaction limit of the dLV model.In the strong-interaction regime, a systematic bifurcation analysis of the transition from stable xed points to chaos under the adiabatic scheme outlines in Supplementary Figure S8 may also provide insights.
Our model can be extended in several directions.Differences in growth rates and carrying capacities between species are expected for natural communities, and revealed when tting time series [74]; one may also consider a less than fully connected interaction network, e.g.sparsity [75][76][77].Such generalizations will however break assumptions-in particular that of independent Gaussian interactions-underlying some of our analytical approximations.Preliminary numerical explorations suggest that existence of chaotic turnover dynamics is robust to these features, but that they may bias some species towards abundance or rarity (Supplementary Figure S5).A systematic investigation is however warranted.Given the proposed connection of our work to plankton ecology, considering a trophically structured ecosystem might be a particularly relevant generalization, as grazers [78] and viruses [79] have been put forward as key ecological actors.On the one hand, we have shown that non-structured competition between many species can produce high-dimensional chaos; on the other hand, the dynamics between functional groups such as bacteria, phyto-and zooplankton, and detrivores could potentially also be chaotic, but in low dimension.It is therefore an interesting question how uctuations across di erent levels of coarse-graining of a community might be intertwined.
To conclude, we have demonstrated the emergence of a chaotic turnover of rare and abundant species in a strongly interacting community with minimal model assumptions.By deriving an explicit focal-species model to capture this complex dynamics, we have identi ed community-level observables and e ective parameters that constrain the variation of the power-law exponent of the species abundance distribution.These insights may prove valuable for interpreting eld data [3,18], as well as for predicting dynamical features of synthetic communities [13].

Appendices A Numerical implementation
For Lotka-Volterra simulations we used a xed time-step Euler scheme with ∆ = 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 xed-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.

B Similarity metrics
The Bray-Curtis similarity index [80] is de ned as where is the relative abundance of species with respect to the joined abundances + .By de nition, BC( , ) = 1 i = , and BC ≈ 0 when, for each , either ≫ or ≫ ; this makes it suitable for communities where abundances span orders of magnitude.
For the similarity graph Figure 1B, we have plotted BC( dom ( ), * ( )), where dom ( ) is the restriction of ( ) in the reference simulation to only the dominant species at time , and * ( ) is the xed point reached from dom ( ) as initial condition, with = 0.
To compare the similarity the AFD of species , ( ), to the species-averaged AFD = ∑ ∕ , we de ne the index where and are the cumulative distribution functions of and , respectively; i.e., the index is based on the Kolmogorov-Smirnov distance [51] of the AFDs.

C Derivation of time-averaged total abundance
Direct summation of Eq. ( 1) over (assuming = 1), and then division on both sides by ( ) = ∑ ( ), yields ) with e as in Eq. ( 6) and ( ) as Eq. ( 12) but without the time average.We denote the time-average operator by Applying it to Eq. ( 17), the left-hand side becomes lim →∞ ( ( ) − (0))∕ , which evaluates to zero on the assumption that no species diverges in abundance.The righthand side contains terms such as [ ∕ e ] and [ ].If the relative uctuations in , e , are small (see Supplementary Figure S9), or these functions are at most weakly correlated to one another, then we obtain, approximately, As the immigration is small compared to the other terms it can be neglected; solving for nds Eq. ( 11).The relative error in between Eq. ( 11) for simulated values of the community-level observables in the right-hand side, and the simulated value of , is typically less than a few percent (see Supplementary Figure S14).

D Selective advantage
The dynamics of the relative abundance = ∕ , with = ∑ , is found by summing and di erentiating Eq. ( 1) as Using Eqs. ( 3), (6), and ( 12) in de ning we can write Eq. ( 20) as The term is responsible for the bias of species against the reference proportion .As a heuristic means of calculating the time-averaged bias, we suppose the 's can be treated independently of the and be replaces by ≈ 1∕ ; then we obtain ≈ − ∕ √ .On this basis, we expect to be indicative of a species' dominance bias.

E Derivation of the stochastic focal-species model from dynamical mean-eld arguments
We write Eq. ( 1) as If we suppose that the abundances { ( )} (or, rather, their statistical properties) are independent of the particular realization [ ] of the interaction matrix, then, for a given realization of { ( )}, based on the properties of sums of Gaussian variables.The time-varying mean and variance of means that, averaged over time, does not necessarily follow a Gaussian distribution.We introduce which are found to exhibit signi cant relative uctuations, with skewed distributions (Figure 7A and B).However, once we shift and scale ( ) into the "e ective noise" we recover (closely) a (0, 1) distribution, for both the set { ( )} 1,…, at any given time , or for the stationary distribution of ( ), at least for typical species (Figure 7C and  D).The empirical distribution of the across all species and times is closely approximated by the stationary distribution ( , ) (Figure 7E).Therefore, we suppose that, despite their uctuations, we can replace ( ) and ( ) with their time-averages and model as a stochastic process where ( ) is a process with stationary distribution (0, 1).The parameter correspondence in Eq. ( 9) follows by = − , = ≈ ∕ √ e , and = , the correlation time of .Note that, up to neglecting a diagonal term of the sum, the e ective noise can be written with ∼ (0, 1), and ( ) = ( )∕ ( ) 2 .Given the chaotic turnover pattern, the latter is expected to perform something like a random walk on the -sphere, with a decorrelation time corresponding to the turnover of dominant species.This timescale is inherited by the e ective noise.More precisely, we compare autocorrelation functions (ACF).The ACF of a function is de ned as species and time (grey), over just species for one random time (green), over all time for the rst/mid/last-ranked species with respect to average abundance (blue/pink/yellow), with (0, 1) (black, dashed) for reference.E The empirical distribution of in Eq. ( 23) over all species and times, compared to the distribution ( , ) assumed for in the focal-species model.F Autocorrelation functions: for every species (grey), rst/mid/last-rank species (blue/pink/yellow)), and the average over all species (black).The left inset compares the ACFs of (green), (black), and the exponential t to the latter (red); the right inset shows the distribution of the parameter in exponential ts to each species ACF. with = − and using the notation Eq. (18).By denition (0) = 1.For each species' e ective noise we compute numerically ( 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 The decay of correlation is well-approximated by the exponential exp(− lag ∕ ), where the parameter ( tted by least squares) represents the noise correlation timescale for a 'typical' species.
The approximately (0, 1) distribution and exponential autocorrelation function of the e ective noise suggest that it can be modelled as an Ornstein-Uhlenbeck process, the only Markov process with these two properties; where ( ) is a Gaussian white noise; ⟨ ⟩ = 0, ⟨ ( ) ( )⟩ = ( − ).The timescale referred to as dom in the main text can be de ned as , the decay time of the exponential t to the ACF of the abundance vector.For a vector-valued function, Eq. ( 29) gives Comparing and , they match very well (inset of Figure 7F) for the reference simulation; as do the associated timescales and for all ( , ) in the chaotic phase (Supplementary Figure S13).This observation motivates identifying of the focal-species model with the turnover timescale dom .Thus, the focal species model and its parameters have been fully speci ed.
The crucial di erence to dynamical mean-eld theory developed in the weak-interaction limit is the fact that ( ) and ( ) are, for strong interactions, determined by a small number ( e ) of dominant species, whose abundances uctuate substantially, and, during their time of co-dominance, have signi cant e ects on each other; i.e., they are conditionally correlated.Therefore, self-consistent determination of the e ective parameters fails, because species cannot be treated as independent realizations of the focal-species model.For example, the self-consistency relation for in Eq. (7b) is = ⟨ ⟩ − 1.In our reference simulation = 0.26, whereas ⟨ ⟩ − 1 = −0.31even has the wrong sign.This discrepancy is due to the neglected inter-species correlations needed for the collective correlation , Eq. ( 12), to exceed the critical value Eq. ( 13) associated with > 0 and boom-bust dynamics.

F Steady-state solution of the focal-species model under the uni ed coloured noise approximation
The uni ed coloured noise approximation [56] assumes overdamped dynamics to replace a process ̇ = ( ) + ( ) , driven by Gaussian correlated noise of correlationtime , with a process driven by white noise.The approximation is exact in the limits → 0 or → ∞.The stationary distribution of the corresponding white-noise process is * ( ) ∝ exp ∫ ( ) d , with and a function of , , and .For Eq. ( 7)

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

S1 Supplementary Figures
As in the    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 ≈ ≈ 8000 re ects a time when some dominant component persisted for an unusually long time.B For a few well-separated time points (one graph each), we show how ( , ) decays over time on a timescale of 200 time units (top panel), and how it uctuates 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 re ecting that one or a few species will eventually reappear in a dominant community that is otherwise di erently composed.6)), increases slowly (but super-logarithmically) with the overall richness .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.where is a single, xed realization of a standard Gaussian random matrix.A For each value of , we initialized separate simulation runs starting at = 1.4,and let their abundances evolve until an attractor was found.For each run, we then changed by small increments = −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 di erent attractor.If a state diverged, the initial abundances for the next value of were set as the most recent non-divergent attractor.Thus, each simulation traced a sequence of attractors from = 1.4 → −0.1, corresponding to a horizontal line in the phase diagram.The colour quality re ects the class of the attractor, and the colour gradation indicates the e ective community size, revealing the following features: First, we nd mostly xed points in the multiple attractor region.This is because, once a xed point is converged to, it is "hold on to" until it vanishes or changes stability.If, instead, every simulation at given , would start from newly sampled initial abundances and interaction matrix, we would nd di erent attractors every time, and the diagram becomes more heterogeneous (compare Main Text Figure 5).Second, clear lines radiate from ( , ) = (1, 0) and delineate sectors characterized by the number of high-abundance species coexisting at a xedpoint.In section S4 we show that an invasion analysis predicts such sectors, but not the right scaling of the lines' slope with e .Third, the jump from xed-point to chaotic attractors occurs along a sharply de ned line.B Stacked abundances of the attractor found in an adiabatic sequence = 1.4 → 0.6 (top panel, right to left) and the reverse 0.6 → 1.4 (bottom panel, left to right) at = 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 = 1.15, which jumps to a 6-species equilibrium by the invasion of three more species at = 1.11; another two species displace one of the previous at = 0.9; and at = 0.72 a sudden jump onto a chaotic attractor occurs.Reversing the adiabatic protocol, the transition from chaos to xed point occurs only at = 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.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 ve 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 dened by Main Text Eq. ( 13).C The ratio ∕ tends towards 1 at the boundary to the equilibrium phase.Note that 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.In one dimension, * ( ) must be constant.Since is a non-negative abundance in our case, we must impose a boundary at = 0 through which probability cannot ow.Therefore * ≡ 0. The solution for * is then

Figure 1 :
Figure1: 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 , 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.

Figure 2 :
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.Di erent lines represent observations at well-separated time points.B Species abundance distribution (SAD, blue histogram) corresponding to the blue rank-abundance plot; overlaid, abundance uctuation distribution (AFD), averaged over all species (black line)with ± 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 uctuations.The one bar missing from the SAD is the e ect of nite species richness, as high-abundance bins only ever contain a couple of species for = 500.The vertical dashed line indicates the immigration level which determines a lower limit to abundances.

Figure 3 :
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 uctuation distribution (AFD) from Figure2(black), and the AFD of the focal-species model (pink): excellent agreement is found for the powerlaw section.The 'uni ed 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.

Figure 4 :
Figure4: Species di erences in dominance.A Example of a long abundance time series for the three species who are ranked rst, 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 results in a proportional decrease in median dominance time.C Distribution of dominance biases against relative dominance rank for a range of : there appears to be convergence towards a non-constant limiting distribution, implying that net species di erences are not due to small-e ects.Note that, by de nition, 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 normalized sum of interaction coe cients, Eq. (4): lower net competition correlates with higher dominance bias.Species in the tails of the distribution are also less 'typical', with typicality quanti ed by the index , Eq. (16), representing the similarity of a species AFD to the species-averaged AFD.Panel A and D are both for = 500.

Figure 5 :
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 and (with 0.01 increment), 30 simulations were made, each with a di erent random initial condition ∼ ( ,2∕) 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 multipleattractor phase (in which cycles and multi-stable xed point are common in addition to chaos) is not sharp, unless probed adiabatically in the way explained in Supplementary FigureS8.The unique xed-point phase has been studied analytically in the weak-interaction regime ( ∼ 1∕ ).When inter-speci c competition is in general stronger than intra-speci c competition, a single species (identity depending on initial condition) dominates, in line with the classical competitive exclusion principle[60].

Figure 6 :
Figure 6: Relations between e ective parameters in the chaotic phase.A Colour legend of the chaotic phase (boundaries from Figure5).Each pair of ( , ) has been mapped to a distinct colour.B The exponent of the power-law section of the AFD for the chaotic dLV model plotted against the analogue 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.C Co-dependence of the e ective parameters , , : the amplitude of growth-rate uctuations approximately equals the absolute value of the negative growth rate (only weakly depending on and ; Supplementary FigureS10); is roughly proportional to the inverse turnover time, but the slope of the relationship depends on and .

Figure 7 :
Figure 7: Statistical properties of the e ective noise.A, B Time series and distribution of rel = ∕ − 1, etc. C, D Histograms of ( ) across allspecies and time (grey), over just species for one random time (green), over all time for the rst/mid/last-ranked species with respect to average abundance (blue/pink/yellow), with (0, 1) (black, dashed) for reference.E The empirical distribution of in Eq. (23) over all species and times, compared to the distribution ( , ) assumed for in the focal-species model.F Autocorrelation functions: for every species (grey), rst/mid/last-rank species (blue/pink/yellow)), and the average over all species (black).The left inset compares the ACFs of (green), (black), and the exponential t to the latter (red); the right inset shows the distribution of the parameter in exponential ts to each species ACF.

Figure S1 :
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 a ect 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 coe cients by (10 −6 ); D a perturbation of the initial abundances by (10 −8 ).E Each type of perturbation leads to completely di erent community composition compared to the reference (measured as Bray-Curtis similarity) after a few hundred time units.

Figure S2 :
Figure S2: Convergence to positive maximum Lyapunov exponent (MLE).AThe dominant nite-time Lyapunov exponent (FTLE) over a few integration time steps ( = 2) uctuates 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 S3 :
Figure S3: Decay of community similarity with time..A The temporal similarity matrix has elements given by the Bray-Curtis similarity between the abundance vectors at two time points, ( , ) = BC( ( ), ( )).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 ≈ ≈ 8000 re ects a time when some dominant component persisted for an unusually long time.B For a few well-separated time points (one graph each), we show how

Figure S4 :
Figure S4: Scaling of e ective community size with richness.The time-average of the e ective community size, e (Main Text Eq. (6)), increases slowly (but super-logarithmically) with the overall richness .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 S5 :
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 ∼ (0, 1).B Non-uniform carrying capacities:∼ LogNorm(0, 0.1).C Sparse interactions: each interaction has a 0.1 chance to be non-zero.D Symmetric bias: = 0.2 correlation between diagonally opposed interaction coe cients.E Predator-prey bias: = −0.3correlation between diagonally opposed interaction coe cients.

Figure S6 :
Figure S6: Pairwise correlations in species abundances.While most of the ( − 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 coe cient with lag time lag ( lag ), and the horizontal axis has the rescaled interaction coe cient = ( − )∕ .The inset show that most zero-lag correlation coe cients 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 (0) for every species ; in order to see if correlations are stringer if we optimize over the delay time, we show in light blue max , lag <200 ( lag ).It is seen that the maximal correlations are around 0.25 in size, and clearly associated with < 0, i.e. a less-than-averagely negative (even positive) e ect of species on .Similarly, the extremal anti-correlations (pink for zero time lag, and light pink optimizing over time lag) are associated with > 0, i.e. a particularly negative e ect of on .

Figure S7 :
Figure S7: Scaling of AFD power-law exponent with and .From simulations, we have extracted the slope of the power-law section of the abundance uctuation distribution (AFD) A For varying , we nd empirically that the exponent depends linearly on the logarithm of species richness, with coe cients that depend on the system's other parameters: = 0 + log , where 0 = 0 ( , , ), = ( , , ).B For varying , the exponent appears to follow = 1 − ∕ log , where = ( , , ).The values of 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 = 1 with the value at = 10 −8 .

Figure S8 :
Figure S8: Phase diagram form adiabatic simulations.Adiabatic simulations allow to track, in a numerically e cient fashion, the attractors of the dynamics as model parameters are changed slowly and continuously.To make the interaction statistics and continuous parameters of the model, we use as interaction matrix ( , ) = +where is a single, xed realization of a standard Gaussian random matrix.A For each value of , we initialized separate simulation runs starting at = 1.4,and let their abundances evolve until an attractor was found.For each run, we then changed by small increments = −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 di erent attractor.If a state diverged, the initial abundances for the next value of were set as the most recent non-divergent attractor.Thus, each simulation traced a sequence of attractors from = 1.4 → −0.1, corresponding to a horizontal line in the phase diagram.The colour quality re ects the class of the attractor, and the colour gradation indicates the e ective community size, revealing the following features: First, we nd mostly xed points in the multiple attractor region.This is because, once a xed point is converged to, it is "hold on to" until it vanishes or changes stability.If, instead, every simulation at given , would start from newly sampled initial abundances and interaction matrix, we would nd di erent attractors every time, and the diagram becomes more heterogeneous (compare Main Text Figure5).Second, clear lines radiate from

Figure S9 :
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 uctuations.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 S10 :
Figure S10: Dependence of the e ective parameters , on , .The empirical, approximate relationship ∝ , found across the range of , in the chaotic phase, has a proportionality constant that depends relatively weakly on and .

Figure S11 :
Figure S11: Collective correlation.Within the bounds to the chaotic phase indicated in FigureS8A, 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 ve times, else the point was omitted (chaos probability was shown in Main Text Figure5).A The collective correlation.B The critical value of the collective correlation as dened by Main Text Eq. (13).C The ratio ∕ tends towards 1 at the boundary to the equilibrium phase.Note that changes continuously across this boundary (FigureS9C).The arrow at the upper end of the colour bar implies the range has been capped for clarity.

Figure S12 :
Figure S12: Power-law exponent in the chaotic phase.A Variation of the AFD power-law exponent across the chaotic phase.Apart from outliers, we nd an exponent larger than one.B To test the accuracy of the focalspecies model in predicting the exponent, we measure the relative error in − 1 (since we expect > 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 S13 :
Figure S13: Comparison of autocorrelation times.We compare the autocorrelation time of the abundance vector and the autocorrelation time of the e ective noise .These two parameters are obtained by the exponential t − ∕ applied to the respective autocorrelation functions.Across the chaotic phase, these to timescale are quantitatively close, for reasons explained in Main TextAppendix E.