# Unveiling abundance-dependent metabolic phenotypes of microbial communities

Natalia E Jiménez, Vicente Acuña, María Paz Cortés, Damien Eveillard, Alejandro Eduardo Maass

► To cite this version:

Natalia E Jiménez, Vicente Acuña, María Paz Cortés, Damien Eveillard, Alejandro Eduardo Maass. Unveiling abundance-dependent metabolic phenotypes of microbial communities. mSystems, In press, 10.1128/msystems.00492-23. hal-04198586

HAL Id: hal-04198586

https://hal.science/hal-04198586v1

Submitted on 7 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Unveiling abundance-dependent metabolic phenotypes of microbial communities

Natalia E. Jiménez,¹,² Vicente Acuña,¹,² María Paz Cortés,¹ Damien Eveillard,³ Alejandro Eduardo Maass¹,²,⁴

**AUTHOR AFFILIATIONS** See affiliation list on p. 21.

**ABSTRACT** Constraint-based modeling has risen as an alternative for characterizing metabolism of communities. Adaptations of flux balance analysis have been proposed to model metabolic interactions, which in most cases consider the maximization of biomass production as their objective. In nature, novel essential functions are not directly related to cell growth force communities to display suboptimal growth rates. These suboptimal states allow a degree of plasticity in their metabolism, thus allowing quick shifts between alternative flux distributions as an initial response to environmental changes. In this work, we introduce the abundance-growth space as a representation of metabolic phenotypes of a community. This space is defined by the composition of a community, represented by its members' relative abundances, and their growth rate. The analysis of this space allows us to pinpoint how critical reactions respond to shifts of the environment, showing where changes in community plasticity occur. Interestingly, it highlights the relevance of the relative abundance of its members in the lost or gain of plasticity. This method is applied to two simple communities that exchange metabolites. A synthetic community of two mutant *Escherichia coli* strains and an environmental bioleaching community composed by *Acidithiobacillus ferrooxidans* Wenelen and *Sulfobacillus thermosulfidooxidans* Cutipay, where only Cutipay consumes organic matter disposed of by the community.

**IMPORTANCE** In nature, organisms live in communities and not as isolated species, and their interactions provide a source of resilience to environmental disturbances. Despite their importance in ecology, human health, and industry, understanding how organisms interact in different environments remains an open question.

In this work, we provide a novel approach that, only using genomic information, studies the metabolic phenotype exhibited by communities, where the exploration of suboptimal growth flux distributions and the composition of a community allows to unveil its capacity to respond to environmental changes, shedding light of the degrees of metabolic plasticity inherent to the community.

**KEYWORDS** constrained-based modeling, microbial communities, metabolic plasticity

In nature, organisms live in communities rather than as isolated species. These communities emerge from interactions (1) and provide a source of resilience to perturbations such as biofilm development (2, 3), flux relation-feeding (4), or community resistance to environmental disturbances (5). The study of these interactions is imperative in critical areas such as ecology (6), health (7), and the biotechnological industry (8). Understanding how organisms interact in different environments is still an open discussion despite their relevance.

The last decade saw the rise of metabolic modeling for formalizing microbial interactions. The metabolic cross-talking between microorganisms notably justifies the

**Editor** Karoline Faust, Katholieke Universiteit Leuven, Leuven, Belgium

Address correspondence to Alejandro Eduardo Maass, amaass@dim.uchile.cl.

Natalia E. Jiménez and Vicente Acuña contributed equally to this article. Author order was determined by who performed the implementation of the proposed method.

The authors declare no conflict of interest.

See the funding table on p. 21.

**Received** 24 May 2023

**Accepted** 21 June 2023

**Published** 5 September 2023

Copyright © 2023 Jiménez et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.
metabolic abstraction for sustaining biogeochemical cycles (9). As a natural following, several approaches have been developed for modeling the metabolism of communities (10) using comprehensive genome-scale descriptions for which each organism is described by its inner biochemical machinery (11–14). This genome-scale description allows developing of several computational strategies to identify essential metabolic interactions (e.g., SteadyCom [15], RedCom [13], community flux balance analysis [cFBA] [16], and MICOM [microbial community] [17]). When applied in specific contexts, it results in insightful advances like the determination of essential interactions in anaerobic digestion for biogas production (18), cancer (19), and the gut microbiota (20, 21).

These approaches characterize metabolic flux distributions and community compositions using an adaptation of flux balance analysis (FBA) (22), where an objective function, usually biomass production, is optimized (13, 15–17, 23, 24). Maximization of biomass is of great interest for growth rate estimation, and it has shown substantial benefits in the biotechnological context that aims at controlling single strains to improve yields for a product of interest (25, 26). However, in environmental communities, organisms should be prepared to face shifting conditions. Hence they increase their fitness by maximizing their growth rate for a subset of possible conditions (27), hence displaying suboptimal growth rates. Moreover, metabolic models do not always account for additional requirements, such as the synthesis of secondary metabolites essential to sustain communities (28) and survival in adverse environments (5). Suboptimal growth solutions allow for resources to be used in these additional functions, as it was observed experimentally for model species in references 29 and 30.

A challenge for studying and predicting metabolic interactions in community models is to propose methods allowing the exploration of alternative metabolic states that a community could display while preserving the practicality of traditional methods like FBA and flux variability analysis (FVA) (31). An idea that can be adapted to community models was proposed initially in reference 32, where metabolic phenotypes of an organism are studied by grouping couples of nutrient fluxes by their incidence in the optimal solution of an FBA for biomass. More recently, an alternative approach is proposed to directly explore the flux space defined by exchange reactions to characterize the so-called metabolic niche (33). Comparison between this space for different organisms allows us to determine environments where a community can thrive. However, these approaches still leave out the role and abundance of each organism in the community. Indeed, the relative abundance of members of a community shape the distribution of resources and synthesis of secondary metabolites (15). Although there are promising ideas regarding the laws governing its probabilistic distribution (34), the prediction of abundances is uncertain at steady state. Thus, the exploration of the effect of community composition is critical for studying metabolic interactions.

In this work, we explore the space of metabolic fluxes of a community, focusing on the distribution of abundances of its organisms and community growth rates. This space is partitioned according to displayed metabolic phenotypes compatible with each abundance-growth point, computed based on flux ranges for each reaction. These ranges are seen as descriptors of the system's plasticity (35, 36) and allow us to distinguish between reactions that present flux plasticity, where their fluxes remain positive (i.e., always active) despite flux variations, from reactions associated with structural plasticity that shows zero flux on some solutions, emphasizing their putative replacement by alternative pathways. Moreover, the proposed framework allows projecting the metabolic phenotype in a lower dimensional space for interpretation, thus providing a practical exploration of metabolic phenotypes in a community characterized by different relative abundances and growing at suboptimal growth rates.

For illustration, the method is applied to a synthetic *E. coli* community with a mandatory cross-feeding of leucine and lysine, and an environmental community composed of *Acidithiobacillus ferrooxidans* Wenelen and *Sulfobacillus thermosulfidooxidans* Cutipay, where only the latter consumes organic matter disposed of by the community.
In both examples, the study of the abundance-growth space allows us to pinpoint how critical reactions respond to shifts in the environment, showing where changes in community plasticity occur. Notably, the method not only allows us to confirm that plasticity increases as the growth rate decreases as expected (29) but also highlights the relevance of the relative abundance of its members in the loss or gain of plasticity. In a more quantitative perspective, in both examples we show that despite the different degrees of plasticity, there are strong couplings between key reactions when reaching high growth rates or attaining very unbalanced communities.

# RESULTS

## Rationale of the abundance-growth space

We aim to understand mechanisms that support community structure and the ecosystem's resilience to perturbations. To achieve this, we propose to study the metabolic phenotypes of a community through the definition of what we call the abundance-growth space.

A community model is constructed using single-organism genome-scale metabolic models, considering each organism as a compartment. Exchanges between organisms occur in the pool compartment. We assume the community is stable over time, i.e., all organisms grow at the same rate to maintain constant relative abundances, which is common in community modeling (13, 15, 16). The requirement of an identical growth rate between all organisms applies not at every time point but averaged over time. Thus, this hypothesis is a reasonable approximation for communities where composition is mainly maintained between spaced time points (15).

In this model, $\mu$ is the community growth rate, and the relative abundance $f_i$ of organism $i$ appears as a factor of the flux bounds for any reaction of its respective organism (Materials and Methods). Thus, variations on the vector of abundances $F = (f_1, ..., f_N)$ and on $\mu$ affect the set of flux distributions that are feasible solutions of the model. Specifically, for a given abundance-growth pair $(F, \mu)$, the *polytope* $P_{F,\mu}$ is defined as the set of all flux distributions satisfying model's constraints. We define the abundance-growth space as the set of all feasible pairs $(F, \mu)$, that is, when $P_{F,\mu} \neq \emptyset$. We believe that this space is a good asset to capture the metabolic phenotype when changing growth rates and abundances of members of a community.

Characterizing $P_{F,\mu}$ for each point, $(F, \mu)$ of the abundance-growth space is computationally hard and complex to describe. Thus, we represent it by the minimum $[\min_{F,\mu}(r)]$ and maximum $[\max_{F,\mu}(r)]$ values of the flux $v_r$ of reaction $r$ in $P_{F,\mu}$. We propose the range between $\min_{F,\mu}(r)$ and $\max_{F,\mu}(r)$ as an indicator of metabolic plasticity as described in references 35 and 36. Indeed, a reaction $r$ has flux plasticity if $\min_{F,\mu}(r) \neq \max_{F,\mu}(r)$. In addition, the reaction has structural plasticity if zero is in between these extreme values, meaning that alternative pathways could replace $r$. To obtain a more quantitative description of $P_{F,\mu}$, one can consider the set of combined feasible fluxes of pairs of selected reactions (quantitative coupling). Computationally, this representation is achieved by defining a *grid* $\mathcal{F}$ of feasible abundance-growth points (Materials and Methods). For an example of the abundance-growth space for a community of two organisms, see Fig. 1B [we only plot $(f_1, \mu)$ since $f_2 = 1 - f_1$].

