Improved flux profiling in genome-scale modeling of human cell metabolism

10 Understanding human cell metabolism through genome -scale flux profiling is of interest to diverse 11 research areas of human health and disease. Metabolic modeling using genome -scale metabolic 12 models (GEMs) has the potential to achieve this, but has been limited by a lack of appropriate input 13 data as model constraints. Here we show that GEM-based flux profiling simulations can be improved 14 with an appropriate input data collection procedure and exo -metabolite exchange flux calculation 15 method, called regression during exponential growth phase (REGP). Our results show that the GEM-16 simulated feasible flux space is constrained to a biologically meaningful region, allowing an exploration 17 of the basic organizing principles of the feasible flux space. These improvements help to fulfil the 18 promise of GEMs as a valuable tool in the study of human metabolism and future development of 19 translational applications. 20

21 Cell metabolism; genome-scale metabolic modeling; flux profiling; feasible flux space 22 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint

23 Metabolism of human cells is a highly complex network of thousands of metabolites and reactions. 24 Alterations in cell metabolism are associated with many complex health conditions such as diabetes, 25 inflammatory diseases, and cancer (1,2). Importantly, the defining feature of metabolism is not the 26 concentrations of metabolites in the cell, but the metabolic fluxes (𝑟) through reactions and pathways 27 (3,4). Intracellular m etabolic fluxes can be experimentally determined through isotope -labeled 28 substrate tracing for a small subset of reactions (5,6), but to systematically profile all fluxes of a cell at 29 the genome-scale, mathematical modeling is necessary. 30 Genome-scale metabolic models (GEMs) is a modeling framework wherein the complete metabolic 31 network of a cell is reconstructed in silico (5,7). GEMs can be used for simulations to calculate the 32 optimal (max or min) fluxes of each reaction, and determine the feasible flux space for the entire 33 metabolic network, using techniques called flux balance analysis (FBA) and flux variability analysis 34 (FVA) (8). FBA and FVA requires a small amount of input data as constraints, typically consisting of 35 measured exchange fluxes ( ± measurement error) of a small number of exo -metabolites, such as 36 glucose, lactate, and amino acids. With these experimentally measured input data, GEM simulations 37 have been shown to be strikingly accurate in microorganisms such as E. coli and S. cerevisiae (9–11). 38 Critically, FBA and FVA assume that cells are in steady -state. Thus, input data for these successful 39 applications of GEM simulations have all been collected during exponential growth. 40 Building on the success of GEM simulations in microbial applications, there is considerable interest in 41 studying human cell metabolism using Human-GEM (3,12,13). The current practice to determine exo-42 metabolite exchange fluxes in human cells is to use the consumption and release (CORE) (4) method 43 (Fig 1A). In this method, exo-metabolite concentrations in the cell culture medi a are measured at a 44 single time point , and e xchange fluxes are then calculated based on the difference between the 45 measured ‘spent’ medi a and the unused (‘fresh’) medi a (Fig 1A) . Thus, CORE -calculated exchange 46 fluxes are not true steady-state exchange fluxes, and the use of these values as constraints for FBA and 47 FVA should be done with caution. 48 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint In this study, we measured exo-metabolite concentrations at multiple time points during exponential 49 growth for a human cell line MCF10A, and calculate d exo-metabolite exchange fluxes by regression 50 during exponential growth phase (REGP ; Fig 1B ). We found that REGP-calculated exchange fluxes 51 (𝑟𝑚,𝑅𝐸𝐺𝑃) were substantially different from CORE-calculated exchange fluxes (𝑟𝑚,𝐶𝑂𝑅𝐸). Using 𝑟𝑚,𝑅𝐸𝐺𝑃 52 as input data for FBA and FVA, we showed that the GEM-simulated feasible flux space was constrained 53 to a more biologically meaningful region, allowing an exploration of the basic organizing principles of 54 the feasible flux space. We anticipate tha t future application of 𝑟𝑚,𝑅𝐸𝐺𝑃 as input data for GEM 55 simulations can rapidly advance our understanding of cell metabolism in diverse applications related 56 to human health and disease. 57

