Pinpointing Cancer Sub-Type Specific Metabolic Tasks Facilitates Identification of Anti-cancer Targets

Introduction

Metabolic reprogramming is one of the hallmarks of tumorigenesis (1) Understanding the metabolic changes in cancer cells provides attractive therapeutic targets, and new strategies for cancer therapy (2). Altered metabolic activity is vital for sustaining uncontrolled proliferation, evasion of growth suppression, resistance of cell death and the metastasis to other areas. Cancer cells have many metabolic characteristics. In order to provide enough energy source for rapid proliferation, cancer cells involve a shift from mitochondrial metabolism toward glycolysis, even in the presence of oxygen, which is called Warburg effect (3).

The metabolic states are not the same in different cancer types or subtypes, even within the same sample of solid tumors. Metabolic reprogramming is largely based on genetic changes and environment (4, 5), and those different reprogramming routes contribute to metabolic heterogeneity. Many researches, including transcriptomics (6), proteomics (7) and metabolomics (8), show metabolic heterogeneity across cancer types. Although in most cancers, nucleotide synthesis and glycolysis are upregulated, others, like oxidative phosphorylation, vary a lot (6). Cancer subtypes also display diverse metabolic phenotypes. For instance, triple negative breast cancers (TNBCs) rely more on glucose and glutamine uptake than ER+ breast cancers (9). Understanding metabolic heterogeneity is of significance for identifying metabolic vulnerabilities susceptible to therapeutic targeting (2).

Taken into consideration of all those metabolic objectives, it is impossible to simply synthesize one key general criterion for all cancer types, even subtypes that originate from the same tissue. In comparison of experimental and model-predicted flux, Schuetz et al. demonstrated that the combination of two or three metabolic objectives and modeling by Flux Balance Analysis (FBA) offers a better performance for the explanation of flux distribution in microbe modeling (10, 11). The theory of evolutionary trade-offs can shed light on the biological significance of multi-task optimization, as cells need to perform multiple tasks that all contribute to fitness in gene expression space (12, 13). Using a multi-objective genome-scale metabolic model is a reasonable approach to improve the accuracy and prediction ability of genome-scale metabolic network modeling in evolutionary terms (14).

In order to understand the heterogeneity of cancer cells, we construct a multi-objective metabolic model (MOMM) using transcriptomics data of different cancer types. Specifically, we studied breast cancer, lung cancer, digestive organ cancer, digestive tract cancer, and reproductive cancer. To evaluate our methodology, we applied MOMM in the prediction of cell viability. We implemented MOMM to predict synthetic lethal target pairs, and many of the top-ranking pairs have been reported as lethal target pairs in previous experimental studies. Our analysis showed that the number of activated synthetic lethal target pairs is a good marker of patient survival time. This study opens a new window for the development of anti-cancer drugs and provides promising treatment plans for cancer.

Materials and Methods

Pre-processing of Data

The transcriptomics data used in this study comes from TCGA (15). Cancer types include breast cancer, lung cancer, digestive organ cancer, digestive tract cancer, and reproductive organ cancer.

(1) Breast cancer contains 1,215 samples, with four PAM50 subtypes and normal types, as well as unknown types. Among them, the basal subtype of breast cancer contains 142 samples, the luminal A subtype contains 423 samples, the luminal B subtype contains 194 samples, the HER2 subtype contains 67 samples, and the normal cell contains 137 samples.

(2) Lung cancer contains 1,129 samples, with two pathological types and normal types. Among them, lung adenocarcinoma contains 517 samples, lung squamous cell carcinoma contains 502 samples, and normal type contains 110 samples.

(3) Digestive organ cancer contains 651 samples, and there are three cancer subtypes and normal types. Among them, liver cancer contains 373 samples, bile duct cancer contains 36 samples, pancreatic cancer contains 179 samples, and normal cell contain 63 samples.

(4) digestive tract cancer contains 444 samples, and there are two cancer subtypes and normal types. Among them, colon cancer contains 288 samples, rectal cancer contains 95 samples, and normal type has 51 samples.

(5) Reproductive organ cancer contains 874 samples, and there are four cancer subtypes and normal types. Among them, ovarian cancer contained 308 samples, endometrial cancer contained 177 samples, cervical cancer contained 305 samples, uterine carcinosarcoma contained 57 samples, and normal cells contained 27 samples.

The genome-wide metabolic network used in the model is the Recon 1 model supplemented by Duarte et al. (16). On top of this, we imposed additional constraints on the intake of 20 amino acids and important metabolites (Supplementary Table 4). The network used in the simulation contains 4,924 metabolic reactions (reversible reactions are disassembled and counted as 2 reactions, so that the flow rate is non-negative, which is convenient for subsequent calculations) and 2,767 metabolites. By comparing the existence of metabolic genes in the metabolic network in the transcriptome data, we obtained a total of 1,442 metabolic-related gene expression profile data for subsequent analysis.

Inference of Metabolic Tasks Using Pareto-Tasks Inference Method

We used Pareto-tasks inference method to interpret the number of metabolic tasks (17). By using the default algorithm, we used three criteria to evaluate the propriate number of tasks. (1) Simplicity: The greater the number of metabolic tasks, the greater the explained variance and the more information it contains, but more complex the model will be. We used the Elbow plot to find the best trade-off task number between information content and model complexity. (2) Interpretability, we used the t-ratio test, the smaller the P-value. The more explanatory the model is for biological significance. (3) Stability: The average error value of the vertex position of the multi-objective model is calculated by the bootstrapping method. The smaller the average error value, the stronger the model stability.

Evaluation of Multi-Objective Metabolic Model

