Speciation in a MacArthur model predicts growth, stability, and adaptation in ecosystem dynamics

Ecosystem dynamics is often considered driven by a coupling of species’ resource consumption and its population size dynamics. Such resource-population dynamics is captured by MacArthur-type models. One biologically relevant feature that would also need to be captured by such models is the introduction of new and different species. Speciation introduces a stochastic component in the otherwise deterministic MacArthur theory. We describe here how speciation can be implemented to yield a model that is consistent with current theory on equilibrium resource-consumer models, but also displays readily observable rank diversity metric changes. The model also reproduces a priority effect. Adding speciation to a MacArthur-style model provides an attractively simple extension to explore the rich dynamics in evolving ecosystems.


Introduction
An ecosystem is a set of species, each of finite population size that interact by competing for finite resources that fuel their growth. A single ecosystem can involve dynamics that occurs over a wide range of length and timescales Azaele et al (2016); Darwin already eloquently referred to this in his "tangled bank" remark. The seemingly universal nature of a species' emergence, adaptation and extinction in such ecosystems, has inspired many to describe the phenomenology of ecosystem dynamics with simple modeling with only a few ingredients that are independent of the specific physical mechanisms at play Nowak (2006); Azaele et al (2016). What is then the simplest quantitative description that displays all the salient dynamical features of evolution? This question has a long list of partial answers Nowak (2006); Tikhonov (2016); Posfai et al (2017), although much work in the field concerns equilibria MacArthur and Wilson (1963); MacArthur (1955); Chesson (1990); Hubbell (2001). Here we show that the already successful variants of MacArthur models can be amended with a simple stochastic mechanism that introduces new species, which allows us to show many of the biologically relevant dynamical features of evolution even when the starting point for the dynamics is a single primordial ancestor. In particular, our MacArthur model variant can describe the growth dynamics of an ecosystem in terms of species richness, its evolution towards a dynamic equilibrium size, and even its adaptation to resource influx changes. The predictions of the model are consistent with other existing modeling on for example equilibrium dynamics, and reasonably in line with common observations on for example resource shock experiments.
2 The speciating MacArthur model Evolutionary dynamics modeling including MacArthur type models usually start with describing population dynamics asṅ(t) = n(t)f (n) with n the set of species population sizes and the dot denotes a time derivative. The crux is that f (n) is not constant but a growth rate determining function that depends on ecosystem features such as population sizes and coupling constants which specify inter-species competition and preying efficiency Wangersky (1978); Cressman and Tao (2014). This general approach is tremendously successful even in capturing quantitative experimental observations of low dimension systems Korolev et al (2011). However, describing the dynamics of larger ecosystems with many evolving species is challenging, as high dimensional systems quickly lose their numerical and analytical tractability, even without incorporating the additional complexity of the evolution of each species.
We address this complexity by adding evolutionary dynamics to multispecies ecosystem models by quantifying the growth rate within the context of the ecosystem properties. It has become customary in recent years to define every species j by a strategy vector s j Posfai et al (2017); Tikhonov and Monasson (2017); Pacciani-Mori et al (2020); Caetano et al (2021) that couples the growth dynamicsṅ(t) = n(t)f (n) to a dynamic resource vector r(t) which represents the amount of available resources at time t. Here, each component i of s j describes which fraction of each resource r i is used by every species j at every time step. The time-dependent growth factor is then naturally captured by the alignment s j · r. We can write for each elementṅ j = n j f j (n, s j · r) while introducing resource time-dependence via a functionṙ = g(n, s), where s is the matrix with strategies s j as its columns. The concept of a strategy vector s j that every species j ∈ {1 . . . k} has in order to harvest resources is central in the model. Each component of s j represents a strategy for a particular resource, the ensemble of which is characterized by a time-dependent vector r = {r 1 , r 2 , . . . , r l } where l is the number of resources. s j quantifies how much of each resource every individual would like to take out of the resource bath. We consider the concept of a resource component r i as extremely general: it can refer to a specific molecule, a chemical energy influx, or even to a certain amount of space available in a habitat. Each species has limited energy and time to spend harvesting. Thus, species need to optimize their foraging behavior by choosing how to be efficient as regards different resource consumption MacArthur and Pianka (1966). To express the interdependence of resources, one simple way is to fix a norm of the strategy vector to an arbitrary quantity. For convenience we choose s j 2 = 1, but other norms and values can be chosen. This choice does have a biological significance Caetano et al (2021). For now, we will focus on speciation with this fixed choice, but in Sec. 7.3 we will come back to it.
We assume that the total demand for resources is proportional to s j and the population size of each species n j . The time dynamics of every resource component r i is then described by: Here β is a timescale, and γ i is the resource replenishment factor of resource i, essentially representing a chemostat Posfai et al (2017). In this simplified linearized resource dynamics, r i is not strictly positive, which is unphysical. When r i < 0 we set it to zero. We find that the time-dependence of the model is very sensitive to the choice of the rate of consumption (and growth). However, the approach of a dynamic equilibrium while adding species to the ecosystem is preserved regardless of the choice of growth rate factor. If the preferred resource intake of the species is similar to the composition of the resource environment r, the growth rate should be maximal; in the case where s j and r are not aligned, the species should perform poorly. A species goes extinct when its population size is below a threshold value. We verified that threshold choice is not important for much of the dynamics observed. We do note that setting the threshold lower trivially increases the total ecosystem size. It is now natural to write for n j that Writing explicitly that s j · r ≡ i=1...l s ij r i makes clear that s ij is the resource utilization coefficient of species j for resource i. α is, again, a time constant; δ sets the population decay rate. Eqns. 1 and 2 are a simplified version of MacArthur equations MacArthur (1970); Chesson (1990);Haygood (2002). However, we interpret the coupling matrix s ij much more specifically: it is essential to see how s j here serves as the definition of species j Posfai et al (2017).

