Changing relative intrinsic growth rates of species alter the stability of species communities

11 It is perplexing when species-rich ecosystems change abruptly and, for example, dominant or 12 economically interesting species populations collapse. Although various aspects of such ecosystem 13 regime shift at tipping points have been studied, little attention has been paid to the possible 14 dependence of community stability on the intrinsic growth rates of their species. Intrinsic growth rates 15 of species can vary, e.g., due to evolution, environmental changes or fluctuations, disturbances, or 16 human influences such as exploitation of certain species. We show theoretically and computationally 17 that under certain conditions changing relative intrinsic growth rates of competing species have a strong 18 effect on community stability. In addition, we are investigating which species characteristics favour 19 such changes of communities from stable to unstable and vice versa.


Introduction
A change in stability of ecological communities can lead to regime shifts or collapses (Beisner et al., 2003;Scheffer and Carpenter, 2003;Sahade et al., 2015;Heymans and Tomczak, 2016;Kareiva and Carranza, 2018).Thus, in the current debate on the potential impacts of environmental changes on natural habitats, such as, e.g., climate change, exploitation of certain species (Bland et al., 2018;Gamelon et al., 2019) or the introduction of neobiota into a habitat (Case, 1995;Mieth and Bork, 2010), it is of vital interest to study their influences on community stability.One way how environmental changes influence communities is via intrinsic growth rates of species.These can be altered e.g., by increased temperatures (Hansen et al., 2019;Bauman et al., 2022), more frequent extreme climatic events (Lloret et al., 2012) such as droughts, increased sedimentation (Sahade et al., 2015), through phenotypic or evolutionary changes (Dakos et al., 2019), or direct support, removal or addition of species' individuals (Selaković et al., 2022) or populations (Stinson et al., 2006;Campbell et al., 2022).
To understand community change, Lotka-Volterra (LV) type models are often used, which, like theoretical ecological models in general, are based on a few features of the ecological system under study that are considered essential (Jorgensen and Fath, 2011).These models are applied for thought experiments to better understand the ecological system and evaluate consequences of their settings.A common approach is to analyze the effects of species characteristics of these LV models on the stability of equilibria that reflect the community composition in terms of species identities and abundances.These characteristics, e.g., the intrinsic growth rates, are assumed to depend on the environment.Many studies of LV models focused on the influence of interaction parameters (May, 1972;Allesina and Tang, 2012;Kessler and Shnerb, 2015;Lischke and Löffler, 2017) and carrying capacity (Liautaud et al., 2019) on community stability.However, the influence of intrinsic growth rates on stability has been investigated in only a few studies.The LV model formalisms used differ in how the intrinsic growth rates determine the dynamics.In the r-formalism (Saavedra et al., 2017;Tabi et al., 2020) equilibria depend on the intrinsic growth rates.In this formalism, it was shown that the stability of the special case of feasible and diagonally and thus globally stable equilibria depends on the intrinsic growth rates (Grilli et al., 2017;Song and Saavedra, 2018).In the Modern Competition Theory (MCT) formalism (Chesson, 2000;Song et al., 2020), it has been mathematically shown that different intrinsic growth rates can change the stability of a community (Strobeck, 1973).However, it is unclear which communities change their equilibrium stability with which changes in intrinsic growth rates.To address this question, we use theoretical and computational analyses to systematically examine how changes in positive intrinsic growth rates affect the general local stability of communities.In addition, we investigate how these effects depend on species number and interaction parameter values.We focus on competition, which is important for plant communities (Fort, 2020), microorganisms (Gause, 1932), or marine faunal communities in kelp holdfasts (Allesina and Levine, 2011), for example.

Model formulation
The Lotka-Volterra competition (LVC) model (Gause, 1932)  summarizes the limitation of species i relative to species j and the effect of species j on species i relative to the effect of species j on itself (Adler et al., 2018).

Equilibria
A solution z of the LVC model has 2 n combinations of surviving (zi ≠ 0) and extinct (zi = 0) species and therefore the model has 2 n equilibria (Goh, 1978), some of which are feasible (i.e., community dynamics with non-negative abundances) and stable (Fried et al., 2017;Lischke and Löffler, 2017).In ecological processes, only feasible equilibria   * ∈ ℝ + are of interest.Assuming an equilibrium z * has n 0 = |N 0 | extinct and n + = |N + | surviving species, where N 0 and N + contain the indices of extinct and surviving species, respectively (Lischke and Löffler, 2017;Serván et al., 2018), the equilibria of the LVC model are ,  = (  ) ,∈ + (surviving).Therefore, these equilibria and their feasibility are independent of the intrinsic growth rates ṟi.

