In Silico Metabolic Modeling Reveals Potential Muscle Metabolite Markers of Insulin Resistance in Newly Diagnosed Diabetic Patients

Background Type 2 diabetes mellitus (T2DM) is a challenging and globally ubiquitous metabolic disease caused by insulin resistance. Skeletal muscle is the major insulin-sensitive tissue that plays a great role in blood sugar homeostasis. Dysfunction of muscle metabolism is implicated in the disturbance of glucose hemostasis and the development of insulin resistance and T2DM. Here, we attempted to nd metabolic dysregulations that are associated with the onset of T2DM. Besides, metabolite markers of T2DM were explored. Methods We reconstructed a human muscle-specic metabolic model and applied it to perform metabolic analysis in newly diagnosed T2DM patients. We investigated the metabolism reprogramming by using two topology-based and constraint-based approach. Moreover, we applied a machine learning method to predict potential metabolite markers of insulin resistance in muscle. Results Our results showed that metabolic alterations have occurred in carbohydrate, fatty acids, lipids, amino acids, and inositol phosphate metabolisms as well as pathways implicated in building extracellular matrix (ECM). Also, dysregulation of coenzyme Q10 metabolism was observed. Moreover, 13 exchange metabolites were predicted as the potential metabolite markers of insulin resistance in skeletal muscle. The eciency of these markers in detecting insulin-resistant muscle was validated using a separate muscle gene expression data from another diabetes-related study. Conclusion In this study, the most updated muscle-specic metabolic model was generated and successfully was validated. This model was used for the investigation of metabolic disturbances at the onset of T2DM. Our results indicated the signicance of ECM metabolites in insulin resistance, and reinforce the role of coenzyme Q10 as a candidate for further research in insulin resistance and T2DM treatment. The model is freely available and can be used for other muscle metabolic studies. We also predicted metabolite markers of insulin resistance in the skeletal muscle, which can be considered for further empirical investigations.


Introduction
Type 2 diabetes mellitus is a complex progressive metabolic disease with a heterogeneous etiology. The disease has a high global prevalence and is among the top 10 causes of death worldwide. According to the latest international diabetes federation report, the global prevalence of diabetes in 2019 is estimated to be 463 million people. In addition, 374 million people are at greater risk of developing T2DM [1]. Therefore, this disease has become a major concern in the eld of global health.
Insulin resistance is the main pathophysiologic feature in T2DM. The key insulin-sensitive tissues include skeletal muscle, liver, and adipose tissues. Among these tissues, skeletal muscle is responsible for the clearance of over 75% of glucose from the bloodstream thus plays a great role in lowering the blood glucose level [2]. Dysfunction of muscle metabolism is associated with disturbance in glucose hemostasis and implicated in the development of obesity, insulin resistance, and T2DM. Therefore, the study of the molecular mechanisms underlying insulin resistance in skeletal muscle will improve our understanding and treatments for T2DM.
In complex and multifactorial diseases such as diabetes, the use of systems biology provides a holistic view to study affected subsystems related to the disease in a cell. Meanwhile, genome-scale metabolic modeling serves as a context for in silico exploration of metabolism through a systems biology approach. Since there is not particularly a one-to-one relationship between transcription and translation, the study of gene expression data alone does not provide an accurate understanding of cellular metabolism while metabolic network simulation gives an in-depth insight into the molecular mechanisms involved in cellular processes. Generic metabolic networks are reconstructed based on the information encoded in the genome. This information includes all possible reactions that can occur in all types of human cells. By having gene expression data from a particular cell or disease state and integrating it into the generic model, active reactions are identi ed and a context-speci c model can be reconstructed.
Using genome-scale metabolic models (GEMs), some efforts have been made to study the molecular mechanisms contributing to muscle insulin resistance, so far. Bordbar et.al. have reconstructed a multitissue type genome-scale metabolic model and integrated gene expression data of obese and type 2 obese gastric bypass patients to study metabolic activity between these states [3]. Nogiec et.al. have applied ux balance analysis to model muscle insulin resistance metabolism in fasted to fed states. In addition, they tried to identify reactions that reproduce key features of insulin resistance by perturbing the muscle metabolic network [4]. Varemo et.al. have reconstructed a myocyte metabolic network and have utilized it to nd reporter metabolites in type 2 diabetes mellitus [5]. However, so far there is not any study demonstrated the metabolic disturbance at the onset of diabetes. Since T2DM is a progressive disease, system-level identi cation of disorders in newly diagnosed T2DM patients who are at the early stages of the disease process can provide signi cant information about the underlying onset disorders in the disease and may somewhat clarify the causes of insulin resistance.
In this study, we attempted to nd muscle metabolism dysregulations at the onset of type 2 diabetes mellitus. We used muscle gene expression data from newly diagnosed diabetic patients that their blood glucose levels were slightly higher than the diabetic diagnostic criterion, and had not taken any medications [6]. We rst reconstructed the functional muscle-speci c metabolic model and then employed it for our analyses. We applied a topology-based approach for the identi cation of reporter metabolites, which are associated with the dysregulated genes. Then using a constraint-based approach, the changes in muscle metabolic capabilities were studied. Moreover, potential metabolite markers related to the muscle insulin resistance were predicted using a machine learning approach. Finally, we validated these markers with gene expression data from a separate study [7]. The overall work ow of this study is shown in Fig. 1.

