Self-organization of ecosystems to exclude half of all potential invaders

I. INTRODUCTION

Biodiverse ecosystems on large spatial scales exhibit such a wide range of reproducible, emergent high-level properties that we might consider them a distinct state of matter, accessible to methods from statistical mechanics. Entropy maximization and field-theoretical methods have been invoked to describe such systems and explain, e.g., widely observed noninteger power laws relating study area to species counts (species richness) .

Such spatiotemporal structures formed by interacting species populations are reproduced in simulations of spatially coupled systems of Lotka-Volterra competition models of the form

(1)

Here $𝑡$ is time, $𝑏 𝑖 \geq 0$ are the time-dependent population biomasses of $𝑆 pool$ species $𝑖$ competing in the region, $𝑠 𝑖$ are population growth rates in absence of competition, and the elements $𝐴 𝑖 𝑗$ of matrix $𝐀$ quantify the strengths of intraspecific ( $𝑖 = 𝑗$ ) and interspecific ( $𝑖 \neq 𝑗$ ) competition. The term in brackets in Eq. approximates the momentary population growth rates of interacting species as a linear function of population biomasses. Rescaling biomass variables, one can set $𝐴 𝑖 𝑖 = 1$ ( $1 \leq 𝑖 \leq 𝑆$ ). The $𝜄 𝑖$ represent spatial coupling (immigration) to neighboring locations. With $𝜄 𝑖 = 0$ , we can distinguish the $𝑆$ species extant ( $𝑏 𝑖 > 0$ ) in an equilibrium and all others ( $𝑏 𝑖 = 0$ ). Equation then has several equilibria, called ecological communities, distinguished by different (index) sets of extant species $ℰ$ , called community compositions.

We consider a weaker interspecific ( $𝑖 \neq 𝑗$ ) than intraspecific competition ( $0 \leq 𝐴 𝑖 𝑗 < 1$ ) and absent or very weak spatial coupling ( $𝜄 𝑖 = 0$ ). Under weak coupling, extant and extinct species can be distinguished using heuristic criteria . Substantial progress was recently made in understanding Eq. when interspecific $𝐴 𝑖 𝑗$ are assigned random values (fixed through time) and species richness $𝑆 = | ℰ |$ is large. We focus on the empirically well-supported case where the $𝐴 𝑖 𝑗$ are independently and identically distributed (i.i.d.) with mean $𝜇$ and variance $𝜎 2$ (though correlated $𝐴 𝑖 𝑗$ and $𝐴 𝑗 𝑖$ are also often considered ). Higher moments play only minor roles in the limit $𝑆 \to \infty$ .

Random matrix theory predicts that for high species richness Eq. becomes an ill-defined problem, a phenomenon called ecological structural instability (, Eq. (18.1)): when

𝑆 > 𝑆 ESI def = (1 - 𝜇) 2 𝜎 2,

(2)

the area in the complex plane covered by the eigenvalues of the matrix ${𝐴 𝑖 𝑗} 𝑖, 𝑗 \in ℰ$ includes zero . Population abundances are then highly sensitive to perturbations such as changes in the $𝑠 𝑖$ or addition or removal of species. This implies more species than $𝑆 ESI$ are unlikely to coexist . For $𝑆 pool < 2 𝑆 ESI$ , the “unique fixed point” phase, populations reach an equilibrium where $𝑆 < 𝑆 ESI$ species coexist, with the rest extinct . At $𝑆 pool = 2 𝑆 ESI$ , where $𝑆 = 𝑆 ESI$ , the system exhibits a phase transition . For $𝑆 pool > 2 𝑆 ESI$ , in the so called “multiple attractors” (MA) phase, system dynamics is only partially understood. There is, with all $𝜄 𝑖 = 0$ , generally no stable equilibrium in this phase, because for each equilibrium there is at least one locally extinct species that could invade (i.e., its linear population growth rate is positive). These invasions are represented by heteroclinic orbits in phase space that each lead from one equilibrium to another . The resulting heteroclinic orbit network contains (as a directed graph) at least one strongly connected component that is a stable attractor of system dynamics. On this attractor, the community transitions between different community compositions at an approximately constant rate on a logarithmic time scale . For fixed small but positive $𝜄 𝑖$ , this network is perturbed, leading to slow switching between different community compositions on a linear time scale.

In the explicitly spatially extended case, so-called metacommunity models , the $𝜄 𝑖$ are dynamically determined by neighboring communities, quantifying the coupling between the local Lotka-Volterra equations . For these models a similar phase transition is observed in simulations as the number of interacting communities increases .

In this paper, we make advances characterizing the statistical mechanics of the MA phase for $𝑆 pool ≫ 𝑆 ESI ≫ 1$ . This is achieved by approximating Eq. by a community assembly model . In this model we construct an ecological community by iterated invasions of random species. Extant species go extinct as a result of competition and a steady state is reached where invasions and extinctions are balanced. We assume, supported by simulations over a wide parameter range, that extinction occurs only to low-biomass species and only after the invasion of a new species. In other words, the changes in community composition and biomasses as a result of a single invasion are small—there are no mass extinctions or sudden changes. The community assembly process, particularly the probability of successful invasion $𝑝 i n v$ , have been the subject of many studies . It is a quantity of practical relevance in conservation ecology, where alien species invasions are considered a problem, and an important characteristic of community structure and dynamics. Our objective is to calculate $𝑝 i n v$ in a steady-state community experiencing continuous species turnover through invasions and extinctions.

