Stability of ecosystems enhanced by species-interaction constraints

Ecosystem stability is a central question both in theoretical and applied biology. Dynamical systems theory can be used to analyze how growth rates, carrying capacities, and patterns of species interactions affect the stability of an ecosystem. The response to increasing complexity has been extensively studied and the general conclusion is that there is a limit. While there is a complexity limit to stability at which global destabilisation occurs, the collapse rarely happens suddenly if a system is fully viable (no species is extinct). In fact, when complexity is successively increased, we ﬁnd that the generic response is to go through multiple single-species extinctions before a global collapse. In this paper we demonstrate this ﬁnding via both numerical simulations and elaborations of theoretical predictions. We explore more biological interaction patterns, and, perhaps most importantly, we show that constrained interaction structures—a constant row sum in the interaction matrix—prevent extinctions from occurring. This makes an ecosystem more robust in terms of allowed complexity, but it also means singles-species extinctions do not precede or signal collapse—a drastically different behavior compared to the generic and commonly assumed case. We further argue that this constrained interaction structure—limiting the total interactions for each species—is biologically plausible.


I. INTRODUCTION
In theoretical studies of ecosystem stability, dynamical systems are often used. These models were initially extended from a few species [1,2] to whole ecosystems by using random matrices to represent the interaction network among species. In particular, May found that increased complexity in terms of number of species and/or interactions typically leads to instability-a surprising result at the time [3]. Since then, more powerful numerical methods have led to additional insights into ecosystem stability by using approaches such as dynamical modeling with higher-order or other forms of species interactions [4][5][6][7][8] as well as topological studies [9][10][11][12]. Although higher-order interactions and features such as adaptive foraging have been shown to sometimes reverse the generic conclusions first found by May [13,14]complexity leads to system-wide instabilities-many studies confirm it [6][7][8].
The simplest and most widely used dynamical model capturing the interaction complexity, and the most closely related to May's result, is the generalized Lotka-Volterra (GLV) Published where the x i are species abundances for species i = 1, 2, . . . , N, and r i and K i are intrinsic growth rates and carrying capacities for each species respectively. A is the N × N interaction matrix and σ a parameter that regulates the standard deviation of the interaction strengths. The diagonal of A is set to zero (A ii = 0), because the first term on the right side of Eq. (1) already captures the self-interactions. The fixed points of Eq. (1) are Until recently, a common approach when analyzing stability of the GLV has been to assume that all species are nonextinct [2,3,[15][16][17][18]. This restricts the analysis to the nontrivial fixed point (called a feasible fixed point) of Eq. (2), and effectively rules out extinctions. To analyze the stability of the model the Jacobian evaluated at a fixed point (J * ), including nonfeasible fixed points, is used. When taking extinctions into account it can be written as Here X * , K, R, and D are diagonal matrices with x * i , r i , 1/K i , and (Ax * ) i on the diagonal respectively. E is a diagonal matrix with ones corresponding to extinct species, E ii = 1 when x * i = 0, and zero otherwise. Note that rows corresponding to extinct species in J > are zero, and equivalently rows for nonextinct species in J 0 . This means the spectrum of J * has two separate parts {λ 1 , . . . , λ n , γ 1 , . . . , γ m } where n and m are the numbers of no-extinct and extinct species respectively.
The Jacobian used by May, J M = (σ A − I ), is similar to J > or the entire Jacobian under the assumption that all species are nonextinct, when r i = K i = 1. Assuming feasibility by using J M , drawing entries of A from a random normal distribution N (μ = 0, 1) with probability c (connectance), and using Wigner's semicircle law [19], May derived a boundary for fixed-point stability which was interpreted as a boundary for ecosystem stability in terms of σ , N, and c. The equation for the boundary is i.e., the point when we expect the Jacobian to have an eigenvalue with positive real part. This result has ever since been extensively debated and its relation to real ecosystems both questioned and defended [20]. Many studies have used it to assess stability of real and simulated ecosystems [17,21,22]. It has also been extended, for example, in terms of modularity in interactions and mean interaction strength [18,23]. The Jacobian introduced by May and its extensions do not include fixed-point abundances. The Jacobian's sign structure (hence the stability) can in some cases be affected by the X * scaling [24]. However, the main qualitative difference is the possibility of extinctions in the dynamical model. Extinctions occurring at lower complexities (smaller values of σ ) than for destabilization, σ M = 1/ √ cN, will shift the stability boundary to larger σ , a feature that has long been overlooked. In our previously published work, we showed that this in fact is the generic behavior of systems with random interaction structures, we updated the complexity limit for destabilization/collapse, and we introduced a limit beyond which systems can no longer be feasible [25]. This picture substantially changes both the prediction of a system's response to perturbations as well as approach to collapse. Figure 1 demonstrates the general behavior of the GLV model with random interactions, as analyzed in [25], when no restrictions on feasibility are imposed. Although for sufficiently large σ there may be several fixed points that are qualitatively equivalent, the figure shows a specific but generic trajectory that is followed when increasing σ . In particular, systems do not switch between fixed points following perturbations in σ or in species abundances. However, randomly chosen initial species abundances for the dynamical system could induce a switch. Taken together, this means that, after a perturbation in σ , the dynamical system can respond in two qualitatively different ways: either a single species goes extinct or the system loses stability.
With the generic behavior of random systems mapped in [25], in this paper we follow up by answering the question, are there any biologically relevant conditions for which extinctions will not occur at lower complexities than the global collapse? To answer this, we begin with an elaboration of our previous derivation of the extinction boundary.