Stochastic speciation
So far, we did not consider any stochasticity: the skeleton of the dynamics is deterministic and embeds a selection for the fittest species Vellend (2010). The novel feature in this work is that we add speciation dynamics in two ways: (1) by adding species and (2) by ensuring the added species can be different from existing ones. We choose to make (i) mutants appear randomly, adding a new species equation to the system in a Monte Carlo way. We shall see that this captures both evolution and invasion. (ii) The strategy vector of the newborn species is stochastically generated Posfai et al (2017); Pacciani-Mori et al (2020); Drake (1990); Serván et al (2018); May (1972). Thus, the community assembly happens sequentially at random times with randomly evolved species-see Fig. 1. We call our speciating MacArthur approach towards ecosystem dynamics "SMA" for brevity.
We can now use s j to specify how a species evolves with a simple Monte Carlo evolutionary model: mutants are generated from existing species from which they differ only in terms of the harvesting strategy vector. The initial species in the ecosystem is defined by uniformly drawing an ancestral strategy vector s 1 , with elements in the range (0, 1), and normalizing it to one in the Euclidean norm. A new species can be spawned for every time step and for every alive species, when a random number, drawn from a standard normal distribution, is larger than ν standard deviations. We can thus define a mutant k + 1 by taking any existing species strategy s j and by adding a noise vector: The noise vector is composed of a random vector drawn from a normal distribution, ψ, weighted by a parameter η quantifying the amplitude of the mutation. Therefore, the noise vector represents a shift in the species' resource utilization composition as a consequence of mutations. Note that by Fig. 1 The phylogenetic tree derived from the emergence of species in the evolution of one realization. Time is expressed in terms of integration time steps. The species lifetime is indicated in color to provide the chronology of emerged species in the tree. Leaf 1 indicates the emergence of the last spawned species, and so on. The red line shows the lineage of the species with which the tree started; it went extinct around t ∼ 4.4 × 10 4 steps.
so "phenotypically" defining our species solely in terms of s j , a natural link to genetic variation within a species is lost. The stochastic arrival of a new species and the resource richness influence the local selection outcome by possibly inducing historical contingency and priority effect in the community assembly Fukami (2015); Almany (2003); Sale (1977). The level of historical dependency on species' arrival hinges on the value of η and the number of different resource types. We assume that a small η represents an infinitesimal evolutionary mutation in a species' survival strategies, originating speciation events in response to the dynamic resource landscape. Indeed, η defines the degree of strategies' divergence from parent species to daughter species. On the contrary, foreign invasions are modeled by considering the arrival of a new species with an entirely new set of characteristics, uncorrelated with the ones already present in the system. We call such case η → ∞ and will be discussed more in detail in Sec. 5.2. As we will see, SMA does show a remarkable ability to reproduce behaviors that can be interpreted in any eco-evolutionary context. We will discuss the interpretation of η in more detail in Sec. 6.

Biological example: biofilms
Despite the sober mathematical formulation, SMA conceptually captures some essential features of ecological systems: in several natural ecosystems, ecological successions are intertwined with populations' adaptation to environmental conditions. Moreover, both native species' evolution and alien species' invasion contribute to determining the community's fate, together with environmental responses or sudden changes. One simple yet effective example is that of microbial communities in biofilm formation, where evolutionary and ecological timescales are comparable Hansen et al (2007); Goyal et al (2022). In subaerial biofilms, such as those that grow in monumental buildings, the substrate of stones is firstly colonized by pioneer airborne microbes, which then leads to further successive stages with the subsequent invasions of other microorganisms Gorbushina (2007). The different substrate characteristics and environmental conditions define the dynamic resource richness, which in SMA is rendered via the types of resources and their influx vector C. Gaylarde (2020); Ariño et al (2010); Caneva et al (2004). Moreover, microorganisms not only can feed on others' metabolic discards, but are also able to evolve rapidly via strains' mutations, which in our model is captured by the signs and mutations in the strategy vectors Gorbushina (2007). As a result, several biofilms quickly show resistance to chemical anti-degradation treatments, often leading to unexpected and new community structures Simões et al (2009);Gorbushina (2007). It is evident that biofilms are much more complex ecosystems than the ones described by SMA. In biofilm formation, both evolutionary speciation and ecological invasion act simultaneously, while in the current version of SMA we consider these processes separately for simplicity. SMA can of course include both effects simultaneously, yet we aim to disentangle the dynamics observed in a simple general framework that captures existing natural systems. Many other interpretations besides biofilms are possible and welcome.

