Species dynamics and interactions via metabolically informed consumer-resource models

Quantifying the strength, sign, and origin of species interactions, along with their dependence on environmental context, is at the heart of prediction and understanding in ecological communities. Pairwise interaction models like Lotka-Volterra provide an important and flexible foundation, but notably absent is an explicit mechanism mediating interactions. Consumer-resource models incorporate mechanism, but describing competitive and mutualistic interactions is more ambiguous. Here, we bridge this gap by modeling a coarse-grained version of a species’ true, cellular metabolism to describe resource consumption via uptake and conversion into biomass, energy, and byproducts. This approach does not require detailed chemical reaction information, but it provides a more explicit description of underlying mechanisms than pairwise interaction or consumer-resource models. Using a model system, we find that when metabolic reactions require two distinct resources we recover Liebig’s Law and multiplicative co-limitation in particular limits of the intracellular reaction rates. In between these limits, we derive a more general phenomenological form for consumer growth rate, and we find corresponding rates of secondary metabolite production, allowing us to model competitive and non-competitive interactions (e.g., facilitation). Using the more general form, we show how secondary metabolite production can support coexistence even when two species compete for a shared resource, and we show how differences in metabolic rates change species’ abundances in equilibrium. Building on these findings, we make the case for incorporating coarse-grained metabolism to update the phenomenology we use to model species interactions.

anaerobic-usually sugar consuming-metabolic lifestyle, and is the primary anaerobic energy-producing reaction for many microorganisms (Gottschalk, 1986). A signature of fermentation is that it results in byproducts such as organic acids, alcohols, and/or gases, which are produced due to the incomplete resource oxidation during the energy producing reactions. When organic acids (e.g., lactate) are produced, they are often used as resources by other microorganisms. In our model, we simplify the true biochemistry (see Supplemental, Fig. S1) by assuming glucose (R 1 ) and phosphate (R 2 ) enter the cell at uptake rates k i , biomass (N 1 ) is generated at rate β M 1 M 2 , and the metabolic byproduct, lactate, is exported from the cell at rate ν 3 . In addition, we assume that any unused glucose and phosphate are exported at rates ν 1 and ν 2 .
production of a byproduct, lactate, which is exported back into the environment (i.e., excretion); therefore, efficiency is 100 emergent property determined by the balance between the biomass and energy production. 101 Here, we assume that cellular metabolism relies on the interaction of sugar and phosphate, producing new biomass and a byproduct (i.e., lactate).
R 1 here is the sugar, R 2 represents a source of phosphorus, and R 3 is the metabolic byproduct lactate. We model this 102 intracellular reaction as occurring at rate β . Next, we assume that the organism grows in a chemostat-type environment, 103 where resource inflow rates are held constant. In addition, we assume that metabolic rates are not limited by any other 104 factors and that resources are not toxic at high concentrations and therefore do not inhibit growth. Given that we consider 105 one phosphate molecule (i.e., we assume 1:1 stoichiometry), this already simplifies the resource requirements relative 106 to the true process. But we will see general principles emerge. 107 We now assume that intracellular processes are always close to equilibrium, so that we can balance fluxes for the 108 internal cellular densities of the two input resources (labeled M 1 and M 2 ). This leads to: for uptake rates k 1 and k 2 , which can depend in an arbitrary way on external resource concentrations R 1 , R 2 , and R 3 , 110 outflow (i.e., export) rates ν 1 , ν 2 , and ν 3 , and biomass production ν bio (Fig. 1). We make a natural assumption that 111 export of each resource is determined by passive excretion: i.e. that ν α = λ α M α for each metabolite, where λ α is a 112 species-and resource-specific constant. This resource excretion, or leakiness, prevents unbounded internal resource 113 accumulation, but it can also mediate species interactions (e.g., cross-feeding) and therefore structure communities 114 (Pfeiffer and Bonhoeffer, 2004;Morris, 2015;Sun et al., 2019). As such, the outflow rates (ν α ) represent use it or lose 115 it internal resource dynamics and therefore intracellular resources that are not used are exported back into the external 116 environment.

