Robust coexistence in competitive ecological communities

Domain: neuroscience
Published: 2026-02-11
Editorial outcome: minor_revision

Abstract
Darwin already recognized that competition is fiercest among conspecifics, a principle that later made intraspecific competition central to ecological theory through concepts such as niche differentiation and limiting similarity. Beyond shaping coexistence, strong intraspecific competition can also stabilize community dynamics by ensuring that populations return to equilibrium after disturbance. Here we investigate a more fundamental question: how intraspecific competition influences the very existence of a steady state (feasibility) in large random ecological communities dominated by competition. We show that, in analogy with classical results on stability, there is a critical level of intraspecific competition above which a feasible steady state is guaranteed to exist. We derive a general expression for the probability of feasibility and prove that, asymptotically (as species number grows), the transition to stability occurs before the transition to feasibility with probability one. Thus, in large competitive communities, any feasible equilibrium is automatically stable. This ordering persists even when many species in the initial pool cannot coexist and extinctions occur: the dynamics prune the community, shifting feasibility and stability thresholds but never reversing their order. These results imply that large competitive communities generically converge to a globally stable equilibrium, making sustained oscillations or chaos unlikely—consistent with experimental observations.

Introduction
The notions of intraspecific competition and population self-regulation play a central role in ecological theory, leading to foundational concepts such as limiting similarity and niche differentiation1–6. They also are crucial for coexistence: for example, May7 showed that ecological dynamics can be stabilized by imposing sufficiently strong intraspecific competition, and that for large systems with random species interactions the stability transition is universal (independent of the probability distribution used to model interactions) and sharp (always achieved at the same critical value of intraspecific competition). Importantly, the existence of an equilibrium (feasibility) is a necessary condition for coexistence. Further, results on stability can dramatically change when accounting for feasibility. The first contribution to this problem can be traced back to a response to May’s work on stability8, showing that in competitive Lotka-Volterra models with random interactions the probability of stability, conditioned on the presence of a feasible equilibrium, does increase, rather than decrease with system size. Since then, significant advances have been gradually accumulating in the literature9–15 (Supplementary Note 3, Sec. 2). However, the progress on the study of feasibility has yet to match that reached for stability7,16–18. Here we show that accounting for feasibility is key to determining the expected dynamics of large competitive systems.

We compute the probability that a system would display a feasible equilibrium given its level of intraspecific competition, size, and the statistics on the random distribution of species interactions. We show that the transition to feasibility is universal, as for stability, but is not sharp—thereby extending, generalizing and refining previous results8,9,12,14,19–22. More importantly, we show that, in large competitive systems, when we increase the strength of intraspecific competition, the transition to stability precedes that to feasibility8,9,20. Therefore, in these large systems feasibility implies stability. This result holds even when the initial pool of species cannot coexist, and extinctions ensue. That is, during ecological dynamics, species extinctions will alter the values for the stability and feasibility thresholds, but not their ordering. Differently from May, who focused on local asymptotic stability of an equilibrium, here we consider the stronger criterion of global stability—robust coexistence in which populations rebound from any perturbation short of extinction.

Overall, these findings have far-reaching implications for the ecological dynamics of large random competitive ecosystems. They rule out coexistence through limit cycles or chaos23,24, and instead indicate that the community dynamics converge to a highly-robust equilibrium, such that (1) any perturbation short of extinction would be buffered, (2) the community resists invasions, and (3) it can be assembled from the ground up25,26.

Results
Stability and feasibility thresholds
We consider the Generalized Lotka-Volterra27–30 (GLV) model for n interacting populations:



where intrinsic growth rates r i are organized in the vector r, and interactions A ij in the matrix A. We say that Eq. (1) has a feasible (i.e., biologically attainable) equilibrium whenever there is a x⋆ > 0 such that x_ i = 0 for all populations, and that such a feasible equilibrium is asymptotically stable if dynamics return to it following a perturbation.
Following May7, we model A as the sum of two matrices, A = αI + B, where B collects the inter-specific interactions (B ij, i≠j). I is the identity random matrix, and α ≥ 0 the intraspecific competition in the system, assumed to be the same for all species (relaxing this assumption does not change the results, see Supplementary Note 3, Sec. 9).
Now we consider increasing α from zero to an arbitrary large positive value, and track the effects of this change on coexistence. We denote α S as the level of intraspecific competition that is sufficient to guarantee the (global) stability of an equilibrium (Methods, and Supplementary Note 1, Sec. 3). Importantly, for any α ≥ α S stability is maintained, i.e., once surpassed this threshold, dynamics always converge to a globally stable equilibrium. Here we ask whether, in analogy to stability, there exist a critical level α F that a) guarantees the existence of a positive equilibrium, and b) further increasing α ≥ α F retains the feasibility of the equilibrium. In particular, we want to probe the relationship between α F and α S.
We start with a toy example (Fig. 1, n = 6, see Supplementary Note 3 and Supplementary Fig. 3 for an example with two populations), illustrating the effect of increasing α on feasibility and stability. We first sample a baseline interaction matrix B, and compute the thresholds for stability, α S, and feasibility, α F (Methods, and Supplementary Note 1, Sec. 3, and Supplementary Note 2, Sec. 4), marked in Fig. 1 with green and yellow triangles, respectively. We then construct three interaction matrices A that differ only in their intraspecific competition values, α I, α II, and α III (Fig. 1). For each system, we integrate the dynamics in Eq. (1) from two different initial conditions (Fig. 1 bottom panels, solid vs. dashed lines), revealing three qualitatively distinct dynamical regimes.
System I (Fig. 1, left) undergoes extinctions and displays alternative stable states, indicating that the coexistence equilibrium is not feasible, and the resulting sub-community equilibrium is not globally stable (indeed, neither stability nor feasibility thresholds have been surpassed, since α I < α S < α F). System II (Fig. 1, middle), also displays extinctions (as there is no feasible equilibrium), but starting from different initial conditions no longer leads to different outcomes, because the equilibrium of the sub-community is globally stable (α S < α II < α F, and see Supplementary Note 1, Sec. 3). Finally, system III (Fig. 1, right), exhibits all populations coexisting stably (α S < α F < α III). From this single realization, we can see that stability is achieved before feasibility. Next, we show that this is the case for large random ecological communities (but see Supplementary Note 2, Sec. 3 for results on small or arbitrary competitive systems).

Feasibility implies robust coexistence almost surely
As long as populations are statistically equivalent (i.e., there is no a priori structure making some populations behave differently from the rest), the distribution of the equilibrium values, y⋆, has always the same mean and correlation structure (Methods). Moreover, when n is large, we can assume the components of y⋆ to be approximately normally distributed9 and uncorrelated (“Methods”, and Supplementary Note 3, Sec. 7). Then, the probability of feasibility can be taken as the product of the probabilities that each component is positive, and can thus be written as the cumulative distribution function of a standard normal distribution, raised to the nth power (Supplementary Note 3, Sec. 7):



Substituting this expression into Eq. (3), we obtain an extremely accurate analytical approximation of the probability of feasibility (Fig. 3, black curves): the probability p F(α) is vanishingly small when α = 0 (Supplementary Note 3, Sec. 3), and increases smoothly to 1 with increasing α. Eq. (4) shows that the details of the distribution of interactions (Fig. 3, linetypes), once accounted for the correlation (Fig. 3, colors), are irrelevant when the system is large. This means that asymptotically, the curve describing the probability of feasibility is universal. Interestingly, exactly as a negative correlation between interaction strengths improves the chance of stability16, it also increases the probability of feasibility21. Contrary to current expectations14, the effect of the pairwise correlation ρ on the probability of feasibility remains significant for large systems (n = 1000 in Fig. 3). Additionally, the transition to feasibility also remains smooth compared to stability (Fig. 3, and Supplementary Fig. 10, see Supplementary Note 3, Sec. 8 for a detailed discussion) for large systems.
Importantly, our calculations show that for large systems the amount of self-regulation required for feasibility is high enough (i.e., α S < α F) to yield global stability (as argued before9,14,20,32,33). Then, a feasible system will be stable with probability 1, and thus p F(α) measures the probability of coexistence for a large, random competitive system. These results have profound implications for dynamics. Competitive GLV systems can give rise to any type of dynamics, including limit cycles (when n≥3) and chaos (n≥4)23,24. Importantly, however, non-equilibrium coexistence still requires the existence of a feasible equilibrium25. But we have seen that feasibility generically implies global stability, and therefore non-equilibrium coexistence is precluded.
When the system is not feasible, extinctions will ensue, pruning the original pool of species to a subsystem with different thresholds for feasibility and stability (Fig. 2c). As such, in principle, any type of dynamics could be observed. Next, we show that, despite the change in values of α S and α F caused by extinctions, the ordering α S < α F is maintained almost surely, provided that n is large. As a result, dynamics still converge to stable equilibria even when feasibility is not attained for the initial species pool.