Eigenvalues of the Jacobian
To determine the stability of equilibria, the eigenvalues of the Jacobian matrix J of the LVC model are needed.The eigenvalues are the roots of the characteristic polynomial of J and the sign of their real parts determines the stability of the LVC model.The intrinsic growth rates can change this sign for  ≥ 3 (Strobeck, 1973;Marcus, 1990).J has the elements ),  =  −      ,  ≠  .
With  * =  |= * the Jacobian at an equilibrium z * and using a similarity transformation (Lischke and Löffler, 2017),  * becomes a block triangular matrix  = (  ++  +0 0  00 ). ++ describes the influence of the n + = |N + | surviving (  * > 0, i ∈ N + ) and  00 the influence of the n 0 = |N 0 | extinct species (  * = 0, i ∈ N 0 ) on the stability, where N + and N 0 are the index sets of surviving and extinct species, respectively (Lischke and Löffler, 2017;Serván et al., 2018).The determinant of  is the product of the determinants of  ++ and  00 and consequently the eigenvalues of  and also of  * are the eigenvalues of  ++ and  00 . 00 is a n 0 x n 0 diagonal matrix with elements   (1 − ∑
For the extinct species, the intrinsic growth rates ṟi ∈ + cannot change the sign of the eigenvalues of eqn.2, i.e., they do not affect the stability of the equilibria.Thus, it is sufficient to consider the matrix  ++ of survivors.For convenience in the following, we write n and J instead of n + and  ++ and omit N + .By normalizing the intrinsic growth rates by ri = ṟi / rc, Because the constant factor rc ∈ + has no influence on the sign of the real part of the eigenvalues of J, only the relative intrinsic growth rates ri influence the stability of the equilibria.The values of the relative intrinsic growth rates are in the (n -1)-unit simplex 1).

Numerical examinations
We studied the feasibility and stability by calculations in the s-and r-space (Figure 1).We used a numerical analysis framework (Figure 1, cf. 2.2.1) to investigate the feasibility (cf.2.2.2) and the effects of changes in the intrinsic growth rate on the stability of the feasible matrices of the s-space (cf.2.2.1-2.2.5).In addition, the results were analysed using different statistical measures in different reference regions of the s-space (cf.2.2.6-2.2.11) and against the number of species (cf.2.2.9).

Randomly generated S matrices
The values smin and smax, determining the range of competition strengths in a community, were uniformly sampled (i.e., randomly sampled from a uniform distribution) in the s-space (Figure 1) within specified regions.The results of the analyses were visualized in the (smin, smax)-plane.For every point s = (smin, smax) the elements sij of the matrix S were uniformly sampled between smin and smax, with one sij = smin and one skl = smax for off-diagonal positions in the matrix.The intra-specific coefficients on the diagonal were set to sii = 1.Every point s = (smin, smax) defines infinitely many matrices S and the generation of S matrices stopped with the first feasible matrix found, i.e., one which defines a feasible solution z * of the LVC model, or when reaching a predefined maximum number  ̃m of trials to avoid endless loops (cf.2.2.2).

Probability of feasibility
Since feasibility of the S matrices is a prerequisite for the investigations of stability change, we examined how likely it is to find a feasible matrix in the s-space (cf. Figure 2 and S1).Testing the feasibility of uniformly sampled matrices S of a point s = (smin, smax) is a Bernoulli experiment, i.e., the number of trials of finding the first feasible S follows a geometric distribution.Thus, if the probability of finding a feasible matrix of a point s = (smin, smax) is pf,s, then the mean number of trials required is   = 1  , ⁄ (Ibe, 2013).Per point s, the sampling rate  ̃, =  ̃  ̃, ⁄ for finding a feasible matrix is an estimator for pf,s.Thereby,  ̃, is the number of sampled S matrices and  ̃ is the found number of feasible ones.To ensure that the estimator  ̃, satisfies a certain accuracy,  ̃, was iteratively increased until p f,s converged, i.e., the last two estimations p f,s,i and p f,s,i+1 differed by less than the absolute abs and relative rel With this probability estimator value p f,s,i+1, the minimum number of trials required to find the first feasible matrix S is  ̃min, = 1  ̃,,+1 ⁄ , i.e.,  ̃min, =  ̃,,+1  ̃,+1 ⁄ (if 0 <  ̃,+1 ).The number  ̃ = ( ̃min, ) was used as upper limit in further simulations to find a feasible S matrix (cf.2.2.1) with the predefined accuracy.This setup was limited by computational power to n = 17 for the chosen 30,000 s points in the chosen region.

