Modeling tissue-specific Drosophila metabolism identifies high sugar diet-induced metabolic dysregulation in muscle at reaction and pathway levels

Domain: biology
Published: 2026-01-19
Editorial outcome: major_revision

Abstract
Individual tissues perform highly specialized metabolic functions to maintain whole-body metabolic homeostasis. Although Drosophila serves as a powerful model for studying human metabolic diseases, modeling tissue-specific metabolism has been limited in this organism. To address this gap, we reconstruct 32 tissue-specific genome-scale metabolic models (GEMs) by integrating a curated Drosophila metabolic network with pseudo-bulk single-nuclei transcriptomics data, revealing distinct metabolic network structures and subsystem coverage across tissues. We validate enriched pathways identified through tissue-specific GEMs, particularly in muscle and fat body, using metabolomics and pathway analysis. Moreover, to demonstrate the utility in disease modeling, we apply muscle-GEM to investigate high sugar diet (HSD)-induced metabolic dysregulation. Constraint-based semi-quantitative flux and sensitivity analyses identify altered NAD(H)-dependent reactions and distributed control of glycolytic flux, including GAPDH. This prediction is further validated through in vivo 13C-glucose isotope tracing study. Notably, decreased glycolytic flux, including GAPDH, is linked to increased redox modifications. Finally, our pathway-level flux analyses identify dysregulation in fructose metabolism. Together, this work establishes a quantitative framework for tissue-specific metabolic modeling in Drosophila, demonstrating its utility for identifying dysregulated reactions and pathways in muscle in response to HSD.

Introduction
Multicellular organisms consist of tissues that perform highly specialized metabolic functions1. As such, analyzing tissue-level metabolism is critical to delineate the complex metabolic interplay among tissues and understand organismal physiology2. However, direct measurements of tissue-specific enzyme levels or fluxes are often unavailable. Instead, mRNA levels have served as a proxy to reconstruct context-relevant or tissue-specific genome-scale metabolic models (GEMs)3–5, enabling quantitative metabolic network analyses. In fact, numerous computational algorithms such as mCADRE6, tINIT7,8, CORDA5, and MERGE4, have been developed and applied for reconstruction of tissue-specific GEMs across complex organisms from C. elegans4 to Homo sapiens6. Yet, tissue-specific metabolic models in Drosophila melanogaster are unavailable, limiting quantitative assessments of metabolic processes at tissue-level in this organism9,10.

The high conservation of metabolic genes and functionally analogous organs between humans and fruit flies, together with the short generation time, lifespan, and availability of numerous genetic and analytical tools for Drosophila, has made flies an excellent organism in which to model human metabolic diseases11,12. Various mechanisms underlying dysregulated metabolism associated with type 2 diabetes13,14, aging15, and cancer16 have been elucidated in flies. Moreover, several generic GEMs, such as FlySilico17, iDrosophila18, and Fruitfly119, have recently been developed, enabling systems-level evaluations of metabolic networks in this organism. However, they lack the resolution needed for tissue-specific metabolic analysis.

While GEMs are powerful computational tools for analyzing metabolic network structures and simulating fluxes, experimental techniques such as metabolomics, 13C isotope tracing, and redox proteomics also serve as valuable complementary and validation approaches for quantitative metabolic analyses. Metabolomics provides a comprehensive snapshot of metabolite profiles20,21, while 13C isotope tracing techniques can estimate specific metabolic pathway activity22,23. Furthermore, redox proteomics reveals the extent of redox modifications on sensitive amino acid residues (e.g., cysteine and methionine), providing insight into the post-translational regulation of enzyme activity24–27. However, such experimental analyses are often underrepresented in GEM-focused metabolic studies, limiting both quantitative metabolic investigation and empirical validation when studying complex metabolic diseases. As such, leveraging both computational modeling and experimental approaches is essential for advancing quantitative understanding of complex metabolic disorders like type 2 diabetes (T2D).

T2D is characterized by insulin resistance and hyperglycemia, posing a significant global health challenge28–30. Drosophila effectively models key metabolic features of human T2D through high-sugar diet (HSD) feeding, which impairs glucose homeostasis, and induces mitochondrial dysfunction and tissue-specific metabolic dysregulation13,14,31,32. Among affected tissues, muscle plays a central role in glucose metabolism and insulin sensitivity. While perturbations in glycolytic and TCA cycle activities have been observed in diabetic muscle29,33,34, many questions remain unanswered. Specifically, the extent to which NAD(H)-dependent reactions are altered, the identity of potential rate-controlling steps within glycolysis, and the additionally dysregulated metabolic pathways beyond central carbon metabolism remain to be elucidated.

In this study, we first reconstructed 32 tissue-specific GEMs for adult Drosophila by integrating a curated Drosophila metabolic network model with pseudo-bulk single-nuclei transcriptomics data. These GEMs enabled us to systematically evaluate similarities and differences in metabolic network structures across individual tissues. After validating GEM-based predictions on enriched pathways using targeted metabolomics and pathway enrichment analyses, we demonstrated the utility of tissue-specific GEMs in modeling human metabolic disease. Specifically, we applied the muscle-GEM to simulate metabolic changes induced by a HSD, a well-established Drosophila model of T2D. Constraint-based flux analyses revealed altered fluxes in NAD(H)-dependent fluxes, including a decreased maximal NADH production capacity, and sensitivity analysis indicated distributed control of glycolytic fluxes across several enzymes, including GAPDH. [U-¹³C]-glucose tracing and redox proteomics further confirmed decreased glycolytic flux and revealed that these changes were associated with increased redox modifications of glycolytic enzymes. Pathway-level flux comparisons further highlighted dysregulated fructose metabolism. Together, this work establishes a quantitative framework for tissue-specific metabolic modeling in Drosophila and demonstrates its utility for identifying high sugar diet-induced metabolic perturbations at reaction and pathway levels.

Results
Reconstruction of 32 tissue-specific genome-scale metabolic models (GEMs) in Drosophila melanogaster
To evaluate tissue-specific metabolism in Drosophila, we designed a strategy to reconstruct 32 tissue-specific GEMs and analyze their metabolic network structures and functions (Fig. 1a and Supplementary Data 1a–e). The reconstructed tissue-specific GEMs revealed variations in their metabolic network structures, reflected in the differing numbers of reactions, metabolites, and genes across tissues (Supplementary Fig. 1a). Among 32 tissue-specific GEMs, fat body and oenocytes, analogous to human adipose tissue and liver, contained the highest number of reactions (n = 5447 and 5015), while most neuronal tissues had the fewest number of reactions (n = 2533 to 3077). A similar trend was observed for the metabolite numbers across tissue-specific GEMs. Additionally, we observed a positive correlation between the number of reactions and metabolites (R2 = 0.92), consistent with the expected increase of metabolite numbers as reactions number increases. The number of genes also varied across models, ranging from 565 to 858, with the germline and fat body containing the higher numbers. Similarly, we observed a positive correlation between the number of reactions and genes (R2 = 0.54) (Supplementary Fig. 1b).

Based on these differences, we hypothesized that the tissues performing similar metabolic functions would show similar metabolic network structures. To test this hypothesis, we compared the metabolic network structures of 32 tissue-specific GEMs by specifically comparing their reaction contents (Fig. 1b). Hierarchical clustering analysis revealed 13 distinct clusters of similar metabolic network structures (Supplementary Fig. 1c and Supplementary Data 1f). As expected, tissues of similar functions were grouped together—muscle and indirect muscle; fat body and oenocytes; hindgut and enteroendocrine cells; seven neuron GEMs, and six glia cells distributed across three different clusters. Additionally, to assess how these GEM-based clusters compared to transcriptomics-based clustering, we performed hierarchical clustering of gene expression data and calculated Jaccard indices to quantify overlap (Supplementary Fig. 1d, e). While we observed high overlap between the two approaches for muscle and neuron clusters (Jaccard index = 1) and fat body/oenocytes (0.75), clustering of glia and gut tissues showed lower overlap between the two approaches (Fig. 1c and Supplementary Data 1g). These results indicate that clustering based on metabolic network structures, which incorporate stoichiometric constraints of biochemical reactions, reveals metabolic distinctions not apparent from transcriptomics alone.

Next, to further evaluate differences in metabolic network structures across tissue-specific GEMs, we analyzed their metabolic subsystem coverage. We found that transport reaction subsystem accounted for the largest proportion of all subsystems, comprising approximately 35 ± 4% of all reactions (Fig. 1d). Within this subsystem, extracellular transporter reactions accounted for approximately 80 ± 6%. Interestingly, while muscle group did not have the most transport reactions compared to other tissue groups, it had the highest fraction of the extracellular transport reactions (f = 0.86 ± 0.02) (Fig. 1e and Supplementary Data 1h). Additionally, we calculated the percent difference in subsystem coverage and identified 33 subsystems with more than 50% deviation from mean coverage across tissues (Fig. 1f and Supplementary Data 1i, j). Notably, fat body and oenocyte GEMs had the highest reaction counts, particularly in beta-oxidation subsystems (n = 36.6 ± 15, p < 0.0001) (Supplementary Fig. 1f), supporting their known function in beta-oxidation of fatty acids11,35.

