Identification of Anticancer Target Combinations to Treat Pancreatic Cancer and Its Associated Cachexia Using Constraint-Based Modeling

1. Introduction

Cancer is fundamentally characterized by profound alterations in cellular physiology, among which metabolic dysregulation plays a central role. Tumor cells undergo extensive metabolic rewiring, driven by the activation of oncogenes, inactivation of tumor suppressor genes, and adaptive modifications in intracellular signaling pathways [1,2,3,4]. This reprogramming supports the heightened energetic and biosynthetic demands of malignant proliferation and facilitates survival under dynamic and often hostile microenvironmental conditions, such as hypoxia [3]. The resulting dependence of specific metabolic pathways introduces selective vulnerabilities that can be strategically targeted for therapeutic intervention [5].

Beyond the primary tumor, advanced cancer frequently gives rise to systemic complications such as cachexia—a multifactorial metabolic disorder marked by the progressive loss of skeletal muscle mass, depletion of adipose tissue, and widespread metabolic dysfunction. Cachexia is particularly prevalent in pancreatic cancer (PC) and significantly compromises patient outcomes by reducing treatment tolerance and overall survival [6]. Among the key molecular mediators of cancer cachexia are inflammatory proteins such as Growth Differentiation Factor-15 (GDF-15) and interleukin-6 (IL-6). GDF-15 contributes to muscle atrophy through both extracellular and intracellular mechanisms: it is secreted by tumor cells—often via exosomes—and actives apoptotic signaling cascades within myocytes [7,8]. IL-6, a central cytokine in inflammatory responses, is strongly associated with systemic hyperinflammation in cachexia, and its inhibition has shown therapeutic promise [10,11]. Studies have further demonstrated that antibody-mediated inhibition of GDF-15 not only suppresses excessive lipid oxidation but also aids in maintaining body weight [8]. On a molecular level, GDF-15 exacerbates muscle wasting by activating the Bcl-2/caspase-3 apoptoic pathway and upregulating Atrogin-1, a muscle-specific E3 ubiquitin ligase involved in proteasomal degradation of muscle proteins [13,14].

Identifying therapeutic targets that can effectively suppress tumor progression while concurrently alleviating cancer-associated cachexia and minimizing off-target effects remains a major clinical and translational challenge. Cell-specific genome-scale metabolic models (GSMMs) have emerged as powerful computational tools for elucidating the metabolic reprogramming that differentiates malignant cells from their healthy counterparts. By integrating high-throughput omics data—such as transcriptomics, proteomics, and metabolomics—these models enable simulation and comparative analysis of metabolic flux distributions under diverse physiological and pathological conditions. This systems-level approach allows for the systematic identification of dysregulated metabolic pathways and enzyme-catalyzed reactions that may represent viable intervention points for anticancer therapy.

Despite substantial progress, most previous studies have focused primarily on the intrinsic metabolic alterations of cancer cells [15,16,17], often overlooking the broader metabolic consequences of therapeutic interventions and the development of cancer-associated cachexia. In this study, we employed cell-specific GSMMs to address three interrelated dimensions of cancer management: the identification of anticancer targets, the discovery of anti-cachexia targets, and the prediction of drug-induced side effects. To identify anticancer targets, a cell mortality index was incorporated to pinpoint metabolic pathways uniquely dysregulated in tumor cells. Potential side effects were anticipated by assessing metabolic perturbation induced in healthy cells following target-specific intervention. Additionally, anti-cachexia targets were identified by characterizing the degradation reactions of key structural muscle proteins—myosin, actin, and titin—which serve as surrogate biomarkers for muscle wasting. This integrative strategy provides a comprehensive framework for understanding the systemic impact of cancer therapies and supports the development of metabolic interventions that concurrently suppress tumor growth and ameliorate cachexia.

2. Results and Discussion

2.1. Cell-Specific Genome-Scale Metabolic Models

The human genome-scale metabolic network (GSMN) Recon3D [18] was employed as a foundation to reconstruct cell-specific genome-scale metabolic models (GSMMs) for identifying potential biomarkers and enzyme targets that inhibit cancer cell proliferation and mitigate muscle degradation. However, Recon3D lacks essential components, including pathways for the synthesis, degradation, and recycling of proinflammatory cytokines (e.g., TNF-α, IFN-γ, IL-6) and key structural proteins such as myosin, actin, and titin. To bridge this gap, a lumped modeling strategy was adopted, wherein protein sequences were used to define the stoichiometry of the missing reactions. These reactions were subsequently integrated into Rencon3D to form an extended GSMN, as illustrated in . This extended network served as a generic template, which was then customized using RNA-sequence expression data from cancerous and cachectic tissues to reconstruct cell-specific GSMMs for downstream metabolic analysis and therapeutic design.

