Metabolic Network-Based Identification and Prioritization of Anticancer Targets Based on Expression Data in Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) is a deadly form of liver cancer with high mortality worldwide. Unfortunately, the large heterogeneity of this disease makes it difficult to develop effective treatment strategies. Cellular network analyses have been employed to study heterogeneity in cancer, and to identify potential therapeutic targets. However, the existing approaches do not consider metabolic growth requirements, i.e., biological network functionality, to rank candidate targets while preventing toxicity to non-cancerous tissues. Here, we developed an algorithm to overcome these issues based on integration of gene expression data, genome-scale metabolic models, network controllability, and dispensability, as well as toxicity analysis. This method thus predicts and ranks potential anticancer non-toxic controlling metabolite and gene targets. Our algorithm encompasses both objective-driven and—independent tasks, and uses network topology to finally rank the predicted therapeutic targets. We employed this algorithm to the analysis of transcriptomic data for 50 HCC patients with both cancerous and non-cancerous samples. We identified several potential targets that would prevent cell growth, including 74 anticancer metabolites, and 3 gene targets (PRKACA, PGS1, and CRLS1). The predicted anticancer metabolites showed good agreement with existing FDA-approved cancer drugs, and the 3 genes were experimentally validated by performing experiments in HepG2 and Hep3B liver cancer cell lines. Our observations indicate that our novel approach successfully identifies therapeutic targets for effective treatment of cancer. This approach may also be applied to any cancer type that has tumor and non-tumor gene or protein expression data.

Introduction

Hepatocellular carcinoma (HCC) is a primary form of liver cancer and one of the main causes of cancer mortality globally (Ferlay et al., 2015). Patient prognosis is usually poor and most result in patient death (Altekruse et al., 2009). This high mortality is partly due to the high tumor heterogeneity of this cancer (Friemel et al., 2015; Benfeitas et al., 2017), which makes it difficult to identify suitable and effective therapeutic targets. It is therefore urgent to devise suitable approaches capable of capturing the high tumor heterogeneity for identifying effective therapeutic targets for this cancer.

Biological network analysis is emerging as an efficient and suitable in silico approach to identify cellular targets for disease treatment incorporating cellular heterogeneity (Barabasi et al., 2011). For instance, network topology features such as centralized protein hubs have been employed in the analysis of protein-protein interaction (PPI) networks to identify potential therapeutic targets (Guney et al., 2016; Lv et al., 2016). A few efforts have applied network controllability to identify minimum sets of driver proteins for controlling PPI networks (Liu et al., 2011; Yuan et al., 2013), or indispensable proteins from a network controllability perspective (Vinayagam et al., 2016). Central and highly connected proteins and genes appear in bottleneck interactions (Wuchty, 2014), and tend to be essential from a lethality perspective (Jeong et al., 2001; Yu et al., 2008; Najafi et al., 2014). However, these approaches have a limited scope in their considered essentiality analyses and do not take into account comprehensive descriptions of biological functionality, or the stoichiometry of the interactions.

Genome-scale metabolic models (GEMs) are whole cell stoichiometric representations of metabolism that take into account network functionality through prediction of a model's capacity to achieve one or more biological tasks (e.g., biomolecule synthesis and growth) (Mardinoglu and Nielsen, 2012, 2015; Mardinoglu et al., 2013). GEMs have been employed to identify essential genes or anticancer metabolites, build tissue-specific cellular characterizations, and explain cell behavior at the metabolic, signaling, gene, and protein level (Folger et al., 2011; Agren et al., 2014; Mardinoglu et al., 2014a,b; Varemo et al., 2015). However, GEM-based methods often do not prioritize between different candidate therapeutic targets, and assessing toxicity to normal tissues is not always possible if the essential nodes are not contemplated by the flux distribution. Further, these approaches often do not allow for ranking between different candidate targets.

Here, we overcome the limitations of current state-of-the-art methods by introducing a network-based prioritization approach to identify and prioritize non-toxic metabolic targets for disease treatment using network controllability, topology analysis, and constraint-based modeling techniques. Our algorithm approach permits analyzing cellular network behavior using any kind of expression data such as microarray, transcriptomics (RNA-seq), or protein expression. Briefly, this algorithm combines GEMs, personalized metabolite-metabolite and reaction-reaction association networks. Based on the analysis of these networks, it determines metabolites/genes whose perturbation has a strong effect on network dynamics (i.e., minimum driver set nodes) (Yuan et al., 2013), indispensable metabolites/genes (Vinayagam et al., 2016), and several network topology parameters (i.e., node centrality analysis). Based on this information, this algorithm determines and ranks anticancer metabolites/genes, excluding those that would also become toxic for non-cancerous tissues.

We employed this algorithm in the analysis of RNA-seq data of hepatocellular carcinoma (HCC), the most prevalent form of liver cancer and third leading cause of cancer-related mortality worldwide (Ferlay et al., 2010). Through integration of transcriptomic data from cancerous and non-cancerous samples and personalized systems biology approaches, we identified cancer-specific targets. These targets were experimentally validated in HepG2 and Hep3B cells, cell lines frequently used as models for liver cancer.

Methods

Algorithm Overview and Requirements