Changeable stability
At each point s = (smin, smax), m = 1000 relative intrinsic growth rates vectors r ∈ n-1 were uniformly sampled from the r-space and the eigenvalues of the LVC model consisting of the matrix S and each r vector were calculated to determine (in)stability by determining numerically the eigenvalues of the Jacobian J (eqn. 3).Different stabilities for different r vectors mark then a stability change with r.Points s = (smin, smax) that change their stability with r will be referred to as changeable stability points below.

Triangle Tn of changeable stability
Changeable stability points occurred only in a triangular region Tn of the s-space, for each number n of species (cf.green lines in Figure 3, 4, and S2-S4).Outside Tn, the examined points were either always unstable or always stable.These triangles Tn are determined by the points (1, 1), (0, smax,u,n), and (0, smax,l,n), where smax,u,n (smax,l,n) is the maximum (minimum) intersection point on the smax-axis (i.e., intercept) of all lines through (1,1) and each s-point (smin,i,n, smax,i,n,) with changeable stability.Describing such a line as linear equation by smax=b smin+a and inserting the point (1,1) results a=1-b, thus smax=b smin+1-b and therefore  = . For a given line through a point (smin,i,n, smax,i,n,), it is . The intercept is at smin=0 and its value is therefore ai.

Adaptive sampling region
Since the research topic of changeble stability is limited to the regions Tn, examinations can be restriced to Tn, to minimize the computational effort.Since the triangular regions Tn (cf.green lines in Figure 3, 4, and S2-S4) change with n, we determined iteratively search regions around Tn. Starting for n = 3 with smax = [0, 40], the first region T3 was chosen arbitrarily large for the changeable stability search.To determine the search region for changeable stability for the current n, the green triangle Tn-1 of the last iteration was used.It was heuristically (based on the pre-studies) extended by multiplying the upper and lower intercepts smax,u,n-1 and smax,l,n-1 (cf.2.2.4) with 1.2 and 0.6, respectively, to ensure no changeable stability points were missed.These search regions are the yellow triangles in figs.3, 4, suppl. figs. S2-S4.This reduced the computational effort required to find enough points of changeable stability and allowed to study changeable stability for up to n = 23 species on a computer cluster.

Stability change probability
As probability with which a point s changes its stability, we calculated the mean of the stability states over all r vectors at that point.The (in)stability of a LVC model with r ∈ n-1 and matrix S was recorded by  = 1 for stability and  = 0 for instability.The mean value  , =

Mean probability of stability change
To assess the influence of the community size n on the stability behavior, we wanted to know the mean probability of a stability change in the region Tn and for any larger region containing Tn.

Mean probability of stability loss
Considering a community that is stable, but of which we know neither r nor s, we calculated the average probability that it will lose its stability if r changes.Such a stable community can be either located in the purely stable region of the s-space or in Tn, where it can be either purely or changeable stable.We looked first at a point s in Tn with changeable stability that is stable for a given r.The probability that this stable point changes to unstable is 1 -pst,s and the mean probability of such a stability loss in Tn is  , = (cf.fig.5).

Fit of axis intercepts of the lines of triangles Tn
To get a general idea of the change of the triangles' Tn size and location in the s-space with increasing community size n, we approximated the intercepts smax,l,n and smax,u,n (cf., 2.2.4) by functions of n.To this aim, we determined these intercepts of of the lower and upper lines of the triangles for each n.
Then we fitted these intercepts with different functions.Intercepts for n = 3 were excluded because the first simulation range was chosen by heuristic pre-studies.We chose the decreasing functions ) and fitted them to the smax,l,n and smax,u,n by a nonlinear least squares fit (Team, 2021), weighted with  =  3 ∑  3 =4,…,23 ⁄ to emphasize higher n.Using the lowest value of the second-order Akaike Information Criterion (AICc) to test the quality of the fit (Barton, 2020), the best fitting function was chosen (cf.suppl.tab.T1, fig. 5 (A)).