120
Using these intracellular concentrations, we can now develop a set of generalized consumer-resource equations using 121 uptake rates and this simplified flux balance analysis (Eq. 4) as the only building blocks. First, we focus on one species 122 and model its uptake rate k α of resource α as C α1 R α , where R α is the external (environmental) concentration of this 123 resource. Applying this assumption, we derive a general set of equations for the consumer and three resources: where ρ α and η α are inflow and outflow rates for each of the three resources, and µ 1 is the mortality rate of the consumer 125 (see Table 1). For the resources (Eqs. 5a-5c), the extracellular resource concentrations are determined by the inflow and 126 outflow rates minus the consumer-density dependent components which include uptake, leakage, and conversion into 127 biomass (Eq. 4). For the consumer (Eq. 5d), the production of new biomass is determined by the density dependent 128 uptake, leakage, and conversion of resources into biomass (Eq. 4) minus the density dependent mortality. Here, we have 129 assumed that uptake, k α , increases linearly with resource concentration, but this could easily be modified to saturate 130 using Monod dynamics by changing k α to (C α1 R α )/(K α1 + R α ), where K, the half-saturation constant, is a species-and 131 resource-specific constant.

132
We now note two limits of the simplified flux balance analysis (Eqs 3-4). The first limit is where β (i.e., the internal 133 metabolic rate) is large relative to the other rates, Where O(1/β ) represents the linear approximation of the higher order terms of the Taylor Series expansion. Note that 135 we have to keep the O(1/β ) terms for M 1 and M 2 , because for at least one of the two (depending on whether k 1 > k 2 or 136 k 1 < k 2 ) the O(1/β ) term vanishes in this limit of fast reaction rate β . Therefore, when β is large relative to the other 137 rates the consumer-resource equations become: Hence, we recover Liebig's law (i.e. the net growth rate of the consumer is min(R 1 C 11 , R 2 C 21 )). Using these equations 139 in numerical simulations, we observe that the consumer abundance saturates ( Fig. 2A), and the final abundance depends 140 on the inflow rate of the more limiting resource (Fig. 2B). As such, the model in this limit behaves similar to classical 141 consumer-resource models. In addition, we generate novel terms for the byproduct production rate, which in this fast 142 reaction rate limit is N 1 min(R 1 C 11 , R 2 C 21 ), and we find corresponding equations for the uptake and export of glucose 143 and phosphate.  We used numerical simulations of our single-species model to determine how the resource landscape changes population dynamics. In our model, N 1 is the population, R 1 and R 2 are the consumed resources (i.e., sugar and phosphate), R 3 is a metabolic byproduct, ρ is the resource inflow rate, and β is the internal metabolic rate of the consumer. A & B: In the limit where β is high relative to other rates (Eq. 7), we recover Liebig's Law of the Minimum. Consumer abundance saturates, but the final abundance depends on the inflow rate, ρ α , of the more limiting resource. When the inflow of R 2 (ρ 2 ) is lower than the inflow of R 1 (ρ 1 ), R 2 becomes depleted. At saturation, R 1 is in excess while R 2 is limiting (A). Changes in the inflow rates alters the final abundance (B), but the final abundance will always be determined by the more limiting resource. When ρ 1 = 100 (blue) the final abundance will increase until ρ 2 = 100 (blue dashed line), and when ρ 1 = 108 (red) the final abundance will increase until ρ 2 = 108 (red dashed line). In contrast, in the limit where β is low relative to other rates (Eq. 9), we recover multiplicative co-limitation. Consumer abundance still saturates, albeit at a lower final abundance given the same resource inflow rates and with a growth lag-phase (C). As consumer abundance increases, both R 1 and R 2 are depleted. In addition, the final abundance increases as inflow rates increase (D). When ρ 1 = 100 (blue) the final abundance will continue to increase past ρ 2 = 100, and when ρ 1 = 108 (red) the final abundance will continue to increase past ρ 2 = 108. Simulations based on Eq. 5 using the following parameters: η α = 0.01, C α1 = 0.01, µ 1 = 0.01, λ α = 0.1, β 1 = 10 (A) or 10 −6 (B), ρ 1 = 100 or 108, and ρ 2 varies between 80 and 120 (A & C: ρ 1 = 100, ρ 2 = 80). Resource concentrations are expressed in cell equivalents (i.e., moles required to generate a cell) and we assume 1:1 resource stoichiometry.
The second limit is where β (i.e., the internal metabolic rate) is small relative to the other rates, so that the consumer-resource equations become: Hence, we recover multiplicative co-limitation by the two resources (i.e., the net growth rate of the consumer is ). Using numerical simulations, we observe that consumer abundance, N 1 saturates as expected but now 148 includes a growth lag-phase (Fig. 2C). The lag-phase (i.e., period with no net population growth: dN dt = 0) occurs while 149 β < µ 1 λ 1 λ 2 C 11 R 1 C 21 R 2 . When R 1 and R 2 accumulate ( Fig. 2C), the denominator increases and the inequality flips which yields 150 population growth. Due to the parameters used in the simulation above, the equation for the lag-phase simplifies to 151 β = 1 R 1 R 2 ; therefore, our simulation demonstrates how the lag-phase is controlled by an interaction between the resource 152 environment and the internal metabolic rates. Finally, the final abundance depends on the inflow rate of both resources 153 ( Fig. 2D). An increase in the inflow rate of either resource will yield a higher final population abundance.

