Diese Website nutzt Technologien wie Cookies, um wesentliche Funktionen der Website zu ermöglichen, sowie für analyse, personalisierung und gezielte werbung. Sie können Ihre Einstellungen jederzeit ändern oder die Standardeinstellungen übernehmen. Sie können dieses Banner schließen, um nur mit essenziellen Cookies fortzufahren. Datenschutz-Bestimmungen

Speichereinstellungen

- Gezielte Werbung
- Personalisierung
- Analyse

- Go Mobile »
- Access by Staats- und Universitaetsbibliothek Bremen

# Ecological communities with Lotka-Volterra dynamics

Guy Bunin

- Department of Physics, Technion-Israel Institute of Technology, Haifa 32000, Israel

PDF
Share

- X
- Facebook
- Mendeley
- LinkedIn
- Reddit
- Sina Weibo

Phys. Rev. E 95, 042414 – Published 28 April, 2017

DOI: https://doi.org/10.1103/PhysRevE.95.042414

274

## Abstract

Ecological communities in heterogeneous environments assemble through the combined effect of species interaction and migration. Understanding the effect of these processes on the community properties is central to ecology. Here we study these processes for a single community subject to migration from a pool of species, with population dynamics described by the generalized Lotka-Volterra equations. We derive exact results for the phase diagram describing the dynamical behaviors, and for the diversity and species abundance distributions. A phase transition is found from a phase where a unique globally attractive fixed point exists to a phase where multiple dynamical attractors exist, leading to history-dependent community properties. The model is shown to possess a symmetry that also establishes a connection with other well-known models.

### Physics Subject Headings (PhySH)

- Ecological population dynamics
- Populations
- Cavity methods

## Article Text

Natural ecological conditions are often spatially heterogeneous, with conditions varying between the different local habitats. Individuals migrate between these habitats and interact locally within each habitat (see, e.g.,
). As a result of these processes, local ecological communities assemble (see Fig.
). These processes can be modeled on many levels; one of the simplest and most popular is where a single community is described explicitly, and the rest of the ecosystem is modeled via a regional pool of species from which individuals of all species can invade
.

FIG. 1.

Communities assembled in different locations from a regional species pool. The model focuses on one such community.

Characterizing ecological communities is a central subject in ecology. Some of the key questions are: Does a given habitat reach a stable composition of species? If so, does it depend on historical factors such as the initial conditions following an environmental change, or the order of colonization by different species? How many species coexist in a community? What is the distribution of abundances (number of individuals) of the species in the community? How much do the abundances fluctuate in time? The purpose of this paper is to address these questions within one popular setting.

To proceed, a model for the dynamics of the populations needs to be specified. Here we choose the well-known generalized Lotka-Volterra equations. As for many models, these equations include many system parameters that encode the interactions between all pairs of species. Since the details of the interactions between all pairs of species are typically not available, and for many purposes not needed, here the system parameters are replaced with random numbers, drawn from distributions characterized by a few model parameters. Ever since the pioneering work of May
, models with random parameters have played an important role in theoretical ecology. However, in contrast to Ref.
, here it is the properties of the species pool that are drawn at random, rather than the community, and the composition of the community results from the dynamics. This approach can be viewed as a way of studying the outcomes of migration and species interactions on the community, allowing us to disentangle these processes from other factors such as evolution.

The resulting communities are thus an outcome of the heterogeneous interactions, with the fates of different species intertwined in complex ways. As such, analytical calculations of their properties are not trivial. Tools from physics of disordered systems (in particular, spin glasses) are ideally suited to address such problems, as has been realized many years ago
. Works
that followed in the footsteps of Ref.
described the population dynamics using the replicator equations, which are commonly used in game theory and other fields
. In describing species interactions, the Lotka-Volterra equations are very popular
but have only been analyzed with these tools in
. However, due to the model assumptions in
, some of the main phenomena discussed here do not show up, including the multiple attractor phase and partial coexistence of species. Recently, related techniques have been used to study other models
, as part of a revival in applying methods from statistical mechanics to understand ecology
. Motivated by this renewed interest, we here attempt to provide derivations that are as self-contained and elementary as possible.