Local sensitivity of stability to changes of r
Generally, the local sensitivity quantifies how an output variable responds to changes in a given (local) value of an input variable or parameter.This measure helps to interpret the effects of the arrangement of stable and unstable r points in r-space on stability changes.Here, the equilibrium of the LVC model is either stable ( = 1) or unstable ( = 0) and hence its local sensitivity to changes in the relative intrinsic growth rate vectors r ∈ n-1 is ∫   Δ −1 Δ −1 |  (Pianosi et al., 2016).This can be aproximated by in stability is for a given r relative to all other r.The value of  depends on the location r ∈ n-1 and the distances d of the r to r with a different stability state and is able to distinguish between different arragements of stable and unstable r vectors (cf.example in suppl.fig.S6).The  were grouped in nand pst,s classes and represented as frequency distributions (fig.6, S5).

Sensitivities for random and clustered arrangements of stable points
To to get an idea about the spatial arrangements in n-1 of stable and unstable r vectors underlying the obtained frequency distributions of sensitivity , we compared these frequency distributions to those generated from given spatial arrangements of stable and unstable points.To this aim, we generated random and clustered (in the center and in the corner of the simplex) spatial patterns of 10,000 points with different probabilities p of two states of red and blue points (suppl.fig.S6).Blue and red were interpreted as stable and unstable, respectively.The sensitivity measure  was applied to these arrangements for all n and p. Comparing the resulting distributions with the distributions of the sensitivities  to a change in r (suppl.fig.S5) with the Kolmogorov Smirnov test statistic (R-function ks.test (Team, 2021)) revealed that the latter were most similar to the clustered distributions for n up to about n = 10 and to the random distributions for n > 10 (suppl.fig.S7).

Theoretical results
The theoretical examinations show that (a) the equilibria z * and their feasibility are independent of the intrinsic growth rates ṟi (eqn. 1) and (b) that the intrinsic growth rates of the extinct species cannot change the sign of the eigenvalues and therefore do not effect the stability of the equilibria (eqn 2).The Jacobian matrix J of the surviving species contains the intrinsic growth rates ṟi which (c) have therefore an influence on the eigenvalues and the stability of the equilibria (eqn.3).This influence on stability results (d) from the relative intrinsic growth rates ri, since the constant factor rc ∈ + has no effect on the sign of the real part of the eigenvalues of J (eqn. 3).

Feasibility
We studied the feasibility by calculations in the s-space.Figure 2 and suppl.fig.S1 show the probabilities pf,s to find a feasible matrix S for points s = (smin, smax) with smin ≤ sij ≤ smax (cf.2.2.2).
Overall, the probability of finding a feasible matrix decreased with increasing the number n of species, but remained high in the bottom triangle close to the diagonal dmin=max (i.e., where smin = smax).For the studied n ≤ 17, each with 30,000 s points in the selected range, always a feasible matrix S was found with probability pf,s ≥ 10 -6 .Therefore, in all subsequent simulations, a maximum of  ̃min, = 10 6 random matrices S were generated to find a feasible one for each s point (cf.2.2.1).

Stability
We used the numerical analysis framework (Figure 1, cf. 2.2) to investigate the effects of changes in the intrinsic growth rate on the stability of the feasible matrices of the s-space.In combination with heuristic pre-studies that identified interesting parameter ranges, the available computing capacity allowed a systematic numerical analysis for species numbers up to n=23.The eigenvalues of the Jacobian J, which depends on the matrix S and the relative intrinsic growth rate vectors rk = (ri)k ∈ n-1, k = 1,…,1000, showed that often the LVC model for the same S matrix was stable for some rk vectors (stable r vectors) and unstable for others (unstable r vectors).In such cases, stability changes between some of the different rk vectors (changeable stability).Points s = (smin, smax) with changeable stability are located in triangles Tn (green lines through the point s = (1, 1) in figs.3, 4, suppl.figs.S2-S4).Figure 3 shows the example of a 9 species community (n = 9) with points that were unstable (left), had changeable stability (center), and were stable (right).

