Ecological communities with Lotka-Volterra dynamics

The species pool includes 𝑆 species that might reach the local community (see Fig. 1). Within the community the abundances 𝑁𝑖 of the species, with 𝑖=1..𝑆, evolve according to the generalized Lotka-Volterra equations:

Here 𝑟𝑖 are the intrinsic growth rates, 𝐾𝑖 are the carrying capacities, and 𝛼𝑖⁢𝑗 for 𝑖≠𝑗 encode interspecies interactions, with positive values representing competition. Here we focus on the effect of interactions and neglect noise [22].

In order to proceed, the parameters 𝑟𝑖,𝐾𝑖 and 𝛼𝑖⁢𝑗 of the species in the pool—as would be measured in the local community conditions—need to be specified. Here the 𝛼𝑖⁢𝑗 are sampled independently, except for possibly a correlation between 𝛼𝑖⁢𝑗 and 𝛼𝑗⁢𝑖, set by the coefficient 𝛾≡𝑐⁢𝑜⁢𝑟⁢𝑟⁢(𝛼𝑖⁢𝑗,𝛼𝑗⁢𝑖) with −1≤𝛾≤1. It is not restricted to the symmetric case, 𝛾=1.

The analytical technique is controlled at large pool sizes 𝑆 with individually weak interspecies interactions, i.e., keeping the asymmetry 𝛾 and the parameters 𝜇≡𝑆⁡⟨𝛼𝑖⁢𝑗⟩ and 𝜎2≡𝑆⁢𝑣⁢𝑎⁢𝑟⁢(𝛼𝑖⁢𝑗) constant. (The mean will be denoted by angular brackets, ⟨..⟩.) It will be instructive and convenient to work with variables 𝑎𝑖⁢𝑗 defined through

with ⟨𝑎𝑖⁢𝑗⟩=0,⟨𝑎2𝑖⁢𝑗⟩=1/𝑆 and ⟨𝑎𝑖⁢𝑗⁢𝑎𝑗⁢𝑖⟩=𝛾/𝑆, so the relations for 𝛼𝑖⁢𝑗 are satisfied. 𝑎𝑖⁢𝑗 can be sampled from any distribution with the required mean and variance. The carrying capacities 𝐾𝑖 are sampled independently of the 𝛼𝑖⁢𝑗. By rescaling 𝑁𝑖→𝑁𝑖/⟨𝐾𝑖⟩ we set ⟨𝐾𝑖⟩=1. Where 𝐾𝑖 are not identical they are taken to be Gaussian with 𝜎2𝐾≡𝑣⁢𝑎⁢𝑟⁢(𝐾𝑖). (Calculations can be carried out for other distributions of 𝐾𝑖 as well. Note that 𝐾𝑖 is allowed to be negative, as in the case of predators that rely on prey for their survival.)

While the analytical theory is exact for large 𝑆, in practice the number of species need not be large; good agreement is found for modest values of 𝑆 (see Figs. 3 and 4 for 𝑆=15 and communities down to 6 species). The rationale behind the scaling in Eq. (2) is to keep the interaction term in Eq. (1), ∑𝑗,(𝑗≠𝑖)𝛼𝑖⁢𝑗⁢𝑁𝑗, finite when 𝑆 is large. This has been employed in many works [4–10,15,18], while others have considered the limit where the distribution of 𝛼𝑖⁢𝑗 is kept constant as 𝑆 grows [19,20]. One notable difference is the fraction of species that coexist in the community, which attains a finite value in the present scaling and is very small if 𝛼𝑖⁢𝑗 is maintained constant. This may serve as a rough guideline for the scaling limit most appropriate for a specific ecological setting, if 𝑆 and the number of species 𝑆* in the community can be estimated.

We study the properties of the community at long times. The dynamics of Eq. (1) can exhibit various different behaviors. Equation (1) might converge to a fixed point that is stable against small perturbations in values of the persistent species (defined as the species for which 𝑁𝑖>0). This fixed point could also be resistant to invasions, i.e., species for which 𝑁𝑖=0 will decay to zero if introduced at a small abundance, [1𝑁𝑖⁢𝑑⁢𝑁𝑖𝑑⁢𝑡]𝑁𝑖→0+<0. In such a case, the community composition will be maintained indefinitely, if migration acts on long enough time scales which allow the community to relax between colonization attempts [19,23–26]. Other scenarios are also possible: the system may settle to a dynamical attractor, in which the abundances fluctuate indefinitely. In this case uninvadable fixed points might not exist [25]. In addition to the stability and uninvadability conditions just described, all species in the community by definition have positive 𝑁𝑖. This is known as feasibility when formulated as a separate condition.