II. METHODS AND RESULTS

The state of our model, approximating Eq. , is characterized by the community composition $ℰ$ . It changes by invasion of species $𝑘 \notin ℰ$ , sampled at random conditional to $𝑘$ being able to invade (explained below). Since we are considering the case $𝑆 pool ≫ 𝑆 ESI = 𝒪 (𝑆), 𝑘$ 's interactions $𝐴 𝑗 𝑘, 𝐴 𝑘 𝑗$ ( $𝑗 \in ℰ$ ) are effectively i.i.d., sampled from a distribution with support $0 \leq 𝐴 𝑗 𝑘 < 1$ as above. The immigration terms $𝜄 𝑖$ in Eq. are disregarded in this approximation. Denoting $ℰ + = {𝑘} \cup ℰ$ , population dynamics of the assembly model between invasions is given by

𝑑 𝑏 𝑖 𝑑 𝑡 = [𝑠 𝑖 - \sum 𝑗 \in ℰ + 𝐴 𝑖 𝑗 𝑏 𝑗] 𝑏 𝑖 (𝑖 \in ℰ +),

(3)

with initial conditions such that the species in $ℰ$ are in equilibrium and $𝑘$ has small, positive abundance ( $0 < 𝑏 𝑘 ≪ 1 / | ℰ |$ ). From there, the system typically reaches a new unique stable equilibrium [found, e.g., by simulating Eq. ]. After this, the set $ℰ$ of extant species is updated accordingly and the process is repeated, leading to a stochastic steady state where invasion and extinctions balance each other on average.

Our calculation makes use of a close relation between two quantities. First, the harvesting resistance of extant species $𝑖$ , defined as the inverse susceptibility $𝑟 𝑖 = (d l n (𝑏 𝑖) / d 𝑠 𝑖) - 1$ , Sec. 14.5. It can be computed as

𝑟 𝑖 = 𝑠 𝑖 - \sum 𝑗 \in ℰ ∖ {𝑖} 𝐴 𝑖 𝑗 𝑏 𝑗,

(4)

with the $𝑏 𝑗$ chosen to satisfy $𝑠 𝑙 - \sum 𝑗 \in ℰ ∖ {𝑖} 𝐴 𝑙 𝑗 𝑏 𝑗 = 0$ for all $𝑙 \in ℰ ∖ {𝑖}$ (see details in the Appendix) , Sec. 14.5. Second, invasion fitness, the linear growth rate of new invaders $𝑘 \notin ℰ$ , which, by comparison of Eqs. and , is given by Eq. with $𝑖 = 𝑘$ . As a general rule, $𝑟 𝑖 > 0$ for all extant species, and $𝑘$ can invade if and only if $𝑟 𝑘 > 0$ . Invasion probability $𝑝 i n v$ is the probability that $𝑟 𝑘 > 0$ for a randomly sampled invader $𝑘$ .

For simplicity, we assume the same value $𝑠$ for all $𝑠 𝑖$ . In this case, $𝑟 𝑖$ and $𝑏 𝑖$ are approximately proportional to each other, such that (, Eq. (17.23))

𝑟 𝑖 𝑏 𝑖 \approx const. \approx (1 - 𝜇) .

(5)

For new, random invaders $𝑖 = 𝑘$ , the $𝐴 𝑘 𝑗$ in Eq. are random. Assuming a sparse matrix $𝐴 𝑖 𝑗$ such that the extant biomasses are weakly correlated, the invasion fitness $𝑟 𝑘$ , Eq. , is therefore for large $𝑆$ normally distributed. Since community turnover leads to slow random changes in invasion fitness, we therefore assume (and verify in the Appendix), that, for a randomly chosen potential invader $𝑘$ that never actually invades, $𝑟 𝑘$ follows an approximate Ornstein-Uhlenbeck (OU) process (, Sec. 14.6) with mean $‾ 𝑟$ and variance $v a r 𝑟$ . We introduce a new, slow time scale $𝑇$ that measures the number of successful invasions in the assembly model and denote the OU relaxation rate on this time scale by $𝜌$ .

Now consider the alternative case where $𝑖$ is one of many ( $𝑆 ≫ 1$ ) extant species. Because $𝑟 𝑖$ is given by the same formula as $𝑟 𝑘$ , Eq. , we assume that between successive invasions the dynamics of $𝑟 𝑖$ follows the same kind of random walk as $𝑟 𝑘$ . At the moment when $𝑖$ invades, $𝑟 𝑖$ follows a normal distribution with mean $‾ 𝑟$ and variance $v a r 𝑟$ , conditional to $𝑟 𝑖 > 0$ . Its extinction occurs when $𝑟 𝑖$ reaches zero.

The relaxation rate $𝜌$ is therefore determined in two different ways: (i) as the relaxation rate $𝜌$ of the underlying OU process and (ii) as the relaxation rate of the sum in Eq. , of which each term performs, by Eq. , a random walk with the same relaxation rate parameter, while invasions and extinctions lead to removal and addition of terms. This gives a self-consistency condition that we can evaluate to determine $𝑝 i n v$ .