Robust coexistence emerges through community dynamics
So far, we have investigated how coexistence is affected by increasing α. However, intraspecific competition models the effects of a species on itself. As such, when species are lost to extinction, the level of intraspecific competition should not be affected: the coefficients on the diagonal do not change when the matrix is pruned of the populations that have gone extinct. Thus, we turn our attention to a pool of n populations where we fix α < α S < α F (analogous to system I in Fig. 1), and consider “top-down assembly” (where all the species in the pool are introduced in the environment at the same time, Fig. 4a). Since we are below the feasibility threshold, some populations go extinct and, as a result, α F and α S changes (both values depend on n; Fig. 2c). We re-compute their values every time an extinction occurs (plotted in Fig. 4b), and find that their initial ordering is maintained throughout the dynamics until the system converges to a stable equilibrium after ecological dynamics have elapsed. Furthermore, we show (Supplementary Note 3, Sec. 12) that the same equilibrium can be achieved via “bottom-up assembly”, through successive invasions26.
To confirm that these observations are a generic feature, we plot the ratio α F/α S when we run the same simulation in the previous example for many draws of random systems with the same characteristics (Fig. 4c). We see that the ratio is kept above one, implying that the relationship α S < α F is maintained throughout the dynamics for all the simulations. This shows that, with high probability, the dynamical outcome of large random competitive systems is always the same: a sub-set of the species pool will coexist robustly at an equilibrium that is characterized by α S < α F < α; the equilibrium is stable, and the last few populations to go extinct (if started in phase I), or all the extinct populations (if started in phase II) cannot re-invade the system.