154
In summary, from this coarse-grained representation, we recover two classic outcomes of consumer-resource theory 155 by taking limits of the internal reaction rate β . We can also generalize these classic limits for intermediate β , in a 156 way that is not currently used in consumer-resource models and falls neither into the category of Liebig's Law nor 157 multiplicative co-limitation. Finally, we find functional forms for the production rate of metabolic byproducts (and 158 excretion of other resources) in each of these limits and in between. Our model therefore demonstrates how we can 159 generalize the functional form of consumer-resource models by considering realistic, simplified intracellular processes. In the two species model, two species interact through competition and facilitation. For each species, we use a simplified metabolic networks to describe resource uptake rates (k αi ), biomass generation rates (β i M αi M αi ), and metabolic byproduct export rates (ν αi ) for species i and resource α. Species N 1 consumes glucose (R 1 ) and phosphate (R 2 ) and produces lactate (R 3 ). Species N 2 competes for phosphate (R 2 ) but also uses the lactate produced by species N 1 . In addition, all internal resources (byproducts and unused resources) are exported at rate ν α (not shown). As such, facilitation is modeled by cross-feeding (i.e., production and consumption of lactate) between species N 1 and N 2 .
The Two Species Model. We now modify the model above to incorporate a second species. Here, both species use R 2 161 (phosphate), but the second species, N 2 , uses a combination of R 2 and R 3 (lactate) to generate new biomass (Fig. 3).

162
While we are using this as a model with both competition (e.g., shared resources) and facilitation (e.g., metabolic 163 cross-feeding) interactions, it also represents the natural cross-feeding interaction between organisms. For example, the 164 model could represent the interactions between the lactate producing (e.g., Bifidobacterium adolescentis) and lactate 165 consuming (e.g., Eubacterium hallii) bacteria found in human and animal digestive systems (Duncan et al., 2004).

171
In our model, productivity (and efficiency) increase due to the second species using the metabolic byproducts of the 172 first (Fig. 3). As such, our model demonstrates how competition (i.e., shared resources) and facilitation (i.e., production 173 and consumption of intermediate resources) mediate species dynamics and coexistence conditions and can be used to 174 understand natural and engineered microbial systems.

175
First, we define two distinct internal metabolic processes, one for each consumer species: Consumer N 2 may also produce a metabolic byproduct, but we have not included such a process here because we are 177 focusing on competition for R 2 and facilitation through the production of R 3 by species N 1 (Fig. 3). Importantly, we 178 now have two internal reaction rates, β 1 and β 2 . Here we focus on how these rates, both relative to each other and also 179 to the other rates in the model, affect species coexistence. This approach allows us to use the coarse-grained metabolic 180 processes as the mechanisms underlying species interactions. Furthermore, it allows us to explore the potential for 181 changes in species coexistence due to metabolic (i.e., reaction rates) and landscape (i.e., inflow rates) factors.

