{'fraction_non_alphanumeric': 0.025082731648616127, 'fraction_numerical': 0.007972322503008424, 'mean_word_length': 5.346778042959428, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling microbial cross-feeding at intermediate scale portrays community dynamics and species coexistence

Social interaction between microbes can be described at many levels of details, ranging from the biochemistry of cell-cell interactions to the ecological dynamics of populations. Choosing the best level to model microbial communities without losing generality remains a challenge. Here we propose to model cross-feeding interactions at an intermediate level between genome-scale metabolic models of individual species and consumer-resource models of ecosystems, which is suitable to empirical data. We applied our method to three published examples of multi-strain Escherichia coli communities with increasing complexity consisting of uni-, bi-, and multi-directional cross-feeding of either substitutable metabolic byproducts or essential nutrients. The intermediate-scale model accurately described empirical data and could quantify exchange rates elusive by other means, such as the byproduct secretions, even for a complex community of 14 amino acid auxotrophs. We used the three models to study each community's limits of robustness to perturbations such as variations in resource supply, antibiotic treatments and invasion by other"cheaters"species. Our analysis provides a foundation to quantify cross-feeding interactions from experimental data, and highlights the importance of metabolic exchanges in the dynamics and stability of microbial communities.


Introduction
Most microorganisms that affect the environments we live in 1 and that impact our health 2 do not live in isolation: they live in complex communities where they interact with other strains and species. The past decade has seen a surge of scientific interest in microbial communities, such as the human microbiome, but most studies remain limited to cataloguing community composition 3 .
Our mechanistic understanding of how biochemical processes occurring inside individual microbial cells command the interactions occurring between cells, and lead to the emergent properties of multi-species communities remains limited 4 .
Microorganisms consume, transform and secrete many kinds of chemicals, including nutrients, metabolic waste products, extracellular enzymes, antibiotics and cell-cell signaling molecules such as quorum sensing autoinducers 5,6 . The chemicals produced by one microbe can impact the behaviors of other microbes by promoting or inhibiting their growth 7 , creating multidirectional feedbacks that drive ecological interactions which may be beneficial or detrimental to the partners involved 8,9 .
If a community is well-characterized and given sufficient data on population dynamics, it should be possible to parameterize the underlying metabolic processes involved in microbemicrobe interactions by fitting mathematical models 10 . Any model can potentially yield insights 11 , but the complexity of most models so far has been either too high for parameterization, or too low to shed light on cellular mechanisms. Microbial processes may be modelled across a range of details: At the low end of this spectrum of details we have population dynamic models such as generalized Lotka-Volterra (gLV) 12 and Consumer-Resource (C-R) models 13 , which treat each organism as a 'black-box' at the cellular level. For example, C-R models assume a linear or Monod dependence of microbial growth on resource uptake kinetics. At the high end of this spectrum, we 5 have detailed single-cell models such as dynamic flux balance analysis (dFBA) 14 and agent-based models 15 that have too many parameters to be parameterizable by experimental data. For example, the linear equations for fluxes obtained from quasi-steady-state assumption of dFBA are highly underdetermined. What is the appropriate level of details to model and constrain microbial processes using data, that may produce not only accurate predictions but also mechanistic insights on microbial communities?
Here we propose a generalizable framework that couples classical ecological models of population and resource dynamics with coarse-grained intra-species metabolic networks. We show that modeling communities at this intermediate scale can accurately quantify metabolic processes from population dynamics data alone. We demonstrate the value of this approach on three engineered communities of Escherichia coli (E. coli) strains with increasing levels of complexity: (1) unilateral acetate-mediated cross-feeding 16 , (2) bilateral amino-acid-mediated cross-feeding between leucine and lysine auxotrophs 17 , and (3) multilateral amino-acid-mediated cross-feeding between 14 distinct amino acid autotrophs 18 . The models report inferred leakage fractions of metabolic byproducts that are difficult to measure directly by experiments, reveal how resource supply and partitioning alter the coexistence and ecological relationships between cross-feeders, and predict the limits of community robustness against external perturbations.

Modeling microbial metabolic processes at an intermediate level is appropriate to fit the
population dynamics data. Inspired by the classical MacArthur's CR models 19 and many followups 13,20-22 , we propose to integrate CR models with a coarse-grained yet mechanistic description of cell metabolism. Metabolic reactions can be broadly classified as catabolic and anabolic, where 6 catabolic reactions break down complex substrates from culture media into smaller metabolic intermediates that can be used to build up biomass components by anabolic reactions. A minimal representation of cell metabolism is a three-layer network composed of growth substrates at the top, metabolic intermediates in the middle, and biomass at the bottom (Fig. 1). Despite its simplicity, this model is flexible enough to describe the transformation of resources into other resources or non-consumable chemicals and biomass, regardless of the specific reactions involved.
Real cells can consume multiple nutritional resources that may be either substitutable or complementary for cell growth. Our model focuses on complementary resources for three reasons: (1) many microorganisms in natural samples are auxotrophs 23 whose growth relies on complementary essential nutrients; (2) minimal medium-popular for cultivating microbial communities in laboratory conditions including the data analyzed in our study-is composed of complementary nutrients; (3) substitutable metabolites can be mathematically lumped into functional groups.
Based on these assumptions, we developed a dynamic modeling framework that contains six kinds of biochemical reactions describing resource consumption, transformation, secretion, utilization for biomass synthesis, and inactivation (Supplementary Equations (S1)-(S6)). Briefly, substrates available in the growth media can be imported into cells. A certain fraction of the imported substrates is then broken down into metabolites, which can either be released back to the surrounding environment or used by cells for biomass production. Released metabolites can be imported by cells in a way similar to externally supplied substrates, except that their uptake may be inhibited by other substitutable substrates that are assumed to be preferentially used. To model the effects of toxic compounds 24 we allow the growth rate of any cell population to be not only governed by a birth-death process that constantly produces and loses cell material due to 7 biosynthetic and maintenance processes respectively, but can be additionally inhibited by accumulation of toxic metabolites in the environment.
The six types of reactions can be translated to differential equations by specifying their kinetic rate expressions. We assumed quasi-steady-state for intracellular substrates and metabolites, as metabolic reactions typically occur at faster time scales compared to ecological dynamics. The time-scale separation thus simplifies our model by excluding intracellular variables, leaving only three types of variables that describe the population density of active cells ( ! , = 1,2, , ⋯ , " ), the extracellular concentrations of substrates ( [ # ], = 1,2, ⋯ , $ ), and the concentrations of metabolic byproducts excreted by cells ([ % ], = 1,2, ⋯ , & ). Assuming a chemostat environment with dilution rate (which reduces to a batch culture when = 0), the differential equations associated with the three state variables are given below where ',# and ',% are the feed medium concentrations of substrate # and metabolite % respectively. !,# )*+,, and !,% )*+,9 represent uptake fluxes of substrates and metabolites respectively, Fitting the model to microbial community data. We applied our framework to published datasets of two two-species communities with increasing level of complexity: a uni-lateral 16 and a bi-lateral 17 cross-feeding between laboratory evolved and engineered strains of E. coli respectively.
Our goal was to manually parameterize the intrinsic metabolic processes relevant for the interactions between the community members, directly from time series data of community composition and experimentally measured metabolite concentrations. The number of metabolites essential for E. coli growth is estimated of the order of hundreds 25 . Therefore, we chose to include in our model as model variables only the metabolites known to mediate interpopulation interactions, together with the most limiting growth substrate. 9 The first community is a well-documented unilateral acetate-mediated cross-feeding polymorphism evolved from a single ancestral lineage of E. coli in laboratory conditions 16 Table 1). The community contains two polymorphic subpopulations (E. coli subspecies) whose metabolism differs in their quantitative ability to uptake and efflux carbon sources: a glucose specialist strain (CV103) which has a faster glucose uptake rate but cannot grow on acetate, and an acetate specialist strain (CV101) which can grow on acetate but has a lower glucose uptake rate. CV103 secretes acetate-a major by-product of its aerobic metabolism-and this way creates a new ecological niche for CV101.  Table 2). The two mutants differ only by single gene deletions in the lysine (DlysA) and leucine (DleuA) biosynthesis pathways. Neither mutant can grow in monoculture, but their coculture can survive by creating a bilateral dependency of two mutants cross-feeding each other missing essential amino acids. Fig. 2G, H show that our model was able to quantitatively recapitulate the growth and nutrient dynamics in both monoculture and coculture conditions. The fitted values of parameters reveal that the maximum growth rate of the lysine auxotroph is over 50% larger than that of the leucine auxotroph ( Fig. 2I), which is consistent with the data showing that the biosynthesis of leucine is more costly than the biosynthesis of lysine 18 . Nonetheless, the parameters also indicate that the mortality rate of the lysine auxotroph (about 20% of its maximum growth rate) is also substantially higher than that of the leucine auxotroph ( Fig. 2J), which qualitatively agrees with cell viability experiments in the monoculture and absence of amino acid supplementation 17 . Since cell mortality rate is determined by the ratio of maintenance rate to nutrient recycling efficiency from dead cells 26 , this finding suggests that the lysine auxotroph has either or both of high maintenance cost and low biomass recovering yield.
Comparison of these two cross-feeding models suggests that resource sharing between natural (CV103 and CV101) and engineered (DlysA and DleuA) cross-feeders can be markedly different. We predicted that the glucose specialist lost 33% carbon in acetate overflow resulting in nearly equal flux values between acetate secretion and glucose uptake, a quantitative relationship that has been observed in a different E. coli strain 27 . By contrast, the engineered interaction between the DlysA and DleuA is much weaker with only 0.3% and 1.4% carbon loss in releasing leucine and lysine respectively. Although the acetate-mediated cross-feeding may have been an incidental finding, the high efflux of acetate could facilitate adaptive co-evolution and accumulation of degenerative mutations 16 .
Metabolic secretion fluxes modulate likelihood of genotypic coexistence. The stable coexistence of different genotypes is a prerequisite for mixed microbial communities. But how strong are the metabolic secretion fluxes necessary to maintain genotypic coexistence in the absence of metabolite supplementation? We leveraged the two cross-feeding models above to address this question by simulating cocultures in chemostats at varied levels of resource supply and partitioning, which independently and synergistically modulate the actual secretion flux values.
We constructed phase diagrams that show how the community composition at steady state has distinct patterns between the two cross-feeding systems (Fig. 3A,B). First, competitive exclusion does not occur when cross-feeding is obligate and bidirectional ( and /,0 are the maximum glucose uptake rates of the glucose and acetate specialists respectively.
When : > ∆ 0 , the glucose specialist releases more acetate than the amount needed to help the acetate specialist overcome its basal growth disadvantage, causing a declining self-balancing capacity of population dynamics and reduced likelihood of coexistence. By contrast, coexistence of the lysine and leucine auxotrophs is only weakly constrained by the resource partitioning level, but ultimately determined by the total amount of resources put into the system (Supplementary Within the region of coexistence, the relative frequency of the acetate specialist increases continuously with the fraction of acetate leaked (Fig. 3A), whereas increasing the fraction of lysine leaked by the leucine auxotroph triggers a discontinuous, abrupt switch from a steady state dominated by the leucine auxotroph to a steady state dominated by the lysine auxotroph (Fig. 3B).
Such abrupt, discontinuous regime shifts are a common feature of microbial communities limited by several essential nutrients 28 . Interestingly, growth of the dominant and rare auxotrophs are always limited by its auxotrophic amino acid and glucose respectively, which suggests an implicit negative feedback loop that maintains their relative abundance ratio before and after the switch: increasing population size of the dominant auxotroph impairs the growth of the rare auxotroph by consuming more glucose but eventually, its own growth is inhibited because a smaller amount of amino acid it needs to grow can be produced by its partner. Taken together, our models show that the likelihood of coexistence can be modulated by varying the metabolic secretion fluxes, but the effect of varying those fluxes depends on the approach used to modulate the system (resource supply or partitioning) and the cross-feeding type (unilateral or bilateral).

Environmental changes to nutrients can reverse the sign of microbial social interactions.
Cross-feeding interactions within a microbial community may be described as social interactions with costs and benefits to the members involved 29,30 . Those costs and benefits may be altered by environmental perturbations that supply or remove the cross-fed metabolites form the environment. The ecological relationship between the glucose and acetate specialists was displayed on a 2-dimensional phase space spanned by the feed medium concentrations of glucose and acetate (Fig.   4A). The entire space is divided into six distinct regions with diverse outcomes, including population collapse, competitive exclusion, and stable coexistence. Notably, it is very difficult to select supplementation resulting in stable coexistence. This is because, as explained above, the inferred value of : (0.33) is much greater than that of ∆ 0 (0.12). The remaining diversity of the phase space structure is primarily driven by the dose-dependent effect of acetate 24 : it serves as a 13 nutrient for the acetate specialist at low concentration but becomes inhibitory to growth of both strains when abundant ( Supplementary Fig. 1). To illustrate this effect, we increased glucose supplementation from P1 to P3 (gray dots in Fig. 4A) in the phase space, which induced higher release of acetate to environment (Fig. 4B, top row) and switch of winners of the coculture competition (Fig. 4B,  therefore, the glucose specialist wins again when the negative inhibitory effect of acetate outweighs its positive nutritional effect on the acetate specialist. Compared to unilateral cross-feeding, new ecological relationships such as mutualism and parasitism emerges in the phase space when cross-feeding is bidirectional (Fig. 4C). The mutualistic relationship was maintained over a broad range of supplied amino acid concentrations, even though amino acid supplementation releases the dependence of one auxotroph on the other and is hence detrimental to mutualism. In the regime of mutualism, glucose is in excess and both strains are limited by the essential amino acids they cannot produce (Fig. 4D, left column). Further addition of amino acids leads to strain dominance, but not necessarily competitive exclusion. The lysine auxtroph was excluded when leucine was provided to release the leucine auxotroph from its growth dependence (Fig. 4D, middle column), whereas adding lysine only reduced the relative abundance of the leucine auxotroph, rather than leading to the loss of its entire population (Fig.   4D, right column).
14 Amino acid supplementation may lead to competitive exclusion or parasitism depending on whether one or both auxotrophs are limited by glucose. When glucose limits both auxotrophs, the leucine auxotroph wins because it has the same growth rate as the lysine auxotroph on glucose but lower death rate (Fig. 2I,J). When only the lysine auxotroph is limited by glucose, the leucine auxotroph can sustain its population by occupying a different niche and growing on leucine released by its competitor. Regardless of the outcome, our results suggest that adding cross-fed nutrients can induce competition between community members that previously interacted mutualistically, and shift positive interactions to negative interactions.

Uncovering complex cross-feeding interactions between 14 amino acid auxotrophs.
Next we demonstrated the utility of our model to study cross-feeding interactions within communities of more than two members. We modeled a community of 14 amino acid auxotrophs engineered from been accumulated to a critical size that allows the actual growth to compensate for mortality and dilution. As the pool size increases, its net growth rate (growth minus mortality) surpasses the dilution rate and recovers its population density, which eventually levels off when the positive and negative forces reach equilibrium. By fitting the population density dynamics, we concomitantly inferred the concentration dynamics of glucose and all amino acids ( Supplementary Fig. 3), which are hidden states (not yet observed) that are relatively costly and inaccurate to measure in experiments.
Cross-feeding network is prone to collapse upon external perturbations. By simulating the 14-auxotroph community model to steady state, we predicted that the initial mixture converges to a stable coexisting subset that contains 4 auxotrophs that are deficient in biosynthesis of isoleucine (DI), lysine (DK), methionine (DM), and threonine (DT) (Fig. 6A). The predicted coexistence state was successfully validated by two independent observations over 50-day serial dilution 18 , a much longer period of time than the duration of the training dataset (7-day serial dilution; Fig. 5D). The We computationally tested how external perturbations, including nutrient downshift, the addition of antibiotics, and invasion of cheating phenotypes (the same auxotrophic dependence but no amino acid leakage) affect the stability of coexistence among the 4 auxotrophs (see Methods).
The 4-strain community was able to cope with these disturbances to a certain extent and remained integrated. Beyond the thresholds, all three perturbation types resulted in community collapse as a result of domino effect (Fig. 6C-E), implying that tightly coupled cooperative communities are fragile and prone to collapse. Since antibiotics inhibit growth of individual strains (targeting consumer nodes in the bipartite network) while cheaters are amino acid sinks (targeting resource nodes in the bipartite network), we identified that DT and methionine as the weakest consumer node (Fig. 6D) and resource node (Fig. 6E) in the bipartite network respectively. Our results suggest that DTàK (secretion of lysine by the threonine auxotroph) and MàDM (uptake of methionine by the methionine auxotroph)-the outgoing links from the two weakest nodes that are also essential to maintain community integrity-are the weakest metabolic fluxes that may set the resistance level of the community to external perturbations 32 . Our model could extract new insights from previously published empirical data. The analysis shows that unidirectional cross-feeding is equivalent to a commensalism and bidirectional cross-feeding is equivalent to a mutualism. As shown by our study (Fig. 4) and previous work 24,29 , the actual relationship between cross-feeders, however, can be diverse in even simple and constant environments (e.g., glucose minimal medium) due to a combination of positive effects of crossfeeding with negative effects of competition and toxicity of cross-fed metabolites, suggesting that the exact outcome cannot be precisely delineated by the cross-feeding type alone. Moreover, mechanistic models can help identify knowledge gaps 43 . For example, recent experiments have demonstrated that the coexistence of two carbon source specialists in the unilateral cross-feeding example is mutualistic in the sense that the consortium is fitter than the individuals 44 . The syntropy 19 can be explained by a null expectation from theoretical ecology models 45 : the glucose specialist provides acetate in an exchange for a service provided by the acetate specialist which scavenges the acetate down to a level at which growth inhibition is insignificant. Although we thoroughly considered the mechanism of resource-service exchange, additional features of our model and/or the use of data-consistent parameter values did not support mutualistic coexistence in any environmental condition we tested (however, competitive coexistence is possible). The discrepancy suggests that our model and even the classical resource-service exchange theory have missed some qualitative or quantitative details that are the key to understanding of syntrophic mechanisms in this specific example.

Discussion
What could we have missed? Since mutualism occurs when the reciprocal benefits associated with cross-feeding outweigh competitive costs 46 , our model should logically predict either or both of lower benefits and higher costs than the null expectation from simpler models. In the classical theory of syntropy, it is typically assumed that leaking chemicals are by-products which are inhibitory to producers but beneficial to consumers 45 . Since acetate was shown to inhibit growth of both cell types (Fig. 2D) and acetate specialist (the consumer) is more sensitive, its population density may be insufficient to reward the glucose specialist to a level that allows benefits higher than costs. On the other hand, costs are potentially similarly high since both cell types are polymorphic and share similar glucose uptake kinetics. We estimated that the relative difference in their maximum growth rates is 12%, which is much smaller than the observed value in experiments (33%) 16 . This quantitative difference may be important considering that the competition is stronger between populations with similar nutrient acquisition strategies. Recently, it was theoretically proposed that controlled metabolic leakages optimize resource allocation and can be beneficial to producers even under nutrient limitation 47 . We speculate that in case where 20 acetate overflow improves, rather than negatively impacts, the growth of producers, the likelihood of forming a mutualistic pair between two cell types would be much higher. Overall, the costbenefit nature of the cross-feeding interaction between polymorphic E. coli strains is more complex than thought before and warrants further research.
So far, the current framework has been applied to well-characterized communities with known chemicals and associated interactions. Can it be applied to infer community structure of complex microbiomes (e.g., human gut microbiome) where most of the metabolic exchanges involved in microbe-microbe interactions are still unknown? Our model has the potential if some technical challenges can be solved. First, direct modeling of a real-world microbiome with hundreds of species would be hurdled by too many unknown kinetic parameters. One way to solve this problem is to simply ignore the rare species 35 . Another-arguably better-approach might be by grouping species composition into functional guilds using unsupervised methods that infer those groups from the data alone 48 , or to use prior knowledge from genomics or taxonomy to create such functional groups. Second, inferring chemical mediators within a community of interacting populations is a nontrivial task. It can be facilitated by prior knowledge such as searching the literature or leveraging systems biology tools such as community-level metabolic network reconstruction 49 . Finally, our model is nonlinear, so that an efficient and robust nonlinear regression approach for parameter estimation is essential. Manual parameter selection is often the only possible approach for small datasets like the experimental systems we analyzed here. Indeed, non-linear optimization algorithms often fail to converge to a realistic set of parameters. Although we chose the manual method to calibrate our models in this proof-of-concept study, manual fitting requires an expert operator and is a time-consuming process, which for now precludes it from being applied to large-scale microbial communities. On the positive side, the process of trial-and- They were mixed with the resident community in varying ratios at the beginning of simulation.
For all three perturbation types, the feed medium glucose concentration is 0.2 wt% in the unperturbed condition and serial dilution was run to steady state at 60 days.