Modeling microbial diversity with metabolic trade-offs 1

11 12 Nature exhibits much higher biodiversity than predicted by theories of 13 competition. One solution for reconciling this “paradox of the plankton” is to 14 imposes metabolic trade-offs, where species need to allocate limited cellular 15 resources into multiple functions. However, two questions exist for metabolic 16 models: first, as many such models have been proposed with diverse 17 assumptions and different results, can we find a universal language to summarize 18 various models into one unified framework? Second, under the pressure of 19 evolution, will there be a single optimal metabolic strategy that finally dominates 20 over all others? In this work, we address these two questions by constructing a 21 generalizable framework to describe the species-environment feedback in 22 chemostat-type resource-competition models. Employing this framework, a 23 fitness landscape based on the strategy-growth rate relationship can be 24 constructed. Species are capable of creating their own fitness landscape by 25 shaping their nutrient environment, which allows for dynamic fitness landscapes 26 and rich ecological behaviors, and is crucial for biodiversity in all the models we 27 examined. A non-invasible strategy corresponds to a species creating a fitness 28 landscape that places itself at the top. Under certain conditions, more than one 29 species is required to complete this task, which leads to evolutionarily stable 30 coexistence. Our approach facilitates quantitative understanding of chemostat 31 experiments, and provides insight into the competitive-exclusion paradox. 32


INTRODUCTION 33
In the natural world, species are in constant competition. So why doesn't the fittest 34 species outcompete the others and become the sole survivor? This question, captured 35 by the "paradox of the plankton" (Hutchinson, 1961), has perplexed community 36 ecologists for nearly a century. On the basis of simple resource-competition models, it 37 has been argued that the number of stably coexisting species cannot exceed the 38 number of resources, leading to the so-called competitive exclusion principle 39 (Armstrong & McGehee, 1980;Hardin, 1960;Levin, 1970 conditions. In a model where species compete for essential resources, different nutrient 60 requirements can produce intrinsically oscillatory or even chaotic dynamics that allows 61 for increased diversity  Recently, a simple model with a trade-off in nutrient uptake, was shown to self-organize 66 to a state of unlimited coexistence . This large 67 variety of models and the richness of possible behaviors raises the question of 68 unification: is there a simple framework that consolidates this diverse group of models 69 into one easily understandable picture? 70 71 A key challenge to producing such a framework is that fitness landscapes are not static. 72 Not only can extrinsic environmental factors fluctuate in space and time (Mustonen & 73 Lässig, 2009), but species can also actively reshape their habitats (Laland,Matthews,& 74 Feldman, 2016; Leibold, 1995; Odling-Smee, Laland, & Feldman, 2003). The feedback 75 loop between species and their environment produces an intrinsically dynamic fitness 76 landscape in which the action of one species can influence the fitness of all species. A 77 profound example is the Great Oxygenation Event, when cyanobacteria created an 78 oxygen-rich atmosphere (Kasting & Siefert, 2002), causing a massive extinction of 79 anaerobic bacteria but also stimulating an explosion of biodiversity (Schirrmeister,de 80 Vos, Antonelli, & Bagheri, 2013). Today, species continue shaping their habitats on all 81 scales: from humans inducing the sixth mass extinction (Ceballos et al., 2015) to 82 microbes consuming nutrients, releasing wastes, and producing toxins (Callahan,83 Fukami, & Fisher, 2014). 84 85 Resource-competition models provide a simple context to explore the interaction 86 between species and their environment (Smith & Waltman, 1995). In such models, 87 species interact only indirectly, via consumption (and sometimes production) of a 88 common pool of nutrients. A steady state can be reached if the species present can 89 shape the nutrient concentration to support a growth rate equal to their dilution or death 90 rate (David Tilman, 1982). Resource-competition models underpin many ecosystem 91 theories including contemporary niche theory as pioneered by MacArthur (MacArthur,  92 1970), popularized by Tilman (David Tilman, 1980, 1982, and extended by Chase and 93 Leibold (Chase & Leibold, 2003). A central component of contemporary niche theory is 94 a graphical approach, generally consisting of three components: zero net growth 95 isoclines (ZNGIs) in nutrient space, an impact vector representing a species' nutrient 96 consumption, and a supply point to described the external resource supply (Koffel,97 Daufresne, Massol, & Klausmeier, 2016). This graphical approach is a powerful and 98 intuitive way of evaluating the outcome of competition and community assembly, but 99 has not been commonly utilized to understand models of coexistence (Letten,Ke,& 100 Fukami, 2017), especially coexistence beyond the limit of competitive exclusion 101 Posfai et al., 2017).