The algorithm presented here requires Matlab R2016a, RAVEN and tINIT (Agren et al., 2013, 2014) and has been implemented on a Dell laptop running Windows 7. It comprises the 4 major steps detailed below: 1. Objective-dependent identification of potential targets; 2. Creation of personalized metabolic and reaction networks; 3. Controllability analysis; and 4. Target prioritization. We provide an explanation for each of these steps below, as well as pseudo-code or original references showing their implementation.

Step 1. Objective-Dependent Identification of Potential Targets

We identify metabolites and genes that potentially serve as therapeutic targets based on antimetabolites and in silico gene silencing of personalized GEMs. The two approaches rely on objective function and metabolic tasks for cancer and non-cancerous tissues to determine the functional outcome, i.e., viable vs. non-viable. Viability is determined for instance by testing whether objective-functions are non-null. Importantly, toxicity testing is done by identifying those antimetabolites/silenced genes that render tumors non-viable, but maintain viability in non-cancerous cells.

Step 1a. Antimetabolites are those structurally similar to endogenous metabolites and that prevent cell growth. Basically, for each metabolite and regardless of compartmentalization, an antimetabolite is identified by removing the reactions where the metabolite serves as substrate. This step was implemented using tINIT as previously indicated (Agren et al., 2014).

Step 1b. The concept of in silico gene silencing is similar to lethality analyses but considers metabolic tasks (Supplementary File 1) in addition to objective functions to assess model functionality, to assess the effect of gene silencing. In order to simulate the effect of gene silencing for each gene on the models, we removed the reaction(s) catalyzed by the gene.

ALGORITHM STEP 1b

Algorithm STEP 1b. In silico gene silencing

Step 1c. The predicted antimetabolites (Step 1a) and silenced genes (Step 1b) are filtered to select only those found in HCC GEMs but not in non-cancerous GEMs.

Step 2. Creation of Personalized Metabolic and Reaction Networks

Personalized GEMs are used to generate directed Metabolite-Metabolite network (MMN, Step 2a), and directed Reaction-Reaction network (RRN, Step 2b). MMN and RRN are simple graphs G(V, E) with V nodes and E edges. In MMN, nodes (metabolites) are connected by an edge if they are involved in the same reaction. In RRN, nodes (reactions) are connected if the product of a reaction serves as substrate by the other (Supplementary Figure 1).

Step 3. Controllability Analysis

We identify minimum driver set nodes (MDS), those responsible for controlling network dynamics upon perturbation (Yuan et al., 2013), as well as indispensable nodes (Vinayagam et al., 2016), from MMN (Step 2a) and RRN (Step 2b). MDS and indispensable nodes (metabolites or genes) are jointly used for subsequent steps, and termed controlling nodes.

Step 3a. We identify MDS in MMN and RRN according to the Popov–Belevitch–Hautus (PBH) rank condition (Sontag, 2013; Yuan et al., 2013).

ALGORITHM STEP 2a

Algorithm STEP 2a. Creation of directed metabolite-metabolite networks

Step 3b. Node dispensability is determined by assessing the maximum geometric multiplicity of the eigenvalue λm upon removing each node.

Step 3c. We excluded those nodes that were simultaneously identified in both HCC and non-cancerous networks in Steps 3a and 3b.

Step 4. Target Prioritization

In the final step, we filter nodes by selecting those with anticancer properties (Step 1, metabolic task dependent) that are also controlling metabolites/genes (Step 3, metabolic task independent). This generates a list of anticancer non-toxic controlling nodes that are then ranked based on their centrality properties.

Step 4a. Filter nodes simultaneously identified as controlling nodes (Step 3) with anticancer properties (Step 1).

Step 4b. We compute topological parameters for MMN and RRN. Specifically, we compute node degree and betweenness centrality using default commands [e.g., Matlab R2016a built-in degree() and centrality() functions].

Step 4c. Based on each node's topological measures, we rank anticancer controlling nodes. Indispensable nodes are ranked based on degree centrality whereas MDS are ranked based on betweenness centrality. This results in a list of nodes (metabolites/genes) with anticancer and controlling properties, sorted according to their topological importance in the networks.

ALGORITHM STEP 2b

Algorithm STEP 2b. Creation of directed reaction-reaction networks

Gene Expression Data and Modeling for HCC Patients

Transcriptomic data (RNA-seq) for 50 HCC and 50 adjacent non-cancerous liver samples were retrieved from NCI's Genome Data Commons (Grossman et al., 2016) as Fragments Per Kilobase of transcript per Million mapped reads (FPKM). Personalized GEMs were built by integrating transcriptomic data with the liver-specific GEMs through tINIT and RAVEN (Agren et al., 2013, 2014) using the reference human GEM HMR2 (Uhlen et al., 2017).

The following thresholds for gene levels were considered: no expression (FPKM < 1), low expression (1 ≤ FPKM < 10), medium expression (10 ≤ FPKM < 50), and high (FPKM ≥ 50). We considered biomass production or ATP consumption as objective functions for HCC and non-cancer GEMs, respectively. Model functionality of HCC and non-cancer GEMs was respectively determined based on 57 and 56 metabolic tasks, as previously determined (Agren et al., 2014). Cell growth is considered as an additional metabolic task for HCC.

HepG2 and Hep3B Cell Line Experiments