Methods
A canonical form for competitive systems
We consider the competitive GLV system in Eq. (1)

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript addresses a theoretical question whether feasible equilibria in large random competitive ecological systems are also dynamically stable. Based on analytical and numerical investigations of random matrix models, the authors concluded that feasibility implied stability in large systems governed by Lotka–Volterra dynamics with random interactions. The central claim is that the transition to dynamical stability consistently precedes the transition to feasibility, thus any feasible equilibrium is necessarily stable. These results provide a plausible explanation for the empirical observations in nature that communities of organisms interacting primarily through competition—such as species at the same trophic level—often exhibit stable equilibria, whereas sustained cycles and chaotic dynamics are more commonly associated with predator-prey systems or structured populations.\nWe found these results mathematically elegant and conceptually important. The analytical framework and insights offered here contribute to bridging ecological theories and observations, and the paper shall stimulate further research on the mechanisms underlying ecological stability. However, we have several concerns about the presentation, generality, and robustness of the results.\n(1) The definition of stability needs to be further clarified. The ecological literature covers a broad set of different definitions of stability (local stability, resistance, recovery, variability, etc.), therefore it is important to clarify the meaning of ‘stability’ in the introduction as early as possible. Reading the paper, we had been thinking about ‘local stability’, but later it turned out to be ‘global stability’ after reading Fig. 1 and line 71. We believe that it is important to clarify in the Introduction which measure of stability will be focused on, why this measure is chosen (from the broad set, as well as between global vs. local stability), how they are defined and quantified, etc. In this regard, the authors may incorporate some general discussion form the supplementary materials explaining local and global stability.\n(2) The paper included many arguments that was presented as general conclusions (e.g., L28, L93, L139, L172), which however – to the best of our understanding – hold only under certain contexts. For instance, in the abstract the statement “feasible systems will necessarily be stable” (line 28) was presented as a strong conclusion, but clearly it holds only under the study context using large competitive systems with random interaction strengths. Even in this context, the argument was supported primarily by asymptotic approximations and simulations, and neither the main text nor the appendix provided formal mathematical proofs establishing that α_S<α_F holds for all random matrices within the model class. Our preliminary simulations revealed counterexamples (i.e., feasible but unstable equilibria) even for n=50. So we recommend to be more cautious with conclusions, perhaps – feasibility implies stability with high probability in large random systems. Similarly, the main text and figure captions included quite a number of strong arguments with “always” “invariably” and “necessarily”. We recommend the authors to check throughout and avoid overstatement of the generality.\n(3) The robustness of the results needs more examination. The paper did not provide sufficient information on the uncertainty involved in the simulations, such as confidence intervals, failure rates, or sensitivity bounds on the probability that feasibility implies stability (e.g., Figs. 2 & 3). Without such information, it is difficult to assess how robust the results are to variation in model settings. In particular, since the analytical solutions were obtained under the assumption of large communities with random interactions, we were wondering how large (cid:1) needs to be for the asymptotic results to reach a certain precision. The bottom-left panel in Fig. 2 provided some insights into this point, but an extended investigation would be useful for understanding the behaviors of the analytical solutions.\n(4) Overall, the paper was quite brief, in terms of both results and discussion. This had made additional difficulties in the understanding of the paper. We recommended the authors to expand the manuscript by adding necessary explanations (e.g., context of analyses and conclusions, as explained above) and expanding the analyses by incorporating some results from the appendix. In particular, the main text presented the solutions mainly for feasibility, but little for the global stability. Is there a solution for the latter? If so, please provide them. If not, explain how were the thresholds obtained for stability in Figs. 1 & 2 (triangles)?\nMoreover, the discussion should be substantially expanded to make connections to the literature. For instance, the recent paper by Hu et al. (2022) illustrated transitions from stable & feasible coexistence, to stable and unfeasible coexistence, and to fluctuating coexistence, as the number of species and average interaction strength increased. How do the findings in the current paper relate to and differ from Hu et al. (2022)? Also, many ecological models can produce scenarios where feasible equilibrium is not stable, such as the simplest two-species model with L-V type competition. How to reconcile the present study with these cases? Again, this needs to clarify the context of the results in this paper and reconcile them with broader models and patterns under different settings in the literature.\nMinor comments:\nLine 42. The statement “the stability of an equilibrium requires its existence (feasibility) in the first place” did not seem correct. As an example, Fig. 1 clearly showed the case (case II) where the community was infeasible but globally stable. Please specify the context where such argument should hold.\nLine 104: vague arguments. What does it mean by ‘several transitions’? Please specify it if it is important, otherwise simplify the text.\nLine 147: “the levels of intraspecific competition remain more or less constant over time” – why? Environmental fluctuations can certainly lead to variation in carrying capacity (K_i) and thus intraspecific competition (as defined here, i.e., α_ii=1/K_i). We felt that this argument argument was an assumption, but presented as a conclusion. This should be clarified.\nLine 206: why ‘proportional to each other’? In the equation, y* also occur in the dominator, making it nonlinear.\nFig. 1: this is a very clear and interesting illustration. We think that such transitions could (though not necessarily always) occur already in L-V models with two competitors – that is, alternative stable states at low α_ii, deterministic exclusion at higher α_ii, and coexistence at very high α_ii. If so, this may help the readers to make the bridge to existing knowledge.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "I found the work presented in ‘Robust coexistence in ecological competitive communities’ very interesting, elegantly presented and quite exciting frankly. The authors demonstrate the significance of self-regulation as a mechanism for community persistence. In particular, the authors use random matrices to generate competitive ecological communities and show that increasing self-regulation (intra-specific competition) first achieves stability (fixed point equilibrium) and then feasibility (positive density of species at equilibrium). The order of stability being attained prior to feasibility (along the self-regulation gradient) implies that feasibility guarantees stability. An important corollary of this work is that competitive communities (that are large and have randomly distributed interaction strengths, a fundamental assumption and caveat) will have stable equilibria rather than non-equilibrium steady states.\nI have one comment which touches upon the basic assumptions of this work. The authors use a single value for self-regulation, alpha, across species. Do I understand this correctly? This scenario would imply that all species share the same strength of self-regulation, while the strength of inter-specific interactions can vary. I wonder what would happen if there was a mean alpha for a matrix, rather than a fixed alpha for all species? How strongly would this impact the insights (obviously given different assumptions about the distribution of alpha’s)? Adding to this, what would happen if these alpha values were somehow correlated to the inter-specific interaction terms? For instance, a weak competitor with respect to other species, would presumably also have weak self-regulation to increase its survival chances. I understand that the authors simplified a very complex problem and I appreciate the value in such an approach. However, since the fundamental ecological parameter studied here is self-regulation, I would like the authors to explore the limitations of their simplifying assumptions regarding this parameter. Ultimately, it is important to understand the generality of their results to relaxing this assumption.\nBesides this and to the extent that I was able to follow the mathematical analyses in the supplement, I found the justification of the results persuasive. Moreover, I liked how the manuscript was written; concisely, clearly and with a linearly progressing narrative. Overall, I really enjoyed reading this study and I would be happy to see it published so as to stimulate a debate in the broader community of interested ecologists.\nMy only regret was the absence of an ecological interpretation for the results, beyond the very obvious statements. For instance, as I mention in a comment below, I see some relation of certain results to niche theory; it would be very insightful to link the two. I believe there’s plenty of scope to expand on such nuances in the manuscript and I feel the authors would be doing ecologist readers a great service to expand on this aspect. Ultimately, given the journal and the style of the manuscript, I view this suggestion more as a ‘nice to have’ than a requirement for publication.\nI have made more specific comments below, mostly for clarification purposes.\nTitle\nWouldn’t ‘Robust coexistence in competitive ecological communities’ sound better?\nAbstract\nThe final sentence is quite strong. The absence of non-equilibrium steady states rests on certain assumptions (large, random interaction matrices). Consider toning down the absoluteness of the finding. Also, the authors could explain/expand on the notion of ‘robust’ coexistence (i.e., that this equilibrium cannot be invaded).\nS1\nl. 53 the authors can clarify the notion of ‘universal’ convergence to feasibility at this point so that it becomes easier to follow this line of reasoning.\nl. 60 ‘resilient to any perturbation sort of extinction’: this statement is not very clear\nS2\nIn fig1, add a brief description of the top panel, to make it immediately obvious that this refers to the interaction matrix entries, with the shade reflecting the magnitude of the parameters.\nl. 83 ‘neither feasible’\nS3\nl.124-127 & eq. (4): it appears that as n increases, only the correlation (rho) will impact the variance and hence the probability of feasibility. If that is the case, I am not sure I understand l. 125, where I perceive the message as being that correlation (besides the mean and variance of the interaction matrix) is irrelevant.\nl.131 e.g., what does P_F=0.5 mean for a community? In other words, what determines whether or not a community will be feasible? Are there other aspects of the interaction matrix distribution (higher moments?) whose consideration will make this outcome more definite? Is feasibility higher dimensional that stability and we require additional dimensions to make its transition sharp?\nl. 142-144 I have a question which I am not sure I can present clearly, but I will try because I think there’s some value in answering it (if possible). We have the following information: for a given n and fixed alpha and rho, we can find the P_F (fig2, eqs. 3&4). If P_F<1, does this imply the community will ‘prune’ itself, to a lower n and attain P_F=1 (l. 143)? It would be interesting to seek a pattern linking these alpha’s to the ‘pruned’ community size n, for which P_F=1. If I understand correctly, it won’t matter whether these certainly feasible communities are the result of ‘pruning’ or bottom-up assembly. Perhaps, there are insights to be gained from a figure linking n to alpha such that P_F=1. Also, linked to this, I was wondering whether the slope of these curves (fig2) changes with n, i.e., does very large n lead to a sharp transition for instance?\nfig2: the correlation in the interaction matrix can be seen as a proxy for specialisation, if I understand correctly. Positively correlated entries imply no trade-offs, while negative correlation indicate trade-offs and hence some kind of specialisation. This can be linked to the notion of the niche, which would seem to agree with the results: niche structure guarantees stability and facilitates feasibility (rho=-1), niche overlaps (rho=1) require much stronger self-regulation for the same outcomes. I think the manuscript would benefit significantly if such links between the mathematics and ecology were more pronounced and greater in number."}, "author_response": "Responses to Reviewers\nManuscript NCOMMS-25-39378\nReviewers 1 & 2\nThis manuscript addresses a theoretical question whether feasible equilibria in large\nrandom competitive ecological systems are also dynamically stable. Based on ana-\nlytical and numerical investigations of random matrix models, the authors concluded\nthat feasibility implied stability in large systems governed by Lotka–Volterra dynamics\nwith random interactions. The central claim is that the transition to dynamical sta-\nbility consistently precedes the transition to feasibility, thus any feasible equilibrium\nis necessarily stable. These results provide a plausible explanation for the empirical\nobservations in nature that communities of organisms interacting primarily through\ncompetition—such as species at the same trophic level—often exhibit stable equilibria,\nwhereas sustained cycles and chaotic dynamics are more commonly associated with\npredator-prey systems or structured populations.\nWe found these results mathematically elegant and conceptually important. The\nanalytical framework and insights offered here contribute to bridging ecological the-\nories and observations, and the paper shall stimulate further research on the mech-\nanisms underlying ecological stability. However, we have several concerns about the\npresentation, generality, and robustness of the results.\nWethankthereviewersforthesuccinctsummaryandtheirpositiveevaluationofourmanuscript,\nas well as for their interest in our work and their helpful suggestions. We have found the comments\nwell-chosen and constructive: addressing them helped us improving the clarity of our results, and\nled us to add a new material to our contribution.\nMajor comments\n1. The definition of stability needs to be further clarified. The ecological literature\ncovers a broad set of different definitions of stability (local stability, resistance,\nrecovery, variability, etc.), therefore it is important to clarify the meaning of\nstability in the introduction as early as possible. Reading the paper, we had\nbeen thinking about local stability, but later it turned out to be global stability\nafter reading Fig. 1 and line 71. We believe that it is important to clarify in\nthe Introduction which measure of stability will be focused on, why this measure\nis chosen (from the broad set, as well as between global versus local stability),\nhow they are defined and quantified, and so on. In this regard, the authors may\nincorporate some general discussion from the supplementary materials explaining\nlocal and global stability.\nin the text. We have now added a sentence in the introduction to make clear that we focus\non global asymptotic stability:\n“Differently from May, who focused on local asymptotic stability of an equilibrium,\nhere we consider the stronger criterion of global stability—robust coexistence in\nwhich populations rebound from any perturbation short of extinction.”\nAdditionally, we have added a section to the methods in the main text where the threshold\nfor global stability is derived.\n2. The paper included many arguments that were presented as general conclusions\n(for example, line 28, line 93, line 139, line 172), which however hold only under\ncertain contexts. For instance, in the abstract the statement “feasible systems\nwill necessarily be stable” (line 28) was presented as a strong conclusion, but\nclearly it holds only under the study context using large competitive systems with\nrandom interaction strengths. Even in this context, the argument was supported\nprimarily by asymptotic approximations and simulations, and neither the main\ntext nor the appendix provided formal mathematical proofs establishing that\nα < α holds for all random matrices within the model class. Our preliminary\nS F\nsimulations revealed counterexamples, that is feasible but unstable equilibria,\neven for n = 50. We recommend being more cautious with conclusions, perhaps\nphrasingthemasfeasibilityimpliesstabilitywithhighprobabilityinlargerandom\nsystems. Similarly, the main text and figure captions included quite a number of\nstrong arguments with terms such as “always”, “invariably” and “necessarily”.\nWe recommend the authors to check throughout and avoid overstatement of the\ngenerality.\nWe agree that the use of definitive words such as “always”, “invariably”, “necessarily”, and\nso on, is inconsistent with the derivations presented in the text, which are based on proba-\nbilities and asymptotics. We have carefully revised the manuscript and changed the wording\nthroughout to reflect the probabilistic nature of our results.\n3. The robustness of the results needs more examination. The paper did not provide\nsufficient information on the uncertainty involved in the simulations, such as\nconfidence intervals, failure rates, or sensitivity bounds on the probability that\nfeasibility implies stability. Without such information, it is difficult to assess how\nrobust the results are to variation in model settings. In particular, since the\nanalytical solutions were obtained under the assumption of large communities\nwith random interactions, we were wondering how large n needs to be for the\nasymptotic results to reach a certain precision. The bottom-left panel in Fig. 2\nprovided some insights into this point, but an extended investigation would be\nuseful for understanding the behaviors of the analytical solutions.\nWe estimate a probability of feasibility, and show that our estimate matches simulations very\naccurately when the number of interacting populations is large. As the reviewers point out,\nthis leaves the problem of how large n must be for the asymptotic results to be valid in\npractice. To provide further insight into this question, we now include a figure (Fig. 2 in\nmanuscript, see Fig. 1 here) in the main text bridging the toy example and the analytical\nresults. This figure (panel d) shows that the probability of observing feasible but unstable\n2\nequilibria decreases rapidly with n, across different correlation values. Note that for n = 15,\nthe probability that a feasible equilibrium is not stable is below 0.01. The simulations in\nthe main text are performed for n = 1000. Additionally, (Fig. 2) we provide an analysis\n(a)\n1\n0\n0 1 α F 3\nIntraspecific competition\nmuirbiliuqE\n(b)\n2\n1\n5 10 15 20\nNumber of popoulations n\nsnoitisnart\nfo\n#\nρ=1\nρ=0\nρ= 1\n−\n(c)\n3\n2\n1\n0\n101 102 103\nNumber of populations n\nα\ndlohserhT\ne\n0.3\n(d)\nα α\nF S\n0\ne e 5 10 15 20\nα<\nα\nP\nS\nF\n(cid:1)\n(cid:0)\nFigure 1: Auxiliary simulations bridging toy example and analytical results. (a) Bi-\nfurcation diagram of equilibrium abundances as a function of intraspecific competition α. Shaded\nregionscorrespondtoparametervalueswherenofeasibleequilibriumexists, whileunshadedregions\nmark feasibility. For a system with four populations (colors), we see that when α = 0 no feasible\nsolution of Ax = (αI+B)x = 1 exists; as α increases, feasible equilibria can appear and disappear\nmultiple times, but for sufficiently large α α feasibility is always guaranteed.\nF\n≥\n(b)–(d) Results\nfrom simulations across many system sizes. For each n, we sample many random interaction ma-\ntrices with entries drawn in correlated pairs from a bivariate normal distribution with mean 0, unit\nvariance, and correlation ρ (colors), and calculate both feasibility and stability thresholds. (b) The\naverage number of transitions in and out of feasibility (y-axis) decreases with n and converges to\none, indicating that for large systems feasibility is excluded below α and guaranteed above it.\nF\n(c) Normalized thresholds α = α /√n 1 and α = α /√n 1 as a function of n. While α\nF F S S S\n− −\nconverges to a constant, α grows without bound. Triangles mark α for n = 103 (d) Probability\nF S\nthat α\nF\n< α\nS\nas a function(cid:101)of n, i.e. the chance t(cid:101)hat a feasible but unstable equilibrium exis(cid:101)ts.\nThis probability decrease(cid:101)s rapidly with system size. (cid:101)\nof accuracy of analytical approximation to empirical cumulative distribution function of α\nF\n(i.e., probability of feasibility) for different choices of correlation (colors), distribution (beta,\ndiscrete, and normal), and number of populations (20, 40, 100, and 200). The analytical\napproximation becomes extremely accurate for n = 200, and the simulations in the main\ntext are performed for n = 1000. We did not attempt to quantify discrepancies found for\nsmall sizes, because our emphasis is squarely on the asymptotic case of many interacting\npopulations.\n4. Overall, the paper was quite brief, in terms of both results and discussion. This\nhad made additional difficulties in the understanding of the paper. We recom-\nbetween the analytical approximation to the cumulative distribution of the feasibility threshold α\nF\nand simulation results, for different numbers of species (n = 20,40,100,200), correlation values\n(colors), and distributions of interaction strengths (beta, discrete, normal). The approximation is\nalready highly accurate by n = 200; in the main text simulations were performed at n = 1000,\nwhere agreement is essentially perfect.\nmend expanding the manuscript by adding necessary explanations and expanding\nthe analyses by incorporating some results from the appendix. In particular, the\nmain text presented the solutions mainly for feasibility, but little for the global\nstability. Is there a solution for the latter? If so, please provide them. If not,\nexplain how were the thresholds obtained for stability in Figs. 1 and 2 (triangles)?\nWe recognize that parts of the paper would have benefitted from expanded explanations. To\naddress this, we have:\n4\n• Added an auxiliary figure (Fig. 1 here) in the main text connecting toy simulations and\nthe mathematical results.\n• Added two subsections to the Methods section: one deriving the global stability thresh-\nolds,andoneshowingtheequivalencebetweenintraspecificcompetitionandthevariance\nof the interaction distribution.\n5. Moreover, the discussion should be substantially expanded to make connections\nto the literature. For instance, the recent paper by Hu et al. (2022) illustrated\ntransitions from stable and feasible coexistence, to stable and unfeasible coex-\nistence, and to fluctuating coexistence, as the number of species and average\ninteraction strength increased. How do the findings in the current paper relate\nto and differ from Hu et al. (2022)?\nThis is an excellent point, which should be included in the manuscript. We have added a\nparagraph to the discussion section detailing how our findings relate to Hu et al. (2022). The\nmain difference between our model and the model used in Hu et al. (2022) is the inclusion of\na constant immigration term in their model, which is absent in ours (we analyze the classic\nGLV model). The addition of constant immigration indeed explains the discrepancy between\nthese results.\nWhen constant immigration is added to all species, as done by Hu et al., extinction of any\nspecies is prevented. Take a system that in the absence of immigration would result in one or\nmorealternativestablestates; immigrationpreventsspeciesfromreachinganyofthesestates,\nwith the result that the community hovers around one of the points for a while, and then\njumps to the other, in an endless sequence. Hu et al. clearly show that without immigration\nthe behavior would be quite different: in their Figure S1, they show simulations in which\nimmigration is stopped at t = 1000. In all cases, the communities collapse, resulting in stable\ndynamics for the remaining species, exactly as expected according to our results.\nTo situate both results in a common framework, and see exactly when results converge—\nwith dynamics changing smoothly in transition from one model to another—we ran further\nsimulations (Fig. 3) where we consider the community in the work by Hu et al.: n = 50\nspecies with interactions sampled from a uniform distribution U[0,2a], where a is the mean\ninteraction strength between species (which also determines the variance of interactions, as\nthese two values increase proportionally). To reproduce their results (Fig. 1 in their paper),\nwe integrate the dynamics for the same initial conditions, and several values of the mean\nstrength from 0.08 to 0.64 (Figure 3, columns). To see how their results are highly dependent\non the immigration parameter, we integrate the same community for several values of the\nimmigration rate D (rows).\nLooking at Fig. 3 we se that when D = 0 (first row, red rectangles), their model collapses to\nthe GLV model studied in our manuscript. Importantly, their parameter a controls the vari-\nanceofinteractions, whichwehaveshown(Methods)isequivalenttovaryingtheintraspecific\ncompetition (inversely proportionally; see Methods). That is, the smaller their parameter a,\nthe larger our parameter for intraspecific competition. Indeed, we observe feasibility (and\nglobal stability) for the smallest value of their parameter a, which would correspond to the\nlargestvalueofintraspecificcompetition(topleftpanel),butfeasibilityislostforlargervalues\n(four remaining panels of the first row). These cases correspond to the results we present in\nSection 5 (robust coexistence emerging even when feasibility is not present). When we depart\nfrom our model, D > 0 (rows second to last), constant immigration prevents extinctions to\ndifferent immigration levels D, columns show increasing mean interaction strength a for n = 50\nspecies with interactions sampled from U[0,2a]. When D = 0 (first row, red rectangles), the\nmodel reduces to the classic GLV studied here: feasibility is present at low a (high α), and breaks\nfor higher a values. In all cases, however, robust coexistence emerges after species extinctions,\nconsistent with our results. For D > 0, immigration prevents extinction and produces persistent\nmetastable dynamics, consistent with Hu et al. (2022).\noccur, making the system bounce between metastable states for large enough values of a and\nsmall enough values of D.\nOverall, we hope the reviewers will find this discussion helpful to understand that the results\npresented here and those in Hu et al. (2022) are not mutually exclusive. Without constant\nimmigration, we recover the results presented in our manuscript; when constant immigration\nis added (thus deviating from the classic GLV model), things become more complicated, and\ndynamics show wide fluctuations reminiscent of heteroclinic cycles, as shown by Hu et al. To\nbriefly discuss these matters, we have added the following text to the discussion:\n“Recently, it has been shown that a dynamical phase characteristic of persistent\nfluctuationscanariseforthehigh-diversity,high-interactionstrengthlimitinbacte-\nrial systems dominated by competitive interactions. Furthermore, this observation\nmirrors dynamical phases arising in a GLV model with dispersal. Crucially, these\nfindings hinge on the presence of constant immigration—without it, irrespective of\ninteraction strength, robust coexistence emerges for large communities (see Fig. S1\ninforresultswhenn = 50). Assuch,ourfindingsarenotincontradictionwiththese\nprevious results—including constant immigration for all species precludes extinc-\ntions, thereby preventing the system from relaxing to the equilibria studied here.\nTherefore, the results presented can be used to form an expectation for isolated\nsystems in which populations can go extinct, or systems in which immigration is\n6\nrare. For open systems in which immigration is widespread, outcomes could depart\nseverely from what is presented here.”\n6. Also, many ecological models can produce scenarios where feasible equilibrium is\nnot stable, such as the simplest two-species Lotka–Volterra competition model.\nHow to reconcile the present study with these cases?\nThis is indeed the case: even when choosing parameters at random, small systems often\ndisplay feasible, unstable equilibria. What we are showing in the manuscript is that, on\nthe other hand, for large competitive systems with random parameters, feasibility almost\nsurely guarantees stability — i.e., the self-regulation needed to achieve a feasible equilibrium\nis large enough to guarantee stability with high probability. Put in another way, to find a\nfeasible, unstable equilibrium we need to either have a very small size, or “fine-tune” some\nof the parameters. This result is central to the manuscript: while similar arguments were\nmade for example by Stone, our calculations show that for large systems the amount of\nself-regulation required for feasibility is typically high enough to yield global stability, and\ntherefore preventing cycles and chaos. We have added a figure (See Fig. 1, panels b, c and\nd – Fig. 2 in the manuscript) to show the asymptotic behavior of our results (see above), to\ngive readers a sense of what is the specific size for which we should expect them to apply.\nAdditionally, we have added a paragraph in the discussion relating our results with other\nlarge scale random competitive models in the literature (see above).\nMinor comments\n• Line 42. The statement “the stability of an equilibrium requires its existence\n(feasibility) in the first place” did not seem correct.\nWesimplifiedthestatementto: “Importantly, theexistence(feasibility)ofanequilibriumisa\nnecessary condition for coexistence.” We clarified this in relation to Fig. 1 in the manuscript,\nnoting that when the feasible equilibrium does not exist, it cannot be stable, although the\ndynamics converge to a reduced stable equilibrium.\n• Line 104: vague arguments. What does it mean by “several transitions”?\nWe have clarified this point and added an auxiliary figure (Fig. 1 panel a) illustrating these\ntransitions explicitly.\n• Line 147: “the levels of intraspecific competition remain more or less constant\nover time” – why? Environmental fluctuations can certainly lead to variation\nin carrying capacity (K ) and thus intraspecific competition (as defined here,\ni\ni.e., α = 1/K ). We felt that this argument argument was an assumption, but\nii i\npresented as a conclusion. This should be clarified.\nIndeed, this assumption was not presented as such in writing, and we have modified the text.\nThe intraspecific competition models the effects of a species on itself. As such, when species\nare lost, the level of intraspecific competition should not be affected: the coefficients on the\ndiagonal do not change when the matrix is pruned. We have clarified that this is indeed a\nreasonable assumption in the main text:\n“So far, we have investigated how coexistence is affected by increasing α. However,\nintraspecific competition models the effects of a species on itself. As such, when\nspecies are lost to extinction, the level of intraspecific competition should not be\n7\naffected: the coefficients on the diagonal do not change when the matrix is pruned\nof the populations that have gone extinct.”\n• Line 206: why “proportional to each other”? In the equation, y⋆ also occur in\nthe dominator, making it nonlinear.\nWe apologize for the confusion. What we mean is that if we were to plot x⋆ against y⋆, we\nwould obtain points laying on a line with intercept zero, and slope 1/(1+mTy⋆) (a positive\nscalar); hence, x⋆ and y⋆ would have the same signs (feasibility is unchanged), and would be\nproportional to each other. We have reworded the paragraph to make this clearer.\n• Fig. 1: This is a very clear and interesting illustration. Could such transitions\noccur already in two-species Lotka–Volterra models?\nThis is a very good point, and we have added a figure (Fig. 4 here) to the SI with the\nequivalent scheme for the simple two species GLV competitive model. In this figure, we also\nshow the bifurcation diagram for the parameter α, showing how the system traverses the\nthree phases: alternative stable states, convergence to a stable equilibrium with extinctions,\nand robust coexistence. We have added this figure to the beginning of the supplementary\ninformation, and reference it in the caption of Figure 1, and the main text.\nReviewer 3\nI found the work presented in ‘Robust coexistence in ecological competitive commu-\nnities’ very interesting, elegantly presented and quite exciting frankly. The authors\ndemonstrate the significance of self-regulation as a mechanism for community per-\nsistence. In particular, the authors use random matrices to generate competitive\necological communities and show that increasing self-regulation (intra-specific compe-\ntition) first achieves stability (fixed point equilibrium) and then feasibility (positive\ndensity of species at equilibrium). The order of stability being attained prior to feasi-\nbility (along the self-regulation gradient) implies that feasibility guarantees stability.\nAn important corollary of this work is that competitive communities (that are large\nand have randomly distributed interaction strengths, a fundamental assumption and\ncaveat) will have stable equilibria rather than non-equilibrium steady states.\nWe thank the Reviewer for acknowledging the relevance of our results in their context, together\nwith their positive evaluation of our work. Their suggestions were very helpful and, by addressing\nthem point by point, the presentation of the manuscript improved and, overall, the contribution\nbecame strengthened.\nMajor comments\n1. I have one comment which touches upon the basic assumptions of this work. The\nauthors use a single value for self-regulation, α, across species. This scenario\nwould imply that all species share the same strength of self-regulation, while the\nstrength of interspecific interactions can vary. I wonder what would happen if\nthere was a mean α for a matrix, rather than a fixed α for all species? How\nstrongly would this impact the insights? Adding to this, what would happen if\nthese α values were somehow correlated to the interspecific interaction terms?\n8\n2\n1.5\n1\n0.5\n0\n0 50 100\nTime\necnadnubA\nPhaseI(α=0.55) PhaseII(α=0.60) PhaseIII(α=0.75)\n1.5 1\n1\n0.5\n0.5\n0\n0 50 100 0 50 100\nTime Time\n2\n1.5\n0.5\n0\n0.5 0.6 0.7\nIntraspecificcompetitionα\necnadnubadnE\nx x\n1 2\nIC1 IC2\nFigure 4: Three dynamical phases in the classic two–species Lotka–Volterra competi-\ntion model. Top row: time series for α α ,α ,α (here, 0.55,0.60,0.75), with two initial\nI II III\n∈ { }\nconditions shown (linetypes; colors denote species). Bottom row: bifurcation diagram of equilib-\nrium abundances versus the intraspecific competition strength α (x–axis), with line type indicating\nthe initial condition. As α increases while a ,a are fixed and r = r = 1, the system transitions\n12 21 1 2\nthrough three regimes separated by the vertical guides: Phase I (low α): interspecific effects\ndominate and priority effects appear—different initial conditions lead to alternative stable out-\ncomes (each species can exclude the other). Phase II (intermediate α): competitive exclusion\nis deterministic—the same species survives from any initial condition (unique attracting boundary\nequilibrium). Phase III (high α): intraspecific competition dominates and a feasible interior\ncoexistence equilibrium exists and attracts all initial conditions. All three panels use the same\ninterspecific coefficients, varying only α across scenarios.\nThe Reviewer suggests a very reasonable extension – allowing intraspecific interactions to\nvary across species. However, it substantially complicates the analysis. For example, the\nDyson (or self-consistent) equation used to characterize spectra and resolvents holds cleanly\nunder independent and identically distributed assumptions, whereas heteroscedastic diagonal\nentries introduce a variance profile that breaks those simplifications. Although there has\nbeen significant progress for stability questions with variance profiles (for instance, Cook et\nal., Electron. J. Probab. 2018), feasibility requires controlling the distribution of equilibrium\ncomponents, which is more delicate.\nWe therefore adopt a controlled relaxation in which diagonal and off–diagonal variances can\ndiffer modestly. Concretely, let\nVar(B ) = σ2, Var(B ) = 1 for i = j, Cor(B ,B ) = ρ,\nii α ij ij ji\n̸\nwith σ2 close to 1. Under this change we derive an analytical expression for the feasibility\nα\n9\nprobability\np (α) = P(α α),\nF F\n≤\nand find that the entire curve p (α) shifts slightly to the right relative to the σ2 = 1 case.\nF α\nTheshiftissmallforσ2 1(seeSISec.3.8). Forglobalstability, resultsfromrandommatrix\nα\n≈\ntheory with variance profiles justify that, in the limit n , the stability threshold α is\nS\n→ ∞\nunchanged by this modification, so the ordering α < α persists asymptotically.\nS F\nWe also test robustness under strong deviations by simulation, allowing very large diagonal\nvariance, up to σ2 = 100. Consistent with the analytical trends, even when relaxing the\nα\nequal self–regulation assumption widely, our main conclusion holds: the stability transition\nprecedesthefeasibilitytransition. Thefigurebelowdisplaysempiricalcumulativedistribution\nfunctions (CDFs) of feasibility and stability thresholds across many random communities.\nSolid curves are the CDFs of α (probability of feasibility), dashed curves are the CDFs of\nF\nα (probability of stability). Color encodes the diagonal variance level σ2. At each α on the\nS α\nhorizontal axis, the vertical coordinate gives p (α) = P(α α) or p (α) = P(α α).\nF F S S\n≤ ≤\nWe use n = 200, correlation ρ = 0, off–diagonal entries drawn from a standard normal, and\nN = 1000 random matrices per σ2. Curves farther to the right correspond to larger typical\nα\nthresholds. Within each color, the dashed curve lies to the left of the solid one, indicating\nα is typically smaller than α . Increasing σ2 shifts both distributions to the right, implying\nS F α\nhigher required α.\nFigure 5: Empirical CDFs of feasibility and stability thresholds under heterogeneous\ndiagonal variance. Solid curves: p (α) = P(α α). Dashed curves: P (α) = P(α α).\nF F S S\n≤ ≤\nColors indicate σ2 0,100 . Parameters: n = 200; off–diagonal entries B for i = j are\nα ij\n∈ { } ̸\ni.i.d. (0,1) with Cor(B ,B ) = 0; diagonal entries B (0,σ2); N = 1000 matrices per σ2.\nij ji ii α α\nN ∼ N\nFor each color, the dashed curve lies to the left of the corresponding solid curve, demonstrating\nthat the stability transition precedes the feasibility transition for each σ level.\nα\n2. My only regret was the absence of an ecological interpretation for the results,\nbeyond the obvious statements. I see some relation to niche theory; it would be\ninsightful to link the two.\nthis deliberately, being aware that we are studying a quite contrived/unrealistic parameteri-\nzation, in which all populations are essentially equivalent. We recognize that natural systems\nwould not abide by these strict requirements, and that therefore the direct application of\nthese results to natural/experimental systems would be premature.\nOur goal is more modest, but in our opinion still essential to make progress in the discipline.\nWe want to provide a baseline expectation for the relationship between feasibility and stabil-\nity; to this end, we analyze simple systems that can be studied mathematically — thereby\nidentifying a few of the main levers regulating the behavior of the community, and deriving\nformulas that can be then contrasted with experimental and observational data.\nAs it is often the case in science, what is interesting is to understand when, how and why\ndo natural systems deviate from these expectations. We can turn to May, who in his 1973\nbook prompted ecologists to “elucidate the devious strategies which make for stability in\nenduring natural systems”—escaping the relationship between complexity and stability he\nhad derived. It is very much in this spirit that we presented a minimal discussion, setting out\nclear expectations for the behavior of random, unrealistic systems, with the hope this will\nhelp elucidate how natural systems behave.\nMinor comments\n• Title: Would not “Robust coexistence in competitive ecological communities”\nsound better?\nYes, we have updated the title accordingly.\n• Abstract: The final sentence is quite strong. The absence of non-equilibrium\nsteady states rests on certain assumptions. Consider toning down the absolute-\nness and expand on robust coexistence.\nWe have restated our assumptions in the abstract two more times, to make sure the reader\nunderstands that the results presented apply to large, random, competitive communities.\nSection 1\n• line 53, clarify the notion of universal convergence to feasibility at this point so\nthat it becomes easier to follow this line of reasoning\nWe have added an intuitive definition of universality:\n“They also are crucial for coexistence: for example, May showed that ecological\ndynamicscanbestabilizedbyimposingsufficientlystrongintraspecificcompetition,\nand that for large systems with random species interactions the stability transition\nisuniversal(independentoftheprobabilitydistributionusedtomodelinteractions)\nand sharp (always achieved at the same critical value of intraspecific competition)”\nWe cannot be more precise with our wording before presenting the more technical results in\nthe following section.\n• line 60 ‘resilient to any perturbation sort of extinction’: this statement is not\nvery clear\nment – that the perturbations we refer to are those involving changes in species abundances,\nnot community composition.\nSection 2\n• In Fig. 1, add a description of the top panel to make it immediately obvious\nthat this refers to the interaction matrix entries, with the shade reflecting the\nmagnitude of the parameters.\nWe added:\n“The three matrices, with the same off-diagonal coefficients and increasingly strong\ndiagonal coefficients, are shown as tile plots above the arrow (darker colors reflect\nlarger values).”\n• line 83 ‘neither feasible’\nWe have updated the text\nSection 3\n• Lines 124–127 and Eq. (4): it appears that as n increases, only the correlation\nρ will impact the variance and hence the probability of feasibility. If that is\nthe case, I am not sure I understand line 125, where I perceive the message as\nbeing that correlation (besides the mean and variance of the interaction matrix)\nis irrelevant.\nApologies for the misunderstanding. The Reviewer is correct that the correlation will indeed\nimpact the variance even when n is large. In the main text (n = 1000), the correlation has a\nclear effect on the probability of feasibility.\nThere was a small mistake in the text, which we have corrected: the mean of the distribu-\ntion (µ) does not affect probability of feasibility (see Methods)—only the variance and the\ncorrelation do. It is true that in Eq. (4) the variance does not appear explicitly, because it\nis absorbed into the choice of the excess intraspecific competition α through rescaling, as we\nnow show in the Methods of the paper. Thus, we set the variance to 1, and focus solely on\nthe effect of α. Naturally, one could have fixed α and varied σ instead. We kept the variance\nfixed and focused on α to keep with the classic formulation introduced by May in his seminal\nwork on stability of ecological systems.\nInsum,correlationρandvarianceσ2 (or,equivalently,intraspecificcompetitionα)arecrucial\nfordeterminingthevarianceoftheequilibriumy(σ2),andthereforecomputingtheprobability\ny\nof feasibility. However, once we fix ρ and σ2 (equivalently, α), the probability of feasibility\nwill be the same regardless of what distribution we use, provided that n is large. We have\nupdated the main text to correct this:\n“Eq. 4 shows that the details of the distribution of interactions (Fig. 3, linetypes),\nonce accounted for the correlation (Fig. 3, colors), are irrelevant when the system\nis large.”\nα.\n• Line 131: for example, what does p = 0.5 mean for a community?\nF\nThisshouldbeinterpretedprobabilistically: whensettingαsuchthatp = 0.5,sampling1000\nF\ninteraction matrices would result in approximately 500 of them having a feasible equilibrium.\nAs such, if one wanted to always find a community with feasible equilibrium, α would need to\nbe infinite. For any smaller value there will always be a (small) chance to find an unfeasible\ncommunity.\n• Are there other aspects of the interaction matrix distribution (higher moments)\nwhose consideration will make this outcome more definite?\nThis is a very interesting point, which has a parallel in the study of stability. For stability,\nprovided that interactions are sampled independently (or pairs of interactions are sampled\nindependently), higher moments of the distribution of the coefficients (that is, skewness,\nkurtosis, etc.) have no effect when the system is large.\nIn Appendix B1 we show that the same holds for feasibility: when we derive the variance of\ny∗ (which controls the probability of feasibility) using the Richardson approximation, we see\nthat higher moments of the distribution do appear in the formulas. However, the coefficients\nassociated with the higher moments vanish as n grows, showing that indeed asymptotically\nonly the mean, variance, and correlation will have an effect. This is key for the universality\nof the results.\nThe study of the higher moments of the distribution was the main reason that led us to con-\nsider non-symmetric distributions, such as the Modified Beta and Skew-Normal distributions.\nProbability density functions for three skew values of each distribution are shown in the first\nrow of Fig. 6.\nInthesecondrowofFig.6weshowtheprobabilityoffeasibilityasafunctionofαfordifferent\nvalues of skewness and different distributions. For small systems, the higher-order terms have\na small but noticeable effect.\n• Is feasibility higher dimensional than stability and do we require additional di-\nmensions to make its transition sharp?\nThis is an interesting question. Simulations and analytical results show that the transition\nfrom low to high probability of stability is very sharp when considering stability. In fact, this\nwas the main result of the famous paper by Robert May on stability:\n“The central feature of the results for large systems is the very sharp transition\nfrom stable to unstable behaviour as the complexity exceeds a critical value” (May,\n1972).\nWhat we observe for feasibility is more nuanced. First of all, for both stability and feasibility\none needs to consider a scaled version of the system — the eigenvalues of an infinitely-sized\nmatrixwoulddivergetoinfinity,andthisiswhyonealwaysanalyzesanappropriatelyrescaled\nmatrix, A/√n. Even when rescaling the entries of the matrix, we find that the threshold\nfor feasibility keeps growing with n (while that for stability converges to a constant, as\nshown in Fig. 2c of the manuscript, Fig 1c here). This is consistent with the fact that the\nthreshold for feasibility is higher than that for stability. The width of the sigmoidal part\n13\n0.4\n0.2\n0\nytilibaborP\nSkew normal Modified beta Discrete\n1\n0.5\n0\n60 100 140\nP\nytilibisaef\nfo\nytilibaborP\nF\n60 100 140 60 100 140\nIntraspecific competition α\nLow Neutral High\nFigure 6: Effect of skewness of interaction distributions on feasibility transitions. Top\nrow: probability density functions of skew-normal. modified beta, and discrete distributions (as in\nthe manuscript) and their skewed versions in both directions (colors). Bottom row: probability of\nfeasibilityp (α)asafunctionofintraspecificcompetitionα fordifferentskewnessvalues. Forsmall\nF\nn, higher moments such as skewness shift the curves slightly, but for large n the effect vanishes (see\nAppendix B1), confirming universality of the feasibility transition.\nof the curve describing the probability of feasibility shrinks as n increases; however, we find\nthat to observe a sharp-enough transition we would need to consider an unrealistically large\nnumberofinteractingpopulations. Forinstance, theEarthisbelievedtosupportontheorder\nof ten million species, and we would need to consider trillions to have an appreciably sharp\ntransition. Hence, the sigmoidal curve describing feasibility remains noticeably sigmoidal for\nany reasonable community size.\n• Lines 142–144: If p < 1, does this imply the community will ‘prune’ itself to a\nF\nlower n and attain p = 1?\nF\nThisisanexcellentpoint, whichwetriedtocapturewhendiscussingtherelationshipbetween\ncommunity assembly and the probability of feasibility (Section 4). In general, p is always\nF\nless than 1 (it would only be equal to 1 for α ). Given that feasibility is necessary for\n→ ∞\ncoexistence, and that the only possible outcome for an unfeasible competitive system is the\nextinctionofoneormorespecies, wenecessarilyhavethatwheneverthesystemisnotfeasible\n(which happens more often when p is low), the system will be pruned, resulting in a smaller\nF\nsystem that is feasible.\nNote that this does not mean that p = 1 — just that the system was pruned to a point at\nF\nwhich feasibility was achieved (α > α ). An example clarifies this subtle point. Suppose that\nF\ninitially n = 1000 and p = 0.01. Chances are that a random system would be unfeasible,\nF\nandthusthe firstspecies wouldbe lost. Nown = 999andp would increaseslightly(because\nF\n14\nof the decrease in n). Again, it would still be improbable to have feasibility (we would still\nhave a probability of about 1 in 100). Another species would be lost, and so on. Note that\nevery time a species is lost, we have a small chance of having achieved feasibility, such that\nafter a large number of extinctions, we could have found a feasible system, even though for\nmatrices with the same size and statistics the probability would still be p 1.\nF\n≪\n• It would be interesting to seek a pattern linking α to the ‘pruned’ community\nsize n for which p = 1.\nF\nThis is a great, but subtle point. A combination of Figures in the main text can answer\nthis question (Figs. 2c and 4b). Figure 4b in the main text shows that after the community\nprunes, α α. That is, the assembly of a community can be seen as the decrease in α due\nF F\n≈\nto extinctions until the system reaches the first community composition with α below the\nF\ncommunity intraspecific competition α. Figure 2c shows what is the tipical value of α for\nF\neach community size. Together, they provide an answer to the pattern the Reviewer is asking\nfor. Typically, the number of species at which assembly will stop will be the one given when\nwe set α = α and look at the x-axis for the corresponding solid curve shown in Fig. 2c.\nF\n• Does very large n lead to a sharp transition?\nThis is a very interesting question. The new Fig. 2 aims to address part of it. Specifically\npanel (c), where we calculate α for many communities of increasing size, and compare it\nF\nwith the results for α . We observe that the latter grows asymptotically with n, to an upper\nS\nbound. On the contrary, the growth of the former is faster, and does not seem to have\nan asymptote for the values of n studied. Indeed, we confirm that this is the case in the\nsupplement (SI, Sec. 3.8, page 48), where we show that:\n“Inthelimitn , thewidthofthedistributiontendstozero, thusthetransition\n→ ∞\nbecomes sharper and sharper, and displaces to the right without reaching a finite\ntransition value in the limit, because the transition occurs at α (the median grows\nF\nunboundedly). It is important to remark that the convergence of α to a step\nF\nfunction centered at infinity is very slow, as Fig. 14 (SI) shows.”\nWe point to this in the main text too:\n“Additionally, the transition to feasibility also remains smooth compared to stabil-\nity (Fig. 3, and SI Fig. 14, see Sec. 3.8 for a detailed discussion) for large systems.”\n• Fig. 2: the correlation in the interaction matrix can be seen as a proxy for\nspecialisation. Positively correlated entries imply no trade-offs, while negative\ncorrelation indicates trade-offs and hence some kind of specialisation. This can\nbe linked to the notion of the niche.\nIndeed, as is often the case in community ecology, our results can be interpreted as a “niche-\nbasedmechanism”: strongintraspecificcompetitionsignalsaccesstoresourcesthatareprivate\nto the species (and hence a separation of niches), while a positive correlation makes the\npopulations more similar to each other compared to the anticorrelated case. We now mention\nthese facts in the discussion.\n“Intraspecific competition models the fact that each species has access to private\nresources or regulating factors that are species-specific; as such, it plays an impor-\ntant role in niche differentiation, and sets a limit to niche similarity. Our results\n15\nshow that, to ensure feasibility, one had to postulate a level of niche differentiation\nthat is sufficient to guarantee stability. This observation reinforces the notion that,\nin competitive systems, feasibility is key for coexistence. This is in analogy with\nconsumer-resource models, as well as in GLV models with symmetric interactions\nproduced by niche overlap or species’ similarity, for which stability is guaranteed,\nand thus feasibility determines coexistence. ”\n16"}, {"round": 2, "reviewer_comments": {"reviewer_1": "This revision substantially improves both rigor and readability. In particular, the revamped treatment around Fig. 2 now clearly articulates the relationships among feasibility and stability thresholds, contrasts the finite-n “multi-window” phenomenon with the large-n single transition, and maps these mathematical results onto ecological regimes (priority effects → deterministic exclusion → stable coexistence). The Discussion’s ecological interpretation is now well aligned with niche theory.\nWe have only one point to consider before publication. The extension to random diagonal terms (heterogeneous intraspecific competition) and the claim that the main qualitative transition persists for large n strengthen ecological realism. However, it remains unclear whether the number of transitions into and out of feasibility converges to one as community size increases when intraspecific competition is heterogeneous. If this “single-window” property is not robust to such heterogeneity (which should be the case in reality), conclusions that rely on a unique α_F may be fragile.\nMinor comments\nL140: “with [increasing] alpha”\nLine 621: “we we have” → “we have”.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "I appreciate the amount of work the authors have put it to thoroughly answer all the question in the first round of reviews. I continue to find the discussion too light on ecological context, though I understand the authors’ perspective. I only have some minor comments, below.\nl. 118-121: interesting how inter-specific competition strength does not impact feasibility. So, the persistence/feasibility of a community depends on the strength of self-regulation of the individual populations, i.e., if self-regulation is ‘strong enough’, a community will persist no matter the inter-specific interaction strengths. As discussed in l. 193-199, this relates to niche differentiation. Can we say that, sufficient (or appropriate) niche differentiation in a competitive community guarantees the stable coexistence of all populations? Could we begin to make connections to other metrics of niche differentiation and potentially empirically test these findings? I would have liked the authors to discuss (speculate even) such aspects of their work.\nl. 150: ‘almost surely guaranteed’. In the context of this modelling framework, is it not guaranteed?\nl. 191: again ‘almost surely’. I realise this may be in response to comment 2 of reviewers 1 & 2 (with which I agreed), but I am not sure that is what they meant. If a_f > a_S holds always (for large n and the given set-up with its specific assumptions), then state it. Otherwise, it causes confusion. The authors do this clearly in l.200-201."}, "author_response": "Responses to Reviewers (round 2)\nManuscript NCOMMS-25-39378\nReviewers 1 & 2\nThis revision substantially improves both rigor and readability. In particular, the\nrevamped treatment around Fig. 2 now clearly articulates the relationships among\nfeasibility and stability thresholds, contrasts the finite-n “multi-window” phenomenon\nwith the large-n single transition, and maps these mathematical results onto ecolog-\nical regimes (priority effects → deterministic exclusion → stable coexistence). The\nDiscussion’s ecological interpretation is now well aligned with niche theory.\nWe thank the reviewers for the careful consideration of our manuscript, and for engaging with\nsuggestions of Reviewer 3 in a constructive way.\nMajor comments\n1. We have only one point to consider before publication. The extension to ran-\ndom diagonal terms (heterogeneous intraspecific competition) and the claim that\nthe main qualitative transition persists for large n strengthen ecological realism.\nHowever, it remains unclear whether the number of transitions into and out of\nfeasibility converges to one as community size increases when intraspecific com-\npetition is heterogeneous. If this “single-window” property is not robust to such\nheterogeneity (which should be the case in reality), conclusions that rely on a\nunique α may be fragile.\nF\nIndeed, without showing that there is one single transition for feasibility (asymptotically) for\nthe case of heterogeneous variance, the argument we made on our previous response to justify\nthat heterogeneous intraspecific competition does not change the results, is incomplete. As\nsuch, we run simulations where the matrix of interactions is sampled as\nA = αI +D(a)+B,\nwhere B (i ̸= j) are i.i.d. random variables with zero mean and variance 1, D(x) stands for\nij\na diagonal matrix with vector x in the diagonal, and a = (a ) is a vector of i.i.d. random\nk\nvariables with the same distribution as B but zero mean and variance σ2.\nij α\nIn the figure below, we plot the average number of feasible intervals for increasing values\nof diagonal variance σ2. We see that increasing such variance only decreases the number\nα\nof transitions, but the average number of transitions converges to one at approximately the\nsame rate. Therefore, asymptotically, the picture is the same: there is only one transition to\nfeasibility.\n1\n1.2\n1.1\n1.0\n0.9\n0.8\n5 10 15 20\nSystem size (n)\nslavretni\nelbisaef\nrebmun\negarevA\nAverage number of feasible intervals vs system size\ns 2\na\n5\n4\n3\n2\n1\nvergence to a single transition to feasibility. Average number of feasible intervals as a\nfunction of system size (n) across different variance levels (σ2). Each line represents a different\nα\ndiagonal variance parameter ranging from σ2 = 0 (purple, deterministic diagonal) to σ2 = 5 (yel-\nα α\nlow, high variance). Results are averaged over 1000 simulations for each (n,σ2) combination. The\nα\nmatrix of interactions is formed by normally distributed off-diagonal elements (without correlation\norskewness),withdiagonalelementsdrawnfromthespecifiedvariancedistribution. Thenumberof\nfeasible intervals represents the count of parameter regions where the system maintains feasibility.\n• L140: “with [increasing] alpha”\nWe have corrected this, thanks.\n• Line 621: “we we have” → “we have”.\nCorrected, thanks again.\nReviewer 3\nI appreciate the amount of work the authors have put it to thoroughly answer all the\nquestion in the first round of reviews. I continue to find the discussion too light on\necological context, though I understand the authors’ perspective. I only have some\nminor comments, below.\nWe thank the Reviewer for their positive evaluation of our work. We appreciate their under-\nstanding of our perspective in relation to the discussion.\nMinor comments\n• l. 118-121: interesting how inter-specific competition strength does not impact\nfeasibility.\nYes, this is correct. Feasibility is not affected by strength of competition. However, the value\nof the equilibrium abundances would decrease with increasing competitive strength. As such,\nif we were to consider stochastic effects(as is the case for real ecosystems), low abundances\nmay lead to extinctions. For the results presented here however, inter-specific competition\nstrength is irrelevant. We have added a sentence to the methods stating\n“Although increasing the strength of competition µ does not affect feasibility, it decreases\nequilibrium abundance values. As such, increasing µ too much will result in very low abun-\ndances, which could lead to extinctions in real systems where stochasticity plays a factor.”\n• So, the persistence/feasibility of a community depends on the strength of self-\nregulation of the individual populations, i.e., if self-regulation is ‘strong enough’,\na community will persist no matter the inter-specific interaction strengths. As\ndiscussed in l. 193-199, this relates to niche differentiation. Can we say that,\nsufficient (or appropriate) niche differentiation in a competitive community guar-\nantees the stable coexistence of all populations?\nThis is also is correct, intraspecific competition is a proxy for niche differentiation. As such,\nthe Reviewer’s statement is consistent with the results of this work.\n• Could we begin to make connections to other metrics of niche differentiation and\npotentially empirically test these findings? I would have liked the authors to\ndiscuss (speculate even) such aspects of their work.\nThis is indeed a very exciting and promising possibility. We have added a sentence in the\ndiscussion stating\n“While direct empirical test of these patterns is out of the scope of this work, it is an exciting\navenue for future research.”\nprovide a null expectation of the behavior of random communities, which we acknowledge to\nbe not representative of communities in the real world. The spirit is similar to that pioneered\nby May.\n“Thatstabilitymayusuallygowithcomplexityinthenaturalworld,butnotusually\ningeneralmathematicalmodelsisnotreallyparadoxical. Innaturewedealnotwith\narbitrary complex systems, but rather with ones selected by a long and intricate\nprocess. The emergent moral is that theoretical work should not try to prove any\ngeneral theorem that “complexity implies stability”, but instead should focus on\nelucidating the very special sorts of complexity, the singular strategies, which may\npromote such mathematically atypical stability.”\nHere, we aim to parallel his work on stability, outlining the typical behavior of the feasibility\nofeasibility\nof large, random, competitive communities.\n• l. 150: ‘almost surely guaranteed’. In the context of this modelling framework,\nis it not guaranteed?\nWe agree that, written as such, the statement is rather technical. We have re-written it as\n“a feasible system will be stable with probability 1”, to avoid confusion.\n• l. 191: again ‘almost surely’. I realise this may be in response to comment 2\nof reviewers 1 & 2 (with which I agreed), but I am not sure that is what they\nmeant. If α > α holds always (for large n and the given set-up with its specific\nF S\nassumptions), then state it. Otherwise, it causes confusion. The authors do this\nclearly in l.200-201.\nWe agree that, since the assumptions of the modeling framework are stated in the same\nsentence, we can use the word ‘always’ as a synonym for “with probability 1”, rather than\n“almost surely”, which is more technical. As such, we have changed the wording to: “the\ntransition to stability always precedes that for feasibility”. Note that the results presented\nhere are only true in probability, in the sense that, for large random matrices, the probability\nto sample a non-stable feasible interaction system is never zero, but vanishingly small, all\nthe statements about frequency in this paper should be understood under this lens. We have\nadded a sentence in the discussion stating this\n“Importantly,alltheresultsinthispaperaretrueinprobability. Thatis,whiletheprobability\nof finding feasible but unstable systems is never exactly zero, it becomes vanishingly small\nwith system size.”\n4"}]

Decision letter