Following this rationale, for studying metabolic phenotypes given by a set of reactions of interest $R'$ on an abundance-growth point $(F, \mu)$, we assign a categorical value to $r \in R'$ according to Table 1. In an exploratory study, $R'$ can be the whole set of reactions, but in a more precise context, for instance when focusing on the community's metabolic interactions, it can be restricted to exchange reactions. Then, for each point of the grid, we compute a categorical vector for reactions in $R'$, generating a partition of the abundance-growth space in zones with identical vectors. Depending on the size of $R'$, the number of zones can be huge. To overcome this, we define a cluster-partition, whose zones are obtained by a clustering algorithm, and we attribute a consensus categorical vector to each cluster-partition (Materials and Methods).
<table><caption>TABLE 1 Categorical values for a reaction *r* and a point (*F*, *μ*) of the abundance-growth space</caption><thead><tr><td>Condition</td><td>Category</td><td>Flux sign</td><td>Flux variability</td><td>Flux requirement</td><td>Level of plasticity</td></tr></thead><tbody><tr><td>min<sub>F,μ</sub>(r) = max<sub>F,μ</sub>(r) = 0</td><td>0</td><td>Zero</td><td>Fixed</td><td>Off</td><td>No plasticity</td></tr><tr><td>min<sub>F,μ</sub>(r) = max<sub>F,μ</sub>(r) &lt; 0</td><td>–</td><td>Negative</td><td>Fixed</td><td>Mandatory</td><td>No plasticity</td></tr><tr><td>0 &lt; min<sub>F,μ</sub>(r) = max<sub>F,μ</sub>(r)</td><td>+</td><td>Positive</td><td>Fixed</td><td>Mandatory</td><td>No plasticity</td></tr><tr><td>min<sub>F,μ</sub>(r) &lt; max<sub>F,μ</sub>(r) &lt; 0</td><td>––</td><td>Negative</td><td>Variable</td><td>Mandatory</td><td>Flux plasticity</td></tr><tr><td>0 &lt; min<sub>F,μ</sub>(r) &lt; max<sub>F,μ</sub>(r)</td><td>++</td><td>Positive</td><td>Variable</td><td>Mandatory</td><td>Flux plasticity</td></tr><tr><td>min<sub>F,μ</sub>(r) &lt; max<sub>F,μ</sub>(r) = 0</td><td>–0</td><td>Negative or zero</td><td>Variable</td><td>Optional</td><td>Structural plasticity</td></tr><tr><td>0 = min<sub>F,μ</sub>(r) &lt; max<sub>F,μ</sub>(r)</td><td>0+</td><td>Positive or zero</td><td>Variable</td><td>Optional</td><td>Structural plasticity</td></tr><tr><td>min<sub>F,μ</sub>(r) &lt; 0 &lt; max<sub>F,μ</sub>(r)</td><td>–+</td><td>Any</td><td>Variable</td><td>Optional</td><td>Structural plasticity</td></tr></tbody></table>

## The abundance-growth space for a synthetic *E. coli* community