To obtain an explicit expression for the relaxation rate given by (ii), we first make an additional simplifying assumption, later relaxed: with two parameters $0 < 𝐼, 𝐶 < 1$ , the i.i.d. distribution of $𝐴 𝑖 𝑗$ for $𝑖 \neq 𝑗$ is such that $𝐴 𝑖 𝑗 = 𝐼$ with probability $𝐶$ and $𝐴 𝑖 𝑗 = 0$ otherwise. We denote by $𝒜 𝑖$ the subset of indices $𝑗$ in $ℰ$ for which $𝐴 𝑖 𝑗 = 𝐼$ . We require that $𝐶$ is sufficiently small such that $| 𝒜 𝑖 |$ follows an approximate Poisson distribution.

We standardize the variance of the OU process driving harvesting resistance by introducing scaled variables $𝑦 𝑖 = 𝑟 𝑖 / (v a r 𝑟) 1 / 2$ and defining $‾ 𝑦 = ‾ 𝑟 / (v a r 𝑟) 1 / 2$ . Hence, $𝑝 i n v = Φ (‾ 𝑦)$ , with $Φ$ denoting the cumulative normal distribution. By Eq. , the scaled harvesting resistances $𝑦 𝑗$ differ from the population biomasses $𝑏 𝑗$ only by a constant factor. The relaxation rate of the sum in Eq. therefore equals the relaxation rate of

𝑅 𝑖 = \sum 𝑗 \in 𝒜 𝑖 𝑦 𝑗,

(6)

where the OU process followed by each $𝑦 𝑗$ from the time it invades until it becomes extinct can be expressed through a Langevin equation

d 𝑦 𝑗 (𝑇) d 𝑇 = - 𝜌 (𝑦 𝑗 (𝑇) - ‾ 𝑦) + \sqrt 2 𝜌 𝜉 𝑗 (𝑇),

(7)

with $𝜉 𝑗 (𝑇)$ denoting uncorrelated unit white noise $[c o r (𝜉 𝑖 (𝑇), 𝜉 𝑗 (𝑇')) = 𝛿 𝑖, 𝑗 𝛿 (𝑇 - 𝑇')]$ . The coefficient of the noise term follows because the unconstrained OU process has unit variance by construction. The noise is uncorrelated due to our assumption of a sparse matrix $𝐴 𝑖 𝑗$ and weak interspecific interactions. In addition, $𝑅 𝑖$ changes through invasions and extinctions of species $𝑗$ . When a species $𝑗$ invades, its scaled harvesting resistance $𝑦 𝑗$ follows a normal distribution conditional to $𝑦 𝑗 > 0$ , with density

𝑃 i n v (𝑦) = 1 \sqrt 2 𝜋 Φ (‾ 𝑦) e x p (- (𝑦 - ‾ 𝑦) 2 2) .

(8)

Invasions of species $𝑗$ with $𝐴 𝑖 𝑗 > 0$ occur at a rate $𝐶$ , adding the index $𝑗$ to $𝒜 𝑖$ and a term $𝑦 𝑗$ to the sum defining $𝑅 𝑖$ . Extinction of a species $𝑗 \in 𝒜 𝑖$ occurs when $𝑦 𝑗$ reaches zero. The index $𝑗$ is then removed from $𝒜 𝑖$ and the process $𝑦 𝑗$ from $𝑅 𝑖$ .

We now evaluate the relaxation rate of $𝑅 𝑖$ using moment equations, denoting expectation values by brackets $⟨ \cdot ⟩$ . With our assumption of negligible correlations, suppressing the index $𝑖$ and abbreviating $𝑍 = | 𝒜 |$ , we can write $⟨ 𝑅 ⟩ = ⟨ 𝑍 ⟩ ⟨ 𝑦 ⟩$ , where the second factor is the expectation of $𝑦 𝑖$ for a randomly chosen extant species. According to above considerations, $⟨ 𝑅 ⟩$ changes over a small time interval $𝛿 𝑇 \geq 0$ to

⟨ 𝑅 (𝑇 + 𝛿 𝑇) ⟩ = ⟨ \sum 𝑗 \in 𝒜 𝑦 𝑗 - 𝛿 𝑇 𝜌 (𝑦 𝑗 - ‾ 𝑦) + 𝛿 𝑇 𝐶 𝑦 i n v ⟩ = ⟨ 𝑍 ⟩ ⟨ 𝑦 ⟩ + [- 𝜌 ⟨ 𝑍 ⟩ (⟨ 𝑦 ⟩ - ‾ 𝑦) + 𝐶 𝑦 i n v] 𝛿 𝑇,

(9)

where $𝑦 i n v$ denotes the mean associated with $𝑃 i n v (𝑦)$ . The sum above is over the extant species at time $𝑇$ . Extinctions of species can be disregarded at lowest order, because shortly before extinction their contribution to $𝑅$ is small. In the steady state the time-dependent term in Eq. vanishes:

- 𝜌 (⟨ 𝑦 ⟩ - ‾ 𝑦) ⟨ 𝑍 ⟩ + 𝐶 𝑦 i n v = 0 .

(10)

Making use of Eq. , we similarly derive in the Appendix the second moment:

⟨ 𝑅 (𝑇) 𝑅 (𝑇 + 𝛿 𝑇) ⟩ = ⟨ 𝑍 ⟩ 2 ⟨ 𝑦 ⟩ 2 + ⟨ 𝑍 ⟩ ⟨ 𝑦 2 ⟩ + ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ ‾ 𝑦 - ⟨ 𝑦 2 ⟩) 𝛿 𝑇 .

(11)

Combining the moment equations, we first evaluate

v a r 𝑅 = ⟨ 𝑅 2 ⟩ - ⟨ 𝑅 ⟩ 2 = ⟨ 𝑍 ⟩ ⟨ 𝑦 2 ⟩ .

