Fuzzy optimization for identifying anti-cancer targets with few side effects in constraint-based models of head and neck cancer

Computer-aided methods can be used to screen potential candidate targets and to reduce the time and cost of drug development. In most of these methods, synthetic lethality is used as a therapeutic criterion to identify drug targets. However, these methods do not consider the side effects during the identification stage. This study developed a fuzzy multi-objective optimization for identifying anti-cancer targets that not only evaluated cancer cell mortality, but also minimized side effects due to treatment. We identified potential anti-cancer enzymes and antimetabolites for the treatment of head and neck cancer (HNC). The identified one- and two-target enzymes were primarily involved in six major pathways, namely, purine and pyrimidine metabolism and the pentose phosphate pathway. Most of the identified targets can be regulated by approved drugs; thus, these drugs are potential candidates for drug repurposing as a treatment for HNC. Furthermore, we identified antimetabolites involved in pathways similar to those identified using a gene-centric approach. Moreover, HMGCR knockdown could not block the growth of HNC cells. However, the two-target combinations of (UMPS, HMGCR) and (CAD, HMGCR) could achieve cell mortality and improve metabolic deviation grades over 22% without reducing the cell viability grade.


Introduction
Head and neck cancers (HNC) are among the most common types of cancer worldwide [1,2]. HNC is characterized by a high locoregional recurrence rate and includes cancers of the oral cavity, oropharynx, hypopharynx and larynx. People who smoke, abuse alcohol or chew betel nuts are at greater risk of developing HNC than people who do not have these habits [3]. The incidence of oral cancer (excluding lip cancer) is high in South and Southeast Asia, Western and Eastern Europe, Latin America, and the Caribbean and Pacific regions. Taiwan ( part of East Asia) has one of the world's highest incidence rates of oral cancers [3]. Most HNC cases are squamous cell carcinoma. Recent advances in molecular biology have enabled the documentation of substantial genetic differences between head and neck squamous cell carcinoma (HNSCC) cells and normal cells, potentially enabling the development of new therapeutics and chemoprevention methods [4].
The efficient identification of drug targets and their side effects is a key challenge in drug discovery and development [5][6][7][8][9][10][11][12][13][14][15]. Drug side effects cause substantial clinical and economic burdens. The existing computational methods used to predict the side effects of drugs assume that similar drugs have comparable properties in terms of chemical and biological characteristics, such as their structures and targets. However, a literature review revealed that few studies have attempted to predict side effects of each candidate target in the early stage of drug discovery processes. In recent years, applications of computational systems biology have enabled the modelling and analysis of metabolic pathways, regulatory networks and signal transduction networks to better understand cellular behaviour [5,6]. Tumour initiation and progression require the metabolic reprogramming of cancer cells. Cancer cells autonomously alter their metabolic flux through various metabolic pathways; such pathways include aerobic glycolysis, reduced oxidative phosphorylation and the increased generation of biosynthetic intermediates that are required for cell growth and proliferation [16][17][18]. Numerous studies have employed metabolic reprogramming to understand cancer-specific metabolic networks in order to facilitate drug discovery for cancer treatment [19][20][21][22][23][24][25][26][27][28][29][30][31][32].
Constraint-based modelling approaches entail the use of context-specific genome-scale metabolic networks to identify anti-cancer targets [19][20][21][22][23][24][25][26][27][28][29][30][32][33][34][35][36][37][38]. However, most computer-aided screening methods apply synthetic lethality as a therapeutic indication and ignore side effects during the identification of drug targets. To the best of our knowledge, few studies have attempted to predict side effects of the identified potential targets in the early stage of drug discovery processes. Metabolic flux distributions can alter when normal cells are perturbed due to dysregulation or exerting drug treatment. Metabolic perturbation can be applied to evaluate the influence of the perturbed cell relative to its normal cell counterpart and derives as a measure of side effect for each identified target. The present study used the metabolic perturbation of normal cells relative to the normal template to evaluate side effects for each identified target. A study introduced a fuzzy multi-objective optimization framework for identifying anti-cancer targets (IACT) in order to identify targets that not only increase the mortality of cancer cells but also minimize side effects due to treatment [39]. In addition, the existence and limitations of the IACT framework were theoretically verified [40]. The IACT framework uses RNA-seq expressions of colon cancer cells and their normal counterparts to reconstruct context-specific genome-scale metabolic networks. However, this framework does not include these RNA-seq expressions in computations for inner optimization problems to obtain uniform flux distributions (UFDs). The present study extended the IACT framework to account for RNA-seq information in inner optimization problems in order to yield uniform flux patterns for treated and perturbed cells; these were then applied to identify anti-cancer targets for treating HNSCC.

