Prediction and screening of potential targets of PLB
The chemical structure information of PLB was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Molecular targets of PLB were predicted based on ligand-based two-dimensional and three-dimensional similarity using the SwissTarget (https://swisstargetprediction.ch/index.php), SEA (https://sea.bkslab.org/), TargetNet (http://targetnet.scbdd.com/calcnet/index/), PharmMapper (https://www.lilab-ecust.cn/pharmmapper/index.html), Comparative Toxicogenomics Database (CTD, https://ctdbase.org/), Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, https://www.tcmsp-e.com/index.php), and HERB (http://herb.ac.cn/) databases. To ensure data reliability and reproducibility, database‑specific filtering thresholds were applied as follows: TCMSP retained targets with oral bioavailability (OB) ≥ 30% and drug‑likeness (DL) ≥ 0.18; SwissTargetPrediction kept entries with a prediction probability ≥ 0.5; PharmMapper selected targets with a normalized fit score ≥ 4.0; SEA included only targets with an E‑value < 0.001 as significantly enriched; TargetNet retained targets with a prediction probability > 0.5; CTD included only targets with curated evidence (inferred from chemical–gene interactions) and an interaction score > 0.3; and HERB retained targets with a literature‑supported score ≥ 0.4. All retrieved protein identifiers from the above databases were standardized to official human HGNC gene symbols using the UniProt database (https://www.uniprot.org/), with species restricted to Homo sapiens. After removing duplicate entries across all sources, a total of 500 unique PLB‑related targets were obtained for subsequent analysis.
Retrieval of prostate cancer-associated genes
Disease-related genes for prostate cancer were obtained from GeneCards (https://www.GeneCardss.org/), DrugBank (https://go.drugbank.com/), CTD, and HERB. Database-specific inclusion criteria were applied to ensure data reliability. GeneCardss retained genes with a relevance score ≥ 0.5, as this threshold captures genes with moderate-to-strong evidence linking them to the disease query. DrugBank included only entries with experimental evidence (e.g., FDA-approved or investigational drugs for prostate cancer) and excluded computationally predicted or theoretical interactions. CTD retained only records with evidence levels classified as "marker" or "mechanism" based on curated chemical-gene-disease interactions. HERB included targets with a literature-supported confidence score ≥ 0.4 to ensure sufficient experimental or text-mining evidence. All gene symbols were unified to the HGNC nomenclature via the UniProt database (https://www.uniprot.org/), with species restricted to Homo sapiens, and duplicates were removed, yielding 1,199 unique prostate cancer–related targets for subsequent analysis.
Overlapping target identification and conflict annotation correction
The 500 PLB targets and 1,199 prostate cancer targets were intersected using standardized HGNC symbols in Venny 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/index.html), yielding 151 overlapping candidate targets. To resolve cross‑database annotation conflicts, a three‑step correction pipeline was implemented: (i) gene alias inconsistencies were unified via the UniProt ID mapping tool; (ii) paralogous redundant genes were excluded using CD‑HIT with a sequence similarity threshold > 0.4; and (iii) contradictory functional annotations were retained only if supported by at least two independent databases, while unique conflicting descriptions from a single source were discarded.
Network construction and hub identification
A protein-protein interaction (PPI) network was constructed by submitting the overlapping targets to the STRING database (version 12.0, https://cn.string-db.org/), with the search restricted to Homo sapiens. Only interactions with a combined confidence score ≥ 0.700 were retained to ensure high reliability, and all nodes lacking any connections were excluded from the network. The constructed PPI network was visualized and analyzed topologically in Cytoscape (version 3.10.2), after which hub genes were determined using the cytoHubba plugin, employing MCC (Maximal Clique Centrality) as the primary ranking algorithm and Degree as an auxiliary cross‑validation metric.
Functional enrichment analysis
Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for the differentially expressed genes were performed using the Bioinformatics Online Analysis Platform (https://www.bioinformatics.com.cn, last accessed on May 4, 2026). The platform employs Fisher's exact test to calculate raw P‑values, and the Benjamini-Hochberg (BH) method was used for multiple-testing correction to control the false discovery rate (FDR). GO terms (biological process, cellular component, and molecular function) and KEGG pathways with an FDR < 0.05 were considered significantly enriched.
Molecular docking
Crystal structures of the selected target proteins, representing active conformations, were obtained from the RCSB Protein Data Bank (https://www.rcsb.org/). Protein preparation was performed in Discovery Studio, including removal of water molecules and heteroatoms, deletion of redundant chains, repair of missing residues, addition of hydrogen atoms at pH 7.4, assignment of Gasteiger charges, and energy minimization. The binding pocket was defined based on the co‑crystallized ligand coordinates, using a grid that encompassed all critical substrate‑binding residues. Molecular docking of PLB was carried out using the SwissDock web server (http://www.swissdock.ch/) via the Attracting Cavities module. Binding free energies were calculated to assess affinity, and the lowest‑energy conformation for each target was selected as the final docking pose. The results of the molecular docking simulations were visualized to validate them.
Molecular dynamics simulation
Molecular dynamics (MD) simulations were carried out using GROMACS 2022.2. The Amber14SB force field was used to describe the protein, with the system solvated in TIP3P water. Plumbagin parameters, including AM1-BCC partial charges and GAFF2 atom types, were generated using Antechamber, followed by topology conversion with ACPYPE and assignment of Joung-Cheatham ion parameters. Each protein-ligand complex was embedded in a truncated dodecahedral simulation box with a minimum protein-to-box distance of 1.2 nm, solvated with TIP3P water molecules, and neutralized by adding 0.15 M Na⁺/Cl⁻. After energy minimization using the steepest descent algorithm (Fmax < 1,000 kJ·mol⁻1·nm⁻1), the system underwent sequential NVT and NPT equilibration for 200 ps each at 298 K. Production simulations were then performed for 200 ns under NPT conditions with a 2 fs integration time step using the Verlet cutoff scheme. Long-range electrostatic interactions were calculated with the particle mesh Ewald (PME) method, while both electrostatic and van der Waals interactions employed a cutoff distance of 1.2 nm. The LINCS algorithm was used to constrain hydrogen-containing bonds. Temperature was maintained at 298 K using the Nosé-Hoover thermostat and pressure at 1 bar using the Parrinello-Rahman barostat. Coordinates were saved every 10 ps for subsequent analyses. Trajectory analysis and visualization were performed using GROMACS utilities, VMD, and PyMOL, whereas MM-PBSA binding free energy calculations were carried out with gmx_MMPBSA where appropriate.
Equilibration assessment of molecular dynamics trajectories and MM-PBSA free energy sampling protocol
For each protein–plumbagin complex, a 200 ns all-atom production molecular dynamics simulation was performed using GROMACS. The 100 ns time point was set as the equilibration boundary: the 0–100 ns segment was designated as the conformational relaxation phase, during which the protein backbone and ligand-binding pocket underwent continuous conformational adjustments; the 100–200 ns segment was identified as the thermodynamically stable plateau, as evidenced by the absence of unidirectional drift in RMSD, radius of gyration (Rg), residue-wise RMSF, buried SASA of the ligand, intermolecular hydrogen bonds, and ligand–receptor interaction energies, which exhibited only minor steady-state fluctuations. All quantitative kinetic parameters and MM/PBSA binding free energies were calculated exclusively from the 100–200 ns equilibrium phase, from which frames were uniformly extracted every 100 ps, yielding 1,000 equilibrium snapshots per system for gmx_MMPBSA input. The first 100 ns of relaxation trajectories were discarded to eliminate interference from conformational drift in free-energy calculations.
Prognostic analysis of target genes in prostate cancer
Prognostic analyses were conducted using prostate cancer datasets from The Cancer Genome Atlas (TCGA). RNA sequencing (STAR-counts) data and corresponding clinical information were obtained from the Genomic Data Commons portal (https://portal.gdc.cancer.gov). Gene expression counts were converted to transcripts per million (TPM) and normalized using log2(TPM + 1). Following the exclusion of samples with incomplete clinical information, 498 cases were included in the analysis. The median expression level of each gene was used to divide patients into high- and low-expression groups. Kaplan-Meier survival analysis with the log-rank test and univariate Cox proportional hazards regression were performed to evaluate overall survival (OS) for SRC and progression-free survival (PFS) for MAPK3, with hazard ratios (HRs) and 95% confidence intervals (CIs) reported. Predictive performance was further assessed by generating time-dependent receiver operating characteristic (ROC) curves at 1, 3, and 5 years. Statistical analyses were carried out in R version 4.0.3, with P < 0.05 considered statistically significant.