Data
We searched transcriptome databases for skeletal muscle data from healthy and T2DM individuals. We found the most comprehensive dataset from the dbGaP database with the accession code phs001068.v1.p1. The gene expression data were obtained from vastus lateralis muscle using highthroughput RNA-Seq technology. The data comprises of 91 healthy and 63 newly diagnosed subjects [6]. This transcriptome data were downloaded and were used for our analyses.
For validation of metabolite markers the muscle gene expression data of healthy and diabetic subjects were taken from the GEO repository with the accession numbers GSE63887 and GSE81965 [7]. This data contains RNA-Seq samples from 24 individuals in four equal following groups: normal glucose tolerant (NGT)/non-obese, NGT/obese, T2DM/non-obese, and T2DM/obese.

Differential gene expression
The Bioconductor R package DESeq2 was used for differential gene expression analysis [8]. A preprocessing of data was applied involving pre-ltering of genes whose expression level was below a minimum cutoff level (< 5 read counts in less than 25 percentages of samples). Between-samples normalization was done according to the DESeq2 manual. DEGs were found based on a negative binomial distribution. Multiple testing correction based on Benjamini-Hochberg procedure was applied and a false discovery rate lower than 0.1 was considered as the differentially expressed criterion.

Muscle-speci c metabolic model reconstruction
The muscle metabolic model was reconstructed based on the HMR2 model as the generic model and gene expression data of healthy muscles [9,10]. The Gene-protein-reaction relationships in HMR2 were corrected using the Metabolic Atlas database [11]. For the preprocessing, low expressed genes (< 5 read counts in less than 25 percentages of samples) were ltered and between-sample normalization according to the DESeq2 normalization method with gene length adjustment and log2 (gene expression + 1) transformation were applied [8]. Moreover, the myocyte biomass reaction from the Bordbar muscle metabolic model was added to our model [3]. E-Flux algorithm was employed to integrate geneexpression data to HMR2 model and to reconstruct the context-speci c metabolic model [12]. To build a fully functional model we removed dead-end reactions. Figure 5 shows the work ow of this section. This model is freely available at https://github.com/Maryamkhn/Muscle_MetNet_CBB.

Reporter metabolite analysis
The Bioconductor R package piano was used for gene set analysis with the identi cation of reporter metabolites [13]. To determine the metabolite gene sets, the metabolite-reaction-gene associations of muscle metabolic model were used and gene IDs were assigned to each related metabolites. Adjusted pvalues with the log2-fold changes obtained from DESeq2 were used here.
Personalized metabolic models and ux variability analysis Personalized GEMs were reconstructed by integrating gene expression data from each individual with the muscle metabolic model. Pre-processed data and E-Flux algorithm were used here. Media conditions were set by body uid metabolites [14]. We used maximizing ux through the production of mitochondrial ATP as the objective function. In addition, to ensure having alive cell, we set the lower bound of biomass reaction to 80% of the maximum potential of the healthy model to produce biomass [15]. To obtain the minimum (min) and maximum (max) possible ux through each reaction in each personalized model, FVA was employed using Cobra Toolbox [16]. The two-sample t-test was applied on the min and max uxes to nd perturbed reactions between T2DM and healthy states. To apply multiple testing correction, the Benjamini-Hochberg procedure was used and reactions with false discovery rate less than 0.1 were regarded as signi cant perturbed reactions. Signi cant reactions were used to get perturbed pathways using ux enrichment analysis in the Cobra toolbox.
Classi cation and feature selection A SVM classi er with a polynomial kernel of order 2 was used for classi cation. The classi er was evaluated by 10-fold strati ed cross-validation and analysis of accuracy, sensitivity, and speci city were reported.
For feature selection, a feature selection method, based on a hybrid of genetic algorithm and support vector machine was employed. In the GA, the feature subset is encoded as the candidate solution on a chromosome-like structure. The population is formed from a set of chromosomes that mutation and crossing over in it occur to generate the next generation. A tness score is calculated for each chromosome representing adaptation of it to the environment and better feature subsets have more chance of reproducing the next generation. Creating the next generations will be continued until a stopping criterion is satis ed.
Here, a binary genetic algorithm was employed that the binary values of 1 and 0 represent the presence or absence of a speci c feature at the particular chromosome, respectively. The chromosome length was set to the number of features. The population size and the maximal number of generations were set to 300 chromosomes and 100 generations, respectively. We used the accuracy of the SVM classi er as the tness score. The genetic algorithm was stopped if the tness score was reached an accuracy higher than 90 percent or the maximum number of generations was reproduced.