(12)

Then we calculate the short-term autocorrelation function

c o r [𝑅 (𝑇), 𝑅 (𝑇 + 𝛿 𝑇)] = c o v [𝑅 (𝑇), 𝑅 (𝑇 + 𝛿 𝑇)] v a r 𝑅

(13)

and from this, considering that for an Ornstein Uhlenbeck process $c o r [𝑅 (𝑇), 𝑅 (𝑇 + 𝛿 𝑇)] = e x p (- ˜ 𝜌 | 𝛿 𝑇 |)$ , the relaxation rate of $𝑅$ as

˜ 𝜌 = - l i m 𝛿 𝑇 \to 0 + d c o r [𝑅 (𝑇), 𝑅 (𝑇 + 𝛿 𝑇)] d 𝛿 𝑇 = - ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ ‾ 𝑦 - ⟨ 𝑦 2 ⟩) ⟨ 𝑍 ⟩ ⟨ 𝑦 2 ⟩ = - 𝜌 (⟨ 𝑦 ⟩ ‾ 𝑦 - ⟨ 𝑦 2 ⟩) ⟨ 𝑦 2 ⟩ .

(14)

Equating $˜ 𝜌$ and $𝜌$ in Eq. yields our self-consistency condition. It simplifies to

‾ 𝑦 = 0, and so 𝑝 i n v = Φ (‾ 𝑦) = 12,

(15)

our main result. Our model ecosystems self-organise to exclude half of potential invaders.

Using the fact that dependence on $⟨ 𝑍 ⟩$ and so $𝐶$ drops out of Eq. , we can now address the case of general distributions of the off-diagonal coefficients $𝐴 𝑖 𝑗$ . Assume first that $𝐴 𝑖 𝑗$ can attain at random not one but two values other than zero. In this case the dynamics of the sum in Eq. through $𝑇$ is given by a linear combination of the dynamics of two sums of the form of Eq. . Since both have the same relaxation rate, Eq. , so has Eq. . One easily extends this argument to arbitrary discrete distributions for $𝐴 𝑖 𝑗$ and then to the limit of continuous distributions.

Following , Sec. 14.6, we compute further properties of the system steady state. Applying the classical mean first passage time formula to the process, Eq. , taking into account that the harvesting resistance of invaders is distributed as in Eq. , the mean time between invasion and extinction of a species is $𝜌 - 1 l n (2)$ . Since species invade the community at a rate one by construction, the steady-state mean number of species in the community is $⟨ 𝑆 ⟩ = 𝜌 - 1 l n (2)$ , so that

𝜌 = l n (2) ⟨ 𝑆 ⟩ \approx 0.693 ⟨ 𝑆 ⟩ .

(16)

A previous attempt at deriving $𝑝 i n v$ , Sec. 14.6 had incorrectly assumed $𝜌 = ⟨ 𝑆 ⟩ - 1$ on heuristic grounds.

To compute the distribution of $𝑦$ in a steady-state community, we construct the steady-state Fokker-Planck equation corresponding to Eq. for $‾ 𝑦 = 0$ :

0 = 𝜌 𝜕 [𝑦 𝑃 (𝑦)] 𝜕 𝑦 + 𝜌 𝜕 2 𝑃 (𝑦) 𝜕 𝑦 2 + \sqrt 2 𝜋 e x p (- 𝑦 2 2) .

(17)

The addition of $𝑃 i n v (𝑦)$ in the last term of Eq. describes invasion at rate one of new species with random harvesting resistances. These are balanced by extinctions, described the absorbing boundary conditions:

𝑃 (0) = 0, 𝜕 𝑃 (𝑦) 𝜕 𝑦 | 𝑦 = 0 = 𝜌 - 1 .

(18)

The solution to this equation is the distribution of $𝑦$ amongst extant species in the community, subject to Lotka-Volterra competition and a continuous turnover of interacting species through extinction and invasions. The distribution can be expressed using a hypergeometric function and is shown in Fig. . It has the moments $⟨ 𝑦 ⟩ = 2 1 / 2 𝜋 - 1 / 2 l n (2) - 1 = 1.1511$ and $⟨ 𝑦 2 ⟩ = 1 + l n (4) - 1 = 1.7213$ , which can be used to evaluate other characteristics of the steady state. From Eq. (17.25) of Ref. , Sec. 14.6, for example, we obtain the approximate mean steady-state species richness

⟨ 𝑆 ⟩ = 𝑆 ESI ⟨ 𝑦 2 ⟩,

(19)

and from Eq. (17.26) the approximate standard deviation of the fitness of new random invaders $𝑘$ ,

SD 𝑟 = (v a r 𝑟 𝑘) 1 / 2 = ⟨ 𝑦 2 ⟩ ⟨ 𝑦 ⟩ \cdot 𝜎 2 𝜇 (1 - 𝜇) 𝑠 .

(20)

FIG. 1.

The parameter-free distribution of standardized harvesting resistance in the steady state of our model.