Insert Figure 3.
The points of changeable stability in Tn had a probability for feasibility higher then 10  S2) increased from close to 0 (red: unstable communities for almost all rk) to 0.5 (white: stable/instable communities for half of the rk each) to close to 1 (blue: stable communities for almost all rk).Consequently, the stability change probability psc,s (cf.2.2.6) is 0 for nearly stable and instable s points, increasing to 1 at pst,s = 0.5.The triangles Tn with the changeable stability overlap with the regions of purely unstable points above and purely stable points below for all examined n (figs.3, suppl.figs.S2-S4).Thus, the points in the s-space represent communities that transition from unstable to changeable stability to stable as smax decreases.Figure 4 and suppl.fig.S2 show that with increasing n, Tn rotates counterclockwise and its area decreases, i.e., its upper and lower bounds smax,u,n and smax,l,n both decrease and approach each other.

Changeable stability depending on species number
Figure 5(A), shows the functions of the number of species (4 ≤ n ≤ 23) that best fit the smax,u,n and smax,l,n intercepts of Tn (suppl.tab.T1, 2.2.9).As the number n of species increases, smax,u,n and smax,l,n decrease and approach each other.This corroborates that the region of Tn shrinks and rotates in the direction of the diagonal dmin=max (i.e., where smin = smax).Therefore, as the number of species increases, feasible and stable or stability changing competitive communities have more and more weaker and similar interactions smin ≤ sij ≤ smax.However, the coefficients a0 of the fits of different functions (suppl.tab.T1) indicate that smax,u,n and smax,l,n possibly do not converge (a0 of smax,u,n > a0 of smax,l,n).Thus, for every n there is possibly a non-empty Tn, which probably also contains purely stable communities.
Since the fit qualities (AIC in suppl.tab.T1) of the functions for smax,l,n were very similar, it is unclear, whether a0 of smax,l,n is zero or positive, i.e., whether the region containing only stable points completely coincides to dmin=max for large n. Figure 5 (B) shows that the mean probability  ,  (cf.2.2.7) of a r induced stability change in Tn ranges between close to 0 and 4%.n,sc, the n dependent part of the mean probability  ,   of stability change in a reference region of area A ref which includes Tn, is of the same order of magnitude and decreases linearily with n.  ,   can be calculated for any reference region by dividing n,sc through A ref (cf.2.2.7).In the region Dn between Tn and the dmin=max diagonal, all s points are purely stable.With the lower boundary of Tn turning in the direction of the diagonal, Dn also becomes smaller.In the shrinking combined region of pure and changeable stability (Dn and Tn), the mean probability  , + of a (changeable or purely) stable point to become unstable starts at 0.25 and decreases to about 0.075 with increasing n (cf.2.2.8).

Sensitivity of changeable stability to intrinsic growth rate
We determined how sensitive the stability of each changeable stability point s = (smin, smax) is to distances from each vector r to all other vectors rk ∈ n-1 (cf.2.2.10).These distances indicate how much r must be changed to move the LVC model to another stability state.We calculated the relative frequency distribution of the local sensitivity measure  (cf..2.2.10) over all s points in combined classes of stability probability pst,s and n (panels in fig.6, suppl.fig.S5).The mean sensitivity (white numbers, fig.6, suppl.fig.S5) increases with n up to 0.49 and is highest for pst,s = 0.5.For pst,s close to 0 and 1,  is very high for a few r and very low for most r (left-skewed).With pst,s at 0.3 and 0.7,  is either low or high (bimodal) and at pst,s at 0.5, the distribution of  is unimodal.For large n, the peaks of the distribution are narrow, while for small n they are broader and more left-skewed.Comparison with the pattern resulting from random and clustered spatial arrangements (cf.example in suppl.fig.S6) indicates that the arrangement of the stable r vectors in the r-space is close to random for large and close to clustered for small n (suppl.fig.S7).The s points with a stability probability pst,s about 0.5 in the center of the triangle Tn (white points in figs.3, 4, suppl.fig.S2) are in average most sensitive to changes in r (white numbers in the panels, fig.6, suppl.fig.S5), and the single relatively narrow peaks imply that the sensitivity  is similar for all r.However, for s points with pst,s at 0.3 (0.7) bordering the center of Tn (light blue and light red points in figs.3, 4, suppl.fig.S2) the average sensitivity is smaller and the bimodality implies that the sensitivity to changes in r is different over the r-space.