The paper is organized as follows: The model is presented in Sec.
, and its phase diagram is described in Sec.
. In Sec.
the model is shown to posses a symmetry by which the properties of its fixed points depend only on certain parameter combinations. In Sec.
the diversity and species abundance distributions are derived, and the transition line to the multiple attractor phase is discussed in Sec.
. The paper concludes with a discussion of the implications of the results to ecology and possible future research directions.

### I. MODEL DEFINITION

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

𝑑⁢𝑁 𝑖 𝑑⁢𝑡=𝑟 𝑖 𝐾 𝑖⁢𝑁 𝑖⁢(𝐾 𝑖−𝑁 𝑖−∑𝑗,(𝑗≠𝑖)𝛼 𝑖⁢𝑗⁢𝑁 𝑗).
(1)

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
.

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

𝛼 𝑖⁢𝑗=𝜇 𝑆+𝜎⁢𝑎 𝑖⁢𝑗,
(2)

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.
and
for 𝑆=1 5 and communities down to 6 species). The rationale behind the scaling in Eq.
is to keep the interaction term in Eq.
, ∑𝑗,(𝑗≠𝑖)𝛼 𝑖⁢𝑗⁢𝑁 𝑗, finite when 𝑆 is large. This has been employed in many works
, while others have considered the limit where the distribution of 𝛼 𝑖⁢𝑗 is kept constant as 𝑆 grows
. 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.

### II. PHASE DIAGRAM FOR THE DYNAMICAL BEHAVIOR

We study the properties of the community at long times. The dynamics of Eq.
can exhibit various different behaviors. Equation
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
. 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
. 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.
and
. For small values of 𝑆 the sharp transitions are broadened, as discussed below.

### III. A SYMMETRY OF THE FIXED POINTS

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

0=𝑛 𝑖⁢(𝜆 𝑖−𝑢⁢𝑛 𝑖−∑𝑗≠𝑖 𝑎 𝑖⁢𝑗⁢𝑛 𝑗+ℎ),
(3)

where 𝑛 𝑖 are the normalized abundances

𝑛 𝑖=𝑁 𝑖/(1 𝑆⁢𝑆∑𝑗=1 𝑁 𝑗)

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

𝑢=1−𝜇/𝑆 𝜎,𝜆 𝑖=𝐾 𝑖−1 𝜎⁡⟨𝑁 𝑖⟩,ℎ=1/⟨𝑁 𝑖⟩−𝜇 𝜎,
(4)

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.
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.
depends on the original parameters only through the combinations 𝑢,𝛾,𝜎 2 𝜆 in Eq.
, 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.
and Fig.
]. In particular, note that the transition between the UFP and MA phases in Fig.
is a straight line. This is due to a symmetry discussed in the next section, as is manifested for 𝜎 𝐾=0.

FIG. 2.

The model exhibits three distinct dynamical behaviors, depending on model parameters. In the unique fixed-point (UFP) phase, a unique, stable fixed point that is resistent to invasion exists for any system. In multiple attractor (MA) phase, multiple dynamical attractors can be found for any system. These may be stable fixed points or other attractors, and the community composition is history dependent. In the third phase abundances grow without bound; here the Lotka-Volterra equations likely break down. The phases are shown for asymmetric interactions ( 𝛾=0) and equal carrying capacities, all set to 𝐾 𝑖=1.

The mapping also establishes a connection to the replicator equations and literature that studies its properties, as the fixed points of Eq.
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
, 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
, ⟨𝛼 𝑖⁢𝑗⟩=1,𝜎=0.

### IV. DIVERSITY AND SPECIES ABUNDANCE

To study the properties of the community, we use a variant of the cavity method
. It is based on the dynamical cavity method
, which does not require 𝛼 𝑖⁢𝑗 to be symmetric but replaces its generating functional formalism by a more elementary derivation, closer in spirit to
. Given the mapping described in Sec.
, it is no wonder that the final equations will be similar to those obtained in the context of the replicator equations
, 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.
, 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.
is extended to include additional small perturbations 𝜉 𝑖 to the 𝜆 𝑖 of each species, later set to zero:

0=𝑛 𝑖⁢(𝜆 𝑖−𝑢⁢𝑛 𝑖−∑𝑗≠𝑖 𝑎 𝑖⁢𝑗⁢𝑛 𝑗+ℎ+𝜉 𝑖).
(5)

The response to the perturbation is defined by

𝑣 𝑖⁢𝑗≡[𝜕𝑛 𝑖/𝜕𝜉 𝑗]𝜉 𝑗=0.
(6)

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

