Spatial ecology of territorial populations

Significance All organisms live in spatial communities. In many cases, such as vegetation or bacterial biofilms, dense surface-bound populations compete for both resources and physical space. How do these territorial interactions impact ecosystem behavior and biodiversity? We study a theoretical model of territorial resource competition with trade-offs and show that many features of real ecosystems emerge naturally, including slow population dynamics that render community composition susceptible to demographic and other noise. We also observe alternate steady states, including the Allee effect in which survival requires a minimum population. Importantly, we demonstrate that biodiversity occurs robustly and can arise in territorial communities simply due to competition for resources.

In order to prevent ν from changing the overall population size or leading to negative abundances, we project ν onto the 28 surface σ νσ = 0 and rescale its magnitude such that min(nσ + νσ) ≥ 0. In Fig. 4B of the main text, we use noise strength 29 ω = 5 × 10 −3 and period ∆t = 5/δ, so ρ = 2000. In this section we show that the the well-mixed model is an explicit limit of the spatial territory model. This allows us to make 32 precise comparisons and isolate the effects of spatial structure. The spatial variation of the nutrient environment depends on 33 the nutrient diffusion time τD ≡ L 2 E/D. When nutrients diffuse across the whole system much more quickly than species 34 consume them, τD → 0 and the spatial dynamics reduce to the well-mixed dynamics (Eq. 2) with corrections of order O(τD).

35
Consider m species competing for p resources. In the region occupied by species σ, the full equation for the concentration of 36 nutrient i is given by [3] 38 The spatial coordinate x runs from 0 to L, but the nutrient equations are simpler under the change of variables x → x− σ <σ n σ , 39 so that each cσi(x) runs from 0 to nσ. The following argument applies to each nutrient independently, so we will suppress the 40 nutrient index i for brevity. We can rewrite Eq. 3 in terms of dimensionless variables u ≡ x/L,ασ ≡ ασ/E,si ≡ Si/S, to where ≡ √ τD. We will also use relative populations pσ = nσ/L, so σ pσ = 1. When is small, we can expand Eq. 4 around 44 the point = 0. We require that Cσ(u) be continuous and differentiable at the boundaries between populations, to all orders in 45 . To meet this constraint, we express the coefficients {Aσ, Bσ} as series in τD, yielding 46 Cσ(u) = S Es ασ 47 and gather terms of like order to express the nutrient concentration as 48 Cσ(u) = k k C (k) σ (u).
[6] 49 It is simpler to use the sum and difference variables F σ . The O(1) (or k = 0) term is which is spatially uniform. The only solution to the boundary conditions C (0) σ+1 is a constant value across all territories. 50 We call this C (0) , and its value will be fixed by the second-order correction.
[8] 56 This gives us an equation at each of the m boundaries. Summing over all m equations, we obtain [9]

58
Now we want to show that the O( ) correction vanishes, i.e. C (1) each territory, the only way to satisfy the boundary conditions of equal C are equal. Then we add up the m equations for C (3) (u) to be differentiable at the boundaries to obtain where C (2) σi (u) is quadratic in u (Eq. 7) and has coefficients determined by higher-order corrections. Recalling that = √ τD, 65 the population dynamics are then given by [12]

67
In units where the nutrient value v = 1, the O(1) term is equivalent to the well-mixed population dynamics (Eq. 2) with a 68 supply vector rescaled by the system size L. Therefore, for small τD, the spatial dynamics reduce to the well-mixed dynamics 69 with a correction of O(τD).

71
In Fig. 2 of the main text, we show that the diffusion time τD is a control parameter for diversity. Why does diversity decrease 72 with increasing τD? Figure S1A and B show the average fraction of surviving species and effective number of species for 500 73 random communities where s = (0.5, 0.5) (cf. Fig. 2, where s = (0.4, 0.6)). Because the nutrient supply is perfectly balanced,

74
there are no oligotrophs to cause extinctions. However, ecosystem diversity still depends on relative abundances. Figure S1C 75 shows the average steady-state populations of the three most abundant species from the communities in Fig. S1A and B.
the Rank-1 species grows to occupy the additional space. As the population of this dominant species grows, the abundances 78 become less even, reducing biodiversity as measured by both rank-abundance and by effective number of species M .