182
Generalizing the approach in the previous section, and again assuming internal flux balance, we can define the fluxes 183 for an individual belonging to species N 1 as: for internal concentrations M α1 and uptake rates k α1 . While for an individual of species N 2 we have: for internal concentrations M α2 and uptake rates k α2 . We also assume that all resources can be excreted from both con-186 sumers, but to simplify the model slightly we will assume that the specific export rates are equal, λ , broadly consistent 187 with passive diffusion across sufficiently similar cell wall types. Based on these assumptions, we can solve for internal 188 equilibrium in both cell types to obtain: where, similarly to the one species case, the function F(x, a, b) =  192 To focus on the effects of internal reaction rates and the resource landscape, we will further simplify our model by 193 making a few assumptions. We will assume that the outflow rates for each resource are the same, so that η i = η, and 194 that the per capita mortality rates for each consumer are equal, so that µ i = µ. We will also again assume that the per 195 capita uptake rate of resource α by species i can be written as C αi R α . We will then determine the effects of internal 196 reaction rate by independently changing β (i.e., the internal metabolic rate) for each species. In addition, we will change 197 the landscape conditions by exploring the inflow rates for each resource ρ α . Given these assumptions, our general two 198 species consumer-resource model is: where ρ α and η α are inflow and outflow rates for each of the three resources, and µ 1 is the mortality rate of the consumer 200 (see Table 1). For the resources (Eqs. 14a-14c), the extracellular resource concentrations are determined by the inflow 201 and outflow rates minus the consumer-density dependent components of each consumer which include uptake, leakage, 202 and conversion into biomass (Eq. 13). For the consumers (Eqs. 14d & 14e), the production of new biomass is determined 203 by the density dependent uptake, leakage, and conversion of resources into biomass (Eq. 13) minus the density dependent 204 mortality.

205
Using numerical simulations, we model the resource and consumer dynamics and determine equilibrium conditions 206 (Figs. 4 & 5). First we consider when the internal metabolic rates, β i , are the same. When internal metabolic rates are 207 both high, species coexist at a density determined by the shared resource inflow rate (i.e., ρ 2 ) until the inflow rate of the 208 unshared resource, R 1 , exceeds the inflow rate of the shared resource, R 2 . When ρ 1 is greater than ρ 2 (i.e., ρ 1 /ρ 2 > 1),

209
species N 1 will outcompete species N 2 for R 2 , and N 2 will become rare (Fig. 4A). In our simulations, species N 2 remains 210 in the system with rare abundance due to resource leakage from species N 1 , but would likely go extinct due to stochastic 211 fluctuations, and so coexistence in this limit is somewhat tenuous. These findings expand the expectations of Liebig's 212 Law to two cross-feeding species and demonstrate both species dynamics and equilibrium abundances across various 213 resource landscapes. In contrast, when internal metabolic rates are low, species relative abundances are determined by 214 their respective required resources based on multiplicative co-limitation and both species now demonstrate growth lag-215 phases. These findings expand the expectations of multiplicative co-limitation to two cross-feeding species. However, 216 since species N 1 will not be resource limited when ρ 1 is greater than ρ 2 , then species N 2 will maintain a higher relative 217 abundance across wider resource landscape (Fig. 4B, see Eq. 9). Together, we find that if internal rates are the same 218 but either high or low compared to the other rates in our model, then our results expand the expectations of Liebig's 219 Law and multiplicative co-limitation to the two-species system with a metabolic dependency. Between these limits, we 220 find a smooth transition between the two scenarios as coexistence conditions become decoupled from classical resource 221 limitation predictions when the internal metabolic rates decrease (Fig. 4C). In addition, we find that coexistence depends 222 on both the internal metabolic rates and the resource inflow rates even when uptake rates are the same. Furthermore, 223 while species N 2 depends on species N 1 for the production of R 3 , external inputs of R 3 (i.e., ρ 3 > 0) alone are not enough 224 to allow species N 2 to outcompete species N 1 but will result in a lower final density (see Supplemental).

