A hierarchical approach employing metabolic and gene expression profiles to identify the pathways that confer cytotoxicity in HepG2 cells
BMC Systems Biology volume 1, Article number: 21 (2007)
Free fatty acids (FFA) and tumor necrosis factor alpha (TNF-α) have been implicated in the pathogenesis of many obesity-related metabolic disorders. When human hepatoblastoma cells (HepG2) were exposed to different types of FFA and TNF-α, saturated fatty acid was found to be cytotoxic and its toxicity was exacerbated by TNF-α. In order to identify the processes associated with the toxicity of saturated FFA and TNF-α, the metabolic and gene expression profiles were measured to characterize the cellular states. A computational model was developed to integrate these disparate data to reveal the underlying pathways and mechanisms involved in saturated fatty acid toxicity.
A hierarchical framework consisting of three stages was developed to identify the processes and genes that regulate the toxicity. First, discriminant analysis identified that fatty acid oxidation and intracellular triglyceride accumulation were the most relevant in differentiating the cytotoxic phenotype. Second, gene set enrichment analysis (GSEA) was applied to the cDNA microarray data to identify the transcriptionally altered pathways and processes. Finally, the genes and gene sets that regulate the metabolic responses identified in step 1 were identified by integrating the expression of the enriched gene sets and the metabolic profiles with a multi-block partial least squares (MBPLS) regression model.
The hierarchical approach suggested potential mechanisms involved in mediating the cytotoxic and cytoprotective pathways, as well as identified novel targets, such as NADH dehydrogenases, aldehyde dehydrogenases 1A1 (ALDH1A1) and endothelial membrane protein 3 (EMP3) as modulator of the toxic phenotypes. These predictions, as well as, some specific targets that were suggested by the analysis were experimentally validated.
Elevated levels of free fatty acids (FFAs) have been implicated in the pathogenesis of many obesity-related metabolic disorders [1–4], such as fatty liver disease and steatohepatitis. Dietary fatty acids produce a variety of metabolic and genetic effects on liver cells. Fatty acids compete with glucose for oxidation at the TCA cycle . Fatty acids also cause changes in the enzyme make-up of the cells by regulating the transcription of enzymes of metabolism. FFAs exert their transcriptional effects by activating transcription factors (TFs) such as sterol receptor element binding protein (SREBP), peroxisome proliferator activated receptors (PPARs) and hepatic nuclear factors (HNFs) . PPARs regulate the expression of proteins involved in fatty acid oxidation, SREBP regulates phospholipid and cholesterol synthesis and HNF affects both lipid and carbohydrate metabolism.
Tumor necrosis factor alpha (TNF-α) is another factor that has been shown to affect the function of hepatocytes in numerous ways. It has been associated with the development of hepatic insulin resistance and hepatocyte cell death [7–10]. TNF-α also activates transcription factors such as nuclear factor kappa B (NFκB) and c-Jun [11, 12]. These transcription factors alter the expression of genes involved in cellular metabolism, cell proliferation and cell death [11, 12]. Hepatocytes are known to be resistant to the cytotoxic action of TNF-α due to prompt upregulation of cytoprotective genes mediated by the activation of NF-κB in response to TNF-α . Therefore, the cytotoxic effect of TNF-α requires a secondary insult, e.g., transcriptional inhibition  or glutathione depletion .
In vivo, under conditions of obesity, hepatocytes are simultaneously exposed to elevated FFAs and TNF-α. The importance of these factors in the pathogenesis of many diseases motivated this study of the physiological, metabolic and genetic effects of the simultaneous exposures to different types of FFAs and TNF-α. Among the many responses to FFA and TNF-α, the mechanism of cell death in response to simultaneous exposure to these factors is not well characterized. Hepatocyte cell death is suggested to play an important role in the development of various hepatic disorders, e.g. in non-alcoholic steatohepatitis (NASH). Previous studies on the toxic effects of different types of FFAs on liver cells have identified that saturated FFAs are much more toxic than unsaturated FFAs [16–19]. These studies have suggested that ER-ROS (reactive oxygen species) stress , mitochondrial alterations  and lysosomal permeabilization  are the major mechanisms in the toxicity of saturated FFAs. Studies on the toxicity of saturated FFAs to other cells have suggested that increased ROS production and ceramide synthesis [20, 21] are the major mechanisms of palmitate-toxicity in those cells. Similarly, increased ceramide and ROS generation have been suggested to play important roles in the toxicity of TNF-α to hepatocytes [10, 22]. However, there has not been any study on the cytotoxic effects of simultaneous exposure to FFAs and TNF-α.
Human hepatoblastoma cells (HepG2) were treated simultaneously with different types of FFAs (saturated, monounsaturated and polyunsaturated) and levels of TNF-α and the cytotoxic, metabolic and genetic responses were studied. Exposing the cells to saturated fatty acid (palmitate) was cytotoxic and the exposure to TNF-α increased this toxicity, whereas the unsaturated FFAs induced increased triglyceride (TG) accumulation and were not cytotoxic. TNF-α alone was not toxic, either alone or in combination with the unsaturated FFAs. Our objective was to develop an approach to elucidate the underlying pathways that confer the cytotoxicity. The hierarchical framework developed to integrate the metabolic and genetic information to identify the genes and biological processes regulating a phenotypic response is shown in Figure 1. The framework consisted of three stages. First, the metabolic changes associated with the cytotoxic phenotype were identified with Fisher's Discriminant Analysis (FDA) . Ketone body (e.g., beta-hydroxybutyrate (BOH) and acetoacetate (AcAc)) release and TG accumulation were identified to be the most important in separating the phenotypes, suggesting their involvement in inducing the cytotoxicity. To identify the signaling and gene regulatory pathways involved in regulating the toxicity, the genomic responses were measured using cDNA microarrays and analyzed with GSEA . GSEA has been applied successfully to identify the signaling, transcriptional regulatory and metabolic pathways involved in the development of type 2 diabetes . GSEA indicated the involvement of mitochondria-related, oxidative stress-related and FFA metabolism pathways in inducing the cytotoxic phenotype. Finally, the gene expression and metabolic profiles were integrated with multi-block partial least squares (MBPLS) regression analysis  to identify the genes most relevant to the metabolic changes which correlated most with the cytotoxicity. Some of the identified genes were experimentally perturbed to validate their predicted roles in regulating the cytotoxicity. This analysis identified that NADH dehydrogenases, aldehyde dehydrogenase 1A1 (ALDH1A1) and endothelial membrane protein 3 (EMP3) played important roles in the toxicity of FFAs.
1. Metabolic changes relevant to cytotoxicity
As shown in Figure 2, the mono- and poly- unsaturated fatty acids, TNF-α and their combinations were found to be non-toxic. However, the saturated fatty acid palmitate was toxic to the cells and TNF-α increased this toxicity. Because the release of LDH is a late event in the process of cell death, we also measured the activation of caspase-3. Exposure to palmitate caused a significant elevation in the caspase-3 activity (Figure 3). Exposure to TNF-α significantly increased the caspase-3 activation only in the palmitate-treated cells, but not in the cells treated with control medium or oleate. These results indicated significant interactive effects of palmitate and TNF-α on the process of cell death. Finally, we also found that cells exposed to palmitate were also TUNEL positive (Figure 4), further corroborating with the evidence of initiation of apoptosis in these cells. Because palmitate was found to be the primary factor in the toxicity and the effects of TNF-α were only secondary to that of palmitate, we further investigated the mechanisms of toxicity of this FFA. To evaluate the metabolic changes responsible for mediating the toxicity, the rates of uptake/release of 27 metabolites, of glucose, fatty acid and amino acid metabolism, were measured and shown in additional data file 1. FDA was applied to identify the metabolic changes responsible for separating the phenotypic response (cytotoxic versus nontoxic as defined by the level of lactate dehydrogenase (LDH) release). As shown in Figure 5, the palmitate-treated samples were separated from the rest in a space defined by the first FDA projection (T1) with the palmitate samples taking on negative values while the rest of the samples taking on positive values on the Y axis. Furthermore, as shown in Figure 6 the production of BOH and acetoacetate, with the most negative coefficients, separated the palmitate samples in the direction of toxicity, while the uptake of cysteine, ornithine and phenylalanine, with the most positive coefficients, separated the samples away from the toxicity. In addition, the palmitate-treated samples were also separated from the rest of the sample in the space defined by the second FDA projection (T2) with the palmitate samples having smaller values than the rest of the samples as shown in Figure 7. Correspondingly, we found that TG and ornithine with the most positive coefficients in this projection separated the palmitate samples away from the toxicity as shown in Figure 8. Acetoacetate and BOH are the products of fatty acid oxidation. Increased fatty acid oxidation is known to increase oxidative stress and cell death [27, 28]. TG accumulation has been found to protect cells from palmitate induced lipoxicity . Thus, FDA analysis confirmed that fatty acid oxidation is positively associated with the palmitate-induced cytotoxicity, while accumulation of TG is negatively associated to cytotoxicity.
2. Functional pathway analysis with GSEA
37 gene sets involved in a variety of cellular processes and organelles which are known to play important roles in the changes caused by FFA and TNF-α were evaluated to identify which of these processes were associated with palmitate-induced cytotoxicity. The gene sets evaluated included those of fatty acid metabolism, cell death, TNF-α signaling, etc., (see Table 1 for a full list). 14 of the 37 gene sets were significantly enriched with nominal p values less than 0.05, and shown in Table 2. The enriched pathways included electron transport chain (ETC), fatty acid metabolism, glycolysis, oxidative phosphorylation, ROS, the pentose phosphate pathway (PPP), cell death, fatty acid beta-oxidation, TCA cycle, fructose metabolism, glutathione, and the ERK 1, ERK 2, MAP kinase pathways. This suggested that the toxicity may be associated with changes in oxidative stress related pathways such as ROS and glutathione, as well as energy generating processes of glycolysis, oxidative phosphorylation and ETC. Importantly, the gene set of sphingolipid (ceramide) metabolism was not selected, indicating that ceramide metabolism may not play an important role in the observed toxicity. Notably, most of the enriched gene sets belong to the mitochondria, for example fatty acid beta-oxidation, electron transport chain and oxidative phosphorylation, depicting a central role of mitochondria in the cytotoxicity. The gene set of mitochondria, which included 229 genes, was also significantly altered in the palmitate-treated cells, with a nominal p value of 0.00189.
3. Integrating the metabolic and the gene expression profiles to identify the pathways relevant to the cytotoxicity
To identify the genes that regulate the metabolic functions most closely associated with cytotoxicity, MBPLS models were developed to predict the metabolic processes based upon the expression data of the gene sets identified by GSEA. Ketogenesis and TG accumulation were identified to be positively and negatively related to LDH release in step 2, respectively. Therefore, two MBPLS models were developed to model BOH and TG, respectively. The MBPLS models contained 14 blocks each, corresponding to the 14 enriched gene sets identified by the GSEA. The importance of individual genes within these functional groups was identified by evaluating the regression coefficients of the genes. In particular, the genes with high positive regression coefficients to ketogenesis and TG accumulation are discussed below as the predicted roles of these genes were evaluated using inhibitors or RNA interference (RNAi).
Increased production of ketone bodies, such as BOH and acetoacetate, has been shown to be associated with oxidative stress [27, 28] and an important role of oxidative stress has been identified in the palmitate toxicity . Therefore, genes with large (positive or negative) regression coefficients to ketone body release could suggest pathways and genes relevant to the cytotoxicity. The genes with the largest positive regression coefficients to B-OH are listed in Table 3. The complete list of regression coefficients can be found in the additional data file 2. Genes relevant to ROS production, fatty acid metabolism, and detoxification of lipid peroxidation products were found to have high regression coefficients. GSTM5 was most positively related to the BOH (and in turn, the cytotoxicity), while other isoforms of GST, e.g., GSTM4 and GSTM1 (listed in additional file 2), had negative regression coefficients. The basal level of GSTM5 in the liver is much lower than those of GSTM1 and GSTM4 . Therefore, in spite of a large positive coefficient of GSTM5, the combined negative effect of GSTM1 and GSTM4 on the cytotoxicity is expected to dominate. Similarly, aldehyde dehydrogenase 1 (ALDH1A1) was found to have a positive relation to cytotoxicity, while ALDH1A3 (listed in additional file 2) had a negative regression coefficient. We found that the basal levels of ALDH1A3 were about 10-5-times lower than that of ALDH1A1 (Table 4). Therefore, the cytotoxic effects of ALDH1A1 would be greater than the cytoprotective effects of ALDH1A3. Indeed, over-expression of ALDH1A3 had no significant effect on the cytotoxicity or the caspase-3 activation by palmitate (not shown). Two NADH dehydrogenases had highly positive regression coefficients, indicating a positive role of these genes in the cytotoxicity. Their combined coefficient was greater than any other gene, indicating that these genes may play very important roles in the toxicity. NADH dehydrogenase complex, or complex I of the electron transport chain, is a main site of superoxide production [30, 31]. In the course of electron transfer some of the activated oxygen is released as superoxides or H2O2 . Their selection suggested involvement of ROS in the toxicity. Indeed, elevated ROS production was observed in cells treated with palmitate (see section 4 b ii).
Intracellular TG was identified to have significantly negative association to the toxicity. In agreement with this, the channeling of palmitate to TG has been shown to protect the cells from its toxic effects . Thus, genes with high positive regression coefficients to intracellular TG accumulation may be potential targets for cytoprotection. Such genes are listed in Table 5. A complete list of regression coefficients is available in the additional data file 3. Endothelial membrane protein 3 (EMP3) was the gene with greatest positive coefficient, indicating that a reduction in the expression of this gene would be cytotoxic and vice-versa. The function of EMP3 is not clearly understood, it is known to be involved in cell growth, differentiation, and apoptosis [32, 33]. Many genes involved in glutathione metabolism had positive regression coefficients. Examples of genes related to glutathione that were selected included glutamate-cysteine ligase, catalytic subunit (GCL-c), glutathione S-transferase M4 (GSTM4), microsomal glutathione S-transferase 1 (MGST1), glutathione S-transferase A2 (GSTA2) and glutathione reductase (GSR). The selection of these genes indicated a cytoprotective role of these genes, and corroborated with the role of oxidative stress in the toxicity. Finally, some genes involved in regulation of glycolysis such as hexokinase 1 (HK1) and glucokinase (hexokinase 4) regulatory protein (GCKR) had large positive regression coefficients, indicating that upregulation of glycolysis may be cytoprotective. Exposure to FFAs is associated with a reduction in cellular energetics, e.g., a decrease in glycolytic enzymes [34, 35] and a reduction in ATP synthesis due to mitochondrial uncoupling . This reduction in cellular ATP levels may play a role in the toxicity  or exacerbate it. Under these conditions, increasing glycolysis may provide an alternative route for ATP synthesis and reduce the toxicity.
4. Experimental validations
Hepatocytes are resistant to the cytotoxic effects of TNF-α due to rapid upregulation of cytoprotective genes, mediated in part by the activation of NF-κB by TNF-α. It is for this reason that most of the previous studies on the cell death caused by TNF-α employ a secondary insult, such as transcriptional inhibition or glutathione depletion [4, 15]. The observation of the dependence of TNF-α toxicity on that of palmitate suggests that the toxicity of palmitate can act as a secondary/additional insult. Because the saturated free fatty acid (palmitate) was found to have the greatest toxicity of all the treatments and the toxicity of TNF-α depended on the effect of palmitate, the subsequent validations were conducted for the palmitate conditions.
a. Validation of GSEA results
a.i. The role of mitochondria
Mitochondria are known to be central to cell death pathways, and most cytotoxic insults are associated with a loss of mitochondrial integrity. Mitochondrial membrane potential is an indicator of mitochondrial integrity and functioning. GSEA suggested involvement of these organelles in the cytotoxicity. We measured the mitochondrial potentials in response to the various treatments, as shown in Figure 9. Palmitate-treatment significantly decreased the mitochondrial potential. The loss of mitochondrial integrity will cause many deleterious effects, including the breakdown of energy generation, increased ROS production, and the release of apoptogenic factors such as cytochrome c from the mitochondria .
a.ii. Lack of involvement of de novo ceramide synthesis
De novo ceramide synthesis has been suggested to play an important role in the toxicity of saturated FFAs to beta cells and cardiomyocytes. However, the ceramide metabolism pathway was not found to be significantly enriched according to GSEA, suggesting that this pathway may not play an important role in the toxicity of the saturated FFAs to hepatoma cells. To test this prediction, we treated the cells with palmitate in the presence of an inhibitor of de novo ceramide synthesis, Fumonisin B1 (FB1). Treatment with this inhibitor reduced the ceramide levels significantly (by about 77%, Table 6). However, no change in the toxicity was observed. Similarly, treatments with higher concentrations of this inhibitor (up to 100 uM) had no effect on the toxicity of the saturated FFA. Therefore, these experiments corroborated the GSEA prediction that ceramide metabolism does not play an important role in the toxicity. These results were also in agreement with a previous study which found that hepatocyte cell death induced by the saturated FFAs is independent of ceramide synthesis [17, 18].
b. Validating hierarchical model predictions
b.i. Important role of NADH dehydrogenase
Exposure to palmitate has been shown to increase the ROS production in CHO cells . We investigated the role of ROS production in the palmitate-toxicity in our case. As shown in Figure 10, the ROS levels were significantly increased in the palmitate-treated cells. In contrast, ROS levels were not increased in the cells treated with the unsaturated FFA, oleate, and treatments with radical scavengers significantly reduced palmitate toxicity . This suggested that increased ROS production plays an important role in palmitate toxicity. We performed further studies with co-supplementation of FFAs. Exposure to 0.4 mM palmitate was also toxic to the cells (about 20% LDH release on day 2), while treatment with 0.4 mM palmitate and 0.3 mM oleate were non-toxic (LDH release < 1%). These results indicated that the production of ROS and LDH release were specific to palmitate treatment. NADH dehydrogenases (mitochondrial complex I) has been suggested as an important source of ROS. The hierarchical model identified that 2 isoforms of NADH dehydrogenases had important roles in the cytotoxicity of the FFA. To experimentally validate their role in the cytotoxicity, we measured the ROS level and LDH release in palmitate cultured cells pretreated with irreversible NADH dehydrogenase inhibitor. Pretreatment with the NADH dehydrogenase inhibitor rotenone (0.5 uM) significantly reduced the ROS levels (Figure 10) as well as the LDH release (Figure 11), confirming the role of NADH dehydrogenase in the observed toxicity as predicted by the model. The cytotoxicity was, however, not completely prevented by the complex I inhibitor, suggesting the involvement of other, ROS-independent, mechanisms in the cytotoxicity.
b.ii. The role of aldehyde dehydrogenases
The analysis of regression coefficients suggested that many ALDH family members had positive association with the toxicity. Among these, ALDH1A1 was found to have the greatest positive regression coefficient to cytotoxicity (Table 3). Exposure to palmitate, but not to oleate, significantly increased the expression of ALDH1A1 (Figure 12). Therefore, we tested the effect of RNA interference (RNAi) of the ALDH1A1 gene on the toxicity. RNAi of ALDH1A1 significantly reduced its expression (Figure 13). The reduction in ALDH1A1 significantly reduced the caspase-3 activity of the cells in response to palmitate alone (Figure 14), confirming the cytotoxic role of ALDH1A1.
b.iii. The role of endothelial membrane protein 3 (EMP3)
EMP3 was identified as a potentially cytoprotective gene. This suggested that knocking down EMP3 would increase the toxicity of palmitate. To test this prediction, EMP3 was knocked down using RNAi (Figure 15). The RNA interference of EMP3 significantly increased the toxicity of palmitate treatment (Figure 16), verifying the predicted cytoprotective role of EMP3.
In this paper, a framework was developed to hierarchically integrate the metabolic and gene expression profiles to identify the genes which play important roles in determining the phenotypic responses. Application of this framework to the FFA and TNF-α induced cytotoxicity in HepG2 cells yielded novel targets to regulate the toxicity. The genes and pathways relevant to the cytotoxicity were identified and experimentally validated. For example, GSEA identified that mitochondrial alterations, but not ceramide synthesis, were associated with the toxicity. An advantage of incorporating GSEA analysis into our hierarchical approach is that incorporating knowledge-based information enhances the signal to noise ratio and the robustness of the analysis, and permits the detection of genes with modest changes . As illustrated by 1) the identification of ceramide metabolism not as an important player in the observed toxicity, which was experimentally validated, and 2) the novel targets identified by the hierarchical approach, suggest that the GSEA pathway analysis can compensate for the limited number of replicates to provide useful information. The former prediction is supported by the observation of no effect on the palmitate-induced toxicity upon inhibition of ceramide synthesis. The MBPLS prediction of the role of NADH dehydrogenase in regulating cytotoxicity was validated with complex I inhibitor studies, while those of ALDH1A1 and EMP3 were confirmed by RNAi studies.
Identification of genetic changes that control the phenotypic responses is an area of active research. This requires identification of genes that regulate the altered metabolic/physiological changes. Among the simpler models to identify such genes are the multivariate linear regression models . We also conducted multivariate linear regression analysis to relate the genetic and metabolic profiles, and found that the linear regression identified a much smaller number of genes relevant to a metabolite production or release which were assigned regression coefficients different from zero. For example, only 25 genes out of the 272 genes had weights different from zero in the linear regression model to predict BOH (see additional data file 4 for the list of genes). None of the important genes identified by MBPLS such as NADH dehydrogenases, glutathione S-transferases, ALDHs or EMP3 were selected. Linear regression model is ill-conditioned when the number of variables (genes) exceeds the number of observations (conditions), resulting in most of the regression coefficients taking on values of zero. MBPLS is a multivariate approach capable of modeling a large number of variables using a relatively small set of observations. It circumvents typical problems associated with the highly correlated and collinear nature of experimental data by projecting the data onto a few independent latent factors. These latent factors simplify the complex and diverse relationships by capturing the variable interactions contained in the original data in a new set of fewer unobserved/latent variables.
In this study, in addition to the generation of important information such as the roles of NADH dehydrogenases and the lack of involvement of de novo ceramide in the toxicity, two novel modulators of the palmitate-toxicity (ALDH1A1 and EMP3) were identified and validated. A notable point is the differences in the response of ALDH1A1 and EMP3 silencing. While silencing ALDH1A1 reduced caspase-3 activity significantly, it did not affect LDH release. On the contrary, silencing EMP3 increased LDH release without any affecting the caspase-3 activity. These differences in the cellular responses are most likely due to the cellular location of the corresponding proteins (cytosolic for ALDH1A1 and cell membrane for EMP3). Exposure to elevated levels of FFA would lead to increased omega oxidation, in which ALDHs play an important role. Omega oxidation can be a source of ROS and toxicity . Under these conditions, ALDH knock-down would reduce the omega oxidation and hence the cytotoxicity. In addition to its cytotoxic roles, ALDH1A1 may have potentially cytoprotective effects mediated by the detoxification of reactive lipid aldehydes . EMP3 is a protein of the peripheral membrane protein 22 (PMP22) family . The proteins of this family have been suggested to play important roles in cell proliferation and apoptosis . However, there has not been any study on the role of EMP3 in fatty acid toxicity. Its effect on the LDH release suggests the possibility that this protein may regulate membrane integrity.
In summary, this paper illustrated how phenotypic, metabolic and genetic profiles can be integrated hierarchically to identify phenotype relevant genes and pathways and novel targets. This approach identified the involvement of ROS generation, altered fatty acid and energy, but not ceramide metabolism in the cytotoxicity. In addition, novel targets such as ALDH1A1 and EMP3 were identified to modulate the toxicity of saturated FFA. Thus, the integration of metabolic and genetic information provides a more comprehensive picture of the perturbations as well as novel targets to regulate cellular responses.
HepG2/C3A cells and Fetal Bovine Serum were purchased from American Type Culture Collection (ATCC, Manassas, VA). Dulbecco's modified Eagle's medium with high glucose and no pyruvate (DMEM), Penicillin-Streptomycin (P/S), phosphate buffered saline (PBS, pH 7.4) and trizol reagent were purchased from Invitrogen (Carlsbad, CA). Fatty acid free bovine serum albumin (BSA) was purchased from MP Biomedicals (Chillicothe, OH). Sodium salts of all the fatty acids (palmitate, oleate and linoleate) were purchased from Sigma Aldrich chemical company (St. Louis, MO). 6-carboxy-2',7'-dichlorodihydrofluorescein diacetate, di(acetoxymethyl ester) (DCFDA dye) was obtained from Molecular Probes (Eugene, OR). Recombinant human TNF-α was from Peprotech (Rocky Hill, NJ).
One million HepG2/C3A cells were seeded into each well of a 6-well culture plate. The cells were cultured in 2 ml of medium containing DMEM supplemented with 10% fetal bovine serum (FBS) and 2% Penicillin-streptomycin (P/S). Cells were incubated at 37°C and in 10% CO2 atmosphere. After the cells reached confluence, the medium was replaced with 2 ml of the chosen medium, either HepG2, or the FFA medium containing 0.7 mM palmitate, oleate or linoleate; or the FFA-TNF-α medium. The fatty acids chosen are the most prevalent in the class of saturated (palmitate), monounsaturated (oleate) and polyunsaturated (linoleate) fatty acids in the plasma. The concentration of fatty acids chosen (0.7 mM) is commonly found in conditions of obesity. The FFAs were dissolved in 4% fatty acid-free BSA. Therefore, in addition to the HepG2 medium control, 4% fatty acid-free BSA in HepG2 medium was used as another control. TNF-α was added from a 100ug/ml stock in deionized water to make the desired final concentrations of either 20 or 100 ng/ml.
The cytotoxicity of the treatments was measured as the fraction of lactate dehydrogenase (LDH) released into the medium. LDH is a cytosolic enzyme which leaks out into the medium on exposure to toxic chemicals. Thus, the fraction of LDH released into the medium provides an accurate estimate of the cytotoxicity. For the LDH measurements, the cells were cultured in different media for 24 h and the supernatant was collected. Cells were washed with phosphate buffered saline (PBS) and kept in 1% triton-X-100 in PBS for another 24 h at 37°C. Cell lysate was then collected, vortexed for 15 seconds and centrifuged at 7000 rpm for 5 minutes. Cytotoxicity detection kit (Roche Applied Science, Indianapolis, IN) was used to measure the LDH release. LDH released was normalized to the total LDH (released + lyzed). The LDH released into the medium was normalized to total LDH, as shown in the equation below
Caspase-3 assay and TUNEL staining
Cells were treated with different FFA in the presence and absence of TNF-α. Caspase-3 activity was measured using a commercially available fluorescence-based assay kit from Biomol, according to manufacturer's protocol. Data were normalized to protein measurements from parallel experiments. For TUNEL staining, cells were cultured in Labtek Chamber Slide System. TUNEL staining was performed using the Dead-End fluorimetric TUNEL system from Promega Biosciences, according to manufacturer's protocol and imaged under a fluorescence microscope (Leica).
Measurement of metabolic uptake and production
The net uptake or production of a species was calculated by the difference in the concentration of the specie in the medium, before and after the treatment. The concentrations of metabolites were measured using enzymatic assays or HPLC. The amino acids were measured by HPLC using the AccQTag method (Waters) according to the manufacturer's protocol. Briefly, the media collected after treatments were de-proteinized by adding 4X acetonitrile (by volume) and keeping on ice for 1 h, after which the samples were centrifuged and the supernatant was derivatized using the reagent in the AccQTag kit as per the instructions given. The samples were analyzed with HPLC on a Waters 2690 separations module, using an AccQTag column (Waters) and fluorescence detection (Waters detector). Intracellular triglycerides were measured by lysing the cells with 1% triton-X-100 for 24 h and measuring the triglycerides by an enzymatic assay kit from Sigma. Glucose, lactate and glycerol were measured using enzymatic analysis kits from Sigma. Free fatty acid half micro kit (Roche Diagnostics) was used to find the concentration of free fatty acids in the medium. Beta-hydroxybutyrate (BOH) released was measured using a kit from Stanbio. Acetoacetate was measured using enzymatic assay .
Fisher's discriminant analysis
FDA identifies the projection axis that maximizes the ratio of the between-group and the within-group variations. Details on the FDA algorithm can be found in . FDA was applied to identify which among the 27 measured metabolic uptake/production contributed to the separation of the different phenotypes (cytotoxic versus nontoxic). Because there were 2 classes (toxic and non-toxic), a single discriminating vector is sufficient to separate them . The relative importance of the metabolites was identified by the projection of the metabolites in the new discriminant vector space.
Cells were cultured in 10 cm tissue culture plates until confluence and then exposed to different treatments for 24 h. RNA was isolated with the Trizol reagent. The gene expression profiles were obtained with the cDNA microarrays. The microarray analyses were conducted at the Van Andel Institute, Grand Rapids, MI. The protocols are available online at . There were two biological replicates for each condition and each replicate was measured with the Cy3 and Cy5 dyes (i.e. there were two technical replicates for each biological replicate). The microarray data has been deposited at the GEO website , with a query number of GSE5441.
Measurement of reactive oxygen species
Confluent cells were treated with the control medium or 0.7 mM palmitate for different time periods. Positive control cells were treated with 2 mM H2O2 for 1 h at 37°C. Measurement of ROS was performed by flow cytometry, using 6-carboxy-2',7'-dichlorodihydrofluorescein diacetate, di(acetoxymethyl ester) (C-2938, Molecular Probes, Eugene, OR) dye. A 5 mM stock solution of the dye was prepared in DMSO, and diluted to the 5 umol/L in DMEM. After the desired treatment times, the cells were exposed to the dye for 30 min, following which the cells were washed and imaged under a fluorescence microscope. Quantitative measurements were made by reading the fluorescence of the cells in a microplate fluorimeter with excitation at 488 nm and emission at 527 nm.
Measurement of ceramide
The cells were treated with 0.7 mM palmitate in the presence or absence of 20 μM Fumonisin B1 (FB1) for 24 h, following which the lipids were extracted by the method of Bligh and Dyer . The C-16 ceramide in the samples were then quantified with LC-MS as per a previously published protocol .
Gene Set Enrichment Analysis (GSEA) of the gene data
GSEA integrates a priori knowledge of a gene's functional role with the expression data to detect concerted expression changes in a set of genes responsible for producing a phenotype. The software GSEA-P, from , was used for the GSEA analysis. Thirty seven gene sets were selected from the molecular signature database, MsigDB  functional gene group c2, as shown in Table 1. These sets included 10 metabolic pathways, 26 signal pathways and 1 cellular component. An enrichment score of a gene set S characterizes whether the set of genes randomly distributed across the list or falls mainly at the bottom or top of the list. The null hypothesis that a gene set S randomly distributes across the ranked gene list was tested with Kolmogorov-Smirnov test, with the statistical significance value estimated through 1000 random permutations of the phenotype label. The gene sets with a high significance of enrichment are considered important in separating the distinct phenotypes.
Integrating the gene expression and metabolic profiles
MBPLS is a hierarchical multivariate analysis method [49, 50], where the variables are divided into different blocks based upon a priori knowledge, for example according to different stages of a process  or different metabolic pathways in a cell . For details of the MBPLS algorithm, refer . In the MBPLS model, the different gene sets identified by GSEA formed the blocks of the MBPLS. This ensured that the blocks were separated according to the functional role of the genes. Block scores were extracted from each block to predict the metabolic uptake/productions. This facilitated the identification of an important block (gene set) to a metabolic flux, and identified the important genes within the block. Important genes sets were identified by evaluating the weights of each block and importance of individual genes were then further identified by evaluating the regression coefficients of the genes within the block. Two latent variables were selected to be extracted from each block based upon the prediction accuracy of the MBPLS model. The N-way toolbox  was applied to conduct the MBPLS modeling.
RNA interference and reverse transfection
Silencer® Validated siRNAs targeting human EMP3 and ALDH1A1 mRNA were purchased from Ambion (Austin, TX). The synthesized oligonucleotides for siRNA of EMP3 are 5'-GUCCCUGAAUCUCUGGUACtt-3' and 5'-GUACCAGAGAUUCAGGGACtc-3', and the synthesized oligonucleotides for siRNA of ALDH1A1 are 5'-GGAACAGUGUGGGUGAAUUtt-3' and 5'-AAUUCACCCACACUGUUCCtg-3'. Reverse transfection of siRNA was performed. In general, siRNA and the transfection reagent, Lipofectamine RNAiMAX (Invitrogen, Carlsbad, CA), were diluted in Opti-MEM (Invitrogen) devoid of any serum and antibiotics and added to the 6-well plates. After 20 minutes, the same numbers of HepG2 cells in antibiotic free media were plated in each well and incubated at 37°C. After 24 h of transfection, the HepG2 cells were cultured for another 24 hours in regular media with other additives, for example palmitate. Then cells were harvested, washed twice with phosphate-buffered saline and lysed.
Overexpression and forward transfection
The ALDH1A3 plasmid, pCMV6-XL4-hALDH1A3, was purchased from Origene (Rockville, MD). Transient transfection was performed according to the Lipofectamine 2000 (Invitrogen) method. In general, the HepG2 cells were seeded in 6-well plates and cultured until reaching 80% confluency. Before transfection, the cells were washed twice with phosphate buffered saline, and medium was replaced with 2 ml of Opti-MEM medium. 1 μg/well of pCMV6-XL4-hALDH1A3 was then mixed with 5ul/well of Lipofectamine 2000 in Opti-MEM, and 20 minutes later, the mixture was added to the wells. After 6–12 hours of transfection, the cells were then cultured in regular media, treated with other additives like palmitate or TNF-α, and harvested after the treatment.
Real-time quantitative RT-PCR analysis
Total RNA was extracted from cells with an RNeasy mini kit and depleted of contaminating DNA with RNase-free DNase (Qiagen, Valencia, CA). Equal amounts of total RNA (1 μg) were reverse-transcribed using an iScript cDNA synthesis kit (Bio-RAD). The first-strand cDNA was used as a template. The primers used for quantitative RT-PCR analyses of human EMP3 (5'-GTGGTCTCAGCCCTTCACA-3' and 5'-ACGTGCAGTCGTACCAGAGA-3'), human ALDH1A1 (5'-AGCCTTCACAGGATCAACAGA-3' and 5'-GTCGGCATCAGCTAACACAA-3'), human ALDH1A3 (5'-GCCCTTTATCTCGGCTCTCT-3' and 5'-CGGTGAAGGCGATCTTGT-3') and human GAPDH (5'-AACTTTGGTATCGTGGAAGGA-3' and 5'-CAGTAGAGGCAGGGATGATGT-3') were synthesized by Operon Biotechnologies, Inc. (Huntsville, AL). RT-PCR was performed in 25-μl reactions using 1/10 of the cDNA produced by reverse transcription, 0.2 μM each primer, 1 X SYBR green supermix from Bio-RAD, and an annealing temperature of 57°C for 40 cycles. Each sample was assayed in three independent RT reactions and triplicate reactions each and normalized to GAPDH expression. Negative controls included the absence of enzyme in the RT reaction and absence of template during PCR. The cycle threshold (C T ) values corresponding to the PCR cycle number at which fluorescence emission in real time reaches a threshold above the base-line emission were determined using MyIQ™ Real-Time PCR Detection System.
- ACAC :
- ALDH :
- BOH :
- EMP3 :
endothelial membrane protein 3
- FFA :
free fatty acid
- GSEA :
gene set enrichment analysis
- GST :
- HNF :
hepatic nuclear factors
- LDH :
- MBPLS :
multi-block partial least squares analysis
- NASH :
- NF-κB :
nuclear factor kappa B
- PPAR :
peroxisome proliferator activated receptors
- ROS :
Reactive oxygen species
- SREBP :
sterol receptor element binding protein
- TNF-α :
tumor necrosis factor alpha
Felber JP, Golay A: Pathways from obesity to diabetes. Int J Obes Relat Metab Disord. 2002, 26 Suppl 2: S39-45. 10.1038/sj.ijo.0802126.
Kobayashi M: Molecular mechanism of insulin resistance. Saishin Igaku. 1998, 53 (6): 1210-1216.
Tilg H: Cytokines and liver diseases. Can J Gastroenterol. 2001, 15 (10): 661-668.
Watada H, Kawamori R: Insulin resistance and NASH. BIO Clinica. 2003, 18 (10): 874-879.
Randle PJ, Garland PB, Newsholme EA, Hales CN: The glucose fatty acid cycle in obesity and maturity onset diabetes mellitus. Ann N Y Acad Sci. 1965, 131 (1): 324-333. 10.1111/j.1749-6632.1965.tb34800.x.
Jump DB: Fatty acid regulation of gene transcription. Crit Rev Clin Lab Sci. 2004, 41 (1): 41-78. 10.1080/10408360490278341.
Cheung AT, Ree D, Kolls JK, Fuselier J, Coy DH, Bryer-Ash M: An in vivo model for elucidation of the mechanism of tumor necrosis factor-alpha (TNF-alpha)-induced insulin resistance: evidence for differential regulation of insulin signaling by TNF-alpha. Endocrinology. 1998, 139 (12): 4928-4935. 10.1210/en.139.12.4928.
Lang CH, Dobrescu C, Bagby GJ: Tumor necrosis factor impairs insulin action on peripheral glucose disposal and hepatic glucose output. Endocrinology. 1992, 130 (1): 43-52. 10.1210/en.130.1.43.
Schwabe RF, Brenner DA: Mechanisms of Liver Injury. I. TNF-alpha-induced liver injury: role of IKK, JNK, and ROS pathways. Am J Physiol Gastrointest Liver Physiol. 2006, 290 (4): G583-9. 10.1152/ajpgi.00422.2005.
Ding WX, Yin XM: Dissection of the multiple mechanisms of TNF-alpha-induced apoptosis in liver injury. J Cell Mol Med. 2004, 8 (4): 445-454. 10.1111/j.1582-4934.2004.tb00469.x.
Heyninck K, Wullaert A, Beyaert R: Nuclear factor-kappa B plays a central role in tumour necrosis factor-mediated liver disease. Biochem Pharmacol. 2003, 66 (8): 1409-1415. 10.1016/S0006-2952(03)00491-X.
Brenner DA: Signal transduction during liver regeneration. J Gastroenterol Hepatol. 1998, 13 Suppl: S93-5.
Tilg H, Diehl AM: Cytokines in alcoholic and nonalcoholic steatohepatitis. N Engl J Med. 2000, 343 (20): 1467-1476. 10.1056/NEJM200011163432007.
Rousset S, Bringuier A, Lardeux B, Feldmann G: Apoptosis induced by tumor necrosis factor a in human hepatoma cell lines. Falk Symposium. 113: 303-313.
Nagai H, Matsumaru K, Feng G, Kaplowitz N: Reduced glutathione depletion causes necrosis and sensitization to tumor necrosis factor-alpha-induced apoptosis in cultured mouse hepatocytes. Hepatology. 2002, 36 (1): 55-64. 10.1053/jhep.2002.33995.
Ji J, Zhang L, Wang P, Mu YM, Zhu XY, Wu YY, Yu H, Zhang B, Chen SM, Sun XZ: Saturated free fatty acid, palmitic acid, induces apoptosis in fetal hepatocytes in culture. Exp Toxicol Pathol. 2005, 56 (6): 369-376. 10.1016/j.etp.2005.02.003.
Wei Y, Wang D, Topczewski F, Pagliassotti MJ: Saturated fatty acids induce endoplasmic reticulum stress and apoptosis independently of ceramide in liver cells. Am J Physiol Endocrinol Metab. 2006, 291 (2): E275-81. 10.1152/ajpendo.00644.2005.
Srivastava S, Chan C: Hydrogen peroxide and hydroxyl radicals mediate palmitate-induced cytotoxicity to hepatoma cells: relation to mitochondrial permeability transition. Free Radic Res. 2007, 41 (1): 38-49. 10.1080/10715760600943900.
Feldstein AE, Werneburg NW, Li Z, Bronk SF, Gores GJ: Bax inhibition protects against free fatty acid-induced lysosomal permeabilization. Am J Physiol Gastrointest Liver Physiol. 2006, 290 (6): G1339-46. 10.1152/ajpgi.00509.2005.
Listenberger LL, Ory DS, Schaffer JE: Palmitate-induced apoptosis can occur through a ceramide-independent pathway. J Biol Chem. 2001, 276 (18): 14890-14895. 10.1074/jbc.M010286200.
Lu ZH, Mu YM, Wang BA, Li XL, Lu JM, Li JY, Pan CY, Yanase T, Nawata H: Saturated free fatty acids, palmitic acid and stearic acid, induce apoptosis by stimulation of ceramide generation in rat testicular Leydig cell. Biochem Biophys Res Commun. 2003, 303 (4): 1002-1007. 10.1016/S0006-291X(03)00449-2.
Gomez EO, Mendoza-Milla C, Ibarra-Sanchez MJ, Ventura-Gallegos JL, Zentella A: Ceramide reproduces late appearance of oxidative stress during TNF-mediated cell death in L929 cells. Biochem Biophys Res Commun. 1996, 228 (2): 505-509. 10.1006/bbrc.1996.1690.
Chan C, Hwang D, Stephanopoulos GN, Yarmush ML, Stephanopoulos G: Application of multivariate analysis to optimize function of cultured hepatocytes. Biotechnol Prog. 2003, 19 (2): 580-598. 10.1021/bp025660h.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005, 102 (43): 15545-15550. 10.1073/pnas.0506580102.
Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstrale M, Laurila E, Houstis N, Daly MJ, Patterson N, Mesirov JP, Golub TR, Tamayo P, Spiegelman B, Lander ES, Hirschhorn JN, Altshuler D, Groop LC: PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003, 34 (3): 267-273. 10.1038/ng1180.
Hwang D, Stephanopoulos G, Chan C: Inverse modeling using multi-block PLS to determine the environmental conditions that provide optimal cellular function. Bioinformatics. 2004, 20 (4): 487-499. 10.1093/bioinformatics/btg433.
Gill HK, Wu GY: Non-alcoholic fatty liver disease and the metabolic syndrome: effects of weight loss and a review of popular diets. Are low carbohydrate diets the answer?. World J Gastroenterol. 2006, 12 (3): 345-353.
Sanyal AJ, Campbell-Sargent C, Mirshahi F, Rizzo WB, Contos MJ, Sterling RK, Luketic VA, Shiffman ML, Clore JN: Nonalcoholic steatohepatitis: association of insulin resistance and mitochondrial abnormalities. Gastroenterology. 2001, 120 (5): 1183-1192. 10.1053/gast.2001.23256.
Listenberger LL, Han X, Lewis SE, Cases S, Farese RV, Ory DS, Schaffer JE: Triglyceride accumulation protects against fatty acid-induced lipotoxicity. Proc Natl Acad Sci U S A. 2003, 100 (6): 3077-3082. 10.1073/pnas.0630588100.
Takahashi Y, Campbell EA, Hirata Y, Takayama T, Listowsky I: A basis for differentiating among the multiple human Mu-glutathione S-transferases and molecular cloning of brain GSTM5. J Biol Chem. 1993, 268 (12): 8893-8898.
Moller IM: PLANT MITOCHONDRIA AND OXIDATIVE STRESS: Electron Transport, NADPH Turnover, and Metabolism of Reactive Oxygen Species. Annu Rev Plant Physiol Plant Mol Biol. 2001, 52: 561-591. 10.1146/annurev.arplant.52.1.561.
Taylor V, Suter U: Epithelial membrane protein-2 and epithelial membrane protein-3: two novel members of the peripheral myelin protein 22 gene family. Gene. 1996, 175 (1-2): 115-120. 10.1016/0378-1119(96)00134-5.
Jetten AM, Suter U: The peripheral myelin protein 22 and epithelial membrane protein family. Prog Nucleic Acid Res Mol Biol. 2000, 64: 97-129.
Dentin R, Benhamed F, Pegorier JP, Foufelle F, Viollet B, Vaulont S, Girard J, Postic C: Polyunsaturated fatty acids suppress glycolytic and lipogenic genes through the inhibition of ChREBP nuclear protein translocation. J Clin Invest. 2005, 115 (10): 2843-2854. 10.1172/JCI25256.
Bolon C, Gauthier C, Simonnet H: Glycolysis inhibition by palmitate in renal cells cultured in a two-chamber system. Am J Physiol. 1997, 273 (5 Pt 1): C1732-8.
Brookes PS: Mitochondrial H(+) leak and ROS generation: an odd couple. Free Radic Biol Med. 2005, 38 (1): 12-23. 10.1016/j.freeradbiomed.2004.10.016.
Parke DV: Mechanisms of chemical toxicity--a unifying hypothesis. Regul Toxicol Pharmacol. 1982, 2 (4): 267-286. 10.1016/0273-2300(82)90001-0.
Kim JS, He L, Lemasters JJ: Mitochondrial permeability transition: a common pathway to necrosis and apoptosis. Biochem Biophys Res Commun. 2003, 304 (3): 463-470. 10.1016/S0006-291X(03)00618-1.
Jia Z, Xu S: Clustering expressed genes on the basis of their association with a quantitative phenotype. Genet Res. 2005, 86 (3): 193-207. 10.1017/S0016672305007822.
Ockner RK, Kaikaus RM, Bass NM: Fatty-acid metabolism and the pathogenesis of hepatocellular carcinoma: review and hypothesis. Hepatology. 1993, 18 (3): 669-676.
Choudhary S, Xiao T, Vergara LA, Srivastava S, Nees D, Piatigorsky J, Ansari NH: Role of aldehyde dehydrogenase isozymes in the defense of rat lens and human lens epithelial cells against oxidative stress. Invest Ophthalmol Vis Sci. 2005, 46 (1): 259-267. 10.1167/iovs.04-0120.
Olsen C: An enzymatic fluorimetric micromethod for the determination of acetoacetate, -hydroxybutyrate, pyruvate and lactate. Clin Chim Acta. 1971, 33 (2): 293-300. 10.1016/0009-8981(71)90486-4.
Belhumeur PN, Hespanha JP, Kriegman DJ: Eigenfaces vs. Fisherfaces: Recognition Using Class Specific Linear Projection. IEEE Trans Pattern Analysis and Machine Intelligence. 1997, 19 (7): 711-720. 10.1109/34.598228.
cDNA microarry protocol at Van Andel Institute. http://www.vai.org/Research/Services/LMT/Protocols.aspx
GEO website . http://www.broad.harvard.edu/gsea/
Bligh EG, Dyer WJ: A rapid method of total lipid extraction and purification. Can J Biochem Physiol. 1959, 37 (8): 911-917.
Yamaguchi M, Miyashita Y, Kumagai Y, Kojo S: Change in liver and plasma ceramides during D-galactosamine-induced acute hepatic injury by LC-MS/MS. Bioorg Med Chem Lett. 2004, 14 (15): 4061-4064. 10.1016/j.bmcl.2004.05.046.
GSEA website. http://www.ncbi.nlm.nih.gov/geo
MacGregor JF, Jaeckle C, Kiparissides C, Koutoudi M: Process Monitoring and Diagnosis by Multi-Block PLS Methods. AIChE Journal. 1994, 40 (5): 826-838. 10.1002/aic.690400509.
Lopes JA, Menezes JC, Westerhuis JA, Smilde AK: Multiblock PLS analysis of an industrial pharmaceutical process. Biotechnol Bioeng. 2002, 80 (4): 419-427. 10.1002/bit.10382.
Andersson CA, Bro R: The N-way Toolbox for MATLAB. Chemometrics & Intelligent Laboratory Systems. 2000, 52 (1): 1-4. 10.1016/S0169-7439(00)00071-X.
This work is supported in part by the National Science Foundation (BES 0222747, BES 0331297, and 0425821), the National Institute of Health (1R01GM079688-01), the Environmental Protection Agency (RD83184701), the MSU Foundation and the Center for Systems Biology and the Whitaker Foundation.
ZL conceived the methodology and performed part of the experiments and wrote the manuscript. SS performed the cell culture, LDH measurement, microarray measurement, caspase, TUNNEL and metabolite measurement, wrote part of the manuscript. SM performed cell culture and extracted RNA for microarray measurement. XY performed ALDH and EMP3 experiment and wrote part of the manuscript. PN, JR, BH prepared, did QC/QC and performed cDNA microarry measurements. CC conceived the study and supervised the experiment and writing of the manuscript. All authors read and approved the final manuscript.
Zheng Li, Shireesh Srivastava contributed equally to this work.
Electronic supplementary material
Additional data file 1: Metabolic profiles. The data is a table of 27 measured metabolic uptake/productions. (XLS 28 KB)
Additional data file 2: BOH regression coefficients. Additional data file 2 is a table listing the regression coefficients of genes to BOH. (XLS 83 KB)
Additional data file 3: TG regression coefficients. Additional data file 3 is a table listing the regression coefficients of genes to TG. (XLS 77 KB)
Additional data file 4: BOH linear regression coefficients. Additional data file 4 is a table listing the regression coefficients to BOH with the linear regression model. (XLS 26 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Li, Z., Srivastava, S., Yang, X. et al. A hierarchical approach employing metabolic and gene expression profiles to identify the pathways that confer cytotoxicity in HepG2 cells. BMC Syst Biol 1, 21 (2007). https://doi.org/10.1186/1752-0509-1-21