Numerical simulations of the assembly model (detailed in the Appendix) reveal the main effects limiting the accuracy of our calculations. First, $𝑟 𝑖$ as given by Eq. cannot exceed $𝑠 𝑖 = 𝑠$ , which constrains the degree to which the $𝑟 𝑖$ and $𝑦$ can be described by an OU process or Eqs. and . This becomes a problem when $SD 𝑟$ , given by Eq. , is of the magnitude of $𝑠$ . The smaller $SD 𝑟$ is compared to $𝑠$ ; the better simulations agree with the analytic predictions (Fig. ). Second, the actual proportion of nonzero off-diagonal entries in ${𝐴 𝑖 𝑗} 𝑖, 𝑗 \in ℰ$ is smaller than $𝐶$ , possibly a reflection of the constraint that ${𝐴 𝑖 𝑗} 𝑖, 𝑗 \in ℰ$ must permit coexistence . This effect disappears when $SD 𝑟 ≪ 𝑠$ [Fig. ]. We therefore extrapolated all our numerical results to $SD 𝑟 \to 0$ [Figs. –]. In this limit, we verified our main result, Eq. , for invasion probability to within $1.8 %$ relative error; $𝑝 i n v$ is particularly robust to parameter changes in simulations [Fig. ]. We further verified Eq. for the relaxation rate of invasion fitness to within $0.6 %$ [Fig. ], Eq. for expected species richness $⟨ 𝑆 ⟩$ to within $1.0 %$ [Fig. ], and Eq. for $SD 𝑟$ itself to within $1.9 %$ [Fig. ].

FIG. 2.

Numerical test of theory. Panels compare five predictions and assumptions of our theory with steady-state outcomes of 105 simulations. Simulation parameters were chosen as $𝐶 = 10 - 2, 10 - 2.5, 10 - 3$ , and $𝐼$ such that $𝑆 ESI$ as given by Eq. equals 300, 350, 400, ..., 2000 [N.B.: $𝜇 = 𝐼 𝐶, 𝜎 2 = 𝐼 2 𝐶 (1 - 𝐶)$ ; see Appendix for details]. The abscissa in each graph represents the standard deviation $SD 𝑟$ of the fitness $𝑟$ of random invaders, as observed in simulations. Solid lines and 95% confidence ranges (gray) correspond to the quadratic regressions used to extrapolate to $SD 𝑟 = 0$ in the large-system limit. Horizontal dashed lines represent our analytic predictions. Panel (a) compares the sampled proportion of nonzero off-diagonal entries in ${𝐴 𝑖 𝑗} 𝑖, 𝑗 \in ℰ$ (“connectance”) with the assumed value $𝐶$ . (b) tests our prediction for invasion probability, Eq. , (c) our prediction for the relaxation rate $𝜌$ , Eq. . Panel (d) compares observed species richness with that predicted by Eq. , with $𝑆 ESI$ given by model parameters. (e) compares the observed standard deviation of invasion fitness, computed as in Eq. , with that predicted by Eq. (using the realized connectance for $𝐶$ rather than the model parameter).

III. DISCUSSION AND CONCLUSIONS

Structure and dynamics of ecosystems result from processes far from thermodynamic equilibrium, and fundamental physical laws help little in illuminating laws of ecology. Yet, just as quantum physics and relativity explain why some particles have spin $1 / 2$ , basic ecological principles imply that $𝑝 i n v = 1 / 2$ . Potentially, this is not the only constant of nature emerging from ecological principles. Other candidates are the diet partitioning exponent of predators , which empirically appears narrowly constrained around 0.54(2), or the ratio by which species richness declines from one level to the next, observed to be close to $1 / 3$ in both models and data .

Except for substantial work seeking universal exponents linking individual body masses to physiological and ecological rates , researchers do not generally consider that universal ecological constants might exist. Most ecological studies are designed to identify dependencies between quantities rather than testing constancy under well defined conditions. Nevertheless, there is striking empirical evidence consistent with our prediction that $50 %$ of attempted invasions succeed.

Invasion ecology distinguishes two processes: “establishment,” the successful invasion of a local ecological community by an introduced species, and “spread,” the subsequent successful invasion into a metacommunity. A metaanalysis of the invasion success of vertebrates transported between Europe and North America found that the mean establishment success for species introduced from Europe to North America was $59.6 \pm 11.6 %$ (standard error) and $52.4 \pm 11.9 %$ in the other direction. For spread, the corresponding figures were $65.0 \pm 16.5 %$ and $54.4 \pm 17.4 %$ , respectively. A more detailed metaanalysis confirmed that for birds invasion probability is $50 \pm 2.6 %$ , but reports $79 \pm 1.7 %$ for mammals . These probabilities are largely unchanged over a factor $106$ in invaded land mass area . One study estimating the invasion probability for 15 plant species into $0.5 \times 2 m$ plots in an old field reports probabilities around $50 %$ for most species (for the equilibrium richness of $\approx 13$ species per plot, we extract a median of $51 %$ and an average of $56 \pm 5 %$ ).

Considering the higher value for mammals, we note that our theory assumes that $𝜏 = c o r (𝐴 𝑖 𝑗, 𝐴 𝑗 𝑖)$ is zero, while interspecific resource competition can lead to $𝜏 > 0$ and predator-prey interactions to $𝜏 < 0$ . It is known that $𝜏 > 0$ decreases $𝑝 i n v$ (, Sec. 17.3). Predator-prey interactions between mammals might have the opposite effect, increasing $𝑝 i n v$ .

Other studies, especially with plants, report lower invasion probabilities. Only $10 %$ of $27009$ plant species introduced into Australia have formed self-sustaining wild populations . However, this figure is dominated by $25360$ ornamental gardening plants , for many of which we can assume that the climate and soil of Australia are unsuitable, violating our assumption that all $𝑠 𝑖$ in Eq. are equal. Our theory must break down when variability in $𝑠 𝑖$ (called environmental filtering ) is comparable or larger in magnitude than variability in the sum in Eq. , given by Eq. (called biotic filtering ). Including environmental filtering into the theory is an important next step.