The relationship between gene expression and survival probability is represented by hazard ration (HR). Genes with HR > 1 means that their low expression corresponds to high survival probability; genes with HR < 1 means that their high expression corresponds to high survival probability. We select genes whose objective function value are within [Obj_min + 0.1*(Obj_wt−Obj_min), Obj_wt] after their knockout to calculate their receiver operating characteristic (ROC) curve, where Obj_wt is the wild-type objective function value before gene knockout simulation, Obj_min is the minimum value of the objective function in all gene simulations, that is, the most affected state.

Integration of Metabolic Tasks and Genome-Scale Metabolic Network Model

The metabolic tasks that predicted by Pareto task inference (ParTI) method are represented by the metabolic gene expression level. Integrating the metabolic tasks into the genome-scales metabolic network is essentially integrating the gene expression to the flux model.

First, we mapped the gene expression level to metabolic enzyme expression level. In the genome-scale metabolic network, the relationship between reactions and genes and proteins is connected by gene-protein-reaction (GPR) relationships, which are Boolean expressions between transcripts, proteins, and the corresponding reactions. Isozymes are represented by OR relationship, and we used the sum of the expression of isozymes as the total expression of the enzyme. If the protein contains multiple subunits, it is represented by AND relationship. We used the minimum value of the subunits to represent the total expression of the enzyme.

${\mathrm{r}}_1\hspace{1em}{\mathrm{g}}_1\text{OR}{\mathrm{g}}_2\mathit{\hspace{1em}\hspace{0.5em}\hspace{1em}}\mathrm{sum}({\mathrm{g}}_1,{\mathrm{g}}_2)$

${\mathrm{r}}_2\hspace{1em}{\mathrm{g}}_3\text{AND}{\mathrm{g}}_4\mathit{\hspace{1em}}\min ({\mathrm{g}}_3,{\mathrm{g}}_4)$

Then, we integrated the enzyme expression to flux by E-Fmin method (18),

$\text{max}\sum _i{w}_i{v}_i(1)$

Subject to,

$\mathrm{S}=0(2)$

$v_L\le \mathrm{\nu }\le {\mathrm{\nu }}_U(3)$

$\mathrm{\epsilon }\le v_{biomass}(4)$

Where w_i is the function of enzyme expression, $w_i=\overline{r_i}-1$, $\overline{r_i}$ is the normalized enzyme expression. v_i is the flux of the i^th reaction. S is the stoichiometric matrix (m× n), m represents the number of metabolites, n represents the number of reactions. ν is the sum of n flux, v_L and ν_U are the lower and upper bound of ν. v_biomass is the flux of biomass, its lower bound ε depends on the real situation. Through the prediction of the model, we can calculate the optimal flux distribution of the corresponding metabolic task.

In silico Simulation of Gene Knockouts

We used minimization of metabolic adjustment (MOMA) to calculate the flux distribution after knockout (19). Essentially, we use this method to find new flux distribution which is the closest to the original flux distribution in new feasible area.

$\text{D}(\text{w},\text{x})=\sqrt{\sum _{i=\mathrm{\hspace{0.33em}1}}^N{\left(w_i-v_i\right)}^2}(5)$

Where w is the original flux distribution, w_i is the ith flux, v ∈ ϕ_j is the flux distribution after knockouts, v_i is the ith flux after gene modification, N is the number of fluxes. D(w,v) represents the Euclidean distance between w and v. We want to find v which let D minimizes. So we turn (5) into (6), i.e., minimize,

$\text{f}\left(\text{x}\right)=L\cdot v+\frac{1}{2}{v}^{T}Qv(6)$

Where L and Q are the linear and quadratic parts of the objective function, the length of L is N, the size of Q is N × N, v^T is the transpose of v. Minimizing (6) is equivalent to minimizing the quadratic formula inside the rood, and the constant term after expansion can also be ignored. Therefore, minimizing D is equivalent to minimizing f(v). If only the optimal solution exists, and the feasible region is not empty after adjustment, we can always find a solution to this optimization problem. The uniqueness of the solution is guaranteed by the convexity of f(v).

Identification of Synthetic Lethality Pairs

We simulated the paired knockouts of all genes in the metabolic network, and used MOMA to calculate the flux distribution afterward. According to “highest single agent” synergy model (20), i.e., synthetic lethality pair represents the combined effect outperforms either its single component. Here, we exclude “transport” and “exchange” related reactions, and only consider intracellular reactions.

We introduce Impact Score (IS) to represent the impact of the cell after knockouts. IS=(obj_max−obj_val)/(obj_max−obj_min), where obj_max represents the maximum of the objective function, which is the original value without gene modification; obj_min represents the minimal objective function value in all the combined knockout simulations; obj_val represents the corresponding objective function value after the specific gene pair knockouts. We choose gene pairs with IS≥ 0.5 in the following analysis.

For survival analysis in synthetic lethality pairs evaluation, we choose the breast cancer GEO data set, with 1,809 patients (21). And we use the subtype classification method derived by Mihály, etc., to category the sample (22).

Results

The Dominant Metabolic Tasks for Cancers

We used multi-objective metabolic model (MOMM) to answer the questions of how many and what are the dominant metabolic tasks for a specific type of cancer. We assume that every single cell should be in a Pareto-optimal situation in order to survive. Here, we used TCGA gene expression data by RNA-seq (data was obtained from UCSC Xena)¹ to build the model. Only metabolism-related expression profile data were used. The metabolism-related genes were defined as the genes that are involved in the metabolic reactions in Recon 1, a commonly used genome-scale metabolic network (16), and at the same time, have expression records in the TCGA data set (Details of data were shown in Supplementary Table 1).