RNA-Seq expression data for pancreatic duct adenocarcinoma cancer (PDAC) cells were obtained from Gene Expression Omnibus (GEO) database of the National Center for Biotechnology Information (NCBI) under accession number GSE183795, based on the study by Yang et al. [19]. This dataset comprises expression profiles from 139 tumor samples and 102 matched adjacent non-tumor tissues, which were used to reconstruct cell-specific genome-scale metabolic models (GSMMs) for both PDAC and its healthy tissue counterparts (HT), with a particular focus on proinflammatory cytokine production. Additionally, gene expression data for cachectic rectus abdominis muscle from PDAC patients were obtained from the study by Narasimhan, et al. [20], including 23 cachectic PDAC samples (PDAC-CX) and 11 non-cancer controls. Following reconstruction protocols described by Cheng et al. [15] and Wang et al. [17], these datasets were used to develop cell-specific GSMMs to characterize muscle degradation involving key structural proteins such as myosin, actin and titin. The four reconstructed cell-specific GSMMs are provided in Supplementary File S1 through S4.

Figure 1. Integration of Protein Synthesis, Degradation, and Recycling into the Human Genome-Scale Metabolic Network (Recon3D) to Generate an Extended Network.

Figure 1. Integration of Protein Synthesis, Degradation, and Recycling into the Human Genome-Scale Metabolic Network (Recon3D) to Generate an Extended Network.

summarizes the number of metabolites (species) and reactions in these reconstructed models. In this context, a ‘species’ refers to a metabolite assigned to one of nine cellular compartments represented in the metabolic model. The four models had 2863 species, and 3741 reactions shared in common, as shown in the overlapping region in . The PDAC model comprised 3669 species, and 5669 reactions, while the PDAC-CX model consisted of 3419 species, and 5543 reactions. Of these, 52 species and 162 reactions were unique to the PDAC model, whereas 87 species and 465 reactions were specific to the PDAC-CX model. The PDAC network retains approximately 62% of the species and 52.9% of the reactions present in the extended Recon3D model. In comparison, PDAC-CX network encompasses about 57.8% of the species and 51.8% of the reactions.

Figure 2. Summary Statistics of Reconstructed Metabolic Models for Pancreatic Duct Adenocarcinoma Cancer (PDAC), Cachectic PDAC (PDAC-CX), and Their Corresponding Healthy Counterparts (PDAC-HT and PDAC-CX-HT).

Figure 2. Summary Statistics of Reconstructed Metabolic Models for Pancreatic Duct Adenocarcinoma Cancer (PDAC), Cachectic PDAC (PDAC-CX), and Their Corresponding Healthy Counterparts (PDAC-HT and PDAC-CX-HT).

2.2. Identified Potential Biomarkers

The identification of potential biomarkers is crucial for inhibiting cancer cell growth and mitigating muscle degradation, particularly in conditions like pancreatic cancer-associated cachexia. To achieve this, a computational framework incorporating Parsimonious Metabolite Flow Variability Analysis (pMFVA) combined with differential expression analysis was employed to identify statistically significant changes in metabolite flow rates (p-value < 0.05), as illustrated in . This analytical approach was applied across five distinct nutritional media to assess the impact of nutrient availability on biomarker detection. A key aspect of our approach was to identify “medium-independent biomarkers” by focusing on metabolites that exhibited consistent directional changes, showing either a complete increase or a complete decrease across a panel of five distinct nutrient media used in simulations. This rigorous selection process helps to overcome nutrient dependency, ensuring the robustness and broader applicability of the identified biomarkers. summarizes the number of detected biomarkers under each condition. In the PDAC model, 524, 507, 483, 536, and 433 biomarkers were identified across the five media, DMEM, HAM, HPLM, RPMI, and VMH, respectively. In contrast, the PDAC-CX model yielded 428, 474, 512, 496, and 504 biomarkers across the same media. These differences highlight the sensitivity of metabolic flux distributions to nutrient availability and suggest the potential for partial or complete overlap in metabolite flow between diseased and healthy states.

Table 1. Number of Identified Biomarkers in PDAC and PDAC-CX Based on Metabolite Flow Rates Across Five Nutritional Media. Abbreviations: CI―complete increase; PI―partial increase, II―inclusive increase; ID―inclusive decrease; PD―partial decrease; CD―complete decrease (as defined in ).

Medium	Type	CI	PI	II	ID	PD	CD	Total
DMEM	PDAC	97	240	12	26	4	145	524
	PDAC-CX	80	134	16	58	24	116	428
HAM	PDAC	132	212	12	20	6	125	507
	PDAC-CX	78	130	20	72	20	154	474
HPLM	PDAC	125	211	17	32	7	91	483
	PDAC-CX	88	175	17	112	26	94	512
RPMI	PDAC	134	222	11	15	6	148	536
	PDAC-CX	95	146	20	89	21	125	496
VMH	PDAC	63	203	10	31	8	118	433
	PDAC-CX	140	171	25	43	29	96	504

Table 1. Number of Identified Biomarkers in PDAC and PDAC-CX Based on Metabolite Flow Rates Across Five Nutritional Media. Abbreviations: CI―complete increase; PI―partial increase, II―inclusive increase; ID―inclusive decrease; PD―partial decrease; CD―complete decrease (as defined in ).

Medium	Type	CI	PI	II	ID	PD	CD	Total
DMEM	PDAC	97	240	12	26	4	145	524
	PDAC-CX	80	134	16	58	24	116	428
HAM	PDAC	132	212	12	20	6	125	507
	PDAC-CX	78	130	20	72	20	154	474
HPLM	PDAC	125	211	17	32	7	91	483
	PDAC-CX	88	175	17	112	26	94	512