Regarding the solutions of Eq. , here we have obtained a comprehensive description of the situation when $𝑆 pool ≫ 𝑆 ESI$ , deep inside the MA phase. The equilibrium species richness $⟨ 𝑆 ⟩$ predicted by Eq. is by a factor $⟨ 𝑦 2 ⟩$ smaller than the richness $𝑆 ESI$ at the phase transition. This is because steady states in species turnover arise already when the system is sufficiently sensitive that the perturbation induced by a new invader leads, on average, to extinction of a rare species, while $𝑆 = 𝑆 ESI$ corresponds to infinite sensitivity [, Eq. (17.24)]. The precise nature of the transition from $𝑆 pool = 2 𝑆 ESI$ to $𝑆 pool ≫ 2 𝑆 ESI$ remains unclear. Simulations indicate the existence of an intermediate parameter range where the community dynamically switches between a few recurring configurations (Clementsian temporal turnover, ). Whether this constitutes a distinct phase in the macroecological limit $𝑆 ESI \to \infty$ remains unclear.

ACKNOWLEDGMENTS

This work forms part of the project “Mechanisms and prediction of large-scale ecological responses to environmental change” funded by the Natural Environment Research Council (NE/T003510/1).

Appendices

APPENDIX

1. Derivation of Eq.

We derive Eq. , making first use of Eq. and then of the fact that $v a r 𝑍 = ⟨ 𝑍 2 ⟩ - ⟨ 𝑍 ⟩ 2 = ⟨ 𝑍 ⟩$ on account of the Poisson distribution of $𝑍$ :

⟨ 𝑅 (𝑇) 𝑅 (𝑇 + 𝛿 𝑇) ⟩ = ⟨ \sum 𝑗, 𝑙 \in 𝒜 𝑦 𝑗 (𝑇) 𝑦 𝑙 (𝑇 + 𝛿 𝑇) + \sum 𝑗 \in 𝒜 𝑦 𝑗 𝐶 𝑦 i n v 𝛿 𝑇 ⟩ = ⟨ \sum 𝑗, 𝑙 \in 𝒜 𝑦 𝑗 [𝑦 𝑙 - 𝜌 (𝑦 𝑙 - ‾ 𝑦) 𝛿 𝑇] ⟩ + ⟨ 𝑦 ⟩ ⟨ 𝑍 ⟩ 𝐶 𝑦 i n v 𝛿 𝑇 = ⟨ \sum 𝑗, 𝑙 \in 𝒜 𝑗 \neq 𝑙 𝑦 𝑗 [𝑦 𝑙 - 𝜌 (𝑦 𝑙 - ‾ 𝑦) 𝛿 𝑇] ⟩ + ⟨ \sum 𝑗 \in 𝒜 𝑦 𝑗 [𝑦 𝑗 - 𝜌 (𝑦 𝑗 - ‾ 𝑦) 𝛿 𝑇] ⟩ + ⟨ 𝑦 ⟩ ⟨ 𝑍 ⟩ 𝐶 𝑦 i n v 𝛿 𝑇 = ⟨ 𝑍 (𝑍 - 1) ⟩ [⟨ 𝑦 ⟩ 2 - 𝜌 (⟨ 𝑦 ⟩ 2 - ⟨ 𝑦 ⟩ ‾ 𝑦) 𝛿 𝑇] + ⟨ 𝑍 ⟩ [⟨ 𝑦 2 ⟩ - 𝜌 (⟨ 𝑦 2 ⟩ - ⟨ 𝑦 ⟩ ‾ 𝑦) 𝛿 𝑇] + ⟨ 𝑦 ⟩ ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ - ‾ 𝑦) ⟨ 𝑍 ⟩ 𝛿 𝑇 = ⟨ 𝑍 ⟩ 2 ⟨ 𝑦 ⟩ 2 + ⟨ 𝑍 ⟩ ⟨ 𝑦 2 ⟩ + (⟨ 𝑍 ⟩ 2 - ⟨ 𝑍 2 ⟩) 𝜌 (⟨ 𝑦 ⟩ 2 - ⟨ 𝑦 ⟩ ‾ 𝑦) 𝛿 𝑇 + ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ 2 - ⟨ 𝑦 2 ⟩) 𝛿 𝑇 = ⟨ 𝑍 ⟩ 2 ⟨ 𝑦 ⟩ 2 + ⟨ 𝑍 ⟩ ⟨ 𝑦 2 ⟩ - ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ 2 - ⟨ 𝑦 ⟩ ‾ 𝑦) 𝛿 𝑇 + ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ 2 - ⟨ 𝑦 2 ⟩) 𝛿 𝑇 = ⟨ 𝑍 ⟩ 2 ⟨ 𝑦 ⟩ 2 + ⟨ 𝑍 ⟩ ⟨ 𝑦 2 ⟩ + ⟨ 𝑍 ⟩ 𝜌 (⟨ 𝑦 ⟩ ‾ 𝑦 - ⟨ 𝑦 2 ⟩) 𝛿 𝑇 .

2. Computation of harvesting resistance

