Network Pharmacology-Guided Identification of Secondary Metabolites as Modulators of HIV Latency via Oncogenic Signaling Pathways.

1. Introduction
Eradication of human immunodeficiency virus type 1 (HIV-1) remains an unsolved global health challenge, primarily due to the persistence of latent viral reservoirs that evade clearance by combination antiretroviral therapy (cART). These reservoirs, predominantly housed within resting CD4+ T cells, lymphoid tissues, gut-associated lymphoid tissue (GALT), and anatomical sanctuaries such as the central nervous system (CNS), persist in a transcriptionally silent state, rendering them refractory to immune surveillance and pharmacological intervention [1,2]. Notably, the CNS poses a unique obstacle due to limited drug penetrance, immune privilege, and the presence of specialized reservoir cell types, necessitating novel therapeutic strategies capable of targeting latency across compartmentalized reservoirs [3]. Current latency-modulating strategies are primarily centred on the “shock and kill” and “block and lock” paradigms. The former seeks to reactivate latent proviruses using latency-reversing agents (LRAs), such as histone deacetylase inhibitors (HDACis), protein kinase C (PKC) agonists, toll-like receptor (TLR) agonists, and bromodomain and extraterminal domain (BET) inhibitors, followed by immune-mediated or cytopathic clearance [4]. In contrast, “block and lock” approaches aim to suppress viral transcription, thereby locking HIV into deep latency [5]. Although clinical trials of LRAs such as romidepsin, bryostatin-1, and JQ1 have demonstrated mechanistic activity, their clinical success has been hindered by insufficient latency reversal, immunotoxicity, and heterogeneity of latent reservoirs [6,7]. This underscores the ongoing need for novel and safer LRAs, particularly those with microbiome-derived origin, which may offer unique therapeutic advantages in HIV-candidiasis co-infection.
Opportunistic infections are a hallmark of HIV disease progression, with Candida albicans emerging as one of the most prevalent fungal pathogens in people living with HIV (PLWH) [8]. C. albicans frequently colonizes mucosal surfaces and not only breaches epithelial barriers but also influences immune homeostasis by modulating host signaling pathways, such as NF-κB, IL-17, and MAPK signaling [9,10]. These immunomodulatory effects may influence the activation state of CD4+ T cells and cytokine environment, contributing to the establishment or reversal of HIV latency. Among the secreted metabolites, candidalysin and farnesol have been characterized for their immunological activity. Candidalysin is a pore-forming peptide toxin secreted during hyphal growth that disrupts epithelial membranes and activates innate immune signaling pathways, including the MAPK and NF-κB pathways. It also induces ligand shedding, primarily via epiregulin and epigen, which activates EGFR signaling and amplifies innate immune responses. Additionally, candidalysin activates the NLRP3 inflammasome in immune cells, such as macrophages and dendritic cells, leading to IL-1β secretion and promoting inflammatory cascades [11]. Farnesol, a sesquiterpene alcohol, functions as a quorum-sensing molecule that regulates fungal morphogenesis and virulence and modulates host immune signaling pathways. Farnesol activates innate immune cells by inducing pro-inflammatory signaling through the MAPK and NF-κB pathways, resulting in cytokine release and oxidative burst responses. At the same time, it has been reported to impair dendritic cell maturation by downregulating the granulocyte-macrophage colony-stimulating factor (GM-CSF) receptor, thereby reducing T-cell activation. Farnesol may also alter immune polarization by promoting inflammatory responses and suppressing Th1-mediated immunity [12]. Despite these known host–pathogen interactions, the potential of C. albicans secondary metabolites to directly target human proteins involved in HIV latency regulation has not been systematically explored, representing an untapped chemical space for latency-modulating therapeutic discovery. Importantly, fungal metabolites from other species have shown potential to modulate HIV latency. For instance, gliotoxin from Aspergillus fumigatus disrupts the 7SK small nuclear ribonucleoprotein (snRNP) complex, thereby releasing positive transcription elongation factor b (P-TEFb) and promoting HIV reactivation [13]. Chaetocin from Chaetomium minutum acts as a Suv39H1 histone methyltransferase inhibitor (HMTi) and has been shown to reverse HIV latency by demethylating repressive histone marks, thereby enhancing transcription at the HIV long terminal repeat (LTR) region [14]. Harziachalasin D, derived from the endophytic fungus Trichoderma harzianum MLJ-4, reverses HIV latency by activating the NF-κB pathway. This activation leads to the phosphorylation and degradation of inhibitor of nuclear factor kappa B alpha (IκBα), facilitating NF-κB nuclear translocation and HIV transcriptional activation [15]. These findings underscore fungal species as a largely unexplored source of bioactive small molecules with potential as LRAs.
Given the immunological relevance of C. albicans in HIV co-infection and the precedent set by other fungal metabolites, we hypothesized that certain C. albicans-derived secondary metabolites could modulate key host pathways involved in HIV-1 latency. To investigate this, we employed an integrative network pharmacology and in silico approach. Network pharmacology provides a holistic framework for characterizing the multi-target effects of bioactive compounds by integrating pharmacokinetic profiling, target prediction, protein–protein interaction (PPI) analysis, and pathway enrichment within the context of systems biology [16,17]. Although emerging in HIV cure research, this approach has begun to elucidate host dependency pathways and potential therapeutic targets involved in viral persistence, as well as the repurposing of drugs for HIV cure [18,19,20]. In this study, C. albicans metabolites were curated from the literature and screened for drug-likeness, oral bioavailability (OB), blood–brain barrier (BBB) permeability, and toxicity. Potential human protein targets were predicted and cross-referenced with the HIV latency-associated proteins. Subsequently, PPI network construction, topological analysis, and functional enrichment analyses were performed to reveal key biological processes and signaling pathways. Finally, molecular docking, molecular dynamics (MD) simulations, and binding free energy (BFE) calculations were performed to evaluate metabolite–target (MET–TAR) binding affinities.
Figure 1 shows the methodological framework used in this study. Our findings highlight that the enrichment of signaling cascades, while predominantly oncogenic, are critical host dependency factors for HIV-1 transcriptional regulation. To the best of our knowledge, this is the first comprehensive investigation on the role of C. albicans secondary metabolites as modulators of HIV-1 latency. These findings suggest a previously underexplored microbial metabolite space with potential relevance to HIV latency modulation.

2. Results

2.1. Drug-likeness and Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) Screening for Therapeutic Metabolites
Six candidate metabolites were identified as described in Section 4. All six metabolites complied with Lipinski’s Rule of Five (Ro5) with no observed violations, supporting their drug-like properties and OB. The results are presented in Table 1. Molecular weights (MWs) ranged from 222.37 g/mol (MET 181) to 411.43 g/mol (MET 119), and all compounds had favorable bioavailability score (BS) of 0.55 and 0.85 (MET 176). The topological polar surface area (TPSA) values ranged from 20.23 to 138.51 Å2. Gastrointestinal (GI) absorption was predicted to be high for five compounds, with MET 119 being the only exception (low absorption). Moreover, all compounds are non-substrates of P-glycoprotein (Pgp). Three metabolites (MET 28, MET 176, and MET 181) were predicted to cross the BBB. Cytochrome P450 3A4 (CYP3A4) inhibition was predicted for MET 28 and MET 119. Solubility predictions showed moderate solubility for most compounds, with MET 34 classified as very soluble and MET 119 as soluble. None of the compounds contain pan-assay interference compounds (PAINS) alerts, and the synthetic accessibility (SA) scores ranged from 2.51 to 3.26. The median lethal dose (LD50) values ranged from extremely low in MET 34 (41 mg/kg, Class II) and MET 176 (48 mg/kg, Class II) to exceptionally high in MET 119 (25,000 mg/kg, Class VI). LD50 values were moderately high for MET 15 (1500 mg/kg, Class IV) and MET 181 (5000 mg/kg, Class V). All the metabolites were inactive with respect to nephrotoxicity, cardiotoxicity, immunotoxicity, and cytotoxicity. However, only one metabolite flagged specific concerns based on high-confidence predictions (confidence score ≥ 0.70). MET 15 was predicted to be neurotoxic (active: 0.88). Other endpoints such as hepatotoxicity, mutagenicity, and respiratory toxicity were either inactive or predicted with low confidence (confidence score < 0.70). The predicted toxicity profiles of the metabolites are illustrated in the radar chart in Supplementary Figure S1.