According to the theory of evolutionary trade-offs between tasks, whole-cell gene expression can fall into a low-dimensional polytope whose vertices are the optimal points in each task alone (12, 13). We therefore inquired whether a low-dimensional polytope (e.g., line, triangle, or tetrahedron) could enclose the metabolism-related gene expression data. We used the ParTI method to fit the metabolism-related transcriptomics data into a low-dimensional polyhedron, in which each vertex, i.e., archetype, represents an optimal phenotype for a single metabolic tasks (17). This approach was used in many biological contexts, including population and single-cell level gene expression (13, 23).

The more the metabolic tasks are, the more accurate the model, but less simplicity. To identify the dominant metabolic tasks that can cover as much information as possible without being overly complicated, we set three criteria: simplicity, interpretability, and stability (see section “Materials and Methods”). Simplicity represents the complexity level of metabolic tasks model and is analyzed by explained variability of different number of archetypes. Interpretability represents the biological significance of the calculated metabolic tasks, i.e., the description of the data of these metabolic tasks. Stability requires that the metabolic tasks are not changed along with the sample size changes, and is calculated by bootstrapping.

By applying these criteria, we found that the combination of three dominant metabolic tasks is the best representation for breast cancer (Figures 1A–C). Although the combination of four metabolic tasks also works comparably, it adds more model complexity and is less stable.

FIGURE 1

Figure 1. Determination of the number of metabolic tasks. (A) Simplicity, the elbow plot of the number of metabolic tasks and the explained variance in breast cancer. The point with the largest distance is the best trade-off point between the amount of information and model complexity. This figure show three-task model is the best. (B) Interpretability, P-value of t-ratio test for different metabolic target numbers of breast cancer. (C) Stability, the average variance of the vertices position of breast cancer. (D–F) Lung cancer.

Applying the same criteria, we found that the combination of three metabolic tasks is the best representation for lung cancer (Figures 1D–F). From the perspective of simplicity and stability, models with three tasks are the optimal modeling direction. However, from an interpretability point of view, models with three or four metabolic tasks are similar. Therefore, we choose three tasks for the follow-up analysis of lung cancer.

We also analyzed the metabolic tasks in digestive organ cancers (include liver cancer, bile duct cancer, and pancreatic cancer), digestive tract cancers (include colon cancer and rectal cancer) and reproductive system cancers (include ovarian cancer, endometrial cancer, cervical cancer, and uterine carcinosarcoma). The models for digestive organ cancers, digestive tract cancers and reproductive system cancers contain four, three, and three dominating metabolic tasks, respectively (Details were given in Supplementary Material).

The Biological Significance of Metabolic Tasks

Based on the estimation of the number of metabolic tasks, we further analyzed the biological significance of metabolic tasks, including their corresponding subtype information and biological functions. We mainly focus on the breast cancer and lung cancer. In terms of the propensity of different subtypes of metabolic tasks, we calculate the degree of enrichment of the subtypes around the metabolic task. We inferred the significance of the metabolic task through the gene sets that enriched around the metabolic task. Through the analysis of the cancer subtype enrichment, we found that different cancer subtypes have different metabolic task tendencies, and their biological significance is also quite different.

By analyzing the PAM50 subtypes of breast cancer, including luminal A subtype, luminal B subtype, basal subtype, HER2 subtype, and normal subtype, we found that luminal B, basal and normal subtype have clear metabolic task tendency (Figure 2).

FIGURE 2

Figure 2. Schematic representation of the three-task of breast cancer and enrichment analysis of breast cancer subtypes. (A) The vertices in the figure represent the best metabolic tasks. The ellipse represents the error value calculated by the bootstrap method. The dots represent patient samples, and the color represents the clinical molecular subtype of PAM50. See the legend at the top right. (B–D) The enrichment distribution relative to the three metabolic tasks of basal, luminal B, and normal type. The horizontal axis represents the Euclidean distance ranking of the data bin (see section “Materials and Methods” for definition) and the metabolic task, and the vertical axis is the degree of enrichment, which is defined as the ratio of the density of the subtype in this bin to the overall density.

Luminal B accounts for about 35% of breast cancer, and is estrogen receptor (ER) positive, and either human epidermal growth factor receptor 2 (HER2) positive or negative (24). We found that luminal B has dominated energy metabolism (Hypergeometric test P-value <1×10^–4). There are only two enriched kyoto encyclopedia of genes and genomes (KEGG) metabolic pathways of Luminal B subtype, namely steroid synthesis and oxidative phosphorylation (Supplementary Table 2). That is, the gene expression levels of these two metabolic pathways are significantly higher than other metabolic tasks. Steroid hormones can exert their mitogenic effects by binding to estrogen, progesterone and androgen receptors (25). Targeting steroid synthesis pathway to control the development of luminal B cancer has been clinically applied, such as the treatment using anastrozole, a competitive inhibitor of androgen synthesis (26). Oxidative phosphorylation is related to energy production and cell metabolism intermediate production. Experiments at the cellular level have shown that luminal B cell line is more active in oxidative phosphorylation, and inhibiting the upstream regulatory pathway of oxidative phosphorylation, i.e., mTOR pathway, can reduce cell respiration, while has no effect on basal subtype (27).