Implementation
We run all ecosystems starting from one species, with α = 0.005, β = 0.01, δ = 0.1. We assume that all resources have an equal influx rate given by γ i = 1. We focus on the case of ν = 3.8 and we consider systems where we vary l. The value for l can be representative of several different community scales. For example, when considering microbial communities, 100 or more different resource types are an appropriate choice Fischbach and Clardy (2007); Fischbach and Sonnenburg (2011);Tikhonov and Monasson (2017). We then varied η within previously defined limits and studied the ecosystem's evolution.

Solver
To integrate the species dynamics in Eq. 2, we employ a fourth-order Runge-Kutta (RK4) method with stochastic elements that can only generate one new species per existing one at every time step. We checked that the RK4 accuracy used in all our calculations does not affect the results. We verified that our custom implementation provides similar performance to the standard MATLAB ode45 solver for the deterministic l = 1, k = 1 case, with α = 0.05, β = 0.01, γ i = 1, δ = 0.2. In the stochastic setting, we loop in every time step over the extant species and compute for each species the RK4 step and add a new species when a number drawn from a standard normal distribution is, in absolute value, larger than ν. After the loop over the species, we perform an Euler forward step for the resources dynamics Eq. 1 where we use the updated values for the species. Then we set all negative values for r i to zero to ensure positivity of the resources. Representative code is provided on Zenodo zen (2022).

Example ecosystem
The SMA model can show a wide range of different dynamics, depending on the (initial) number of species k and total number of resources l and other parameters. The phenomenology of speciation embedded in the SMA can however already be observed for starting evolution with the most stringent starting condition of one species, that is, k = 1. Note that capturing the emergence of an ecosystem with interacting species from a single primordial reproducing entity is one explicit aim of the current modeling approach. Evolving an ecosystem from a single species is, for example, not possible in the classical MacArthur, Lotka-Volterra (LV) or replicator equation contexts.
In this example we choose an l = 5 resource space and set r i (0) = 10 for all i as initial resource amount available. Due to the presence of multiple species, the vector r(t) will be time-dependent and may not always be aligned with a particular species vector s j . We use α = 0.005, β = 0.01, γ i = 1, δ = 0.1, ν = 3.7 and η = 1 and evolve the system for up to 10 5 time steps of size h = 0.1. Note that the value of ν is intrinsically linked to the choice of h because h also sets the frequency at which new species are generated; we come back to this point in Sec. 4.2. We consider a species extinct if its size is smaller than 0.1n start , where n start is the initial size of the population; this threshold effectively captures the role of fluctuations in small populations Reichenbach et al (2006); Parker and Kamenev (2009);Huang et al (2015). The specific value for the extinction threshold does not affect the essence of the evolutionary dynamics of r and n. For repeated independent ecosystem calculations, we generate a new s 1 for every iteration. Fig. 2 shows the typical ecosystem evolution initiated from one species k = 1 starting at size n 1 = 10. Due to the stochastic nature of species emergence and extinction, every realization of ecosystem dynamics is different. However, several important qualitative features reproduce and are visible in this and any example: (i ) The initial species size oscillates in time until a viable new species has emerged; in sync, the resource dynamics is also oscillatory for the smallest component of s as all the other resources get quenched to zero Huang et al (2017). Note that this excludes the case when an s i is strictly zero, which is possible but rare. This is such that the resulting behavior dynamically balances the resource usage with the resource influx. (ii ) The emergence of new species affects the timescale of periodic oscillations; also new species can make older species go extinct. (iii ) Later in the evolution, the population fluctuations shift in frequency and decay in amplitude and multiple resources become utilized. The interpretation of these three trends is clear: the randomly selected initial species favors the survival of one resource, for which s i is the smallest. After this transient, the dynamics follow the l = 1, k = 1 system which is pseudo-LV in character and allows for periodic orbits of fixed frequency. In this phase, the possibility of the random emergence of new species is consequential: the emergence of new species that are ηψ different from their parents will suppress the dominant role of the first species and limit its overuse of other resources, thus making the remaining resources emerge again as they are always continuously replenished at rate γ i . We will make these statements more quantitative in the next sections.

