Stochastic competitive exclusion leads to a cascade of species extinctions
SStochastic competitive exclusion leads to a cascade of species extinctions
Jos´e A. Capit´an a,c, ∗ , Sara Cuenda b , David Alonso c a Department of Applied Mathematics, Technical University of Madrid, Madrid, Spain b Departamento de An´alisis Econ´omico: Econom´ıa Cuantitativa, Universidad Aut´onoma de Madrid, Madrid, Spain c Center for Advanced Studies (CEAB-CSIC), Blanes, Catalunya, Spain
Abstract
Community ecology has traditionally relied on the competitive exclusion principle, a piece of common wisdom in conceptualframeworks developed to describe species assemblages. Key concepts in community ecology, such as limiting similarity and nichepartitioning, are based on competitive exclusion. However, this classical paradigm in ecology relies on implications derived fromsimple, deterministic models. Here we show how the predictions of a symmetric, deterministic model about the way extinctionsproceed can be utterly different from the results derived from the same model when ecological drift (demographic stochasticity)is explicitly considered. Using analytical approximations to the steady-state conditional probabilities for assemblages with twoand three species, we demonstrate that stochastic competitive exclusion leads to a cascade of extinctions, whereas the symmetric,deterministic model predicts a multiple collapse of species. To test the robustness of our results, we have studied the effect ofenvironmental stochasticity and relaxed the species symmetry assumption. Our conclusions highlight the crucial role of stochasticitywhen deriving reliable theoretical predictions for species community assembly.
Keywords:
Competitive exclusion, Ecological drift, Continuous-time Markov processes
1. Introduction
Ecological communities are shaped from the complex inter-play of four fundamental processes (Vellend, 2010): selection,in the form of species interactions that favor certain speciesagainst others; speciation, leading to the appearance of newspecies, better adapted to the environment; dispersal, whichpermits spatial propagation of individuals; and ecological drift,a demographic variability in species population numbers due tothe stochastic processes that take place. Ecological drift, in par-ticular, has a prevalent role in modern theoretical frameworks incommunity ecology (Black and McKane, 2012). Accordingly,current approaches reveal the need for process-based, stochas-tic models that help to understand how ecological communi-ties are assembled and their interaction with environmental fac-tors (Wisz et al., 2013).Classical community ecology, however, has mainly relied ondeterministic community models (see Roughgarden (1979) andreferences therein), most of them based on Lotka-Volterra dy-namics, although alternatives have been proposed (Schoener,1974a). There is a long-standing research focus on communityassembly models, in which communities are built up throughspecies invasions, and most of them rely on deterministic ap-proaches (Post and Pimm, 1983; Law and Morton, 1993, 1996;Capit´an et al., 2009; Capit´an and Cuesta, 2011; Capit´an et al.,2011). On the other side, there have been strong theoreti-cal efforts to describe community assemblages in stochastic ∗ Corresponding author.
Email addresses: [email protected] (Jos´e A. Capit´an), [email protected] (Sara Cuenda), [email protected] (David Alonso) terms (Hubbell, 2001; Alonso et al., 2008; Rosindell et al.,2011). In certain situations, the results and conclusions derivedfrom deterministic models have been shown to be quite differ-ent in the presence of stochasticity (Bolker and Grenfell, 1995;Alonso et al., 2007; Haegeman and Loreau, 2011; Bonachelaet al., 2012; Wang et al., 2012).One of the contexts where the differences between deter-ministic and stochastic approaches become apparent is relatedto theoretical formulations of the competitive exclusion princi-ple (Volterra, 1926; Gause, 1934; Hardin, 1960). This principleconstitutes a fundamental pillar of community ecology and be-longs to the traditional body of ecological theory. It provides auseful theoretical framework to explore how complex speciesassemblages persist over time. Important concepts such asadaptation to shared niches (Roughgarden, 1979), species lim-iting similarity (MacArthur and Levins, 1967; Roughgarden,1974) or niche partitioning (Pielou, 1977; Schoener, 1974b)all are immediate derivations of the principle. Classical ap-proaches predict the maximum degree of species similarity thatpermit species stable coexistence (MacArthur, 1969, 1970).However, theoretical predictions for limiting similarity of-ten rely on deterministic community models (see MacArthur(1968); Levin (1970); Haigh and Maynard-Smith (1972); Ches-son (1990) and Appendix A for a discussion on competitive ex-clusion based on deterministic approaches), and the relevanceof stochasticity, in the form of ecological drift, to species co-existence has remained almost unexplored (with the exceptionof Turelli (1980)). The relationship between limiting similar-ity and environmental stochasticity has been studied more thor-oughly (May and MacArthur, 1972; Turelli, 1978, 1981).
Preprint submitted to Journal of Theoretical Biology August 15, 2016 a r X i v : . [ q - b i o . P E ] A ug ecently, we focused on the influence of ecological drifton the similarity of coexisting species via the competitive ex-clusion principle (Capit´an et al., 2015). In that contributionwe showed that, in the presence of ecological drift, the max-imum degree of similarity that ensures stable coexistence canbe significantly lowered when compared to the correspond-ing limits to similarity derived from deterministic models.If similarity is interpreted in terms of an interspecific com-petitive overlap (MacArthur and Levins, 1967; Roughgarden,1974), stochasticity displaces the deterministic threshold to-wards lower values of the competitive overlap (Capit´an et al.,2015). Thus, when stochasticity is considered, the extinc-tion phenomena caused by competitive exclusion takes place atlower values of the competitive overlap (i.e., species have to bemore dissimilar to stably coexist in the presence of ecologicaldrift).Ecological drift becomes a key process determining speciescoexistence in aspects other than the maximum similarity ofco-occurring species. Beyond a more restrictive threshold incompetition induced by ecological drift (which was the mainresult of Capit´an et al. (2015)), we here analyze the influence ofdemographic stochasticity on the extinction mechanism itself,which in principle can lead to either sequential or grouped ex-tinctions as competition strength increases. For that purpose,we considered a deterministic, Lotka-Volterra model and itsstochastic counterpart, both of which treat species interactionssymmetrically. Whereas the deterministic model predicts themultiple extinction of all the species in the community but oneas competition crosses over a certain threshold, in the presenceof demographic stochasticity extinctions proceed progressively,in the form of a cascade, as competition increases. The only dif-ference between both approaches is the explicit considerationof ecological drift in the dynamics. In order to derive our con-clusions, we developed convenient analytical approximations tothe steady-state configurations of the stochastic system for sim-ple species assemblages formed by two or three species. Suchapproximations help us to partition the set of feasible popu-lation numbers into regions associated to coexistence, or theextinction of one, two, or three species. The steady-state proba-bilities, when aggregated over those regions, unveil the extinc-tion cascade phenomenon. Our main result reveals overlappingwindows in competitive strength, at low values related to con-figurations where the coexistence of three species is the mostprobable state, intermediate ranges where it is more likely toobserve two-species assemblages, and large competition valuesfor which the most probable state is formed by one species ornone. We also studied the transition to the deterministic modelwhen demographic stochasticity tends to zero, and our resultsreveal an abrupt transition to situations compatible with smallstochasticity.To test the robustness of our conclusions, we replaced demo-graphic stochasticity by environmental stochasticity and con-firmed that, although the extinction phenomena are qualita-tively different, the extinction cascade persists. We also relaxedthe assumption of symmetry to assess the effect of stochastic-ity on deterministic models that not only predict multiple ex-tinctions, as in the fully symmetric scenario, but also lead to both progressive and grouped extinctions for fixed competi-tive strengths. When stochasticity comes into play, however,the stochastic cascade persists and the expected extinction se-quence is qualitatively different from its deterministic counter-part. Thus, the predictions of both models are significantly dif-ferent in generic, non-symmetric scenarios for species interac-tions.The paper is organized as follows: in Section 2 we describeboth the deterministic and the stochastic frameworks, the latterbased on the formulation of Haegeman and Loreau (2011), andshow that the deterministic, symmetric model predicts a multi-ple species extinction. In Section 3 we start by presenting theanalytical approximations for a two-species stochastic commu-nity model, and we then extend the procedure to a three-speciescommunity. These approximations help us to obtain analyticalformulae for the critical points of the steady-state, joint proba-bility distribution of the community. Formulae for saddle pointsare then used to properly define aggregated probabilities of co-existence, or one-, two-, and three-species extinction, whichreveal themselves the sequential decline of species driven byecological drift. After studying the small stochasticity limit andtesting the robustness of our results, we conclude the paper withseveral implications and prospects (Section 4).
2. Model description
For the sake of simplicity, in this contribution we will focuson the symmetric version of the deterministic Lotka-Volterracompetitive dynamics (see Appendix A),˙ x i = rx i (cid:32) − x i + ρ (cid:80) j (cid:44) i x j K (cid:33) , i = , . . . , S , (1)where x i stands for the population density of species i (spaceis implicitly assumed) and model parameters are uniform andspecies-independent. Here r stands for an intrinsic, species-independent growth rate, ρ measures interspecific competition, K represents a carrying capacity, and S is the species richness ofthe community. The dynamics has an interior equilibrium point, ˆx = ( ˆ x , . . . , ˆ x ), where ˆ x = K / (1 − ρ + ρ S ), which is globally sta-ble if and only if ρ < ρ >
1, all the species become extinct except forjust one of them (see Appendix B). As a consequence, compet-itive exclusion in the symmetric, deterministic model impliesthe joint extinction of S − n i at2ime t , n ( t ) = ( n ( t ) , . . . , n S ( t )). Contrary to the deterministiccase, which focuses on population densities x i = n i / Ω , Ω be-ing a meaningful measure (area, volume) of the system size,here discrete population numbers are considered. The elemen-tary processes that define the stochastic dynamics (local birthsand deaths, immigration, and competition) are characterized byprobability rates that, in the deterministic limit, yield the Lotka-Volterra equations (1). As in Haegeman and Loreau (2011), wechoose the following probability rates to model elementary pro-cesses:1. Local births (deaths) of species i occur at a density-independent rate r + n i ( r − n i ). We adopt the notation r = r + − r − to represent the net growth rate in the absence ofcompetitors.2. Immigration of a new individual of species i takes place ata rate µ . Although the deterministic model (1) does not in-clude immigration, dispersal is an important process driv-ing community assembly (Vellend, 2010). In addition, im-migration is key for the stochastic process to reach a non-trivial steady-state. We consider here the low-immigrationregime, in which the deterministic limit is expected to re-cover results close to those yielded by Eq. (1), see Capit´anet al. (2015).3. Intraspecific competition occurs at a density-dependentrate rn i / K , where K represents a species carrying capac-ity.4. Interspecific competition between species i and j ( i (cid:44) j )takes place at a probability r ρ n i n j / K per unit time (it isalso a density-dependent rate).Population vectors n ( t ) belong to the configuration space N S ,where N = { , , , . . . } . The elementary processes listed abovedefine a birth-death-immigration stochastic process in contin-uous time, and the probability P ( n , t ) of finding a populationvector n ( t ) at time t is determined by the master equation, ∂ P ( n , t ) ∂ t = S (cid:88) i = (cid:8) q + i ( n − e i ) P ( n − e i , t ) + q − i ( n + e i ) P ( n + e i , t ) − [ q + i ( n ) + q − i ( n )] P ( n , t ) (cid:9) . (2)Here q + i ( n ) = r + n i + µ is the overall birth rate for species i , q − i ( n ) = r − n i + rn i ( n i + ρ (cid:80) j (cid:44) i n j ) / K is the overall death rate,and e i = (0 , . . . , , . . . ,
0) is the i -th vector of the canonical basisof R S . Notice the correspondence of these rates with the termsarising in the deterministic model (1) for µ =
0. The steady-state probability distribution is obtained by solving the cou-pled recurrence equation given by the condition ∂ P ( n , t ) /∂ t = S = S >
1, in the absence of competition ( ρ =
0) populations are uncor-related, and the joint probability distribution factors as a prod-uct of marginal probabilities, which reduces this case to a one-dimensional problem. The fully neutral model ( ρ =
1) for S > S ≥ ρ >
0, though. Haegeman and Loreau (2011) devised approx-imations for this case, but we will follow here a different ap-proach to find analytical formulae for the critical points of thesteady-state distribution in the case of small-sized communities.
3. Results
In Capit´an et al. (2015) we showed numerically that thesteady-state distribution for two-species communities presentsa maximum at an interior point of the configuration space, aswell as two boundary maxima with population vectors of theform ( n ,
0) and (0 , n ). This implies that the discrete probabilitydistribution, when extended to be real-valued (by, for instance,cubic spline interpolation), must exhibit by continuity two sad-dle points located in between the coexistence maximum andthe two boundary maxima. These saddle points can be used toconveniently partition the configuration space into regions as-sociated to coexistence, the extinction of one species, or the ex-tinction of two species. In this section we develop analytical ap-proximations that help us to obtain good estimates for the criti-cal points of the joint probability distribution. We use the S = To estimate the location of the critical points of the joint dis-tribution P ( n , n ) we use that the conditions ∂ P ( n , n ) /∂ n = ∂ P ( n , n ) /∂ n = ∂∂ n P ( n | n ) = ∂∂ n (cid:32) P ( n , n ) P ( n ) (cid:33) = ,∂∂ n P ( n | n ) = ∂∂ n (cid:32) P ( n , n ) P ( n ) (cid:33) = , (3) P ( n | n ) being the probability that species 1 has n individualsconditioned to species 2 having n individuals. This means thatcritical points of the joint distribution P ( n , n ) can also be ob-tained through conditional probabilities. By fixing n , we justneed to evaluate the derivative of P ( n | n ) along the n direc-tion. The same applies under the change n ↔ n . We now approximate P ( n | n ) by T ( n | n ) as follows. For S =
2, the steady-state distribution satisfies a two-term recur-3ence in population numbers n and n , namely0 = q + ( n − , n ) P ( n − , n ) + q − ( n + , n ) P ( n + , n ) + q + ( n , n − P ( n , n − + q − ( n , n + P ( n , n + − (cid:2) q + ( n , n ) + q − ( n , n ) (cid:3) P ( n , n ) − (cid:2) q + ( n , n ) + q − ( n , n ) (cid:3) P ( n , n ) . (4)Notice that we are interested in approximating the conditionalprobability P ( n | n ) where n is fixed. Hence, in the approxi-mation we ignore the terms in Eq. (4) that involve variation of n , and assume that n acts as a fixed parameter in the remainingterms. Thus, the approximated conditional probability T ( n | n )satisfies0 = q + ( n , n − T ( n − | n ) + q − ( n , n + T ( n + | n ) − (cid:2) q + ( n , n ) + q − ( n , n ) (cid:3) T ( n | n ) . (5)This expression fulfills a detailed balance condition (Karlin andTaylor, 1975), which yields an approximate one-term recur-rence formula in n , T ( n | n ) = q + ( n , n − q − ( n , n ) T ( n − | n ) = K [ µ + r + ( n − n [ Kr − + r ( n + ρ n )] T ( n − | n ) , (6)and leads to an explicit solution in terms of hypergeometricfunctions, as in Haegeman and Loreau (2011). A symmetricrecurrence holds for T ( n | n ). Our next step is to approximate the partial derivative alongthe n direction by the backwards discrete difference (compara-ble results are obtained with the forward difference), ∂∂ n P ( n | n ) ≈ ∂∂ n T ( n | n ) ≈ T ( n | n ) − T ( n − | n ) . (7)Similarly, ∂∂ n P ( n | n ) ≈ T ( n | n ) − T ( n − | n ) . (8)Eq. (6) allows us to write the system for the critical points ofthe two-dimensional joint-probability surface, ∂ P ( n | n ) /∂ n = ∂ P ( n | n ) /∂ n =
0, as n + ρ n n − Kn + K α = , n + ρ n n − Kn + K α = , (9)where α = ( r + − µ ) / r . Solving the quadratic system yieldsthe following estimates for the interior critical points: M = ( m + , m + ), which can be either a local maximum or a saddlepoint depending on the value of ρ (see below), and the localminimum M = ( m − , m − ), where m ± = K + ρ ) ± (cid:114) − α (1 + ρ ) K , (10) and two (symmetrical) saddle points Q = ( s + , s − ) and Q = ( s − , s + ), where s ± = K ± (cid:115) − α K (1 − ρ ) . (11)To test the goodness of these analytical approximations, crit-ical points of the exact probability distribution are determinednumerically by extending the discrete distribution P ( n , n ) tobe a real-valued function using cubic spline interpolation bothin the n and n directions. We then solve numerically the sys-tem ∂ P ( n , n ) /∂ n = ∂ P ( n , n ) /∂ n =
0, using analytical pre-dictions as initial guesses for iterative root finding. The Hessianmatrix decides whether a given point is a local maximum, a lo-cal minimum or a saddle point. We find that the critical pointsare in excellent agreement with the analytical formulae above(see Fig. 1a, in which we plot the coordinates of each criticalpoint calculated both analytically and numerically).For ρ < ρ c , with ρ c = − α K = − r + − µ ) Kr , (12)the analysis of the Hessian matrix reveals that M is a (coexis-tence) local maximum. For ρ = ρ c , saddle points Q and Q coincide with M and, for ρ > ρ c , M transforms into a saddlepoint and the former saddle points Q and Q (cf. Eq. (11)) nolonger exist. Moreover, when ρ > K / α −
1, all interior criticalpoints (10) become complex and the only persistent maxima arethose located at the boundary.There are four critical points at the boundary, which can beapproximated using Eq. (6) for n =
0, resulting in two localmaxima [( b + , , b + )] and two local minima [( b − , , b − )],where b ± = K ± (cid:114) − α K . (13)Observe that the non-zero coordinates of the local maxima(minima) coincide with those of the interior maximum (mini-mum) for ρ = The method proposed in the previous subsection can be fullyextended to the case of a three-species community. Criticalpoints are obtained by approximating conditional probabilities P ( n | n , n ) and taking discrete derivatives with respect to thefirst argument. As before, we consider the three-term recur-rence relation that fulfills the joint distribution P ( n , n , n ) andignore the terms that involve variation in population numbers n and n . Under this approximation, the steady-state conditionreduces to0 = q + ( n − , n , n ) T ( n − | n , n ) + q − ( n + , n , n ) T ( n + | n , n ) − (cid:2) q + ( n , n , n ) + q − ( n , n , n ) (cid:3) T ( n | n , n ) . (14)Due to detailed balance, the approximate conditional probabil-ities T ( n | n , n ) satisfy the one-term recurrence relation K (cid:2) µ + r + ( n − (cid:3) T ( n − | n , n ) = n (cid:2) Kr − + rn + r ρ ( n + n ) (cid:3) T ( n | n , n ) . (15)4 Competitive overlap, ρ t + t - s + s - m + m - cd Competitive overlap, ρ P opu l a t i on nu m be r , n t + t - s + s - m + m - ab Figure 1: (a) Coordinates of the critical points for S = ρ . They were calculated using the exact steady-state distribution (symbols) as well as theapproximations given by Eqs. (10) and (11) (lines). Once the coexistence maximum M = ( m + , m + ) and the two saddle points Q i = ( s ± , s ∓ ) ( i = ,
2) coalesce, M becomes a saddle point (there is no longer a coexistence maximum). Model parameters are r + = r − = µ = K =
25. (b) Coordinates of the critical pointsfor a three-species community (symbols), compared with theoretical approximations (17) and (19). When the two saddle points coincide, the former maximumbecomes a saddle point. Parameter values are the same as in panel (a). Panels (c)-(d) check the accuracy of our approximations for other set of parameter values,namely r + = r − = . µ = .
1, and K =
25. As before, (c) represents the coordinates of the critical points for S = S = Repeating for S = n + ρ n ( n + n ) − Kn + K α = , n + ρ n ( n + n ) − Kn + K α = , n + ρ n ( n + n ) − Kn + K α = , (16)which yields 8 interior critical points, 6 of which are saddlepoints, and the remaining two points are a local minimum and,as in S =
2, a point that is a maximum or a saddle point depend-ing on ρ . Explicit expressions for their coordinates are givenbelow. For the sake of comparison we have also calculated crit-ical points using the exact joint distribution P ( n , n , n ), seeFig. 1b.The coordinates for the interior critical points are: on the onehand, M = ( m + , m + , m + ) and M = ( m − , m − , m − ), where m ± = K + ρ ) ± (cid:114) − α (1 + ρ ) K . (17)Both solutions turn out to be complex when ρ > (cid:18) K α − (cid:19) . (18) On the other hand, six interior critical points appear at points Q i , i = , . . . ,
6, where Q = ( t + , t + , s − ) and Q , Q are ob-tained as the cyclic permutations of the entries of Q , whereas Q = ( t − , t − , s + ) and the entries of Q and Q are the cyclicpermutations of that of Q , with t ± = K + ρ ) ± (cid:115) − α (1 + ρ ) K (1 − ρ ) , s ± = K ± (cid:115) − α (1 + ρ ) K (1 − ρ ) . (19)Both s ± and t ± are real whenever ρ ≤ ρ c , where ρ c = − r + − µ ) / Kr + r + − µ ) / Kr . (20)In spite that Eqs. (17) and (19) have been obtained using an ap-proximate form for conditional probabilities, the numerical cal-culation of the critical points using the actual distribution is invery good agreement with these approximations (see Fig. 1b).For ρ < ρ c , M is classified as a local maximum. At ρ = ρ c all saddle points and M coalesce in a single point. For ρ > ρ c ,however, M becomes a saddle point, as can be checked numer-ically with the Hessian matrix.5 igure 2: Partitioning of the configuration space N for S = ρ = . ρ = . ρ = .
45 and ρ = .
6; remaining parametersare r + = r − = µ =
1, and K =
50. The smaller coordinates of the two saddle points (used to draw white lines) are used to partition the space into regionsassociated to coexistence (central square), one-species extinction (the two rectangles which contain the boundary maxima) and two-species extinction (the smallsquare that contains the origin). As competition increases, the coexistence maximum approaches to the origin, and saddle points become closer to the maximum.
Boundary maxima are of the form ( n , n , n , ,
0) —andtheir cyclic permutations. The non-zero coordinates of the for-mer are equal to that of the coexistence maximum M obtainedfor 2 species, see Eq. (10); the non-zero entries of the latter arethe same as the boundary maxima for S =
2, see Eq. (13).
For S = s − ) that closestto the boundary (cf. Eq. (11)), the partitioning results as:1. 0 ≤ n i < s − for i = ,
2. This square is associated with fullextinction.2. 0 ≤ n < s − , n > s − or 0 ≤ n < s − , n > s − . Thesetwo rectangles are associated with the extinction of onespecies, since configurations are close to extinction in theform ( n ,
0) or (0 , n ).3. n i > s − for i = ,
2. The rest of the configuration spaceputs together states that can be associated to coexistence.Note that, strictly speaking, there will be configurations whereboth n , n > S = s − and t − in Eq. (19)). Note also that s − > t − . Taking this fact into account, saddle points divide the configuration space N into several regions, which have beendepicted in Fig. 3:1. 0 ≤ n i < t − for i = , ,
3. This cube is associated with theextinction of the three species (since the origin belongs tothis region).2. 0 ≤ n i < t − for i = , n > t − (and the two remainingcombinations). These three parallelepipeds are associatedwith the extinction of two species, because boundary max-ima of the form ( n , ,
0) —and its cyclic permutations—are situated inside this volume, as well as other populationconfigurations close to two-species extinctions.3. n i > s − for i = , ,
3. This cube contains the interiormaximum and is therefore associated with the coexistenceof the three species.4. The remaining volume of the configuration space, whereboundary maxima of the form ( n , n ,
0) —and its cyclicpermutations— are located, is associated with the extinc-tion of one species.The configuration-state partitioning slightly differs from thegeneral case when saddle points have coalesced. If ρ > ρ c , thepoint M = ( m + , m + , m + ) is classified as the only saddle point.Based on its coordinates, the configuration space is partitionedas follows:1. 0 ≤ n i < m + for i = , ,
3. This cube is associated withfull extinction.2. 0 ≤ n i < m + for i = , n > m + (and the two re-maining combinations). These three parallelepipeds areassociated with the extinction of two species.3. n i > m + for i = , ,
3. This cube is associated with coex-istence configurations.4. The remaining volume of the configuration space is relatedto the extinction of one species.The same partitioning applies for the S = igure 3: Partition of the configuration space N for S = N into four regions according to sad-dle points: black circles denote points Q , Q and Q , whereas white squaresrepresent saddle points Q , Q and Q . The latter are used to define regions forcomplete (red) and two-species extinctions (blue), and the former determine thecoexistence volume (green) and the one-species extinction region (yellow). The(green) cube has been displaced vertically to facilitate visualization. Saddle points of the joint probability distribution have al-lowed us to establish a natural partitioning into regions associ-ated to coexistence (containing the coexistence maximum) andto the extinction on one, two, or three species (containing thecorresponding boundary maxima), see Fig. 3. Over these re-gions we aggregate the joint distribution P ( n , n , n ), calcu-lated numerically as described in Appendix C, to define theoverall probability of coexistence, P c = (cid:88) n , n , n > s − P ( n , n , n ) , (21)the probability of three-extinct species configurations, P = (cid:88) n , n , n < t − P ( n , n , n ) , (22) and the probability of two-extinct species, which by symmetryover cyclic permutations of ( n , n , n ) can be expressed as P = (cid:88) n , n < t − (cid:88) n > t − P ( n , n , n ) . (23)The probability of one-extinct species, P , is obtained from thenormalization condition P c + P + P + P =
1. Fig. 4 showsthese aggregated probabilities as a function of ρ for two sets ofmodel parameters. In the first case, the coexistence probabilityis almost one for low values of the competitive overlap and, as ρ increases, at some point the probability declines rapidly. Atthe same time, the probabilities of one and two extinct speciesbegin to increase. Note that, once the threshold ρ c has beencrossed over, the most probable state consists of a single, extantspecies, and the probability of coexistence becomes negligible.In the second case, corresponding to a larger value of themortality rate r − , the threshold ρ c at which the coexistence max-imum M transforms into a saddle point becomes smaller. Theprobability of coexistence rapidly declines as ρ increases and,in addition, there is a non-negligible probability of completeextinction. Remarkably, Fig. 4b shows that, for smaller valuesof the carrying capacity, coexistence is not the only possiblestate even at ρ =
0. This puts a practical limit to the maximumnumber of coexisting species which does have a deterministiccounterpart —recall that, due to global stability, the determinis-tic model permits the packing of an arbitrary number of speciesfor ρ < S − ρ increases (Appendix B), Fig. 4 showsthat ecological drift induces a sequential cascade of extinctions,in which states with a larger number of extinct species are moreprone to be observed as the competitive overlap increases.An important remark is on purpose here. The cascade ofextinction we have just described has nothing to do with thedegree of synchronicity in which extinctions take place alongtime, i.e., the term “cascade” does not refer here to a sequen-tial extinction in time. In particular, the symmetric, determin-istic model leads in general to asynchronous extinctions. Thestochastic phenomenon analyzed in this contribution refers tothe progressive extinctions that occur as competitive strengthincreases. In the absence of stochasticity, the deterministic model pre-dicts the extinction of S − ρ = Competitive overlap, ρ P r obab ili t y coexistence1 extinct2 extinct3 extinct a b Figure 4: (a) Stochastic extinction cascade. The panel shows the dependence of aggregated probabilities over the regions determined by the critical points as afunction of competitive overlap ρ . Parameter values are S = r + = r − = K =
25, and µ =
1. Although initially the three species coexist, at intermediatevalues of ρ the most probable state is formed by only two extant species. For large competition values the most likely configuration comprises a single extantspecies. A vertical, dashed line marks the value ρ c (cf. Eq. (20)) at which the probability of coexistence is negligible. (b) Same as panel a but for a higher intrinsicmortality rate, r − =
35. In this case, the aggregated probability of complete extinction is non-negligible even for ρ =
0. Aggregated probability for one-extinctspecies configurations starts declining and, at the same time, the two-extinct species configurations become more likely as ρ increases. Close to ρ = ρ c . The dot-dashed lineshows the value given by Eq. (18), at which the two critical points M and M no longer exist. Finally, a dotted, vertical line marks the value ( ρ = K / α −
1, seeSection 3.1) at which the boundary maxima of the form ( n , n ,
0) —and permutations— no longer exist. two scenarios as stochasticity decreases. To do so, we havestudied the limit of small stochasticity, in which the determin-istic model is to be recovered. As shown below, the transitionto the small-stochasticity scenario is abrupt, hence incorporat-ing demographic stochasticity to community models should bestrongly considered.Since fluctuations are expected to decrease as population sizeincreases (see Appendix D), the small-stochasticity limit isequivalent to the limit of large population sizes, so we have re-peated the analysis by increasing the carrying capacity at fixed ρ . We have quantified the intensity of demographic noise by thecoefficient of variation of population abundances, ν = σ n / (cid:104) n (cid:105) , σ n being the standard deviation of population numbers and (cid:104) n (cid:105) the average value. As shown in Appendix D, when K (cid:29) σ n ∼ K / and (cid:104) n (cid:105) ∼ K , so ν tends to zero in the limit of largepopulation sizes. From the numerical point of view, to get closeto the deterministic scenario we would have to choose a car-rying capacity value such that the average (cid:104) n (cid:105) is large enoughcompared to the variability in populations. In practice, we willassume that the system is close to a low-noise regime when theactual coefficient of variation, obtained though the joint prob-ability distribution of the stochastic model calculated numeri-cally, is close to that obtained by a Gaussian approximation ofthe joint distribution valid in the limit K (cid:29) K , the probabil-ity of coexistence is almost equal to one (as in the deterministicscenario). However, as ν increases (i.e., the carrying capacity K decreases), the probability that one species becomes extinctgrows sharply and, at the same time, the probability of coex-istence starts declining. As the variability in population sizesaugments, overlapping windows for the likelihood of progres-sive extinctions arise. The transition to situations where lowdemographic stochasticity operates is, therefore, abrupt.Note that the cascade obtained in Fig. 5b,d as a function ofnoise can be immediately translated into a cascade in carryingcapacity (for fixed ρ ). This reinforces our conclusion, since theextinction cascade phenomenon also occurs when other modelparameter ( K ) varies. It is presumably the relative balance be-tween ρ and K that determines the subspace of the parameterspace for which extinctions start appearing. In order to test the robustness of our main result against dif-ferent sources of noise, we have replaced demographic stochas-ticity by environmental stochasticity. We have introduced vari-ability in model parameters so that, to keep the scenario as sim-ple as possible, the competitive overlap ρ in Eq. (1) is replacedby ρ + ξ ( t ), where ξ ( t ) stands for a noise term with zero mean.The deterministic dynamics transforms into a Langevin equa-tion with multiplicative noise,˙ x i = rx i (cid:32) − x i + ρ (cid:80) j (cid:44) i x j K (cid:33) − rx i K (cid:18) (cid:88) j (cid:44) i x j (cid:19) ξ ( t ) , (24)8 Carrying capacity, K C oe ff i c i en t o f v a r i a t i on , ν Coefficient of variation, ν P r obab ili t y c d Carrying capacity, K C oe ff i c i en t o f v a r i a t i on , ν analyticalnumerical Coefficient of variation, ν P r obab ili t y coexsistence1 extinct2 extinct3 extinct a b Figure 5: Cascade persistence for variable levels of demographic stochas-ticity. Panel (a) shows the coefficient of variation of population abundances, ν = σ n / (cid:104) n (cid:105) (diamonds), as a function of increasing carrying capacity, K , for S = r + = r − = µ = . ρ = .
5. For large K , the coeffi-cient of variation tends to the analytical approximation (D.8) derived in Ap-pendix D. Simultaneously, the probability of coexistence (alternative verticalaxis, circles) becomes closer to 1. (b) Extinction cascade as a function of demo-graphic stochasticity. At ν ∼ .
2, the probability that one species goes extinctabruptly increases, and the coexistence probability starts declining. As ν in-creases, higher-order extinctions become more likely. Panels (c)-(d) are equiv-alent to (a)-(b) but for S = r + = r − = µ = . ρ = .
1. Whenthe Gaussian approximation is valid, coexistence probability is almost equalto 1. After an abrupt decline, sequential extinctions occur for higher levels ofdemographic noise. which can be rewritten as z (cid:48) i = z i (cid:18) − z i − ρ (cid:88) j (cid:44) i z j (cid:19) − z i (cid:18) (cid:88) j (cid:44) i z j (cid:19) ξ ( t ) (25)by re-scaling species densities as z i = x i / K and time as t (cid:48) = rt ( z (cid:48) i stands for the derivative with respect to the scaled time t (cid:48) ).We choose the noise as ξ ( t ) = ρκη ( t ), where η ( t ) is a Brownianmotion. The noise ξ has been scaled by ρ in order to avoid thatthe overall competitive strength, ρ + ξ ( t ), becomes negative.Fig. 6 shows the persistence of the extinction cascade whenonly environmental stochasticity is considered. This points tothe robustness of our results: as for demographic stochastic-ity, environmental stochasticity also alters the predictions of thedeterministic dynamics. There are, however, qualitative differ-ences between the predictions yielded by the model when de-mographic or environmental stochasticity come into play. First,the range in competition on which the cascade takes place iswider for demographic noise. It seems that, in the presenceof environmental noise, the range of the cascade can increase moderately when a larger number of species are to be packed(Fig. 6b). More importantly, a second difference arises: no fullextinction is possible in the case of environmental noise. Thisis a peculiarity, not altered by the noise, of a generic competi-tive, Lotka-Volterra dynamics (cf. Eq. (A.2) in Appendix A),for which it can be easily shown that the full extinction equi-librium ( ˆ x i = i = , . . . , S ) is unstable. The Langevin equa-tion, therefore, can not reproduce configurations associated tothe full extinction of the community, contrary to what is ob-served for demographic stochasticity (Fig. 4b). In this section we test the robustness of our results by relax-ing the assumption of species symmetry in model parameters. Itcan be argued that the effect of demographic stochasticity sim-ply consists of breaking the symmetry between species. In ageneric, non-symmetric, deterministic scenario, one could ex-pect progressive extinctions even in the absence of ecologicaldrift. Therefore, the role of stochasticity would be simply tore-establish a deterministic scenario where one-by-one speciesextinctions occur. In this subsection we discuss the implicationsof relaxing the symmetry assumption to determine the true roleof demographic stochasticity in a generic case.In order to address these questions, here we consider two ex-amples of fully non-symmetric, three-species competitive dy-namics, ˙ x i = rx i − K i (cid:88) j = ρ i j x j , (26)where carrying capacities and interspecific competitivestrengths are species-dependent. We set, without loss of gen-erality, off-diagonal competition values as ρ = ρ = ρ , ρ = ρ = ρ + δ , ρ = ρ = ρ + δ , and ρ ii = i = , , δ = . δ = . K = K =
16, and K =
20. After analyzing the stability ofall the equilibrium points of the deterministic system (see de-tails in Appendix E) we find that, although being fully non-symmetric, this model predicts a two-species, grouped extinc-tion when the threshold ρ = . ρ > .
4, only equilibria with a single extant species are stable,whereas no other scenario but three-species coexistence is sta-ble for 0 ≤ ρ < .
4. By continuity of the eigenvalues of theJacobian matrix on model parameters, there are multiple non-symmetric, deterministic scenarios that exhibit the transitionfrom full coexistence to a single-extant-species state. There-fore, grouped extinctions are not a peculiarity of the fully sym-metric scenario.In the second example (Fig. 7b), δ = . δ = . K = K = K =
15. The complete stability analysis draws theconclusion that non-symmetric models yield to multistability.In Appendix E it is shown that, for ρ < .
3, only the interiorequilibrium point is asymptotically stable; for 0 . < ρ < . / (1 + ρ ) , / (1 + ρ ) ,
0) is stable, butfor 0 . < ρ <
1, the former boundary point remains stable aswell as the single-species extinction equilibrium (0 , ,
15) (see9 .8 0.9 1
Competitive overlap, ρ coexistence1-4 extinct5-8 extinct9 extinct P r obab ili t y coexistence1 extinct2 extinct a b Figure 6: Cascade of extinctions in the case of environmental stochasticity. The Langevin equation (25) has been numerically integrated for (a) S = S =
10 species. At the end of the simulated time span, species i is regarded as extinct if the corresponding density z i is at least 1% smaller than the maximumdensity z m = max { z , . . . , z S } . Probabilities are calculated here by averaging over 400 realizations starting from random initial conditions. Here we take κ = . ξ . To ease visualization, in panel (b) we have aggregated together the probabilities from one to four extinctions, as wellas the probabilities of observing from five to eight extinctions. Consistently with the deterministic model, the full extinction state is never observed even in thepresence of noise. also Fig. 7b). For ρ >
1, however, the three boundary equi-libria with a single, extant species are the only ones that re-main stable. As a consequence, there is a range in competition(0 . < ρ <
1) where configurations formed by a single extantspecies or by two coexisting species co-occur. Depending onthe initial condition, the dynamics can end up in one of the twoattractors. The basin of attraction of each equilibrium point willdetermine how frequently one or two species go extinct whencompetition surpasses the value ρ = .
7. This analysis is out ofthe scope of this contribution, though. Importantly, the extinc-tion sequence in non-symmetric scenarios can depend on initialconditions and is not fully determined in principle.Multiple extinctions are not precluded in general even whenthe symmetry between species is broken, as the first exampleshows. Multistability ranges could also lead to grouped ex-tinctions in deterministic scenarios. In a Lotka-Volterra com-munity model with S species there are 2 S attractors. Increas-ing complexity would likely lead to additional overlapping re-gions in competition where multiple stable attractors co-occurand species can decline together. Moreover, the determinis-tic extinction can be ambiguous. It seems difficult to estab-lish the conditions under which a general non-symmetric modelwill produce a well-defined extinction sequence. The variabil-ity introduced by idiosyncratic, species-dependent carrying ca-pacities, growth rates, or intra- and interspecific strengths, maycause the extinction sequence to be analytically unpredictablefor species-rich communities.Remarkably, the extinction sequence predicted by non-symmetric, deterministic models has nothing to do with thecascade observed when demographic stochasticity comes intoplay. We have calculated the probabilities of coexistenceand one-, two- or three-species extinctions for both stochastic parametrizations of the two non-symmetric models consideredin this subsection. To aggregate joint probabilities, we obtainednumerically the critical points by spline interpolation of the ex-act joint distribution. Now the six saddle points have asymmet-ric entries, but three of them exhibit a coordinate close to theboundary, and the remaining three saddle points present two co-ordinates near the axes. Thus, the partition of the configurationspace is conceptually equivalent to that of Fig. 3. Not surpris-ingly, as Fig. 7 evidences, the stochastic model predicts a se-quential cascade of extinctions. The threshold at which extinc-tions start occurring displaces towards smaller values of ρ , lead-ing to narrow ranges of effective stochastic coexistence. Themost likely state in the presence of stochasticity is not necessar-ily the same as in a deterministic scenario, the predictions beingutterly different in terms of the extinction sequence. Therefore,regardless of the inherent lack of symmetries that deterministicdynamics may have, demographic stochasticity can influencesignificantly the way in which extinctions take place.
4. Discussion
In this contribution, we have analyzed the extinction phe-nomenon for a symmetric, Lotka-Volterra competitive systemformed by S species. In particular, we have focused on thedifferences between the deterministic system and its stochasticcounterpart. Our main result is related to the way in whichspecies extinction proceeds: on the one hand, in the deter-ministic system, S − Competitive overlap, ρ P r obab ili t y coexistence1 extinct2 extinct3 extinct a b2 extinct Interior 1 extinct2 extinctInteriorInterior Figure 7: Non-symmetric scenario described in Sec. 3.6. In both panels, r + = . r − = . µ = .
1. (a) In this case, the deterministic model predicts amultiple extinction of two species when the threshold ρ = . K = K = K = δ = . δ = .
05. The stochastic scenario, however, yields a sequential cascade for whichextinctions are spread out along the axis of competitive overlap (symbols). (b) Here carrying capacities are uniform ( K = K = K =
15) and δ = . δ = . overlap increases. Therefore, stochasticity is responsible of aprogressive sequence of species extinction, a phenomenon thatis absent in the deterministic system. Our analyses are basedon a birth-death-immigration stochastic dynamics that was an-alyzed in deep by Haegeman and Loreau (2011) and later usedby Capit´an et al. (2015) to unveil the existence of a more restric-tive threshold in competition when demographic stochasticity isexplicitly considered. In addition to lowering the threshold incompetition at which extinctions begin to occur, ecological driftalso changes drastically the way those extinctions take place.In order to evidence this difference, we have developed con-venient analytical approximations to the critical points of thejoint probability for S = S = ρ , whereas one-species and two-species extinction most likely interchangeamong each other for intermediate values of ρ . The presenceof multiple modes in the joint probability distribution is akin tothe presence of multistability in a deterministic system. How-ever, the symmetric, deterministic model is characterized by theabsence of multiple stable states for ρ <
1. Hence the extinctioncascade described in this work is an entirely new effect causedonly by ecological drift. In addition, we have evidenced that thetransition to the deterministic model is sharp when the intensityof demographic stochasticity tends to zero. Moreover, we haveproven that environmental stochasticity also leads to a cascade of extinctions, although the ranges in competition where ex-tinctions take place are smaller and, more importantly, the fullextinction of the community is not possible in this case.In this work, we have implemented two types of stochas-ticity: demographic noise (ecological drift) and environmen-tal stochasticity. These are two typical sources of noise thatrepresent, respectively: (i) the variability in discrete popula-tion numbers as a consequence of stochastic births and deaths,or (ii) the stochastic variability in model parameters that canbe ascribed to changing environmental conditions. Althoughthey are very different implementations of noise, the main re-sult of this manuscript (the stochastic cascade of extinctions) iscommon to both of them, with some qualitative differences. Inthe case of demographic stochasticity, we do not impose a par-ticular form for the noise distribution since it directly emergesfrom the inherent stochastic dynamics of discrete populationswhose individuals undergo a number of elementary processes(in principle, the noise distribution would follow from the mas-ter equation). Other ways to implement demographic noisehave been discussed in the literature (Bonachela et al., 2012),and they would plausibly lead to mechanisms similar to thosefound here.It can be argued that the role of stochasticity in the fullysymmetric system reduces to break species symmetry and yieldto progressive species extinctions, a scenario that can arise innon-symmetric, deterministic approaches. We have illustratedwith examples that grouped extinctions are not exclusive of afully symmetric situation. Moreover, predicting the extinctionsequence in non-symmetric, deterministic cases is difficult be-cause multiple stable equilibria can co-occur in ranges of com-petition. Besides, although the extinction sequence were com-pletely determined, the cascade in the presence of stochasticity11an be totally different from that predicted by the determinis-tic model. We believe that the examples analyzed in this con-tribution show up the key role of stochasticity in communityassembly.It would be interesting to empirically test the stochastic ex-tinction cascade phenomenon. In principle, a plausible wayto conduct the experiment would involve simple protist micro-cosms where species compete for a shared resource (see, forexample, Violle et al. (2010) and references therein). Lower-ing the amount of resource could be associated to a decreasein the carrying capacity, and we have shown that the extinctioncascade mechanism is expected to arise as long as the carryingcapacity (resource availability) is reduced, see Fig. 5. For theexperiment to reproduce model settings, an individual (immi-grant) coming from the species pool should be inserted in theexperimental community at certain times. From a time serieslisting species identities at certain sampling times, one couldestimate the probabilities for coexistence and for the extinctionof one or more species, and test whether extinctions in empiri-cal systems tend to proceed sequentially or not.Our work has two important implications: first, we have de-veloped analytical approximations for conditional probabilitiesin the cases S = S =
3. As we have shown, these func-tions work well at least around the critical points of the jointprobability. Presumably, the techniques proposed here mightbe extended to approximate the joint probability distribution it-self. This approach, however, has to be performed carefully.The simplest way to approximate the three-species joint distri-bution according to our methodology is to set P ( n , n , n ) ≈ T ( n | n , n ) T ( n | n ) P ( n ), where P ( n ) is the marginal, one-species probability distribution, which can be expressed ana-lytically in terms of hypergeometric functions (Haegeman andLoreau, 2011). Apparently, the approximated distribution lacksof an important property of the exact joint distribution: it is notconserved under cyclic permutations of its arguments. There-fore, it is necessary to devise appropriate combinations of con-ditional probabilities that preserve the symmetry of the jointdistribution under cyclic permutations of its arguments. Thisresearch direction, together with a generalization for communi-ties formed by more than three species, could be worth pursuingand compared with other approximations to the joint probabil-ity, such as those developed by Haegeman and Loreau (2011).The second implication of this contribution is that we remarkthe importance of explicitly considering ecological drift in the-oretical frameworks in community ecology. Natural processesare intrinsically stochastic, because changes in population num-bers are discrete, so ecological communities are more reliablymodeled using stochastic community models, even at regimeswhere their deterministic limits are not expected to fail. In cer-tain situations, the use of deterministic dynamics in communityecology could lead to utterly different predictions, as we haveshown. Traditionally ecology has relied on these kind of mod-els, and currently considerable theoretical progress is made onthe basis of deterministic approaches. However, ecological driftmay play a determinant role that could not be captured by de-terministic formulations. We hope that the present work can in-spire new contributions in the future that highlight the distinct role of ecological drift in species community models.
5. Acknowledgements
We thank the constructive criticisms and comments of theEditor and two anonymous reviewers, who helped to improvesignificantly the original manuscript. This work was funded byproject SITES (CGL2012-39964, DA, JAC) and the Ram´on yCajal Fellowship program (DA). JAC also acknowledges partialfinancial support from the Department of Applied Mathematics(Technical University of Madrid).
Appendix A. Deterministic competitive exclusion
Gause’s competitive exclusion principle is usually stated as“two species competing for a single resource cannot coexist”.Strictly speaking, the competitive exclusion principle was firstformulated by Volterra (1926) as a mathematical proposition.Albeit a mathematical proof of the principle, based on dynam-ical systems theory, can be found in the book by Hofbauer andSigmund (1998), and is essentially the same as that of Volterra(1926) original paper, we here provide a purely algebraic alter-native demonstration (see Roughgarden (1979) as well).Assume that the densities of S species living on R resourcesvary in time according to the dynamics˙ x i = x i (cid:32) R (cid:88) j = γ i j y j − d i (cid:33) , i = , . . . , S , (A.1)where Γ = ( γ i j ) is a S × R matrix with non-negative entries, y j , j = , . . . , R , are amounts of R resources, which are alsoassumed to depend on species densities, and d i , i = , . . . , S are the rates of decline of species when all resources are zero.We now show that if S > R and species densities reach a well-defined steady state, then in the long run S − R species willgo extinct until the number of species equals the number ofresources.Imposing the condition ˙ x i = t and assuming thatall species densities are positive, we get the non-homogeneouslinear system Γ y = d , with y = ( y i ) ∈ R R , and d = ( d i ) ∈ R S .Given that matrix Γ has S rows and only R < S columns, thesystem will be incompatible except for some specific vectors d belonging to the image Im( Γ ) of matrix Γ . As a consequence,to find an equilibrium solution some species densities must goto zero. For those displaced species, the corresponding rows ofmatrix Γ can be removed until rank( Γ ) rows remain. If the rankequals the number R of resources, the system turns out to becompatible and determinate, and R species will stably coexist.Note that it is the rank of matrix Γ rather than the number ofresources itself that induces competitive exclusion.Therefore, competition for shared resources imposes a limitto the maximum number of species that can stably coexist.This result has a counterpart in the stability of Lotka-Volterraequations derived from the MacArthur’s consumer-resourcemodel (MacArthur, 1970; Chesson, 1990), i.e., the particularcase of Eq. (A.1) in which per-capita resource growth rates ˙ y i / y i x i = r i x i (cid:32) − x i + (cid:80) j (cid:44) i ρ i j x j K i (cid:33) , i = , . . . , S . (A.2)Here r i is interpreted as the intrinsic growth rate of species i , K i as a carrying capacity, and ρ i j measures interspecific compe-tition strength between species i and j relative to intraspecificcompetition. As shown by Chesson (1990), if S = R = rank( Γ ),then exists a unique solution of the system x i + (cid:80) j (cid:44) i ρ i j x j = K i , i = , . . . , S . The resulting equilibrium point will be interior ifevery species density remains strictly positive. Moreover, Ches-son (1990) demonstrated that if there exists a unique interiorequilibrium point for the dynamics (A.2), it will be globallystable if and only if S = R = rank( Γ ). Therefore, if competitiveexclusion does not operate, a unique globally stable equilibriumpoint is reached and, conversely, if the Lotka-Volterra equationspresent an interior, globally stable equilibrium point, any pos-itive initial condition will make the dynamics (A.2) convergeto the equilibrium point. This automatically ensures that noneof the S species is driven to extinction by competitive exclu-sion (MacArthur, 1970; Takeuchi, 1996). Appendix B. Stability of the symmetric, deterministicmodel for ρ ≥ In this Appendix we perform a stability analysis of the equi-librium points of the symmetric, deterministic dynamics in thecompetition regime ρ ≥
1. As we mentioned before, when ρ < ρ =
1. In this case, thesystem is also stable and the initial condition determines theattractor which the dynamics converges to. Any equilibriumpoint is such that its densities satisfy S (cid:88) j = x j = K . (B.1)Let J ( t ) = (cid:80) Sj = x j ( t ). Then, by summing up the equations ofthe system (1) for ρ = dJdt = rJ (cid:18) − JK (cid:19) , (B.2)which can be integrated and yields J ( t ) = (cid:34) K + (cid:32) J − K (cid:33) e − rt (cid:35) − , (B.3) J being the initial condition for J ( t ), J = J (0) = (cid:80) Sj = x j (0).Therefore, the dynamics of each species turns out to be decou-pled, dx i dt = rx i (cid:32) − J ( t ) K (cid:33) . (B.4) The equilibrium point to which (B.4) converges is determinedby the initial condition vector x = ( x (0) , . . . , x S (0)). For twodistinct species i and j , (B.4) implies that ˙ x i / ˙ x j = x i / x j . In-tegration yields x i ( t ) / x j ( t ) = x i (0) / x j (0), which means that theproportions of population densities are conserved along the dy-namics, whose orbits are reduced to straight lines starting fromthe initial condition x along the direction determined by thevector x . Therefore, the final equilibrium point is given by theintersection of the hyperplane (cid:80) Sj = x j = K and the line thatlinks the initial condition point, x , and the origin. Any of these(infinite) equilibrium points will be stable, provided that the ini-tial densities satisfy x i (0) ≥ i .The case ρ > S equilibrium points with posi-tive or zero densities. Without loss of generality, for 0 ≤ n ≤ S any equilibrium point will be of the form x n = (cid:32) ˆ x n u S − n n (cid:33) , (B.5)where subscript n indicates that n densities are strictly equalto 0, and u n = (1 , . . . , T is a vector with n entries equal to1. Any equilibrium with n zero densities can be written as apermutation of (B.5), without altering the subsequent stabilityanalysis. In addition, the non-zero entries of x n are the solutionsof the system(1 − ρ ) x i + ρ S − n (cid:88) j = x j = K , i = , . . . , S − n . (B.6)This system admits a single solution for which the S − n specieshave equal densities, x i = ˆ x n , whereˆ x n = K − ρ + ρ ( S − n ) . (B.7)Our demonstration reduces to show that, if ρ >
1, all the non-trivial equilibria act as repellors except for n =
1, i.e., whenonly a single species survives. To this purpose we evaluate theeigenvalue spectra of the Jacobian matrix of the system (1). TheJacobian matrix J can be expressed in a block form as J = (cid:32) J J J J (cid:33) , (B.8)where J = − r ˆ x n K (cid:104) (1 − ρ ) I S − n + ρ u S − n u T S − n (cid:105) , (B.9) J = − r ρ ˆ x n K u S − n u T n , (B.10) J = n T S − n , (B.11) J = r (1 − ρ ) ˆ x n K I n , (B.12)and I n is the n × n identity matrix and n = (0 , . . . , T is the zerovector with n entries. Without loss of generality, eigenvectorscan be written as v = (cid:32) a S − n b n (cid:33) , (B.13)13here a S − n and b n are column vectors with dimensions S − n and n , respectively. First let us assume that b n = n . Thespectral problem reduces to J a S − n = λ a S − n . (B.14)Two different solutions arise: if u T S − n a S − n =
0, i.e., a S − n is or-thogonal to u S − n , then the eigenvalue is λ = − r (1 − ρ ) ˆ x n K (B.15)with algebraic multiplicity S − n −
1. However, if a S − n ∝ u S − n we find λ = − r as an eigenvalue with algebraic multiplicityequal to 1.On the other hand, if b n (cid:44) n we have to solve J a S − n + J b n = λ a S − n , J b n = λ b n . (B.16)Since J is proportional to I n , we find λ = r (1 − ρ ) ˆ x n K , (B.17)with n eigenvectors b n = e ( i ) n ( i = , . . . , n ), e ( i ) n being the i -thvector of the canonical basis of R n . The algebraic multiplicityassociated to (B.17) is equal to n . Substituting these results intothe first equation of (B.16) yields (cid:104) (1 − ρ ) I S − n + ρ u S − n u T S − n (cid:105) a S − n + ρ u S − n = − (1 − ρ ) a S − n . (B.18)Given the structure of the matrices involved, we look for a so-lution of the form a S − n = α u S − n , which implies a non-trivialvalue α = − ρ/ [2(1 − ρ ) + ρ S ].Two eigenvalues determine the asymptotic stability of all theequilibrium points: λ = − r (1 − ρ ) ˆ x n / K and λ = − λ = r (1 − ρ ) ˆ x n / K (the third eigenvalue, λ = − r , is always negative).If ρ > ≤ n < S − λ > S − n − >
0. Therefore, anyequilibrium with less than S − n = S − λ is notan eigenvalue anymore and the two other eigenvalues remain: λ = r (1 − ρ ) ˆ x n / K < S −
1) and λ = − r < S boundary equilibria witha single extant species are asymptotically stable.Since the trivial equilibrium point (associated to complete ex-tinction) is obviously unstable, we deduce that any orbit start-ing from an interior initial condition will be repelled if it getsclose to any equilibrium, except when the equilibrium point isformed by a single extant species. If the orbit enters the basinof attraction of any of those S equilibria, it will end up in itasymptotically. This implies the extinction at a time of S − ρ > Appendix C. Numerical calculation of the steady-stateprobability distribution
In order to compute numerically the stationary joint dis-tribution, we limit the infinite configuration space to the set Ξ ≡ { , , . . . , n max } S by choosing n max large enough so that theprobability of finding a population number equal to n max is neg-ligible (we choose n max as the integer part of 2 K , which fulfillsthe requirement).The stationary distribution is obtained by solving the embed-ded Markov chain associated to the continuous-time Markovprocess (Karlin and Taylor, 1975). The transition matrix of theembedded Markov chain is defined by the transition probabili-ties Pr { n → n ± e i } = q ± i ( n ) Λ ( n ) , (C.1)where Λ ( n ) = (cid:80) Si = [ q + i ( n ) + q − i ( n )]. The remaining transi-tions have zero probability. Elementary events (overall birthsand deaths) take place after exponential times, so that the timelapsed to the next event is drawn from a random variable τ withcumulative distribution Pr( τ ≤ t ) = − e − Λ ( n ) t . Once the steady-state distribution ϕ = ( ϕ ( n )) of the embedded Markov chain —i.e., the left-eigenvector of the transition matrix with eigenvalue1— has been determined, according to the mean time spent bythe process at state n , the probability of finding the continuous-time Markov process at state n turns out to be (Cinlar, 1975) P ( n ) = ϕ ( n ) Λ ( n ) − (cid:80) m ∈ Ξ ϕ ( m ) Λ ( m ) − . (C.2) Appendix D. Coefficient of variation of population abun-dances in the limit of large carrying capacity
In this section we derive an analytical expression for thecoefficient of variation of population abundances in the smallstochasticity limit. In the case of demographic stochasticity,low variability levels can be obtained for large population sizesor, equivalently, in the limit of large carrying capacity. We buildon the Gaussian approximation for the joint probability distri-bution, which is valid in the limit K (cid:29) S .The Gaussian approximation can be obtained as the solutionto the Fokker-Planck equation deduced from the master equa-tion (2). We do not reproduce its derivation here; it can befound at the Supplemental Information of Capit´an et al. (2015).Under this approximation, the joint probability distribution isexpressed as Π ( n ) = Z exp {− ( n − ˆ x u S ) T Q ( n − ˆ x u S ) } , (D.1)where Z is an appropriate normalization factor,ˆ x = K − ρ + ρ S ) + (cid:114) + µ (1 − ρ + ρ S ) rK , (D.2)and the covariance matrix Q − is given by Q − = b a (cid:18) I S − ca + cS u S u T S (cid:19) . (D.3)In terms of model parameters, a = r ˆ x (1 − ρ ) / K + µ/ ˆ x , b = r + ˆ x + µ ), and c = r ρ ˆ x / K . In the large carrying capacity limit,14he average population abundance is expressed through a seriesexpansion on powers of K as (cid:104) n (cid:105) = ˆ x = K − ρ + ρ S + µ r + O ( K − ) . (D.4)Similarly, series expansions give a = r (1 − ρ )1 − ρ + ρ S + O ( K − ) , b = r + K − ρ + ρ S + µ (cid:32) r + r + (cid:33) + O ( K − ) , c = r ρ − ρ + ρ S + O ( K − ) . (D.5)Inserting these expressions into (D.3) yields, up to order K , theapproximation Q − = r + Kr (1 − ρ ) (cid:32) I S − ρ − ρ + ρ S u S u T S (cid:33) + O ( K ) . (D.6)Therefore the standard deviation of population abundance, σ n ,can be obtained as the square root of diagonal elements of ma-trix Q − , σ n = (cid:115) r + (1 − ρ + ρ S ) r (1 − ρ )(1 − ρ + ρ S ) K / + O ( K − / ) , (D.7)which (not surprisingly) scales with K as K / . Finally, thecoefficient of variation of population abundances is expressed,in the limit K (cid:29)
1, as ν = σ n (cid:104) n (cid:105) = (cid:115) r + (1 − ρ + ρ S )(1 − ρ + ρ S ) r (1 − ρ ) K − / . (D.8)Strictly speaking, the deterministic scenario ( ν =
0) is onlyachieved in the limit K → ∞ . However, low stochasticityregimes can be assessed using Eq. (D.8): if the actual coeffi-cient of variation is close to that yielded by the Gaussian ap-proximation, both of which are small for large K , then extinc-tion configurations are precluded and the variability of popula-tions with respect to the mean value is small. We adopt, as apractical definition for low stochasticity, the parameter combi-nations for which the actual coefficient of variation is close tothe approximation given by Eq. (D.8). Appendix E. Stability analysis for two deterministic mod-els with non-symmetric competition
Here we analyze the stability of the equilibrium points of thetwo non-symmetric, three-species competitive dynamics of theform (26) considered in the main text, for which the interactionmatrix is written as R = ( ρ i j ) = ρ ρ + δ ρ ρ + δ ρ + δ ρ + δ , (E.1) δ and δ being two positive numbers. The sign of equilibriumdensities and the corresponding eigenvalues of the Jacobian ma-trix determine the ranges of ρ for which the system is stable. Weimpose the condition ρ ≥ r >
0, in both cases thefull extinction equilibrium ˆx = (0 , ,
0) is unstable.In the first example, δ = . δ = .
05. Species carryingcapacities are non-uniform: K = K =
16, and K = ˆx = R − K , with K = ( K , K , K ) T , has three positive densities andis asymptotically stable if and only if 0 ≤ ρ < . λ of the stability (Jacobian) matrix are expressed as λ = λ (cid:48) r , i.e., they are scaled by the growth rate r > ˆx = (cid:32) − ρ )1 − ρ , − ρ )1 − ρ , (cid:33) , (E.2)are positive if and only if 0 ≤ ρ < . ρ > .
5. Thescaled eigenvalues λ (cid:48) are given by λ (cid:48) ∈ (cid:40) − , − ρ + ρ − ρ ) , − (5 − ρ )(2 − ρ )10(1 − ρ ) (cid:41) . (E.3)The condition λ (cid:48) i < i = , , . < ρ < .
5. Therefore, this point is never feasible and stable at thesame time.(b) For ˆx = (cid:32) − ρ )99 − ρ − ρ , , − ρ )99 − ρ − ρ (cid:33) , (E.4)we require positivity for the first and third entries, whichyields 0 ≤ ρ < . ρ > .
9. On the other hand, theeigenvalues in this case are λ (cid:48) ∈ (cid:40) − , − ρ − ρ (9 − ρ )(11 + ρ ) , − − ρ )(19 − ρ )(9 − ρ )(11 + ρ ) (cid:41) . (E.5)Stability implies 0 . < ρ < .
9, which is incompatiblewith the feasibility condition.(c) The third equilibrium with a single extinct species is ˆx = (cid:32) , − ρ )399 − ρ − ρ , − ρ )399 − ρ − ρ (cid:33) . (E.6)This point is feasible if and only if 0 ≤ ρ < .
75 or ρ > . λ (cid:48) ∈ (cid:40) − , − − ρ )(6 − ρ )(19 − ρ )(21 + ρ ) , − ρ − ρ (19 − ρ )(21 + ρ ) (cid:41) . (E.7)Since the system of inequalities λ (cid:48) i < i = , ,
3) turnsout to be incompatible, this point is unstable for all valuesof ρ .15d) ˆx = (40 , , {− , (2 − ρ ) / , − ρ ) / } . This point is stable for ρ > . ˆx = (0 , , λ (cid:48) ∈ {− , (5 − ρ ) / , − ρ ) / } . Theequilibrium point is stable if and only if ρ > . ˆx = (0 , , λ (cid:48) ∈ {− , (19 − ρ ) / , − ρ ) / } .Asymptotic stability is achieved for ρ > . ≤ ρ < .
4, the only stable point isthe coexistence equilibrium. However, for ρ > .
4, only two-extinct species equilibria remain asymptotically stable. Sincethe eigenvalues are continuous functions of model parameters,close to this example we can find multiple non-symmetric sys-tems that exhibit a grouped, two-species extinction as ρ in-creases.The second example shows that multiple stable equilibria canco-occur when interactions are chosen non-symmetrically. Inthis case, we have taken δ = . δ = . K = K = K =
15. The densities of the interior equilibrium point areexpressed as x = − ρ )(7 − ρ ) / D ( ρ ) , x = − ρ )(11 − ρ ) / D ( ρ ) , x = − ρ )(1 − ρ ) / D ( ρ ) , (E.8)where D ( ρ ) = − ρ − ρ + ρ . The feasibilityanalysis of the equilibrium point yields, for ρ ≥
0, the ranges0 ≤ ρ < . . < ρ < . ρ > .
1. The eigenvaluesof the stability matrix can be fully calculated, although theirexpressions are too cumbersome to be reproduced here. It iseasy to check that the three eigenvalues are negative if and onlyif 0 ≤ ρ < .
3. Therefore, this equilibrium point is interior andasymptotically stable if and only if 0 ≤ ρ < . ˆx = (cid:16) + ρ , + ρ , (cid:17) : the scaled eigenvalues λ (cid:48) are (cid:110) − , − ρ + ρ ) , − + ρ + ρ (cid:111) . Therefore, this point is stable for 0 . <ρ < ˆx = (cid:16) + ρ , , + ρ (cid:17) : λ (cid:48) ∈ (cid:110) − , − ρ + ρ ) , − + ρ + ρ (cid:111) . The stabilityconditions form an unfeasible problem, so this point turnsout to be unstable for any ρ .(c) ˆx = (cid:16) , + ρ , + ρ (cid:17) : λ (cid:48) ∈ (cid:110) − , − ρ + ρ , − + ρ + ρ (cid:111) . Again,this point is unstable for all values of competitive overlap.(d) ˆx = (15 , , (cid:110) − , − ρ, − ρ (cid:111) . This point is stable for ρ > ˆx = (0 , , λ (cid:48) ∈ (cid:110) − , − ρ , − ρ (cid:111) . The equilibriumpoint is stable if and only if ρ > ˆx = (0 , , λ (cid:48) ∈ (cid:110) − , − ρ , − ρ (cid:111) . Stability is attainedfor ρ > .