Given the various possible dynamical behaviors, we first review the phase diagram of the model's dynamical behavior before going into the details of the calculations. Depending on the parameters 𝜇,𝜎,𝜎𝐾, and 𝛾, the model exhibits three distinct dynamical phases, separated by sharp transitions at large 𝑆. In the unique-fixed-point (UFP) phase, any given system admits a unique, stable fixed point that is resistant to invasion.

In the multiple attractor (MA) phase, multiple dynamical attractors generally exist so that the community composition depends on assembly history, for example, the initial conditions or the order of species' invasions. These attractors may be stable fixed points or other dynamical attractors. Finally, in the third phase the abundances grow without bound; here the description in terms of Lotka-Volterra equations eventually breaks down and this regime will not be further discussed.

The analytical technique is based on assuming that an uninvadable fixed point exists, calculating its properties and checking the validity of this assumption self-consistently. Since in the UFP phase all uninvadable fixed points are also stable, the analytical predictions are exact there. The calculation of the transition lines is described in Secs. IV and V. For small values of 𝑆 the sharp transitions are broadened, as discussed below.

We begin with a change of variables that reveals an underlying symmetry of properties of the fixed points. Fixed points 𝑑⁢𝑁𝑖/𝑑⁢𝑡=0 of Eq. (1) require that for all 𝑖, 𝑁𝑖⁢(𝐾𝑖−𝑁𝑖−∑𝑗≠𝑖𝛼𝑖⁢𝑗⁢𝑁𝑗)=0. By using the definition of 𝑎𝑖⁢𝑗 and rearranging this becomes

where 𝑛𝑖 are the normalized abundances

so that 1𝑆⁢∑𝑆𝑖=1𝑛𝑖=1, and

so that ⟨𝜆𝑖⟩=0 and 𝜎2𝜆≡⟨𝜆2𝑖⟩=𝜎2𝐾⁡/(𝜎⁡⟨𝑁𝑖⟩)2. In this language the problem becomes as follows: Given the input parameters 𝑢,𝛾, and 𝜎2𝜆, find a value of ℎ such that Eq. (3) holds for all species; ∑𝑆𝑖=1𝑛𝑖=𝑆; species for which 𝑛𝑖=0 cannot invade; and the persistent species ( 𝑛𝑖>0) are stable against small perturbations in 𝑛𝑖. Thus ℎ is a result of the calculation from which ⟨𝑁𝑖⟩ can be obtained. At large 𝑆, 𝑢≃1/𝜎 to lowest order in 1/𝑆, which is used throughout the paper in comparisons with Lotka-Volterra simulations.

As Eq. (3) depends on the original parameters only through the combinations 𝑢,𝛾,𝜎2𝜆 in Eq. (4), new nontrivial symmetries of the original problem are revealed. For example, if all carrying capacities are identical ( 𝐾𝑖=1), then 𝜎2𝜆=0 and the problem depends on all the original parameters through 𝑢 and 𝛾. Therefore at at large 𝑆 all properties of the normalized abundances 𝑛𝑖 do not depend on 𝜇 [see Fig. 3(a) and Fig. 4(b)]. In particular, note that the transition between the UFP and MA phases in Fig. 2 is a straight line. This is due to a symmetry discussed in the next section, as is manifested for 𝜎𝐾=0.

FIG. 2.

The mapping also establishes a connection to the replicator equations and literature that studies its properties, as the fixed points of Eq. (3) are precisely those of the replicator equations 𝑑⁢𝑛𝑖/𝑑⁢𝑡=𝑛𝑖⁢(𝜆𝑖−𝑢⁢𝑛𝑖−∑𝑗≠𝑖𝑎𝑖⁢𝑗⁢𝑛𝑗+ℎ). This mapping only holds for the fixed-point properties. It is not the well-known mapping of the full dynamics [7,12], in which the change of variables mixes them in a way that generates statistical dependencies between interaction strengths, even when they do not exist in the original Lotka-Volterra equations.