Basal subtype is also known as triple-negative breast cancer, the overlap of the two subtypes is about 70–80% (28). Basal subtype accounts for about 10–20% of breast cancer with poor survival rate (29). Basal subtype tends to accomplish rapid proliferation (Hypergeometric test P-value < 1 × 10^–4). It enriches six KEGG metabolic pathways, including folic acid synthesis, oxidative phosphorylation, glyoxylic acid and dicarboxylic acid metabolism, etc (Supplementary Table 2). Folic acid is a carbon donor for one-carbon metabolism. Folic acid supports the production of NADPH, nucleotide production and methylation. Cancer cells up-regulate the folic acid metabolism pathway that is related to DNA production and cell growth (30). Basal subtype also enriches oxidative phosphorylation pathway. Clinical studies have shown that targeting oxidative phosphorylation provides an effective way to treat basal subtype (31, 32). Glyoxylic acid and dicarboxylic acid pathways are highly expressed in basal breast cancer (33).

Normal cells tend to fulfill the metabolic tasks other than cancer tasks (Hypergeometric test P-value < 1 × 10^–4). Interestingly, there are 60 gene sets enriched for normal cells metabolic tasks, far exceeding luminal B and basal. The metabolic tasks of normal cells should be maintaining metabolic homeostasis, such as the balance of redox potential and the stable rhythm of the circadian clock.

We performed the same analysis for lung cancer (Figure 3). There are two pathological types and the normal type of lung cancer samples. The pathological types are lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC).

FIGURE 3

Figure 3. Schematic representation of the three-task of lung cancer and enrichment analysis of lung cancer subtypes. (A) The vertices in the figure represent the best metabolic tasks. The ellipse represents the error value calculated by the bootstrap method. The dots represent patient samples, and the color represents the clinical molecular subtype of PAM50. See the legend at the top right. (B,C) The enrichment distribution relative to the three metabolic tasks of LUSC and normal type. The horizontal axis represents the Euclidean distance ranking of the data bin (see section “Materials and Methods” for definition) and the metabolic task, and the vertical axis is the degree of enrichment, which is defined as the ratio of the density of the subtype in this bin to the overall density.

Through the enrichment analysis of the KEGG metabolic pathway for lung cancer metabolic tasks, we found that LUSC cancer cells tend to complete the metabolic task 3 (Hypergeometric test P-value <1×10^–4), but with no significantly enriched KEGG metabolic pathways. However, the metabolism task 1, which normal cells tend to accomplish (Hypergeometric test P-value < 1 × 10^–4), has 36 KEGG metabolic pathways significantly enriched. As no enriched KEGG metabolic pathways were found in LUSC, in order to understand the biological significance of the metabolic task of cancer cells, we alternatively conducted enrichment analysis of the REACTOME metabolic pathways and found that the metabolic tasks of LUSC enriched in six pathways, including basic immunoglobulin interactions, vascular wall cell surface interactions, and amino acids synthesis and transamination, glycolysis, Gastrin-CREB signaling pathways through PKC and MAPK, and purine nucleoside monophosphate biosynthesis, which as a whole, are related to biomass production and energy production. Normal cellular metabolism tasks are enriched in 33 REACTOME pathways, and are related to the maintenance of tissue function and homeostasis. No highly enriched pathways are found for LUAD.

The reason that there are more metabolic pathways enriched in normal cells than cancer cells is largely because normal cells need to maintain function and body homeostasis, while cancer cells do not have clear physiological functions, but only to meet the monolithic requirements of proliferation and migration.

Multi-Objective Metabolic Model Outperforms Classical Objective Function in Gene Knockout Simulations

We use MOMM to simulate gene knockout to predict cell state. In order to verify the predictive power of the model, we compare the MOMM prediction results with the maximization of biomass production model.

We first use Kaplan–Meier online analyzer² to analyze the effect of gene expression on the survival probability of breast cancer patients. This online tool collected the relationship between the expression levels of 22,277 genes in 1,809 patients and the prognosis of cancer patients (34). We selected genes that have effect on cells after knockout, and compared the ROC curves of the two models. The selection criteria of genes and the calculation method of ROC curve are shown in section “Materials and Methods.”

We used breast cancer basal and luminal B subtype specific MOMM to simulate the impact of all metabolic gene knockout on cells, and the ratio of the deviation degree of the objective function value was used to express the degree of the gene knockout effect. We found that our model performed better than the maximize biomass model (Figures 4A,B).

FIGURE 4

Figure 4. Comparison between MOMM and the maximization of biomass model. (A,B) The survival comparison of basal and luminal B, the green line represents MOMM, and the red line represents the biomass model. The shade of the color represents the degree of influence on the objective function, the darker the color, the stronger the influence, and the 10% influence degree of the original function is used as the threshold. (C–H) Comparison of the accuracy of MOMM and the biomass model on the CRISPR-Cas9 knockout gene database DepMap.

Subsequently, we also used the CRISPR-Cas9 gene knockout database of cell line data to further compare MOMM and biomass model. We obtained the Boolean relationship matrix between cell lines and gene essentiality from the DepMap database³. The essentiality means that after knocking out with CRISPR-Cas9, it has a significant impact on the adaptability of the cell. We define that for the same type of cell line, if a gene is necessary in at least one cell line, this gene is also necessary for this type of cancer. See the Supplementary Table 3 for the cell type classification of cell lines.

We used the model to predict the basal and luminal B subtype, lung adenocarcinoma in lung cancer, liver cancer in digestive organ cancer, colon cancer in digestive tract cancer, and ovarian cancer in reproductive organ cancer. We found that MOMM outperforms the maximization of biomass production in all these cancer types (Figures 4C–H).

Exploration of Potential Synthetic Lethality Gene Pairs of Breast Cancer Subtypes

Synthetic lethality refers to the phenomenon that when two genes are interfered at the same time, it will cause cell death, while single ones will not be lethal (35). The interference includes loss of function mutation, RNA interference (RNAi), and drug therapy.

