Cellular Genome-Scale Metabolic Modeling Identifies New Potential Drug Targets Against Hepatocellular Carcinoma.

Genome-scale metabolic modeling (GEM) is one of the key approaches to unpack cancer metabolism and for discovery of new drug targets. In this study, we report the Transcriptional Regulated Flux Balance Analysis-CORE (TRFBA-), an algorithm for GEM using key growth-correlated reactions using hepatocellular carcinoma (HCC), an important global health burden, as a case study. We generated a HepG2 cell-specific GEM by integrating this cell line transcriptomic data with a generic human metabolic model to forecast potential drug targets for HCC. A total of 108 essential genes for growth were predicted by the TRFBA-CORE. These genes were enriched for metabolic pathways involved in cholesterol, sterol, and steroid biosynthesis. Furthermore, we silenced a predicted essential gene, 11-beta dehydrogenase hydroxysteroid type 2 (HSD11B2), in HepG2 cells resulting in a reduction in cell viability. To further identify novel potential drug targets in HCC, we examined the effect of nine drugs targeting the essential genes, and observed that most drugs inhibited the growth of HepG2 cells. Some of these drugs in this model performed better than Sorafenib, the first-line therapeutic against HCC. A HepG2 cell-specific GEM highlights sterol metabolism to be essential for cell growth. HSD11B2 downregulation results in lower cell growth. Most of the compounds, selected by drug repurposing approach, show a significant inhibitory effect on cell growth in a wide range of concentrations. These findings offer new molecular leads for drug discovery for hepatic cancer while also illustrating the importance of GEM and drug repurposing in cancer therapeutics innovation.


a from this c
ll type and a generic human metabolic model.We next assessed the strength of the generated model in predicting cancer essential genes and further examined the knockdown effect of one of the predicted essential genes on the reduction of HepG2 cell growth.Moreover, taking a drug repurposing approach, we selected several approved drugs targeting the predicted essential genes identi ed by our algorithm, and examined their treatment effect on HepG2 cell viability in vitro.We speci cally compared the inhibitory effect of these compounds in different concentrations to Sorafenib, a well-known approved drug for advanced HCC.


Results And Discussion


Metabolic modeling of HCC

Here, we employed TRFBA-CORE, a benchmark-driven algorithm developed by an extensive evaluation of other cancer genome-scale metabolic modeling methods, to nd potential drug targets for HCC 27 .We also reconstructed two cell-speci c GEMs using mCADRE and CORDA algorithms to compare their predictive capabilities with TRFBA-CORE in simulating a set of de ned metabolic tasks for the liver 12,19,26 .CORDA was able to predict a wider range of liver metabolic tasks (Figure 1).However, the enrichment of cell-speci c models with essential genes for HepG2 showed that only GEM generated by TRFBA-CORE was able to signi cantly predict essential genes (Enrichment p values = 0.005, 0.098, and 0.086, for TRFBA-CORE, CORDA, and mCADRE, respectively).

The use of growt -correlated key reactions and the FASTCORE minimal modeling approach may explain the ability of TRFBA-CORE to signi cantly predict HepG2 essential genes.The size of generated models further shows that the TRFBA-CORE model has the minimum number of reactions required to render the network ux consistent (Figure 2).


Identi cation of potential drug targets for HCC

Of 139 essential genes predicted by TRFBA-CORE, 31 genes (22%) disrupted de ned metabolic functions in energy metabolism and redox equilibrium, or biomass production, predominantly in the OXPHOS pathway (Supplementary table S1).Therefore, we excluded these genes due to their potential adverse effects on healthy cells, and performed a gene-set enrichment analysis using the remaining 108 predicted essential genes to better understand their functional impacts.Interestingly gene ontology (GO) biological processes showed an enrichment of processes involved in cholesterol, sterols and steroids biosynthesis (Supplementary table S2).Furthermore, KEGG (Supplementary table S3) and Reactome (Supplementary table S4) biochemical pathways involved in lipid and lipoprotein metabolism, biosynthesis, cholesterol and steroids regulation, and carbohydrate and amino acids metabolisms had a sig