II. DERIVING THE EXTINCTION BOUNDARY
As seen in Fig. 1 species extinction occur for smaller sigma than May's stability boundary. To show that this behavior is generic we need to calculate the expected value of σ at the first extinction, σ f , to compare with the stability boundary in Eq. (5). The boundary σ f is the extinction boundary introduced in [25], then as a first-order approximation. Here, we go even further and derive an exact expression.
Assuming no extinction has occurred, the fixed-point equation is linear: where we set r i = K i = 1 for simplicity (r i = 1 and K i = 1, demonstrated in [25]). Reexpressing the inverse (I − σ A) −1 as a von Neumann series expansion gives For guaranteed convergence of the von Neumann expansion ||σ A|| ∞ < 1 (or equivalent norm). This strict condition does not hold for σ A in the σ ranges of interest here, although there are less strict criteria of convergence and rigorous numerical investigations confirm the sum typically converges until in the vicinity of the stability boundary, σ M . From Eq. (7) it follows that x * i are sums of products of the random variables A i j , and from now on we treat them as stochastic variables, X * i . For simplicity we start by assuming that the mean of interaction strengths, μ, is zero. The general case with μ = 0 is treated below. Assuming independence between the terms in the outer sum (scalar product with the vector 1) in Eq. (7), the X * i will converge to the normal distribution for sufficiently large cN in accordance with the central limit theorem. Then the variance of the X * i can be found using the rules for mean and variance of product and sum distributions for independent stochastic variables, for example Finally the outer sum in Eq. (7) gives the total variance of X * i , i.e., . It follows that for μ = 0 the species abundances have a normal distribution with mean 1, variance 0 at σ = 0, and diverging variance at May's stability limit Eq. (5). The divergence is however artificial since when increasing σ one of the x * i will become zero, breaking the assumptions in Eq. (6), marking the first extinction event. Statistically we can express this event in terms of order statistics where the distribution of the qth smallest species abundance of X * i [with distributions f (x)] can be expressed as The first extinction boundary, σ f , is the σ for which Eq. (10) has mean zero.
To account for systems with nonzero mean interaction strengths, μ = 0, the interaction matrix can be written as 1), and M is a matrix with μ at the same positions as the nonzero entries of A. Then for (I − σ A μ=0 + σ M )x * = 1 and c = 1, Mx * is approximately a constant vector x T OT 1. Each of these sums exclude the diagonal element but is approximately correct when the variance of the species abundances is small. Therefore FIG. 2. Extinctions or destabilization first? The top plot shows simulation averages of the mean value (with one-standard-deviation error bars) of the variation in interaction strength, σ , at which first extinction occurs. The simulation averages are from "Random" systems with interaction strengths from a normal distribution A i j ∼ N (0, 1), predator/prey systems where sgn(A i j ) = − sgn(A ji ) (ρ = −2/π ), mutualistic/competitive systems with sgn(A i j ) = sgn(A ji ) (ρ = 2/π ), and random systems with different means of the normal distribution A i j ∼ N (μ = ±0.1, 1). All simulations were implemented with r i = k i = 1. The dashed lines are theoretical first extinction boundaries, σ f . The dash-dotted lines are May's stability boundaries, σ M , color-coded according to the legend (some colors are "missing" since they coincide; for example all first extinction boundaries coincide and collapse boundaries for systems with varying interaction strength mean, for these only "Random" is visible). May's stability boundary is given by [11]. Note the spread in σ M compared to σ f , which is equal for all systems. The bottom panels show an example simulation of each type, with the same color coding. The vertical lines indicate σ f , σ M , and the actual loss of stability σ c (operational collapse definition), in order of increasing σ . Three phases are marked by the different shades: strict stability (SS) before σ f , extinction continuum (EC), a phase of single-species extinctions, and collapse (C) where no stable nearby fixed-point exists. Note that σ M always falls in the extinction continuum for every type of system. and we conclude that μ = 0 only adds to a multiplicative scaling factor, that does not affect the first extinction event. This result is verified in Fig. 2.
The above analysis can be extended to interaction-webs that capture more biologically inspired interaction structures 062405-3 found in the literature, for example predator/prey, sgn(A i j ) = − sgn(A ji ) [15] and mutualistic/competitive, sgn(A i j ) = sgn(A ji ) [11]. The symmetric or antisymmetric property of the interaction matrix can be modeled by a parameter ρ, where ρ = 1 gives a symmetric matrix and ρ = −1 an antisymmetric one. The above predator/prey and mutualistic/competitive interaction structures then correspond to values of ρ = −2/π and ρ = 2/π respectively [15]. Even for these structures, entries generally do not match when producing elements of A p in Eq. (7). Thus, the outer sum can be assumed to be a sum of independent variables. This implies no substantial change in the species abundance distributions, so we expect the first extinction event, as in the case with μ = 0, to be little affected by these structures. A comparison between first extinction and the stability boundary shown in Fig. 2 confirms this and our assumption of independence.