HepG2 and Hep3B, human hepatocellular carcinoma cell lines frequently used as models for liver cancer (Qiu et al., 2015), were cultured in RPMI-1640 Medium (R2405, Sigma-Aldrich) and Minimum Essential Medium Eagle (M4655, Sigma-Aldrich), respectively. Both cell lines were supplemented with 10% fetal bovine serum (FBS, F2442, Sigma-Aldrich) and incubated in 5% CO₂ humidity at 37°C. To confirm the effect of silencing of CRLS1, PRKACA and PGS1 genes on cancer cell growth and viability, we used RNA interference (RNAi) to inhibit gene expression. Cells were infected by predesigned targeted Silencer^® siRNAs for each gene (Thermo Fisher Scientific, clone IDs shown in Supplementary File 2), at 35 nM by using Lipofectamine^® RNAiMAX (13778075, Life Technologies). Cells incubated in medium with non-targeting negative control siRNA at 35 nM (4390843, Life Technologies) were assigned as control.

ALGORITHM STEP 3b

Algorithm STEP 3b. Determination of node dispensaiblity

Total RNA was isolated with RNeasy mini Kit (15596026, Thermo Fisher Scientific) after treatment with siRNA for 24 h. The expression profiles of key genes (CRLS1, PRKACA and PGS1) in HepG2 and Hep3B cells separately were measured and analyzed via quantitative real-time PCR with iTaq Universal SYBR Green One-Step Kit (1725151, Bio-Rad), using anchored oligo (dT) primer based on CFX96™ detection system (Bio-Rad). GAPDH was set as the internal control for normalization. Primer sequences are listed in Supplementary File 2. Cell Counting Kit-8 (CCK-8, Sigma-Aldrich) was used to detect variation in cell proliferation after interference by siRNA for 24 h.

All experiments were performed strictly according to the manufacturer's instructions and repeated in triplicate.

Statistics

Statistic comparison between siRNA-targeted gene expression and non-negative controls was performed using Welch's t-test, corrected for Benjamini & Hochberg.

Results

Algorithm

GEMs include metabolites (x), reversible and irreversible reactions (r), and associated genes (g). Reactions are represented as a stoichiometric matrix (S) of size m × n, where rows (m) in S represent metabolites and columns (n) correspond to reactions in the model (Supplementary Figure 1). Each entry is a stoichiometric coefficient s_ij of the metabolites that participate in a reaction. s_ij ≥ 1 and s_ij ≤ −1 respectively indicate that metabolite i is a product or reactant in reaction j, whereas s_ij = 0 indicates no involvement in the reaction. Reaction reversibility is discriminated by the rev vector, where 0 and 1 respectively indicate an irreversible or reversible reaction. Currency metabolites (x_curr) such as cofactors, coenzymes and H₂O are found in a number of reactions and considered by the GEMs, but are disregarded when creating metabolite-metabolite and reaction-reaction association networks(Khosraviani et al., 2016).

Our algorithm (Figure 1, Supplementary Figure 2) predicts potential metabolite and gene targets that display anticancer activity but are non-toxic to non-cancerous tissues. To do so, it integrates expression data (transcriptomics, microarray or protein levels) into GEMs (Methods, Algorithm steps 1b–2b) to generate cancer and non-cancer specific models. After network manipulation and analysis, it identifies and prioritizes anticancer metabolites and gene targets (Algorithm step 3b). These were implemented in Matlab, and require RAVEN and tINIT (Agren et al., 2013, 2014). Our hypothesis follows the determination of metabolite or genes that display anticancer and network-controlling properties, and sorting them according to their topology parameters. The highest ranking (most central) nodes are potentially better targets due to their pivotal and central role in network dynamics, and are excluded if identified in non-cancerous tissues to prevent toxicity. Several filtering steps are taken to exclude false positives: candidate anticancer targets must have network controlling properties and cannot be found in non-cancerous networks, i.e., they must be anticancer non-toxic controlling nodes.

FIGURE 1

Figure 1. Outline of network-based identification and prioritization of metabolite and gene anticancer targets.

This algorithm consists of 4 main steps, applied to both cancer and non-cancer personalized GEMs (details in section Methods), constructed using expression data and reference models (Uhlen et al., 2017). Step 1 is independent of steps 2 and 3, and all are integrated at Step 4.

(1) Objective-dependent identification of potential targets

a. Antimetabolite analysis

b. In silico gene silencing

c. Toxicity analysis

(2) Creation of personalized metabolite and reaction networks

a. Creating metabolite-metabolite networks

b. Creating reaction-reaction networks

(3) Controllability analysis

a. Identification of Minimum Driver Set Nodes

b. Dispensability analysis

c. Toxicity analysis

(4) Target prioritization

a. Determining anticancer non-toxic controlling nodes

b. Computation of centrality parameters

c. Target prioritization

All steps are described in detail according to a toy model (Supplementary Figures 1, 2) and pseudo-code is described in section Methods.

Determination of Antimetabolites, Anticancer Genes, and Networks Based on Personalized GEMs in HCC

As a test case, we sought to identify non-toxic cancer targets using 50 HCC and 50 non-cancerous transcriptomic samples from TCGA. This test case employs liver-specific GEMs (Uhlen et al., 2017), which encompass 7,780 reactions and 2,857 metabolites controlled by 2,892 genes. We built functional personalized GEMs for HCC and non-cancerous samples (Supplementary File 3) where model functionality was determined based on the metabolic tasks and the objective functions.

