Alternative steady states in random ecological networks
AAlternative steady states in random ecological networks
Yael Fried, Nadav M. Shnerb, and David A. Kessler
Department of Physics, Bar-Ilan University, Ramat-Gan IL52900, Israel
In many natural situations one observes a local system with many competing species which is coupledby weak immigration to a regional species pool. The dynamics of such a system is dominated byits stable and uninvadable (SU) states. When the competition matrix is random, the number ofSUs depends on the average value of its entries and the variance. Here we consider the problem inthe limit of weak competition and large variance. Using a yes/no interaction model, we show thatthe number of SUs corresponds to the number of maximum cliques in a network close to its fullyconnected limit. The number of SUs grows exponentially with the number of species in this limit,unless the network is completely asymmetric. In the asymmetric limit the number of SUs is O (1).Numerical simulations suggest that these results are valid for models with continuous distributionof competition terms. I. INTRODUCTION
The richness of ecological communities poses a pro-longed theoretical challenge. Focusing on guilds of manyspecies competing for a common resource (and neglect-ing, for the moment, processes like predation or mutual-ism) the main problems are two. First, the competitiveexclusion principle [1, 2] suggests that the result of com-petition for a single limiting resource is the extinctionof all except the fittest speices, and that in the presenceof a few limiting resources the equilibrium number ofspecies is smaller than or equal to the number of resources[3]. Second, even if the number of limiting resources islarge, May [4] pointed out that if the niche overlap be-tween species is substantial the chance of a system of N species to admit a stable equilibrium decreases expo-nentially with N . May’s result is based on the spectralproperties of random stability matrices. Practically, itimplies that to achieve a stable coexistence of more than6-8 species one has to fine-tune the competition parame-ters in an unrealistic manner.However, in many situations the (inter and intra) spe-cific dynamics takes place on local patches, which arecoupled by migration to each other or to a regional pool.Accordingly, many ecological models are focused on thedynamics of a local patch, putting aside the global sta-bility problem. A mainland-island model [5, 6] is thesimplest scenario considered in this context: a set of lo-cal populations of different species are competing witheach other and the island is exposed to weak migrationfrom a static pool of N species on the mainland. Thestructure of the community on the island reflects a bal-ance between local extinctions and colonization by im-migrants from the mainland. Extinctions may be eitherdeterministic, due to the pressure that a species suffersfrom its competitors (or from the local environment), orstochastic, caused by the random nature of the birth-death process [7, 8] possibly superimposed on the effectof environmental variations [9, 10].In a recent work [11], Kessler and Shnerb suggested aclassification of the qualitative features of the communityon the island. Four different “phases” were identified. I. Full coexistence : If the interspecific competitionis weak (say, if different species use essentially dif-ferent resources) any species in the mainland mayinvade the island and establish a finite population,so all N species are present on the island. Localextinctions still occur but if the local populationsin steady state are large, these events are rare andtransient. Technically speaking, the deterministic(i.e., noise-free) model allows for a stable fixed pointwith all the N species coexisting.II. Partial coexistence : As the competition amongspecies grows, the species’ abundances decay as theyfeel more pressure from other species. Since thecompetition matrix is heterogenous, some speciesfeel more pressure than others, and the determinis-tic model eventually allows for a stable fixed pointfor a finite subset containing S of the N species. Theother N − S species on the mainland cannot invadethe island, i.e., their growth rate at low densities onthe island is negative.III. Disordered : When the competition increases evenfurther, and the competition matrix is not symmet-ric (meaning that species 1 may put a lot of pressureon species 2 but species 2 puts much less pressureon species 1, say), the system may not have an at-tractive fixed point at all or, even if it have one,its basin of attraction will be very narrow. In thepresence of noise, the system fails to converge toan equilibrium state and instead it shows intermit-tent behavior with many long excursions that reflecthigh-dimensional chaotic/periodic trajectories.IV.
Alternative steady states : Finally, when thecompetition terms are large, there will be a num-ber of different subsets of the N species which areboth stable and uninvadable. For example, if theinterspecific competition is extremely large and theisland is colonized by a single species, all other main-land species cannot invade, so one needs to wait fora rare stochastic extinction in order to see speciesturnover. If the competition is not so strong, subsetsof more than one species play the same role: species a r X i v : . [ q - b i o . P E ] M a r within the subset interact only weakly so they maylive happily together, but other species cannot in-vade.The aim of this paper is to discuss this last phase,which is characterized by strong competition and al-ternative steady states. The immediate motivation forthis discussion comes from a recent paper by Fisher andMehta [12], who suggested that in this phase the dynam-ics of the island exhibits a glass transition: with weaknoise/immigration the system is trapped for most of thetime in one of the SUs, while when stochastic effects arestrong it behaves like a “liquid” and its dynamics is closerto the disordered phase discussed above. In [12], a versionof the symmetric competition model with strong interac-tion was mapped into a well-known physical model forglassy behavior, the random energy model [13].Technically, the appearance of a glass transition in therandom energy model is related to the exponential in-crease of the number of local minima with the systemsize (here, the number of mainland species, N ). Whenboth the energy and the entropy increase linearly withthe system size, a glass transition appears at finite tem-perature (level of noise). Therefore, it is natural to inves-tigate the scaling of the number of SUs on an island withthe number of species. In fact, this problem is consideredby ecologists for many years [14, 15].Recently, we have studied this problem and derived afew exact results [16]. Using a version of the model wecall the Binomial model, we have mapped the problemof counting SUs to that of finding the number of max-imal cliques in a random network. We showed that ina particular parameter regime the number of SUs is not exponential; it goes like N ln( N ) for symmetric networks,and like N/ ln / ( N ) for (fully) asymmetric networks.Here we are going to analyze the very same model in adifferent parameter regime, which includes the case wherethe competition is weak and the heterogeneity is strong(this is the case considered in [12] and [17], see discus-sion). We will show that in this regime the number ofSUs indeed increases exponentially with N , if the systemis not fully asymmetric. On the other hand, for the asym-metric system (see the definitions below) the number ofSUs in this regime is order one. We also show how tomake a connection between this weak competition regimeand the results obtained for strong competition in [16],and provide some intuitive argument.This paper is organized as follows. In the next sectionwe summarized the results of [16]; we present the general-ized competitive Lotka-Volterra model (GCLV) and oursimplified, Binomial model, and show how to map SUsto maximal cliques, arriving at the formula of Bollob´asand Erd¨os [18]. In the next two sections we present ourmain result, the number of SUs, as calculated from thisformula, for the symmetric and the asymmetric case. Wealso present numerical computations showing that the re-sults of the Binomial model also describe the qualitativebehavior of the more realistic Gamma model. Finally we discuss the works of Refs. [12] and [17] and the relevanceof the results presented herein and in these papers torealistic ecological networks. II. THE MODEL
To get oriented, we start with a system of two compet-ing species without noise and immigration. The GCLVreads dx dt = x − x ( x + ˜ c , x ) (1) dx dt = x − x ( x + ˜ c , x ) , where x i is the abundance of each of the species. Thissystem is characterized by the competition matrix (cid:18) c , ˜ c , (cid:19) . where the intraspecific density dependence (a decreasein the growth rate with abundance, manifested in thediagonal term) was taken to be one and is not part of thematrix. The stress put upon species 1 by species 2 is ˜ c , and the stress put upon 2 by 1 is ˜ c , . ρ ≡ ˜ c , + ˜ c , isa rough measure for the niche overlap, or total strengthof competition in the system. κ ≡ ˜ c , − ˜ c , measuresthe heterogeneity of the competition matrix, i.e., it tellsus how much the species differ from each other in theirresponse to an increase in the density of a competitor. Weconsider a model as symmetric if ˜ c i,j = ˜ c j,i for any pairof species, and as asymmetric if there is no correlationbetween ˜ c i,j and ˜ c j,i .A steady solution for (1) in which both x and x arenon-negative is called a feasible solution (we cannot allownegative densities). If both densities are positive andthe solution is stable, we called it a coexistence solution.Such a solution for (1) exists as long as ρ < − | κ | ,meaning that, for a given level of niche overlap ρ thesystem becomes less stable as the heterogeneity grows.This basic logic holds also in more diverse systems.For a system of many competing species the GCLV is: dx i dt = x i − x i − (cid:88) j (cid:54) = i ˜ c i,j x i x j . (2)The mean of the terms of the competition matrix, C ≡ N ( N − (cid:88) i,j ˜ c i,j , reflects the overall strength of the competition in the sys-tem. The variance of these entries, ˜ σ , is the simplestmeasure for its heterogeneity. To emphasize these prop-erties we will factor out the average from the competitionmatrix, so the GCLV takes the form, dx i dt = x i − x i x i + C N (cid:88) j (cid:54) = i c i,j x j , (3)where c i,j = 1.May’s analysis [4] of the complexity-stability problemis based on the observation that a linear stability fea-tures of a feasible solution of (3) yields an N × N ran-dom matrix which is similar to the interaction matrix.For such a state to be stable all the eigenvalues of thismatrix should be negative. However, a random matrixwith ( −
1) on the diagonal and off diagonal terms withmean zero and variance C ˜ σ has its eigenvalues between − C ˜ σ √ N and − − C ˜ σ √ N , so a feasible solution for (3)is almost surely unstable when N , the number of species,is above N c = 1 / ( Cσ ) . The applicability of this argu-ment to purely competitive systems requires some morediscussion, since the main problem in these systems is toensure feasibility [11, 19], but the main insight turns outto be valid here as well.In this paper, as in [16], we are interested in the fea-tures of the system way above this “May limit”, i.e., when N (cid:29) N c and the system supports alternative steadystates. What we would like to know is how many sta-ble and uninvadable (SU) subsets of the N species exist,i.e., how many S -subsets of the N species satisfy the fol-lowing two conditions:1. Stability and feasibility:
Eq. (3), when limited to aspecific size S subset, S , yields a time independentsolution for which ¯ x i > S , where ¯ x i is the equilibrium density of the i -thspecies in the subcommunity.2. Uninvadability:
Eq. (3), when applied to all absent N − S species and linearized around the fixed point x i = ¯ x i for i ∈ S and x i = 0 for i (cid:54)∈ S , yieldsnegative growth rates ˙ x i /x i for all i (cid:54)∈ S .We are interested in the SU enumeration problem for arandom matrix, so we would like to draw the c i,j s from auniform, positive semi-definite, distribution with a meanone and a given variance [11, 12]. For our numerics wehave used the Gamma distribution for this purpose, anddenote this as the Gamma model. A c i,j matrix (forsimplicity the examples are given for the symmetric case)may look like, .
95 1 .
63 0 . .
95 0 0 .
48 0 . .
63 0 .
48 0 1 . .
96 0 .
97 1 .
12 0 . . To map this model to the maximum clique problem, wetreat an alternative model, the Binomial (yes/no) model,where all the elements of the c i,j matrix (in the asymmet-ric case; the pair c i,j = c j,i in the symmetric case) eitherare strictly zero (with probability p ) or (with probabil-ity 1 − p ) are equal to a finite constant C · A , so theinteraction matrix ˜ c i,j takes the form, say, C A AA A A A . The Gamma and the Binomial model have the samecompetition strength, C , if A = 11 − p . (4)The variance of the matrix elements of the Binomialmodel is given by, ˜ σ = C p − p . (5)For the symmetric model, if C is large enough, eachpair of species i and j is either non-interfering, ˜ c i,j =˜ c j,i = 0, or mutually exclusive. Accordingly, as explainedin [16], the SU problem has a geometrical interpretation.For a graph in which each node represent a species andeach pair of non-interfering species is connected by anedge, a stable state corresponds to a subset S of nodessuch that the induced subgraph is complete. For thisstable state to be uninvadable any vertex that is not apart of the clique is required to have at least one mutu-ally exclusive species in the clique, i.e., that the cliqueis maximal such that it cannot be extended by includingany other connected vertex. Accordingly, for large C thenumber of SUs is equal to the number of maximal cliquesof the corresponding graph.In [16] we showed that, as long as p is O (1), the growthof the number of maximal cliques, SU ( N ) with N is not exponential, and in fact for a symmetric system it growsas SU ( N ) ∼ N ζ ( p ) ln( N ) , (6)where ζ ( p ) = 1 / [2 ln(1 /p )] . Clearly, the expression (6)must fail when the value of p is close to one, since p → /p ) → SU ( N ) → ∞ . On the otherhand when p = 1 we reach an extreme stabilizationand the system has only one stable uninvadable state, SU ( N ) = 1. In order to clarify the behavior of the systemin this limit, in the next two sections we will find a for-mula for SU ( N ) in the limit N → ∞ , where p = 1 − α/N and α = O (1). III. THE SYMMETRIC CASE
In this section we consider the symmetric version of thebinomial model. Every pair of species is noninterfering( c i,j = c j,i = 0) with probability p and have symmetriccompetition ( c i,j = c j,i = A ) with probability 1 − p . Weassume that C is large such that no pair of competingspecies is allowed on the island (if they are interacting,then they are mutually exclusive), and a species cannotinvade the island in the presence of one of its competi-tors. Accordingly, each maximal clique of the network isa stable and uninvadable state.To get the basic intuition for the results we derive inthis section, let us consider the case p = 1, i.e., all speciesare noninteracting. Clearly, in this case there is only onemaximal clique - the one with all the N species.Now let us break the link between, say, species 1 and2, so c , = c , = A . The number of maximal cliques isnow two: species 1 and all the species 3 ..N and the set2 ..N . Breaking the next link (without loss of generality,between 3 and 4) doubles the number of maximal cliquesand so on. Hence, the number of maximal cliques growsexponentially as p decreases, until we start to break morethan one link per node. Since there are O ( N ) links inthe system, this will happen when the number of brokenlinks is O ( N ). Accordingly, one expects that, when thedeviation of p from one is O (1 /N ), the number of SUswill be exponentially large in N .Bollob´as & Erd¨os [18] showed that the number of max-imal cliques of size S in a random graph of N nodes is given by SU ( N, S ) = (cid:18) NS (cid:19) p S ( S − / (1 − p S ) N − S . (7)In [16], we performed the sum over S , giving the behaviorof SU ( N ), Eq. (6), when p is not too close to unity.Now let us find the leading asymptotic behavior of thesum (7) when 1 − p is small, of order O (1 /N ). Moreprecisely, we define p ≡ − αN , S ≡ N β, (8)where α and β are both O (1). With these definitions, (7)reads SU ( N, β ) = (cid:18)
NN β (cid:19) (cid:16) − αN (cid:17) Nβ ( Nβ − / (cid:16) − (1 − αN ) Nβ (cid:17) N − Nβ . (9)In the large N limit, the expression (9) may be writtenas, SU ( N, β ) ≈ (cid:18) NN β (cid:19) e − Nαβ / e N (1 − β ) ln(1 − e − αβ ) (10)Since for β ∼ O (1), both N β and N (1 − β )are large, wecan approximate the combinatorial factor using Stirling’sformula, giving (cid:18) NN β (cid:19) ≈ e − N [(1 − β ) ln(1 − β ) − β ln( β )] (cid:112) πN β (1 − β ) . (11)Accordingly, SU ( N, β ) ≈ e Nf ( β ) (cid:112) πN β (1 − β , (12)where F ( β ) = − (1 − β ) ln(1 − β ) − β ln( β ) (13) − αβ − β ) ln(1 − e − αβ ) . The total number of SUs is the sum over S of SU ( N, S ), which is translated to an integral over β of(12). This integral may be approximated via Laplace’smethod, as F ( β ) has a maximum in the range 0 < β < β , which de-pends on α and satisfies F (cid:48) ( β ) = − αβ + α (1 − β ) e αβ − − β ) − ln( β ) − ln(1 − e − αβ ) = 0 . The graph of β ( α ) is depicted in Fig. 1. Then, toleading order, the total number of SUs, SU ( N ), is givento leading order by SU ( N ) ≡ N (cid:90) SU ( N, β ) dβ = e NF ( β ) N (cid:112) | F (cid:48)(cid:48) ( β ) | β (1 − β ) , (15)so that indeed the number of SUs increasing exponen-tially with N in this parameter range. Here we will havean interest only in the controlling factor, so we focus on F ( β ). For general α , this needs to computed numer-ically, with the results shown in Fig. 2. We see that F ( β ) rises from 0 at α = 0, reaches a maximum andthen decreases slowly for large α . The behavior at largeand small α is accessible to analysis. For small α , we seefrom Fig. 1 that β is close to unity and thus to leadingorder in α , 1 − β we have F (cid:48) ( β ) ≈ ln(1 − β ) − ln α (16)so that β ≈ − α and F ( β ) ≈ α/
2. For large α , since β is small but αβ is large, the dominant balance of termsfor large α is − αβ + αe αβ − ≈ ln( β ) . (17)We can exactly solve this equation using an auxiliaryvariable r , α ≡ r r ln r , in terms of which β = r − r , ascan be directly verified by substituting into the equation.This implicit approximate solution is correct to order1 /α for large α . If we try to produce an explicit solutionfrom this, we run into correction terms like log(log( α ))and the convergence is super-slow. Nevertheless, a simplerough approximation is β ≈ ln αα , (18) FIG. 1. Solutions for Eq. (14). The full, black line is theexact numerical solution, the dashed blue is 1 − α and thedotted red line depicts ln( α ) /α .FIG. 2. F ( β ( α )), or ln( SU ) /N , vs. α for 0 < α < N , where thecoefficient of the exponent is between 0 and 0 . up to corrections of ln(ln r ) /α , which is correct to betterthan 6% for α >
10. We can now approximate F ( β ) forlarge α , where the αβ / F ( β ) ≈ r − r ln r ≈ ln α α (19)Our result connects directly with our previous result,Eq. (6), when α is O (1). Remember the relation between p and α , as α becomes large, p moves away from theregion close to unity, and soln SU ( N ) ≈ ln N (1 − p )2(1 − p ) ≈ ln N /p ) (20)as expected.In the opposite limit, when α is very small, say, α = γ/N (so p = 1 − γ/N ), β = 1 − α , and the exponen-tial term of (15), exp( γ/ γ = 0 and growsexponentially with γ , as discussed above. ln ( p ) l n ( nu m be r o f S U s ) FIG. 3. ln( SU ) vs. ln( p ), as obtained from a numerical sum-mation of the Bollob´as-Erd¨os formula for the symmetric case,for N = 20 (red) N = 40 (green) and N = 80 (blue). Thegrowth of the number of SUs with N becomes exponentialwhen 1 − p is O (1 /N ), as expected. When 1 − p is O (1 /N ),there is a drop towards one SU with all the N species. In thisregime the number of SU s is independent of N . These three regimes are depicted in Figure 3. As op-posed to the case where p is O (1) (or α is O ( N )) consid-ered in [16], where the growth is N ln N type, when α is O (1) the growth is exponential while if α is O (1 /N ) thenumber of SUs is close to one. In figure 4 we show thatthe same qualitative behavior is observed in the corre-sponding symmetric Gamma model [16]. THE ASYMMETRIC NETWORK
Unlike the symmetric case where c i,j = c j,i , in anasymmetric system c i,j and c j,i are drawn independentlyfrom a given distribution. In this section we consider theBinomial model in this case.The strong competition phase of the asymmetric Bino-mial interaction model allows for three types of relation-ships between species i and j . As in the symmetric case,it may happen that c i,j = c j,i = 0, so the species are non-interfering, and c i,j = c j,i = A , meaning that for largeenough A the two species are mutually exclusive. Theasymmetry allows for a third, dominance , relationship atlarge A : if c i,j = A, c j,i = 0, species i may invade j butthe opposite process is forbidden. Accordingly, j may bea member of a maximal clique only if another species inthis clique is uninvadable by i .In the asymmetric Binomial model, we define ˜ p to bethe chance that a single entry of the interaction ma-trix is zero. The argument presented at the beginning ofthe last section here yields a completely different answer.Starting from ˜ p = 1 and breaking one link (say, between1 and 2), implies that 1 can invade 2 but 2 cannot invade1, so the number of maximum cliques remains one. Thiswill be the case until we hit both links between two spe-cific species, and again this will happen only when 1 − ˜ p is O ( N ). Accordingly, at the asymmetric case we expectthat the number of SUs for p close to one will be O (1).By extending this argument, one can develop some in- FIG. 4. Number of SUs (averaged over 500 samples) vs. σ γ as obtained from an exact enumeration of all the stable anduninvadable combinations of species in the Gamma model for N = 20. Red points are actual results, each obtained from ex-amination for stability and uninvadability of all the 2 com-binations of species, the dashed line is just to guide the eye.The Gamma model is described by Eq. (3) with C = 1 andwhere each pair of numbers c i,j = c j,i is picked independentlyfrom a Gamma distribution with mean one and standard de-viation σ γ . While the number of SUs is smaller than theirnumber in the corresponding Binomial model (we have sug-gested in [16] that the Binomial model yields an upper boundfor the Gamma) we still observe the growth in the number ofSUs when σ γ ∼ √ N (around 4-5), then it drops towards thefull coexistence phase. tuition for the generic case which is neither symmetricnor asymmetric. In general one may expect that thestress species 1 suffers from 2 is not exactly the same asthe stress species 2 suffers from 1, but that they are re-lated to each other. For example, if there is some nicheoverlap between species 1 and 2, but the niche of 1 iswider than the niche of 2, one expects c , < c , , buttheir values are correlated.Naively, one may guess that symmetry is a “fragile”property, so any deviation from perfect symmetry willsend the system to the equivalence class of the asymmet-ric model. However, our argument allows us to realizethat the opposite is true. As long as the system allows fora finite fraction of symmetric links, the breaking of eachof them doubles the number of maximal cliques when thesystem is very close to its complete graph limit. Accord-ingly, as long as there is some symmetry in the problem( c i,j is positively correlated with c j,i ) one should expectan exponentially large number of SUs when α is O (1),although the coefficient of N in the exponent falls alongwith the degree of correlation. As we shall see below, theresult in the purely asymmetric case reflects a “miracu-lous” cancelation of terms, so this turns out to be thefragile case.In Ref. [16], we extended the Bollob´as-Erd¨os formulato the asymmetric case, showing that the number of SUs satisfies, SU ( N, S ) = (cid:18) NS (cid:19) ˜ p S ( S − (1 − ˜ p S ) N − S . (21)Interestingly, the only difference between (21) and (7) isthe factor of 2 in the second term, reflecting the fact thatfor a collection of S species to be noninterfering one needsall the S ( S − / c i,j s to be zero in the symmetric case,while in the asymmetric case c i,j and c j,i are picked atrandom so the number of independent links is doubled.As we shall see, this innocent looking modification hashighly nontrivial consequences.Implementing the method used for the symmetric case,one finds, SU ( N, S ) = e NF as ( β ) √ πN , (22)where, F as ( β ) = − (1 − β ) ln(1 − β ) − β ln( β ) (23) − αβ + (1 − β ) ln(1 − e − αβ ) . As before, this pair of formulas appear to suggest that,as long as both α and β are O (1), the number of SUs isexponential in N . However, we shall see that in this case F as ( β ) = 0 and the actual large-N asymptotic turns outto be non-exponential.The equation for β now reads F (cid:48) as ( β ) = − αβ − α ( − β ) − e αβ + ln(1 − β ) (24) − ln( β ) − ln(1 − e − αβ ) = 0 . Surprisingly, one can find an exact solution to thisequation, β = W ( α ) α , (25)where W is the Lambert W function, defined by W ( x ) exp[ W ( x )] = x . Plugging this into F , we find that F as ( β ) vanishes identically . This implies that the firstcontribution from the Laplace integral is O (1) (insteadof being exponential in N ) so we should repeat the ex-ercise from its starting point, keeping all the O (1) terms(omitting only O (1 /N ) and other small terms). F as in the controlling factor of (22) then takes the form, F as = − (1 − β + 12 N ) ln(1 − β ) − ( β + 12 N ) ln( β ) (26) − αβ + αβN − α β N + (1 − β ) ln(1 − e − αβ − βα N ) . To continue, let us write F as as, F as = F (0) as + F (1) as N (27)where F as is given in (23) and F (1) as = α β α β e αβ − − α β e αβ − − αβ + 12 log(1 − β ) + log( β )2 . (28)Evidently, the main contribution in the large N limit stillcomes from β given in (25) and, e F (1) as ( β ) = 12 (cid:18) log (cid:18) − W ( α ) α (cid:19) − W ( α ) (cid:19) , (29) meaning that there is no exponential growth of the num-ber of maximal cliques with N .Now we can implement the Laplace integral schemeto (22) (the sum over S is converted to an integral over N dβ ) to obtain, SU ( N ) = N (cid:90) dβ e NF as ( β ) √ πN = N √ π e F (1) as ( β ) (cid:90) ∞−∞ e − N | ( F (0) as ) (cid:48)(cid:48) | β | ( β − β ) = e F (1) as ( β ) (cid:112) | ( F as ) (cid:48)(cid:48) | , (30)so SU ( N ) = αW ( α ) + W ( α ) (31)In the limit where α is O ( N ) considered in [16] thewidth of the Gaussian in the integration of (30) is 1 /N ,meaning that only a single large term in the sum of SU ( N, S ) over S contributes (see Fig. 5). In such acase there is no contribution from the integration aroundthe maximum, and the number of cliques is SU ( N ) ∼ e F (1) as ( β ) √ πN ∼ N ln / ( N ) (32)as shown in [16].The behavior of the number of SUs in different regimesof the Binomial model is depicted in Figure 5, and theresults of the corresponding Gamma model are shown inFigure 6. IV. DISCUSSION
In this paper we have studied the number of stable anduninvadable states in an ecological network. We assumeda local community which is coupled to a regional speciespool. In the local community the particular level of com-petition between any given pair of species was drawn atrandom.Many empirical ecological networks (in particular foodwebs [20] and networks with mutualistic interactions [21])were shown to admit a nontrivial structure (like modu-larity or nestedness) so their description as random net-works is problematic. Still, we believe that the analysispresented here is relevant to various aspects of the gen-eral problem. First, there are less evidence, as far as weknow, for a general structure in systems of competingspecies (see, e.g., [22]). Second, even if the mainland in-teractions are structured, there is no a priori reason toassume that this is the case on the island. A third point N u m b e r o f S U s FIG. 5. ln( SU ) vs. p , as obtained from a numerical summa-tion of the Bollob´as-Erd¨os formula for the asymmetric case,for N = 1000. In general the number of SUs decays with p ,with no exponential peak close to the fully connected limit.This general trend is superimposed on oscillations in the re-gion where p is O (1), since in this regime there is only oneinteger value of maximal clique sizes that dominate the sum,as explained in the text. To demonstrate that, thick pointswere added to mark the p values for which S = β N is 2 , , , (which is complementary to the second) is that, when theinteractions are inferred from empirical studies of localcommunities, one would like to understand what aspectsof these interactions are the result of the restriction of aregional system with a (possibly) different structure toits SUs.The model considered here is characterized by threeparameters: the number of species in the regional pool N , the mean value of the off-diagonal entries of the com-petition matrix C and the parameter that reflects theheterogeneity of the competition terms, ˜ σ . In the Bi-nomial model ˜ σ = C p/ (1 − p ), so in the limit when p = 1 − α/N , ˜ σ ∼ C Nα .
In the works of Mehta and Fisher [12] and Bunin [17]
FIG. 6. Number of SUs (averaged over 500 samples) vs. σ γ asobtained from an exact enumeration of all the stable and un-invadable combinations of species in the asymmetric Gammamodel for N = 20. Red points are actual results, each ob-tained from examination for stability and uninvadability ofall the 2 combinations of species, the dashed black line isjust to guide the eye. The asymmetric Gamma model is de-scribed by Eq. (3) with C = 1 and where each c i,j is pickedindependently from a Gamma distribution with mean one andstandard deviation σ γ . the average value of an (off diagonal) interaction matrixterm and the variance of these terms both are taken tobe of order 1 /N . This parameter regime is right on theborder of the regime defined by May’s stability criteriamentioned above. Translating this to the notations of ourpaper, one has C ∼ µ/N , say (when µ is O (1)), hence σ = µ /N α . So the regime of parameters covered bythe p = 1 − α/N limit of the Binomial model (with α order one) includes the regime considered in [12, 17] as aspecial case.The main outcome of the analysis presented here andin [16] is that the number of SUs grows exponentiallywith N if α is O (1), and the matrix in not purely asym-metric. If p is order one, or in the case of an asymmetric competition terms, the growth is subexponential, rang-ing from N ln( N ) dependency to sublinear growth. Thisimplies that, as long as an exponential number of SUs isrequired for a glass transition (as suggested by the anal-ogy with the random energy model presented in [12]),such a transition occurs only in the regime of very weakcompetition and very large systems.When ecologists consider high-diversity assemblagesand try to understand the forces that shape their struc-ture, they usually have in mind systems like tropicaltrees [23], coral reef [24] or plankton [25]. In these casesthe level of niche overlap between species is evidentlyquite high, as all these species are using the same set ofa few key resources, more or less in the same manner.Accordingly, one should expect these systems to be inthe regime where the interaction terms of the competi-tion matrix are O (1) (see, e.g. the recent study [26]),where the number of species in an SU scales logarithmi-cally with N [16], the number of SUs is subexponential,and there is no glass transition.To the best of our understanding, the parameterregime considered in [12, 17] and here corresponds to acompletely different scenario. This is the case of a com-munity with many species but with strong niche parti-tioning (say, many bird species with different beak size,eating different kinds of food) that still have weak compe-tition between species (due to some overlap in the typeof food they are eating, weak nest site competition ordue to predation by a common predator). Most ecolo-gist feel that the coexistence of many different species insuch a scenario needs no explanation (since the main is-sue they consider is the competitive exclusion principle)but in fact there is still a theoretical problem, namelyMay’s complexity-diversity relationship, meaning thateven a community with very weak interactions will col-lapse when the number of species is large. Here we haveshown that in this case one should expect to see a lo-cal community with O ( N ) species ( β is order one), andpossibly some kind of a glass transition. The relevance ofthis theoretical framework to empirical systems appearsto be an open problem. [1] G. F. Gause, The struggle for existence (Williams andWilkins, Baltimore, 1934).[2] G. Hardin, Science , 1292 (1960).[3] D. Tilman,
Resource competition and community struc-ture (Princeton Univ. Press, 1982).[4] R. M. May, Nature , 413 (1972).[5] R. H. MacArthur,
The theory of island biogeography ,Vol. 1 (Princeton University Press, 1967).[6] J. Losos and R. Ricklefs, Ecology , 2806 (2010).[7] D. A. Kessler and N. M. Shnerb, Journal of StatisticalPhysics , 861 (2007).[8] O. Ovaskainen and B. Meerson, Trends in ecology & evo-lution , 643 (2010).[9] R. Lande, S. Engen, and B.-E. Saether, Stochastic popu- lation dynamics in ecology and conservation (Oxford Uni-versity Press, 2003).[10] M. Danino, N. M. Shnerb, S. Azaele, W. E. Kunin, andD. A. Kessler, Journal of Theoretical Biology , 155(2016).[11] D. A. Kessler and N. M. Shnerb, Physical Review E ,042705 (2015).[12] C. K. Fisher and P. Mehta, Proceedings of the NationalAcademy of Sciences , 13111 (2014).[13] B. Derrida, Physical Review B , 2613 (1981).[14] M. Gilpin and T. Case, Nature , 40 (1976).[15] H. Lischke and T. Lffler, TheoreticalPopulation Biology, in press (2017),http://dx.doi.org/10.1016/j.tpb.2017.02.001. [16] Y. Fried, D. A. Kessler, and N. M. Shnerb, ScientificReports (2016).[17] G. Bunin, arXiv preprint arXiv:1607.04734 (2016).[18] B. Bollob´as and P. Erd¨os, Mathematical Proceedings ofthe Cambridge Philosophical Society , 419 (1976).[19] I. D. Rozdilsky and L. Stone, Ecology Letters , 397(2001).[20] B. Drossel, P. G. Higgs, and A. J. McKane, Journal ofTheoretical Biology , 91 (2001).[21] S. Suweis, F. Simini, J. R. Banavar, and A. Maritan,Nature , 449 (2013).[22] I. Volkov, J. R. Banavar, S. P. Hubbell, and A. Maritan,Proceedings of the National Academy of Sciences ,13854 (2009). [23] H. Ter Steege, N. C. Pitman, D. Sabatier, C. Baraloto,R. P. Salom˜ao, J. E. Guevara, O. L. Phillips, C. V.Castilho, W. E. Magnusson, J.-F. Molino, et al. , Science , 1243092 (2013).[24] S. R. Connolly, M. A. MacNeil, M. J. Caley, N. Knowl-ton, E. Cripps, M. Hisano, L. M. Thibaut, B. D. Bhat-tacharya, L. Benedetti-Cecchi, R. E. Brainard, et al. ,Proceedings of the National Academy of Sciences ,8524 (2014).[25] M. Stomp, J. Huisman, G. G. Mittelbach, E. Litchman,and C. A. Klausmeier, Ecology , 2096 (2011).[26] F. Carrara, A. Giometto, M. Seymour, A. Rinaldo, andF. Altermatt, Ecology96