58 Exo-metabolite exchange fluxes at steady-state 59 We measured the exo-metabolite concentrations of exponentially-growing MCF10A cells at five time 60 points during the exponential growth steady-state (Supplemental Table 1), and used both the CORE 61 (Fig 1A) and REGP (Fig 1B) methods to calculate the exchange fluxes of exo-metabolites, 𝑟𝑚. Fig 1C 62 shows the comparison between 𝑟𝑚,𝑅𝐸𝐺𝑃 and 𝑟𝑚,𝐶𝑂𝑅𝐸 in the MCF10A cell lines, as well as the 𝑟𝑚,𝐶𝑂𝑅𝐸 63 of 11 cell lines of the NCI -60 panel that were previously considered reliable (3,4,14). By convention, 64 consumption of metabolites (e.g. glucose) is represented by a negative flux, and release of metabolites 65 (e.g. lactate) by a positive flux. As expected, 𝑟𝑚,𝐶𝑂𝑅𝐸 were comparable between MCF10A cells and the 66 11 cell lines of the NCI-60 panel (Fig 1C). However, the 𝑟𝑚,𝑅𝐸𝐺𝑃 and 𝑟𝑚,𝐶𝑂𝑅𝐸 in MCF10A cells, based on 67 the same raw metabolite measurements and cell count data, were substantially different. For several 68 exo-metabolites, for example glutamate and glycine, 𝑟𝑚,𝐶𝑂𝑅𝐸 values were positive, indicating that cells 69 were releasing these metabolites into the culture media; while 𝑟𝑚,𝑅𝐸𝐺𝑃 values were negative, 70 indicating that cells were consuming these metabolites as nutrients. As the CORE method 71 encompasses both lag phase and exponential growth phase (Fig 1A) , whereas the REGP method 72 calculates the exchange flux during exponential growth only (Fig 1B), this reflects that cell metabolism 73 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint differs between lag phase and exponential growth (Fig 2). Our results indicate that glutamate is 74 released by the cells during lag phase, and consumed during exponential growth (Fig 2B). Similarly, 75 consumption of glycine differs between lag phase and exponential growth (Fig 2C). 76 For other nutrients, such as glucose, glutamine, and other amino acids, i n general we observed that 77 𝑟𝑚,𝑅𝐸𝐺𝑃 for nutrient consumption were larger (that is, more negative) than 𝑟𝑚,𝑅𝐸𝐺𝑃, consistent with a 78 higher nutrient consumption rate during exponential growth compared to the lag phase (Fig 1C). For 79 the metabolic waste lactate, we observed that the exchange flux for lactate release was smaller (that 80 is, less positive) by the REGP method (Fig 1C), suggesting that lactate production is elevated during the 81 lag phase, and reduced during exponential growth. 82 Simulating steady-state cell growth 83 A common way to benchmark FBA- and FVA-based GEM simulations is to estimate the cell growth rate 84 (3,15), which can be easily validated experimentally. To do this, we first constructed cell line-specific 85 GEMs by tINIT (16) using cell-line specific transcriptomics data, generated in-house for MCF10A cells 86 (Supplemental Table 2; GEO accession GSE293588) and mined from the Cancer Cell Line Encyclopedia 87 (17) for the 11 cell lines of the NCI-60 panel. We then specified the metabolites that are present in the 88 cell culture media (Ham’s medium), followed by constraining the exo-metabolite exchange fluxes (Fig 89 1C). For MCF10As, either 𝑟𝑚,𝐶𝑂𝑅𝐸 or 𝑟𝑚,𝑅𝐸𝐺𝑃 were used; for the 11 cell lines of the NCI-60 panel, only 90 the available 𝑟𝑚,𝐶𝑂𝑅𝐸 were used (see Fig 3). Based on these input data as model constraints, the range 91 of feasible in silico growth rates were simulated by maximizing and minimizing biomass production. 92 We found that the 𝑟𝑚,𝐶𝑂𝑅𝐸-constrained MCF10A model was infeasible (Fig 3A), meaning that the in 93 silico cell was unable to ‘grow’ with the CORE -calculated metabolite uptake and secretion rates. In 94 contrast, the 𝑟𝑚,𝑅𝐸𝐺𝑃-constrained MCF10A model was feasible (Fig 3A). Critically, the experimentally 95 measured growth rate fell within the GEM-simulated solution space (Fig 3A) , indicating that GEM -96 simulations are physiologically relevant when using 𝑟𝑚,𝑅𝐸𝐺𝑃 as constraints, but not with 𝑟𝑚,𝐶𝑂𝑅𝐸. 97 Similar to the 𝑟𝑚,𝐶𝑂𝑅𝐸-constrained MCF10A model, most of the 𝑟𝑚,𝐶𝑂𝑅𝐸-constrained models of the 98 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint NCI-60 panel cell lines were also either infeasible, or do not contain the experimentally measured 99 growth rate within the simulated solution space (Fig 3A). 100 As a sensitivity analysis, we included a flexibilization factor ranging from 0-20% for all 𝑟𝑚 used as model 101 constraints, and found that t hat this did not play a role in determining model feasibility or the 102 physiological relevance of the simulations (Fig 3A-B). The 𝑟𝑚,𝐶𝑂𝑅𝐸-constrained MCF10A model only 103 produced physiologically relevant simulations when the flexibilization factor was increased to 700%; 104 even at this point, the 𝑟𝑚,𝐶𝑂𝑅𝐸-constrained model of the SR cell line still performed poorly (Fig 3C). 105 Organization of the feasible flux space 106 We then took the 𝑟𝑚,𝑅𝐸𝐺𝑃-constrained MCF10A model as described above, and added a constraint of 107 the biomass production reaction with the experimentally measured growth rate , to produce a 108 constrained GEM of the MCF10A cell line that makes use of all available data. We used this model to 109 explore the feasible flux space of the entire metabolic network of the cell . Following FVA for every 110 metabolic reaction, we calculated the fractional representation of different metabolic subsystems in a 111 sliding window of 200 reactions, ordered by increasing flux variability (Fig 4; Supplemental Table 3). 112 This analysis showed that metabolic subsystems related to fatty acid metabolism , including for 113 example fatty acid biosynthesis pathways and beta oxidation pathways, exhibited low variability in 114 reaction flux (Fig 4). In contrast, amino acid metabolism (AAM) and most central carbon metabolism 115 (CCM) pathways showed intermediate to high levels of variability (Fig 4) , even though the exo -116 metabolite exchange fluxes used as model constraints were all related to CCM (glucose, lactate) and 117 AAM, consistent with previous observations (13,18). For nucleotide metabolism, we observed two 118 distinct regions in this analysis , one with intermediate variability and another with high variability. 119 Finally, we found that reactions in sphingolipid and steroid metabolism, as well as miscellaneous 120 reactions such as pool reactions and artificial reactions necessary for model simulations, exhibited high 121 flux variability (Fig 4). 122 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint

123 Understanding human cell metabolism at the systems level is of critical interest in many areas of health 124 and medicine. GEM-based simulations have shown very promising applications in microorganisms (9–125 11), but obtaining the necessary input data in human cells has proven to be difficult. A number of 126 methodologies have been developed to leverage transcriptomics data as model constraints (19–23), 127 with mixed results (24), likely because metabolic fluxes are poorly reflected by the abundance of 128 (transcripts of) enzymes in a cell. More recently, exo -metabolite exchange fluxes have been used 129 (3,9,13,25), based on a comparison of exo-metabolite concentrations between the ‘spent’ and ‘fresh’ 130 media (4,13). The caveat of this method is that it violates the steady-state assumption of FBA and FVA, 131 and thus should be used with caution. To address this limitation, here we determined exo-metabolite 132 exchange fluxes by collecting multiple samples exclusively during exponential growth phase (Fig 1B-C). 133 Our results indicated a substantial difference in exchange fluxes in different phases of cell growth (Fig 134 2), underscoring the importance of distinguishing between growth phases when studying cell 135 metabolism. With the exponential growth -phase exchange fluxes as model constraints, GEM 136 simulations were biologically meaningful, with the measured cell growth rate falling within the 137 simulated solution space (Fig 3). This allowed us to explore the entire metabolic network of the cell 138 with a physiologically relevant flux profile, revealing a distinct organization of the feasible flux space 139 by metabolic subsystems (Fig 4). 140 Previously, cell-specific GEMs constrained by the exchange fluxes of glucose, lactate, and threonine (all 141 calculated by the CORE method), were shown to predict the cell growth rate to a reasonable degree 142 of agreement with experimentally measured cell growth rates (3). However, incorporation of 143 additional (CORE-calculated) exo-metabolite exchange fluxes, with a flexibilization factor of up to 20%, 144 lead to a large number of infeasible models (Fig 2A-B); suggesting underlying problems with the model 145 constraints and the biological relevance of the simulations. With the REGP method, model simulations 146 remained feasible with a larger number of measured exo-metabolite exchange fluxes, and the feasible 147 flux space was constrained to a biologically meaningful region (Fig 2). We found that the maximum in 148 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint silico growth rate exceeded the experimentally determined growth rate by approximately 2-fold (Fig 149 3A-B). This suggests that a portion of the consumed nutrients is diverted into non-growth related 150 metabolic tasks, consistent with the notion that, perhaps with the exception of fast -growth cancer 151 cells, human cells do not operate to solely maximize growth (24,26). 152 Our results show that, even though the input constraints of our model were all related to central 153 carbon metabolism (glucose, lactate) or amino acid metabolism (see Fig 1C), there is nevertheless an 154 intermediate level of flux variability in these subsystems (Fig 4). This is in line with previous work 155 showing that these subsystems do not operate at full capacity in growing cells (18,27). In contrast, we 156 found that reactions in fatty acid metabolism exhibit ed low flux variability, while reactions in 157 sphingolipid metabolism and steroid metabolism exhibited high flux variability, likely reflecting the 158 degrees of connectivity (i.e. pathway branching) in these subsystems (3,28). 159 Our findings demonstrate that constraining GEMs with exo-metabolite exchange fluxes calculated by 160 the REGP method allows for accurate model simulations . While this approach demands more 161 resources for exo-metabolite measurements, we believe that it is crucial to profile the metabolic fluxes 162 of human cells at the genome-scale, which can lead to a better understanding of the metabolic process 163 in healthy cells and the identification of potential metabolic targets in diseases. 164 165