Human Muscle-speci c Metabolic Model Reconstruction
To reconstruct the muscle-speci c metabolic network, the generic metabolic model and gene expression data of muscle were needed. For model reconstruction and analysis of metabolic disturbance in T2DM, we looked for the most comprehensive gene expression data from the muscle of healthy and diabetic individuals who had the following conditions: 1) newly diagnosed with T2DM. 2) Taking no medication for diabetes. 3) Do not have any other disease affecting the analysis. This resulted in data with the accession code phs001068.v1.p1 in the dbGaP database. Diabetic patients in this sample have been newly diagnosed as T2DM patients, and in most of them, the fasting blood glucose concentration was around 7 mmol/l, thus they were at the onset of diabetes. Table 1 shows subjects characteristics of the samples used in this study. Investigation of metabolism reprogramming in these patients who are at the onset of T2DM has signi cant value in identifying the underlying disorders in T2DM. We used gene expression data of healthy individuals in this sample for the reconstruction of the muscle-speci c metabolic model. The muscle-speci c metabolic model was reconstructed by integrating muscle gene expression data from healthy group into the Human Metabolic Reaction 2 (HMR2) [10] as the generic metabolic model.
Muscle-speci c biomass reaction was added to the model in order to have a functional and alive muscle model. The model comprises of 3585 metabolites, 5740 reactions and, 3627 genes. Table 2 shows the characteristics of the muscle-speci c metabolic model in comparison with HMR2. For quality controls and validation, we rst checked for the biomass production to assess that all the required reactions to have alive model are satis ed. Then, 56 metabolic tasks known to occur in all types of human cells were checked; these metabolic functions can be found in [17]. Moreover, for muscle-speci c validation, like Bordbar et.al. [3], we tested the ability of the model to produce ATP from several different sources like as glucose, fatty acids, glycogen branched-chain amino acids, and ketone bodies. The model successfully passed all these tests. The reconstructed muscle-speci c metabolic model was employed for further analyses. Reporter metabolites associated with the onset of diabetes were identi ed by topology-based analysis In this analysis, the reconstructed muscle-speci c metabolic network was employed as the scaffold to explore the effects of transcriptional reprogramming in the context of metabolism. Using this approach, it is possible to identify so-called reporter metabolites that are those metabolites associated with the most signi cant transcriptional changes [18]. Gene sets associated with each metabolite can be identi ed from a metabolic network through gene-reaction-metabolite relation. Here, the Piano R package was used to perform reporter metabolites analysis. Results demonstrated that from 3310 metabolites having associated genes, 275 ones were identi ed as reporter metabolites with adjusted p-value < 0.05 in at least one directional/non-directional class. Most of the reporter metabolites were affected by down-regulated genes. Also, most of them were implicated in mitochondrial functions such as the production of mitochondrial ATP, activated fatty acids, pyruvate, ubiquinol, ubiquinone, AKG, succinyl-CoA, NADH, Acetyl-CoA, lipoamide, 5,10-methenyl-THF, and 5,10-methylene-THF. In addition, some metabolites including ceramide, tyrosine, and tryptophan were affected by up-regulated genes. The reporter metabolites were mainly involved in glycolysis, tricarboxylic acid, fatty acids, and amino acids metabolism. The complete list of the reporter metabolites can be found in Additional File1 Table S1.
Constraint-based analysis of metabolism at the onset of T2DM reveals reprogramming of pathways related to extracellular matrix We attempted to capture cell metabolism at the system-level by using constraint-based metabolic modeling. In this approach, the metabolic model is converted to the mathematical format, and several constraints (e.g. adjusting lower and upper bound of ux in each reaction) were imposed on the model. Combining high-throughput transcriptional data and constraint-based models provides an opportunity to infer cellular metabolism at the relevant physiological state like diabetes [19]. Flux variability analysis (FVA) is a well-known constraint-based approach, which can provide allowable ux states in the metabolic model through linear programming based strategy [20].
In order to study metabolic capability changes in muscle between healthy and newly diagnosed diabetic patients, personalized GEMs were reconstructed and FVA was applied. Signi cant perturbed reactions between the two states identi ed using the two-sample t-test. This analysis resulted in perturbed reactions involving in several pathways such as glucose, amino acid, and lipid metabolisms. The perturbations in the glucose to glucose-6-phosphate conversion reaction and the glucose exchange reaction were also observed. Many reactions contributing to fatty acids and lipid metabolism were dysregulated. Also, changes in uxes of several reactions implicated in inositol phosphate, chondroitin, and heparin sulfate metabolism were observed. To nd associated pathways with these perturbed reactions, ux enrichment analysis was employed. These signi cantly perturbed pathways are shown in Table 3.

