Integrative network toxicology, transcriptomic, and molecular docking approaches to elucidate the toxicity and mechanisms of bisphenol A in stroke.

In the present study, stroke-related transcriptomic data were retrieved from the public GEO database ( https://www.ncbi.nlm.nih.gov/geo ), with a focus on the GSE16561 dataset associated with stroke. This dataset contains gene expression profiles derived from peripheral whole blood RNA samples, including 39 cases from patients with ischemic stroke and 24 samples from healthy individuals serving as controls. Additionally, BPA target genes were collected using three databases: ChEMBL ( http://www.ebi.ac.uk/chembl/ ), DrugBank ( https://go.drugbank.com/ ), and SwissTargetPrediction ( http://www.swisstargetprediction.ch/ ). Before conducting the analysis, all datasets were subjected to strict preprocessing steps, such as eliminating duplicate entries and performing quality control checks, to guarantee the data’s integrity and dependability for subsequent bioinformatics analyses. The ProTox-3.0 platform ( https://tox.charite.de ) was employed to predict the toxicity classification, toxicological endpoints, and putative molecular targets of BPA, thereby evaluating its potential relevance to stroke. A toxicity radar plot and target interaction network for BPA were generated, and a comprehensive toxicological report was downloaded for further interpretation. After applying Z-score normalization to the transcriptomic data from the GSE16561 dataset, we utilized the WGCNA R package to perform WGCNA and build a weighted gene co-expression network. A correlation matrix was computed, with a soft-thresholding power (β) chosen according to the scale-free topology standard. The optimal β value was established by assessing the association between scale-free fit and mean connectivity. Sample clustering dendrograms were constructed to visualize gene expression patterns before and after module detection. Gene modules were identified through hierarchical clustering based on the Topological Overlap Matrix (TOM) and delineated using the dynamic tree cut algorithm. Each module was given a distinct color for easy identification. Subsequently, module eigengenes were correlated with stroke phenotype, and a trait heatmap was created to pinpoint the module with the strongest association with stroke, which was then selected for further analyses. For PCA analysis, the normalized and filtered gene expression matrix was used as input. PCA was performed using the standardize gene expression values (mean = 0, standard deviation = 1) to avoid bias from genes with high absolute expression levels. The first two principal components (PC1 and PC2) were extracted, and their respective variance explanation rates were calculated. A scatter plot was generated using the ggplot2 package, where each point represents a sample, and samples were colored by group (Disease vs. Control) to visualize the separation pattern between groups. Additionally, a scree plot was constructed to display the variance explanation rates of the top 10 principal components, evaluating the overall distribution of transcriptional variation across samples. After preprocessing the data with R software (version 4.2.1), the resulting normalized gene expression matrix was fed into limma. A linear model was built via the lmFit function, and empirical Bayes moderation was implemented through the eBayes function to calculate moderated t-statistics and p-values for each gene. For managing multiple testing, p-values were adjusted via the Benjamini-Hochberg method to regulate the false discovery rate (FDR). DEGs were visualized using the ggplot2 package. A volcano plot was produced with -log10(P-value) (P-value  0.263) were marked red, and downregulated genes (logFC < -0.263) were shaded blue. Additionally, violin plots were created to illustrate the distribution of expression levels for top DEGs across sample groups. To evaluate sample clustering patterns based on DEGs, hierarchical clustering was conducted with Euclidean distance and complete linkage. Following Z-score normalization of DEG expression values, a heatmap was generated using the pheatmap package, with rows denoting genes and columns corresponding to samples, and color intensity indicating expression levels to highlight distinct patterns between stroke and control groups. The DEGs, WGCNA disease-related module genes, and predicted BPA target genes were collated, and their intersections, representing the potential targets of BPA in stroke, were displayed using a Venn diagram. The screened cross-targets were imported into the STRING database (version 12.0; https://string-db.org/ ) with Homo sapiens selected as the species [ 30 ]. To ensure biologically meaningful interactions, the minimum confidence score was set to greater than 0.700, and discrete nodes were hidden to remove low-confidence or unconnected proteins. The resulting PPI data, including nodes and edges, were exported and visualized in Cytoscape (version 3.10.3) for topological analysis [ 31 ]. Hub genes were identified using the CytoHubba plugin. Three classical algorithms were applied, including maximal clique centrality (MCC), closeness, and degree, to independently rank the connectivity of all nodes. The top ten genes from each method were extracted. To obtain reliable core targets, the three gene lists were intersected, and the intersection was visualized on the Microbiome online tools ( http://www.bioinformatics.com.cn/ ). These intersecting genes were defined as the core targets involved in BPA-mediated neurotoxicity and were used for subsequent enrichment and molecular validation. To explore the functions of the core target genes and the associated signaling pathways, GO and KEGG enrichment analyses were performed using the DAVID database, with human species as the focus. GO analysis was conducted across biological processes (BP), molecular functions (MF) and cellular components (CC) to clarify the roles of core target genes in pathways related to energy metabolism and cellular respiration. KEGG analysis was used to identify signaling pathways enriched for these targets, with emphasis on those associated with metabolism, oxidative stress and neurological function. KEGG analysis was used to identify signaling pathways enriched for these targets, with emphasis on those associated with metabolism, oxidative stress and neurological function. The parameters for enrichment analysis: gene set size between 5 and 5000, a P value threshold of < 0.05, false discovery rate (FDR) < 0.2, and multiple hypothesis testing correction using the Benjamini & Hochberg method. The 3D crystal structures of all target proteins were retrieved from the RCSB PDB database ( https://www.rcsb.org/ ) using their corresponding gene symbols. During protein structure preparation, all water molecules and original ligands were removed to prevent potential interference with the initial docking conformation and to focus on direct interactions between BPA and the protein receptor. Polar hydrogen atoms were added, and Kollman charges were assigned using AutoDockTools. The 3D structure of BPA was obtained from the PubChem database, and Gasteiger charges, non-polar hydrogens, and rotatable bonds were assigned using AutoDockTools. A global docking approach was employed without predefined active sites. The docking grid box was set to encompass the entire protein receptor to explore all potential binding modes. Docking results were ranked by binding free energy, and representative binding conformations were selected for further analysis. For these conformations, the ligand-binding pocket, defined as residues within 5 Å of the ligand, was characterized, including calculation of the Kyte-Doolittle hydrophobicity index to assess the physicochemical properties of the binding microenvironment. The stability of binding between the small molecule and the protein receptor is inversely related to binding energy—lower binding energy indicates stronger affinity. A binding free energy below − 5 kcal/mol was taken as an indicator of good binding activity [ 32 ]. Ultimately, the conformation with the lowest binding energy was chosen, and PyMOL software was used to generate a 3D binding mode diagram, which further illustrates how BPA binds to the target protein.

Toxicity predictions from the ProTox3.0 database indicated that BPA has an acute oral lethal dose (LD₅₀) of 4950 mg/kg, classifying it as a low-toxicity compound (Class 5) (Fig.  1 A). However, BPA exhibited higher-than-average activity scores for estrogen receptor α, estrogen receptor ligand-binding domains, and mitochondrial membrane potential, suggesting potential endocrine and mitochondrial toxicity (Fig.  1 B). Network analysis of predicted toxicity–target interactions further revealed that BPA may affect several critical biological systems, including the blood–brain barrier (BBB), cytochrome P450 enzymes (e.g., CYP2C19), and estrogen receptor–related signaling pathways (Fig.  1 C). Consistent with these findings, the toxicity model predicted BPA to be active in estrogen receptor signaling (probability 1.0), mitochondrial membrane potential disruption (probability 1.0), cytochrome P450–mediated metabolism (probability 0.68–0.77), ecotoxicity (probability 0.63), and BBB penetration (probability 0.53) (Fig.  1 D). Fig. 1 Toxicological profile of BPA based on ProTox3.0 analysis. ( A ) The BPA toxicity panel presents its molecular structure, chemical properties and the key indicators for toxicity prediction. ( B ) Radar map for toxicity prediction. Blue represents the predicted data for BPA, and orange represents the average probability distribution of known reference molecules. ( C ) The toxicity target interaction network illustrates predicted targets and pathways, where red nodes indicate probabilities ≥ 0.7 and pink nodes indicate probabilities<0.7. ( D ) The toxicity prediction report summarizes the classification results, target names, abbreviations, prediction outcomes and associated probabilities Toxicological profile of BPA based on ProTox3.0 analysis. ( A ) The BPA toxicity panel presents its molecular structure, chemical properties and the key indicators for toxicity prediction. ( B ) Radar map for toxicity prediction. Blue represents the predicted data for BPA, and orange represents the average probability distribution of known reference molecules. ( C ) The toxicity target interaction network illustrates predicted targets and pathways, where red nodes indicate probabilities ≥ 0.7 and pink nodes indicate probabilities<0.7. ( D ) The toxicity prediction report summarizes the classification results, target names, abbreviations, prediction outcomes and associated probabilities Based on the stroke-related transcriptome dataset GSE16561 from the GEO database, we constructed a gene co-expression network using WGCNA. Sample clustering revealed that all samples fell within a reasonable range, with no outliers detected, indicating high data quality (Fig.  2 A). The soft-thresholding power (β) was determined according to the scale-free topology criterion, with β = 19 selected when the scale-free fitting index (R²) exceeded 0.80 while maintaining moderate mean connectivity (Fig.  2 B). Dynamic tree cutting was then applied to identify distinct co-expression modules, followed by correlation analysis with clinical traits. The results demonstrated that several modules were significantly associated with stroke-related features. Specifically, the ME2 module showed a strong positive correlation with age, whereas the ME3 module was negatively correlated with stroke type (Fig.  2 C–D). The ME3 module, which contains 1161 module genes, was selected for subsequent analysis. Fig. 2 Identification of key gene modules by WGCNA. ( A ) The sample clustering dendrogram illustrates the relationships among all samples and assists in identifying potential outliers. ( B ) The soft-threshold screening plot shows the scale-free topology fit index across threshold values on the left and the corresponding changes in average connectivity on the right. ( C ) The gene module clustering dendrogram displays the hierarchical relationships among modules. ( D ) The module–trait heatmap presents the associations between gene modules and clinical traits. Red indicates positive correlation, blue indicates negative correlation, numerical values represent correlation coefficients and the corresponding P values are shown in parentheses Identification of key gene modules by WGCNA. ( A ) The sample clustering dendrogram illustrates the relationships among all samples and assists in identifying potential outliers. ( B ) The soft-threshold screening plot shows the scale-free topology fit index across threshold values on the left and the corresponding changes in average connectivity on the right. ( C ) The gene module clustering dendrogram displays the hierarchical relationships among modules. ( D ) The module–trait heatmap presents the associations between gene modules and clinical traits. Red indicates positive correlation, blue indicates negative correlation, numerical values represent correlation coefficients and the corresponding P values are shown in parentheses PCA was performed to assess the overall expression patterns between normal and stroke samples. In the two-dimensional space defined by the first two principal components, PC1 and PC2 explained 51% and 15% of the total variance, respectively. Although partial overlap was observed along PC1, the two groups displayed a distinct separation trend along the PC2 axis (Fig.  3 A). Contribution analysis of the top ten principal components confirmed that PC1 accounted for the majority of the total variance (> 50%), indicating it was the main source of variation (Fig.  3 B). Differential expression analysis identified 1,801 significantly DEGs, including 689 upregulated and 1,112 downregulated genes (|logFC| > 0.263, FDR < 0.01) (Fig.  3 C). The consistent distribution of normalized expression values across samples demonstrated reliable data preprocessing and good comparability between groups (Fig.  3 D). Bidirectional hierarchical clustering based on all DEGs revealed a clear transcriptomic distinction between stroke and control samples, highlighting substantial gene expression heterogeneity associated with stroke (Fig.  3 E). Similarly, a heatmap of the top 20 most significantly dysregulated genes further confirmed distinct clustering patterns between the two groups (Fig.  3 F). Fig. 3 Results of PCA and differential expression analysis. ( A ) PCA of transcriptomes from the disease and control groups, with the control group shown in blue and the disease group in orange. ( B ) Bar chart illustrating the contributions of the top ten principal components. ( C ) Volcano plot of differential expression, where blue indicates down-regulated genes, gray indicates non-significant genes and red indicates up-regulated genes; key differentially expressed genes are labeled. ( D ) Gene expression profiles of all samples, with blue representing the control group and orange representing the disease group. ( E ) Hierarchical clustering heatmap of DEGs, in which each cell reflects the expression level of one gene in one sample; red denotes high expression and blue denotes low expression. ( F ) Heatmap of the top 20 significantly DEGs Results of PCA and differential expression analysis. ( A ) PCA of transcriptomes from the disease and control groups, with the control group shown in blue and the disease group in orange. ( B ) Bar chart illustrating the contributions of the top ten principal components. ( C ) Volcano plot of differential expression, where blue indicates down-regulated genes, gray indicates non-significant genes and red indicates up-regulated genes; key differentially expressed genes are labeled. ( D ) Gene expression profiles of all samples, with blue representing the control group and orange representing the disease group. ( E ) Hierarchical clustering heatmap of DEGs, in which each cell reflects the expression level of one gene in one sample; red denotes high expression and blue denotes low expression. ( F ) Heatmap of the top 20 significantly DEGs A total of 689 BPA-associated target genes were retrieved from the ChEMBL and SwissTargetPrediction databases. Intersection analysis between these targets, stroke-related DEGs, and WGCNA-identified module genes yielded 22 overlapping candidate genes (Fig.  4 A). These genes were then imported into the STRING database for PPI analysis, applying a minimum interaction confidence threshold of 0.700. Twelve genes with strong interaction relationships were retained and used to construct a PPI network, which was visualized in Cytoscape (Fig.  4 B). To identify key nodes within the network, three topological parameters (MCC, closeness, and degree) were calculated to rank gene importance (Fig.  4 C–E). Integrating the results from the three ranking methods, ten hub genes consistently ranked among the top across all metrics: NDUFB7, NDUFA11, NDUFB11, NDUFB8, NDUFS8, NDUFA12, NDUFS5, NDUFA8, NDUFV1, and NDUFA13 (Fig.  4 F). Fig. 4 Intersection analysis of DEGs, WGCNA module genes, and BPA targets for identification of key genes. ( A ) Venn diagram showing overlapping genes. ( B ) PPI network of the potential targets, where nodes represent proteins and edges represent their interactions. ( C – E ) Ranking of network nodes based on ( C ) MCC, ( D ) closeness, and ( E ) degree centrality. ( F ) Venn diagram of top10-ranked hub genes across three metrics Intersection analysis of DEGs, WGCNA module genes, and BPA targets for identification of key genes. ( A ) Venn diagram showing overlapping genes. ( B ) PPI network of the potential targets, where nodes represent proteins and edges represent their interactions. ( C – E ) Ranking of network nodes based on ( C ) MCC, ( D ) closeness, and ( E ) degree centrality. ( F ) Venn diagram of top10-ranked hub genes across three metrics To elucidate the potential biological functions of the ten core targets and their associated signaling pathways, GO and KEGG enrichment analyses were performed using the DAVID database, with Homo sapiens selected as the reference species. GO analysis identified twelve significantly enriched terms, including six BP, four CC, and two MF. In the BP category, the core genes were primarily involved in mitochondrial energy metabolism, encompassing ATP synthesis coupled to proton transport, aerobic respiration, and electron transport processes (Fig.  5 A). CC enrichment indicated a predominant localization within energy-related mitochondrial structures, such as respiratory chain complex I and the mitochondrial inner membrane (Fig.  5 B). MF analysis further highlighted NADH dehydrogenase activity and 4Fe–4 S cluster binding as key functional features (Fig.  5 C). KEGG pathway enrichment revealed fourteen significantly associated pathways, predominantly linked to core energy metabolism processes such as oxidative phosphorylation. Additional pathways were related to intracellular signaling, metabolic disorders, oxidative stress, and several neurodegenerative diseases (Fig.  5 D). Collectively, these findings suggest that BPA may contribute to metabolic dysregulation and neuronal injury by disturbing mitochondrial energy metabolism, enhancing oxidative stress, and impairing signal transduction. We therefore hypothesize that BPA exposure disrupts oxygen utilization and cellular energy supply in brain tissue, promotes reactive oxygen species generation, and interferes with DNA replication, ultimately leading to neural dysfunction and paralysis. Fig. 5 Functional enrichment analysis of potential core genes. ( A – C ) GO enrichment results for ( A ) BP, ( B ) CC, and ( C ) MF. ( D ) KEGG pathway enrichment analysis identified fourteen enriched pathways. Circle size reflects the number of genes within each pathway, and circle color indicates enrichment significance. Darker colors represent higher –log 10 (p-value) values and therefore stronger enrichment Functional enrichment analysis of potential core genes. ( A – C ) GO enrichment results for ( A ) BP, ( B ) CC, and ( C ) MF. ( D ) KEGG pathway enrichment analysis identified fourteen enriched pathways. Circle size reflects the number of genes within each pathway, and circle color indicates enrichment significance. Darker colors represent higher –log 10 (p-value) values and therefore stronger enrichment Molecular docking analysis was conducted to further validate the interactions between BPA and ten potential stroke-related target proteins. The results demonstrated that BPA could spontaneously bind to all target proteins, with binding energies below − 5.0 kcal/mol, suggesting stable interactions (Fig.  6 ; Table  1 ). Among these, BPA exhibited the strongest binding affinity with NDUFA11 (–7.7 kcal/mol), with Tyr9 identified as the key interacting residue; this interaction was stabilized by a combination of hydrogen bonding and pronounced van der Waals contacts between the BPA aromatic rings and surrounding hydrophobic residues. NDUFV1 exhibited the second-strongest interaction (–7.4 kcal/mol), primarily involving Ser323, accompanied by substantial hydrophobic contributions to binding stabilization. BPA also showed notable binding to NDUFS8 and NDUFA12, with binding energies of − 6.5 kcal/mol and key residues Asn202 and Trp110, respectively. The hydrogen bonds and van der Waals contacts jointly contributed to complex stability. Similarly, BPA bound to NDUFA13 and NDUFB7 with binding energies of − 6.4 kcal/mol, primarily through Val117, Glu105, and Glu101, forming a stable hydrophobic van der Waals interaction network around the aromatic rings. The remaining targets (NDUFA8, NDUFS5, NDUFB8, and NDUFB11) displayed binding energies ranging from − 5.3 to − 6.3 kcal/mol, with major interacting residues including Arg47, Ala52, Trp40, Asp87, Arg86, and Gln83, and these interactions were also maintained through hydrogen bonding and extensive van der Waals contacts. Additional hydrogen-bonding interactions involving residues such as Arg113 and Tyr112 further supported the stability and specificity of BPA–protein binding. Overall, BPA formed at least one hydrogen bond with each target protein, complemented by systematic hydrophobic and van der Waals interactions, which collectively constitute the dominant forces stabilizing the predicted complexes. Fig. 6 Molecular docking of BPA with ten core target proteins. ( A ) BPA–NDUFB7; ( B ) BPA–NDUFA11; ( C ) BPA–NDUFB11; ( D ) BPA–NDUFB8; ( E ) BPA–NDUFS8; ( F ) BPA–NDUFA12; ( G ) BPA–NDUFS5; ( H ) BPA–NDUFA8; ( I ) BPA–NDUFV1; ( J ) BPA–NDUFA13. In the visualizations, the protein backbone is shown in pink, BPA is highlighted in red, and key amino acid residues are depicted in turquoise. Stable hydrogen bonds between BPA and the protein are indicated by yellow dashed lines, with bond lengths labeled Molecular docking of BPA with ten core target proteins. ( A ) BPA–NDUFB7; ( B ) BPA–NDUFA11; ( C ) BPA–NDUFB11; ( D ) BPA–NDUFB8; ( E ) BPA–NDUFS8; ( F ) BPA–NDUFA12; ( G ) BPA–NDUFS5; ( H ) BPA–NDUFA8; ( I ) BPA–NDUFV1; ( J ) BPA–NDUFA13. In the visualizations, the protein backbone is shown in pink, BPA is highlighted in red, and key amino acid residues are depicted in turquoise. Stable hydrogen bonds between BPA and the protein are indicated by yellow dashed lines, with bond lengths labeled Table 1 Molecular docking analysis of BPA binding affinities with core targets Gene_name Uniprot_id PDB ID Protein_name Binding energy (kcal/mol) NDUFB7 P17568 9CWT: Entity 42 NADH dehydrogenase [ubiquinone] 1 beta subcomplex subunit 7 -6.4 NDUFA11 Q86Y39 9CWT: Entity 20 NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 11 -7.7 NDUFB11 Q9NX14 9CWT: Entity 28 NADH dehydrogenase [ubiquinone] 1 beta subcomplex subunit 11, mitochondrial -5.3 NDUFB8 O95169 9CWT: Entity 26 NADH dehydrogenase [ubiquinone] 1 beta subcomplex subunit 8, mitochondrial -5.8 NDUFS8 O00217 9CWT: Entity 2 NADH dehydrogenase [ubiquinone] iron-sulfur protein 8, mitochondrial -6.5 NDUFA12 Q9UI09 9CWT: Entity 13 NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 12 -6.5 NDUFS5 O43920 9CWT: Entity 31 NADH dehydrogenase [ubiquinone] iron-sulfur protein 5 -6.2 NDUFA8 P51970 9CWT: Entity 41 NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 8 -6.3 NDUFV1 P49821 9CWT: Entity 1 NADH dehydrogenase [ubiquinone] flavoprotein 1, mitochondrial -7.4 NDUFA13 Q9P0J0 9CWT: Entity 21 NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 13 -6.4 Molecular docking analysis of BPA binding affinities with core targets

Our comprehensive multi-omics analysis systematically uncovered the molecular connections between BPA exposure and stroke pathogenesis, identifying ten core target genes. Molecular docking further confirmed strong binding affinities between BPA and these targets, with the most stable interactions observed for NDUFA11 (–7.7 kcal/mol) and NDUFV1 (–7.4 kcal/mol). These targets were significantly enriched in key stroke-related pathways, including oxidative phosphorylation, neurodegenerative disease, and metabolic signaling pathways. The convergence of pathway enrichment and strong binding affinity provides a mechanistic framework suggesting that environmental BPA exposure may contribute to cerebrovascular injury by disrupting mitochondrial function and energy metabolism. Although additional experimental validation is required, our findings (1) reveal novel molecular links between BPA exposure and stroke, (2) identify priority molecular targets for future mechanistic and therapeutic studies, and (3) emphasize the need to re-evaluate the safety of BPA in food-contact materials given its potential to impair brain energy metabolism and mitochondrial integrity. Overall, this study provides new insights into the environmental pathogenesis of stroke and demonstrates the value of computational multi-omics approaches in advancing neurotoxicology research.

Although BPA has long been recognized for its neurotoxic potential as an environmental endocrine disruptor, its specific molecular mechanisms in acute neurological disorders such as stroke, which has high mortality and disability rates, remain poorly understood. In this study, we integrated network toxicology, multi-omics bioinformatics, and molecular docking analyses to systematically investigate BPA’s effects in stroke for the first time. Our findings suggest that BPA may contribute to stroke pathology by directly targeting the structure and function of mitochondrial respiratory chain Complex I. Through this approach, we identified ten core molecular targets and mapped a mechanistic pathway linking BPA exposure to cerebrovascular injury, providing a novel theoretical framework for understanding how environmental factors influence stroke development. Our multi-omics analysis and cross-validation revealed that the toxic effects of BPA converge on Complex I, a central hub of mitochondrial energy metabolism. The ten identified core targets showed distinct functional roles distributed across multiple structural modules of Complex I: β subcomplexes responsible for electron transfer (NDUFB7, NDUFB8, NDUFB11), α subcomplexes maintaining structural stability (NDUFA11, NDUFA13), iron–sulfur proteins mediating electron transport (NDUFS5, NDUFS8), and FAD-bound flavoproteins involved in oxidation–reduction reactions (NDUFV1) [ 33 – 35 ]. This broad distribution indicates that BPA’s disruption of Complex I likely impairs the enzyme’s overall functional integrity rather than affecting a single component. Given the energy-demanding nature of the brain, Complex I dysfunction can trigger a cascade of damaging effects during stroke. Reduced ATP synthesis efficiency compromises neuronal membrane potential and ion homeostasis, while BPA’s inhibition of respiratory chain enzymes further worsens mitochondrial dysfunction, leading to neuronal swelling and necrosis [ 36 – 38 ]. Simultaneously, impaired electron transport promotes excessive production of ROS, causing oxidative stress and activating mitochondria-mediated apoptosis [ 39 , 40 ]. Collectively, our results demonstrate that BPA induces mitochondrial membrane potential collapse, decreases ATP generation, and elevates ROS levels. These findings extend existing toxicological evidence to the context of stroke, suggesting that BPA may aggravate cerebral energy failure and oxidative injury by concurrently disrupting multiple subunits of Complex I. The results of the functional enrichment analyses further reinforced and expanded our understanding of the mechanisms described above. GO analysis revealed that the core targets were significantly enriched in energy metabolism–related biological processes, including mitochondrial ATP synthesis coupled with electron transport and proton transmembrane transport. These findings are consistent with previous reports showing that BPA disrupts mitochondrial structure and function, and even exposure at reference doses can alter the expression of genes regulating energy metabolism [ 41 ]. KEGG pathway analysis indicated that the identified targets were involved not only in the oxidative phosphorylation pathway but also in pathways commonly associated with neurodegenerative diseases such as Alzheimer’s disease, Parkinson’s disease, and amyotrophic lateral sclerosis. Epidemiological studies have likewise reported that exposure to bisphenolic compounds is linked to an increased risk of neurodegenerative disorders [ 42 ]. The proposed mechanisms underlying these effects include neurotransmitter disruption, excitotoxicity, neuroinflammation, BBB impairment, neuronal apoptosis, intracellular calcium imbalance, excessive ROS production, lactate dehydrogenase release, axonal degeneration, microglial DNA damage, astrocyte proliferation, and reduced myelination [ 43 , 44 ]. Together, these findings suggest that stroke and neurodegenerative diseases may share a common molecular foundation centered on mitochondrial dysfunction. In addition, the core targets were enriched in pathways associated with metabolic disorders, including nonalcoholic fatty liver disease and diabetic cardiomyopathy. This implies that BPA may indirectly compromise cerebrovascular health by promoting systemic metabolic disturbances. Consistent with this notion, previous studies have shown that bisphenolic compounds can modulate inflammatory pathways and lipid metabolism, contributing to metabolic syndrome [ 45 ]. Moreover, BPA exposure has been significantly correlated with a higher risk of abdominal obesity [ 46 ]. Notably, BPA is metabolized by uridine diphosphate–glucuronosyltransferase (UGT) enzymes, and genetic variation in UGT1A9 has been associated with altered BPA metabolism and increased susceptibility to BPA-induced mitochondrial dysfunction [ 47 ]. This finding provides a molecular basis for understanding individual differences in the biological response to BPA exposure. Molecular docking analysis indicated that BPA binds strongly to NDUFA11 and NDUFV1, suggesting a direct interference with Complex I structure and electron transport. This aligns with increasing biological evidence showing that BPA disrupts mitochondrial function in neuronal and other tissues. In vivo studies have demonstrated that BPA exposure reduces mitochondrial membrane potential, impairs oxidative phosphorylation and decreases ATP production, changes accompanied by excessive mitochondrial ROS generation and oxidative damage [ 24 , 48 ]. These effects are further supported by neuronal cell experiments showing disrupted cristae structure, mitochondrial fragmentation and inhibited mitophagy, leading to accumulation of dysfunctional mitochondria and enhanced oxidative stress [ 49 ]. Under ischemic–hypoxic conditions relevant to stroke, neurons are highly dependent on intact Complex I to sustain energy metabolism [ 50 ]. BPA-induced destabilization of NDUFA11 and NDUFV1 may therefore act as an additional metabolic stressor, further weakening electron transport and accelerating ATP depletion. Several animal studies have reported that BPA aggravates ischemia–reperfusion injury, increases oxidative burden and amplifies neuronal loss, indicating that BPA functions as a secondary insult that exacerbates mitochondrial vulnerability during ischemic events [ 51 , 52 ]. Consistent with our findings, the predicted reduction in ATP synthesis and increased ROS accumulation caused by BPA-mediated Complex I inhibition provide a mechanistic explanation for the intensified mitochondrial dysfunction observed in stroke models. Recent studies showed that molecular docking applications for BPA have considerably broadened beyond classical endocrine receptors. While many studies still examine BPA–estrogen receptor and other nuclear receptor interactions, a growing body of dock-and-network literature has targeted metabolic regulators (e.g., PPAR family members and AKT1), signaling kinases and growth factors (e.g., AKT1, SRC, IGF1), matrix-remodeling enzymes (e.g., MMP9), and a range of disease-relevant proteins identified by network toxicology approaches [ 53 , 54 ]. For example, recent network-docking studies have reported stable BPA docking to AKT1, PPARG/PPARA and MMP9 in the context of diabetic cardiomyopathy (binding energies − 5.3 to − 7.5 kcal/mol) [ 55 ]. Similarly, integrated bioinformatics and docking analyses in endometriosis and renal cell carcinoma research prioritized hub proteins such as HSP90AA1, AKT1, SRC, IGF1, CHRM3 and GABBR1 as putative BPA targets, with Vina scores commonly below − 6.5 kcal/mol for the top hits [ 56 , 57 ]. These examples demonstrate that BPA has been computationally evaluated against a wide spectrum of protein families relevant to metabolism, inflammation, cell-cycle control and extracellular matrix regulation. Expanding the catalog of putative BPA–protein interactions beyond classical hormone receptors and underlining the chemical’s pleiotropic interaction potential. Despite this diversification, docking efforts explicitly directed at mitochondrial respiratory chain proteins — and Complex I subunits in particular — remain limited in the published literature. Therefore, our demonstration of sub-7 kcal/mol binding of BPA to NDUFA11 and NDUFV1 addresses a relative gap by proposing a direct mitochondrial target mechanism that complements the better-characterized receptor-mediated pathways. This study has several limitations. (1) The main conclusions are derived from computational simulations and bioinformatics predictions, which have not yet been experimentally validated. (2) Although the transcriptomic data were obtained directly from human brain tissue, the samples comprised multiple cell types, making it impossible to distinguish target expression patterns in specific populations such as neurons or glial cells. (3) The analysis primarily focused on the direct effects of BPA on mitochondria, without fully considering the indirect effects mediated by nuclear receptors such as estrogen receptors or potential multi-target synergistic mechanisms. (4) This study focused exclusively on BPA and did not include structurally related analogues such as BPS or BPF. Although these substitutes are increasingly used, current evidence regarding their exposure levels, toxicokinetic behavior, and mechanistic toxicity remains limited, making it unclear whether their biological effects resemble or diverge from those of BPA. Future studies should therefore aim to experimentally verify the regulatory effects of BPA on the identified core targets and assess its impact on mitochondrial function in both in vivo and in vitro models. In addition, single-cell sequence or spatial transcriptomic techniques could be employed to delineate the cell type–specific expression and functional characteristics of these targets within brain tissue. Finally, expanding the analytical framework to include bisphenol analogues and integrating multi-omics data may help to systematically elucidate the full spectrum of direct and indirect pathways through which bisphenols contribute to neuronal injury. In summary, this study moves beyond traditional toxicological observations to systematically elucidate the potential molecular link between BPA exposure and stroke, centering on mitochondrial energy metabolism revealed through multi-omics analysis. This work provides a novel molecular perspective and experimental basis for understanding how environmental pollutants contribute to cerebrovascular diseases and identifies several potential therapeutic targets. These findings establish a theoretical foundation for developing strategies to prevent and mitigate BPA-induced neurotoxicity, including approaches to preserve mitochondrial function and identify biomarkers of susceptibility. Based on these insights, we recommend a more comprehensive evaluation of BPA’s environmental safety, with particular emphasis on the potential cerebrovascular risks associated with chronic low-dose exposure.

Bisphenol A (BPA), a typical industrial compound, is widely applied in manufacturing polycarbonate plastics and epoxy resins—materials that are further processed into daily necessities ranging from food storage containers and infant feeding bottles to feminine care products, thermal printing paper, and household electronic devices [ 1 ]. As a highly prevalent synthetic endocrine disruptor, BPA can either mimic the activity of endogenous hormones or act as their antagonist, ultimately disrupting the homeostasis of the endocrine regulatory system [ 2 ]. Human exposure to BPA occurs through ingestion, dermal absorption, and inhalation [ 3 ]. Ingestion is the primary route, as BPA migrates from food-contact materials into food and beverages. Polyethylene terephthalate (PET) bottled water contains an average of 0.0055 ± 0.0022 mg/kg BPA [ 4 ], and storage at higher temperatures increases migration from 5.19 ng/L to 11.30 ng/L [ 5 ]. Other packaged drinking water also contributes to exposure [ 6 ]. BPA is common in food substrates, with levels reaching 0.089 mg/kg in canned shrimp [ 7 ] and appearing in 36% of non-canned foods at an average of 4.73 µg/kg [ 8 ], supporting food as a major exposure medium [ 9 ]. Skin absorption represents another pathway, as thermal paper used for store receipts contains substantial BPA that can transfer through contact [ 10 ]. Inhalation also contributes, BPA present in indoor air can be taken up through the respiratory tract, especially in confined spaces where BPA-containing products are used. Water environments form an additional source, with concentrations up to 1465 ng/L in wastewater and up to 3113 ng/L in surface water [ 11 ]. Research findings indicate that excessive BPA exposure can impair multiple physiological systems, including the reproductive, endocrine, cardiovascular, cerebrovascular, and nervous systems. Of particular concern is BPA-mediated neurotoxicity, which has garnered growing research interest due to its close link to the development and progression of both cerebrovascular disorders and neurodegenerative conditions—posing a significant risk to human well-being [ 12 ]. However, current research is mostly focused on chronic neurotoxicity, neurodevelopmental abnormalities, or neurodegenerative diseases, while the potential mechanisms of its involvement in acute central nervous system injury events, especially in acute neurological emergencies characterized by tissue necrosis such as stroke, remain underexplored. Stroke refers to an acute neurological condition triggered by interrupted blood flow to the brain. It features sudden onset, quick progression, and frequently unfavorable outcomes, and stands as a major factor behind adult disability and death across the globe [ 13 ]. As stated by the World Health Organization, over ten million individuals suffer from stroke annually, with approximately one-third of patients dying within one year of onset [ 14 ]. Survivors frequently experience long-term cognitive and motor impairments [ 15 , 16 ]. The core pathological mechanisms of stroke involve a coordinated series of events, including energy metabolism disruption, mitochondrial dysfunction, inflammatory cascades, and programmed cell death [ 17 , 18 ]. Notably, mitochondrial dysfunction stands out as a key feature: research indicates that in stroke models, mitochondria undergo critical functional impairments, including collapsed membrane potential, disrupted electron transport chain activity, and depleted ATP levels. These abnormalities work together to trigger apoptotic signaling cascades, eventually resulting in neuronal loss [ 19 ]. Stroke is the result of the synergistic effect of multiple factors, among which modifiable factors dominate, including underlying cardiovascular diseases and unhealthy lifestyles. Effective intervention targeting these factors could significantly reduce the risk of onset. Additionally, emerging research over the past few years has highlighted that environmental exposure factors—particularly long-term contact with neurotoxic substances—might contribute significantly to the onset and progression of stroke [ 20 ]. As a well-recognized environmental neurotoxicant, BPA has been shown to disrupt mitochondrial function [ 21 ]; such disruption includes direct changes to mitochondrial membrane permeability [ 22 ], suppression of respiratory chain key complex activity, stimulation of excessive reactive oxygen species (ROS) generation, and initiation of apoptotic or necrosis-like cellular processes [ 23 , 24 ]. However, despite BPA having clearly been established to possess neurotoxicity, systematic toxicological studies and integrative molecular mechanism analyses targeting BPA remain relatively scarce [ 25 ]. BPA also remains the most widely used bisphenol, and biomonitoring studies consistently show that BPA is detected at higher levels than its analogues such as BPS and BPF in human samples [ 26 , 27 ]. Although BPS and BPF are increasingly used as BPA substitutes, available data indicate that their toxicokinetic and mechanistic evidence remain less comprehensive compared with BPA [ 28 ]. Furthermore, traditional neurotoxicity assessment methods have been unable to fully elucidate the mechanisms by which environmental factors act in complex disease models. With the development of computational toxicology, systems biology, and big data mining technologies, an integrative research model has gradually become popular: that is, based on network toxicology [ 29 ], combining high-throughput transcriptomic data with molecular simulation techniques to construct a network of interactions between toxicants, targets, pathways, and diseases. This study aims to clarify how BPA contributes to acute injury of the central nervous system, with a particular focus on stroke. We integrate toxicological prediction, analysis of stroke-related transcriptome datasets and multi-level computational validation to identify the molecular events underlying BPA-related neurotoxicity. A weighted gene co-expression network will be constructed using Gene Expression Omnibus (GEO) stroke datasets, and principal component analysis (PCA) will be applied to distinguish disease-related expression patterns. Key molecules and pathways will be identified through protein–protein interaction (PPI) networks and Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses. The interaction between BPA and core targets will be further assessed through molecular docking. Together, these approaches will help define the specific mechanisms through which BPA influences stroke pathology and will provide a basis for developing targeted intervention strategies.