To illustrate this method, we first consider a synthetic cooperative community composed of two *E. coli* strains auxotrophic for leucine (eco_L) and lysine (eco_K) (File S1, https://github.com/mathomics/ecosystem), studied in reference 37. Consumption and production of these amino acids were constrained based on their reported fluxes (Materials and Methods). In this community, cross-feeding of leucine and lysine is critical for survival (Fig. 1).

We obtain an abundance-growth space whose shape resembles a concave triangle (Fig. 1B). Its dimensions are influenced by changes in the production capacity of leucine and lysine, as well as its external supplementation (Fig. S1), as discussed below. This shape shows that at high growth rates, abundances need to be balanced ($f_{eco_K} \sim f_{eco_L}$), which is consistent with reference 37. On the other hand, communities with a dominant strain can only grow at lower rates, which was attributed to the less abundant strain being unable to satisfy amino acid requirements from the other (37). Indeed, at a fixed abundance of eco_K (eco_L), the minimum flux value for the production of leucine (lysine) reaches the upper bound in the model, making it unfeasible for greater growth of the community (Fig. S2).

From the 5,068 reactions of the model, a cluster-partition is performed. Eight clusters allow us to distinguish relevant changes in this space (Materials and Methods and Fig. S3a). We distinguish 116 reactions whose consensus values change between clusters, 1,880 that are unable to carry any flux in the abundance-growth space, and 2,371 presenting fixed non-zero flux values (File S2, https://github.com/mathomics/ecosystem). Blocked reactions are associated with the exchange of metabolites absent in the environment of the simulated community, as well as consequences of gene knockout performed on both strains. Fixed reactions are associated with biomass production requirements, such as NAD and tetrahydrofolate synthesis and transport of ions present in biomass. Interestingly, only leucine consumption by eco_L exhibits a fixed flux value in all the abundance-growth space (status –), while eco_K displays a small range in which consumption of lysine is observed (status –).

From the 116 changing reactions, 98 appear only in shifts between clusters of unbalanced communities at very low growth rates (clusters 2, 5, 1, and 4). Among the remaining 18 reactions, we determine that 10 are able to describe all the transitions between clusters, since redundant reactions (such as transporters of O₂ and H₂O) can be removed. Changes of these 10 reactions show that the community is set toward cell growth by breaking down glucose, rather than amino acid synthesis and consumption. Indeed, these shifts correspond to indirect measurements of how active glucose catabolism is, hence displaying exchanges for oxygen, carbon dioxide, and water, for the community. Between clusters 7 and 8, shifts are observed in reactions depicting an active metabolism toward biomass synthesis: ATP synthesis (from $–$ to $+$ [Fig. 1C]) and ribose-5-phosphate isomerase (RPI), which reverses its flux toward nucleotide synthesis.

In addition, to get a deterministic description of the abundance-growth space, we compute the partition given by the qualitative vectors associated with these 10
<figure><img src="image_2.png" /><figcaption>FIG 1 Analysis of the abundance-growth space for a synthetic *Escherichia coli* community. (A) A synthetic cooperative microbial community of two *E. coli* mutants is constructed by knocking out key genes for leucine (*eco_L*) and lysine (*eco_K*) production. Bounds on exchange reactions for each strain are constrained based on their abundances (depicted as intervals). (B) Cluster-partition ($k = 8$) of the abundance-growth space for the presented community with no leucine or lysine supplementation. (C) Table of reactions in which qualitative states change between neighboring clusters. Clusters are organized according to their location in the abundance-growth space from low to high growth rate and from $f_{eco_K}$ 0 to 1. Lysine and leucine exchange reactions (gray rows) do not exhibit changes among the analyzed clusters. EX: exchange reactions, h2o: water, o2: oxygen, lys: lysine, leu: leucine, ATPS4rpp: ATP synthase, RPI: ribose-5-phosphate isomerase.</figcaption></figure>

reactions. This yields nine areas (Fig. S4) revealing in more detail the shifts observed in the cluster partition. Particularly, shifts forcing ATP synthesis and water exchange by eco_L (Fig. 1A).

Finally, since states for both leucine and lysine exchanges remain constant along the abundance-growth space, we study the degree of quantitative coupling between these reactions. For this, we plot flux values for different couples of these reactions in six points of the space (Fig. 2A). In particular, Fig. 2B shows the consumption of lysine and production of leucine by eco_K. Here, flux ranges are affected by both composition and growth in such a way that when eco_K is more abundant and requires higher growth rates, lysine consumption by this bacterium tends to be fixed, but the production of leucine can vary on a wide interval. On the contrary, leucine production by eco_K tends to be fixed when this bacterium is less abundant at higher growth rates.

Also, in Fig. 2C, consumption of lysine by eco_K vs consumption of leucine by eco_L are plotted, showing the following assymetry: eco_L requires a fixed flux of leucine that increases at higher growth; on the other hand, the flux of lysine consumption by eco_K is more flexible, living in an interval which decreases with its abundance. Fixed values on leucine consumption reflect its exclusive requirement for biomass synthesis, while lysine
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**FIG 2** Quantitative flux coupling analysis for the exchange of lysine and leucine in the synthetic *E. coli* community. (A) Selected points in the abundance-growth space. (B) Feasible fluxes for the consumption of lysine and production of leucine by eco_K at the different points marked in panel A. (C) Feasible fluxes for the consumption of lysine by eco_K and consumption of leucine by eco_L. Fluxes are in units of (mmol/g DWorg/h).</figcaption></figure>

can be transformed into 1,2-diaminepentane through the lysine decarboxylase reaction and this metabolite can be exported.

Although this has not been reported (37), by simulating the supplementation of leucine and lysine, this method can give insights into how metabolism is affected by shifts in its environment. When both lysine and leucine are supplemented (Fig. 3) at a flux of 0.01 (mmol/g DWcom/h), we get that a cluster-partition of 10 sets accurately represent qualitative states (Fig. S3B).

We observe an area of low growth rates, where the community can opt to not synthesize leucine and lysine (clusters 7, 9, 8, and 10). By increasing the community growth rate, both bacteria are required to produce leucine and lysine. First, in clusters 4 and 6, the key reaction for leucine synthesis (eco_K_IPPS) changes from 0+ to ++. But its supplementation for the other bacterium (eco_K_EX_leu_L_e) is only required when eco_K is less abundant (cluster 6) (Fig. 3). In Fig. S5, a partition of the abundance-growth space for the previous reactions highlights these phenomena. For lysine synthesis, an analogous phenomena is observed, where its key reaction (eco_L_DAPDC) shifts from 0+ to ++. From this point upwards (clusters 5, 1, 2, and 3), both bacteria are required to
A

<figure><img src="image_6.png" /></figure>

B

<table><caption>FIG 3 Analysis of the abundance-growth space for a synthetic community of *E. coli* supplemented with leucine and lysine. (A) Cluster-partition ($k = 10$) of the abundance-growth space. (B) Key reactions associated with changes in their states between neighboring clusters together with their qualitative values along the clusters. Clusters are organized according to their location in the abundance-growth space from low to high growth rates and from $f_{eco_K}^0$ to 1. EX: exchange reactions, h2o: water, o2: oxygen, lys: lysine, leu: leucine, ATPS4rpp: ATP synthase, IPPS: 2-isopropylmalate synthase, RPI: ribose-5-phosphate isomerase, DAPDC: diaminopimelate decarboxylase</caption><thead><tr><th></th><th>c7</th><th>c9</th><th>c8</th><th>c10</th><th>c6</th><th>c4</th><th>c5</th><th>c1</th><th>c2</th><th>c3</th></tr></thead><tbody><tr><td>EX_h2o_e</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>++</td><td>++</td></tr><tr><td>EX_o2_e</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>--</td><td>--</td></tr><tr><td>eco_K_ATPS4rpp</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>++</td><td>++</td></tr><tr><td>eco_K_EX_h2o_e</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>++</td><td>++</td></tr><tr><td>eco_K_EX_leu_L_e</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>++</td><td>+</td><td>++</td><td>++</td><td>++</td></tr><tr><td>eco_K_IPPS</td><td>0+</td><td>0+</td><td>0+</td><td>0+</td><td>++</td><td>++</td><td>++</td><td>++</td><td>++</td><td>++</td></tr><tr><td>eco_K_RPI</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>--</td></tr><tr><td>eco_L_ATPS4rpp</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>++</td><td>++</td></tr><tr><td>eco_L_DAPDC</td><td>0+</td><td>0+</td><td>0+</td><td>0+</td><td>++</td><td>+</td><td>+</td><td>++</td><td>++</td><td>++</td></tr><tr><td>eco_L_EX_h2o_e</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>++</td><td>++</td></tr><tr><td>eco_L_EX_lys_L_e</td><td>+</td><td>+</td><td>+</td><td>-</td><td>-</td><td>+</td><td>+</td><td>+</td><td>++</td><td>++</td></tr><tr><td>eco_L_EX_o2_e</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>--</td><td>--</td></tr><tr><td>eco_L_RPI</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>--</td></tr></tbody></table>

supplement lysine and leucine, respectively (++ state). No more changes are observed with respect to the scenario without the supplementation of lysine and leucine.

The quantitative analysis of Fig. 2 in this scenario only reflects that there is more plasticity due to amino acid supplementation (Fig. S6), measured as an increase in the range of minimum and maximum fluxes.

## The abundance-growth space for a bioleaching community

To test the method in an environmental community, we consider a simple microbial consortium between *A. ferrooxidans* Wenelen and *S. thermosulfidooxidans* Cutipay (see File S3 and File S4 in https://github.com/mathomics/ecosystem and Materials and Methods), obtaining a community competing for inorganic energy sources (Fe(II) and thiosulfate). Both bacteria are able to fix carbon available as CO₂ but only Cutipay can consume organic matter made available by the degradation of community biomass, which is modeled as a pseudo-reaction transforming a fraction $\alpha$ of the biomass produced by the community (Materials and Methods).

We consider the case study where $\alpha = 0.2$ and a substrate level composed of energetically equivalent amounts of Fe(II) and thiosulfate. The obtained abundance-growth space has a convex shape such that the maximum community growth rate increases when Wenelen decreases, showing that Cutipay is more efficient in producing biomass than Wenelen. Since both bacteria have been observed together in equilibrium (38), it seems that they are not maximizing the growth rate of the community in these conditions.

Considering $R'$ as the set of 100 exchange reactions, a cluster-partition is computed using 20 clusters (Fig. 4B) (Materials and Methods and Fig. S7 to justify the number of clusters).

We find that 23 reactions change between neighboring clusters (Fig. 4B). The 27 remaining reactions are blocked, coupled with their biomass production requirements, or exhibited structural plasticity on each grid point. The curation of these 23 reactions allowed to remove those that are coupled in all clusters, such as consumption of Fe(II) and production of Fe(III), pairs of exchange reactions for H₂O and H, and SO₄ and H. Additionally, reactions in which changes are only attributed to a few clusters were
A

<figure><img src="image_7.png" /></figure>

B

<figure><img src="image_8.png" /></figure>

C

<table><caption>FIG 4 Analysis of the abundance-growth space for the bioleaching community. (A) A community composed of *A. ferrooxidans* Wenelen (wen) and *S. thermosulfidooxidans* Cutipay (cut) is studied in an environment with Fe(II) and thiosulfate in the presence of organic matter. A pseudo-reaction represents organic matter (OM) availability as a fraction ($α$, here $α = 0.2$) of the total community biomass. (B) Twenty clusters are computed for the 100 exchange reactions of the community model. (C) Key reactions associated with qualitative changes between neighboring clusters define six areas for the behavior of this community. A status is associated with a cluster when over 80% of the points in the cluster exhibit a certain state. If no consensus is reached, an empty cell is presented. Clusters are organized according to their location in the abundance-growth space from low to high growth rate and from $f_{wen}$ 0 to 1. OM: organic matter, EX: exchange reactions, h: hydrogen (proton), tsul: thiosulfate, co2: carbon dioxide, fe2: Fe(II), glc: organic matter, CW-P: Cutipay Wenelen plasticity.</caption><thead><tr><th rowspan="2"></th><th colspan="3">High plasticity</th><th colspan="3">Proton plasticity</th><th colspan="3">Wenelen plasticity</th><th colspan="2">CW-P</th><th colspan="3">Cutipay plasticity</th><th colspan="3">Low plasticity</th></tr><tr><th>c3</th><th>c2</th><th>c1</th><th>c4</th><th>c11</th><th>c12</th><th>c10</th><th>c20</th><th>c19</th><th>c16</th><th>c14</th><th>c13</th><th>c9</th><th>c8</th><th>c7</th><th>c15</th><th>c17</th><th>c18</th><th>c5</th><th>c6</th></tr></thead><tbody><tr><td>EX_f2_e</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>EX_tsul_e</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>EX_h_e</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>cut_EX_f2_e</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>-0</td></tr><tr><td>cut_EX_tsul_e</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>--</td><td>--</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>-0</td><td>-0</td></tr><tr><td>cut_EX_h_e</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>--</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>--</td><td>--</td><td>+</td></tr><tr><td>cut_EX_glc_e</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>cut_EX_co2_e</td><td>--</td><td>+</td><td>+</td><td>--</td><td>--</td><td>--</td><td>+</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>wen_EX_f2_e</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>-0</td><td>-0</td><td>--</td><td>--</td></tr><tr><td>wen_EX_tsul_e</td><td>-0</td><td>-0</td><td>--</td><td>-0</td><td>-0</td><td>--</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>-0</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>wen_EX_h_e</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr></tbody></table>

removed, notably the exchange of H₂, which production is optional only in cluster 1. After this curation, 11 critical reactions associated with iron, sulfur, and carbon metabolism were obtained. Five of them belong to Cutipay, three to Wenelen, and three to the microbial community as a whole (Fig. 4C).

Six sections are distinguished (Fig. 4B and C). A high plasticity section is defined at lower growth rates, where most exchanges are optional (–0 and –+) and independent of abundances. As growth requirements increase, oxidation of both inorganic energy sources by the community becomes mandatory (–), as well as the consumption of
sulfur and iron, which indicates a loosening of global structural plasticity. Next, we define a proton plasticity section (clusters 10, 11, and 12), where the community needs to consume Fe(II) and thiosulfate, preserving structural plasticity (−+) for the exchange of protons (H). This indirectly measures the balance between Fe(II) and thiosulfate oxidation, since the first requires protons while the second produces them (39). Thus, in this section, the community can still opt not to affect the pH of its environment.

In the Wenelen plasticity section, all reactions associated with this bacterium exhibit structural plasticity. Here Cutipay shifts from no plasticity for Fe(II) oxidation in cluster 20 (where it is abundant) to structural plasticity in cluster 16. Also, we observe a less active Fe(II) oxidation (cut_EX_h_e, from - to +). Analogously, in the Cutipay plasticity section its associated reactions display structural plasticity, while Wenelen loses it for thiosulfate (cluster 13) and iron (cluster 9) consumption when approaching $f_{wen} = 1$. Cluster 7 shows an increase in carbon metabolism structural plasticity, where it no longer requires CO₂ consumption to support its growth (cut_EX_co2_e).

In between the above two sections stands the Cutipay-Wenelen plasticity section, where high plasticity is observed for both organisms. It is a desirable area where organisms can adapt to shifts in their environment, maintaining high growth rates. Finally, in the low plasticity section, characterized by the highest growth rates, all reactions lose plasticity, and the high energy requirements force Cutipay to consume organic matter (represented by equivalent carbon units in the form of glucose [Materials and Methods]).

Eleven critical exchange reactions were selected from the previous clustering analysis, producing a partition of 40 areas, among which an organic matter consumption area is clearly defined for Cutipay as well as a much specific description of the low plasticity section (Fig. S8), thus confirming in a deterministic way previous observations regarding plasticity in the abundance-growth space.

To quantitatively illustrate previous results, explore competition for resources (Fig. 5B; Fig. S9), and use of different substrates to support cell growth (Fig. 5C; Fig. S10), we computed the feasible solution space for critical pairs of reactions associated with the consumption of energy sources at specific points of the abundance-growth space. As expected, there is an overall reduction in plasticity while moving into the growth-oriented area, where exchange fluxes for energy sources are perfectly coupled. This result is evident in the points (0.388 and 0.095), where for a given consumption of thiosulfate, there is a unique consumption of iron (Fig. 5B). Plasticity reduction is also observed when moving into extreme compositions of the community, where being more abundant results in having less specific consumption rates. Notably, when the community becomes balanced, the feasible solution space of critical pairs of reactions concentrates more at higher growth rates. This is less evident just from the abundance-growth space observation.

## Impact of organic matter availability on a bioleaching community

Organic compounds are detrimental for some autotrophic bioleaching bacteria (40), toxic for chemolithotrophic (41), having inhibitory effects in iron oxidation (42), thus favoring the abundance of heterotrophs or facultative heterotrophs in bioleaching communities (43).

Scenarios with increasing amounts of organic matter availability were simulated by moving $\alpha$, and partitions of the abundance-growth space were computed for essential bioleaching reactions (Fig. 6). By increasing $\alpha$, a direct effect in the enlargement of the space with the maximum growth rate of the community and an increase of the higher structural plasticity areas are observed (Fig. 6). In Wenelen, this only happens until the maximum uptake of inorganic sources is reached when $\alpha = 0.6$. Interestingly, in Cutipay, the enlargement of the space is given by the mandatory consumption of organic matter (−−), and it exhibits a shift where CO₂ consumption is no longer required to support its growth (gray area in Fig. 6A). Moreover, new areas appeared near the maximum community growth section, where this bacterium is forced to break down organic matter,
<figure><img src="image_9.png" /><figcaption>FIG 5 Quantitative flux coupling analysis for competition for thiosulfate and distribution of resources by *A. ferrooxidans* in the bioleaching community. (A) Selected points in the abundance-growth space. (B) Feasible fluxes for thiosulfate exchanges (EX_tsul_e) for *A. ferrooxidans* (wen) and *S. thermosulfidooxidans* (cut) at the different points marked in panel A. (C) Feasible fluxes of Fe(II) and thiosulfate exchanges (EX.fe2_e and EX_tsul_e) for *A. ferrooxidans* at different points marked in panel A. Fluxes are in units of (mmol/g DW<sub>org</sub>/h).</figcaption></figure>

hence producing CO₂ to survive (Fig. 6). This shows that organic matter is beneficial for achieving higher growth rates and that Cutipay does not present a preference for organic matter over inorganic sources.

This result, at higher growth rates, initially contradicts literature regarding the behavior of bioleaching communities in the presence of organic matter, which is characterized by an increased abundance of heterotrophs and facultative heterotrophs, which favor the consumption of carbon sources for energy production (43) and will be explored by imposing additional objectives. This could mean that these communities grow at low rates.
<figure><img src="image_10.png" /><figcaption>A</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

**FIG 6** Effects of organic matter availability in the bioleaching community. A partition of the abundance-growth space is computed for key reactions of the bioleaching community when increasing organic matter availability (α). α represents the ratio of the biomass generated by the community that could be degraded to serve as an additional carbon source. The reference case is highlighted with thicker edges; dashed lines show the maximum community growth rate for α = 0. (A) Qualitative states of reactions of *S. thermosulfidooxidans* Cutipay (cut). (B) Qualitative states of reactions of *A. ferrooxidans* Wenelen (wen). Partition tables are ordered according to the location of different zones in the abundance-growth space from lower to higher growth rates, from *f*<sub>wen</sub> 0 to 1, and their occurrence in different scenarios from left to right. EX: exchange reactions, fe2: Fe(II), tsul: thiosulfate, h: hydrogen (proton).

**Impact of substrate composition on the bioleaching community**

Different stages of the bioleaching process are characterized by changes in chemical and physical properties, which shape the distribution of resources, metabolic states, and community compositions (43). In particular, uniform distribution of substrates is not
always guaranteed, and nutrient bioavailability could pose a bottleneck for bioleaching efficiency.

Consistently, in our model, changing the relative availability of Fe(II) and thiosulfate has drastic effects on the behavior of the community (Fig. 7). Iron-predominant scenarios present new metabolic states where iron availability overpasses the effect of preference for thiosulfate as an energy source, which has been observed up to this point in these analyses. It is worth noticing that sulfur compounds have higher energy yields than Fe(II) since they have more electrons available (44). Although, *A. ferrooxidans* is believed to exhibit a preference for iron (45), evidence shows the simultaneous expression of genes

<figure><img src="image_14.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**FIG 7** Effect of relative substrate availability in the bioleaching community. A partition of the abundance-growth space defined by qualitative states for key reactions of the bioleaching community in scenarios with changing ratios between both energy sources. Two partition tables for each bacteria are displayed for the qualitative states in different scenarios. (A) Partition tables of reactions for *S. thermosulfidooxidans* Cutipay. (B) Partition tables of reactions for *A. ferrooxidans* Wenelen. The reference case is highlighted with thicker edges. Partition tables are ordered according to the location of different zones in the abundance-growth space from lower to higher growth rates, from $f_{wen}$ 0 to 1 and from their occurrence in different scenarios from left to right.</figcaption></figure>
for the consumption of iron and sulfur (46), being solubilization of sulfur compounds is key to consuming both energy sources (47). This result is consistent with our analysis, where a preference for Fe(II) oxidation is only displayed when its availability surpasses thiosulfate.

## Impact of considering alternative objectives

For testing a more realistic representation of a bioleaching community, a parsimonious FBA, where minimization of the sum of fluxes is performed as a proxy for energetic efficiency. This has surged as a realistic alternative to retrieve flux distributions (29, 48). This analysis allows a 10%, 20%, and 50% deviation from the optimal solution to maintain a certain degree of plasticity. The results (Fig. 8) show a significant change with the one observed in Fig. S8.

A

<figure><img src="image_16.png" /></figure>

B

<figure><img src="image_17.png" /><figcaption>**FIG 8** Effect on qualitative states when optimizing energetic efficiency in the bioleaching community. A partition of the abundance-growth space defined by qualitative states for key reactions of the bioleaching community in scenarios where 10%, 20%, and 50% of deviation from the optimal solution for energy efficiency is allowed. (A) Two partition tables are displayed for *S. thermosulfidooxidans* Cutipay (cut) for the qualitative states presented in different scenarios (B). Partition table of reactions for *A. ferrooxidans* Wenelen (wen). Partition tables are ordered according to the location of different zones in the abundance-growth space from lower to higher growth rates, from $f_{wen}$ 0 to 1, and from their occurrence in different scenarios from left to right.</figcaption></figure>
Regarding energy source preferences, it is remarkable how thiosulfate becomes mandatory at significantly lower growth rates as we are close to the optimal energy efficiency value. However, the degrees of plasticity in Fe(II) consumption remain unaltered concerning the reference case in both bacteria. This observation is valid in all areas of the space, except for Cutipay in $f_{wen} \ge 0.8$, where this bacteria retains structural plasticity for both inorganic sources. Interestingly, organic matter consumption by this bacteria is required at lower growth rates, which is in agreement with the literature, where the presence of organic matter is beneficial for heterotrophs that use this energy source for biomass production (43, 49, 50).

DISCUSSION

Environmental communities are subjected to constant shifts that affect their survival. Adaptation to such changes includes gene mutations (51), gene transfer, and even gene loss in collaborative communities (52), as well as the development of regulatory processes that result in efficient resource distribution such as flux distributions that result on increased metabolic plasticity (30, 53). This raises the long-standing discussion of whether the community **opts** to live at suboptimal growth rates to increase the capacity to use alternative metabolic pathways and if it is possible to predict their metabolic phenotypes at a given condition, as proposed in references 29 and 30.

Constraint-based modeling for single organisms has historically based most of their flux predictions on the maximization of biomass production (54). This assumption yields promising results for organisms well characterized in laboratory settings for being efficient in the growth and synthesis of a product of interest (55, 56). In contrast, organisms in nature are exposed to constant shifts in their environment, which implies the production of secondary metabolites that provide favorable conditions for survival, deviating resources to the detriment of growth maximization (5, 29). Additionally, communities interact in a complex manner (51, 52), where the abundance of each member has implications in the growth rate of the community.

Several optimization strategies have been developed to characterize interactions between organisms (12–15). From a more ecological perspective, the notion of metabolic niche was proposed to characterize all possible environments in which an organism could live in suboptimal conditions (33). This niche concept uses the feasible flux space to study how organisms interact with their environment, represented by few exchange reactions, depicting new relationships between genotype and metabolic phenotypes (33). Even though these methods consider all feasible metabolic fluxes, they do not fully explore emergent properties arising from organism relative abundances.

In this work, we introduce the abundance-growth space defined by the composition of a community and its growth rate, as an approach to characterize the metabolic phenotype attained by a community and the plasticity of its metabolism. More precisely, we describe, on each point of this space, which reactions are required to carry flux (mandatory vs optional flux) and whether their function can be partially replaced by alternative pathways (variable vs fixed flux). Since these characteristics vary across the abundance-growth space, partitioning this space into zones of similar metabolic profiles gives a characterization of the most relevant phenotypes presented and their relation with the composition of the community.

In addition, the partition of the abundance-growth space allows us to pinpoint critical reactions involved in significant changes in interactions of the community with their environment and establishes, in a precise manner, where critical changes in community plasticity occur. Notably, in both examples developed here, this method allows to unveil how plasticity changes in the abundance-growth space, not only confirming that plasticity increases as the growth rate decreases as expected (29) but also highlighting the relevance of the relative abundance of its members in the loss/gain of plasticity and in the flux coupling of some reactions. Additionally, in these examples only few reactions are found to explain most of the metabolic shifts and many reactions become fixed. Particularly, reactions related to inorganic energy sources exchange and organic
matter consumption in the bioleaching community, and reactions for leucine and lysine exchange and energetic metabolism in the synthetic E. coli community are critical for community plasticity. Certainly, this number can increase as the environment becomes more complex.

However, the characterization of reactions solely based on their qualitative values does not allow to recover the degree of coupling between reactions. This is complemented by performing a quantitative flux coupling analysis, which illustrates how the coupling between reactions changes along the abundance-growth space. In both examples, we show that despite the different degrees of plasticity, there are strong couplings when reaching high growth rates or attaining very unbalanced communities.

The examples presented show the versatility of the method to study changes in the metabolic phenotype in communities with different interactions. On the one hand, those that supply each other with metabolites (cross-feeding) and on the other hand, communities where their complexity lies in variable capacities of its members to consume energy sources, such as the bioleaching community where autotroph and heterotroph bacteria coexist. Interestingly, the abundance-growth space shows a very different behavior in different cases.

An obvious observation of the method has to do with its scalability when the size of the community grow, since the number of grid points can increase exponentially. It is worth mentioning that, both the computation of the grid points in the abundance-growth space and the categorical values of the reactions associated to each feasible grid point admit a parallel implementation, which makes this method usable for communities of few organisms. However, for bigger communities, the method it is still versatile. For instance, one could compute the abundance-growth space for selected combinations of two or three organisms, fixing the abundance of the remaining members of the community, at a reasonably computing cost and still keeping a suitable visualization, or one could reduce this analysis to a fixed number of interesting points of the abundance-growth space.

In any case, we think this method provides a valuable tool for studying, in a detailed manner, different metabolic scenarios where a small community decides to live given a certain environmental condition and to reveal the critical drivers for such functioning, considering the relevance of their relative abundances. In addition, we believe that it could provide a useful approach for the design of synthetic communities of high impact in biotechnology and medicine.

# MATERIALS AND METHODS

## Modeling microbial community metabolism

An expansion of constraint-based modeling for microbial communities was proposed by Koch et al. (13). In summary, N organisms belonging to a community are considered. For each organism $i \in \{1, ..., N\}$, a single model comprised a set of $M_i$ metabolites differentiated by organism and compartment and a set $R_i$ of all reactions of organism i. These reactions include all internal reactions, transport between compartments, and exchange of metabolites with their media. A biomass reaction $r_{biom_i}$ is also included in $R_i$, whose stoichiometric values on their substrates correspond to the amounts (in millimoles) of metabolites present in 1 g of dry weight of biomass of the organism i (g DWi).

If $v_r$ denotes the flux over reaction $r$ expressed in (mmol/gDWi/h), then the following constraints exist for organism i:

$$ LB_r \le v_r \le UB_r, \text{ for all } r \in R_i, \\ \sum_{r \in R_i} S_{mr}v_r = 0, \text{ for all } m \in M_i. $$
The constants LB<sub>r</sub> and UB<sub>r</sub> are the lower and upper bounds for the flux over a reaction *r* and are used to control the amount and direction of fluxes in the system.

Given *N* single-organism models, a community is represented as a single compart-
mentalized system as previously established in references 16 and 57. In this model,
each metabolic network is incorporated as a meta-compartment with an additional pool
compartment that represents the shared media where metabolites can be exchanged
among community members and the exterior.

A set of pool metabolites M<sub>pool</sub> is created containing metabolites present in the extracellular compartment of at least one organism. Consequently, any exchange reaction of organism *i* is considered a transport reaction between the extracellular compartment of *i* and the pool. Finally, a new set R<sup>EX</sup> of exchange reactions is defined to control external conditions for all metabolites in the pool.

In this community model, each member *i* makes up a fraction *f*<sub>*i*</sub> of the total amount of community biomass (with Σ<sup>*N*</sup><sub>*i* = 1</sub> *f*<sub>*i*</sub> = 1). These fractions are relevant to compare the feasible fluxes on reactions of different organisms accurately.

Indeed, in this community model, fluxes are expressed in a single unit (mmol/g DW<sub>com</sub>/h), where g DW<sub>com</sub> is a gram of dry weight of the community biomass. Since bounds on reaction fluxes of organism *i* were originally expressed in (mmol/g DW<sub>i</sub>/h) in the single-organism model, they must be recomputed in the community model, which is done by weighting by *f*<sub>*i*</sub> each bound previously defined on a flux reaction of organism *i* (15).

Given this unit conversion, the flux v<sub>biom_i</sub> of the biomass reaction of organism i is not expressing the growth μ<sub>i</sub> of organism i, but the grams of organism i produced by each gram of the community per hour. Hence, if we consider a state of balanced growth, where the fraction f<sub>i</sub> of each organism i is maintained over time, then we should impose that all organisms are growing at the same rate: μ<sub>1</sub> = ... = μ<sub>N</sub> = μ<sub>com</sub>, which is equivalent to imposing: v<sub>biom_i</sub> = f<sub>i</sub> μ<sub>com</sub> for all i ∈ {1, ..., N} (15).

The balanced growth assumptions indicate that the community has reached a stable state that allows the organisms to maintain their relative abundances over time, which is a common assumption in community modeling (13, 15, 16). If this constrain is not assumed, the faster-growing organism will ultimately displace all other organisms in the community. The requirement of identical growth rate between all organisms applies not necessarily at every time point but averaged over time, allowing small variations in the relative abundances. Thus, this hypothesis is a reasonable approximation for communities where composition is mainly maintained between spaced time points (15).

Considering all these observations, the set of constraints on the fluxes v<sub>r</sub> of the
community are the following:

$$
\begin{align*}
\sum_{r \in R_i} S_{mr} v_r &= 0, \quad m \in M_v, i \in 1, \dots, N, \\
V_{m\_pool\_exchange} &= \sum_{i=1}^{N} V_{m\_extracell\_i \to pool}, \quad m \in M_{pool} \\
f_i L B_r \le v_r \le f_i U B_r, \quad r \in R_i, i \in \{1, \dots, N\} \\
LB_{m\_pool\_exchange} \le V_{m\_pool\_exchange} \le UB_{m\_pool\_exchange}, \quad m \in M_{pool} \\
\mu_{com} &= \sum_{i=1}^{N} V_{bio\_i} \\
v_{biom\_i} &= f_i \mu_{com}, \quad i \in \{1, \dots, N\}.
\end{align*}
$$

If the fractions *f*<sub>*i*</sub> are fixed (as in cFBA [16]) or if the community growth rate *μ*<sub>com</sub> is
fixed (as in the case of SteadyCom [15] and RedCom [13]), then all constraints are linear.

Eventually, any other interaction between the organisms or additional information that can be expressed as a linear constraint on the fluxes can be easily added to the model. For instance, the model can directly include fixing values to reaction fluxes or imposing coupling between reactions. In the case of the bioleaching community, a
pseudo-reaction models the capacity of one organism to use degraded organic matter as an additional energy source.

## Analysis of the abundance-growth space

The definition and analysis of the abundance-growth space are performed as follows.

**Definition of abundance-growth space.** Considering N organisms in the community, an abundance-growth pair is given by $(F, \mu) = (f_1, ..., f_N, \mu) \in \mathbb{R}^{N+1}$, where $f_1, ..., f_N$ are the relative abundances of the organisms (with $\sum_{i=1}^N f_i = 1$) and $\mu$ is the growth rate of the community (15). Each pair $(F, \mu)$ defines a polytope $P_{F,\mu}$ of flux distributions satisfying all constraints of the model for the abundance F and growth rate $\mu$. We say that the pair $(F, \mu)$ is feasible if given F and $\mu$, at least one flux distribution satisfies all constraints of the community model, that is, $P_{F,\mu} \neq \emptyset$. We define the abundance-growth space as the set of all feasible abundance-growth pairs. We use the term "space" to remark that each point $(F, \mu)$ is a representation of the set of $P_{F,\mu}$ of feasible fluxes. It is important to note that, since the constraints related to the flux of biomass of each organism are not linear when both F and $\mu$ are variables, then the set of feasible abundance-growth points is not necessarily convex, as shown in the synthetic E. coli community example. However, for fixed values of the vector F or $\mu$ they become linear, and so, the convexity is assured. This means that convexity is assured in horizontal and vertical sections of the abundance-growth space.

**Defining a grid of feasible abundance-growth points.** A set $\mathcal{F}$, corresponding to a discretization of the abundance-growth space, is defined as follows. Specifically, for a given resolution parameter $\ell$, we define a discretization of the possible relative abundances values by considering the set $\overline{[0, 1]}$ corresponding to $\ell$ equidistant values in the interval $[0, 1]$. Thus, the set $\mathcal{A}$ of relative abundance points $(f_1, ..., f_N)$ that we consider is defined by:

$$ \mathcal{A} = \{(f_1, ..., f_N) \in \overline{[0, 1]}^N \mid \sum_{i=1}^N f_i = 1\}. $$

A second discretization is done for the community growth rate $\mu$ by defining the set $\overline{[0, \text{MAX}_{\text{biomass}}]}$ corresponding to $\ell$ equidistant values in the interval $[0, \text{MAX}_{\text{biomass}}]$, where $\text{MAX}_{\text{biomass}}$ is the maximum biomass computed in the model for all abundance points in $\mathcal{A}$. The set of abundance-growth points considered in the analysis is defined by:

$$ \mathcal{A} \times \mathcal{G} = \{(f_1, ..., f_N, \mu) \in \overline{[0, 1]}^N \times \overline{[0, \text{MAX}_{\text{biomass}}]} \mid \sum_{i=1}^N f_i = 1\}. $$

Finally, the set of points that define the grid $\mathcal{F}$ are those that are feasible. That is, the points $(f_1, ..., f_N, \mu) \in \mathcal{A} \times \mathcal{G}$ such that the community can reach a growth $\mu$ for the given relative abundances $(f_1, ..., f_N)$. Note that, since $\sum_{i=1}^N f_i = 1$, then an abundance of one organism can be omitted in an abundance-growth point, which is especially useful in the case when $N = 2$, since this implies that abundance-growth points $(f_1, \mu)$ can be depicted in a two-dimensional plot.

**Computing qualitative vectors for the abundance-growth points.** Each point $(F, \mu)$ of the abundance-growth space is a representation of all the flux distributions in its associated polytope $P_{F,\mu}$. Given a point $(F, \mu) \in \mathcal{F}$ in the abundance-growth space, the range of feasible fluxes of a reaction $r$ is exactly the interval $[\min_{F,\mu}(r), \max_{F,\mu}(r)]$, where $\min(r) = \min\{v_r \mid \text{for all } v \in P_{F,\mu}\}$, $\max(r) = \max\{v_r \mid \text{for all } v \in P_{F,\mu}\}$. Both amounts are easily determined in the community model by running flux variability analysis (FVA) (31) fixing F and $\mu$ in the constraint-based model. Thus, the interval $[\min_{F,\mu}(r), \max_{F,\mu}(r)]$ is computed for all reactions $r$ in each point $(F, \mu)$ of the grid
F. Computation of these amounts can be done in parallel or by efficiently implementing consecutive LP formulations with the same feasible space (58). The range [min<sub>F,μ</sub>(r), max<sub>F,μ</sub>(r)] is used to give a categorical value to each reaction *r* and abundance-growth point (*F*, *μ*) (Table 1). These qualitative values are encapsulated in the vector *X*(*F*, *μ*) ∈ ℂ<sup>|R|</sup> of dimension equal to the number of analyzed reactions in the community model, where ℂ is the set of categories. In the applications of this method, a subset *R'* ⊂ *R* of selected reactions will be considered, and thus we will focus our analysis on the restriction of *X*(*F*, *μ*) to the coordinates associated with *R'*, say *X*′(*F*, *μ*) ∈ ℂ<sup>|R'|</sup>.

Partition and cluster-partition of the abundance-growth space. As stated previously, the vectors X'(F, μ) ∈ ℂ^{|R|} are qualitative descriptions of the set of reactions on each point (F, μ) of the grid F. We define two ways of partitioning this grid. First, using the categorical vectors X'(F, μ), a partition of the grid can be defined by stating that two feasible abundance-growth points in the grid (F₁, μ₁) and (F₂, μ₂) are equivalent if X'(F₁, μ₁) = X'(F₂, μ₂). Thus, each set of the partition can be considered a zone where selected reactions in the model have identical categorical descriptions. It is called the partition of the abundance-growth space defined by the selected reactions in R'.

Depending on the size of the set of selected reactions $R'$, the number of zones in this partition can be huge. To address this issue, we propose to use a coarser partition based on the degree of similarity of the categorical vectors. More precisely, we perform a classification of the vectors $X'(F, \mu)$ for each $(F, \mu) \in \mathcal{F}$ and $r \in R'$ using a hierarchical clustering algorithm with a Jaccard distance, where the number of clusters ($k$) is defined a priori. In other words, we can identify zones of the grid where reactions have a similar categorical behavior according to the clustering method used. We call this partition the cluster-partition of the abundance-growth space (defined by the selected reactions in $R'$). Alternatively, this method can be implemented using other clustering methods.

Additionally, we want to summarize the qualitative state for each reaction on a given zone defined by a cluster partition. To achieve this, for each reaction, given a cluster, we assign a qualitative value to each reaction if such a qualitative state is present in over 80% of the points of that cluster.

*Determination of the number of clusters for analysis.* Determination of the number of clusters to be used is performed as follows: for a given number *k* of clusters, each categorical vector *X*'(*F*, *μ*) is compared to its assigned cluster-partition descriptor. For each reaction, a score value between 0 and 1 represents the fraction of points consistent with the cluster-partition descriptor. This analysis is performed with different values of *k* clusters to retrieve the minimum number of clusters where all analyzed reactions are represented correctly in at least 80% of points of the grid.

*Quantitative flux coupling analysis.* Given an abundance-growth point, the space of feasible fluxes for two given reactions *r*₁ and *r*₂ is computed by defining a homogeneous grid of 50 values between the minimum and maximum value of the flux of *r*₁ computed by FVA. Then, for each flux value of *r*₁ in this grid, an FVA is performed for reaction *r*₂ to compute the minimum and maximum flux values of *r*₂. The interval defined by these values is the feasible fluxes of *r*₂ for the given flux value of *r*₁. Normalization by *f*ᵢ is performed to retrieve flux values in units of (mmol/g DWᵢ/h) for each analyzed reaction.

Analysis of a synthetic *E. coli* community

A community defined by two strains of *Escherichia coli* auxotrophic for leucine and lysine was selected from reference 37. The *E. coli* genome-scale model iAF1260 (59) was modified by constraining reactions associated to gene knock-outs made in reference 37: diaminopimelate decarboxylase (DAPDC) for the lysine auxotrophic strain (eco_K) and 2-isopropylmalate synthase (IPPS) for the leucine auxotrophic strain (eco_L). Glucose availability was set up to be 10 (mmol/g DW<sub>org</sub>/h), amino acid production was set up to be 0.2 (mmol/g DW<sub>org</sub>/h) for both leucine and lysine, and consumption for both amino acids for both strains was constrained according to reported values of production and
pseudo-reactions representing the degradation of community biomass and its transformation into carbon equivalent units of glucose:

$$ \text{Biomass}_{\text{community}} \rightarrow \alpha \text{Biomass}_{\text{degraded}} + (1 - \alpha) \text{Biomass}_{\text{net}} $$

$$ \text{Biomass}_{\text{degraded}} \rightarrow 6.169 \text{glc}^{\text{pool}} $$

where parameter $\alpha$ represents that the death rate of the organisms is equivalent to a fraction $\alpha$ of the community growth rate. The choice of glucose as a unique representation of organic matter is an arbitrary simplification in the absence of more information about how the organisms are degraded. This was done because glucose can be easily assimilated by the metabolic network (File S3, https://github.com/mathomics/ecosystem), as well as that it can be incorporated to synthesize nucleotides, amino acids, as well as other simpler organic acids found in bioleaching environments, such as acetate.

The set of 1,717 reactions of this community model comprises 1,668 reactions from the original single models and 49 reactions that appear in the construction of the community model. The set of 1,611 metabolites of the community model comprises 1,564 metabolites from single models and 47 metabolites from the community construction.

*Definition of parameters and scenarios for simulations of bioleaching communities.* Different external and internal conditions are explored to determine how these scenarios change the distribution of qualitative states in the abundance-growth space. External conditions are represented by available resources, both organic and inorganic, for the community. Meanwhile, internal conditions correspond to alternative requirements for the obtained flux distributions, mainly being energetically efficient.

Modeling the availability of organic matter is achieved by setting the parameter $\alpha$. The availability of inorganic sources is modeled, defining the available fluxes of Fe(II) and thiosulfate in the community. It is achieved by adjusting the lower bounds for the community exchange reactions of both sources.

Thus, we define the values $\text{Max}_{\text{Fe(II)}} = 150$ (mmol/g $\text{DW}_{\text{com}}$/h) and $\text{Max}_{\text{tsul}} = 10$ (mmol/g $\text{DW}_{\text{com}}$/h) as reference values, which correspond to the uptakes of Fe(II) and thiosulfate that produce the maximum level of biomass on both organisms according to their single models.

Any convex combination of these amounts generates an equivalent amount of inorganic sources. Formally, a relative combination of $\lambda$ Fe(II) and $1 - \lambda$ thiosulfate means that there is a flux of $\lambda \text{Max}_{\text{Fe(II)}}$ of Fe(II) and a flux of $(1 - \lambda) \text{Max}_{\text{tsul}}$ of thiosulfate available to the community. With these definitions, the following scenarios are defined:

1. A reference case is defined as having a low organic availability ($\alpha = 0.2$), and a relative combination of 0.5 Fe(II) and 0.5 thiosulfate ($\lambda = 0.5$).

2. Analysis of the effects of organic matter is performed by changing $\alpha$ in the reference case scenario ($\lambda = 0.5$). The analyzed values of $\alpha$ are 0, 0.2, 0.4, 0.6, and 0.8.

3. Analysis of the effect of substrate composition is performed by variations in the value of $\lambda$ in the reference case ($\alpha = 0.2$). The analyzed values of $\lambda$ are 0.25, 0.4, 0.5, 0.6, and 0.75.

4. Analysis of alternative objective functions is performed by analyzing the reference case ($\alpha = 0.2, \lambda = 0.5$) when additional functions are required. Specifically, energetic efficiency is represented as an additional constraint which states that the sum of fluxes of all reactions is less than a factor of the minimum sum of fluxes. The analyzed values for this factor are 1.1, 1.2, and 1.5.

ACKNOWLEDGMENTS

This work was supported by the Center for Mathematical Modeling (CMM) grant FB210005, which is a Basal fund for centers of excellence from ANID-Chile, Center for Genome Regulation, which is the Millennium Institute Project ICN2021 044 supported
by the Millennium Scientific Initiative of the Ministry of Economy, Development and Tourism (Chile), Grant Exploración no. 13220002, and the National Laboratory for High Performance Computing (NLHPC) at CMM.

## AUTHOR AFFILIATIONS

¹Center for Mathematical Modeling, University of Chile, Santiago, Chile

²Center for Genome Regulation, Millennium Institute, University of Chile, Santiago, Chile

³Nantes Université, Ecole Centrale Nantes, CNRS, Nantes, France

⁴Department of Mathematical Engineering, University of Chile, Santiago, Chile

## AUTHOR ORCIDs

Natalia E. Jiménez ID http://orcid.org/0000-0003-2084-1256

Alejandro Eduardo Maass ID http://orcid.org/0000-0002-7038-4527

## FUNDING

<table><thead><tr><th>Funder</th><th>Grant(s)</th><th>Author(s)</th></tr></thead><tbody><tr><td>Agencia Nacional de Investigación y Desarrollo (ANID)</td><td>FB2100005</td><td>Alejandro Eduardo Maass<br>Natalia E. Jiménez<br>Maria Paz Cortés<br>Vicente Acuña</td></tr><tr><td>Ministerio de Economía, Fomento y Turismo, Chile (Ministerio de Economía, Fomento y Turismo)</td><td>ICN2021 044</td><td>Vicente Acuña<br>Natalia E. Jiménez<br>Alejandro Eduardo Maass<br>Maria Paz Cortés</td></tr></tbody></table>

## AUTHOR CONTRIBUTIONS

Natalia E. Jiménez, Conceptualization, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing | Vicente Acuña, Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review and editing | María Paz Cortés, Data curation, Investigation, Validation | Damien Eveillard, Investigation, Writing – original draft, Writing – review and editing | Alejandro Eduardo Maass, Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – original draft, Writing – review and editing

## DATA AVAILABILITY

This method was implemented in Python using COBRApy for constructing community models, FBA, and FVA simulations (75), and it is available at Github.

## ADDITIONAL FILES

The following material is available online.

### Supplemental Material

**Fig. S1 (mSystems00492-23_S0001.pdf).** Effect of modifying production capacity and supplementation of leucine and lysine in the abundance-growth space of the *E. coli* synthetic community.

**Fig. S2 (mSystems00492-23_S0002.pdf).** Partition of the abundance-growth space determined by exchange of leucine and lysine in a synthetic *E. coli* community.

**Fig. S3 (mSystems00492-23_S0003.pdf).** Determination of the number of clusters for a synthetic *E. coli* community.
Fig. S4 (mSystems00492-23_S0004.pdf). Partition of the abundance-growth space determined by key reactions in a synthetic *E. coli* community.

Fig. S5 (mSystems00492-23_S0005.pdf). Effect of modifying leucine and lysine supplementation in leucine production and export by eco_K.

Fig. S6 (mSystems00492-23_S0006.pdf). Quantitative flux coupling analysis for exchanges of lysine and leucine in the synthetic *E. coli* supplemented with these amino acids.

Fig. S7 (mSystems00492-23_S0007.pdf). Effect of the variation of the number of clusters in the bioleaching community.

Fig. S8 (mSystems00492-23_S0008.pdf). Partition of the abundance-growth space determined by 11 selected key bioleaching reactions describes relevant shifts in the bioleaching community.

Fig. S9 (mSystems00492-23_S0009.pdf). Quantitative flux coupling analysis for iron(II) competition in the bioleaching community.

Fig. S10 (mSystems00492-23_S00010.pdf). Quantitative flux coupling analysis for the distribution of resources by *S. thermosulfidooxidans* in the bioleaching community.

## REFERENCES

1. Madsen JS, Sørensen SJ, Burmølle M. 2018. Bacterial social interactions and the emergence of community-intrinsic properties. Curr Opin Microbiol 42:104–109. https://doi.org/10.1016/j.mib.2017.11.018

2. Røder HL, Herschend J, Russel J, Andersen MF, Madsen JS, Sørensen SJ, Burmølle M. 2018. Enhanced bacterial mutualism through an evolved biofilm phenotype. ISME J 12:2608–2618. https://doi.org/10.1038/s41396-018-0165-2

3. Raghupathi PK, Liu W, Sabbe K, Houf K, Burmølle M, Sørensen SJ. 2018. Synergistic interactions within a multispecies biofilm enhance individual species protection against grazing by a pelagic protozoan. Front. Microbiol 8:2649. https://doi.org/10.3389/fmicb.2017.02649

4. Pande S, Merker H, Bohl K, Reichelt M, Schuster S, de Figueiredo LF, Kaleta C, Kost C. 2014. Fitness and stability of obligate cross-feeding interactions that emerge upon gene loss in bacteria. ISME J 8:953–962. https://doi.org/10.1038/ismej.2013.211

5. Shade A, Peter H, Allison SD, Baho DL, Berga M, Bürgmann H, Huber DH, Langenheder S, Lennon JT, Martiny JBH, Matulich KL, Schmidt TM, Handelsman J. 2012. Fundamentals of microbial community resistance and resilience. Front Microbiol 3:417. https://doi.org/10.3389/fmicb.2012.00417

6. Tara Ocean Foundation, Abreu A, Bourgois E, Gristwood A, Troublé R, Tara Oceans, Acinas SG, Bork P, Boss E, Bowler C, Budinich M, Chaffron S, de Vargas C, Delmont TO, Eveillard D, Guidi L, Iudicone D, Kandels S, Morlon H, Lombard F, Pepperkok R, Karlusich JJP, Piganeau G, Régimbeau A, Sommeria-Klein G, Stemmann L, Sullivan MB, Sunagawa S, Wincker P, Zablocki O, European Molecular Biology Laboratory (EMBL), Arendt D, Bilic J, Finn R, Heard E, Rouse B, Vamathevan J, European Marine Biological Resource Centre - European Research Infrastructure Consortium (EMBR-C-ERIC), Casotti R, Cancio I, Cunliffe M, Kervella AE, Kooistra WHCF, Obst M, Pade N, Power DM, Santi I, Tsagaraki TM, Vanaverbeke J. 2022. Priorities for ocean microbiome research. Nat Microbiol 7:937–947. https://doi.org/10.1038/s41564-022-01145-5

7. Nicholson JK, Holmes E, Kinross J, Burcelin R, Gibson G, Jia W, Pettersson S. 2012. Host-gut microbiota metabolic interactions. Science 336:1262–1267. https://doi.org/10.1126/science.1223813

8. Lindemann SR, Bernstein HC, Song H-S, Fredrickson JK, Fields MW, Shou W, Johnson DR, Beliaev AS. 2016. Engineering microbial consortia for controllable outputs. ISME J 10:2077–2084. https://doi.org/10.1038/ismej.2016.26

9. Falkowski PG, Fenchel T, Delong EF. 2008. The microbial engines that drive earth's biogeochemical cycles. Science 320:1034–1039. https://doi.org/10.1126/science.1153213

10. Zomorrodi AR, Segrè D. 2016. Synthetic ecology of microbes: mathematical models and applications. J Mol Biol 428:837–861. https://doi.org/10.1016/j.jmb.2015.10.019

11. Gottstein W, Olivier BG, Bruggeman FJ, Teusink B. 2016. Constraint-based stoichiometric modelling from single organisms to microbial

communities. J R Soc Interface 13:124. https://doi.org/10.1098/rsif.2016.0627

12. Zelezniak A, Andrejev S, Ponomarova O, Mende DR, Bork P, Patil KR. 2015. Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proc Natl Acad Sci USA 112:6449–6454. https://doi.org/10.1073/pnas.1421834112

13. Koch S, Kohrs F, Lahmann P, Bissinger T, Wendschuh S, Benndorf D, Reichl U, Klamt S. 2019. Redcom: a strategy for reduced metabolic modeling of complex microbial communities and its application for analyzing experimental datasets from anaerobic digestion. PLOS Comput Biol 15:e1006759. https://doi.org/10.1371/journal.pcbi.1006759

14. Belcour A, Frioux C, Aite M, Bretaudieu A, Hildebrand F, Siegel A. 2020. Metage2Metabo, microbiota-scale metabolic complementarity for the identification of key species. Elife 9:e61968. https://doi.org/10.7554/eLife.61968

15. Chan SHJ, Simons MN, Maranas CD. 2017. Steadycom: predicting microbial abundances while ensuring community stability. PLOS Comput Biol 13:e1005539. https://doi.org/10.1371/journal.pcbi.1005539

16. Khandelwal RA, Olivier BG, Röling WFM, Teusink B, Bruggeman FJ. 2013. Community flux balance analysis for microbial consortia at balanced growth. PLOS ONE 8:e64567. https://doi.org/10.1371/journal.pone.0064567

17. Diener C, Gibbons SM, Resendis-Antonio O. 2020. MICOM: metagenome-scale modeling to infer metabolic interactions in the gut microbiota. mSystems 5:e00606-19. https://doi.org/10.1128/mSystems.00606-19

18. Basile A, Campanaro S, Kovalovszki A, Zampieri G, Rossi A, Angelidaki I, Valle G, Treu L. 2020. Revealing metabolic mechanisms of interaction in the anaerobic digestion microbiome by flux balance analysis. Metab Eng 62:138–149. https://doi.org/10.1016/j.ymben.2020.08.013

19. Damiani C, Maspero D, Di Filippo M, Colombo R, Pescini D, Graudenzi A, Westerhoff HV, Alberghina L, Vanoni M, Mauri G. 2019. Integration of single-cell RNA-Seq data into population models to characterize cancer metabolism. PLoS Comput Biol 15:e1006733. https://doi.org/10.1371/journal.pcbi.1006733

20. Henson MA, Phalak P. 2018. Suboptimal community growth mediated through metabolite crossfeeding promotes species diversity in the gut microbiota. PLOS Comput Biol 14:e1006558. https://doi.org/10.1371/journal.pcbi.1006558

21. Hirmas B, Gasaly N, Orellana G, Vega-Sagardía M, Saa P, Gotteland M, Garrido D. 2022. Metabolic modeling and bidirectional culturing of two gut microbes reveal cross-feeding interactions and protective effects on intestinal cells. mSystems 7:e0064622. https://doi.org/10.1128/mSystems.00646-22

22. Orth JD, Thiele I, Palsson BØ. 2010. What is flux balance analysis? Nat Biotechnol 28 (3):245–248. 10.1038/nbt.1614.

23. Budinich M, Bourdon J, Larhlimi A, Eveillard D. 2017. A multi-objective constraint-based approach for modeling genome-scale microbial
ecosystems. PLOS ONE 12:e0171744. https://doi.org/10.1371/journal.pone.0171744

24. Zomorrodi AR, Maranas CD. 2012. Optcom: A multi-level optimization framework for the metabolic modeling and analysis of microbial communities. PLOS Comput Biol 8:e1002363. https://doi.org/10.1371/journal.pcbi.1002363

25. Yoshikawa K, Toya Y, Shimizu H. 2017. Metabolic engineering of *Synechocystis* sp. PCC 6803 for enhanced ethanol production based on flux balance analysis. Bioprocess Biosyst Eng 40 (5):791–796. 10.1007/s00449-017-1744-8.

26. Alper H, Jin Y-S, Moxley JF, Stephanopoulos G. 2005. Identifying gene targets for the metabolic engineering of lycopene biosynthesis in *Escherichia coli*. Metab Eng 7:155–164. https://doi.org/10.1016/j.ymben.2004.12.003

27. Bruggeman FJ, Planqué R, Molenaar D, Teusink B. 2020. Searching for principles of microbial physiology. FEMS Microbiol Rev 44:821–844. https://doi.org/10.1093/femsre/fuaa034

28. Mas A, Jamshidi S, Lagadeuc Y, Eveillard D, Vandenkoornhuyse P. 2016. Beyond the black queen hypothesis. ISME J 10:2085–2091. https://doi.org/10.1038/ismej.2016.22

29. Schuetz R, Zamboni N, Zampieri M, Heinemann M, Sauer U. 2012. Multidimensional optimality of microbial metabolism. Science 336:601–604. https://doi.org/10.1126/science.1216882

30. Fischer E, Sauer U. 2005. Large-scale *in vivo* flux analysis shows rigidity and suboptimal performance of *Bacillus subtilis* metabolism. Nat Genet 37:636–640. https://doi.org/10.1038/ng1555

31. Mahadevan R, Schilling CH. 2003. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng 5:264–276. https://doi.org/10.1016/j.ymben.2003.09.002

32. Edwards JS, Ramakrishna R, Palsson BO. 2002. Characterizing the metabolic phenotype: a phenotype phase plane analysis. Biotechnol Bioeng 77:27–36. https://doi.org/10.1002/bit.10047

33. Régimbeau A, Budinich M, Larhlimi A, Pierella Karlusich JJ, Aumont O, Memory L, Bowler C, Eveillard D. 2022. Contribution of genome-scale metabolic modelling to niche theory. Ecol Lett 25:1352–1364. https://doi.org/10.1111/ele.13954

34. Grilli J. 2020. Macroeological laws describe variation and diversity in microbial communities. Nat Commun 11:4743. https://doi.org/10.1038/s41467-020-18529-y

35. Almaas E, Oltvai ZN, Barabási AL. 2005. The activity reaction core and plasticity of metabolic networks. PLOS Comput Biol 1:e68. https://doi.org/10.1371/journal.pcbi.0010068

36. Samal A. 2008. Conservation of high-flux backbone in alternate optimal and near-optimal flux distributions of metabolic networks. Syst Synth Biol 2:83–93. https://doi.org/10.1007/s11693-009-9025-8

37. Zhang X, Reed JL, Fong SS. 2014. Adaptive evolution of synthetic cooperating communities improves growth performance. PLoS ONE 9:e108297. https://doi.org/10.1371/journal.pone.0108297

38. Latorre M, Cortés MP, Travisany D, Di Genova A, Budinich M, Reyes-Jara A, Hödar C, González M, Parada P, Bobadilla-Fazzini RA, Cambiazo V, Maass A. 2016. The bioleaching potential of a bacterial consortium. Bioresour Technol 218:659–666. https://doi.org/10.1016/j.biortech.2016.07.012

39. Rohwerder T, Sand W. 2007. Oxidation of inorganic sulfur compounds in acidophilic prokaryotes. Eng. Life Sci. 7:301–309. https://doi.org/10.1002/elsc.200720204

40. Fang D, Zhou LX. 2006. Effect of sludge dissolved organic matter on oxidation of ferrous iron and sulfur by *Acidithiobacillus ferrooxidans* and *Acidithiobacillus thiooxidans*. Water Air Soil Pollut 171:81–94. https://doi.org/10.1007/s11270-005-9014-9

41. Nancucheo I, Johnson DB. 2010. Production of glycolic acid by chemolithotrophic iron-and sulfur-oxidizing bacteria and its role in delineating and sustaining acidophilic sulfide mineral-oxidizing consortia. Appl Environ Microbiol 76:461–467. https://doi.org/10.1128/AEM.01832-09

42. Zhou H, Zhang R, Hu P, Zeng W, Xie Y, Wu C, Qiu G. 2008. Isolation and characterization of *Ferroplasma thermophilum* sp. nov., a novel extremely acidophilic, moderately thermophilic archaeon and its role in bioleaching of chalcopyrite. J Appl Microbiol. 105:591–601. https://doi.org/10.1111/j.1365-2672.2008.03807.x

43. Li Q, Tian Y, Fu X, Yin H, Zhou Z, Liang Y, Qiu G, Liu J, Liu H, Liang Y, Shen L, Cong J, Liu X. 2011. The community dynamics of major bioleaching

microorganisms during chalcopyrite leaching under the effect of organics. Curr Microbiol 63:164–172. https://doi.org/10.1007/s00284-011-9960-y

44. Johnson D. 2018. The evolution, current status, and future prospects of using biotechnologies in the mineral extraction and metal recovery sectors. Minerals 8:343. https://doi.org/10.3390/min8080343

45. Ponce JS, Moinier D, Byrne D, Amouric A, Bonnefoy V. 2012. *Acidithiobacillus ferrooxidans* oxidizes ferrous iron before sulfur likely through transcriptional regulation by the global redox responding RegBA signal transducing system. Hydrometallurgy 127–128:187–194. https://doi.org/10.1016/j.hydromet.2012.07.016

46. Kucera J, Bouchal P, Lochman J, Potesil D, Janiczek O, Zdrahal Z, Mandl M. 2013. Ferrous iron oxidation by sulfur-oxidizing *Acidithiobacillus ferrooxidans* and analysis of the process at the levels of transcription and protein synthesis. Antonie Van Leeuwenhoek 103:905–919. https://doi.org/10.1007/s10482-012-9872-2

47. Inaba Y, Kernan T, West AC, Banta S. 2021. Dispersion of sulfur creates a valuable new growth medium formulation that enables earlier sulfur oxidation in relation to iron oxidation in *Acidithiobacillus ferrooxidans* cultures. Biotechnol Bioeng 118:3225–3238. https://doi.org/10.1002/bit.27847

48. Lewis NE, Hixson KK, Conrad TM, Lerman JA, Charusanti P, Polpitiya AD, Adkins JN, Schramm G, Purvine SO, Lopez-Ferrer D, Weitz KK, Eils R, König R, Smith RD, Palsson BØ. 2010. Omic data from evolved *E. coli* are consistent with computed optimal growth from genome-scale models. Mol Syst Biol 6:390. https://doi.org/10.1038/msb.2010.47

49. Saavedra A, Aguirre P, Gentina JC. 2020. Biooxidation of iron by *Acidithiobacillus ferrooxidans* in the presence of D-galactose: understanding its influence on the production of EPS and cell tolerance to high concentrations of iron. Front Microbiol 11:759. https://doi.org/10.3389/fmicb.2020.00759

50. Liu H, Yin H, Dai Y, Dai Z, Liu Y, Li Q, Jiang H, Liu X. 2011. The co-culture of *Acidithiobacillus ferrooxidans* and *Acidiphilium acidiphilum* enhances the growth, iron oxidation, and CO₂ fixation. Arch Microbiol 193:857–866. https://doi.org/10.1007/s00203-011-0723-8

51. Van Valen L. 1977. The red queen. Am Nat 111:809–810. https://doi.org/10.1086/283213

52. Morris JJ, Lenski RE, Zinser ER. 2012. The black queen hypothesis: evolution of dependencies through adaptive gene loss. MBio 3:e00036–12. https://doi.org/10.1128/mBio.00036-12

53. Zampieri M, Enke T, Chubukov V, Ricci V, Piddock L, Sauer U. 2017. Metabolic constraints on the evolution of antibiotic resistance. Mol Syst Biol 13:3. https://doi.org/10.15252/msb.20167028

54. Feist AM, Palsson BO. 2010. The biomass objective function. Curr Opin Microbiol 13:344–349. https://doi.org/10.1016/j.mib.2010.03.003

55. Yang D, Park SY, Park YS, Eun H, Lee SY. 2020. Metabolic engineering of *Escherichia coli* for natural product biosynthesis. Trends Biotechnol. 38:745–765. https://doi.org/10.1016/j.tibtech.2019.11.007

56. Priola JJ, Calzadilla N, Baumann M, Borth N, Tate CG, Betenbaugh MJ. 2016. High-throughput screening and selection of mammalian cells for enhanced protein production. Biotechnol J 11:853–865. https://doi.org/10.1002/biot.201500579

57. Frioux C, Singh D, Korcsmaros T, Hildebrand F. 2020. From bag-of-genes to bag-of-genomes: metabolic modelling of communities in the era of metagenome-assembled genomes. Comput Struct Biotechnol J 18:1722–1734. https://doi.org/10.1016/j.csbj.2020.06.028

58. Gudmundsson S, Thiele I. 2010. Computationally efficient flux variability analysis. BMC Bioinformatics 11:1–3. https://doi.org/10.1186/1471-2105-11-489

59. Feist AM, Henry CS, Reed JL, Krummenacker M, Joyce AR, Karp PD, Broadbelt LJ, Hatzimanikatis V, Palsson BØ. 2007. A genome-scale metabolic reconstruction for *Escherichia coli* K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Mol Syst Biol 3:121. https://doi.org/10.1038/msb4100155

60. Travisany D, Di Genova A, Sepúlveda A, Bobadilla-Fazzini RA, Parada P, Maass A. 2012. Draft genome sequence of the *Sulfobacillus thermosulfidooxidans* cutipay strain, an indigenous bacterium isolated from a naturally extreme mining environment in northern Chile. J Bacteriol 194:6327–6328. https://doi.org/10.1128/JB.01622-12

61. Travisany D, Di Genova A, Sepúlveda A, Bobadilla-Fazzini RA, Parada P, Maass A. 2012. Data from "draft genome sequence of the Sulfobacillus
Thermosulfidooxidans Cutipay strain, an indigenous bacterium isolated from a naturally extreme mining environment in northern Chile. *J Bacteriol* 194:6327–6328. https://doi.org/10.1128/JB.01622-12

62. Latorre M, Ehrenfeld N, Cortés MP, Travisany D, Budinich M, Aravena A, González M, Bobadilla-Fazzini RA, Parada P, Maass A. 2016. Global transcriptional responses of *Acidithiobacillus ferrooxidans* wenelen under different sulfide minerals. *Bioresour Technol* 200:29–34. https://doi.org/10.1016/j.biortech.2015.09.110

63. Aite Méziane, Chevallier M, Frioux C, Trottier C, Got J, Cortés MP, Mendoza SN, Carrier G, Dameron O, Guillaudeux N, Latorre M, Loira N, Markov GV, Maass A, Siegel A. 2018. "Traceability, reproducibility and wiki-exploration for "à-la-carte" reconstructions of genome-scale metabolic models". PLOS Comput Biol 14:e1006146. https://doi.org/10.1371/journal.pcbi.1006146

64. Aite M, Chevallier M, Frioux C, Trottier C, Got J, Cortés MP, Mendoza SN, Carrier G, Dameron O, Guillaudeux N. 2018. Supplemental material for Aite et al. PLOS Comput Biol. https://hub.docker.com/r/dyliss/aureme-img.

65. Loira N, Zhukova A, Sherman DJ. 2015. Pantograph: a template-based method for genome-scale metabolic model reconstruction. *J Bioinform Comput Biol* 13:1550006. https://doi.org/10.1142/S0219720015500067

66. Campodonico MA, Vaisman D, Castro JF, Razmilic V, Mercado F, Andrews BA, Feist AM, Asenjo JA. 2016. *Acidithiobacillus ferrooxidans's* comprehensive model driven analysis of the electron transfer metabolism and synthetic strain design for biomining applications. *Metab Eng Commun* 3:84–96. https://doi.org/10.1016/j.meteno.2016.03.003

67. Valdés J, Pedroso I, Quatrini R, Dodson RJ, Tettelin H, Blake R, Eisen JA, Holmes DS. 2008. Data from "Acidithiobacillus ferrooxidans metabolism: from genome sequence to industrial applications". *GenBank*. https://www.ncbi.nlm.nih.gov/assembly/GCF_000021485.1/.

68. Quatrini R, Appia-Ayme C, Denis Y, Jedlicki E, Holmes DS, Bonnefoy V. 2009. Extending the models for iron and sulfur oxidation in the extreme acidophile *Acidithiobacillus ferrooxidans*. BMC Genomics 10:1–19. https://doi.org/10.1186/1471-2164-10-394

69. Mangold S, Valdés J, Holmes DS, Dopson M. 2011. Sulfur metabolism in the extreme acidophile *Acidithiobacillus caldus*. *Front Microbiol* 2:17. https://doi.org/10.3389/fmicb.2011.00017

70. Chen L, Ren Y, Lin J, Liu X, Pang X, Lin J, Lebedev N. 2012. *Acidithiobacillus caldus* sulfur oxidation model based on Transcriptome analysis between the wild type and sulfur oxygenase reductase defective mutant. PLOS ONE 7:e39470. https://doi.org/10.1371/journal.pone.0039470

71. Justice NB, Norman A, Brown CT, Singh A, Thomas BC, Banfield JF. 2014. Comparison of environmental and isolate Sulfobacillus genomes reveals diverse carbon, sulfur, nitrogen, and hydrogen metabolisms. BMC Genomics 15:1107. https://doi.org/10.1186/1471-2164-15-1107

72. Lieven Christian, Beber ME, Olivier BG, Bergmann FT, Ataman M, Babaei P, Bartell JA, Blank LM, Chauhan S, Correia K, Diener C, Dräger A, Ebert BE, Ediringhe JN, Faria JP, Feist AM, Fengos G, Fleming RMT, García-Jiménez B, Hatzimanikatis V, van Helvoirt W, Henry CS, Hermjakob H, Herrgård MJ, Kaafarani A, Kim HU, King Z, Klamt S, Klipp E, Koehorst JJ, König M, Lakshmanan M, Lee D-Y, Lee SY, Lee S, Lewis NE, Liu F, Ma H, Machado D, Mahadevan R, Maia P, Mardinoglu A, Medlock GL, Monk JM, Nielsen J, Nielsen LK, Nogales J, Nookaew I, Palsson BO, Papin JA, Patil KR, Poolman M, Price ND, Resendis-Antonio O, Richelle A, Rocha I, Sánchez BJ, Schaap PJ, Sheriff RSM, Shoaie S, Sonnenschein N, Teusink B, Vilaça P, Vik JO, Wodke JAH, Xavier JC, Yuan Q, Zakhartsev M, Zhang C. 2020. MEMOTE for standardized genome-scale metabolic model testing. Nat Biotechnol 38:272–276. https://doi.org/10.1038/s41587-020-0477-4

73. Lieven C, Beber ME, Olivier BG, Bergmann FT, Ataman M, Babaei P, Bartell JA, Blank LM, Chauhan S. Correia K. 2020. Supplemental material for Lieven et al. Available from: https://github.com/opencobra/memote

74. Nemati M, Harrison STL, Hansford GS, Webb C. 1998. Biological oxidation of ferrous sulfate by *Thiobacillus ferrooxidans*: a review on the kinetic aspects. *Biochem Eng J* 1:171–190. https://doi.org/10.1016/S1369-703X(98)00006-0

75. Ebrahim A, Lerman JA, Palsson BO, Hyduke DR. 2013. Cobrapy: constraints-based reconstruction and analysis for python. BMC Syst Biol 7:74. https://doi.org/10.1186/1752-0509-7-74