edicted essential genes.
holesterol is a key molecule in cell membrane structure, and a prerequisite for essential hormones and bile acids.The liver is one of the main organs for the synthesis and metabolism of endogenous cholesterol in the body 28,29 .The growth of HCC tumors is dependent on cholesterol biosynthesis, and the use of statins (inhibitors of HMG-CoA reductase) has been suggested as a candidate therapy for HCC.

Although this enzyme is a key regulatory enzyme in cholesterol synthesis, several other enzymes in the cholesterol biosynthesis pathway have been targeted to prevent or reduce cell proliferation in various cancers 28 .Furthermore, previous ndings suggest that the mevalonate pathway plays a role in the development of HCC [30][31][32] .Of the 18 essential genes predicted for cholesterol metabolism, 17 metabolic genes (excluding sterol O-acyltransferase 1 (SOAT1)) has been previously identi ed to be involved in the production of potential antimetabolites reducing the growth of HCC tumors 28 .Moreover, glucocorticoids, such as dexamethasone, are a class of steroids widely use

as antiin ammatory agents.Ma et al. showed that
he expression of two essential enzymes that regulate endogenous glucocorticoids and glucogenesis was altered in malignant hepatocytes, so that the expression of 11-beta dehydrogenase hydroxysteroid type 1 (HSD11B1) and 2 (HSD11B2) in mouse and human tumor cells were up and down-regulated, respectively.Moreover, the expression ratio of HSD11B1:HSD11B2 was related to the prognosis and survival of patients with HCC, so that patients with a high ratio of HSD11B1: HSD11B2 had a longer survival time than others 33,34 .Also, it has been shown that HSD11B1 expression usually leads to a decrease in cell proliferation, while HSD11B2 expression is involved in an increase in cell proliferation.Several studies showed that the expression of these two genes varies in different tumors and can provide a suitable microenvironment for tumor growth 35,36 .Therefore, HSD11B2, which was the only predicted essential gene in the steroid biosynthesis pathway, was selected for further investigation.

We compared the relative expression level of HSD11B in HepG2 and three non-cancerous cell lines, namely HH, HHSEC and LX-2 (Figure 3a) and observed a relatively higher expression of HSD11B2 in HepG2 cells.Therefore, we silenced HSD11B2 in HepG2 cells (Figure 3b) to examine its knockdown effect on the cell viability.Interestingly, we were able to observe a modest but signi cant reduction in cell viability after 48 (p-value = 0.021) and 72 (p-value = 6.5E-4) hr (Figure 4) Glucocorticoids are essential hormones in the body, regulating the nuclear expression of approximately 5% of human genes.One of their function is to increase glucose synthesis by increasing gluconeogenesis.Ma et al. showed that glycogenesis is impaired in HCC, and that the abnormal expression of two genes, HSD11B1 and HSD11B2, affects some key genes in the glucose pathway.Although the mechanism by which malignant hepatocytes target these two genes for malignant progression is not clear, targeting HSD11B2 in a mouse model reduced tumor mass 34 .Furthermore, by comparing the RNA-Seq expression data in GEPIA web server, we observed that HSD11B1 and HSD11B2 were respectively down-and up-regulated in HCC samples (P <0.05).The expression of these two genes showed a negative correlation in healthy and cancerous samples from Atlas of Cancer Genome (TCGA, Spearman R=-0.18, p=3×10 -4 ), while, the correlation in healthy samples was positive (Spearman R = 0.36, p = 0.011).Therefore, given the lower expression of HSD11B1 in HCC cells, HSD11B2 may be a potential target gene in HCC to regulate the HSD11B1: HSD11B2 expression ratio 37,38 .


Drug repositioning for HCC