We used MOMM to perform pairwise knockout simulations for all reactions, to calculate the value of the objective function after double knockout, and used the highest single agent method to determine whether it is a valid combination (see section “Materials and Methods”). For luminal B subtype, there are 9,854 synthetical lethality pairs (SLPs). For basal subtype, there are 11,198 SLPs. Tables 1, 2 list the top 20 SLPs for basal and luminal B, respectively.

TABLE 1

Table 1. The top 20 SLPs of basal breast cancer.

TABLE 2

Table 2. The top 20 SLPs of luminal B breast cancer.

We found that the top 1 and 2 SLPs in both basal and luminal B are the same, one is glutathione oxidoreductase (also named as glutathione reductase, GR) and cytochrome c oxidase, mitochondrial complex IV (COX), the other is L-lactate dehydrogenase (LDH) and ubiquinol-6 cytochrome c reductase, mitochondrial complex III (UQCR).

Glutathione reductase is an antioxidant that catalyzes the reduction of glutathione disulfide to glutathione, and its high expression is related to the resistance of cancer cells to oxidative stress (36). COX is the last enzyme in the mitochondrial respiratory chain. It is the oxygen receptor of the respiratory chain and catalyzes the reduction of oxygen to water. Anthracyclines such as daunomycin (DAU) and doxorubicin (DOX) are breast cancer chemotherapy drugs, and their side effects are that they can cause cardiomyopathy during long-term treatment (37, 38). Both DAU and DOX can inhibit COX and promote the production of ROS, and the inhibitory effect is related to the dose of the drug (39). Flavonoids are inhibitors of glutathione reductase (40). Studies have shown that flavonoid inhibitors can reduce the cardiotoxic side effects of anthracyclines (41). Breast cancer cell line experiments have shown that the combined use of flavonoid drug quercetin and anthracycline DOX has better anti-tumor effect on highly aggressive breast cancer cells than the effect of single use, and weakens the side effects of DOX on non-tumor cells (42). Our results provide another theoretical basis for the combined use of flavonoids and anthracyclines in cancer treatment.

Lactate dehydrogenase catalyzes the conversion of lactic acid to pyruvate, and is usually highly expressed in breast cancer (27). LDH is located in the mitochondrial matrix. UQCR is overexpressed in breast cancer, and its knockdown can reduce the aggressiveness of breast cancer (43). Jeong et al. found that after inhibiting LDH in Chinese hamster ovary (CHO) cells, the respiration rate of the cells increased and they were more sensitive to UQCR inhibitor Antimycin A. After inhibiting both LDH and UQCR, cell activity decreased more than UQCR alone was inhibited (44).

We also counted the frequency of metabolic reactions in the SLPs, and found that mitochondrial energy metabolism-related enzymes occur most frequently, including ATP synthase (four protons for one ATP), acetyl-CoA carboxylase, cytochrome c oxidase, NADH dehydrogenase, ubiquinol-6 cytochrome c reductase, pyruvate dehydrogenase, acetate CoA-transferase. This suggests the importance of mitochondrial-related enzymes in combination medication. In fact, the inhibition of mitochondria already has some applications in tumor treatment. Wang et al. combined treatments targeting mitochondria and radiotherapy to reduce the growth of multidrug resistant tumors without significant systemic toxicity (45). Inhibition of mitochondria can also overcome dose toxicity. In addition to reduce cardiotoxicity of anthracyclines, it can also reduce the systemic toxicity of glycolysis inhibitors. Cheng et al. used mitochondrial inhibitors and glycolysis inhibitor 2-deoxy-D-glucose treatment of basal subtype of human breast cancer xenograft model, and found that the tumor size was significantly reduced, and the kidney, liver, and heart had no obvious morphological changes (46).

Other lethal combinations for which there is no direct experimental or clinical evidence may provide new guidelines for combination drug design. Most metabolic pathways also play important roles in normal cells, so the design of combination drugs may reduce the dosage of single drugs and avoid serious side effects.

We also compared the ranking of SLPs in single-target models. In most cases, the two targets in one top-ranking SLP pair do not rank high at the same time in single-target analysis. For example, combination of GR and COX of the Basal subtype is the top 1 SLP, while GR ranks 214th and COX ranks 4th in Basal single-target analysis. For the Luminal B subtype, combination of GR and COX are also the top 1 SLP, while GR ranks 139th and COX ranks 5th in the single-target analysis of Luminal B.

The Number of Activated Synthetic Lethality Gene Pairs as a Marker for Relapse-Free Survival

We analyzed the correlation between SLPs and breast cancer patient survival to evaluate the model. Patients data come from Györffy et al. (21). It should be pointed out that we did not select TCGA data that used to infer metabolic tasks to avoid circular reasoning. The data set contains 230 basal patients and 265 luminal B patients. We define activated SLPs as SLPs with low expression of gene A and high expression of gene B. We calculated the number of activated SLPs in patients with different subtypes. When the number of activated SLPs in the patient is higher than three-quarters of the number of activated SLPs in all patients of this subtype, we define that the patient has a high level of activated SLPs. When the patient’s activated SLPs is less than a quarter of all patients of this subtype, we define that the patient has a low level of activated SLPs.

We found that the number of activated SLPs can be used as an indicator of relapse free survival (RFS). Patients with higher survival probability have less activated SLPs, possibly due to less severe disorders of the patients’ metabolic network. While patients with lower survival probability have more activated SLPs (Basal: HR = 2.54, P-value = 1.13 × 10^–3; Luminal B: HR = 2.12, P-value = 1.15 × 10^–2; Figure 5).

FIGURE 5

Figure 5. The survival curve of the high and low SLP level. (A) Basal subtype; (B) Luminal B subtype.