𝑢⁢𝑛 𝑖=𝜆 𝑖−∑𝑗≠𝑖 𝑎 𝑖⁢𝑗⁢𝑛 𝑗+ℎ−𝑎 𝑖⁢0⁢𝑛 0.

This is the same as Eq.
, 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

𝑛 𝑗=𝑛 𝑗\0−𝑛 0⁢∑𝑘 𝑣 𝑗⁢𝑘⁢𝑎 𝑘⁢0.

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

𝑛+0≡𝜆 0−∑𝑗 𝑎 0⁢𝑗⁢𝑛 𝑗\0+ℎ+𝜉 0 𝑢−∑𝑗,𝑘 𝑣 𝑗⁢𝑘⁢𝑎 0⁢𝑗⁢𝑎 𝑘⁢0.
(7)

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
). The sum over the 𝑗=𝑘 terms in ∑𝑗,𝑘 𝑣 𝑗⁢𝑘⁢𝑎 0⁢𝑗⁢𝑎 𝑘⁢0 gives

⟨∑𝑗 𝑣 𝑗⁢𝑗⁡𝑎 0⁢𝑗⁢𝑎 𝑗⁢0⟩=∑𝑗 𝑣 𝑗⁢𝑗⁡⟨𝑎 0⁢𝑗⁢𝑎 𝑗⁢0⟩=𝛾 𝑆⁢∑𝑗 𝑣 𝑗⁢𝑗,

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.
, the term 𝜆 0−∑𝑗 𝑎 0⁢𝑗⁢𝑛 𝑗\0+ℎ has mean ℎ and variance 𝜎 2 𝜆+∑𝑗⟨𝑎 2 0⁢𝑗⟩⁢⟨𝑛 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.,
), and so the numerator is Gaussian, ℎ+𝜉 0+√𝑞+𝜎 2 𝜆 0⁢𝑧 with 𝑃⁡(𝑧)=1√2⁢𝜋⁢𝑒−𝑧 2/2. Setting 𝜉 0=0, Eq.
becomes

𝑛+0=1 𝑢−𝛾⁢𝑣⁢(ℎ+√𝑞+𝜎 2 𝜆⁢𝑧).
(8)

From the Lotka-Volterra equations, Eq.
, 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=m a x⁡(0,𝑛+0) with 𝑛+0 given in Eq.
. 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,

𝑛=m a x⁡⎛⎜ ⎜ ⎜⎝0,ℎ+√𝑞+𝜎 2 𝜆⁢𝑧 𝑢−𝛾⁢𝑣⎞⎟ ⎟ ⎟⎠.
(9)

The distribution of 𝑛 is therefore a truncated Gaussian (as was obtained in other models
).

It remains to find the values of 𝑞,𝑣,ℎ,𝜙. Using Eq.
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⁢𝜋⁢𝑒−𝑧 2 2⁢(𝑧+Δ)𝑘 𝑑 𝑧, these three relations read

𝜙=𝑤 0⁡(Δ),𝑢−𝛾⁢𝑣=√𝑞+𝜎 2 𝜆⁢𝑤 1⁡(Δ),(𝑢−𝛾⁢𝑣)2=(1+𝜎 2 𝜆/𝑞)⁢𝑤 2⁡(Δ),
(10)

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

𝑣=⟨𝑣 0 0⟩=𝜙⁢1 𝑢−𝛾⁢𝑣.
(11)

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

ℎ=𝑞⁢[𝑢−𝑣⁢(1+𝛾+𝜎 2 𝜆/𝑞)].
(12)

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

These equations can be solved numerically by evaluating 𝑞=𝑤 2/𝑤 2 1,𝑣=𝑤 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.
and ⟨𝑁 2 𝑖⟩=𝑞⁢⟨𝑁 𝑖⟩2=𝑞/(𝜎⁢ℎ+𝜇)2. The species abundance of 𝑁 𝑖 is a truncated Gaussian from Eq.
, fully characterized by ⟨𝑁 𝑖⟩ and ⟨𝑁 2 𝑖⟩ [see Fig.
]. It also shows the fraction of persistent species and the moments ⟨𝑁 𝑖⟩ and ⟨𝑁 2 𝑖⟩.

FIG. 3.