Main Phenomenology
The first significant result from SMA is the naturally bounded ecosystem it produces in both size and structure while we neither fix the (maximum) number of hosted species in the community nor the maximum population size of the individual species; only the influx of resources γ is bounded. In modeling, this is traditionally captured with logistic growth models and/or to restrict oneself to probing the dynamics of an ecosystem with a fixed number of species MacArthur (1970); Chesson (1990); Haygood (2002) (2017); Grilli et al (2017). SMA embeds size limitations naturally, as we observe that for enough simulation time, the number of species grows towards a long-term stationary value-see below. This maximum number of coexisting species is solely determined by the distribution of strategies and resource availability. The bounded growth feature allows us to explore the long-term species abundance distribution (SAD) Hubbell (2001);McGill et al (2007), with the knowledge that a system will maintain, on average, a constant number of competitors and, as we will see, a finite global size. Note that during equilibrium size, the model allows for and will randomly let species emerge and go extinct; the ceiling represents a dynamic equilibrium. Note that in much of the dynamics explored, ν mostly sets the speciation rate for the system evolution and is

Growth towards equilibrium dynamics
Rather than considering each species or resource type separately, we gain insight into the system evolution as a whole by considering the total number of individuals (N ), and the total amount of resources (R). Interest in the total abundances dynamics for similar trophic species is seldom suggested and, to our knowledge, rarely explored Posfai et al (2017). Yet, empirical evidence shows that the aggregate biomass could provide valid information on the stability and composition of a community Tilman et al (1997); Doak et al (1998). Thus, studying the evolution of total abundances allows us to explore the system behavior in greater depth from a new perspective and, simultaneously, reduces the variables involved. Resultant equilibrium dynamics for ecosystem averages for different l are shown in Fig. 3a. The inset shows how the equilibrium is reached for a particular example setting of l = 2. We find a family of fixed points for the (R, N ) dynamics that are all on a line defined by This total abundance dynamics can be understood by simply considering separately the sum of the resources R and the sum of the species N . In this way, the dynamically evolving dimension of the system reduces to a two-dimensional problem Here, we emphasize that the number of living species, k, is a function of time: the equilibrium is dynamic in nature. The stationary solution (R * , N * ) of this system should solve the equations but note again that the elements that make up R * and N * do not have to be stationary. Inspired by Fig. 2, we now assume r i is approximately constant in i, meaning that the mean abundance per resource does not vary too much per resource, we can pull r i = R * l out of the sum, resulting indeed in the fraction It turns out that this equation predicts the slope of the line in the (R, N ) phase space on which all the attractors of the total abundances dynamics lie, as is shown in Fig. 3a. Note however that the pictures shown are for a constant γ, and numerical results seem to indicate that the assumption r i ≈ R * l becomes less valid when γ i = γ j . Using the results from Posfai et al (2017), we can understand why this line has such predictive power. In the deterministic version of our model, i.e. without speciation, any number of species can coexist. That is, as long as the geometric conditions introduced in Posfai et al (2017) on the strategy vectors s j and the replenishment γ are met. When these conditions are met, the system converges to a fixed point where all r i attain the same value. Hence, what we observe is that every time a new species is introduced (or an old one removed), the dynamics converges to a new fixed point that is indistinguishable from the old one in the R-N dynamics. We conclude that the stochasticity in our system always results in an ecosystem where the necessary geometric conditions for the coexistence of many species are met. To be precise, in Posfai et al (2017) results were obtained for a nonlinear version of our model, such as Eq. 8, with normalization in L 1 instead of the Euclidean norm. This norm changes some of the details, see Sec. 7.3. Also, the nonlinear version of the model changes the slope of Eq. 7 somewhat, while making it valid under more general types of γ.

Transient scaling with ν
A second feature in SMA is that the speciation threshold, ν, induces a timescale τ for the evolution. At every time step, each species has a probability p m to mutate. This probability is given by the tail (≥ ν) of the standard normal distribution, i.e. p m = erfc(ν/ √ 2). Therefore, in absence of extinction, we expect the average number of species S a to grow as ∼ exp(erfc(ν/ √ 2)t). However, it should be noted that ν is intrinsically entangled with the choice of the time step h: indeed, ν defines the probability of an alive species generating a new species in the h time unit. In our simulations, we always kept h constant at the value 0.1 which ensures the stability of the solver. Because of this choice, we must take S a ∼ exp(erfc(ν/ √ 2)t/h). For ν to set such a timescale, the average behavior of a dynamic observable obtained for different values of ν, when plotted as a function of τ = exp(erfc(ν/ √ 2)t/h), should collapse into the same master curve. We explored this possibility by considering as observable the time-dependent average of living species S a over 150 realizations, keeping η and l fixed but setting ν = 3.8, 3.9 and 4. We rescaled the time axis of each curve associated with a different value of ν to τ . For small l and high values of η, we observe that ν does induce a timescale for the overall evolution: indeed, the S a curves associated with the different ν tend to collapse towards a single curve-see Fig. 4a,b in which we display results for l = 2 and l = 100 respectively and η → ∞. It follows that for low ν the evolution is faster, while high values of ν slow the community's formation.
It is also clear from the figures that τ is a valid timescale for small communities. When the number of species grows, extinction becomes important and the curves for S a deviate weakly from the timescale τ . Quantitatively, when the average number of extinctions S e > 1 we find that the rescaling becomes less accurate. Curiously, the rescaling works very well for all S e .