2.2. Target Prediction and HIV Latency Relevance Mapping of Candida albicans Secondary Metabolites
A total of 369 non-redundant targets were identified for the six selected Candida albicans metabolites, spanning diverse functional classes, including enzymes, kinases, lyases, and proteases (Supplementary Figure S2). The predicted targets were 39, 127, 115, 49, 124, and 97 for MET 15, MET 28, MET 34, MET 119, MET 176, and MET 181, respectively. To assess their relevance in HIV latency, a curated list of 2510 human genes associated with HIV latency was obtained from the GeneCards database using the keyword “HIV latency.” A cross-comparison of the predicted targets with the latency-associated gene set using Venny 2.1.0 revealed 176 overlapping genes (Figure 2A). The PPI network comprised 176 nodes and 2464 edges, exhibiting a significantly enriched interaction pattern (PPI enrichment p < 1.0 × 10−16), far above what would be expected for a random set of proteins (expected edges = 1009). The average node degree was 28 and the average local clustering coefficient was 0.54 (Figure 2B).

2.3. Construction and Topological Analysis of MET–TAR Network
The MET–TAR network comprised 182 nodes (six metabolites and 176 target proteins) and 263 edges, representing the predicted interactions between each metabolite and its associated target(s) (Figure 3). Topological analysis was performed using the NetworkAnalyzer plugin v4.5.0 in Cytoscape v3.10.3 to assess the degree centrality (DC), betweenness centrality (BC), and closeness centrality (CC) of individual nodes, which reflects the number of direct connections of a node within the network, the relative extent to which a node is positioned on the shortest paths between other nodes, and the relative proximity of a node to all other nodes in the network, respectively. MET 28 recorded the highest DC of 69, followed by MET 34, MET 176, MET 181, MET 119, and MET 15 with DC of 65, 49, 42, 24, and 14, respectively. Similarly, MET 28 showed the highest BC (0.454) and CC (0.452), followed by MET 34 (BC = 0.427, CC = 0.444), MET 176 (BC = 0.305, CC = 0.411), MET 181 (BC = 0.229, CC = 0.399), MET 119 (BC = 0.093, CC = 0.369), and MET 15 (BC = 0.040, CC = 0.355) (Table 2). Among the protein targets, PTGS1 and PTGS2 both had a DC of 5, followed by AR with a DC of 4. Other highly connected proteins included GSK3B, CASP9, HDAC6, PIK3CA, TLR9, and HDAC4 with a DC of 3, as well as CDK2, PARP1, CDK4, MIF, EGFR, AKT1, HDAC8, ESR1, CDK5, SIRT1, RELA, JAK1, JAK2, PIK3CD, and MDM2 with a DC of two (Supplementary Table S2).

2.4. PPI Network Analysis and Hub Target Identification
The PPI network was imported into Cytoscape v3.10.3 for advanced visualization and topological analyses. The network had a diameter of 4, density of 0.164, characteristic path length of 1.988, and centralization value of 0.542. Notably, 174 nodes were part of a single connected component, suggesting robust molecular interplay among the targets. To prioritize the key regulatory targets within the network, CytoHubba v0.1 was used to calculate the DC score for each node (Figure 4A). The nodes were visualized using a red–yellow color gradient corresponding to a decreasing DC score and highlighted in a focused subnetwork (Figure 4B). The top 20 hub genes identified were AKT1, TNF, EGFR, STAT3, MAPK3, ESR1, GSK3B, MTOR, PPARG, PTGS2, MMP9, CCND1, ERBB2, MDM2, JAK2, MAPK1, PIK3CA, SIRT1, RELA, and PARP1. Furthermore, the co-expression relationships among these genes were analyzed using ProteomeHD on STRING database v12.0, with edge coloration in the network indicating the confidence level of gene–gene associations. Darker edges represent higher confidence scores, suggesting robust co-regulation or functional interdependence among several hub genes involved in HIV latency (Figure 4C). A detailed summary of the hub genes, including their Ensembl identifiers, target families, DC scores, GeneCards inferred functionality scores (GIFTS), and HIV latency relevance scores, is provided in Table 3.

2.5. Functional Annotation and Pathway Enrichment Analysis of Overlapping Targets
The top 30 enriched gene ontology (GO) terms in the biological process (BP), cellular component (CC), and molecular function (MF) categories, as well as the Kyoto encyclopedia of genes and genomes (KEGG) pathways, are illustrated in Figure 5A–C and Figure 6, respectively.

2.5.1. GO Term Enrichment Analysis
A total of 1000 GO BP terms were identified. Among these, response to organonitrogen compound (GO:0010243), cellular response to oxygen-containing compound (GO:1901701), and response to nitrogen compound (GO:1901698) emerged as the top three enriched processes, with false discovery rate (FDR) values as low as 1.72 × 10−45, 1.26 × 10−55, and 7.81 × 10−45 respectively, and fold enrichment (FE) scores > 8.0. Notably, these terms also had high gene counts (71–85 genes). Although regulation of biological quality (GO:0065008) and response to chemical (GO:0042221) had slightly lower FE (<4.0), they were associated with substantially larger gene sets (121 and 140 genes, respectively) (Figure 5A and Supplementary Table S3). CC category yielded 317 enriched terms. The top-ranking terms included caveola (GO:0005901), plasma membrane raft (GO:0044853), and integral component of synaptic membrane (GO:0099699), with FE scores of 14.21, 10.43, and 8.89, respectively. However, these terms included relatively small gene sets (10–11 genes). In contrast, broader components, such as vesicle (GO:0031982), nuclear lumen (GO:0031981), and nucleoplasm (GO:0005654), overlapped with high gene sets of 75, 73, and 71, albeit at lower FE scores of 2.17, 1.90, and 2.00, respectively (Figure 5B and Supplementary Table S4). Among the 546 MF terms, the most significantly enriched were various kinase-related activities, including protein tyrosine kinase activity (GO:0004713), transcription coregulator binding (GO:0001221), protein serine/threonine/tyrosine kinase activity (GO:0004712), and protein serine kinase activity (GO:0106310), with FE scores of 18.70, 15.08, 13.75, and 11.00, respectively. These high-FE terms were complemented by broad-spectrum functions, such as catalytic activity acting on a protein (GO:0140096), small molecule binding (GO:0036094), anion binding (GO:0043168), and nucleotide binding (GO:0000166), with lower FE scores of 4.16, 3.77, 3.83, and 4.02, but high gene counts of 83, 80, 78, and 74, respectively (Figure 5C and Supplementary Table S5).

2.5.2. KEGG Pathway Enrichment Analysis
Of the 217 KEGG pathways identified, VEGF signaling pathway (hsa04370) emerged as the most overrepresented, with an FE of 41.63, overlapping with 19 genes. Other strongly enriched pathways included cancer-related pathways, including non-small cell lung cancer (hsa05223; FE = 37.70), pancreatic cancer (hsa05212; FE = 35.72), acute myeloid leukemia (hsa05221; FE = 34.73), EGFR tyrosine kinase inhibitor resistance (hsa01521; FE = 34.36), prostate cancer (hsa05215; FE = 33.32), and central carbon metabolism in cancer (hsa05230; FE = 33.24). In contrast, broader signaling mechanisms, such as microRNAs in cancer (hsa05206; 30 genes), PI3K-Akt signaling (hsa04151; 42 genes), and pathways in cancer (hsa05200; 57 genes), demonstrated extensive gene participation but slightly lower FE scores of 24.09, 15.34, and 13.90, respectively (Figure 6B and Supplementary Table S6). A Sankey diagram was generated to visualize the interplay between prioritized genes and pathways (Figure 6A). The visualization revealed functional pleiotropy for several hub genes, including PIK3CA, MAPK3, MAPK1, AKT1, EGFR, MTOR, RELA, and STAT3, which mapped to multiple key pathways, such as pathways in cancer (hsa05200), PI3K-Akt signaling pathway (hsa04151), proteoglycans in cancer (hsa05205), microRNAs in cancer (hsa05206), and prostate cancer (hsa05215). Interestingly, TNF and PARP1, despite being high-ranking genes, did not map to any of the top 30 enriched KEGG pathways, suggesting a potentially distinct mechanism of action outside this oncogenic-rich canonical pathway enrichment. Furthermore, based on their extremely low FDR values, high FE scores, and dense number of associated overlapping and hub genes, pathways in cancer (hsa05200) and PI3K-Akt signaling pathway (hsa04151) emerged as the most significantly enriched pathways in our analysis. Notably, both pathways have established roles in the regulation of HIV latency (Figure 7). The enriched KEGG pathways from our analysis were systematically compared with those in published literature on the mechanisms underlying HIV latency. We identified ten pathways with established or emerging evidence of involvement in HIV latency (Table 4).