Diversity and species abundance distribution. 𝜙 is the fraction of persistent ( 𝑁 𝑖>0) species. Solid and dashed lines show analytical predictions, exact in the UFP phase, left of the vertical dotted lines in (a). Parameter values are 𝛾=0 and 𝜎 2 𝐾=0, and in (b) 𝜇=4 and 𝜎=1.1.

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

FIG. 4.

Moments of the species abundance distribution. Symbols and parameter values are the same as in Fig.
.

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

(2) The analytical predictions are indeed exact for large 𝑆 in the UFP phase, left of the vertical dotted line in Fig.
and Fig.
. 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 𝑆=1 5, 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.
].

The transition to the unbounded growth phase, shown in Fig.
, 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.
, with crosses indicating the position for the transition found in simulations.

FIG. 5.

(a) Phase diagram for 𝜎 𝐾=0 and different 𝛾. The 𝛾=0 lines correspond to Fig.
in the main text. Crosses mark the transition to diverging solutions, found numerically. (b) Numerical check of the phase boundary between UFP and MA phases, the fraction of systems for which there are multiple solutions, at 𝜇=4 and 𝛾=1. At large 𝑆 this fraction jumps sharply at the phase transition. The dashed line marks the analytically calculated transition point, 𝜎=1/√2.

### V. STABILITY AND MULTIPLE ATTRACTORS

As in similar models
, the transition to the MA phase proceeds via the loss of stability of the single fixed point which was described in the previous section
. 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.
, 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
. One finds that

⟨(𝛿⁢𝑛)2⟩⟨𝜉 2⟩=𝜙(𝑢−𝛾⁢𝑣)2−𝜙.
(13)

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
). The transition lines for different values of 𝛾 are shown in Fig.
. The transition lines for different 𝜎 𝐾>0 are presented in Fig.
. 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.
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.
, for different initial conditions on the same system. The fraction of systems that have multiple fixed points is plotted in Fig.
. For large 𝑆 this fraction jumps abruptly at the predicted position. The transition is broadened for finite values of 𝑆 and is quite broad even for 𝑆=1 0 0, 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
, and is also similar to transitions described in
.

FIG. 6.

Phase boundary between UFP and MA phases for different values of 𝜎 𝐾 at 𝛾=0. Solid lines correspond to the 𝛾=0 phase boundaries in Fig.
.

The calculation in Eq.
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.
around ⃗𝑁*, one finds that the matrix 𝑀*𝑖⁢𝑗=𝛼*𝑖⁢𝑗⁢𝑟*𝑖⁢𝑁*𝑖/𝐾*𝑖 must be positive definite, which does not generally follow from 𝛼* being positive definite
. 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).

### VI. DISCUSSION

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
. This may happen for dynamical attractors
, or if species that go below some abundance cutoff are removed from the community, as seems inevitable in any realistic situation
. We see this phenomenon in numerical simulations (see Appendix
). 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
. 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.
and
).

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
. 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
, especially if the species are very similar in their properties) and the full spatial effects (as in meta-community models
). It has also been suggested
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 𝑆=1 0 0.

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
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
. In models of this class that use Lotka-Volterra dynamics, the matrix 𝛼 is constructed from an underlying niche space
. 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
. These models may exhibit a transition from a single fixed-point regime to one where a subspace of marginally stable fixed points exists
; this is different from the MA phase, where multiple fully stable fixed points can appear [see Fig.
].

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

### ACKNOWLEDGMENTS

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.

### Appendices

### APPENDIX A: TRANSITION TO MULTIPLE ATTRACTORS

The stability of 𝛼* can be analyzed through the changes in 𝑛 𝑖 when a perturbation ⃗𝜉 is applied. Consider an on-site perturbation ⃗𝜉=ɛ⁢⃗𝜂 to all species with 𝑛 𝑖>0, where the 𝜂 𝑖 are independent random variables, which are also independent of the 𝑛 𝑖 and 𝑎 𝑖⁢𝑗 and with ⟨𝜂 𝑖⟩=0. ɛ is a parameter, and we wish to calculate 𝑑⁢𝑛+0/𝑑⁢ɛ. We see that it diverges at some point, signaling the loss of the stability of the unique fixed point and a transition to a phase with multiple fixed points.

The average over realizations of the existing variables, 𝑛 𝑗\0 and 𝜂 𝑖>0, is denoted by ⟨..⟩; the average also over 𝜂 0 is denoted by ‾‾‾‾‾(..). As in the other cavity derivations,