Adaptation
The resource alignment interpretation of SMA clearly gives it many physically meaningful links to real world ecosystem dynamics. Empirically, a promising constraint is to provide a time-varying resource influx by introducing γ i (t), sometimes also called a "pulse" experiment Tilman (1987); Hiltunen et al (2015). We demonstrate in what follows that SMA shows adaptation under such conditions. Additionally, we explore how the magnitude of changes in subsequent generations as characterized by η affects adaptation dynamics.

Rank abundance
To demonstrate the effects of a "pulse", we focus on the results obtained for a system characterized solely by invasion events, and l = 100. We chose to make the resource shock occur at t c = 8 × 10 4 steps, and we doubled the length of the simulation to provide enough time for the system to respond to the perturbation. For t < t c , the influx rates are γ i = 1 for all i; when t ≥ t c , the new i γ i is three times that before the perturbation. We chose to distribute 75% of the new γ among only 25% of the resources. This abrupt change in resource influx induces adaptation dynamics by the ecosystem. Solving SMA with time-dependent resource influx over several realizations at previously defined α, β, δ, ν, we observe that ecosystems are able to recover from such a resource shock: when the perturbation occurs, there is an initial stage, after which R(t) and N (t) gradually restore their limit-cycles (not shown). However, for the linear model used here, the position of the attractor in the (R, N ) phase space changes according to the new resource influx rates. The center of the oscillations do not lie on the line with slope l i γi δβ α anymore. On the contrary, given the unevenness of the new resource influx vector, the correct slope seems now proportional to the average influxes of the resource types that are not fully depleted, which in general are the ones associated with the highest resource influx. Fig. 5 (a) The number of living species, Sa, for a single realization of a system with l = 100, ν = 3.8 and η → ∞. i) From t ≈ 4.5 × 10 4 steps on, the curve reaches a plateau; ii) The resource shock perturbs the system at tc = 8 × 10 4 steps; iii) After the shock, a second plateau is reached. (b) Rank-abundance plot for the realization in (a). The curves show the trend at every 1000 time steps from t = 8 × 10 4 to t = 9.1 × 10 4 . From t = 9.1 × 10 4 to the end of the simulation, corresponding to when the Sa reach the second plateau in (a), the Rank-abundance curves are displayed every 14000 time steps. The color scheme follows the colors on the left. The solid red line is the curve when the shock occurs. The arrows indicate the time immediately after the perturbation.
We quantify the pulse response by probing species occurrence. Interestingly, the number of coexisting species is strongly affected by the resource shock. In the absence of perturbations, the number of living species hosted in the system spontaneously grows until it reaches an average maximum value in time-see Fig. 5a. When the perturbation occurs, the increase in available resources initially encourages the system to welcome new species, resulting in a sharp peak in the number of coexisting species. Subsequently, the living species curve decays with a characteristic timescale until it reaches a substantially lower new stationary value-see Fig. 5a. The new rank abundance distribution corresponds to having fewer species that are all large in population size.
Remarkably, the shock also influences the SAD-see Fig. 5b. Before the resource shock occurs, the rank-abundance plot, also known as Whittaker plot Magurran (2013), displays a curve that gradually collapses towards lognormal-like behavior in the tail. Such behavior resembles that observed in empirical data Sugihara (1980) May (1975), although a few methodological aspects that give rise to such distribution are still debated Magurran (2013); May (1975). After the perturbation, the curve still preserves its characteristic shape. However, its slope gets steeper in time, and the curve reaches a new asymptotic behavior characterized by less species evenness. Curiously, for systems with l = 2, the rank-abundance trend shows a strongly uneven species distribution that is often associated with harsh environments or early stages of successions Magurran (2013)

The role of η
The noise amplitude η has two biologically different limiting cases. For η → 0, the community's evolution proceeds via infinitesimal steps, with all the new species occupying the same niche. On the contrary, we can imagine a scenario in which a foreign species invades the community from an external pool. In this case, we assume that the foreign species evolved from a different ancestral species. Thus, we define its strategy by drawing a new ancestral one and adding a noise vector with the maximum noise amplitude η = 1. This approach preserves the biological interpretation of the ancestral strategy and the mutations, originating from two different distributions: uniform and standard normal. For simplicity, we will refer to this scenario with the term η → ∞, although mathematically we do not explore the limit of η → ∞.