III. DO EXTINCTIONS HAPPEN BEFORE COLLAPSE?
From our analysis and Fig. 2 we reconfirm that extinctions are generic in the GLV with random interactions, as well as for more biologically inspired systems. Thus, in contrast to previous studies that assume only one boundary between stable and unstable, we find two boundaries: one marking the first extinction event, σ f , and a second, σ c , beyond which there are no stable fixed points for a system without a substantial loss of species. Our operational collapse definition in present simulations is the absence of any stable solution unless there is an average loss of ten species or more, with the averaging and rounding being across multiple initial conditions (σ c in the figures). A rigorous explanation of the splitting of the σ -parameter space of ecological systems into three stability phases strict stability (SS), extinction continuum (EC), and collapse (C), as seen in Fig. 1, can be found in [25]. Our previous analysis in [25], however, lacks the exact expression for the variance of species abundances [Eq. (8)].
The extinction continuum--a phase between the firstspecies extinction and collapse that encompasses May's prediction for collapse-is a region where a system can remain stable if perturbed through single-species extinctions. Thus, the EC shows how ecosystems can respond to perturbations without collapsing. The new collapse boundary-after the succession of extinctions-cannot be obtained directly from order statistics as σ f . Although, we can use order statistics to estimate the number of viable species, n, for a certain σ , and thus in combination with σ M (n) (which is accurate when all species are assumed nonextinct [23]), the collapse of the system can be predicted, as outlined in [25]. Consequently, an estimation for actual collapse is at the σ where the predicted amount of nonextinct species n reaches σ = 1/ √ cn. Having now developed a comprehensive framework for assessing stability of these complex systems, we are ready to answer the main question of this paper: Are there biologically relevant structures that substantially shift the stability boundaries explained above, possibly in such a way that σ M < σ f , meaning single-species extinctions do not occur before collapse? As seen in Fig. 2 and noted in previous studies [15,23], sign symmetry (ρ = 2/π ) shifts σ M to smaller values of σ . Although some previous studies have found that more symmetric interaction matrices can be stabilizing, we find that symmetry acts to destabilize GLV-type models [26].
Despite the shift in σ M for sign-symmetric matrices, since σ f is almost unchanged (seen in Fig. 2) the boundaries do not cross. In the case of symmetric matrices (ρ = 1), which have the smallest σ M , symmetry will violate the assumption of independence in the outer sum in Eq. (7), shifting σ f . However, this shift of σ f is in the same direction as σ M , thus keeping separation between the boundaries. A shift of boundaries in the same direction also occurs for large values of |μ|.