Moreover, we performed a metabolic task analysis to determine the extent to which tissue-specific GEMs could perform 219 pre-defined metabolic tasks (Supplementary Data 1k)7,36. Tissue-specific GEMs could pass an average of 40 ± 10 tasks, with germline cells completing the most metabolic tasks (n = 66) and polar follicle cells pass the least (n = 17) (Supplementary Fig. 1g, h). Additionally, consistent with the known role of fat body in performing gluconeogenesis and trehalose synthesis37, we confirmed that fat body-GEM could synthesize trehalose from various substrates, including alanine, pyruvate, glutamine, and glycerol (Supplementary Data 1l). Furthermore, after categorizing metabolic tasks into seven metabolic systems, we performed Fisher’s exact test and found significant tissue-specific associations with the specific metabolic systems: fat body with carbohydrate/amino acid metabolism, and CNS glia/gustatory neurons with lipid metabolism (Supplementary Fig. 1i and Supplementary Data 1m, n). Taken together, our reconstructed tissue-specific GEMs provide a systems-level framework to evaluate metabolic differences and similarities across tissues in Drosophila, enabling quantitative comparisons of tissue-specific metabolic network structures and metabolic capabilities.

Validation of tissue-specific GEMs through regional metabolomics and pathway analysis
Validating the predictions from tissue-specific GEMs is essential to establish confidence in GEM-based analyses. Although gene essentiality analysis is commonly used for validation38,39, in vivo growth rate data were unavailable for the tissues we investigated. Thus, we performed targeted metabolomics and pathway enrichment analysis, while comparing these results to GEM-based pathway analysis. We hypothesized that enriched pathways identified through metabolomics would also be represented in the corresponding tissue-specific GEMs. To test the hypothesis, we first profiled 303 polar metabolites in four dissected Drosophila regions—head (containing neuronal and glial cells), thorax (containing muscle cells), gut (containing hindgut/enteroendocrine cells), and abdomen (containing fat body, oenocytes, and other tissues) (Fig. 2a and Supplementary Data 2a).

Principal component analysis and hierarchical clustering revealed distinct metabolite profiles among these regions (Fig. 2b, c). Specifically, we identified 10 significantly enriched metabolites in thorax, 13 in head, 20 in gut, and 1 in abdomen (Log FC > 2 and p-value < 0.05) (Fig. 2d and Supplementary Data 2b). Qualitatively, these enriched metabolites aligned with known tissue functions. For instance, in thorax, energy-related metabolites such as adenosine, adenylosuccinate, adenine, and AMP were enriched, consistent with muscle’s high demand on energy metabolism compared to other tissues11 (Supplementary Fig. 2a). In heads, 4-aminobutyrate (GABA), N-acetyl-L-aspartate, GMP, and ascorbates, were enriched. GABA is a major inhibitory neurotransmitter, and N-acetyl-L-aspartate is often used as a marker for neuronal health40,41. In gut, including Malpighian tubules (MT), dietary-related metabolites, such as pyridoxine, betaine, and uric acids, were enriched, consistent with their dietary origin and uric acids generated from MT’s purine excretion42–45. In abdominal carcass, containing fat body and oenocytes, nucleotide and one-carbon metabolites were enriched, such as xanthosine and hypoxanthine, methionine s-adenosyl-l-methioninamine and 5-methyl-THF.

Next, we performed KEGG pathway enrichment analysis using the hypergeometric test based on the enriched metabolite sets46 (Fig. 2e and Supplementary Data 2c). The most significantly enriched pathway was purine metabolism (KEGG ID: dme00230) for thorax, pyrimidine metabolism (dme00240) for gut, and one carbon pool by folate (dme00670) for head and abdomen (Supplementary Fig. 2b). We also identified six pathways commonly enriched across all regions, such as TCA cycle, glyoxylate, purine, and pyrimidine metabolism (Supplementary Fig. 2c). Moreover, several region-specific enriched pathways were identified, such as starch/sucrose metabolism, histidine metabolism, and fatty acid degradation in thorax.

Next, we compared these metabolomics-derived enriched pathways to those predicted from subsystem coverage analysis of tissue-specific GEMs, using weighted Jaccard index47,48 (Fig. 2f and Supplementary Note 1). In brief, the index quantifies the degree of pathway overlap between the two datasets. We expected higher index values for tissues associated with specific regions (e.g., muscle with thorax49 and fat body/oenocytes with abdomen50). Indeed, muscle-GEM showed the highest index value with thorax (index: 0.8 ± 0.1), compared to those for other regions (head: 0.48 ± 0.2, abdomen: 0.44 ± 0.1, gut: 0) (Fig. 2g and Supplementary Data 2d–f). The overlapped pathways included butanoate, glycolysis, lipoic, and starch metabolism. Additionally, fat body-GEM also showed high index value with abdomen, with overlapped pathways including cysteine, nicotinamide, and pantothenate metabolism (Fig. 2h). However, pathway overlaps were less distinct for neuron/glia cells and gut groups, as no uniquely enriched pathways were identified for these regions in our datasets (Supplementary Fig. 2d). Furthermore, when we repeated this analysis using gene expression data, we observed significantly lower Jaccard index values across all tissues, with no distinct trends (Supplementary Fig. 2e and Supplementary Data 2g, h). This indicated that GEM-predicted enriched pathways more accurately reflect metabolomics-derived pathways, compared to gene expression-derived enriched pathways. Taken together, these results validate GEM-predicted pathway analysis, particularly for muscle and fat body.

Constraint-based flux analyses predict perturbations in NAD(H)-dependent reactions in muscle under high sugar diet
Among the reconstructed tissue-specific GEMs, we further applied the muscle-GEM to evaluate how muscle metabolism could be rewired in response to high sugar diet (HSD), which induces T2D-like phenotypes in Drosophila14,31,51. We first defined a HSD-muscle-GEM by constraining reaction rates to simulate metabolic phenotypes observed in type 2 diabetic muscle based on the literature evidence28,29,33,34,52–54 (Fig. 3a and Supplementary Data 3a, b). These constraints included glucose uptake rates, GAPDH, and several reactions in TCA cycle53,55. Consistent with our applied constraints, flux variability analysis with sampling (FVA-sampling) showed a reduction in glucose uptake rate, along with decreased fluxes through GAPDH, OGDH, and SDH (Supplementary Fig. 3a, b, Supplementary Note 2, and Supplementary Data 3c, d). The decrease in model-generated glucose uptake rate was confirmed experimentally by measuring glucose uptake rates in w1118 male flies fed with HSD (Supplementary Fig. 3c).

Next, to identify perturbed reactions beyond those we directly constrained, we systematically compared differential fluxes obtained from flux balance analysis (FBA), FVA-sampling, and parsimonious FBA (pFBA) (Supplementary Note 3). In brief, we performed pFBA to complement FVA-sampling analysis, which could generate non-physiologically high fluxes in reactions such as NAD(H)-dependent cycling reactions and generate fluxes near boundaries due to under-constrained nature of network (Supplementary Fig. 3d–g and Supplementary Data 3e, f). By comparing these analyses, we identified 77 reactions showing consistently decreased fluxes and 18 with increased fluxes in the HSD-muscle-GEM relative to the unconstrained model (Supplementary Fig. 3h and Supplementary Data 3g). Among the reactions with decreased fluxes, many were associated with transport, aminoacyl-tRNA biosynthesis, and glycolysis, with associated genes including Dic1, Cd98hc, Gapdh1/2, Aldo, Eno, Pfk, and aminoacyl-tRNA synthetases such as AsnRS and AspRS (Fig. 3b). In contrast, reactions with increased fluxes were mainly linked to transport, aromatic amino acid biosynthesis, and TCA/glyoxylate metabolism, with associated genes including Hrm, Cg6231, Got1/2, Kdn, Scsa1, and Mdh2.

As perturbed NADH metabolism has also been implicated in T2D56,57, we further investigated individual NAD(H) dependent reactions to identify specific reactions showing altered fluxes. By comparing differential fluxes of 135 active NAD(H) dependent reactions across three different flux simulations, we identified 23 reactions showing consistently decreased fluxes and 3 reactions with increased flux (Supplementary Fig. 3i, j and Supplementary Data 3h–j). Associated genes with decreased fluxes included Dhfr, Ldh, Gs, and Gapdh1/2, while those with increased fluxes included Mdh2 and CG5955. Profiling differential fluxes of NAD(H)-dependent reactions further highlighted a variety of significantly perturbed reactions and their magnitudes (Fig. 3c). Next, to further investigate whether the network’s maximum capacity of NADH production was perturbed, we introduced an artificial NADH demand reaction and maximized its flux (Fig. 3d). This analysis revealed approximately a 27% reduction in NADH maximum production capacity in the HSD-muscle GEM, accompanied by altered contributions of individual NADH-producing reactions (Fig. 3e and Supplementary Fig. 3k). Although this change does not directly predict the cellular NAD⁺/NADH ratio, it suggests impaired NADH turnover, consistent with our experimental observation of a decreased NAD⁺/NADH ratio in thoracic muscle of w1118 flies fed with HSD (Fig. 3f). Together, our constraint-based semi-quantitative flux analyses predicted perturbations in specific subsystems, reactions, and genes in muscle under HSD, highlighting disrupted NAD(H)-dependent reactions characterized by altered flux distributions and decreased NADH production capacity.