RPMI	PDAC	134	222	11	15	6	148	536
	PDAC-CX	95	146	20	89	21	125	496
VMH	PDAC	63	203	10	31	8	118	433
	PDAC-CX	140	171	25	43	29	96	504

Figure 3. Workflow for Biomarker Identification Using Differential Expressions Analysis and Parsimonious Metabolite Flow Variability Analysis.

Figure 3. Workflow for Biomarker Identification Using Differential Expressions Analysis and Parsimonious Metabolite Flow Variability Analysis.

The findings demonstrate that the choice of nutritional medium significantly influences biomarker detection in the PDAC model. To refine biomarker selection, we focused on metabolites that exhibited consistent directional changes, showing either complete increase or complete decrease across all five media, as illustrated in (K). This strategy allowed for the identification of medium-independent biomarkers. presents these medium-independent biomarkers, highlighting those with flow rates that consistently increase or decrease across all media.

Figure 4. Nutritional Medium-Independent PDAC Biomarkers Identified from Intracellular Metabolite Flow Rates.

Figure 4. Nutritional Medium-Independent PDAC Biomarkers Identified from Intracellular Metabolite Flow Rates.

For example, proinflammatory cytokines such as IL-6 and TNF-α, key components of the proinflammatory biosynthetic pathway, consistently show elevated flow rates in the cancer state, reflecting their established roles in tumor-associated inflammation and pancreatic cancer progression [21,22]. Metabolites involved in glycolysis and redox balance, including pyruvate, NADH, NAD⁺, and lactate, also exhibit consistently increased fluxes. These changes align with the Warburg effect and alter redox regulation commonly observed in pancreatic cancer cells [23,24]. Ubiquinol-10, located in the inner mitochondrial membrane, facilitates electron transport and ATP production through oxidative phosphorylation. Correspondingly, TCA cycle intermediates such as oxalate and oxaloacetate show sustained increases in flux, supporting enhanced nucleotide and lipid biosynthesis. Additionally, oxaloacetate production via GOT1 contributes to maintaining NAD⁺ regeneration, which is critical for sustaining metabolic and redox homeostasis in cancer cells [21].

Conversely, twelve metabolites exhibited consistent decreases in PDAC, indicative of disrupted biosynthesis, compromised NADPH-dependent redox defense. Notably, oxalosuccinate depletion in cancer cells likely results from the diversion of TCA cycle intermediates toward anabolic processes and redox balancing. Fatty acyl-CoA derivatives such as palmitoyl-CoA and trans-hexadec-2-enoyl-CoA also decrease, reflecting altered lipid metabolism. Thymidine levels decline, suggesting perturbations in nucleotide metabolism. Meanwhile, uptake fluxes of phosphate, L-glutamate, and ubiquinol increase consistently, as illustrated in , underscoring their importance in supporting cancer cell metabolism. Additionally, secretion of 2-deoxyadenosine is reduced in cancer cells, which may influence nucleotide salvage pathways [25,26].

Using similar analytical procedures, the PDAC-CX model identified a set of medium-independent biomarkers, as shown in , comprising 28 metabolites with consistently increased fluxes and 14 metabolites with consistent decreases across all nutritional media. Notably, peptides and degradation products derived from titin, myosin and actin were markedly elevated, indicating enhanced muscle protein catabolism. These degradation products were subsequently secreted into the extracellular space, as illustrated in , aligning with the muscle wasting characteristic of cancer cachexia. Additionally, increased uptake fluxes were observed for L-alanine, glycine, L-glutamate, phosphate and hydrogen peroxide. Phosphate is critical for the biosynthesis of nucleic acids and membrane phospholipids, its elevated uptake in proliferative disease state such as cancer supports increased nucleotide synthesis and membrane production, thereby promoting anabolic growth. In contrast, increased uptake of hydrogen peroxide contributes to redox imbalance. While low to moderate levels of hydrogen peroxide act as signaling molecules to promote cellular adaptation and survival, excessive accumulation induces oxidation stress, leading to macromolecular damage and potentially triggering cell death.

Figure 5. Uptake Reactions and Secretion Reactions Identified as Nutritional Medium-Independent PDAC Biomarkers.

Figure 5. Uptake Reactions and Secretion Reactions Identified as Nutritional Medium-Independent PDAC Biomarkers.

Figure 6. Nutritional Medium-Independent PDAC-CX Biomarkers Identified from Intracellular Metabolite Flow Rates.

Figure 6. Nutritional Medium-Independent PDAC-CX Biomarkers Identified from Intracellular Metabolite Flow Rates.

Figure 7. Uptake and Secretion Reactions Identified as Nutritional Medium-Independent PDAC-CX Biomarkers.

Figure 7. Uptake and Secretion Reactions Identified as Nutritional Medium-Independent PDAC-CX Biomarkers.

2.3. Enzyme Targets Predicted Using Constraint-Based Modeling