Without loss of generality, we show that Eq. evaluates to harvesting resistance $𝑟 1 = (d l n (𝑏 1) / d 𝑠 1) - 1$ if species $𝑖 = 1$ is extant, assuming the other $𝑆 - 1$ extant species have indices $2, . . ., 𝑆$ .

Let $𝐀$ and $𝐬$ be, respectively, the interaction matrix and the vector of linear growth rates $𝑠 𝑖$ entering Eq. , but restricted to the $𝑆$ extant species. Define $𝐌 = 𝐀 - 1$ . We break both $𝐀$ and $𝐌$ into four blocks to separate the entries for the species with index 1 from all other species, assigned a mnemonic index “o”

𝐀 = (𝐴 11 𝐀 1 o 𝐀 o 1 𝐀 o o),

(A1)

and, using a well-known block-wise inversion formula ,

𝐌 = (𝑀 11 𝐌 1 o 𝐌 o 1 𝐌 o o)

(A2)

= (𝑀 11 - 𝑀 11 𝐀 1 o 𝐀 - 1 o o - 𝐀 - 1 o o 𝐀 o 1 𝑀 11 𝐀 - 1 o o + 𝐀 - 1 o o 𝐀 o 1 𝑀 11 𝐀 1 o 𝐀 - 1 o o),

(A3)

where $𝑀 11 = (𝐴 11 - 𝐀 1 o 𝐀 - 1 o o 𝐀 o 1) - 1$ . Similarly, we break up

𝐬 = (𝑠 1 𝐬 o) .

(A4)

By Eq. of the main text, $𝐛 = 𝐀 - 1 𝐬 = 𝐌 𝐬$ is the vector of the equilibrium biomasses $𝑏 𝑖$ of the extant species. In particular, $𝑏 1 = 𝑀 11 𝑠 1 + 𝐌 1 o 𝐬 o$ . Hence, by the definition of harvesting resistance,

𝑟 1 = (d l n 𝑏 1 d 𝑠 1) - 1 = 𝑏 1 d 𝑏 1 / d 𝑠 1 = 𝑀 11 𝑠 1 + 𝐌 1 o 𝑠 o 𝑀 11 .

(A5)

By Eq. , $𝐌 1 o = - 𝑀 11 𝐀 1 o 𝐀 - 1 o o$ . Inserting this above, $𝑀 11$ cancels out and we obtain

𝑟 1 = 𝑠 1 - 𝐀 1 o 𝐀 - 1 o o 𝑠 o .

(A6)

Defining $𝐛' o = 𝐀 - 1 o o 𝐬 o$ , this becomes

𝑟 1 = 𝑠 1 - 𝐀 1 o 𝐛' o,

(A7)

which is Eq. of the main text for $𝑖 = 1$ , with $𝐛' o$ defined as prescribed below this equation.

Hence, we have shown that Eq. of the main text provides the value of harvesting resistance in line with its definition. The fact that some entries of the abbreviation $𝐛' o = 𝐀 - 1 o o 𝐬 o$ entering this equation may be negative does not invalidate this derivation and does not imply that negative population biomasses actually arise in our calculation.

3. General simulation method

With a direct simulation of the model, computation time increases with system size approximately as $𝑆 4 ESI$ : by a factor $𝑆 2 ESI$ because of the vector-matrix operation in Eq. , by a factor $𝑆 ESI$ because the time to reach equilibrium grows approximately as $𝑆 ESI$ , and by another factor $𝑆 ESI$ because turning over an existing community of size $𝑂 (𝑆 ESI),$ so to generate statistically independent data, requires $𝑂 (𝑆 ESI)$ invasions.

To allow evaluation of the model for values of $𝑆 ESI$ up to 2000, we therefore computed the new equilibrium of Eq. after invasion of a species by directly solving the linear system of equilibrium equations,

0 = 𝑠 𝑖 - \sum 𝑗 \in ℰ + 𝐴 𝑖 𝑗 𝑏 𝑗 (𝑖 \in ℰ +),

(A8)

rather than through numerical integration of the model. When one or more species in the solution of Eq. had negative population biomass $𝑏 𝑗$ , we considered the one with lowest $𝑏 𝑗$ extinct, removed it from $ℰ$ , and solved Eq. once more, repeating until all populations biomasses were nonnegative. In the simulation, we counted as iterations of the model both cases where a randomly sampled species could invade and where it could not.

The above procedure is plausible when only a single species attains negative biomass after an invasion. When there are several candidates for extinction, one might be concerned that the order of species removal affects the outcome. We note, however, that (i) in the model steady state only a single species goes extinct after each successful invasion on average and fluctuations around this mean are small; (ii) indirect interactions between species are of similar magnitude as direct interactions , implying that extinct species are not usually directly interacting with the added species or with each other; (iii) when a species $𝑖$ goes extinct, $| 𝑏 𝑖 |$ is generally small both before and after the species addition step (especially for large $𝑆 ESI$ ), implying that the effect of its removal on other species in $ℰ +$ remains small. Because of this, concerns about the order of species removal have little bearing in practice.

To achieve the desired efficiency gain from directly solving Eq. , one must keep in mind that the direct numerical solution of this equation requires time of order $𝒪 (𝑆 3 ESI)$ , implying that this would not fundamentally improve the algorithm compared to direct simulation. Instead, we therefore kept track of both the matrix ${𝐴 𝑖 𝑗} 𝑖, 𝑗 \in ℰ$ and its inverse between evaluations of Eq. , using again the block-wise inversion formula to account for addition or removal of a single species in $ℰ$ .