Drug repurposing is an attractive approach providing cancer patients with faster and cost-effective potential therapeutics with well-characterized safety pro les 39,40 .Here, we used the DrugBank database to evaluate potential compounds targeting the predicted metabolic essential genes 41 .Thus, a set of metabolic drugs targeting one or more of the predicted genes was obtained, among which we selected 9 compounds by reviewing the literature and applying the following criteria; rst, metabolic drug has not been used in the treatment of HCC, and secondly, has been evaluated in at least one of the other cancers, or the target gene was suggested as a therapeutic target for HCC (Table 1).We next examined the treatment effect of these compounds on HepG2 cell viability in 5 different concentrations.Interestingly, most of the drugs, except Terbina ne, were able to signi cantly reduce the cell viability (Figure 5).In addition, some of these drugs showed a comparatively higher inhibitory effect compared to Sorafenib, the well-known approved drug to treat advanced HCC.For instance, Clofarabine, Alanosine, and Pralatrexate were more effective than Sorafenib at lower concentrations (100 and 1 nm).Alanosine targeted two genes adenylosuccinate synthetase like 1 (ADSSL1) and adenylosuccinate synthetase (ADSS), which TRFBA-CORE predicted as synthetic lethal.These two genes encode two isozymes, adenylosuccinate synthase 1 and 2, respectively, catalyzing a key reaction in purine nucleotide biosynthesis.The effect of Alanosine on inhibiting the growth of other cancerous cells has also been demonstrated 42 .Targeting AMP, another critical precursor by Alanosine, has also been shown to be effective at the nucleotide levels 43 .

AMP is produced from the I

by ADSS and ADSSL1, and an
essential precursor in the production of DNA and RNA 44 .Synthetic lethality occurs in a cell when perturbations of two genes combined are lethal, while a genetic defect in either gene alone is harmless.In this case, a genetic change, such as a defect in one tumor suppressor gene causes another gene to become essential for cell viability, and therefore, offering a therapeutic approach by selective targeting of cancer cells 45 .TRFBA-CORE predicted a total of 25 synthetic lethal pairs, nine of which had previously been identi ed in the study of Folger et al. 46 .Interestingly, none of these predicted gene pairs disrupted the de ned metabolic functions in Recon 1.

Among the novel predicted synthetic lethal pairs, hydroxyl methylglutaryl-CoA synthase (HMGCS2), which is involved in cholesterol metabolism, contained eight gene pairs all of which were involved in amino acid metabolism (Supplementary table S5).Moreover, carnitine O-acyl transferase (CRAT), which facilitates the release of acetyl-CoA mitochondria into the cytosol 47 , had three gene pairs that were involved in TCA and citrate transport from mitochondria to the cytosol.The role of HMGCS2 and CRAT in several cancers has been already shown 47,48 .

In addition to Alanosine, Clofarabine, one of the most successful drugs against leukemia 49 , targets the large subunit of ribonucleoside-diphosphate reductase (RRM1).Notably, few other drugs targeting RRM1

have previously been suggested to inhibit HCC cell growth and other liver diseases 50,51 .RRM1 and RRM2 are responsible for ribonucleotide reductase expression, which acts as a rate-limiting enzyme in dNTP production.Also, Pralatrexate targets dihydrofolate reductase (DHFR) and thymidylate synthetase (TYMS), bot of which were predicted to be essential genes.Pralatrexate is an approved drug in the treatment of peripheral T cell lymphoma, and its promising outcome has been shown in some other cancer cells 52,53 .However, the effect of Pralatrexate on growth was relatively constant at different concentrations (Figure 5).TYMS plays a vital role in dTMP synthesis, which is signi cantly upregulated in tumor cells 54 .Tetrahydrofolate (THF), the biologica ly active form of folic acid, is an essential metabolite synthetized by DHFR and utilized by TYMS to produce dTMP, a critical compound in DNA replication 55,56 .

Interestingly, Tri uridine had a larger inhib tory effect on cell growth than Sorafenib at a low concentration of 1 µm.Tri uridine/ tipiracil combination therapy has been shown to have a bene cial effect on survival in those with colon cancer 57 .Furthermore, several other drugs targeting TYMS have been suggested for the treatment of HCC 58,59 .