79
Why does the most abundant species reap all the benefits of increased system size? Consider the example in Fig. S1D 80 and E, which show the steady-state nutrient concentrations for the same community of species at L = 20 and L = 30, 81 respectively. When τD > 1, the species σ with the lowest Rσ = i Si/ασi will dominate because it has the lowest overall 82 nutrient requirements. (σ will drive other species to extinction if it is an oligotroph, but this is not necessary to become the 83 dominant species.) Once nσ is large enough, the cσi(x) ≈ Rσi in the "bulk" of species σ's territory. Recall that Rσi is the 84 steady-state nutrient concentration of an isolated population, so ci = Rσi ∀ i gives zero net growth locally. Mathematically, 85 this is expressed as [13] Finally, it is worth noting that at high τD, some sets of strategies lead to a computationally stiff system. In Fig. 2C and 98 D and Fig. S1A and B, we analyzed each set of strategies across the full range of τD. When Mathematica detected a stiff

108
Throughout the text, we have characterized steady-state diversity using the fraction of surviving species and the effective 109 number of species. However, ecologists often employ rank abundance curves to visualize the full distribution of abundances in 110 an ecosystem. Fig. S2 shows steady-state rank abundance curves for unbalanced and balanced nutrient supplies across a range 111 of τD. Each point is an average over many random communities. For both nutrient supplies, the well-mixed system (τD = 0) 112 has a shallow curve, indicating that all species have similar abundances. (This is for the case of equal initial abundances.

113
Because the well-mixed system has a degenerate manifold of steady states, less even initial conditions typically yield less 114 even distributions.) When the nutrient supply is unbalanced (left), even rapid nutrient mixing (τD = 0.01) allows the most 115 abundant species to occupy ≈ 70% of space. As a result, the other species' shares of abundances decrease precipitously with 116 rank. Increasing τD exacerbates this effect, with the rank-1 species controlling > 90% of space at τD = 1600. It should be 117 noted that most steady states with an unbalanced nutrient supply have fewer than ten surviving species, so many abundances 118 contributing to the high rank averages are zero. A balanced nutrient supply leads to higher diversity (right), but the population 119 of the rank-1 species still grows from 50 − 85% as the nutrient mixing time increases.  Figure S3 shows an example for (m, p) = (3, 2). Fig. S3A shows the 124 directions of the eigenvectors vi in population space. The first eigenvector v1 is the "stiff" mode of changing n1, and has a large 125 eigenvalue that does not scale with τD ( Fig. S3B). v2 is the "soft mode" of the balance between n2 and n3. This is the slow 126 manifold, with |λ2| ∝ τD (Fig. S3C). As τD → 0 in the well-mixed limit, v2 becomes the degenerate manifold of well-mixed 127 steady states. v3 corresponds to perturbations that change the total population size L. Because the dynamics do not allow n to leave the surface σ nσ = L, this is a new system so λ3 is not relevant to the dynamics.

129
Below, we show that the slow manifold is a general feature of the spatial model. To do so, we calculate the Jacobian, analyze 130 the stability of the well-mixed dynamics, and treat the spatial dynamics as a perturbation of the well-mixed eigenvalues.

131
A. Calculating the Jacobian. We found in Section 2 that we can express the population dynamics as a series in τD by expanding 132 c(x): σi (x)dx and we have truncated the series at O(τD).

135
At steady state, 137 which implies that in region σ, where Qσi satisfies Eq. 16 says that at steady state, the nutrient concentrations cσi = δ/E at O(1) for all i and σ, with spatial corrections of O(τD) which depend on both the nutrient and the region. (Note that the total population is fixed, so δ = S in units where v = 1. Thus this is equivalent to cσi = S/E.) This simplifies calculation of the matrix elements of the Jacobian J: where gσγ( n * ) is a complicated function containing all the spatial interactions. This expression holds for the diagonal elements 142 σ = γ as well. Now we can express J as the well-mixed Jacobian with an O(τD) correction: where we have used the fact that L = σ nσ = δ/E. In Eq. 21, J0 ≡ −( δ E )DM is the Jacobian of the well-mixed system,  To find the number of zero eigenvalues, it suffices to find the number of linearly independent solutions to the equation Writing this in terms of matrix elements, we have where yi ≡ (α1ix1 + α2ix2 + ... + αmixm)/Si. This is invertible for small , so the characteristic polynomial from Eq. 25 can be expressed as The small eigenvalues c can be found by dropping terms of O( ) and solving for c: where c is independent of . Thus the matrix A + B has m − p eigenvalues of the form λ = c, corresponding the m − p 174 solutions of Eq. 28.

