Abstract The toxic side effects of acetyl tributyl citrate (ATBC) on humans are concerning, but studies related to its effects on osteoarthritis (OA) are lacking. Therefore, this study aimed to explore the potential targets and mechanisms of action of ATBC in OA through network toxicology. We obtained ATBC-related targets from the ChEMBL, Swiss Target Prediction, and STITCH databases and OA-related targets from the GeneCards, DisGeNET, and OMIM databases and identified overlapping targets. Core targets (key molecules in the progression of diseases) were determined via the STRING database and Cytoscape software, followed by further Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to determine potential mechanisms in depth. Moreover, a gene interaction and competing endogenous RNA (ceRNA) network for the core targets was constructed. Additionally, the expression levels of the core targets were preliminarily validated using single-cell data from the GEO database. Furthermore, in-depth validation of the core targets was carried out through molecular docking and molecular dynamics simulations. A total of 132 overlapping targets between ATBC and OA were identified, and six core targets (TP53, EZH2, HDAC1, HDAC2, SIRT1, and SMARCA4) were further screened. The results of the enrichment analysis revealed that the core pathways related to the effect of ATBC on OA mainly involved key signaling cascades, including the thyroid hormone signaling pathway, the Notch signaling pathway, and cellular senescence. Single-cell analysis revealed that the core target is expressed in different cell subpopulations. Molecular docking and molecular dynamics simulation results indicate that there is a stable binding interaction between ATBC and the core target. This study provides a theoretical foundation for the molecular mechanisms of OA triggered by ATBC, highlighting the value of network toxicology in assessing the toxicity of emerging environmental pollutants. However, further clinical and experimental investigations are needed to validate these findings. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-11178-5. Keywords: ATBC, Osteoarthritis, Network toxicology, Molecular Docking, Core targets Subject terms: Drug discovery, Molecular biology Introduction The production and utilization of industrial chemicals have increased dramatically over the past few decades, leading to the emergence of many new types of environmental pollutants. Among them, plastics have become an indispensable part of the daily lives of people because of their versatility, durability, and utility^[38]1. However, substances in plastics may also have toxic effects on the environment and organisms. These substances may disrupt the endocrine system of living organisms^[39]2triggering an inflammatory response and affecting or even exacerbating symptoms of osteoarthritis (OA). Plasticizers are widely used to soften plastics and are mainly classified into phthalates and nonphthalates^[40]3. Acetyl tributyl citrate (ATBC) has been widely used as an alternative plasticizer in many fields, such as personal care products, medical equipment, and laboratory tubing, and has replaced phthalates^[41]4. ATBC is a safe plasticizer at typical concentrations, but with the prolonged exposure of humans to ATBC in the environment, its hazards to humans need to be evaluated. The high-mobility characteristics of ATBC also indicate that it may be hazardous to humans^[42]5. ATBC has a migration value of 79.8 mg/kg, which greatly exceeds the migration value of phthalates, leading to leakage into the environment from some products that contain it^[43]6which may be associated with the lack of strong chemical bonds between ATBC and polymers in plastics^[44]7. Therefore, ATBC may not be very safe. For example, ATBC might be a potential environmental obesogen with the risk of interfering with metabolism and inducing fatty liver disease^[45]8. Exposure to 20 mg/kg/d ATBC exacerbates disorders of glucose and lipid metabolism and causes cognitive deficits in type 2 diabetes mellitus (T2DM) mice^[46]9; at the cellular level, ATBC can significantly inhibit the proliferation and migration of osteoclasts and disrupt bone metabolism homeostasis at low doses^[47]7. Based on the effects of ATBC on specific levels of estrogen, androgen, and thyroid hormones and considering that estrogen has a significant protective effect on OA, we hypothesized that ATBC may indirectly promote OA by disrupting its normal function. OA is a disease of the joint, with lesions in the articular cartilage, subchondral bone, and synovial membrane, and eventually progresses to joint failure^[48]10. OA not only affects the health of individuals but also places an incredible burden on the healthcare system. The annual direct medical cost of OA in the United States between 2008 and 2011 was about $339.7 billion^[49]11. Early prevention of OA is essential because current treatment options do not prevent its progression. In OA, risk factors can be categorized into those that act on the susceptibility level of individuals and those that modify the biomechanical stability of the joints. Studies have shown a significant increase in the incidence of OA in populations exposed to certain plastic contaminants in the environment, such as perfluorooctanoates and perfluorooctane sulfonates^[50]12. People are exposed to plastic-containing pollutants for prolonged periods, but few studies have investigated whether plastic pollutants in the environment increase the risk of OA. Plasticizers are an integral part of plastics, and research on the role of ATBCs in OA is lacking. Therefore, the mechanism underlying the toxic effects of ATBC on OA needs to be elucidated to support the early prevention and subsequent diagnosis and treatment of OA. In the last few decades, several industrial chemicals have been synthesized and used, and many new environmental pollutants have emerged, the toxicological characteristics of which are largely unknown^[51]13. Conventional toxicological methods are technically and instrumentally inadequate for fully assessing the complex biological effects of these contaminants. Therefore, toxicological research methods need to be developed to assess the human health risks of many uncharacterized environmental pollutants effectively. As an emerging scientific field, the core concept of network toxicology involves the use of multiple databases to structure networks of relationships between compounds, toxicities, and targets^[52]4. Computer modeling and bioinformatics are used to construct relational networks and convert complex mechanisms into graphical models. Network toxicology can help us predict the toxicity of chemical substances and molecular mechanisms and facilitate the screening of potential toxic substances because of its simple operation, accurate results, and wide application^[53]14. Therefore, this study aimed to explore the potential targets and mechanisms of action by which ATBC affects OA through network toxicology, providing theoretical support for the further diagnosis and treatment of OA related to ATBC exposure. Materials and methods The overall research design flowchart is illustrated in Fig. [54]1. The specific process is as follows: the candidate targets of ATBC were retrieved from ChEMBL, SwissTargetPrediction, and STITCH databases, while OA-related targets were obtained from GeneCards, DisGeNET, and OMIM databases. A total of 132 overlapping genes were identified as potential targets. These were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis and to PPI network construction. Core targets (n = 6) were identified using cytoHubba. Further analysis was conducted using gene-gene and competing endogenous RNA (ceRNA) networks, as well as single-cell RNA sequencing data. The core targets were also evaluated through molecular docking and molecular dynamics simulations. Fig. 1. [55]Fig. 1 [56]Open in a new tab The flowchart of target identification and analysis strategy for ATBC affecting OA. Collection of ATBC targets The standard structure and SMILES symbol of ATBC were determined via the query “acetyl tributyl citrate” in the PubChem database ([57]https://pubchem.ncbi.nlm.nih.gov/, PubChem release 2021.10.14). Following this search, by using the keyword “acetyl tributyl citrate” and restricting the species to “Homo sapiens”, potential ATBC targets were obtained from the ChEMBL database ([58]https://www.ebi.ac.uk/chembl/, ChEMBL ID: CHEMBL1904556). Moreover, the SMILE symbol ATBC was uploaded to the Swiss Target Prediction database ([59]http://www.swisstargetprediction.ch/) and the STITCH database ([60]http://stitch.embl.de/) to identify unnoticed potential targets. We subsequently integrated the targets from different sources, removed duplicate data, and verified structural consistency. Next, the UniProt database ([61]https://www.uniprot.org/) was used to normalize the names of the obtained targets: the original target list obtained from the ChEMBL, SwissTargetPrediction, and STITCH databases (including Gene Symbols, Entrez IDs, and other heterogeneous identifiers) was uniformly mapped using the batch retrieval module of UniProt. For targets that cannot be automatically matched, manual verification is employed, and targets that do not pass UniProt validation are removed. Ultimately, an extensive target library for ATBC was established, encompassing all consolidated and verified targets^[62]7. Selection of the OA-related target network The OA-related targets were retrieved from the following databases: GeneCards ([63]https://www.genecards.org/), DisGeNET ([64]https://www.disgenet.org/), and OMIM ([65]https://www.omim.org/). Next, we integrated the targets from different sources, removed duplicate data, and established a dedicated OA target library. Furthermore, a Venn diagram was generated to identify potential targets shared between ATBC and OA, and the overlapping section designed as potential targets through which ATBC affects OA^[66]15. Common target functions and pathway enrichment analysis To explore the biological functions of common targets influenced by ATBC in the process of OA, the R software (v 4.2.1) was used with the “clusterProfiler” (v 4.4.4) package to perform enrichment analyses on the common potential targets, the ID conversion library was org.Hs.eg.db, and the P value correction method was the Benjamini-Hochberg method. We conducted a comprehensive GO analysis, including biological process (BP), cellular component (CC), and molecular function (MF) terms. We also identified the pathways related to OA affected by ATBC through KEGG analysis^[67]16. Finally, the results were presented in the form of bar plots and bubble plots using the “ggplot2” (v 3.4.4) package. The bar chart integrates information on enrichment categories, significance levels, and enriched entries, whereas the bubble chart integrates information on gene ratios, significance levels, and gene counts^[68]17with significance indicated by adj. P < 0.05. Construction of the PPI network and core target selection The potential target genes of ATBC affecting OA were input into the STRING database ([69]https://string-db.org/). The species was limited to “Homo sapiens”, and the “minimum required interaction score” was set to “high confidence > 0.7” for analysis. The data produced by the STRING database were subsequently imported into Cytoscape software (v 3.9.1) for visualization via the default parameters. Next, the cytoHubba plugin was used to identify core targets based on four common algorithms: MCC, MNC, Degree, and Closeness. Finally, a Venn diagram was generated to determine the overlap of the top 10 ranked core targets from the four algorithms, and chromosomal localization analysis of the core targets was conducted utilizing the “circlize” (v 0.4.15) package. Core target function and pathway enrichment analysis To elucidate the potential mechanisms involved in the process by which ATBC affects OA, we performed GO and KEGG enrichment analyses of the core targets. Construction of gene interaction networks The GeneMANIA database ([70]https://genemania.org/) was utilized to pinpoint potential genes associated with the core targets, and a gene-gene interaction network was constructed^[71]18. The starBase database ([72]https://rnasysu.com/encori/) was used to set “CLIP-Data ≥ 4” for analysis. First, miRNAs for core targets were predicted; then, the qualified miRNAs were used as input to predict the corresponding lncRNAs. Next, the organized and standardized prediction results were imported into Cytoscape software to establish the mRNA-miRNA-lncRNA ceRNA regulatory network. Single-cell data validation Download the [73]GSE176308 and [74]GSE283080 datasets from the GEO database ([75]https://www.ncbi.nlm.nih.gov/geo/) for validation of core targets at the single-cell level. [76]GSE176308 includes 3 samples, derived from synovial fibroblasts, thus only cell subtypes are classified. The following preprocessing is performed on [77]GSE176308 (set the threshold according to the original article): individual sample data are merged, and low-quality or double cells are filtered by excluding cells with more than 7,100 expressed genes and more than 15% mitochondrial gene expression; principal component analysis (PCA) is conducted on highly variable genes, and cells are clustered through Uniform Manifold Approximation and Projection (UMAP) analysis with a resolution of 0.2^[78]19. [79]GSE283080 includes 4 samples (2 fibrous (fibro) and 2 inflammatory (infla)). The following preprocessing is performed on [80]GSE283080 (set the threshold according to the original article): individual sample data are merged, and low-quality or double cells are filtered by excluding cells with more than 2,500 expressed genes and more than 5% mitochondrial gene expression; PCA analysis is conducted on highly variable genes, and cells are clustered through UMAP analysis with a resolution of 0.5^[81]20. Marker genes for different cell subpopulations are presented using bubble plots, and the expression differences of core targets in different OA subtypes are analyzed for different cell subpopulations. Different datasets are classified according to OA subtypes: [82]GSE176308 samples are divided into early and endstage OA patients, while [83]GSE283080 samples are divided into fibro and infla OA patients. All analyses are conducted using the “Seurat” (v 4.0) package. Molecular Docking verification The core target protein structure was downloaded from the RCSB PDB database ([84]https://www.rcsb.org/), and the ATBC structure was obtained from PubChem. Next, PyMOL software (v 2.6) was employed to eliminate water molecules and the original ligand for protein preprocessing, which was saved as a PDB file. The online molecular docking website CB-Dock2 ([85]https://cadd.labshare.cn/cb-dock2/php/index.php) was used for online docking. Finally, PyMOL software was used to compare and visualize the mode of binding of the core target protein with ATBC. Molecular dynamics simulation verification Molecular dynamics simulations were performed via Gromacs software (v 2022.3). Small molecules were pretreated via AmberTools22 to add a GAFF force field to the small molecules, Gaussian16W was used for hydrogenation and calculation of restrained electrostatic potential potential operation, and the potential data were added to the molecular dynamics system topology file. The simulation conditions are carried out at static temperature (300 K) and atmospheric pressure (1 bar), the force field is aber99sb-ildn, the solvent is a water molecule (TIP3P water model), and the total charge of the simulated system is neutralized by the addition of an appropriate amount of Na^+^[86]21. The molecular dynamics simulation system was first subjected to energy minimization using the steepest descent method, followed by 100,000 steps of constant number of particles, volume, and temperature (NVT) equilibrium and constant-pressure, constant-temperature (NPT) equilibrium, respectively. Finally, a free molecular dynamics simulation was run with 5,000,000 steps of 2 fs duration for a total of 100 ns. After the simulation was complete, the trajectories were analyzed via software to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), radius of gyration (Rg), free energy of binding (MM/GBSA), free energy landscape (FEL) and other data. Results Identification of the targets affected by ATBC in OA A total of 362 ATBC-related targets were acquired by integrating the ChEMBL, Swiss Target Prediction, and STITCH databases; 5,290 OA-related targets were obtained by integrating the GeneCards, DisGeNET, and OMIM databases. Subsequently, 132 common cross-targets shared between ATBC and OA were discerned (Fig. [87]2, Supplemental Table [88]1), and these overlapping targets were potential candidates for the effect of ATBC on OA. Fig. 2. Fig. 2 [89]Open in a new tab Venn diagram of the targets of ATBC and OA. GO and KEGG enrichment analyses of common potential targets Enrichment analysis was conducted on 132 common potential targets, resulting in 1,390 significant GO terms, including 1,194 BP terms, 62 CC terms, and 134 MF terms. The GO terms were systematically ranked on P value, and the top 10 terms exhibiting the most notable statistical disparities in each category the 10 most significant terms in each category were selected for visualization (Fig. [90]3A and B). Within the BP classification, the targets were enriched mainly in histone modification (adj. P = 7.77 × 10^–40), peptidyl-lysine modification (adj. P = 3.37 × 10^–23), and histone methylation (adj. P = 1.66 × 10^–19); within the CC classification, the targets were enriched mainly in chromosomes and telomeric regions (adj. P = 3.47 × 10^–13), transcription regulator complexes (adj. P = 2.40 × 10^–12), and histone deacetylase complexes (adj. P = 4.12 × 10^–12), and within the MF classification, the targets were enriched mainly in histone deacetylase activity (adj. P = 9.88 × 10^–17), protein lysine deacetylase activity (adj. P = 9.88 × 10^–17), and transcription coregulator activity (adj. P = 1.40 × 10^–16). Fig. 3. [91]Fig. 3 [92]Open in a new tab GO and KEGG enrichment analyses of common potential targets. (A, C) Bar plot of the results of the GO and KEGG enrichment analyses. (B, D) Bubble diagram of the results of the GO and KEGG enrichment analyses (bar plot: different colors represent different categories, the x-axis indicates significance, and the y-axis represents enrichment result terms; bubble diagram: the color of the bubbles represents significance, the size represents quantity, the x-axis indicates the gene ratio, and the y-axis represents enrichment result terms). Moreover, 99 signaling pathways with significant differences were identified, and the 20 entries with the most significant differences were chosen for visualization (Fig. [93]3C and D). The results revealed that the common potential targets of ATBC and OA were associated with various pathways, such as the FoxO signaling pathway (adj. P = 7.89 × 10^–7), microRNAs in cancer (adj. P = 7.89 × 10^–7), and the thyroid hormone signaling pathway (adj. P = 1.57 × 10^–6). Interaction network of common potential targets and acquisition of core targets A PPI network with 109 nodes and 547 edges was generated utilizing the STRING database (Fig. [94]4A). The top 10 potential targets ranked by MCC (Fig. [95]4B), MNC (Fig. [96]4C), Degree (Fig. [97]4D), and Closeness (Fig. [98]4E) were generated via the cytoHubba plugin. Finally, six overlapping targets (TP53, EZH2, HDAC1, HDAC2, SIRT1, and SMARCA4) were identified by intersecting results from the four algorithms (Fig. [99]5; see Table [100]1 for relevant information). Fig. 4. [101]Fig. 4 [102]Open in a new tab (A) The PPI network of common potential targets. Core target network based on the MCC (B), MNC (C), Degree (D) and Closeness (E) parameters. Fig. 5. [103]Fig. 5 [104]Open in a new tab (A) Venn diagram of the four algorithms. (B) Chromosomal localization map of core targets. Table 1. Information related to core targets. Core targets Full name Chromosome PDB ID Vina score TP53 Tumor protein P53 chr17 8SVI −4.9 EZH2 Enhancer of zeste homolog 2 chr7 5U5T −6.3 HDAC1 Histone deacetylase 1 chr1 7Z1K −4.5 HDAC2 Histone deacetylase 2 chr6 7KBG −4.9 SIRT1 Sirtuin 1 chr10 1YC5 −6.6 SMARCA4 SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily A member 4/brahma-related gene 1 chr19 7TAB −4.4 [105]Open in a new tab GO and KEGG enrichment analysis of core targets Enrichment analysis was conducted on six core targets, and the results revealed 779 significant GO terms, comprising 688 BP terms, 30 CC terms, and 61 MF terms. As shown in Figs. [106]6A and B, within the BP classification, the targets were enriched mainly in circadian rhythm (adj. P = 7.98 × 10^–7), heterochromatin organization (adj. P = 1.72 × 10^–6), and histone deacetylation (adj. P = 1.87 × 10^–6); within the CC classification, the targets were enriched mainly in heterochromatin (adj. P = 2.14 × 10^–5), the SWI/SNF superfamily type complex (adj. P = 2.14 × 10^–5), and the transcription regulator complex (adj. P = 8.02 × 10^–4). Within the MF classification, the targets were enriched mainly in promoter-specific chromatin binding (adj. P = 1.89 × 10^–10), p53 binding (adj. P = 6.44 × 10^–8), and nicotinamide adenine dinucleotide (NAD)-dependent histone deacetylase activity (adj. P = 1.82 × 10^–7). Fig. 6. [107]Fig. 6 [108]Open in a new tab GO and KEGG enrichment analyses of core targets. (A, C) Bar diagram of the results of the GO and KEGG enrichment analyses. (B, D) Bubble diagram of the GO and KEGG enrichment analyses. Moreover, 17 signaling pathways with significant differences were enriched. As shown in Fig. [109]6C and D, the core targets of ATBC and OA were related mainly to signaling pathways involving microRNAs in cancer (adj. P = 2.76 × 10^–5), longevity-regulating pathways in multiple species (adj. P = 2.35 × 10^–4), the thyroid hormone signaling pathway (adj. P = 7.17 × 10^–4), the cell cycle (adj. P = 7.17 × 10^–4), the Notch signaling pathway (adj. P = 4.69 × 10^–3), and cellular senescence (adj. P = 2.29 × 10^–2). Gene-gene interaction network The gene-gene interaction network of the core targets is shown in Fig. [110]7, with the central hub nodes symbolizing these genes, surrounded by 20 nodes (SMYD1, CDK4, EED, HSPA9, SMARCB1, NFKBIA, RERE, SUZ12, MAP2K1, MAP1LC3B, NR2F2, RUNX2, MCM10, PPARD, EIF3D, RBBP4, NCOA1, HIF1A, CCAR2, and RUNX1T1), representing genes significantly associated with them. The size of the node represents the degree of correlation, with larger nodes indicating greater correlations. The network highlights the top seven related functions associated with these genes, such as chromatin binding, chromatin organization involved in the regulation of transcription, and protein deacetylase activity, among others. Fig. 7. [111]Fig. 7 [112]Open in a new tab Gene-gene interaction network. CeRNA network The ceRNA regulatory network of mRNAs-miRNAs-lncRNAs for the core targets is shown in Fig. [113]8. Through miRNA prediction, 26 miRNAs (hsa-miR-218-5p, hsa-miR-3163, hsa-miR-27b-3p, hsa-miR-513a-5p, hsa-miR-27a-3p, hsa-miR-363-3p, hsa-miR-92a-3p, hsa-miR-33a-5p, hsa-miR-32-5p, hsa-miR-217, hsa-miR-25-3p, hsa-miR-137, hsa-miR-1914-3p, hsa-miR-5194, hsa-miR-3126-5p, hsa-miR-125b-5p, hsa-miR-125a-5p, hsa-miR-18b-5p, hsa-miR-874-3p, hsa-miR-494-3p, hsa-miR-128-3p, hsa-miR-2355-5p, hsa-miR-624-5p, hsa-miR-556-3p, hsa-miR-22-3p, and hsa-miR-942-5p) were found to be associated with six core targets. Subsequently, lncRNA prediction was performed on the identified miRNAs, resulting in the discovery of 55 lncRNAs that may associate with the miRNAs. Fig. 8. [114]Fig. 8 [115]Open in a new tab A ceRNA network of mRNAs-miRNAs-lncRNAs was constructed based on core targets; the orange circles represent mRNAs, the green circles represent miRNAs, and the blue circles represent lncRNAs. Single-cell data validation After quality control, the dataset [116]GSE176308 was annotated to obtain 7 subtypes of fibroblasts, defined as subtypes C0-C6. The expression of core targets was analyzed in the dataset, and the results showed that TP53, EZH2, HDAC1, HDAC2, SIRT1, and SMARCA4 were expressed in different cell subpopulations of early and endstage OA patients (Fig. [117]9, Supplementary Fig. [118]1 A, B). After quality control, the dataset [119]GSE283080 was annotated to obtain 8 cell types, and the expression of core targets was analyzed in the dataset, revealing that TP53, EZH2, HDAC1, HDAC2, SIRT1, and SMARCA4 were expressed in different cell subpopulations of fibro and infla OA patients (Fig. [120]10, Supplementary Fig. [121]1 C, D). Fig. 9. [122]Fig. 9 [123]Open in a new tab [124]GSE176308 dataset single-cell validation. (A) Cell annotation UMAP plot. (B) Bubble chart of marker genes for different cell subpopulations. (C) UMAP plot of core targets validated in early OA patients. (D) UMAP plot of core targets validated in endstage OA patients. Fig. 10. [125]Fig. 10 [126]Open in a new tab [127]GSE283080 dataset single-cell validation. (A) Cell annotation UMAP plot. (B) Bubble chart of marker genes for different cell subpopulations. (C) UMAP plot of core targets validated in fibro OA patients. (D) UMAP plot of core targets validated in infla OA patients. Molecular Docking of ATBC with core target proteins The interaction between ATBC and six core target proteins was investigated through molecular docking technology (Fig. [128]11). The findings revealed that the binding energy of the six target proteins ranged from − 4.4 to − 6.1 kcal/mol: TP53 (–5.2 kcal/mol, 1548Ser, 1554Ser, 1577Gln exist hydrogen bonds, TYR1552 exists π-bonds), EZH2 (–6.1 kcal/mol, 313Arg, 100Gln exist hydrogen bonds, 103Trp exists π-bonds), HDAC1 (–4.4 kcal/mol, 3510Gln, 3541Ser exist hydrogen bonds), HDAC2 (–4.8 kcal/mol, 206Phe, 179His exist hydrogen bonds), SIRT1 (–5.7 kcal/mol, 100Ile, 98Gln exist hydrogen bonds), and SMARCA4 (–4.4 kcal/mol,1540Asn exists hydrogen bonds). All six core target proteins exhibited strong binding with ATBC, with binding energies < − 4.25 kcal/mol, which indicated that ATBC can spontaneously bind to each core target^[129]22 and implied that they play a key role in the molecular mechanism by which ATBC affects OA. Fig. 11. [130]Fig. 11 [131]Open in a new tab Molecular docking of ATBC with its core targets and the binding energy of the six complexes: (A) TP53, − 5.2 kcal/mol, (B) EZH2, − 6.1 kcal/mol, (C) HDAC1, − 4.4 kcal/mol, (D) HDAC2, − 4.8 kcal/mol, (E) SIRT1, − 5.7 kcal/mol, and (F) SMARCA4, − 4.4 kcal/mol. Molecular dynamics simulation verification Molecular dynamics simulations were performed to validate the molecular docking simulation results of ATBC with the core targets. The RMSD curves of all the complexes reached equilibrium after 20 ns and remained relatively stable thereafter (Fig. [132]12). RMSF analysis revealed key flexible regions in each core target protein, and the RMSF values of most residues fluctuated in the range of 0.05–0.20 nm, indicating structural stability (Fig. [133]13). The corresponding SASA value gradually decreases with time, indicating that the system tends toward a tighter conformation during the simulation (Fig. [134]14). The Rg values of the six complexes remained relatively balanced, indicating that the overall shape of the system remained relatively stable throughout the simulation without significant unfolding or contraction (Fig. [135]15). After the complex system stabilized, we calculated the average binding free energies of the six groups of proteins with ATBC via the MM/GBSA method. As shown in Table [136]2, the average free energies of TP53, EZH2, HDAC1, HDAC2, SIRT1, and SMARCA4 binding to ATBC were − 26.10 kcal/mol, − 43.66 kcal/mol, − 26.89 kcal/mol, − 22.55 kcal/mol, − 43.67 kcal/mol and − 30.25 kcal/mol, respectively; the results show that all six core targets have good binding energies and high binding affinities with ATBC, suggesting that the binding of the six core targets to ATBC is thermodynamically favorable, driven mainly by hydrophobic effects and electrostatic interactions, further validating the molecular docking results. Fig. 12. [137]Fig. 12 [138]Open in a new tab The RMSD value of the six complexes provides insights into their overall stability and conformational changes: (A) TP53-ATBC, (B) EZH2-ATBC, (C) HDAC1-ATBC, (D) HDAC2-ATBC, (E) SIRT1-ATBC, and (F) SMARCA4-ATBC. Fig. 13. [139]Fig. 13 [140]Open in a new tab The RMSF value of the six complexes highlights the flexibility of residues: (A) TP53-ATBC, (B) EZH2-ATBC, (C) HDAC1-ATBC, (D) HDAC2-ATBC, (E) SIRT1-ATBC, and (F) SMARCA4-ATBC. Fig. 14. [141]Fig. 14 [142]Open in a new tab The SASA value tracks the exposure of the six complexes to the solvent: (A) TP53-ATBC, (B) EZH2-ATBC, (C) HDAC1-ATBC, (D) HDAC2-ATBC, (E) SIRT1-ATBC, and (F) SMARCA4-ATBC. Fig. 15. [143]Fig. 15 [144]Open in a new tab The Rg value examines the compactness of the six complexes: (A) TP53-ATBC, (B) EZH2-ATBC, (C) HDAC1-ATBC, (D) HDAC2-ATBC, (E) SIRT1-ATBC, and (F) SMARCA4-ATBC. Table 2. MM-GBSA parameters for the six protein-ligand complexes. Energy TP53-ATBC EZH2-ATBC HDAC1-ATBC HDAC2-ATBC SIRT1-ATBC SMARCA4-ATBC ΔVDWAALS −34.65 ± 1.39 −54.05 ± 2.69 −33.27 ± 2.08 −32.18 ± 0.32 −55.04 ± 0.68 −38.86 ± 0.08 ΔE[elec] −18.32 ± 0.57 −31.21 ± 0.20 −13.32 ± 1.96 −8.92 ± 0.19 −12.89 ± 0.75 −3.45 ± 1.25 ΔE[GB] 32.50 ± 1.91 49.02 ± 0.52 23.99 ± 0.07 22.96 ± 0.34 32.19 ± 1.57 17.54 ± 0.18 ΔE[surf] −5.63 ± 0.08 −7.42 ± 0.00 −4.29 ± 0.10 −4.42 ± 0.01 −7.94 ± 0.07 −5.47 ± 0.04 ΔG[gas] −52.97 ± 1.51 −85.26 ± 2.70 −46.59 ± 2.86 −41.10 ± 0.37 −67.92 ± 1.01 −42.31 ± 1.25 ΔG[solvation] 26.87 ± 1.91 41.60 ± 0.52 19.71 ± 0.12 18.55 ± 0.34 24.25 ± 1.57 12.06 ± 0.18 ΔTotal −26.10 ± 2.43 −43.66 ± 2.75 −26.89 ± 2.86 −22.55 ± 0.50 −43.67 ± 1.86 −30.25 ± 1.27 [145]Open in a new tab ΔVDWAALS: Van der Waals energy; ΔE[elec]: electrostatic energy; ΔE[GB]: polar solvation energy; ΔE[surf]: non-polar solvation energy; ΔG[gas]: molecular mechanical term energy (meteorological chemical energy) = ΔVDWAALS + ΔE[elec]; ΔG[solvation]: solvation energy = ΔE[GB] + ΔE[surf]; ΔTOTAL: comprehensive total energy = ΔG[gas] + ΔG[solvation]. To provide rich visual information, as well as a deeper understanding of the dynamic behavior and stability of the complexes, we plotted the RMSD values of the small molecules docked to the protein receptor complexes and the Rg values derived from the PCA of the PC1 and PC2 axes, combined with the Gibbs relative free energy along the Z-axis, to create a three-dimensional FEL (Fig. [146]16), which was plotted in combination with the Gibbs relative free energy along the Z-axis. The region with a blue hue in the plot prominently marks the steady-state conformations of the complexes, which are clustered at lower energies within the minimum free energy region. The smoothness and depth of this region reflect the strength of the complex stability. If the interaction between the protein and the ligand is weak or unstable, the FEL map exhibits multiple rough-surface minimum energy clusters. However, in our study, the FEL plots of the small molecule-protein receptor docking complexes in each system exhibited nearly single and sharp minimum energy clusters. This finding indicates that the complexes maintain a highly stable state during the simulation. These results are consistent with the RMSD curves of the complexes, further confirming that all the complexes have excellent stability. This stability is crucial for the study of the potential risk effect of ATBC on OA, as it implies that the complexes can maintain their structural and functional integrity under physiological conditions, thus exerting their biological activities more efficiently. Fig. 16. [147]Fig. 16 [148]Open in a new tab Free energy landscape, the blue and green regions represent low-energy conformational states, suggesting stable configurations reveals the six complexes remain highly stable throughout the simulation: (A) TP53, (B) EZH2, (C) HDAC1, (D) HDAC2, (E) SIRT1, and (F) SMARCA4. Discussion In this study, network toxicology methods were used, integrating data from the ChEMBL, Swiss Target Prediction, STITCH, GeneCards, DisGeNET, and OMIM databases, and 132 common potential targets associated with ATBC exposure and OA were systematically identified. A PPI network was constructed via the STRING database and Cytoscape software, which identified six core targets that play key roles via the impact of ATBC on OA, including TP53, EZH2, HDAC1, HDAC2, SIRT1, and SMARCA4, providing valuable insights for a deeper understanding of their potential mechanisms. TP53 (tumor protein P53), also known as p53, has been implicated in cell cycle arrest and apoptosis during disease development and is a key regulator of cellular activity. Targeted modulation of the STAG1/TP53/P21 signaling pathway may be a potential therapeutic strategy to delay cartilage tissue degeneration in OA patients^[149]23. The mechanism of action involves significant upregulation of TP53 expression levels in the pathological state of OA, and available studies suggest that this phenomenon may be related to the promotion of IL-38 expression upregulation by TP53, which mediates anti-inflammatory effects^[150]24. Some studies have also shown that the secretory phenotype implicated in OA cartilage senescence can be attenuated by inhibiting the STAG1/TP53 pathway^[151]23. Additionally, the up-regulation of TP53 and IL-38 and the activation of IL-36R reduce inflammation and improve chondrocyte apoptosis in OA^[152]24. EZH2 (enhancer of zeste homologue 2) is a histone methyltransferase that catalyzes H3K27me3, thereby repressing gene transcription. EZH2 is highly expressed in OA chondrocyes and silences the Wnt-antagonist SFRP1 via H3K27me3. In turn, H3K27me3 activates the Wnt/β-catening pathway, accelerating cartilage matrix degradation^[153]25. EZH2 enhances bone metabolism and promotes the development of OA by regulating related genes in multiple ways^[154]26. EZH2 disorders accelerate cell proliferation and prolong cell survival, which may contribute to cancer development^[155]27. EZH2 can also regulate E2F1-dependent apoptosis through epigenetic regulation of Bim expression, thus inhibiting apoptosis^[156]28. HDAC1 (histone deacetylase 1) belongs to the family of histone deacetylases related to the regulation of histone acetylation^[157]29. HDAC1 expression is elevated in OA, and HDAC1 inhibits miR-146a expression in the synovium and exacerbates cartilage damage^[158]30. HDAC1/EEF2 regulates the activation of the NF-κB and MAPK pathways^[159]31. In double knockout mice, deletion of the HDAC1/2 gene causes sustained DNA damage, along with senescence and loss of podocytes, indicating that HDAC1/2 plays a key role in developmental processes^[160]32. By regulating HDAC1 expression, nuclear factor of activated T-cell 2 modulates the mesenchymal (MES) transition of glioma stem cells (GSCs)^[161]33. HDAC2 (histone deacetylase 2) is very similar to HDAC1 and plays a pro-inflammatory role in the development of OA^[162]34. A study revealed elevated HDAC2 gene expression in OA chondrocytes compared to controls^[163]35. This upregulation also applies to HDAC1 at the protein level. HDAC2 represses cartilage-specific genes in human chondrocytes^[164]36,[165]37. Inhibition of HDAC2 expression promotes chondrogenic matrix expression in chondrocytes and delays OA progression^[166]38. SIRT1 (sirtuin 1) acts as a NAD-dependent enzyme that participates in the regulation of various biological processes, such as cellular metabolism, senescence, and inflammatory responses. SIRT1 expression is reduced in human chondrocytes samples isolated from the knees of OA patients undergoing total knee arthroplasty. In addition, SIRT1 expression was reduced in mouse knee OA, and OA progressed more rapidly in mice with specific knockout of SIRT1 in cartilage^[167]39. SIRT1 reduces the expression of IL-1β-induced inflammatory factors by inhibiting the NF-κB signaling pathway, thus attenuating inflammatory responses^[168]39. Taken together, the evidence from these studies suggests that changes in the expression level of SIRT1 participate in the molecular regulatory network of OA pathological processes as key regulators. SMARCA4 (brahma-related gene 1) commonly regulates gene expression and affects cell differentiation and proliferation. The expression of miR-155 is notably increased in OA chondrocytes, and SMARCA4 is considered to be one of its target genes^[169]40. SMARCA4 plays a key role in the pathological process of OA as a stage-specific epigenetic regulatory hub during chondrocyte senescence. Pathological analyses confirmed that OA-mesenchymal stromal cells are the main effector cells involved in the formation of a pro-inflammatory microenvironment in cartilage tissues, whereas SMARCA4 maintains the cellular homeostatic balance by regulating the molecular cascade of molecular responses involved in the transformation of osteoarthritic chondrocytes to OA-MSCs. Notably, functional deletion of SMARCA4 can trigger multiple pathological effects: accelerating the senescence process of MES stromal cells, triggering the release of inflammatory factor cascades, and ultimately leading to cartilage matrix degradation and increased synovial inflammatory responses^[170]41. In our study, the KEGG analysis of the core targets revealed that pathways associated with Longevity regulating pathway - multiple species, Thyroid hormone signaling pathway, and the Notch signaling pathway, were significantly enriched. Longevity regulating pathway - multiple species pathway, plays important roles in various organisms, and our core target, SIRT1, is a very important part of this pathway. In OA, Longevity regulating pathway may influence disease progression by modulating cellular metabolism and stress responses. For example, the mTOR signaling pathway can regulate cell growth, metabolism, and proliferation and is related to chronic inflammatory and fibrotic diseases^[171]42. Thyroid hormones act intracellularly through their receptors to regulate many biological processes, such as cell differentiation, metabolism, and inflammatory responses. Studies suggest that thyroid hormones may facilitate the development of OA and that increased availability of 3,5,3′-L-triiodothyronine may be detrimental to joint maintenance and articular cartilage homeostasis^[172]43. Additionally, deiodinase iodothyronine type 2 (DIO2) is one of many OA susceptibility genes; changes in DIO2 activity in articular cartilage might upset the balance in vivo by promoting the differentiation of hypertrophic chondrocytes and other undesirable events that ultimately lead to the onset or progression of OA^[173]44. Abnormal thyroid hormone receptor α (THRα) signaling in osteoblasts may affect microangiogenesis in a portion of subchondral bone, and local inhibition of THRα may be a target for the treatment of OA^[174]45. Notch is a single-channel transmembrane cell surface receptor that plays an important role in cell fate allocation by regulating various developmental systems, such as differentiation and apoptosis during embryogenesis, as well as neurogenesis and hematopoiesis^[175]46. The role of the Notch pathway in chondrocyte extracellular matrix (ECM) homeostasis, inflammation, phenotypic changes, and hypoxia regulation in OA^[176]47 has been well studied. The role of Notch signaling in OA disease progression appears to be complex; it is repressive in the initial phase but stimulatory in the terminal phase during chondrocyte differentiation and endochondral ossification^[177]48. Moreover, RBPjκ-dependent Notch signaling in postnatal articular chondrocytes is required for articular cartilage and joint maintenance^[178]49. Additionally, aberrant transmission of hedgehog signaling in OA may be limited by activation of the Notch pathway^[179]50. Some studies suggest that circ_0104873 promotes OA by targeting the miR-875-5p/NOTCH3 axis^[180]51. Non-coding RNAs, such as miRNAs and lncRNAs, can directly interact and function cooperatively in gene regulation. Non-coding RNAs strongly affect the pathogenesis of OA^[181]52. They have both positive and negative effects on joints; most studies have shown that miRNAs negatively regulate target gene expression by binding to the 3′-untranslated region of target genes, and a variety of miRNAs have been shown to play specific roles in OA progression. MicroRNA-27a attenuates articular cartilage degradation and IL-1β-induced inflammatory responses through the TLR4/NF-κB signaling pathway in articular chondrocytes^[182]53. An analysis of OA synovial tissue revealed six differentially expressed genes identified by LASSO analysis, including the SCRG1 gene, which might be regulated by hsa-miR-363-3p; thus, hsa-miR-363-3p may be a new therapeutic target for OA^[183]54. Moreover, hsa-miR-92a-3p can promote chondrocyte proliferation and differentiation by targeting the WNT5A gene^[184]55. MiRNAs play both positive and negative regulatory roles in OA, and the mechanisms underlying their effects are highly complex and require further in-depth studies. Like miRNAs, lncRNAs play a key role in post-OA regulation, and the lncRNA KCNQ1OT1 can reduce osteoarthritic chondrocyte dysfunction through the miR-218-5p/PIK3C2A axis^[185]56. In human chondrocytes, the overexpression of MALAT1 can inhibit cell viability and promote cartilage ECM degradation by targeting miR145^[186]57. The overexpression of GAS5 increases the levels of chondrocyte catabolism and apoptotic activity^[187]58. Moreover, lncRNAs trigger their function by activating their transcription, transcription elongation and processing, or by the lncRNA molecules themselves^[188]59. In this study, we predicted and characterized 32 miRNAs and 55 lncRNAs with 6 core targets. We constructed an mRNA-miRNA-lncRNA regulatory network and identified a close relationship between them. This regulatory network may represent a potential regulatory mechanism of ATBC in OA and needs to be further investigated. Single-cell analysis has provided novel insights into the heterogeneity of cellular activity in OA. Based on this, we conducted a brief single-cell analysis. Results from the dataset [189]GSE283080 revealed that the core targets are distributed within fibroblasts. Further data analysis indicated that the core targets are expressed across different fibroblast subtypes. Fibroblasts, particularly synovial fibroblasts, play a critical and complex role in the pathological progression of OA. One of their functions is to stimulate synovial hyperplasia and angiogenesis. Synovial thickening is one of the significant pathological features of OA^[190]10holding important diagnostic value for OA onset. Furthermore, pain associated with synovitis is also a key factor linked to OA. This pain can emerge in the early stages when radiographic signs of cartilage damage are not yet significant. Synovial fibroblast hyperplasia is a key characteristic of this process and constitutes the primary source for the increased pro-inflammatory cytokines within the OA joint and the intensified infiltration of activated immune cells. Studies demonstrate that activated fibroblasts are a major source of matrix metalloproteinases (MMPs) and aggrecanases (ADAMTS) in the joint fluid^[191]60. The degradation of type II collagen by MMPs and the cleavage of the core protein of aggrecan by ADAMTS lead to the destruction and loss of the cartilage matrix, representing one of the core mechanisms of cartilage damage in OA. In this study, the molecular docking results confirmed the important role of the core target in the influence of ATBC on OA, as they stably bind to the active site of ATBC (with binding energies lower than − 4.25 kcal/mol). To further validate its stability, we conducted molecular dynamics simulations and systematically calculated the following indicators: RMSD, RMSF, Rg, SASA, MM/GBSA, and FEL, which collectively confirmed the structural stability of the aforementioned complex in a physiological environment. Analyzing from the perspective of molecular mechanisms, this strong stability may enhance the in vivo residence time of the aforementioned complex, thereby continuously regulating the function of the target protein and affecting the progression of OA, providing a key theoretical basis for elucidating the molecular pathological mechanism of ATBC’s influence on OA^[192]7. Although there is no direct evidence conclusively proving a direct effect of ATBC on OA, we hypothesized that ATBC has a deleterious effect on OA through a rigorous screening and prediction process using network toxicology approach combined with an analysis of real-world exposure scenarios. We argue that any factor that can lead to an adverse effect of ATBC on OA needs to be emphasized to minimize unforeseen consequences. Limitations Although the results of this study are novel and encouraging, it has several limitations. First, the findings were primarily based on network toxicology analysis, bioinformatics analysis, and computational molecular docking; direct experimental evidence confirming the relationship between ATBC and increased OA risk is lacking. Second, the research results may be biased due to related issues such as source variability, data quality, and algorithmic bias. Although this study lacked in vivo and in vitro experimental validation, it provides a valuable methodological framework for evaluating the safety of ATBC, offering certain reference values for future research. In the future, we aim to perform standardized cellular and animal experiments to validate the core targets and signaling pathways determined by network analysis and perform long-term, multicenter, high-quality epidemiological studies to monitor and evaluate the evolving relationship between exposure to ATBC and OA. Conclusion To summarize, this research systematically investigated the relationship between ATBC and OA for the first time using network toxicology, bioinformatics, molecular docking, and molecular dynamics simulation methods. These results indicate that ATBC may be closely associated with OA through multiple targets and signaling pathways. This study also promoted the utilization of emerging network toxicology methods related to the assessment of the toxicity of environmental pollutants. Electronic supplementary material Below is the link to the electronic supplementary material. [193]Supplementary Material 1^ (14.3KB, xlsx) [194]Supplementary Material 2^ (840.6KB, docx) Acknowledgements