We then identified 374 antimetabolites out of 2,857 metabolites that inhibited growth of HCC GEMs, but passed all metabolic tasks in non-cancerous GEMs, i.e., non-toxic cancer-specific antimetabolites (Steps 1a and 1c, Supplementary File 4). Most antimetabolites are involved in carnitine shuttle, fatty acid activation, and acyl-CoA hydrolysis, as well as metabolism of cholesterol, glycerolipid, glycan, purine, pyrimidine, nucleotide, amino acid, and proteins (Figure 2A).

FIGURE 2

Figure 2. General MMN and RRN node features in HCC. (A) Number of antimetabolites identified per pathway. Note that a predicted antimetabolite may be found in multiple pathways. (B) Cumulative degree distribution for MMN (purple) and RRN (magenta), compared to randomly generated models (Erdos–Rényi, dots). RRNs show qualitatively similar properties (results not shown). (C) Mean degree for indispensable, neutral, and dispensable metabolites in noncancer and HCC MMNs. RRNs show qualitatively similar properties (results not shown). (D) Number of genes leading to lethality identified per pathway. Note that a predicted synthetic lethal gene may be found in multiple pathways.

Based on in silico gene silencing, we found that 2,204 genes did not fulfill all 56 metabolic tasks in non-cancerous models. We identified 283 genes that inhibited the growth in at least one HCC, and 8 genes that inhibit growth in all HCC models, but induced no change in non-cancer GEMs (Steps 1b and 1c, Figure 2D, Supplementary File 4).

We also generated MMN and RRN for each personalized GEM (Steps 2a and 2b, respectively). The generated MMN and RRN showed scale free (power-law) degree distribution, and the cumulative degree distribution of the networks were compared with random networks generated by Erdos–Rényi model (Figure 2B).

Identification of Non-toxic Anticancer Controlling Metabolites Using Metabolite-Metabolite Network Controllability

We identified cancer-exclusive controlling metabolites in all GEMs (Step 3, Supplementary File 4). Specifically, we found 201 MDS metabolites, most of which are involved in exchange and transporting reactions. Most amino acids, (deoxy)nucleotides, fatty acids and cholesterol, and some sugars including mannose, glucose, fructose, and galactose, are MDS in HCC models. In turn, we identify 578 indispensable metabolites, most of which either start metabolic pathways or are important metabolites in the networks, such as malonyl-CoA, palmitoyl-CoA, glucose, mannose, glucose-6-phosphate, HMG-CoA, (deoxy)nucleotides, DNA, acetyl-CoA, amino acids, L-carnitine, fatty acids, and coenzymes. Indispensable metabolites present higher degree in MMN (Supplementary File 4), in comparison with other node types in both non-cancer and cancer networks (Figure 2C). Because of this, the indispensable nodes are ranked based on node degree parameters (Step 4c). After removing metabolites simultaneously found in HCC and non-cancerous networks, we identify 142 cancer-specific controlling metabolites which were used in the subsequent steps (Supplementary File 4). The 374 antimetabolites found in Steps 1a and 1c were filtered based on the 142 cancer-specific controlling metabolites. This results in 74 non-toxic anticancer metabolites with network controlling properties, potential targets for antimetabolite treatment (Step 4a). Among the 74 anticancer controlling metabolites, 54 are indispensable and 20 are MDS (Supplementary File 4). The anticancer indispensable controlling metabolites shows substantially higher degree (Supplementary File 4). For instance, the top ranking indispensable antimetabolites show high degree (leucine 79, isoleucine 78, UMP 66, malonyl-CoA 34, palmitoyl-CoA 16) when compared with mean metabolite degree = 6. The 75 anticancer controlling metabolites are mostly involved in metabolism of DNA and nucleotides, amino acids, and fatty acids.

Importantly, while several nodes present very high centrality, they may not be suitable therapeutic targets because their silencing either does not lead to lethality, or leads to lethality in both cancer and non-cancer networks (i.e., toxicity, Supplementary Files 4, 5).

Overall, our network-based approaches highlight several potential metabolic targets for cancer-specific treatment.

Identification of Gene Targets Exclusive to HCC Using Reaction-Reaction Network and Experimental Validation

Silencing of controlling genes in RRN shows specific features in HCC and non-cancerous networks (Figure 3). Among the eight genes that inhibit growth in all HCC but not in non-cancerous models, we identify three genes that are controlling nodes in HCC RRN but not in non-cancerous RRN (Step 3, Supplementary File 5): PRKACA (protein kinase cAMP-activated catalytic subunit alpha), PGS1 (phosphatidylglycerophosphate synthase 1), CRLS1 (cardiolipin synthase 1). Since they are indispensable nodes, we then ranked them based on degree parameter as PRKACA, PGS1, and CRLS1 (Supplementary File 5), from most to least central genes.

FIGURE 3

Figure 3. Pivotal genes in HCC (Left) and non-cancerous (Right) RRN. Node size is proportional to degree centrality, and controlling anticancer genes are highlighted (red means higher number of patients where gene silencing leads to lethality and where a gene is a controlling node).