𝑛+0=1 ˆ 𝑢⁢(𝜆 0−∑𝑗 𝑎 0⁢𝑗⁢𝑛 𝑗\0+ℎ+ɛ⁢𝜂 0),

where ˆ 𝑢≡𝑢−𝛾⁢𝑣. Differentiating with regard to ɛ, and denoting 𝑦 𝑖\0≡𝑑⁢𝑛 𝑖\0/𝑑⁢ɛ and 𝑦+0≡𝑑⁢𝑛+0/𝑑⁢ɛ,

𝑦+0=1 ˆ 𝑢⁢(−∑𝑗 𝑎 0⁢𝑗⁢𝑦 𝑗\0+𝜂 0).

If 𝑛 0>0, then 𝑑⁢𝑛 0 𝑑⁢ɛ=𝑑⁢𝑛+0 𝑑⁢ɛ=𝑦+0; otherwise 𝑑⁢𝑛 0 𝑑⁢ɛ=0. Taking the square of this equation,

ˆ 𝑢 2⁢(𝑦+0)2=∑𝑗,𝑗′𝑎 0⁢𝑗⁢𝑎 0⁢𝑗′⁢𝑦 𝑗′⁢\0⁢𝑦 𝑗\0−2⁢𝜂 0⁢∑𝑗 𝑎 0⁢𝑗⁢𝑦 𝑗\0+𝜂 2 0.

The term ∑𝑗,𝑗′𝑎 0⁢𝑗⁢𝑎 0⁢𝑗′⁢⟨𝑦 𝑗′⁢\0⁢𝑦 𝑗\0⟩𝜂 𝑖>0 is self-averaging and equal to 𝑆−1⁢∑𝑗≥1⟨𝑦 2 𝑗\0⟩𝜂 𝑖>0=𝜙⁢⟨𝑦 2 𝑗\0⟩+, where ⟨..⟩+ denotes the average only over the community (i.e., 𝑛 𝑗\0>0) variables. To see this, note that 𝑦 2 𝑖\0 is 𝑂⁡(𝑆 0) and with positive mean, while 𝑦 𝑗\0 are weakly correlated for different 𝑗, and with ⟨𝑑⁢𝑛 𝑗\0 𝑑⁢ɛ⟩𝜂 𝑖>0=0 (as ⃗𝜂 𝑖>0→−⃗𝜂 𝑖>0 leads to 𝑑⁢𝑛 𝑗\0 𝑑⁢ɛ→−𝑑⁢𝑛 𝑗\0 𝑑⁢ɛ). Therefore ∑𝑗≠𝑗′𝑎 0⁢𝑗⁢𝑎 0⁢𝑗′⁢⟨𝑦 𝑗′⁢\0⁢𝑦 𝑗\0⟩𝜂 𝑖>0=𝑆⁢𝑂⁡(𝑆−3/2), while ∑𝑗 𝑎 2 0⁢𝑗⁢⟨𝑦 2 𝑗′⁢\0⟩𝜂 𝑖>0=𝑂⁡(1). Thus

ˆ 𝑢 2⁢(𝑦+0)2=𝜙⁢⟨𝑦 2 𝑗\0⟩+−2⁢𝜂 0⁢∑𝑗 𝑎 0⁢𝑗⁢𝑦 𝑗\0+𝜂 2 0.

Averaging over 𝜂 0, which is independent from both 𝑎 0⁢𝑗 and 𝑦 𝑗\0,

ˆ 𝑢 2⁢‾‾‾‾‾‾(𝑦+0)2=𝜙⁢⟨𝑦 2 𝑗\0⟩++‾‾‾𝜂 2 0,

we see that ‾‾‾‾‾‾(𝑦+0)2 is independent of 𝑛+0. Thus ‾‾‾‾‾‾(𝑦+0)2 is precisely the average over 𝑑⁢𝑛 0 𝑑⁢ɛ in all cases where 𝑛 0>0. If 𝑛 0>0, once it has been added to the community it is no different from all the other community variables, so we must have ‾‾‾‾‾‾(𝑦+0)2=⟨𝑦 2 𝑗\0⟩+. Also, ⟨𝑦 2 𝑗\0⟩+=⟨𝑦 2 𝑗⟩+ because the differences between 𝑦 𝑗\0 and 𝑦 𝑗 are 𝑂⁡(𝑆−1/2). Together this gives