IV. THE ROW-SUM CONSTRAINT
From our discussion so far we can glean two intriguing features. First, different aspects of the interaction matrix are involved in the two boundaries. The first extinction boundary is dependent on sums and multiplications of the entries of A, while the stability boundary depends on its spectrum. These aspects can interact but need not. Second, the way to significantly shift the first extinction boundary is to correlate the entries of the outer sum in Eq. (7), thereby decreasing the variance in Eq. (8) of the minimum distribution Eq. (10). Based on these insights, our strategy to answer the question of whether there are biologically relevant structures that significantly alter the EC is to search for features that correlate the entries of the sum in Eq. (7) but that leave the spectrum unaffected or with less negative real parts for smaller values of σ .
From Eq. (7) the most obvious departure from independence is to constrain the row sums of the interaction matrix to some constant. Such a constraint leaves the spectrum unaffected. A global constraint of this kind, akin to energy or momentum constraints in physics, is indeed biologically reasonable. Physical resources, energy, time, and space are all degrees of freedom that constrain biological systems and affect the interaction patterns of species. Anyone of these may result in a row-sum constraint, where species balance negative and positive interactions with their ecosystem neighbors. This connection is further discussed below.
To investigate row-sum constraints and their dynamical consequences, we use an interaction matrix constructed as A = (1 − ξ )A c + ξ A 0 , with A c a random matrix with rows shifted to make the sums constant j A c;i j = B according to Here A c;i j is an entry of the matrix A c and n i the number of nonzero entries in row i. A 0 is a random matrix A 0 ∼ N (0, 1) with entries at the same positions as A c (only disturbing the existing interactions), and we set B = 0 for simplicity. The parameter ξ regulates the row-sum variance of A. Written in this form, the row-sum constraint can be enforced to different degrees by varying ξ . The distribution of the nonzero entries of A will slightly change with this construction N (0, 1) → N (0, (1 − ξ ) 2 + ξ 2 ). This in turn shifts the stability boundary to σ M = 1/ cN (1 − ξ ) 2 + cNξ 2 .
This construction of A together with the conclusion that μ = 0 does not affect σ f , allows us to handle these types of  7). The terms with A c in the rightmost position disappear in the outer sum because the row sums are zero. This leads to X * i normally distributed (central limit theorem) with mean 1 but with variance according to This reduces to the solution without constraint when ξ = 1. Averages of the first extinction event from dynamics with row-sum constrained interaction matrices are shown in Fig. 3 together with example simulations. As expected, σ f occurs at larger σ as the constraint is more strictly enforced. For sufficiently constrained systems σ f converges to σ M . This is seen in Eq. (13). Small ξ reduces the variance of the species abundances unless the denominator compensates by approaching zero, which is exactly when σ → σ M = 1/ cn(1 − ξ ) 2 + cnξ 2 . For weakly constrained systems the variance of species abundances, σ + , is large enough for extinctions to occur before it diverges at σ M , meaning collapse is always preceded by an extinction continuum. For significantly constrained systems, when σ f converges to σ M , the approach to collapse is highly variable.
We notice two types of collapses in these cases. Either the system collapses at σ M before any extinction, or some species abundances suffer a sharp decrease as the variance diverges at σ M , leading to an extinction continuum and eventual collapse as for less constrained systems. These two collapse types are shown in Fig. 3. For the type-2 collapse, the drop in abundance is sharp enough (and not predicted by |Re(λ) max | as in [27]) to merit being called collapse in and of itself.
For systems both with or without the row-sum constraint (except in cases with collapse type 1) an estimation for actual collapse is at the σ where the amount of nonextinct species n reaches σ = 1/ √ cn. Both with and without the row-sum constraint, the nonrandom extinctions in the EC that precede collapse result in the same small changes in structure: a slight increase in interaction strength mean, and weak correlations in the interaction matrix, positive for interaction strengths within column (A ji , A ki ) and negative for those within row (A i j , A ik ). The only salient difference is that the rate of extinctions when increasing σ tends to be higher for highly constrained systems (see near the first-extinction boundary for collapse type 2 in Fig. 3).