We tested the anticancer properties of these three genes using siRNA in HepG2 and Hep3B cells (Figure 4). Cell viability was tested 24 h after gene silencing using the CCK-8 assay (see Methods). As negative controls, cells were exposed to a non-targeting siRNA. Our observations indicate that PRKACA silencing led to 22–24% decrease in cell growth in both cell lines (P < 0.05). CRLS1 silencing led to 30% lower growth in HepG2 (P < 0.05) but no change in Hep3B. Finally, PGS1 silencing produced the most substantial growth change, decreasing growth in Hep3B by 35% (Q < 0.05).

FIGURE 4

Figure 4. siRNA inhibition of CRLS1, PGS1, and PRKACA expression in Hep3B and HepG2 liver cell lines leads to decreased growth after 24 h in comparison with controls. Comparisons performed using Welch's T-test. See methods for details.

Discussion

Network analyses have been extensively used in understanding cellular response to disease and predicting targets for treatment. Due to their comprehensive and integrative nature, networks successfully and efficiently capture multidimensional biological heterogeneity (Mardinoglu and Nielsen, 2012, 2015; Mardinoglu et al., 2013; Najafi et al., 2014; Benfeitas et al., 2017). Such cases include topology (Asgari et al., 2013), controllability (Liu et al., 2011; Yuan et al., 2013; Wuchty, 2014; Kanhaiya et al., 2017), dispensability (Vinayagam et al., 2016), and genome-scale metabolic models (Mardinoglu and Nielsen, 2012, 2015; Mardinoglu et al., 2013; Agren et al., 2014) in prediction of disease targets. For instance, a recent analysis shows that a small number of proteins essential for survival are controlled by a small set of driver nodes (Kanhaiya et al., 2017). Network controllability analyses have also highlighted that highly connected proteins are critical for cell survival (Jeong et al., 2001; Jonsson and Bates, 2006). In turn, others have used metabolic networks to identify candidate metabolite analogs and genes with anticancer properties (Agren et al., 2014). While these efforts have undoubtedly provided great contributions to clarify the underlying disease biological phenomena, and identified pivotal components for controlling network dynamics, such methods have specific shortcomings. For instance, topology-based approaches often do not take into account interaction stoichiometry, and essentiality analyses do not consider multiple functionality properties necessary for growth, e.g., energetic and redox balance, maintenance of central metabolism (see Agren et al., 2014 and references therein).

Here, we developed an algorithm that integrates expression data with both objective-dependent and objective-independent approaches to identify and rank metabolites and genes for cancer treatment. Objective dependent approaches such as flux balance analysis, employed here, base their predictions on satisfying specific biological phenomena. However, our method overcomes some of the disadvantages in FBA-based approaches which, unlike our method, do not permit ranking multiple candidate targets, and assessing potential toxic targets is not always possible if the essential nodes are not contemplated by the flux distribution. Objective independent approaches, including network topology, controllability, and dispensability analysis employed here, permit classifying and ranking between nodes. In our algorithms, candidate targets are filtered at multiple steps based on their toxicity, network controllability and anticancer properties. Thus, and unlike typical network topology analysis, our proposed method incorporates biological functionality expressed as metabolic requirements, and reaction stoichiometry. To our knowledge, this is the first such case where prioritization and classification of potential metabolite/gene targets is determined while considering toxicity to non-cancerous tissues and biological functionality.

We employed this algorithm for the analysis of hepatocellular carcinoma, a deadly form of liver cancer (Ferlay et al., 2010). This transcriptomic dataset comprised 50 tumor and 50 non-cancerous samples, retrieved form TCGA (Weinstein et al., 2013). We identify 142 metabolites, and 3 genes as potential therapeutic targets for non-toxic treatment of this liver cancer.

The non-toxic anticancer metabolites identified in here tend to be involved in metabolism of DNA and nucleotides, amino acids, and fatty acids. Several observations support utilization of analogs of these metabolites or targeting these biological processes for cancer treatment. For instance, several pyrimidine or purine analogs are among FDA-approved drugs for cancer therapy in DrugBank (Law et al., 2014), such as Gemcitabine, Fluorouracil, Clofarabine, Azacitidine, Floxuridine, Cladribine, Raltitrexed, Tioguanine, Azacitidine. All amino acids analogs were found as antimetabolites but most of them were filtered based on cancer controlling metabolites as toxic targets based on controllability analysis. Among all amino acids, glutamine, leucine, and isoleucine were introduced as nontoxic anticancer metabolite targets. Several studies show the critical role of amino acids such as glutamine in cancer (Xie et al., 2012; Bhutia et al., 2015; Xiao et al., 2016; Fung and Chan, 2017; Maddocks et al., 2017). For instance, leucine and isoleucine supplementation promote tumorigenesis in rat bladder cancer (Xie et al., 2012), and leucine deprivation inhibits breast cancer proliferation in humans (Xiao et al., 2016).

Other metabolites such as fatty acids also display a crucial role in cancer. For instance, malonyl-CoA serves as substrate for fatty acid synthase, and inhibition of this enzyme with the synthetic inhibitor C75 (Kuhajda et al., 2000) significantly decreases growth in HepG2 cells (Gao et al., 2006). Our previous work identified L-carnitine and metabolites involved in L-carnitine biosynthesis as non-toxic anticancer HCC targets using protein levels for 6 HCC samples (Agren et al., 2014). Our algorithm presented here improves on this work by showing that L-carnitine is a non-toxic anticancer target for HCC treatment (Step 2, Supplementary File 3) using a different technology and sample, and network topology analyses indicate that other metabolites may have better network controlling properties than L-carnitine.

