Computer-Aided Design for Identifying Anticancer Targets in Genome-Scale Metabolic Models of Colon Cancer

Abstract

Simple Summary

Discovery of anticancer targets with minimal side effects is a major challenge in drug discovery and development. This study developed a fuzzy optimization framework for identifying anticancer targets. The framework was applied to identify not only gene regulator targets but also metabolite- and reaction-centric targets. The computational results show that the combination of a carbon metabolism target and any one-target gene that participates in the sphingolipid, glycerophospholipid, nucleotide, cholesterol biosynthesis, or pentose phosphate pathways is more effective for treatment than one-target inhibition is, and a two-target combination of 5-FU and folate supplement can improve cell viability, reduce metabolic deviation, and reduce side effects of normal cells.

Abstract

The efficient discovery of anticancer targets with minimal side effects is a major challenge in drug discovery and development. Early prediction of side effects is key for reducing development costs, increasing drug efficacy, and increasing drug safety. This study developed a fuzzy optimization framework for Identifying AntiCancer Targets (IACT) using constraint-based models. Four objectives were established to evaluate the mortality of treated cancer cells and to minimize side effects causing toxicity-induced tumorigenesis on normal cells and smaller metabolic perturbations. Fuzzy set theory was applied to evaluate potential side effects and investigate the magnitude of metabolic deviations in perturbed cells compared with their normal counterparts. The framework was applied to identify not only gene regulator targets but also metabolite- and reaction-centric targets. A nested hybrid differential evolution algorithm with a hierarchical fitness function was applied to solve multilevel IACT problems. The results show that the combination of a carbon metabolism target and any one-target gene that participates in the sphingolipid, glycerophospholipid, nucleotide, cholesterol biosynthesis, or pentose phosphate pathways is more effective for treatment than one-target inhibition is. A clinical antimetabolite drug 5-fluorouracil (5-FU) has been used to inhibit synthesis of deoxythymidine-5$${}^{\prime }$$-triphosphate for treatment of colorectal cancer. The computational results reveal that a two-target combination of 5-FU and a folate supplement can improve cell viability, reduce metabolic deviation, and reduce side effects of normal cells.

1. Introduction

The process of drug discovery and development is challenging as well as cost and time intensive. Recent progress in omics fields (e.g., genomics [1], transcriptomics [2], proteomics [3], metabolomics [4], and fluxomics [5]) can promote the development of technology for discovering new drug targets [6]. Advancements in high-throughput data acquisition have been combined with systems biology approaches to increase the time and cost effectiveness of drug target discovery through computer-aided simulation techniques. Metabolism is the primary biological mechanism linking genotype with phenotype; studying metabolism facilitates understanding of cell physiology and disease phenotypes caused by metabolic dysregulation [7,8,9,10,11]. Genome-scale metabolic networks (GSMNs) relate metabolites and reactions and represent the full set of intracellular metabolic processes curated using knowledge of cellular functions from the literature. GSMNs combined with constraint-based approaches are leading approaches for developing methods to simulate cell behavior, such as flux balance analysis (FBA) [12,13].

Systems biology, a holistic approach to the study of biological systems, involves the modeling and analysis of metabolic pathways, regulatory networks, and signal transduction networks to understand cellular behavior. Cellular metabolism is often altered during disease; therefore, metabolic analysis can facilitate drug discovery. During tumor development, the metabolic processes of cancer cells are altered to sustain their uncontrolled proliferation. Due to progress in research in the last decade, metabolic reprogramming has become a common focus of cancer studies [14,15]. Numerous studies have employed metabolic rewiring to understand cancer-specific metabolic networks and to predict anticancer targets that could impair tumor growth or viability [9,10,16,17,18,19,20,21].

Several publications have applied cancer-specific GSMNs to identify anticancer targets [7,9,16,21,22,23,24,25,26,27,28,29]. The integration of omics data, cancer-specific GSMNs, and FBA has recently utilized the heterogeneity of metabolic patterns to discover biomarkers of cancers [30]. Mapping these tissue-specific metabolisms in GSMNs provides deeper insight into the metabolic basis of physiological and pathological processes. Current context-specific-model-building algorithms can be broadly categorized into flux-dependent methods and pruning methods [16,22,23,28,31,32,33,34,35,36,37]. Flux-dependent methods identify an optimal GSMN and include the maximum number of high-confidence reactions as supported by substantial experimental data. By contrast, pruning methods start with a core set of reactions obtained through literature reviews or experimental data and proceed by removing reactions from the reconstruction while maintaining functionality in a core reaction set.

Colorectal cancer (CRC) is a worldwide health burden and it is the third most common type of cancer and the fourth most common cause of cancer-related death [38]. An estimated 51,020 deaths from CRC were reported in 2019 in the United States [38]. In 2017, the Taiwan Cancer Registry reported that CRC was the most frequently diagnosed cancer in men and the second most frequently diagnosed cancer in women [39]. In our previous study, the human metabolic network Recon 2.2 [40] was incorporated with the Human Protein Atlas (HPA) [3] to reconstruct GSMNs for cancerous colorectal tissue and its healthy counterpart. An oncogene inference optimization algorithm was introduced to integrate both models to predict which dysregulated genes cause tumorigenesis [41]. The algorithm was also applied to identify oncogenes for head-and-neck cancer [42] and non-small-cell lung carcinoma [43]. This study introduced a tri-level optimization framework (TLOP) for identifying anticancer targets (IACT) for treatment of cancers. The IACT framework, extended from the oncogene inference optimization algorithm [41,42], identified not only gene regulator targets but also metabolite- and reaction-centric targets. Fuzzy set theory was also applied to investigate cell viability, metabolic deviation, and side effects after target treatment. CRC was used as a case study to illustrate the performance of IACT.

2. Materials and Methods

2.1. Reconstruction of Tissue-Specific GSMNs

This study combined a human metabolic network (Recon 3D) [12] with RNA-Seq expression data from The Cancer Genome Atlas (TCGA) database [44] to reconstruct population-based tissue-specific GSMNs for CRC and healthy counterpart tissue. The RNA-Seq data of 41 healthy colorectal samples with FPKM-UQ normalized expression value and 478 colon adenocarcinoma samples with different TCGA barcode were downloaded from the TCGA database and normalized using quantile normalization; then the mean, confidence interval, and coefficient of dispersion were calculated for each gene. This information was used to evaluate supportive genes and classify them into four groups (high, medium, low, and not detected). Based on the four gene groups and Recon 3D, the CORDA algorithm was used to reconstruct GSMNs for cancer and healthy tissue. We developed a systems biology program to automatically build stoichiometric and gene-protein-reaction models in the General Algebraic Modeling System (GAMS) to discover the anticancer targets with few side effects. The workflow of the reconstruction of the GSMM for CRC is shown in Figure 1.

2.2. Optimization Framework for Target Identification

We aimed to identify anticancer targets that not only are lethal to cancer cells but also minimize the side effects of toxicity-induced tumorigenesis for normal cells and have reduced metabolic perturbation. We established a TLOP to mimic a wet-lab experiment for identifying targets. The flowchart of the in silico experiment is displayed in Figure 2. The optimal design concept of the IACT framework is described in Table 1. The TLOP is a hierarchical optimization problem with four objectives and subject to inner optimization problems describing the characteristics of cancer cells for targeting treatment and metabolic perturbation of normal cells caused by treatment. The first objective is to evaluate the mortality of cancer cells, a common criterion for discovering target problems [18,26,45,46,47], that the biomass growth rate of cancer cells for treatment (denoted as treated cells) has a value as small as possible. Anticancer targets may cause toxicity-induced tumorigenesis in normal cells and lead to harmful metabolic perturbations (referred to as metabolic perturbation). Therefore, the metabolic perturbation is defined as the normal cells accompanied with treatment (referred to as perturbed cells) to alter their metabolic flux distributions. The second objective is to obtain superior cell viability of perturbed cells, that is equivalent to maximize the ATP production rate and minimize the cell growth rate. We defined two types of metabolic deviations for perturbed cells to evaluate a side-effect grade. They are the differences of flux distributions of perturbed cells from cancer template and normal template, respectively. The third objective is to keep the metabolic deviations of perturbed cells as dissimilar as possible to the cancer template, and the fourth objective is to keep the metabolic deviations of perturbed cells as similar as possible to the normal template. The four objectives are formulated based on the fuzzy set theory [48,49,50,51,52], and detailed in Supplementary Materials (Figure S1). The constraint-based models for cancer and normal cells in inner optimization problems are expressed as follows:

(1)$\begin{array}{c}\mathrm{Inner}\mathrm{optimization}\mathrm{problems}:\hfill \\ \left\{\begin{array}{c}\mathrm{Treatment}\mathrm{of}\mathrm{cancer}\mathrm{cells}:\hfill \\ \begin{array}{c}\left\{\begin{array}{c}\mathrm{FBA}\mathrm{problem}\hfill \\ \underset{{\mathbf{v}}_{f/b}}{max}ob{j}_{CA}\equiv \left(w_{ATP}{v}_{ATP}+w_{biomass}{v}_{biomass}\right)\hfill \\ \mathrm{subject}\mathrm{to}\hfill \\ {\mathbf{N}}^{\mathbf{CA}}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=\mathbf{0}\hfill \\ v_{f/b,i}^{LB,TR}\le v_{f/b,i}\le v_{f/b,i}^{UB,TR},z_i\in {\Omega }^{TR}\hfill \\ v_{f/b,j}^{LB}\le v_{f/b,j}\le v_{f/b,j}^{UB},z_j\notin {\Omega }^{TR}\hfill \end{array}\right.\hfill \end{array}\hfill \\ \begin{array}{c}\left\{\begin{array}{c}\mathrm{UFD}\mathrm{problem}\hfill \\ \underset{{\mathbf{v}}_{f/b}}{min}\sum _{i\in {\Omega }^{Int}}(v_{f,i}^2+v_{b,i}^2)\hfill \\ \mathrm{subject}\mathrm{to}\hfill \\ {\mathbf{N}}^{\mathbf{CA}}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=\mathbf{0}\hfill \\ v_{f/b,i}^{LB,TR}\le v_{f/b,i}\le v_{f/b,i}^{UB,TR},z_i\in {\Omega }^{TR}\hfill \\ v_{f/b,j}^{LB}\le v_{f/b,j}\le v_{f/b,j}^{UB},z_j\notin {\Omega }^{TR}\hfill \\ obj\ge ob{j}_{CA}^{*}\hfill \end{array}\right.\hfill \end{array}\hfill \\ \mathrm{Perturbation}\mathrm{of}\mathrm{normal}\mathrm{cells}:\hfill \\ \begin{array}{c}\left\{\begin{array}{c}\mathrm{FBA}\mathrm{problem}\hfill \\ \underset{{\mathbf{v}}_{f/b}}{max}\mathrm{ob}{j}_{BL}\equiv \left(w_{ATP}{v}_{ATP}+w_{biomass}{v}_{biomass}\right)\hfill \\ \mathrm{subject}\mathrm{to}\hfill \\ {\mathbf{N}}^{\mathbf{BL}}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=\mathbf{0}\hfill \\ v_{f/b,i}^{LB,TR}\le v_{f/b,i}\le v_{f/b,i}^{UB,TR},z_i\in {\Omega }^{TR}\hfill \\ v_{f/b,j}^{LB}\le v_{f/b,j}\le v_{f/b,j}^{UB},z_j\notin {\Omega }^{TR}\hfill \end{array}\right.\hfill \end{array}\hfill \\ \begin{array}{c}\left\{\begin{array}{c}\mathrm{UFD}\mathrm{problem}\hfill \\ \underset{{\mathbf{v}}_{f/b}}{min}\sum _{i\in {\Omega }^{Int}}(v_{f,i}^2+v_{b,i}^2)\hfill \\ \mathrm{subject}\mathrm{to}\hfill \\ {\mathbf{N}}^{\mathbf{BL}}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=\mathbf{0}\hfill \\ v_{f/b,i}^{LB,TR}\le v_{f/b,i}\le v_{f/b,i}^{UB,TR},z_i\in {\Omega }^{TR}\hfill \\ v_{f/b,j}^{LB}\le v_{f/b,j}\le v_{f/b,j}^{UB},z_j\notin {\Omega }^{TR}\hfill \\ obj\ge ob{j}_{BL}^{*}\hfill \end{array}\right.\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}$

where ${\mathbf{v}}_{f/b}$ is the forward/backward flux vector of reversible reactions; the stoichiometric matrices ${\mathbf{N}}^{\mathbf{CA}}$ and ${\mathbf{N}}^{\mathbf{BL}}$ for tissue-specific cancer and normal cells, respectively, can be reconstructed by using Recon 3D [12] with the TCGA [44] or HPA [3] databases; $v_{f/b,i}^{LB}$ and $v_{f/b,i}^{UB}$ are the positive lower and upper bounds of the ith backward/forward flux, respectively; the integer vector $\mathbf{z}$ is used to determine mutated enzymes; and $ob{j}_{CA/BL}^{*}$ is the maximum cellular objective obtained from FBA. The aim of the IACT framework is to determine modulated reactions for metabolite-centric and reaction-centric approaches as well as for the gene-centric approach. The approaches are dependent on the restrictions for the lower and upper bounds $v_{f/b,i}^{LB,TR}$ and $v_{f/b,i}^{UB,TR}$ of the ith modulated reactions in the inner optimization problem. The restrictions on the bounds are discussed in Supplementary Materials (Figure S2).

2.3. Hierarchical Fitness in Outer Optimization

The IACT framework (Table 1) is expressed as a fuzzy multiobjective optimization problem (i.e., a hierarchical optimization problem). A bilevel optimization problem is a simple hierarchical optimization problem that converts the inner optimization problem into constraints in its outer-level problem by using duality theory. However, the inner problems in Equation (1) include two loops that are difficult to convert to constraints in the outer problem. The nested hybrid differential evolution (NHDE) algorithm was applied to predict oncogenes of various cancers [41,42]. This study extends the NHDE to solve the IACT problem. The computational procedures are presented in Supplementary Materials (Figure S3) and the implementation code in the GAMS (General Algebraic Modeling System) modeling language can be downloaded from https://chopin.ccu.edu.tw/?link=38a8SPuu6p (accessed on 1 July 2021). The outer optimization problem consists of three fuzzy goals and one crisp goal as shown in detail in Supplementary Materials (Figure S1). In fuzzy set theory [52], fuzzy objectives can be attributed membership functions to convert the objectives into decision criteria and thus convert the fuzzy optimization problem into a maximizing decision-making problem [50].

A linear membership grade was applied to normalize each objective between zero and one; thus, the multiobjective functions were converted to a hierarchical fitness function for evaluating the fitness of the NHDE algorithm. The definition and computation for each membership grade are described in Supplementary Materials (Figure S1), and the hierarchical fitness function is defined as follows:

(2)${\eta }_D=\frac{{\eta }_{CV}+min({\eta }_{CV},{\eta }_{DV})}{2}$

where the first priority grade ${\eta }_{CV}$ in the hierarchical fitness is the membership grade for the cell growth rates of the treated and perturbed cells. This grade is computed by the mean-min evaluation for both cells as follows:

(3)${\eta }_{CV}=\frac{\frac{{\eta }_{CV}^{TR}+{\eta }_{CV}^{PB}}{2}+min({\eta }_{CV}^{TR},{\eta }_{CV}^{PB})}{2}$

where the membership grade ${\eta }_{CV}^{TR}$ is the measure of the mortality of cancer cells in response to treatment, and the membership grade ${\eta }_{CV}^{PB}$ is an evaluation of cell viability for normal cells perturbed by treatment. This evaluation assesses the minimization of cell growth and maximization of ATP production for the perturbed cells and is computed by the mean-min evaluation of both grades as follows:

(4)${\eta }_{CV}^{PB}=\frac{\frac{{\eta }_{biomass}^{PB}+{\eta }_{ATP}^{PB}}{2}+min({\eta }_{biomass}^{PB},{\eta }_{ATP}^{PB})}{2}$

where ${\eta }_{biomass}^{PB}$ and ${\eta }_{ATP}^{PB}$ are the membership grades for the cell growth and ATP production of perturbed cells, respectively.

The second priority grade ${\eta }_{DV}$ in the hierarchical fitness function is used to evaluate metabolic deviations of the perturbed cells (the third and fourth goals of the framework) by using the following mean-min evaluation:

(5)${\eta }_{DV}=\frac{\frac{{\eta }_{S_F}+{\eta }_{S_M}+{\eta }_{L_F}+{\eta }_{L_M}+{\eta }_v+{\eta }_M}{6}+min({\eta }_{S_F},{\eta }_{S_M},{\eta }_{L_F},{\eta }_{L_M},{\eta }_v,{\eta }_M)}{2}$

where the membership grades ${\eta }_{S_F},{\eta }_{S_M},{\eta }_{L_F},$ and ${\eta }_{L_M}$ of the third goal are used to evaluate the maximization of the differences in the flux patterns of perturbed cells compared with templates generated from cancer and normal cells. The membership grades ${\eta }_v$ and ${\eta }_M$ of the fourth goal are used to measure the maximization of the flux and metabolite-flow similarities between perturbed cells and their normal counterparts. A higher membership grade for metabolic deviation implies a smaller metabolic perturbation due to treatment. The membership grade of side effects ${\eta }_{SE}$ is defined as follows:

(6)${\eta }_{SE}=\frac{\frac{{\eta }_{CV}^{PB}+{\eta }_{DV}}{2}+min({\eta }_{CV}^{PB},{\eta }_{DV})}{2}$

This membership grade was used to evaluate the metabolic deviation of the perturbed cells from templates generated from cancer and normal cells.

2.4. Factor Analysis

Log${}_2$ fold changes of all metabolite-flows for each perturbation (denoted as PB) to the normal (BL) state, and the template at cancer (CA) and normal states were described in Supplementary Materials (Figure S1), and thus,

(7)$\left\{\begin{array}{c}\mathrm{Perturbation}:L_{M,m}^a={log}_2\left(\frac{r_m^{P{B}_a}}{r_m^{BL}}\right),m\in {\Omega }^m,a\in {\Omega }^a\hfill \\ \mathrm{Template}:L_{M,m}^T={log}_2\left(\frac{r_m^{CA}}{r_m^{BL}}\right)\hfill \end{array}\right.$

where ${\Omega }^m$ is the set of metabolites in a GSMN, ${\Omega }^a$ is the set of identified anticancer targets, and the metabolite-flow rates $r_m^{PB/CA/BL}$ of the mth metabolite in each perturbation, cancer, and normal states were computed as follows:

(8)$r_m^{PB/CA/BL}=\sum _{i\in {\Omega }^c}\left(\sum _{N_{ij}>0,j}{N}_{ij}{v}_{f,j}^{PB/CA/BL}-\sum _{N_{ij}<0,j}{N}_{ij}{v}_{b,j}^{PB/CA/BL}\right),m\in {\Omega }^m$

${\Omega }^c$ is the set of metabolites located in different compartments of the cell. The expression enclosed in brackets in Equation (8) indicates the synthesis rate of the ith metabolite calculated by summing the influxes of the forward and backward reactions. Iterated principal factor analysis by using an orthogonal quartimax rotation method in SAS software (https://www.sas.com/, accessed on 1 July 2020) was used to analyze the ${log}_2$ fold changes $(L_{M,m}^T,L_{M,m}^a;a\in {\Omega }^a)$ of the template and all perturbations.

3. Results and Discussion

3.1. Reconstruction of Healthy and Cancerous Models

The GSMN of Recon 3D consisted of 5835 species, 10,600 reactions, and 2248 associated genes. We retrieved RNA-Seq data of 41 healthy colorectal samples with FPKM-UQ normalized expression value and 478 colon adenocarcinoma samples with different TCGA barcode from the TCGA database to reconstruct tissue-specific GSMNs for healthy (HT) and cancer (CA) states. The CORDA algorithm was used to reconstruct the HT and CA colorectal models. The HT model comprised 3742 species, 6023 reactions, and 1934 genes, and the CA model comprised 4402 species, 7027 reactions, and 1920 genes. Both models were merged into the basal (BL) or normal model that was the union set of the HT and CA models and included 4541 species, 7453 reactions, and 1986 genes. The numbers listed in the overlapping regions of Figure 3 denote the number of identical species, reactions, and genes for the HT and CA models.

3.2. Gene-Centric Approach

The NHDE algorithm [41,42] was first applied to solve the IACT problem by using a gene-centric approach to identify optimal targets. The NHDE algorithm determined a set of one-target genes having the highest hierarchical fitness among 1934 candidates and the best 20 target genes, as shown in Table 2. The determined 17 out of 20 target genes are the same genes predicted by using the NCI-60 cancer cell lines [16]. This study identified three new targets genes CRLS1, PGS1 and ADSS2 for treatment. We downloaded the dataset of cancer cell lines from the Cancer Dependency Map (DepMap, https://depmap.org/portal/, accessed on 1 July 2020), and 51 colon cancer cell lines from the dataset (2021Q1 version) were collected. Surveying the dataset, we observed that most of the target genes could cause cell death for a high percentage of colon cancer cell lines except for EBP, LSS, and NSDHL. These genes participate in the cholesterol biosynthesis III pathway. Using the STRING database (https://string-db.org/, accessed on 1 July 2020) and a Markov clustering method, we classified these gene-encoded enzymes into four classes of protein-protein interaction (PPI) (Figure 4A) that participate in the sphingolipid, glycerophospholipid, nucleotide, cholesterol biosynthesis, and pentose phosphate pathways.

On the basis of this computation, we found that cancer cells died (cell growth rate $\le {10}^{-10}$) for each one-target gene treatment and that although the growth of normal cells was nearly zero, the ATP production rate was 63% of the maximum level. The cell viability grade ${\eta }_{CV}^{PB}$ for each perturbed cell was 0.625. Thus, the cell viability grade for treated and perturbed cells reached 0.719. However, the metabolic deviation grades ${\eta }_{DV}$ of these genes were greater than 0.65, except for that of PGS1, PTDSS1, and DHODH, which were less than 0.448. Therefore, side effect grades ${\eta }_{SE}$ for PGS1, PTDSS1, and DHODH were lower than those of other genes. A lower ${\eta }_{SE}$ implies more side effects, that is, the normal cells have a higher chance of tumorigenesis or significant metabolic perturbation due to the treatment.

The NHDE is a genetic algorithm that can obtain and rank targets with higher grade. We used two groups of candidates in the algorithm to identify multiple targets for reducing computational burden. The first candidate group includes 20 identified targets in Table 2, and the second group includes the other candidate genes. We performed a series of computations to obtain a set of two-target combinations and determine their optimal grades as shown in Figure 5. For each combination, cancer cells died and the ATP production rate of the perturbed cells was greater than 95%. Thus, cell viability grades were greater than 0.944 for all treatments except for combinations including PCK1 (Figure 5). Most metabolic deviation grades ${\eta }_{DV}$ improved by approximately 5% compared with similar single-target treatments. All combinations with PCK1 achieved cell viability grades of at least 0.67. These results demonstrate that cell viability grade can be reduced to improve metabolic deviation and reduce side effects. We used the identified enzyme-encoding genes to investigate the PPI network (Figure 4B) by using the STRING database. The PPI network has five classes of interaction that are similar to interactions of the one-target enzymes. We observed that treatment predictions were superior for two-target combinations that had one target involved in glycosaminoglycan metabolism or central carbon metabolism, and the other target involved in any of the cholesterol, sphingolipid, glycerophospholipid, or nucleotide pathways. The gene-encoded enzyme PTDSS1 combined with any one target in central carbon metabolism does not have better performance than one-target inhibitors did (Figure 5). We found the two-target combination of PTDSS1 and PTDSS2 increase metabolic deviation grade to 0.674 (or side effect grade of 0.676), but decrease cell viability grade to 0.762. We also identified a three-target combination (PTDSS1, PTDSS2, and ENO1) with cell variability grade of 0.953, metabolic deviation grade of 0.689, and side effect grade of 0.751. The side effect grade for the three-target combination has been a significant improvement, and the performance of treatment is nearly identical to those of two-targets shown in Figure 5.

The metabolic pathway modulated by the identified genes is shown in Figure 6. These identified genes could be modulated by many drugs from DrugBank [53]. Of the drugs retrieved from DrugBank, 26 act on DHODH and 20 act on HMG-CoA reductase (HMGCR) (Table 2). DHODH catalyzes the oxidation of dihydroorotate (dhor-S) to orotate (orot) by using ubiquinone as an electron acceptor. Orotate is then catalyzed by uridine monophosphate synthetase (UMPS) to generate uridine monophosphate (ump), which is essential for the de novo production of pyrimidines for RNA and DNA replication. UMPS was also identified by the IACT framework as shown in Table 2. Both DHODH and UMPS inhibition may be effective for treating CRC [54,55]. Some literature has suggested that DHODH can be used to treat other diseases such as small-cell lung cancer [56], acute myeloid leukemia [57], and autoimmune diseases [58]. On the basis of the computation, we observed that DHODH achieved a lower side effect grade (${\eta }_{SE}=0.492$), implying that normal cells are more likely to undergo tumorigenesis or have significant metabolic perturbation due to the treatment.

DHODH was acted on by 26 drugs surveyed from DrugBank. Table 2 lists the number of available drugs for the identified one-target inhibitors. These drugs were used to investigate the grades of the adverse events by using a SIDER survey (http://sideeffects.embl.de/, accessed on 1 July 2020) and the ADDReSS (http://www.bio-add.org/ADReCS/, accessed on 1 July 2020) databases. The National Cancer Institute Common Terminology Criteria for Adverse Events (CTCAE) provide unique clinical descriptions of severity for adverse events (AE) from mild to death graded on a scale from 1 to 5. The average grade for each drug determined by using the CTCAE is shown in Supplementary Materials (Table S1). Of the 26 drugs acting on DHODH, 3 (leflunomide, atovaquone, and teriflunomide) have AE grades. The overall average AE grade was at most 10.04 (Table 2). This trend was consistent with the smaller computed ${\eta }_{SE}$ implying higher flux and metabolite-flow perturbations. Eight drugs (Supplementary Materials, Table S2) acting on HMGCR had an overall average AE grade of 8.77. The computed ${\eta }_{SE}$ was 0.637 for this target.

Two-target inhibitors, DHODH and PCK1, have been investigated as metabolic therapeutic targets for the treatment of CRC metastatic progression [55]. The IACT framework was applied to investigate the performance of the two-target enzymes in inhibiting treatments and to obtain a superior side effect grade (${\eta }_{SE}=0.634$), as shown in Figure 5. This side effect grade improved by approximately 29% compared with that of one-target DHODH. However, the PCK1 treatment was unable to inhibit cancer cell proliferation according to the computations. PCK1 combined with the other genes such as UMPS was also evaluated and found to have better grades (Figure 5).

As discussed previously, PCK1 could improve reaction synergism and catalyze the reversible decarboxylation and phosphorylation of oxaloacetate (oaa) to produce phosphoenolpyruvate (pep), as shown in Figure 6. The metabolite pep is also produced by the conversion of 2-phosphoglycerate (2pg) catalyzed by $\alpha$-enolase (ENO1). From Figure 5, we observed that two-target inhibition of DHODH and ENO1 had higher cell viability and metabolic deviation grades than the two-target inhibition of DHODH and PCK1 did. Combinations targeting ENO1 and the one-target genes in Table 2 also have higher grades. The identified target enzymes (e.g., GAPDH, PGK1, BPGM, and RPIA) involved in central carbon metabolism had similar results.

HMGCR, encoded by HMGCR, is a rate-controlling enzyme of the mevalonate pathway producing cholesterol and other isoprenoids (Figure 6). A total of 20 drugs targeting the enzyme HMGCR were retrieved from DrugBank. Cholesterol-lowering drugs targeting HMGCR are widely available and collectively known as statins. HMGCR has been investigated as an anticancer target for the clinical treatment of CRC, breast cancer, and other cancers [59,60,61,62]. Computations revealed that the one-target inhibitor could cause cancer cell toxicity and had a membership grade of 0.719 (Table 2). Two-target combinations targeting HMGCR and an enzyme participating in the central carbon metabolism (Figure 5) achieved higher cell viability grades (${\eta }_{CV}>0.94$) than did HMGCR one-target treatments. Moreover, three enzymes (MVK, MVD, and PMVK) participating in the mevalonate pathway could block cancer cell growth, and treatments targeting these enzymes could achieve results approaching that of those targeting HMGCR.

3.3. Metabolite-Centric Approach

The IACT framework was also used for a metabolite-centric approach for identifying antimetabolites. One-target antimetabolites and two-target antimetabolites for treatment of CRC were determined as shown in Table 3 and Figure 7, respectively. HMGCR catalyzes the conversion of hydroxymethylglutaryl coenzyme A (hmgccoa) to mevalonate (mev-R)-a necessary step in the biosynthesis of cholesterol (Figure 6). A cell viability grade greater than 0.71 and metabolic deviation grade greater than 0.66 could be achieved with a side effect grade of 0.63 by blocking either metabolite (hmgcoa or mev-R) (Table 3).

The metabolites orot, dtmp, and ump are nucleotides participating in the pyrimidine biosynthetic pathway of DNA synthesis. If these metabolites were inhibited, grades of 0.7, 0.66, and 0.63, were achieved for cell viability, metabolic deviation, and side effects, respectively. The results were nearly identical to those obtained using the gene-centric approaches (Table 2 and Figure 6). We also determined other one-target metabolites achieving satisfactory grades for cellular viability and side effects (Table 3). These metabolites included glu-L and ser-L in amino acids, pe-hs in a class of glycerophospholipids, and crm-hs and sphmyln-hs in the ceramide and sphingolipid families, respectively. Serine is a precursor for numerous other metabolites, including sphmyln-hs and folate (fol), and is the principal donor of one-carbon fragments in biosynthesis. Serine is also important in metabolism in that it participates in the biosynthesis of purines and pyrimidines. In this study, these antimetabolites were identified by the IACT framework, as displayed with gray symbols in Figure 6.

Although DHODH inhibition could cause cancer cell toxicity, worse side effects were predicted. Targeting DHODH and either PCK1 or ENO1 could improve efficacy. DHODH inhibition blocks orotate (orot) production, whereas PCK1 and ENO1 produce pep. Two-target antimetabolite inhibition achieved grades of 0.936, 0.695, and 0.751 for cell viability, metabolic deviation, and side effects, respectively (Figure 7). These grades were nearly equal to those of the two-target enzyme inhibition for DHODH and ENO1 but superior than those for DHODH and PCK1. Statins, also known as HMGCR inhibitors, are drugs that block the conversion of mev-R to hmgcoa and thereby inhibiting cholesterol synthesis. We found that the two-target combinations of hmgcoa or mev-R with a metabolite in the central carbon metabolism such as g3p, 3pg, 2pg, pep, or pyr could be applied to treat CRC satisfactorily (Figure 7). From the data displayed in Figure 7, we identified various two-target antimetabolite combinations and observed that the highest grades were for treatments targeting glu-L with central carbon metabolism metabolites such as g3p, 3pg, 2pg, pep, and pyr.

To verify the performance of the IACT framework, we evaluated cell viability and side effect grades of perturbation for normal cells treated with a clinical anticancer drug, 5-fluorouracil (5-FU). 5-FU is an antimetabolite drug commonly used for the treatment of CRC [63]. The drug inhibits the biosynthesis of dtmp (Figure 6) by acting on a key enzyme, thymidylate synthetase (TYMS). We used IACT for dtmp synthesis inhibition and computed a cell viability grade of 0.719 and a metabolic deviation grade of 0.665. We also evaluated two-target metabolite combination of dtmp and fol and found that inhibiting dtmp and enhancing fol synthesis achieved a cell viability grade of 0.978 but reduced the metabolic deviation grade to 0.597 (Figure 7). However, the side effect grade increased to 0.678. The reduction of ${\eta }_{DV}$ was consistent with observations in review articles [63,64]. To increase the anticancer activity of 5-FU, a modulation cotreatment with leucovorin to enhance intracellular levels of folate has been used [64,65,66]. This strategy has been demonstrated to increase the in vitro and in vivo toxicity of 5-FU for numerous cancer cell lines [67,68,69].

Although folate is a substrate for DNA synthesis in cell cycle S phase, the combination of folate and 5-FU results in synthetic anti-cancer effects due to the shift in the dump/dtmp metabolism and sensitization of 5-FU cytotoxicity. Such metabolism reprogramming can be predicted from our model. In clinical trial for CRC patients, similar results can also be observed especially in 5-FU continuous infusion which enhancing thymidylate synthase inhibition. The combination regimen of 5-FU/folinic acid increased anti-tumor effects without increasing toxicities while 5-FU was continuously infused. In EORTC40952 clinical trial for untreated metastatic CRC patients, Köhne et al. [70] showed significantly longer progression free survival in 5-FU continuous infusion plus leucovorin group (5.6 months vs. continuous 5-FU alone 4.1 months vs. bolus 5-FU plus leucovorin 4.0 months, $p=0.029$) and similar stomatitis was observed between continuous 5-FU plus leucovorin (5%) and continuous 5-FU alone (3%), while especially high in 5-FU bolus plus leucovorin group (11%).

3.4. Factor Loading of Identified Targets

The ${log}_2$ fold changes of metabolite-flows were used to form a $1478\times 201$ matrix that excluded unchanged flows and included the identified anticancer target genes, identified antimetabolites, and the template. The matrix data were analyzed through factor analysis to obtain nine factors and their associated factor loadings as listed in Supplementary Materials (Table S3). The key concept of factor analysis is that multiple observed variables have similar response patterns because they are all associated with a latent (i.e., not directly measured) variable. According to the analysis, the first and second factor loadings (Figure 8) were 64.19% and 19.32%, respectively. Therefore, these two factors could explain the key flux alterations for perturbations. The first group around Factor 1 displayed in the red circle of Figure 8 contained 117 anticancer target genes and antimetabolites from the two-target combinations and three-target combination of (PTDSS1, PTDSS2, ENO1). The flux alternations for genes in the first group (Supplementary Materials, Table S3) are very similar, because the values of cell viability grades, metabolic deviation grades, and side effect grades (Figure 5 and Figure 7) of all genes in the group are close to each other. The mean fold changes of metabolite-flows for the first group are displayed in the green boxes of Figure 6A. We observed that the metabolite-flows in the central carbon pathway were reduced in comparison with others. By contrast, the template increased as indicated in the yellow boxes of Figure 6A. The key metabolic reprogramming in cancer cells is rapid glucose metabolism to pyruvate, which is then largely converted to lactate. Glutamine is then replenished by the TCA cycle for proliferation requirement [14,15]. From Figure 6A, we observed that the flux alterations of the template coincide with the metabolic reprogramming of cancer progression. Moreover, the fold changes of metabolite-flows for the perturbed cells were different from the template.

The second group around Factor 2 contained 42 anticancer target genes and antimetabolites. Half of these were from two-target combinations. These two-target combinations involved the enzyme PCK1. The mean fold changes of metabolite-flow for the second group (Figure 6B) are little different from the results for the first group. Furthermore, the factor loadings could explain relationship with the membership grades for the identified targets. The factor loadings for PTDSS1, (PTDSS1, PTDSS2), and (PTDSS1, PTDSS2, ENO1) are shown by red dots in Figure 8. The identified targets have different cell viability, metabolic deviation, and side effect grades. The grades of the three-target combinations were very close to those of the others in the first group. Moreover, the grades for HMGCR (Table 2) are different to those (Figure 5) of the two-target combination (HMGCR, ENO1), so that they are separated into two groups in Figure 8. This result indicates that the ${log}_2$ fold changes of metabolite-flows for HMGCR can be discriminated from those of (HMGCR, ENO1).

4. Conclusions

A TLOP was designed to identify anticancer targets for the treatment of colon cancer. This IACT framework was applied to not only determine gene regulator drug targets but also discover metabolite- and reaction-centric targets. For the gene-centric approach, we determined one-target-gene-encoding enzymes that participate in the sphingolipid, glycerophospholipid, nucleotide, cholesterol biosynthesis, and pentose phosphate pathways. For two-target treatments, combinations of any of the one-target inhibitors and an enzyme (e.g., GAPDH, PGK1, ENO1, RPIA, BPGM, or PCK1) in the central carbon metabolism achieved better performance than one-target inhibitors did. For the metabolite-centric approach, 10 one-target metabolites (two fatty acids, three nucleotides, two amino acids, two sphingolipids, and one glycerophospholipid) were determined. These targets had nearly identical cell viability and side effects compared with those of one-target enzymes. To examine the performance of the IACT framework, cellular viability and side effects due to metabolic perturbation of normal cells after treatment with a clinical antimetabolite drug (5-FU) were investigated using the model. The computational results were consistent with 5-FU inhibiting dtmp synthesis to block cancer cell proliferation and with folate supplementation improving cell viability, slightly reducing metabolic deviation, and reducing side effects in comparison with the one-target 5-FU treatment [64,65,66]. These computational results are supported by in vitro experimental observations.

Author Contributions

Conceptualization, F.-S.W.; Data curation, C.-T.C., T.-Y.W. and P.-R.C.; Formal analysis, C.-T.C., T.-Y.W. and P.-R.C.; Funding acquisition, C.-Y.F.H. and F.-S.W.; Investigation, J.-M.L. and Y.-R.H.; Methodology, W.-H.W.; Project administration, F.-S.W.; Resources, P.M.-H.C. and C.-Y.F.H.; Software, W.-H.W. and F.-S.W.; Supervision, F.-S.W.; Validation, J.-M.L., P.M.-H.C. and Y.-R.H.; Visualization, C.-T.C., T.-Y.W. and P.-R.C.; Writing—original draft, F.-S.W.; Writing—review & editing, W.-H.W. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Ministry of Science and Technology of Taiwan (Grant MOST109-2320-B-030-007 to JML, MOST109-2320-B-075-003 to PMHC, MOST109-2320-B-010-026 to CYFH, MOST109-2320-B-037-032 and MOST109-2320-B-037-012 to YRH, and MOST109-2320-B-194-003 to FSW).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data generated or analyzed during this study are included in this published article.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

Symbol	Enzyme
ADSL	Adenylosuccinate lyase
ADSS2	Adenylosuccinate synthetase isozyme 2
BPGM	Bisphosphoglycerate mutase
CAD	CAD protein
CRLS1	Cardiolipin synthase (CMP-forming)
DHODH	Dihydroorotate dehydrogenase (quinone)
EBP	3-$$\beta$$-hydroxysteroid-$$\delta$$(8),$$\delta$$(7)-isomerase
ENO1	$$\alpha$$-enolase
FDFT1	Squalene synthase
G6PD	Glucose-6-phosphate 1-dehydrogenase
GAPDH	Glyceraldehyde-3-phosphate dehydrogenase
H6PD	GDH/6PGL endoplasmic bifunctional protein
HMGCR	3-hydroxy-3-methylglutaryl-coenzyme A reductase
KDSR	3-ketodihydrosphingosine reductase
LSS	Lanosterol synthase
MVD	Diphosphomevalonate decarboxylase
MVK	Mevalonate kinase
NSDHL	Sterol-4-alpha-carboxylate 3-dehydrogenase
PCK1	Phosphoenolpyruvate carboxykinase
PGK1	Phosphoglycerate kinase 1
PGS1	CDP-diacylglycerol-glycerol-3-phosphate 3-phosphatidyltransferase
PMVK	Phosphomevalonate kinase
PTDSS1	Phosphatidylserine synthase 1
RPIA	Ribose-5-phosphate isomerase
SPTLC1/2/3	A complex of serine palmitoyltransferase
SQLE	Squalene monooxygenase
UMPS	Uridine 5$${}^{\prime }$$-monophosphate synthase
Symbol	Metabolite
13dpg	3-Phospho-D-Glyceroyl Phosphate
2pg	2-Phospho-D-Glycerate
3dsphgn	(2S)-1-Hydroxy-3-Oxooctadecan-2-Aminium
3pg	3-Phospho-D-Glycerate
5dpmev	(R)-5-Diphosphomevalonate
5mthf	5-Methyltetrahydrofolate
5pmev	(R)-5-Phosphonatomevalonate
6pgc	6-Phospho-D-Gluconate
6pgl	6-Phospho-D-Glucono-1,5-Lactone
accoa	Acetyl Coenzyme A
ach	Acetylcholine
akg	2-Oxoglutarate
amet	S-Adenosyl-L-Methionine
amp	Adenosine Monophosphate
betald	Betaine Aldehyde
cdpdag-hs	Cytidine-5$${}^{\prime }$$-Diphosphate-Diacylglycerol
chol	Choline
cit	Citrate
clpn-hs	Cardiolipin
crm-hs	N-Acylsphingosine
ctp	Cytidine-5$${}^{\prime }$$-Triphosphate
dhap	Dihydroxyacetone Phosphate
dhcrm-hs	Dihydroceramide
dhf	7,8-Dihydrofolate
dhor-S	(S)-Dihydroorotate
dmgly	N,N-Dimethylglycine
dtdp	Deoxythymidine-5$${}^{\prime }$$-Diphosphate
dtmp	Deoxythymidine-5$${}^{\prime }$$-Phosphate
dttp	Deoxythymidine-5$${}^{\prime }$$-Triphosphate
dudp	Deoxyuridine-5$${}^{\prime }$$-Diphosphate
dump	Deoxyuridine-5$${}^{\prime }$$-Monophosphate
dutp	Deoxyuridine-5$${}^{\prime }$$-Triphosphate
f6p	D-Fructose 6-Phosphate
fdp	D-Fructose 1,6-Bisphosphate
fol	Folate
g1p	D-Glucose 1-Phosphate
g3p	Glyceraldehyde 3-Phosphate
g3pc	Glycerophosphocholine
g6p	D-Glucose 6-Phosphate
glc-D	D-glucose
gln-L	L-Glutamine
glu-L	L-Glutamate
gly	Glycine
gmp	Guanosine-5$${}^{\prime }$$-Monophosphate
hmgcoa	Hydroxymethylglutaryl Coenzyme A
icit	Isocitric Acid
imp	Inosine-5$${}^{\prime }$$-Monophosphate
ipdp	Isopentenyl Diphosphate
maloca	Malonyl Coenzyme A
met-L	L-Methionine
mev-R	(R)-Mevalonate
mlthf	5,10-Methylenetetrahydrofolate
oaa	Oxaloacetate
orot	Orotate
pa-hs	Phosphatidate
pail-hs	1-Phosphatidyl-1D-Myo-Inositol
pchol-hs	Phosphatidylcholine
pe-hs	Phosphatidylethanolamine
pep	Phosphoenolpyruvate
pglyc-hs	Phosphatidylglycerol
pgp-hs	Phosphatidyl Glycerol Phosphate
pmtcoa	Palmitoyl Coenzyme A
prpp	5-Phospho-Alpha-D-Ribose 1-Diphosphate
ps-hs	Phosphatidylserine
pyr	Pyruvate
r5p	Alpha-D-Ribose 5-Phosphate
ser-L	L-Serine
sphmyln-hs	Sphingomyelin
succoa	Succinyl Coenzyme A
thf	5,6,7,8-Tetrahydrofolate
udp	Uridine Diphosphate
ump	Uridine-5$${}^{\prime }$$-Monophosphate

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Roadmap of the reconstruction of genome-scale metabolic models. Workflows for reconstructing genome-scale metabolic models for CRC and its corresponding healthy tissue. (A) Download RNA-Seq data of CRC from the TCGA database. (B–G) A series statistical analyses of the download RNA-Seq data. (H) Download general human metabolic network model (Recon3D) from Bigg Models or VMH database. (I) Integrate the Recon3D model with the RNA-Seq data to classify all enzyme-encoded genes into four classes. (J) Retrieve gene-protein-reaction associations from the Recon3D model. (K) Compute the confidence score for each reaction based on the gene-protein-reaction associations, and classify all reactions into four groups. (L) Reconstruct tissue-specific metabolic models for healthy and cancer cells using the CORDA algorithm. (M) Create GAMS codes for tissue-specific metabolic models using the SBP platform.

Figure 2. Work flowchart for identifying anticancer target framework. (A) Tissue-specific genome-scale metabolic models of cancerous (CA) and normal (BL) cells were reconstructed through biological data. (B) Flux distribution patterns for cancer tissue can be provided from clinical data if available; otherwise the template can be computed through FBA and UFD problem without considered dysregulated restriction. (C) Flux distribution patterns for normal tissue can be provided from clinical data if available; otherwise the template can be computed through FBA and UFD problem without considered dysregulated restriction. (D) A set of anticancer targets are identified by the nest hybrid differential algorithm (NHDE), and provided to compute the flux distributions for each cancer treatment. (E) The same targets are provided to compute the perturbated flux distributions of normal cell during treatment. (F) Using cancer cell growth rate, cell mortality is evaluated. (G) Using membership function, cancer template and perturbated fluxes are used to compute dissimilarity grade. (H) Cell viability of perturbed cell is computed using ATP synthesis and cell growth rate. (I) Using membership function, normal template and perturbated fluxes are used to compute similarity grade. (J) The four-objective grades are used to evaluate fitness for making decision in the NHDE algorithm. (K) The next anticancer targets are generated in the NHDE algorithm if the fitness is unsatifactory, and repeat the procedures. (L) The optimal targets are obtained if the fitness is satisfactory.

Figure 3. Comparison of metabolic network data between HT and CA models. Statistics of cancer (CA) and healthy (HT) reconstructed metabolic models. The basal (BL) model is the union set of the CA and HT models.

Figure 4. Protein-protein interactions. Protein-protein interactions of identified (A) one-target anticancer genes and (B) two-target combinations. MCL clustering in the STRING database was applied to classify one-target enzymes into four classes and two-target enzymes to five classes. The first class contained nine enzymes in terpenoid backbone biosynthesis, the second class included five enzymes in metabolism of nucleotides, the third class had four enzymes in sphingolipid metabolism, and the fourth class included two enzymes in glycerophospholipid biosynthetic pathway. For the one-target case, RPIA in the pentose phosphate pathway was categorized in the first class, but that in the two-target case was considered to participate in glycosaminoglycan metabolism or central carbon metabolism.

Figure 5. Membership grades for two-target combinations. Membership grades of cell viability ([Forumla omitted. See PDF.]), metabolic deviation ([Forumla omitted. See PDF.]), and side effect ([Forumla omitted. See PDF.]) for two-target combinations of anticancer enzymes for colon cancer treatment. Cancer cell cytotoxicity was observed for each treatment. Therefore, the cell viability grade represents the cell growth viability of normal cells during treatment. Higher metabolic deviation grades indicate that the flux pattern of the perturbed cells was more dissimilar to the cancer template and more similar to the normal template. Higher side effect grades indicate fewer side effects. The numbers of drugs identified from DrugBank acting on each first target are shown in Table 2. The numbers of drugs acting on second targets are listed in brackets as follows: GAPDH (7), PGK1 (6), ENO1 (6), RPIA (1), BPGM (Not Available), and PCK1 (6). Error bar around each estimate was obtained through ten repeated executions.

Figure 6. Mean fold changes of metabolite-flows for perturbed cases and the template. The [Forumla omitted. See PDF.] fold changes of metabolite-flows for the identified targets in (A) the first group and (B) the second group. The identified gene-encoding enzymes are represented in red. The identified antimetabolites are presented in gray. The values in yellow boxes denote [Forumla omitted. See PDF.] fold changes of metabolite-flows [Forumla omitted. See PDF.] for the template, whereas the green box displays [Forumla omitted. See PDF.] for each perturbation. A positive value indicates that the metabolite-flow of cancer cells or the perturbation is higher than basal. A negative value represents lower flow.

Figure 7. Membership grades for two-target combinations of antimetabolites. Membership grades of cell viability ([Forumla omitted. See PDF.]), metabolic deviation ([Forumla omitted. See PDF.]), and side effect ([Forumla omitted. See PDF.]) for two-target combinations of antimetabolites for colon cancer treatment. Cancer cell cytotoxicity was observed for each treatment. Inhibition of two-target antimetabolites (except folate enhancement) improved membership grades compared with one-target counterparts. Glutamate combined with the metabolites in the central carbon metabolism such as g3p, 3pg, 2pg, pep, and pyr had the highest grades. Error bar around each estimate was obtained through ten repeated executions.

Figure 8. Membership grades for two-target combinations. First and second factor loadings of the identified anticancer target genes and antimetabolites. The template is the green dot. The first group around Factor 1 in the red circle contained 116 two-target combinations and a three combination of (PTDSS1, PTDSS2, ENO1) comprising anticancer target genes and antimetabolites. The second group around Factor 2 contained 42 anticancer target genes and antimetabolites. Half of these were two-target combinations. These two-target combinations included the enzyme PCK1.

Table 1

Optimization framework for IACT by evaluating the performance of identified targets according to four objectives.

Objectives in the outer optimization problem

1. The first objective is to measure the mortality of treated cancer cells.

2. The second objective is to obtain superior cell viability of normal cells in cancer treatment.

3. The third objective is to keep the metabolic deviation of perturbed cells as dissimilar as possible to the cancer template.

4. The fourth objective is to keep the metabolic deviation of perturbed cells as similar as possible to the normal template.

subject to the constraint-based models in the inner optimization problems

1. FBA and UFD problems for treated cancer cells.

2. FBA and UFD problems for perturbation of normal cells due to cancer treatment.

Table 2

Top 20 one-target genes obtained using the IACT framework. SPTLC1/2/3 is a complex of serine palmitoyltransferase constructed by SPTLC1, SPTLC2, and SPTLC3. Other genes each encode for a single enzyme, as shown in the abbreviations section. The symbol (–) denotes that data are unavailable.

Gene	$${\mathit{\eta }}_{\mathbf{CV}}{}^{†}$$	$${\mathit{\eta }}_{\mathbf{DV}}{}^{‡}$$	$${\mathit{\eta }}_{\mathbf{SE}}{}^{\S }$$	N/D $${}^{♭}$$	nDrugs $${}^{♮}$$	Ave. AE $${}^{♯}$$	Pathway $${}^{\mathbf{\P }}$$
HMGCR	0.719	0.675	0.637	50/51	20	8.77	Cholesterol Biosynthesis, Statin Pathway, Mevalonate Pathway
MVK	0.719	0.669	0.636	50/51	1	–	Mevalonate Pathway, Regulation of Cholesterol Biosynthesis By SREBP
MVD	0.719	0.669	0.636	48/51	–	–	Mevalonate Pathway, Regulation of Cholesterol Biosynthesis By SREBP
PMVK	0.719	0.669	0.636	23/51	–	–	Mevalonate Pathway, Regulation of Cholesterol Biosynthesis By SREBP
SQLE	0.719	0.672	0.637	3/51	4	6.73	Cholesterol Biosynthesis III, Statin Pathway
FDFT1	0.719	0.674	0.637	8/51	1	–	Cholesterol Biosynthesis III, Statin Pathway
EBP	0.719	0.674	0.637	0/51	1	8.94	Cholesterol Biosynthesis III
LSS	0.719	0.672	0.637	0/51	2	–	Cholesterol Biosynthesis III
NSDHL	0.719	0.674	0.637	0/51	1	–	Cholesterol Biosynthesis III
SPTLC1/2/3	0.719	0.670	0.636	37/51	2	–	Sphingolipid Metabolism
KDSR	0.719	0.670	0.636	6/51	–	–	Sphingolipid Metabolism
CRLS1	0.719	0.669	0.636	38/51	–	–	Glycerophospholipid Biosynthesis
PGS1	0.719	0.448	0.492	50/51	–	–	Glycerophospholipid Biosynthesis
PTDSS1	1.0	0.303	0.477	8/51	1	–	Glycerophospholipid Biosynthesis
ADSL	0.719	0.669	0.636	49/51	–	–	Metabolism of Nucleotides, Purine Metabolism
ADSS2	0.719	0.669	0.636	37/51	3	–	Metabolism of Nucleotides, Purine Metabolism
UMPS	0.719	0.655	0.632	29/51	2	9.82	Metabolism of Nucleotides, Pyrimidine Biosynthesis
DHODH	0.719	0.448	0.492	19/51	26	10.04	Metabolism of Nucleotides, Pyrimidine Biosynthesis
CAD	0.719	0.669	0.636	27/51	3	8.9	Metabolism of Nucleotides, Pyrimidine Biosynthesis
RPIA	0.719	0.675	0.637	6/51	1	–	Pentose Phosphate Pathway

$${}^{†}$$ Cell viability grade as evaluated from cancer cell treatment and the perturbations of normal cells due to treatment; ‡ Metabolic deviation grade indicating the perturbance of the cellular flux patterns as measured by dissimilarity to the cancer template and similarity to the basal template; $${}^{\S }$$ Side effect grade. A higher $${\eta }_{SE}$$ indicates fewer predicted side effects; $${}^{♭}$$ The cell death number (N) divided by the total number of colon cancer cells (D) used for the test from DepMap; $${}^{♮}$$ The number of drugs retrieved from DrugBank that modulate each gene; $${}^{♯}$$ Average grade of adverse events for drugs acting on an identified gene; $${}^{\P }$$ Accessed from GeneCards.

Table 3

Membership grades of cell viability ($${\eta }_{CV}$$), metabolic deviation ($${\eta }_{DV}$$), and side effect ($${\eta }_{SE}$$) for the top 10 one-target antimetabolites determined by the IACT framework.

Metabolite	Symbol	$${\mathit{\eta }}_{\mathbf{CV}}$$	$${\mathit{\eta }}_{\mathbf{DV}}$$	$${\mathit{\eta }}_{\mathbf{SE}}$$	Subclass
Hydroxymethylglutaryl Coenzyme A	hmgcoa	0.712	0.668	0.636	Fatty acyl thioesters
Mevalonate	mev-R	0.719	0.674	0.637	Fatty acids and conjugates
Orotate	orotsumpplementation	0.718	0.673	0.637	Pyrimidines and pyrimidine derivatives
Deoxythymidine-5$${}^{\prime }$$-Phosphate	dtmp	0.719	0.665	0.635	Pyrimidine deoxyribonucleotides
Uridine-5$${}^{\prime }$$-Monophosphate	ump	0.701	0.67	0.626	Pyrimidine ribonucleoside monophosphates
Phosphatidylethanolamine	pe-hs	0.795	0.691	0.7	Glycerophospholipids
L-Glutamate	glu-L	0.783	0.663	0.675	Amino acid
L-Serine	ser-L	0.789	0.656	0.672	Amino acid
N-Acylsphingosine	crm-hs	0.75	0.664	0.665	The parent compounds of the ceramide family
Sphingomyelin	sphmyln-hs	0.75	0.655	0.658	Phosphosphingolipids, fatty acyl group

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/biology10111115/s1, Figure S1: The mathematical formulation of the outer optimization problem in IACT framework, Figure S2: The restrictions on the bounds of the inner optimization problem, Figure S3: Introduction to the Nested Hybrid Differential Evolution (NHDE) algorithm, Table S1: Average AE grades for drugs acting on DHODH, Table S2: Average AE grades for drugs acting on HMGCR, Table S3: The $${log}_2$$ fold changes of metabolite-flows for the template, identified anticancer target genes, and identified antimetabolites.

References

1. Morozova, O.; Marra, M.A. Applications of next-generation sequencing technologies in functional genomics. Genomics; 2008; 92, pp. 255-264. [DOI: https://dx.doi.org/10.1016/j.ygeno.2008.07.001]

2. Barrett, T.; Wilhite, S.E.; Ledoux, P.; Evangelista, C.; Kim, I.F.; Tomashevsky, M.; Marshall, K.A.; Phillippy, K.H.; Sherman, P.M.; Holko, M. et al. NCBI GEO: Archive for functional genomics data sets–update. Nucleic Acids Res.; 2013; 41, pp. D991-D995. [DOI: https://dx.doi.org/10.1093/nar/gks1193]

3. Uhlen, M.; Oksvold, P.; Fagerberg, L.; Lundberg, E.; Jonasson, K.; Forsberg, M.; Zwahlen, M.; Kampf, C.; Wester, K.; Hober, S. et al. Towards a knowledge-based human protein atlas. Nat. Biotechnol.; 2010; 28, pp. 1248-50. [DOI: https://dx.doi.org/10.1038/nbt1210-1248] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21139605]

4. Wishart, D.S.; Feunang, Y.D.; Marcu, A.; Guo, A.C.; Liang, K.; Vázquez-Fresno, R.; Sajed, T.; Johnson, D.; Li, C.; Karu, N. et al. HMDB 4.0: The human metabolome database for 2018. Nucleic Acids Res.; 2018; 46, pp. D608-D617. [DOI: https://dx.doi.org/10.1093/nar/gkx1089]

5. Winter, G.; Krömer, J.O. Fluxomics-connecting ’omics analysis and phenotypes. Environ. Microbiol.; 2013; 15, pp. 1901-16. [DOI: https://dx.doi.org/10.1111/1462-2920.12064] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23279205]

6. Paananen, J.; Fortino, V. An omics perspective on drug target discovery platforms. Briefings Bioinform.; 2020; 21, pp. 1937-1953. [DOI: https://dx.doi.org/10.1093/bib/bbz122]

7. Agren, R.; Mardinoglu, A.; Asplund, A.; Kampf, C.; Uhlen, M.; Nielsen, J. Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling. Mol. Syst. Biol.; 2014; 10, 721. [DOI: https://dx.doi.org/10.1002/msb.145122]

8. Brul, S.; Angione, C. Human systems biology and metabolic modelling: A review–from disease metabolism to precision medicine. Biomed Res. Int.; 2019; 2019, 8304260.

9. Kell, D.B.; Goodacre, R. Metabolomics and systems pharmacology: Why and how to model the human metabolic network for drug discovery. Drug Discov Today; 2014; 19, pp. 171-182. [DOI: https://dx.doi.org/10.1016/j.drudis.2013.07.014]

10. Yizhak, K.; Chaneton, B.; Gottlieb, E.; Ruppin, E. Modeling cancer metabolism on a genome scale. Mol. Syst. Biol.; 2015; 11, 817. [DOI: https://dx.doi.org/10.15252/msb.20145307]

11. Zielinski, D.C.; Jamshidi, N.; Corbett, A.J.; Bordbar, A.; Thomas, A.; Palsson, B.O. Systems biology analysis of drivers underlying hallmarks of cancer cell metabolism. Sci. Rep.; 2017; 7, pp. 41241-41241. [DOI: https://dx.doi.org/10.1038/srep41241] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28120890]

12. Brunk, E.; Sahoo, S.; Zielinski, D.C.; Altunkaya, A.; Dräger, A.; Mih, N.; Gatto, F.; Nilsson, A.; Preciat Gonzalez, G.A.; Aurich, M.K. et al. Recon3D enables a three-dimensional view of gene variation in human metabolism. Nat. Biotechnol.; 2018; 36, pp. 272-281. [DOI: https://dx.doi.org/10.1038/nbt.4072] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29457794]

13. Heirendt, L.; Arreckx, S.; Pfau, T.; Mendoza, S.N.; Richelle, A.; Heinken, A.; Haraldsdóttir, H.S.; Wachowiak, J.; Keating, S.M.; Vlasov, V. et al. Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat. Protoc.; 2019; 14, pp. 639-702. [DOI: https://dx.doi.org/10.1038/s41596-018-0098-2]

14. Hanahan, D.; Weinberg, R. Hallmarks of cancer: The next generation. Cell; 2011; 144, pp. 646-74. [DOI: https://dx.doi.org/10.1016/j.cell.2011.02.013]

15. Pavlova, N.N.; Thompson, C.B. The emerging hallmarks of cancer metabolism. Cell Metab.; 2016; 23, pp. 27-47. [DOI: https://dx.doi.org/10.1016/j.cmet.2015.12.006]

16. Folger, O.; Jerby, L.; Frezza, C.; Gottlieb, E.; Ruppin, E.; Shlomi, T. Predicting selective drug targets in cancer through metabolic networks. Mol. Syst. Biol.; 2011; 7, 501. [DOI: https://dx.doi.org/10.1038/msb.2011.35] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21694718]

17. Ghaffari, P.; Mardinoglu, A.; Nielsen, J. Cancer metabolism: A modeling perspective. Front. Physiol.; 2015; 6, pp. 382-382. [DOI: https://dx.doi.org/10.3389/fphys.2015.00382]

18. Jerby, L.; Ruppin, E. Predicting drug targets and biomarkers of cancer via genome-scale metabolic modeling. Clin. Cancer Res.; 2012; 18, 5572. [DOI: https://dx.doi.org/10.1158/1078-0432.CCR-12-1856]

19. Lewis, N.E.; Abdel-Haleem, A.M. The evolution of genome-scale models of cancer metabolism. Front. Physiol.; 2013; 4, pp. 237-237. [DOI: https://dx.doi.org/10.3389/fphys.2013.00237] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24027532]

20. Nilsson, A.; Nielsen, J. Genome scale metabolic modeling of cancer. Metab. Eng.; 2017; 43, pp. 103-112. [DOI: https://dx.doi.org/10.1016/j.ymben.2016.10.022]

21. Robinson, J.L.; Nielsen, J. Anticancer drug discovery through genome-scale metabolic modeling. Curr. Opin. Syst. Biol.; 2017; 4, pp. 1-8. [DOI: https://dx.doi.org/10.1016/j.coisb.2017.05.007]

22. Agren, R.; Bordel, S.; Mardinoglu, A.; Pornputtapong, N.; Nookaew, I.; Nielsen, J. Reconstruction of genome-scale active metabolic networks for 69 human cell types and 16 cancer types using INIT. PLoS Comput. Biol.; 2012; 8, e1002518. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1002518] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22615553]

23. Bintener, T.; Pacheco, M.P.; Sauter, T. Towards the routine use of in silico screenings for drug discovery using metabolic modelling. Biochem. Soc. Trans.; 2020; 48, pp. 955-969. [DOI: https://dx.doi.org/10.1042/BST20190867]

24. Frezza, C.; Zheng, L.; Folger, O.; Rajagopalan, K.N.; MacKenzie, E.D.; Jerby, L.; Micaroni, M.; Chaneton, B.; Adam, J.; Hedley, A. et al. Haem oxygenase is synthetically lethal with the tumour suppressor fumarate hydratase. Nature; 2011; 477, pp. 225-228. [DOI: https://dx.doi.org/10.1038/nature10363]

25. Li, L.; Zhou, X.; Ching, W.K.; Wang, P. Predicting enzyme targets for cancer drugs by profiling human metabolic reactions in NCI-60 cell lines. BMC Bioinf.; 2010; 11, 501. [DOI: https://dx.doi.org/10.1186/1471-2105-11-501] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20932284]

26. Pacheco, M.P.; Bintener, T.; Ternes, D.; Kulms, D.; Haan, S.; Letellier, E.; Sauter, T. Identifying and targeting cancer-specific metabolism with network-based drug target prediction. EBioMedicine; 2019; 43, pp. 98-106. [DOI: https://dx.doi.org/10.1016/j.ebiom.2019.04.046]

27. Resendis-Antonio, O.; Checa, A.; Encarnación, S. Modeling core metabolism in cancer cells: Surveying the topology underlying the warburg effect. PLoS ONE; 2010; 5, e12383. [DOI: https://dx.doi.org/10.1371/journal.pone.0012383]

28. Turanli, B.; Zhang, C.; Kim, W.; Benfeitas, R.; Uhlen, M.; Arga, K.Y.; Mardinoglu, A. Discovery of therapeutic agents for prostate cancer using genome-scale metabolic modeling and drug repositioning. EBioMedicine; 2019; 42, pp. 386-396. [DOI: https://dx.doi.org/10.1016/j.ebiom.2019.03.009]

29. Yizhak, K.; Gaude, E.; Le Dévédec, S.; Waldman, Y.Y.; Stein, G.Y.; van de Water, B.; Frezza, C.; Ruppin, E. Phenotype-based cell-specific metabolic modeling reveals metabolic liabilities of cancer. eLife; 2014; 3, e03641. [DOI: https://dx.doi.org/10.7554/eLife.03641] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25415239]

30. Jalili, M.; Scharm, M.; Wolkenhauer, O.; Damaghi, M.; Salehzadeh-Yazdi, A. Exploring the metabolic heterogeneity of cancers: A benchmark study of context-specific models. J. Pers. Med.; 2021; 11, 496. [DOI: https://dx.doi.org/10.3390/jpm11060496]

31. Jerby, L.; Shlomi, T.; Ruppin, E. Computational reconstruction of tissue-specific metabolic models: Application to human liver metabolism. Mol. Syst. Biol.; 2010; 6, 401. [DOI: https://dx.doi.org/10.1038/msb.2010.56]

32. Machado, D.; Herrgård, M. Systematic evaluation of methods for integration of transcriptomic data into constraint-based models of metabolism. PLoS Comput. Biol.; 2014; 10, e1003580. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1003580] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24762745]

33. Opdam, S.; Richelle, A.; Kellman, B.; Li, S.; Zielinski, D.C.; Lewis, N.E. A systematic evaluation of methods for tailoring genome-scale metabolic models. Cell Syst.; 2017; 4, pp. 318-329.e6. [DOI: https://dx.doi.org/10.1016/j.cels.2017.01.010]

34. Schultz, A.; Qutub, A.A. Reconstruction of tissue-specific metabolic networks using CORDA. PLoS Comput. Biol.; 2016; 12, e1004808. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1004808]

35. Vlassis, N.; Pacheco, M.P.; Sauter, T. Fast reconstruction of compact context-specific metabolic network models. PLoS Comput. Biol.; 2014; 10, pp. e1003424-e1003424. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1003424] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24453953]

36. Wang, Y.; Eddy, J.A.; Price, N.D. Reconstruction of genome-scale metabolic models for 126 human tissues using mCADRE. BMC Syst. Biol.; 2012; 6, 153. [DOI: https://dx.doi.org/10.1186/1752-0509-6-153]

37. Wu, M.; Chan, C. Human metabolic network: Reconstruction, simulation, and applications in systems biology. Metabolites; 2012; 2, pp. 242-53. [DOI: https://dx.doi.org/10.3390/metabo2010242]

38. Siegel, R.L.; Miller, K.D.; Jemal, A. Cancer statistics, 2019. CA A Cancer J. Clin.; 2019; 69, pp. 7-34. [DOI: https://dx.doi.org/10.3322/caac.21551]

39. Taiwan Health Promotion Administration. Taiwan Cancer Registry, Cancer Statistics: Incidence and Mortality Rates for the Top 10 Cancer in Taiwan, 2008–2014. 2018; Available online: http://tcr.cph.ntu.edu.tw/ (accessed on 24 November 2018).

40. Swainston, N.; Smallbone, K.; Hefzi, H.; Dobson, P.D.; Brewer, J.; Hanscho, M.; Zielinski, D.C.; Ang, K.S.; Gardiner, N.J.; Gutierrez, J.M. et al. Recon 2.2: From reconstruction to model of human metabolism. Metabolomics; 2016; 12, 109. [DOI: https://dx.doi.org/10.1007/s11306-016-1051-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27358602]

41. Wang, F.S.; Wu, W.H.; Hsiu, W.S.; Liu, Y.J.; Chuang, K.W. Genome-scale metabolic modeling with protein expressions of normal and cancerous colorectal tissues for oncogene inference. Metabolites; 2019; 10, 16. [DOI: https://dx.doi.org/10.3390/metabo10010016] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31881674]

42. Wu, W.H.; Li, F.Y.; Shu, Y.C.; Lai, J.M.; Chang, P.M.H.; Huang, C.Y.F.; Wang, F.S. Oncogene inference optimization using constraint-based modelling incorporated with protein expression in normal and tumour tissues. R. Soc. Open Sci.; 2020; 7, 191241. [DOI: https://dx.doi.org/10.1098/rsos.191241]

43. Wang, Y.T.; Lin, M.R.; Chen, W.C.; Wu, W.H.; Wang, F.S. Optimization of a modeling platform to predict oncogenes from genome-scale metabolic networks of non-small-cell lung cancers. FEBS Open Bio; 2021; 11, pp. 2078-2094. [DOI: https://dx.doi.org/10.1002/2211-5463.13231] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34137202]

44. National Cancer Institute of U.S. The Cancer Genome Atlas Program. Department of Health and Human Services. 2020; Available online: https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga (accessed on 1 July 2020).

45. Gatto, F.; Miess, H.; Schulze, A.; Nielsen, J. Flux balance analysis predicts essential genes in clear cell renal cell carcinoma metabolism. Sci. Rep.; 2015; 5, 10738. [DOI: https://dx.doi.org/10.1038/srep10738]

46. Megchelenbrink, W.; Katzir, R.; Lu, X.; Ruppin, E.; Notebaart, R.A. Synthetic dosage lethality in the human metabolic network is highly predictive of tumor growth and cancer patient survival. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 12217-12222. [DOI: https://dx.doi.org/10.1073/pnas.1508573112] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26371301]

47. Pratapa, A.; Balachandran, S.; Raman, K. Fast-SL: An efficient algorithm to identify synthetic lethal sets in metabolic networks. Bioinformatics; 2015; 31, pp. 3299-3305. [DOI: https://dx.doi.org/10.1093/bioinformatics/btv352]

48. Chen, Y.F.; Wang, F.S. Crisp and fuzzy optimization of a fed-batch fermentation for ethanol production. Ind. Eng. Chem. Res.; 2003; 42, pp. 6843-6850. [DOI: https://dx.doi.org/10.1021/ie0210107]

49. Huang, H.J.; Wang, F.S. Fuzzy decision-making design of chemical plant using mixed-integer hybrid differential evolution. Comput. Chem. Eng.; 2002; 26, pp. 1649-1660. [DOI: https://dx.doi.org/10.1016/S0098-1354(02)00098-4]

50. Hsu, K.C.; Wang, F.S. Fuzzy optimization for detecting enzyme targets of human uric acid metabolism. Bioinformatics; 2013; 29, pp. 3191-3198. [DOI: https://dx.doi.org/10.1093/bioinformatics/btt564]

51. Wang, F.S.; Jing, C.H.; Tsao, G.T. Fuzzy-decision-making problems of fuel ethanol production using a genetically engineered yeast. Ind. Eng. Chem. Res.; 1998; 37, pp. 3434-3443. [DOI: https://dx.doi.org/10.1021/ie970736d]

52. Zimmermann, H.J. Fuzzy set theory. WIREs Comp. Stat.; 2010; 2, pp. 317-332. [DOI: https://dx.doi.org/10.1002/wics.82]

53. Wishart, D.S.; Feunang, Y.D.; Guo, A.C.; Lo, E.J.; Marcu, A.; Grant, J.R.; Sajed, T.; Johnson, D.; Li, C.; Sayeeda, Z. et al. DrugBank 5.0: A major update to the DrugBank database for 2018. Nucleic Acids Res.; 2018; 46, pp. D1074-D1082. [DOI: https://dx.doi.org/10.1093/nar/gkx1037]

54. Satoh, K.; Yachida, S.; Sugimoto, M.; Oshima, M.; Nakagawa, T.; Akamoto, S.; Tabata, S.; Saitoh, K.; Kato, K.; Sato, S. et al. Global metabolic reprogramming of colorectal cancer occurs at adenoma stage and is induced by MYC. Proc. Natl. Acad. Sci. USA; 2017; 114, pp. E7697-E7706. [DOI: https://dx.doi.org/10.1073/pnas.1710366114] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28847964]

55. Yamaguchi, N.; Weinberg, E.M.; Nguyen, A.; Liberti, M.V.; Goodarzi, H.; Janjigian, Y.Y.; Paty, P.B.; Saltz, L.B.; Kingham, T.P.; Loo, J.M. et al. PCK1 and DHODH drive colorectal cancer liver metastatic colonization and hypoxic growth by promoting nucleotide synthesis. eLife; 2019; 8, e52135. [DOI: https://dx.doi.org/10.7554/eLife.52135]

56. Li, L.; Ng, S.R.; Colón, C.I.; Drapkin, B.J.; Hsu, P.P.; Li, Z.; Nabel, C.S.; Lewis, C.A.; Romero, R.; Mercer, K.L. et al. Identification of DHODH as a therapeutic target in small cell lung cancer. Sci. Transl. Med.; 2019; 11, [DOI: https://dx.doi.org/10.1126/scitranslmed.aaw7852]

57. Sykes, D.B. The emergence of dihydroorotate dehydrogenase (DHODH) as a therapeutic target in acute myeloid leukemia. Expert Opin. Ther. Targets; 2018; 22, pp. 893-898. [DOI: https://dx.doi.org/10.1080/14728222.2018.1536748] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30318938]

58. Abdel-Magid, A.F. Use of dihydroorotate dehydrogenase inhibitors for treatment of autoimmune diseases and cancer. ACS Med. Chem. Lett.; 2020; 11, pp. 2072-2074. [DOI: https://dx.doi.org/10.1021/acsmedchemlett.0c00466]

59. Clendening, J.W.; Penn, L.Z. Targeting tumor cell metabolism with statins. Oncogene; 2012; 31, pp. 4967-78. [DOI: https://dx.doi.org/10.1038/onc.2012.6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22310279]

60. Gray, R.T.; Loughrey, M.B.; Bankhead, P.; Cardwell, C.R.; McQuaid, S.; O’Neill, R.F.; Arthur, K.; Bingham, V.; McGready, C.; Gavin, A.T. et al. Statin use, candidate mevalonate pathway biomarkers, and colon cancer survival in a population-based cohort study. Br. J. Cancer; 2017; 116, pp. 1652-1659. [DOI: https://dx.doi.org/10.1038/bjc.2017.139] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28524155]

61. Lipkin, S.M.; Chao, E.C.; Moreno, V.; Rozek, L.S.; Rennert, H.; Pinchev, M.; Dizon, D.; Rennert, G.; Kopelovich, L.; Gruber, S.B. Genetic variation in 3-hydroxy-3-methylglutaryl CoA reductase modifies the chemopreventive activity of statins for colorectal cancer. Cancer Prev. Res.; 2010; 3, pp. 597-603. [DOI: https://dx.doi.org/10.1158/1940-6207.CAPR-10-0007]

62. Tian, S.; Siu, F.M.; Lok, C.N.; Fung, Y.M.E.; Che, C.M. Anticancer auranofin engages 3-hydroxy-3-methylglutaryl-coenzyme A reductase (HMGCR) as a target. Met. Integr. Biometal Sci.; 2019; 11, pp. 1925-1936. [DOI: https://dx.doi.org/10.1039/C9MT00185A]

63. Longley, D.B.; Harkin, D.P.; Johnston, P.G. 5-fluorouracil: Mechanisms of action and clinical strategies. Nat. Rev. Cancer; 2003; 3, pp. 330-338. [DOI: https://dx.doi.org/10.1038/nrc1074] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12724731]

64. Ishiguro, L.; Yang, M.; Sohn, K.J.; Streutker, C.J.; Grin, A.; Croxford, R.; Kim, Y.I. Folic acid supplementation adversely affects chemosensitivity of colon cancer cells to 5-fluorouracil. Nutr. Cancer; 2016; 68, pp. 780-90. [DOI: https://dx.doi.org/10.1080/01635581.2016.1170168] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27175995]

65. Handali, S.; Moghimipour, E.; Rezaei, M.; Ramezani, Z.; Kouchak, M.; Amini, M.; Angali, K.A.; Saremy, S.; Dorkoosh, F.A. A novel 5-Fluorouracil targeted delivery to colon cancer using folic acid conjugated liposomes. Biomed. Pharmacother.; 2018; 108, pp. 1259-1273. [DOI: https://dx.doi.org/10.1016/j.biopha.2018.09.128]

66. Matherly, L.H.; Czajkowski, C.A.; Muench, S.P.; Psiakis, J.T. Role for cytosolic folate-binding proteins in the compartmentation of endogenous tetrahydrofolates and the 5-formyl tetrahydrofolate-mediated enhancement of 5-fluoro-2^′-deoxyuridine antitumor activity in vitro. Cancer Res.; 1990; 50, pp. 3262-3269.

67. Focaccetti, C.; Bruno, A.; Magnani, E.; Bartolini, D.; Principi, E.; Dallaglio, K.; Bucci, E.O.; Finzi, G.; Sessa, F.; Noonan, D.M. et al. Effects of 5-Fluorouracil on Morphology, cell cycle, proliferation, apoptosis, autophagy and ROS production in endothelial cells and cardiomyocytes. PLoS ONE; 2015; 10, e0115686. [DOI: https://dx.doi.org/10.1371/journal.pone.0115686]

68. Hou, Z.y.; Tong, X.p.; Peng, Y.b.; Zhang, B.k.; Yan, M. Broad targeting of triptolide to resistance and sensitization for cancer therapy. Biomed. Pharmacother.; 2018; 104, pp. 771-780. [DOI: https://dx.doi.org/10.1016/j.biopha.2018.05.088]

69. Zeng, H.; Zhu, X.; Tian, Q.; Yan, Y.; Zhang, L.; Yan, M.; Li, R.; Li, X.; Wang, G.; Ma, J. et al. In vivo antitumor effects of carboxymethyl chitosan-conjugated triptolide after oral administration. Drug Deliv.; 2020; 27, pp. 848-854. [DOI: https://dx.doi.org/10.1080/10717544.2020.1770370] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32508161]

70. Köhne, C.H.; Wils, J.; Lorenz, M.; Schöffski, P.; Voigtmann, R.; Bokemeyer, C.; Lutz, M.; Kleeberg, C.; Ridwelski, K.; Souchon, R. et al. Randomized phase III study of high-dose fluorouracil given as a weekly 24-hour infusion with or without leucovorin versus bolus fluorouracil plus leucovorin in advanced colorectal cancer: European organization of Research and Treatment of Cancer Gastrointestinal Group Study 40952. J. Clin. Oncol.; 2003; 21, pp. 3721-3728.