V. DISCUSSION
It is natural to ask if any biological or physical constraints in real ecosystems could lead to a row-sum constraint. In ecosystem interaction webs, there are abiotic constraints in terms of energy, space, time, and other resources. Part of such constraints are translated into biotic constraints such as cost of camouflage, prey abundance, foraging time, and many others. Because of both abiotic and biotic constraints, species in general need to make tradeoffs when interacting with their fellow species and environment [28][29][30]. Our row-sum constraintwhere positive and negative interactions are balanced-can be interpreted as representing these trade-offs. A species investing energy in foraging thereby gains more resources (large positive interaction with its prey) while at the same time exposing itself to predators (larger negative interaction with its predators). Another type of tradeoff, in case of a nonzero row-sum constraint in, for example, mutualistic interaction webs, is that between specialist and generalist. Species can either spread their efforts on many weak interactions or interact strongly with a few other species. These interpretations make the row-sum constraint not only biologically plausible, but even likely.
One of the strongest and most pervasive constraints across consumer and resources is Damuth's law: the rate of energy consumption of each species is approximately constant across species [31]. This law leads to the number density of individuals within a species decreasing as body mass, M, to the −3/4 power (M −3/4 ) and to biomass density increasing as (M 1/4 ). That is, species with small body sizes (which consume less per individual) compensate with higher numbers to consume the total amount of resources per time as a species with a large body size. If a row-sum constraint is imposed on interaction strengths (i.e., consumption rates) in the interaction matrix, then this would immediately lead to a correspondence between Dmauth's law and a row-sum constraint, providing a compelling ecological reasoning as to why this constraint might exist and systems might exhibit the properties that we have explained in our current study. Indeed, previous theoretic studies have shown that stability is enhanced when species abundances and body sizes scale according to Damuth's law [12] but that this scaling alone cannot stabilize a community without also having negative intraspecific interactions [32]. If the row-sum constraint is found to be a valid representation of Damuth's law in future research, this would support and enlarge the stabilizing role of the energy scaling of ecosystems.
The likely existence of constraints and their large impact on the first extinction boundary and width of the extinction continuum can lead to radically different behaviors of ecosystems under pressure. A highly relevant example of this is loss of habitat. In case of habitat loss, organisms are forced together more tightly, often increasing interaction strengths that correspond to an increase in σ . From this perspective, the EC mimics a sequential loss of habitat and demonstrates that such changes do not inexorably lead to collapses but can be preceded and anticipated by singlespecies extinctions. However, the stricter the constraints are, the higher the frequency of extinctions is under continued habitat loss and the shorter the period of extinctions is before collapse.
If constraints are such that the EC is present, the singlespecies extinctions are not equivalent to a random draw. Consequently, interaction matrices of systems in the EC are no longer completely random. This makes an interesting connection to assembly models where immigration or speciation together with extinctions are used to form ecosystems [33][34][35]. The closest resemblance is to assembly models implementing immigration of species with random interspecific interaction strengths since systems in the EC are subsets of initial random communities. Remarkably, in such systems the increase in structure is very slight. Interaction strengths are still normally distributed with a slight shift to larger mean [25]. This finding agrees with random immigration assembly models that find such communities indistinguishable from random ones based on analysis of their eigenvalue spectrum [34]. However, the extinctions notably do introduce small correlations, positive for (A ji , A ki ) and negative for (A i j , A ik ) [36]. On the other hand when speciation is used to assemble, communities develop a richer structure with more clustering, and subsequently more cascades in extinctions [34]. This might imply that for highly niched systems the EC would have a more clustered succession of extinctions. In additio,n the fact that communities formed in assembly models saturate at a certain diversity, with turnover if continued migration/speciation occurs, makes it plausible that many systems reside in the EC, although it is hard to tell how constraints would alter this conclusion.
Notwithstanding all insights from GLV models as used in this paper and elsewhere, it is important to put the results in context. Modeling choices are always made, many of which can influence a system's properties' effects on stability. As mentioned above, symmetry and antisymmetry in species interactions are examples of such properties. In GLV type models with direct and (in terms of mechanism) unspecified interactions, symmetry is destabilizing while antisymmetry acts to stabilize [23,37]. In contrast, models with an explicit mechanism-such as consumers interacting indirectly with a resource species-have demonstrated that the influence of symmetry can be stabilizing [26]. Overall, comparisons across studies reveal that the impacts of symmetry and antisymmetry on stability are not yet resolved. Or, to phrase it differently, the impacts clearly depend on the details in the models, much like the influence of diversity [38][39][40]. It thus remains unexplored and left for future research what the effects of constraints may be in mechanistically explicit models and models with other types of functional responses.
In summary, the behavior of Lotka-Volterra dynamics including extinctions is generic but can be altered with certain types of constraints on the interaction matrix that break the assumption of independence of the random entries without affecting the spectrum. For such systems the approach to collapse is highly variable and can be without extinctions such that the collapse is located approximately at the classical stability boundary. Or the system can have a phase of single-species extinctions that locates the actual collapse at larger values of σ . It is therefore vital to know what type of constraints exist for a system because they qualitatively change the system's response to species abundance perturbations or changes in structure, represented here by the standard deviation of interaction strengths.