2.6. Topological Analysis of the Metabolite–Target–Pathway (MET–TAR–PATH) Network
To investigate the integrative relationships among Candida albicans metabolites, their associated biological targets, and enriched signaling pathways, a MET–TAR–PATH network was constructed using Cytoscape v3.10.3. The merged network comprised 56 nodes, including six metabolites, 20 hub gene targets, and 30 KEGG pathways interconnected by 301 edges (Figure 8). Topological analysis involved DC calculations to identify several key nodes. Among the metabolites, MET 28 exhibited the most extensive interactions, mapping to 13 hub targets. MET 34, MET 176, MET 119, MET 181, and MET 15 mapped to eight, six, five, three, and three hub genes, respectively. Among the hub targets, PIK3CA exhibited the highest DC (32), followed by MAPK3 (29), MAPK1 (29), AKT1 (29), EGFR (21), RELA (20), and MTOR (19). This makes PIK3CA the “central bridge” of the tripartite network. In contrast, less connected targets included SIRT1 (3), PPARG (2), PARP1 (2), and TNF (1). Among the enriched pathways, pathways in cancer (hsa05200) mapped to 17 hub targets, human cytomegalovirus infection (hsa05163), proteoglycans in cancer (hsa05205), PI3K-Akt signaling (hsa04151), microRNAs in cancer (hsa05206), and prostate cancer (hsa05215) each mapped to 12 hub targets. The lowest hub target-expressed pathways were cAMP signaling pathway (hsa04024), Rap1 signaling pathway (hsa04015), sphingolipid signaling pathway (hsa04071), and VEGF signaling pathway (hsa04370), each with five connections. (Supplementary Table S8). Collectively, these findings underscore the potential multitarget therapeutic relevance of C. albicans metabolites in addressing HIV latency through the coordinated modulation of interconnected signaling networks that govern both cell cycle/oncogenic processes and viral latency.

2.7. Binding Affinity Calculations
To evaluate the binding affinity of the prioritized Candida albicans secondary metabolites with key HIV latency-related targets, molecular docking simulations were performed for 16 hub proteins identified through integrated MET–TAR–PATH analyses. The docking scores (binding energies in kcal/mol) are summarized in Supplementary Table S7 and are visualized as a heatmap (Figure 9A). Docking scores below −5.0 kcal/mol are generally indicative of favorable interactions and stable binding conformations. Across all six metabolites (MET 15, MET 28, MET 34, MET 119, MET 176, and MET 181), the binding scores ranged from −11.6 to −4.1 kcal/mol, confirming appreciable interactions with multiple hubs. Notably, MET 15 consistently exhibited strong binding, particularly with MMP9 (−11.6 kcal/mol), MAPK1 (−10.9 kcal/mol), ERBB2 (−10.4 kcal/mol), and both MAPK3 and PIK3CA (−10.3 kcal/mol), in some cases approaching those of the reference ligands. MET 119 also demonstrated broad and favorable activity, with binding energies ≤ −9.0 kcal/mol against MMP9, ERBB2, MAPK3, JAK2, MAPK1, GSK3B, MTOR, and PIK3CA. MET 28 displayed selective high-affinity interactions, particularly with ERBB2 (−9.7 kcal/mol) and MAPK3 (−9.0 kcal/mol), whereas MET 181 exhibited moderate but consistent interactions across several targets (−8.1 to −4.6 kcal/mol). In contrast, MET 34 and MET 176 generally showed weaker docking scores (−7.6 to −4.1 kcal/mol). The reference ligands displayed the expected high binding affinities, including SCH772984 for MAPK3 (−14.3 kcal/mol), RLY2608 for PIK3CA (−13.7 kcal/mol), VTX-11e for MAPK1 (−12.2 kcal/mol), Torin2 for MTOR (−11.6 kcal/mol), SYR127063 for ERBB2 (−11.4 kcal/mol), and BE4 for MMP9 (−10.6 kcal/mol). Although the reference ligands consistently outperformed the metabolites, except for MMP9 (−10.6 kcal/mol), several metabolites, particularly MET15, MET28, and MET119, achieved competitive docking scores. Collectively, MET 15 and MET 119 are multi-target binders because of their engagement with several high-priority latency-related hubs, whereas MET 28 and MET 181 retain selective potential. In contrast, MET 34 and MET 176 show weaker binding affinities. Hub targets with higher binding affinities and topological importance in the MET–TAR–PATH network were considered as key targets. The comparative binding poses of the selected complexes are shown in Figure 9B.

2.8. Protein–Ligand Interaction Analysis
To delineate the molecular underpinnings of ligand binding, detailed interaction profiling was performed (Figure 10). In the MAPK1MET 15 complex, hydrophobic interactions were dominant, involving ILE-31, ALA-35, TYR-36, VAL-39, ALA-52, and LEU-156. In comparison, VTX-11e shared hydrophobic contacts (ILE-31, VAL-39, ALA-52, ILE-56, and LEU-156) but extended the binding profile with two hydrogen bonds to MET-108 and a π-cation interaction with LYS-54 (Figure 10A). For PIK3CA binding, MET 15 engaged the hydrophobic residues PHE-933, LYS-937, PHE-941, TYR-943, ILE-1015, and ILE-1018, supported by hydrogen bonding to LYS-937 and π-stacking with PHE-998. RLY2608 displayed a more complex profile with hydrophobic contacts (VAL-948, TYR-1017, and ILE-1018), hydrogen bonds (LEU-907 and ASP-1014), π-stacking (PHE-998 and TYR-1017), a π-cation interaction with LYS-937, and halogen bonds (GLN-805, ILE-906, and GLU-1008) (Figure 10B). For MAPK3 binding, MET 15 formed hydrophobic contacts with ILE-48, TYR-53, VAL-56, LYS-71, ILE-73, GLU-88, and ASP-184, complemented by a hydrogen bond with ASP-184, π-stacking with TYR-53, and a π-cation interaction with LYS-71. SCH772984, in contrast, exhibited a more extensive network comprising hydrophobic interactions (TYR-53, ALA-69, LYS-71, TYR-81, GLN-122, and LEU-173), multiple hydrogen bonds (LYS-71, MET-125, and LYS-131), π-stacking with TYR-81, and salt bridges with GLU-88 and ASP-184 (Figure 10C). For EGFR binding, MET 15 interacted through hydrophobic contacts with LEU-23, LYS-50, and THR-95 and established hydrogen bonds with THR-95, MET-98, and GLY-101. Erlotinib displayed an expanded hydrophobic network (LEU-23, LYS-50, LEU-93, THR-95, and LEU-97) while retaining a hydrogen bond with MET-98, suggesting a broader binding stability (Figure 10D). Furthermore, in the MTORMET 15 complex, there were extensive hydrophobic contacts involving LEU-801, GLU-806, TYR-841, ILE-853, ILE-972, ASP-973, and TRP-1045, with two additional π-cation interactions contributed by LYS-803. Torin 2 displayed a comparable but distinct profile, involving hydrophobic residues (LEU-801, ILE-853, THR-861, ILE-972), hydrogen bonds with GLU-806 and VAL-856, and two π-stacking interactions with TRP-855 (Figure 10E). Finally, for AKT1 binding, MET 15 engaged in hydrophobic interactions with PHE-21, VAL-24, LEU-41, ILE-46, and PHE-298, supported by a π-stacking interaction with PHE-21. In contrast, EX4 formed a more diverse interaction profile, comprising hydrophobic contacts (LEU-16, VAL-24, GLU-138), hydrogen bonds (GLU-88, ALA-90), and a salt bridge with GLU-94 (Figure 10F).

2.9. Time-Resolved Thermodynamic Assessment: Structural Conformational Changes and BFE Thermodynamic Profiling