Discussion
In our MCT formulation of the LVC model, a community is represented by the combination of the competition matrix S and the intrinsic growth rate vector r.The matrix S determines the composition of species at equilibrium and, together with r, the trajectory leading to that equilibrium.A change in the stability of the equilibrium implies a change in the composition of the community and of its dynamics.Gain of stability sets in motion the assembly of a new community, with species provided, e.g., by immigration.Loss of stability means that a community that is on its way to or already in a stable equilibrium moves away from it and toward another stable equilibrium with different species (Chesson, 2018).Such a loss has different consequences in real systems ranging from subtle changes to regime shifts (Scheffer and Carpenter, 2003).It can mean a replacement of few by other species up to the complete exchange of a community, the loss of few species up to mass extinction (Reyer et al., 2015).
Our results show that community stability can be altered by changes in the relative intrinsic growth rates of species.However, this only happens within a certain range of competitive strengths.
Changes in relative intrinsic growth rates, like other aspects of community dynamics, can be caused by many processes or drivers.Such drivers affects species differently, can trigger a change in a community from stable to unstable or vice versa, i.e., Consequently, even a change in the r of a single species can change the stability of an entire community.includetemporal or spatial environmental changes (e.g., of climate (Louthan and Morris, 2021;Usinowicz and Levine, 2021), nutrients (Griffiths et al., 2015), sediments (Sahade et al., 2015), pathogens (Mordecai, 2011;Metz et al., 2012)), human activities (e.g., climate change, selective harvesting (Turkalo et al., 2017)).Even if the environment remains the same, the relative intrinsic growth rates r change when in a community species go extinct or new species are added by introduction, invasion, or evolution (Munson and Lauenroth, 2009;Gioria et al., 2018;Linders et al., 2019).
As long as the equilibrium of the community remains stable, i.e., a tipping point has not yet been reached, a changing r does not affect the composition of the community at equilibrium which therefore does not provide an early warning sign of a change in stability (Boerlijst et al., 2013).On the other hand, if the community is not in equilibrium, a change in r will affect species abundances and dynamics, which could be used as an early warning sign.For example, different recovery rates under different r after perturbations can be such an early warning sign (Munson and Lauenroth, 2009;Gioria et al., 2018;Brondizio et al., 2019;Linders et al., 2019).At the tipping point the stability changes abruptly, but the inertia of a real system can cause it to take some time to reach a new stable equilibrium.(Lloret et al., 2012).
In detail, our results show that positive intrinsic growth rates of species do not affect species composition, abundances and community feasibility (Chesson, 2018), but that they can affect community stability when there are at least three species.Thereby, the intrinsic growth rates of extinct species have no influence on the stability of equilibria (Chesson, 2018).However, the relative intrinsic growth rates of the surviving species influence the stability and thus their coexistence.Interestingly, in the permanence theory, which takes a different view on stable coexistence, the intrinsic growth rates of species in an equilibrium do not matter (Chesson, 2018).Studies using the LVC model in the rformalism have also shown that intrinsic growth rates influence community stability, but also the feasibility (Saavedra et al., 2017;Saavedra et al., 2020;Song et al., 2021;Selaković et al., 2022).Our results show that it is the relative values of intrinsic growth rates that affect the stability of a community.This is consistent with results of the r-formalism LVC model that it is sufficient to use normalized r vectors as a measure of structural stability (Grilli et al., 2017).Any influence that alters relative intrinsic growth rates and affects species differently, can trigger a change in a community from stable to unstable or vice versa, i.e., a tipping point is crossed.Consequently, even a change in the r of a single species can change the stability of an entire community.
Our approach to determine the sensitivity has a similarity to the structural stability approach in that it combines feasibility and stability depending on parameter ranges (Rohr et al., 2014;Grilli et al., 2017).However, it is less restrictive as it considers arbitrary competition matrices and thus community compositions and applies the general Lyapunov stability criterion (Logofet, 2005).The local sensitivity used represents the risk that a change in r will move a community to a different stability state, emphasising small changes in r that are assumed to be more likely.We found that the mean local sensitivity increases with community size and is highest in the center part of Tn, where also stability change probability psc,s is highest.The sensitivity of a community is not homogenously distributed in r-space (Grilli et al., 2017) and therefore depends on the position of the community in r-space.That means, local sensitivity can be similar but also highly different for communities with the same interaction matrices but different r vectors.Overall, there is a non-negligible risk that already a small environmental change affecting r can trigger a stability change.
Each point s = (smin, smax) represents an infinite number of competition matrices S, one of which was selected that allowed a feasible community.For all species numbers studied, the points s with changeable stability are in a region Tn between and overlapping with the regions of completely unstable and completely stable points.As smax in Tn decreases, there are more and more relative intrinsic growth rates for which the communities are stable and at the same time the probability of a stability change first increases and then decreases again.The region Tn shrinks as species numbers n increase and approaches the diagonal dmin=max, implying that the competition strengths sij of the communities become smaller and more similar.The mean probability  ,  of a community uniformly sampled in Tn to change its (in)stability is less than 4% independent of n.The mean probability  ,   of a stability change in a reference region (which includes Tn,), decreases with n and also with increasing the reference region.Thus, the larger the reference region, the less frequently a changeable stability community is found by random sampling.However, the mean probability  , + of a given stable community -that is necessarily either purely or changeable stable with a stable r vector -to become unstable is between 25% and 8% and decreases with n (cf.fig. 5 (B)).This risk of stability loss is not negligible, but species rich communitities are less susceptible to it than species poor ones (De Boeck et al., 2018).Additionally, a change in intrinsic growth rates is much more likely to cause a stable community to lose its stability than a community randomly selected in an arbitrary region to change its stability.Because stable communities are in restricted regions (Tn and Dn) whose sizes decrease with n, they are difficult to hit by random sampling in any larger region, in particular for large communities, which contributes to the stability-diversity debate (Gardner and Ashby, 1970;May, 1972;Goh and Jennings, 1977;McCann, 2000).Additionally, the regions with (changeable) stability are located where the probability of feasibility is high.This is in accordance to and results of the structural stability approach (May, 1972), where feasibility implies stability, that "the larger the system is, the smaller is the set of conditions leading to coexistence" (Grilli et al., 2017).But the fits of the intercepts of the changeable stability regions suggest that feasible and stable randomly assembled matrices exist even for high n, similarly to previous findings (Serván et al., 2018).This provides a new aspect to the theoretical explanation for the existence of species-rich competing communities (Hutchinson, 1961;Chesson, 2000;Li and Chesson, 2016;Saavedra et al., 2017;Schreiber et al., 2023).Stone, L. (2018).The feasibility and stability of large complex biological networks: a random matrix approach.Scientific Reports 8(1), 8246.doi: 10.1038/s41598-018-26486-2.Strobeck, C. (1973).N Species Competition.Ecology 54(3), 650-654. doi: 10.2307/1935355. Tabi, A., Pennekamp, F., Altermatt, F., Alther, R., Fronhofer, E.A., Horgan, K., et al. (2020).Species multidimensional effects explain idiosyncratic responses of communities to environmental change.Nature Ecology & Evolution 4(8), 1036-1043.doi: 10.1038/s41559-020-1206-6. Team, R.C. (2021)."R: A language and environment for statistical computing.R Foundation for Statistical Computing".).Turkalo, A.K., Wrege, P.H., and Wittemyer, G. (2017)