Potential Metabolite Markers Prediction
The genome-scale model of human metabolism provides an opportunity to predict potential bio uid markers associated with the diseases. The overall approach involves identifying exchange metabolites that their ux range in disease state is shifted compared to healthy ones. The successful performance of this method has already been evaluated and validated [21]. Here, we applied this method to nd exchange metabolite markers of insulin resistance in muscle. Two generic metabolic models of healthy and diabetic muscles were reconstructed using average gene expression data of each group and E-Flux method; Then, FVA analysis was applied to them. The exchange reactions in which the ux interval was shifted relative to the normal state were selected (Fig. 2). This approach resulted in 32 exchange reactions. In order to evaluate the discriminative capability of these reactions at the individual level, a machine learning approach was employed. To do this, personalized metabolic models of healthy and diabetic individuals were reconstructed and FVA analysis was applied. The generated minimum and maximum uxes of these reactions based on FVA analysis were used as the features for classi cation with support vector machine (SVM). This evaluation resulted in 72.22% accuracy, 81.11% speci city, and 60.32 sensitivity. We found that these exchange metabolites altered at the individuals' level.
To achieve a fewer and better feature subset with improvement in classi cation accuracy, a wrapper feature selection method was employed that comprises a combination of genetic algorithm (GA) and SVM as the feature selector and classi er, respectively. Several subsets of features with which SVM classi er can discriminate diabetic patients from healthy ones with approximately 90 percentages accuracy were found. This GA-SVM procedure was iterated 100 times resulting in 100 features subsets; then, the observation frequency of each feature in these iterations calculated. These frequencies were considered as indicative of the importance of features and features with at least 80% frequency regarded as top-ranked features. The top-ranked features were related to 13 reactions exchanging glucose, alanine, aspartate, sodium, hyaluronate, galactose, GQ1b, acetaldehyde, deoxyribose, 4-OH-Estradiol, hypoxanthine, retinoate, and methylglyoxal. Subsequently, the SVM classi cation performance with the top-ranked features was assessed. Our analysis revealed that using these top-ranked features improved classi cation accuracy to 81.04% with 73.02% sensitivity and 86.67% speci city. Since, the classi er performance will vary depending on which samples are assigned to the training set and which ones to the test set, we repeated the 10-fold cross-validation 100 times and evaluated the SVM performance. The boxplot of this evaluation is shown in Fig. 3.
For further validation, we assessed the performance of these top-ranked features with data obtained from another study [7]. This data comprises of RNA-Seq samples from vastus lateralis muscle of 24 participants, which is divided into four following groups: normal glucose tolerant (NGT)/non-obese, NGT/obese, T2DM/non-obese, and T2DM/obese. Preprocessing of the read counts was applied. Then personalized metabolic models were reconstructed and FVA analysis was performed. Associated uxes of top-ranked features were obtained and were used as the features for SVM classi er.
We divided these people into NGT and diabetic groups regardless of the obesity state and the generated uxes of marker reactions in these 24 participants were used for validation with SVM. This analysis resulted in 49.46% accuracy, 42.25% speci city, and 56.67% sensitivity, which is remarkably low, probably due to the presence of normoglycemic obese people in the healthy group. Although the blood glucose level in these individuals is in the normal range, high levels of fasting blood insulin level (60 ± 33 pmol/l in NGT/obese vs 25 ± 7 pmol/l in NGT/non-obese) indicates insulin resistance in this group, which is compensated by more insulin secretion. Thus, we decided to remove data from NGT/obese individuals and re-evaluated markers in this set. The result signi cantly altered with the 78.50% accuracy, 70.17% speci city, and 82.67% sensitivity. Therefore, the metabolites markers are related to insulin-resistant metabolism and can successfully discriminate insulin-resistant individuals from healthy ones.