In addition, the stricter the threshold for judging the number of activated SLPs, the more obvious the predictive effect on survival (Figure 6). HR > 1 means that patients are more likely to have cancer recurrence. The larger the HR is, the more dangerous situation the patients are in. When we decrease the threshold for distinguishing the number of activated SLPs, the value of HR will be lower, which proves the predictive effect of our model. If we choose 10% as threshold, i.e., high-level SLPs are the top 10% and low-level SLPs are the bottom 10% of all patients, HR is 3.99 (P-value = 8.19 × 10^–4). HR is 2.53 (P-value = 1.13 × 10^–3) when 20% is the threshold; HR is 2.40 (P-value = 2.77 × 10^–4) when 30% is the threshold; and HR is 1.81 (P-value = 1.64 × 10^–3) when 50% is the threshold.

FIGURE 6

Figure 6. The determination of the relationship between the threshold of activated SLP levels and survival prediction in basal subtype. (A–D) The thresholds for the number of high-level SLPs are the top 10, 20, 30, and 50%, and the thresholds for the number of low-level SLPs are the bottom 10, 20, 30, and 50% of all patients.

Discussion

Due to the lack of fluxomics data, we used transcriptomics data instead to defer the metabolic tasks. Although transcriptomics data and fluxomics data are positively correlated, there are discrepancies due to the presence of post-transcriptional regulation and translation process regulation (47, 48). More experimental fluxomics data in the future will help to improve the current analysis.

Metabolic heterogeneity not only exists among different cancer subtypes, but also among patients within the same cancer subtype. When patient-specific transcriptomics data are available, personalized drug targets and treatment plans can be predicted. Along with the rapid development of single cell technologies, heterogeneity of cancer cells within one cancer tissue of a patient can be analyzed and treatment plans considering the diversity and spatial arrangement of the cancer cells can be given. Moreover, the tumor-associated immune cells also have metabolic reprogramming during tumor progression to facilitate the escape of tumor cells from immune surveillance (49, 50). Our metabolic tasks analysis method may also be used to uncover the metabolic regulators in the tumor-associated immune cells in future studies.

Metabolic reprogramming is a key feature and vulnerable point of cancer cells. Predictive mathematical modeling can incorporate multiple levels of data to enable a better understanding of metabolic reprogramming in cancer cells and uncover therapeutic targets (51–53).

The heterogeneity of cancer cell metabolism due to different metabolic tasks increases the complexity of research and cancer drug discovery. In this study, we used ParTI method to analyze the metabolic tasks of breast cancer, liver cancer, digestive organ cancer, digestive tract cancer and reproductive system cancer. We found that different cancer types have different metabolic task tendencies with different biological implications. The metabolic tasks of cancer cells are more “uniform” compared to normal cells. The metabolic tasks of different sub-types of the same cancer are also different. For example, the luminal B subtype of breast cancer is more prone to energy metabolism, and the basal subtype is more prone to rapid proliferation. The anticancer target prediction by MOMM outperforms the conventional biomass maximization objective as evaluated by patient survival analysis.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: http://cancergenome.nih.gov/.

Author Contributions

SG and LL conceived and designed the project, developed the model, analyzed the data, and wrote the manuscript. ZD and HX contributed analytic tools. LL supervised the whole project. All authors contributed to the article and approved the submitted version.

Funding

This study was supported in part by the Ministry of Science and Technology (2016YFA0502303), National Natural Science Foundation of China (21633001), and the Chinese Academy of Medical Sciences (2021-I2M-5-014).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Prof. Minghua Deng (Peking University, Beijing, China), all members in LL’s group for their helpful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmed.2022.872024/full#supplementary-material

Footnotes

1. ^ http://xena.ucsc.edu/
2. ^ http://kmplot.com/
3. ^ https://score.depmap.sanger.ac.uk/

References

1. Ward PS, Thompson CB. Metabolic reprogramming: a cancer hallmark even warburg did not anticipate. Cancer Cell. (2012) 21:297–308. doi: 10.1016/j.ccr.2012.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Vander Heiden MG. Targeting cancer metabolism: a therapeutic window opens. Nat Rev Drug Discov. (2011) 10:671–84. doi: 10.1038/nrd3504

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Warburg O. On the origin of cancer cells. Science. (1956) 123:309–14.

Google Scholar

4. DeBerardinis RJ, Chandel NS. Fundamentals of cancer metabolism. Sci Adv. (2016) 2:e1600200. doi: 10.1126/sciadv.1600200

PubMed Abstract | CrossRef Full Text | Google Scholar

