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.


Introduction
A central goal in ecology is to understand and predict the dynamics in communities of interacting species (Holt 1977;Loreau 2010;Vellend 2010Vellend , 2016. Mathematical models allow us to generate and test theoretical predictions, and the development of such models leads to a hierarchy of challenges. The first challenge is determining an appropriate functional form describing species dynamics. A range of functional forms with increasing complexity has been used and each has strengths and weaknesses which are often context dependent (Holling 1959;Abrams 1982;DeAngelis et al. 1989;Murdoch et al. 2003;Mougi and Kondoh 2012). Second, we need to parameterize these equations, for example by quantifying the strength and sign of species interactions in a given environmental context. While some attempts have been made to parameterize natural systems, accurately fitting interaction strengths remains challenging in both empirical and theoretical work (Schoener 1983;Tilman 1987;Ives et al. 2003;Marino et al. 2014;Carrara et al. 2015;Terry et al. 2017; Barner et al. 2018). Finally, we may wish to determine how changes in the environmental context will modify species interactions, dynamics, and even coexistence. Integrating each of these goals will lead to the development of robust models which can predict the dynamics of complex communities, even when the environmental landscape changes within and across ecosystems.
There are two primary approaches to mathematically model the dynamics of interacting species. The Lotka-Volterra equations characterize species interactions in terms of the net, direct effect of one population on another's growth rate (Lotka 1932;Volterra 1926, see Supplemental). However, because these models lack an explicit description of the mechanisms mediating interactions, it is often difficult to translate the inferred interactions in different environmental contexts (Abrams 1983;Grilli et al. 2017). For example, if two species compete, it is often because they consume common resources (Gause and Witt 1935;MacArthur 1970;Schoener 1983). But, these models assume that resource dynamics can be safely ignored because resource dynamics are faster than consumer dynamics (MacArthur 1970). This exposes an important context-dependence of Lotka-Volterra type equations: the strength and even the sign of a pairwise interaction may depend on what resources are present (Xiao et al. 2017). An alternate approach, consumer-resource equations, characterizes competitive interactions as the explicit result of shared resource consumption (Tilman et al. 1982;Grover 1990;Litchman 2003;Abrams 2009, see Supplemental). These models produce species interactions as an emergent property dependent on shared resource consumption. Assuming we can infer or otherwise estimate consumer preferences, one critical aspect of the environmental context is now explicitly characterized, via resource input rates. As such, species dynamics across resource landscapes can be understood better than in the case of Lotka-Volterra, where the effect of the environment is implicit (Tilman 1977;Grover 1990Grover , 2011. The more explicit mechanism in a consumer-resource model gives an advantage, but there is a cost. Lotka-Volterra equations are extremely flexible and can straightforwardly incorporate a mixture of antagonistic and mutualistic interactions simply by altering the signs of entries in the community matrix (Mougi and Kondoh 2012;Allesina and Tang 2012). But for consumer-resource models we have to be careful about how those mechanisms are formulated. Consumption can take a variety of forms depending, for example, on whether resources are substitutable or essential (Tilman 1980), and mutualistic interactions can occur via resource production and exchange (i.e., cross-feeding) (Loreau 2001;Freilich et al. 2011;Zelezniak et al. 2015;Sun et al. 2019). The latter in particular is underexplored relative to consumption (Butler and O'Dwyer 2018), and an open question is to what extent the precise form of exchange might affect community dynamics and species coexistence.
So how do we retain the advantages of consumerresource models, but also incorporate a mixture of species interactions? Current consumer-resource models are largely agnostic to what happens inside cells or organisms. For many systems, this approach may be valid especially when the resources (i.e., prey) belong to higher trophic levels, self-regulate, and/or closely match the stoichiometric requirements of the consumer (Sterner and Elser 2002;Cherif and Loreau 2007;Hall 2009). However, this may not be valid for microorganisms. First, most microorganisms consume abiotic resources which do not self-regulate. Second, any single resource generally does not satisfy the full stoichiometric requirements of metabolism. For example, heterotrophic microorganisms require organic carbon, but they still require nitrogen, phosphorus, and other resources to grow and reproduce. Because growth depends on multiple resources, population dynamics may depend on a single limiting resource (e.g., Liebig's Law of the Minimum: von Liebig and Gregory 1842; Odum 1959) or an interaction between resources (e.g., multiplicative colimitation: Harpole et al. 2011). Likewise, the consumption and transformation of resources depends on how cells produce biomass, energy, and metabolic byproducts, and species interactions may therefore depend on metabolic rates and the release of metabolic byproducts.
Here, we propose that going one level deeper into cellular metabolism will allow us to generalize consumer-resource models in a meaningful way and give them the same flexibility to describe multiple interaction types as models of direct species interactions. First, we introduce a general metabolically informed consumer-resource model for a single species. In this model, we include internal cellular dynamics via a simplified metabolic network, which requires less knowledge of a particular species' idiosyncrasies but still captures the major metabolic events transforming resources. Using this model, we explore population dynamics in variable resource environments. Next, we expand the approach to incorporate a second species. Here, both species share a common resource, but the second species also uses a metabolic byproduct produced by the first species. As such, our model includes competition for the shared resource and facilitation via metabolic cross-feeding. We then explore how the resource landscape and internal metabolic rates interact to demonstrate how competition and facilitation mediate species dynamics and coexistence conditions under different resource use scenarios. In summary, our findings demonstrate how a metabolically informed consumer-resource model for a single species can produce both Liebig's Law and multiplicative co-limitation in particular limits of the metabolic rates. Furthermore, when a second species is included we expand the expectations of Liebig's Law and multiplica-tive co-limitation to include cross-feeding and demonstrate both species dynamics and equilibrium abundances across resource landscapes.

Methods and results
The metabolically informed consumer-resource model We know a substantial amount about the internal physiology of microorganisms, and there have been large advancements in the development of flux balance models which use biochemistry and genomics to describe (to some level of approximation) every reaction that occurs within a cell (Kauffman et al. 2003;Orth et al. 2010). Likewise, models based on the dynamic energy budget theory provide a complete mass and energy balance approach for molecular processes, cellular physiology, and population growth (Kooijman 2001;Nisbet et al. 2008). More recently, flux balance models have been applied in the context of entire microbial communities and their interactions (Embree et al. 2015;Zomorrodi and Segrė 2016;Pacheco et al. 2019). However, we propose that including a full description of a metabolic network may not be required to develop a useful ecological model. Simplified metabolic models have been used to understand metabolic partitioning (Kempes et al. 2012) and trophic tradeoffs (Chakraborty et al. 2017). Furthermore, models that add some intracellular dynamics, such as Droop's model which accounts for cell quotas and intracellular resource storage (Droop 1974), have been used to better understand how consumers contend with nutrient limitations and how trait trade-offs underlie coexistence (Cherif and Loreau 2010;Litchman et al. 2015). Here, we model the internal dynamics of a microorganism via simplified metabolic networks, which require less knowledge of the particular species' idiosyncrasies but still capture the major metabolic events transforming resources. Our prototypical metabolic model is based on a basic fermentation reaction, homolactic fermentation, which uses glucose and phosphate and produces lactate ( Fig. 1 & S1). In this reaction, glucose and phosphate are used to form new biomass, and glucose is used to produce chemical energy via fermentation. The energetic component results in the production of a byproduct, lactate, which is exported back into the environment (i.e., excretion); therefore, efficiency is emergent property determined by the balance between the biomass and energy production.
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 intracellular reaction as occurring at rate β. Next, we assume that the organism grows in a chemostat-type environment, where resource inflow rates are held constant. In addition, we assume that metabolic rates are not limited by any other factors and that resources are not toxic at high concentrations and therefore do not inhibit growth. Given that we consider one phosphate molecule (i.e., we assume 1:1 stoichiometry), this already simplifies the resource requirements relative to the true process. But we will see general principles emerge. We now assume that intracellular processes are always close to equilibrium, so that we can balance fluxes for the 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 , outflow (i.e., export) rates ν 1 , ν 2 , and ν 3 , and biomass production ν bio (Fig. 1). We make a natural assumption that export of each resource is determined by passive excretion: i.e., that ν α = λ α M α for each metabolite, where λ α is a species-and resource-specific constant. This resource excretion, or leakiness, prevents unbounded internal resource accumulation, but it can also mediate species interactions (e.g., cross-feeding) and therefore structure communities (Pfeiffer and Bonhoeffer 2004;Morris 2015;Sun et al. 2019). As such, the outflow rates (ν α ) represent use it or lose it internal resource dynamics and therefore intracellular resources that are not used are exported back into the external environment. We can solve this system of polynomial equations (2) for internal resource concentrations to obtain: while biomass production is given by: and where F (β, k 1 , k 2 ) = k 1 + k 2 + λ 1 λ 2 β 2 − 4k 1 k 2 depends on the uptake, export, and reaction rates.  Fig. 1 Conceptual model. Our metabolically informed consumerresource model uses fermentation as a prototype. is an anaerobicusually 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 Using these intracellular concentrations, we can now develop a set of generalized consumer-resource equations using uptake rates and this simplified flux balance analysis (4) as the only building blocks. First, we focus on one species and model its uptake rate k α of resource α as C α1 R α , where R α is the external (environmental) concentration of this 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 (see Table 1). For the resources (5a-c), the extracellular resource concentrations are determined by the inflow and outflow rates minus the consumer-density dependent components which include uptake, leakage, and conversion into biomass (4). For the consumer (5d), the production of new biomass is determined by the density dependent uptake, leakage, and conversion of resources into biomass (4) minus the density dependent mortality. Here, we have assumed that uptake, k α , increases linearly with resource concentration, but this could easily be modified to saturate using Monod dynamics by changing k α to (C α1 R α )/(K α1 + R α ), where K, the half-saturation constant, is a species-and resource-specific constant. We now note two limits of the simplified flux balance analysis (3-4). The first limit is where β (i.e., the internal 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 we have to keep the O(1/β) terms for M 1 and M 2 , because Note: resource concentrations are given in cell equivalents -moles required to generated one cell for at least one of the two (depending on whether k 1 > k 2 or k 1 < k 2 ) the O(1/β) term vanishes in this limit of fast reaction rate β. Therefore, when β is large relative to the other 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 in numerical simulations, we observe that the consumer abundance saturates (Fig. 2a), and the final abundance depends on the inflow rate of the more limiting resource (Fig. 2b). As such, the model in this limit behaves similar to classical consumer-resource models. In addition, we generate novel terms for the byproduct production rate, which in this fast 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 and phosphate.
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 βC 11 R 1 C 21 R 1 λ 1 λ 2 ). Using numerical simulations, we observe that consumer abundance, N 1 saturates as expected but now includes a growth lag-phase (Fig. 2c). The lag-phase (i.e., period with no net population growth: dN dt = 0) occurs while β < μ 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  Multiplicative Colimitation a b c d Fig. 2 Single species simulation. 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 (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 (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 yields population growth. Due to the parameters used in the simulation above, the equation for the lag-phase simplifies to β = 1 R 1 R 2 ; therefore, our simulation demonstrates how the lag-phase is controlled by an interaction between the resource environment and the internal metabolic rates. Finally, the final abundance depends on the inflow rate of both resources (Fig. 2d). An increase in the inflow rate of either resource will yield a higher final population abundance.
In summary, from this coarse-grained representation, we recover two classic outcomes of consumer-resource theory by taking limits of the internal reaction rate β. We can also generalize these classic limits for intermediate β, in a way that is not currently used in consumer-resource models and falls neither into the category of Liebig's Law nor multiplicative co-limitation. Finally, we find functional forms for the production rate of metabolic byproducts (and excretion of other resources) in each of these limits and in between. Our model therefore demonstrates how we can generalize the functional form of consumer-resource models by considering realistic, simplified intracellular processes.

The two species model
We now modify the model above to incorporate a second species. Here, both species use R 2 (phosphate), but the second species, N 2 , uses a combination of R 2 and R 3 (lactate) to generate new biomass (Fig. 3). While we are using this as a model with both competition (e.g., shared resources) and facilitation (e.g., metabolic cross-feeding) interactions, it also represents the natural cross-feeding interaction between organisms. For example, the model could represent the interactions between the lactate producing (e.g., Bifidobacterium adolescentis) and lactate consuming (e.g., Eubacterium hallii) bacteria found in human and animal digestive systems (Duncan et al. 2004). Likewise, other examples of these interactions would include the exchange of vitamin precursors and amino acids by auxotrophic organisms (Carini et al. 2014;Embree et al. 2015). These metabolic cross-feeding interactions are common in microbial systems (Pfeiffer and Bonhoeffer 2004;Morris et al. 2012;Mee et al. 2014;Morris 2015;Tasoff et al. 2015;Zelezniak et al. 2015;Sun et al. 2019) and have industrial applications (Jiao et al. 2012). Cross-feeding has been shown to promote species diversity and increase net ecosystem production (Morris et al. 2012;Pande et al. 2014). In our model, productivity (and efficiency) increase due to the second species using the metabolic byproducts of the first (Fig. 3). As such, our model demonstrates how competition (i.e., shared resources) and facilitation (i.e., production and consumption of intermediate resources) mediate species dynamics and coexistence conditions and can be used to understand natural and engineered microbial systems.
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 focusing on competition for R 2 and facilitation through the production of R 3 by species N 1 (Fig. 3). Importantly, we now have two internal reaction rates, β 1 and β 2 . Here we focus on how these rates, both relative to each other and also to the other rates in the model, affect species coexistence. This approach allows us to use the coarsegrained metabolic processes as the mechanisms underlying species interactions. Furthermore, it allows us to explore the potential for changes in species coexistence due to metabolic (i.e., reaction rates) and landscape (i.e., inflow rates) factors.  Fig. 3 Two species model. 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 Generalizing the approach in the previous section, and again assuming internal flux balance, we can define the fluxes 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 consumers, but to simplify the model slightly we will assume that the specific export rates are equal, λ, broadly consistent with passive diffusion across sufficiently similar cell wall types. Based on these assumptions, we can solve for internal equilibrium in both cell types to obtain: where, similarly to the one species case, the function Finally, we can put all of this together to generate a set of equations for the dynamics of both populations and the external concentrations of all three resources. To focus on the effects of internal reaction rates and the resource landscape, we will further simplify our model by making a few assumptions. We will assume that the outflow rates for each resource are the same, so that η i = η, and that the per capita mortality rates for each consumer are equal, so that μ i = μ. We will also again assume that the per capita uptake rate of resource α by species i can be written as C αi R α . We will then determine the effects of internal reaction rate by independently changing β (i.e., the internal metabolic rate) for each species. In addition, we will change the landscape conditions by exploring the inflow rates for each resource ρ α . Given these assumptions, our general two 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 (see Table 1). For the resources (14a-c), the extracellular resource concentrations are determined by 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. In this simulation, the resource inflow ratio needs to be greater than around ∼ 1.5 before N 1 outcompetes N 2 and N 2 becomes rare. Between these metabolic rates, there is a smooth transition (c). When both resources have equal inflow rates, the final abundance of N 1 decreases and N 2 increases as β drops below 1. Simulations based on Eq. 14a-e using the following parameters: η α = 0.01, C αi = 0.01, μ i = 0.01, λ α = 0.1, β i = 10 (high) or 10 −6 (low), ρ 2 = 80 or 100, and ρ 1 /ρ 2 varies between 0.5 and 1.  5 Two species with different metabolic rates. We simulated species dynamics in the two-species model where the species have different internal metabolic rates, β i . In this model, one of the species has high and the other low metabolic rates compared to the other rates 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 ratios (ρ 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 Eq. 14a-e 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 left & b left: ρ 1 = 64, ρ 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 inflow and outflow rates minus the consumer-density dependent components of each consumer which include uptake, leakage, and conversion into biomass (13). For the consumers (14d & e), the production of new biomass is determined by the density dependent uptake, leakage, and conversion of resources into biomass (13) minus the density dependent mortality.
Using numerical simulations, we model the resource and consumer dynamics and determine equilibrium conditions (Figs. 4 & 5). First we consider when the internal metabolic rates, β i , are the same. When internal metabolic rates are both high, species coexist at a density determined by the shared resource inflow rate (i.e., ρ 2 ) until the inflow rate of the 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), 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 in the system with rare abundance due to resource leakage from species N 1 , but would likely go extinct due to stochastic fluctuations, and so coexistence in this limit is somewhat tenuous. These findings expand the expectations of Liebig's Law to two cross-feeding species and demonstrate both species dynamics and equilibrium abundances across various resource landscapes. In contrast, when internal metabolic rates are low, species relative abundances are determined by their respective required resources based on multiplicative co-limitation and both species now demonstrate growth lag-phases. These findings expand the expectations of multiplicative co-limitation to two cross-feeding species. However, since species N 1 will not be resource limited when ρ 1 is greater than ρ 2 , then species N 2 will maintain a higher relative abundance across wider resource landscape (Fig. 4b, see Eq. 9). Together, we find that if internal rates are the same but either high or low compared to the other rates in our model, then our results expand the expectations of Liebig's Law and multiplicative co-limitation to the two-species system with a metabolic dependency. Between these limits, we find a smooth transition between the two scenarios as coexistence conditions become decoupled from classical resource limitation predictions when the internal metabolic rates decrease (Fig. 4c). In addition, we find that coexistence depends on both the internal metabolic rates and the resource inflow rates even when uptake rates are the same. Furthermore, 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 to allow species N 2 to outcompete species N 1 but will result in a lower final density (see Supplemental).
Finally, we consider the dynamics and equilibria when the internal metabolic rates, β i , are different. We find that, when the rates differ the outcome depends on which species has the higher metabolic rate. When the byproduct producer, N 1 , has the higher rate, then the results are similar to when both species have high internal metabolic rates (Fig. 5a). We do note, however, two important differences: 1) species N 2 alone exhibits a growth lag-phase, and 2) both R 1 and R 2 are depleted as the species reach an equilibrium. However, when species N 2 has the higher rate the coexistence conditions and high relative abundances for both species are greatly expanded. In fact, we find coexistence with moderate abundances along all inflow rates tested and the final abundances of both species are determined only by the shared 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 the depleted resources. Between these limits, we find that the transition from classical resource limitation prediction to complete coexistence depends on both metabolic rates decreasing (Fig. 5c). These findings highlight how variation in both the internal metabolic reaction rates and in the environmental conditions can influence species interactions and change coexistence expectations.

Discussion
Here, we reformulated the classic consumer-resource model by independently including resource uptake, internal metabolic rates, and byproduct export. We found that when internal metabolic dynamics are included -in addition to uptake -two common models of resource limitation (Liebig's Law and multiplicative co-limitation) appear in particular limits of the internal reaction rates. In this model, we make predictions for the functional form of the production of metabolic byproducts which can be used in metabolic exchanges, and we balance the requirements for growth, energy, storage, and export. Because our model includes an independent term for internal metabolic rates, it provides a flexible and more general approach to understand how the resource landscape and other physical variables affect species dynamics in microbial communities. We further expand our model to include a second species which uses the metabolic byproduct of the first species. In this metabolically-informed approach, both species interactions (competition and facilitation) and resource use efficiency are emergent properties of the system. In addition, we show how internal metabolic reaction rates and the resource landscape determine species dynamics and equilibria. We find that the growth dynamics change when resources are metabolically limiting which depends on the internal metabolic rates. Therefore our model shows how metabolic rates and the resource landscape change the interactions between cross-feeding species which can shift the contributions of facilitation and competition at equilibrium. We further show that changes in the environmental landscape and/or metabolic rates can change resource and consumer dynamics in ways that alter speciesspecies and species-resource correlations (see resource insets in Figs. 4 & 5); therefore, correlations may not accurately reflect interactions (see: Barner et al. 2018).
Classic formulations for pairwise interactions and consumer-resource dynamics have each led to insights regarding species coexistence, community stability, population self-regulation (Schoener 1983;Barabás et al. 2017;Allesina and Tang 2012;Leibold and McPeek 2006). Lotka-Volterra (and related) equations provide a flexible approach to modeling a range of interactions between species but are unable to generalize across environmental variation because they do not provide an unambiguous way to include the resource landscape. While it may be possible to modify per capita growth rates to be a function of environmental conditions (e.g., temperature) and the resource landscape, how to incorporate these effects in the species interaction rates is ambiguous, at best. However, these models allow positive and negative species interactions to be explored straightforwardly, with far-reaching implications for coexistence and stability. Consumer-resource models explicitly include the interaction between the resource landscape and consumers, but at the expense of introducing more explicit mechanism, and therefore more choices in the way interactions are implemented (O'Dwyer 2018). Including positive interactions through the production of resources has led to new predictions regarding the stability of communities O'Dwyer 2018, 2019), indicating that incorporating resource exchange (i.e., cross-feeding) may be important for understanding the dynamics of natural microbial communities. However, we do not know how sensitive these results are to the precise way consumption and exchange are formulated. Here, we incorporated metabolism more explicitly into consumer-resource dynamics and it allowed us to explain a broader range of community dynamics and natural phenomena, with less ambiguity in the functional form of interactions, and these metabolic rates may also reveal how environmental conditions that can alter metabolic rates (e.g., temperature) will contribute to species dynamics.
Our model used a simplified metabolic network to describe the internal physiology of cells and this allowed us to determine how uptake, growth, and export contribute to the dynamics and coexistence of interacting species. Specifically, we were able to incorporate competition via shared resource consumption and facilitation via the production and consumption of a metabolic byproduct (i.e., cross-feeding). Our model accomplishes this by using a consumer-resource framework with multiple inputs and outputs and by explicitly incorporating internal resource dynamics (see Figs. 1 & 3). Because we incorporate multiple resources as required inputs, there are similarities between our model and the synthesizing unit concept used in the dynamic energy budget theory (Kooijman 1998(Kooijman , 2001. For example, both are able to produce limiting resource dynamics based on Liebig's Law. However, our model was also able to produce multiplicative co-limitation, and it allowed us to explore metabolic cross-feeding by providing a mechanism for the production of cellular byproducts. We accomplish this by allowing internal resources, including luxury uptake and metabolic byproducts, to leak out of cells, which establishes strong conceptual similarities to Droop's model (Droop 1974). Because our model uses this internal flux balance approach, efficiency becomes an emergent property rather than a fixed term. This contrasts with approaches that include fixed efficiency terms that account for the internal processes without providing physiological mechanism (Moore et al. 1993(Moore et al. , 2005Schimel and Weintraub 2003). Likewise, our model differs from others that include explicit resource transfer pathways and corresponding efficiencies between species (Wiegert and Owen 1971;DeAngelis et al. 1975;De Ruiter et al. 1995;Moore et al. 2005) because the byproducts and internal resource excesses are returned to the environment as 'public goods' rather than transferred directly. As such, our model uniquely allows us to incorporate a mixture of species interactions in addition to exploring how each population responds to a complex resource landscape because we include a simplified metabolic network.
Our model allows us to address scenarios that are unique to microbial communities -metabolic cross feeding and internal cellular metabolic rates -and this provided novel insights. First, by providing a mechanism for the production and subsequent consumption of metabolic byproducts (i.e., metabolic cross-feeding) our model provides a way to explore the role of facilitation interactions in a consumerresource framework. Cross-feeding is a widely observed phenomenon in microbial communities (Morris et al. 2012;Morris 2015;D'Souza et al. 2018). However, while crossfeeding appears to be ubiquitous, there are very few hypotheses or theoretical prediction for how it will shape the dynamics of interacting species. Our model provides a novel approach to explore this phenomenon and how it contributes to the diversity, complexity, and dynamics of microbial communities. In our model, the coexistence of species N 2 is dependent on the production of a metabolic byproduct, R 3 , from species N 1 . Assuming there is no external input of the byproduct, the coexistence of N 2 is dependent on N 1 . However, due to the way interactions are formulated, external inputs of byproduct alone are not enough to let N 2 outcompete N 1 (see Supplemental). As such, our model provides a new tool to explore the benefits and consequences of cross-feeding in microbial communities. Second, by separating the uptake and internal metabolic rates, our model is able to produce unique species dynamics and it recapitulates two common models of resource limitation (Liebig's Law and multiplicative co-limitation). One major finding based on our model is that internal metabolic rates can have a strong impact on the dynamics of interacting species. Specifically, changes in the internal metabolic rates can alter species coexistence (see Figs. 4 & 5). There are multiple mechanisms that could modify these metabolic rates, including temperature dependence, cofactors, and redox status (Price and Sowers 2004;Glass and Orphan 2012;Falkowski et al. 2016;Ramírez-Flandes et al. 2019;Russell and Cook 1995). Therefore, our model demonstrates how environmental conditions can radically change the relative abundances of interacting species even when resource inflow rates are fixed. Together, these two components of the metabolically-informed consumer--resource model highlight the potential to address questions and scenarios that may be unique to microbial communities.
Our model provides a mechanistic approach to include diverse species interaction in a consumer-resource framework and generates novel insights, but it also provides a solid foundation to explore new questions. For example, we have added mechanisms to describe both competition and facilitation based species interactions. These interactions are not set parameters as in Lotka-Volterra and therefore may be more realistic in terms of modeling natural systems. However, because we have proposed mechanisms for these interactions, some of our conclusions may not be robust due to the precise way we have implemented the specific interactions. This provides an opportunity to explore how changing the implementation yields additional insights and tests new hypotheses. Likewise, we assumed that uptake and export rates increase linearly with resource concentration, but these parts of our equations could easily be modified to incorporate Monod dynamics (as noted above). These simplifications provide a useful starting point from which other hypotheses can be generated and tested. For example, our simplified metabolic models may not capture important idiosyncrasies regarding how metabolism changes along environmental gradients. Perhaps a multi-species genomescale metabolic model is needed to better describe the metabolic transformations and exchanges between species (Kauffman et al. 2003;Orth et al. 2010;Embree et al. 2015). However, we believe that our simplified version is more tractable, especially when complete genomic information is lacking, and similar approaches have been successfully implemented to describe natural systems (Kempes et al. 2012;Litchman et al. 2015;Chakraborty et al. 2017). The aim of our approach was to find a balance that allows us to more mechanistically describe species interactions while still being tractable for describing natural systems and explaining observed variation.
We propose that when extended more broadly, this approach will lead to more mechanistic predictions for the role of positive interactions along stress gradients (Callaway and Walker 1997;Brooker and Callaghan 1998;Koffel et al. 2018), the possibility of species interactions to stabilize or de-stabilize communities (Butler and O'Dwyer 2018;Allesina and Tang 2012), and the mechanisms underlying biodiversity-ecosystem function relationships (Duffy et al. 2007;Flynn et al. 2011;Cardinale et al. 2012). Species interactions often change across environmental gradients (Brooker and Callaghan 1998;He et al. 2013) and a shift that promotes facilitation may increase ecosystem function (e.g., productivity; see: Cardinale et al. 2002). Likewise, increases in nutrient availability has been shown to favor competition over facilitation (Koffel et al. 2018) and can decrease the importance of mutualisms (Keller and Lau 2018). Our model provides a mechanistic approach to understand these shifts based on the resource landscape and internal metabolic rates. In addition, our model demonstrates that even when competition is favored over facilitation, two species can still coexist. Other studies suggest that positive species interactions, such as facilitation mediated through cross-feeding, can promote coexistence via negative frequency dependence (Morris 2015). Our model demonstrates how even when competition is favored between cross-feeding species, the poorer competitor can remain as a rare member of the community due to the leakiness of the cross-feeder interaction; although, we predict that stochastic fluctuations could lead to exclusion. As such, our model can be used to better understand why some systems favor facilitation versus competition and how facilitation promotes coexistence and potentially the coevolution between species. In short, we propose that metabolicallyinformed consumer-resource dynamics will provide a platform to explore the consequences of cooperative and competitive interactions across environmental contexts.