Sensitivity analysis reveals distributed control of glycolytic flux
Given that NAD(H)-dependent reaction fluxes were perturbed in the HSD-muscle GEM and GAPDH functions as one of the key NAD⁺-dependent glycolytic enzymes, we examined whether GAPDH could serve as a rate controlling step in glycolysis under HSD condition. To test this, we performed a sensitivity analysis (See “Methods”). In brief, we systematically decreased the baseline flux of each glycolytic reaction, defined as the median flux from FVA-sampling, and quantified the total pyruvate consumption flux as a proxy for glycolytic output flux (Fig. 3g and Supplementary Fig. 3l). At relatively larger perturbation (e.g., 20 to 50%), most of the lower glycolytic enzymes showed positive sensitivity coefficients, suggesting that decreases in these fluxes led to reductions in total pyruvate consumption (Fig. 3h and Supplementary Data 3k). In contrast, at smaller perturbations (e.g., 5 and 10%), several enzymes, including Hex-a/b, Pgk, and Gapdh1/2, showed strong negative sensitivity coefficients, suggesting that decreases in these fluxes may trigger compensatory increases in total pyruvate consumption through alternative reactions. Indeed, evaluation of individual pyruvate-consuming reactions revealed diverse responses to perturbations. Specifically, while alanine transaminase (CG1640) consistently showed decreased flux in response to high perturbation (e.g., 50%), monocarboxylate transporter (MCT1), or other reactions showed increased fluxes in response to varying perturbations, reflecting the nonlinear and compensatory nature of metabolic network (Supplementary Fig. 3m). When comparing absolute magnitudes of sensitivity coefficients across all perturbations, we found Aldo exhibited the highest average sensitivity coefficient (0.13), with Pgk (0.12), Gapdh1/2 (0.12), and Hex-a/b (0.11) showing comparable sensitivity (Fig. 3h). Altogether, these results suggest that glycolytic flux is not solely controlled by GAPDH as initially hypothesized, but is distributed among several other enzymes and varies with perturbation magnitude under this condition.

Model-predicted decreases in glycolytic flux, including GAPDH, validated through 13C-glucose tracing
To validate the predicted fluxes in glycolysis, particularly involving GAPDH, we performed in vivo 13C-glucose isotopic tracing experiments. After feeding w1118 male flies with a HSD containing uniformly labeled 13C-glucose tracer for five days, we dissected the thoracic muscle and evaluated both the metabolite intensities and labeling patterns of downstream metabolites using mass spectrometry (Fig. 4a). As the [U-13C]glucose is metabolized in cells, glycolytic intermediates will have M+6 or M+3 mass isotopomers, and TCA cycle metabolites will show M+2 mass isotopomers after the first cycle (Fig. 4b). Based on our sensitivity analysis, where we observed strong control around Aldo, Pgk, and Gapdh, we hypothesized that a bottleneck step would result in decreased fractional labeling of downstream isotopomers, accompanied by accumulation of upstream substrates and decreased levels of downstream products.

Indeed, we observed that upper glycolytic intermediates increased by approximately 23%, including F6P, F16BP, and glyceraldehyde 3-phosphate, a substrate of GAPDH (Supplementary Fig. 4a, b and Supplementary Data 4a, b). Conversely, lower glycolytic intermediates, such as 3PG and lactate, decreased by 27% and 14%, respectively. Moreover, the fractional labeling of lower glycolytic intermediates was significantly decreased compared to that of upper glycolytic intermediates (Fig. 4c). Specifically, the M+3 fractional labeling of 1,3-BPG, the product GAPDH, was nearly 90 % lower in response to HSD, indicating a significantly reduced GAPDH activity. In TCA cycle, the abundances of citrate, succinate, fumarate, malate, and glutamate, increased by 40%, 34%, 14%, 97%, and 24%, respectively (Supplementary Fig. 4c, d and Supplementary Data 4c). Moreover, their M+2 fractional labeling decreased relative to M+6 glucose, suggesting a decreased contribution of glucose-derived carbon to the TCA cycle (Fig. 4d).

Furthermore, to investigate the extent to which our model-predicted fluxes align with experimental measurements, we performed a correlation analysis by mapping the flux of each reaction in the HSD-muscle-GEM to the fractional labeling of its corresponding metabolite product. Notably, we observed a strong correlation between the fractional labeling data and model-predicted relative fluxes within glycolysis (ρ = 0.81, p = 0.03), supporting the accuracy of model predictions for this pathway (Fig. 4e). In contrast, correlation in the TCA cycle was weaker (ρ = 0.12), potentially due to the contributions from other carbon sources (e.g., fatty acids58) accounting the predicted fluxes in the TCA cycle (Supplementary Fig. 4e). Taken together, these results confirm the model-predicted flux changes, specifically for glycolysis, and reveal GAPDH as one of the key regulatory steps in muscle in response to HSD.

Model-predicted decreases in glycolytic flux correlate with increased redox modification of glycolytic enzymes
Since redox modification, including oxidation of cysteine residue, can directly alter enzyme activity such as GAPDH24,26,59,60, we were wondering whether the decreased glycolytic fluxes, particularly GAPDH, were linked to redox modifications (Fig. 5a). Given the increased oxidative stress observed in diabetic muscle53,61,62, we hypothesized that redox modification at cysteine or methionine residues would also be elevated in HSD-fed fly muscle. Indeed, PCA of redox proteomics revealed greater variability in peptide oxidation profiles in HSD samples compared to NSD controls (Fig. 5b). Among 1049 detected peptides, 189 peptides exhibited significantly increased redox modifications (|Log FC| > 0.5 and adjusted p-value < 0.05) (Fig. 5c, Supplementary Fig. 5a and Supplementary Data 5a). Next, based on these significantly oxidized peptides, we performed KEGG over-representation analysis63 to identify pathways enriched for redox-modified enzymes (Fig. 5d and Supplementary Data 5b). Consistent to prior studies, the most affected pathways were oxidative phosphorylation (dme00190), a well-known site of mitochondrial dysfunction observed in diabetic muscle64,65. Within this pathway, we observed increased redox modifications in peptides, including succinate dehydrogenase, NADH dehydrogenase, V-ATPase, citrate synthase, and isocitrate dehydrogenase (Supplementary Fig. 5b). Interestingly, glycolysis was also significantly affected, with increased redox modifications observed in peptides, including Ald1, Gapdh1/2, Pyk, Pgi, Pfk, Pgm1, Eno, Fdh, and Pgk (Fig. 5e).

To further evaluate the relationship between the extent of redox modification and predicted fluxes in glycolysis, we performed a correlation analysis and found a significant negative correlation (r = −0.78, p = 0.025) (Fig. 5f). Specifically, Pyk exhibited the highest level of redox modification and corresponded to the largest flux decrease, whereas Pgm1 showed minimal redox modification and the smallest flux change (Supplementary Fig. 5c, d). Notably, protein levels of most glycolytic enzymes remained unchanged, except for Pyk, suggesting that redox modifications, rather than enzyme abundance, more strongly associate with changes in glycolytic flux (Supplementary Fig. 5e, f).

Moreover, since our flux analyses and isotope tracing results revealed GAPDH as one of the key perturbed steps in glycolysis, we further evaluated this enzyme. In Drosophila, two isoforms of GAPDH are present, both closely related to human GAPDH and known to perform partially redundant functions in glycolysis (Supplementary Fig. 5g)66. While the enzyme levels of both GAPDH1 and GAPDH2 remained unchanged (Supplementary Fig. 5h), we observed significant redox modifications at Met-40 residue in GAPDH1 and at Cys-130, Met-127, 141 and 172 in GAPDH2 (referenced to Uniprot: P07486 for GAPDH1 and M9PJN8 for GAPDH2) (Fig. 5g and Supplementary Fig. 5i, j). Structural modeling using AlphaFold further predicted that Met-40 is in near to the NAD+ binding site (~9 Å), suggesting that redox modification at this residue could indirectly affect NAD+ binding and impair activity (Supplementary Fig. 5k). Additionally, to further investigate the functional role of GAPDH in muscle, we used a muscle-specific Mhc-Gal4 driver combined with a temperature-sensitive Gal80ts system to knock down GAPDH1 in fly muscle. Indeed, muscle-specific downregulation of GADPH1 led to a significant decline in climbing ability in male flies fed with HSD starting at day 10, suggesting a physiologically important role for GAPDH1 in maintaining muscle function (Fig. 5h and Supplementary Fig. 5l). In summary, these results indicate a strong association between increased redox modifications and decreased glycolytic fluxes, identifying the redox modification of GAPDH and its knockdown causing climbing defects.

