Evolution of prudent predation in complex food webs

Abstract Prudent predators catch sufficient prey to sustain their populations but not as much as to undermine their populations’ survival. The idea that predators evolve to be prudent has been dismissed in the 1970s, but the arguments invoked then are untenable in the light of modern evolution theory. The evolution of prudent predation has repeatedly been demonstrated in two‐species predator–prey metacommunity models. However, the vigorous population fluctuations that these models predict are not widely observed. Here we show that in complex model food webs prudent predation evolves as a result of consumer‐mediated (‘apparent’) competitive exclusion of resources, which disadvantages aggressive consumers and does not generate such fluctuations. We make testable predictions for empirical signatures of this mechanism and its outcomes. Then we discuss how these predictions are borne out across freshwater, marine and terrestrial ecosystems. Demonstrating explanatory power of evolved prudent predation well beyond the question of predator–prey coexistence, the predicted signatures explain unexpected declines of invasive alien species, the shape of stock–recruitment relations of fish, and the clearance rates of pelagic consumers across the latitudinal gradient and 15 orders of magnitude in body mass. Specific research to further test this theory is proposed.

We explain here in detail why Eqs. (3) and (4) are plausible approximations for sampling the interaction strengths of new species entering our model community.
Under quite general conditions it is possible to approximate the dependence of attack rates on the traits of consumers and resources in the form (Rossberg et al., 2010;Nagelkerke & Rossberg, 2014;Rossberg, 2013, Ch. 8): with D denoting the dimensionality of trophic niche space and v D vulnerabilityand foraging traits of resources and consumers, respectively. These trophic traits can be computed as functions of observable biological traits (Nagelkerke & Rossberg, 2014). A similar representation has been proposed by Rohr et al. (2010). The constant a 0 has dimensions of attack rates and σ k = ±1. There is some ambiguity in how to choose a 0 , σ k and the functions mapping observed traits to trophic traits. However, when imposing a condition that the mean of (v  (Rossberg, 2013, Ch. 8) and the choice of a 0 . With the mean of (v (j) 0 ) 2 thus minimised, we shall approximate v (j) 0 = 0. For large D and sufficient statistical independence of the components of v (j) and f (k) (Rossberg, 2013, Ch. 11), one can approximate the sum in Eq. (S1) for randomly sampled consumer-resource pairs (j, k) by a normal distribution. Denoting the mean of this normal distribution by µ, its variance by σ 2 , and defining a k = a 0 exp(f D is known to be slow (Rossberg et al., 2006;Bersier & Kehrli, 2008;Eklöf & Stouffer, 2016)-a median of 25 times slower in an analysis of Rossberg et al. (2006). It shall here be disregarded.
Mutations of any observable biological traits will affect several foraging traits f (k) 0 , . . . , f D . The question whether this increases of decreases short-term fitness (Goodnight et al., 2008) in a given community depends not only on all traits f (k) 0 , . . . , f (k) D of the focal consumer k but also on the sets of resources and competitors in the community. Even when a mutation leads to an increase in short-term fitness, the change in f (k) 0 associated with this mutation might be positive or negative, because, provided niche space dimensionality D is not too low, the associated change in f (k) 0 is just one of many random changes in foraging traits resulting from this mutation. Mutants arriving at the focal patch from a source patch may therefore have f (k) 0 values that can be higher or lower than the f (k) 0 of the propagule that founded the population in the source patch. Because, all else equal, smaller f (k) 0 correspond to consumers that, overall, forage less efficiently than consumers with larger f (k) 0 , and low efficiency is mechanically easier to achieve than high efficiency, one must plausibly assume that degeneration of traits through mutations (Pomiankowski et al., 1991) leads to a decay of f (k) 0 on average unless this is counteracted by selection pressure. Recalling that a k = a 0 exp(f (k) 0 − µ), this leads to Eq. (3). We assume that the relevant species pools are large and diverse, such that different patches have in effect statistically independent, typically non-overlapping species compositions. The random variables ξ jk in Eq. (4) are therefore sampled anew as a propagule arrives at the focal patch, independent of a consumer's interactions with the residents of its source patch. Only the inheritance of a k must be accounted of.
As a caveat, we note that in reality vulnerability traits do not cover the D-dimensional trophic traits space evenly, e.g. because these traits carry phylogenetic signal (related species have similar consumers, Bersier & Kehrli 2008). Then foraging traits other than f (k) 0 might contribute to long-term fitness as well. For simplicity, we disregard this complication in our model.

Appendix S2 Derivation of the sub-models of the deconstructed formulation
We provide the rationale and outline the derivation of the four criteria Eqs. (13-17) driving invasions and extirpations in the deconstructed model formulation.
The invasibility criterion, Eq. (13), predicts invasibility when disregarding the presence of all but the focal consumer in the dimensionless full model, Eq. (2). Formally, it is obtained by computing the equilibrium state of Eq. (2) for S C = 1 and B C k = 0 (with k = 1), which is B R j = 1 for 1 ≤ j ≤ S R , and then extracting the condition that, by Eq. (2b), this equilibrium is unstable such that the consumer can invade: The condition for the overexploitation of resource j during the expansion phase of an invading consumer k, Eq. (15), is obtained by analysing the dimensionless full model, Eq. (2), for the case of only one consumer and one resource: S C = 1, S R = 1 (with j = k = 1). We consider again the situation where the consumer is initially absent B C k = 0 and the resource at equilibrium B R j = 1, dB R j /dt = 0. Then the consumer invades at low abundance. To estimate the minimum of B R j attained during the consumer invasion, i.e. during the transient before a new equilibrium is reached, we approximate dynamics by disregarding the density dependence of resource production expressed by the term −B R j in Eq. (2a). This approximation is justified because we are interested in situations where B R j falls below M min 1. It reduces the model to the classical Lotka-Volterra predator-prey equations Evaluating the conservation law known for this system (Lotka, 1920) for the initial conditions (Rossberg, 2013, Sec. 20.3.3). Since we are interested in situations where the minimum is deep (B R j < M min 1), this condition can be approximated as ln(B R j ) = −H jk . It follows that B R j falls below M min during consumer k's invasion if ln(M min ) > −H jk , which is equivalent to Eq. (15).
The conditions for consumer-mediated competitive exclusion, for exploitative competitive exclusion and for Pyrrhic competition all derive directly from exact equilibrium solutions of the dynamic model. The general multispecies model, Eq. (2), is well studied (MacArthur, 1970(MacArthur, , 1972Case & Casten, 1979;Chesson, 1990). To write down its equilibrium solution, let H be the matrix with entries H jk and define the competition matrix as the matrix with entries Denote by s the vector of intrinsic consumer growth rates with R k = SR j=1 H jk defined as in the main text. The vector b C of consumer population biomasses B C j at equilibrium is then given by and that of resource population biomasses B R j by In the case of only one consumer (S C = 1, k = 1), the biomass of the resource j is therefore B R j = 1 − H jk (C kk ) −1 s k . The resource with the lowest biomass is that with the largest H jk , i.e., the main resource of k. Its biomass is negative, implying resource extinction (Holt, 1977), if The criterion for consumer-mediated competitive exclusion, Eq. (16), spells out this condition. For the two-consumer (S c = 2) problem, we have, with k = 1 and l = 2, Combining Eqs. (S5) and (S8), we find that (for or equivalently Our criterion of exploitative competitive exclusion, Eq. (14) spells out this condition. Now, assume that Eq. (S10) and the corresponding condition with l's and k's role reversed both fail to be satisfied. This alone does not guarantee coexistence of all species. Combining Eqs. (S5), (S6) and (S8), one can see that the equilibrium abundance of resource B R i is predicted to be negative if This can be re-arranged to and our condition for Pyrrhic competition, Eq. (17), spells out this inequality. We now outline how these conditions can efficiently be evaluated for large S R and S C . The most time-consuming step is the computation of C in Eq. (S3), as (for practical purposes) the number of operations this requires increases as O(S 2 C S R ) with system size. All remaining calculations can be done using just O(S 2 C ) or O(S C S R ) operations. Denote, for any square matrix A, by diag(A) the vector formed by its diagonal elements, and by Diag(v), for any vector v, the diagonal matrix with v on the diagonal. We can evaluate the S C × S C matrix Φ with entries Φ kl given by the left-hand side of Eq. (S9) as (S13) To test for extirpations, set the diagonal of Φ to exactly zero to remove small numerical errors. Extirpation of consumer k by our (simplified) criterion follows if row k of Φ constrains negative elements. The S C ×S C matrix D with entries D kl = C kk C ll −C 2 kl , containing the determinants of all two-consumer competition problems (the denominators in Eqs. (S8), (S11)), can be computed as with • denoting elementwise multiplication. After finding for each consumer k the index m(k) of its main resource, one can constructed the S C × S C matrix M with entries Using this, we obtain the S C × S C matrix ∆ with entries given by the difference between left-and right-hand side of Eq. (S12) for the main resource of each consumer k as To test for extirpations, set the diagonal of ∆ to exactly zero to remove small numerical errors. Extirpation of the main resource of consumer k by our (simplified) criterion follows if row k of ∆ contains negative elements. By striking a new balance between code complexity, speed, and accuracy in the multi-objective optimisation problem of finding fast, simple and accurate models, our deconstructed formulation carves out emergent properties (sensu Rossberg, 2007) of the full model, Eq. (2), e.g., those shown in Figs. 2 and 3.

Appendix S3 Evolutionary steady-state condition including mutation bias
We derive the steady state condition for base attack rate, Eq. (6).
To understand the effect of mutation bias, we invoke the Price equation (Price, 1972). It predicts that the expected rate of evolutionary change of a trait q of a species is given by with f (q) denoting the invasion fitness (for a given environment) of lineages of type q, and the last term representing the mutation bias (the mean inherent rate of change of traits). For trait values q * corresponding to evolutionary steady states, both sides of Eq. (S17) must evaluate to zero. Following Page & Nowak (2002), we expand f (q) to first order at q = q * . Combined with the population-dynamical equilibrium condition f (q * ) = 0, this leads to 0 = f (q * ) var q + Eq, or equivalently This condition generalises the conventional criterion for evolutionary singular strategies, f (q * ) = 0, to situations with mutation bias.
To apply this result to our model, we first note that our way of implementing evolution of base attack rate in our model through Eq. (3) makes the scheme formally analogous to how individual-based models describe evolution of a trait in a population. In this analogy, (i) coexisting populations of consumers in our model correspond to individuals in conventional models, (ii) mutations accumulated between successive invasion events along lineages in our model correspond to mutations between successive generations in conventional models, and (iii) the community of consumers at the focal patch in our model corresponds to the evolving population in conventional models.
While this formal analogy does not represent the process our community assembly model actually describes, it allows us to use Eq. (S18), with q = ln a, to obtain the steady-state condition for mean logarithmic base attack rate in this model. With the assumed statistical equivalence of all patches in the metacommunity (mean-field approximation), such a steady state can arise only if the evolution of the base attack rates of all consumers in the metacommunity has reached a steady state. The steady state condition for the community mean of ln a therefore implies a steady-state condition for the evolution of ln a k for all species k in the metacommunity (see also Appendix S6).
To evaluate Eq. (S18) with q = ln a, we set where L * is the mean lifetime of populations in the community. The standing mutational variance var q = var(ln a) is obtained from the distribution of a over the simulation steady state. We approximate steady-state invasion fitness, i.e. the mean intrinsic rate of increase (f (q) > 0) or decrease (f (q) < 0) of the number of populations of type q in the simulation steady state, as f (q) ≈ ln[R(a)]/L(a), where R(a) and L(a) are as defined in the main text. With a * representing the geometric mean of a over the simulation steady state, such that ln a * is the arithmetic mean of ln a, we expect that R(a * ) = 1. This leads to Putting Eqs. (S19) and (S20) into Eq. (S18) and multiplying both sides with L * gives Applying the identity ln(x) = log 10 (x)/ log 10 [exp(1)] yields Eq. (6).

Appendix S4 Mechanisms determining 'birth rate'
Within the mean-field approximation, the 'birth rate' function b(a) is in our model the rate at which resident populations with base attack rate a give rise to successful invaders into the model community via Eq. (3). To derive an analytic approximation of b(a), we must account for three model elements, which are common to both model formulations: 1. The mutation step, Eq.
(3), determining the new consumer's base attack rate from that of the resident.
2. The sampling of the new consumer's attack rates according to Eq. (4), and the test whether it can invade.
3. The fact that time is measured in numbers of successful consumer invasions.
Crucial is the probability of successful invasion in 2. We begin with an analysis of this element, adding subsequently considerations of 1 and 3. We first consider the deconstructed model formulation. On one hand, competitive exclusion by a resident consumer according to Eq. (14) of the algorithm always implies an inability to invade according to Eq. (13), so that (for S C > 0) only Eq. (14) needs to be considered. On the other hand, Eq. (14) can be read as just a correction of the invasibility criterion, Eq. (13) to account for the presence of competitors. To see this, re-arrange Eq. (14) as The term in square brackets represents the population biomass ( Because there is no mechanism active in the model that would favours values of the expression in square brackets above that are particularly close to zero (see also Fig. S1), most of the variation among the addends in the sum over j in Eq. (S22) is due to the log-normal distribution of the invader's attack rates H jk . The presence of competitors merely moderates the effect of this variation. It can therefore be approximated by substituting the square bracket by a suitable constant 0 < β < 1: the fitting parameter introduced in the main text.
The sum over j in Eq. (S22) can then be written as α 0 a k β SR j=1 e σξ jk . The distribution of the sum in this last expression is, for a given number of resources S R , often well approximated by a single log-normal distribution with suitable choices for mean µ SR ≈ σ √ 2 ln S R and standard deviation σ SR ≈ σ/ √ 1 + 2 ln S R of the logarithm (Rossberg et al., 2011). (We estimated µ SR and σ SR numerically form 10,000 samples of log-normal sums, which is more accurate.) From this log-normal approximation, the invasion probability for species with given base attack rate a k is obtained as with Φ(x) denoting the cumulative standard normal distribution function. For the full model, the same functional form as in Eq. (S23) can be chosen based on the same rationale: compared to the variation in link strengths, the variation among the population biomasses of resident resources is small. Denote by P * inv (a) the probability that the "offspring" of a resident species with base attack rate a can invade the community. The log-normal approximation for the sum in α 0 a k β SR j=1 e σξ jk used above combines seamlessly with the log-normal distribution of a k resulting from the mutation of base attack rate a l of the 'parent' population l as per Eq. (3). We can therefore obtain P * inv (a) from Eq. (S23) by correcting µ * SR = µ SR + ln γ 0 and σ * SR = [σ 2 SR + (ln γ 1 ) 2 ] 1/2 to account for mutational variance and bias. Hence Because we measure time in units of consumer invasions, and both variants of our model attempt consumer invasions until one succeeds, the probability for "offspring" of resident consumer l to invade in a given time step is P * inv (a l )/ SC k=1 P * inv (a k ) (guaranteeing that the probability for offspring of some consumer k to invade evaluates to 1). Since species richness and the distribution of the a l fluctuate somewhat through time, we calculated the 'birth rate' in Fig. 3c,g as the average of this probability for a given base attack rate a over the model steady state: We derive Eq. (10), which predicts a consumer's intrinsic growth rate after serial extinction in the limit of high base attack rate and species richness. Note first that, because during serial extinction resources are successively removed in decreasing order of the consumer's attack rate (and also in the simplified model, Box 2), the distribution of attack rates after serial extinction is the same as before, except for being truncated from above at the point where Eq. (16) gets violated. In situations where the sums in Eq. (16) are not dominated by just a few resources, the central limit theorem can be invoked and the sums approximated by their expectation values, which then permits analytic computation of the truncation threshold H * and other properties of the end state.
The calculations simplify by first approximating the relevant section of the upper tail of the log-normal attack-rate distribution, Eq. (4), by a Pareto distribution, which can be derived in the limit of high resource richness S R (Rossberg et al., 2011;Rossberg, 2013). By this approximation, the consumer has on average Z resources with H jk larger than some "observation threshold" H 0 , and for these with ν = σ −1 √ 2 ln S R . Empirically, typical values for ν are in the range 0.5 to 0.6 (Rossberg et al., 2011;Rossberg, 2013). Values ν ≥ 1 would correspond to extreme omnivory where the proportional contribution of each resource species to a consumer's diet scales as 1/S R , i.e. no resource makes a sizeable contribution to the diet. We are unaware of such a situation occurring in nature, and therefore assume 0 < ν < 1 in this study.
For a given observation threshold H 0 one can define the consumer link density Z, i.e. the mean over all consumers of the number of resources with scaled attack rate H jk above the threshold H 0 . By choosing Z, we can control the typical strengths H max of the strongest attack rate before serial extinction, specifically the exp(−1)-quantile of the distribution of max j H jk . In the limit of large Z, this leads to the condition It goes without saying that H max is proportional to base attack rate a k and can therefore be use as a proxy for the latter.
With this preparation, we can now take expectation values on both sides of Eq. (16) for the case of truncation of the link-strength distribution at H * , the largest value for which Eq. (16) is not violated. This leads to a condition which can be written as where p(x) = −(d/dx)P [H jk ≤ x] is the probability density of the untruncated attack rate distribution. Evaluation of the integrals after inserting Eq. (S26) leads to and, after inserting Eq. (S28) and taking the limit of low observation threshold (H 0 → 0),  Figure S2: Dependence of intrinsic growth rate term C = 1− j H jk on base attack rate in the course of repeated serial extinction and resource turnover. Panel (a) shows geometric means of C over 10 6 replicated runs of the model of Box 2 over 10 4 iterations, panel (b) arithmetic means. For high base attack rates a k (dark lines), both geometric and arithmetic means approach the same value ≈ 0.86 (indicating a near-deterministic outcome) after the first iteration of consumer-mediated competitive exclusion, largely independent of base attack rate, as predicted by the analytic theory. The value is different from the analytic prediction 1 − σ −1 √ 2 log S R ≈ 0.14 valid for large because S R , because S R = 224 is not sufficiently large.
This equation can be solved for H * , yielding The expected intrinsic growth rate of the consumer after serial extinction equals the left-hand sides of Eqs. (S29)-(S32). When putting Eq. (S33) into the left-hand side of Eq. (S32) it simplifies considerably, leading to the final result H jk therefore has contributions from many small terms, justifying our application of the central limit theorem to approximate the sums entering Eq. (16) by their expectation values. Figure S2b qualitatively confirms this result.
Interestingly, above considerations imply that, despite having the same niche width in terms of the spread σ of the log-normal attack-rate distribution, invaders with higher base attack rate will have more diverse diets post Impact than those with lower attack rates. This might explain why invasive alien consumers are often found to be 'generalists'.
Box S1 Algorithm of the evolutionary metapopulation model.
The model state is given by N patches which are either empty or occupied by a population with base attack rate ai (1 ≤ i ≤ N ). The model is simulated as follows: 1. Occupy a proportion p of patches with populations with identical initial base attack rates ai.
2. Select an occupied source patch l for dispersal. Sample the base attack rate a k of a propagule according to Eq. (3).
3. Sample a target patch j.
4. If patch j is occupied: (a) If aj < a k , replace the new population of patch j with one that has base attack rate a k , otherwise do nothing.
5. If patch j is not occupied: (a) With invasion probability P inv (a k ), establish in patch j a new population with base attack rates a k and then remove the population from another occupied patch m, sampled at random from all occupied patches with probability proportional to 1/L(am). P inv (a) is our approximation of invasion probability for the deconstructed community model, Eq. (S23) with β = 0.45 and S R = 224 (corresponding to the mean equilibrium richness in Fig. 2b), and L(a) the polynomial fit to mean population lifetime in Fig. 3h (log 10 L = −0.04105026(log 10 a) 2 −0.78404937 log 10 a−0.77341520).

Continue from
Step 2 for a predetermined number of iterations.
The values of γ0, γ1, and σ are as in Tab. 1.
The algorithm can be reformulated in such a way that only a list of the ai value of occupied patches i is kept in memory. In each iteration, Step 4a is then executed with probability p and otherwise Step 5a. When invasion is successful in Step 5a, the new a k value is stored in the memory location where am was previously stored. This formulation permits us to take the limit p → 0 while keeping the number of occupied patches pN fixed.

Appendix S6 The limited impact of cheaters
Cheaters exploit benefits offered by more altruistic conspecifics to their advantage, thus potentially counteracting the evolution of altruism. To obtain a bound on the impact of cheaters on prudent predation, we devised a simple evolutionary metapopulation model. For a single species, this model explicitly describes the kind of evolutionary process that the evolutionary community model, introduced in the main text, implicitly describes for many species using a mean-field approximation. By comparing the evolutionary steady states reached by the two models for the special case where cheating is absent, the plausibility of the mean-field approximation can be tested.
Specifically, the evolutionary metapopulation model describes a landscape of N patches that are either occupied by the focal species or not. The population occupying patch i (1 ≤ i ≤ N ) has an associated base attack rate a i . We assume that cheating occurs if a population of the focal species disperses to a patch that is already occupied, and the propagule's base attack rate is larger than that of the resident in that patch. The propagule then replaces the resident population. This model disregards that conspecific propagules will not only differ in their base attack rates from residents, but also in other foraging traits (Appendix Appendix S1), and therefore have, on average, a reduced likelihood of encountering suitable resources and establishing themselves. Our metapopulation model is therefore biased to overestimate the likelihood of cheating. We shall see that the predicted impact of cheating remains limited despite this.
Contrasting conventional stochastic patch occupancy models in the tradition of Levins (1969), patch occupancy p, i.e., the proportion of occupied patches, is a parameter in our model. The reason is evidence that species richness both at patch level (α) and at landscape level (γ) is regulated through ecological structural stability limits (O'Sullivan et al., 2019), which our metapopulation model cannot explicitly represent. Mean occupancy is uniquely determined by α and γ as p = α/γ. By fixing p we represent these limits implicitly.
The model is detailed in Box S1. Crucially, the dependencies on a of invasion probability and of the mean lifetime of populations are chosen to reproduce those of the deconstructed formulation of our food-web model (Fig. 3). We chose pN = 1000 over a range of p values, evaluated the algorithm over 4 · 10 7 iterations, and sampled base attack rates from the last 3/4 of each run to characterise the steady state (which was reached after less than a 10 th of iterations).
In the limit p → 0, where cheating does not occur, the metapopulation model should recover the mean base attack rate in the evolutionary steady state of our food-web model. Indeed, the metapopulation model attained a steady state with mean logarithmic base attack rate log 10 a = −5.14, close to the value of −4.96 obtained with the deconstructed model formulation. We also obtained an approximately normal distribution of log 10 a in the steady state of the metapopulation model similar to that in Fig. 3e. These results provide evidence that our reconstruction of the fitness landscape in Fig. 3 does indeed represent the fitness landscape experienced by an evolving metapopulation of a single species.
As shown in Fig. S3, log 10 a increases linearly with p for low p. An occupancy of p = 0.3, for example, leads to an approximate 3-fold increase in geometric mean base attack rates. Hence, cheating makes consumers somewhat less prudent, but does not fundamentally undermine the evolution of prudent predation.  Figure S3: The impact of cheaters on evolutionary stable base attack rate a. Simulation results from the metapopulation model described in Box S1. The higher the occupancy p of patches by the metapopulation, the larger the probability that occupied patches are overtaken by invading cheaters with higher base attack rates. This effect increases steady state base attack rates but does not prevent a steady state from being reached.

Appendix S7 Steepness and basic reproduction number
We derive the relation between basic reproduction number and the steepness of stock recruitment relations given in Eq. (12).
Consider first the following caricature model of a fish stock with standing stock biomass SSB that feeds on a single resource: The parameter F denotes the fishing mortality rate, otherwise model structure and parameterization are as in Eq. (1). If one assumes, for simplicity, that (i) all mature individuals have the same body mass m, (ii) recruits are produced instantaneously, and (iii) the parameter ρ is dominated by natural mortality rather than respiration, then recruitment is given by the first term on the right-hand side of Eq. (S35b): In the second step we eliminated B R by solving Eq. (S35a) with dB R /dt = 0 for B R > 0. Stock-recruitment relations of this quadratic form are frequently used in fisheries science and named after Schaefer (1954).
The basic reproduction number R is defined as recruitment per mature individual (of which there are SSB/m) in units of ρ, in the limit SSB → 0, which evaluates to Hence Eq. (S38) implies Eq. (12). We now verify that Eq. (12) remains valid if one generalises Eq. (S35) to a situation with multiple resources. We assume that the fish stock is initially fully established at SSB 0 , such that resources that would not withstand its consumption have been extirpated. By Eq. (S6), the biomass of each resource is then a linear function of consumer biomass, here SSB. With the linear functional response of Lotka-Volterra models, this implies m Rec = (c 1 − c 2 SSB) SSB (S40) with some positive constants c 1 and c 2 . As above, we can evaluate yielding c 1 = ρR. Furthermore, recruitment balances mortality for an unfish stock with SSB = SSB 0 . So m Rec(SSB 0 ) = ρ SSB 0 , which implies With these values for c 1 , c 2 , plugging Eq. (S40) into the definition of steepness, Eq. (11), yields again Eq. (12).  Figure S4: Activation frequency of consumer extirpation mechanisms. The classification relates to different steps in the deconstructed model formulation (Box 1). Exploitative competition refers to Step 5; Pyrrhic competition to failure to meet the invasibility condition by a consumer losing its main resource, or by the competitor causing this, in Steps 6, 7; Bust after boom refers to Step 3d; and Invasibility criterion to failure to satisfy Eq. (13) at any other point in the algorithm. Extirpations through Pyrrhic competition are very rare, and those through bust after boom contribute just a few percent of cases.