⟨𝛿⁢𝑛 2⟩⟨𝜉 2⟩=𝜙⁢⟨𝛿⁢𝑛 2⟩+⟨𝜉 2⟩=𝜙⁢⟨𝑦 2 𝑗\0⟩+‾‾‾𝜂 2 0=𝜙 ˆ 𝑢 2−𝜙.

This diverges when

(𝑢−𝛾⁢𝑣)2=𝜙=𝑣⁢(𝑢−𝛾⁢𝑣)⇒𝑢 𝑐≡(1+𝛾)⁢𝑣.

Thus when 𝑢<𝑢 𝑐 the cavity solution becomes unstable.

For 𝜎 𝜆=0, Δ 𝑐=𝑞⁢[𝑢 𝑐−𝑣 𝑐⁢(1+𝛾)]=0, so 𝑣 𝑐⁢(𝑢 𝑐−𝛾⁢𝑣 𝑐)=∫∞0 𝐷⁢𝑧=1/2, and (𝑢 𝑐−𝛾⁢𝑣 𝑐)2=∫∞−Δ 𝐷⁢𝑧⁢(𝑧+Δ)2=1/2 and so together 𝑣 𝑐=1/√2; hence

𝑢 𝑐=(1+𝛾)/√2(for 𝜎 𝜆=0).

### APPENDIX B: NUMERICAL SIMULATIONS

To numerically find persistent solutions, the variables 𝛼 𝑖⁢𝑗 and 𝐾 𝑖 are first sampled. 𝛼 𝑖⁢𝑗 are sampled from a normal distribution throughout. A uniform distribution was checked to give identical results at large 𝑆. The Lotka-Volterra dynamics, Eq.
, are then integrated, using a Runge-Kutta 45 solver, from random initial conditions sampled uniformly on [0,1]. All species that go below an abundance cutoff 𝑁 𝑖<1 0−1 4 are removed from the community ( 𝑁 𝑖 set to zero). The solver is terminated when an equilibrium solution is found, in which for every 𝑖 either 𝑑⁢𝑁 𝑖/𝑑⁢𝑡 is small or 𝑁 𝑖<1 0−1 4. Solutions that do not terminate are stopped after a long time ( 𝑇=1 0 7) and all variables with 𝑁 𝑖>1 0−1 4 are considered part of the community. The solution is checked against invasion of the pool species not present. As some species are removed during the dynamics due to the abundance cutoff, it is possible that they would be able to invade later. If any such species are found, the dynamics are run from the end point of the first simulation with additional small abundance ( 𝑁 𝑖=1 0−1 0) to the species that may invade. This process is repeated until an uninvadable solution is reached or after ten iterations. In the UFP phase the resulting community is always found to be uninvadable and is usually reached on the first run of the dynamics. For a given system all initial conditions give the same final community. In the MA phase and for asymmetric interactions (specifically 𝛾=0), this process did not always converge to an uninvadable solution after ten iterations and then was stopped (see Fig.
). All numerical results shown in the paper show only minor differences when plotted after the first run, compared to iterations of the invasion process.

To test if a system has more than one fixed point, as shown in Fig.
, the same system (same 𝛼 𝑖⁢𝑗 and 𝐾 𝑖) is run starting from different initial conditions.

FIG. 7.

Fraction of species that may invade. Data plotted with dashed lines show the results after multiple invasion attempts.

## References (40)