Figure 1:
Overview over computational analysis of stability changes by the relative intrinsic growth rates r, demonstrated for the example n = 3.The analysis was repeated for different numbers n of species.For uniformly sampled points s in the (smin, smax)-space (short s-space), matrices S with smin ≤ sij ≤ smax were uniformly sampled.A point s stands for an infinite number of communities, each of which is characterised by its competition matrix S. The grey tones of the matrices and the tree sizes exemplify competition strengths.The sampling was repeated until a matrix S was feasible, i.e., all species had positive abundances at the equilibrium, or the maximal number of trials was reached.This matrix S was combined with m different relative intrinsic growth rate vectors r uniformly sampled from the unit simplex = 1} (r-space) (red and blue points on grey triangle in the unit cube).The stability of the Lotka-Volterra competition model with matrix S and each sampled relative intrinsic growth rate vector r was determined by calculating the eigenvalues of its Jacobian matrix (blue points in the unit simplex are stable and red ones are unstable).Finally, three stability measures were calculated from the stability behaviour in the unit simplex for the matrices S of points s = (smin, smax).

Figure 2:
Examples of probability pf,s of finding a feasible matrix S of s = (smin, smax) points with smin ≤ sij ≤ smax for a selected number of species n by a Bernoulli experiment (cf.2.2.2).In the chosen range, for all simulated n, 30,000 s points were uniformly sampled.The probability to find a feasible matrix decreased with increasing the number n of species, but remained high in the bottom triangle Dn (cf.2.2.8).