5. DelNero P, Hopkins BD, Cantley LC, Fischbach C. Cancer metabolism gets physical. Sci Transl Med. (2018) 10:1–8. doi: 10.1126/scitranslmed.aaq1011

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Hu J, Locasale JW, Bielas JH, O’Sullivan J, Sheahan K, Cantley LC, et al. Heterogeneity of tumor-induced gene expression changes in the human metabolic network. Nat Biotechnol. (2013) 31:522–9. doi: 10.1038/nbt.2530

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Jimenez CR, Zhang H, Kinsinger CR, Nice EC. The cancer proteomic landscape and the HUPO cancer proteome project. Clin Proteomics. (2018) 15:1–7. doi: 10.1186/s12014-018-9180-6

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kaushik AK, DeBerardinis RJ. Applications of metabolomics to study cancer metabolism. Biochim Biophys Acta Rev Cancer. (2018) 1870:2–14. doi: 10.1016/j.physbeh.2017.03.040

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ogrodzinski MP, Bernard JJ, Lunt SY. Deciphering metabolic rewiring in breast cancer subtypes. Transl Res. (2017) 189:105–22. doi: 10.1016/j.trsl.2017.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Schuetz R, Kuepfer L, Sauer U. Systematic evaluation of objective functions for predicting intracellular fluxes in Escherichia coli. Mol Syst Biol. (2007) 3:119. doi: 10.1038/msb4100162

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Schuetz R, Zamboni N, Zampieri M, Heinemann M, Sauer U. Multidimensional optimality of microbial metabolism. Science. (2012) 336:601–4. doi: 10.1126/science.1216882

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Shoval O, Sheftel H, Shinar G, Hart Y, Ramote O, Mayo A, et al. Evolutionary trade-offs, Pareto optimality, and the geometry of phenotype space. Science. (2012) 336:1157–60. doi: 10.1126/science.1217405

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Korem Y, Szekely P, Hart Y, Sheftel H, Hausser J, Mayo A, et al. Geometry of the gene expression space of individual cells. PLoS Comput Biol. (2015) 11:e1004224. doi: 10.1371/journal.pcbi.1004224

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Dai Z, Yang S, Xu L, Hu H, Liao K, Wang J, et al. Identification of cancer–associated metabolic vulnerabilities by modeling multi-objective optimality in metabolism. Cell Commun Signal. (2019) 17:124–38. doi: 10.1186/s12964-019-0439-y

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours. Nature. (2012) 490:61–70. doi: 10.1038/nature11412

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Duarte N, Becker S, Jamshidi N, Thiele I, Mo ML, Vo TD, et al. Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proc Natl Acad Sci USA. (2007) 104:1777–82. doi: 10.1073/pnas.0610772104

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hart Y, Sheftel H, Hausser J, Szekely P, Ben-Moshe NB, Korem Y, et al. Inferring biological tasks using Pareto analysis of high-dimensional data. Nat Methods. (2015) 12:233–5. doi: 10.1038/nmeth.3254

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Song HS, Reifman J, Wallqvist A. Prediction of metabolic flux distribution from gene expression data based on the flux minimization principle. PLoS One. (2014) 9:e112524. doi: 10.1371/journal.pone.0112524

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Segrè D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA. (2002) 99:15112–7. doi: 10.1073/pnas.232349399

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Berenbaum MC. What is synergy? Pharmacol Rev. (1989) 41:93–141.

Google Scholar

21. Györffy B, Lanczky A, Eklund AC, Denkert C, Budczies J, Li Q, et al. An online survival analysis tool to rapidly assess the effect of 22,277 genes on breast cancer prognosis using microarray data of 1,809 patients. Breast Cancer Res Treat. (2010) 123:725–31. doi: 10.1007/s10549-009-0674-9

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Mihály Z, Gyõrffy B. Improving pathological assessment of breast cancer by employing array-based transcriptome analysis. Microarrays. (2013) 2:228–42. doi: 10.3390/microarrays2030228

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Hausser J, Szekely P, Bar N, Zimmer A, Sheftel H, Caldas C, et al. Tumor diversity and the trade-off between universal cancer tasks. Nat Commun. (2019) 10:5423. doi: 10.1038/s41467-019-13195-1

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Cejalvo JM, De Dueñas EM, Galván P, García-Recio S, Gasión OB, Paré L, et al. Intrinsic subtypes and gene expression profiles in primary and metastatic breast cancer. Cancer Res. (2017) 77:2213–21. doi: 10.1158/0008-5472.CAN-16-2717

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Africander D, Storbeck KH. Steroid metabolism in breast cancer: where are we and what are we missing? Mol Cell Endocrinol. (2016) 466:86–97. doi: 10.1016/j.mce.2017.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Dunbier AK, Anderson H, Ghazoui Z, Salter J, Parker JS, Perou CM, et al. Association between breast cancer subtypes and response to neoadjuvant anastrozole. Steroids. (2011) 76:736–40. doi: 10.1016/j.steroids.2011.02.025

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Pelicano H, Zhang W, Liu J, Hammoudi N, Dai J, Xu RH, et al. Mitochondrial dysfunction in some triple-negative breast cancer cell lines: role of mTOR pathway and therapeutic potential. Breast Cancer Res. (2014) 16:434. doi: 10.1186/s13058-014-0434-6

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kathryn JC, Sireesha VG, Stanley L. Triple negative breast cancer cell lines: one tool in the search for better treatment of triple negative breast cancer. Breast Dis. (2012) 32:35–48. doi: 10.3233/BD-2010-0307.Triple

CrossRef Full Text | Google Scholar

29. Neophytou C, Boutsikos P, Papageorgis P. Molecular mechanisms and emerging therapeutic targets of triple-negative breast cancer metastasis. Front Oncol. (2018) 8:31. doi: 10.3389/fonc.2018.00031

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wang YP, Lei QY. Perspectives of reprogramming breast cancer metabolism. Adv Exp Med Biol. (2017) 1026:217–32. doi: 10.1007/978-981-10-6020-5_10

CrossRef Full Text | Google Scholar

31. Kanchi M, Hirpara J, Sachaphibulkij K, Tan T, Dietzel H, Lim L, et al. Targeting mitochondrial function for the treatment of triple negative breast cancer: development of a small molecule inhibitor against mitochondrial STAT3. In: Proceedings of the 2018 San Antonio Breast Cancer Symposium. (Philadelphia, PA: American Association for Cancer Research (AACR)) (2019). p. 6–18–20.

Google Scholar

32. Meric-Bernstam F, Evans K, Zheng X, Su X, Yuca E, Scott S, et al. Oxidative phosphorylation as a target in triple negative breast cancer therapy. In: Proceedings of the American Association for Cancer Research Annual Meeting 2017. (Philadelphia, PA: American Association for Cancer Research (AACR)) (2017). p. 4970–4970.