175
Therefore the Jacobian J0 + τDG, has m − p eigenvalues |λ| ∝ τD, just as we found via numerical stability analysis. The 176 slow manifold is a general feature of our spatial model whenever there is a diverse community with more species than nutrients 177 (m > p). 178 6. Alternative steady states 179 Figure 5 of the main text shows that spatial structure can lead to multistability and the Allee effect. Figure S4  one species drives the other to extinction. In the brown regions labeled "Coexistence", the two species specialize in different 185 resources, and there is a single steady state in which they coexist. In the tan region, both species specialize in Nutrient 2 but 186 have strategies similar to the nutrient supply vector, and we observe the Allee effect: Species 2 comprises half or more of the 187 system if its initial population passes a threshold, but goes extinct otherwise. Finally, in the beige region labeled "Bistable 188 coexistence", there are two steady states where both species coexist. Typically, one of these steady states will have relatively 189 equal populations, but not the other, (see Fig. 5), so a small initial variation can change the final community composition 190 substantially.

191
The fraction of parameter space in which alternative steady states occur depends on the nutrient environment. Figure S4B 192 and D show that making the nutrient supply more symmetric shrinks the regions of multistability and the Allee effect. On the 193 other hand, Fig. S4C and D show that increasing the diffusion time τD substantially increases the size of these regions. Thus

194
τD not only controls the diversity of steady states (Fig. 2, main text and Fig. S1), but also determines whether alternative 195 outcomes are possible.
196 7. Unequal enzyme budgets Figure 6 of the main text demonstrates that in spatial ecosystems diversity beyond the competitive-exclusion limit persists for 198 unequal enzyme budgets. In fact, diversity peaks for intermediate variance in the enzyme budget E. Figure S5 shows that 199 this is due to an asymmetric effect on oligotroph versus non-oligotroph strategies. In Fig. S5A, δE * is the decrease in E at 200 which an oligotroph stops driving others extinct, or the increase in E at which a non-oligotroph strategy begins driving others 201 extinct. In other words, δE * is how much the enzyme budget must change for an oligotroph to behave as a non-oligotroph, or 202 vice versa. The further a strategy is from the oligotroph region, the more E must increase for that species to drive others 203 extinct like an oligotroph. In the oligotroph region, however, a very small decrease in E allows the oligotroph to coexist with 204 many other species. As a result, the random enzyme budgets in Fig. 6 often suppress the oligotroph effect without allowing 205 other species to become dominant.

206
Why are oligotrophs so sensitive to changes in enzyme budget? Recall that oligotrophs are defined by Rσ < p. Figure S5B 207 shows Rσ as a function of strategy. It is suggestive that in the oligotroph region, Rσ is very close to the threshold p, whereas 208 strategies away from the oligotroph region have Rσ p. For any particular species and community, the critical budget change 209 δE * is not simply the change in E necessary to cross p, but some more complicated function of τD and the other species present.

210
However, Fig. S5C shows that δE * increases monotonically with Rσ. In particular, it rapidly increases as Rσ crosses the 211 oligotroph threshold before leveling off. At lower τD there is more variation (Fig. S5D), but the separate scales for oligotrophs 212 and non-oligotrophs remain apparent.

213
8. Lattice territory model 214 We have shown that spatial territories reduce biodiversity in a one-dimensional model. Does this remains true in two spatial "territory" occupied by a single species σ. A species never leaves its territory, but can exchange nutrients with its neighbors on 226 the lattice (see Fig. S6A for an illustration). Within each territory, the nutrients are assumed to be well-mixed. The population 227 size nσ in each territory does not influence the spatial arrangement of territories on the lattice, but does influence the flux of 228 nutrients between territories. The total amount of nutrient i in territory σ, Nσi, obeys 230 where cσi = Nσi/nσ is the concentration of nutrient i in territory σ. The sum over σ includes all z neighbors of species σ on 231 the lattice. Note that the flux of each nutrient into a territory from the external supply scales with the territory size nσ, as in 232 our continuous 1D model. If nutrient levels equilibrate much more quickly than populations change, we can again assume a 233 separation of timescales (∂N/∂t = 0) and obtain equations for the steady-state nutrient concentrations: Note that if D = 0, cσi = Si/ασi = Rσi, which is the nutrient environment of an isolated species in the 1D continuous model.