2.9.1. Structural Conformational Changes of Hub Targets
Backbone root mean square deviation (RMSD) and per-residue root mean square fluctuation (RMSF) analyses were used to quantify the conformational stability of MAPK1, PIK3CA, MAPK3, EGFR, MTOR, and AKT1 in their apo, MET 15-bound, and reference inhibitor–bound states over the 200 ns MD simulations (Figure 11). The average RMSD values for MAPK1, MAPK1MET 15, and MAPK1VTX−11e were 1.47 Å, 1.56 Å, and 1.53 Å, respectively. All three systems equilibrated rapidly within the initial simulation window (~40 ns) and remained stable thereafter. MAPK1MET 15 exhibited a compact RMSD distribution comparable to that of MAPK1VTX−11e. The average RMSF values for MAPK1, MAPK1MET 15, and MAPK1VTX−11e were 0.94 Å, 0.82 Å, and 0.83 Å, respectively. The RMSF profiles of the bound systems showed reduced fluctuations across most residues relative to the apo system, particularly outside the terminal regions (Figure 11A). The average RMSD values for PIK3CA, PIK3CAMET 15, and PIK3CARLY2608 was 2.08 Å, 2.08 Å, and 2.12 Å, respectively. The RMSD trajectories of all the complexes overlapped closely throughout the simulation. The average RMSF values of PIK3CA, PIK3CAMET 15, and PIK3CARLY2608 were 1.27 Å, 1.16 Å, and 1.25 Å, respectively. Moderate fluctuations were localized to the loop regions, with PIK3CAMET 15 exhibiting slightly restricted residue mobility compared to the apo and PIK3CARLY2608 systems (Figure 11B). The average RMSD of MAPK3, MAPK3MET 15, and MAPK3SCH772984 was 2.57 Å, 1.78 Å, and 2.18 Å, respectively. MAPK3MET 15 showed consistently lower RMSD trajectories than both the apo and MAPK3SCH772984 systems throughout the simulation. The average RMSF values for MAPK3, MAPK3MET 15, and MAPK3SCH772984 were 1.01 Å, 1.01 Å, and 1.13 Å, respectively. The RMSF profiles were generally uniform among the three systems, with apo maintaining slightly restrained flexibility across the catalytic domain (Figure 11C). The average RMSD values for EGFR, EGFRMET15, and EGFRErlotinib were 3.08 Å, 2.93 Å, and 5.42 Å, respectively. The average RMSF values for EGFR, EGFRMET15, and EGFRErlotinib were 1.72 Å, 1.40 Å, and 2.93 Å, respectively. EGFRErlotinib exhibited pronounced RMSD drift and elevated RMSF fluctuations, particularly toward the C-terminal region. In contrast, EGFRMET15 showed substantially lower RMSD and reduced RMSF fluctuations, closely mimicking the stability of the apo. MET 15 demonstrated superior conformational restraint compared to the clinical standard (Figure 11D). The average RMSD values for MTOR, MTORMET 15, and MTORTorin2 were 3.65 Å, 3.87 Å, and 2.98 Å, respectively. MTORTorin2 maintained the lowest RMSD trajectories, while MTORMET 15 showed moderate deviation but stable trajectories without abrupt conformational transitions, except for the first 40 ns. The average RMSF values for MTOR, MTORMET 15, and MTORTorin2 were 1.63 Å, 1.77 Å, and 1.75 Å, respectively. Comparable fluctuation patterns were observed across the three systems, with localized peaks corresponding to the flexible loop regions (Figure 11E). The average RMSD values for AKT1, AKT1MET 15, and AKT1EX4 were 2.02 Å, 2.08 Å, and 2.31 Å, respectively. AKT1MET 15 and apo exhibited nearly identical RMSD trajectories and remained stable throughout the simulation. In contrast, AKT1EX4 exhibited profound increased RMSD drifts, particularly from 180 to 200 ns. The average RMSF values for AKT1, AKT1MET 15, and AKT1EX4 were 0.93 Å, 1.07 Å, and 1.12 Å, respectively. AKT1EX4 showed elevated RMSF peaks, particularly at the C-terminus region, compared to AKT1MET 15 and apo (Figure 11F).

2.9.2. BFE Thermodynamic Profile of Hub Targets
BFE calculations quantified the thermodynamic stability of MET 15 in complex with MAPK1, PIK3CA, MAPK3, EGFR, MTOR, and AKT1 compared to their respective reference inhibitors (Table 5). ΔGbind was decomposed into ΔGgas and ΔGsolv. For MAPK1 binding, MET 15 exhibited a ΔGbind of −33.66 ± 0.16 kcal/mol, driven by favorable ΔEvdW (−36.95 ± 0.11 kcal/mol) and ΔEelec (−43.08 ± 0.19 kcal/mol) interactions. VTX-11e showed a stronger ΔGbind (−56.08 ± 0.18 kcal/mol), associated with a more favorable ΔGgas (−96.90 ± 0.35 kcal/mol) despite a significant ∆Gsolv (40.82 ± 0.21 kcal/mol). In the PIK3CA system, MET 15 achieved ΔGbind of −39.03 ± 0.14 kcal/mol, with dominant contributions from ΔEvdW (−50.54 ± 0.08 kcal/mol). In contrast, RLY2608 displayed a substantially stronger ΔGbind of −69.32 ± 0.10 kcal/mol, reflecting the enhanced ΔGgas (−110.69 ± 0.14 kcal/mol). For MAPK3 binding, MET 15 produced a ΔGbind of −25.61 ± 0.14 kcal/mol. Although both ΔEvdW and ΔEelec were favorable, the ΔGbind was weaker than that of SCH772984, which showed a markedly stronger ΔGbind of −88.34 ± 0.14 kcal/mol. In the EGFR complex, MET 15 exhibited a ΔGbind of −31.45 ± 0.10 kcal/mol, primarily supported by ΔEvdW (−44.44 ± 0.08 kcal/mol). Erlotinib demonstrated stronger ΔGbind (−41.90 ± 0.13 kcal/mol), accompanied by higher ΔEelec (−21.40 ± 0.18 kcal/mol) and increased ∆Gsolv (31.51 ± 0.15 kcal/mol). For MTOR binding, MET 15 recorded a weaker ΔGbind of −21.01 ± 0.16 kcal/mol, with moderate ΔGgas (−45.43 ± 0.26 kcal/mol). Torin2 exhibited a stronger ΔGbind of −34.94 ± 0.11 kcal/mol, consistent with stronger ΔEvdW (−44.73 ± 0.10 kcal/mol) and ΔEelec (−28.16 ± 0.24 kcal/mol). In the AKT1 system, MET 15 displayed the weakest binding among the MET 15 complexes, with ΔGbind of −19.72 ± 0.11 kcal/mol. In contrast, EX4 showed a substantially more favorable ΔGbind of −51.41 ± 0.13 kcal/mol, dominated by a very large ΔEelec (−279.81 ± 0.55 kcal/mol) offset by a high ∆Gsolv (284.07 ± 0.56 kcal/mol).