At the gene level, only 3 genes pass all the filtering criteria, and show anticancer non-toxic controlling properties: protein kinase cAMP-activated catalytic subunit alpha (PRKACA), phosphatidylglycerophosphate synthase 1 (PGS1), and cardiolipin synthase 1 (CRLS1). We experimentally confirmed that the 3 genes display anticancer properties using siRNA in HepG2 and Hep3B liver cancer cell lines. We observe that depletion of PGS1 and PRKACA greatly decreases growth in both HepG2 and Hep3B, whereas depletion of CRLS1 decreases expression in HepG2. These genes have been associated with alterations in proliferation, malignancy and apoptosis in several cancers. The cAMP-dependent protein kinase A has crucial signaling roles in cancer, promotes resistance to therapy and its inhibition blocks tumor invasion (Sabbisetti et al., 2005; Honeyman et al., 2014; Moody et al., 2015). PGS1 metabolizes phosphatidylglycerol in lipid metabolism, is involved in cellular transformation, and its upregulated in metastases (Hirsch et al., 2010; Hartung et al., 2017). Finally, cardiolipin synthase 1 is involved in phosphatidylglycerol metabolism, and alterations in cardiolipin metabolism are associated with proliferation and cell death (Schug and Gottlieb, 2009; Sapandowski et al., 2015).

Overall, these observations indicate that our algorithm successfully identifies several metabolites and genes that may be targeted for cancer treatment. The observations in both MMNs and RRNs show the necessity for combining both topology- and objective-dependent analyses in identification of suitable therapeutic targets. These approaches may be employed together with other omics technologies, to identify and develop targeted therapeutic strategies in other cancers.

Author Contributions

GB and AM conceptualization. GB methodology. GB, RB, and MA formal analysis. GB and RB data curation. EE and MNK experimental validation. GB and RB visualization. GB, RB, EE, and MNK writing initial draft. All authors final draft. AM supervision.

Conflict of Interest Statement

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.

Acknowledgments

This work was financially supported by the Knut and Alice Wallenberg Foundation.

Supplementary Material

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

References

Agren, R., Liu, L., Shoaie, S., Vongsangnak, W., Nookaew, I., and Nielsen, J. (2013). The RAVEN toolbox and its use for generating a genome-scale metabolic model for Penicillium chrysogenum. PLoS Comput. Biol. 9:e1002980. doi: 10.1371/journal.pcbi.1002980

PubMed Abstract | CrossRef Full Text | Google Scholar

Agren, R., Mardinoglu, A., Asplund, A., Kampf, C., Uhlen, M., and Nielsen, J. (2014). Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling. Mol. Syst. Biol. 10:721. doi: 10.1002/msb.145122

PubMed Abstract | CrossRef Full Text | Google Scholar

Altekruse, S. F., Mcglynn, K. A., and Reichman, M. E. (2009). Hepatocellular carcinoma incidence, mortality, and survival trends in the United States from 1975 to 2005. J. Clin. Oncol. 27, 1485–1491. doi: 10.1200/JCO.2008.20.7753

PubMed Abstract | CrossRef Full Text | Google Scholar

Asgari, Y., Salehzadeh-Yazdi, A., Schreiber, F., and Masoudi-Nejad, A. (2013). Controllability in cancer metabolic networks according to drug targets as driver nodes. PLoS ONE 8:e79397. doi: 10.1371/journal.pone.0079397

PubMed Abstract | CrossRef Full Text | Google Scholar

Barabasi, A. L., Gulbahce, N., and Loscalzo, J. (2011). Network medicine: a network-based approach to human disease. Nat. Rev. Genet. 12, 56–68. doi: 10.1038/nrg2918

PubMed Abstract | CrossRef Full Text | Google Scholar

Benfeitas, R., Uhlen, M., Nielsen, J., and Mardinoglu, A. (2017). New challenges to study heterogeneity in cancer redox metabolism. Front. Cell Dev. Biol. 5:65. doi: 10.3389/fcell.2017.00065

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhutia, Y. D., Babu, E., Ramachandran, S., and Ganapathy, V. (2015). Amino Acid transporters in cancer and their relevance to “glutamine addiction”: novel targets for the design of a new class of anticancer drugs. Cancer Res. 75, 1782–1788. doi: 10.1158/0008-5472.CAN-14-3745

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., Shin, H. R., Bray, F., Forman, D., Mathers, C., and Parkin, D. M. (2010). Estimates of worldwide burden of cancer in 2008: GLOBOCAN 2008. Int. J. Cancer 127, 2893–2917 doi: 10.1002/ijc.25516

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., Soerjomataram, I., Dikshit, R., Eser, S., Mathers, C., Rebelo, M., et al. (2015). Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int. J. Cancer 136, E359–386. doi: 10.1002/ijc.29210

PubMed Abstract | CrossRef Full Text | Google Scholar

Folger, O., Jerby, L., Frezza, C., Gottlieb, E., Ruppin, E., and Shlomi, T. (2011). Predicting selective drug targets in cancer through metabolic networks. Mol. Syst. Biol. 7:501. doi: 10.1038/msb.2011.35

PubMed Abstract | CrossRef Full Text | Google Scholar

