Identification of Anticancer Enzymes and Biomarkers for Hepatocellular Carcinoma through Constraint-Based Modeling

Hepatocellular carcinoma (HCC) results in the abnormal regulation of cellular metabolic pathways. Constraint-based modeling approaches can be utilized to dissect metabolic reprogramming, enabling the identification of biomarkers and anticancer targets for diagnosis and treatment. In this study, two genome-scale metabolic models (GSMMs) were reconstructed by employing RNA sequencing expression patterns of hepatocellular carcinoma (HCC) and their healthy counterparts. An anticancer target discovery (ACTD) framework was integrated with the two models to identify HCC targets for anticancer treatment. The ACTD framework encompassed four fuzzy objectives to assess both the suppression of cancer cell growth and the minimization of side effects during treatment. The composition of a nutrient may significantly affect target identification. Within the ACTD framework, ten distinct nutrient media were utilized to assess nutrient uptake for identifying potential anticancer enzymes. The findings revealed the successful identification of target enzymes within the cholesterol biosynthetic pathway using a cholesterol-free cell culture medium. Conversely, target enzymes in the cholesterol biosynthetic pathway were not identified when the nutrient uptake included a cholesterol component. Moreover, the enzymes PGS1 and CRL1 were detected in all ten nutrient media. Additionally, the ACTD framework comprises dual-group representations of target combinations, pairing a single-target enzyme with an additional nutrient uptake reaction. Additionally, the enzymes PGS1 and CRL1 were identified across the ten-nutrient media. Furthermore, the ACTD framework encompasses two-group representations of target combinations involving the pairing of a single-target enzyme with an additional nutrient uptake reaction. Computational analysis unveiled that cell viability for all dual-target combinations exceeded that of their respective single-target enzymes. Consequently, integrating a target enzyme while adjusting an additional exchange reaction could efficiently mitigate cell proliferation rates and ATP production in the treated cancer cells. Nevertheless, most dual-target combinations led to lower side effects in contrast to their single-target counterparts. Additionally, differential expression of metabolites between cancer cells and their healthy counterparts were assessed via parsimonious flux variability analysis employing the GSMMs to pinpoint potential biomarkers. The variabilities of the fluxes and metabolite flow rates in cancer and healthy cells were classified into seven categories. Accordingly, two secretions and thirteen uptakes (including eight essential amino acids and two conditionally essential amino acids) were identified as potential biomarkers. The findings of this study indicated that cancer cells exhibit a higher uptake of amino acids compared with their healthy counterparts.