Tissue-specific genome-scale metabolic models
We combined a human metabolic network (Recon3D) [41] with RNA-seq expression data obtained from The Cancer Genome Atlas (TCGA) database [42] to reconstruct population-based tissue-specific genomescale metabolic models (GSMMs) for HNSCC and healthy counterpart tissues. RNA-seq data of 44 healthy tissue samples with FPKM-UQ normalized expression levels and 500 cancer tissue samples were downloaded from the TCGA database [42] and normalized through quantile normalization; the mean, confidence interval and coefficient of dispersion were then calculated for each gene. These data were used to evaluate supportive genes and classify them into four groups (high, medium, low and not detected). The four gene groups were used to compute the confidence score for each reaction based on the gene-protein-reaction (GPR) associations in Recon3D [41] and classify all reactions into four groups (high, medium, negative and other confidence). The CORDA algorithm [43] used a cost optimization reaction dependency assessment for each confidence reaction to build concise, but not royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 minimalistic, tissue-specific metabolic models for HNSCC and healthy counterpart cells. We developed an integrated toolbox that comprised the aforementioned model reconstruction procedures and to automatically build stoichiometric and GPR models in a computer language of the general algebraic modelling system (https://www.gams.com/) to discover anti-cancer targets with few side effects.
Genome-scale reconstructions provide a mechanistic link between genotype and phenotype. Stoichiometric models represent the material balance of metabolites in cells. GPR associations are typically implemented as Boolean rules, which link metabolic reactions in the stoichiometric models to gene-encoded enzymes in cells. A GPR model generally comprises one-to-one enzymes that catalyse a reaction, and promiscuous enzymes that modulate multiple reactions. Therefore, a reaction can be catalysed by one-to-one enzymes or by isozymes. A simple enzyme is encoded by a single gene, whereas a complex enzyme is encoded by at least two genes. In this study, a simplified network (figure 1) was used to describe the GPR associations. In this network, reactions r 2 and r 3 are catalysed by enzyme E 2 or E 3 (figure 1a); thus, the regulation of both enzymes is identical. We deleted the redundant enzymes in the GSMMs to represent the reduced GPR association (figure 1b), thereby avoiding surplus computations in the evolutionary optimization procedures.

Anti-cancer target discovery problem
During discovery and development, drugs generally fail in the clinical trials for two reasons: they are less efficacious or are demonstrated to have serious side effects during initial trials. This study presents an anti-cancer target discovery (ACTD) platform for identifying potential candidate targets with few side effects. The platform integrates RNA-seq expressions for cancer cells and their healthy counterparts into the IACT framework [39,40] to yield uniform flux patterns for treated and perturbed cells. We reaction gene association enzyme association reduced association (c) Boolean rules of GPR models. The enzyme E 3 is a redundant enzyme that is identical to E 2 and catalyses the same reactions: r 2 and r 3 . The reaction r 3 is catalysed by the isozymes E 2 or E 4 . E 1 is encoded by G 1 to regulate r 1 , and the promiscuous enzymes E 2 and E 4 are encoded by G 2 and a complex of G 4 and G 5 , respectively.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 formulated this extended platform as a fuzzy hierarchical optimization problem. We considered four fuzzy goals in the outer optimization problem in the ACTD platform; we established these goals to account for cancer (CA) cell mortality, healthy (HT) cell viability, perturbation effects and side effects. We directly used fluxes and metabolite flow rates to evaluate the fuzzy objectives. This approach differs from that of the IACT framework, which uses logarithmic fold changes that may have numerical inaccuracies if a value is close to zero. Moreover, in the ACTD platform, the goal targeting CA cell mortality is equivalent to the minimization of not only the cell growth rate but also ATP production rate. The four goals are described as follows: The first goal is to measure the mortality of CA cells for treatment: Fuzzy minimization of the cell growth rate and ATP production rate of treated CA cells: The second goal is to evaluate cell viability of HT cells due to perturbation: Fuzzy minimization of the cell growth rate of perturbed HT cells: g min d,z v PB biomass % 0 Fuzzy maximization of the ATP production rate of pertubed HT cells: The third goal is to measure flux pattern of perturbed HT cells relative to the CA template: Fuzzy dissimilarity of fluxes and metabolite flow rates of perturbed HT cells relative to the CA template: The fourth goal is to measure flux pattern of perturbed HT cells relative to the HT template: Here g min and g max denote fuzzy minimization and maximization, respectively. Moreover, g dissimilarity denotes fuzzy dissimilarity and is used to evaluate the disparity between the fluxes v PB j and metabolite flow rates r PB m for perturbed HT cells relative to the CA template. A significant disparity indicates that the flux distributions and metabolite flow rates of the perturbed HT cells differ considerably from those of the CA template, implying that the perturbation of the HT cells may not lead to tumorigenesis due to treatment. In addition, g similarity denotes fuzzy similarity and is equivalent to a fuzzy equal operator [44]; it is used to evaluate metabolic deviations between the fluxes and metabolite flow rates of the perturbed HT cells relative to the HT template. Smaller metabolic deviations imply lower side effects for a treatment.
The pioneering fuzzy optimization is introduced by Bellman & Zadeh [45] and reveals that fuzzy objectives and/or the constraints constitute classes of alternatives whose boundaries are not sharply defined, and discriminates from a traditional (crisp) optimization. Accordingly, fuzzy minimization (maximization) indicates that the minimum (maximum) objective is less (higher) than its aspiration as possible, and not substantially to fulfil the optimal result completely. Fuzzy similarity is different to a least square error estimation used in a crisp optimization to account for how closeness of the predicted values to the standard level. Furthermore, a crisp optimization method is unable to evaluate a dissimilarity criterion. The present study introduces fuzzy dissimilarity to measure how discrepancy between the predicted values and their standard levels.
Regarding the inner optimization problem, we used constraint-based models of CA and HT cells to calculate the metabolic distributions of treated CA cells and perturbed HT cells, respectively. We reconstructed the tissue-specific GSMMs by using the RNA-seq expression levels of the CA and HT cells to classify metabolic reactions into four groups: high, medium, negative and other confidence, respectively. The weighting factors for the four reaction classes were used to solve the UFD problem royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 to derive a uniform distribution. The inner optimization problem is described as follows: Perturbed HT model: FBA problem: UFD problem: , ð2:2Þ where v f/b is the forward-backward flux vector of reactions; N CA and N HT are stoichiometric matrices for tissue-specific CA and HT cells, respectively, and are reconstructed from those presented in the previous section; v LB f=b,j and v UB f=b,j are the positive lower bound (LB) and upper bound (UB) of the jth forwardbackward flux, respectively; and Ω TR is the set of reactions for CA and HT models. The weighting factors c CA k and c HT k depended on the classified groups and can be assigned as follows: For a high confidence reaction, we can set the smallest weighting factor to obtain a higher flux value from the UFD problem. v LB,TR f=b,i and v UB,TR f=b,i denote the LB and UB of the regulated forward-backward fluxes depended on gene-or metabolite-centric approach for activation. The regulation bounds for the gene-centric approach can be expressed as follows: where v basal f,i and v basal b,i are the basal value of the ith forward-backward flux obtained from XA and HT templates; Ω IZ is the set of reactions regulated by isozymes determined using the GPR associations, and δ is the modulation parameter determined by a nested hybrid differential evolution (NHDE) algorithm [39,40,[46][47][48]. A reaction catalysed by isozymes remains around its basal level; thus, we set the flux ratio ε to 0.03 in this study to restrict the flux value. Metabolite-centric regulators modulate the synthesis reactions of the active metabolites. The LBs and UBs of modulated reactions for the ith active metabolite are restricted as follows: Regulated bounds for the z i th active metabolite: Up-regulation: where N ij is the stoichiometric coefficient of the ith metabolite and the jth reaction.

Maximizing decision-making problem
The ACTD problem expressed in equations (2.1) and (2.2) can be transformed into a maximizing decision-making (MDM) problem by using fuzzy set theory to derive Pareto solutions (figure 2a). Fuzzy minimization and maximization objectives can be assigned membership grades according to their corresponding one-side line membership functions as follows: where FV represents the fluxes and metabolite flow rates for the treated CA and perturbed HT cells, which are computed from equation (2.2); moreover, the LBs and UBs can be provided by a user in advance. These bounds can be obtained from clinical experimental data (if available). Fuzzy dissimilarity is a complement of fuzzy similarity; therefore, the relationship of both membership functions is convertible (figure 2b). A fuzzy similarity grade can be assigned using a two-sided line membership function as follows: Left-hand side membership function: Right-hand side membership function: royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 where ST represents the standard level of CA or HT cells and can be obtained from the corresponding templates of the cells. Therefore, the metabolic deviation grades can be derived from the two-sided line membership function as follows: The metabolic deviation grade ranges between 0 and 1. A value of 0 indicates that an objective has not been met at all, whereas a value of 1 indicates a completely fulfilled objective.
In the ACTD problem, the first goal is to consider the mortality of the treated CA cells; the membership grades h TR biomass and h TR ATP can be derived using equation (2.6). We define a mean-min operator to evaluate the cell mortality grade h TR as follows:  2) and then used to compute the metabolite flow rate r m of the mth metabolite as follows: where Ω c denotes the set of species located in different compartments and N ij denotes stoichiometric coefficients. Both fluxes and metabolite flow rates can be used to compute two-sided membership grades, as indicated in equation ( Þ for metabolic alterations for perturbed cells from the CA and HT templates, respectively. Accordingly, the overall metabolic deviation grade can be derived as follows: The metabolic deviation grade ranges between 0 and 1. A higher metabolic deviation grade implies a smaller degree of metabolic perturbation due to treatment, suggesting that the treatment has fewer side effects. The ACTD problem in equations (2.1) and (2.2) can be transformed into an MDM problem through the aforementioned membership functions. The optimality and limitation of the transformation have been proved in Wang et al. [40]. The MDM problem is expressed as follows: subject to inner optimization problems 1: FBA and UFD problems for treated CA cells 2: FBA and UFD problems for perturbed HT cells;  Figure 3. The comparison of metabolic network data between the HT and CA models was reconstructed through the CORDA algorithm. Statistics of the reconstructed metabolic models for CA and HT cells are also indicated. The numbers listed in the overlapping regions denote the number of identical species, reactions, genes, enzymes and feasible enzymes in the HT and CA models. Feasible enzymes were separated from redundant enzymes through GPR association reduction. where the decision objective η D is a hierarchical criterion. This criterion states that the mortality grade η TR for treated CA cells is the first priority in the decision-making problem; moreover, if either the cell viability or metabolic deviation grade is less than the mortality grade, then one of the lowest grades in the set of {η TR , η CV , η MD } is the second priority for decision-making. The decision-making problem in equation (2.13) is a mixed-integer optimization problem that involves linear and quadratic programming problems in its inner loop. It is a high-dimensional NP-hard problem that cannot be solved by available commercial software [49,50]. We employed the NHDE algorithm to solve the MDM problem. The NHDE algorithm is a parallel direct search procedure, which is an extended version based on hybrid differential evolution [51]. The detailed computational procedures are provided in the electronic supplementary material [52].

Results and discussion
We reconstructed tissue-specific GSMMs for CA and HT cells, that are provided in the electronic supplementary material [52], and figure 3 illustrates the corresponding numbers of species, reactions, genes and encoded enzymes. As indicated by the overlapping regions in figure 3, both models shared numerous similarities in terms of species, reactions, genes and enzymes. Feasible enzymes were separated from redundant enzymes through the reduced GPR associations and were used as candidate variables in equation (2.13) to solve the MDM problem, thus avoiding trivial optimization searches. The iMAT algorithm [53] was also used to reconstruct tissue-specific GSMMs for CA and HT cells, and the corresponding numbers of species, reactions, genes and encoded enzymes are listed in the electronic supplementary material [52] for comparison. The CORDA algorithm [43] is not to keep the reconstructed model as concise as possible. As a result, these numbers are greater than those reconstructed by the iMAT algorithm.

Identifying anti-cancer target genes
We first applied the NHDE algorithm [39,40,[46][47][48] to solve the MDM problem in order to identify optimal anti-cancer target genes; we determined a set of 21 one-target genes with the highest hierarchical fitness among 1104 feasible candidates (table 2). GSMMs for CA and HT cells reconstructed by the iMAT algorithm were also used for the inner optimization problem in equation (2.2) to solve the MDM problem. The optimal anti-cancer target genes were identical to those of the reconstructed models using the CORDA algorithm (see the electronic supplementary material [52]). The cell mortality and cell viability grades obtained using the models reconstructed by CORDA were less than those of the models reconstructed by iMAT, but achieved higher metabolic deviation grade. A survey of the Cancer Dependency Map (DepMap, https://depmap.org/portal/) revealed that most of the identified target genes were compatible to HNC cell lines from DepMap and that these genes could engender a high percentage of cell death when knocked out, except for SLC7A6, SPTLC3 and SLC2A13. The computation results (table 2) revealed that targeting all identified genes for treatment could lead to complete cell death; however, the ATP production rates differed between the targets. Thus, targeting most of these genes for treatment, except for DHFR, SLC2A13, SLC5A3 and CEPT1, could lead to cell mortality grade of 0.668. Similarly, targeting most of the identified genes for the treatment of CA cells and perturbed HT cells could lead to cell viability grades of 0.751. Among the identified one-target genes (table 2), we observed that RIPA knockout yielded the highest metabolic deviation grade (0.682). The metabolic deviation grade is used to evaluate metabolic perturbation for fluxes and metabolite flow rates from the HT and CA templates; thus, a higher metabolic deviation grade implies a smaller metabolic perturbation due to treatment; that is, the treatment would have fewer side effects.
The protein-protein interaction (PPI) networks of the identified 21 one-target genes were investigated using STRING (https://string-db.org/), as presented in figure 4. Markov clustering in STRING was used to classify the PPI network into five classes, which are presented as colour nodes in figure 4. The first class (red nodes) comprised 12 enzymes involved in purine and pyrimidine metabolism and in the pentose phosphate pathway. The enzymes in the second to fifth classes were determined to be royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 involved in the glycerophospholipid biosynthetic pathway, sphingolipid metabolism, NRF2 pathway and basigin interactions. We also performed a brute-force enumerative search to evaluate 1104 feasible enzymes individually to validate the computational results presented in table 2. However, targeting genes such as HMGCR and FDFT1 involved in the statin-induced mevalonate pathway could not eliminate CA cell growth. These results are not consistent with those derived using the IACT framework for IACTs for colon cancer treatment [35].
The ACTD platform can be used to identify multi-target combinations for treating HNSCC. An enumerative search is a time-consuming process because it requires the use of more than 500 000 combinations. Accordingly, we applied the NHDE algorithm by using two groups of candidates to identify multiple targets. The first candidate group comprised the 21 targets earlier (table 2), and the second group comprised other candidate genes from the feasible enzymes. This strategy substantially reduced the computational time and reduced the search space to only 22 743 (21 × 1083) possible combinations between the two candidate groups. We performed a series of computations to obtain a set of two-target combinations. Numerous two-target combinations were identified, and table 3 presents the optimal two-target combinations with cell viability grades greater than 0.75; other Table 3. Top 30 two-target combinations obtained using the ACTD framework. TR and MD denote the cell mortality grade for treated CA cells and the overall metabolic deviation grade for perturbed HT cells to the CA and HT templates, respectively. 'N/D' represents the cell death number (N) divided by the total number of head and neck CA cells (D) used for the test from DepMap. 'no. drugs' denotes the number of drugs retrieved from DrugBank that modulate each gene. UQCR: commonly referred to as complex III. Cytochrome is a heme-containing subunit of the cytochrome b-c1 complex. COX: Complex of cytochrome c oxidase.  (2,9) royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 combinations are listed in the electronic supplementary material [52]. We observed that the combination of phosphoglycerate kinase 1 (PGK1) with each target enzyme (CAD, PAICS and DHODH) in the pyrimidine metabolic pathway improved the cell viability grade by more than 22% (η CV > 0.92) but reduced the metabolic deviation grade by approximately 6% (η MD < 0.526) compared with the use of the single-target counterparts of these enzymes (table 2). Among the targets, ATIC was determined to royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 participate in purine metabolism, and its combination with glutathione reductase (GSR) could increase the cell viability grade by more than 26% and the metabolic deviation grade by 23% (η MD = 0.636).
The cytochrome complex (abbreviated as CYC-COX) contains multiple ubiquinol-cytochrome c oxidoreductase and cytochrome c oxidase subunits. It is part of the mitochondrial respiratory chain and plays a key role in oxidative phosphorylation. We could not eliminate CA cell growth through CYC-COX downregulation only. However, the combination of CYC-COX with any of the 21 singletarget enzymes (table 1) could achieve cell viability and metabolic deviation grades of more than 0.86 and 0.61, respectively. This metabolic deviation grade was noted to be higher than those derived for the PGK1 combinations, indicating that CYC-COX combination treatments would have fewer side effects. Additionally, HMGCR knockdown could not block CA cell growth. However, the two-target combinations of (UMPS, HMGCR) and (CAD, HMGCR) could increase the metabolic deviation grade by more than 22% without reducing the cell viability grade.
We composed a set of the perturbed metabolite flow rates derived from the identified targets to perform factor analysis to investigate their flux patterns. The set comprised 2380 metabolite flow rates and 88 cases, including the CA template, HT template and perturbed cases (see the electronic supplementary material [52]). We conducted an iterated principal factor analysis by using the orthogonal varimax rotation method in SAS software (https://www.sas.com/) to analyse the data; the analysis yielded three factors and their corresponding loadings (see electronic supplementary material [52]). 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. Our analysis revealed that the first-, second-and third-factor loadings (figure 5) were 55.35%, 35.15% and 9.50%, respectively. For each factor, we collected target variables with a loading of more than 0.8 and classified them into three groups (figure 5). The first group comprised the HT template and 41 one-target and two-target combinations with a first factor loading of greater than 0.8. An investigation using GeneCards (https://geneanalytics.genecards.org/) revealed that these target genes mainly participate in the following metabolic pathways: lipid and lipoprotein metabolism, purine and pyrimidine metabolism, and glycosaminoglycan metabolism. The second group comprised 20 two- Table 4. Membership grades for cancer cell mortality (η TR ), cell viability grade (η CV ) and metabolic deviation (η MD ) for the top 10 antimetabolites of one-and two-target combinations. The symbol ↑ indicates upregulation. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 target combinations; specifically, they comprised combinations of the CYC-COX complex with any one target involved in lipid and lipoprotein metabolism, purine and pyrimidine metabolism, or glycosaminoglycan metabolism. The third group comprised two transporters, namely SLC2A13 and SLC5A3. Targets with a factor loading of less than 0.8 were classified into the 'other' group; this group included the CA template and 23 targets. In addition, the flux distributions for all treatments, CA templates and HT templates are presented as a heatmap in the electronic supplementary material [52] to illustrate differences between the treatments. We observed that the flux distributions for each treatment were similar within groups but differed slightly between groups.

Identifying antimetabolites
A metabolite-centric approach in the ACTD framework was also used to discover one-and two-target antimetabolites for treating HNSCC. Table 4 presents some of the optimal targets, and the remaining targets are listed in the electronic supplementary material [52]. The one-target antimetabolites orotate (orot), orotidine 5'-phosphate (orot5p), uridine-5'-monophosphate (ump), N-carbamoyl-L-aspartate (cbasp) and carbamoyl phosphate (cbp) are precursors of nucleotides participating in the pyrimidine biosynthetic pathway of DNA synthesis. The inhibition of these metabolites could lead to cell mortality, cell viability and metabolic deviation grades greater than 0.668, 0.751 and 0.537, respectively. Figure 6 illustrates the synthesis of orot from (S)-dihydroorotate (dhor_S) and quinone  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220633 (q10), which is catalysed by the enzyme DHODH ( figure 6a). Moreover, the enzyme UMPS catalyses two reactions that convert orot to orot5p and subsequently to ump (figure 6a). The enzyme CAD catalyses three reactions (figure 6b,c) that convert dhor_S to cbasp and then to cbp. Moreover, CAD catalyses the synthesis of L-glutamic acid (glu_L) from L-glutamine (gln_L). The identified antimetabolite treatments were noted to be nearly identical to those identified using the gene-centric approaches (tables 2 and 4). These results indicate that the two-target combinations could improve the metabolic deviation grade by more than 13% compared with one-target treatments; however, increases in the cell viability grade were smaller.

Conclusion
Computer-aided drug discovery technology has developed rapidly over recent decades. Computational approaches can help screen potential candidate targets to reduce the time and cost of drug development. Current constraint-based modelling methods can discover potential candidate targets for treatments inhibiting CA cell metabolism. Most computer-aided methods use synthetic lethality as a therapeutic criterion for identifying drug targets; however, they do not attempt to predict side effects of each candidate target in the early stage of drug discovery processes. This study developed an ACTD platform and considered four fuzzy goals for treatment: high cell mortality for CA cells, high cell viability for HT cells, minimal treatment effects on HT cells and low side effects for HT cells. Moreover, RNA-seq expression data were applied to reconstruct tissue-specific GSMMs for CA and HT cells and to improve flux balance models for inner optimization problems to obtain uniform distributions.
We applied the ACTD platform to identify one-target and two-target combinations of anti-cancer enzymes and antimetabolites for HNSCC treatment. For the gene-centric approach, the one-and twotarget combinations mostly comprised targets participating in purine and pyrimidine metabolism, the pentose phosphate pathway, the glycerophospholipid biosynthetic pathway, sphingolipid metabolism, the NRF2 pathway and basigin interactions. By using DrugBank (https://go.drugbank.com/), we found that many of the identified targets could be modulated by approved drugs; thus, these drugs are potential candidates for repurposing against HNSCC. For the metabolite approach, we identified antimetabolite that was nearly identical to those obtained using the gene-centric approach. HMGCR knockdown could not block HNSCC growth. However, the two-target combinations of (UMPS, HMGCR) and (CAD, HMGCR) could improve the metabolic deviation grade by more than 22% without reducing the cell viability grade.  [52]. The data are provided in the electronic supplementary material [54].