The case of η → 0
As discussed, by decreasing η, the realizations start to depend on their initial conditions. When η is infinitesimal, ecosystem dynamics are mainly determined by the ancestor features. However, such development of a neutral community at a species level MacArthur and Wilson (1967) is made dynamic by a priority effect. Simply put, a small η is likely to lead to a successful species only if its ancestor was also successful. For small values of l, each realization still defines limit cycles around a fixed point. By increasing l, the dynamics becomes aperiodic. When l is small, as expected the family of fixed points lie along the line defined by the ratio l γi δβ α . Also for l = 100, the dynamics, even though aperiodic, is still contained in a region of the phase space close to the line-see Fig. 6. We conclude that tuning η allows us to apply SMA to both evolutionary and invasion-type dynamics. The so embedded co-occurrence of both selection and priority effect mirrors empirical evidence and theoretical hypotheses suggesting that stochasticity and determinism in community assembly work hand in hand Chase and Myers (2011)

The case of η → ∞
For systems invaded by foreign species, the dynamics of several different realizations of one system exhibit the same attractor and qualitative behavior-see Fig. 3. From said figure it is clear that at long time scales, the dynamics settles on quasiperiodic orbits around points on the line of fixed points defined by Eq. 4, and higher values of l result in higher values of R and N . Varying the parameters α, β, γ, and δ gives similar results, only changing the slope of the line. Moreover, for high η, in the range [0.5, 1], the results are similar to those obtained for η → ∞. In this limit of η, especially the large l limit is interesting, because at small l, the ecosystem quickly selects the best adjusted strategies, all the others going extinct. For large l, species' strategies are constantly evolving towards an existing optimum that is however statistically unlikely to achieve, leading to slow dynamics. The effect of introducing new species is now also determined by their time of arrival, which now defines their competitiveness, inducing a priority effect Fukami (2015); Almany (2003); Sale (1977) that we will see is the dominant driver of dynamics in the case η → ∞.

Priority effect for η → ∞
In the simulations we see that species that spawn at the end of the simulation have a lower chance of surviving than at the start of the simulation. This can be interpreted as a priority effect. In order to quantify this, we count in a given time interval the number of species that went extinct immediately, i.e. decay exponentially from the initial population to the extinction threshold. Dividing this number by the total number of spawned species in the same interval gives us an immediate extinction probability P ext . When the number of spawned species is large enough, we can divide the list of species up into bins and calculate P ext for each bin separately which indicates how P ext changes over time. Fig. 7 shows a clear trend for an example ecosystem which indicates that there indeed is a priority effect.
6 Interpreting the role of η 6.1 Noise, priority effects and historical contingency We observed that the SMA dynamics solely depends on the number of different resources if the community assembly emerges from adding essentially random species. This scenario comes about in the limit of large η, thus with more significant differentiation between parent and child species, resulting eventually in always the same type of community structure and species' distribution.The late-time community that emerges in this limit is commonly referred to as the climax community and indicates the final ecological succession of the community formation Morin (2009);Weiher and Keddy (1995). The species, specifically the strategies selected to survive in the climax community, are considered resistant to the invasion of new strategies' variants, which might be considered a type of "priority effect". On the other limit, we observe that the community assembly is affected by the history of the species' arrival for low values of η. The community is thus historically contingent Morin (2009); Belyea and Lancaster (1999) ;Schröder et al (2005), an observation that is coherent with the literature. The smaller values of η result in more minor divergences between parent and child species. The community is locally neutral because the species belong to similar trophic levels and have comparable survival probabilities. Thus, the community assembly is likely susceptible to the species' arrival history and the pioneer species' features Fukami et al (2007).

Resource-noise interactions
What is the influence of evolutionary noise when there is not much room to be different in resource space? What is the role of noise when there are many resources? Clearly, l and η span a phase space of dynamics that we briefly explore for the biologically relevant dynamics it can describe.
For most realizations of systems with η → 0 and l = 2, the number of living species S a per realization displays a sigmoidal behavior in time-see Fig. 8a. This results from the resource-consumer feedback: the total consumers increase by consuming resources before reducing again due to excessive competition and resource lack. However, in this scenario, the species are almost equivalent in their competitiveness, and they can only go extinct if they arrive in the community at an unfavorable time. It seems reasonable to assume that when Fig. 8 (a) The number of living species, Sa, for 5 realizations uniformly drawn from a 150 pool of a system with l = 2, ν = 3.8 and η → 0. (b) Same as (a) for a system with l = 100, ν = 3.8 and η → ∞. The simulation length of the system in (b) is 8 × 10 4 steps due to memory load issues. For both the system in (a) and (b) the Sa curves reaches a long time plateau whose value are different for each realization. l = 2, only a limited pool of strategies can survive in what we can call a "harsh" environment. Consequently, the majority of invading species do not present the necessary conditions and quickly go extinct. Successful invasions become rare events, and growth of S a stalls, reaching an equilibrium.
Similarly, also systems with η → ∞ and l = 100 show a long-term stationary value for S a -see Fig. 8b. Here the strong species selection provides additional resource-consumer feedback by choosing the most convenient strategies, leading the less fit species to extinction. Regardless of the interpretation, it is notable that again we observe that long term stability is achieved under very different settings. Note also the much larger number of species that manages to coexist in the limit η → ∞ and l = 100 than when l is small. On the contrary, for the case with η → 0, l = 100 the system does not present any long-term stationary value for the number of living species within the available computation time (data not shown). In this limit, the species take a long time to adjust their size according to the availability of the many resources types. In the limit of η → ∞ and l = 2, many realizations of the system hint at the existence of a long time plateau value of S a . Also here, it is reasonable to assume that when l = 2, a limited pool of strategies can survive in such a harsh environment, but apparently invasion in harsh environments is notably different in its dynamics than evolutionary speciation.