225
Finally, we consider the dynamics and equilibria when the internal metabolic rates, β i , are different. We find that, 226 when the rates differ the outcome depends on which species has the higher metabolic rate. When the byproduct producer,

227
N 1 , has the higher rate, then the results are similar to when both species have high internal metabolic rates (Fig. 5A). 228 We do note, however, two important differences: 1) species N 2 alone exhibits a growth lag-phase, and 2) both R 1 and  Fig. 4 Two Species with the Same Metabolic Rates. We simulated species dynamics in the two-species model with equal internal metabolic rates, β i . In this model, two species compete for a shared resource (R 2 ), but N 1 consumes R 1 and produces R 3 as a metabolic byproduct which is consumed by N 2 . When metabolic rates are high, species abundances saturate and species can coexist provided the appropriate resource conditions (A). When the shared resource inflow rate (ρ 2 ) is higher than the inflow rate of R 1 (ρ 1 ), then the two species coexist at a relatively high abundance because N 1 is more limited by its exclusive resource (R 1 ). However, when ρ 1 is higher than ρ 2 (i.e., ρ 1 /ρ 2 > 1), then N 1 can outcompete N 2 for the shared resource (R 2 ) and the final abundance of N 2 will be reduced. Even when the shared resource inflow rate is increased (ρ 2 = 100 vs. 80), the final abundance of N 2 does not increase when ρ 1 /ρ 2 > 1. When metabolic rates are low, species still coexist, but the range of coexistence is greatly expanded (B). While the final abundance of N 2 still decreases as the resource inflow ratio (ρ 1 /ρ 2 ) approaches 1, it does not decrease at the same rate and the final abundance remains relatively high even when the resource inflow ratio is > 1.
conditions and high relative abundances for both species are greatly expanded. In fact, we find coexistence with moderate 231 abundances along all inflow rates tested and the final abundances of both species are determined only by the shared 232 resource, R 2 (Fig. 5B). In addition, we find that both species exhibit growth lag-phases and that R 2 and R 3 are now 233 the depleted resources. Between these limits, we find that the transition from classical resource limitation prediction 234 to complete coexistence depends on both metabolic rates decreasing (Fig. 5C). These findings highlight how variation 235 in both the internal metabolic reaction rates and in the environmental conditions can influence species interactions and 236 change coexistence expectations.