236
In the case of an extinction, with nσ = 0, we note that cσi = σ c σ i /z, so nutrients still flow through the territories of extinct 237 species. Because the nutrients are well-mixed within a territory, the population dynamics are simply given by [31] 239 The total system population size σ nσ = nT remains fixed so long as S = δ. As in the continuous 1D model, we assume [32] 253 Defining x ≡ α n D + 2, we can rearrange these equations to yield [33]

255
In general, the solution to this system of m equations is quite complicated. We proceed by finding the ratios of the 256 concentrations between adjacent sites, then studying the limit m 1, which yields an exponential form for c(σ). From Eq. 33, [34] 259 1/(x − 1) < 1, so the last site has a lower concentration than the penultimate site, as expected. We can find the subsequent 260 concentration ratios iteratively: [35] 262 The proportionality constant f (k) is defined iteratively: [36]

264
The recurrence relation Eq. 35 does not hold for c(2)/c(1) because of the source at σ = 1. Using Eq. 33 and Eq. 35, we find 265 the following solution for the nutrient profile: [37] 267 We seek only the functional dependence of c(σ) on σ, so we can leave the boundary term c(1) unevaluated. Equation 37 is 268 exact, but we are interested in the decay of c(σ) near the source, in the limit of many sites. Then m − l >> 1 for every term 269 in the product, and we can focus on the behavior of f (k) as k becomes large. Note that f (k) can be written as a continued 270 fraction: [38]

272
Consider the limit F ≡ lim k→∞ f (k). An infinite continued fraction converges if the sequence of convergents approaches a limit.
273 f (k) satisfies this property for x ≥ 2 (i.e. αn/D ≥ 0), so the limit F exists. Then we can solve the following equation for F : [39] 275 We choose the smaller of the two solutions on the physical grounds that F < 1, so each site has a smaller concentration than 276 the last. Then if m σ, we can approximate Eq. 37 as We have achieved our goal of calculating the decay length λ. However, we are interested in comparing the lattice territory 281 model to the 1D continuous model, where the entire system is diffusively coupled. We expect that the lattice model will be most [44]

294
For the rest of Section 8 and in Fig. S6-S8, Eq. 43 and Eq. 44 define τD when we use it with respect to a particular lattice 295 territory model. When referring to the continuous model, we use the main text's definition with L = nT , so τD = n 2 T E/D. Figure S6B shows the outcome of 16 species competing for two resources. The left subplot uses the continuous model, and the 298 right subplot uses the lattice model. In the well-mixed case, all 16 species would coexist. We see that both spatial models 299 allow coexistence beyond the competitive-exclusion limit, but much of the biodiversity is lost relative to the well-mixed case.

300
Indeed, even the identity of the surviving species is similar between the two 1D models, although the outcome in the lattice 301 model is slightly more diverse. Figure S6C shows how the effective number of species M at steady state depends on the 302 nutrient diffusion time τD, for the same strategies and nutrient supply as A. In both models, increasing τD reduces biodiversity 303 by making the abundances of survivors less equal. Finally, Fig. S6D shows the outcome of a competition with the same 304 nutrient supply and strategies as A, except that the oligotroph strategy α4 has been moved outside the oligotroph region. In 305 both models, the number of species coexisting at steady state increases dramatically. Thus the 1D lattice territory model 306 recapitulates three essential features of the 1D continuous model:

307
• Spatial structure reduces biodiversity relative to the well-mixed case, but allows coexistence beyond competitive exclusion.

308
• τD acts as a control parameter for biodiversity by decreasing the evenness of abundances.

309
• The presence of oligotroph strategies leads to low-diversity steady states.

310
Moreover, the lattice territory model generalizes readily to 2D, allowing us to ask how dimensionality might affect these features 311 of the model.