The Anticancer Target Discovery (ACTD) platform proposed by Wang and Zhang [17], as illustrated in , was employed to identify enzyme targets for treating PDAC and PDAC-CX across five distinct nutritional media. The predicted targets are detailed in Supplementary File S1. In this framework, cell viability (CV) and metabolic deviation (MD) were used as dual decision-making criteria to evaluate the therapeutic efficacy of each enzyme target. For PDAC, five objective functions were simultaneously minimized to determine CV, with membership grade ranging from 0 and 1. The highest priority was assigned to abolishing biomass production, while additional objectives included the reduction of ATP synthesis and proinflammatory cytokine production (IL6, IFN-γ, and TNF-α).

For PDAC-CX, the primary objective was to minimize myosin degradation, followed by reduction in actin and titin degradation―reflecting muscle preservation goals in cachexia management. Each candidate target was further evaluated for its impact on healthy tissue (denoted as PB, perturbed by treatment), by examining alterations in flux distributions. The MD metric―normalized between 0 and 1 using a linear membership function―quantified the similarity between PB and healthy tissue (HT), and the dissimilarity between PB and the corresponding diseased models. This metric served as proxy for predicting the potential side effects associated with each enzyme target.

Using the ACTD platform, twenty-seven single-enzyme targets for PDAC treatment were identified across five distinct nutritional media, as detailed in Supplementary File S5. Several targets resulted in zero cell viability (CV) in certain media, indicating complete treatment failure under those specific conditions, However, other targets were able to inhibit cancer growth. These findings suggest that the efficacy of the identified targets is dependent on the nutritional environment. Notably, none of these single-enzyme targets were effective against the PDAC-CX model. Conversely, thirty single-enzyme targets identified for PDAC-CX treatment failed to produce therapeutic effects in PDAC, highlighting the context-specific nature of these interventions. To address this limitation, the ACTD platform was extended to a combinatorial targeting strategy aimed at identifying enzyme target pairs effective against both PDAC and PDAC-CX. Some of these combinations, also listed in Supplementary File S5, yielded zero CV but remained sensitive to medium conditions. Among them, seven target combinations emerged as robust candidates exhibiting medium-independent efficacy. The average values of cell viability (CV) and metabolic deviation (MD) values for these medium-independent combinations in treating both PDAC and PDAC-CX are summarized in .

Figure 8. Framework for Anticancer Target Discovery Using Constraint-Based Modeling and Fuzzy Multiobjective Optimization. (A-E) Reconstruction of cell-specific GSMMs and gene–protein–reaction associations, as described in . (F) Formulation of a parsimonious multiobjective flux balance analysis (pMOFBA) problem for treated CA cells. (G) Formulation of a pMOFBA for perturbed HT cells. (H) Calculation of metabolic flux distributions in treated CA cells for each anticancer candidate by solving the pMOFBA model. (I) Calculation of metabolic flux distributions in perturbed HT cells for each anticancer candidate via pMOFBA. (J) Derivation of the metabolic template of CA cells from clinical data (if available) or baseline pMOFBA without target inregulation. (K) Derivation of the metabolic template of HT cells from clinical data (if available) or baseline pMOFBA. (L) Transformation of fuzzy multiobjective functions into a fuzzy decision score (η_D) using fuzzy set theory. (M) Evaluation of each anticancer candidate’s fitness based on η_D to guide target selection. (N) Generation of new anticancer candidates using a nested hybrid differential evolution algorithm if decision criteria are unmet, with iteration through Steps (F) to (N). (O) Identification of optimal anticancer targets when the decision criterion is satisfied.

Figure 8. Framework for Anticancer Target Discovery Using Constraint-Based Modeling and Fuzzy Multiobjective Optimization. (A-E) Reconstruction of cell-specific GSMMs and gene–protein–reaction associations, as described in . (F) Formulation of a parsimonious multiobjective flux balance analysis (pMOFBA) problem for treated CA cells. (G) Formulation of a pMOFBA for perturbed HT cells. (H) Calculation of metabolic flux distributions in treated CA cells for each anticancer candidate by solving the pMOFBA model. (I) Calculation of metabolic flux distributions in perturbed HT cells for each anticancer candidate via pMOFBA. (J) Derivation of the metabolic template of CA cells from clinical data (if available) or baseline pMOFBA without target inregulation. (K) Derivation of the metabolic template of HT cells from clinical data (if available) or baseline pMOFBA. (L) Transformation of fuzzy multiobjective functions into a fuzzy decision score (η_D) using fuzzy set theory. (M) Evaluation of each anticancer candidate’s fitness based on η_D to guide target selection. (N) Generation of new anticancer candidates using a nested hybrid differential evolution algorithm if decision criteria are unmet, with iteration through Steps (F) to (N). (O) Identification of optimal anticancer targets when the decision criterion is satisfied.

The results reveal that the identified combinations comprise eight enzyme targets, including five individually encoded enzymes and one enzyme complex―namely the RNF20–RNF40 complex. Among these, four enzymes (SLC29A2, RNF20–RNF40, CRLS1, and SGMS1) are knocked out, while two enzymes (CERK and PIKFYVE) are upregulated to achieve therapeutic efficacy. In the PDAC model, the average cell viability (CV) ranged from 0.66 to 0.98, and the metabolic deviation (MD) from 0.32 to 0.43. In the PDAC-CX model, CV ranged from 0.86 to 1.0 and MD from 0.25 to 0.33. Higher values of CV and MD reflect increased therapeutic efficacy and reduced potential side effects. Collectively, the eight identified target combinations represent effective therapeutic strategies for concurrently treating both PDAC and PDAC-CX cells. Furthermore, the CV and MD metrics provide a quantitative basis for selecting optimal trade-off between efficacy and safety.