166 Cell culture 167 MCF10A cells were purchased from ATCC (Cat# CRL-10317). Cells were cultured in DMEM/F12 (Cat# 168 11320033, Gibco) supplemented with MEGM Mammary Epithelial Cell Growth Medium SingleQuots 169 Kit (Cat# CC-4136, Lonza) without GA-1000, along with 0.1 µg/mL cholera toxin (Cat# BML-G117, Enzo 170 Life Sciences) and 100 U/mL penicillin -streptomycin (Cat# 15140, Gibco). Cells were tested for 171 mycoplasma contamination routinely. 172 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint Cell proliferation assays 173 Absolute cell counts were measured at 22, 26, 30, 46, 50, and 54 h after seeding using the CyQUANT 174 Cell Proliferation Assay kit (Cat# C7026, ThermoFisher Scientific) with a CLARIOstar Plus plate reader 175 (BMG LABTECH). Cell proliferation at all other time points were measured by Incucyte ZOOM (Essen 176 Bioscience), then converted to absolute cell counts based on the corresponding cell counts from the 177 CyQUANT measurements. 178 Biomass determination 179 Cells were harvested with 0.05% trypsin-EDTA (Cat# 25300054, Gibco) and counted using 0.4% trypan 180 blue (Cat# 15250061, Gibco) in a TC20 Automated Cell Counter (BioRad). The cell suspension was 181 transferred into pre -weighed Eppendorf tubes and pelleted by centrifugation at 200 g for 5 mins. 182 Pellets were dried in a microwave at 360 W for 20 mins, and desiccated in a desiccator for >3 days. 183 Exo-metabolite measurements 184 Sampling for exo-metabolites was done during cellular exponential growth phase between 22 -30 h 185 after seeding by collection of culture supernatant. Glucose and lactate concentrations were quantified 186 as described before (29), using an HPLC (Shimadzu) with an Aminex HPX-87H column (Cat# 1250140, 187 BioRad) at 65 °C and an IR detector. The column was eluted with 5 mM H2SO4 at a flow rate of 0.6 188 mL/min for 26 min. Amino acids were quantified as described before (30), with the aTRAQ Kit (Cat# 189 4442673, AB Sciex) using a Nexera UHPLC system (Shimadzu) coupled to a Qtrap 6500+ system (AB 190 Sciex) with a BEH C18 column (150 x 2.1 mm, 1.7 m) (Cat# 186002353, Waters) at 50 °C. A gradient 191 elution of water (eluent A) and methanol ( eluent B), both containing 0.1% formic acid and 0.01 % 192 heptafluorobutyric acid, were used as the mobile phases with a constant flow of 1 mL/min . The 193 following MS parameters were used: Curtain Gas: 50; Collision Gas: Medium; IonSpray Voltage: 5500 194 V; Temperature: 500 °C; Ion Source Gas 1: 60; Ion Source Gas 2: 50. The gradient method was: 2% B 195 from 0–2.5 min, 2% to 40% B from 2.5–3.9 min, held at 40% B until 4.2 min, 40% to 90% B from 4.2–196 6.0 min, held at 90% B until 6.1 min, 90% to 2% B from 6.1–8.0 min. Data acquisition and processing 197 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint were performed using Analyst and MultiQuant 3.0.3 software. Following data QC, exo -metabolite 198 exchange fluxes (𝑟𝑚) were calculated by CORE (4) or REGP (see below). 199 Exo-metabolite exchange flux (𝑟𝑚) calculation by REGP 200 Exo-metabolite exchange fluxes (𝑟𝑚) were determined by calculating the ratio between spent media 201 and unused media samples, then normalizing to the known metabolite concentrations of the cell 202 culture media (DMEM/F12). Exo-metabolite concentrations were normalized for culture volume and 203 cell dry weight , and a linear model was fitted to regress the concentrations against time. 𝑟𝑚,𝑅𝐸𝐺𝑃 is 204 then taken as the slope of the linear model, i.e. the derivative of the fitted metabolite concentration 205 [𝑚̂] with respect to time 𝑡 (Fig 1B). Goodness of fit was determined by R 2, with an arbitrary cutoff of 206 0.7. 207 RNA sequencing 208 RNA was sampled at 27 h after seeding. RNA was extracted using the Quick-RNA Microprep Kit (Cat# 209 R1051, Zymo Research) according to the manufacturer’s protocol. Libraries were prepared from 300 210 ng RNA with the KAPA RNA HyperPrep Kit with RiboErase (HMR) (Cat# KK8502, Kapa Biosystems). RNA 211 fragmentation was performed for a desired library insert size of 200-300 bp by fragmentation for 6 min 212 at 94 °C. Library concentrations were determined using the Qubit HS kit (Cat# Q32854, Invitrogen). 213 Library size distribution was determined using a High Sensititivy DNA analysis (Cat# 5067-4626, 214 Agilent) on a Bioanalyzer 2100 (Agilent). Libraries with an average between 300 -400 bp were loaded 215 on the NextSeq 2000 system (Illumina). Reads were quality controlled, mapped to the human genome 216 hg38, and counted by Seq2science (31), available at https://github.com/vanheeringen-217 lab/seq2science. 218 Genome-scale metabolic modeling 219 The consensus Human-GEM, Human1 v1.12.0 (3), was used for all procedures detailed below. Human-220 GEM is a ‘generic’ model which contains all observed metabolites and reactions in human cells. For 221 each cell line, contextualized models were constructed using tINIT (16), where the generic Human -222 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint GEM is pruned in a cell line-specific manner based on whether or not the (transcript of) an enzyme is 223 expressed. We used transcriptomics data measured in-house (MCF10As) or mined from CCLE (17), with 224 an arbitrary expression level cutoff of 1 TPM . Each cell line-specific model was first constrained with 225 components of the growth media (without specifying the exchange fluxes), and then constrained by 226 the measured exchange fluxes with a flexibilization factor ranging from 0% to 700%. For simulation of 227 growth rate, the range of feasible in silico growth rate was simulated by sequentially minimizing and 228 maximizing the biomass reaction, as implemented in the COBRA toolbox (32). The MCF10A-specific 229 model was further constrained with the measured growth rate to perform FVA for all reactions. The 230 flux variability for each reaction, i.e. max(flux) – min(flux), is sorted from lowest to highest. A sliding 231 window of 200 reactions with a step size of 10 reactions was used to assess the variation in reaction 232 flux across different subsystems. Z-scores were calculated for the fraction of reactions per subsystem 233 in a given window. Exchange reactions, transport reactions, and subsystems with less than 5 reactions, 234 were excluded from the visualization of this analysis in Fig 4. All simulations were performed using 235 MATLAB 2023a (MathWorks, Inc.) with Gurobi solver v10.0.1 (Gurobi Optimizer). 236 Data and code availability 237 Raw RNAseq data are available a t GEO, accession GSE293588. All other data and code used in this 238 paper are available in the GitHub repository (https://github.com/Radboud-YuLab/FluxProfilingREGP). 239 Procedures related to metabolic modeling are implemented in MATLAB (2023a). Numerical analyses 240 and graphics are done in R v4.2.3. 241

242 The Yu lab is supported by grants from the Dutch Cancer Society and the Radboud -Western 243 Collaboration Fund. We thank Niky Thijssen for technical support during the experimental work. 244 Author contributions 245 Conceptualization: RY . Data curation: CH, XJ. Formal analysis: CH, XJ. Funding acquisition: YC, RY . 246 Supervision: YC, RY . Writing – original draft: CH. Writing – review and editing: XJ, YC, RY . 247 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint Competing interests 248 None. 249 250 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint Figures and figure legends 251 252 Figure 1. Exo-metabolite exchange fluxes. (A-B), schematic overview of exo-metabolite exchange flux 253 calculation by the CORE and REGP methods. [𝑚], exo-metabolite concentration; A, area under the 254 curve; [𝑚̂], linear-regression-fitted exo-metabolite concentration (see Methods section) . (C), exo-255 metabolite exchange fluxes (𝑟𝑚) for 11 cell lines in the NCI-60 panel calculated by the CORE method, 256 and for the MCF10A cell line calculated by either the CORE or the REGP method. 257 258 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint 259 Figure 2. Exo-metabolite exchange fluxes in different growth phases. (A), growth curve of MCF10A cells 260 showing a distinct lag phase (blue box) and an exponential growth phase (red box and red line). (B-C), 261 𝑟𝐺𝑙𝑢 and 𝑟𝐺𝑙𝑦 showing distinct metabolic profiles during the lag phase and the exponential growth 262 phase. The REGP method is used to calculate the exchange fluxes during the exponential phase (solid 263 red line). Exchange fluxes in the lag phase is estimated by connecting a straight line from the fresh 264 media sample at t=0 h, to the projected exo -metabolite concentration at t=15 h based on REGP 265 calculations (blue dashed line). 266 267 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint 268 Figure 3. FVA simulation of cell growth rates. (A-C), each cell line-specific GEM model was constrained 269 with the corresponding cell line -specific exo-metabolite exchange fluxes ( 𝑟𝑚) with a flexibilization 270 factor of 0% (A), 20% (B), and 700% (C). Black bars, FVA-simulated minimum and maximum in silico 271 cell growth rate for each cell line. Red dots, experimentally measured growth rate for each cell line. 272 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint 273 Figure 4 . Organization of the feasible flux space of MCF10A cells. The MCF10A -specific GEM was 274 constrained with 𝑟𝑚,𝑅𝐸𝐺𝑃 and the measured growth rate, both with a flexibilization factor of 20%. FVA 275 was performed to determine the flux variability (i.e. the feasible solution space) for every metabolic 276 reaction. The fraction of each metabolic subsystem was calculated in a sliding window of 200 reactions 277 of increasing flux variability, followed by a z-score transformation to facilitate comparison. For ease of 278 visualization, z-scores of > 3 or < -3 were set to 3 or -3, respectively. 279 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint

280 1. Hanahan D, Weinberg RA. Hallmarks of Cancer: The Next Generation. Cell. 2011 Mar 4;144(5):646–281 74. 282 2. Ghesquière B, Wong BW, Kuchnio A, Carmeliet P . Metabolism of stromal and immune cells in health 283 and disease. Nature. 2014 Jul 10;511(7508):167–76. 284 3. Robinson JL, Kocabaş P , Wang H, Cholley PE, Cook D, Nilsson A, et al. An atlas of human metabolism. 285 Science Signaling. 2020 Mar 24;13(624):eaaz1482. 286 4. Jain M, Nilsson R, Sharma S, Madhusudhan N, Kitami T, Souza AL, et al. Metabolite Profiling 287 Identifies a Key Role for Glycine in Rapid Cancer Cell Proliferation. Science. 2012 May 288 25;336(6084):1040–4. 289 5. Nielsen J. Systems Biology of Metabolism. Annu Rev Biochem. 2017 Jun 20;86:245–75. 290 6. Jang C, Chen L, Rabinowitz JD. Metabolomics and Isotope Tracing. Cell. 2018 May 3;173(4):822–37. 291 7. Cook DJ, Nielsen J. Genome -scale metabolic models applied to human health and disease. Wiley 292 Interdiscip Rev Syst Biol Med. 2017 Nov;9(6). 293 8. Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nat Biotechnol. 2010 Mar;28(3):245–8. 294 9. Bordel S, Agren R, Nielsen J. Sampling the Solution Space in Genome -Scale Metabolic Networks 295 Reveals Transcriptional Regulation in Key Enzymes. PLoS Comput Biol. 2010 Jul 15;6(7):e1000859. 296 10. Ishchuk OP , Domenzain I, Sánchez BJ, Muñiz-Paredes F , Martínez JL, Nielsen J, et al. Genome-297 scale modeling drives 70 -fold improvement of intracellular heme production in Saccharomyces 298 cerevisiae. Proc Natl Acad Sci U S A. 2022 Jul 26;119(30):e2108245119. 299 11. Edwards JS, Ibarra RU, Palsson BO. In silico predictions of Escherichia coli metabolic capabilities 300 are consistent with experimental data. Nat Biotechnol. 2001 Feb;19(2):125–30. 301 12. Herrgård MJ, Fong SS, Palsson BØ. Identification of Genome-Scale Metabolic Network Models 302 Using Experimentally Measured Flux Profiles. PLOS Computational Biology. 2006 Jul 7;2(7):e72. 303 13. Zielinski DC, Jamshidi N, Corbett AJ, Bordbar A, Thomas A, Palsson BO. Systems biology analysis 304 of drivers underlying hallmarks of cancer cell metabolism. Sci Rep. 2017 Jan 25;7(1):41241. 305 14. Nilsson A, Haanstra JR, Teusink B, Nielsen J. Metabolite Depletion Affects Flux Profiling of Cell 306 Lines. Trends Biochem Sci. 2018 Jun;43(6):395–7. 307 15. Sánchez BJ, Zhang C, Nilsson A, Lahtvee P , Kerkhoven EJ, Nielsen J. Improving the phenotype 308 predictions of a yeast genome-scale metabolic model by incorporating enzymatic constraints. 309 Molecular Systems Biology. 2017 Aug;13(8):935. 310 16. Agren R, Mardinoglu A, Asplund A, Kampf C, Uhlen M, Nielsen J. Identification of anticancer 311 drugs for hepatocellular carcinoma through personalized genome -scale metabolic modeling. 312 Molecular Systems Biology. 2014 Mar;10(3):721. 313 17. Ghandi M, Huang FW , Jané-Valbuena J, Kryukov GV, Lo CC, McDonald ER, et al. Next-generation 314 characterization of the Cancer Cell Line Encyclopedia. Nature. 2019 May;569(7757):503–8. 315 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint 18. Yu R, Campbell K, Pereira R, Björkeroth J, Qi Q, Vorontsov E, et al. Nitrogen limitation reveals 316 large reserves in metabolic and translational capacities of yeast. Nat Commun. 2020 Apr 317 20;11(1):1881. 318 19. Becker SA, Palsson BO. Context-specific metabolic networks are consistent with experiments. 319 PLoS Comput Biol. 2008 May 16;4(5):e1000082. 320 20. Zur H, Ruppin E, Shlomi T. iMAT: an integrative metabolic analysis tool. Bioinformatics. 2010 321 Dec 15;26(24):3140–2. 322 21. Colijn C, Brandes A, Zucker J, Lun DS, Weiner B, Farhat MR, et al. Interpreting Expression Data 323 with Metabolic Flux Models: Predicting Mycobacterium tuberculosis Mycolic Acid Production. PLOS 324 Computational Biology. 2009 Aug 28;5(8):e1000489. 325 22. Lee D, Smallbone K, Dunn WB, Murabito E, Winder CL, Kell DB, et al. Improving metabolic flux 326 predictions using absolute gene expression data. BMC Systems Biology. 2012 Jun 19;6(1):73. 327 23. Navid A, Almaas E. Genome -level transcription data of Yersinia pestis analyzed with a New 328 metabolic constraint-based approach. BMC Systems Biology. 2012 Dec 6;6(1):150. 329 24. Machado D, Herrgård M. Systematic evaluation of methods for integration of transcriptomic 330 data into constraint-based models of metabolism. PLoS Comput Biol. 2014 Apr;10(4):e1003580. 331 25. Aurich MK, Fleming RMT, Thiele I. A systems approach reveals distinct metabolic strategies 332 among the NCI-60 cancer cell lines. PLOS Computational Biology. 2017 Aug 14;13(8):e1005698. 333 26. Bordbar A, Monk JM, King ZA, Palsson BO. Constraint -based models predict metabolic and 334 associated cellular functions. Nat Rev Genet. 2014 Feb;15(2):107–20. 335 27. Xia J, Sánchez BJ, Chen Y , Campbell K, Kasvandik S, Nielsen J. Proteome allocations change 336 linearly with the specific growth rate of Saccharomyces cerevisiae under glucose limitation. Nat 337 Commun. 2022 May 20;13(1):2819. 338 28. Metabolic Atlas [Internet]. [cited 2025 Apr 4]. Available from: https://metabolicatlas.org/ 339 29. Fletcher E, Feizi A, Bisschops MMM, Hallström BM, Khoomrung S, Siewers V, et al. Evolutionary 340 engineering reveals divergent paths when yeast is adapted to different acidic environments. Metab 341 Eng. 2017 Jan;39:19–28. 342 30. Yu R, Vorontsov E, Sihlbom C, Nielsen J. Quantifying absolute gene expression profiles reveals 343 distinct regulation of central carbon metabolism genes in yeast. Elife. 2021 Mar 15;10:e65722. 344 31. van der Sande M, Frölich S, Schäfers T, Smits JGA, Snabel RR, Rinzema S, et al. Seq2science: an 345 end-to-end workflow for functional genomics analysis. PeerJ. 2023;11:e16380. 346 32. Heirendt L, Arreckx S, Pfau T, Mendoza SN, Richelle A, Heinken A, et al. Creation and analysis 347 of biochemical constraint -based models using the COBRA Toolbox v.3.0. Nat Protoc. 2019 348 Mar;14(3):639–702. 349 350 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted May 15, 2025. ; https://doi.org/10.1101/2025.05.13.653488doi: bioRxiv preprint