103
In this work, we utilize and extend the graphical tools of resource-competition theory to 104 relate and unify multiple models for microbial diversity, emphasizing the consequences 105 of species creating their own environment. The nutrient environment shaped by one 106 species through growth and consumption may be inviting or prohibiting to another 107 species. We represent this species-environment interplay via an intuitive geometric 108 visualization of a "rule of invasion". Under various model assumptions, the species-109 environment feedback allows intransitivity of fitness, in which there is no strict 110 competition hierarchy and therefore no single best species or group of species 111 (Soliveres et al., 2015). We demonstrate how such intransitivity can lead to rich 112 ecosystem dynamics, including mutual invasion, multistability, and oscillations, and how 113 all of these behaviors can be simply related via our graphical representation combined 114 with the rule of invasion. 115 116 We extend our investigation of coexistence to encompass evolution. As species evolve 117 to adapt to their environment, an ongoing threat to diversity is that mutation/selection 118 may produce a supreme winner that takes over the habitat. To quantify the impact of 119 evolution on biodiversity, we focus on models with metabolic trade-offs: with limited 120 cellular resources, the growth rate of cells cannot increase without bound. Rather 121 evolution optimizes over cells' internal resource allocation strategies (S. Goyal In an idealized model of a chemostat (Fig 1A), types of nutrients are supplied at rate 178 and concentrations $%&&'( = ( +, supply , 3, supply , … 5, supply ), meanwhile cells and medium 179 are diluted at the same rate . However, the environment that directly impacts cells is 180 the nutrient concentration inside the chemostat, = + , 3 , … 5 , which influences the 181 intracellular metabolite concentrations and the growth rate of each species. 182 Accordingly, the biomass density of each species in the chemostat obeys: 183 (1) The concentration > of the i-th nutrient is a variable, influenced by its rate of 184 consumption by cells. For a single species with import rate > per cell volume, the 185 concentration > satisfies 186 where r is a constant representing the biomass per cell volume. (If the volume of the 187 chemostat is BCDEF$GHG and the total volume of cells is BD''$ , the import flux of the i-th 188 nutrient BD''$ • > implies a rate of change of concentration inside cells of > and a 189 corresponding rate of change of the concentration in the chemostat of BD''$ / BCDEF$GHG • 190 In this manuscript, we define as the "nutrient environment", and all possible values of 194 constitute the "nutrient space". Within a cell, the concentration of metabolites is 195 influenced by intake rate, and influences the growth rate. Different metabolic models 196 assume different forms for such influences, and we use , , to represent the 197 rate of change of : 198 = , , . (  Fig 1D). At steady state, the relation = 0 (Eq. (1) ) requires 215 the growth rate to be exactly equal to the dilution rate (assuming nonzero cell density). 216 Therefore, the contour in nutrient space satisfying = indicates all possible 217 environments that could support the steady state of the species (red curve in Fig 1D). 218 219 2. The flux-balance curve and the supply line reflect how cells shape the nutrient 220 environment for a given supply condition. At steady state, nutrient influx, dilution, and 221 consumption need to be balanced such that > = 0 in Eq. (2). Flux balance can be 222 expressed in two ways, depending on whether the unknown is the nutrient environment 223 or the supply condition: First, given a supply condition, different values of cell density 224 lead to different steady-state nutrient concentrations (Eq. (S5)), constituting a one-225 dimensional "flux-balance curve" in nutrient space (purple, cyan, and blue curves in Fig  226  1D). Alternatively, given a specified steady-state nutrient concentration, different values 227 of cell density lead to a straight line in the space of supply conditions, which we call 228 the "supply line" (see Methods for details). Despite the fact that the supply space and 229 the nutrient space are distinct, they share the same units of concentration in each 230 dimension. Therefore, for ease of visualize we typically show supply lines along with 231 other features in the nutrient space (Eq. (S6), black dashed line in Fig 1D). 232 233 3. The steady-state nutrient environment ss created by the species is the intersection of 234 the growth contour and the flux-balance curve (Fig 1D, Fig 1D), as long as 247 they fall on the same supply line, can lead to identical steady-state nutrient 248 environments (red dot in Fig 1D). Therefore, these supply conditions will lead to 249 identical nutrient limitations and thus the same cellular response. This phenomenon is 250 distinct from ratio-sensing, in which cells process dissimilar environmental inputs into result that the steady-state nutrient concentration and thus the nutrient limitation 253 depends on the relative supply of different nutrients is consistent with our experimental 254 observation (Fig 1B and C) that decreasing the nitrogen or carbon supply from a P-255 limited condition induces a similar proteome allocation as that induced by an N/C-limited 256 supply with P in excess.

258
Interestingly, changing the dilution rate alone has the potential to switch the limiting 259 nutrient. As shown in Fig 1E,  to produce the blue flux-balance curve, this curve intersects with the yellow, orange, 261 and deep red growth contours on the horizontal, horizontal, and vertical edges, 262 respectively. Therefore, as the dilution rate and growth rate increase, the species will 263 transition from nutrient b-limited to nutrient a-limited growth, even though the supply 264 concentrations $%&&'( are kept unchanged.

266
Metabolic trade-off and strategies 267 As demonstrated by our R/P ratio measurements of E. coli, microorganisms allocate 268 their limited internal resources according to the nutrient environment they perceive. In 269 our models, We use M to denote the fraction of internal resources allocated to the -th 270 cellular function, with = ( + , 3 … ) representing a metabolic strategy. An exact 271 metabolic trade-off is assumed, such that M M = 1. For example, Figure 2A shows a 272 simple metabolic model with two substitutable nutrients a and b, such as glucose and 273 galactose (see Methods for details), that contribute linearly to biomass increase. Since a 274 substantial investment of protein and energy is required for nutrient intake, the model 275 assumes a trade-off between the allocation of internal resources to import either 276 nutrient. Specifically, a fraction P of resources is allocated to import a and a fraction Q 277 (= 1 − P ) to import b. All values of P from 0 to 1 define a continuous spectrum of 278 metabolic strategies. How then shall we evaluate these strategies given that a single 279 species adopting any one of these strategies will grow at exactly the same rate as 280 dilution in a steady-state chemostat? 281 282 Rule of invasion 283 We use the outcome of competition between species to evaluate metabolic strategies, 284 assuming each species adopts a given strategy. In particular, we focus on invasion: the 285 introduction of a small number of an "invader" species to a steady-state chemostat 286 already occupied by an "indigenous" species.

288
In the graphical representation of species-environment interaction, the outcome of an 289 invasion can be summed up by a simple geometric rule, as demonstrated in Fig 2B and  290 C. The growth contour of the invader (species Red) separate the nutrient space into 291 two regions: an "invasion zone" where the invader grows faster than dilution (green-292 colored region in Fig 2B and C), and "no-invasion zone" where the invader has a growth 293 rate lower than dilution. If the steady-state environment constructed by the indigenous 294 species (species Blue) is located within the invasion zone of the invader, the invader will 295 initially grow faster than dilution. Therefore, the invader will expand its population and 296 the invasion will be successful ( Fig 2B). By contrast, if the steady-state nutrient 297 environment created by the indigenous species lies outside of the invasion zone, the 298 invasion will be unsuccessful ( Fig 2C, Figure 2D shows an example of mutual invasibility. With a supply condition different 305 from those in Fig 2B and C, while the steady-state environment created by Blue is 306 located within the invasion zone of Red, the steady-state environment created by Red is 307 also located within the invasion zone of Blue. According to the rule of invasion, each 308 species can therefore invade the steady-state environment created by the other. In the 309 face of such successful invasions, the only possible stable nutrient environment for this 310 system is at the intersection of two growth contours, where the two species can coexist.

312
This mutual invasion can be readily understood within a "fitness-landscape" picture. 313 Given an environment, we define the fitness landscape as the relation between the 314 instantaneous growth rate and the metabolic strategy of any invader (Eq. (S8)-(S9), see 315 Methods for details). Different environments give rise to different fitness landscapes. In 316 the steady-state environment created by Red ( P = 0.6), strategies with smaller P have 317 higher fitness ( Fig 2D, upper inset, red curve). In the steady-state environment created 318 by Blue ( P = 0.2), strategies with larger P have higher fitness ( Fig 2D, upper inset, 319 blue curve). Therefore, each species creates an environment that is more suitable for its 320 competitor, which leads to coexistence.

322
For the environment co-created by Blue and Red ( Fig 2D, purple dot), the fitness 323 landscape becomes flat ( Fig 2D, upper inset, purple curve): species with any metabolic 324 strategy will grow at the same rate as dilution in this environment. Therefore, in this 325 system, once a pair of species with a mutual-invasion relationship construct the nutrient 326 environment together, all species become neutral and can coexist indefinitely (see 327 Methods for details). This graphical approach to mutual invasion and the flat fitness 328 landscape provide an intuitive representation of species self-organizing to a state of 329 unlimited coexistence beyond competitive exclusion, as first reported by Posfai et al. 330 . 331 332

Rock-paper-scissor invasion loop and oscillation 333
Resource-competition models focusing on various aspects of cellular metabolism vary 334 in their assumptions regarding ( , , ), ( , , ), and ( , , ), and can lead to 335 diverse results for community structure and coexistence. However, the above general 336 "rule of invasion" allows us to treat these divergent resource-competition models in a 337 unified framework. In the following example, we utilized a metabolic model slightly 338 different from that in Fig 2, to show that a dynamic fitness landscape is indispensable 339 for coexistence.

341
In the metabolic model shown in Fig 3A, three substitutable nutrients, a, b, and c, 342 contribute additively to cell growth. In this three-dimensional nutrient space, the growth 343 contour for each species is a two-dimensional surface ( Fig 3B). In addition to requiring 344 enzymes to import the raw forms of these nutrients as in the model of Fig 2A, enzymes  345 are also required to convert the imported raw materials into biomass. In this model, a 346 six-element is required to describe the metabolic strategy, and there is the possibility 347 of "mismatches" in the fraction of internal resources allocated to import and to convert 348 the same nutrient. Such mismatches can produce a "rock-paper-scissor" invasion loop 349 ( Fig S1A): In the environment created by species 1, species 2 has a higher fitness but 350 not species 3; therefore species 2 can invade species 1 and establish its own 351 environment; however, this environment lies within the invasion zone of species 3 ( Fig  352  3B) but not of species 1, therefore species 3 subsequently invades; then in turn, species 353 3 create an environment where species 1 has the highest fitness. Such a loop of 354 invasions leads to oscillatory population dynamics ( Fig 3C, upper panel), with an ever-355 changing fitness landscape ( Fig 3C, lower panel).

357
Oscillation and even chaos in resource-competition model have been demonstrated by 358 Huisman et al. , and shown to allow dynamical coexistence 359 beyond competitive exclusion. The simple model presented here illustrates how 360 oscillation can be understood as a loop of invasion creating an ever-changing fitness 361 landscape. 362 363

Multi-stability, chain of invasion, and non-invasible strategies 364
When species create environments that are more favorable for their competitors, nutrients, so that a resource allocation strategy is fully characterized by the fraction of 371 resources P allocated to import nutrient a. The growth rate is taken to be the minimum 372 of the two input rates (Odum & Barrett, 1971). As shown in Fig 4B, two species, Red 373 and Blue, each creates a nutrient environment outside of the invasion zone of each 374 other. According to the rule of invasion, neither can be invaded by the other. Therefore, 375 the steady state of the community depends on initial conditions -whichever species 376 occupies the chemostat first will dominate indefinitely.

378
From the perspective of the strategy-growth relationship ( Fig 4B, inset), species Red 379 ( P = 0.65) creates a fitness landscape where small P is disfavored. Symmetrically, 380 species Blue ( P = 0.35) creates a fitness landscape where large P is disfavored. 381 However, neither Red nor Blue sits on the top of the fitness landscape each one creates 382 ( Fig 4C). In the fitness landscape created by Blue, a slightly larger P (green diamond in 383 Fig 4C) has the highest growth rate. Consequently, species adopting the Green strategy 384 can invade Blue. Nevertheless, species Green is not on the top of its own fitness 385 landscape as an even larger P (yellow diamond in Fig 4C) maximizes the growth rate 386 in the environment created by Green. A series of replacements by the fastest-growing 387 species in the environment created by the former species creates a chain of invasion.

389
It is worth noting that in this model after four steps of replacement, bistability appears. 390 The species with P marked by Deep Purple, which is reached by the chain of invasion 391 going from Blue, to Green, to Yellow, to Deep Green, cannot invade the original species 392 Blue. A similar relationship holds between Cyan and Red. This phenomenon highlights 393 the difference between ecological stability and evolutionary stability: Ecologically, as 394 both Blue and Deep Purple create a fitness landscape where the other species grows 395 slower than dilution, they constitute a bistable system. However, evolutionarily, mutants 396 with slightly larger P can invade Blue, bringing the system towards Deep Purple until 397 bistability collapses. 398 399 In this model, with symmetric parameters, the only evolutionarily stable strategy is P = 400 0.5 (black diamond in Fig 4C). This is the only strategy that locates itself on the top of 401 the fitness landscape it creates, and therefore cannot be invaded by any other species. 402 This simple model demonstrates a general definition of evolutionarily stable (aka 403 optimal or non-invasible) strategies: those strategies that create a fitness landscape 404 which places themselves on the top (Eq. (S10)).

406
A nutrient environment defines a fitness landscape, and the steady-state nutrient 407 environment created by species is influenced by supply condition, dilution rate, and cell 408 metabolism. Therefore, different chemostat parameters and different metabolic models 409 lead to different optimal strategies. In the following, we described a generally applicable 410 protocol for obtaining the non-invasible strategies, using the metabolic model in Fig  First, under a nutrient environment , the maximal growth rate EHW (background 413 color in Fig 4D) and the corresponding resource allocation strategy EHW can be 414 obtained analytically or via numerical search through the strategy space (Eq. (S11)). 415 EHW and EHW are independent of the chemostat parameters $%&&'( and . 416 Second, the maximal growth contour for dilution rate is defined as all nutrient 417 environments that support a maximal growth rate of (Eq. (S12)). Different 418 maximizing strategies EHW exist at different points of the maximal growth contour, as 419 shown by the colors of the curve in Fig 4D. By definition, the maximal growth contour 420 envelops the growth contour of any single strategy, and nutrient environments on the 421 maximal growth contour are outside of the invasion zone of any strategy. Therefore, if a 422 species is able to create a steady-state environment on the maximal growth contour, it 423 cannot be invaded. Finally, different $%&&'( form different maximal flux-balance curves 424 (Eq. (S14)), which intersect with the maximal growth contour at one point F&G . Species 425 EHW ( F&G ) that adopt the maximizing strategy at F&G create the environment F&G , and 426 are therefore immune to invasion. Under different $%&&'( , different species become 427 non-invasible (orange, green, and blue growth contours in Fig 4D), and the supply lines 428 emanating from different points on the maximal growth contour indicate the supply 429 conditions for which the corresponding strategies are evolutionarily stable.

431
Evolutionarily stable coexistence 432 Given and $%&&'( , the maximal growth contour and the maximal flux-balance curve 433 are unique, therefore there is only one F&G . Does the uniqueness of F&G imply a single 434 evolutionarily stable species? Or is coexistence still possible even in the face of 435 evolution? In a recent work , this question was addressed by 436 modeling a population of microbes competing for steadily supplied resources. Though 437 in-silico evolution and network analysis, the authors found that multiple species with 438 distinct metabolic strategies can coexist as evolutionarily-stable co-optimal consortia, 439 which no other species can invade. 440 441 Using a simplified version of Taillefumier et al.'s model (Fig 5A), we employ the 442 graphical approach to help identify the requirements for such evolutionarily-stable 443 coexistence and the role of each species in supporting the consortium. In this model, at 444 the cost of producing the necessary enzymes, cells are not only able to import external 445 nutrients, but can also convert any one of the internal nutrients into any other. 446 Meanwhile, nutrients passively diffuse in and out of the cell. The internal concentrations 447 of nutrient a and nutrient b are both essential for cell growth (see Methods for detail). 448 Therefore, metabolic trade-offs in this system have four elements: the fraction of internal 449 resources allocated to import nutrient a ( P ) or nutrient b ( Q ) and/or convert one 450 nutrient into another ( PQ converts internal b into a, and QP converts internal a into b). 451 Each species is defined by its internal resource allocation strategy = ( P , Q , PQ , QP ). 452 453 Following the general protocol described in the previous section, we first identified the 454 maximal growth rates EHW and the corresponding strategy or strategies EHW at 455 each point in the nutrient space, and generated maximal growth contours for different 456 dilution rates (Fig 5B). The maximal growth contours are not smoothly continuous, nor 457 are the corresponding strategies. In nutrient space, three distinct sectors of maximizing 458 strategies appear (Fig 5B, Fig S3A): When nutrient a is very low compared to b, the 459 maximizing strategy is a "b-a converter" which imports b and converts it into a (blue 460 sector, only Q and PQ are non-zero). Symmetrically, when a is comparatively high, the 461 optimal strategy is a "a-b converter" (green sector, only P and QP are non-zero). 462 Otherwise, the maximizing strategy is an "importer" which imports both nutrients without 463 conversion (red sector, only P and Q are non-zero). On the border between sectors, 464 the maximal growth contour has a discontinuous slope.

466
Evolutionarily stable coexistence occurs at these discontinuous points. If an 467 environment point X is located in a continuous region of the maximal growth contour, 468 only one maximizing strategy EHW X exists for that environment (maximizing 469 strategies along the maximal growth contour are indicated by colored squares in Fig  470  5C). Supply conditions that make EHW X the optimal strategy (i.e. allow EHW X to 471 create the steady-state environment X ) constitute the supply line for X and EHW X . 472 However, at the discontinuous points of the maximal growth contour, where two classes 473 of strategies meet, two different strategies are both maximizing. For example, at the 474 purple dot in Fig 5C a  Blue and species Red both create nutrient environments that lie within the maximal 482 growth contour (blue and red dots, Fig 5C), and are thus subject to invasion by other 483 species. Nevertheless, the species-specific growth contours of Blue and of Red 484 intersect at the purple point on the maximal growth contour. Therefore, only when Blue 485 and Red coexist can they co-create an environment on the maximal growth contour, 486 and thus resistant to invasion from any other species. Indeed, when we simulate 487 multiple species with different maximizing strategies under the supply condition 488 indicated by the open black circle, species Blue and species Red are the only two that 489 survive (Fig 5C, inset).

491
The optimal coexistence of species Blue and species Red can be understood intuitively 492 from the dynamic fitness landscape. Given a nutrient environment, the relation between 493 P and growth rate of importer (red curve) or a-b converter (green curve), and that 494 between Q and growth rate of b-a converter (green curve) constitute the fitness 495 landscape of species adopting different possible maximizing strategies (Fig 5D). In the 496 environment created by species Blue (blue dot in Fig 5C), not only will some importers 497 grow faster than Blue, species Blue (strategy marked by blue diamond) is not even on 498 the fitness peak of its own class ( Fig 5D, upper box). Similarly, in the environment 499 created by species Red, the strategy of Red is not at the top of the fitness landscape 500 (Fig 5D, middle box). By contrast, in the environment co-created by species Blue and 501 Red (purple dot in Fig 5C), their strategies are at the top of the fitness landscapes of 502 their own classes and at equal height. For all supply conditions in the gray region, 503 species Blue and species Red jointly drive the nutrient concentrations to the 504 discontinuous point of the optimal growth contour, and thereby achieve evolutionarily 505 stable coexistence. 506 507 Species creating a new nutrient dimension 508 As discussed in the introduction, one possible solution to the paradox of plankton is the 509 creation of new nutrient "dimensions" by species secreting metabolites that can be To explore the possibilities of evolutionarily stable coexistence when species create 517 new nutrients, we used a simplified model to represent multi-step energy generation 518 with a dual-role intermediate metabolite (Fig 6A). A single chemical energy source S is 519 supplied into the chemostat. The pathway for processing S consists of four relevant 520 reactions driven by designated enzymes: External S can be imported and converted into 521 intermediate Methods for details.) 532 533 The metabolic strategy in this model has four components: = ( YZ[+ , YZ[3 , DW& , 534 \E& ). When we examine the maximizing strategies and maximal growth rates in the 535 nutrient space, three distinctive classes of strategy emerge (Fig 6B). When S is 536 abundant and I ext is low, the maximizing strategies have only two non-zero components, 537 YZ[+ and DW& (Fig S3B), meaning this class of species only imports S for the first 538 energy-generating reaction then exports intermediate as waste. Therefore, we call 539 strategies in this class "polluters" (blue section in Fig 6B, Fig S3C). When I ext is high 540 while S is low, the maximizing strategies have two different non-zero components, YZ[3 541 and \E& (Fig S3B), meaning this class of species solely relies on I ext as its energy 542 source. We call these strategies "cleaners" as they clean up the I ext in the environment, 543 which is detrimental to the polluters (green section in Fig 6B, Fig S3C). When there are 544 comparable amounts of S and I ext present, a third class of maximizing strategies 545 appears: these cells neither export nor import intermediates, but rather allocate all 546 enzyme budget to YZ[+ and YZ[3 to carry out both energy-producing reactions. We 547 call species in this class "generalists" (red section in Fig 6B, Fig S3C).

549
As shown in Fig 6B, on the borders between classes of strategies in nutrient space, the 550 maximal growth contours turn discontinuously. These points of discontinuity, as in the 551 previous section, are nutrient environments corresponding to evolutionarily stable 552 coexistence of species from distinct metabolic classes. The classes of optimally 553 coexisting species change with dilution rate. When the dilution rate is low ( = 0.4, Fig  554  6C), at the discontinuous point of the maximal growth contour, the corresponding two 555 maximizing strategies are one polluter (species Blue) and one cleaner (species Green). 556 Their supply lines span a gray region where both species Blue and species Green are 557 required to create a steady-state environment on the maximal growth contour. As we 558 are only supplying the system with S, the supply condition always lies on the x-axis of 559 concentration space. For the supply condition shown by the black open circle in Fig 6C,  560 polluter Blue creates a nutrient environment (blue dot) far from the maximal growth 561 contour. When the cleaner Green is added to the system, not only does the biomass of 562 Blue increase (inset), but also the steady-state nutrient environment moves to the 563 discontinuous point of the maximal growth contour (cyan dot), where both Blue and 564 Green occupies the peaks of their fitness landscapes (Fig 6D). This result is consistent 565 with the long-term evolution experiment of E. coli and also intuitive: polluter Blue and 566 cleaner Green form a mutually beneficial relationship by, respectively, providing 567 nutrients and cleaning up waste for each other, thereby reaching an optimal cooperative 568 coexistence.

570
A quite different coexistence occurs at higher dilution rate ( = 0.6, Fig 6E). Growth 571 contours at this dilution rate show two turning points, but neither are between the 572 polluter and the cleaner class. One discontinuous point is between the cleaner class 573 (green squares) and the generalist class (red squares), but the gray region spanned by 574 the corresponding supply lines does not cover the x-axis and so does not represent an 575 attainable coexistence when only S is supplied. The other discontinuous point is 576 between the generalist class and the polluter class (blue squares). The gray region 577 spanned by the supply lines of the corresponding two maximizing strategies of 578 generalist class (species Red) and polluter class (species Blue) does cover the x-axis. 579 Therefore, a supply condition with only S within the gray region (e.g., the black open 580 circle) leads to the optimal coexistence of generalist Red and polluter Blue on the 581 discontinuous point (purple dot), despite the fact that they do not directly benefit each 582 other. Indeed, when the generalist Red is added to a system with polluter Blue and a 583 cleaner Green, the cleaner Green goes extinct and the biomass of the polluter Blue 584 decreases (inset). Nevertheless, the steady-state nutrient environment is moved from a 585 cyan dot lying inside the maximal growth contour to the purple dot lying on the maximal 586 growth contour. In the environment of the cyan dot created by cleaner Green and 587 polluter Blue, Blue is not on the top of the fitness landscape of the polluter class (Fig 6F,  588 upper box). By contrast, for the fitness landscape created by polluter Blue and 589 generalist Red (Fig 6F, lower box), despite being lower in biomass, Blue occupies the 590 top of the landscape. Therefore, the optimal coexistence of this polluter and this 591 generalist does not arise from direct cooperation, but rather from collaborating to defeat 592 other competitors. 593

DISCUSSION 594
The phrase "survival of the fittest" used to describe natural selection can be applied 595 both to competition within species and to species competing in the same environment. 596 One doctrine governing resource competition among species is the competitive 597 exclusion principle: theoretically, there should be no more surviving species than the 598 number of resources. However, the enormous diversity of coexisting species in the 599 natural world seems to contradict the competitive exclusion principle. This so-called 600 "paradox of the plankton" has stimulated many theoretical models of resource 601 competition, each with its own assumptions and different conclusions. In this work, we 602 examined a range of models for metabolic competition among microbes within a unified 603 framework, using the species-environment feedback as an organizing principle and the 604 geometric "rule of invasion" and dynamic fitness landscapes as common tools. Under 605 this unified framework, it becomes apparent how metabolic tradeoffs promote diversity 606 by allowing a dynamic fitness landscape without a fittest peak. 1980, 1982), and the growth contours in our work reduce to the zero-net growth 617 isoclines (ZNGI) introduced by this school in the particular case of two growth-promoting 618 resources. Our framework, nonetheless, differs in several aspects: First, we focused on 619 metabolic models with trade-offs, for which there are not only different ZNGIs for 620 different species but a continuous family of growth contours, and the envelope of all 621 growth contours is the maximal growth contour. Given these resource allocation trade-622 offs, the growth contours of any pair of species must intersect, clearly demonstrating 623 why metabolic trade-offs prevent a single species from unconditional dominance. 624 Moreover, the definition of growth contours is not limited by the number of resources, 625 nor constrained by whether external chemical concentrations contribute positively or 626 negatively to growth, making the approach suitable to address more realistic metabolic 627 models. Second, the introduction of the flux-balance curve, in additional to the supply 628 line, makes it easier to determine the species-specific environment for a given supply 629 condition, which is particularly useful for determining the non-invasible environment as 630 the intersection between maximal growth contour and the maximal flux-balance curve. 631 In brief, our new graphical approach is well suited to our goals of understanding 632 coexistence from the perspective of species-environment feedback, demonstrating how 633 the fitness landscape is changed by the species present, and identifying evolutionarily 634 stable strategies.

636
Our work is not aimed at adding another solution to the paradox of the plankton. Rather 637 we provide a graphic tool to unify several approaches, and suggest how different 638 proposed solutions to the paradox can emerge intrinsically from competition for 639 resources. When nutrients are substitutable, resource competition among species with 640 metabolic trade-offs has been shown to lead to emergent neutrality , 641 as the fitness landscape is made flat by the competing species. Other commonly 642 invoked solutions to the paradox of the plankton are extrinsic temporal and/or spatial 643 heterogeneity. In this work, we showed that both types of heterogeneity can also 644 emerge intrinsically from species-environment feedback. When the rule of invasion 645 allows non-transitive loops, oscillations and chaos can occur, which have been shown 646 to allow coexistence beyond competitive exclusion . 647 In addition, when multiple nutrients are all essential, the ability of each species to create 648 an environment that favors itself allows for the spontaneous emergence of spatial 649 heterogeneity in an extended system (Fig S2).

651
Yet, even if fixed species can coexist in an ecosystem, will coexistence survive the 652 ceaseless process of mutation and adaptation? Our approach provides a general 653 protocol to determine non-invasible/evolutionarily stable metabolic strategies, which we Another advantage of our graphical approach, besides providing an intuitive picture of 670 species competition, is that it can help understand and control nutrient limitation in 671 chemostat experiments. The capacity of species to shape their own environment, even 672 in a system as simple as a chemostat, presents challenges to controlling which nutrient 673 or nutrients are limiting. By traditional definition, if increasing a certain nutrient leads to 674 an increase of a cell's growth rate, that nutrient is considered "limiting". However, 675 growth rate is invariant in a chemostat, being set experimentally by the dilution rate, so 676 inferring nutrient limitation requires special attention. For example, if one sees the same 677 cellular responses under different nutrient supplies, what can one conclude? Cells may 678 be creating the same nutrient environment out of different supply conditions (cf. Fig 1D) Total RNA and total protein measurements 713 The method for total RNA and protein measurements is described in (Li et al., 2018).
>,_ is the concentration of the -th nutrient within the medium of the -th chemostat. All possible constitute the "nutrient space". In modeling population dynamics in a chemostat, multiple assumptions need to be 719 made concerning how cells sense the environment, import nutrients, export metabolites, 720 utilize resources, and grow in biomass. Different assumptions result in different 721 metabolic models. Some metabolic models focus on trade-offs in resource allocation, as 722 the amount of resources "owned" by a cell, including proteins and energy, is limited. 723 Cells need to allocate these limited resources into different cellular functions, such as 724 metabolism, gene expression, reproduction, motility, maintenance, etc. We use M,a to 725 represent the fraction of resources allocated to the -th cellular function of species , 726 with a = ( +,a , 3,a … ) representing the resource allocation strategy of species . For 727 simplicity, we assume each species has a fixed resource allocation strategy.

equations for a single species in a chemostat 730
In a chemostat with nutrient supply $%&&'( , dilution rate and a single species with 731 fixed strategy a and intracellular metabolite concentration a , the cell biomass density 732 a and the chemostat nutrient concentrations are generally described by the following 733 equations: 734 a = a • , a , a − , (S1) = • $%&&'( − − a / • , a , a .

(S2)
In considering the details of cellular metabolism, one may choose to incorporate the 735 dynamics of intracellular metabolites that originate from nutrient import and influence 736 cell growth. We make the assumption that the biomass concentration , e.g. protein 737 concentration, is constant for cells under all growth conditions. Thus, an increase of 738 total cell mass induces a linear increase of total cell volume. a is the cell mass per 739 volume in the chemostat, and is the cell mass per volume within a cell. For a 740 chemostat-to-cell flux of mass J, the concentration of the metabolite in chemostat 741 decrease by / mnopqrsPs while the concentration in cell increase by / mott . As a result, 742 the metabolites imported into cells are enriched by a factor of , and metabolites 743 secreted by cells are diluted by 1/ . Also, all metabolites are diluted by cellular growth, 744 which is generally a slow process compared to metabolic reactions and can be ignored 745 in most cases. We use a function ( ( , , ), , ) to represent the rate of change of 746 intracellular metabolite a : 747 a = , a , a , a , a .

Species-specific steady state 751
In the steady state of chemostat, Eqs. (S1)-(S3) should be all equal to zero. 752 For intracellular metabolites, as , a , a , a , a = 0 as a result of Eq. (S3)=0, 753 given an environment , the steady state of a can be expressed as a function of : 754 a * = v+ ( , a ). 755 Growth contour (GC): From the perspective of the environment influencing species, at 756 each constant environment, the steady-state growth rate is fully determined by : 757 * , a = , v+ ( , a ) . If the biomass of a species is non-zero ( ≠ 0), Eq. (S1) 758 requires * = . In the -dimensional nutrient space, this requirement defines a ( − 1)-759 dimensional surface, constituted by all environments that support an equal-to-dilution 760 growth rate. This surface reduces to the zero-growth isoclines in contemporary niche 761 theory when the nutrient space is two-dimensional and the growth rate solely relies on 762 monotonically, but is not necessarily limited by the nutrient dimension or the form of 763 the growth function. For convenience, we name this surface the "growth contour" (GC) 764 for species : 765 a ≔ , v+ , a , a = }. (S4) An example of growth contours is shown in Fig 1D.

767
Flux-balance curve (FB): Eq. (S2) describes how species act on the environment. In 768 steady state, the influx, out-flux, and consumption by species should be balanced for 769 each nutrient, which enables calculation of the biomass density-to-dilution ratio for every 770 :  In nutrient space, the steady-state environment ( a,ss ) with non-zero biomass of species 776 will be located at the intersection of the growth contour and the flux-balance curve. 777 This environment is constructed by the species via its consumption of nutrients. If a,ss 778 exists, this species can survive in the chemostat. Otherwise, this species will be washed 779 out by dilution even without competition from other species. For the following 780 discussion, we only consider species that can survive when alone in a chemostat.

782
Supply line (SL): The flux-balance curve is determined by the supply condition $%&&'( . In 783 many cases, it is helpful to derive the supply conditions that enable a species to 784 construct a steady-state environment a,ss . All possible values of $%&&'( that can 785 produce a given a,ss , form a straight line in the space of supply concentrations, 786 described by: 787 • a,ss , v+ a,ss , a , a + a,ss }, with varying non-negative values of a / . An example of a supply line is shown in 788 Fig 1D.

790
Dynamic equations for multiple species in a chemostat 791 In nutrient competition models, multiple species ( = 1 … ) each with biomass density 792 a compete for resources. They have species-specific growth rates , a , a and 793 import rates , a , a , yet all experience the same nutrient environment . Therefore, 794 Eq. (S1) and Eq. (S3) remain the same for each species, while the rate of change of 795 chemostat nutrient concentrations is influenced by the summed action of all species: 796

Multiple species steady state 798
Multiple species, even if each alone can survive in chemostat, do not generally coexist 799 when competing together. For a system starting with different species, the stable 800 steady state contains * (1 ≤ * ≤ ) species with non-zero biomass. We define these 801 * surviving species as a stable species set * , and mark the steady-state 802 environment created by this set as {a * },ss . If > 1, according to Eq. (S1), {a * },ss must be 803 located at the common intersection of growth contours formed by every species in { * }.

805
Invasion 806 Invasion is defined as introducing a small number of invaders (with biomass density 807  In evaluating invasion by a species with strategy a of any environment , we make 813 two assumptions: 814 1. The biomass of the invader is so small that it does not disturb the environment at the 815 time of introduction. 816 2. There is a separation of timescales such that the concentrations of intracellular 817 metabolites reach equilibrium instantaneously at the time of introduction of the invader, 818 therefore Eq. (S3) is always equal to zero and a = v+ ( , a ) holds. 819 Therefore, we define the "invasion growth rate" of a species with strategy a 820 introduced into environment as: 821 Invasion zone: By definition, the growth contour of the invader inv divides the nutrient 822 space into two regions: an "invasion zone" that includes all environments where the 823 invader has an invasion growth rate higher than dilution, and "no-invasion zone" where 824 the invader has an invasion growth rate lower than dilution. If the steady- Two examples of this rule of invasion are presented in Fig 2A and 2B. 831 If the growth rate monotonically increases with the concentration of each nutrient, it can 832 be proven that the invasion zone is always above the growth contour of the invader (an 833 environment " "above" the growth contour \•Ž is defined as ∃ X ∈ \•Ž , s. t. >," ≥ 834 >,X ∀ ). If the growth rate is not monotonically increasing with nutrient concentrations, 835 identifying the invasion zone requires more model-specific analysis. 836 837 Fitness landscape 838 We quantified the fitness landscapes in the chemostat via the relationship between 839 metabolic strategies and the invasion growth rates of an invader adopting strategy 840 in a given nutrient environment . maximal growth rate of constitute the "maximal growth contour": 863 EHW ≔ X EHW X = }.
(S12) EHW is generally formed by many species, with each species adopting the maximizing 864 strategies EHW X corresponding to one environment X on the maximal growth 865 contour. 866 EHW is outside of the invasion zone for any species . (Otherwise, if a species could 867 invade an environment X on EHW , \•Ž a X ∈ EHW > , this would directly violate 868 the requirement by Eqs. (S11) and (S12) that max • \•Ž a | X = .) Therefore, the 869 necessary and sufficient condition for a set of species to be evolutionarily stable, is to 870 construct a steady-state environment on the maximal growth contour: 871 a * •žŸ ,$$ ∈ EHW . (S13) Therefore, a strategy belonging to the non-invasible set must be a maximizing strategy. 872 An example of maximal growth contour is shown in Fig 4D.  873  874 3. Non-invasible strategy: Nevertheless, adopting one of the maximizing strategies 875 along the maximal growth contour does not guarantee that a species will satisfy Eq.

876
(S13) and become non-invasible, as a maximizing strategy for environment + may end 877 up constructing a different environment 3 . To identify a non-invasible species for supply 878 condition $%&&'( , the flux-balance condition needs to be considered, with the strategies 879 maximized at each environment. This requirement forms a "maximal flux-balance curve" 880 in the nutrient space:  When there are discontinuous points on the maximal growth contour, there can be 892 "gaps" in the nutrient supply space, where no single strategy on the maximal growth 893 contour satisfies Eq. (S14) . Under this condition, more than one strategy are required 894 to co-create an environment on the conjunctions of discontinuous points of the maximal 895 growth contour. Therefore, discontinuous points of the maximal growth contour permit 896 evolutionarily stable coexistence, where * F&G contains more than one species. Two 897 examples of such discontinuities and coexistence are shown in Fig 5 and  When two nutrients are both essential for growth, such as nitrogen and phosphorus, 910 and both require a substantial allocation of resources for import, the system can be 911 abstractly modeled as shown in Fig 4A. In this metabolic model, we assume an exact 912 trade-off between the allocation of limited resources to import nutrient a or nutrient b. 913 The fraction of resources allocated to import nutrient a is represented by P , thus 914 leaving a fraction Q = 1 − P to import nutrient b. The import rate of nutrient is 915 assumed to follow the Monod equation as a function of nutrient concentration, and is 916 proportional to > : 917 (S16) Import of both nutrients is required for cell growth: 918 = • min( P ( ), Q ( )). (S17) For this model, for simplicity we do not explicitly consider intracellular metabolites. 919 Rather, import directly determines growth. 920 In this model, a "species" is defined by its value of P . 921 Nutrient limitation can be clearly quantified in this system: if P > Q ( ), the system is 922 limited by nutrient b; if P < Q ( ), the system is limited by nutrient a.

924
A species with the following parameters was used to generate Fig 1D and  In Fig 1E, to demonstrate how dilution rates may switch the nutrient, we used the same 931 supply condition as the blue condition in Fig 1D (  species Red has P = 0.65. In Fig 4C, we started with species Blue and species Red.
We then generated the fitness landscape for each species at the steady-state 938 environment it constructed, then chose the strategy that maximized invasion growth rate 939 for this fitness landscape to generate a new species, and iterated this process five 940 times. The species Black has P = 0.5. 941 In generating Fig 4D, Fig 2B- 3. Metabolic model with substitutable nutrients that require assimilation 962 In cells, the assimilation of imported raw material, such sugars, into biomass such as 963 proteins, takes multiple steps and enzymes and consumes a considerable amount of 964 energy. When the resources allocated to nutrient assimilation are considered, a cell's 965 strategy becomes more complex. A mathematical model involving three substitutable 966 nutrients a, b, c that need assimilation is shown in Fig 3A, with >+ represents the 967 fraction of resources allocated to importing nutrient i into internal metabolite, and >3 968 represents the fraction of resources allocated to assimilate the internal i into biomass. In 969 this model, the import rate has a similar form to the previous two models, 970 (S19) If two nutrients are both essential for growth, and a cell is able to convert one nutrient 987 into another albeit at a certain cost, as shown in Fig 5A,  To implement trade-offs, the sum of elements of = ( P , Q , PQ , QP ) is taken to be 994 equal to 1.

995
In this metabolic model, cells internalize nutrient a and nutrient b from the chemostat to 996 supply internal concentration of nutrients, P,\•GD³•H' and Q,\•GD³•H' . Meanwhile, the 997 internal nutrients can be converted into each other. Nutrients also diffuse in and out of 998 the cell passively with rate . Cell growth requires both internal nutrients, and depletes 999 them in a fixed proportion. 1000 In this model, the growth rate of a cell is taken to be: 1001 The net import rate, including passive diffusion, is: Cells N-limit P-limit, N 10-fold dilution P-limit, N 5-fold dilution P-limit, N 2-fold dilution P-limit C-limit P-limit, C 10-fold dilution P-limit, C 5-fold dilution P-limit, C 2-fold dilution P-limit 0.  concentrations (dots) for the three species in a three-dimensional nutrient space. Black 1336 curves with arrows show the system's limit-cycle trajectory. 1337 1338 C. The top panel shows the time course of species biomass for the limit cycle in (B). 1339 The bottom panel shows how the fitness landscape changes with time over one period 1340 of the oscillation. 1341 1342 1343 for import of two essential nutrients, with the lower of the two import rates determining 1348 growth rate. Species Red and species Blue allocate resources differently (indicated by 1349 parameter a a , see Methods Black, which places itself on the peak of its fitness landscape. The same procedure is 1362 also performed starting with species Red. 1363 1364