Notably, while Sorafenib showed a relatively stronger inhibitory effect on cell viability at 10 µm, thioconazole, an imidazole antifungal drug, showed a larger effect compared to the other compounds at higher concentrations (25 µm and 50 µm, Figure 5).One of the drug targets of thioconazole is lenosterol 14-alpha-demethylase (CYP51A1), which has been linked to the size of HCC tumors in the carriers of the hepatitis C virus 60 .Recently, the effect of thioconazole on inhibiting the growth of cancer cells has been shown 61 .Moreover, at these concentrations, the inhibitory effect of Alendronic acid, Alanosine, Ibandronat and Atovaquone, on growth was similar to that of Sorafenib.Alendronic acid and Ibandronate both inhibit farnesyl pyrophosphate synthase (FDPS).Pamidronic acid and zoldronic acid, both of which are FDPS inhibitors, have also shown antiproliferative effects in HCC 62,63 .Moreover, the e cacy of Atovaquone, an antimalarial drug, has recently been dem nstrated in various cancers 64,65 .Atovaquone inhibits the activity of dihydroorotate dehydrogenase (DHODH), which is involved in de novo pyrimidine synthesis, resulting in a forced pause of DNA and RNA synthesis in cancer cells 66-68 .


Materials And Methods


General human model set-up

The consistent Recon 1 general model (i.e. after removal of all blocked reactions incapable of carrying a non-zero ux) w

t-speci c models 13,23,46 .
t has been shown that alanine and glutamate are secreted (but not uptaken) by cancer cells; therefore, the corresponding exchange reactions were not constrained in the input model 69,70 .TRFBA-CORE algorithm was used to generate a functional cell-speci c model from HepG2 expression data for HepG2 cells (taken from

e Gene Expression Omnibus
GEO) public repository under the accession number GSE38124) 71,72 .The generated cellspeci c GEM was then constrained with simulated MEM + 10% FBS medium 46,73-75 to be consistent with in-vitro experiments.


Liver metabolic functions

A set of 256 previously de ned metabolic functions (such as the synthesis of fatty acids, amino acids, cholesterol, and bile acids), which are physiologically relevant to healthy hepatocytes, were used to validate the generated cell-speci c GEM 9,12 .We kept only a set of 197 tasks that could be performed by consistent Recon 1, and GEM performance was reported as the fraction of successfully passed metabolic functions.To further evaluate the performance of the generated cell-speci c model, two additional cellspeci c GEMs were generated using CORDA and mCADRE algorithms.CORDA 26 uses a dependency assessment to minimize the ux of cost-consuming reactions, and only the non-key reactions required for the nal network activity are considered.Similarly, mCADRE reconstructs context-speci c modes by using gene expression and network topology to score all reactions in a model 19,27 .


Predict

n of essential genes

Gene dep
ndence scores calculated by the DEMETER2 algorithm for the HepG2 cell line were obtained from the DepMap database 76 .The DEMETER2 algorithm analyzes the RNA interference (RNAi) screens, reduces the off-target effects of RNAi and corrects for confounding factors such as batch effects 76,77 .

Here, the genes with a DEMETER2 dependence score below zero were considered as being essential.

Finally, a hypergeometric enri hment test was used to evaluate the over-representation of essential genes present in each of the generated GEMs.

To account for possible side-effects of the predicted essential genes, we used a set of 56 metabolic tasks occurring in all cell types 28 ;therefore, we excluded genes that disrupted any of the tasks involved in energy metabolism and redox balance, or biomass production 78 .Double ge e deletion simulation was also performed to identify synthetic lethal genes.A synergy score was de ned as KO AB / min (KO A , KO B );

where KO A , KO B , and KO AB are simulated growth after the repression of gene A, gene B, and simultaneous suppression of genes A and B, respectively 46 .Fin

ly, double ge
es with a synergy score greater than 0.5

were selected as synthetic lethal 46 .Similar to single gene deletion simulations, double gene deletions that perturbed any of the above-mentioned metabolic functions were excluded.

Predicted essential genes were further used to perform gene ontology (GO) biological processes, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome biochemical pathways enrichment analysis using Enrichr.Adjusted P values and combined scores (log P value × z-score) were reported 79 .


Cell culture