Generalizations of SMA
SMA allows for many further generalizations MacArthur (1970); Chesson (1990);Haygood (2002). In Eq. 1 we have only considered that species consume resources; however, they may also provide resources-see Sec. 7.1. Oxygenic photosynthesis Knoll and Nowak (2017) is one example; on a different scale also gut microbes provide natural resources for each other Faust and Raes (2012); Vet et al (2018). Besides, predator dynamics can be introduced by adding another predator coupling matrix term i m ij n i , which can have both positive and negative elements, when species j is a predator or prey respectively. Growth rates can be made an explicit function of r i , preserving much of the dynamics presented here but adding more biologically relevant constraints. Several other choices can be modified, such as making the speciation rate or η a function of n. Making r i < 0 for some i can account for stressors. We discuss some of these aspects in this section.

Negative strategies
Enabling the components of the strategies s j to have also negative values is a natural choice or ecological dynamical modeling, as it allows for species to contribute to the resources of other species. We find that allowing for the sign change of s j results in a significantly different transient. For example, when one samples s j from the full normal distribution, we observe that for infinitesimal mutations, η → 0, and small numbers of resources, l, no negative strategy appears or survives. For high l and η values, on the contrary, a portion of negative strategies survives. As a result, the total abundance dynamics presents aperiodic behavior, as the negative strategies work as an additional resource influx rate, with an intrinsic stochastic nature given by the random arrival of species with such features-see Fig. 9.

Consumption and growth rates as Monod functions
Until now, we assumed that the resource consumption rate β was constant, but it can be reasonable to assume that β is a function of the resource availability, so β = β(r i ). The consumption rate then depends on the species opportunity to find and consume resources. It turns out that also for SMA, this consumption rate function is an important factor that determines certain characteristics of the dynamics. To explore the role of the consumption rate function, one relevant choice for β(r i ) is the Monod function Monod (1949): β max ri K+ri , where K (K > 0) defines the half-saturation constant, that is, the resource availability that is present in the system when the consumption rate reaches half-speed, β = β max /2. The Monod function usually describes bacterial communities' growth dependency on substrate concentration outside the lag phase Liu (2020). However, we employ it here to express the intuitive concept that the consumption rate will vary depending on the substrate concentration, assuming β(0) = 0, and saturating over a certain level of resource availability β max . More generally, K could be different for each resource; for simplicity, we set it equal for all the resource types.
Similar to many other resource-consumer models, we further assume that resource uptake by the species is, up to a constant, equal to the depletion of the resources Posfai et al (2017); Pacciani-Mori et al (2020). Consequently, α(r i ) in Eq. 2 of the main text also becomes a Monod function, differing from β(r i ) only on the proportionality constant α max : α(r i ) = α max ri K+ri . Summarizing, it follows that the equation for the system dynamics now becomes Concretely, we explore the scenario in which η → ∞ and varied l and K. The choice of K affects both the resource depletion and the species growth times-scales. The dynamics in (R, N ) space does not display limit-cycles as observed for the linearized resource dynamics defined by the systems described in Eq. 1 and 2 (see Fig. 2). However, in the Monod-version of SMA, the attractors still lie along a vertical line, now defined by N * = αmax i γi δβmax , obtained by the total abundances dynamics stationary solution. This can be observed in Fig. 10a for l = 2 and Fig. 11a for l = 100.
The long term stationary state reached is also reproduced by the Monod version of SMA. In Fig. 10b and Fig. 11b for both cases l = 2, 100, we see that independently of K, all the systems reach a long time stationary value. In fact, the S a associated to l = 2 and η → ∞ reach a plateau. This might result from the fact that when a Monod function describes both consumption and growth rates, the systems with low ν are still computable, while in the linearized version of the model, a low ν value leads to a too large number of species for the system to remain computationally tractable (but would presumably otherwise reach equilibrium). We see that, a long time plateau is reached for many parameter choices and functional implementations of SMA, strongly suggesting that this feature is a robust property of the SMA model proposed. Of course, the obtained plateau values do depend on the choice of model details.

