Stability criteria for the consumption and exchange of essential resources

Models of consumer effects on a shared resource environment have helped clarify how the interplay of consumer traits and resource supply impact stable coexistence. Recent models generalize this picture to include the exchange of resources alongside resource competition. These models exemplify the fact that although consumers shape the resource environment, the outcome of consumer interactions is context-dependent: such models can have either stable or unstable equilibria, depending on the resource supply. However, these recent models focus on a simplified version of microbial metabolism where the depletion of resources always leads to consumer growth. Here, we model an arbitrarily large system of consumers governed by Liebig’s law, where species require and deplete multiple resources, but each consumer’s growth rate is only limited by a single one of these resources. Resources that are taken up but not incorporated into new biomass are leaked back into the environment, possibly transformed by intracellular reactions, thereby tying the mismatch between depletion and growth to cross-feeding. For this set of dynamics, we show that feasible equilibria can be either stable or unstable, again depending on the resource environment. We identify special consumption and production networks which protect the community from instability when resources are scarce. Using simulations, we demonstrate that the qualitative stability patterns derived analytically apply to a broader class of network structures and resource inflow profiles, including cases where multiple species coexist on only one externally supplied resource. Our stability criteria bear some resemblance to classic stability results for pairwise interactions, but also demonstrate how environmental context can shape coexistence patterns when resource limitation and exchange are modeled directly.


Model Form and Simplifications
We consider the consumer-resource model with parameters as defined in the main text. Throughout our simulations and analytical theory, we generally set the ii values all to some small fraction and set the remaining ij = 1 for i = j, 32 meaning that the diagonal consumer-resource pairs are limiting. We focus on this case mainly for convenience, as it simplifies our numerical and mathematical analysis. 34 In addition, there can be added complexity in this model when consumers do not require and therefore deplete all of the resources in the system. In this case, the sum over theP ji parame- 36 ters will not actually be over all of the k indices, because only those resources that are actually consumed will lead to cross-feeding. Similarly, consumer growth rates are actually dependent on 38 the total number of resources that each consumer depletes. When consumers don't deplete all resources, this number can vary between consumers. Throughout our analysis, we measure con- 40 sumer abundances in units that eliminate this additional multiplicative factor in consumer growth rates, simplifying the model. 42

Stability Criteria
In this section, we derive the sufficient stability criteria described in the main text. We show that, 44 if a fixed point exists in the model dynamics and if our stability criteria are satisfied, the fixed point is guaranteed to be stable. We describe how these stability criteria apply to our example 46 consumption and production structures. Last, we discuss under what conditions our sufficient stability criteria are also necessary. system in Eq. (2), so we evaluate the Jacobian J of Eq. (2) at the fixed point R and N to 52 determine its stability. If the eigenvalues of the Jacobian all have negative real part, then the equilibrium is stable to small perturbations. Let's compute the Jacobian of this system: In matrix notation, the Jacobian is where C diag (and similarly diag ) is the diagonal matrix with entries from the diagonal of C (re-56 spectively, ). For convenience in the calculations ahead, we will define the short-hand notations A and B for the upper-left and upper-right blocks of J, respectively (so we have A = . Note that we abuse notation here by including the abundances R d in this definition of B, which is 60 slightly different from the definition in the main text. However, in the cases we focus on analytically, this new B is directly proportional to the B in the main text. In the following sections, we will derive 62 criteria for the stability of the Jacobian under some additional assumptions. 64 First, we show that 0 is not an eigenvalue of J, so that we can invert the matrix −λI (where λ is an eigenvalue of J) when applying block determinant rules. Suppose 0 is an eigenvalue of J; then J 66 is singular and there exists a non-zero vector, with first S components denoted by v, and second S components denoted by w, such that