Google Scholar

33. Chen WY, Wu F, You ZY, Zhang ZM, Guo YL, Zhong LX. Analyzing the differentially expressed genes and pathway cross-talk in aggressive breast cancer. J Obstet Gynaecol Res. (2015) 41:132–40. doi: 10.1111/jog.12495

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Nagy Á, Lánczky A, Menyhárt O, Gyorffy B. Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets. Sci Rep. (2018) 8:9227. doi: 10.1038/s41598-018-27521-y

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Gao S, Lai L. Synthetic lethality in drug development: the dawn is coming. Future Med Chem. (2018) 10:2129–32. doi: 10.4155/fmc-2018-0227

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Freitas M, Baldeiras I, Proença T, Alves V, Mota-Pinto A, Sarmento-Ribeiro A. Oxidative stress adaptation in aggressive prostate cancer may be counteracted by the reduction of glutathione reductase. FEBS Open Bio. (2012) 2:119–28. doi: 10.1016/j.fob.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Lovitt CJ, Shelper TB, Avery VM. Doxorubicin resistance in breast cancer cells is mediated by extracellular matrix proteins. BMC Cancer. (2018) 18:41. doi: 10.1186/s12885-017-3953-6

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Guo J, Zhou J, Ying X, Men Y, Li RJ, Zhang Y, et al. Effects of stealth liposomal daunorubicin plus tamoxifen on the breast cancer and cancer stem cells. J Pharm Pharm Sci. (2010) 13:136–51. doi: 10.18433/j3p88z

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Srinivasan S, Avadhani NG. Cytochrome c oxidase dysfunction in oxidative stress. Free Radic Biol Med. (2012) 53:1252–63. doi: 10.1161/CIRCULATIONAHA.110.956839

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Paulíková H, Berczeliová E. The effect of quercetin and galangin on glutathione reductase. Biomed Pap Med Fac Univ Palacky Olomouc Czech Repub. (2005) 149:497–500. doi: 10.5507/bp.2005.089

CrossRef Full Text | Google Scholar

41. Bast A, Kaiserová H, Den Hartog GJM, Haenen GRMM, Van Der Vijgh WJF. Protectors against doxorubicin-induced cardiotoxicity: flavonoids. Cell Biol Toxicol. (2007) 23:39–47. doi: 10.1007/s10565-006-0139-4

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Staedler D, Idrizi E, Kenzaoui BH, Juillerat-Jeanneret L. Drug combinations with quercetin: doxorubicin plus quercetin in human breast cancer cells. Cancer Chemother Pharmacol. (2011) 68:1161–72. doi: 10.1007/s00280-011-1596-x

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Owens KM, Kulawiec M, Desouki MM, Vanniarajan A, Singh KK. Impaired OXPHOS complex III in breast cancer. PLoS One. (2011) 6:e23846. doi: 10.1371/journal.pone.0023846

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Jeong DW, Cho IT, Kim TS, Bae GW, Kim IH, Kim IY. Effects of lactate dehydrogenase suppression and glycerol-3-phosphate dehydrogenase overexpression on cellular metabolism. Mol Cell Biochem. (2006) 284:1–8. doi: 10.1007/s11010-005-9004-7

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang H, Gao Z, Liu X, Agarwal P, Zhao S, Conroy DW, et al. Targeted production of reactive oxygen species in mitochondria to overcome cancer drug resistance. Nat Commun. (2018) 9:562. doi: 10.1038/s41467-018-02915-8

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Cheng G, Zielonka J, Dranka BP, McAllister D, Mackinnon AC, Joseph J, et al. Mitochondria-targeted drugs synergize with 2-deoxyglucose to trigger breast cancer cell death. Cancer Res. (2012) 72:2634–44. doi: 10.1158/0008-5472.CAN-11-3928

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Maier T, Güell M, Serrano L. Correlation of mRNA and protein in complex biological samples. FEBS Lett. (2009) 583:3966–73. doi: 10.1016/j.febslet.2009.10.036

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Vogel C, De Sousa Abreu R, Ko D, Le SY, Shapiro BA, Burns SC, et al. Sequence signatures and mRNA concentration can explain two-thirds of protein abundance variation in a human cell line. Mol Syst Biol. (2010) 6:400. doi: 10.1038/msb.2010.59

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Xu Y, He L, Fu Q, Hu J. Metabolic reprogramming in the tumor microenvironment with immunocytes and immune checkpoints. Front Oncol. (2021) 11:4775. doi: 10.3389/FONC.2021.759015/BIBTEX

CrossRef Full Text | Google Scholar

50. Biswas SK. Metabolic reprogramming of immune cells in cancer progression. Immunity. (2015) 43:435–49. doi: 10.1016/j.immuni.2015.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Roy M, Finley SD. Computational model predicts the effects of targeting cellular metabolism in pancreatic cancer. Front Physiol. (2017) 8:217. doi: 10.3389/FPHYS.2017.00217/BIBTEX

CrossRef Full Text | Google Scholar

52. Gorick CM, Saucerman JJ, Price RJ. Computational model of brain endothelial cell signaling pathways predicts therapeutic targets for cerebral pathologies. J Mol Cell Cardiol. (2022) 164:17–28. doi: 10.1016/J.YJMCC.2021.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Nobile MS, Votta G, Palorini R, Spolaor S, De Vitto H, Cazzaniga P, et al. Fuzzy modeling and global optimization to predict novel therapeutic targets in cancer cells. Bioinformatics. (2020) 36:2181–8. doi: 10.1093/BIOINFORMATICS/BTZ868

PubMed Abstract | CrossRef Full Text | Google Scholar