HepG2 (hepatocellula

carcinoma) derived fr
m the liver tissue of a 15-year-old Caucasian male who had a well-differentiated hepatocellular carcinoma were purchased from ATCC (ATCC® HB8065™).After thawing, the cells were plated in T-75 asks and cultured at 37 °C in Minimum Essential Medium (MEM), containing Fetal Bovine Serum (FBS) 10%, L-glutamine 2 mM, sodium pyruvate 1 mM and non-essential amino acids 1X (NEAA) in a humidi ed atmosphere of 5% CO2.


Gene expression

Total RNA of HepG2 cells was extracted by the RNeasy Mini kit (Qiagen, Valencia, CA) and rst-strand cDNA was synthesized using the High Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scienti c) according to the manufacturer's instructions.HSD11B2 expression in HepG2, human hepatocytes (HH), Human Hepatic Sinusoidal Endot elial Cells (HHSEC), and Human Hepatic Stellate Cells (LX-2) (Sciencell, Carlsbad, CA) was measured using TaqMan probes and Master Mix (Thermo Fisher Scienti c) using quantitative real-time PCR (qPCR) according to the manufacturer's protocol.Expression levels were determined using the 2 -ΔΔCT method and normalized to β-actin in triplicates 80,81 .


Gene silencing

HepG2 were transfected with small interfering RNA (siRNA) or negative control (scrambled) siRNA (Ambion-Lif

Technologies
by Thermo Fisher Scienti c, Rockford, IL) with Lipofectamine Transfection Reagent (Life Technologies, Carlsbad, CA) according to the manufacturer's instructions.The cells were rst cultured in 6-well plates (3×105 cells/well) for one day, and transfected with 50 nM siRNA oligos in a serum-free medium at 50-70% con uency.After 48 h, the RNA was extracted by the RNeasy Mini Kit (Qiagen, Hilden, Germany), and after cDNA synthesis, t HSD11B2 knockdown e ciency was

n rmed.All
xperiments were performed in triplicates.


Cell viability assay

To observe the effect of gene knockdown on cell growth, cell viability was investigated by CellTiter 96® AQueous One Solution Assay (Promega, Madison, WI, USA).This assay is based on the ability of living cells to reduce the tetrazolium, 3-(4, 5-dimethylthiazol-2-yl) -5-(3-carboxymethoxyphenyl) -2-(4sulfophenyl) -2H-tetrazolium (MTS), to formazan.For this purpose, at 24, 48, nd 72 hours after knockdown, the culture medium was changed, and 20% of the volume of the culture medium containing cells, MTS solution was added to the wells containing transfected cells.After incubation for 1 hour, the absorbance of the samples was measured at 490 nm in a microplate reader.All experiments were performed at least three times, with ve repetitions each time.
o evaluate the treatment effect of candidate drugs on cell viability, HepG2 cells were cultured at a density of 5000 cells/well in a 96-well plate, and after incubating for overnight, the cells were treated with each drug at 5 different concentrations (100 nm, 1 µm, 10 µm, 25 µm and 50 µm in 0.5% DMSO).Cells treated with 0.5% DMSO were used as contro .All experiments were performed in triplicates.After 72 and 96 hours, cell viability was measured as described above.


4.

This study has at least two limitations: a) the use of a single cell type, namely HepG2, and results may be driven by the speci c genetic background

f these cells
deriving from a single donor.For example, they are homozygotes for the PNPLA3 I148M genetic variant, the strongest common inborn genetic risk factor predispos ng to hepatocellular carcinoma 3,82,83 ; b) the e cacy of the compounds should be con rmed in other cell types and in murine models of HCC.


Conclusion

we used our recently developed algorithm, TRFBA-CORE, and by integrating HepG2 expression data and a general human model, generated a speci c cell model for HepG2.Many of the reactions present in this minimal HepG2 model overlapped with model content of the two models generated by mCADRE and CORDA.However, the use of growth cor elation reactions by TRFBA-CORE improved the capability of the model in predicting essential HepG2 genes.

We further con rmed that the knocking down of one predicted essential gene, HSD11B2, signi cantly reduced HepG2 cell growth.Previous ndings suggest that the expression of this gene is directly related to the growth o