3. Discussion
This study investigated the therapeutic potential of Candida albicans secondary metabolites in modulating HIV pathogenesis through latency dynamics, extending beyond well-known molecules like candidalysin and farnesol to identify novel bioactive scaffolds. Although C. albicans has been extensively characterized as an opportunistic pathogen in HIV/AIDS, its metabolites remain largely underexplored as therapeutic agents, despite the emerging interest in fungal metabolite profiling [21,22,23,24,25]. For example, the quorum-sensing secondary metabolite farnesol, which C. albicans utilizes in vivo to regulate virulence via immunomodulation, has attracted pharmacological attention for its antimicrobial, anticancer, anti-obesity, anti-inflammatory, and antioxidant properties [26,27]. Notably, Kadian et al. [28] demonstrated that trans, trans-farnesol (TF) alleviated cognitive deficits, oxidative stress, mitochondrial dysfunction, and neuroinflammation in ICV-STZ-induced Alzheimer’s-like rat model. These neuroprotective effects were attributed to TF’s ability to restore antioxidant balance, suppress acetylcholinesterase and TNF-α levels, and normalize mitochondrial enzyme activities, including NADH dehydrogenase and succinate dehydrogenase [28].
In the present study, we employed computational screening, network pharmacology, and molecular modeling to identify C. albicans metabolites and their host targets that are implicated in HIV latency regulation. Six candidate metabolites were predicted to exhibit enhanced intracellular permeability and a high synthetic feasibility. Among them, three metabolites (MET 28, MET 176, and MET 181) demonstrated BBB permeability, which is pharmacologically advantageous given the role of CNS reservoirs in HIV persistence [2]. However, two metabolites (MET 28 and MET 119) may present drug–drug interaction (DDI) liabilities when co-administered with antiretroviral therapies metabolized by CYP3A4. MET 119 presents an “Absorption Paradox.” It has low GI absorption but is also soluble. This suggests that MET 119 may be better suited for intravenous (IV) administration rather than oral administration, which is common for intensive HIV reservoir-purging therapies. Toxicological profiling revealed generally favorable safety characteristics, with minimal organ and systemic risks observed, except for mild neurotoxicity alert associated with MET 15. Therefore, targeted delivery or structural optimization strategies can help reduce its toxicity while maintaining efficacy. Overall, the ADMET properties indicated acceptable solubility, drug-likeness, and manageable safety margins for most metabolites. Notably, BBB permeability and neurotoxicity predictions enabled the identification of candidate compounds with the potential to safely penetrate anatomical sanctuaries while maintaining systemic activity against peripheral HIV reservoirs. Apart from farnesol (MET 181), the remaining compounds were isolated from methanolic extracts of C. albicans cultures and possess diverse structural scaffolds that merit further exploration. The methanolic extract was found to exhibit broad-spectrum antibacterial potential [21]. These findings suggest that the metabolite repertoire of C. albicans could be significantly expanded through solvent variation and optimized culture conditions, consistent with the one strain, many compounds (OSMAC) approach successfully employed in endophytic fungi for the discovery of novel bioactive molecules [29]. Furthermore, high-throughput metabolomics strategies could further unlock the “dark metabolome” of C. albicans for HIV research. The metabolites were mapped to 369 human targets, of which 176 overlapped with HIV latency-associated genes. These overlapping targets represent potential mediators through which C. albicans metabolites influence viral persistence. Network topology analysis revealed dense interconnectivity, reflecting the coordinated regulation among latency-related targets. This “small-world network” implies that information (or signals) can spread very quickly across the network, which is biologically significant for HIV latency regulation. The comparatively high DC values observed for these metabolites suggest biological significance, as their extensive interactions with multiple target proteins imply a potential role in modulating host–hijacked HIV–signaling dynamics and network stability.
The 20 hub genes identified in the PPI network represent central regulators of cellular signaling and transcriptional control that are intricately linked to HIV latency maintenance and reactivation. They include kinases (AKT1, EGFR, MAPK3, GSK3B, MTOR, ERBB2, JAK2, MAPK1, and PIK3CA), enzymes (PTGS2, MMP9, MDM2, and PARP1), transcription factors (STAT3 and RELA), nuclear receptors (ESR1 and PPARG), epigenetic regulators (SIRT1), and signaling molecules (TNF and CCND1). Their high DC underscores their pivotal role as network regulators that can potentially govern HIV latency modulation. The polypharmacological interaction profile of the metabolites further supports the hypothesis that C. albicans secondary metabolites may exert synergistic effects by simultaneously targeting multiple latency-associated proteins. Among the kinases, AKT1 and MTOR serve as pivotal mediators of cellular survival and metabolic homeostasis, influencing both HIV transcription and reservoir viability. Modulation of the AKT/MTOR axis has been shown to alter latent-reservoir stability and reactivation dynamics. Similarly, MAPK1 and MAPK3 (ERK2 and ERK1) regulate transcriptional programs essential for proviral activation, as ERK-mediated signaling modulates activator protein 1 (AP-1) and positive transcription elongation factor b (P-TEFb) activities, thereby affecting chromatin accessibility and elongation of integrated proviruses [30]. Growth factor–related kinases, including EGFR, ERBB2, PIK3CA, and GSK3B, converge on PI3K/AKT signaling nodes that define T-cell activation thresholds and NF-κB signaling balance, which are key determinants of HIV proviral transcription [30,31,32]. Similarly, the JAK/STAT axis influences cytokine-mediated control of latency, and inhibition of JAK2 represses viral reactivation, whereas hyperactivation of STAT3 supports immune dysfunction and the persistence of latent infection [33,34]. Enzymatic hubs, such as MMP9, PTGS2, MDM2, and PARP1, integrate extracellular remodeling, inflammatory signaling, and DNA damage responses. Notably, PARP1 functions as a transcriptional co-activator of NF-κB and interacts with heat shock factor 1 (HSF1) to promote latent proviral reactivation and viral gene expression [35,36]. MDM2, a negative regulator of p53, modulates cell survival and may indirectly sustain latent reservoirs [37]. Transcriptional and epigenetic regulators RELA (NF-κB p65) and SIRT1 directly orchestrate proviral latency control, with SIRT1 deacetylating RELA at Lys310 to dampen NF-κB transcriptional activity and maintain latency [38]. Conversely, TNF serves as a potent latency-reversing cytokine by inducing NF-κB–mediated proviral transcription [39]. Nuclear receptors ESR1 and PPARG modulate host transcriptional and inflammatory pathways that shape HIV transcriptional competence [40,41]. Finally, CCND1, which governs the G1/S phase transition via CDK4/6 activation, remains low in quiescent CD4+ T cells, the major latent HIV reservoir, thereby reinforcing a non-permissive state for viral transcription [42,43].
GO enrichment highlights functional themes that mechanistically link Candida albicans metabolites to HIV latency control. Within the BP category, the enrichment of responses to oxygen- and nitrogen-containing compounds underscores the pivotal role of oxidative stress in NF-κB–dependent HIV reactivation [44,45]. CC enrichment of plasma membrane microdomains, specifically caveolae and lipid rafts, is particularly notable, as these dynamic regions orchestrate signaling assemblies, such as EGFR/PI3K/AKT and MAPK/ERK complexes, that drive transcriptional reprogramming essential for proviral activation and reservoir maintenance [46]. Moreover, HIV-1 exploits lipid rafts during viral entry, assembly, and budding, thereby co-opting the host membrane architecture for persistence [47]. In the MF category, significant enrichment for protein tyrosine kinase and serine/threonine kinase activities reflects the dominance of MAPK and PI3K/AKT/MTOR signaling hubs identified within the PPI network. ERK1/2 (MAPK3/MAPK1) signaling regulates Tat-dependent HIV transcription [48], while MTOR signaling modulates metabolic and translational landscapes that determine reservoir inducibility [49]. Moreover, the enrichment of transcription co-regulator binding and catalytic activity further corroborates the roles of SIRT1, PARP1, and NF-κB as critical regulators of proviral chromatin and transcriptional competence. SIRT1 deacetylates both Tat and RELA, thereby fine-tuning their transactivation potential [38,50], whereas PARP1 contributes to transcriptional regulation and latency dynamics through its co-activator functions [35]. NF-κB also cooperates with AP-1 downstream of T-cell receptor (TCR) signaling to trigger latent proviral reactivation [51].
KEGG pathway enrichment analysis highlights that Candida albicans secondary metabolites are predicted to target highly conserved growth and pro-survival signaling circuits that overlap extensively with mechanisms governing HIV latency and reactivation. Among the enriched pathways, the VEGF signaling pathway (hsa04370) was the most overrepresented, underscoring the convergence of angiogenic and viral regulatory networks. VEGF-mediated activation of the PI3K/AKT and MAPK cascades is known to feed into the NF-κB and HIF-1α pathways, both of which regulate HIV-1 transcriptional activation and latent reservoir maintenance [52]. The enrichment of several cancer-related pathways, including non-small cell lung cancer (hsa05223), pancreatic cancer (hsa05212), acute myeloid leukemia (hsa05221), and prostate cancer (hsa05215), further emphasizes the mechanistic parallels between tumorigenic signaling and HIV persistence [53]. These networks converge on MAPK/ERK and PI3K/AKT signaling hubs that fine-tune T-cell activation thresholds, chromatin accessibility, and transcriptional elongation of latent proviruses [54]. Consistent with this, therapeutic modulation of the PI3K/AKT axis has been proposed as a strategy to destabilize latent HIV reservoirs [30]. Ten of the enriched pathways identified in this study have previously been linked to HIV latency regulation, reactivation, or immune evasion. The PI3K–Akt signaling pathway (hsa04151) and pathways in cancer (hsa05200) exhibit the highest gene participation, reinforcing the interconnectedness of oncogenic and viral persistence mechanisms [34,38,54,55]. The Ras signaling (hsa04014) and focal adhesion (hsa04510) pathways interact with the MAPK/ERK cascade to promote transcriptional activation and intercellular viral transmission. Focal adhesions facilitate LFA-1–ICAM–mediated viral spread and contribute to reservoir establishment [56], whereas Ras signaling activates the Raf/MEK/ERK cascade, enhancing HIV transcription through RBF-2–dependent regulation of the LTR [54]. PD-L1 expression and PD-1 checkpoint signaling (hsa05235) also featured prominently, reflecting their central role in maintaining latency via T-cell exhaustion and functional anergy. The engagement of PD-1 suppresses TCR signaling and restricts proviral reactivation, whereas PD-1 blockade enhances latency reversal ex vivo [57,58]. The HIF-1 signaling pathway (hsa04066), downstream of PI3K/AKT and MAPK activation, modulates HIV transcription through metabolic reprogramming under hypoxic conditions: HIF-1α enhances LTR activity, while HIF-2α suppresses viral transcription within hypoxic lymphoid niches [59,60]. Calcium (hsa04020) and cAMP (hsa04024) signaling pathways are canonical regulators of latency reversal. Calcium influx activates the calcineurin–NFAT axis, while elevated cAMP stimulates the PKA/CREB cascade, both of which enhance proviral transcription and have been leveraged in “shock and kill” latency-reversal strategies [61,62]. The enrichment of microRNAs in cancer (hsa05206) and sphingolipid signaling (hsa04071) further underscores the epigenetic and lipid-mediated modulation of HIV persistence, respectively. MicroRNAs such as miR-28, miR-125b, and miR-150 suppress HIV transcripts in resting CD4+ T cells [63], while sphingolipid metabolites such as sphingosine-1-phosphate (S1P) activate NF-κB to promote reactivation, and analogs such as FTY720 suppress viral replication [64,65]. Sphingolipids, including sphingosine and sphinganine, are essential structural components of cellular membranes and key intermediates in lipid metabolism pathways. Together with cholesterol, they contribute to the formation of lipid raft microdomains that organize membrane signaling and host–pathogen interactions [66]. Sphingolipid metabolism plays an important role in viral infection by influencing viral attachment, membrane fusion, intracellular trafficking and budding. Several viruses, including HIV, influenza virus, rhinovirus, and SARS-CoV-2, exploit sphingolipid-rich membrane domains to facilitate cellular entry and replication [67]. Mechanistically, sphingosine and sphinganine alter membrane lipid organization and lipid raft stability, affecting virus–receptor interactions and membrane fusion. In addition, sphingosine serves as a precursor for S1P, a signaling lipid that regulates pathways involved in endocytosis, inflammation, and cytoskeletal dynamics, which can influence viral entry and replication [68]. Collectively, these findings suggest that C. albicans metabolites may influence HIV latency through multilayered mechanisms, spanning transcriptional regulation, metabolic control, and immune modulation, highlighting their polypharmacological potential for targeting cellular programs that sustain viral persistence. The enrichment of these cancer-related pathways reflects the signaling convergence within multifunctional host pathways that HIV exploits rather than disease-specific oncogenic mechanisms.
The integrated MET–TAR–PATH network delineates a central cell growth and survival signaling module anchored by PIK3CA, MAPK1, MAPK3, AKT1, EGFR, and MTOR, linking Candida albicans metabolites to multiple HIV latency-relevant pathways. The prominence of these hubs is mechanistically significant, as PI3K/AKT/MTOR and MAPK/ERK cascades are master regulators of cellular metabolism, transcriptional reprogramming, and survival processes that collectively govern proviral maintenance. While these hub proteins are often described as oncogenic in cancer contexts, they function as key metabolic and signaling regulators in quiescent T-cells. Their modulation in this setting is not aimed at oncogenesis but rather at controlling pathways that influence HIV latency, effectively “flipping the switch” of the latent provirus. Metabolites such as MET 28 appear to exert polypharmacological, multi-node modulation capable of altering the latency–reactivation equilibrium. Lower-degree nodes, including SIRT1, PPARG, PARP1, and TNF, likely represent context-dependent regulatory axes encompassing chromatin remodeling, DNA repair, and cytokine-mediated inflammatory responses, processes that complement but are less central to the canonical KEGG-defined latency framework. The strong multitarget binding affinities observed for MET 15, MET 28, and MET 119 further reinforce their potential as broad-spectrum modulators of latency. Their stable engagement with MAPK1/3, PIK3CA, AKT1, MTOR, and EGFR implies coordinated regulation of both MAPK/ERK and PI3K/AKT/MTOR signaling—two convergent axes that dictate the transcriptional and survival landscape of HIV-infected cells. Sustained PI3K/AKT/MTOR activation maintains latency by promoting cellular survival and repressing viral transcription, whereas MAPK/ERK signaling mediates the transcriptional reactivation of latent proviruses [30,48,49]. Thus, the engagement of these kinases by C. albicans metabolites could destabilize latency-supportive signaling, shifting the balance toward proviral reactivation or immune clearance. Notably, the predominance of hydrophobic and π–stacking interactions observed in molecular docking mirrors those exploited by the optimized kinase inhibitors, suggesting favorable druggability and bioavailability profiles. Collectively, the polypharmacological binding patterns of key metabolites, particularly MET 15 and MET 119, underscore their potential to modulate multiple latency-supportive networks concurrently, thereby providing a rational scaffold for the development of next-generation latency-modulating therapeutics.
MD simulations and MM/GBSA analyses provided thermodynamic support for the network-level predictions by demonstrating that MET 15 forms stable complexes with multiple kinase hubs without inducing large-scale structural perturbations. In silico studies of kinase–ligand systems have shown that the preservation of native conformational dynamics is a desirable feature for compounds intended to modulate host signaling pathways rather than fully abrogating kinase function [69,70]. The dominance of van der Waals contributions in MET 15 binding, accompanied by moderate electrostatic interactions, is consistent with the binding modes reported for several ATP-competitive kinase inhibitors, in which hydrophobic packing within conserved kinase pockets represents a major energetic driver. This stable binding attests to the excellent shape complementarity of MET 15 with the hydrophobic pockets of these kinases. Importantly, MET 15 exhibited moderate but consistent binding across multiple hubs, supporting its classification as a polypharmacological modulator rather than a high-affinity single-target inhibitor. From a biological perspective, the simultaneous engagement of the PI3K/AKT/MTOR and MAPK/ERK axes is notable, given their documented roles in regulating cellular survival and inducible HIV transcription. Latency is increasingly understood as a network-stabilized state rather than a single-pathway phenomenon, suggesting that moderate multi-node perturbation may be more effective than selective inhibition [38,54]. The key hub genes (PIK3CA, MAPK1/3, EGFR, MTOR, and AKT1) are widely recognized signaling hubs and may appear as high-degree nodes in PPI networks due to study and annotation bias [71,72]. To address this challenge, their prioritization in this study was not based solely on network centrality but also on their intersection with HIV latency-associated host factors, convergence across enriched KEGG pathways, and subsequent structural validation through molecular docking and MD simulation [20]. Collectively, the systems-level coherence strengthens the biological plausibility of fungal secondary metabolites as unconventional but mechanistically useful scaffolds for host-directed HIV latency intervention.
The candidate metabolites identified in this study represent structurally diverse chemical classes, including nitrogen-rich heterocycles, quinoline derivatives, sulfonamides, fatty acids, and terpenoid alcohols. Many compounds containing these pharmacophores have been associated with antiviral activity or with modulation of host signaling pathways that regulate viral latency. For example, MET 15 contains a fused nitrogen-rich heterocyclic system combining pyridazine and pyrrolo motifs. Pyrrolo-pyrimidine analogs have demonstrated antiviral activity against flaviviruses, such as Zika and dengue viruses [73]. In addition, fused nitrogen heterocycles share structural similarities with pharmacophores found in bromodomain and extraterminal domain (BET) inhibitors, which influence HIV transcription through pathways involving BRD4 and P-TEFb [74]. MET 28, a quinoline derivative, belongs to a class of compounds with broad antiviral properties, including activity against the dengue virus, SARS-CoV-2, and HIV [75]. Quinoline-based compounds, such as AV6 derivatives, have also been shown to reactivate latent HIV in cellular models and to synergize with histone deacetylase (HDAC) inhibitors, thereby enhancing viral transcription without inducing excessive T-cell activation [76]. MET 34 contains a pyrrolidine ring and tertiary amide functional groups, structural motifs commonly found in enzyme-binding ligands. Pyrrolidine-containing compounds have been reported to interact with viral enzymes, including proteases and polymerases [77]. MET 119 combines pyrimidine and sulfonamide moieties. Pyrimidine derivatives can interfere with nucleic acid synthesis [78], while sulfonamide-containing molecules interact with epigenetic enzymes such as HDACs [79]. MET 176 is a long-chain unsaturated fatty acid, a class of molecules that can inhibit HIV entry or fusion with host cells [80]. Fatty acids may also modulate host metabolic signaling pathways through nuclear receptors, such as peroxisome proliferator-activated receptors (PPARs), which may influence viral persistence [81]. Finally, MET 181, a sesquiterpene alcohol, has been reported to modulate host inflammatory signaling pathways, such as NF-κB, which plays a key role in HIV transcriptional activation [82]. Collectively, the structural features and previously reported biological activities of these compounds provide a plausible basis for their predicted interactions with host targets involved in HIV latency, thereby supporting their prioritization for further experimental investigation. The identification of latency-modulating metabolites from Candida albicans should not be interpreted as implying that fungal colonization or infection provides therapeutic benefits for PLWH. Rather, consistent with established natural-product drug discovery approaches, microorganisms, including opportunistic and/or pathogenic species, produce chemically diverse secondary metabolites for pharmacological development following isolation, structural optimization, and experimental validation [83,84].