Pathway-level flux analysis reveals dysregulated fructose metabolism
Next, to identify pathway-level perturbations beyond glycolysis under HSD conditions, we evaluated differential pathway fluxes obtained from FVA-sampling and pFBA analyses (Fig. 6a and “Methods”). In brief, we defined pathway flux as the average of non-zero flux magnitudes within each pathway. Using this approach, we identified 77 significantly perturbed pathways, with 33 pathways overlapping between flux sampling and pFBA analyses (Fig. 6b and Supplementary Data 6a). Consistent with our prior findings, we observed significant decrease in pathway fluxes through glycolysis, oxidative phosphorylation, and the TCA cycle, with fructose showing the most pronounced decrease in the HSD muscle compared to control. Among the top perturbed pathways with increased pathway fluxes, we found fatty acid metabolism (e.g., β-oxidation, desaturation, elongation, and omega-3/6

Methods
Reconstruction of tissue-specific GEMs
We reconstructed 32 tissue-specific GEMs for Drosophila melanogaster by integrating a curated Drosophila metabolic network model (FlySilico17) with pseudo-bulk single-nuclei transcriptomics data. The transcriptomics data were obtained from a previously published dataset (Supplementary Data 1a). We used the mCADRE algorithm6 to reconstruct tissue-specific GEMs, which utilizes gene expression data to prune the generic metabolic network while maintaining the metabolic tasks. We defined the metabolic tasks based on the literature (Supplementary Data 1k). The reconstruction process was performed using the COBRA Toolbox in MATLAB.

Metabolomics analysis
We performed targeted metabolomics on four dissected Drosophila regions: head, thorax, gut, and abdomen. Flies were flash-frozen in liquid nitrogen, and tissues were dissected on dry ice. Metabolites were extracted using a cold methanol/water/chloroform extraction method. The extracts were analyzed using LC-MS/MS. Metabolite identification and quantification were performed using the MAVEN software. The data were normalized to the total protein content of each sample. Statistical analysis was performed using one-way ANOVA followed by Benjamini–Hochberg correction.

Constraint-based flux analysis
We performed constraint-based flux analysis using the COBRA Toolbox. We defined the HSD-muscle-GEM by constraining the glucose uptake rate and the fluxes of key glycolytic and TCA cycle reactions based on the experimental data. We used FBA, pFBA, and FVA-sampling to estimate the reaction fluxes. The sensitivity analysis was performed by systematically perturbing the baseline flux of each glycolytic reaction and quantifying the change in the total pyruvate consumption flux. The sensitivity coefficient (S) was calculated as the ratio of the relative change in the output flux to the relative change in the input flux.

13C-glucose isotope tracing
We performed in vivo 13C-glucose isotope tracing by feeding w1118 male flies with a HSD containing [U-13C]glucose for five days. Thoracic muscles were dissected, and metabolites were extracted and analyzed using LC-MS/MS. The mass isotopomer distribution (MID) was calculated for each metabolite. The fractional labeling was calculated as the ratio of the labeled isotopomers to the total pool of the metabolite.

Redox proteomics
We performed switch-tag redox proteomics to identify oxidized peptides. Proteins were extracted from thoracic muscles, and oxidized cysteine/methionine residues were labeled with a switch-tag reagent. The labeled peptides were analyzed using LC-MS/MS. The intensity of each peptide was normalized to the total protein intensity. Differentially oxidized peptides were identified using a threshold of |Log FC| > 0.5 and adjusted p-value < 0.05.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Noteworthy results:\n• Moon and colleagues describe an approach to apply genome scale metabolic models (GEMs) to individual tissues using\ntranscriptomic data from each tissue. They then integrate the muscle GEM with results of a labelled glucose tracer study to\nprofile reaction rates for glycolysis and the TCA cycle, identifying GAPDH as a rate limiting enzyme. In a separate study,\nthey integrate the muscle GEM with metabolomics data from flies fed either normal or high sugar diets, identifying sorbitol\nand ketone body metabolism as dysregulated by high sugar diet.\nSignificance and Novelty:\n• The authors make a valid point that it is crucial to understand the contributions of individual organs to systemic metabolism\nand efforts towards this aim are likely to be welcomed by the community. However, it is not made clear how the results of this\napproach differ to that which would be produced by traditional gene set enrichment. Clustering of the raw pseudobulk\ntranscriptomic data presented in figure 1b is remarkably similar to the clustering of tissues in their final results in figure 3f.\nThus, it is unclear how flux variability analysis adds to our understanding of metabolic activities beyond simply comparing\nenzyme expression. A similar concern exists for the results of the 13C flux analysis. These experiments are interesting,\nhowever, it seems that conclusions are largely apparent from looking at the raw fractional enrichment data.\n• The conclusion relating to GAPDH is not novel, as GAPDH is a well known rate limiting enzyme (see for example", "reviewer_pmid": "n25009227). Conclusions relating to ketone body metabolism are interesting, however, the provided results would benefit\ngreatly from further analysis of the mechanism by which this occurs, as was reported for oxidation of GAPDH.\nDoes the work support the conclusions and claims, or is additional evidence needed?\n• The message of many figures is unclear, especially early in the manuscript. Specifically figures 2b, 3b-e.\n• Figure 4 is interesting, however, the significance of this experiment, and how it relates to the conclusions of the paper\n(GAPDH, ketone metabolism) are unclear. The conclusion that metabolite abundance is reflected by the expression of\nrelated mRNA is not overly surprising. The authors do not quantify the accuracy of their predictions of metabolite flux match\nand so it would be interesting to see how well their predictions correlate with metabolite abundance (across all metabolites,\nrather than only for specific examples).\nIs the methodology sound? Does the work meet the expected standards in your field?\n• I have several concerns relating to the prediction of enzyme activity from mRNA abundance. Protein turnover rates are\nhighly protein-specific and likely vary between tissues. For example, the turnover of beta-oxidation proteins is greater in\nwhite than brown adipose tissue given increased exposure of brown adipose proteins to reactive oxygen species and\noxidative damage. This is not accounted for in the authors estimation of protein abundances, as data from all tissues is\nimputed from a whole organism reference dataset. This seems at odds with the motivation to increase the tissue-specificity\nof analyses.\n• Enzymatic activities of proteins are not directly linked to their abundances. For example, the activity of glycogen synthase\nis regulated exquisitely by protein phosphorylation and allostery. Many metabolic enzymes are regulated in this way, given\nthat they are activated by dynamic signalling processes (often phosphorylation) and changes in intracellular substrate levels,\nhowever, the approach taken here does not account for this. The authors argue that they confirm “previous reported\ndifferences in the estimated rates of metabolic processes” between tissues, however, as both sets of conclusions appear to\nbe based on very similar approaches, the results are not truly independent.", "reviewer_2": "The authors analyze tissue-specific metabolism in Drosophila melanogaster using metabolic network modeling,\nmetabolomics, and carbon tracing. They first construct metabolic network models for multiple tissues by integrating a Fly\nmetabolic model with tissue transcriptomics, obtained by aggregating data single cell RNA-sequencing. These networks are\nthen converted to quantitative models using estimated variables of enzyme kinetics together with flux balance analysis\n(FBA) and flux variability analysis (FVA). They then perform an experimental case study of high (HSD) vs normal (NSD)\nsugar diet in muscle. They show that glyceraldehyde 3-phosphate dehydrogenase (GAPDH) is very likely a rate limiting step\nof glycolysis in HSD and they show evidence that this is probably associated with an increased NADH/NAD ratio.\nGiven that the fly is an important model organism for studies of metabolism and metabolic diseases, bringing tissue-level\nresolution to the metabolic network analysis of this organism is significant and timely. However, I have concerns about the\nmodeling methods and the overall robustness of the models presented. Enzymatic variables are loosely defined and mixed\nwith FBA using an arbitrary objective function. It is not clear if these models are valid or relevant (in part because flux data is\nnot shared in supplementary tables). Furthermore, it is challenging to assert that the modeling results are validated by the\nHSD case study, as it appears that the experimental carbon tracing (and some other measurements such as NADH/NAD+),\nrather than the modeling, predominantly drives the discovery around GAPDH. Therefore, the contribution of modeling to the\nfindings is not convincing. Combined with a general deficiency in the quality and thoroughness of the methods section the\npaper as presented is not a robust advance.\nMajor Comments\n1) Vmax values are used to constrain FBA. Vmax is a multiplication of kcat and enzyme abundance, both of which are rough\napproximations. Thus Vmax values themselves are likely to have poor precision. Indeed, the analysis in the paper shows that\ngood predictions of just the enzyme abundance are only within an order of magnitude of experimental values. Such\nvariance raises doubts about the robustness of a quantitative model, particularly when multiple imprecise estimates are used\nwithin the same network. This could lead to the accumulation or amplification of errors, potentially skewing results\nsignificantly. Multiple specific examples would be needed to justify the usefulness of these constraints (such as predictions\nthat are true with these constraints but wrong without them).\nAlso related to this issue, some reported Vmax values are extremely large (more than a million). I have never seen such high\nnumbers for flux, especially with the unit of mmol/gdW/h, and I was not able to find anything close with a search. To address\nthese concerns, I recommend the authors provide detailed calculations for some of these Vmax values, starting from the\ninitial data to the final output. Furthermore, can the authors comment on what these high numbers mean for material\nconversion rates in a whole fly (can be calculated using the mass of a fly or its relevant organ to get an overall rate per hour,\nsecond, or day)? Are they realistic?\n2) The objective function of FBA adds up ATP, NAD and NADP demand fluxes together. What is the rationale behind this\nformulation? And how do these components interact within the model? Does the model sometimes pick one over the others\nand make some of these fluxes zero? Flux tables would be useful to provide answers.\nHere, ATP demand is a classical objective function and quantifies energy production capacity. But what the authors call\nNAD and NADP demand are reactions that convert NAD to NADH (or NADP to NADPH) freely, and therefore, do not\nquantify electron demand or capacity to extract reducing power from nutrients. Instead, they might inadvertently diminish the\nimportance of these fluxes, analogous to how a reversed ATP reaction (ADP + P -> ATP) would make the electron transport\nchain redundant. This raises questions about the effectiveness and accuracy of the simulations. Further explanations and\njustification are needed for this approach.\n3) The authors seem to relate high flux variability to metabolite concentrations, and they use the flux variability in a given\npathway to validate the enrichment of metabolites to that pathway in metabolomics data. However, it is difficult to link steady\nstate flux values to metabolite abundance in general. Also, variance of flux is an indicator of uncertainty; how would it be\nrelated to metabolite abundance?\nParticularly concerning is the formulation of pathway flux index (PFI). This formula involves multiplying the number of\nreactions in a pathway by a sum that itself multiplies the mean flux of each reaction by its flux range. Typically, one would\nexpect a normalization by the number of reactions rather than a straightforward multiplication, as this would prevent the\nindex from disproportionately escalating with the length of the pathway. Moreover, incorporating flux variability (vmax-vmin)\ninto the index suggests that the PFI might reflect the underlying uncertainty of the flux estimates more than the actual\npathway activity. Contrary to this, Figure 6a implies that PFI is low for pathways carrying low flux and vice versa; this is not\nnecessarily true based on the mathematical formulation. Even if this were true, we again have the problem with trying to link\nflux with metabolite abundance. The lack of a clear theoretical basis for the PFI formula suggests potential overfitting of the\nmetabolomics data.\nThe authors should first show the relationship between flux variability and metabolite enrichment or concentration by doing\nstatistical analysis, such as randomization. If there is indeed a significant relationship, then they could comment on why this\nwould be the case.\n4) The manuscript suffers from significant issues with clarity, most notably within the Methods section but also in other parts\nof the text. These issues affect the overall quality and robustness of the paper.\nMinor Comments\n1) Ln 124, Figs 1b and S1b: Mentioned clusters are not clear in these figures. S1b does not have a dendrogram to show\nclusters.\n2) Ln 133: \"The analysis revealed that germline cells, male accessory glands, and the fat body were able to complete\ngreater number of metabolic tasks than other tissues\"\nWhat is the significance of this and how do we know these chosen tasks truly represent the metabolic capabilities of a cell?\nIsn't this a biased analysis?\n3) Fig S1e: Color legend missing.\n4) Fig S1f: I believe y-axis is indicating a fraction, not number of transporters?\n5) Fig S2c: what is rho for this plot and how does it compare with S2b? It seems from the plot that there is no real\nimprovement in predicting protein abundance using transformed data instead of original mRNA. It is possible that I am not\ngetting it because S2d is not clear (please see below)\n6) Fig S2d: what is plotted as a histogram is not clear from legend or text.\n7) Fig. S2j: Please explain gene ratio in the legend.\n8) Ln 162: \"found that nearly 92 % of the maximum rates were within the default maximum values set in the Fruitfly1, with the\nmedian Vmax of 23.5 mmol/g-DCW/hr\"\nIt is not significant that predicted values are consistent with the default maxima. The default setting of upper and lower\nreaction boundaries (nearly all of which are 1000 and -1000, respectively), just indicate an arbitrarily high flux which you\ndon't expect to reach.\n9) Ln 178: \"certain tissues with higher Vmax were also capable of completing more metabolic tasks than those with lower\nVmax.\"\nWhat is the significance or biological relevance of this?\n10) In multiple figures, the font is too small to even read by zooming. The most notable example is Fig. S3.\n11) Lns 214-216: \"Further, for the protein degradation pathway, muscle cells showed higher flux variability, an observation\nthat agrees with the established function of muscle as a major site of amino acid breakdown.\"\nWhy should flux variability indicate activity? Doesn't it just mean we are uncertain about the flux level? Also, the exact range\nof the flux is important for interpretation. If zero is within the range, then we cannot even conclude if the reaction is predicted\nto carry flux. I would think reactions that are predicted to be active with high confidence are those with high mean flux and\nlow variability.\n12) Ln 226, Figure 3f: It is not clear how data points are grouped together. There are no apparent clusters.\n13) Ln 244: \"GABA serves as an inhibitor of neurotransmitters\"\nI think this is not accurate. GABA is an inhibitory neurotransmitter itself.\n14) Fig 5i: Not clear how FBA is done and what is shown. There are several problems here: (i) UB and LB of some reactions\nseem to be the standard value, 1000. So, does this mean vmax constraints on muscle are not used here? (ii) I think GAPD in other figures is labeled GAP in this one. (iii) Total flux into GAP is more than the flux out of it. Where does excess GAP go? (iv) What is imposed as a constraint and what is obtained? Please give examples to clarify the method, and provide the data. (v) Flux\nvariability is not shown on the HSD map. The authors can show the numbers as in the NSD pathway figure and still keep the\ncolors to indicate ratios.\n15) Ln 360-366: The metabolites discussed in the context of a specific butanoate metabolism pathway, e.g. acetoacetate,\nare connected to other pathways of the metabolic network, too. So how do we know that the change in proportions are due to\nenzymatic activities in this particular pathway?\n16) Ln 521: \"T is the mRNA nucleotides per cell.\"\nHow is this obtained and how can it be way less than 1 (1E-15)?\n17) In methods, equations are not numbered, making it difficult to follow and refer to.\n18) In equations of Lns 524 and 526, do mRNA and [mRNA] mean the same thing?\n19) Ln 535: \"Mw = 110 and mtot = 150 × 10−12 were used\"\nPlease indicate the unit of mtot. Also I think 1.5E-10 would be the correct form in scientific notation.\n20) Ln 539: PaxDb appears as abbreviation before this line.\n21) Ln 533, ATP demand equation: When this equation is used, what happens to the original maintenance requirement\nreaction?\n22) Ln 597: \"The objective function with (cid:1)(cid:2) of 1/3 is defined as\"\ncT is a vector, so this statement or the equation needs to be modified.\n23) Equation after Ln 609: What is the rationale for using this constraint?\n24) Ln 623: is alpha always taken as 0.9 in all analyses or is it varied as the equations suggest?\n25) Ln 633-639:\"The perturbation effect was defined by the summed reactions rates in glycolysis, normalized to the change\nof maximum reaction rates of perturbed enzyme rate.\"\nThis paragraph is not clear and this is the most confusing sentence. Is the quantity that is summed \"mean reaction rate\"? It is\nnot clear what this is normalized with. Please provide an example or specific equations?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "Overview:\nThis manuscript leverages single-nuclei gene expression data to infer tissue-specific genome-scale metabolic models in\nDrosophila. These are constraint-based models that enable prediction of the flux through all metabolic reactions for a given\ntissue. Additionally, kinetic information (maximum reaction rates, estimated from enzyme abundances and kcat values) is\nused to further constrain the predicted flux values. The authors compare features of the resulting metabolic networks (network size\nand structure, ability to perform defined metabolic tasks, max reaction rates, flux values). The models are used to predict the\nmetabolic features in the context of a high sugar diet, and metabolite profiling and heavy-labeled isotope tracing are used to\nqualitatively test some model predictions.\nOverall, the methods are well laid out, appropriate and rigorous. The results presented are clear. However, a primary issue with the paper is regarding significance/impact – it is not clear how the many results presented could be leveraged (by the authors or by others in the field). Specific comments are below.\nMajor issues:\n- A significant weakness is that the paper does not provide a compelling argument that the insights uncovered through\nmetabolic modeling are useful. Thus, it is not clear what the significance and potential impact of this work could be. The\ndiscussion mentions cancer cachexia as a future application, but in general, the implications of the metabolic differences\nbetween tissues are not presented.\n- The title emphasizes the application of studying the effects of high sugar diet; however, most of the paper focuses on\ncomparing the models (p. 6-10; figs 1-4). Only the last part is about high sugar diet (p. 11-13). And related to comment 1\nabove, the question of how we can use their results remains.\n- A central finding for predicting metabolic features with high sugar diet is that GAPDH is a bottleneck reaction (** note, in\nmultiple places, “GADPH” is incorrectly used). This is not a novel result, as GAPDH is widely known to be a limiting reaction\nin many conditions. This further emphasizes the need to explain why the reader should be excited about the results\npresented.\n- The authors should add a paragraph to the Discussion about limitations and assumptions that affect their findings. For\nexample, using estimates of protein abundances and limited availability of kcat values are expected to influence the results\nand subsequent conclusions.\nAdditional issues:\n- The “clusters” shown in Figure 1b are not distinguishable. A more quantitative/rigorous definition of the clusters is needed,\nas many appear to overlap.\n- The paper presents that tissues containing enzymes with high Vmax are able to complete more metabolic tasks (fig. 2b).\nThere appear to be several counter examples to this that should be discussed. Do the points corresponding to low Vmax\nrepresent a low ability to perform tasks? Is this meaningful?\n- A finding from the flux variability analysis is there are “similarities and differences across tissues” (line 225). This is a very\nvague sentence that is not insightful or actionable. Again, it aligns with the major critique that the results are not presented in\na way that highlights their utility.\n- It would be helpful to the reader to have a table that presents the major model predictions and whether each aligns with\nexperimental evidence (either from this paper or from published works). This would better support the authors’ statements\nthat model predictions are consistent with experimental findings.", "reviewer_5": "The manuscript by Sun Jin Moon et al. represents a comprehensive study on the use of tissue-specific GEMs to study the\neffect of a high glucose diet in drosophila. This is quite an impressive study in scope, and will be a good contribution and\nhighly cited study. I think the authors should feel proud of such comprehensive work. However, I think there are still some\nvery important items that need to be addressed regarding the model construction and validation process, along with\nstatistical analyses, which the authors should address, given that others in our community will be using these models, and\nthe models themselves are critical to the claims of the paper. Also there a many critical details that are glossed over, along\nwith community standards in tissue-specific model construction that need to be addressed. However, I hope to see these\nconcerns addressed since I’m excited about this fine work!\nMajor Comments:\n- We applaud the effort to further refine the drosophila reconstruction. These efforts for improving network curation are\nvaluable to the community.\n- Numerous algorithms exist for generating tissue specific models, and a number of decisions need to be made in\nprocessing the data. Such decisions can have a huge impact on the content that is obtained in tissue-specific models, and\nthus any interpretations of model analyses. Why was tINIT selected, as opposed to mCADRE, CORDA, or other high-quality\nmethods? What validations were done to believe that the models obtained by tINIT were valid? Often, community best-\npractices validate models by using gene knockouts, correlations to experimental metabolic fluxes (though this is confounded\nby the models used for fitting C-13 data), cell culture growth rates compared to model growth rates, other experimental\nvalidations, or comparison of model functions to comprehensive lists of known metabolic functions (not just citing a few\npapers that suggest a metabolic function could be active in a couple models, as done here for muscle metabolism). This is\ncritical since other studies (e.g., Opdam, Cell Systems, 2017; Machado, PLoS Comp Bio, 2014) showed some methods\nwere quite poor at making reliable models.\n- tINIT, as its name suggests, requires defined metabolic tasks. How were these defined? The methods just said they used a\nlist of 57 tasks in building the models. What were these tasks, how were they defined, and what is the support for each task?\n(see, e.g., Thiele, Nat Biotech, 2013; Blais Nat Commun, 2017; Richelle, PLoS Comp Bio, 2019). They later mentioned\nusing 256 human metabolism tasks to test model functions. How many of these tasks are relevant to drosophila, as opposed\nto just human?\n- Many tissue-specific model generation methods require one to define which genes are expressed in a given tissue. How\nwas this performed? Was it done by global or local thresholding? (Richelle, PLoS Comp Bio, 2019)\n- In figure 1c only 37/256 tasks are shown. It is stated that these were significant. How was significance defined and\ndetermined?\n- Figure 1d is unclear, and just reporting a Z-score provides limited insight. What % of tasks were functional for the different\ntissues? Do these conform with known tissue and cell-type specific tasks better than permuting tissue labels and comparing\ntasks with known tissue-specific functions? On a related topic, it seems that the authors provide some literature support for a\ncouple tissue-specific functions, but can they define tissue-specific functions prior to the analysis so that they can do a\nproper assessment of true/false positives/negatives to make sure their models are better than randomly-generated models?\n- How much does the Vmax for specific reactions vary across tissues? Do you see, e.g., glycolytic and lactogenic Vmax\nvalues to be higher for more metabolically active cells, like neurons?\n- In setting up simulations for FBA and related approaches, one must consider several things, including the objective\nfunction. Along these lines, multiple studies have benchmarked these, but I don’t see them referenced. For example, in lines\n183-186 the authors explain that all flux analyses in the paper are being conducted by setting the objective with a DM\nreaction that takes into consideration ATP, NAD and NADPH production. The authors state that this is a common metabolic\nobjective across tissues. However, they are not encountering other metabolic activities that the cell is subjected to and is\nrepresented by the biomass function. In this regard the authors should cite previous work that have carried a similar\napproach and also justify the addition of the NAD+ to the objective function as they only explained the incorporation of\nNADPH and ATP to the objective. This could be easily done by referencing the data and rationale in papers such as\nSavinell and Palsson, J Theor Bio, 1992; Schinn, Metabolic Engineering, 2021; Schuetz, Mol Syst Bio, 2007\n- A lot of flux spans are reported from flux variability analysis. These are merely descriptions of the models and parameters\nthat were inferred. I only see one single citation mentioning how the span on follicle and germ cells could be meaningful. Is\nthere any biological support for all of the differences discussed between lines 194-229?\n- Supplementary Fig 1 needs to incorporate labels that indicates the meaning of each color. Specifically panel A and E.\n- Supplementary Fig 3 needs to incorporate labels that indicates the meaning of each color. Specifically panel B, H and I.\n- In the metabolic profiling section starting with line 231, along with figures 4 and S4, the authors have generated a fantastic\ndataset and provided a mildly qualitative analysis, but with no statistical tests to see if there is a significant overlap between\nthe metabolomic data and their flux analyses. However, in theory, metabolic task analysis, which they did for figure 1,\nshould allow one to see, for each tissue, which metabolic tasks are more active (i.e., which metabolites could increase in\nactivity). I believe Sanjay Nigam’s lab at UCSD did such an analysis that he presented 3 or 4 years ago. Thus, it would provide\nmuch greater confidence in their tissue-specific models, if the authors here did a quantitative analysis comparing the metabolic\ndata to the tissue-specific metabolic task activities.\n- Line 342, the authors say there’s a decrease in PFI. What was the magnitude and statistical significance of the decreases?\nBordel, et al. PLoS Comp Bio, 2010 and Bar Evan, PNAS, 2010 from the Nielsen and Milo labs, respectively, presented\nways to quantify significance of changes in pathway fluxes using different sampling approaches. We’ve personally seen in\nmany instances increases and decreases in flux are not significant if using empirical p-values, but the approach is valuable\nfor finding the ones that really are.\n- The PFV subtracts the max and min fluxes for the reactions in a pathway, from what I understand. Shouldn’t this consider\nthe absolute value, since sometimes reactions are written backwards from the way they’re used, and thus would contribute a\nnegative flux value. On a related note, some reactions will be heavily skewed when doing sampling. Did you try using the\nmode instead for the mean flux for PFI?\nMinor Comments:\n- Line 136: Figure incorrectly referenced. Correct reference is Supplementary fig. 1e, instead of Supplementary fig. 1d.\n- Line 137: Is not clear what the authors are stating. I suggest rephrasing this sentence for better understanding.\n- Figure 2a should be incorporated as an equation in the manuscript and being removed from Figure 2.\n- As in Figure 2a, Supplementary Figure 2 equations in a and e should be addressed as equations rather than a panel in the\nfigure.\n- Figure legend of Supplementary Fig 4 panels h and I are wrongly referenced. I suggest reorganizing the figure.\n- Supplementary Fig 5, heatmaps corresponding to panels b and c should be increased in resolution for better readability.\n- Line 289: Add “Glyceraldehyde 3-phosphate” initials for better localization in the bar plot of Supplementary Fig.5a.", "reviewer_6": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Noteworthy results:\n• Moon and colleagues describe an approach to apply genome scale metabolic models\n(GEMs) to individual tissues using transcriptomic data from each tissue. They then\nintegrate the muscle GEM with results of a labelled glucose tracer study to\nprofile reaction rates for glycolysis and the TCA cycle, identifying GAPDH as a rate limiting\nenzyme. In a separate study, they integrate the muscle GEM with metabolomics data\nfrom flies fed either normal or high sugar diets, identifying sorbitol and ketone body\nmetabolism as dysregulated by high sugar diet.\nSignificance and Novelty:\n• “The authors make a valid point that it is crucial to understand the contributions of\nindividual organs to systemic metabolism and efforts towards this aim are likely to be\nwelcomed by the community. However, not made clear how the results of this approach\ndiffer to that which would be produced by traditional gene set enrichment. Clustering of\nthe raw pseudobulk transcriptomic data presented in figure 1b is remarkably similar to\nthe clustering of tissues in their final results in figure 3f. Thus, it is unclear how flux\nvariability analysis adds to our understanding of metabolic activities beyond simply\ncomparing enzyme expression.”\nWe thank the reviewer for the valuable comment. We agree that transcriptomic\nclustering or gene set enrichment analysis can also reveal metabolic similarities and\ndifferences across individual tissues. In the revised manuscript, we included\nadditional analyses by directly comparing clustering derived from either gene\nexpression and tissue-specific GEMs.\nThese comparisons have been added to the revised manuscript (Lines 126–133,\nFig. 1c, Supplementary Fig. 8–9). In brief, while we observed strong agreement for\ncertain tissues such as muscle, neurons, and fat body/oenocyte, agreement was\nlower for glial cells and gut tissues. These results suggest that GEM-based\nclustering, which incorporates stoichiometric constraints of biochemical reactions,\nreveals metabolic features not apparent from transcriptomics alone.\nAdditionally, we agree that the observed similarity between the clustering based on\nGEMs and the flux variability analysis-based clustering (previous Fig. 3f) may raise\nquestions about the added value of FVA-based analysis. However, this similarity is in\nfact expected, as the flux variability analyses were performed using the\nreconstructed tissue-specific GEMs. We acknowledge that we did not sufficiently\nclarify this point in the original submission. To improve clarity and focus the scope of\nthe current study, we have removed the tissue-wide FVA-based clustering results\nfrom the revised manuscript (original figure 3).\nAdditionally, we would like to emphasize that unlike traditional gene set enrichment\nanalysis, our tissue-specific GEMs enabled:\n Comparison of metabolic network structures across tissues, including reaction\ncontents, metabolite composition, subsystem coverage, and predicted\nmetabolic task performance (Fig 1, and supplementary fig1-6, and 10 - 13)\n Simulation of pseudo steady-state flux distributions (Fig 3)\n Identification of rate-controlling steps (e.g., GAPDH) via perturbation-based\nsensitivity analysis (Fig3h)\n Discovery of pathway-level metabolic perturbations using the pathway flux\nvariability index (PFVI) (Fig 6)\nWe recognize that these features were not clearly emphasized in the original\nmanuscript and have revised the text accordingly (Line 105 – 158).\n“A similar concern exists for the results of the 13C flux analysis. These experiments are\ninteresting, however, it seems that conclusions are largely apparent from looking at the\nraw fractional enrichment data.”\nWe thank the reviewer for this important point. In the revised manuscript, we\nclarified that the primary purpose of the ¹³C-glucose tracing experiments was to\nexperimentally validate model-predicted changes in glycolytic fluxes derived from\nthe muscle-GEM (line 252 – 283). To directly link model predictions with\nexperimental data, we included a correlation analysis between predicted fluxes\nand the fractional labeling of corresponding downstream metabolites (Fig 4e).\nWhile we agree that fractional enrichment itself can indicate glycolytic activities,\nour GEM-based analyses provided additional layers of metabolic analysis,\nincluding the identification of altered NADH-dependent reactions (Fig 3c-f,h) and\nadditional perturbed pathways beyond glycolysis (Fig 6a-c).\n• “The conclusion relating to GAPDH is not novel, as GAPDH is a well known rate\nlimiting enzyme (see for example"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have performed considerable additional experimentation in the revision process. This has addressed many/most of the comments raised. The manuscript is much improved.", "reviewer_2": "GENERAL COMMENTS\nThe authors made significant changes in their predictive simulations with the fly tissue GEMs and in the organization of the manuscript. The methodology has been simplified by replacing the previous modeling approach, which relied on kinetic parameters, with regular FBA. Despite this simplification, the predictive power is maintained (or possibly improved), consistent with earlier criticism that the added complexity from inaccurate kinetic parameters was not necessarily beneficial. The methodology is now clearer, aided by better manuscript organization and more straightforward language, in addition to the simplifications in modeling.\nHowever, some methodological concerns remain. A major issue is the treatment of reaction boundaries. As in many GEMs, including the human model used here, reaction flux bounds are uniformly set to +/-1000. This arbitrary value is not biologically meaningful but is meant to represent an effectively unbounded range, under the assumption that biologically relevant fluxes (e.g., those obtained with reasonable nutrient uptake rates and predicting reasonable biomass production rates) will not approach these limits. The authors’ approach, however, treats these arbitrary values as literal and perturbs them in some of their tests (see comments below), which raises technical concerns.\nIn addition, there are some remaining issues related to quality and clarity. Detailed comments are provided below for another round of revision.\nMAJOR COMMENTS\n1) The rationale for applying perturbations by altering reaction boundaries is unclear. A 20% reduction in the upper flux bound (or increase in lower bound) is arbitrary and does not affect all reactions equally. For example, if Reaction 1 and Reaction 2 have fluxes of 700 and 850 mmol/gDW/hr, lowering the upper bound from 1000 to 800 will have no impact on Reaction 1, but will constrain Reaction 2, as its reference flux exceeds the new boundary. This inconsistency undermines the interpretation of such perturbations. A more principled approach would be to apply perturbations directly to the flux values themselves; for example, by imposing a constraint that forces a reaction to carry no more than 80% of its unperturbed flux. While I acknowledge that the authors rely on flux sampling and do not define a single reference flux, it would be more appropriate to apply perturbations based on representative statistics of the sampled distribution (e.g., the mode or median), rather than manipulating arbitrary boundary values. This would avoid the conceptual problem of treating +/-1000 as meaningful limits and ensure a more uniform and interpretable perturbation strategy.\n2) Related to the issue above, some of the flux sampling distributions shown in Supplementary Figure 49 appear to cluster near the flux boundaries. This raises three concerning possibilities:\n(i) The default flux bounds (+/-1000), which are conventionally intended to approximate infinity, may be too restrictive relative to the feasible solution space. In this case, increasing the bounds (e.g., to +/-10000) might be appropriate to avoid artificial truncation of the distributions.\n(ii) The nutrient uptake rates may not be well constrained, allowing unrealistically high fluxes. As a result, key fluxes, particularly for the uptake of limiting nutrients, may saturate at the default boundary of 1000, leading to biologically implausible solutions.\n(iii) The high-magnitude fluxes may reflect the presence of internal flux loops. For example, one reaction (A + NADH <-> C + NAD) may carry a flux of +1000, while a parallel reaction (A + NADPH <-> C + NADP) carries -990, yielding a net flux of 10 while recycling redox cofactors through an internal loop. This scenario could be evaluated by performing parsimonious FBA (pFBA), which minimizes the total flux and suppresses such loops by favoring net conversions over unnecessary cycling.\nMore broadly, it would be informative to compare the predictive power of simple pFBA solutions to that of flux samples.\nGiven the potential ambiguity in sampled distributions, especially when they include artificial extremes or loops, pFBA could offer a more interpretable baseline for certain predictions.\n3) The rationale for using flux variability range as a multiplier in the Pathway Flux Variability Index (PFVI) is not well justified. The interpretation of this index appears problematic, particularly in cases where pathways are essential and consistently active. In such cases, reactions may exhibit narrow flux ranges that do not include zero (reflecting obligatory activity) but the PFVI will be suppressed due to the low variability. According to the current interpretation, this could misleadingly suggest low pathway activity.\nConversely, as illustrated in Supplementary Figure 50a, a reaction with a broad flux range that spans zero (e.g., the red curve) may be inactive in many feasible solutions, yet it receives a high PFVI score due to its variability. This undermines the biological interpretation of PFVI as a proxy for pathway activity.\nThis issue should be addressed either by adjusting the PFVI formulation or clarifying its intended meaning. One might also question whether similar qualitative predictions could be made without using variability at all; e.g., by relying solely on representative flux modes (or pFBA fluxes, as discussed in the previous comment). This would avoid conflating variability with activity and may yield more interpretable results.\n4) Despite the significant improvements, several figures and supplementary tables still suffer from quality and clarity issues. For example, Supplementary Figure 12 has unreadable column labels; Supplementary Data 2a includes #VALUE! errors in the Excel file; Supplementary Figure 15 and 38, as well as Supplementary Data 2b, also require attention. Please see additional issues detailed in the minor comments below.\nMINOR COMMENTS\n1) Figure 2f (pathway overlap between GEMs and metabolomics): It appears the authors use KEGG pathways as an intermediary to match GEM activity with metabolomics data. Why is this done indirectly through KEGG? Moreover, the method used to map GEM reactions to KEGG pathways is not clearly explained.\n2) Figure 3b: The numbers and parenthetical expressions in the panel are not explained in the legend. Please clarify.\n3) Page 9:\n“Specifically, glucose uptake rate was reduced to about 55%, while fluxes through glyceraldehyde-3-phosphate\ndehydrogenase…”\nThis sounds like a circular argument. If glucose uptake was manually reduced by constraining the upper bound, then observing a reduction in glycolytic flux may simply reflect that constraint. It may still be a valid observation if the flux remains well below the new boundary, but this should be clarified (See Major Comments 1 and 2).\n4) Page 9:\n“Overall, TCA cycle fluxes were decreased by approximately 50% in the HSD-muscle-GEM, consistent with decreased TCA\nflux observed in diabetic muscle.”\nThis does not appear consistent with Figure 3b, which shows an increase in upper TCA cycle flux. Additionally, please explain where excess α-ketoglutarate goes and where the additional malate comes from, as the alternate pathways consuming or producing these intermediates must account for any internal flux imbalance in the cycle.\n5) Page 10:\n“Comparison of these four reactions revealed approximately 2-fold increase in NADH generating fluxes and a 2-fold decrease in NAD+ regenerating fluxes (Fig. 3f and Supplementary Fig. 27) ... confirming the perturbed NADH metabolism.”\nThis reasoning is problematic. By mass balance, any NADH generation must be matched by NAD+ regeneration elsewhere in the network. How, then, is this imbalance quantified and connected to the NADH/NAD+ ratio? Are the authors referring to net electron flow (e.g., from NADH to other carriers)? Furthermore, focusing on just a couple reactions is difficult to justify when the rest of the network may compensate.\nTo assess differences in the network’s capacity to extract electrons, the authors could add a demand reaction (NADH -> NAD; in mass and charge balanced form and not in reverse direction as in the original submission) and maximize its flux across conditions.\n6) Figure 4c and the like: Normalizing everything with NSD eliminates valuable information. Why not normalize both NSD and HSD with their respective Glc uptake?\n7) Supp. Fig. 31: How can downstream metabolites show higher fractional labeling than glucose, the uniformly labeled carbon source used in the experiment?\n8) Fig. 6g: Black arrow means Flux is ND. What does it mean for a model prediction? Please clarify.\n9) Page 15: \"Consistent with this, our flux analysis predicted muscle could catabolize fructose and showed a significant decrease in the activity of Hex-A, a Drosophila hexokinase capable of phosphorylating fructose to fructose-6-phosphate75 in response to HSD. This aligned with the accumulation of upstream metabolic intermediates, such as sorbitol and trehalose.\"\nAn alternative hypothesis is that metabolite levels (rather than enzyme abundance) drive flux in high-sugar conditions. In that case, low enzyme expression (inferred from mRNA) does not necessarily indicate low flux. This complicates interpretation of perturbed states. Do the authors have data to rule out this alternative?\n10) Supp. Data 1c: What is the type and unit of the data presented in this table?\n11) Page 18, last paragraph: It is not clear what was achieved with getINITModel2().\n12) Pages 25.26: \"This formulation was used here as a general proxy for basal cell maintenance, including protein turnover, membrane lipid remodeling, and nucleotide recycling. These processes remain active and physiologically relevant to tissues like muscle70,95,96. Thus, rather than using energy and redox demands as sole objective functions, which could bias the solution space or do not recapitulate other core metabolic activities, we adopted this biomass objective function to capture broader metabolic requirements. We confirmed the flux through this reaction was approximately 1.1 mmol/gDCW/hr, which is orders of magnitude lower than the major glycolytic and TCA cycle fluxes, indicating that the objective function serves as a minimal biosynthetic demand (Supplementary Fig. 25, Supplementary Data 3h).”\nI do not understand the reasoning in flux comparisons. Biomass contains pools such as TAG. For example, fatty acid synthesis required for biomass may translate into a very large energetic demand, depending on the biomass composition.\n13) Page 28: “For visualization of the predicted structures we utilized ChimeraX (PMID: 37774136).”\nThere appears to be a formatting issue with this reference.\n14) Page 30: “To evaluate the statistical significance of pathway-level flux changes between NSD and HSD conditions, we performed a two-sample Z-test using flux distributions obtained from sampling analysis, following an approach described in a previous study105.”\nHowever, the bibliography does not contain 105 references. Please correct this.\n15) Page 30: “where (cid:3)̅ represents the mean pathway flux”\nIt is unclear how this average is computed. For example, in glycolysis, flux doubles mid-pathway due to carbon splitting (from 6C to 3C intermediates). How is this handled when averaging across a pathway?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "Overview:\nThis manuscript leverages single-nuclei gene expression data to infer tissue-specific genome-scale metabolic models in Drosophila. These are constraint-based models that enable prediction of the flux through all metabolic reactions for a given tissue. The authors first briefly compare features of the resulting metabolic networks (network size and structure, ability to perform defined metabolic tasks, representation of metabolic subsystems). The muscle-specific model is used to explore the effects in the context of a high sugar diet, and metabolite profiling and heavy-labeled isotope tracing are used to validate some model findings.\nOverall, the revised manuscript is much improved. The paper is more streamlined and better demonstrates the utility of tissue-specific GEMs. The methods are well laid out and explained. The results presented are clear, and the title and introduction align with the main message of the paper. Only a few smaller issues remain.\nSpecific comments:\n- The authors have used experimental methods to confirm several model predictions. This is laudable and nicely bridges communities using computational and experimental methods to study metabolism. However, the authors have not sufficiently emphasized what computational modeling using GEMs uniquely offers above and beyond experimental approaches. This should be better highlighted. This is an important point to address.\n- More justification is needed for the perturbation-based sensitivity analysis. In the first submission, the upper flux bound was decreased by 10%. Here, 20% is used. Either seems arbitrary and should be justified why that perturbation is meaningful. I could not find any support in the paper cited (ref. 99).\n- The number of supplementary figures grew by nearly 10-fold, from 6 to 56! It is the responsibility of the authors to provide a well-curated story – this task has not been met if there are 56 supplementary figures. I believe the authors must reduce this number to only include those that are absolutely needed to help tell the story (rather than adding figures that are “nice to have”). This is a critical point to address.\nEvaluating the response to the points raised by Reviewer 5:\nIn my reading, most of the comments from Reviewer 5 were addressed appropriately. I believe major points related to comparing to experimental data/observations were useful. Additionally, the changes to the manuscript to remove flux distributions and Vmax calculations made some of the reviewer points no longer relevant.\nList of metabolic tasks. This response is still unclear. The authors should mention all essential tasks (57 from human metabolic task list and 219 metabolic functions) in a single place within the methods. As it stands, these two lists are described in different locations in the manuscript. Furthermore, it is not clear whether the 57 is part of the 219. This requires greater clarity from the authors.\nExperimental validation. The authors can go further with their validation, as comparing enriched pathways (data) with subsystem coverage (model) does not consider if the reactions present carry flux or not, it just checks if the subsystem is there or not. I believe it is important to perform the same model predicted enrichment after flux values are calculated. It is not clear that the predictions and FEA would agree.", "reviewer_1_attached": "Reviewer #1 (Remarks to the Author):\nNoteworthy results:\n• Moon and colleagues describe an approach to apply genome scale metabolic models\n(GEMs) to individual tissues using transcriptomic data from each tissue. They then\nintegrate the muscle GEM with results of a labelled glucose tracer study to\nprofile reaction rates for glycolysis and the TCA cycle, identifying GAPDH as a rate limiting\nenzyme. In a separate study, they integrate the muscle GEM with metabolomics data\nfrom flies fed either normal or high sugar diets, identifying sorbitol and ketone body\nmetabolism as dysregulated by high sugar diet.\nSignificance and Novelty:\n• “The authors make a valid point that it is crucial to understand the contributions of\nindividual organs to systemic metabolism and efforts towards this aim are likely to be\nwelcomed by the community. However, not made clear how the results of this approach\ndiffer to that which would be produced by traditional gene set enrichment. Clustering of\nthe raw pseudobulk transcriptomic data presented in figure 1b is remarkably similar to\nthe clustering of tissues in their final results in figure 3f. Thus, it is unclear how flux\nvariability analysis adds to our understanding of metabolic activities beyond simply\ncomparing enzyme expression.”\nWe thank the reviewer for the valuable comment. We agree that transcriptomic\nclustering or gene set enrichment analysis can also reveal metabolic similarities and\ndifferences across individual tissues. In the revised manuscript, we included\nadditional analyses by directly comparing clustering derived from either gene\nexpression and tissue-specific GEMs.\nThese comparisons have been added to the revised manuscript (Lines 126–133,\nFig. 1c, Supplementary Fig. 8–9). In brief, while we observed strong agreement for\ncertain tissues such as muscle, neurons, and fat body/oenocyte, agreement was\nlower for glial cells and gut tissues. These results suggest that GEM-based\nclustering, which incorporates stoichiometric constraints of biochemical reactions,\nreveals metabolic features not apparent from transcriptomics alone.\nAdditionally, we agree that the observed similarity between the clustering based on\nGEMs and the flux variability analysis-based clustering (previous Fig. 3f) may raise\nquestions about the added value of FVA-based analysis. However, this similarity is in\nfact expected, as the flux variability analyses were performed using the\nreconstructed tissue-specific GEMs. We acknowledge that we did not sufficiently\nclarify this point in the original submission. To improve clarity and focus the scope of\nthe current study, we have removed the tissue-wide FVA-based clustering results\nfrom the revised manuscript (original figure 3).\nAdditionally, we would like to emphasize that unlike traditional gene set enrichment\nanalysis, our tissue-specific GEMs enabled:\n Comparison of metabolic network structures across tissues, including reaction\ncontents, metabolite composition, subsystem coverage, and predicted\nmetabolic task performance (Fig 1, and supplementary fig1-6, and 10 - 13)\n Simulation of pseudo steady-state flux distributions (Fig 3)\n Identification of rate-controlling steps (e.g., GAPDH) via perturbation-based\nsensitivity analysis (Fig3h)\n Discovery of pathway-level metabolic perturbations using the pathway flux\nvariability index (PFVI) (Fig 6)\nWe recognize that these features were not clearly emphasized in the original\nmanuscript and have revised the text accordingly (Line 105 – 158).\n“A similar concern exists for the results of the 13C flux analysis. These experiments are\ninteresting, however, it seems that conclusions are largely apparent from looking at the\nraw fractional enrichment data.”\nWe thank the reviewer for this important point. In the revised manuscript, we\nclarified that the primary purpose of the ¹³C-glucose tracing experiments was to\nexperimentally validate model-predicted changes in glycolytic fluxes derived from\nthe muscle-GEM (line 252 – 283). To directly link model predictions with\nexperimental data, we included a correlation analysis between predicted fluxes\nand the fractional labeling of corresponding downstream metabolites (Fig 4e).\nWhile we agree that fractional enrichment itself can indicate glycolytic activities,\nour GEM-based analyses provided additional layers of metabolic analysis,\nincluding the identification of altered NADH-dependent reactions (Fig 3c-f,h) and\nadditional perturbed pathways beyond glycolysis (Fig 6a-c).\n• “The conclusion relating to GAPDH is not novel, as GAPDH is a well known rate\nlimiting enzyme (see for example"}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors addressed all of my previous comments and I agree that the manuscript has improved further. Below are two minor issues with the updated sections:\n1. NADH -> NAD conversion reaction: I could not find this in supplementary tables. The authors should make sure this reaction is mass-balanced using a proton, so that the observed effect is not due to proton consumption (but not charge-balanced, since the reaction represents electron demand).\n2. Methods, \"Quantification of pathway overlap using weighted Jaccard index\": This section is rewritten but is still unclear. In particular, the explanation on wB(i) (normalized weight for pathway i) is confusing. For example, in the equation (cid:4)(cid:5)((cid:6)) = 1/|(cid:5)|, there is no i on the right hand side. What does this mean? Also, why is B in absolute value brackets if it is a number of enriched pathways, which should always be positive?", "reviewer_4": "The authors have addressed my previous comments. I particularly appreciate the changes made to address concerns related to the perturbation analysis (which another reviewer also raised).\nI do not have any further comments to address."}, "author_response": "We thank the reviewer for pointing out this important detail. In fact, we used the NADH demand reaction that was mass balanced by including a proton, represented as NADH -> NAD+ + H+. Specifically, the objective function of NADH demand reactions across compartments was defined as:\nNADH[c] + NADH[m] + NADH[p] → NAD+[c] + H+[c] + NAD+[m] + H+[m] + NAD+[p] + H+[p]\n, where [c], [m], and [p] refer to the cytosol, mitochondria, and peroxisome, respectively.\nWe have now clarified this in the Methods section (“Evaluation of maximum NADH production capacity”) and updated the reaction equation by including the proton in Fig. 3d. The corresponding optimal flux values of NADH -> NAD + H is shown as bar graph in Figure 3e and Source Data Fig 3e.\n\nWe thank the reviewer for this valuable comment. We agree that this method required more detail and now expanded the method (Methods: Quantification of pathway overlap using weighted Jaccard index). Specifically, we provided details of 𝑤 (𝑖) as follows:\n𝐵\n𝑤 (𝑖) refers to the normalized weights assigned to GEM-based enriched pathways. Since\n𝐵\nthe enriched pathways identified from GEM do not have statistical significance values,\neach pathway i in set B is assigned an equal weight, calculated as 𝑤 (𝑖) = 1/𝑛 , where\n𝐵 𝐵\n𝑛 is the number of enriched pathways in that tissue-specific GEM (set B). For any\n𝐵\npathway i that is not part of set B (e.g., 𝑖 ∈ 𝐴∪𝐵 and 𝑖 ∉ 𝐵), 𝑤 (i) = 0. The sum of all\n𝐵\n𝑤 (i) equals 1.\n𝐵\nWe originally used |B| as standard set notation for the number of elements in set B, but\nrecognize that this may be misleading and have replaced it with the clearer notation 𝑛 .\n𝐵\nWe believe this update resolves the ambiguity and improves clarity."}]

Decision letter