The exact expression for 𝑢, 𝑢=(1−⟨𝛼𝑖⁢𝑗⟩)/𝜎, can be used in a different scaling limit where ⟨𝛼𝑖⁢𝑗⟩ is kept constant, instead of ⟨𝛼𝑖⁢𝑗⟩=𝜇/𝑆. Here 𝑢 can be interpreted as the angle in ⟨𝛼𝑖⁢𝑗⟩,𝜎 plane from the Hubbell point [15,20], ⟨𝛼𝑖⁢𝑗⟩=1,𝜎=0.

To study the properties of the community, we use a variant of the cavity method [4–6,27,28]. It is based on the dynamical cavity method [5,6], which does not require 𝛼𝑖⁢𝑗 to be symmetric but replaces its generating functional formalism by a more elementary derivation, closer in spirit to [27]. Given the mapping described in Sec. III, it is no wonder that the final equations will be similar to those obtained in the context of the replicator equations [4,6], and equivalent for 𝜎𝜆=0. For completeness the derivation provided is self-contained and does not require prior knowledge of spin-glass techniques.

The argument proceeds by adding a new species along with newly sampled interactions with the existing system, and comparing the properties of the solution with 𝑆 species to that with 𝑆+1 species, requiring that the new species has the same properties as the rest. Starting from Eq. (3), assume that the abundances 𝑛𝑖\0 of the species in the pool 𝑖=1..𝑆 are known. Introduce a new species with interactions {𝑎0⁢𝑖,𝑎𝑖⁢0}𝑖=1..𝑆 and 𝜆0. For the purposes of the derivation, Eq. (3) is extended to include additional small perturbations 𝜉𝑖 to the 𝜆𝑖 of each species, later set to zero:

The response to the perturbation is defined by

Once the new species is introduced, it might invade and its final abundance will be 𝑛0>0, or else 𝑛0=0. The effect it has on the species 𝑖≥1 is [29]

This is the same as Eq. (5), with 𝜉𝑖=−𝑎𝑖⁢0⁢𝑛0. For large 𝑆 each 𝑎𝑖⁢0⁢𝑛0 is small, as 𝑛0 does not scale with 𝑆 (its mean is 1 by the normalization), and 𝑎𝑖⁢0 scales as 1/√𝑆. Therefore linear response can be used, 𝑛𝑗=𝑛𝑗\0−∑𝑘𝑣𝑗⁢𝑘⁢𝜉𝑘, giving

If 𝑛0>0 we substitute this equation into 0=𝜆0−𝑢⁢𝑛0−∑𝑗𝑎0⁢𝑗⁢𝑛𝑗+ℎ+𝜉0 and rearrange to find that 𝑛0=𝑛+0, where

The denominator of this equation will be a finite number with negligible fluctuations. To see this, note that 𝑣𝑖⁢𝑖=𝑂⁡(𝑆0), while 𝑣𝑖⁢𝑗, which is mediated by the 𝑎𝑖⁢𝑗 interactions, is expected to be 𝑣𝑖⁢𝑗=𝑂⁡(𝑆−1/2) (as can later be verified [30]). The sum over the 𝑗=𝑘 terms in ∑𝑗,𝑘𝑣𝑗⁢𝑘⁢𝑎0⁢𝑗⁢𝑎𝑘⁢0 gives

with 𝑂⁡(𝑆−1/2) fluctuations, while the sum over the 𝑗≠𝑘 terms is 𝑂⁡(𝑆−1/2). Together, up to 𝑂⁡(𝑆−1/2) fluctuations, the denominator is equal to 𝑢−𝛾⁢𝑣 with 𝑣≡⟨𝑣𝑗⁢𝑗⟩. All in all, the indirect effects of the existing community on the new species changes the denominator from 𝑢 to 𝑢−𝛾⁢𝑣.

Turning to the numerator of Eq. (7), the term 𝜆0−∑𝑗𝑎0⁢𝑗⁢𝑛𝑗\0+ℎ has mean ℎ and variance 𝜎2𝜆+∑𝑗⟨𝑎20⁢𝑗⟩⁢⟨𝑛2⟩=𝜎2𝜆+𝑞, where 𝑞≡⟨𝑛2⟩. This follows from the distributions of 𝜆0 and 𝑎0⁢𝑗 (all independent from each other by construction). As a sum of many weakly correlated terms, −∑𝑗𝑎0⁢𝑗⁢𝑛𝑗\0 is Gaussian (see, e.g., [31]), and so the numerator is Gaussian, ℎ+𝜉0+√𝑞+𝜎2𝜆0⁢𝑧 with 𝑃⁡(𝑧)=1√2⁢𝜋⁢𝑒−𝑧2/2. Setting 𝜉0=0, Eq. (7) becomes