Figure 3:
Example for n = 9 of stability probability pst,s over 1000 different relative intrinsic growth rate vectors rk (0 means always unstable and 1 always stable) of feasible LVC model at points s = (smin, smax).Center: The found points with changeable stability are in the triangle T9 which is delimited by the green lines through the point s = (1, 1) and through the points with the maximum and minimum slope (larger red and blue dot).The values smax,u,n and smax,l,n are the intercepts of the green lines at the smin = 0 axis.The search range (yellow delimited triangle) was chosen heuristically slightly larger than the triangle T8 to ensure that points with changeable stability were not overlooked.Left: purely unstable and right: purely stable s points.In T9 the region with the changeable stability overlaps with the regions of purely unstable points above and purely stable points below.For all simulated n, cf.suppl.figs.S2-S4.
Figure 4: Points s = (smin, smax) with changeable stability for selected n, with the stability probability pst,s over 1000 different relative intrinsic growth rate vectors rk.All points with changeable stability were in the triangles Tn which are delimited by the green lines through the point s = (1, 1) and through the points with the maximum and minimum slope (two larger dots per panel).The values smax,u,n and smax,l,n are the intercepts of the green lines at the smin = 0 axis.The search range (yellow delimited triangle) was heuristically chosen slightly larger than the triangle Tn-1 to ensure that points with changeable stability were not overlooked.For n = 3, the simulation range was heuristically chosen large.Note, the scale for smax differs in the panels.For all simulated n, cf.suppl.fig.S2. + is the mean probability of stability loss, i.e., that a point s in a stable state becomes unstable.
1, … , ,  ≠  with the squared Euclidian distance   2014), which emphasizes the effect of short distance changes.Normalizing with the maximum possible local sensitivity 1, … , ,  ≠ , where all r have a stability state other than r, the local sensitivity measure is   =  ≤ 1.The measure   indicates how likely a change

Figure 6 :
Figure 6: Examples of relative frequency distributions of local sensitivities  for all intrinsic growth rate r vectors of all changeable stability s = (smin, smax) points, for selected species numbers n = 3, 5, 10, 20 (rows) and stability probabilities pst,s (columns are the means of classes of width 0.1).The white numbers are the means of .For all n, cf.suppl.fig.S5.
) ,  = , … , ,   ,   ,   ∈ ℝ + ,   ∈ ℝ is a classical nonlinear theoretical model that is often used to study community ecological research questions.It describes the changing state xi of a population characteristic, such as abundance.Ecological processes are abstracted by three parameters, the intrinsic growth rates ṟi, competition strengths ḏij, and carrying capacities Ḵi.For our study, it is important to separate ṟi from the other parameters.We transformed the LVC model so that besides ṟi only one parameter remains.Normalizing the competition strengths and carrying capacities with the intra-specific competition ḏii by dij = ḏij / ḏii and Ki = Ḵi / ḏii (which are the equilibria of the species if they do not interact), then substituting zi = xi / Ki, (the proportion of abundance to equilibrium abundance of species if they do not interact) and eij = Kj / Ki = Ḵj ḏii / Within the triangle Tn we calculated the mean probability of a stability change over all kn examined s points as In an arbitrary reference region with area A ref which contains Tn, the sample size is A ref  and the expected value of the number of changeable stability points is  ,     , which is the same as  , . Then the expected value of the number of changeable stability points in Tn is  ,      with area    = 1 2 ( max,, −  max,, ) of Tn, sample density , and consequently sample size    .  because changeable stability only happens in Tn.Therefore, the mean probability of a stability change over all kn examined s points in A ref is  ,   = ,  and inverse linearly on A ref .In the evaluation (cf.fig.5), we used  , =  ,     , i.e., the n dependent part, from which  ,   can be calculated for any reference region by dividing through A ref .
With kn (kn,st)the number of all (always stable) s points in Tn, the fraction of the changeable stability points s in Tn is kn,sc / kn and that of changeable stability and always stable points s is (kn,sc + kn,st) / kn.Additionally, the triangle Dn below Tn down to the diagonal dmin=max = (smin = smax) (i.e., where smin = smax) which according to the simulations contains only purely stable s points has the )is the fraction of changeable stability points of all stable points.The mean probability that a stable point s loses its stability is then -4 (compare fig. 2 with 4, and suppl.fig.S1 with S2).The stability probability pst,s for points with changeable stability (figs.3, 4, suppl.fig.