Sufficient conditions for stability
Since (1 − θ) N d diag C diag is a diagonal matrix with positive entries by assumption, it is invertible and we know that v = 0. Therefore, B w = 0. But we notice that −B −1 is precisely the matrix that 70 defines N in Eq.
(3) of the main text. So B is also invertible by our assumption that N is welldefined, and we know that w = 0. This contradicts our initial assumption that ( v, w) is non-zero. 72 Thus, we conclude that 0 is not an eigenvalue of J.
Now, we apply block determinant rules to the matrix J − λI and compute the characteristic poly-74 nomial for J: As we have ruled out λ = 0, this implies that the eigenvalues of J are defined by At this point we concentrate on a simplified version of the model, in which R = r 1 and N = n 1, where 1 is a vector of all ones. We also assume that C ii = C d and ii = d for all i. Now the relevant 78 matrices in the characteristic polynomial simplify to 80 We aim to show that, given the two stability criteria described in the main text, the eigenvalues of J must have strictly negative real parts. We use proof by contradiction: Assume that there is an 82 eigenvalue λ with Re(λ) > 0. Then we will demonstrate that the matrix A + B has all eigenvalues with strictly negative real part. Eq. (S5) shows that λ are also the eigenvalues of A + B , so this 84 results in a contradiction. Therefore, we will conclude that Re(λ) < 0.
First, notice that the stability criteria directly imply that B, and hence B , given the assumption 86 Re(λ) > 0, is negative definite. Next, observe that A can be expressed as The positive factor n will not affect our conclusions, so we neglect it for simplicity. By the stability 88 criteria in the main text, the matrices I − P T and [C 1] d − C are individually symmetric. In fact, these matrices are (weakly) diagonally dominant, and so the Gershgorin circle theorem implies 90 that they are individually positive semi-definite [1]. In particular, each matrix has an eigenvector 1 with corresponding eigenvalue 0; then assuming P and C are invertible, the remaining S − 1 92 eigenvalues of the above matrices are strictly positive. For these conclusions, we use the fact that P is column stochastic (i.e. P T 1 = 1). However, the product (I − P T )([C 1] d − C) is not generally 94 symmetric, and so not positive semi-definite.
We will now show that despite this non-symmetry, A shares the same eigenvalues as a negative 96 definite matrix, and thus has all real negative eigenvalues. Before proceeding, it is necessary to deal with a small technical detail. We have seen that [C 1] d − C has one zero eigenvalue, which 98 renders this matrix uninvertible. This will complicate our calculations, so we are motivated to express A as where c is an arbitrary small constant. For c sufficiently small, the matrix c (I − P T ) + 1 r B = B is negative definite, because B is negative definite and the eigenvalues are a continuous function of 102 the matrix entries (this can be shown formally using Weyl's inequality). But now [C 1] d − C − cI is negative definite, and in particular invertible. Let V ΛV T be the eigendecomposition of The first term of Q, in the second line above is congruent to −(I − P T ), which is a negative semidefinite matrix. Then Sylvester's law of inertia implies that this term is also negative semi-definite 108 [1]. The remaining two terms are similar to negative definite matrices (B and B , respectively), implying that they are also negative definite. Thus, Q is a sum of negative definite and negative 110 semi-definite matrices. This implies that Q is negative definite, and therefore Q and A + B have all real negative eigenvalues. This contradicts the assumption that Re(λ) > 0, so we have proved 112 that the criteria in the main text are sufficient to ensure the stability of J. diagonal coefficients, then the system is more likely to be stable. Because the diagonal coefficients of B are proportional to the diagonal consumption coefficients, this relationship helps to explain 120 why C d is a key bifurcation parameter in the model. In Fig D and E, we plot the predictions for this sufficient Gershgorin condition in addition to the predictions from the second stability condition we 122 stated in the main text. 124 We now show that the example consumption and production structures which we described in the main text automatically generate symmetric B matrices, satisfying our first sufficient stability 126 criterion. We consider the two terms of the matrix

Stability criteria for the example consumption and production patterns
separately. 128 In the tradeoff parameterization, P ij = 1 S−1 when i = j while P ii = 0. C is a symmetric matrix with row (or column) sums Since C is symmetric, C ij = C ji and P T C is symmetric 132 as well. Therefore, B is symmetric since it is the sum of symmetric matrices, and the first stability condition is satisfied. Circulant matrices commute with one another, and P T − θP T ˜ C d =

134
C d P T − θP T ˜ is the product of two symmetric circulant matrices, and therefore is symmetric itself. Overall, B is symmetric. 136 Now, we consider the symmetric circulant matrix case. Since C and P T are both symmetric and circulant, we have that (P T C) T = C T P = CP T = P T C. In other words P T C is symmetric. 138 Similarly, the second term in the B matrix is also a product of two symmetric circulant matrices. Therefore, B is symmetric and the first condition is always satisfied. It is possible for the second 140 condition to be false, so we need to explicitly check this condition, as we do in the main text.

142
When there is no cross-feeding, P ij =P ij = 0 for all i, j. In this section, we assume that j =i C ij = (S − 1)C 0 for analytical tractability, and we also continue to assume that R = r 1, 144 N = n 1, C ii = C d and ii = for all i. The eigenvalues λ of the Jacobian J satisfy the characteristic polynomial is not an eigenvalue of J. Therefore, we can safely apply block determinant rules to find all the eigenvalues. 150 We find that (S13) so we can write down the eigenvalues of J in terms of the eigenvalues of C. If ω is an eigenvalue 152 of C, then the roots of the polynomial are eigenvalues of J. So, we want to know when the roots of this complex quadratic have negative 154 real part. We can derive constraints on ω under which λ is guaranteed to have negative real part by using the quadratic formula and the identity Re √ x + iy = 1 √ 2 x 2 + y 2 + x, or just by 156 appealing to the Routh-Hurwitz stability criteria. In either case, if hold, then Re(λ) < 0. The first constraint in (S15) is always satisfied from the definitions. The 158 second constraint is the interesting one. First of all, if Re(ω) ≤ 0 (ie. if −C is unstable), then the Jacobian J is unstable, since the right hand side of (S15) is non-negative. So, we can already 160 rule out consumption structures (C matrices) that are unstable by themselves. In addition to the stability of C, (S15) shows that the magnitude of Im(ω) must be small compared the magnitude 162 of Re(ω). If Im(ω) = 0, which is true for a generic consumption matrix, then by sending n → 0, we will eventually get instability. These results directly mirror those presented in the main text, where 164 low consumer abundances lead to instability. Moreover, the criteria we have derived in this simpler case unify our two stability criteria in the more general case. To see this, let's interpret our more general stability criteria in the non-cross-feeding parameterization. The first symmetry criterion in the general case forces the imaginary part of the eigenvalues of C to be zero (Im(ω) = 0) and 168 the second stability criterion ensures that the matrix C is itself stable. In the non-cross-feeding case, however, we see that these stability criteria are overly restrictive. It is possible to have 170 stability without necessarily being precisely symmetric (and therefore forcing Im(ω) = 0), as long as the real parts of the eigenvalues of the matrix C are large enough. These results mirror our 172 conjecture in the main text that if the matrix B has small imaginary parts, then it will be protected from stability at low consumer abundances (see also Fig I). In addition, the bound in (S15) predicts 174 that, for larger values of r, the system will become unstable at larger values of n, which we then observe for our more complex model in Fig   In the previous section, we showed how the eigenvalues of J are directly related to the eigenvalues of C when there is no cross-feeding. Now, we prove a similar result in a different context. We 180 show that the eigenvalues of J can be written in terms of the eigenvalues of the matrices A and B provided that these matrices are simultaneously diagonalizable. In particular, the two consumption 182 and production structures which we present in the main text (the tradeoff and circulant cases) give rise to A and B matrices which are simultaneously diagonalizable. 184 In the proof of our sufficient stability criteria, we showed that each non-zero eigenvalue λ of J satisfies which means that the matrix λA + n rC d B − λ 2 I has an eigenvalue of 0. Equivalently, there exists an vector v such that for each eigenvalue λ of J. When A and B are simultaneously diagonalizable, they share the same eigenspaces, so we can find a vector v that is an eigenvector of both A and B. Let λ A 190 denote the corresponding eigenvalue of A and λ B the corresponding eigenvalue of B. Then, we can re-write the above vector equation in terms of these eigenvalues and solve for λ to get that  Fig J, we plot the predictions of Eq. S18 as well as the spectrum of J and the agreement is exact. 200 We have shown in two different simplified scenarios that J is stable if and only if B is stable provided that B is symmetric, which suggests that this may be true in general. We cannot prove 202 this statement in its full generality. In this section, we first provide additional numerical evidence that it is true. We consider a more general parameterization of the Jacobian:

Is J stable if and only if B is stable when B is symmetric?
where A is an arbitrary negative definite matrix, B is an arbitrary symmetric matrix and n is some positive constant. Note that these new definitions of A, B and n slightly abuse the notation we have 206 developed because they absorb additional constants (such as and r), but this choice does not affect any of the qualitative results. This more general parameterization of J includes the Jacobian 208 we have analyzed, but it removes the relationship between A and B that exists in our model. Instead, A and B are arbitrary matrices, as long as they are negative definite and symmetric 210 respectively. In Fig K, we plot the maximum eigenvalue of B against the real part of the eigenvalue of J with largest real part for A and B matrices with randomly sampled entries and for a few 212 different system sizes S. We vary the diagonal entries of B to create stable and unstable matrices (see the R code on Github at https://github.com/theogibbs/essential-stability-criteria for details). 214 Across many replicates, the sign of the leading eigenvalue of B predicts exactly the sign of the leading eigenvalue of J even in this more general parameterization of the Jacobian. 216 We have already proven that, if B is stable, then J is stable when B is symmetric. One way to prove the reverse implication is to show that, if B is unstable, then so is J. We will actually show 218 that, if B is unstable with an odd number of positive eigenvalues and no repeated eigenvalues, then J is unstable. Since B has no repeated eigenvalues, it admits an eigendecomposition. Therefore, 220 there exists a matrix Q such that B = Q −1 ΛQ where Λ is a diagonal matrix containing the eigenvalues of B. Because B is symmetric, we also know that Q is orthogonal so that Q −1 = Q T . Now, 222 let's set so that we can define the matrix J 1 is similar to J, so they share the same eigenvalues. Let's denote the eigenvalues of B by λ 1 > λ 2 > · · · > λ k > 0 > λ k+1 > · · · > λ S for some odd integer k ∈ {1, . . . , S} and let's define 226 the matrix |Λ| as the entry-wise absolute value of the matrix Λ. Let's also define the matrix T as a 2S × 2S block matrix with S × S size blocks. The off-diagonal blocks of T are all zeroes, the upper 228 left-hand block is the identity and the lower right-hand block is a diagonal matrix with −1 in the first k columns and 1 in the remaining columns. Then, we can write the matrix J 1 as where we have defined J 2 . Last, let's define We have already shown that J 3 is stable because it satisfies our two stability criteria. In par-232 ticular, Q|Λ|Q −1 is symmetric because Q is orthogonal: It is also similar to the matrix −|Λ| so it has negative eigenvalues and is nega- is odd. From this, we can conclude that J has at least one positive real eigenvalue. Since J is a real matrix, its complex eigenvalues come in complex conjugate pairs. The determinant is the 238 product of the eigenvalues, so these conjugate pairs do not change its sign, regardless of whether or not they have positive or negative real part, since they enter the determinant only through their 240 modulus. Since J 1 and J 2 have opposite signs, we know that at least one of the real eigenvalues of J 1 must be positive. If this were not true, then J 1 and J 2 would have determinants with the same 242 sign. So, because J is similar to J 1 , J is unstable as well. has only one positive eigenvalue because it does not have any degenerate eigenvalues. From the above argument, we know that J is also unstable. We cannot analytically rule out the possibility

Intuition for why small consumer abundances cause instability
As the consumer abundances are reduced in our model, the spectrum of the Jacobian begins to 252 have many eigenvalues with non-zero imaginary parts. These complex eigenvalues are separated into two clouds above and below the real axis (see Fig M). The imaginary parts of the eigenvalues 254 of B control the width of these clouds, while the consumer abundances determine where they are centered. If B has purely real eigenvalues, then these two clouds have no width, so regardless 256 of how small the consumer abundances become, no eigenvalues will cross over the imaginary axis. If instead B has eigenvalues with non-zero imaginary part, then the clouds in the spectrum 258 of J will have non-zero width, and some of these eigenvalues will cross the imaginary axis at a small value of n, creating an unstable fixed point. This description also helps to explain why the 260 interaction network B need not be precisely symmetric to still promote, but not guarantee, stability.
As the imaginary parts of the eigenvalues of B become smaller, so too does the width of the 262 eigenvalue clouds in the spectrum of J. As a result, networks which are not precisely symmetric but do have eigenvalues with small imaginary parts are better protected from instability at low 264 consumer abundance than those with larger imaginary parts.

266
Throughout the stability analysis in the previous section, we assumed that a fixed point R and N existed. We now evaluate when such a fixed point exists, and describe how our feasibility criteria 268 relate to our stability criteria.

General feasibility criteria
are all positive for the specified resource inflow ρ. We can solve for the resource abundances im-274 mediately because the consumer dynamics involve only one resource. We get that 276 We now show that the abundances N = n 1 is a fixed point for the example consumption and production structures that we consider in the main text. We need the vector

Feasibility for the example consumption and production patterns
to have all positive components. In each of our parameterizations, the C, P andP matrices are symmetric and have constant row and column sums. Specifically,

Variable resource abundances
In Fig 5 of the main text, we showed that our qualitative stability conclusions still applied to a model 324 where resource inflows were not chosen to ensure that all consumers had the same abundance. We still set η i = η, however, so all resource abundances were equal (as in a chemostat with 326 a constant dilution rate). In Fig H, we plot the analogous probabilities as in G, except we now sample the washout rates η i from a uniform distribution. As a result, the resource abundances 328 now vary between resources, but we observe the same qualitative transition into instability at low resource supply rates.

Matrix Parameterizations
In this section, we summarize the different matrix parameterizations we used to generate Fig   332 G which test our hypothesis that the imaginary parts of the eigenvalues of B control stability. All the code used to run simulations and generate the figures is available on GitHub at 334 https://github.com/theogibbs/essential-stability-criteria. In each case, we set the diagonal values of the consumption coefficients to be C d . We also set the diagonal values of the P matrix to be 336 zero. In Fig B, we visualize all of the consumption and production matrices we did not display in the main text.

Consumption matrices
Tradeoff: C is a symmetric matrix with row (or column) sums given by j C ij = C d + (S − 1)C 0 340 for each i. To produce this parameterization numerically, we generate a random symmetric matrix from the underlying uniform distribution. Then, we find the closest matrix with constant row sums 342 using the R package Spbsampling [2].
Circulant: First, we sample a random vector of length S from the underlying uniform distribution. 344 Then, we assign the matrix elements in the permuted fashion described in the main text. Last, we make the matrix symmetric by taking its Hermitian part.

346
Unstructured: We randomly sample all the off-diagonal entries from the underlying uniform distribution.

348
Banded: We randomly sample the off-diagonal entries displaced from the diagonal by one index from the underlying uniform distribution. Then, we constrain each row to have a row sum of C 0 . 350 Therefore, the matrix B when P is constant has only real eigenvalues but is not symmetric.
Correlated: We specify a parameter p which controls the degree of correlation between each off-352 diagonal pair (C ij , C ji ). For each pair, we sample one random value from the underlying uniform distribution, which we will call M ij . Then, C ij = pM ij + (1 − p)C ij and C ji = pM ij + (1 − p)C ji 354 whereC ij andC ji are new samples from the underlying uniform distribution. When p = 1, C is the symmetric matrix given by M , while for p = 0, it is fully random. 356 Non-symmetric Tradeoff: We sample all the off-diagonal elements of C fully randomly from the underlying uniform distribution, and then we constrain all the row sums to be given by (S − 1)C 0 + C d .
Sparse: We randomly select S/3 off-diagonal entries from each row to be non-zero. Then, we 360 sample these non-zero entries from the underlying uniform distribution.       We plot the probability of finding a feasible and stable fixed point in 5 replicates across a range of C d values and resource inflows ρ. We use three different possible resource inflow profiles (the rows). The constant inflow has all resources supplied equally, the one resource inflow supplies only the first resource and the random inflow has all resources supplied at rates sampled from a uniform distribution. In all cases, we ensure that the total resource supply to be given by ρ. We enforce that the fixed point is realized under the Liebig's law dynamics, where each consumer grows on the most limiting nutrient of all the resources. The columns are the different consumption matrix parameterizations, as shown in Fig B. Fig B also shows the corresponding production matrix used to generate the heatmap. As we described in the main text, the matrix parameterizations which generate B matrices with eigenvalues that have non-zero imaginary parts show a transition to instability at low resource inflow. In contrast, parameterizations whose B matrices have purely real eigenvalues do not show this transition. In particular, systems with larger diagonal consumption C d tend to be stable at lower levels of resource inflow ρ. In each of these cases, the resource inflow profile does not significantly change the probability of feasible and stable fixed points. Parameters: S = 10, θ = 0.5, ii = 0.05, ij = 1 for i = j and consumption coefficients sampled from uniform distributions on [0.5, 1.5] before the constraints are imposed.  (A) We plot the smallest value of C d at which the system becomes unstable or unfeasible averaged over 5 replicates as a function of the parameter p. In the correlated C matrix case, p determines the correlation between the pair (C ij , C ji ) by the formula C ji = pC ij + (1 − p)C ij where C ij and C ji are uniform random variables on [0, 1]. We take P as in the constant case. The different colors are different values of the consumer abundance n. As the consumer abundances decrease, C d must be large to maintain stability for uncorrelated matrices (ie. small p). However, when p is near 1, and hence the off-diagonal pairs (C ij , C ji ) are highly correlated, n no longer affects the numerically determined value of C d . Because the transition between these two regimes is smooth, even approximate symmetry in the C matrix promotes stability. (B) For the same range of p values and choices of consumer abundance n, we plot the average of the absolute value of the imaginary part of the eigenvalues of the matrix B = −C + P T C. This quantity does not depend on n, but decreases to zero as p goes to 1 when the matrix B has purely real eigenvalues. Parameters: S = 15, r = 1, θ = 0.5, ii = 0.9 and ij = 1 for i = j.    panel (B)). We then multiply these random samples by the equilibrium abundances and add this quantity to the equilibrium abundances to get the initial conditions. Parameters: S = 15, C d = 5, θ = 0.5, ii = 0.05, ij = 1 for i = j and η = 1.

Production matrices
The C matrices are sampled according to the sparse matrix parameterization with an underlying uniform distribution on [0.5, 1.5]. The matrix P is given by the constant parameterization. The resource inflow was ρ = 10 −4 1. These parameters are the same as those used to generate the unstable dynamics in Fig 1 of   When the consumer abundances are small, these two eigenvalue bulks appear as two clouds of eigenvalues with non-zero imaginary parts for both parameterizations. For the tradeoff parameterization, the clouds have no width, because the underlying B matrix has purely real eigenvalues. For the random parameterization, the clouds have non-zero width because B has eigenvalues with non-zero imaginary part. As the consumer abundances become small, the non-zero width of these clouds means that an eigenvalue complex conjugate pair eventually crosses the imaginary axis at a value n > 0. Parameters: S = 100, C d = 15, θ = 0.5, = 0.05, ij = 1 for i = j and r = 1. In each of the parameterizations, the consumption coefficients are sampled from a uniform distribution on [0, 2] before the constraints are imposed. In the tradeoff parameterization, the production matrix is constant, while in the random parameterization, the production coefficients are sampled from a uniform distribution on [0, 1] before the constraints are imposed.