Friemel, J., Rechsteiner, M., Frick, L., Bohm, F., Struckmann, K., Egger, M., et al. (2015). Intratumor heterogeneity in hepatocellular carcinoma. Clin. Cancer Res. 21, 1951–1961 doi: 10.1158/1078-0432.CCR-14-0122

PubMed Abstract | CrossRef Full Text | Google Scholar

Fung, M. K. L., and Chan, G. C. (2017). Drug-induced amino acid deprivation as strategy for cancer therapy. J. Hematol. Oncol. 10:144. doi: 10.1186/s13045-017-0509-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Lin, L. P., Zhu, C. H., Chen, Y., Hou, Y. T., and Ding, J. (2006). Growth arrest induced by C75, A fatty acid synthase inhibitor, was partially modulated by p38 MAPK but not by p53 in human hepatocellular carcinoma. Cancer Biol. Ther. 5, 978–985. doi: 10.4161/cbt.5.8.2883

CrossRef Full Text | Google Scholar

Grossman, R. L., Heath, A. P., Ferretti, V., Varmus, H. E., Lowy, D. R., Kibbe, W. A., et al. (2016). Toward a shared vision for cancer genomic data. N. Engl. J. Med. 375, 1109–1112. doi: 10.1056/NEJMp1607591

PubMed Abstract | CrossRef Full Text | Google Scholar

Guney, E., Menche, J., Vidal, M., and Barabasi, A. L. (2016). Network-based in silico drug efficacy screening. Nat. Commun. 7:10331. doi: 10.1038/ncomms10331

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartung, F., Wang, Y., Aronow, B., and Weber, G. F. (2017). A core program of gene expression characterizes cancer metastases. Oncotarget 8, 102161–102175. doi: 10.18632/oncotarget.22240

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirsch, H. A., Iliopoulos, D., Joshi, A., Zhang, Y., Jaeger, S. A., Bulyk, M., et al. (2010). A transcriptional signature and common gene networks link cancer with lipid metabolism and diverse human diseases. Cancer Cell 17, 348–361. doi: 10.1016/j.ccr.2010.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Honeyman, J. N., Simon, E. P., Robine, N., Chiaroni-Clarke, R., Darcy, D. G., Lim, I. I., et al. (2014). Detection of a recurrent DNAJB1-PRKACA chimeric transcript in fibrolamellar hepatocellular carcinoma. Science 343, 1010–1014. doi: 10.1126/science.1249484

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeong, H., Mason, S. P., Barabasi, A. L., and Oltvai, Z. N. (2001). Lethality and centrality in protein networks. Nature 411, 41–42. doi: 10.1038/35075138

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonsson, P. F., and Bates, P. A. (2006). Global topological features of cancer proteins in the human interactome. Bioinformatics 22, 2291–2297. doi: 10.1093/bioinformatics/btl390

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanhaiya, K., Czeizler, E., Gratie, C., and Petre, I. (2017). Controlling directed protein interaction networks in cancer. Sci. Rep. 7:10327. doi: 10.1038/s41598-017-10491-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Khosraviani, M., Saheb Zamani, M., and Bidkhori, G. (2016). FogLight: an efficient matrix-based approach to construct metabolic pathways by search space reduction. Bioinformatics 32, 398–408. doi: 10.1093/bioinformatics/btv578

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhajda, F. P., Pizer, E. S., Li, J. N., Mani, N. S., Frehywot, G. L., and Townsend, C. A. (2000). Synthesis and antitumor activity of an inhibitor of fatty acid synthase. Proc. Natl. Acad. Sci. U.S.A. 97, 3450–3454. doi: 10.1073/pnas.97.7.3450

PubMed Abstract | CrossRef Full Text | Google Scholar

Law, V., Knox, C., Djoumbou, Y., Jewison, T., Guo, A. C., Liu, Y., et al. (2014). DrugBank 4.0: shedding new light on drug metabolism. Nucleic Acids Res. 42, D1091–D1097. doi: 10.1093/nar/gkt1068

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y. Y., Slotine, J. J., and Barabasi, A. L. (2011). Controllability of complex networks. Nature 473, 167–173. doi: 10.1038/nature10011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, W., Xu, Y., Guo, Y., Yu, Z., Feng, G., Liu, P., et al. (2016). The drug target genes show higher evolutionary conservation than non-target genes. Oncotarget 7, 4961–4971. doi: 10.18632/oncotarget.6755

PubMed Abstract | CrossRef Full Text | Google Scholar

Maddocks, O. D. K., Athineos, D., Cheung, E. C., Lee, P., Zhang, T., Van Den Broek, N. J. F., et al. (2017). Modulating the therapeutic response of tumours to dietary serine and glycine starvation. Nature 544, 372–376. doi: 10.1038/nature22056

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardinoglu, A., Agren, R., Kampf, C., Asplund, A., Uhlen, M., and Nielsen, J. (2014a). Genome-scale metabolic modelling of hepatocytes reveals serine deficiency in patients with non-alcoholic fatty liver disease. Nat. Commun. 5:3083. doi: 10.1038/ncomms4083

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardinoglu, A., Gatto, F., and Nielsen, J. (2013). Genome-scale modeling of human metabolism - a systems biology approach. Biotechnol. J. 8, 985–996. doi: 10.1002/biot.201200275

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardinoglu, A., and Nielsen, J. (2012). Systems medicine and metabolic modelling. J. Intern. Med. 271, 142–154. doi: 10.1111/j.1365-2796.2011.02493.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardinoglu, A., and Nielsen, J. (2015). New paradigms for metabolic modeling of human cells. Curr. Opin. Biotechnol. 34, 91–97. doi: 10.1016/j.copbio.2014.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardinoglu, A., Uhlen, M., and Nielsen, J. (2014b). The use of Genome-scale metabolic models for drug target and biomarker identification. N. Biotechnol. 31, S49–S49. doi: 10.1158/1078-0432.CCR-12-1856