1. M. A. Leibold et al., Ecology Lett. 7, 601 (2004).
2. R. H. MacArthur and E. O. Wilson, Theory of Island Biogeography (Princeton University Press, Princeton, NJ, 2015).
3. R. M. May, Nature (London) 238, 413 (1972).
4. S. Diederich and M. Opper, Phys. Rev. A 39, 4333 (1989).
5. H. Rieger, J. Phys. A: Math. Gen. 22, 3447 (1989).
6. M. Opper and S. Diederich, Phys. Rev. Lett. 69, 1616 (1992).
7. K. Tokita, Phys. Rev. Lett. 93, 178102 (2004).
8. K. Tokita, Ecological Informatics 1, 315 (2006).
9. Y. Yoshino, T. Galla, and K. Tokita, Phys. Rev. E 78, 031924 (2008).
10. Y. Yoshino, T. Galla, and K. Tokita, J. Stat. Mech.: Theory Exp. (2007) P09003.
11. T. Obuchi, Y. Kabashima, and K. Tokita, Phys. Rev. E 94, 022312 (2016).
12. J. Hofbauer and K. Sigmund, Evolutionary Games and Population Dynamics (Cambridge University Press, Cambridge, UK, 1998).
13. R. May and A. McLean, Theoretical Ecology: Principles and Applications (Oxford University Press, Oxford, 2007).
14. R. V. Solé and J. Bascompte, Self-Organization in Complex Ecosystems (Princeton University Press, Princeton, NJ, 2006).
15. B. Dickens, C. K. Fisher, and P. Mehta, Phys. Rev. E 94, 022423 (2016).
16. M. Tikhonov and R. Monasson, Phys. Rev. Lett. 118, 048103 (2017).
17. T. Biancalani, L. DeVille, and N. Goldenfeld, Phys. Rev. E 91, 052107 (2015).
18. C. K. Fisher and P. Mehta, Proc. Natl. Acad. Sci. USA 111, 13111 (2014).
19. Y. Fried, D. A. Kessler, and N. M. Shnerb, Sci. Rep. 6, 35648 (2016).
20. D. A. Kessler and N. M. Shnerb, Phys. Rev. E 91, 042705 (2015).
21. S. Azaele, S. Suweis, J. Grilli, I. Volkov, J. R. Banavar, and A. Maritan, Rev. Mod. Phys. 88, 035003 (2016).
22. See, for example, [18, 21], and references therein.
23. J. A. Drake, J. Theor. Biol. 147, 213 (1990).
24. T. J. Case, Proc. Natl. Acad. Sci. USA 87, 9610 (1990).
25. R. Law and R. D. Morton, Ecology 77, 762 (1996).
26. U. Bastolla, M. Lässig, S. C. Manrubia, and A. Valleriani, J. Theor. Biol. 235, 521 (2005).
27. M. Mézard, G. Parisi, and M. A. Virasoro, Europhys. Lett. 1, 77 (1986).
28. A. Crisanti, H. Horner, and H.-J. Sommers, Z. Phys. B: Condens. Matter 92, 257 (1993).
29. This assumes that 𝑣 𝑖⁢𝑗=0 if 𝑛 𝑖=0. A small fraction, of order 1/√𝑆, of the species with 𝑛 𝑗=0 may acquire a positive abundance of order 1/√𝑆, but this effect is negligible.
30. From Eq. (6), the change 𝛿⁢𝑛 𝑖 in response to a perturbation vector ⃗𝜉 is 𝛿⁢𝑛 𝑖=∑𝑗 𝑣 𝑖⁢𝑗⁢𝜉 𝑗. If the elements of ⃗𝜉 are sampled independently, then ⟨(𝛿⁢𝑛)2⟩/⟨𝜉 2⟩=𝑣 2+(𝑆−1)⁢⟨𝑣 2 𝑖⁢𝑗⟩. As long as this is finite [see Eq. (13)], 𝑣 𝑖⁢𝑗 scales as 1/√𝑆.
31. H. Nishimori, Statistical Physics of Spin Glasses and Information Processing: An Introduction (Clarendon Press, Oxford, 2001).
32. J. R. L. De Almeida and D. J. Thouless, J. Phys. A: Math. Gen. 11, 983 (1978).
33. H. Sompolinsky, A. Crisanti, and H.-J. Sommers, Phys. Rev. Lett. 61, 259 (1988).
34. A. Berman and D. Hershkowitz, SIAM J. Algebra. Discr. Meth. 4, 377 (1983).
35. J. A. Capitán, J. A. Cuesta, and J. Bascompte, Phys. Rev. Lett. 103, 168101 (2009).
36. R. V. Solé, D. Alonso, and A. McKane, Philos. Trans. R. Soc., B 357, 667 (2002).
37. R. MacArthur and R. Levins, Am. Nat. 101, 377 (1967).
38. M. Scheffer and E. H. van Nes, Proc. Natl. Acad. Sci. USA 103, 6230 (2006).
39. R. Mac Arthur, Proc. Natl. Acad. Sci. USA 64, 1369 (1969).
40. M. Mézard, G. Parisi, and M. Virasoro, Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications (World Scientific, Singapore, 1987).