Introduction
The liver is a pivotal metabolic organ that processes and transforms nutrients, stores energy, and synthesizes essential molecules.Hepatocytes, constituting approximately 85% of the liver's total mass, perform the majority of these metabolic processes [1].Hepatocellular carcinoma (HCC) is a liver disease wherein numerous metabolic processes are dysregulated, fostering tumorigenesis [2].Globally, HCC ranks as the third-leading cause of cancer-related mortality [3,4].The dysregulation of liver metabolism alters cellular metabolic pathways and is associated with cancer cells that appear different from normal cells [5].Similar to other cancer cells, HCC cells undergo metabolic changes, including aerobic glycolysis, diminished oxidative phosphorylation, and heightened biosynthetic intermediate generation to sustain rapid growth and proliferation [6][7][8][9][10].
Numerous studies have leveraged metabolic characteristics and reprogramming to comprehend cancer-specific metabolic networks, identify potential anticancer targets, and, thus, inhibit tumor growth or viability [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].Many drug discovery strategies, as noted in studies [11][12][13][14][15][16][17][18][19], have employed cancer cell mortality as a metric to gauge the effectiveness of the identified targets.However, the drug development process poses significant challenges due to its high expenses and prolonged timelines.A major hurdle is the efficient identification of both drug targets and potential side effects.Traditionally, side effects manifest primarily during pre-clinical or clinical trials, allowing for assessment of their severity.However, there is a notable lack of research focusing on predicting candidate target side effects in the early stages.In this study, we propose the introduction of a side effect index for target discovery procedures.This index aims to assess deviations between metabolic patterns of perturbed healthy cells in comparison to the respective cancer and healthy templates.
Tissue-specific genome-scale metabolic models (GSMMs) have been reconstructed through constraint-based modeling approaches to identify anticancer targets for treatment [13,15,[17][18][19][20][21][22][23][24][25][26][27][28] and predict biomarkers for diagnosis [16,[29][30][31].Synthetic lethality frequently serves as a therapeutic index in screening potential targets [25,28].Notably, few studies have attempted to predict the side effects of the identified targets during the early stages of drug discovery.Considering that the dysregulation of human cells can lead to metabolic perturbations or alterations in healthy cells, the present study used the cell mortality index to identify anticancer targets while accounting for metabolic perturbations in healthy cells to predict side effects during target discovery.However, such approaches lack prediction of metabolic perturbations for each identified target.Fuzzy similarity and fuzzy dissimilarity have been introduced as indicators of metabolic perturbations for use in predicting the side effects associated with the identified targets [24].Moreover, parsimonious flux variability analysis [32] was applied to evaluate the differential expression of metabolites between cancer cells and their healthy cells to identify potential biomarkers.

Weighting Factors in the UFD Problem
We extracted RNA-Seq expression data from 50 healthy liver samples with FPKM-UQ normalized data and 371 hepatocellular carcinoma samples, each with a distinct TCGA barcode, from the TCGA database [33].These data were used to reconstruct tissue-specific GSMMs for healthy (HT) and cancer (CA) states using the CORDA algorithm [34].The HT model comprised 4154 metabolites, 6652 reactions, and 2143 gene-encoding enzymes, whereas the CA model comprised 3982 metabolites, 6561 reactions, and 2107 gene-encoding enzymes.The GSMMs for the CA and HT models, stored in Excel format, are available in Supplementary File S1.
The Warburg effect, a common occurrence in the metabolic reprogramming of cancer cells, which involves an increase in the rate of glucose uptake and preferential production of lactate, even in the presence of oxygen [35], was observed, leading to the hypothesis that it could increase or decrease the production rate of metabolites [35][36][37].We examined the trends in flux changes based on the Warburg hypothesis found in the literature [6,8,20,37] and compared them with the computational results obtained for the pFBA problems using Equation (1).As expressed by Equation (2), quartile weighting factors were employed to solve the pFBA problems and calculate the flux distributions for both the CA and HT models.This approach led to a remarkable overlap of approximately 63.9% when compared with the flux change trends derived from the Warburg hypothesis.Additionally, pFBA evaluations using equal weightings (i.e., w CA/HT j = 1 for all reactions) resulted in a similarity of 59.1% with the trends observed in the Warburg hypothesis.

Single Anticancer Targets
We used Dulbecco's Modified Eagle Medium (DMEM), exhibiting 54 uptake reactions that were configured as reversible and exchangeable reactions, as illustrated in Supplementary File S1.Secretion reactions in the GSMM were configured as irreversible reactions.The NHDE algorithm [38] was applied to solve the MDM problem, to identify anticancer target genes.The algorithm can be used to identify all conceivable single-target anticancer enzymes and suppress biomass growth in HCC.Multiple iterations resulted in the identification of 22 single-target enzymes (refer to Table 1) from a pool of 2034 candidate enzymes.Downregulating these enzymes can suppress the biomass growth rate of cancer cells from 4% to 100%.The identified enzymes-ADSL, ADSS2, GK, CTSA*, and GMPR2-in Table 1 represent duplicate enzymes in GSMMs.CTSA* is a complex enzyme consisting of the four genes CTSA, GLB1, NEU1, and PSAP.The ACTD framework was employed to identify overall single-target genes and examine the protein-protein interaction (PPI) network based on the STRING (https://string-db.org/, accessed on 21 August 2023) and GeneCards (https://www.genecards.org/,accessed on 21 August 2023) databases.Subsequently, the network was divided into four distinct classes, as depicted in Figure 1A.The first class comprised eleven genes involved in cholesterol biosynthesis; the second class comprised five genes involved in nucleotide metabolism; the third class comprised five genes involved in glycerophospholipid biosynthesis; and the fourth class comprised four genes involved in sphingolipid metabolism.These targets enhanced the cancer cell elimination rate from 4% to 100% and reduced the maximum ATP synthesis rate of 38 mmol/gDW/h to 6.618 mmol/gDW/h for TR cells.The corresponding cell viability grade (η TR CV ) values ranged from 0.224 to 0.869.The metabolic deviation grade ( η TP MD for a treatment, an indicator of side effects, can be calculated based on the fuzzy similarity and fuzzy dissimilarity of the metabolic patterns for TR cells and PH cells relative to the metabolic templates of the HT and CA cells.The metabolic deviation grade values ranged from 0.269 to 0.479.As indicated by the data in Table 1, the highest cell viability grade ( ) of 0.896 was obtained for the respective treatment for the knockout of cytidine/uridine monophosphate As indicated by the data in Table 1, the highest cell viability grade (η TR CV ) of 0.896 was obtained for the respective treatment for the knockout of cytidine/uridine monophosphate kinase 1 (CMPK1) or downregulation of dCMP deaminase (DCTD).The metabolic deviation grade was 0.472 for the CMPK1 treatment and 0.451 for DCTD, respectively.Thus, CMPK1 induces fewer side effects than DCTD.The enzyme CMPK1 is involved in nucleotide metabolism, specifically in the phosphorylation of cytidine monophosphate (CMP) to form cytidine diphosphate (CDP).DCTD catalyzes the conversion of deoxycytidine monophosphate (dCMP) into deoxyuridine monophosphate (dUMP) by removing the amino group (NH2) from the cytosine base of dCMP.Both enzymes play a crucial role in nucleotide metabolism and are indirectly associated with the Gemcitabine pathway.Gemcitabine is a chemotherapeutic drug commonly used to treat various cancers.We examined DrugBank [39] to determine the number of drugs approved for the treatment of human diseases using the enzymes identified in Table 1.Nine drugs acting on CMPK1 and one on DCTD emerged as potential candidates for drug repurposing and treating HCC.

Various Nutrient Media
Ham's medium was employed as the nutrient source, with 68 uptake reactions (see Supplementary File S1) employed to identify anticancer targets for treating the CA model.Ten target genes were identified using the ACTD framework, as illustrated in Table 1.The computational results indicated that the cell viability grades obtained using Ham's medium were nearly equivalent to those derived from DMEM.However, the metabolic deviation grades with Ham's medium were slightly higher than those with DMEM.
The PPI network of the identified genes is depicted in Figure 1B.A comparison of Figure 1A,B reveals that enzymes associated with nucleotide metabolism, glycerophospholipid biosynthesis, sphingolipid metabolism, and the transport of vitamins and nucleotides reduced the biomass growth rate of treated cancer cells in Ham's medium.However, enzymes linked to cholesterol biosynthesis did not suppress the biomass growth (Table 1) and, thus, were not associated with the identified enzymes (Figure 1B).
The biomass growth suppression by an enzyme involved in cholesterol biosynthesis depends on whether the medium includes cholesterol uptake.To explore this further, we extended our analysis to three different media: VMH medium, human plasma-like medium (HPLM), and Roswell Park Memorial Institute (RPMI) medium.The uptake reactions in each medium are illustrated in Supplementary File S1.The VMH medium, sourced from the VMH database (https://www.vmh.life/#homeaccessed on 10 April 2023), was utilized to simulate human nutrient uptake.The ACTD framework identified potential single-target enzymes for the three media, with Table 1 summarizing the identified enzymes for each medium.The PPI network of the identified enzymes for these three media is illustrated in Figure 1C,D.Notably, the enzymes involved in cholesterol biosynthesis were unidentifiable using the VMH medium.
A comparison of nutrient uptake reactions in five different nutrient media revealed the absence of the cholesterol uptake reaction in DMEM, HPLM, and RPMI.We used five additional media to assess nutrient uptake: DMEM, HPLM, and RPMI, which exhibited the cholesterol uptake reaction and were denoted as DMEM+cholesterol, HPLM+cholesterol, and RPMI+cholesterol, respectively.In contrast, Ham and VMH were employed without the cholesterol uptake reaction, denoted as Ham−cholesterol and VMH−cholesterol, respectively.The ACTD framework was applied to identify single anticancer targets and investigate the relationships between biomass growth and different nutrient components across these nutrient media.The ACTD framework identified all potential single-target enzymes (Supplementary File S2) for the ten nutrient media.Table 2 lists the enzymes identified for each medium.The PPI networks of the identified enzymes for these media are depicted in Figure 2A-D.Our analysis revealed that anticancer enzymes involved in cholesterol biosynthesis were identifiable when cholesterol uptake did not occur in a medium, as illustrated in Figures 1A,D,E and 2B,C.Conversely, when the cholesterol uptake reaction occurred in the medium, the anticancer genes involved in cholesterol biosynthesis became unidentifiable, as illustrated in Figures 1B,C and 2A,D,E.Table 2. Cell viability (CV) grades and metabolic deviation (MD) grades of identified enzymes in various nutrient media with or without cholesterol uptake."--" indicates that the enzyme was unidentifiable by the nutrient medium used."*" indicates that CTSA is the representative of duplicate enzymes (CTSA, GLB1, NEU1, and PSAP).

Combination of Targets with Exchange Reactions
We incorporated a two-group representation within the ACTD framework to identify potential target combinations for treating HCC.Some of the single-target genes listed in Tables 1 and 2 could not reduce the biomass growth rate of TR cells satisfactorily.Therefore, a two-group approach was employed to identify target combinations, pairing a single-target gene from Table 1 with an additional exchange reaction.Numerous combinations were identified using various nutrient media, with or without the cholesterol uptake reaction, as detailed in Supplementary File S3.Table 3 presents select combinations using DMEM, Ham's, and VMH media from Supplementary File S3.Computational results revealed that the cell viability grades for all two-target combinations exceeded those for their corresponding one-target enzymes.Hence, combining a target enzyme with the modulation of an additional exchange reaction effectively reduced the biomass growth rate and ATP production rate in the treated cancer cells.Moreover, most two-target combinations exhibited an enhanced metabolic deviation grade relative to their corresponding single-target enzymes.
The cholesterol uptake reaction is pivotal in both Ham's and VMH media and renders single-target enzymes involved in the cholesterol pathway inadequate for treating HCC.However, incorporating one of these enzymes in an additional exchange reaction in these media can suppress the biomass growth of TR cells, thereby achieving a cell viability grade of 0.842 in Ham's medium and 0.980 in the VMH medium, as illustrated in Table 3.Thus, enabling an additional exchange reaction using a target enzyme can substantially enhance the therapeutic outcomes while minimizing the metabolic disturbances.As illustrated in Supplementary File S3, some combinations require at least two exchange reactions to improve the therapeutic outcomes, for example, knockout of NSDHL combined with two additional exchange reactions, as illustrated in Table 3.
Table 3. Combinations of various single-target enzymes with additional exchange reactions in a nutrient medium."↓" indicates that an enzyme is downregulated or an exchange reaction is to be taken up."∆" denotes that an enzyme is knocked out or an exchange reaction is to be blocked."↑" indicates that an enzyme is upregulated or an exchange reaction is to be secreted.

Biomarker Identification
We evaluated the flux distributions of the pFVA problem for the CA and HT models using ten nutrient media.Each nutrient medium yielded corresponding metabolite differential expression patterns between the CA and HT models.These patterns were used to calculate the mean of the distributions and determine the minimum and maximum values of secretion, uptake fluxes, and metabolite flow rates for the CA and HT models.The metabolites with differential expression for the CA and HT models with a log2 fold change greater than two are illustrated in Supplementary Files S4 and S5.Subsequently, potential biomarkers were identified using categories of "complete increase" or "complete decrease".However, some identified biomarkers depended on the nutrient medium used.We compiled the identified biomarkers shared across ten nutrient media, as illustrated in Figures 3 and 4. In the cancer state illustrated in Supplementary File S4, the secretion fluxes of 3 ′ -S-hydroxy-pravastatin-tetranor and acetoacetate were zero, indicating a significant decrease relative to those in the healthy state (Figure 3).The uptake of 3 ′ -S-hydroxypravastatin significantly decreased, whereas that of the others increased.According to our findings, two secretions and thirteen uptakes could serve as potential biomarkers.The uptakes included eight essential amino acids and two conditional essential amino acids, indicating that cancer cells exhibit higher uptake rates toward amino acids compared with healthy cells.

Biomarker Identification
We evaluated the flux distributions of the pFVA problem for the CA and HT models using ten nutrient media.Each nutrient medium yielded corresponding metabolite differential expression patterns between the CA and HT models.These patterns were used to calculate the mean of the distributions and determine the minimum and maximum values of secretion, uptake fluxes, and metabolite flow rates for the CA and HT models.The metabolites with differential expression for the CA and HT models with a log2 fold change greater than two are illustrated in Supplementary Files S4 and S5.Subsequently, potential biomarkers were identified using categories of "complete increase" or "complete decrease".However, some identified biomarkers depended on the nutrient medium used.We compiled the identified biomarkers shared across ten nutrient media, as illustrated in Figures 3 and 4. In the cancer state illustrated in Supplementary File S4, the secretion fluxes of 3′-S-hydroxy-pravastatin-tetranor and acetoacetate were zero, indicating a significant decrease relative to those in the healthy state (Figure 3).The uptake of 3′-S-hydroxy-pravastatin significantly decreased, whereas that of the others increased.According to our findings, two secretions and thirteen uptakes could serve as potential biomarkers.The uptakes included eight essential amino acids and two conditional essential amino acids, indicating that cancer cells exhibit higher uptake rates toward amino acids compared with healthy cells.with those identified in Figure 3. Additionally, 34 metabolites exhibited elevated levels, with thymidine, choline, and six amino acids among them exhibiting patterns similar to the uptake reactions depicted in Figure 3. Computational results indicated that the log2 fold changes of N-acylsphingosine and D-glucosyl-N-acylsphingosine between CA and HT cells were completely decreased using HPLM regardless of the occurrence or absence of the cholesterol uptake reaction; however, the levels of the other enzymes were elevated.Thus, N-acylsphingosine and D-glucosyl-N-acylsphingosine were influenced by the specific medium employed and, thus, cannot be considered candidate biomarkers.

Figure 4.
Log2 fold changes of metabolite flow rates between CA and HT cells.The minimum and maximum synthesis rates for each metabolite reaction in CA and HT cells were entirely increased or decreased, and the identified biomarkers were shared for ten nutrient media, except Nacylsphingosine and D-glucosyl-N-acylsphingosine."(+)" indicates that a nutrient medium included additional cholesterol uptake."(−)" indicates that a nutrient medium excluded cholesterol uptake.
Figure 4 illustrates the log2-fold changes in metabolite flow rates between CA and HT cells.The metabolite flow rate, as defined in Equation ( 4), corresponds to the overall synthesis rate of the mth metabolite in the GSMM.As illustrated in the classification of Figure 5, six metabolites exhibited a complete decrease.Among these metabolites, 3 ′ -Shydroxy-pravastatin-tetranor, 3 ′ -S-hydroxy-pravastatin, and acetoacetate were consistent with those identified in Figure 3. Additionally, 34 metabolites exhibited elevated levels, with thymidine, choline, and six amino acids among them exhibiting patterns similar to the uptake reactions depicted in Figure 3. Computational results indicated that the log2 fold changes of N-acylsphingosine and D-glucosyl-N-acylsphingosine between CA and HT cells were completely decreased using HPLM regardless of the occurrence or absence of the cholesterol uptake reaction; however, the levels of the other enzymes were elevated.Thus, N-acylsphingosine and D-glucosyl-N-acylsphingosine were influenced by the specific medium employed and, thus, cannot be considered candidate biomarkers.and  ⁄ denote the maximum and minimum levels for each flux and metabolite flow rate in the cancer and healthy state.(I) Groups in "complete decrease" or "complete increase" are identified as potential biomarkers.(J) The reconstructed GSMMs are used for the ACTD framework to identify the anticancer targets discussed in Figure 6.

Parsimonious Flux Balance Analysis and Flux Variability Analysis
The metabolic processes of HCC and its healthy counterpart were modeled through a constraint-based modeling approach to formulate it as a parsimonious flux balance analysis (pFBA) problem, which is an extension of flux balance analysis (FBA).In this modeling technique, the most efficient and minimal flux distribution through a metabolic network is determined to achieve a given objective.Notably, pFBA encompasses an FBA problem and a uniform flux distribution (UFD) problem, as follows: where N CA/HT represents the stochiometric matrices for the cancer (CA) and healthy (HT) models reconstructed based on a generic human network (Recon3D) [40] using RNA-Seq expression data from The Cancer Genome Atlas (TCGA) database [33].The parameters v LB f /b,i and v UB f /b,i represent the positive lower bound and positive upper bound of the ith forward flux and ith backward flux, respectively.The row vector c T,CA/HT contains raw weights indicating each reaction's contribution to the objective function for the CA and HT models.The weighting factors w CA/HT j in the UFD problem depend on confidence scores derived from gene-protein-reaction (GPR) associations.obj* is the maximal objective value obtained from the FBA problem, and the parameter ζ is assigned a value of 1 for pFBA.All reactions are classified using RNA-Seq expression of CA cells and HT cells and the GPR association of Recon3D into four types of confidence groups and can be assigned as follows: A small weighting factor indicates that a higher flux value in the UFD problem can be achieved when the corresponding reaction has a high confidence score in the GPR association.
The optimal fluxes for the pFBA problem are obtained with the fixed parameter ζ = 1, a value that can generally vary between 0 and 1.To account for this variability, we introduced parsimonious flux variability analysis (pFVA), which involves evaluating the metabolite differential expressions between CA and HT cells.The pFVA formulation is expressed as follows: where r m is the mth metabolite flow rate, and is expressed as follows: where Ω c is a set of metabolites located in various compartments of the cell and N m i ,j is the stoichiometric coefficient of the mth metabolite at the ith compartment participating in the jth reaction in a GSMM.The overall synthesis rates for the mth metabolite in various compartments are summed to yield the mth metabolite flow rate.Figure 5 illustrates a workflow using pFVA problems to calculate distributions for fluxes and metabolite flow rates toward achieving potential biomarkers.This workflow comprises two frameworks: one for cell-specific metabolic network reconstruction and another for biomarker identification.The initial framework involves employing reconstruction methods like CORDA [34] or iMAT [41] to reconstruct cell-specific GSMMs for CA and HT cells, as depicted in Figure 5A-D.The procedures for reconstruction are detailed in the literature [21][22][23][24].Both models are used to establish the pFVA problems to yield distributions for the flux values and metabolite flow rates for CA and HT cells (Figure 5E,F).The computational data are utilized to calculate the corresponding normal or Gaussian distributions for cancer and healthy cells.
As depicted in Figure 5G, a volcano plot illustrates the relationship between log2 fold change and their p-values.Biomarkers are commonly identified using log2 fold change and a p-value less than 0.05.However, this approach may yield numerous biomarkers, and determining clinical relevance becomes challenging when they fall within the overlap region of the distributions, making it difficult to discern their association with cancer or a healthy state.Utilizing these distributions, we established the minimum and maximum levels for each flux and metabolite flow rate in both cancer and healthy states.In this study, the flux variances shown in Figure 5G can be subsequently categorized into seven distinct groups, as shown in Figure 5H.A potential biomarker can be identified through comparative analysis: a flux value falling into the "complete increase" or "complete decrease" category indicates a significant difference in flux and metabolite flow rate between the cancer and healthy states, with no overlap in variability.

Anticancer Target Discovery Problem
In this study, we established an anticancer target discovery (ACTD) framework aimed at identifying potential anticancer targets, as depicted in Figure 6.This framework bears a resemblance to the antiviral target discovery platform.Tissue-specific genome-scale metabolic models (GSMMs) were reconstructed for HCC cells (referred to as CA cells) and their healthy counterparts (referred to as HT cells), as demonstrated in Steps A-D of Figure 5. Subsequently, both models were employed to formulate the pFBA problem, expressed in Equation ( 1), to assess flux distributions resulting from target treatment and the ensuing perturbations.
Conventional screening strategies apply the death of cancer cells as an evaluation criterion to discover targets.The ACTD framework not only reduces the growth rate of cancer cells but also minimizes the metabolic deviation of the metabolic patterns (Steps D and E) for treated cancer cells (referred to as TR cells) and perturbed healthy cells (referred to as PH cells) in comparison to their healthy counterparts.Metabolic patterns for TR and PH cells are represented as distributions of fluxes and metabolite flow rates resulting from target modulation.The metabolic templates (Steps F and G) of the CA and HT cells are defined as standard levels of fluxes and metabolite flow rates in their corresponding metabolic networks.Although these templates are typically obtained from clinical data, such a genome-scale standard is not currently available.In this study, a pFBA problem was formulated to compute the templates for CA and HT cells, respectively.The metabolic patterns were computed through Steps B and C for each candidate target and then compared against the metabolic templates to identify optimal targets.Subsequently, the metabolic patterns and templates were used to assess multi-objective fuzzy membership functions (Step H) and then transformed into a decision-making problem to maximize decision fitness (η D ).The fitness value for each anticancer candidate was used as a basis to identify the candidates that achieved a satisfactory grade (Step I).If the decision criterion was not satisfied, a subsequent set of anticancer candidates was generated using a nested hybrid differential evolution algorithm (Step J).The anticancer target was considered optimal if it had a good fit (Step K).

Fuzzy Multi-Objective Hierarchical Optimization Problem
The ACTD framework was formulated as a fuzzy multi-objective hierarchical optimization (FMHO) problem, as illustrated in Figure 7A, encompassing four fuzzy objectives.We defined a cell viability grade (η TR CV ) for TR cells as a criterion for evaluating cell mortality in cancer cells after treatment using a target.The cell viability grade for TR cells was determined through the fuzzy minimization of the biomass growth rate and ATP production rate.In contrast, a cell viability grade (η PH CV ) for PH cells was defined as a metric to evaluate the fuzzy minimization of the biomass growth rate and maximization of the ATP synthesis rate for PH cells due to treatment.More than 5000 fluxes and metabolite flow rates were acquired from the corresponding pFBA problems involving TR and PH cells.These values, denoted as metabolic patterns, were then contrasted with their templates to assess the disparities in metabolic distributions between the two models.We incorporated a novel concept of fuzzy dissimilarity to assess differences in the metabolic patterns of the TR and PH cells relative to the metabolic template of the CA cells.Additionally, the degree of similarity between the metabolic patterns of the TR and PH cells and the metabolic template of the HT cells was evaluated by employing fuzzy similarity.Fuzzy dissimilarity and similarity were determined using two-sided membership functions, which yielded a metabolic deviation grade (η TP MD ) that can be used as a metric for evaluating the side effects of treatment.The procedures for computation are detailed in the literature [24] and Supplementary File S6.
Molecules 2024, 29, 2594 15 of functions (Step H) and then transformed into a decision-making problem to maximi decision fitness (ηD).The fitness value for each anticancer candidate was used as a basis identify the candidates that achieved a satisfactory grade (Step I).If the decision criterio was not satisfied, a subsequent set of anticancer candidates was generated using a nest hybrid differential evolution algorithm (Step J).The anticancer target was considered o timal if it had a good fit (Step K).

Fuzzy Multi-Objective Hierarchical Optimization Problem
The ACTD framework was formulated as a fuzzy multi-objective hierarchical op mization (FMHO) problem, as illustrated in Figure 7A, encompassing four fuzzy obje tives.We defined a cell viability grade ( ) for TR cells as a criterion for evaluating c mortality in cancer cells after treatment using a target.The cell viability grade for TR ce was determined through the fuzzy minimization of the biomass growth rate and ATP pr duction rate.In contrast, a cell viability grade ( ) for PH cells was defined as a metric evaluate the fuzzy minimization of the biomass growth rate and maximization of the AT synthesis rate for PH cells due to treatment.More than 5000 fluxes and metabolite flo rates were acquired from the corresponding pFBA problems involving TR and PH cel These values, denoted as metabolic patterns, were then contrasted with their templates assess the disparities in metabolic distributions between the two models.We incorporat a novel concept of fuzzy dissimilarity to assess differences in the metabolic patterns of t TR and PH cells relative to the metabolic template of the CA cells.Additionally, the degr of similarity between the metabolic patterns of the TR and PH cells and the metabo template of the HT cells was evaluated by employing fuzzy similarity.Fuzzy dissimilari and similarity were determined using two-sided membership functions, which yielded metabolic deviation grade ( ) that can be used as a metric for evaluating the side effec of treatment.The procedures for computation are detailed in the literature [24] and Su plementary File S6.

Computational Procedures
Based on fuzzy set theory, the FMHO problem illustrated in Figure 7A can be transformed into a maximizing decision-making (MDM) problem, as depicted in Figure 7B.The MDM problem can be solved using the nested hybrid differential evolution (NHDE) algorithm [38].The optimality conditions of this transformation were proven in a previous study [38].According to the optimality conditions, a Pareto solution to the FMHO problem can be derived based on the transformed MDM problem.MDM problems are challenging optimization problems that cannot be solved directly using commercially available software [42,43].The high-dimensional, bilevel, and mixed-integer linear characteristics of these problems make them NP-hard.
In this study, we applied the NHDE algorithm, which is a stochastic and parallel direct search algorithm based on procedures in hybrid differential evolution (HDE) [44], an extension of the original differential evolution algorithm [45].The computational procedures of the NHDE algorithm are detailed in Table 4 and Supplementary File S6.The initialization process involved randomly generating a population of Np individuals (z i ) to cover the search uniformly.A structure array was used to represent the order number of a candidate target, its corresponding regulation, and the regulated strength parameter.Additionally, a two-group strategy was established to represent candidate targets and identify combinations of candidate genes and nutrient uptakes, as illustrated in Figure 8. Table 4. Computational procedures of the nested hybrid differential evolution algorithm to solve the maximizing decision-making problem.

1.
Represent and initialize for a population of target individuals (z 0 ) i = uniformInt z min , z max , i = 1, . . ., N p Each individual is represented as a structure array generated by a mixed-integer random array between z min and z max with a uniform distribution 2.
Use mutation with a rounding operation to generate a new population of target individuals (ẑ Use crossover operation to reproduce the next target individuals Use restriction operation to restrain a individual within its bounds Use selection and evaluation to determine the surviving individuals at the next generation where Ψ z i is a positive penalty to discipline the target individual (z i ) if its corresponding solution is infeasible.(c) Use one-to-one competition between the parent individual and its corresponding offspring, and select the best individual from the population of individuals.

6.
Use migration operation to prevent clustering around the best individual in a subsequent iteration, achieving a premature solution. 7.
Repeat steps 2 to 6 initialization process involved randomly generating a population of Np individuals (zi) to cover the entire search space uniformly.A structure array was used to represent the order number of a candidate target, its corresponding regulation, and the regulated strength parameter.Additionally, a two-group strategy was established to represent candidate targets and identify combinations of candidate genes and nutrient uptakes, as illustrated in Figure 8.

Conclusions
HCC, the predominant subtype of primary liver cancer, is a leading cause of cancerrelated mortality globally.Metabolic dysregulation is increasingly recognized as a significant risk factor for HCC, playing a crucial role in tumorigenesis and progression.Similar to other tumors, HCC tumors are marked by the neoplastic transformation of normal hepatocytes through substantial metabolic rewiring.This paper introduces an ACTD framework to identify potential biomarkers and anticancer targets for HCC treatment.The framework involves reconstructing tissue-specific GSMMs for both cancer cells and their healthy counterparts and employing them to formulate parsimonious flux variability analysis problems for biomarker identification.Furthermore, this study employed an FMHO framework to identify promising anticancer targets for HCC treatment.
Quartile weighting factors were employed in the pFBA problems for both models instead of equal weightings to calculate the flux patterns of the CA and HT models.Biomarkers are commonly identified using both flux patterns to evaluate log2 fold change and a p-value less than 0.05.This approach makes it difficult to discern the overlapping region of the distributions.We used confidence limits to categorize the flux patterns into seven groups and then to identify completely decreased or increased biomarkers.According to this classification, two secretions and thirteen uptakes were identified as potential biomarkers.Among these uptakes, eight were essential amino acids and two were conditional essential amino acids, indicating higher amino acid uptake rates in cancer cells relative to healthy cells.
The nutrient composition can significantly influence target identification.Within the ACTD framework, we employed ten different media to explore this effect and identify corresponding anticancer enzymes.The findings indicated that the identifiability of an enzyme participating in the cholesterol biosynthetic pathway may vary when it is used in conjunction with a cholesterol uptake medium.Computational results indicated that enzymes in the cholesterol biosynthetic pathway can be identified when a cholesterol uptake reaction does not occur in the medium used.Moreover, enzymes in the cholesterol biosynthetic pathway were identifiable in a cholesterol-free cell culture medium.A

Figure 2 .
Figure 2. Protein-protein interactions of anticancer target enzymes using various nutrient media with or without cholesterol uptake.(A) DMEM+cholesterol: DMEM included additional cholesterol uptake.(B) HAM-cholesterol: HAM excluded cholesterol uptake from the medium.(C) VMH-cholesterol: VMH excluded cholesterol uptake from the medium.(D) HPLM+cholesterol: HPLM included additional cholesterol uptake.(E) RPMI+cholesterol: RPMI included additional cholesterol uptake.

Figure 2 .
Figure 2. Protein-protein interactions of anticancer target enzymes using various nutrient media with or without cholesterol uptake.(A) DMEM+cholesterol: DMEM included additional cholesterol uptake.(B) HAM−cholesterol: HAM excluded cholesterol uptake from the medium.(C) VMH−cholesterol: VMH excluded cholesterol uptake from the medium.(D) HPLM+cholesterol: HPLM included additional cholesterol uptake.(E) RPMI+cholesterol: RPMI included additional cholesterol uptake.

Figure 3 .
Figure 3. Log2 fold changes of exchange reactions between CA and HT cells.The minimum and maximum flux values for each exchange reaction in CA and HT cells were entirely increased or decreased, and the identified biomarkers were shared for ten nutrient media."(+)" indicates that a nutrient medium included additional cholesterol uptake."(−)" indicates that a nutrient medium excluded cholesterol uptake.

Figure 3 .
Figure 3. Log2 fold changes of exchange reactions between CA and HT cells.The minimum and maximum flux values for each exchange reaction in CA and HT cells were entirely increased or decreased, and the identified biomarkers were shared for ten nutrient media."(+)" indicates that a nutrient medium included additional cholesterol uptake."(−)" indicates that a nutrient medium excluded cholesterol uptake.

Figure 5 .
Figure 5. Workflow for identifying biomarkers using differential expression patterns between CA and HT models through pFVA problems.(A) Use of a generic human genome-scale metabolic network to reconstruct tissue-specific metabolic models.(B) Utilization of biological databases to specify the context of nutrients and pathways.(C) Incorporation of RNA-sequence expression patterns of cancer cells and their healthy counterparts for reconstructing tissue-specific metabolic models.(D) Reconstruction of tissue-specific GSMMs and gene-protein-reaction models for CA cells and HT cells.(E) Formulation of the pFVA problem for the CA model to compute a set of flux distributions.(F) Formulation of the pFVA problem for the HT model to compute a set of flux distributions.(G) Computation of log2(CA/HT) with respect to the p-values.(H) Normal or Gaussian distributions are categorized into seven distinct groups.The red curve indicates the distributions in the cancer state, and the green curve indicates that of the healthy state. ⁄and  ⁄ denote the maximum and minimum levels for each flux and metabolite flow rate in the cancer and healthy state.(I) Groups in "complete decrease" or "complete increase" are identified as potential biomarkers.(J) The reconstructed GSMMs are used for the ACTD framework to identify the anticancer targets discussed in Figure6.

Figure 5 .
Figure 5. Workflow for identifying biomarkers using differential expression patterns between CA and HT models through pFVA problems.(A) Use of a generic human genome-scale metabolic network to reconstruct tissue-specific metabolic models.(B) Utilization of biological databases to specify the context of nutrients and pathways.(C) Incorporation of RNA-sequence expression patterns of cancer cells and their healthy counterparts for reconstructing tissue-specific metabolic models.(D) Reconstruction of tissue-specific GSMMs and gene-protein-reaction models for CA cells and HT cells.(E) Formulation of the pFVA problem for the CA model to compute a set of flux distributions.(F) Formulation of the pFVA problem for the HT model to compute a set of flux distributions.(G) Computation of log2(CA/HT) with respect to the p-values.(H) Normal or Gaussian distributions are categorized into seven distinct groups.The red curve indicates the distributions in the cancer state, and the green curve indicates that of the healthy state.r max CA/HT and r min CA/HT denote the maximum and minimum levels for each flux and metabolite flow rate in the cancer and healthy state.(I) Groups in "complete decrease" or "complete increase" are identified as potential biomarkers.(J) The reconstructed GSMMs are used for the ACTD framework to identify the anticancer targets discussed in Figure 6.

Figure 6 .
Figure 6.Framework for identifying the anticancer target discovery platform.(A) Reconstruction of tissue-specific GSMMs and gene-protein-reaction models for CA cells and HT cells.(B) Formulation of a parsimonious flux balance analysis problem for treating CA cells.(C) Formulation of a parsimonious flux balance analysis problem for perturbed HT cells.(D) Obtaining the metabolic pattern of treated CA cells for each anticancer candidate target by conducting the pFBA problem.(E) Obtaining the metabolic pattern of perturbed HT cells for each anticancer candidate target by conducting the pFBA problem.(F) Deriving the metabolic template of CA cells from clinical data (if available) or computed by conducting the pFBA problem without considering anticancer target regulation.(G) Deriving the metabolic template of HT cells from clinical data (if available) or computed by conducting the pFBA problem without considering anticancer target regulation.(H) Use of fuzzy set theory to convert fuzzy multiobjective functions into a fuzzy decision (ηD) in a maximizing decision-making problem.(I) Evaluation of the fitnesses of all anticancer candidates to make a decision.(J)Generating the subsequent set of anticancer candidates by using a nested hybrid differential evolution algorithm if the decision criterion is unsatisfied and then repeating Steps (B,J) until the maximum iteration assigned by the user is achieved.(K) The optimal targets are obtained if the fitness is satisfactory.

Figure 6 .
Figure 6.Framework for identifying the anticancer target discovery platform.(A) Reconstruction of tissue-specific GSMMs and gene-protein-reaction models for CA cells and HT cells.(B) Formulation of a parsimonious flux balance analysis problem for treating CA cells.(C) Formulation of a parsimonious flux balance analysis problem for perturbed HT cells.(D) Obtaining the metabolic pattern of treated CA cells for each anticancer candidate target by conducting the pFBA problem.(E) Obtaining the metabolic pattern of perturbed HT cells for each anticancer candidate target by conducting the pFBA problem.(F) Deriving the metabolic template of CA cells from clinical data (if available) or computed by conducting the pFBA problem without considering anticancer target regulation.(G) Deriving the metabolic template of HT cells from clinical data (if available) or computed by conducting the pFBA problem without considering anticancer target regulation.(H) Use of fuzzy set theory to convert fuzzy multiobjective functions into a fuzzy decision (η D ) in a maximizing decision-making problem.(I) Evaluation of the fitnesses of all anticancer candidates to make a decision.(J) Generating the subsequent set of anticancer candidates by using a nested hybrid differential evolution algorithm if the decision criterion is unsatisfied and then repeating Steps (B,J) until the maximum iteration assigned by the user is achieved.(K) The optimal targets are obtained if the fitness is satisfactory.

Figure 7 .
Figure 7.Transformation of a fuzzy multi-objective hierarchical optimization problem into a ma imizing decision-making problem using fuzzy membership functions.The lower bound (LB), upp bound (UB), and standard value (ST) are provided by a user.These values can be obtained fro clinical data (if available) or estimated from CA and HT templates.One-sided linear membersh functions are used to evaluate the fuzzy minimization and maximization (dashed and dot-dash lines, respectively).Two-sided linear membership functions are used to evaluate the fuzzy similar and fuzzy dissimilarity (red and green lines, respectively).

Figure 7 .
Figure 7. Transformation of a fuzzy multi-objective hierarchical optimization problem into a maximizing decision-making problem using fuzzy membership functions.The lower bound (LB), upper bound (UB), and standard value (ST) are provided by a user.These values can be obtained from clinical data (if available) or estimated from CA and HT templates.One-sided linear membership functions are used to evaluate the fuzzy minimization and maximization (dashed and dot-dashed lines, respectively).Two-sided linear membership functions are used to evaluate the fuzzy similarity and fuzzy dissimilarity (red and green lines, respectively).

( a )
For each target individual, solve the inner pFBA problem for the treated CA model and perturbed HT model, respectively.(b) Compute the fitness for each target individual (z i ) of the maximizing decision-making problem fitnes s

Figure 8 .
Figure 8. Representation of target individuals comprising candidate genes and nutrient uptakes.(A) Each two-group target is represented by a structure array that comprises a combination of order numbers for each two-group target, a regulated strength parameter (δ), and its corresponding regulation (up, down, or knockout).(B) The order numbers are mapped to the corresponding genes and nutrient uptakes.(C) The candidate genes are used to restrict the lower and upper bounds of the fluxes through the GPR association, and the candidate uptakes are also limited by their lower and upper bounds.

Table 4 .NHDE 1 .
Computational procedures of the nested hybrid differential evolution algorithm to solve the maximizing decision-making problem.Represent and initialize for a population of target individuals ( ) = uniformInt( , ),  = 1, … ,  Each individual is represented as a structure array generated by a mixed-integer random array between z min and z max with a uniform distribution 2. Use mutation with a rounding operation to generate a new population of target individuals

Figure 8 .
Figure 8. Representation of target individuals comprising candidate genes and nutrient uptakes.(A) Each two-group target is represented by a structure array that comprises a combination of order numbers for each two-group target, a regulated strength parameter (δ), and its corresponding regulation (up, down, or knockout).(B) The order numbers are mapped to the corresponding genes and nutrient uptakes.(C) The candidate genes are used to restrict the lower and upper bounds of the fluxes through the GPR association, and the candidate uptakes are also limited by their lower and upper bounds.

Table 1 .
Cell viability (CV) grades and metabolic deviation (MD) grades of identified enzymes in various nutrient media.D/T was obtained from the DepMap portal (https://depmap.org/portal/,accessed on 21 August 2023) and is defined as the ratio of the cell death number (N) and the total number (T) of colon cancer cell lines used in the experimental test."No.Drugs" denotes the number of drugs retrieved from DrugBank (https://go.drugbank.com/,accessed on 21 August 2023) that modulate each enzyme."--" indicates that the enzyme was unidentifiable by the nutrient medium used."*" indicates that CTSA is the representative of the duplicate enzymes (CTSA, GLB1, NEU1, and PSAP).
; and the fourth class comprised four genes involved in sphingolipid metabolism.These targets enhanced the cancer cell elimination rate from 4% to 100% and reduced the maximum ATP synthesis rate of 38 mmol/gDW/h to 6.618 mmol/gDW/h for TR cells.The corresponding cell viability grade ( ) values ranged from 0.224 to 0.869.The metabolic deviation grade ( ) for a treatment, an indicator of side effects, can be calculated based on the fuzzy similarity and fuzzy dissimilarity of the metabolic patterns for TR cells and PH cells relative to the metabolic templates of the HT and CA cells.The metabolic deviation grade values ranged from 0.269 to 0.479. biosynthesis

Table 2 .
Cell viability (CV) grades and metabolic deviation (MD) grades of identified enzymes in various nutrient media with or without cholesterol uptake."--" indicates that the enzyme was unidentifiable by the nutrient medium used."*" indicates that CTSA is the representative of duplicate enzymes (CTSA, GLB1, NEU1, and PSAP).