SLC29A2 (solute carrier family 29 member 2), also known as ENT2, is a bidirectional transporter that regulates intracellular homeostasis by mediating nucleosides and nucleobase uptake and efflux. Its knockout disrupts nucleotide salvage pathways, impairs cellular energy balance, and inhibits cancer cell proliferation. In the PDAC-CX model, the RNF20–RNF40 complex functions as a key E3 ubiquitin ligase responsible for ubiquitinating protein-L-lysine residues―including those in myosin, actin, and titin. Knockout of this complex may impair lysine and ubiquitin recycling, thereby promoting enhanced protein degradation (). Consequently, dual targeting of SLC29A2 and the RNF20–RNF40 complex represents a promising therapeutic strategy for PDAC and PDAC-CX cells.

Knockout of the SGMS1 (sphingomyelin synthase 1) gene leads to significant metabolic disruption by impairing sphingolipid metabolism, thereby reducing cancer cell viability and inhibiting tumor progression. SGMS1 catalyzes the synthesis of sphingomyelin by transferring a phosphocholine group from phosphatidylcholine to N-acylsphingosine. Similarly, knockout of the CRLS1 (cardiolipin synthase 1) gene impairs mitochondrial lipid metabolism by blocking cardiolipin biosynthesis, which involves the conversion of cytidine-5′-diphosphate-diacylglycerol and phosphatidylglycerol. Both sphingomyelin and cardiolipin are critical lipid components essential for biomass synthesis in the PDAC model. Consequently, knockout of either SGMS1 or CRLS1 abolishes biomass production. Therefore, co-targeting SGMS1 or CRLS1 in combination with the RNF20-RNF40 complex constitutes a promising therapeutic strategy for the simultaneous treatment of PDAC and PDAC-CX cells.

Upregulation of CERK (ceramide kinase) or PIKFYVE (phosphoinositide kinase) alone is insufficient to prevent muscle degradation in the PDAC-CX model across all five nutritional media, as shown in Supplementary File S5. CERK catalyzes the phosphorylation ceramide to generate ceramide-1-phosphate (C1P). Deletion or downregulation of CERK results in ceramide accumulation, promoting pro-catabolic and pro-apoptotic signaling in muscle cells. This effect is particularly detrimental in PDAC-associated cachexia, where inflammation, energy imbalance, and lipid metabolism are already dysregulated. Although CERK upregulation may exert protective effects by shifting sphingolipid balance toward C1P, the current dataset provides no direct evidence supporting this outcome. Nevertheless, previous studies have associated CERK deficiency with muscle atrophy, underscoring its therapeutic potential for preserving muscle mass under cachectic conditions [27,28,29].

Similarly, upregulation of PIKFYVE may modulate muscle metabolism through its involvement in insulin sensitivity and glucose uptake. While most existing studies emphasize the detrimental effects of PIKFYVE deficiency―such as impairs insulin signaling and glucose homeostasis [30]―its upregulation could potentially enhance these metabolic pathways, thereby contributing muscle preservation. However, direct evidence supporting this effect remains limited. Consequently, combined the upregulation of CERK or PIKFYVE with the knockout of SLC29A2, SGMS1 or CRLS1 may constitute a more effective therapeutic strategy for concurrently targeting both PDAC and PDAC-CX cells across all five media conditions, as illustrated in

Figure 9. Average Cell Viability (CV) and Metabolic Deviation (MD) of Medium-Independent Target Combinations for the Treatment of PDAC and PDAC-CX Cells.

Figure 9. Average Cell Viability (CV) and Metabolic Deviation (MD) of Medium-Independent Target Combinations for the Treatment of PDAC and PDAC-CX Cells.

3. Methods

3.1. Extension of the Human Genome-Scale Metabolic Network

illustrates the workflow used to reconstruct cell-specific genome-scale metabolic networks (GSMNs) for PDAC and PDAC-CX cells, based on transcriptomic data. RNA-Seq data for PDAC cells were obtained from Gene Expression Omnibus (GEO) database at the National Center for Biotechnology Information (accession number GSE183795), as reported by Yang et al. [19]. Transcriptomic data for PDAC-CX cells were sourced from the study by Narasimhan et al. [20]. These datasets were processed using established reconstruction protocols described in previous studies [15,16,17] to generate the respective cell-specific GSMNs. However, the standard human genome-scale metabolic network, Recon3D [18], lacks essential pathways for the synthesis, degradation, and recycling of key proinflammatory cytokines (e.g., TNF-α, IFN-γ, IL-6) and structural proteins (e.g., myosin, actin, and titin). To address this limitation, a lumped-reaction approach was applied, wherein protein sequences were used to define and construct the missing reactions. These reactions were subsequently incorporated into Rencon3D to create an extended GSMN, as depicted in .