From the Lotka-Volterra equations, Eq. (1), it follows that if 𝑛+0>0, the solution 𝑛0=0 is not stable against invasion ([1𝑁𝑖⁢𝑑⁢𝑁𝑖𝑑⁢𝑡]𝑁𝑖→0+ >0), so 𝑛0=𝑛+0. This is where the resistance to invasion enters. Together, 𝑛0=max⁡(0,𝑛+0) with 𝑛+0 given in Eq. (8). But once species “0” has been added to the system it is in no way different from the other species, so we may drop the subscript 0 to obtain the species abundance distribution of all species,

The distribution of 𝑛 is therefore a truncated Gaussian (as was obtained in other models [6,7]).

It remains to find the values of 𝑞,𝑣,ℎ,𝜙. Using Eq. (9) for 𝑛, the relations 1=⟨𝑛⟩,𝑞=⟨𝑛2⟩,𝜙=⟨Θ+⁡(𝑛)⟩ can be used. Here Θ+⁡(𝑛) is the Heaviside function with Θ+⁡(0)=0, i.e., Θ+⁡(𝑛)=0 if 𝑛<0 and 1 otherwise. Note that 𝜙 is the fraction of persistent species. Denoting 𝑤𝑘⁡(Δ)≡∫∞−Δ1√2⁢𝜋⁢𝑒−𝑧22⁢(𝑧+Δ)𝑘𝑑𝑧, these three relations read

where Δ≡ℎ⁡/√𝑞+𝜎2𝜆0. A fourth equation is obtained by differentiating Eq. (7) with respect to 𝜉0: if 𝑛0>0 it gives 𝑣00=1/(𝑢−𝛾⁢𝑣) and otherwise 𝑣00=0. Together,

This completes the set of four coupled equations for the unknowns 𝑞,𝑣,ℎ,𝜙. Using the identity 𝑤2⁡(Δ)=𝑤0⁡(Δ)+Δ·𝑤1⁡(Δ) and the definition of Δ, we also have

These equations were first derived, for 𝜎2𝜆=0, in the context of the replicator equations in [4,6].

These equations can be solved numerically by evaluating 𝑞=𝑤2/𝑤21,𝑣=𝑤0/(𝑤1⁢√𝑞+𝜎2𝜆0) and 𝑢=𝛾⁢𝑣+√𝑞+𝜎2𝜆0⁢𝑤1 as functions of Δ,𝛾 and 𝜎2𝜆. One then finds that at any fixed values of 𝜎2𝜆,𝛾, and for all 𝑢>0, the function 𝑢 is an increasing function of Δ. Therefore 𝑢⁡(Δ) is one-to-one and can be numerically inverted to obtain Δ⁡(𝑢), allowing one to calculate the values of 𝑣,𝑞,ℎ as a function of the input parameters 𝑢,𝛾, and 𝜎2𝜆. The explicit expressions for 𝑞,𝑣,𝑢,ℎ as functions of Δ,𝛾, and 𝜎2𝜆, along with the uniqueness of the inversion 𝑢⁡(Δ), guarantee that the solution as a function of 𝑢,𝛾, and 𝜎2𝜆 is unique.

Returning to the Lotka-Volterra variables 𝑁𝑖=⟨𝑁𝑖⟩⁢𝑛𝑖, one has ⟨𝑁𝑖⟩=1/(𝜎⁢ℎ+𝜇) from Eq. (4) and ⟨𝑁2𝑖⟩=𝑞⁢⟨𝑁𝑖⟩2=𝑞/(𝜎⁢ℎ+𝜇)2. The species abundance of 𝑁𝑖 is a truncated Gaussian from Eq. (9), fully characterized by ⟨𝑁𝑖⟩ and ⟨𝑁2𝑖⟩ [see Fig. 3(b)]. It also shows the fraction of persistent species and the moments ⟨𝑁𝑖⟩ and ⟨𝑁2𝑖⟩.

FIG. 3.

To simulate the model, a realization of the interaction matrix 𝛼𝑖⁢𝑗 and the carrying capacities 𝐾𝑖 is sampled, and the dynamics in Eq. (1) are then run. The details are given in Appendix B. The predictions for the diversity and abundance distributions are shown in Figs. 3 and 4, and compared to simulations with 𝑆=200 and 𝑆=15. The following points are noteworthy:

FIG. 4.

(1) At large 𝑆 (here, 𝑆=200), data collapses for 𝜇=2,4 for the quantities in Fig. 3(a) and Fig. 4(b), which are properties of the normalized 𝑛𝑖 alone, in accordance with the predictions of Sec. III.

(2) The analytical predictions are indeed exact for large 𝑆 in the UFP phase, left of the vertical dotted line in Fig. 3(a) and Fig. 4. Note that the analytical predictions are also a good approximation beyond this point, well into the MA phase.

(3) The analytical predictions fit quite well also for 𝑆=15, and communities with 𝑆* as low as 6 species. The transition to the unbounded growth phase is broadened at low 𝑆, leading to higher values and large fluctuations in ⟨𝑁2𝑖⟩ [see Fig. 4(b)].

The transition to the unbounded growth phase, shown in Fig. 2, is marked by the divergence of ⟨𝑁𝑖⟩. As ⟨𝑁𝑖⟩=1/(𝜎⁢ℎ+𝜇), the boundary with the unbounded growth can be located analytically through the condition 𝜎⁢ℎ+𝜇=0. The analytical expression for ℎ is exact in the UFP phase and approximate in the MA phase, so the prediction for this phase boundary will accordingly be exact when it limits the UFP phase and approximate when it limits the MA. The phase diagram for different values of 𝛾 is shown in Fig. 5(a), with crosses indicating the position for the transition found in simulations.

FIG. 5.

As in similar models [4,6], the transition to the MA phase proceeds via the loss of stability of the single fixed point which was described in the previous section [32,33]. The transition can thus be located by calculating the stability of this fixed point to perturbations. Consider the change of the normalized abundances ⃗𝑛 in response to a perturbation ⃗𝜉 defined in Eq. (5), with the 𝜉𝑖's sampled independently from each other. (The average includes 𝛿⁢𝑛𝑖=0 for species outside the community.) The derivation is similar to the one in the previous section and to standard techniques. It is given in Appendix A. One finds that

This diverges at the boundary between the UFP and MA phases, 𝜙=(𝑢−𝛾⁢𝑣)2, indicating the loss of stability of the fixed point. For 𝜎𝐾=0 this line lies along 𝜎=√2/(1+𝛾) for all 𝜇>0 (see Appendix A). The transition lines for different values of 𝛾 are shown in Fig. 5. The transition lines for different 𝜎𝐾>0 are presented in Fig. 6. These figures demonstrate that lower symmetry (reducing 𝛾) and higher variability in the carrying capacities (increasing 𝜎𝐾) both delay the transition from the UFP to the MA phase. Beyond this line, the right-hand side of Eq. (13) is negative, which cannot be correct. This means that the assumption of a single uninvadable fixed point is not longer consistent. To observe this transition numerically, one checks whether the same final community is obtained in repeated runs of the Lotka-Volterra dynamics, Eq. (1), for different initial conditions on the same system. The fraction of systems that have multiple fixed points is plotted in Fig. 5(b). For large 𝑆 this fraction jumps abruptly at the predicted position. The transition is broadened for finite values of 𝑆 and is quite broad even for 𝑆=100, with a single fixed point found even well inside the MA phase. The transition between the UFP and MA phases is closely related to those found for the replicator equations [4,6,9,10,18], and is also similar to transitions described in [15,20].

FIG. 6.

The calculation in Eq. (13) is of the response to changes in the carrying capacities 𝐾𝑖. This is equivalent to the requirement that 𝛼*, the interaction matrix reduced to the community species 𝑁𝑖>0, is positive definite. An a priori different requirement is for the fixed point ⃗𝑁* to be dynamically stable, returning to ⃗𝑁* when abundances are initialized close to it. By expanding Eq. (1) around ⃗𝑁*, one finds that the matrix 𝑀*𝑖⁢𝑗=𝛼*𝑖⁢𝑗⁢𝑟*𝑖⁢𝑁*𝑖/𝐾*𝑖 must be positive definite, which does not generally follow from 𝛼* being positive definite [34]. However, in the UFP phase we always find that 𝑀* is positive definite if ⃗𝑟* is sampled independently from 𝛼* and ⃗𝐾*. This was tested for various distributions (including identical values, exponential, power-law, and uniform distributions).