237
Discussion 238 Here, we reformulated the classic consumer-resource model by independently including resource uptake, internal metabolic 239 rates, and byproduct export. We found that when internal metabolic dynamics are included -in addition to uptake -two 240 common models of resource limitation (Liebig's Law and multiplicative co-limitation) appear in particular limits of the 241 internal reaction rates. In this model, we make predictions for the functional form of the production of metabolic byprod-242 ucts which can be used in metabolic exchanges, and we balance the requirements for growth, energy, storage, and export.   in the model. The species compete for R 2 , and species N 1 consumes R 1 and produces R 3 which is used by species N 2 . When N 1 has the higher metabolic rate, the species dynamics are similar and the coexistence conditions are identical to the model where both species have high metabolic rates (A). However, there are now two important differences. First, N 2 (gray line) has a growth lag-phase. Second, both R 1 and R 2 are now depleted as the population abundance saturates. When N 2 has the higher metabolic rates, the coexistence conditions are drastically different (B). N 1 and N 2 have growth lag-phases, and R 2 and R 3 are now the depleted resources. In addition, both species coexist at a final abundance determined only by the shared resource inflow rate (ρ 2 ) across all resource inflow rates (ρ 1 /ρ 2 ) tested. When we independently change β 1 and β 2 , the final abundances depend on both values (C). As β 1 decreases from 1 to 10 −4 (blocks), the final abundance of N 1 decreases and N 2 increases for a given resource inflow ratio. In addition, as B 2 increases from 10 −4 to 1 (within each block), the final abundance of N 1 decreases more and N 2 increases for a given inflow ratio. Simulations based on Eqs. 14a-14e using the following parameters: η α = 0.01, C αi = 0.01, µ i = 0.01, λ α = 0.1 and β i = 10 (high) or 10 −6 (low), ρ 2 = 80 or 100, and ρ 1 /ρ 2 varies between 0.5 and 1.7 (A le f t & B le f t : ρ 1 = 64, ρ 2 = 80).
a range of interactions between species but are unable to generalize across environmental variation because they do not 259 provide an unambiguous way to include the resource landscape. While it may be possible to modify per capita growth 260 rates to be a function of environmental conditions (e.g., temperature) and the resource landscape, how to incorporate 261 these effects in the species interaction rates is ambiguous, at best. However, these models allow positive and nega-262 tive species interactions to be explored straightforwardly, with far-reaching implications for coexistence and stability.

263
Consumer-resource models explicitly include the interaction between the resource landscape and consumers, but at the 264 expense of introducing more explicit mechanism, and therefore more choices in the way interactions are implemented metabolic rates may also reveal how environmental conditions that can alter metabolic rates (e.g., temperature) will 272 contribute to species dynamics.

273
Our model used a simplified metabolic network to describe the internal physiology of cells and this allowed us to 274 determine how uptake, growth, and export contribute to the dynamics and coexistence of interacting species. Specif-275 ically, we were able to incorporate competition via shared resource consumption and facilitation via the production 276 and consumption of a metabolic byproduct (i.e., cross-feeding). Our model accomplishes this by using a consumer-277 resource framework with multiple inputs and outputs and by explicitly incorporating internal resource dynamics (see both are able to produce limiting resource dynamics based on Liebig's Law. However, our model was also able to pro-281 duce multiplicative co-limitation, and it allowed us to explore metabolic cross-feeding by providing a mechanism for 282 the production of cellular byproducts. We accomplish this by allowing internal resources, including luxury uptake and 283 metabolic byproducts, to leak out of cells, which establishes strong conceptual similarities to Droop's model (Droop, 284 1974). Because our model uses this internal flux balance approach, efficiency becomes an emergent property rather than 285 a fixed term. This contrasts with approaches that include fixed efficiency terms that account for the internal processes 286 without providing physiological mechanism (Moore et al., 1993;Schimel and Weintraub, 2003;Moore et al., 2005).
287 Likewise, our model differs from others that include explicit resource transfer pathways and corresponding efficiencies 288 between species (Wiegert and Owen, 1971;DeAngelis et al., 1975;De Ruiter et al., 1995;Moore et al., 2005) because 289 the byproducts and internal resource excesses are returned to the environment as 'public goods' rather than transferred 290 directly. As such, our model uniquely allows us to incorporate a mixture of species interactions in addition to exploring 291 how each population responds to a complex resource landscape because we include a simplified metabolic network.

292
Our model allows us to address scenarios that are unique to microbial communities -metabolic cross feeding and in-293 ternal cellular metabolic rates -and this provided novel insights. First, by providing a mechanism for the production and 294 subsequent consumption of metabolic byproducts (i.e., metabolic cross-feeding) our model provides a way to explore 295 the role of facilitation interactions in a consumer-resource framework. Cross-feeding is a widely observed phenomenon 296 in microbial communities (Morris et al., 2012;Morris, 2015;D'Souza et al., 2018). However, while cross-feeding ap-297 pears to be ubiquitous, there are very few hypotheses or theoretical prediction for how it will shape the dynamics of 298 interacting species. Our model provides a novel approach to explore this phenomenon and how it contributes to the 299 diversity, complexity, and dynamics of microbial communities. In our model, the coexistence of species N 2 is dependent 300 on the production of a metabolic byproduct, R 3 , from species N 1 . Assuming there is no external input of the byprod-301 uct, the coexistence of N 2 is dependent on N 1 . However, due to the way interactions are formulated, external inputs of 302 byproduct alone are not enough to let N 2 outcompete N 1 (see Supplemental). As such, our model provides a new tool 303 to explore the benefits and consequences of cross-feeding in microbial communities. Second, by separating the uptake 304 and internal metabolic rates, our model is able to produce unique species dynamics and it recapitulates two common 305 models of resource limitation (Liebig's Law and multiplicative co-limitation). One major finding based on our model is 306 that internal metabolic rates can have a strong impact on the dynamics of interacting species. Specifically, changes in the