The geometry of evolved community matrix spectra

Random matrix theory has been applied to food web stability for decades, implying elliptical eigenvalue spectra and that large food webs should be unstable. Here we allow feasible food webs to self-assemble within an evolutionary process, using simple Lotka–Volterra equations and several elementary interaction types. We show that, as complex food webs evolve under \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${10^5}$$\end{document}105 invasion attempts, the community matrix spectra become bi-modal, rather than falling onto elliptical geometries. Our results raise questions as to the applicability of random matrix theory to the analysis of food web steady states.

: Vulnerability of resident species in three types of food webs. a, Cumulative probability distribution of residence times. b, Cumulative probability distribution of extinction event size.
Resident times (Fig. S1a) and fractional extinction event sizes (Fig. S1b) of three different food web types. Resident times in all three food webs fall off approximately like ∝ exp −bt 1/c , though less accurately for large t. For the treelike food web we observe b ≈ 11 and c ≈ 9, whereas for the two food webs with network loops we observe b ≈ 12.5 and c ≈ 8. S3 Analytical spectrum of food webs with two species The only feasible food web with two species is that of one producer and one species consuming the producer. The steady states of this food web are for the producer and consumer, respectively. Here β = β 21 and η = η 21 . Inserting this in Eq. (6) and diagonalising yields the eigenvalues The food web is stable if the real parts of all eigenvalues are negative. The real part of Re(λ − ) is always negative. Re(λ + ) is only non-negative if the square root is real and equal to or greater than kα2 2βη . After some basic algebraic operations this criterion reduces to the following restriction on the species parameters Eq. (S3) is further restricted by feasibility of the food web, which requires the steady states of Eq. (S1) to be positive. S * 1 is always positive, but S * 2 is only positive if βη > α2 1−α1/k . Inserting this in Eq. (S3) now yields This condition can never be satisfied and consequently a feasible food web with two species will always be stable.
The eigenvalues are purely real if the square root of Eq. (S3) is bigger than or equal to zero. Again using some basic algebra this criterion transforms to This second order polynomial is zero when and negative in the interval between the two roots. (βη) − is always negative, since kγ > 0 for all α 1 < k. This root is not physically meaningful. Accordingly, the eigenvalues are purely real if and only if βη ≤ (βη) + .
S4 Omnivorous eigenvalue spectra of species richness 2-10 Figure S3: Complex eigenvalue spectra of evolved omnivorous food webs with β = 0.75. Each panel represents the two-dimensional histogram in the complex plane. Species richness as labelled in panels. Note that the colour scale is logarithmic, with green marking the areas with largest likelihood of eigenvalues. Left row corresponds to the omnivorous food webs in Fig. 3. As discussed in the main text, changing the invasion mechanics does not affect the spectrum notably.
Raw moment of n i Table S2: The first three raw moments of the Binomial distribution, n i ∼ B(N, p i ). Adapted from [2] S5 Deriving the connectivity of omnivorous food webs In this section we derive Eq. (4), that is the connectivity of omnivorous food webs as a function of number of species in the food web, N , given no extinctions occur. The connectivity of a food web is the total number of interaction links in the food web, multiplied by two (since every interaction appears twice in the community matrix), and divided by the total number of off-diagonal elements, N (N − 1). Here, we disregard the diagonal of the community matrix since it represents self-regulation and is set to −d in the random matrix [1].
There are three types of species in an omnivorous food web in terms of number of interaction links: primary producers, with links to all other primary producers in the food web, consumers with one resource (1 link), and consumers with two resources (2 links). In addition, any of the three can be the resource of another species. To obtain the average number of interaction links stemming from each species type, we need to consider all possible configurations of the food web.
Starting with the primary producers, we compute the average number of producer links, L p Here, n p denotes the number of producers in the food web, and n p (n p − 1) is the number of producer links in a given food web. Since all food webs start with a producer, the position of the first producer is fixed and only n p − 1 producers can be rearranged among N − 1 species. For the same reason, the probability that an invasive species is a primary producer, 1/3, is raised only to the power of n p − 1. 1 − 1 3 = 2/3 is the probability that an invasive species is not a producer, i.e. is a consumer with one or two resources. We also need to distinguish the consumer types when we consider the different configurations of food webs with N > 2. However, this is already accounted for by the probability 2/3. The probability of a specific configuration of consumers is (1/3) N −np . The fact that any consumer can be replaced by the other consumer type without affecting the configuration of producers adds a factor 2 N −np .
We can rewrite Eq. (S7) as and use that ⟨n Then we compute the average number of links of the consumer species, following mostly the same pattern as Eq. (S7). First we consider consumers of one resource, L c1 . Because the probability that the second species in the food web is a one-resource consumer is twice the probability that any other species is a one-resource consumer, we treat the two cases separately (S10) The first sum represents the case where the second species to be added to the food web is a consumer of one resource, whereas the second sum represents the case where the second species is a producer. The factor in front of each sum is the probability of each case, respectively. In both sums we sum over 2n c1 since each one-resource consumer contributes with one interaction, and every interaction appears twice in the community matrix. In the first sum the upper limit is N − 1. Here, the position of one of the one-resource consumers is fixed, hence only n c1 − 1 one-resource consumers can be rearranged among N − 2 species. In the second sum the upper limit is N − 2 because the two first species in the food web are producers, therefore n c1 one-resource consumers can be rearranged among N − 2 species.
1/3 is the probability that an invasive species is a one-resource consumer in food webs of N ≥ 2, and is therefore raised to n c − 1 in the first sum and n c in the second. As in Eq. (S7), we do not need to explicitly consider the different configurations of producers and two-resource consumers. We rewrite Eq. (S10) as (S12) Here, we have used Tab. S2 to replace the sums. Lastly, we compute the average number of links stemming from the consumer of two resources, L c2 Where the sum only goes to N − 2 since two-resource consumers are only added to food webs of N ≥ 2. Finally, we combine Eq. (S9) and Eqs. (S12)-(S13) and divide by N (N − 1) to obtain the connectivity S6 The effect of β on eigenvalue spectra there is a peak around 0 on the imaginary axis, as expected since most spectra contain a significant fraction of purely real eigenvalues. As β increases, so does the width of the imaginary distribution, whereas the overall shape remains the same. nTotal is a global variable that keeps track of the total number of species in the food web. addAttempt is the iteration during which a given species invaded the food web. density, decay and level represent S i , α i and l i , respectively. resources and consumers contain interaction parameters w.r.t. the species' resources and consumers, respectively. That is, resources contains β i1 η i1 , β i2 η i2 , ...β in η in , and consumers contains η 1i , η 2i , ...η ni . All η ij and η ji representing interactions that are not present in the food web are set to zero. isProducer is a Boolean variable that is true if the given species is a primary producer and false otherwise. Lastly, the class Species contains the function computeDerivative() that computes the derivative of the given species according to Eq. (2). The class Producer contains the additional global variable nProducer which is the total number of primary producers in the food web. Furthermore, growth represents k i and computeDerivative() computes the derivative according to Eq. (1).
The class FoodWeb is shown in Fig. S6b which is mainly created to keep track of indices during the integration. feasible and stable are true if the food web during a given invasion is feasible or linearly stable, respectively, and false otherwise. prevIteration is an integer that describes how the food web behaved after the previous invasion, and is used to optimise detection of non-convergent food webs. If the food web converged to the steady stated after the previous invasion, prevIteration is set to one. If the food web did not reach the steady states within the time limit, prevIteration is set to two (linearly stable food webs) or nine (unstable food webs). Finally, prevExtinction represents the previous species to go extinct.