Specifically, if $𝐀'$ is the matrix obtained from $𝐀$ by the addition of row $𝐯$ and column $𝐮$ , i.e.,

𝐀' = (𝐀 𝐮 𝐯 T 1),

(A9)

then

𝐀' - 1 = (𝐀 - 1 + 𝑐 - 1 (𝐀 - 1 𝐮) (𝐯 T 𝐀 - 1) - 𝑐 - 1 𝐀 - 1 𝐮 - 𝑐 - 1 𝐯 T 𝐀 - 1 𝑐 - 1)

(A10)

with $𝑐 = 1 - 𝐯 T 𝐀 - 1 𝐮$ , which allowed us to compute $(𝐀') - 1$ from $𝐀 - 1$ and vice versa with $𝒪 (𝑆 2 ESI)$ operations. To ensure against loss of accuracy, we computed the inverse of ${𝐴 𝑖 𝑗} 𝑖, 𝑗 \in ℰ$ directly every 1000 iterations and when $| 𝑐 | < 10 - 7$ in a species removal step.

4. Simulation parameters

We evaluated the above assembly model for communities with values of connectance $𝐶 = 10 - 2, 10 - 2.5, 10 - 3$ and interaction strength

𝐼 = [𝑆 ESI 𝐶 (1 - 𝐶)] - 1 / 2

(A11)

for $𝑆 ESI = 300$ , 350, 400, ..., 2000, running the model over $200 \times 𝑆 ESI$ iterations, starting with a single species.

5. Ornstein-Uhlenbeck test and determination of relaxation rate $𝜌$

Of the quantities displayed in Fig. of the main text, only the relaxation rate $𝜌$ is difficult to extract from simulations. To obtain $𝜌$ , we simulated the invasion fitness of 200 test species (which never actually invade) throughout a simulation.

To avoid artifacts due to small numbers of interaction partners or due to the fact that invasion fitness cannot exceed $𝑠 𝑖$ ( $= 1$ ) in our model, which might arise for low $𝐶$ , we artificially fixed the connectance of test species (i.e., the probability for these species to be suppressed by a given resident species) at $40 / 𝑆 ESI$ .

For each test species, we sampled the sum of the biomasses $𝑏 𝑖$ of the suppressing resident species (hereafter called “suppression”) after every $𝑆 ESI / 10$ iterations and discarded the first 100 samples as burn in. To circumvent the problem of determining the expected mean suppression, we paired the first 100 test species up with the other 100 test species and computed time series for the difference of suppression of each pair, as this difference has mean zero and the same autocorrelation function as both components.

To these 100 time series, we then applied a cosine-bell taper and computed their periodograms. We then averaged the 100 periodograms to estimate the power spectrum. To the lowest 200 nonzero frequencies $𝑓$ of this power spectrum we obtained the least-square fit to a Cauchy profile (the Fourier transform of the exponential autocovariance function of an Ornstein-Uhlenbeck process),

𝑎 (𝜌') 2 + (2 𝜋 𝑓) 2,

(A12)

with fitting parameters $𝑎, 𝜌' > 0$ .

As a test for our assumption that invasion fitness (and so suppression) follows an Ornstein-Uhlenbeck process, we verified that the Cauchy profile, Eq. , described the observed power spectra well (see example in Fig. ). Remaining small deviations from the Cauchy profile are likely attributable to numerical limitations, including the finite duration of the simulated time series.

FIG. 3.

Test for consistency of invasion fitness (or suppression) with Ornstein-Uhlenbeck process. Top: average of 100 periodograms (circles) compared with a fitted Cauchy profile for a community with $𝐶 = 10 - 2$ and $𝑆 ESI = 2000$ , illustrating the good fit. Bottom: ratio of averaged periodogram to Cauchy profile. In principle, the ratio should be one at each frequency, with residual variance close to $1 / 100 = 0.01$ (since periodograms have a coefficient of variation of one). We find that there is only a small increase with frequency—over the frequency range in which the Cauchy profile declines by a factor 180, the ratio increases by just $7 %$ —and residual variance is 0.0097, demonstrating the good fit.

From the fitted $𝜌'$ we computed $𝜌$ by taking account of the subsampling of the time series and, dividing by the observed value of $𝑝 i n v$ , of the fact that $𝑇$ in the main text counts successful invasions only.

Self-organization of ecosystems to exclude half of all potential invaders

Abstract

Physics Subject Headings (PhySH)

Article Text

I. INTRODUCTION

II. METHODS AND RESULTS

III. DISCUSSION AND CONCLUSIONS

ACKNOWLEDGMENTS

Appendices

APPENDIX

1. Derivation of Eq.

2. Computation of harvesting resistance

3. General simulation method

4. Simulation parameters

5. Ornstein-Uhlenbeck test and determination of relaxation rate $𝜌$

References (34)

Self-organization of ecosystems to exclude half of all potential invaders

Abstract

Physics Subject Headings (PhySH)

Article Text

I. INTRODUCTION

II. METHODS AND RESULTS

III. DISCUSSION AND CONCLUSIONS

ACKNOWLEDGMENTS

Appendices

APPENDIX

1. Derivation of Eq. (11)

2. Computation of harvesting resistance

3. General simulation method

4. Simulation parameters

5. Ornstein-Uhlenbeck test and determination of relaxation rate 𝜌

References (34)

1. Derivation of Eq.

5. Ornstein-Uhlenbeck test and determination of relaxation rate $𝜌$