Abstract Colorectal cancer (CRC) remains a leading cause of global cancer mortality, underscoring the need for novel therapeutic strategies. This study used a systems pharmacology approach integrated with molecular docking and molecular dynamics (MD) simulations to evaluate the potential of repurposing terfenadine and domperidone for inhibition of apoptotic gene associations in CRC. Network pharmacology analysis identified 4 principal targets—SLC6A4 (5I6X), DRD2 (7DFP), HTR2A (6WGT), and EGFR (6LUD)—involved in the apoptotic regulatory network. Molecular docking studies demonstrated high binding affinities of both terfenadine and domperidone against all selected targets (−7.1 to −11.5 kcal/mol), with the strongest interaction observed with DRD2, where both compounds exhibited a binding affinity of −11.5 kcal/mol. Detailed interaction profiling revealed critical hydrogen bonding and hydrophobic interactions stabilizing the drug-target complexes. Molecular dynamics simulations over a 100 ns timescale confirmed the structural stability and conformational fidelity of the docked complexes, evidenced by low root mean square deviation values and consistent hydrogen bond occupancy. Furthermore, post-MD simulation study supports the stable score landscape and stability of complex. In conclusion, this integrative computational analysis highlights terfenadine and domperidone as promising candidates capable of modulating key apoptotic pathways in CRC. The findings provide a strong rationale for subsequent in vitro and in vivo studies to validate their therapeutic potential and facilitate clinical translation in CRC management. Keywords: Colorectal cancer, domperidone, drug repurposing, terfenadine, molecular docking Introduction Colorectal cancer (CRC) represents one of the most prevalent malignancies worldwide, ranking third in incidence and second in cancer-related mortality. Recent global estimates report over 1.9 million new cases and nearly 935 000 deaths in 2024 alone, reflecting the urgent need for more effective therapeutic strategies. Colorectal cancer is the third most prevalent cancer globally, responsible for about 10% of all cancer diagnoses and the second leading cause of cancer deaths.^ [39]1 Despite advances in surgery, chemotherapy, and targeted therapy, the prognosis for advanced-stage CRC remains poor, with a 5-year survival rate below 15%. These epidemiological challenges highlight an urgent need for novel therapeutic strategies, including drug repurposing approaches that can offer faster clinical translation.^[40]2 The pathophysiology of CRC involves a complex, multistep process characterized by the accumulation of genetic and epigenetic alterations.^ [41]3 Key mutations typically affect tumor suppressor genes such as APC, TP53, and SMAD4, as well as oncogenes like KRAS. Dysregulation of signaling pathways—including Wnt/β-catenin, PI3K/AKT, MAPK (mitogen-activated protein kinase), and TGF-β (transforming growth factor beta)pathways—promotes uncontrolled cell proliferation, resistance to apoptosis, angiogenesis, and metastasis.^[42]4,[43]5 Chronic inflammation, aberrant cellular communication, and tumor microenvironment remodeling further contribute to CRC initiation and progression. Importantly, defects in apoptotic mechanisms allow tumor cells to evade programmed cell death, a hallmark feature of malignancy and a critical therapeutic target.^ [44]6 Treatment options for CRC at present consist of surgical resection, chemotherapy (such as 5-fluorouracil, oxaliplatin, irinotecan), targeted therapy (anti-EGFR (epidermal growth factor receptor) antibodies such as cetuximab, anti-VEGF (anti-vascular endothelial growth factor) agents such as bevacizumab), and immunotherapies (such as checkpoint inhibitors for microsatellite instability-high tumors). While these therapies have optimized survival for earlier stages, they are not effective in metastatic or chemoresistant CRC. These include drug resistance, extreme toxicity, tumor heterogeneity, and failure to efficiently target the principal oncogenic pathways. These have the consequence of resulting in treatment failure and relapse of disease.^[45]7,[46]8 In addition, the time-consuming and expensive nature of developing new anticancer drugs makes alternative strategies more desirable. Drug repurposing has emerged as a promising strategy to accelerate the development of novel cancer therapies by identifying new uses for existing drugs with established safety profiles.^ [47]9 In this context, terfenadine, an antihistamine, and domperidone, a dopamine D2 receptor antagonist, have attracted attention for their potential anticancer properties. Terfenadine has been shown to induce oxidative stress, mitochondrial dysfunction, and caspase-dependent apoptosis in various cancer models by disrupting intracellular calcium homeostasis and suppressing pro-survival signaling pathways.^[48]10,[49]11 Meanwhile, domperidone, through D2 receptor antagonism, may inhibit oncogenic cascades such as AKT/PI3K and MAPK pathways, facilitating apoptosis and impairing tumor cell proliferation. Furthermore, it induces apoptosis through STAT3 signaling pathway.^[50]12,[51]13 In CRC models, targeting neurotransmitter-associated receptors like DRD2 has been associated with reduced tumor growth and enhanced apoptotic response.^ [52]14 This study used a systems pharmacology approach to identify potential apoptotic gene targets associated with terfenadine and domperidone in CRC. Key targets were prioritized using network pharmacology analysis, followed by molecular docking and molecular dynamics (MD) simulations to evaluate drug-target interactions and binding stability. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis highlighted key signaling pathways, suggesting that modulation of neurotransmission, hormonal signaling, and intercellular communication may underlie the therapeutic effects of these repurposed drugs. Material and Methods Screening for molecular properties, bioactivity and toxicity profiling of terfenadine and domperidone The 3-dimensional (3D) structure of terfenadine and domperidone is obtained from the PubChem database ([53]http://pubchem.ncbi.nlm.nih.gov/) by searching for information about these drugs and comprehensive literature evaluations of cancer characteristics. With MoleSoft LLC ([54]http://molsoft.com/), we used a SMILE query to determine each molecular characteristics and drug similarity scores of terfenadine and domperidone. ProTox v3.0 ([55]https://tox.charite.de/protox3/) was used for toxicity study. ProTox v3.0 is a computer models based on information gathered from in vitro or in vivo investigations to forecast the toxicity and metabolic stability.^ [56]15 PerMM server ([57]https://permm.phar.umich.edu/) was used to predict the cell membrane permeability.^ [58]16 Drug-associated targets The sdf file format of terfenadine and domperidone was uploaded to SwissTargetPrediction database ([59]http://www.swisstargetprediction.ch/). The Homo sapiens option was selected from the “Select a species” option bar and then targets associated with the drugs were predicted. All the predicted targets were saved in Excel sheet to further analysis.^ [60]17 Retrieval of CRC genes The GeneCard database ([61]https://www.genecards.org/) provided the target gene for CRC. We used “colorectal cancer” as a prompt. The pseudogene, duplicate gene, and uncategorized gene were excluded from the analysis. About 49 001 genes were selected for further analysis. Target overlap and protein-protein interaction network construction Common targets between drug compounds and CRC were identified using a Venny 2.1 ([62]https://bioinfogp.cnb.csic.es/tools/venny/). The intersected targets were then used to construct a compound-target-disease interaction network using Cytoscape v3.10.2. In the network, nodes represent compounds or proteins, and edges represent interactions.^ [63]18 Furthermore, the overlapping targets were submitted to the STRING database v12.0 ([64]https://string-db.org/) to analyze their protein-protein interaction (PPI), node and edge characteristic. Initial screening was conducted with a minimum interaction score of 0.4. This initial screening provided a large less accurate network. We also did not get accurate node and edge from this analysis. After that, we increased the confidence score up to 0.7 to get more accurate and closed network. The resulting network was visualized in Cytoscape v3.10.2, and hub genes were identified based on topological parameters such as degree centrality using the CytoHubba plugin.^ [65]19 Gene ontology and KEGG pathway enrichment analysis A prominent method for analyzing genomic data, especially large-scale transcriptome data, is gene ontology (GO) analysis. Biological process (BP), cellular component (CC), and molecular function (MF) were the 3 types of potential targets that were examined for GO functional enrichment. KEGG and GO enrichment analysis were conducted in the ShinyGO v0.82 database ([66]https://bioinformatics.sdstate.edu/go/) to improve understanding of the roles of cross-targets and signaling pathways involved in the common top 10 targets of terfenadine, domperidone, and CRC. ShinyGO is a web-based tool designed for gene enrichment analysis and visualization, allowing users to perform GO and KEGG pathway analysis with interactive graphical outputs. It supports comprehensive functional annotation, hierarchical clustering, and network-based interpretation of gene sets across multiple species.^[67]20,[68]21 Survival analysis Using a standard processing pipeline, GEPIA v2 ([69]https://gepia2.cancer-pku.cn/) is a web service that analyzes the RNA sequencing expression data of 8587 normal samples and 9736 cancers from the TCGA and GTEx projects. Using GEPIA v2 cancer genomics servers, the predictive relevance of each hub gene was evaluated, with a focus on overall survival (OS) for CRC. Colorectal cancer patients were divided into 2 groups, high and low expression groups, and compared using the Kaplan-Meier survival table. The log rank p-value for survival and the risk ratio (hazard ratio (HR); 95% confidence interval) were calculated using log rank p 0.05 as the threshold for statistical significance.^ [70]22 Molecular docking Ligands, including terfenadine and domperidone, were obtained from the PubChem database ([71]https://pubchem.ncbi.nlm.nih.gov/) in SDF format. These structures were converted to PDBQT format using OpenBabel integrated within PyRx v0.8. Prior to docking, score minimization was conducted using the Merck Molecular Force Field (MMFF94) via Chimera v1.16. The 3-dimensional crystal structure of the target protein was retrieved from the RCSB Protein Data Bank ([72]https://www.rcsb.org/).^ [73]23 The 3-dimensional crystal structures of the target proteins—SLC6A4 (5I6X),^ [74]24 DRD2 (7DFP),^ [75]25 HTR2A (6WGT),^ [76]26 and EGFR (6LUD)^ [77]27 — were retrieved from the RCSB Protein Data Bank. The downloaded structures were prepared for docking by removing water molecules, heteroatoms, and bound ligands using BIOVIA Discovery Studio Visualizer 2021 ([78]https://discover.3ds.com/discovery-studio-visualizer-download). The cleaned structures were then converted from PDB to PDBQT format using plugin OpenBabel in PyRx v0.8.^ [79]28 Docking was conducted using PyRx v0.8 with an exhaustiveness parameter set to 8 for thorough sampling of the binding site. The specific grid box dimensions and centers were defined as follows: * For DRD2 (7DFP): the grid center was set at (center_x = −87.5631, center_y = −15.9276, center_z = 234.1955). * For SLC6A4 (5I6X): the grid center was set at (center_x = −35.9027, center_y = −22.0573, center_z = 2.809). * For HTR2A (6WGT): the grid center was set at (center_x = 24.814, center_y = 43.0225, center_z = 30.0521). * For EGFR (6LUD): the grid center was set at (center_x = −57.7471, center_y = −8.0019, center_z = −24.695). The docking of ligands to the target protein was carried out in PyRx v0.8, and resulting protein-ligand complexes were visualized and analyzed using BIOVIA Discovery Studio Visualizer 2021.^ [80]29 Physiochemical properties and validation of proteins The physiochemical properties of the selected proteins—7DFP ([81]P14416), 5I6X ([82]P31645), 6WGT ([83]P28223), and 6LUD ([84]P00533)—were analyzed using the ExpasyProtParam web server ([85]https://web.expasy.org/protparam). Protein sequences were retrieved from the UniProtKB database and input into ProtParam to compute parameters including molecular weight, theoretical isoelectric point (pI), extinction coefficient (assuming all cysteine residues are reduced), instability index, aliphatic index, Grand Average of Hydropathy (GRAVY), and estimated half-life in Escherichia coli. The resulting data provided insights into the structural and stability characteristics of each protein.^ [86]30 Docking protocol validation To validate the molecular docking protocol, we performed re-docking of the native ligands co-crystallized with the target proteins obtained from the RCSB Protein Data Bank. Each native ligand was separated from its receptor using BIOVIA Discovery Studio Visualizer 2021. The native ligand was then re-docked into the prepared binding site using PyRx v0.8 with the same grid parameters and docking settings applied to domperidone and terfenadine. The root mean square deviation (RMSD) between the predicted and experimentally observed binding poses of the native ligands was calculated using PyMOL to assess the reliability of the docking protocol. An RMSD value of less than 2.0 Å was considered acceptable, indicating that the docking protocol could accurately reproduce the crystallographic binding orientation.^ [87]31 MD simulation The dynamic behavior of the DRD2 (7DFP) protein in complex with terfenadine and domperidone was investigated through 100 ns all-atom MD simulations using the Desmond simulation package (Schrödinger, v2021-4). The prepared protein-ligand complexes were imported into the system builder module. Each complex was embedded in an orthorhombic box with a 10 Å buffer distance and solvated using the TIP3P water model. The system was neutralized by adding counterions and a physiological salt concentration of 0.15 M NaCl was maintained to mimic in vivo conditions. Score minimization was performed using the OPLS_2005 force field to relieve any steric hindrances. The system was equilibrated under NPT ensemble conditions at 300 K temperature and 1 atm pressure.^ [88]32 A trajectory snapshots were recorded every 100 ps. Separate simulations were conducted for the DRD2-terfenadine and DRD2-domperidone complexes to capture individual dynamic behaviors. Post-MD simulation analysis Trajectory analyses were carried out using the Simulation Interaction Diagram and in-built analysis modules in Maestro. The RMSD and root mean square fluctuation (RMSF) were used to evaluate overall protein-ligand stability and residue-specific flexibility, respectively. The Radius of Gyration (rGyr) was analyzed to monitor the compactness of the protein structure over time. Solvent accessible surface area (SASA), molecular surface area (MolSA), polar surface area (PSA), and intra-molecular hydrogen bonds (intraHB) were measured to assess structural compactness, surface exposure, and hydration effects.^[89]33 Principal component analysis, free score landscape, and probability density function analysis For deeper insight into the conformational landscape, principal component analysis (PCA) was performed on Cα atoms of the DRD2 backbone. The first 2 eigenvectors (PC1 and PC2) were plotted to identify dominant motion patterns. Free score landscape (FEL) was constructed using the PC1 and PC2 coordinates to visualize low-score conformational states and assess structural stability. The probability density function (PDF) analysis was conducted to evaluate the frequency distribution and transitions between metastable states during simulation.^ [90]34 Binding free energy calculation The binding free energy of the DRD2-terfenadine and DRD2-domperidone complexes was computed using molecular mechanics/generalized born surface area (MM/GBSA) methodology through the Prime MM/GBSA module in Schrödinger. A total of 100 frames were extracted from each trajectory and used for MM/GBSA calculation. The ΔG_bind was calculated as follows: [MATH: ΔG=Gcom plex−(Gprotein+Gligand) :MATH] Where each component includes the MM score, solvation free score, and surface area score contributions. These energies offered insight into the comparative binding affinity of the ligands with DRD2 post-dynamics.^[91]35,[92]36 Results Pharmacokinetic profiling, bioactivity prediction, and toxicity assessment of terfenadine and domperidone To evaluate the therapeutic suitability of terfenadine and domperidone for CRC, a comprehensive screening of their pharmacokinetic, physicochemical, and toxicity properties was conducted ([93]Supplementary file, Tables S1-S4). The drug-likeness scores calculated through ADME profiling tools indicated favorable properties, with domperidone and terfenadine scoring 1.26 and 1.61, respectively. Oral acute toxicity predictions revealed LD[50] values of 5000 mg/kg for domperidone and 715 mg/kg for terfenadine, suggesting a comparatively wider safety margin for domperidone. In addition to basic pharmacological characteristics, membrane permeability and interaction properties were computationally simulated to assess drug absorption and bioavailability ([94]Supplementary file, Table S1). These simulations supported the ability of both drugs to effectively interact with membrane components, essential for cellular uptake and target engagement in a therapeutic context. Furthermore, target prediction analysis ([95]Supplementary, Figure S1) identified a broad spectrum of interactive targets for both drugs, reinforcing their potential for modulating multiple CRC-associated pathways. The convergence of these predictions with targets identified in the PPI and compound-target networks strengthens the pharmacological relevance of terfenadine and domperidone in CRC intervention. Identification of shared therapeutic targets in CRC A Venn diagram analysis was conducted to identify common molecular targets among terfenadine, domperidone, and CRC-associated genes ([96]Figure 1A). Terfenadine was associated with 66 targets, while domperidone had 101 predicted targets. In contrast, the CRC gene pool encompassed 49 001 targets. The intersection revealed 35 common genes, representing a small but significant subset (0.1%) shared across all 3 datasets. Notably, terfenadine and domperidone shared only one direct target, emphasizing their unique target profiles. The overlapping genes provide a focused starting point for investigating the potential of these drugs in modulating CRC-related pathways. Figure 1. [97]Common targets of repurposed drugs terfenadine, domperidone and colorectal cancer (A), Protein-protein interaction analysis of target genes using string database (B), Top 10 hub genes in network string interaction ranked by degree method (C), and drug compound-target gene network construction (D) [98]Open in a new tab Common targets of repurposed drugs terfenadine, domperidone and colorectal cancer (A), Protein-protein interaction analysis of target genes using string database (B), Top 10 hub genes in network string interaction ranked by degree method (C), and drug compound-target gene network construction (D). PPI network analysis To examine the functional connectivity of the 35 shared targets ([99]Supplementary file, Table S5), PPI network was constructed using the STRING database ([100]Figure 1B). The network comprised 34 nodes and 109 edges, indicating a highly interconnected structure. The average node degree was 6.41, and the average local clustering coefficient was calculated as 0.546. Compared with the expected number of edges (13), the actual network showed a substantial enrichment, confirmed by a PPI enrichment P value of <1.0 × 10^−16. This high level of statistical significance suggests that the interactions among these proteins are nonrandom and biologically meaningful, supporting their potential collective involvement in CRC pathophysiology. Hub gene prioritization through network topology A gene interaction network was further mapped to prioritize key regulators among the 35 targets ([101]Figure 1C). Top 10 best genes were identified through maximal clique centrality (MCC) analysis ([102]Table 1). Several hub genes, including DRD2, HTR2A, OPRM1, EGFR, and ADRA2A, displayed extensive connectivity, suggesting their critical roles in the overall network architecture. Interestingly, genes related to neurotransmission and receptor signaling such as HTR2C, SLC6A2, and DRD2 formed interconnected modules, hinting at neuromodulatory mechanisms potentially influencing cancer progression. Table 1. Best 10 targets identified from network pharmacology maximal clique centrality (MCC) analysis along with their rank score. Rank Name of targets Score 1 SLC6A4 17 2 DRD2 16 3 HTR2A 13 3 SLC6A2 13 5 HTR2C 12 6 ADRA2A 11 7 HTR1A 10 8 OPRM1 9 9 ADRA2C 8 9 EGFR 8 [103]Open in a new tab Drug-target network construction To validate drug-specific associations, a bipartite drug-target interaction network was constructed for terfenadine and domperidone ([104]Figure 1D) using the Cytoscape software. In the network, there were 70 edges and 37 nodes. Domperidone and terfenadine had degree values of 35 and 2, respectively. Terfenadine showed preferential binding to targets such as CHRM4, HTR2A, DRD2, and SLC6A2, whereas domperidone was associated with HTR1A, DRD3, JAK2, and OPRK1 among others. Both compounds shared common targets including EGFR, PRKCD, JAK2, and HTR2C, highlighting convergence on cancer-related signaling pathways such as cell proliferation, apoptosis, and neurotransmitter regulation. These findings suggest that both drugs may exert anticancer effects through modulation of overlapping but distinct molecular circuits. KEGG pathway enrichment reveals neurological and signaling pathway associations The KEGG pathway analysis of 35 shared targets revealed significant enrichment in neurotransmission-related pathways, notably neuroactive ligand-receptor interaction and gap junctions. Additional pathways, including estrogen, calcium, and cAMP signaling, were linked to cancer progression and apoptosis ([105]Figure 2A). Notably, the neuroactive ligand-receptor interaction pathway showed the highest significance (−log10FDR ≈ 8) and involved the largest number of overlapping genes (n = 7), emphasizing the central role of neurotransmission-related mechanisms in the disease-drug interaction landscape. These findings suggest that terfenadine and domperidone may modulate CRC via neuronal, hormonal, and communication pathways. Figure 2. [106]KEGG Enrichment pathway (A & B). Graphical comparison of biological process fold-enrichment scores of gene ontology terms (A) and KEGG pathway categories (B), with N of Genes represented color-coded dots and FDR values indicated by the text. [107]Open in a new tab KEGG enrichment pathway (A). Biological process fold enrichment analysis of gene ontology (B). GO enrichment highlights neurotransmission and secretory pathways The GO enrichment analysis of the 35 shared target genes revealed significant involvement in neurotransmission and secretory processes. Monoamine transport and catecholamine transport were the most enriched BPs, showing high fold enrichment and significance (−log[10]FDR > 11), suggesting strong modulation of neurotransmitter dynamics ([108]Figure 2B). Additional enriched terms included G protein-coupled receptor (GPCR) signaling, signal release, and regulation of secretion, all critical for cell communication and tumor progression. The hierarchical clustering ([109]Supplementary file, Figure S3) grouped related synaptic processes such as trans-synaptic signaling, chemical synaptic transmission, and cell-cell signaling, indicating coordinated regulation of neuronal and secretory pathways. GO CC analysis reveals synaptic localization Cellular component enrichment analysis identified significant enrichment of gene targets in synapse-associated regions. The most prominent terms included the integral and intrinsic components of the presynaptic and postsynaptic membranes, with high fold enrichment and strong statistical significance (−log[10]FDR > 8) ([110]Figure 3A). Additional enriched components included synaptic membrane, membrane raft, and plasma membrane, indicating that the identified genes are predominantly localized at neuronal interfaces involved in signal transmission. Hierarchical clustering further highlighted strong correlations among synaptic and membrane-related compartments ([111]Supplementary file, Figure S4 and [112]S5). Figure 3. [113]compare enrichment of gene ontology of cellular component vs molecular function in table [114]Open in a new tab Cellular component enrichment analysis of gene ontology (A). Molecular function enrichment analysis of gene ontology (B). GO MF analysis highlights neurotransmitter and receptor binding activity The GO enrichment analysis for MFs revealed a strong association with neurotransmitter receptor activity, particularly involving G protein-coupled serotonin and adrenergic receptors. Top enriched terms included Gq/11-coupled serotonin receptor activity, alpha-adrenergic receptor activity, and amine binding, indicating the involvement of neurotransmission-related mechanisms ([115]Figure 3B). Hierarchical clustering showed a tight grouping of GPCR-related functions reinforcing the role of these receptors in the therapeutic potential of terfenadine and domperidone ([116]Figure 5S, supplementary file). High fold enrichment values and significant p-values (−log[10]FDR > 7) further validate the robustness of these associations. The findings suggest that the compounds may modulate CRC progression via interaction with neurotransmitter receptors, particularly through GPCR-mediated signaling pathways. KEGG-gap junction The KEGG Gap Junction pathway diagram illustrates the intricate network of proteins and signaling molecules involved in intercellular communication via gap junctions, composed primarily of connexin subunits forming hemichannels and full channels ([117]Figure 4). These channels enable the direct exchange of ions, second messengers, and small metabolites between adjacent cells, playing a critical role in maintaining tissue homeostasis, embryonic development, and electrical/metabolic coupling. The highlighted proteins—RTK, DRD2, and HTR2—indicate key involvement of receptor tyrosine kinases and neuroactive ligand-receptor interactions (notably dopaminergic and serotonergic pathways), suggesting modulation of gap junction activity by external stimuli such as growth factors and neurotransmitters. The pathway links to multiple downstream cascades, including MAPK, PKA/PKC signaling, and calcium signaling, affecting processes like apoptosis, differentiation, and synaptic transmission. The diagram underscores the regulatory complexity and cross-talk between membrane-bound receptors, intracellular kinases, and cytoskeletal elements (eg, tubulins), reflecting the centrality of gap junctions in cellular communication networks. Figure 4. [118]Cytokine Receptor Interaction - RTK, MAP Kin- ERK1/2, and AP1 activated by RTK. GPCRs activated by RTK. [119]Open in a new tab KEGG pathway analysis of the gap junction signaling network highlighting RTK, DRD2, and HTR2 interactions. Survival analysis The disease-free survival of the hub genes in CRC was examined using the Kaplan-Meier survival plot in order to determine whether or not hub genes positively impacted patient prognosis. [120]Figure 5 presents Kaplan-Meier survival curves comparing high and low DRD2 expression in CRC patients. In [121]Figure 5A, the OS rates between the low DRD2 group (blue line) and high DRD2 group (red line) were compared. The survival curves indicated that patients in the high DRD2 group had slightly poorer OS compared with the low DRD2 group, although the difference was not statistically significant (Logrank P = .24). The HR for the high DRD2 group was 1.6, with a corresponding p-value of 0.24, suggesting a trend toward increased risk but without statistical significance. The number of patients in the high and low groups was 35 and 36, respectively. Figure 5. [122]Represents overall survival (35 patients) and disease free survival (35 patients: 1.1) of DRD2 gene target. [123]Open in a new tab Represents overall survival (A) and diseases free survival of DRD2 gene target (B). In [124]Figure 5B, the disease-free survival rates were similarly compared between the 2 groups. The curves showed that disease-free survival was comparable between the high and low DRD2 expression groups, with no significant differences observed (Logrank P = .82). The HR for the high DRD2 group was 1.1, and the p-value for the HR was 0.81, indicating a negligible association between DRD2 expression and disease-free survival. The number of patients remained consistent with 35 in the high group and 36 in the low group. Analysis of physiochemical properties and validation of proteins The physiochemical properties of the 4 proteins analyzed from UniProtKB entries [125]P14416 (7DFP), [126]P31645 (5I6X), [127]P28223 (6WGT), and [128]P00533 (6LUD) revealed notable differences in their structure and stability ([129]Table 2). The number of amino acids ranged from 443 in 7DFP to 1210 in 6LUD, with corresponding increases in molecular weight, from 50 619.43 to 134 277.40 Da. The theoretical isoelectric points (pI) varied significantly, with 7DFP exhibiting a basic pI of 9.55, while 5I6X and 6LUD showed more acidic values of 5.89 and 6.26, respectively. All proteins demonstrated long estimated half-lives in E. coli, exceeding 10 hours. The extinction coefficients under reducing conditions ranged from 0.932 M^−1 cm^−1 in 6LUD to 2.143 M^−1 cm^−1 in 5I6X, indicating variability in aromatic residue content. The instability index suggested that 7DFP might be the least stable (53.08), while 5I6X appeared more stable (33.13). Aliphatic indices were relatively high across the proteins, reflecting potential thermostability, with values exceeding 98. GRAVY (Grand Average of Hydropathy) scores ranged from slightly hydrophilic in 6LUD (−0.316) to more hydrophobic in 5I6X (0.422), indicating differences in solubility and interaction with aqueous environments. Table 2. Physiochemical properties of proteins. Parameters Proteins and UniProtKB 7DFP: [130]P14416 5I6X: [131]P31645 6WGT: [132]P28223 6LUD: [133]P00533 Number of amino acids 443 630 471 1210 Molecular weight 50619.43 70324.86 52603.17 134277.40 Formula C[2293]H[3628]N[620]O[617] S[27] C[3278]H[4949]N[781]O[879] S[30] C[2367]H[3742]N[602]O[688] S[30] C[5875]H[9284]N[1646]O[1786] S[85] Theoretical pI 9.55 5.89 7.83 6.26 Ext. coefficient (M^-1cm^-1) 1.261 2.143 1.128 0.932 Estimated half-life (E. coli) >10 hours >10 hours >10 hours >10 hours Instability index 53.08 33.13 39.77 44.59 Aliphatic index 98.78 100.44 101.83 GRAVY 0.060 0.422 0.221 -0.316 [134]Open in a new tab Molecular docking In this study, we used molecular docking to examine the molecular interaction of terfenadine and domperidone against SLC6A4 (5I6X), DRD2 (7DFP), HTR2A (6WGT), and EGFR (6LUD) in cancer cell and the anti-apoptotic properties ([135]Table 3). The molecular docking analysis revealed that both terfenadine and domperidone demonstrated strong binding affinities against the selected protein targets. For DRD2 (PDB ID: 7DFP), both drugs exhibited identical binding scores of −11.5 kcal/mol. Terfenadine formed a key hydrogen bond with ILE184 (2.6 Å), whereas domperidone interacted with CYS182 at 2 distances 2.4 Å and 3.3 Å, indicating stable binding at the active site. In the case of HTR2A (PDB ID: 6WGT), domperidone exhibited superior binding score (−11.3 kcal/mol) compared with terfenadine (−9.0 kcal/mol), with domperidone forming 2 strong hydrogen bonds with ASN343 (2.3 Å) and SER131 (2.5 Å), whereas terfenadine interacted with CYS227 (1.9 Å) and GLU223 (2.7 Å) ([136]Figure 6A and [137]B). For EGRF (PDB ID: 6LUD), binding scores were relatively weaker, with terfenadine showing −7.1 kcal/mol and domperidone −8.4 kcal/mol. Domperidone formed multiple hydrogen bonds with MET793 (2.1 Å and 2.3 Å), GLN791 (3.2 Å), and ASP800 (1.9 Å), indicating more extensive interactions compared with terfenadine ([138]Figure 7A and [139]B). Against SLC6A4 (PDB ID: 5I6X), terfenadine showed a binding score of −11.0 kcal/mol via a strong hydrogen bond with GLU493 (2.3 Å), while domperidone displayed a slightly lower affinity (−9.7 kcal/mol) interacting with ASP98 (2.2 Å) ([140]Figure 7C and [141]D). Overall, domperidone consistently demonstrated either comparable or stronger binding interactions across most targets, suggesting a slightly higher potential as a therapeutic candidate for targeting apoptotic genes in CRC. Table 3. Binding score and hydrogen bond interaction of terfenadine and domperidone against selected 4 different cancer targets (SLC6A4, DRD2, HTR2A, and EGFR). PDB ID Terfenadine (binding score kcal/mol) H-bind interaction (distance) Domperidone (binding score kcal/mol) H-bond interaction (distance) 7DFP -11.5 ILE184 (2.6Å) -11.5 CYS182 (2.4 & 3.3 Å) 5I6X -11.0 GLU493 (2.3 Å) -9.7 ASP98 (2.2 Å) 6WGT -9.0 CYS227 (1.9 Å), GLU223 (2.7) -11.3 ASN343 (2.3 Å), SER131 (2.5 Å) 6LUD -7.1 ASN841 (2.4 Å), ASP855 (2.0 Å) -8.4 MET793 (2.1 & 2.3 Å), GLN791 (3.2 Å), ASP800 (1.9 Å) [142]Open in a new tab Figure 6. Molecule interaction between 6WGT and terfenadine A, domperidone B: 3D top and 2D stick models for each. [143]Open in a new tab Molecular interaction of terfenadine (A) and domperidone (B) against 6WGT protein. Both surface view model and stick model are presented. Figure 7. The image depicts molecular interactions between terfenadine (A) and domperidone (B) against the 6LUD protein, and between terfenadine (C) and domperidone (D) against the 5I6X protein. [144]Open in a new tab Molecular interaction of terfenadine (A) and domperidone (B) against 6LUD protein. (C) and (D) represent the interaction of terfenadine and domperidone against 5I6X respectively. Docking protocol validation result To ensure the reliability of the docking protocol, native ligand re-docking was carried out for each target protein. The RMSD values between the docked and crystallographic poses of the native ligands were as follows: * DRD2 (7DFP) – 1.62 Å * HTR2A (6WGT) – 1.48 Å * EGFR (6LUD) – 1.75 Å * SLC6A4 (5I6X) – 1.33 Å All values were below the 2.0 Å threshold, confirming the accuracy and reproducibility of the docking methodology. These results validate the reliability of the docking protocol applied to domperidone and terfenadine, supporting the credibility of the predicted binding poses and interaction profiles. MD simulation Among 4 targets, DRD2 (7DFP)-drug complexes showed strong binding interaction in molecular docking study. Therefore, we consider DRD2-drug complexes for MD simulation and post-MD simulation analysis. Root mean square deviation Molecular dynamics simulations were performed for 100 ns to evaluate the stability and dynamic behavior of the DRD2 (7DFP) protein in complex with domperidone and terfenadine. The RMSD plots ([145]Figure 8A and [146]B) revealed that both complexes achieved relative stability throughout the simulation. The 7DFP-domperidone complex ([147]Figure 8A) exhibited an average backbone RMSD fluctuation between 1.5−2.5 Å, indicating stable protein conformation. The ligand RMSD remained within 1.0−1.8 Å, suggesting that domperidone maintained a consistent binding pose. In contrast, the 7DFP-terfenadine complex ([148]Figure 8B) showed slightly higher backbone fluctuations, ranging from 1.8 to 3.2 Å, while the ligand RMSD gradually increased, peaking around 4.0 Å, suggesting more pronounced ligand flexibility or partial repositioning within the binding pocket. Figure 8. [149]RMSD and RMSF of protein-ligand complexes, including domperidone and terfenadine, over time and residue index. [150]Open in a new tab RMSD of 7DFP-domperidone complex (A), 7DFP-terfenadine complex (B), RMSF of 7DFP-domperidone complex (C), 7DFP-terfenadine complex (D). Root mean square fluctuation The RMSF profiles ([151]Figure 8C and [152]D) provided insights into the residue-wise flexibility of the protein backbone in both complexes. In the 7DFP-domperidone complex ([153]Figure 8C), higher fluctuations were observed at the N- and C-terminal regions and a few loop segments, while most residues exhibited limited flexibility below 2.0 Å, indicating overall rigidity in the binding site. Similar fluctuation patterns were observed in the 7DFP-terfenadine complex ([154]Figure 8D), though several regions, especially around residue index ~50 and ~270, exhibited slightly elevated RMSF values, reflecting localized flexibility. These RMSF results closely align with the amino acid interaction data obtained from molecular docking. Key interacting residues such as ASP114, PHE390, and CYS182—identified through hydrogen bond and hydrophobic interaction analysis—showed minimal fluctuations during the simulation, suggesting that these residues maintain consistent contacts with the ligands. The reduced RMSF values at these positions indicate structural stability and sustained interaction between the ligand and the binding pocket. Conversely, residues distant from the binding interface showed higher RMSF, consistent with their peripheral or flexible loop locations, and did not significantly contribute to ligand stabilization. [155]Figure 9A and [156]B, on the contrary, shows the RMSF of the ligands themselves — domperidone and terfenadine — fitted to the protein. These trajectories depict how much each atom of the ligand fluctuates during the simulation, indicating ligand conformational stability within the binding pocket. The ligand RMSF analysis revealed that domperidone ([157]Figure 9A) exhibited relatively lower atomic fluctuations when fit on the binding pocket compared with the protein, indicating stable binding, whereas terfenadine ([158]Figure 9B) showed higher overall fluctuations, particularly when fit on the protein, suggesting greater conformational flexibility within the binding pocket. Secondary structure elements derived from the simulation ([159]Supplementary file, Figure S6 and [160]S7) indicated that the 7DFP structure was predominantly composed of α-helices and β-strands, contributing to its structural integrity. This correlation between docking and dynamics reinforces the notion that domperidone and terfenadine anchor firmly within the active site via persistent interactions with structurally stable residues, thereby supporting their potential efficacy as inhibitors of DRD2-mediated anti-apoptotic signaling in CRC. Figure 9. [161]RMSF plots showing protein and ligand fit for domperidone and terfenadine [162]Open in a new tab RMSF plot of domperidone (A) and terfenadine (B). Protein-ligand compactness The MD simulation analysis of the physicochemical properties for the 7DFP complexes with domperidone ([163]Figure 10A) and terfenadine ([164]Figure 10B) provided key insights into the stability and compactness of the ligand-protein complexes over the 100 ns trajectory. For the domperidone complex, the rGyr remained stable around ~5.6 Å, suggesting a consistent compactness of the ligand. Intra-molecular hydrogen bonding (intraHB) remained stable at 1, indicating sustained internal interactions within the ligand. The MolSA and SASA values fluctuated slightly but maintained averages around 392 Å^2 and 400 Å^2, respectively, indicative of stable exposure to the solvent. The PSA also showed minimal fluctuations around 135 Å^2, further supporting the structural stability of domperidone within the binding pocket ([165]Figure 10A). Figure 10. [166]compare the time courses of rGyr, intraHB, MolSA, SASA, and PSA of domperidone (A) and terfenadine (B). [167]Open in a new tab rGyr, intraHB, MolSA, SASA, PSA of domperidone (A) and terfenadine (B). In contrast, terfenadine exhibited higher variability in its properties. The rGyr fluctuated more prominently, ranging between 5.7 and 6.6 Å, reflecting greater conformational flexibility ([168]Figure 10B). IntraHB dropped to zero at several points, suggesting a loss or rearrangement of internal hydrogen bonds during the simulation. The MolSA and SASA values displayed more significant variations, peaking at 482.6 Å^2 and 386.7 Å^2, respectively, around 15.5 ns, indicating transient expansions in surface exposure. Notably, PSA dropped sharply to 72.2 Å^2 at the same time point, reflecting changes in polar surface characteristics, which may suggest a shift in binding orientation or solvation state. Domperidone exhibited more stable physicochemical behavior throughout the simulation, while terfenadine showed signs of greater flexibility and transient instability within the binding site. Protein-ligand interaction The molecular docking analysis demonstrated the interaction profiles of 7DFP with domperidone and terfenadine. Domperidone established key hydrogen bond interactions with residues such as ASP114, PHE390, and CYS182, contributing to stable binding. The histogram in [169]Figure 11A and [170]B quantifies the interaction fractions, highlighting dominant contributions from ASP114, PHE390, and surrounding hydrophobic residues. [171]Figure 11C and [172]D depicts the binding interactions and interaction histogram for terfenadine. Terfenadine exhibited strong hydrogen bonding with ASP114 and CYS182, as well as hydrophobic contacts with PHE390, HIS393, and TRP100. The histogram in panel [173]Figure 11D showed a high interaction fraction primarily at ASP114, PHE390, and neighboring residues, indicating stable and persistent binding throughout simulation. Both compounds demonstrated strong and consistent binding profiles with DRD2, supporting their potential as effective inhibitors targeting apoptotic pathways in CRC. Figure 11. [174]The image depicts molecular diagrams and histograms illustrating the interactions between the dopamine receptor D2 (7DFP) and the antihistamine drugs domperidone and terfenadine. [175]Open in a new tab Molecular interaction diagrams and interaction histograms of DRD2 (7DFP) with domperidone and terfenadine. Domperidone formed key hydrogen bonds with ASP114, PHE390, and CYS182, with interaction fractions shown in the histogram (A, B). Terfenadine interacted strongly with ASP114, CYS182, PHE390, and HIS393, with corresponding interaction fractions (C, D). Post-MD simulation analysis The post-MD simulation analyses of the domperidone-7DFP and terfenadine-7DFP complexes were conducted using PCA, FEL evaluations (both 2-dimensional (2D) and 3D), and PDF analyses. These techniques offered critical insights into the dynamic behavior, conformational flexibility, and stability of the complexes following MD simulations. FEL analysis The 2D FEL of the domperidone-7DFP complex revealed 2 major low-score basins, indicating the presence of well-defined and energetically favorable conformations. This suggested that the complex remained relatively stable during the simulation, exploring only a few conformational states. The corresponding 3D FEL supported this observation by displaying smoother valleys and a limited number of local minima, which reflected lower energetic variability and a more consistent structural configuration throughout the simulation ([176]Figure 12). In contrast, the FEL analysis of the terfenadine-7DFP complex exhibited a rugged score landscape with several low-score regions. The 2D FEL showed multiple scattered basins, while the 3D FEL presented sharp score peaks and deep valleys ([177]Figure 13). These features indicated that the complex sampled a broader range of conformations and underwent significant structural transitions during the MD simulation. This pointed to higher conformational flexibility and the presence of multiple metastable states. Figure 12. [178]free energy landscape, 3D free energy landscape, probability density function analysis, PCA projection, PCA principal components, domperidone-7DFP complex, 2D color plot. [179]Open in a new tab PCA, FEL, and PDF of domperidone-7DFP complex. Figure 13. [180]the protein structure resolved by X-ray diffraction of terfenadine-7DFP complex using pca [181]Open in a new tab PCA, FEL, and PDF of terfenadine-7DFP complex. Principal component analysis The PCA projection of the domperidone-7DFP complex showed a moderately clustered distribution of data points, implying that the complex underwent restricted movements during the simulation. The major variance in motion was primarily accounted for by the first 2 principal components (PC1 and PC2), indicating that the system followed a limited number of dominant motion pathways ([182]Figure 12). This aligned well with the FEL results, further supporting the structural rigidity and stability of the complex. The PCA plot for the terfenadine-7DFP complex displayed a more dispersed distribution of conformations across the principal components. The data points spanned a wider range in 3D space, reflecting greater structural flexibility and broader conformational sampling. This dispersion confirmed that terfenadine induced dynamic motions in the 7DFP receptor, consistent with the observations from the FEL analysis ([183]Figure 13). PDF analysis The PDF plot for the domperidone-7DPF complex revealed relatively sharp and narrow peaks for both PC1 and PC2, indicating that the complex predominantly resided in 1 or 2 conformational states. The density distributions suggested a high probability of maintaining a particular structural configuration, which reinforced the notion of conformational stability and limited fluctuations ([184]Figure 12). Similarly, the PDF analysis of the terfenadine-7DPF complex demonstrated broader and more complex peaks for PC1 and PC2. These features indicated a wider range of conformational transitions and a more flexible structural behavior during the simulation ([185]Figure 13). The presence of multiple peaks suggested that the complex existed in several intermediate states with comparable probabilities. MM/GBSA score calculation In [186]Figure 14A, the box plot revealed that the overall binding free energy (ΔG_bind) was predominantly negative, indicating a favorable binding interaction. The mean ΔG_bind value was −74.97 kcal/mol. Among the decomposed components, the van der Waals interaction (ΔG_bind_Vander) contributed the most negative value (−49.95 kcal/mol), followed by Coulombic interactions (−20.01 kcal/mol) and lipophilic contributions (−29.82 kcal/mol). Hydrogen bonding and covalent interactions showed minor contributions with mean values around −2.07 and −5.07 kcal/mol respectively. The relatively tight distribution for ΔG_bind_Vander suggested its consistency across simulation frames, whereas other terms exhibited higher variability. Figure 14. [187]Molecular dynamics simulation study of the complex between drug and receptor in energy decomposition. [188]Open in a new tab Presented the MM/GBSA-based post-MD binding free energy decomposition for the domperidone-7DFP complex. The analysis comprised 3 parts: a box plot of score components (A), time evolution of the total binding free energy (B), and frequency distributions of individual energy terms (C). [189]Figure 14B depicted the ΔG_bind time series across 100 ns of MD simulation. The ΔG_bind fluctuated stably around the mean value (−74.97 kcal/mol), as indicated by the red dashed line. Despite short-term fluctuations, the absence of drift confirmed the energetic stability of the complex throughout the trajectory. In [190]Figure 14C, histogram distributions of ΔG_bind and its decomposed terms provided further statistical insights. The ΔG_bind distribution exhibited a unimodal Gaussian-like shape centered near −74.97 kcal/mol. Similarly, Coulombic, covalent, van der Waals, hydrogen bonding, and lipophilic components all followed approximately normal distributions, indicating consistent sampling across the trajectory. The tightest distribution was observed for ΔG_bind_Vander, corroborating its dominant and stable contribution to binding. In [191]Figure 15A, the box plot highlighted the mean binding free energy (ΔG_bind) to be approximately −75.70 kcal/mol, suggesting a strong and energetically favorable interaction between terfenadine and 7DFP. The van der Waals interaction (ΔG_bind_Vander) was again the most substantial contributor, with mean energy close to −49.83 kcal/mol. Lipophilic interactions also played a significant role (−30.06 kcal/mol), followed by Coulombic (−19.57 kcal/mol), covalent (−4.92 kcal/mol), and hydrogen bonding interactions (−1.97 kcal/mol). All components showed stable distributions with limited outliers, indicating consistency across the MD trajectory. Figure 15. [192]Visualize the MM/GBSA score decomposition for the terfenadine–7DFP complex after MD simulation, including a box plot of energy components, ΔG_bind over time, and histograms of each component's frequency. [193]Open in a new tab Illustrated the MM/GBSA score decomposition analysis for the terfenadine–7DFP complex post molecular dynamics (MD) simulation. The figure comprised 3 panels: (A) box plot of binding score components, (B) ΔG_bind fluctuations over time, and (C) histograms of each energy component. [194]Figure 15B showed the time evolution of the total ΔG_bind over 100 ns. The binding free energy remained stably centered around −75.70 kcal/mol (indicated by the red dashed line), with periodic fluctuations that did not exhibit any upward or downward drift. This stability further confirmed the reliability of the interaction throughout the simulation. All distributions approximated normal Gaussian profiles, affirming statistical reliability and uniform sampling across the trajectory ([195]Figure 15C). The ΔG_bind values clustered tightly around the mean, reinforcing the notion of a well-converged binding profile. Among the components, van der Waals and lipophilic interactions showed the narrowest spread, consistent with their major contribution and stability. Discussion In this study, we explored the repurposing potential of terfenadine and domperidone against apoptotic targets relevant to CRC using an integrated computational approach. Molecular docking revealed that both compounds exhibited strong binding affinities across key proteins, with domperidone generally demonstrating either equivalent or superior interactions compared with terfenadine. A systems pharmacology approach was used to identify key apoptotic targets for CRC. Network pharmacology analysis revealed SLC6A4, DRD2, HTR2A, and EGFR as central hubs associated with apoptotic regulation, tumor survival, and metastasis. These targets were prioritized based on topological parameters such as degree, betweenness, and closeness centrality, highlighting their critical role in CRC pathogenesis. To elucidate the functional relevance of the 35 putative target genes shared between terfenadine, domperidone, and CRC, KEGG pathway enrichment analysis was performed ([196]Figure 2A). The results highlighted several significantly enriched pathways, with the top hits including gap junction, serotonergic synapse, synaptic vesicle cycle, and neuroactive ligand-receptor interaction. Additional enriched pathways included estrogen signaling, calcium signaling, CGMP-PKG (cyclic GMP-dependent protein kinase) signaling, Rap1 signaling, and cAMP signaling, all of which are known to play roles in cancer progression, cellular communication, and apoptotic regulation. The fold enrichment for these pathways ranged broadly, with gap junctions exhibiting the highest enrichment score, suggesting its prominent involvement in modulating intercellular communication and tumor suppression. The hierarchical clustering tree ([197]Supplementary file, Table S6 and [198]Figure S2) further illustrated the relationships among the enriched pathways based on shared gene content. Pathways such as estrogen signaling, cAMP signaling, and Rap1 signaling clustered closely, reflecting functional similarities in second messenger and hormone-related signaling mechanisms. Similarly, neuroactive ligand-receptor interaction, serotonergic synapse, and synaptic vesicle cycle were grouped, underscoring a convergence on neurotransmission and receptor-ligand regulatory systems. The positioning of gap junction and calcium signaling in close proximity also supports their interconnected roles in intracellular signal propagation and apoptosis. These results suggest that terfenadine and domperidone may influence CRC progression through modulation of neuronal signaling, hormone response, and intercellular communication pathways. These mechanistic insights provide a rationale for targeting neurotransmitter systems in CRC therapy and validate the multi-targeted potential of the repurposed drugs. The dopamine receptor D2 (DRD2), a GPCR primarily known for its role in the central nervous system, has emerged as a significant regulator in various cancers, including CRC. In CRC, DRD2 has been implicated in promoting tumor cell proliferation, survival, and metastasis.^ [199]37 CRC is characterized by a complex interplay of genetic mutations, epigenetic alterations, and dysregulated signaling pathways such as Wnt/β-catenin, PI3K/AKT, MAPK, and TGF-β, which collectively contribute to uncontrolled cell proliferation, evasion of apoptosis, angiogenesis, and metastasis. The loss of apoptotic control is a central hallmark of CRC, often mediated by overactivation of survival-promoting receptors and kinases.^ [200]38 In this study, we identified 4 key targets—DRD2, HTR2A (5-Hydroxytryptamine Receptor 2A), EGFR (Epidermal Growth Factor Receptor), and SLC6A4 (Serotonin Transporter)—which are deeply involved in these pathological mechanisms. Domperidone, a dopamine D2 receptor antagonist, binds strongly to DRD2, disrupting its G-protein-coupled signaling. This inhibition attenuates downstream PI3K/AKT and MAPK pathways, thereby restoring apoptotic sensitivity in CRC cells. Terfenadine, an antihistaminic agent, was shown to interact with HTR2A, a serotonin receptor implicated in mitogenic and anti-apoptotic signaling. By antagonizing HTR2A, terfenadine may reduce serotonin-mediated proliferation and enhance caspase-dependent apoptosis. In addition, both drugs showed interaction with EGFR, a receptor tyrosine kinase commonly overexpressed in CRC, leading to inhibition of RAS-RAF-MEK-ERK signaling—a pathway that drives tumorigenesis and chemoresistance.^ [201]39 The SLC6A4 transporter, which regulates serotonin uptake and modulates tumor microenvironment signaling, was another critical node targeted in this study; terfenadine’s interaction with SLC6A4 may hinder serotonin availability, further impairing survival pathways. By targeting these receptors and transporters, terfenadine and domperidone modulate multiple converging pathways involved in CRC pathogenesis, particularly those affecting apoptosis, cell cycle regulation, and neuroendocrine signaling. The integrative analyses demonstrated that both compounds can stably engage these targets under physiological conditions. Thus, this mechanistic overview supports the repositioning potential of domperidone and terfenadine as multi-target modulators capable of disrupting CRC progression through synergistic effects on key oncogenic and apoptotic circuits. The proteins showed distinct physiochemical properties, with variations in size, pI, and stability. While all had long half-lives in E. coli, 7DFP appeared less stable based on its higher instability index. Differences in GRAVY scores and aliphatic indices suggested varied solubility and potential thermostability among the proteins. Molecular docking results showed strong binding affinities of terfenadine and domperidone against these apoptotic targets. Against DRD2 (7DFP), both drugs achieved identical binding score of −11.5 kcal/mol. Notably, domperidone formed stable hydrogen bonds with CYS182 at 2 distances (2.4 Å and 3.3 Å), while terfenadine interacted with ILE184 (2.6 Å), suggesting strong and stable anchoring within the receptor’s active site. Binding analysis against SLC6A4 (5I6X) showed terfenadine had a slightly better affinity (−11.0 kcal/mol) than domperidone (−9.7 kcal/mol); however, domperidone maintained strong bonding with ASP98 (2.2 Å), a residue critical for SLC6A4 anti-apoptotic activity. Interestingly, domperidone outperformed terfenadine against HTR2A (6WGT), binding at −11.3 kcal/mol compared with −9.0 kcal/mol for terfenadine, indicating a higher potential to modulate HTR2A-mediated apoptosis. For EGFR (6LUD), domperidone again showed stronger interactions (−8.4 kcal/mol vs −7.1 kcal/mol) and established multiple hydrogen bonds, suggesting enhanced stabilization of p53-mediated tumor suppressor functions. Molecular dynamics simulations further validated the docking results by demonstrating that the complexes, particularly those with domperidone, exhibited minimal fluctuations in RMSD, stable hydrogen bond profiles, and compact rGyr values throughout the 100 ns simulation period. These findings indicated that both drugs, especially domperidone, could maintain stable and strong interactions within the active sites under physiological conditions. The MM/GBSA results indicated that the binding of terfenadine and domperidone with their respective protein was predominantly driven by van der Waals and hydrophobic forces, complemented by electrostatic and minor hydrogen bond contributions. Importantly, previously published in vitro studies also support the repositioning potential of these drugs. Domperidone has been reported to induce apoptosis in cancer cell lines by modulating dopamine receptor pathways and affecting downstream AKT signaling.^ [202]12 Similarly, studies have shown that terfenadine can inhibit cancer cell proliferation and promote apoptosis via oxidative stress and mitochondrial pathway modulation. Specifically, DRD2 antagonism has been linked to reduced cancer cell viability and increased apoptotic markers such as cleaved caspase-3 and Bax in CRC models, aligning with our target engagement findings.^ [203]10 BCL2 inhibitors have consistently shown synergy with standard chemotherapy in CRC cells by promoting intrinsic apoptotic pathways, while p53 activation remains a cornerstone for inducing tumor cell death.^ [204]40 Collectively, the integrated docking and MD simulation analyses provide strong evidence that domperidone possesses a superior binding profile across multiple apoptotic targets compared with terfenadine. This multi-target interaction potential could translate into better therapeutic outcomes by simultaneously modulating several apoptotic pathways critical in CRC progression. Thus, domperidone emerges as a particularly promising candidate for further preclinical investigation in CRC therapeutics. Limitation Although the in silico results are encouraging, this work is entirely computational in nature and does not have experimental or clinical confirmation. In vitro or in vivo investigations to confirm the hypothesized correlations are lacking, which limits the findings’ immediate translational impact. Building on the promising computational results, future research should focus on experimental validation of terfenadine and domperidone in CRC models. In vitro assays, such as apoptosis induction, cell viability, and target gene expression studies, should be conducted to confirm their mechanistic roles. In addition, in vivo studies are needed to assess therapeutic efficacy, pharmacodynamics, and toxicity profiles. Exploring combination therapies with existing chemotherapeutic agents could further enhance anticancer effects. Investigating the broader applicability of this repurposing strategy across different cancer types may also uncover new therapeutic opportunities. Conclusion This study demonstrated the potential of repurposing terfenadine and domperidone as therapeutic agents targeting apoptotic gene associations in CRC through a comprehensive systems pharmacology approach integrated with molecular docking and dynamics simulations. Network pharmacology analysis identified top targets involved in CRC as critical apoptotic regulators, with molecular docking results showing strong binding affinities, particularly against DRD2 at −11.5 kcal/mol for both drugs. The stability and reliability of these interactions were further substantiated by MD simulations, revealing sustained complex stability and favorable conformational dynamics over a 100 ns trajectory. These findings suggest that terfenadine and domperidone can effectively modulate apoptotic signaling pathways relevant to CRC progression. This multi-targeted repositioning strategy not only accelerates drug development timelines but also opens new avenues for enhancing therapeutic outcomes in CRC treatment. Supplemental Material sj-docx-1-bbi-10.1177_11779322251365019 – Supplemental material for Repurposing terfenadine and domperidone for inhibition of apoptotic gene association in colorectal cancer: A system pharmacology approach integrated with molecular docking, MD simulations, and post-MD simulation analysis [205]sj-docx-1-bbi-10.1177_11779322251365019.docx^ (711.2KB, docx) Supplemental material, sj-docx-1-bbi-10.1177_11779322251365019 for Repurposing terfenadine and domperidone for inhibition of apoptotic gene association in colorectal cancer: A system pharmacology approach integrated with molecular docking, MD simulations, and post-MD simulation analysis by Pushpaveni C, Hemavathi S, Santosh Prasad Chaudhary Kurmi, Biswa Ranjan Patra, V Angelin Esther, Chandrajeet Kumar Yadav, Mahalakshmi Suresha Biradar and Shankar Thapa in Bioinformatics and Biology Insights Acknowledgments