Study Limitations
This study offers an integrative systems-level perspective on the potential of Candida albicans metabolites to modulate HIV-1 latency; however, several limitations should be acknowledged. These findings are based on computational predictions from network pharmacology, molecular docking, and dynamics simulations, which, while robust, remain theoretical and require experimental confirmation. In silico approaches cannot fully capture the multifactorial complexity of HIV latency regulation, including chromatin remodeling, transcriptional feedback, and immune microenvironmental influences. PPI network outcomes are constrained by existing database annotations and potential knowledge biases. Additionally, pathway enrichment analyses are inherently constrained by database structures and annotation biases. The apparent over-representation of molecular targets within the KEGG pathways in cancer map (hsa05200), for example, largely reflects the aggregation of extensively characterized oncogenic signaling cascades rather than a direct or exclusive association with cancer-specific pathogenic drivers. Many of these pathways, such as PI3K/AKT/MTOR, MAPK/ERK, and EGFR signaling, function as multifunctional host regulatory hubs that are equally central to immune signaling, cellular survival, and HIV-1 latency maintenance. Consequently, pathway enrichment in cancer-associated maps should be interpreted as evidence of signaling convergence and network centrality, rather than disease specificity. Furthermore, the pharmacokinetic behavior, bioavailability, and cytotoxicity profiles of the identified metabolites remain uncharacterized. Despite these limitations, the multi-target and pathway-convergent insights generated in this study provide a rational foundation for future experimental and translational exploration of fungal-derived latency-modulating scaffolds.