CrossRef Full Text | Google Scholar

Moody, S. E., Schinzel, A. C., Singh, S., Izzo, F., Strickland, M. R., Luo, L., et al. (2015). PRKACA mediates resistance to HER2-targeted therapy in breast cancer cells and restores anti-apoptotic signaling. Oncogene 34, 2061–2071. doi: 10.1038/onc.2014.153

PubMed Abstract | CrossRef Full Text | Google Scholar

Najafi, A., Bidkhori, G., Bozorgmehr, J. H., Koch, I., and Masoudi-Nejad, A. (2014). Genome scale modeling in systems biology: algorithms and resources. Curr. Genomics 15, 130–159. doi: 10.2174/1389202915666140319002221

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, G. H., Xie, X., Xu, F., Shi, X., Wang, Y., and Deng, L. (2015). Distinctive pharmacological differences between liver cancer cell lines HepG2 and Hep3B. Cytotechnology 67, 1–12. doi: 10.1007/s10616-014-9761-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabbisetti, V. S., Chirugupati, S., Thomas, S., Vaidya, K. S., Reardon, D., Chiriva-Internati, M., et al. (2005). Calcitonin increases invasiveness of prostate cancer cells: role for cyclic AMP-dependent protein kinase A in calcitonin action. Int. J. Cancer 117, 551–560. doi: 10.1002/ijc.21158

PubMed Abstract | CrossRef Full Text | Google Scholar

Sapandowski, A., Stope, M., Evert, K., Evert, M., Zimmermann, U., Peter, D., et al. (2015). Cardiolipin composition correlates with prostate cancer cell proliferation. Mol. Cell. Biochem. 410, 175–185. doi: 10.1007/s11010-015-2549-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Schug, Z. T., and Gottlieb, E. (2009). Cardiolipin acts as a mitochondrial signalling platform to launch apoptosis. Biochim. Biophys. Acta 1788, 2022–2031. doi: 10.1016/j.bbamem.2009.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sontag, E. D. (2013). Mathematical Control Theory: Deterministic Finite Dimensional Systems. New York, NY: Springer Science & Business Media.

Google Scholar

Uhlen, M., Zhang, C., Lee, S., Sjostedt, E., Fagerberg, L., Bidkhori, G., et al. (2017). A pathology atlas of the human cancer transcriptome. Science 357:eaan2507. doi: 10.1126/science.aan2507

PubMed Abstract | CrossRef Full Text | Google Scholar

Varemo, L., Scheele, C., Broholm, C., Mardinoglu, A., Kampf, C., Asplund, A., et al. (2015). Proteome- and transcriptome-driven reconstruction of the human myocyte metabolic network and its use for identification of markers for diabetes. Cell Rep. 11, 921–933. doi: 10.1016/j.celrep.2015.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinayagam, A., Gibson, T. E., Lee, H. J., Yilmazel, B., Roesel, C., Hu, Y., et al. (2016). Controllability analysis of the directed human protein interaction network identifies disease genes and drug targets. Proc. Natl. Acad. Sci. U.S.A. 113, 4976–4981. doi: 10.1073/pnas.1603992113

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinstein, J. N., Collisson, E. A., Mills, G. B., Shaw, K. R. M., Ozenberger, B. A., Ellrott, K., et al. (2013). The Cancer genome atlas pan-cancer analysis project. Nat. Genet. 45, 1113–1120. doi: 10.1038/ng.2764

PubMed Abstract | CrossRef Full Text | Google Scholar

Wuchty, S. (2014). Controllability in protein interaction networks. Proc. Natl. Acad. Sci. U.S.A. 111, 7156–7160. doi: 10.1073/pnas.1311231111

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, F., Wang, C., Yin, H., Yu, J., Chen, S., Fang, J., et al. (2016). Leucine deprivation inhibits proliferation and induces apoptosis of human breast cancer cells via fatty acid synthase. Oncotarget 7, 63679–63689. doi: 10.18632/oncotarget.11626

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, X. L., Wei, M., Yunoki, T., Kakehashi, A., Yamano, S., Kato, M., et al. (2012). Long-term treatment with L-isoleucine or L-leucine in AIN-93G diet has promoting effects on rat bladder carcinogenesis. Food Chem. Toxicol. 50, 3934–3940. doi: 10.1016/j.fct.2012.07.063

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, H., Braun, P., Yildirim, M. A., Lemmens, I., Venkatesan, K., Sahalie, J., et al. (2008). High-quality binary protein interaction map of the yeast interactome network. Science 322, 104–110. doi: 10.1126/science.1158684

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, Z., Zhao, C., Di, Z., Wang, W. X., and Lai, Y. C. (2013). Exact controllability of complex networks. Nat. Commun. 4:2447. doi: 10.1038/ncomms3447

PubMed Abstract | CrossRef Full Text | Google Scholar