Intratumoral heterogeneity and metabolic cross-feeding in a three-dimensional model of breast cancer: an in silico perspective.

Intratumoral heterogeneity and metabolic cross-feeding in a three-dimensional model of breast cancer: an in silico perspective. | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Intratumoral heterogeneity and metabolic cross-feeding in a three-dimensional model of breast cancer: an in silico perspective. Osbaldo Resendis-Antonio, Jorge Arellano-Villavicencio, Aarón Vázquez-Jiménez, and 4 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4864972/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Today, intra-tumoral composition is a relevant factor associated with the progression and aggresivity of cancer. Despite it suggests a metabolic interdependence among the subpopulations inside the tumor, a detailed map of how this interdependence contributes to malignant phenotype is still lacking. To address this question, we developed a systems biology approach integrating single-cell RNASeq and genome-scale metabolic reconstruction for mapping the metabolic cross-feeding among the subpopulations previously identified in spheroids of MCF7 breast cancer. By calibrating our model with expression profiles and the experimental growth rate, we concluded that the Reverse Warburg effect emerges as a mechanism to optimize the community growth. Furthermore, through an in silico analysis, we identified lactate, alpha-ketoglutarate, and some amino acids as key metabolites whose disponibility alters the growth rate of the spheroid. Altogether, this work contributes with a strategy for assessing how space and intra-tumoral heterogeneity influence the metabolic robustness of cancer, issues that claim computational strategies to move toward the design of optimized treatments. Biological sciences/Systems biology/Computer modelling Biological sciences/Cancer/Cancer metabolism Biological sciences/Systems biology/Biochemical networks Cancer Metabolism Tumor microenvironment Heterogeneity MCF-7 Reverse Warburg Effect Lactate Systems Biology. Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 I - Introduction The Warburg effect, described by Otto Warburg in 1956, has been one of the most critical metabolic markers in the study of cancer, where tumor cells oxidize high rates of glucose as an energy source to its biotransformation into lactate even in the presence of oxygen 1 , acidifying the tumor microenvironment, favoring cell transformation and invasion into other tissues 2 . However, some tumors, such as in breast cancer, adopt different strategies to obtain energy, such as the incorporation of glutamine 3 to adapt to alterations in their microenvironment and compensate for their energy demand or oxidizing fatty acids whose metabolism is closely related to survival and resistance to chemotherapeutic treatment 4 . In recent years, it has been considered that a metabolic association between the tumor and adjacent cells may ensure tumor survival. This metabolic coupling has been experimentally validated in cancer-associated fibroblasts (CAF) and cancer cells 5,6 . This relationship has been hypothesized to be induced by reactive oxygen species (ROS) or acidification of the tumor microenvironment. Among its consequences, the incentive changes in the microenvironment as a consequence of tumor growth induce an elevated glycolytic activity in CAFs. The behavior benefits tumor cells that could dispose of the lactate released into the extracellular medium in response to unfavorable conditions such as hypoxia or changes in pH. This association has been described as a reverse Warburg effect, a feature accepted as a critical metabolic adaptation in the progression of pathology 5,6 . The transient nature of this association and the lack of information on the conditions and metabolic regulatory pathways involved in its establishment add significant complexity to its study. The complexity of studying cancer is attributed to the inter-tumoral (between tumors) and intra-tumoral (within the tumor) heterogeneity, which makes it challenging to classify possible metabolic phenotypes or subphenotypes between different types of cancer and cells of the same tumor 7 . The progress of omics technologies, such as single-cell transcriptome sequencing, has opened a new window to understanding tumor phenotypes and heterogeneity from gene expression data 8 . Accordingly, the implementation of single-cell RNAseq techniques (scRNA-Seq) in an in vitro 3D culture of luminal invasive breast cancer subtype A (MCF-7) and a functional enrichment analysis determined the presence of three cell subpopulations coexisting within the same spheroid with characteristics of proliferation, invasion, and a third phenotype of transition or metabolic reservoir 9 . Some questions arise based solely on the expression data, including the metabolic interdependence between these coexisting subpopulations and how we can modulate metabolism to interfere with tumor cell progression and survival. Enrichingic pathways from sequencing data have significantly contributed to the characterization of specific functions in tumor cells and non-tumor tissue, providing relevant information for prognosis and the development of targeted therapies 10 . However, this information does not determine whether functional heterogeneity is supported by metabolism or how cells may couple to promote tumor survival. The cooperative metabolic landscape that may be present inside a solid tumor and the complexity attributed to its cellular heterogeneity opens the door to implementing comprehensive strategies where mathematical formalisms and bioinformatics tools could provide approximations with the capacity to be experimentally corroborated. The use of computational tools for generating genome-scale metabolic reconstructions (GEMS) from expression data has been suggested 11 ; this has benefited from improvements in single-cell and massive sequencing techniques to obtain better-calibrated tissue-specific reconstructions with the ability to correlate with experimentally measured metabolites 12 , 13 , 14 . Applications in the area of systems biology and metabolic engineering have been handy in simulating the specific production of a metabolite or stimulating cell growth by implementing mathematical formalisms modeling metabolism assuming a steady-state system; this strategy does not consider the quantification of metabolites in the medium or related enzyme kinetic information 15 . The metabolism modeling through flux balance analysis (FBA) has had great validity when implemented to successfully predict the growth rates of some organisms or the maximization of the production of a metabolite of interest, generating a certain degree of validation for metabolic research 16 . Until now, FBA-based models for the study of cancer have been limited to the representation of a single specific biological trait and the modeling of a single cell or tumor. However, in recent years, progress has been made in developing computational models that use GEMs to study communities 17 . Tools such as MICOM (microbial community) have been implemented and validated in the modeling of complex communities, allowing the analysis of the interactions between individuals and extracellular microenvironments with greater precision 18 . Despite this approach being developed for microbial communities, it is not limited, given that it considers the metabolic pathways without being species-specific. Broadening the use of these tools to study cancer offers a complementary approach, as they allow the easy modification of specific conditions, such as the presence of nutrients in a medium and alterations in the presence of genes or specific metabolic reactions. This ability to efficiently generate hypotheses facilitates the identification of therapeutic targets and the creation of approaches that seek to respond to concrete phenomena. This project proposes a strategy to model in silico metabolism from MCF7 spheroids obtained from scRNA-seq data. Algorithms were implemented using the expression profile to identify the presence of enzymes and metabolite transporters and to create tissue-specific metabolic reconstructions of the cell community, constituted by invasive, proliferative, and reservoir cells. The community modeling reflected the metabolic cooperation between the subpopulations by providing the exchange fluxes between each subpopulation, evidencing the presence of the intratumoral reverse Warburg effect and the implementation of alternative carbon sources by the tumor cells in response to the oxygen gradients present due to tumor growth. II - Results 2.1 Community metabolic modeling strategy design. To model the metabolism of a cell community, we developed a strategy that involved creating GEMs based on the scRNA-seq data from a spheroid culture of MCF-7. This approach enabled us to generate three reconstructions representing distinct cell subpopulations: Invasive, Reservoir, and Proliferating (refer to Fig 1a). Next, to address the challenge of constructing a community model from these metabolic reconstructions, we proposed the implementation of tools capable of integrating the reconstructions by specifying a growth medium and the relative proportions of each subpopulation. This setup was designed to simulate two temporal stages of cell spheroid growth (days 6 and 19) to elucidate the impact of cell distribution, subpopulation proportions on community growth, and the exchange fluxes with the microenvironment, as illustrated in Fig 1b. 2.2 Specific Sub-population genome-scale metabolic reconstruction from scRNAseq To survey the metabolic cross-feeding among the cancer subpopulation in MCF7, we accomplished a genome-scale metabolic reconstruction starting from the scRNA-Seq data associated with the three identifiable subpopulations at two times of the growth process: 6 and 19 days of cultured 9 . As originally discovered by the authors, these subpopulations carry out different biological functionalities (proliferative, invasive, and reservoir), and their interactive crosstalk among them presumably is an essential factor that stabilizes the spheroid. We preserved the original annotation made by the authors. Before proceeding with the metabolic reconstructions, an imputation method for scRNA-seq data was applied to the count matrix to denoise the data and correct the 'dropout events.' This minimizes the noise characteristic of scRNA-seq technology 19 . Overall, the imputation contributed to improving the subsequent metabolic reconstruction procedure. After imputing the count matrix, we proceeded with the reconstruction by applying the Cost Optimization Reaction Dependency Assessment (CORDA) method to the gene expression profiles of each subpopulation. Despite the variety of techniques to accomplish the reconstruction, we selected CORDA due to its ability to create non-minimalistic metabolic models with greater fidelity in representing biological characteristics 9,13,20 . As a result, we created three metabolic reconstructions, one for each cell subpopulation described inside the spheroid. As expected, each reconstruction has distinctive metabolic capabilities reflecting the feasible space of functionalities in each subpopulation; see Supplementary Table 2. To ensure the quality of these reconstructions and their nearness to represent actual metabolic processes, we evaluated each network through the Metabolic Modeling Test (MEMOTE) (see methods), a computational platform to assess the quality of genome-scale metabolic reconstructions through a variety of parameters such as the consistency of the stoichiometric matrix, the evaluation of biomass production, the Gene-Protein-Reaction association, among others criteria. The quality report obtained with MEMOTE indicated that the reconstructions have a subtotal consistency score of 95%, suggesting that they present a proper stoichiometry matrix, mass balance, charge balance, and adequate connectivity between metabolites in the reconstructions, see supplementary information section 4. 2.3 Optimization and Community Modeling Once we obtained and proved the quality of each reconstruction, we proceeded to simulate the metabolic capacities when coexisting as a system computationally. We tested the community for each time, at 6 and 19 days of MCF-7 culture (days that were analyzed in the experimental design of Muciño et al. 2020) and constrained by medium (see Fig 1b). This model allowed us to tackle two fundamental questions: postulate the shared metabolites between the subpopulations (cross-feeding metabolic activity) and evaluate how their communication affects metabolic stability and biomass production in the entire community. All reconstructions were integrated and modeled in MICOM 18 . MICOM considers the relative abundance of each subpopulation and the metabolites in the extracellular medium to maximize a community biomass function. Computationally, MICOM combines Flux Balance Analysis by integrating linear optimization with an L2 regularization constraint to estimate the biomass required to keep the biomass for each population (see methods - Community modeling - MICOM). To reproduce the experimental conditions, we defined the extracellular medium in terms of the metabolic composition of culture media used to grow the MCF7 spheroids in the original experiment experimentally, see supplementary Table 1. To simplify our analysis and avoid diffusion processes, we considered that all the media resources are available for each subpopulation at both simulation times (6 and 19 days of progression of the spheroid). Under these assumptions, our metabolic simulation indicated that the community growth rate is 0.026 mmol/cell/h, corresponding to a doubling time of approximately ~ 26 hr on both days; details on duplication time calculations see Supplementary Information section 8. Consistently with our in silico estimations, this double time is nearly associated with the experimental measurements reported for the MCF-7 cell line between 20 and 24 hr doubling time 21 . At this point, our model was a computational platform to create testable hypotheses around the exchanged metabolites with the medium and the metabolic cross-feeding between the subpopulations. This interaction can be represented graphically as groups corresponding to the subpopulations with their respective exchange interactions between them and with the extracellular medium, see Fig 2a. This computational platform can guide design experiments, prove or generate biological hypotheses, and explore strategies for searching metabolic targets for potential therapeutic interventions. 2.4 Effect of pseudo-spatial distribution of subpopulations in MCF-7 spheroid In the previous section, we simulated the metabolic interaction among the three subpopulations by neglecting the diffusion process inside the tumor spheroids. However, as revealed by spatial RNAseq, the cellular distribution and arrangement of metabolites in different tumor strata (such as oxygen and carbon sources) significantly influence the intratumoral profiles of metabolic heterogeneity 22 . To emulate this process and quantitatively assess the effect of metabolic gradients on tumor malignancy and stability, we extend our computational genome-scale metabolic modeling by including layers of gradients over oxygen, a metabolite whose concentrations inside the tumor depend exclusively on the size of the spheroid. In the rest of this work, we will assume that oxygen is a metabolite consumed solely by the environment, a key metabolite in establishing intratumoral cellular heterogeneity. We selected three spatial layers of oxygen concentration in the inner microenvironment to emulate the oxygen gradient. The simulation starts from the outside to the center of the spheroid; the first level emulates cell normoxia, the second level represents an intermediary level of oxygen abundance, and the last level indicates a hypoxia condition. Given the uncertainty around the precise location of the cell subpopulations in the spheroid, we modeled the metabolic interaction over six different spatial configurations, each defined by different spatial arrays of subpopulation and oxygen availability (see Fig 2b). In principle, all these configurations are feasible, at least in an in silico model; however, we focused on those configurations that maximize the growth rate of all the communities to survey the possible metabolic communication between the subpopulations. Hence, we calculated the community growth rate for each spatial distribution, considering the abundance of the subpopulations previously reported at 6 or 19 days. In Fig. 2c and 2d, we can observe the community's growth rate; the results show the cellular organization that the spheroid could have in both days of study when taking into account the best growth rate of the community. Consequently, the spatial configuration denoted by condition E in Fig 2c is the ideal combination to optimize the doubling time of the spheroid at day 6. According to intuition, this configuration suggests that the invasive subpopulation is located in the center of the spheroid, the reservoir at a medium level, and the proliferative subpopulation on the outside layers. This optimized distribution suggests that the invasive cells are better adapted in the center of the spheroids, a place with relatively low uptake of oxygen. In contrast, when changing the corresponding population abundances for day 19 of growth, the spatial configuration E (see Fig 2d) proved the least optimal. Our simulations suggested that placing the invasive cells on the outside enhances community growth (see Fig 2d condition A ). In summary, we conclude that configurations E6 and A19 supply metabolic advantages to optimize the community growth rate at a systemic level. To analyze the metabolic phenotypes that emerged from these configurations, the following sections will: 1) evaluate the metabolic exchange between the environment and the tumor, as well as the metabolic cross-feeding among populations, and 2) identify the metabolic pathways involved in each subpopulation, based on enzymatic reactions classified through the Virtual Metabolic Human (VMH) database 23 . 2.5 Metabolic mapping of cancer spheroids: secretion, consumption, and cooperativity among sub-populations. Having simulated the growth rate at both times of the spheroid, we explored how the metabolic exchange between the tumor and the extracellular milieu is influenced by the cellular organization of the spheroid at both times of spheroid progression (see methods section) To simplify data presentation, we have generated a Boolean representation of metabolic excretion or consumption between the subpopulations and the environment. The metabolites can be classified into three categories: those secreted and produced simultaneously and alternate between production and excretion as time goes on. For example, we noted that most of the amino acids and carbon sources, such as galactose and glutamine, are consumed by each subpopulation from the media at both times (see Fig 3a). Meanwhile, our modeling suggests that metabolites such as carbon dioxide, ammonia, and succinate are excreted from the community to the media. Furthermore, we identified a third class of metabolites whose production consumption changed as time varied. This last category includes 2-oxoglutarate (AKG), acetoacetate, fumarate, L-lactate, and amino acids such as glycine and L-alanine. In terms of the metabolic cross-feeding among subpopulations, results showed that all three subpopulations have the metabolites present in the medium; however, some metabolites were not initially found in the culture medium and whose origin is a metabolic by-product of one of the subpopulations. Figure 3a indicates all the metabolites produced by one subpopulation and consumed by another. Among these metabolites, we identified primarily essential amino acids like L-isoleucine, L-leucine, L-lysine, L-methionine, L-phenylalanine, L-threonine, L-tryptophan, L-valine, non-essential like L-asparagine, L-histidine, L-serine, conditional non-essential like L-arginine, L-cysteine, L-glutamine, L-tyrosine and carbohydrates such as galactose. Our model postulated that these metabolites maintain constant consumption regardless of the growth day or subpopulation type. On the other hand, some amino acids do not follow this behavior. For example, glycine, a metabolite consumed by the invasive population, can be taken from the medium and reservoir cells (see Fig 3a). However, proliferative subpopulations can also contribute to this amino acid on day 6. On the other side, we observed that specific metabolites are produced exclusively by one type of subpopulation. For instance, the addition of alanine to the medium is regulated by the invasive cells on both days of the study, similar to proline by the proliferative cells. As Fig 3b shows, our analysis postulates an intricate metabolic mapping underlying an interdependence between the cancer subpopulations in MCF7 spheroids. The production and consumption of metabolites are crosslinking to meet the metabolic requirements of the growing community. Fig 3b shows the origin of the metabolite from its respective subpopulation and where it is required. Such is the case of alpha-ketoglutarate, fumarate, and lactate, whose consumption and secretion profile changes at both times of the spheroid. The panels of Figure 3b and 3c show the cross-feeding metabolic profile of the community on days 6 and 19. On day six, the proliferative cells secrete alpha-ketoglutarate and lactate, a characteristic marker of the Warburg effect. In contrast, the invasive cells dispose of this metabolite by incorporating it into their pathways. This effect is also represented on day 19 and can be seen as a hallmark between proliferative and invasive subpopulations. Altogether, the cooperative metabolic profile between the three subpopulations creates a complex scenario of metabolic inner/outer communication at different stages of growth. Among these possibilities, some metabolites behave depending on the day of growth. For instance, AKG is secreted to the medium on day six by proliferating cells. On day 19, the invasive cells show this behavior, exchanging roles between excretion and consumption only between these two subpopulations. By analyzing a portion of the top 20 of the approximately ~350 exchange reactions and ~2000 community reactions, we may be ruling out the possibility that AKG-like switches are also present with metabolites whose exchange fluxes are minimal. These scenarios could add a touch of complexity to the design of metabolic targets due to the metabolic adaptation and reconfiguration of a spheroid over time. We estimated the metabolic pathways with the highest activity in each subpopulation to assess metabolic adaptations. Then, based on the reconstruction of metabolic pathways, reactions with non-zero flux units (in any directionality) were isolated to determine each pathway's importance in each subpopulation. A common feature is that all three subpopulations have activity in their main pathways of central metabolism, such as glycolysis, the TCA, and the pentose phosphate pathway (see supplementary Fig 1). However, there were cases, such as lipid metabolism (Fatty acid oxidation and Fatty acid synthesis), where a higher count of positive-flux reactions was observed in proliferative and reservoir cells than in invasive cells. On day 19, the invasive cells had more active reactions due to changes in spheroid abundance and size (see supplementary Fig 1). Despite this analysis supplying a proxy of the metabolic activity associated with each subpopulation, the reactions present within each pathway usually could not be the same, or if so, the behavior of the directionality of the metabolic fluxes may vary in those reactions with reversible capacities. The description of directionality changes in reactions provided a more comprehensive approach to tracking how each subpopulation uptakes, transforms, and utilizes the metabolites in the medium for energy generation and biomass production. 2.6 Metabolic activity inside the populations Once the central metabolic pathways were identified, they were used to analyze individual behavior and determine the existence of metabolic cross-feeding between subpopulations. Among these, we identified energy pathways such as glycolysis, the Krebs cycle, oxidative phosphorylation, the pentose phosphate pathway, and some other metabolic pathways that might be of interest (see Fig. 4 and Supplementary Fig 2). The simulation is not only limited to providing fluxes of exchange within the medium but also fluxes of intracellular reactions of each subpopulation. This information allowed us to generate a descriptive picture of the behavior of each cell type—for example, the glycolytic activity of the community. By having galactose as the main carbohydrate according to the original experimental design, we consequently obtained a high activity of the enzyme galactosidase in most of the subpopulations on both days of study except for the invasive cells that show lower fluxes of activity in this enzyme on day 6 of growth (see Fig 4b)—suggesting that this metabolite is not as essential for the subpopulation in early periods of growth. The information obtained by simulation suggests that the invasive subpopulation does have glycolytic activity, as seen in Fig 4a, as metabolic fluxes are seen in most of the enzymes involved in glycolysis. However, particular details such as lactate dehydrogenase (LDH) enzyme activity do not reflect a directionality similar to the proliferative cells that initially convert pyruvate to lactate but instead that the invasive cells are in a reversal process of this action and not by carrying out the oxidation of galactose to pyruvate by a classical pathway. Reaction directionality reversal behaviors further solidify the association theory between subpopulations by incorporating metabolites secreted into the medium by neighboring cells. The origin of metabolites by alternative pathways, such as anaplerotic pathways, often converges to the TCA for generating reducing power and oxidative phosphorylation. Given the results shown in Fig 4c, we have concluded that full participation of TCA is not necessary to satisfy the energy demand of the cells. Although all reconstructions have full TCA, the optimization suggests a preference for certain enzymes that compensate for the objective function of the simulation. For example, given the directionality of the reaction, the model suggests that succinate coenzyme A ligase, one of the critical enzymes of the TCA, consumes succinate to produce succinate and GTP from succinyl-CoA. These reactions could be preceded by the activity shown by alpha-ketoglutarate dehydrogenase, which provides the Succinyl-CoA originating from the presence of alpha-ketoglutarate (see Fig 4c). Both enzymes are shared among the three subpopulations. However, particular cases, such as invasive cells on day 19 and proliferating cells on day 6, activate another section of TCA, such as succinate dehydrogenase, whose reversibility converts fumarate that might be taken up from the extracellular medium when secreted by reservoir cells to succinate (see Fig 4c). The partial recovery of pathways, such as TCA, suggests that the availability of metabolites such as pyruvate for energy may have different origins. From galactose oxidation or incorporation of lactate from the microenvironment to the involvement of anaplerotic pathways such as alanine incorporation, of which L-alanine transaminase activity was recorded in the invasive and reservoir subpopulations (see Fig 4d). The simulation gave us a broader picture of what alternative energy pathways the cells might use. The simulation also facilitates the analysis of alternative anaplerotic pathways. To identify alternative pathways, we separated already described features of tumor cell metabolism, such as the incorporation and activity of enzymes like mitochondrial glutaminase and glutamate dehydrogenase enzymes, which coordinate for generating alpha-ketoglutarate in TCA, see Fig 4d. This set of results allowed us to characterize each subpopulation and assign specific metabolic characteristics. According to Muciño et al., 2024, the cells have a glycolytic activity. However, our results complement this assertion by providing the variations of each subpopulation to obtain energy. These differences could be explained by the capacity of each cell to adapt its metabolism to the conditions of its microenvironment and the cooperativity between individuals in the community. Under this criteria, we explored the adaptive capacities of the community in response to stress in its microenvironment. Directly, oxygen was taken as an essential element in the community that could quickly stimulate a change in the subpopulations. As we argued in the next section, the action of different oxygen zones makes evident the importance of some amino acids for the survival of the whole community. 2.7 Metabolic response of cancerous subpopulations under gradients of oxygen. The stress induced in the community through the oxygen supply of each subpopulation could be an essential factor for the appearance of heterogeneous profiles in the spheroid; this is because one of the variables that allowed us to establish cooperativity among the subpopulations was the presence of a simulated gradient in the community model. To determine the response and capabilities of the community to an adverse scenario, he designed a modeling approach that could represent the decrease in oxygen availability from normoxia to anoxia. Each simulation point reflects the behavior of the subpopulations in the face of oxygen deprivation (see Fig 5). The response was reflected in a greater or lesser preference for selected metabolites as possible energy sources, such as galactose, pyruvate, L-lactate, and AKG, as well as amino acids such as glycine, L-alanine, L-glutamine and L-serine. An example is glutamine consumption, an amino acid associated with the fundamentals of tumor metabolism whose role has been essential as an alternative carbon source for tumor cells to ensure survival under unfavorable conditions 3 . This metabolite proved crucial for the three subpopulations on both study days. The simulation showed that the proliferating subpopulation had the most glutamine in the medium on day six. As oxygen availability decreases, the invading cells begin to increase their consumption. The opposite is the case on day 19, where the invading cells now have most of the glutamine, and this behavior does not appear to be affected in response to the oxygen gradient (see Fig 5a and 5b) . The cells can incorporate metabolites in the medium, such as pyruvate, to obtain energy under normoxic conditions or in response to gradients. Proliferating cells dispose of pyruvate on day 6 (see supplementary Fig 3a), and as oxygen decreases, this behavior is maintained, as well as in the reservoirs, but in a lower proportion. By day 19, under normoxic conditions, the entire community had a low amount of this metabolite, which could be interpreted as a low preference or essentiality. However, as disposition decreases, the invasive subpopulation begins to mobilize pyruvate into its interior (see supplementary Fig 3b). The effects differed significantly from those observed when consuming the primary carbon source, galactose. An almost linear gradient accompanied the decrease in oxygen on both days (see Supplementary Fig 3c and 3d) as a response to oxygen availability and possibly by incorporating alternative sources. Another type of response observed in the simulation was the cooperativity between subpopulations; in normoxic conditions, there was already an association between proliferating and invasive cells with lactate mobility. However, on day six, the decrease in oxygen increased the transport of lactate into the invasive cells from the proliferating cells and the reservoir. Stages close to zero activated lactate production in the reservoir contribute to the lactate demand in the community (see Fig 5c). This behavior is the first finding to resemble the reverse Warburg effect as an adaptive response to gradients. The simulation showed how the three subpopulations may respond to changes in their microenvironment. This feature was maintained on day 19 (see Fig 5d); however, the exchange interactions diminished as oxygenate approached zero. A further metabolite that showed a peculiar behavior in response to gradients was AKG, which had already been observed to interact between spheroid subpopulations under normoxic conditions, originating in proliferating cells towards invasive cells. However, when hypoxia is approached, the community adapts in such a way that the directionality of this metabolite is inverted, generating the hypothesis that invasive cells have a better capacity to adapt to gradients and provide AKG to proliferating cells to compensate for tumor survival (see Fig 5e). In the same way as with lactate, by day 19, the spheroid organization established an exchange of lactate from invasive to proliferative from normoxia, and oxygen deprivation creates a decay in the interactions, see Fig 5f. This behavior suggests that metabolic cross-feeding between proliferative and invasive cells is established in the first days of growth; day 6 shows a cellular organization that responds better to oxygen deprivation stress than day 19 since more invasive cells do not have cellular allies to help their survival. 2.8 Knock-Out in silico of Metabolic Enzymes To identify those reactions that could be essential for community growth, a Knock-Out (K.O.) essentiality analysis of metabolic reactions was implemented, including exchange reactions with the medium and individual intracellular reactions (see methods section), using growth rate values to determine the impact of K.O., a list of reactions whose absence directly affected the entire community was isolated. This suggested that there is a dependence on metabolites from the medium, especially amino acids, as the most important for survival (see Fig 4b). The relationship was maintained on both days of growth, with the only difference being that Serine might be essential for the day 19 spheroid. To evaluate if the metabolic dependence of the amino acids on the spheroid is robust in space, we repeated the simulations by mimicking the inner zone of the spheroids through the gradient of oxygen described in the previous section. Hence, we separately accomplished the KO analysis by considering three feasible scenarios of oxygen: normoxia, mid-level, and hypoxia. The essentially analysis obtained from each oxygen environment did not show significant changes between days 6 and 19. These results postulated that in this time scale, the metabolic phenotype of the entire community mainly depends on the amino acids and galactose as a carbohydrate present in the medium located in the external environment (see Fig 5). Despite the diversity of components in the extracellular medium, such as pyruvate and some vitamins, the presence or absence of these components did not show importance in sustaining growth. Furthermore, no intracellular reactions appeared as essential in any of the three sub-populations, suggesting that the environment is the most probable metabolic mechanism to alter its metabolic tumor phenotype. Our in silico analysis highlights the importance of external amino acids for tumor growth. Despite some metabolites being essential to growth at any time and oxygen concentration (such as D-galactose, L-glutamine, and L-histidine), others have an interment relevance in the community growth rate, see Figure 5g. For instance, this is the case for L-valine, L-serine, and L-choline, whose relevance to growth rate depends on time and oxygen availability. Overall, these results open an avenue toward an in silico platform to construct testable experimental hypotheses in the design of therapeutic strategies, mainly focused on controlling the growth rate for 3D tumor spheroids. III - Discussion In vitro, 3D tumor spheroids are a perfect system to elucidate how intratumoral heterogeneity shapes the metabolic organization and sustains the malignant phenotype. As we have described here, our in silico analysis, which integrates single-cell RNASeq on MCF7, is a powerful tool. It allows us to map the complex cross-feeding between the cellular components of the MCF-7 spheroid and predict how metabolic exchange profiles change as a function of the relative cell abundances and time. This predictive power instills confidence in the model's ability to anticipate metabolic changes. Among the most remarkable metabolic exchanges that could be identified with the computational model was the secretion and consumption of lactate. This metabolite is used as an energy source by cells in oxygen-restricted conditions, such as invasive cells inside the spheroid. From a diagnostic point of view, lactate has been listed as a biomarker of poor prognosis in cancer patients. This oncometabolite has been positioned as an essential marker due to its involvement as a tumor promoter and facilitator of immune evasion 24 . Our results broaden the picture of lactate involvement beyond the tumor microenvironment regulation and as a potential energy source inside the inner sub-populations. In recent years, the reverse Warburg effect has increased relevance by postulating the importance of lactate as an additional carbon source to produce biomass and not as a product, as described in the classical Warburg effect when cancer cells are near other neighbor cells. For instance, this is the case of the association between Cancer Associated Fibroblast (CAFs) and tumor cells by ensuring cancer survival through an exchange of lactate as an additional carbon source 6 . Although the association has been experimentally induced between the tumor and its stroma, we propose a similar scenario among cells of the same tumor (among the subpopulations) as a possible metabolic adaptation. In agreement with some in vitro experiments, promoting lactic acidosis in the microenvironment by tumor cells can induce metabolic turnover in adjacent cells by promoting the shift from aerobic glycolysis to an oxidative phenotype 25 . The conditions simulated for invasive cells on day six could be equivalent to an experiment in which low oxygen levels and acidic conditions favor lactate consumption in breast cancer cell lines 26 . This behavior has also been observed in other tumor lines, such as lung cancer 26,27 . These results allow us to question whether this behavior may be a hallmark of glycolytic tumors or a cellular adaptation in response to stimuli from their microenvironment. The organization proposed by the simulation for the spheroid at day 6 is valid when considering that proliferative cells condition the acidic microenvironment, and low oxygen levels promote the transformation of cells in the center of the spheroid to an invasive phenotype 28 . The metabolic connection among the two subpopulations creates a principle of intratumoral cooperativity to ensure the entire community's survival. In the simulations of gradual oxygen deprivation in the microenvironment, the response of the invasive cells was to increase the rate of lactate consumption, compensated by the proliferative cells to the point that the reservoir subpopulation began to participate by secreting lactate to the medium. See Fig 5c. This indicates that this exchange is enhanced in response to hypoxia. On the other hand, the proposal on day 19 suggests that this association is already established, and oxygen gradients decrease the consumption rate of metabolites such as lactate, see Fig 5d. This may be related to the relative abundance of each subpopulation on day 19 because, on this day, there are more invasive cells than proliferative cells, which decreases the availability of lactate in the medium; the incorporation of pyruvate and high rates of glutamine consumption compensates this demand. All this suggests that the computational model of day six growth would have greater validity for future analyses where the interest is to evaluate the community's response to oxygen deprivation. On the other hand, the day 19 proposal contemplates advanced tumor scenarios, where approximately 50% of the community is invasive, and its metabolism is based on the consumption of L-glutamine (Fig 5b), galactose (supplementary Fig. 3d), serine (supplementary Fig. 3f) and glycine (supplementary Fig 3h). Nevertheless, being a computational simulation, these tools allowed us to monitor alternative metabolic pathways that could be interesting. Such as the dynamics of AKG exchange from proliferative cells to invasive cells on day 6. The role of AKG in cancer metabolism has been somewhat controversial. Meanwhile, its potential has been observed as an enhancer of radiotherapies by sensitizing triple-negative breast cancer tumor cells 29 , antagonizing tumor progression in P53-deficient malignancies 30 , and as an ammonium carrier through glutamate synthesis via reversible glutamate dehydrogenase, reducing urea production 31 , however, in tumors such as colon cancer AKG is associated with carcinogenesis 32 . Although glutamine is a semi-essential amino acid for tumor cells, the AKG skeleton is more soluble than glutamine, and its metabolic participation is usually similar, compensating for the metabolic demands of invasive and proliferative cells. AKG can be directly incorporated into the TCA as Succinyl-CoA through the action of alpha-ketoglutarate dehydrogenase to obtain energy via OXPHO 33 , which would correspond to their participation within the metabolic synergy between proliferative and invasive as an alternative source of energy or as metabolic compensation since having a metabolite available in its most accessible form would be energetically more efficient. The ability of AKG to diffuse freely through voltage-dependent transporters facilitates its permeability and, through oxoglutarate/malate transporters, without difficulty, ensures its bioavailability in mitochondria for incorporation into TCA 34 . Interestingly, a survey of the transport fluxes of AKG shows that its origin is from proliferative to invasive cells. However, the effect of hypoxia created a switch in this behavior, reversing the roles. The proposal is that in normoxic conditions, the proliferative cells provide AKG and lactate as metabolic sustenance to the invasive cells; this establishes a picture of metabolic compensation since this stage approximates the standard growth conditions of a spheroid. Fig. 6 reflects the proposed interaction between two cell subpopulations with the most significant contribution to the sustenance of the community, which points to the presence of the Reverse Warburg effect in cells of the same tumor type and a connection between the AKG that compensates for the energy rate of the cells in the center of the spheroid. The reversal of the dynamics, on the other hand, is associated with an increase in the rate of lactate production. Since it is maintained as a response to hypoxia and observed as a feature on day 19, this condition could be related to a marker of poor prognosis. Further study is crucial to elucidate the pathways involved in such associations and to explore how they can be intervened by selecting metabolic targets. The simulation suggests a connection between the reservoir and proliferative cells by exchanging fumarate and L-proline and providing glycine and Adenosine Triphosphate (ATP) in the medium available for the invading cells, see Fig 6. Despite the metabolic importance of these metabolites, the third subpopulation plays a crucial role in regulating tumor survival. The constant secretion of ATP, specifically in breast cancer, has been described to be involved in invasion processes to other tissues where ATP in the microenvironment could induce hypoxia-inducible factor (HIF1a), triggering epithelial-mesenchymal transition processes even under normoxic conditions 35 . This makes this subpopulation a possible key piece in promoting invasive processes. On another subject, in tumor processes, the deficiency of the activity of enzymes such as fumarate hydratase (FH), responsible for the conversion of fumarate to malate, produces an intracellular accumulation of fumarate and in cases such as renal cancer, the accumulation of this metabolite inhibits the activity of prolyl hydroxylases, which leads to an accumulation of HIF1a 36 . The importance of the exchange of this metabolite lies beyond the regulation of HIF1a, beyond favoring glycolytic activity over oxidative phosphorylation and regulating the sensitization of reactive oxygen species 37 by orchestrating the antioxidant response via NERF2 when fumarate inhibits its regulator 38 . In summary, the main objective of this project was to develop a strategy that could model the metabolic cross-talk among the MCF-7 cell communities. The modeling implementation tests have also allowed us to deduce and create a proposal for interaction between subpopulations of a spheroid. The proposal designed for the association of the three subpopulations corresponds to the growth period of day 6. In summary, the results suggested the existence of metabolic cooperativity between proliferators and invaders and the possibility that the third reservoir subpopulation fulfills the functions of providing energy in the form of ATP or fumarate. However, the reservoir subpopulation could be preparing the necessary conditions to promote the transformation of invasive cells towards metastatic processes, and the over-stimulation of proliferative cells with fumarate could be an initial marker of transformation toward invasive phenotypes. Finally, the approximations performed with computational simulation are created from the information obtained from sequencing data. However, proteomic and metabolomic data are needed to ensure better results and calibrate the reconstructions. This provides the presence or absence of metabolites that could generate false positives. The proposed methodology to perform these in silico analyses was preliminarily validated by bibliographic information that could qualitatively coincide with what is already known about breast cancer metabolism. Lastly, the in silico study of the inner metabolic phenotype in cancer provides a framework for designing testable experimental hypotheses in oncological studies. With the advent of single-cell technologies (such as sc-RNASeq and single-cell metabolome), these computational methods will be essential to integrate data and generate hypotheses. Considering cancer as a heterogeneous and communal structure, we envision that these methods will be necessary to move toward the design of treatments at a personalized level. Overall, these efforts allowed us to conclude that cancer is more than one phenotype; instead, the presence of subpopulations contributes with an additional degree of complexity whose inclusion should be addressed to improve the wellness of the patients. IV - Conclusions The information offered by computational models has had a growing acceptance in recent years, demonstrating the ability to represent complex biological phenomena such as cancer metabolism and its adaptations to variations in its microenvironment. Our model was able to include tissue-specific metabolic reconstructions from transcriptomic data of MCF-7 spheroids and create a community model that allows glimpses of possible pathways of metabolic interconnections between cell subpopulations, as well as the presence of the Reverse Warburg effect and other alternate pathways for energetic survival. These types of models offer approaches that could become comparable, being a great tool to suggest more targeted designs in the case of the study of cancer metabolism. The calibration of these proposals could become more accurate by having more specific metabolic targets to identify the cornerstone of cancer metabolic regulation. V - Methods 5.1 Data The data for the study was retrieved from the Gene Expression Omnibus (GEO) under the accession number GSE145633 9 . The dataset consisted of two temporal points, day six and day 19, for MCF7 spheroids cultivated in L-15 medium (ATCC 30-2008). The researchers associated the three identified clusters with their respective sample labels: Cluster A was linked to the invasive phenotype, Cluster B to the reservoir, and Cluster C to the proliferative phenotype; the cell labeling was conserved, and we took the author’s classification. 5.2 Data Imputation - sc-PHENIX To recover some of the hidden gene expression due to the noisy nature of single-cell technology, we applied an imputation specialized for scRNA-seq data. We imputed the data using sc-PHENIX (single-cell phenotype recovery via non-linear expression imputation), a recent pipeline with classification performance capabilities currently reported in the literature 19 . The input matrix consisted of 363 cells and 23,923 genes. The sc-PHENIX method uses the count matrix and a UMAP space 9 (https://github.com/resendislab/sc-PHENIX). Imputation was done to recover the missing gene expression from scRNA-seq and to improve differential gene expression. We used the same UMAP space work from the original work of dataset 9 . The parameters for sc-PHENIX were t=5, decay=5, metric='euclidean', and knn=20. 5.3 Reconstruction of Tissue-Specific Metabolic Network - CORDA Based on the imputed data for each subpopulation, we constructed a genome-scale metabolic model (GEM) using the Cost Optimization Reaction Dependency Assessment algorithm (CORDA) 13 , identified at 6 and 19 days of progression. Using the single-cell RNAseq data from an MCF7 culture, we set up the metabolic activation for every cluster that triggers diverse cellular processes. To start the reconstruction, CORDA requires as input a confidence score for each reaction of the base model. We implemented an agglomerative cluster to get the confidence level for each gene in the data, fixing the number of clusters to 3. N-dimensional values were mapped considering samples in each cluster, and the Euclidean distance was computed to associate genes with similar expressions. The confidence value range was from 1 to 3, relating low, medium, and high confidence scores. Additionally, we set the maximum confidence level to genes differentially overexpressed and the minimum confidence value for genes differentially underexpressed in each cluster; we considered differentially expressed genes for those with |fold change| ≥ 2 and p-value ≤ 0.05. The base model for each genome-scale metabolic reconstruction was taken from Recon 2.2; It integrates 7785 reactions, 5324 metabolites, and 1675 genes 39 . The non-mapped reactions from data were assigned with a confidence value of 0 (unknown confidence). We restricted the range to [-1000,1000] to be closer to the biological phenomena. Finally, given the experimental setup 9 , we modified the upper and lower bound for the exchange reaction relating to the growth medium (Supplementary Table 1). To ensure an exchange of metabolites from the medium to the cytoplasm, we set the confidence level for the medium exchange reaction to 3. All codes were implemented in Python 3, which can be reproduced at https://github.com/resendislab/corda. 5.4 Analysis of Metabolic Reconstruction Quality To ensure the reproducibility of results in creating and utilizing genome-level metabolic reconstructions, we utilized MEMOTE (Metabolic Model Testing) 40 . MEMOTE serves as a tool for evaluating the quality of these reconstructions, offering a comprehensive report and an overall assessment. It scrutinizes aspects such as reaction stoichiometry, alignment with databases like BIGG or KEGG, and incorporating a biomass function. In this particular study, all reconstructions evaluated obtained a subtotal score of 94% due to suitable stoichiometry and reaction balance. However, the report shows a 68% total score due to inconsistencies in the IDs with which the MEMOTE database makes the comparisons. However, the simulation of communities with the implemented tool is fine for modeling purposes. 5.5 Community Modeling - MICOM The community modeling was conducted by implementing MICOM (Microbial Community) version 0.34.1, a software tool enabling the creation of cellular community models and the analysis of their metabolic dynamics. In this way, community models were generated from models in SBML (Systems Biology Markup Language) format with the capability to investigate exchange fluxes among different cellular populations 18 To perform community modeling on days 6 and 19, various factors were taken into account, including the growth medium specified in Supplementary Table 1 and Table 2, the three GEMs generated using the CORDA algorithm, and the relative abundance of each cellular subpopulation. On day 6, the relative abundances were as follows: Invasive - 0.0977, Reservoir - 0.3910, Proliferative - 0.5113. On day 19, the relative abundances were: Invasive - 0.5739, Reservoir - 0.2739, Proliferative - 0.1522. 5.6 Essentiality analysis by Knock-out of metabolic reactions This analysis integrates the community generated through MICOM and the functions provided by the COBRA toolbox 41 . In this process, a total count of the reactions in the entire community was performed, and a reaction-by-reaction Knock-out (K.O.) evaluation was carried out. Optimization focuses on the objective function of the community, which is related to community growth. The community in question consists of a specific number of reactions, and a growth value is assigned to each reaction in response to the alteration of these reactions. Those reactions that, when inhibited in silico, interrupt spheroid growth when the biomass value reaches zero or decreases below its optimal value for growth are considered essential. 5.7 Oxygen gradients: Transport of amino acids and other carbon sources In the metabolic reconstructions, individual metabolic reactions have their fluxes restricted by defined upper and lower ranges, which correspond to the directionality of the reactions. To simulate the response of the community to global changes in oxygen availability, the lower limit of the oxygen transport reaction was gradually adjusted, optimizing community growth in each scenario. This adjustment was carried out gradually, decreasing from 0.16 (normoxic conditions) to 0.01 unit decrements until reaching a value of zero, which corresponds to anoxia. The objective behind this modification was to evaluate the community's behavior and determine the existence of metabolic cooperation to guarantee the whole community's survival in the face of this stress induced by the lack of oxygen. 5.8 Fluxes visualization The output generated by the community model is presented in CSV format, which contains a list of reactions along with the fluxes associated with each subpopulation. Total fluxes resulting from the simulation for day 6 and total fluxes for the day 19 simulation (see supplementary_file1) are included. In the visualization of reactions related to metabolite exchange with the environment, a filter was applied that identifies reactions whose characteristic corresponds to "EX_metabolite_m".This indicates that such reactions are related to population and growth medium exchange processes. Exchange reactions with the medium on day six and exchange reactions on day 19 (see supplementary_file1). To visualize the metabolic fluxes of internal reactions, we manually filtered to select reactions involved in glycolytic metabolism, tricarboxylic acids, and galactose metabolism, Fig 4, and some as pentose phosphate pathway and fatty acid oxidation (see supplementary Fig 3). According to the VMH database 23 , metabolites were classified, and reactions were grouped according to their metabolic pathway. This visualization was performed by implementing package ggalluvial from the ggplot2 library of R (version 4.2.2). All other figures were created with the Python matplotlib package (version 3.7.1); the information needed to create them can be found in detail in the GitHub link. 5.9 Availability Code and Data The data are available through Gene Expression Omnibus (GEO) with accession number GSE145633, obtained from Muciño et al. 9 . The source code and simulation results supporting the findings of this study are publicly available at (https://github.com/resendislab/Modeling_Heterogeneity_Cancer_Metabolism_MICOM) Declarations VI. Acknowledgements O.R.-A. thanks to the financial support from CONAHCyT (Grant Ciencia de Frontera 2019, FORDECYT-PRONACES/425859/2020), PAPIIT-UNAM (IA202720), and an internal grant from the National Institute of Genomic Medicine (INMEGEN, México). J. E. A-V. is a doctoral student from Programa de Doctorado en Ciencias Bioquímicas, Universidad Nacional Autónoma de México (UNAM) and received a fellowship from CONAHCYT - CVU 702604. This paper is part of the doctoral thesis and the requirements to obtain a Doctor of Science degree in Science from J.E.A-V. C.P-M is a doctoral student from Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional Autónoma de México (UNAM) and received a fellowship from CONAHCYT, México - CVU 855825. J.J.O-V's work was supported by CONAHCYT (Grant Ciencia de Frontera 2019, FORDECYT-PRONACES/425859/2020) and UNAM Posdoctoral Program DGAPA (POSDOC). We would also like to thank the Instituto Nacional de Medicina Genómica (INMEGEN), where the project was carried out as part of the student program. Finally, we would like to thank all the authors who contributed directly to the development of this project to improve it and the members of the Human Systems Biology Lab for their support and comments that improved this work over time. VII. Author Contributions O.R-A and J.E.A-V design the project. J.E.A-V performed the simulations, community modeling, quality testing and essentiality analysis. O.R-A. supervised the research. A.V-J. created the genome-scale metabolic reconstructions. C.P.M. Performed the data imputation. J.E. A-V., O.R-A., A. V-J., J.J. O-V., C. P-M., H. P-G., and A. T-P. drafted and discussed the main manuscript. All authors reviewed the manuscript. VIII. Competing Interests The authors declare no competing interests. Additional information Correspondence and requests for materials should be addressed to Osbaldo Resendis-Antonio. References Warburg, O. On the origin of cancer cells. Science 123 , 309–314 (1956). Vaupel, P., Schmidberger, H. & Mayer, A. The Warburg effect: essential part of metabolic reprogramming and central contributor to cancer progression. Int. J. Radiat. Biol. 95 , 912–919 (2019). Li, S. et al. Glutamine metabolism in breast cancer and possible therapeutic targets. Biochem. Pharmacol. 210 , 115464 (2023). Wang, T. et al. JAK/STAT3-Regulated Fatty Acid β-Oxidation Is Critical for Breast Cancer Stem Cell Self-Renewal and Chemoresistance. Cell Metab. 27 , 136–150.e5 (2018). Fu, Y. et al. The reverse Warburg effect is likely to be an Achilles’ heel of cancer that can be exploited for cancer therapy. Oncotarget 8 , 57813–57825 (2017). Liang, L. et al. ‘Reverse Warburg effect’ of cancer‑associated fibroblasts (Review). Int. J. Oncol. 60 , (2022). Roulot, A. et al. Tumoral heterogeneity of breast cancer. Ann. Biol. Clin. 74 , 653–660 (2016). Wu, S. Z. et al. A single-cell and spatially resolved atlas of human breast cancers. Nat. Genet. 53 , 1334–1347 (2021). Muciño-Olmos, E. A. et al. Unveiling functional heterogeneity in breast cancer multicellular tumor spheroids through single-cell RNA-seq. Sci. Rep. 10 , 12728 (2020). Bartoschek, M. et al. Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing. Nat. Commun. 9 , 1–13 (2018). Bordbar, A., Monk, J. M., King, Z. A. & Palsson, B. O. Constraint-based models predict metabolic and associated cellular functions. Nat. Rev. Genet. 15 , 107–120 (2014). Wang, Y., Eddy, J. A. & Price, N. D. Reconstruction of genome-scale metabolic models for 126 human tissues using mCADRE. BMC Syst. Biol. 6 , 153 (2012). Schultz, A. & Qutub, A. A. Reconstruction of Tissue-Specific Metabolic Networks Using CORDA. PLoS Comput. Biol. 12 , e1004808 (2016). Huang, Y. et al. Characterizing cancer metabolism from bulk and single-cell RNA-seq data using METAFlux. Nat. Commun. 14 , 4883 (2023). Orth, J. D., Thiele, I. & Palsson, B. Ø. What is flux balance analysis? Nat. Biotechnol. 28 , 245–248 (2010). Edwards, J. S., Ibarra, R. U. & Palsson, B. O. In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nat. Biotechnol. 19 , 125–130 (2001). Computational modeling of metabolism in microbial communities on a genome-scale. Current Opinion in Systems Biology 26 , 46–57 (2021). Diener, C., Gibbons, S. M. & Resendis-Antonio, O. MICOM: Metagenome-Scale Modeling To Infer Metabolic Interactions in the Gut Microbiota. mSystems 5 , (2020). Padron-Manrique, C. et al. Diffusion on PCA-UMAP manifold captures a well-balance of local, global, and continuum structure to denoise single-cell RNA sequencing data. bioRxiv 2022.06.09.495525 (2022) doi:10.1101/2022.06.09.495525. Angione, C. Human Systems Biology and Metabolic Modelling: A Review—From Disease Metabolism to Precision Medicine. Biomed Res. Int. 2019 , (2019). Jensen, B. L., Skouv, J., Lundholt, B. K. & Lykkesfeldt, A. E. Differential regulation of specific genes in MCF-7 and the ICI 182780-resistant cell line MCF-7/182R-6. Br. J. Cancer 79 , 386–392 (1999). Ahmed, R. et al. Single-Cell RNA Sequencing with Spatial Transcriptomics of Cancer Tissues. Int. J. Mol. Sci. 23 , (2022). Brunk, E. et al. Recon3D enables a three-dimensional view of gene variation in human metabolism. Nat. Biotechnol. 36 , 272–281 (2018). de la Cruz-López, K. G., Castro-Muñoz, L. J., Reyes-Hernández, D. O., García-Carrancá, A. & Manzo-Merino, J. Lactate in the Regulation of Tumor Microenvironment and Therapeutic Approaches. Front. Oncol. 9 , 492088 (2019). Prado-Garcia, H., Campa-Higareda, A. & Romero-Garcia, S. Lactic Acidosis in the Presence of Glucose Diminishes Warburg Effect in Lung Adenocarcinoma Cells. Front. Oncol. 10 , 807 (2020). Kennedy, K. M. et al. Catabolism of exogenous lactate reveals it as a legitimate metabolic substrate in breast cancer. PLoS One 8 , e75154 (2013). Carlos-Reyes, A., Romero-Garcia, S. & Prado-Garcia, H. Metabolic Responses of Lung Adenocarcinoma Cells to Survive under Stressful Conditions Associated with Tumor Microenvironment. Metabolites 14 , (2024). Zhang, Y. et al. Hypoxia in Breast Cancer-Scientific Translation to Therapeutic and Diagnostic Clinical Applications. Front. Oncol. 11 , 652266 (2021). Tan, H. et al. Ketoglutaric acid can reprogram the immunophenotype of triple-negative breast cancer after radiotherapy and improve the therapeutic effect of anti-PD-L1. J. Transl. Med. 21 , 1–16 (2023). Morris, J. P., 4th et al. α-Ketoglutarate links p53 to cell fate during tumour suppression. Nature 573 , 595–599 (2019). McLain, A. L., Szweda, P. A. & Szweda, L. I. α-Ketoglutarate dehydrogenase: a mitochondrial redox sensor. Free Radic. Res. 45 , 29–36 (2011). Rzeski, W. et al. Alpha-ketoglutarate (AKG) inhibits proliferation of colon adenocarcinoma cells in normoxic conditions. Scand. J. Gastroenterol. 47 , 565–571 (2012). Xiao, D. et al. The glutamine-alpha-ketoglutarate (AKG) metabolism and its nutritional implications. Amino Acids 48 , 2067–2080 (2016). Zdzisińska, B., Żurek, A. & Kandefer-Szerszeń, M. Alpha-Ketoglutarate as a Molecule with Pleiotropic Activity: Well-Known and Novel Possibilities of Therapeutic Use. Arch. Immunol. Ther. Exp. 65 , 21–36 (2017). Yang, H. et al. Extracellular ATP promotes breast cancer invasion and epithelial-mesenchymal transition via hypoxia-inducible factor 2α signaling. Cancer Sci. 110 , 2456–2470 (2019). Isaacs, J. S. et al. HIF overexpression correlates with biallelic loss of fumarate hydratase in renal cancer: novel role of fumarate in regulation of HIF stability. Cancer Cell 8 , 143–153 (2005). Sudarshan, S. et al. Fumarate hydratase deficiency in renal cancer induces glycolytic addiction and hypoxia-inducible transcription factor 1alpha stabilization by glucose-dependent generation of reactive oxygen species. Mol. Cell. Biol. 29 , 4080–4090 (2009). Alderson, N. L. et al. S-(2-Succinyl)cysteine: a novel chemical modification of tissue proteins by a Krebs cycle intermediate. Arch. Biochem. Biophys. 450 , 1–8 (2006). Swainston, N. et al. Recon 2.2: From reconstruction to model of human metabolism. Metabolomics 12 , 109 (2016). Lieven, C. et al. MEMOTE for standardized genome-scale metabolic model testing. Nat. Biotechnol. 38 , 272–276 (2020). Heirendt, L. et al. Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat. Protoc. 14 , 639–702 (2019). Additional Declarations There is NO Competing Interest. Supplementary Files SupplementaryInformation.docx SFigure1.tif SFigure2.tif SFigure3.tif Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4864972","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":338265789,"identity":"ba9d240a-690b-4e64-bbd3-a1719d813123","order_by":0,"name":"Osbaldo Resendis-Antonio","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA/klEQVRIiWNgGAWjYPACZh4G9sYGBJc4LTwHSdTCwCCRgMLFDfilDx/+8HGHtYz8zMdtj25U1MqZM7BffFzAsC2xAYcWyb60NMmZZ9J5DG4nthvnnDlubNnAU2w8g+G2MS5bDM7wmDHzth3mMZBObJPObTuWuOEAT5o0D8NtOdxa+D9//gvUIj/zIFDLv2P1MC08eGxhkGYEamG4wQjU0lCTYHCA/RheWyR72Mwke9uAfjkDdFjOsQOGGw7zMBvPMMDtF34e5scffrZZ28u3H38mnVNTJ29wvP3h44KK2zhDDB0cBkWrAdDBRKoHgjogZn9AvPpRMApGwSgYCQAAk6hUdAbnrEkAAAAASUVORK5CYII=","orcid":"https://orcid.org/0000-0001-5220-541X","institution":"INMEGEN","correspondingAuthor":true,"prefix":"","firstName":"Osbaldo","middleName":"","lastName":"Resendis-Antonio","suffix":""},{"id":338265790,"identity":"27bf1350-ecb6-4cf7-a36c-de498199eb99","order_by":1,"name":"Jorge Arellano-Villavicencio","email":"","orcid":"https://orcid.org/0009-0009-1034-7722","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Jorge","middleName":"","lastName":"Arellano-Villavicencio","suffix":""},{"id":338265791,"identity":"9b3d906e-28c4-449e-ae9b-3c6c9a7dab20","order_by":2,"name":"Aarón Vázquez-Jiménez","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Aarón","middleName":"","lastName":"Vázquez-Jiménez","suffix":""},{"id":338265792,"identity":"cb0ddbf7-e966-4206-b127-22160e97a68f","order_by":3,"name":"Juan Oropeza-Valdéz","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Juan","middleName":"","lastName":"Oropeza-Valdéz","suffix":""},{"id":338265793,"identity":"1570d258-4e4e-4d15-8636-ee10836d6a67","order_by":4,"name":"Cristian Padrón-Manrique","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Cristian","middleName":"","lastName":"Padrón-Manrique","suffix":""},{"id":338265794,"identity":"2e091a96-5c84-4ba8-8602-94218fb016c6","order_by":5,"name":"Heriberto Prado-García","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Heriberto","middleName":"","lastName":"Prado-García","suffix":""},{"id":338265795,"identity":"082e42e2-fa32-4a98-ac6c-4c309e6dd246","order_by":6,"name":"Armando Tovar","email":"","orcid":"","institution":"Instituto Nacional de Ciencias Médicas y Nutrición Salvador Zubirán","correspondingAuthor":false,"prefix":"","firstName":"Armando","middleName":"","lastName":"Tovar","suffix":""}],"badges":[],"createdAt":"2024-08-06 02:35:21","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4864972/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4864972/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":62836985,"identity":"3df2e5d0-4cfd-4ad0-af54-ef0ae2a7c501","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1104227,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eWorkflow for community modeling from expression data.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea) Genome-scale Metabolic Reconstructions Generation\u003c/strong\u003e. The analysis process started with obtaining scRNA-seq expression data from MCF-7 spheroids. Subsequently, using the CORDA algorithm, the presence of metabolic enzymes was inferred based on the expression data. The three cellular profiles previously described by Muciño et al. were generated, and finally, the three GEMs corresponding to each subpopulation were obtained. To ensure the quality of the reconstructions, they underwent analysis using the MEMOTE software.\u003cstrong\u003e b) Community modeling and metabolic coupling\u003c/strong\u003e. Part of the process involved gathering data about the community, including the abundance of each cellular subpopulation in the spheroid during the two days of analysis, the type of culture medium in which they grew, and a specific GEM for each phenotype. The community was constructed using the MICOM tool, which was used to optimize growth. As a final result, a report detailing the metabolic fluxes of each reaction in the individual reconstructions, as well as their interaction with the extracellular environment and population and individual growth, was obtained.\u003c/p\u003e","description":"","filename":"Figure1.png","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/e39fe4f1e2d9f0e3f2377f43.png"},{"id":62837454,"identity":"86752bf0-77eb-4f09-a3a4-06066d276292","added_by":"auto","created_at":"2024-08-20 05:38:28","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":683843,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSpatial distribution simulation with oxygen gradients.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea) Graphical representation of metabolite exchange between each cell group in the community modeling\u003c/strong\u003e. The Venn diagram represents a count of metabolites exchanged between them and those shared by the three collectively in the central part and, in its section, the number of metabolites they secrete to the extracellular medium.\u003cstrong\u003eb) Proposed oxygen gradients for a three-dimensional model.\u003c/strong\u003e 6 combinations for three subpopulations consider that each stratum change goes from its normoxia approach to levels close to hypoxia. \u003cstrong\u003ec) and d)\u003c/strong\u003e \u003cstrong\u003eCommunity Growth Rates\u003c/strong\u003e. The community rate measurement for each subpopulation shows that combination E for day 6 is the highest and combination A for day 19 is the one suggested by the simulation. In both graphs, the community growth on the Y-axis in flow units and on the X-axis the groups of combinations of possible cellular organization scenarios.\u003c/p\u003e","description":"","filename":"figure2.png","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/f376ad861e3971d9e56c3df1.png"},{"id":62837453,"identity":"2d9c0f8b-8303-4cdf-96f7-52be1b37406c","added_by":"auto","created_at":"2024-08-20 05:38:28","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1976087,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMetabolite exchange dynamics and Intracellular metabolism\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea) The secretion and consumption behavior of the subpopulations\u003c/strong\u003e. The first 20 with the highest secretion flow (red) and the first 20 with the highest consumption (blue) were taken. According to their initials, \u003cem\u003eI\u003c/em\u003e corresponds to invasive cells, \u003cem\u003eP\u003c/em\u003e to proliferative cells, and \u003cem\u003eR\u003c/em\u003e to reservoir cells on days 6 and 19. A blank space indicates that this metabolite does not participate in the exchange with the medium by the corresponding cell. \u003cstrong\u003eb) and c) Cross-feeding among subpopulations.\u003c/strong\u003e The exchange suggests the presence of metabolites from specific subpopulations to any of the other 2. The color represents the origin of each metabolite, which is proliferative (green), invasive (red), and reservoir (blue). Following the lines from left to right, we can appreciate the metabolite's origin (secretion and cellular destination (consumption).\u003c/p\u003e","description":"","filename":"Figure3.png","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/39aeb98609468388b0880438.png"},{"id":62836987,"identity":"21fe8867-664b-4ae4-94da-8840f576e6c9","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":1472070,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eEnzymes with metabolic activity.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe metabolic pathways were grouped according to the VMH database (see methods section), in which each reaction is assigned a subsystem such as a) glycolytic activity or glycolysis, b) galactose metabolism, c) tricarboxylic acid pathway or Krebs cycle and d) glutamate metabolism. Only reactions with activity in at least one of the three subpopulations were considered for visualization. To represent the directionality, a color identifier (reddish) corresponds to the directionality of the reaction in its classical pathway described in the database, and (blue) corresponds to a reversal of the directionality of the reaction when observed. The blanks do not exclude the presence of enzymes in the reconstruction; they only reflect whether there was activity. The units of the FBA modeling correspond to the metabolic biotransformation rate, mmol/cell/hour, reflected in the size of the spheres for each reaction.\u003c/p\u003e","description":"","filename":"Figure4.png","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/2146ede7930d30d975dcb9a5.png"},{"id":62837452,"identity":"e61fb644-69ba-48be-a34f-56f3951f9e6c","added_by":"auto","created_at":"2024-08-20 05:38:28","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1955181,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCommunity response to oxygen depletion and KO in silico.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFrom right to left, the gradient from normoxia to hypoxia is represented, where oxygen fluxes are decreasing, and each of the 16 marked points reflects the behavior of the subpopulations about the consumption or production of a metabolite. Values above zero represent secretion and those below zero consumption. Red lines correspond to invading cells, green lines to reservoirs, and blue lines to proliferative cells. Each point on the X-axis represents a different community simulation or scenario. Since this is a deterministic scenario, we do not have any statistics since there are no variations in the repetitions of the simulations. To visualize the behavior of the community with the disposition of the elements of the medium, we selected metabolites that could be used as an energy source on different study days: a) glutamine, c) lactate, e) alpha-ketoglutarate on day six, and b) glutamine, d) lactate, f) alpha-ketoglutarate on day 19. g) Essential metabolites in silico KO. The essentiality analysis generated a list of 17 metabolites whose absence in the medium proved lethal for the community's growth. According to the growth day, red (day 6) and blue (day 19) were recorded depending on the state of the community, such as normoxia (upper part) or hypoxia (lower part), and a midpoint between the two.\u003c/p\u003e","description":"","filename":"Figure5.png","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/4216b3686a4314465609f028.png"},{"id":62836990,"identity":"0299b001-4283-46c4-9f26-ea58e62d90ff","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":507240,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMetabolic interconnection between subpopulations (cross-feeding).\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe schematic representation of the interaction between the three cell subpopulations of MCF-7 spheroid. Day 6 interaction is presented where metabolism modeling suggested a connection between invasive (reddish) and proliferating (green) cells as metabolic sustenance. The primary exchange metabolites are lactate and AKG, which could be used as energy sources by the action of lactate dehydrogenase (LDH) and alpha-ketoglutarate dehydrogenase (AKGD). The lactate released into the extracellular medium (M) could induce lactic acidosis that inhibits glycolytic activity (bibliographic information), which supports the cells' use of alternative carbon sources. On the other hand, the third subpopulation (blue) or the reservoir shares high energy metabolites such as ATP, which, in addition to being energetic, could have secondary functions such as promoting epithelial-mesenchymal transition (EMT) via the HIF1a pathway. Similar to the effect of fumarate in regulating the inhibition of the HIF1a inhibitor prolyl hydroxylated (PHD). Both metabolites theoretically increase HIF1a expression (red dates) and promote transformation from a proliferating to an invasive spheroid, as observed by day 19 abundance.\u003c/p\u003e","description":"","filename":"Figure6.png","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/18b5d4df70e20816e5b7c612.png"},{"id":62837816,"identity":"eaddd240-8ff3-4871-ba5c-c38144d93397","added_by":"auto","created_at":"2024-08-20 05:46:37","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":9471946,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/baf91e4c-eec9-4594-a3d7-c5f4a8b987b9.pdf"},{"id":62836989,"identity":"caeb25e3-6c5b-48d7-9825-60f5ae5ffc27","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"docx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":2286685,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryInformation.docx","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/7770cf9f6e7094246c24b63e.docx"},{"id":62836994,"identity":"73bbaff4-eeb3-47ca-ab1a-af3343ab967a","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"tif","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":9359196,"visible":true,"origin":"","legend":"","description":"","filename":"SFigure1.tif","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/55a779a0d9eb7c46c264b48c.tif"},{"id":62836993,"identity":"dfd0d04d-60d9-4092-94d6-427defc42869","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"tif","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":13292164,"visible":true,"origin":"","legend":"","description":"","filename":"SFigure2.tif","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/21d8699326a8b4af3fad9bef.tif"},{"id":62836995,"identity":"610b4d57-2acc-4f9a-bb60-737c2c0bf8dd","added_by":"auto","created_at":"2024-08-20 05:30:28","extension":"tif","order_by":10,"title":"","display":"","copyAsset":false,"role":"supplement","size":19831764,"visible":true,"origin":"","legend":"","description":"","filename":"SFigure3.tif","url":"https://assets-eu.researchsquare.com/files/rs-4864972/v1/f130abbc29a4049e4377f535.tif"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Intratumoral heterogeneity and metabolic cross-feeding in a three-dimensional model of breast cancer: an in silico perspective.","fulltext":[{"header":"I - Introduction","content":"\u003cp\u003eThe Warburg effect, described by Otto Warburg in 1956, has been one of the most critical metabolic markers in the study of cancer, where tumor cells oxidize high rates of glucose as an energy source to its biotransformation into lactate even in the presence of oxygen\u003csup\u003e1\u003c/sup\u003e, acidifying the tumor microenvironment, favoring cell transformation and invasion into other tissues\u003csup\u003e2\u003c/sup\u003e. However, some tumors, such as in breast cancer, adopt different strategies to obtain energy, such as the incorporation of glutamine\u003csup\u003e3\u003c/sup\u003e to adapt to alterations in their microenvironment and compensate for their energy demand or oxidizing fatty acids whose metabolism is closely related to survival and resistance to chemotherapeutic treatment\u003csup\u003e4\u003c/sup\u003e. In recent years, it has been considered that a metabolic association between the tumor and adjacent cells may ensure tumor survival. This metabolic coupling has been experimentally validated in cancer-associated fibroblasts (CAF) and cancer cells\u003csup\u003e5,6\u003c/sup\u003e. This relationship has been hypothesized to be induced by reactive oxygen species (ROS) or acidification of the tumor microenvironment. Among its consequences, the incentive changes in the microenvironment as a consequence of tumor growth induce an elevated glycolytic activity in CAFs. The behavior benefits tumor cells that could dispose of the lactate released into the extracellular medium in response to unfavorable conditions such as hypoxia or changes in pH. This association has been described as a reverse Warburg effect, a feature accepted as a critical metabolic adaptation in the progression of pathology\u003csup\u003e5,6\u003c/sup\u003e. The transient nature of this association and the lack of information on the conditions and metabolic regulatory pathways involved in its establishment add significant complexity to its study.\u003c/p\u003e\n\u003cp\u003eThe complexity of studying cancer is attributed to the inter-tumoral (between tumors) and intra-tumoral (within the tumor) heterogeneity, which makes it challenging to classify possible metabolic phenotypes or subphenotypes between different types of cancer and cells of the same tumor\u003csup\u003e7\u003c/sup\u003e. The progress of omics technologies, such as single-cell transcriptome sequencing, has opened a new window to understanding tumor phenotypes and heterogeneity from gene expression data\u003csup\u003e8\u003c/sup\u003e. Accordingly, the implementation of single-cell RNAseq techniques (scRNA-Seq) in an in vitro 3D culture of luminal invasive breast cancer subtype A (MCF-7) and a functional enrichment analysis determined the presence of three cell subpopulations coexisting within the same spheroid with characteristics of proliferation, invasion, and a third phenotype of transition or metabolic reservoir \u003csup\u003e9\u003c/sup\u003e. Some questions arise based solely on the expression data, including the metabolic interdependence between these coexisting subpopulations and how we can modulate metabolism to interfere with tumor cell progression and survival.\u003c/p\u003e\n\u003cp\u003eEnrichingic pathways from sequencing data have significantly contributed to the characterization of specific functions in tumor cells and non-tumor tissue, providing relevant information for prognosis and the development of targeted therapies\u003csup\u003e10\u003c/sup\u003e. However, this information does not determine whether functional heterogeneity is supported by metabolism or how cells may couple to promote tumor survival. The cooperative metabolic landscape that may be present inside a solid tumor and the complexity attributed to its cellular heterogeneity opens the door to implementing comprehensive strategies where mathematical formalisms and bioinformatics tools could provide approximations with the capacity to be experimentally corroborated.\u003c/p\u003e\n\u003cp\u003eThe use of computational tools for generating genome-scale metabolic reconstructions (GEMS) from expression data has been suggested \u003csup\u003e11\u003c/sup\u003e; this has benefited from improvements in single-cell and massive sequencing techniques to obtain better-calibrated tissue-specific reconstructions with the ability to correlate with experimentally measured metabolites \u003csup\u003e12\u003c/sup\u003e\u003csup\u003e,\u003c/sup\u003e\u003csup\u003e13\u003c/sup\u003e\u003csup\u003e,\u003c/sup\u003e\u003csup\u003e14\u003c/sup\u003e. Applications in the area of systems biology and metabolic engineering have been handy in simulating the specific production of a metabolite or stimulating cell growth by implementing mathematical formalisms modeling metabolism assuming a steady-state system; this strategy does not consider the quantification of metabolites in the medium or related enzyme kinetic information\u003csup\u003e15\u003c/sup\u003e. The metabolism modeling through flux balance analysis (FBA) has had great validity when implemented to successfully predict the growth rates of some organisms or the maximization of the production of a metabolite of interest, generating a certain degree of validation for metabolic research\u003csup\u003e16\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eUntil now, FBA-based models for the study of cancer have been limited to the representation of a single specific biological trait and the modeling of a single cell or tumor. However, in recent years, progress has been made in developing computational models that use GEMs to study communities \u003csup\u003e17\u003c/sup\u003e. Tools such as MICOM (microbial community) have been implemented and validated in the modeling of complex communities, allowing the analysis of the interactions between individuals and extracellular microenvironments with greater precision \u003csup\u003e18\u003c/sup\u003e. Despite this approach being developed for microbial communities, it is not limited, given that it considers the metabolic pathways without being species-specific. Broadening the use of these tools to study cancer offers a complementary approach, as they allow the easy modification of specific conditions, such as the presence of nutrients in a medium and alterations in the presence of genes or specific metabolic reactions. This ability to efficiently generate hypotheses facilitates the identification of therapeutic targets and the creation of approaches that seek to respond to concrete phenomena.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;This project proposes a strategy to model in silico metabolism from MCF7 spheroids obtained from scRNA-seq data. Algorithms were implemented using the expression profile to identify the presence of enzymes and metabolite transporters and to create tissue-specific metabolic reconstructions of the cell community, constituted by invasive, proliferative, and reservoir cells. The community modeling reflected the metabolic cooperation between the subpopulations by providing the exchange fluxes between each subpopulation, evidencing the presence of the intratumoral reverse Warburg effect and the implementation of alternative carbon sources by the tumor cells in response to the oxygen gradients present due to tumor growth.\u003c/p\u003e"},{"header":"II - Results","content":"\u003ch2\u003e2.1 Community metabolic modeling strategy design.\u003c/h2\u003e\n\u003cp\u003eTo model the metabolism of a cell community, we developed a strategy that involved creating GEMs based on the scRNA-seq data from a spheroid culture of MCF-7. This approach enabled us to generate three reconstructions representing distinct cell subpopulations: Invasive, Reservoir, and Proliferating (refer to Fig 1a). Next, to address the challenge of constructing a community model from these metabolic reconstructions, we proposed the implementation of tools capable of integrating the reconstructions by specifying a growth medium and the relative proportions of each subpopulation. This setup was designed to simulate two temporal stages of cell spheroid growth (days 6 and 19) to elucidate the impact of cell distribution, subpopulation proportions on community growth, and the exchange fluxes with the microenvironment, as illustrated in Fig 1b.\u003c/p\u003e\n\u003ch2\u003e2.2 Specific Sub-population genome-scale metabolic reconstruction from scRNAseq\u003c/h2\u003e\n\u003cp\u003eTo survey the metabolic cross-feeding among the cancer subpopulation in MCF7, we accomplished a genome-scale metabolic reconstruction starting from the scRNA-Seq data associated with the three identifiable subpopulations at two times of the growth process: 6 and 19 days of cultured \u003csup\u003e9\u003c/sup\u003e. As originally discovered by the authors, these subpopulations carry out different biological functionalities (proliferative, invasive, and reservoir), and their interactive crosstalk among them presumably is an essential factor that stabilizes the spheroid. We preserved the original annotation made by the authors. Before proceeding with the metabolic reconstructions, an imputation method for scRNA-seq data was applied to the count matrix to denoise the data and correct the 'dropout events.' This minimizes the noise characteristic of scRNA-seq technology \u003csup\u003e19\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eOverall, the imputation contributed to improving the subsequent metabolic reconstruction procedure. After imputing the count matrix, we proceeded with the reconstruction by applying the Cost Optimization Reaction Dependency Assessment (CORDA) method to the gene expression profiles of each subpopulation. Despite the variety of techniques to accomplish the reconstruction, we selected CORDA due to its ability to create non-minimalistic metabolic models with greater fidelity in representing biological characteristics \u003csup\u003e9,13,20\u003c/sup\u003e. As a result, we created three metabolic reconstructions, one for each cell subpopulation described inside the spheroid. As expected, each reconstruction has distinctive metabolic capabilities reflecting the feasible space of functionalities in each subpopulation; see Supplementary Table 2. To ensure the quality of these reconstructions and their nearness to represent actual metabolic processes, we evaluated each network through the Metabolic Modeling Test (MEMOTE) (see methods), a computational platform to assess the quality of genome-scale metabolic reconstructions through a variety of parameters such as the consistency of the stoichiometric matrix, the evaluation of biomass production, the Gene-Protein-Reaction association, among others criteria. The quality report obtained with MEMOTE indicated that the reconstructions have a subtotal consistency score of 95%, suggesting that they present a proper stoichiometry matrix, mass balance, charge balance, and adequate connectivity between metabolites in the reconstructions, see supplementary information section 4.\u003c/p\u003e\n\u003ch2\u003e2.3 Optimization and Community Modeling\u003c/h2\u003e\n\u003cp\u003eOnce we obtained and proved the quality of each reconstruction, we proceeded to simulate the metabolic capacities when coexisting as a system computationally. We tested the community for each time, at 6 and 19 days of MCF-7 culture (days that were analyzed in the experimental design of Muci\u0026ntilde;o et al. 2020) and constrained by medium (see Fig 1b). This model allowed us to tackle two fundamental questions: postulate the shared metabolites between the subpopulations (cross-feeding metabolic activity) and evaluate how their communication affects metabolic stability and biomass production in the entire community. All reconstructions were integrated and modeled in MICOM \u003csup\u003e18\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eMICOM considers the relative abundance of each subpopulation and the metabolites in the extracellular medium to maximize a community biomass function. Computationally, MICOM combines Flux Balance Analysis by integrating linear optimization with an L2 regularization constraint to estimate the biomass required to keep the biomass for each population (see methods - Community modeling - MICOM). To reproduce the experimental conditions, we defined the extracellular medium in terms of the metabolic composition of culture media used to grow the MCF7 spheroids in the original experiment experimentally, see supplementary Table 1. To simplify our analysis and avoid diffusion processes, we considered that all the media resources are available for each subpopulation at both simulation times (6 and 19 days of progression of the spheroid). Under these assumptions, our metabolic simulation indicated that the community growth rate is 0.026 mmol/cell/h, corresponding to a doubling time of approximately ~ 26 hr on both days; details on duplication time calculations see Supplementary Information section 8. Consistently with our in silico estimations, this double time is nearly associated with the experimental measurements reported for the MCF-7 cell line between 20 and 24 hr doubling time \u003csup\u003e21\u003c/sup\u003e. At this point, our model was a computational platform to create testable hypotheses around the exchanged metabolites with the medium and the metabolic cross-feeding between the subpopulations. This interaction can be represented graphically as groups corresponding to the subpopulations with their respective exchange interactions between them and with the extracellular medium, see Fig 2a. This computational platform can guide design experiments, prove or generate biological hypotheses, and explore strategies for searching metabolic targets for potential therapeutic interventions.\u003c/p\u003e\n\u003ch2\u003e2.4 Effect of pseudo-spatial distribution of subpopulations in MCF-7 spheroid\u003c/h2\u003e\n\u003cp\u003eIn the previous section, we simulated the metabolic interaction among the three subpopulations by neglecting the diffusion process inside the tumor spheroids. However, as revealed by spatial RNAseq, the cellular distribution and arrangement of metabolites in different tumor strata (such as oxygen and carbon sources) significantly influence the intratumoral profiles of metabolic heterogeneity \u003csup\u003e22\u003c/sup\u003e. To emulate this process and quantitatively assess the effect of metabolic gradients on tumor malignancy and stability, we extend our computational genome-scale metabolic modeling by including layers of gradients over oxygen, a metabolite whose concentrations inside the tumor depend exclusively on the size of the spheroid. In the rest of this work, we will assume that oxygen is a metabolite consumed solely by the environment, a key metabolite in establishing intratumoral cellular heterogeneity.\u003c/p\u003e\n\u003cp\u003eWe selected three spatial layers of oxygen concentration in the inner microenvironment to emulate the oxygen gradient. The simulation starts from the outside to the center of the spheroid; the first level emulates cell normoxia, the second level represents an intermediary level of oxygen abundance, and the last level indicates a hypoxia condition. Given the uncertainty around the precise location of the cell subpopulations in the spheroid, we modeled the metabolic interaction over six different spatial configurations, each defined by different spatial arrays of subpopulation and oxygen availability (see Fig 2b). In principle, all these configurations are feasible, at least in an \u003cem\u003ein silico\u003c/em\u003e model; however, we focused on those configurations that maximize the growth rate of all the communities to survey the possible metabolic communication between the subpopulations. Hence, we calculated the community growth rate for each spatial distribution, considering the abundance of the subpopulations previously reported at 6 or 19 days. In Fig. 2c and 2d, we can observe the community's growth rate; the results show the cellular organization that the spheroid could have in both days of study when taking into account the best growth rate of the community.\u003c/p\u003e\n\u003cp\u003eConsequently, the spatial configuration denoted by condition \u003cem\u003eE\u003c/em\u003e in Fig 2c is the ideal combination to optimize the doubling time of the spheroid at day 6. According to intuition, this configuration suggests that the invasive subpopulation is located in the center of the spheroid, the reservoir at a medium level, and the proliferative subpopulation on the outside layers. This optimized distribution suggests that the invasive cells are better adapted in the center of the spheroids, a place with relatively low uptake of oxygen. In contrast, when changing the corresponding population abundances for day 19 of growth, the spatial configuration \u003cem\u003eE\u003c/em\u003e (see Fig 2d) proved the least optimal. Our simulations suggested that placing the invasive cells on the outside enhances community growth (see Fig 2d condition\u003cem\u003e A\u003c/em\u003e). In summary, we conclude that configurations E6 and A19 supply metabolic advantages to optimize the community growth rate at a systemic level.\u003c/p\u003e\n\u003cp\u003eTo analyze the metabolic phenotypes that emerged from these configurations, the following sections will: 1) evaluate the metabolic exchange between the environment and the tumor, as well as the metabolic cross-feeding among populations, and 2) identify the metabolic pathways involved in each subpopulation, based on enzymatic reactions classified through the Virtual Metabolic Human (VMH) database\u003csup\u003e23\u003c/sup\u003e.\u003c/p\u003e\n\u003ch2\u003e2.5 Metabolic mapping of cancer spheroids: secretion, consumption, and cooperativity among sub-populations.\u003c/h2\u003e\n\u003cp\u003eHaving simulated the growth rate at both times of the spheroid, we explored how the metabolic exchange between the tumor and the extracellular milieu is influenced by the cellular organization of the spheroid at both times of spheroid progression (see methods section)\u003c/p\u003e\n\u003cp\u003eTo simplify data presentation, we have generated a Boolean representation of metabolic excretion or consumption between the subpopulations and the environment. The metabolites can be classified into three categories: those secreted and produced simultaneously and alternate between production and excretion as time goes on. For example, we noted that most of the amino acids and carbon sources, such as galactose and glutamine, are consumed by each subpopulation from the media at both times (see Fig 3a). Meanwhile, our modeling suggests that metabolites such as carbon dioxide, ammonia, and succinate are excreted from the community to the media. Furthermore, we identified a third class of metabolites whose production consumption changed as time varied. This last category includes 2-oxoglutarate (AKG), acetoacetate, fumarate, L-lactate, and amino acids such as glycine and L-alanine.\u003c/p\u003e\n\u003cp\u003eIn terms of the metabolic cross-feeding among subpopulations, results showed that all three subpopulations have the metabolites present in the medium; however, some metabolites were not initially found in the culture medium and whose origin is a metabolic by-product of one of the subpopulations. Figure 3a indicates all the metabolites produced by one subpopulation and consumed by another. Among these metabolites, we identified primarily essential amino acids like L-isoleucine, L-leucine, L-lysine, L-methionine, L-phenylalanine, L-threonine, L-tryptophan, L-valine, non-essential like L-asparagine, L-histidine, L-serine, conditional non-essential like L-arginine, L-cysteine, L-glutamine, L-tyrosine and carbohydrates such as galactose. Our model postulated that these metabolites maintain constant consumption regardless of the growth day or subpopulation type. On the other hand, some amino acids do not follow this behavior. For example, glycine, a metabolite consumed by the invasive population, can be taken from the medium and reservoir cells (see Fig 3a). However, proliferative subpopulations can also contribute to this amino acid on day 6. On the other side, we observed that specific metabolites are produced exclusively by one type of subpopulation. For instance, the addition of alanine to the medium is regulated by the invasive cells on both days of the study, similar to proline by the proliferative cells.\u003c/p\u003e\n\u003cp\u003eAs Fig 3b shows, our analysis postulates an intricate metabolic mapping underlying an interdependence between the cancer subpopulations in MCF7 spheroids. The production and consumption of metabolites are crosslinking to meet the metabolic requirements of the growing community. Fig 3b shows the origin of the metabolite from its respective subpopulation and where it is required. Such is the case of alpha-ketoglutarate, fumarate, and lactate, whose consumption and secretion profile changes at both times of the spheroid. The panels of Figure 3b and 3c show the cross-feeding metabolic profile of the community on days 6 and 19. On day six, the proliferative cells secrete alpha-ketoglutarate and lactate, a characteristic marker of the Warburg effect.\u003c/p\u003e\n\u003cp\u003eIn contrast, the invasive cells dispose of this metabolite by incorporating it into their pathways. This effect is also represented on day 19 and can be seen as a hallmark between proliferative and invasive subpopulations. Altogether, the cooperative metabolic profile between the three subpopulations creates a complex scenario of metabolic inner/outer communication at different stages of growth. Among these possibilities, some metabolites behave depending on the day of growth. For instance, AKG is secreted to the medium on day six by proliferating cells. On day 19, the invasive cells show this behavior, exchanging roles between excretion and consumption only between these two subpopulations. By analyzing a portion of the top 20 of the approximately ~350 exchange reactions and ~2000 community reactions, we may be ruling out the possibility that AKG-like switches are also present with metabolites whose exchange fluxes are minimal. These scenarios could add a touch of complexity to the design of metabolic targets due to the metabolic adaptation and reconfiguration of a spheroid over time.\u003c/p\u003e\n\u003cp\u003eWe estimated the metabolic pathways with the highest activity in each subpopulation to assess metabolic adaptations. Then, based on the reconstruction of metabolic pathways, reactions with non-zero flux units (in any directionality) were isolated to determine each pathway's importance in each subpopulation. A common feature is that all three subpopulations have activity in their main pathways of central metabolism, such as glycolysis, the TCA, and the pentose phosphate pathway (see supplementary Fig 1). However, there were cases, such as lipid metabolism (Fatty acid oxidation and Fatty acid synthesis), where a higher count of positive-flux reactions was observed in proliferative and reservoir cells than in invasive cells. On day 19, the invasive cells had more active reactions due to changes in spheroid abundance and size (see supplementary Fig 1).\u003c/p\u003e\n\u003cp\u003eDespite this analysis supplying a proxy of the metabolic activity associated with each subpopulation, the reactions present within each pathway usually could not be the same, or if so, the behavior of the directionality of the metabolic fluxes may vary in those reactions with reversible capacities. The description of directionality changes in reactions provided a more comprehensive approach to tracking how each subpopulation uptakes, transforms, and utilizes the metabolites in the medium for energy generation and biomass production.\u003c/p\u003e\n\u003ch2\u003e2.6 Metabolic activity inside the populations\u003c/h2\u003e\n\u003cp\u003eOnce the central metabolic pathways were identified, they were used to analyze individual behavior and determine the existence of metabolic cross-feeding between subpopulations. Among these, we identified energy pathways such as glycolysis, the Krebs cycle, oxidative phosphorylation, the pentose phosphate pathway, and some other metabolic pathways that might be of interest (see Fig. 4 and Supplementary Fig 2). The simulation is not only limited to providing fluxes of exchange within the medium but also fluxes of intracellular reactions of each subpopulation. This information allowed us to generate a descriptive picture of the behavior of each cell type\u0026mdash;for example, the glycolytic activity of the community. By having galactose as the main carbohydrate according to the original experimental design, we consequently obtained a high activity of the enzyme galactosidase in most of the subpopulations on both days of study except for the invasive cells that show lower fluxes of activity in this enzyme on day 6 of growth (see Fig 4b)\u0026mdash;suggesting that this metabolite is not as essential for the subpopulation in early periods of growth.\u003c/p\u003e\n\u003cp\u003eThe information obtained by simulation suggests that the invasive subpopulation does have glycolytic activity, as seen in Fig 4a, as metabolic fluxes are seen in most of the enzymes involved in glycolysis. However, particular details such as lactate dehydrogenase (LDH) enzyme activity do not reflect a directionality similar to the proliferative cells that initially convert pyruvate to lactate but instead that the invasive cells are in a reversal process of this action and not by carrying out the oxidation of galactose to pyruvate by a classical pathway. Reaction directionality reversal behaviors further solidify the association theory between subpopulations by incorporating metabolites secreted into the medium by neighboring cells.\u003c/p\u003e\n\u003cp\u003eThe origin of metabolites by alternative pathways, such as anaplerotic pathways, often converges to the TCA for generating reducing power and oxidative phosphorylation. Given the results shown in Fig 4c, we have concluded that full participation of TCA is not necessary to satisfy the energy demand of the cells. Although all reconstructions have full TCA, the optimization suggests a preference for certain enzymes that compensate for the objective function of the simulation. For example, given the directionality of the reaction, the model suggests that succinate coenzyme A ligase, one of the critical enzymes of the TCA, consumes succinate to produce succinate and GTP from succinyl-CoA. These reactions could be preceded by the activity shown by alpha-ketoglutarate dehydrogenase, which provides the Succinyl-CoA originating from the presence of alpha-ketoglutarate (see Fig 4c). Both enzymes are shared among the three subpopulations. However, particular cases, such as invasive cells on day 19 and proliferating cells on day 6, activate another section of TCA, such as succinate dehydrogenase, whose reversibility converts fumarate that might be taken up from the extracellular medium when secreted by reservoir cells to succinate (see Fig 4c).\u003c/p\u003e\n\u003cp\u003eThe partial recovery of pathways, such as TCA, suggests that the availability of metabolites such as pyruvate for energy may have different origins. From galactose oxidation or incorporation of lactate from the microenvironment to the involvement of anaplerotic pathways such as alanine incorporation, of which L-alanine transaminase activity was recorded in the invasive and reservoir subpopulations (see Fig 4d). The simulation gave us a broader picture of what alternative energy pathways the cells might use. The simulation also facilitates the analysis of alternative anaplerotic pathways. To identify alternative pathways, we separated already described features of tumor cell metabolism, such as the incorporation and activity of enzymes like mitochondrial glutaminase and glutamate dehydrogenase enzymes, which coordinate for generating alpha-ketoglutarate in TCA, see Fig 4d.\u003c/p\u003e\n\u003cp\u003eThis set of results allowed us to characterize each subpopulation and assign specific metabolic characteristics. According to Muci\u0026ntilde;o et al., 2024, the cells have a glycolytic activity. However, our results complement this assertion by providing the variations of each subpopulation to obtain energy. These differences could be explained by the capacity of each cell to adapt its metabolism to the conditions of its microenvironment and the cooperativity between individuals in the community. Under this criteria, we explored the adaptive capacities of the community in response to stress in its microenvironment. Directly, oxygen was taken as an essential element in the community that could quickly stimulate a change in the subpopulations. As we argued in the next section, the action of different oxygen zones makes evident the importance of some amino acids for the survival of the whole community.\u003c/p\u003e\n\u003ch3\u003e\u0026nbsp;\u003c/h3\u003e\n\u003ch2\u003e2.7 Metabolic response of cancerous subpopulations under gradients of oxygen.\u003c/h2\u003e\n\u003cp\u003eThe stress induced in the community through the oxygen supply of each subpopulation could be an essential factor for the appearance of heterogeneous profiles in the spheroid; this is because one of the variables that allowed us to establish cooperativity among the subpopulations was the presence of a simulated gradient in the community model. To determine the response and capabilities of the community to an adverse scenario, he designed a modeling approach that could represent the decrease in oxygen availability from normoxia to anoxia. Each simulation point reflects the behavior of the subpopulations in the face of oxygen deprivation (see Fig 5). The response was reflected in a greater or lesser preference for selected metabolites as possible energy sources, such as galactose, pyruvate, L-lactate, and AKG, as well as amino acids such as glycine, L-alanine, L-glutamine and L-serine.\u003c/p\u003e\n\u003cp\u003eAn example is glutamine consumption, an amino acid associated with the fundamentals of tumor metabolism whose role has been essential as an alternative carbon source for tumor cells to ensure survival under unfavorable conditions \u003csup\u003e3\u003c/sup\u003e. This metabolite proved crucial for the three subpopulations on both study days. The simulation showed that the proliferating subpopulation had the most glutamine in the medium on day six. As oxygen availability decreases, the invading cells begin to increase their consumption. The opposite is the case on day 19, where the invading cells now have most of the glutamine, and this behavior does not appear to be affected in response to the oxygen gradient (see Fig 5a and 5b)\u003cem\u003e.\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe cells can incorporate metabolites in the medium, such as pyruvate, to obtain energy under normoxic conditions or in response to gradients. Proliferating cells dispose of pyruvate on day 6 (see supplementary Fig 3a), and as oxygen decreases, this behavior is maintained, as well as in the reservoirs, but in a lower proportion. By day 19, under normoxic conditions, the entire community had a low amount of this metabolite, which could be interpreted as a low preference or essentiality. However, as disposition decreases, the invasive subpopulation begins to mobilize pyruvate into its interior (see supplementary Fig 3b). The effects differed significantly from those observed when consuming the primary carbon source, galactose. An almost linear gradient accompanied the decrease in oxygen on both days (see Supplementary Fig 3c and 3d) as a response to oxygen availability and possibly by incorporating alternative sources.\u003c/p\u003e\n\u003cp\u003eAnother type of response observed in the simulation was the cooperativity between subpopulations; in normoxic conditions, there was already an association between proliferating and invasive cells with lactate mobility. However, on day six, the decrease in oxygen increased the transport of lactate into the invasive cells from the proliferating cells and the reservoir. Stages close to zero activated lactate production in the reservoir contribute to the lactate demand in the community (see Fig 5c). This behavior is the first finding to resemble the reverse Warburg effect as an adaptive response to gradients. The simulation showed how the three subpopulations may respond to changes in their microenvironment. This feature was maintained on day 19 (see Fig 5d); however, the exchange interactions diminished as oxygenate approached zero.\u003c/p\u003e\n\u003cp\u003eA further metabolite that showed a peculiar behavior in response to gradients was AKG, which had already been observed to interact between spheroid subpopulations under normoxic conditions, originating in proliferating cells towards invasive cells. However, when hypoxia is approached, the community adapts in such a way that the directionality of this metabolite is inverted, generating the hypothesis that invasive cells have a better capacity to adapt to gradients and provide AKG to proliferating cells to compensate for tumor survival (see Fig 5e).\u003c/p\u003e\n\u003cp\u003eIn the same way as with lactate, by day 19, the spheroid organization established an exchange of lactate from invasive to proliferative from normoxia, and oxygen deprivation creates a decay in the interactions, see Fig 5f. This behavior suggests that metabolic cross-feeding between proliferative and invasive cells is established in the first days of growth; day 6 shows a cellular organization that responds better to oxygen deprivation stress than day 19 since more invasive cells do not have cellular allies to help their survival.\u003c/p\u003e\n\u003ch3\u003e\u0026nbsp;\u003c/h3\u003e\n\u003ch2\u003e2.8 Knock-Out \u003cem\u003ein silico\u003c/em\u003e of Metabolic Enzymes\u003c/h2\u003e\n\u003cp\u003eTo identify those reactions that could be essential for community growth, a Knock-Out (K.O.) essentiality analysis of metabolic reactions was implemented, including exchange reactions with the medium and individual intracellular reactions (see methods section), using growth rate values to determine the impact of K.O., a list of reactions whose absence directly affected the entire community was isolated. This suggested that there is a dependence on metabolites from the medium, especially amino acids, as the most important for survival (see Fig 4b). The relationship was maintained on both days of growth, with the only difference being that Serine might be essential for the day 19 spheroid.\u003c/p\u003e\n\u003cp\u003eTo evaluate if the metabolic dependence of the amino acids on the spheroid is robust in space, we repeated the simulations by mimicking the inner zone of the spheroids through the gradient of oxygen described in the previous section. Hence, we separately accomplished the KO analysis by considering three feasible scenarios of oxygen: normoxia, mid-level, and hypoxia. The essentially analysis obtained from each oxygen environment did not show significant changes between days 6 and 19. These results postulated that in this time scale, the metabolic phenotype of the entire community mainly depends on the amino acids and galactose as a carbohydrate present in the medium located in the external environment (see Fig 5).\u003c/p\u003e\n\u003cp\u003eDespite the diversity of components in the extracellular medium, such as pyruvate and some vitamins, the presence or absence of these components did not show importance in sustaining growth. Furthermore, no intracellular reactions appeared as essential in any of the three sub-populations, suggesting that the environment is the most probable metabolic mechanism to alter its metabolic tumor phenotype.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eOur in silico analysis highlights the importance of external amino acids for tumor growth. Despite some metabolites being essential to growth at any time and oxygen concentration (such as D-galactose, L-glutamine, and L-histidine), others have an interment relevance in the community growth rate, see Figure 5g. For instance, this is the case for L-valine, L-serine, and L-choline, whose relevance to growth rate depends on time and oxygen availability. Overall, these results open an avenue toward an in silico platform to construct testable experimental hypotheses in the design of therapeutic strategies, mainly focused on controlling the growth rate for 3D tumor spheroids.\u003c/p\u003e"},{"header":"III - Discussion","content":"\u003cp\u003e\u003cem\u003eIn vitro,\u003c/em\u003e 3D tumor spheroids are a perfect system to elucidate how intratumoral heterogeneity shapes the metabolic organization and sustains the malignant phenotype. As we have described here, our in silico analysis, which integrates single-cell RNASeq on MCF7, is a powerful tool. It allows us to map the complex cross-feeding between the cellular components of the MCF-7 spheroid and predict how metabolic exchange profiles change as a function of the relative cell abundances and time. This predictive power instills confidence in the model's ability to anticipate metabolic changes.\u003c/p\u003e\n\u003cp\u003eAmong the most remarkable metabolic exchanges that could be identified with the computational model was the secretion and consumption of lactate. This metabolite is used as an energy source by cells in oxygen-restricted conditions, such as invasive cells inside the spheroid. From a diagnostic point of view, lactate has been listed as a biomarker of poor prognosis in cancer patients. This oncometabolite has been positioned as an essential marker due to its involvement as a tumor promoter and facilitator of immune evasion \u003csup\u003e24\u003c/sup\u003e. Our results broaden the picture of lactate involvement beyond the tumor microenvironment regulation and as a potential energy source inside the inner sub-populations. In recent years, the reverse Warburg effect has increased relevance by postulating the importance of lactate as an additional carbon source to produce biomass and not as a product, as described in the classical Warburg effect when cancer cells are near other neighbor cells. For instance, this is the case of the association between Cancer Associated Fibroblast (CAFs) and tumor cells by ensuring cancer survival through an exchange of lactate as an additional carbon source \u003csup\u003e6\u003c/sup\u003e. Although the association has been experimentally induced between the tumor and its stroma, we propose a similar scenario among cells of the same tumor (among the subpopulations) as a possible metabolic adaptation.\u003c/p\u003e\n\u003cp\u003eIn agreement with some in vitro experiments, promoting lactic acidosis in the microenvironment by tumor cells can induce metabolic turnover in adjacent cells by promoting the shift from aerobic glycolysis to an oxidative phenotype \u003csup\u003e25\u003c/sup\u003e. The conditions simulated for invasive cells on day six could be equivalent to an experiment in which low oxygen levels and acidic conditions favor lactate consumption in breast cancer cell lines \u003csup\u003e26\u003c/sup\u003e. This behavior has also been observed in other tumor lines, such as lung cancer \u003csup\u003e26,27\u003c/sup\u003e. These results allow us to question whether this behavior may be a hallmark of glycolytic tumors or a cellular adaptation in response to stimuli from their microenvironment. The organization proposed by the simulation for the spheroid at day 6 is valid when considering that proliferative cells condition the acidic microenvironment, and low oxygen levels promote the transformation of cells in the center of the spheroid to an invasive phenotype \u003csup\u003e28\u003c/sup\u003e. The metabolic connection among the two subpopulations creates a principle of intratumoral cooperativity to ensure the entire community's survival.\u003c/p\u003e\n\u003cp\u003eIn the simulations of gradual oxygen deprivation in the microenvironment, the response of the invasive cells was to increase the rate of lactate consumption, compensated by the proliferative cells to the point that the reservoir subpopulation began to participate by secreting lactate to the medium. See Fig 5c. This indicates that this exchange is enhanced in response to hypoxia. On the other hand, the proposal on day 19 suggests that this association is already established, and oxygen gradients decrease the consumption rate of metabolites such as lactate, see Fig 5d. This may be related to the relative abundance of each subpopulation on day 19 because, on this day, there are more invasive cells than proliferative cells, which decreases the availability of lactate in the medium; the incorporation of pyruvate and high rates of glutamine consumption compensates this demand. All this suggests that the computational model of day six growth would have greater validity for future analyses where the interest is to evaluate the community's response to oxygen deprivation. On the other hand, the day 19 proposal contemplates advanced tumor scenarios, where approximately 50% of the community is invasive, and its metabolism is based on the consumption of L-glutamine (Fig 5b), galactose (supplementary Fig. 3d), serine (supplementary Fig. 3f) and glycine (supplementary Fig 3h).\u003c/p\u003e\n\u003cp\u003eNevertheless, being a computational simulation, these tools allowed us to monitor alternative metabolic pathways that could be interesting. Such as the dynamics of AKG exchange from proliferative cells to invasive cells on day 6. The role of AKG in cancer metabolism has been somewhat controversial. Meanwhile, its potential has been observed as an enhancer of radiotherapies by sensitizing triple-negative breast cancer tumor cells \u003csup\u003e29\u003c/sup\u003e, antagonizing tumor progression in P53-deficient malignancies \u003csup\u003e30\u003c/sup\u003e, and as an ammonium carrier through glutamate synthesis via reversible glutamate dehydrogenase, reducing urea production \u003csup\u003e31\u003c/sup\u003e, however, in tumors such as colon cancer AKG is associated with carcinogenesis \u003csup\u003e32\u003c/sup\u003e. Although glutamine is a semi-essential amino acid for tumor cells, the AKG skeleton is more soluble than glutamine, and its metabolic participation is usually similar, compensating for the metabolic demands of invasive and proliferative cells. AKG can be directly incorporated into the TCA as Succinyl-CoA through the action of alpha-ketoglutarate dehydrogenase to obtain energy via OXPHO \u003csup\u003e33\u003c/sup\u003e, which would correspond to their participation within the metabolic synergy between proliferative and invasive as an alternative source of energy or as metabolic compensation since having a metabolite available in its most accessible form would be energetically more efficient. The ability of AKG to diffuse freely through voltage-dependent transporters facilitates its permeability and, through oxoglutarate/malate transporters, without difficulty, ensures its bioavailability in mitochondria for incorporation into TCA \u003csup\u003e34\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eInterestingly, a survey of the transport fluxes of AKG shows that its origin is from proliferative to invasive cells. However, the effect of hypoxia created a switch in this behavior, reversing the roles. The proposal is that in normoxic conditions, the proliferative cells provide AKG and lactate as metabolic sustenance to the invasive cells; this establishes a picture of metabolic compensation since this stage approximates the standard growth conditions of a spheroid. Fig. 6 reflects the proposed interaction between two cell subpopulations with the most significant contribution to the sustenance of the community, which points to the presence of the Reverse Warburg effect in cells of the same tumor type and a connection between the AKG that compensates for the energy rate of the cells in the center of the spheroid. The reversal of the dynamics, on the other hand, is associated with an increase in the rate of lactate production. Since it is maintained as a response to hypoxia and observed as a feature on day 19, this condition could be related to a marker of poor prognosis. Further study is crucial to elucidate the pathways involved in such associations and to explore how they can be intervened by selecting metabolic targets.\u003c/p\u003e\n\u003cp\u003eThe simulation suggests a connection between the reservoir and proliferative cells by exchanging fumarate and L-proline and providing glycine and Adenosine Triphosphate (ATP) in the medium available for the invading cells, see Fig 6. Despite the metabolic importance of these metabolites, the third subpopulation plays a crucial role in regulating tumor survival. The constant secretion of ATP, specifically in breast cancer, has been described to be involved in invasion processes to other tissues where ATP in the microenvironment could induce hypoxia-inducible factor (HIF1a), triggering epithelial-mesenchymal transition processes even under normoxic conditions \u003csup\u003e35\u003c/sup\u003e. This makes this subpopulation a possible key piece in promoting invasive processes. On another subject, in tumor processes, the deficiency of the activity of enzymes such as fumarate hydratase (FH), responsible for the conversion of fumarate to malate, produces an intracellular accumulation of fumarate and in cases such as renal cancer, the accumulation of this metabolite inhibits the activity of prolyl hydroxylases, which leads to an accumulation of HIF1a \u003csup\u003e36\u003c/sup\u003e. The importance of the exchange of this metabolite lies beyond the regulation of HIF1a, beyond favoring glycolytic activity over oxidative phosphorylation and regulating the sensitization of reactive oxygen species \u003csup\u003e37\u003c/sup\u003e by orchestrating the antioxidant response via NERF2 when fumarate inhibits its regulator \u003csup\u003e38\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eIn summary, the main objective of this project was to develop a strategy that could model the metabolic cross-talk among the MCF-7 cell communities. The modeling implementation tests have also allowed us to deduce and create a proposal for interaction between subpopulations of a spheroid. The proposal designed for the association of the three subpopulations corresponds to the growth period of day 6. In summary, the results suggested the existence of metabolic cooperativity between proliferators and invaders and the possibility that the third reservoir subpopulation fulfills the functions of providing energy in the form of ATP or fumarate. However, the reservoir subpopulation could be preparing the necessary conditions to promote the transformation of invasive cells towards metastatic processes, and the over-stimulation of proliferative cells with fumarate could be an initial marker of transformation toward invasive phenotypes.\u003c/p\u003e\n\u003ch3\u003e\u0026nbsp;\u003c/h3\u003e\n\u003cp\u003eFinally, the approximations performed with computational simulation are created from the information obtained from sequencing data. However, proteomic and metabolomic data are needed to ensure better results and calibrate the reconstructions. This provides the presence or absence of metabolites that could generate false positives. The proposed methodology to perform these in silico analyses was preliminarily validated by bibliographic information that could qualitatively coincide with what is already known about breast cancer metabolism.\u003c/p\u003e\n\u003cp\u003eLastly, the in silico study of the inner metabolic phenotype in cancer provides a framework for designing testable experimental hypotheses in oncological studies. With the advent of single-cell technologies (such as sc-RNASeq and single-cell metabolome), these computational methods will be essential to integrate data and generate hypotheses. Considering cancer as a heterogeneous and communal structure, we envision that these methods will be necessary to move toward the design of treatments at a personalized level. Overall, these efforts allowed us to conclude that cancer is more than one phenotype; instead, the presence of subpopulations contributes with an additional degree of complexity whose inclusion should be addressed to improve the wellness of the patients.\u003c/p\u003e"},{"header":"IV - Conclusions","content":"\u003cp\u003eThe information offered by computational models has had a growing acceptance in recent years, demonstrating the ability to represent complex biological phenomena such as cancer metabolism and its adaptations to variations in its microenvironment. Our model was able to include tissue-specific metabolic reconstructions from transcriptomic data of MCF-7 spheroids and create a community model that allows glimpses of possible pathways of metabolic interconnections between cell subpopulations, as well as the presence of the Reverse Warburg effect and other alternate pathways for energetic survival. These types of models offer approaches that could become comparable, being a great tool to suggest more targeted designs in the case of the study of cancer metabolism. The calibration of these proposals could become more accurate by having more specific metabolic targets to identify the cornerstone of cancer metabolic regulation.\u003c/p\u003e"},{"header":"V - Methods","content":"\u003ch2\u003e5.1 Data\u003c/h2\u003e\n\u003cp\u003eThe data for the study was retrieved from the Gene Expression Omnibus (GEO) under the accession number GSE145633 \u003csup\u003e9\u003c/sup\u003e. The dataset consisted of two temporal points, day six and day 19, for MCF7 spheroids cultivated in L-15 medium (ATCC 30-2008). The researchers associated the three identified clusters with their respective sample labels: Cluster A was linked to the invasive phenotype, Cluster B to the reservoir, and Cluster C to the proliferative phenotype; the cell labeling was conserved, and we took the author\u0026rsquo;s classification.\u003c/p\u003e\n\u003ch2\u003e5.2 Data Imputation - sc-PHENIX\u003c/h2\u003e\n\u003cp\u003eTo recover some of the hidden gene expression due to the noisy nature of single-cell technology, we applied an imputation specialized for scRNA-seq data. We imputed the data using sc-PHENIX (single-cell phenotype recovery via non-linear expression imputation), a recent pipeline with classification performance capabilities currently reported in the literature \u003csup\u003e19\u003c/sup\u003e. The input matrix consisted of 363 cells and 23,923 genes. The sc-PHENIX method uses the count matrix and a UMAP space \u003csup\u003e9\u003c/sup\u003e (https://github.com/resendislab/sc-PHENIX). Imputation was done to recover the missing gene expression from scRNA-seq and to improve differential gene expression. We used the same UMAP space work from the original work of dataset \u003csup\u003e9\u003c/sup\u003e. The parameters for sc-PHENIX were t=5, decay=5, metric='euclidean', and knn=20.\u003c/p\u003e\n\u003ch2\u003e5.3 Reconstruction of Tissue-Specific Metabolic Network - CORDA\u003c/h2\u003e\n\u003cp\u003eBased on the imputed data for each subpopulation, we constructed a genome-scale metabolic model (GEM) using the Cost Optimization Reaction Dependency Assessment algorithm (CORDA) \u003csup\u003e13\u003c/sup\u003e, identified at 6 and 19 days of progression. Using the single-cell RNAseq data from an MCF7 culture, we set up the metabolic activation for every cluster that triggers diverse cellular processes. To start the reconstruction, CORDA requires as input a confidence score for each reaction of the base model. We implemented an agglomerative cluster to get the confidence level for each gene in the data, fixing the number of clusters to 3. N-dimensional values were mapped considering samples in each cluster, and the Euclidean distance was computed to associate genes with similar expressions. The confidence value range was from 1 to 3, relating low, medium, and high confidence scores. Additionally, we set the maximum confidence level to genes differentially overexpressed and the minimum confidence value for genes differentially underexpressed in each cluster; we considered differentially expressed genes for those with |fold change| \u0026ge; 2 and p-value \u0026le; 0.05.\u003c/p\u003e\n\u003cp\u003eThe base model for each genome-scale metabolic reconstruction was taken from Recon 2.2; It integrates 7785 reactions, 5324 metabolites, and 1675 genes \u003csup\u003e39\u003c/sup\u003e. The non-mapped reactions from data were assigned with a confidence value of 0 (unknown confidence). We restricted the range to [-1000,1000] to be closer to the biological phenomena. Finally, given the experimental setup \u003csup\u003e9\u003c/sup\u003e, we modified the upper and lower bound for the exchange reaction relating to the growth medium (Supplementary Table 1). To ensure an exchange of metabolites from the medium to the cytoplasm, we set the confidence level for the medium exchange reaction to 3. All codes were implemented in Python 3, which can be reproduced at https://github.com/resendislab/corda.\u003c/p\u003e\n\u003ch2\u003e5.4 Analysis of Metabolic Reconstruction Quality\u003c/h2\u003e\n\u003cp\u003eTo ensure the reproducibility of results in creating and utilizing genome-level metabolic reconstructions, we utilized MEMOTE (Metabolic Model Testing) \u003csup\u003e40\u003c/sup\u003e. MEMOTE serves as a tool for evaluating the quality of these reconstructions, offering a comprehensive report and an overall assessment. It scrutinizes aspects such as reaction stoichiometry, alignment with databases like BIGG or KEGG, and incorporating a biomass function. In this particular study, all reconstructions evaluated obtained a subtotal score of 94% due to suitable stoichiometry and reaction balance. However, the report shows a 68% total score due to inconsistencies in the IDs with which the MEMOTE database makes the comparisons. However, the simulation of communities with the implemented tool is fine for modeling purposes.\u003c/p\u003e\n\u003ch2\u003e5.5 Community Modeling - MICOM\u003c/h2\u003e\n\u003cp\u003eThe community modeling was conducted by implementing MICOM (Microbial Community) version 0.34.1, a software tool enabling the creation of cellular community models and the analysis of their metabolic dynamics. In this way, community models were generated from models in SBML (Systems Biology Markup Language) format with the capability to investigate exchange fluxes among different cellular populations \u003csup\u003e18\u003c/sup\u003e To perform community modeling on days 6 and 19, various factors were taken into account, including the growth medium specified in Supplementary Table 1 and Table 2, the three GEMs generated using the CORDA algorithm, and the relative abundance of each cellular subpopulation. On day 6, the relative abundances were as follows: Invasive - 0.0977, Reservoir - 0.3910, Proliferative - 0.5113. On day 19, the relative abundances were: Invasive - 0.5739, Reservoir - 0.2739, Proliferative - 0.1522.\u003c/p\u003e\n\u003ch2\u003e5.6 Essentiality analysis by Knock-out of metabolic reactions\u003c/h2\u003e\n\u003cp\u003eThis analysis integrates the community generated through MICOM and the functions provided by the COBRA toolbox \u003csup\u003e41\u003c/sup\u003e. In this process, a total count of the reactions in the entire community was performed, and a reaction-by-reaction Knock-out (K.O.) evaluation was carried out. Optimization focuses on the objective function of the community, which is related to community growth. The community in question consists of a specific number of reactions, and a growth value is assigned to each reaction in response to the alteration of these reactions. Those reactions that, when inhibited in silico, interrupt spheroid growth when the biomass value reaches zero or decreases below its optimal value for growth are considered essential.\u003c/p\u003e\n\u003ch2\u003e5.7 Oxygen gradients: Transport of amino acids and other carbon sources\u003c/h2\u003e\n\u003cp\u003eIn the metabolic reconstructions, individual metabolic reactions have their fluxes restricted by defined upper and lower ranges, which correspond to the directionality of the reactions. To simulate the response of the community to global changes in oxygen availability, the lower limit of the oxygen transport reaction was gradually adjusted, optimizing community growth in each scenario. This adjustment was carried out gradually, decreasing from 0.16 (normoxic conditions) to 0.01 unit decrements until reaching a value of zero, which corresponds to anoxia. The objective behind this modification was to evaluate the community's behavior and determine the existence of metabolic cooperation to guarantee the whole community's survival in the face of this stress induced by the lack of oxygen.\u003c/p\u003e\n\u003ch2\u003e5.8 Fluxes visualization\u003c/h2\u003e\n\u003cp\u003eThe output generated by the community model is presented in CSV format, which contains a list of reactions along with the fluxes associated with each subpopulation. Total fluxes resulting from the simulation for day 6 and total fluxes for the day 19 simulation (see supplementary_file1) are included. In the visualization of reactions related to metabolite exchange with the environment, a filter was applied that identifies reactions whose characteristic corresponds to \"EX_metabolite_m\".This indicates that such reactions are related to population and growth medium exchange processes. Exchange reactions with the medium on day six and exchange reactions on day 19 (see supplementary_file1). To visualize the metabolic fluxes of internal reactions, we manually filtered to select reactions involved in glycolytic metabolism, tricarboxylic acids, and galactose metabolism, Fig 4, and some as pentose phosphate pathway and fatty acid oxidation (see supplementary Fig 3). According to the VMH database \u003csup\u003e23\u003c/sup\u003e\u003cem\u003e, \u003c/em\u003emetabolites were classified, and reactions were grouped according to their metabolic pathway. This visualization was performed by implementing package ggalluvial from the ggplot2 library of R (version 4.2.2). All other figures were created with the Python matplotlib package (version 3.7.1); the information needed to create them can be found in detail in the GitHub link.\u003c/p\u003e\n\u003ch2\u003e5.9 Availability Code and Data\u003c/h2\u003e\n\u003cp\u003eThe data are available through Gene Expression Omnibus (GEO) with accession number GSE145633, obtained from Muci\u0026ntilde;o et al. \u003csup\u003e9\u003c/sup\u003e. The source code and simulation results supporting the findings of this study are publicly available at\u003c/p\u003e\n\u003cp\u003e(https://github.com/resendislab/Modeling_Heterogeneity_Cancer_Metabolism_MICOM)\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003e\u003cstrong\u003eVI. Acknowledgements\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eO.R.-A. thanks to the financial support from CONAHCyT (Grant Ciencia de Frontera 2019, FORDECYT-PRONACES/425859/2020), PAPIIT-UNAM (IA202720), and an internal grant from the National Institute of Genomic Medicine (INMEGEN, M\u0026eacute;xico).\u0026nbsp;J. E. A-V. is a doctoral student from Programa de Doctorado en Ciencias Bioqu\u0026iacute;micas, Universidad Nacional Aut\u0026oacute;noma de M\u0026eacute;xico (UNAM) and received a fellowship from CONAHCYT - CVU 702604.\u0026nbsp;This paper is part of the doctoral thesis and the requirements to obtain a Doctor of Science degree in Science from J.E.A-V.\u0026nbsp;C.P-M is a doctoral student from Programa de Doctorado en Ciencias Biom\u0026eacute;dicas, Universidad Nacional Aut\u0026oacute;noma de M\u0026eacute;xico (UNAM) and received a fellowship from CONAHCYT, M\u0026eacute;xico - CVU 855825.\u0026nbsp;J.J.O-V\u0026apos;s work was supported by CONAHCYT (Grant Ciencia de Frontera 2019, FORDECYT-PRONACES/425859/2020) and UNAM Posdoctoral Program DGAPA (POSDOC). We would also like to thank the Instituto Nacional de Medicina Gen\u0026oacute;mica (INMEGEN), where the project was carried out as part of the student program.\u0026nbsp;Finally, we would like to thank all the authors who contributed directly to the development of this project to improve it and the members of the Human Systems Biology Lab for their support and comments that improved this work over time.\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003e\u0026nbsp;VII. Author Contributions\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eO.R-A and J.E.A-V design the project. J.E.A-V \u0026nbsp;performed the simulations, community modeling, quality testing and essentiality analysis. O.R-A. supervised the research. A.V-J. created the genome-scale metabolic reconstructions. C.P.M. Performed the data imputation. J.E. A-V., O.R-A., A. V-J., J.J. O-V., C. P-M., H. P-G., and A. T-P. drafted and discussed the main manuscript. All authors reviewed the manuscript.\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003eVIII. Competing Interests\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAdditional information\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eCorrespondence and requests for materials should be addressed to Osbaldo Resendis-Antonio.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eWarburg, O. On the origin of cancer cells. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e123\u003c/strong\u003e, 309\u0026ndash;314 (1956).\u003c/li\u003e\n\u003cli\u003eVaupel, P., Schmidberger, H. \u0026amp; Mayer, A. The Warburg effect: essential part of metabolic reprogramming and central contributor to cancer progression. \u003cem\u003eInt. J. Radiat. Biol.\u003c/em\u003e \u003cstrong\u003e95\u003c/strong\u003e, 912\u0026ndash;919 (2019).\u003c/li\u003e\n\u003cli\u003eLi, S. \u003cem\u003eet al.\u003c/em\u003e Glutamine metabolism in breast cancer and possible therapeutic targets. \u003cem\u003eBiochem. Pharmacol.\u003c/em\u003e \u003cstrong\u003e210\u003c/strong\u003e, 115464 (2023).\u003c/li\u003e\n\u003cli\u003eWang, T. \u003cem\u003eet al.\u003c/em\u003e JAK/STAT3-Regulated Fatty Acid \u0026beta;-Oxidation Is Critical for Breast Cancer Stem Cell Self-Renewal and Chemoresistance. \u003cem\u003eCell Metab.\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 136\u0026ndash;150.e5 (2018).\u003c/li\u003e\n\u003cli\u003eFu, Y. \u003cem\u003eet al.\u003c/em\u003e The reverse Warburg effect is likely to be an Achilles\u0026rsquo; heel of cancer that can be exploited for cancer therapy. \u003cem\u003eOncotarget\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 57813\u0026ndash;57825 (2017).\u003c/li\u003e\n\u003cli\u003eLiang, L. \u003cem\u003eet al.\u003c/em\u003e \u0026lsquo;Reverse Warburg effect\u0026rsquo; of cancer‑associated fibroblasts (Review). \u003cem\u003eInt. J. Oncol.\u003c/em\u003e \u003cstrong\u003e60\u003c/strong\u003e, (2022).\u003c/li\u003e\n\u003cli\u003eRoulot, A. \u003cem\u003eet al.\u003c/em\u003e Tumoral heterogeneity of breast cancer. \u003cem\u003eAnn. Biol. Clin.\u003c/em\u003e \u003cstrong\u003e74\u003c/strong\u003e, 653\u0026ndash;660 (2016).\u003c/li\u003e\n\u003cli\u003eWu, S. Z. \u003cem\u003eet al.\u003c/em\u003e A single-cell and spatially resolved atlas of human breast cancers. \u003cem\u003eNat. Genet.\u003c/em\u003e \u003cstrong\u003e53\u003c/strong\u003e, 1334\u0026ndash;1347 (2021).\u003c/li\u003e\n\u003cli\u003eMuci\u0026ntilde;o-Olmos, E. A. \u003cem\u003eet al.\u003c/em\u003e Unveiling functional heterogeneity in breast cancer multicellular tumor spheroids through single-cell RNA-seq. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 12728 (2020).\u003c/li\u003e\n\u003cli\u003eBartoschek, M. \u003cem\u003eet al.\u003c/em\u003e Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing. \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 1\u0026ndash;13 (2018).\u003c/li\u003e\n\u003cli\u003eBordbar, A., Monk, J. M., King, Z. A. \u0026amp; Palsson, B. O. Constraint-based models predict metabolic and associated cellular functions. \u003cem\u003eNat. Rev. Genet.\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 107\u0026ndash;120 (2014).\u003c/li\u003e\n\u003cli\u003eWang, Y., Eddy, J. A. \u0026amp; Price, N. D. Reconstruction of genome-scale metabolic models for 126 human tissues using mCADRE. \u003cem\u003eBMC Syst. Biol.\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 153 (2012).\u003c/li\u003e\n\u003cli\u003eSchultz, A. \u0026amp; Qutub, A. A. Reconstruction of Tissue-Specific Metabolic Networks Using CORDA. \u003cem\u003ePLoS Comput. Biol.\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, e1004808 (2016).\u003c/li\u003e\n\u003cli\u003eHuang, Y. \u003cem\u003eet al.\u003c/em\u003e Characterizing cancer metabolism from bulk and single-cell RNA-seq data using METAFlux. \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 4883 (2023).\u003c/li\u003e\n\u003cli\u003eOrth, J. D., Thiele, I. \u0026amp; Palsson, B. \u0026Oslash;. What is flux balance analysis? \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 245\u0026ndash;248 (2010).\u003c/li\u003e\n\u003cli\u003eEdwards, J. S., Ibarra, R. U. \u0026amp; Palsson, B. O. In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 125\u0026ndash;130 (2001).\u003c/li\u003e\n\u003cli\u003eComputational modeling of metabolism in microbial communities on a genome-scale. \u003cem\u003eCurrent Opinion in Systems Biology\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 46\u0026ndash;57 (2021).\u003c/li\u003e\n\u003cli\u003eDiener, C., Gibbons, S. M. \u0026amp; Resendis-Antonio, O. MICOM: Metagenome-Scale Modeling To Infer Metabolic Interactions in the Gut Microbiota. \u003cem\u003emSystems\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, (2020).\u003c/li\u003e\n\u003cli\u003ePadron-Manrique, C. \u003cem\u003eet al.\u003c/em\u003e Diffusion on PCA-UMAP manifold captures a well-balance of local, global, and continuum structure to denoise single-cell RNA sequencing data. \u003cem\u003ebioRxiv\u003c/em\u003e 2022.06.09.495525 (2022) doi:10.1101/2022.06.09.495525.\u003c/li\u003e\n\u003cli\u003eAngione, C. Human Systems Biology and Metabolic Modelling: A Review\u0026mdash;From Disease Metabolism to Precision Medicine. \u003cem\u003eBiomed Res. Int.\u003c/em\u003e \u003cstrong\u003e2019\u003c/strong\u003e, (2019).\u003c/li\u003e\n\u003cli\u003eJensen, B. L., Skouv, J., Lundholt, B. K. \u0026amp; Lykkesfeldt, A. E. Differential regulation of specific genes in MCF-7 and the ICI 182780-resistant cell line MCF-7/182R-6. \u003cem\u003eBr. J. Cancer\u003c/em\u003e \u003cstrong\u003e79\u003c/strong\u003e, 386\u0026ndash;392 (1999).\u003c/li\u003e\n\u003cli\u003eAhmed, R. \u003cem\u003eet al.\u003c/em\u003e Single-Cell RNA Sequencing with Spatial Transcriptomics of Cancer Tissues. \u003cem\u003eInt. J. Mol. Sci.\u003c/em\u003e \u003cstrong\u003e23\u003c/strong\u003e, (2022).\u003c/li\u003e\n\u003cli\u003eBrunk, E. \u003cem\u003eet al.\u003c/em\u003e Recon3D enables a three-dimensional view of gene variation in human metabolism. \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cstrong\u003e36\u003c/strong\u003e, 272\u0026ndash;281 (2018).\u003c/li\u003e\n\u003cli\u003ede la Cruz-L\u0026oacute;pez, K. G., Castro-Mu\u0026ntilde;oz, L. J., Reyes-Hern\u0026aacute;ndez, D. O., Garc\u0026iacute;a-Carranc\u0026aacute;, A. \u0026amp; Manzo-Merino, J. Lactate in the Regulation of Tumor Microenvironment and Therapeutic Approaches. \u003cem\u003eFront. Oncol.\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 492088 (2019).\u003c/li\u003e\n\u003cli\u003ePrado-Garcia, H., Campa-Higareda, A. \u0026amp; Romero-Garcia, S. Lactic Acidosis in the Presence of Glucose Diminishes Warburg Effect in Lung Adenocarcinoma Cells. \u003cem\u003eFront. Oncol.\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 807 (2020).\u003c/li\u003e\n\u003cli\u003eKennedy, K. M. \u003cem\u003eet al.\u003c/em\u003e Catabolism of exogenous lactate reveals it as a legitimate metabolic substrate in breast cancer. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, e75154 (2013).\u003c/li\u003e\n\u003cli\u003eCarlos-Reyes, A., Romero-Garcia, S. \u0026amp; Prado-Garcia, H. Metabolic Responses of Lung Adenocarcinoma Cells to Survive under Stressful Conditions Associated with Tumor Microenvironment. \u003cem\u003eMetabolites\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, (2024).\u003c/li\u003e\n\u003cli\u003eZhang, Y. \u003cem\u003eet al.\u003c/em\u003e Hypoxia in Breast Cancer-Scientific Translation to Therapeutic and Diagnostic Clinical Applications. \u003cem\u003eFront. Oncol.\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 652266 (2021).\u003c/li\u003e\n\u003cli\u003eTan, H. \u003cem\u003eet al.\u003c/em\u003e Ketoglutaric acid can reprogram the immunophenotype of triple-negative breast cancer after radiotherapy and improve the therapeutic effect of anti-PD-L1. \u003cem\u003eJ. Transl. Med.\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, 1\u0026ndash;16 (2023).\u003c/li\u003e\n\u003cli\u003eMorris, J. P., 4th \u003cem\u003eet al.\u003c/em\u003e \u0026alpha;-Ketoglutarate links p53 to cell fate during tumour suppression. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e573\u003c/strong\u003e, 595\u0026ndash;599 (2019).\u003c/li\u003e\n\u003cli\u003eMcLain, A. L., Szweda, P. A. \u0026amp; Szweda, L. I. \u0026alpha;-Ketoglutarate dehydrogenase: a mitochondrial redox sensor. \u003cem\u003eFree Radic. Res.\u003c/em\u003e \u003cstrong\u003e45\u003c/strong\u003e, 29\u0026ndash;36 (2011).\u003c/li\u003e\n\u003cli\u003eRzeski, W. \u003cem\u003eet al.\u003c/em\u003e Alpha-ketoglutarate (AKG) inhibits proliferation of colon adenocarcinoma cells in normoxic conditions. \u003cem\u003eScand. J. Gastroenterol.\u003c/em\u003e \u003cstrong\u003e47\u003c/strong\u003e, 565\u0026ndash;571 (2012).\u003c/li\u003e\n\u003cli\u003eXiao, D. \u003cem\u003eet al.\u003c/em\u003e The glutamine-alpha-ketoglutarate (AKG) metabolism and its nutritional implications. \u003cem\u003eAmino Acids\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, 2067\u0026ndash;2080 (2016).\u003c/li\u003e\n\u003cli\u003eZdzisińska, B., Żurek, A. \u0026amp; Kandefer-Szerszeń, M. Alpha-Ketoglutarate as a Molecule with Pleiotropic Activity: Well-Known and Novel Possibilities of Therapeutic Use. \u003cem\u003eArch. Immunol. Ther. Exp.\u003c/em\u003e \u003cstrong\u003e65\u003c/strong\u003e, 21\u0026ndash;36 (2017).\u003c/li\u003e\n\u003cli\u003eYang, H. \u003cem\u003eet al.\u003c/em\u003e Extracellular ATP promotes breast cancer invasion and epithelial-mesenchymal transition via hypoxia-inducible factor 2\u0026alpha; signaling. \u003cem\u003eCancer Sci.\u003c/em\u003e \u003cstrong\u003e110\u003c/strong\u003e, 2456\u0026ndash;2470 (2019).\u003c/li\u003e\n\u003cli\u003eIsaacs, J. S. \u003cem\u003eet al.\u003c/em\u003e HIF overexpression correlates with biallelic loss of fumarate hydratase in renal cancer: novel role of fumarate in regulation of HIF stability. \u003cem\u003eCancer Cell\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 143\u0026ndash;153 (2005).\u003c/li\u003e\n\u003cli\u003eSudarshan, S. \u003cem\u003eet al.\u003c/em\u003e Fumarate hydratase deficiency in renal cancer induces glycolytic addiction and hypoxia-inducible transcription factor 1alpha stabilization by glucose-dependent generation of reactive oxygen species. \u003cem\u003eMol. Cell. Biol.\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 4080\u0026ndash;4090 (2009).\u003c/li\u003e\n\u003cli\u003eAlderson, N. L. \u003cem\u003eet al.\u003c/em\u003e S-(2-Succinyl)cysteine: a novel chemical modification of tissue proteins by a Krebs cycle intermediate. \u003cem\u003eArch. Biochem. Biophys.\u003c/em\u003e \u003cstrong\u003e450\u003c/strong\u003e, 1\u0026ndash;8 (2006).\u003c/li\u003e\n\u003cli\u003eSwainston, N. \u003cem\u003eet al.\u003c/em\u003e Recon 2.2: From reconstruction to model of human metabolism. \u003cem\u003eMetabolomics\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 109 (2016).\u003c/li\u003e\n\u003cli\u003eLieven, C. \u003cem\u003eet al.\u003c/em\u003e MEMOTE for standardized genome-scale metabolic model testing. \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cstrong\u003e38\u003c/strong\u003e, 272\u0026ndash;276 (2020).\u003c/li\u003e\n\u003cli\u003eHeirendt, L. \u003cem\u003eet al.\u003c/em\u003e Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. \u003cem\u003eNat. Protoc.\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 639\u0026ndash;702 (2019).\u003c/li\u003e\n\u003c/ol\u003e\n\u003cp\u003e\u003cstrong\u003e\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Cancer Metabolism, Tumor microenvironment, Heterogeneity, MCF-7, Reverse Warburg Effect, Lactate, Systems Biology.","lastPublishedDoi":"10.21203/rs.3.rs-4864972/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4864972/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"Today, intra-tumoral composition is a relevant factor associated with the progression and aggresivity of cancer. Despite it suggests a metabolic interdependence among the subpopulations inside the tumor, a detailed map of how this interdependence contributes to malignant phenotype is still lacking. To address this question, we developed a systems biology approach integrating single-cell RNASeq and genome-scale metabolic reconstruction for mapping the metabolic cross-feeding among the subpopulations previously identified in spheroids of MCF7 breast cancer. By calibrating our model with expression profiles and the experimental growth rate, we concluded that the Reverse Warburg effect emerges as a mechanism to optimize the community growth. Furthermore, through an in silico analysis, we identified lactate, alpha-ketoglutarate, and some amino acids as key metabolites whose disponibility alters the growth rate of the spheroid. Altogether, this work contributes with a strategy for assessing how space and intra-tumoral heterogeneity influence the metabolic robustness of cancer, issues that claim computational strategies to move toward the design of optimized treatments.","manuscriptTitle":"Intratumoral heterogeneity and metabolic cross-feeding in a three-dimensional model of breast cancer: an in silico perspective.","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-08-20 05:30:23","doi":"10.21203/rs.3.rs-4864972/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"5f19e28f-d805-49c2-bb8a-2b9d917bc95f","owner":[],"postedDate":"August 20th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":35839980,"name":"Biological sciences/Systems biology/Computer modelling"},{"id":35839981,"name":"Biological sciences/Cancer/Cancer metabolism"},{"id":35839982,"name":"Biological sciences/Systems biology/Biochemical networks"}],"tags":[],"updatedAt":"2024-08-20T05:30:23+00:00","versionOfRecord":[],"versionCreatedAt":"2024-08-20 05:30:23","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4864972","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4864972","identity":"rs-4864972","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}