The role of the resource norm
To study the dynamics along the vertical line of attractors, as observed in Fig. 11a, we have to make a connection with the results in Caetano et al (2021). Therefore, let us assume that the strategy vectors are normalized in L p in stead of just L 2 . Furthermore, let us consider two extreme options for the structure of s j : (i) all species specialize into consuming one resource and (ii) all species have an identical strategy vector s * with nonzero components.
For notational convenience, we write g(r i ) for the Monod function and drop the max subscript in α and β. This choice will also highlight that the following results do not depend explicitly on the Monod function, but work for any choice of g that is monotone and passes through 0.
From Posfai et al (2017), we know the strategy s * in L 1 , so inspired by their proof, we make the following computation. From Eq. 8b we find that, in equilibrium, we must have We can rewrite this as We can read this equation as the inner product of the known vector (s * ) p with the unknown vector (s * ) 1−p αg(r * ) δ . As we have just one equation to define this l-dimensional vector, we have an l − 1 dimensional solution space. To close the equation, we pick the solution that minimizes the L p norm. By definition, i (s * i ) p = 1, so the most straightforward solution is for every component to be unity: For p > 1, it is also the unique positive solution that minimizes the L p norm. From Eq. 8a, we deduce that in equilibrium we have the equality so when we fill in the value of N * and the result from Eq. 11, we find that This is identical to the strategy found in Caetano et al (2021).
In the simulations, such as Fig. 11, we noted that the dynamics leads to the minimization of R * , so let us investigate which strategy minimizes R * . When we invert g(r * i ) in Eq. 12, we can find R * by summing over the individual r * i and find which is in line with the numerical results. From this equation we learn two things. First, when we choose our normalization in L 1 , all r * i will become independent of γ, explaining the results from Posfai et al (2017). Second, when p < 1, we see that all r * i become smaller when we choose specialization as the preferred strategy. This is again in line with the results from Caetano et al (2021).

Limitations
After studying several extensions, we must look at the limitations of our approach. We noted that the deterministic version of our model can have an arbitrary large number of coexisting species in absence of an extinction Fig. 11 (R, N ) phase space of the entire evolution of exemplifying realizations (different colors), one for each system with ν = 3.8, η → ∞, l = 100 and different values of K(arrows). Both consumption and species growth rates are described by a Monod function. Here, αmax = 0.5 and βmax = 1. We notice a fast evolution towards the attractor. (b) Examples of living species curve, one for each system with different K, as described in (a). All the systems reach a long time stationary value.
threshold Posfai et al (2017). Unfortunately, this is not a generic property of the model. Indeed, when we perturb the constraint that all species have their strategy normalized to the same value, i.e. by letting the value of the normalization depend on the species, we observe that the complex communities collapse back towards a number of species that is in line with the competitive exclusion principle. It is however very well possible that in, for example, plankton communities, many species have similar functional traits Borics et al (2021); Graco-Roza et al (2021). Therefore, our model might not be the full description of a complex ecosystem, but it can give us insight in the dynamics in a single 'cluster' of similar traits.

Conclusions
We showed that adding speciation in a MacArthur model adds a host of dynamical features reminiscent of commonly observed evolutionary biology. Even when one starts the dynamics with a single species, by introducing new species of slightly different type into an existing ecosystem, we observe equilibration to a maximum number of species on a time scale that is a simple function of the spawning rate. We observe that the system self-maximizes the number of coexisting species, reaching a long-term stationary value. The stationary behavior represents a dynamic equilibrium as attractor in (R, N ) space. Parameters that set the stochastic strength allow the model to explore both invasive and evolutionary dynamics; the size of the resource pool affects the dynamical ability to converge to a niche community Fisher and Mehta (2014). Community aggregate behavior is also stable under perturbations: the system adjusts its features after a resource influx shock; rank abundance plots are in line with commonly observed features. Much phenomenology is robust when different choices are made for resource consumption rates and other model features and some analytical features of the model are consistent with literature results. The perspective embedded in SMA and the range of biologically relevant phenomena it produces offer a flexible interpretation of the term "species" that gives a simple computational tool for a more quantitative understanding of evolution and ecology. One significant open question for this framework is whether the completely random speciation introduced here can also result in the emergence of "clusters" of similar species vectors Maynard et al (2018). Such clustering of species that consume resources in a complementary way would be a computationally tractable representation of a true "tangled bank".

Declarations
• Funding -The authors thank Wageningen University for supporting this work, and declare that no other funds, grants, or other external support was received for the preparation of this manuscript. • Competing interests -The authors have no relevant financial or non-financial interests to disclose. • Availability of data, code and materials -Representative code used to generate the results will be deposited in a public repository. Data and materials used for the study are available. • Author contributions -All authors developed variants of the numerical codes used, performed simulations and interpreted data. E.B. and C.H. derived the main mathematical results, while J.A.D. conceived of the study. All authors contributed to writing the manuscript.