4. Materials and Methods

4.1. Collection of Candida albicans Secondary Metabolites
A comprehensive library of secondary metabolites derived from Candida albicans was compiled through an extensive literature review conducted between 3–10 June 2025. Searches were performed across major scientific databases, including PubMed, Scopus, and Google Scholar, using relevant keywords such as “Candida albicans,” “secondary metabolites,” “natural products,” and “bioactive compounds.” Only metabolites characterized using gas chromatography-mass spectrometry (GC-MS) or nuclear magnetic resonance (NMR) metabolomics were included to ensure reliable chemical identification [21,22,23,24,25]. A total of 185 unique secondary metabolites were identified. Their chemical structures were retrieved in simplified molecular input line entry system (SMILES) format from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/, accessed on 10 June 2025).

4.2. Evaluation of ADMET Properties
To prioritize pharmacologically viable latency-modulating candidates, all 185 curated Candida albicans metabolites (MET 1–MET 185) were subjected to ADMET evaluations (Supplementary Table S1). Drug-likeness and pharmacokinetic properties were predicted using SwissADME (http://www.swissadme.ch/, accessed on 11 June 2025). Toxicity assessment was performed using the ProTox-3.0 platform (https://tox.charite.de/protox3/, accessed on 11 June 2025). The screening workflow employed a rigorous multiparametric filter to ensure the selection of compounds with favourable pharmacokinetic and safety profiles. First, compounds were retained if they met the acceptable MW range (200 ≤ MW ≤ 500 g/mol) for OB. Compliance with Lipinski’s Ro5 was strictly enforced, retaining only those compounds with zero violations. Solubility was assessed using the estimated solubility (ESOL) LogS model with a minimum threshold of –6.0, which is indicative of at least moderate aqueous solubility. Lipophilicity was evaluated using the implicit LOGP (iLOGP) model, with a cutoff value of ≤ 5 to indicate acceptable membrane permeability and bioavailability. Additionally, compounds with a TPSA ≤ 140 Å2 were considered likely to possess favourable permeability characteristics. GI absorption and BBB permeability are essential criteria for prioritizing compounds capable of targeting both peripheral and central HIV reservoirs. Furthermore, to minimize the risk of nonspecific interactions, metabolites containing PAINS were excluded. Additionally, SA scores were considered, and only compounds with SA < 6 were retained, reflecting the reasonable feasibility of the chemical synthesis [85]. Toxicity profiling included the prediction of LD50, toxicity class assignment, organ-specific toxicities (hepatotoxicity, neurotoxicity, nephrotoxicity, respiratory toxicity, cardiotoxicity, and carcinogenicity), and toxicity endpoints such as immunotoxicity, mutagenicity, and cytotoxicity. Compounds are categorized into six toxicity classes based on their predicted LD50 values: Classes I (highly toxic), Class II (moderately toxic), Class III (toxic), Class IV (harmful), Class V (may be harmful), and Class VI (non-toxic) [86].

4.3. Human Target Prediction of Candida albicans Secondary Metabolites
The six secondary metabolites that passed ADMET screening were subjected to target prediction using SwissTargetPrediction (http://www.swisstargetprediction.ch/, accessed on 13 June 2025), with the species parameter set to “Homo sapiens” and a probability threshold > 0. This platform predicts the most probable macromolecular targets of bioactive compounds based on two-dimensional (2D) and three-dimensional (3D) molecular similarity to known ligands with experimentally validated binding data [87]. Its application is particularly relevant to this study, given the limited experimental data available for fungal metabolites in the context of HIV-1 latency. TargetNet (probability ≥ 0.7) was employed to enhance the prediction coverage. TargetNet uses quantitative structure-activity relationship (QSAR) models and machine learning (ML) algorithms to predict high-confidence targets, making it suitable for underexplored natural products. The union of predicted targets from both platforms were cross-referenced with the human proteome using the UniProt Knowledgebase (UniProtKB) (https://www.uniprot.org/) to ensure proper annotation.

4.4. Collection of Candida albicans-HIV Latency-Associated Genes
Human genes annotated in HIV latency were retrieved from the GeneCards database (https://www.genecards.org/) using the keyword “HIV latency,” accessed on 14 June 2025. Only Genes with a relevance score > 0 were compiled. GeneCards is a comprehensive bioinformatics resource that integrates genomic, transcriptomic, proteomic, and disease-related data to provide high-confidence gene function associations. Its built-in relevance scoring system facilitates the prioritization of genes that are most strongly linked to HIV latency mechanisms [88]. To identify potential interactions between Candida albicans and HIV-latency-related human pathways, overlapping genes were identified using Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/, accessed on 14 June 2025). The analysis compared C. albicans-related gene targets (derived from prior target mapping) with the HIV latency-associated gene set obtained from the GeneCards database. Overlapping genes were annotated and standardized using UniProtKB (https://www.uniprot.org/) to ensure consistent nomenclature for downstream analyses.

4.5. MET-TAR Network Construction
To elucidate the multitarget interactions of Candida albicans secondary metabolites with HIV latency-associated host proteins, a MET–TAR network was constructed using Cytoscape v3.10.3 (https://cytoscape.org/, accessed on 15 June 2025). In this network, nodes represent secondary metabolites or the overlapping protein targets, and edges represent their interactions. The topological characteristics of the network were analyzed using the NetworkAnalyzer tool v4.5.0 in Cytoscape to compute metrics such as DC, BC, and CC [89].

4.6. PPI Network Analyses
To investigate the functional interactions among HIV-latency-associated host targets modulated by Candida albicans metabolites, a PPI network was constructed. Overlapping genes between the predicted targets of the secondary metabolites and HIV latency-related gene targets were retrieved from the STRING database v12.0 (https://string-db.org/, accessed on 16 June 2025), with the search restricted to “Homo sapiens.” A minimum interaction confidence score of 0.400 was applied, and associations with FDR > 5% were excluded to enhance network reliability. The resulting interaction network was visualized and analyzed in Cytoscape v3.10.3 [89]. To prioritize key regulatory genes, hub nodes were identified based on DC using the CytoHubba plugin v0.1 [90]. DC quantifies the number of direct connections each node has within the network. Nodes with higher degree values were considered more influential, as they represent highly connected proteins capable of exerting broad control over latency-associated signaling. The top 20 hub genes with the highest degree scores were selected as candidate targets for therapeutic exploration in the present study.

4.7. GO and KEGG Pathway Annotations
To comprehensively explore the biological roles and signaling pathways associated with HIV latency-related targets modulated by Candida albicans metabolites, functional enrichment analysis was conducted using ShinyGO v0.82 (https://bioinformatics.sdstate.edu/go/, accessed on 17 June 2025) [91]. Analysis was performed for overlapping genes between the predicted metabolite targets and HIV latency-associated host genes. All analyses were carried out under the species parameter “Homo sapiens.” GO analysis categorized functional annotations into BP, CC, and MF. KEGG enrichment analysis was used to identify the signaling pathways that are potentially modulated by these genes. To ensure robustness, the results were filtered using an FDR of < 0.05. Enriched terms were ranked by FDR-adjusted p-values and FE, and the top 30 terms from each GO category and KEGG pathway dataset were visualized as bubble plots, representing the gene count and statistical significance (−log10[FDR]).

4.8. Construction of MET–TAR–PATH Network
To delineate the multiscale regulatory interactions between Candida albicans metabolites, host target genes, and HIV latency-associated signaling pathways, a composite MET–TAR–PATH network was constructed using Cytoscape v3.10.3. This network was generated by integrating the MET–TAR associations with target–pathway (TAR–PATH) relationships using the six candidate secondary metabolites of C. albicans, the top 20 hub genes, and the top 30 KEGG pathways. In the network, nodes correspond to either metabolites, target genes, or pathways, and edges represent the predicted interactions or functional associations. Topological characteristics were evaluated using the NetworkAnalyzer plugin v4.5.0 in Cytoscape, with DC as the primary metric to quantify node connectivity and identify key hubs [89].

4.9. Molecular Docking Calculations
Molecular docking calculations were performed to evaluate the binding potential of the six candidate Candida albicans secondary metabolites against 11 key targets with topological importance from MET–TAR–PATH analysis. 3D structures of the target proteins were obtained from the RCSB Protein Data Bank (PDB) (https://www.rcsb.org/, accessed on 20 June 2025). The selection criteria included: (i) X-ray crystal structures resolved at ≤ 3.0 Å, (ii) origin from Homo sapiens, (iii) presence of co-crystallized ligands, and (iv) selection of the highest-resolution structure when multiple entries existed. RELA and CCND1 structures lack co-crystallized ligands; therefore, their binding sites were predicted using CASTp 3.0 (http://sts.bioe.uic.edu/castp/, accessed on 20 June 2025), and validated using published structural evidence [92]. 3D ligand structures were retrieved from the PubChem database, energy-minimized using the General AMBER Force Field (GAFF) v1.19, assigned Gasteiger partial charges, and saved in PDBQT format using UCSF Chimera v1.19. Protein structures were also prepared using UCSF Chimera v1.19, which included the removal of water molecules, separation of non-essential heteroatoms, addition of polar hydrogens, and modeling of missing side chains using Modeller v10.5 [93]. Docking simulations were performed using AutoDock Vina v1.1.2 [94]. For protocol validation, the redocking of native ligands was performed to assess docking accuracy by calculating the RMSD values between the docked and experimental poses (Supplementary Figure S3). Binding sites were defined based on the coordinates of the co-crystallized ligands or the predicted active pockets of the apo structure. Docked poses were analyzed using the protein-ligand interaction profiler (PLIP), which details non-covalent interactions [95].

4.10. MD Simulations
To assess the structural stability and dynamic behavior of the top protein–ligand complexes, MD simulations were performed using the AMBER18 package (University of California, San Francisco [UCSF], USA) with the PMEMD.CUDA module for GPU-accelerated computation [96]. This simulation approach enables the evaluation of atomic motion, intermolecular interactions, and conformational changes over time. The ligands were parameterized using ANTECHAMBER with RESP-generated charges and GAFF. The protein structures were preprocessed using pdb4amber, and each complex was solvated in a TIP3P water box with a 12 Å buffer and neutralized with Na+ or Cl− counterions [97]. Energy minimization was first conducted using restrained minimization (500 kcal/mol/Å2 restraint) for 2500 steps, followed by unrestrained minimization for 200 steps. Heating was carried out from 0 to 310 K over 50 ps under the NVT ensemble using a Langevin thermostat, with 10 kcal/mol/Å2 positional restraints on the solute. This was followed by pressure equilibration under the NPT ensemble at 1 atm using the Berendsen barostat. SHAKE constraints were applied to all hydrogen-containing bonds, and production runs were performed for 200 ns with a 2 fs time step under periodic boundary conditions. Post-simulation analysis was performed using the CPPTRAJ module to compute the RMSD and RMSF. Trajectories were visualized and inspected using UCSF Chimera v1.19 and PyMOL v3.1.4, and plots were generated using Origin 2018 software [20,98].

4.11. BFE Calculations
BFE calculations were performed using the molecular mechanics/Generalized Born surface area (MM/GBSA) method [99]. This approach was implemented using the MMPBSA.py module in the AMBER18 software suite, which couples molecular mechanics with an implicit solvent model to provide dynamic trajectory-based estimations of binding affinities that are superior to static docking scores [100,101]. For each protein–ligand complex, 50,000 frames were uniformly sampled from 200 ns MD trajectories to ensure robust averaging. The BFE (ΔGbind) was computed using the following thermodynamic cycle:
The total energy is further decomposed into the following components:
The gas-phase energy (Egas) was contributed by the internal (Eint), van der Waals (Evdw), and electrostatic (Eele) terms. The solvation free energy (Gsol) was contributed by the polar term (GGB) from the Generalized Born (GB) model and the nonpolar term (GSA) from the solvent-accessible surface area (SA) model, with the surface tension constant (γ) set at 0.0072 kcal/mol/Å2. The entropic contribution (T∆S) was calculated using normal mode analysis.

5. Conclusions
This study provides the first integrative systems-level insight into how Candida albicans secondary metabolites may influence HIV-1 latency through multitarget engagement of interconnected signaling networks. Through network pharmacology, molecular docking, and dynamics interaction analyses, the metabolites MET 15, MET 28, and MET 119 were identified as promising polypharmacological modulators targeting key latency-associated hubs, such as PIK3CA, MAPK1, MAPK3, AKT1, EGFR, and MTOR. The enrichment of KEGG pathways, including PI3K/AKT, MAPK, PD-1/PD-L1 checkpoint, and HIF-1 signaling, underscores the mechanistic plausibility of these metabolites in modulating cellular survival and transcriptional programs that sustain proviral persistence. By uncovering the potential intersections between fungal metabolism and HIV latency regulation, we identified MET 15 as a top-ranked metabolite and an unexplored reservoir of latency-modulating scaffolds with potential translation in HIV cure research.
The dual nature of C. albicans as both an opportunistic pathogen and a reservoir of bioactive secondary metabolites presents a unique paradox for antiviral drug development. Our findings suggest that metabolites such as MET 181 (farnesol) are not merely passive byproducts of fungal metabolism but may actively shape the host’s signaling landscape. By modulating key pathways, including the PI3K/AKT/MTOR and MAPK/ERK axes, these compounds could potentially be leveraged as either “block and lock” agents to reinforce HIV latency or as LRAs to support “shock and kill” strategies, depending on their specific agonistic or antagonistic activities.
Future studies should experimentally validate the identified metabolites using latency models such as J-Lat or primary CD4+ T-cell systems to confirm their effects on proviral transcription. Mechanistic assays integrating transcriptomic and proteomic profiling are recommended to delineate the signaling nodes through which these metabolites influence the PI3K/AKT, MAPK, and NF-κB pathways. In addition, structure–activity relationship (SAR) studies and rational analog design can optimize their potency, selectivity, and pharmacokinetic properties. Finally, combining metabolomics with host–pathogen interaction datasets may clarify whether these fungal metabolites act synergistically with known LRAs or immune modulators.