S8 Varying the sampling distributions of α and η
In the main text we draw decay and interaction rates from uniform distributions (see Tab. 1). Here, we run the simulation for omnivorous food webs, drawing α and η from Gaussian and exponential distributions, respectively. The

S9 Holling type-II response
The simulation is extended to allow consumption rates following Holling type-II functional response [3] r(η ki , where S i is the density of the resource species and η ki is the link specific interaction strength between the resource i and consumer k as defined in Results and Discussion. The parameter h controls the significance of the type-II response. When h = 0 Eq. (S14) reduces to type-I functional response as used in Eqs. (1)-(2). With consumption rates following Holling type-II functional response, the system of equations analogous to Eqs. (1)-(2) is no longer linear in S i . Eq. (5) is therefore not applicable and neither are the convergence criteria described in Materials and Methods. A simpler algorithm is therefore employed in order to compute the eigenvalue spectra here. The food web is now considered to be converged if and the densities that satisfy Eq. (S15) are plugged in as S(t) in Eq. (6). The overall procedure is shown in Algorithm 2. Since the type-II response community matrix eigenvalues are computed only from the food webs that satisfy Eq. (S15), we only obtain the stable eigenvalues. Starting from h = 0, we run the simulation for h = 0.1, 0.2 and 0.5. The effect on the consumption rate can be seen in Fig. S9 for two values of η. In interactions of low interaction strength Figure S10: Complex eigenvalues spectra of evolved omnivorous type-II response food webs with h = 0.5.. Each panel represents the two-dimensional histogram in the complex plane. Species richness as labelled in panels and β = 0.75.
the non-linear effects are negligible for most h, whereas for high interaction strengths the non-linearity is pronounced for all h. As we increase h the average species richness decreases, and for h = 0.5 there are no convergent food webs of species richness higher than 7, though there are non-convergent food webs of up to 9 species. Fig. S10 shows the stable part of eigenspectra of food webs with species richness 2-7. The type-II spectra resemble the type-I spectra from Fig. S3, but are more skewed towards positive real values. It seems intuitive that interactions following type-II response would act as a dampening, thereby stabilising the food webs and allowing for food webs of higher species richness. To get some insight in why this is seemingly not the case, we study the spectrum of a food web with two species, analogous to Sec. S3. The steady states of this food web are i.e. similar to the steady states of type-I response, with β ′ replacing β. From these we see that feasibility now requires Since β ′ ≤ β for all h and α 2 , we see that it requires a higher η to satisfy the feasibility criterion of Eq. (S18), compared that of type-I response.
Then we study the stability of this food web. First we compute the community matrix according to Eq. (6) Since C 22 = 0 and C 12 C 21 < 0, the quadratic formula for the eigenvalues reduces to and we see that both eigenvalues are always stable when C 11 < 0, because the square root is always smaller or equal to C 11 . After some algebra we find the following criterion for stability The RHS is always larger than 1 (for h ≤ 1), which is the upper limit on η. Food webs of species richness 2 are therefore also always stable when we introduce Holling type-II response, given that h ≤ 1. However, they might be unstable for h > 1.
We only use h < 1 in our simulations, and Eq. (S24) is always satisfied with our choice of parameters. Yet the upper limit on η decreases with h, meaning the system somehow becomes "less stable" with h. It therefore seems reasonable that larger food webs with type-II response can indeed be less stable than their type-I counterparts.
Figure S11: Distribution of eigenvalues along the real axis Lastly, we plot the distributions of type-II eigenvalues along the real axis in Fig. S11. Despite a low number of eigenvalues from the simulation with h = 0.5, all distributions appear to follow approximately the same bi-modal form as in the case of h = 0. An exception is the pronounced peak near the lower limit for species richness 2 and h = 0.2. The origin of this peak is not investigated further. S10 Matrix elements of the community matrix From Eq. (6) we have the following mathematical expression for the community matrix C ij ≡ ∂ ∂S jṠ i (S(t)) S=S * .
Here S m , m ∈ {1, . . . , n} represents a resource species of S k and S p , p ∈ {n 1 + 1, . . . , n} represents a consumer species of S k . Since Eq. (2) is linear in S k , we have that C kk = 0 for all k ∈ {n 1 + 1, . . . , n}. All matrix elements representing interactions that are absent in the present food web yield 0.