The relationship between the different constraints on a community—feasibility, uninvadability, and stability—is a central theme in theoretical ecology. The phase diagram of the present model gives an interesting take on the subject. In the UFP phase we find that all feasible and uninvadable fixed points are also automatically stable. Moreover, in this phase they are unique and thus globally stable. Upon transitioning to the MA phase, both these properties break at once: there are multiple fixed points and some of them are unstable. In this phase asymmetric models (e.g., 𝛾=0) may fall into dynamical attractors that are not fixed points. In addition, they might never become uninvadable; instead, invasions can trigger jumps between a number of possible communities [25,26,35]. This may happen for dynamical attractors [25], or if species that go below some abundance cutoff are removed from the community, as seems inevitable in any realistic situation [25,26,35]. We see this phenomenon in numerical simulations (see Appendix B). It is perhaps surprising that all these changes happen upon crossing a single transition.

The analytical technique as employed here only uses feasibility and resistance to invasion, without accounting for stability. Stability may be an important factor affecting the structure of communities [3]. Nonetheless, some of the key results, in particular, the diversity (number of coexisting species) and the abundance distributions, are very reasonably described by the theory even in the MA phase where some fixed points are unstable (see Figs. 3 and 4).

The species abundance distribution 𝑃⁡(𝑁) is Gaussian, truncated at 𝑁=0. This analytical prediction is confirmed by simulations. Abundance distributions are a subject of ongoing research; commonly cited distributions are log-normal or power law rather than Gaussian [13,21,36]. As the community assembly studied here accounts for only some of the processes acting in nature, this result is instructive in disentangling the effect of different processes. Factors that may significantly change the abundance distribution are the introduction of noise (as emphasized by neutral theory [21,36], especially if the species are very similar in their properties) and the full spatial effects (as in meta-community models [1]). It has also been suggested [7] that this 𝑃⁡(𝑁) may appear similar to observed distributions in some parameter regimes.

The applicability of the predictions to small values of 𝑆 is important when studying ecosystems with tens of species (or groups of similar species). Here we find that diversity and species abundance distributions are very well approximated by the large 𝑆 results. The transitions between the different phases are broadened, and in particular, unique fixed points as in the UFP phase are found for model parameters that would be deep in the MA phase. This is a significant effect even for 𝑆=100.

The model studied here should be viewed as a simple null model that may be extended or modified to tackle different questions or specific biological settings. The theoretical framework may be readily applied to a range of such settings. The functional forms of the single-species response and interspecies interactions may be replaced with more realistic, biologically motivated forms. Sparse interactions, where 𝛼𝑖⁢𝑗 vanishes for some of the interaction pairs, can also be studied. In this case the theory is almost unchanged once the definitions of 𝜇,𝜎 are modified to 𝜇=𝐶⁡⟨𝛼𝑖⁢𝑗⟩ and 𝜎2=𝐶⁢𝑣⁢𝑎⁢𝑟⁢(𝛼𝑖⁢𝑗), where 𝐶 is the average number of links of a given species. This straightforward extension requires 𝐶 to be large and with small variability between species, in which case the derivations above generalize trivially. In fact, this case is already covered by the present framework since 𝛼𝑖⁢𝑗 was allowed take any distribution, in particular, one where some of the 𝛼𝑖⁢𝑗 values vanish. If 𝐶 is small, new phenomena may appear [11] that can be approached with variants of the present cavity method.

Models of resource competition play an important role in ecology, yet theoretical techniques related to those used here have only begun to be applied [10,16]. In models of this class that use Lotka-Volterra dynamics, the matrix 𝛼 is constructed from an underlying niche space [17,37,38]. As a result, the properties of 𝛼 are different from those taken in the present work. Under certain conditions on this construction, all fixed points are either stable or marginally stable [39]. These models may exhibit a transition from a single fixed-point regime to one where a subspace of marginally stable fixed points exists [17,37,38]; this is different from the MA phase, where multiple fully stable fixed points can appear [see Fig. 5(b)].

Other future research directions include full spatial or meta-community models [1]. Finally, it would be interesting to study the MA phase in more depth. This may require more powerful and involved theoretical tools [33,40].

It is a pleasure to thank M. Barbier, J. Friedman, J. Gore, M. Kardar, D. Kessler, P. Mehta, D. Rothman, and M. Tikhonov for valuable discussions. The support of the Pappalardo Fellowship in Physics is gratefully acknowledged.