312
D. 2D lattices. How does the extension of the lattice territory model into two spatial dimensions affect biodiversity? Figure   313 S7A shows the outcome of competition for the same strategies and supply as Fig. S6B and C, but with competitors arranged 314 on square and hexagonal lattices. These situations differ from each other and the 1D case in their connectivity. As an example, 315 consider Species 6. In 1D, Species 6 directly exchanges nutrients with two neighbors (σ ∈ {5, 7}). However, the number of 316 neighbors increases to four on the square lattice (σ ∈ {2, 5, 7, 10}) and to six on the hexagonal lattice (σ ∈ {2, 3, 5, 7, 10, 11}).

317
In spite of these differences, the impact of spatial structure on biodiversity remains consistent: diversity is lost compared to the 318 well-mixed model, but exceeds the competitive-exclusion limit. (There are intriguing differences in the number, strategies, and 319 spatial arrangements of surviving species, which will be explored in future studies.) 320 Figure S7B shows how the 2D communities in A are affected by the nutrient diffusion time τD. Again, increasing τD reduces 321 the effective number of species by rendering abundances less even. Figure S7C shows the outcomes of competition for the same 322 community as Fig. S7A and B, except with the oligotroph stategy α4 moved outside the oligotroph region. Again, diversity 323 increases dramatically in the absence of oligotrophs.

324
In Fig. S8 we show that these patterns are not particular to the strategies in Fig. S6 and S7, but are generic features of the 325 lattice model. Figure S8A shows the mean fraction of species surviving at steady state for 300 random sets of 16 strategies, 326 across all three geometries, with and without oligotrophs. For every lattice, spatial structure reduces diversity relative to the 327 well-mixed case, and removing oligotrophs allows many more species to coexist. Indeed, the oligotroph effect seems to be 328 amplified for 2D lattices. Figure S8B shows the effective number of species at steady state as a function of τD, averaged over 329 many random communities. τD controls the evenness of abundances in all cases, just as it does in the 1D continuous model (cf.

331
There are many interesting features of the lattice model which deserve further investigation. For now, however, we have 332 demonstrated that it captures the essential biodiversity patterns of the 1D continuous model, and shown that these generalize 333 to two spatial dimensions. This suggests that our conclusions from the main text -namely, that spatial structure can reduce 334 biodiversity in a resource-competition model with metabolic trade-offs -are not restricted to 1D, nor are they particularly 335 sensitive to model details.  , across a range of τ D . Each curve is averaged over 400 randomly-selected sets of ten initial strategies, and each simulation begins with equal initial populations. Right: Steady-state rank-abundance curves for the balanced nutrient supply s = (0.5, 0.5), across a range of τ D . Each curve is averaged over 500 randomly-selected sets of ten initial strategies, and each simulation begins with equal initial populations.  . In the spatial model, diversity beyond the competitive-exclusion limit persists for unequal enzyme budgets. This is because small changes in the total enzyme budget E turn oligotrophs into non-oligotrophs, but the reverse requires large changes in E. (A) The change in total enzyme budget necessary to change final diversity depends on the strategy α1. For oligotrophs, δE * is the decrease in E required for two more species to survive compared to the corresponding case with equal enzyme budgets. For non-oligotroph strategies, δE * is the increase in E required for two fewer species to survive compared to the case of equal enzyme budgets. For each α1, δE * is averaged over 100 random sets of nine other species without oligotrophs (exp log δE * ± SD , where SD is the standard deviation of log δE * ). (B) Rσ = i S i α σi as a function of strategy ασi. Oligotroph strategies are very close to Rσ = p, while typical non-oligotroph strategies are far from Rσ = p. (C) Critical budget change as a function of Rσ. δE * increases rapidly outside the oligotroph region Rσ < p before leveling off. (D) The same as C, but for smaller τ D . δE * has more variation, but the critical magnitude of δE* is still much smaller for oligotrophs than for non-oligotrophs.  state for three lattice geometries. A population is considered extinct if nσ/n T < 10 −6 (mean ± SD for 300 sets of 16 strategies with s = (0.4, 0.6)). We use n T = 1 in each geometry, yielding τ D = 16 for the 1D lattice and τ D = 1 for the square and hexagonal lattices. (B) Effective number of species at steady state for three lattice geometries, with oligotroph strategies permitted (mean ± SD for 300 sets of 16 strategies at each τ D , with s = (0.4, 0.6)).