Discussions
Here, we attempted to unravel the metabolism reprogramming in skeletal muscle of newly diagnosed diabetic patients. Since diabetes is a metabolic disease, our goal was to identify the metabolic disorders that occur at the onset of the disease. We used the gene expression data from the most comprehensive human skeletal muscle transcriptome study to date. The diabetic sample consists of participants who have newly been diagnosed with diabetes (with an average blood glucose concentration of ~ 7.2 mmol/L) and had not taken any medications. For understanding metabolism alterations in T2DM, we rst reconstructed our functional muscle-speci c metabolic model because we could not utilize Bordbar [3] or Nogiec [4] muscle metabolic models due to the being small-scale. Varemo model [5] has been the most complete muscle metabolic model to date, but because of the wrong gene-protein-reaction (GPR) association and lack of biomass production reaction, it is not suitable for simulation. Thus, we reconstructed the most updated skeletal muscle metabolic model. The model quality was successfully validated for several pathways, and metabolites productions as well as biomass production. This model is freely available and can be used for other muscle metabolic analysis.
We employed the reconstructed model for investigation of dysmetabolism in T2DM by using topologybased and constraint-based analyses. By topology-based analysis, we determined those metabolites that are involved in reactions in which the associated genes are signi cantly dysregulated. We also applied constraint-based analysis to quantitatively compare metabolic capabilities between healthy and newly diagnosed diabetic patients. In these analyses ubiquinone and ubiquinol the oxidized and reduced forms of coenzyme Q10, respectively, were identi ed as reporter metabolites. These metabolites also participate in some perturbed reactions. Coenzyme Q10, a critical component of oxidative phosphorylation, is produced or consumed in the electron transport of the mitochondrial respiratory chain and possesses antioxidant and anti-in ammatory properties. Pieces of evidence have revealed that supplementations of coenzyme Q10 in T2DM patients, preserve mitochondrial function, reduce oxidative stress, and improve glucose tolerance [22,23]. In addition, the administration of coenzyme Q10 in patients with prediabetes has alleviated the progression from prediabetes to diabetes [24]. Therefore, this metabolite should receive more attention in the treatment of diabetes.
We also found that metabolic alterations occurred in carbohydrates, fatty acids, lipids, and amino acids metabolisms. Dysregulation of Branched-chain amino acids (BCAAs) metabolism results in the serine phosphorylation of insulin receptor substrates and subsequent uncoupling of insulin signaling [25]. Moreover, perturbations in metabolism of inositol phosphate, keratin, chondroitin, and heparan sulfate were observed. Inositol mediates insulin signal transduction, associated with glucose uptake and plays an important role in oxidative stress and in ammation. Administration of inositol supplements improves glucose metabolism and insulin resistance [26]. Chondroitin sulfate, keratin sulfate, heparan sulfate, and hyaluronic acid are glycosaminoglycans (GAGs). GAGs play a vital role in cell physiology including cell signaling, proliferation, and cell adhesion. In T2DM, alteration in GAGs structures and functions can occur [27]. The insulin-sensitizing and anti-diabetic impacts of some GAGs have been reported [28,29]. Perturbations in GAGs related pathways and sphingolipid metabolism imply the role of the ECM in insulin resistance, which involved in the regulation of insulin action. Our analyses suggested that it seems muscle dysmetabolism disrupts the abundance of metabolites involved in the process of sensing insulin and transmitting insulin signal into the cell. Decreased insulin sensitivity results in lower expression of insulin-responsive genes, reduced glucose uptake, and consequently changes in energy, glucose, lipids, and amino acids metabolisms (Fig. 4).
As the nal analysis, potential metabolite markers were predicted. For this purpose, we rst identi ed exchange metabolites that their ux interval was shifted in comparison with healthy ones. Then, a wrapper feature selection method applied to nd potential metabolite markers. This approach led to the identi cation of 13 exchange reactions that could discriminate healthy individuals from T2DM patients with 81% accuracy. We validated these markers using a separate gene expression data from another study that has investigated the gene expression pattern in the muscle of obese and non-obese of healthy and diabetic individuals. Their result has shown that transcriptional reprogramming in obesity is similar to that occur in T2DM [7]. Here, we found that using all data from this study, including normoglycemic obese individuals for validation, the accuracy of our proposed marker was notably low (~ 50%). We thought that this low accuracy may be due to the metabolic similarities between obese and diabetic individuals. In fact, as noted in the original data article, obese and diabetic individuals have shown similar gene expression patterns [7] that can lead to similar metabolisms. Moreover, obese individuals had high levels of fasting insulin levels, which demonstrate the insulin resistance in this group. To test this issue, we removed obese healthy individuals from the validation set and examined the classi cation result, which leads to the improvement in accuracy to 78.50%. Therefore, this analysis con rmed both the original data article claim about the similarity of the gene expression pattern between obese and diabetic individuals and the appropriate e ciency of our proposed markers in identifying insulin-resistant individuals. These markers represent some of the insulin-resistance associated abnormalities represented in exchanging metabolites. Important metabolites such as methylglyoxal, hyaluronan, retinoic acid (vitamin A), sodium, alanine, and aspartate are present in the predicted markers. Notably, this method also successfully identi ed glucose as one of the markers. Methylglyoxal, a glycolytic by-product, is a toxic and highly reactive compound involved in cellular dysregulations. This compound modi es nucleotides, proteins, and lipids producing advanced glycation end products (AGEs) which also contribute to the diabetes complications. Methylglyoxal is associated with oxidative stress, cellular in ammation, and age-related disease such as diabetes [30]. Several studies have revealed the impact of methylglyoxal on insulin signaling pathways and insulin resistance [31,32] and recently this metabolite has introduced as an emerging marker for T2DM diagnosis [30]. Hyaluronan is an anionic GAG metabolite implicated in several functions like as cell signaling, proliferation and migration, and angiogenesis. This metabolite also contributes to the in ammation and pathogenesis of T2DM [33,34]. Studies have demonstrated that hyaluronan increases in the serum and skeletal muscle of T2DM subjects [35,36]. Several analyses have shown the possible roles of vitamin A in glucose metabolism, and the progression of insulin resistance [37][38][39]. Association of purine metabolites such as xanthine and hypoxanthine with the risk of T2DM incidence and complications has been reported [40,41]. Also, change in serum concentration of amino acids [42], and sodium [43] is associated with obesity and T2DM. In addition, the implication of gangliosides in insulin resistance has been shown [44][45][46]. Here, GQ1b ganglioside was predicted in the top-ranked metabolites list that can be considered for future analysis. We also checked metabolomicsbased studies of T2DM and the Human Metabolome Database for these metabolites [47]. We found that glucose, hypoxanthine, alanine, aspartate, galactose, hyaluronate, and methylglyoxal levels have been reported to be associated with T2DM [40,47,48]. These metabolite markers can be used for further empirical investigation to verify their prognostic and diagnostic values in insulin resistance.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The dataset analyzed during the current study is available in the dbGaP database with the accession code phs001068.v1.p1. This data is accessible through the repository's data access request procedures (https://dbgap.ncbi.nlm.nih.gov/). The validation datasets are available at GEO repository with the accession numbers GSE63887 and GSE81965 (https://www.ncbi.nlm.nih.gov/geo/). The muscle metabolic model reconstructed in this article is available at https://github.com/Maryamkhn/Muscle_MetNet_CBB

Competing interests
The authors declare no con icts of interest.

Funding
Not applicable.
Author's Contributions MK: Conceptualization, design, data gathering, implementation, analysis, interpretation, writing. AMB-M: Supervision, design, interpretation of data, review & editing. KK: Supervision, interpretation of data, review & editing. AM: Supervision, review. All authors have read and approved the nal manuscript.