Identifying anti-growth factors for human cancer cell lines through genome-scale metabolic modeling

Human cancer cell lines are used as important model systems to study molecular mechanisms associated with tumor growth, hereunder how genomic and biological heterogeneity found in primary tumors affect cellular phenotypes. We reconstructed Genome scale metabolic models (GEMs) for eleven cell lines based on RNA-Seq data and validated the functionality of these models with data from metabolite profiling. We used cell line-specific GEMs to analyze the differences in the metabolism of cancer cell lines, and to explore the heterogeneous expression of the metabolic subsystems. Furthermore, we predicted 85 antimetabolites that can inhibit growth of, or even kill, any of the cell lines, while at the same time not being toxic for 83 different healthy human cell types. 60 of these antimetabolites were found to inhibit growth in all cell lines. Finally, we experimentally validated one of the predicted antimetabolites using two cell lines with different phenotypic origins, and found that it is effective in inhibiting the growth of these cell lines. Using immunohistochemistry, we also showed high or moderate expression levels of proteins targeted by the validated antimetabolite. Identified anti-growth factors for inhibition of cell growth may provide leads for the development of efficient cancer treatment strategies.

H uman cancer cell lines are widely used model systems for studying cellular mechanisms underlying cancer by analyzing their perturbation-response patterns in simplified experimental conditions 1 . Cell lines are derived from human tumors of diverse tissue origin, and have adapted to growth in vitro for prolonged periods. Comparative analysis of cell lines in combination with system-wide profiling techniques may disclose cross cell-type commonalities and variations of biological processes. This knowledge can be used for understanding cancer metabolism, identifying anticancer drugs, in vitro evaluation of predicted drug targets and outlining mechanisms of action of therapeutic agents 2, 3 .
Advances in omics technology have enabled simultaneous measurement of molecular components interacting within complex interconnected networks and hereby provided insights into cellular functions and phenotypic states of the cells. However, analyzing genome-wide data and in particular gaining new biological knowledge, is a nontrivial effort 4 . Reconstruction of genome scale metabolic models (GEMs) can assist in this by integrating omics data in multiple layers to understand the effects of local interactions in the context of the whole network [5][6][7] . This makes GEMs a potential tool for analyzing omics data in health and disease states, and for identifying the underlying cellular mechanisms in the occurrence of complex diseases. To date, several generic human GEMs have been generated [8][9][10][11] and these models have been employed for reconstruction of context-specific GEMs for healthy and cancerous cell-types. GEMs have been developed for studying the metabolic alterations between healthy and cancerous tissues 12,13 as well as within cancerous tissues 14 . Moreover, personalized cancer GEMs has been reconstructed and used to capture common and specific metabolic shifts across the cancer patients and to identify selective anticancer drugs 15 .
Preserved proliferative signaling and avoidance of growth suppressors have been identified as hallmarks of cancer 16 . Adapting the cellular metabolism (e.g. increased synthesis of nucleotides, lipids and proteins) to accommodate tumor expansion is a critical feature of cancer, and has been leveraged for development of novel anticancer drugs 17,18 . In the present study, we used the concept of antimetabolites, which are structural analogues of endogenous metabolites. Antimetabolites subvert cellular processes by serving as inhibitors of all enzymes involved in metabolizing the associated endogeneous metabolite, and hereby dramatically affecting metabolic functions 19,20 . An important characteristic of antimetabolites is their potential to simultaneously inhibit multiple enzymes, and thereby reduce the growth of proliferating cells more effectively 21 . Antimetabolites are among the most common anticancer drugs since the discovery of aminopterin, an effective drug in remission of leukemia. Examples of antimetabolites are antifolates (e.g. Methotrexate), antipyrimidines (e.g. Cytarabine, 5-Fluorouracil) and antipurines (e.g. . In this study, we first followed a systematic approach and analyzed the global mRNA expression pattern of the 20,314 protein coding genes in eleven human cancer cell lines 22 . Secondly, we reconstructed functional cell line-specific GEMs (CL-GEMs) for these eleven cell lines using mRNA expression levels (RNA-Seq) together with our tINIT (task-driven Integrative Network Inference for Tissue) algorithm 15 and Human Metabolic Reaction database (HMR)2 11 ( Figure 1A). We included known metabolic functions of the cell lines during the reconstruction process, and validated functionality of the models based on metabolite consumption and release (CORE) profiles of the cancer cell lines 1 . Thirdly, we compared CL-GEMs by using a statistical multicomparison method as well as here defined heterogeneity scale. We also compared metabolic subsystems in the CL-GEMs based on the distribution of mRNA expression patterns of associated genes across the cell lines. Fourthly, we screened all metabolites in the HMR2 to identify essential metabolites for the growth of cell lines and predicted potential antimetabolites that can inhibit or kill their growth. We defined essential metabolites as metabolites for which their analogues (antimetabolites) may disrupt the metabolic network and lead to cell death. Each antimetabolite was also evaluated for its in silico toxicity by employing GEMs for 83 human healthy cell-types. Next, we experimentally tested the effect of an L-carnitine analogue, one of the potential antimetabolite that was predicted to be effective on all cell lines. Significant difference in proliferation of cultured cell lines in presence or absence of selected analogue, confirmed our computational prediction. Finally, we presented high or moderate protein staining level of the CPT1A and CPT1B, which were targeted by the L-carnitine analogue in cell lines using the antibodies generated in the Human Protein Atlas (HPA).

Results
Transcriptomics data of Human Cancer Cell lines. RNA-Seq have been employed to identify the tissue-specific expression of genes across the major human tissues and combined with the antibody based profiling of the same tissues [23][24][25][26][27] . Previously, we also performed a comparative transcriptomics analysis based on the quantitative RNA profiling of protein-coding genes across eleven human cancer cell lines derived from distinct cellular phenotypes 22 . Fragment per kilobase of exon model per million mapped read (FPKM) values were calculated for each gene, and expression levels were categorized into high (FPKM . 50), medium (FPKM . 20) and low (FPKM . 1), and FPKM , 1 was used to define genes as not detected. Based on this cut-off the number of protein coding genes expressed ranged from 11,000 to 12,500 across the eleven cell lines. Approximately 43% of the expressed genes for each cell line were found in the category of high and medium RNA expression levels and the average transcript abundance score in a log2 scale ranged from 4.90 for the epidermoid carcinoma cell line (A-431) to 5.30 for the metastatic breast adenocarcinoma cell line (MCF-7) ( Figure 1B). The transcript profiling of cell lines was used to study the cross-cell line distribution of detected protein coding genes ( Figure 1C). In total, 15,292 out of 20,314 protein coding genes were expressed in at least one of the cell lines, from which nearly 58% (n 5 8836) were commonly detectable transcripts shared by all cell lines and the remaining 42% (n 5 6456) distributed in different combination of cell lines, with the highest contribution of cell specific genes (n 5 1285). Through in-depth analysis of the data, we found that there is a connection between the mRNA expression levels of detected genes and their presence across cell lines. We found that genes detected in all eleven cell lines had higher transcript expression level compared to cell specific genes, with an average of 9 fold change (Supplementary Dataset 1).
Reconstruction of GEMs for human cancer cell lines. A generic genome-scale metabolic model of cancer has been developed to capture the common metabolic functions between cancer types 28 . A core of 197 highly expressed metabolic genes common to more than 90% of the cell lines in the NCI-60 cell lines database have been employed to reconstruct a generic cancer metabolic network, and the resulting model was used to predict anti-cancer drug targets 28 . Here, we reconstructed CL-GEMs for eleven human cancer cell lines based on cell line specific RNA-Seq data, HMR2 and our task-driven reconstruction approach. Lists of 56 metabolic tasks that are known to occur in all cell types were used during the reconstruction process. These metabolic tasks can be categorized as processes for provision of energy and redox, biosynthesis of metabolites, substrate utilization and internal conversion processes (Supplementary Dataset 2). Growth was also added as a metabolic task to ensure that the cells can proliferate. The CL-GEMs were reconstructed such that the models should be able to perform all the defined tasks and at the same time maintain the consistency with the transcriptomics data. The eleven CL-GEMs are presented in the Human Metabolic Atlas portal (http://www.metabolicatlas.org/) in the Systems Biology Markup Language format.
The resulting CL-GEMs contained 5,297 to 5,584 reactions associated with 2,193 to 2,328 genes and involving 4,209 to 4,432 metabolites (Table 1 . We also analyzed the transcript expression pattern of CL-GEMs and observed an average expression level of 34 to 58 FPKM for all genes incorporated into the metabolic models, approximately 4 fold higher than the expression level of cell linespecific genes (Supplementary Dataset 3), which shows that metabolism represents a core function in the cell lines.
We used CORE metabolite profiles of cancer cell lines 1 to validate the functionality of four CL-GEMs (PC-3, A-549, MCF-7 and U-251) that were common in both studies. Constraining and optimizing the CL-GEMs based on CORE profiles, we observed that all models obtained a feasible solution with flux distribution consistent with published experimental data (Supplementary Dataset 5).
Heterogeneity of CL-GEMs. Cancer cells reprogram their metabolic network to support cell proliferation by synthesizing all the cellular components and the required Gibbs free energy 16 . Despite the shared ability of transformed cells to induce growth stimulating adaptations towards uniformity of metabolism, significant marks of diversity and complexity in metabolic signatures of neoplasms has been recognized in recent years 29,30 . In order to understand the heterogeneity of metabolic networks in cancer, we used the reconstructed CL-GEMs and analyzed the transformability of models based on differences in their constituent parameters: genes, metabolites and reaction. We  Figure 3A). Higher heterogeneity degree of genes compared to reactions indicates the preference of parameter effect to dimension effect, considering the fact that the number of  . A fall in heterogeneity degree, approximately 0.07 degree based on genes and reactions, comparing healthy cell-types to cell lines revealed a tendency towards the uniformity of metabolic networks in cell lines. However, a significant transformation was not observed and the models maintained their heterogeneity. Interestingly, the heterogeneity degree of CL-GEMs based on metabolites was stable, in contrast to genes and reactions, and this indicated the preserved capability of altered metabolic networks to transform the metabolites similar to healthy cell-types. We also checked the number of the flux carrying reactions in each model (Table 1) and found that the models show higher heterogeneity based on the flux carrying reactions even though the number of the reactions in the models decreased by ,2 fold ( Figure 3B and Supplementary Dataset 6). This analysis clearly demonstrated the difference in the metabolism of the cancer cell lines.
Next, using the transcriptome expression profile of genes included into the CL-GEMs, we investigated the global shifts between the cell lines. Kruskal-Wallis one-way analysis of variance 31 was employed to compare the distribution of mRNA expression pattern between models and calculated P-value 5 0.0029 rejected the hypothesis of similarity of distributions at significance level of a 5 0.01. Furthermore, we performed pair wise comparisons of CL-GEMs using Fisher's LSD test 32 with a confidence level of 0.05 (Supplementary Dataset 7). We found that Hep-G2 is the cell line with highest rate of significant difference compared to the other cell lines, and all models except Lung carcinoma cell line (A-549) and CACO-2 cell lines had a significant difference compared to the Hep-G2 model. The A-549 and CACO-2 models were also significantly different from U-251 MG and U-2 OS models.
Thus, comparing GEMs through constituent parameters (genes, metabolites and reactions) and also using mRNA expression pattern disclosed a considerable difference between metabolic models for the eleven cell lines. This indicated the necessity of developing fine-tuned CL-GEMs to capture the metabolic alterations that are specific to each cancer cell line.
Heterogeneous expression of metabolic subsystems. We analyzed the differences in mRNA expression pattern of the individual pathways across the eleven cell lines. We calculated the fraction of cell lines in which each metabolic subsystem had an average high transcript expression (H) and low transcript expression (L), and transferred the results into the (H 2 L) and (H 1 L) coordinates to observe the behavior of subsystems across cell lines ( Figure 4A). We found highly, lowly and heterogeneously expressed subsystems in all eleven cell lines, and this allowed us to focus on subsystems that can be used for identifying generic or cell line-specific targets. As expected, biosynthesis pathways that are essential for biomass production and cell division displayed high expression levels in most of the cell lines, from which fatty acid biosynthesis, b-oxidation of fatty acids and cholesterol biosynthesis had the highest fraction of the H. Along with these pathways, ROS detoxification, folate metabolism and aminoacyl-tRNA biosynthesis displayed frequent high expression pattern in many of cell lines.
In contrast, lipoic acid metabolism and glucocorticoid biosynthesis revealed a relatively high fraction of L across all cell lines, which followed by retinol metabolism and androgen metabolism. Lipoic acid has been shown to reduce the proliferation/viability and to increase apoptosis of tumor cells 33 . Moreover, it has been reported that deficiency of lipoic acid synthetase results in defected mitochondrial energy metabolism and glycine cleavage mechanism as reflected by lactic acidosis and elevation of glycine 34,35 . Glucocorticoids have been used to treat wide range of cancers based on their substantial effect on progression of cell cycle and apoptosis through glucocorticoid receptor mediated mechanisms [36][37][38] . Retinoids are used for cancer treatment because of their ability to promote cell differentiation and also to change the gene expression pattern of tumor cells making them more vulnerable to other therapies [39][40][41] . We next studied the inter-cell expressional variation of subsystems to explore the heterogeneous behavior of metabolic pathways across the different cell lines. Bile acid recycling, fatty acid desaturation, ascorbate and aldarate metabolism, formation and hydrolysis of cholesterol esters, and vitamin D metabolism showed higher heterogeneity among all subsystems in HMR2 ( Figure 4B). In CACO-2, cervical epithelial adenocarcinoma cell line (HeLa) and Hep-G2, genes involved in bile acid recycling pathway mainly had low mRNA expression level, whereas in A-431 and RT-4, these genes had high expression levels. Almost all genes associated with reactions in fatty acid desaturation pathways were highly expressed in CACO-2 and Hep-G2 without showing a considerable high or low express- ion pattern for other cell lines. Formation and hydrolysis of cholesterol esters pathway was completely dominated by highly expressed genes for A-431 and A-549 different from other cell lines. Vitamin D metabolism revealed a different genes expression pattern for U-251 MG with low transcript level for all associated genes which followed by RT-4 and embryonal kidney cell line (HEK 293) with more than 70% of genes with low expression level. These analyses suggested that the activities of heterogeneous pathways are not influenced only by variability of environmental conditions, but also by the specific genetic repertoire of each individual cell line.
Anti-growth Factors for Cell Lines. We used reconstructed CL-GEMs to predict the potential antigrowth factors through analyzing the robustness of metabolic networks against the induced perturbations. The metabolite-centric approach based on the concept of antimetabolites was employed to introduce the growth inhibiting perturbations. By testing the essentiality of metabolites present in the HMR2 on growth of cancer cell lines, we predicted 138 antimetabolites that can kill the growth on any of eleven studied cell lines and 106 of them were found to be effective in all cell lines. However, using these antimetabolites for the treatment of the cancer may have toxic effects to healthy cell-types in the human body. Therefore, we performed in silico toxicity test using the previously published GEMs for 83 different cell types in the human body. We checked if the use of the antimetabolites would potentially disrupt the metabolic tasks in energy production.
We finally predicted 85 potential antimetabolites that can inhibit growth of any cell lines studied here, while not being toxic to healthy cells (Supplementary Dataset 8). We found that 70.6% (n 5 60) were effective in all cell lines, 11.8% (n 5 10) were effective in two cell lines, 17.6% (n 5 15) were effective in just one cell line, and interestingly not any antimetabolite was detected for 3 to 10 combinations of cell lines (Supplementary Dataset 9). Although ,70% of antimetabolites were predicted to be effective in all cell lines, the fact that 30% were cell-line specific indicated the importance of using cell line specific models rather than depending on one generic cancer model. Furthermore, we categorized 85 predicted antimetabolites based on their relevant subsystems in HMR2 and calculated the average mRNA expression level of the subsystems across cell lines ( Figure 5A, Supplementary Dataset 10). It was observed that nucleotide metabolism, amino acids metabolism and cholesterol biosynthesis contain the highest fraction of antimetabolites, more than 68%. This was followed by sphingolipid, glycerophospholipid and folate metabolism. In general, pathways containing higher fraction of predicted antimetabolites showed high average mRNA expression levels in all cell lines and the exceptions were closely connected to one or more pathways with high average expression levels. For example, the carnitine shuttle had relatively low average expression level, but deficiency in this pathway results in deficiency of the b-oxidation pathway that has one of highest average RNA expression levels within all cell lines.
Following the prediction of L-carnitine as an antimetabolite and considering the connection of the L-carnitine shuttle to the mitochondrial b-oxidation, the analogue of L-carnitine were proposed as a potential antigrowth factor targeting the growth in all eleven cell lines ( Figure 5B). L-carnitine is an essential metabolite that plays an important role in fatty acids metabolism and energy production. In humans, carnitine homoeostasis is preserved by dietary absorption and also through a multi-step de novo biosynthesis 42 . L-carnitine can be synthesized in liver and kidney using essential amino acids lysine and methionine as primary precursors to form trimethyllysine (TML) 43,44 . Notably, both lysine and methionine were also detected as antimetabolites through our analysis.
b-Oxidation was found to have high average mRNA expression levels in all eleven cell lines. This pathway consists of recurrent series of reactions that result in cyclical shortening of fatty acids and production of acetyl CoA, NADH and FADH. Generated NADH and FADH are used for ATP production through electron transport chain, while acetyl CoA enters the citric acid cycle (TCA) [44][45][46] . Carnitine shuttle acts as the carrier for the long chain fatty acid acyl groups and transports them from the cytosol to the mitochondria, where b-oxidation occurs. L-carnitine has an important role in transfer of the acetyl-CoA produced by peroxisomal b-oxidation to the mitochondria for oxidation. Other functions performed by L-carnitine include storing energy as acetylcarnitine, regulating the ratio of acyl-CoA to CoA and exerting the poorly metabolized acyl groups as carnitine esters to adjust their toxic effects [42][43][44]47,48 .
The use of L-carnitine analog inhibits the carnitine palmitoyltransferase (CPT) 1 and 2 enzymes that catalyze the L-carnitine and fatty acyl CoAs. We used perhexiline malate salt (perhexiline) that inhibits CPT1 and partly CPT2 to mimic the effect of the potential L-carnitine antimetabolite. To confirm our prediction about the use of L-carnitine analogue as potential antigrowth factor for all cell lines, we tested the effect of Perhexiline on the proliferation of prostate carcinoma cell line (PC-3) and A-431. Perhexiline was used to imitate the L-carnitine mechanism of effect since mitochondrial enzyme CPT1, as a member of carnitine acyltransferases family of enzymes, is responsible for translocation of conjugated L-carnitine and long chain FAs from cytosol to the mitochondria ( Figure 5B). We treated two cell lines with four different concentrations of Perhexiline (2, 4, 8 and 24 mM) and investigated the viability of cells after 24 and 48 hours with eight replicates (Figure 6A-B). For lower concentration of Perhexiline (4 and 6 mM), we repeated our experiments with 16 replicates for 48 hours (Figure 6C-D). Significant reduction (t-test, p-value 0.05) in cell lines viability was observed with Perhexiline concentrations more than 2 mM for PC-3 and 4 mM for A-431. We also analyzed protein expression level of CPT1 and CPT2 using the antibodies generated in the HPA project 49 and observed that CPT1 and CPT2 had strong or moderate protein expression levels. Protein expression levels as detected by using immunoshistochemistry are exemplified in a subset of cell lines in Figure 7.

Discussion
Recent advances in understanding cancer cell biology and continuous development of high throughput analytical methods have revealed that cancer is remarkably complex and heterogeneous. Increased understanding of this genomic heterogeneity between different tumors may open a new window towards the development and optimization of therapeutic methods for specific cancer types. GEMs provide a mechanistic description of relationships between genes and metabolic functions, and they can therefore be used for the interpretation of large experimental datasets. In this study we reconstructed CL-GEMs for eleven human cancer cell lines based on RNA expression data and used the generated models to understand the heterogeneity of metabolic networks in cancer cell lines. Our analysis revealed notable heterogeneity across CL-GEMs, where each Comparing CL-GEMs with GEMs for 83 healthy cell-types, we found a slight tendency towards the uniformity of metabolic networks in CL-GEMs and the heterogeneity of models preserved to a large extent. This finding indicates inadequacy of reconstructing a generic cancer model to capture metabolic alterations in transformed cells.
We analyzed the RNA expression pattern of individual pathways to identify the transcriptome based distribution of subsystems across CL-GEMs. We found that genes involved in fatty acid biosynthesis and cholesterol biosynthesis had highest average mRNA expression level in most of the cell lines, which shows the importance of biosynthesis pathways in transformed cells. On the other hand, lipoic acid metabolism and glucocorticoid biosynthesis revealed lowest average expression level among subsystems consistent with previously published studies indicating their positive correlation to reduced cell proliferation/viability and increased apoptosis of tumor cells 33 . Through in-depth analysis of expressional variation of subsystems across CL-GEMs, we found that bile acid recycling, fatty acid desaturation, ascorbate and aldarate metabolism, formation and hydrolysis of cholesterol esters, and vitamin D metabolism are the pathways with highest expressional heterogeneity among all subsystems in HMR2. Subsystems with high expressional variation may be potential targets for further studies to find cell line specific targets.
Metabolite essentiality, defined as metabolites necessary for cell growth and survivor, has previously been used 50 to identify potential drug targets for an opportunistic pathogen using GEMs. We used this concept to explore essential metabolites in CL-GEMs by perturbing the metabolic networks and found metabolites contributing to growth inhibiting perturbations. The structural analogues of 85 predicted essential metabolites, which passed the in silico toxicity test based on previously generated GEMs for 83 healthy human cell types, proposed as potential antimetabolites. Analogues of L-carnitine were proposed as potential antimetabolites, due to their key role in L-carnitine shuttle pathway, and the contribution of L-carnitine shuttle to the mitochondrial b-oxidation and consequently to the fatty acid synthesis. Two cell lines with different phenotypic origins, PC-3 and A-431 were selected for in vitro testing of our prediction Plasticity is one of the key factors to maintain fitness of transformed cells and b-oxidation of fatty acids may prepare some of this plasticity by supplying metabolic intermediates for cell proliferation, suppressing pro-apoptotic subsystems, eradicating potentially cytotoxic lipids, and empowering the production of NADPH and ATP when required [51][52][53][54][55] . However, b-oxidation of fatty acids cannot be discerned as an independent metabolic pathway in cancer metabolism and there is a challenge to unify the essential role of b-oxidation pathway in tumor cells with the fact that active fatty acid synthesis is necessary for cancer cells division and proliferation.
Despite promising performance of GEMs in identification of biomarkers and drug targets for cancer, there are cautions and limitations regarding this approach. Unlike bacteria, reconstructing and testing the functionality of GEMs is more complex when it comes to cells and tissues. Constructing new generation of unified GEMs that integrate metabolism with other cellular processes, especially signaling and regulation, remains as a nontrivial challenge. Different assumptions to estimate the connection between transcription, translation and flux rate challenges the accuracy and compatibility of models generated by incorporating omics data into GEMs.
In conclusion, we used CL-GEMs to investigate the differences in the metabolism of cell lines and explore potential antigrowth factors, which may impede the growth of cancer cells, and experimentally verified the inhibitory effect of the L-carnitine analogue, one of the predicted antimetabolite. The study of cancer metabolism using GEMs has revealed new therapeutic opportunities, as well as more profound understanding of metabolic reprogramming of cancer cells. Over the past decades we learned that targeting single pathways or enzymes is not enough to cure cancer and agents proposed here as antigrowth factors may be combined with other therapeutic agents and methods to get desired results. The approach presented here may also be extended to study the combinatory effect of potential and standard therapeutic agents to find more effective cancer specific therapies, and may present new possibilities towards personalized medicine.
Methods GEM reconstruction. The tINIT algorithm 15 implemented in the RAVEN toolbox 56 was used for the reconstruction of eleven CL-GEMs. 56 Metabolic functions known to occur in all cell lines as well as growth were used during the reconstruction process of the GEMs (Supplementary Dataset 2). During the reconstruction of the model, we used the content of Ham's media as minimal media required for cell growth. The content of the Ham's media is provided in Supplementary dataset 2.
Heterogeneity of CL-GEMs. A new measure of heterogeneity was used in our study to capture the divergence between metabolic networks based on their constituent parameters: genes, reactions and metabolites. For each model, heterogeneity degree formulated as: Where d h ij is the heterogeneity degree of model i based on parameter j, h i j is the average Hamming distance of model i to all other models based on parameter j, H i j is maximum Hamming distance of model i based on parameter j in comparison to integrated vector of parameter j (Iv j ). In this formulation, n is the number of models (83 for healthy cell-types and 11 for cancer cell lines), m is the number of the constituent parameters (here parameter are genes, reaction and metabolites, m 5 3), and Iv j is a unique union of corresponding parameter across all models. For example, I.v 1 represent a union vector of genes present at least in one of the models.
To compare the CL-GEMs based on RNA expression levels, we built a vector of expected transcriptome expression values of genes linked to metabolic subsystems for each model, and merged these vectors into n 3 m unified ratio matrix. Here, n is number of subsystems and m is the number of the models. We used this matrix to perform multi-comparison analysis employing Kruskal-Wallis non-parametric test with significance level of 0.01. Finally, we performed pair wise analysis of models using post hoc Fisher's LSD test with confidence level of 0.05.
Heterogeneous expression of metabolic pathways. To identify the heterogeneous expression of metabolic subsystems across the cell lines, we calculated the ratios of high and low RNA expression level for genes associated to reactions inside the CL-GEMs in logarithmic scale. For each pathway, we calculated the average ratio of high expression level and low expression level and integrated the results into two different n 3 m ratio matrices, high matrix (HM) and low matrix (LM), respectively. Here, n is number of pathways and m is the number of models. We used a subsystem-specific metrics, HM-LM and HM 1 LM, to compare the expression pattern of pathways across the cell lines.
In silico toxicity test. Inhibition of essential metabolites in cancer cell lines can also disrupt the metabolism of healthy cells, similar to chemotherapeutic drugs. To avoid this side effect, we used reconstructed GEMs for 83 healthy cell types and investigated the response of these metabolic models against inhibition of essential metabolites. We excluded the essential metabolites that disrupted energy metabolism and redox balance in healthy cells. Our approach was based on the central role of energy producing tasks for all cell types, and also high likelihood of neoplastic shifts in energy metabolism of transformed cells.