Each proinflammatory cytokine or skeletal muscle protein is modeled as a polymerization reaction based on its amino acid sequence, retrieved from the UniProt database (https://www.uniprot.org/). The generalized mass balance for protein polymerization is expressed as:

$${\sum }_{A_j}{S}_{A_j}{A}_j+S_{H_{2}O}{H}_{2}O+S_{ATP}ATP\to S_{Protein}〈{\sum }_j{A}_j〉+S_{ADP}ADP+S_{P_i}{P}_i+S_{H^{+}}{H}^{+}$$

(1)

Here, $〈{\sum }_j{A}_j〉$ denotes the synthesized protein composed of amino acids $A_j$. The reaction involves water $\left(H_{2}O\right)$, adenosine triphosphate $\left(ATP\right)$, adenosine diphosphate $\left(ADP\right)$, inorganic phosphate ${(P}_i)$, and proton $\left(H^{+}\right)$. The stoichiometric coefficients (in mmole per gram of protein) are calculated as follows:

$$S_{A_j}={\displaystyle \frac{{10}^3{M}_{A_j}}{M_{Protein}}}$$

(2)

$$S_{H_{2}O}={\displaystyle \frac{{10}^3\left(k_{ATP}-1\right)\left({\sum }_{A_j}{M}_{A_j}-1\right)}{M_{Protein}}}$$

(3)

$$S_{ATP}=S_{ADP}=S_{P_i}=S_{H_{2}O}={\displaystyle \frac{{10}^3{k}_{ATP}\left({\sum }_{A_j}{M}_{A_j}-1\right)}{M_{Protein}}}$$

(4)

$$S_{ATP}=S_{ADP}=S_{P_i}=S_{H_{2}O}={\displaystyle \frac{{10}^3{k}_{ATP}\left({\sum }_{A_j}{M}_{A_j}-1\right)}{M_{Protein}}}$$

(5)

$$S_{Protein}={\displaystyle \frac{1000}{M_{Protein}}}$$

(6)

In these equations, $M_{protein}$ is the molecular weight of the synthesized protein, $M_{A_j}$ is the molar mass of each amino acid $A_j$, and $k_{ATP}$ denotes the number of ATP molecules required per peptide bond, typically ranging from 4 to 4.3. [31]. Detailed amino acid composition and molar contributions are provided in references [32,33].

Finally, the polymerization and degradation reactions, along with amino acid recycling pathways for each cytokine and structural proteins, were incorporated into the extended Recon3D model. Using this augmented network, cell-specific GSMNs for PDAC and PDAC-CX were reconstructed following the established procedures described in references [15,16,17].

3.2. Parsimonious Multi-Objective Flux Balance Analysis

Parsimonious Multi-Objective Flux Balance Analysis (pMOFBA) improves upon classical flux balance analysis by introducing an additional optimization criterion: the minimization of the total fluxes across the metabolic network. This method is formulated as a two-stage optimization framework that integrates Multi-Objective Flux Balance Analysis (MOFBA) with a Uniform Flux Distribution (UFD) strategy:

$$\begin{array}{cc}{1}^{st}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}:\mathrm{M}\mathrm{O}\mathrm{F}\mathrm{B}\mathrm{A}& 2^{nd}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}:\mathrm{U}\mathrm{F}\mathrm{D}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{m}\\ \left\{\begin{array}{c}\begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\max }{z}_1=v_{biomass}\\ \begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\max }{z}_2=v_{ATP}\\ \underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\max }{z}_{P_k}=v_{P_k},\forall k\in \left\{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{s}\right\}\end{array}\\ \begin{array}{c}\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\\ \mathbf{N}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=0\\ \begin{array}{c}0\le {\mathbf{v}}_f^{LB}\le {\mathbf{v}}_f\le {\mathbf{v}}_f^{UB} \\ 0\le {\mathbf{v}}_b^{LB}\le {\mathbf{v}}_b\le {\mathbf{v}}_b^{UB} \end{array}\end{array}\end{array}\\ GPR association \end{array}\right.& \left\{\begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\min }\sum _{i=1}^n{w}_i\left(v_{f,i}+v_{b,i}\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\\ \begin{array}{c}\mathbf{N}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=0\\ 0\le {\mathbf{v}}_f^{LB}\le {\mathbf{v}}_f\le {\mathbf{v}}_f^{UB} \\ \begin{array}{c}0\le {\mathbf{v}}_b^{LB}\le {\mathbf{v}}_b\le {\mathbf{v}}_b^{UB} \\ \mathrm{G}\mathrm{P}\mathrm{R} \mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ z_1\ge z_1^{*},z_2\le z_2^{*}, z_{P_i}\le z_{P_i}^{*} \end{array}\end{array}\end{array}\right.\end{array}$$

(7)

In this framework, the first stage (MOFBA) determines an optimal flux distribution $\left({\mathbf{v}}_f^{*},{\mathbf{v}}_b^{*}\right)$ that simultaneously maximizes multiple biologically relevant objectives―such as biomass production $\left(z_1\right)$, ATP generation $\left(z_2\right)$, and related-protein synthesis $\left(z_{P_k}\right)$. The second stage (UFD) then minimizes the total weighted flux through the network to enforce a parsimonious (biologically efficient) flux profile, while ensuring that the performance levels obtained in the first stage $\left(z_1^{*},z_2^{*},z_{P_k}^{*}\right)$ are maintained.

Here, ${\mathbf{v}}_f$ and backward ${\mathbf{v}}_b$ represent the forward and backward flux vectors, respectively. The stoichiometric matrix $\mathbf{N}\in {\mathbb{R}}^{\mathit{m}\times \mathit{n}}$ defines the metabolic network, where $m$ is the number of metabolites and $n$ the number of reactions. Fluxes are constrained by lower $\left({\mathbf{v}}_f^{LB} \mathrm{a}\mathrm{n}\mathrm{d} {\mathbf{v}}_b^{LB}\right)$ and upper $\left({\mathbf{v}}_f^{UB} \mathrm{a}\mathrm{n}\mathrm{d} {\mathbf{v}}_b^{UB}\right)$ bounds. Gene-Protein-Reaction (GPR) associations impose regulatory constraints based on gene expression and enzyme availability, enabling integration with transcriptomic data. These constraints are consistently applied to both diseased-specific GSMNs and the corresponding healthy models

In the UFD objective, the weights $w_i$reflect the level of transcriptomic confidence associated with each reaction and are assigned based on a quartile-based classification scheme:

$$w_i=\left\{\begin{array}{c}1/4, i\in \left\{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\right\} \\ 1/2,i\in \left\{\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{u}\mathrm{m}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\right\} \\ \begin{array}{c}3/4,i\in \left\{\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\right\} \\ 1,i\in \left\{\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{o}\mathrm{r}\mathrm{n}\mathrm{o}\mathrm{n}-\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}-\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right\}\end{array}\end{array}\right.$$

(8)

This weighting scheme biases the solution toward flux distributions that prioritize reactions with strong transcriptomic support, thereby promoting biologically meaningful and resource-efficient metabolic activity. Consequently, pMOFAB improves both the interpretability and predictive accuracy of reconstructed metabolic models by aligning flux predictions with gene expression profiles and cellular resource constraints.

3.3. Parsimonious Metabolite Flow Variability Analysis

Parsimonious Metabolite Flow Variability Analysis (pMFVA), which incorporates multi-objective functions, is formulated as follows:

$$\left\{\begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b,\varsigma \in \left[\mathrm{0,1}\right]}{\max /\min }{r}_{m_k}, m_k\in {\mathsf{\Omega }}^m\\ \begin{array}{cc}{1}^{st} \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}:\mathrm{M}\mathrm{O}\mathrm{F}\mathrm{B}\mathrm{A}& 2^{nd} \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}:\mathrm{U}\mathrm{F}\mathrm{D}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{m}\\ \left\{\begin{array}{c}\begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\max }{z}_1=v_{biomass}\\ \begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\max }{z}_2=v_{ATP}\\ \underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\max }{z}_{P_k}=v_{P_k},\forall k\in \left\{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{s}\right\}\end{array}\\ \begin{array}{c}\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\\ \mathbf{N}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=0\\ \begin{array}{c}0\le {\mathbf{v}}_f^{LB}\le {\mathbf{v}}_f\le {\mathbf{v}}_f^{UB} \\ 0\le {\mathbf{v}}_b^{LB}\le {\mathbf{v}}_b\le {\mathbf{v}}_b^{UB} \end{array}\end{array}\end{array}\\ GPR association \end{array}\right.& \left\{\begin{array}{c}\underset{{\mathbf{v}}_f,{\mathbf{v}}_b}{\min }\sum _{i=1}^n{w}_i\left(v_{f,i}+v_{b,i}\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\\ \begin{array}{c}\mathbf{N}\left({\mathbf{v}}_f-{\mathbf{v}}_b\right)=0\\ 0\le {\mathbf{v}}_f^{LB}\le {\mathbf{v}}_f\le {\mathbf{v}}_f^{UB} \\ \begin{array}{c}0\le {\mathbf{v}}_b^{LB}\le {\mathbf{v}}_b\le {\mathbf{v}}_b^{UB} \\ \mathrm{G}\mathrm{P}\mathrm{R} \mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ z_1\ge {\zeta z}_1^{*},z_2\le {\zeta z}_2^{*}, z_{P_i}\le {\zeta z}_{P_i}^{*}\end{array}\end{array}\end{array}\right.\end{array}\end{array}\right.$$

(9)

Here, $\left(z_1^{*},z_2^{*},z_{P_k}^{*}\right)$denote the optimal objective values obtained the MOFBA solution at $\varsigma =1$, representing the fully optimal metabolic state. The parameter $\varsigma \in \left[\mathrm{0,1}\right]$ serves as an outer-level decision variable in pMFVA framework, enabling evaluation under optimal $\left(\varsigma =1\right)$ and suboptimal $\left(\varsigma <1\right)$ conditions. The goal of pMFVA is to determine the maximum and minimum feasible flow rates of key metabolites $r_{m_k}$ across the entire range of $\varsigma$, while also generating a set of flux distributions associated with varying levels of metabolic performance. This facilitates a comprehensive assessment of metabolic flexibility and robustness under different physiological or perturbation scenarios.

After solving the MOFBA problem, the UFD problem is iteratively solved for various values of $\varsigma \in \left[\mathrm{0,1}\right)$, generating a spectrum of metabolite flow rates defined as:

$$r_{m_k}\left(\varsigma \right)={\sum }_{N_{ij>0,j=1}}^n{N}_{ij}{v}_{f.j}^{\mathrm{*}}\left(\varsigma \right)-{\sum }_{N_{ij}<0,j=1}^n{N}_{ij}{v}_{b,j}^{\mathrm{*}}\left(\varsigma \right),m_k\in {\mathsf{\Omega }}^m,\varsigma \in \left[\mathrm{0,1}\right]$$

(10)

The overall variability range for each metabolite flow is then computed as:

$$\left[r_{m_k}^{max}, r_{m_k}^{min} \right]=\left[\underset{\varsigma }{\min }\left\{r_{m_k}\left(\varsigma \right)\right\},\underset{\varsigma }{\max }\left\{r_{m_k}\left(\varsigma \right)\right\}\right]$$

(11)

This framework enables a biologically meaningful and parsimonious characterization of metabolite flow variability. By applying Eqs.(9) and (10), two set of metabolite flow rates―corresponding to diseased and healthy conditions―can be obtained. Comparative analysis of these profiles allows for the identification of differential expressed metabolic flows, which may serve as potential biomarkers, as illustrated in .

3.4. Anticancer Target Discovery Framework

represents the Anticancer Target Discovery (ACTD) framework, as proposed by Wang and Zhang [17]. This framework is formulated as a fuzzy multi-objective hierarchical optimization problem designed to identify potential therapeutic targets for the simultaneous treatment of PDAC and PDAC-CX, as illustrated in . As shown in (A), fuzzy minimization and maximization objectives are used to evaluate treatment-associated cell mortality in diseased cells and cell survival in healthy cells. These objectives are assessed using one-sided linear membership functions, represented by dashed and dot-dashed lines, respectively. The resulting membership values are used to define cell viability, denoted as ${\eta }_{CV}\in \left[\mathrm{0,1}\right]$, which quantifies the therapeutic effect in terms of efficacy and selectivity.

Figure 10. Transformation of a Fuzzy Multi-Objective Hierarchical Optimization Problem into a Maximizing Decision-Making Problem Using Fuzzy Membership Functions. The lower bound (LB), upper bound (UB), and standard value (ST) are defined by the user, based either on available clinical data or estimated from the metabolic templates of cancer (CA) and healthy (HT) cells.

Figure 10. Transformation of a Fuzzy Multi-Objective Hierarchical Optimization Problem into a Maximizing Decision-Making Problem Using Fuzzy Membership Functions. The lower bound (LB), upper bound (UB), and standard value (ST) are defined by the user, based either on available clinical data or estimated from the metabolic templates of cancer (CA) and healthy (HT) cells.

Healthy cells are perturbed by the application of a therapeutic target, leading to altered flux distributions. These perturbed flux states—referred to as PH (perturbed healthy) cells—are used to evaluate metabolic deviations in comparison to both cancer (CA) cells and unperturbed healthy tissue (HT) cells. To quantify these deviations, fuzzy dissimilarity and similarity measures are applied to assess these flux differences. Two-sided linear membership functions are used for this purpose, with green lines representing fuzzy dissimilarity and red lines representing fuzzy similarity. These functions are used to define metabolic deviation, denoted as ${\eta }_{MD}\in \left[\mathrm{0,1}\right]$, which serves as an indicator for estimating potential side effects of a given target. Based on these membership function definitions, the fuzzy multi-objective hierarchical optimization problem is reformulated as a maximizing decision-making problem. This problem is then solved using the nested hybrid differential evolution to identify candidate therapeutic targets that balance efficacy and safety.

4. Conclusions

This study presents an integrative systems biology framework for identification of therapeutic targets capable of simultaneously treating pancreatic ductal adenocarcinoma (PDAC) and its associated cachexia (PDAC-CX). By extending the Recon3D network to incorporate pathways for protein synthesis and degradation, and by constructing cell-specific genome-scale metabolic models informed by transcriptomic data, the framework effectively captures critical metabolic alterations underlying both tumor progression and muscle wasting.

The application of parsimonious metabolite flow variability analysis across multiple nutritional environments enabled the identification of robust, medium-independent biomarkers. These biomarkers revealed key aspects of metabolic reprogramming in both PDAC and cachectic muscle cells, including disruptions in redox homeostasis, lipid metabolism, nucleotide biosynthesis, and inflammatory signaling pathways.

Using the anticancer target discovery platform, enzyme targets were systematically evaluated through a fuzzy multi-objective optimization framework. While single-enzyme interventions exhibited context-specific and nutrient-dependent effects, combinatorial targeting strategies revealed eight enzyme target pairs with consistent, medium-independent efficacy. Notably, combinations involving knockouts of SLC29A2, SGMS1, CRLS1, and the RNF20–RNF40 complex, alongside upregulation of CERK and PIKFYVE, demonstrated strong therapeutic potential by achieving favorable trade-offs between efficacy (cell viability) and safety (metabolic deviation).

This study introduces a novel computational framework for dual-target discovery in cancer and cachexia, facilitating the rational design of metabolism-based therapies with enhanced clinical relevance. The findings highlight the necessity of accounting for systemic effects and nutrient variability in anticancer drug development, thereby advancing translation strategies aimed at simultaneously suppressing tumor growth and alleviating cancer-associated muscle wasting.