We used computer simulation analysis to identify AMD-related genes. First, we identified the target organs for amiodarone toxicity using the Protox database (https://tox.charite.de/protox3/) (Banerjee et al., 2018; Banerjee et al., 2024). Then, We obtained its canonical molecular structure and SMILE formula from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (Kim et al., 2024). To identify potential therapeutic targets for AMD, we performed computational target prediction for the active compounds. This process involved querying each compound against three independent databases: ChEMBL (https://www.ebi.ac.uk/chembl/) (Davies et al., 2015; Zdrazil et al., 2024), Similarity Ensemble Approach (SEA) (https://sea.bkslab.org/) (Keiser et al., 2007), and SwissTargetPrediction platform (https://www.swisstargetprediction.ch/) (Daina et al., 2019; Daina and Zoete, 2024). After deduplication, a total of 280 high-confidence AMD-related targets were obtained.

We obtained five PF gene expression datasets (GSE17978, GSE24206, GSE32537, GSE40839, and GSE47460) from the GEO database to explore PF-related genes (Edgar et al., 2002; Barrett et al., 2013). All datasets were retrieved from the GEO database in MINiML format, which contained comprehensive platform, sample, and series metadata. Raw data were normalized using the preprocess Core package. Probe IDs were converted to gene symbols according to corresponding platform annotation files. Batch effect correction was performed using the limma package for samples derived from different batches within the same dataset. When integrating data across multiple datasets or different platforms, only gene symbols common to all datasets were retained. Each dataset or platform was designated as a distinct batch, and batch effects were uniformly adjusted. The effectiveness of batch correction was evaluated by principal component analysis (PCA) conducted before and after normalization. The first four datasets were used as a training cohort. After batch effect correction and standardization, differential expression analysis was performed on this cohort to identify Differentially expressed genes (DEGs). Simultaneously, Weighted Gene Co-Expression Network Analysis (WGCNA) was applied to identify gene modules highly associated with the PF phenotype. The module with the strongest association (MEred) was selected. The union of DEG genes and WGCNA (MEred) genes yielded 233 PF-related genes. All data processing and statistical analysis were performed in the R 4.5.1 environment.

To determine the core targets of AIPF, we intersected 280 AMD-related targets with 233 PF-related genes to derive putative AIPF hub genes. We imported the target gene list into the STRING database (version 11.5) to construct a protein-protein interaction (PPI) network, using Homo sapiens as the species and selecting a medium confidence threshold (>0.4). We exported TSV-formatted data containing nodes, edges, and interaction confidence scores for subsequent network visualization and topological analysis. We subsequently employed the clusterProfiler package in R 4.5.1 to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment studies on these hub genes, considering entries with a corrected p < 0.05 as statistically significant.

To further explore and screen hub AIPF genes, we constructed a ML approach combined with SHAP analysis. 113 ML prediction models were constructed based on transcriptome data, and their performance was validated through independent cohorts using Receiver Operating Characteristic curve (ROC). To determine the optimal predictive model, we conducted a systematic evaluation and comparison based on the model predictive performance and its alignment with the target cohort’s demographic and clinical characteristics. Ultimately, the best algorithm was selected through stratified cross-validation and independent validation on a test set, ensuring the model possesses optimal generalization capability and clinical applicability. The best-performing prediction model was selected for further analysis. SHAP interpretability analysis was applied to quantify the contribution and direction of each gene′s characteristics in the model prediction based on its SHAP value, thereby identifying and ranking hub genes for AIPF classification.

We performed molecular docking simulations to clarify the binding sites between AMD and its target proteins. We obtained the protein structure files of the AIPF hub genes from the PDB database and downloaded the 3D molecular file of AMD from PubChem. Molecular docking was performed using AutoDock Vina, and the data was converted to PDBQT format using MGLTools or OpenBabel. The center coordinates of the binding pocket and the search box size (x, y, z lengths) were defined in the receptor PDBQT, a configuration file (config.txt) was generated, and the exhaustiveness parameter was set. The docking results were evaluated using binding energy (ΔG, kcal/mol); a binding energy below −5.0 kcal/mol indicated strong binding activity. Finally, the binding sites and intermolecular interaction patterns were visualized and analyzed using PyMOL.

Molecular dynamics simulations were performed using the GROMACS 2022, the system was constructed using the AMBER99SB force field and placed in a periodic boundary box, followed by energy minimization to eliminate steric clashes. Subsequently, the system was thoroughly equilibrated successively under NVT and NPT ensembles, employing the V-rescale method for temperature coupling and the Parrinello-Rahman method for pressure coupling, until the temperature, density, and energy of the system stabilized. Following equilibration, a production simulation was conducted for a 50 ns with an integration time step of 2 fs. Short-range interactions were treated using a cutoff distance, while long-range electrostatic interactions were calculated using the Particle Mesh Ewald method. Hydrogen bonds were constrained using the LINCS algorithm. Trajectory data were saved at defined intervals for subsequent structural and dynamical analyses. We employed the built-in “g_sham” script and the “xpm2txt.py” script within GROMACS 2022 to calculate Gibbs free energy based on the RMSD and Rg values of stable complexes, and generated 2D Gibbs free energy contour plots.

We used AMD-induced BEAS-2B cell model for in vitro experiments. The CCK-8 assay was employed to evaluate the proliferation-inhibitory effect of AMD on BEAS-2B cells, thereby determining its optimal intervention concentration and providing a critical dose reference for subsequent mechanism studies. Cells were resuspended in serum-free medium. For the migration assay, cells were seeded at a density of 1 × 10⁶ cells/mL into the upper chamber of Transwell filters (8 μm pore size, Corning). The lower chamber was filled with medium containing 10% FBS to induce migration. For the invasion assay, matrix gel (50 μL/well, Servicebio) was pre-coated onto the lower chamber, and cells were seeded at a density of 5 × 10⁵ cells/mL. After 24–48 h at 37 °C, remove non-migrated cells from the upper chamber. Fix with 4% paraformaldehyde and stain with 0.1% crystal violet. Count cells in the lower chamber across 5 randomly selected fields. Migration/invasion rate = (experimental group cell count/control group cell count) × 100%. Seed cells (5 × 10⁵ cells/well) into 6-well plates. Upon 90% confluence, create linear scratches using sterile pipette tips. Wash debris with PBS and switch to 1% FBS medium to inhibit proliferation. Photograph the fixed sites at 0 h and 24 h. Healing rate = Initial area − 24 h  area / Initial area × 100 % .

BEAS-2B cells were cultured in DMEM +10% FBS and treated with 10 μM AMD for 48 h to establish an epithelial AIPF model. Total RNA was extracted using TRIzol; concentration and purity were assessed by NanoDrop. One microgram of RNA was reverse-transcribed into cDNA. mRNA levels of CTSK, ADORA3, FLT3, and AGER were quantified by SYBR Green qPCR with gene-specific primers. Protein expression was evaluated by Western blot: proteins were extracted in RIPA buffer, quantified by BCA assay, separated by SDS-PAGE, transferred to PVDF membranes, and blocked with 5% non-fat milk. Membranes were incubated overnight at 4 °C with primary antibodies against CTSK, ADORA3, AGER, and FLT3, then with HRP-conjugated secondary antibodies. Subcellular localization and expression were assessed by immunofluorescence: cells on coverslips were fixed, permeabilized, blocked, incubated overnight at 4 °C with the same primary antibodies, then with fluorochrome-conjugated secondaries; nuclei were counterstained with DAPI. Confocal images were acquired and fluorescence intensity quantified in ImageJ. Relative mRNA and protein levels were calculated using the 2^−ΔΔCT method with GAPDH as the internal control.

BEAS-2B cells, STR correctly identified, were purchased from Wuhan Pronosai, catalog number: CL-0016, specification: 1 × 10⁶ cells/T25 culture flask, growth medium: Ham’s F-12K + 10%FBS+1%P/S. The normal human lung epithelial cells BEAS-2B used in this study were purchased from Shanghai Zhongqiao Xinzhou Biotechnology Co., Ltd., catalog number ZQ0381. Anti-Cathepsin K Rabbit pAb was purchased from Servicebio (product number: GB111276). Anti-Adenosine A3 Receptor/A3AR Rabbit pAb was purchased from Servicebio, accession number: GB112420; Anti-RAGE Rabbit pAb was purchased from Servicebio, product number: GB11278. Anti-FLT3 Rabbit pAb was purchased from Servicebio (product number: GB111358). Anti-GAPDH Rabbit pAb was purchased from Servicebio, accession number: GB11002. All primers sequences are listed as follows: GAPDH -F: 5′-CCT​GCA​CCA​CCA​ACT​GCT​TA-3′; GAPDH-R: 5′-GGC​CAT​CCA​CAG​TCT​TCT​GAG-3′; CTSK-F: 5′-CAG​CAG​ACG​CAG​ACA​TTA​CG-3′, CTSK-R: 5′-GAT​GCC​AGG​GTA​CAG​CTC​TC-3′; ADORA3-F: 5′-TGG​TCG​TCG​TCT​TCT​TCA​TC-3′, ADORA3-R: 5′-CAA​AGC​CAC​AAG​GAA​GAC​CT-3′; FLT3-F: 5′-TGC​AAG​CAA​GGA​CGT​GAA​G-3′, FLT3-R: 5′-GAC​TTG​CCA​GGA​AGG​TGA​TG-3′; AGER-F: 5′-GGA​AGG​AAC​GTC​TCA​ACC​AA-3′, AGER-R: 5′-GGT​CTG​GGC​ATT​TAT​CCT​CA-3′.

This study first conducted computer simulation analysis to determine the target sites of AMD in the human body. The definitive molecular structure of AMD (SMILES: CCCCC1 = C(C2 = CC = CC = C2O1)C(=O)C3 = CC(=C(C(=C3)I)OCCN(CC)CC)I) was sourced from the PubChem, as shown in Figures 1A,B. By consulting the Protox database, we found that the predicted LD50 of AMD is 3,000 mg/kg, and the predicted Toxicity Class is 5/6 (1/6 is the most toxic). Table 1 presents the Toxicity Model Report, demonstrating that ADM exhibits the most significant pulmonary toxicity, mostly through immunotoxicity. Subsequently, parallel analysis was conducted utilizing three target prediction platforms: the ChEMBL, SEA, and the SwissTargetPrediction platform.

The prediction results from these three platforms are summarized in Figure 1C. The ChEMBL database identified 189 distinct targets, the SwissTargetPrediction identified 101, and the SEA technique identified 43. Remarkably, only 43 targets were identified in common by those three platforms, illustrating the complementarity of the various computational approaches. By consolidating all unique targets, a total of 280 putative high-confidence targets of AMD were identified, establishing the foundational dataset for further analysis.

This study obtained five separate gene expression datasets from the GEO database (GSE17978, GSE24206, GSE32537, GSE40839, and GSE47460) to uncover hub genes linked to PF pathogenesis. The initial four datasets (GSE17978, GSE24206, GSE32537, and GSE40839) were amalgamated as training sets for preliminary gene screening, whereas GSE47460 was preserved as a subsequent validation set (Table 2). Batch correction was performed to mitigate the confounding effects of technical variation arising from different experimental runs, thereby ensuring that the downstream analyses more accurately reflect underlying biological differences. First, we normalized and corrected for batch effects in the training set to eliminate abiotic variation caused by differences in experimental platforms and techniques (Figures 2A,B). Subsequently, differential expression analysis was performed on the corrected data using R 4.5.1. A total of 109 DEGs were identified, including both upregulated and downregulated genes. The results of the differential analysis were visualized using a volcano plot (Figure 2C). A heat map (Figure 2D) was also created, showing the expression patterns of the 50 most significantly DEGs across each sample. The volcano plot and heat map visually demonstrate the differences in gene expression profiles between the Treat and Control groups. We can observe that in the Treat group, genes such as MGP, COL14A1, and CFH are highly expressed, while genes including NFKBIA, MAP3K8, and THBS1 are downregulated.

Then, we utilized WGCNA to find the gene clusters most closely connected with the PF phenotype. Through the construction of a scale-free co-expression network, we identified several gene modules, each denoted by a distinct color. The examination of module-trait relationships indicated that the red module (MEred) exhibited the strongest association with PF disease state (correlation coefficient = 0.75, p = 2 × 10⁻⁵⁸), as illustrated in Figures 3A–C. Scatter plot depicting a strong positive correlation (cor = 0.86, p = 5.7 × 10⁻⁴⁶) between module membership and gene significance for the treatment group, Each point represents a gene within the red module (Figure 3D). Consequently, the 153 genes within this module were designated as PF-related module genes and utilized for further investigation. We integrated the genes identified by the two methods described above and took the union of DEGs and WGCNA (MEred) genes. As shown in Figure 3E, we obtained 233 PF-related genes.

To clarify the pathogenic mechanisms of AIPF, we compiled gene sets associated with AMD exposure and PF from public databases and performed enrichment analysis. This analysis revealed that 280 genes were significantly associated with AMD, while 233 genes were significantly associated with PF. Notably, eight intersecting genes (MMP1, ADORA3, ABCB1, AGER, SLC6A4, CTSK, SPHK1, and FLT3) were shared by both AMD and PF, suggesting that AMD may induce PF through these hub genes. In the constructed PPI network, several key hub nodes with high connectivity were identified, including CTSK, AGER, FLT3, and MMP9. These molecules occupy central positions in the network topology, suggesting their potential involvement in critical biological processes. To identify key mechanisms underlying AIPF, we performed GO and KEGG analyses. As shown in Figure 4, several biological processes closely related to inflammatory responses, ECM remodeling, collagen metabolism, and the MAPK signaling pathway were significantly enriched.

Functional enrichment analysis of DGEs identified critical processes in AIPF. GO analysis highlighted significant enrichment in critical biological processes, including “positive regulation of p38MAPK cascade” “collagen catabolic process” “microglial cell activation” and “leukocyte activation involved in inflammatory response”. KEGG pathway analysis further established the central role of the “p38MAPK cascade” alongside “ECM disassembly” “positive regulation of NF-κB signaling” and “IL-1 production”. The coordinated enrichment of these pathways indicates a mechanistic framework where AMD triggers a pro-fibrotic response through p38 MAPK/NF-κB-driven inflammation and ECM remodeling.

We employed multiple ML algorithms to construct AIPF diagnostic and predictive models, evaluating their performance on independent cohorts. As shown in Figure 5A, all models were trained on our primary dataset and validated on cohorts GSE47460A and GSE47460B. Results show that the glmBoost + GBM model demonstrated outstanding and stable predictive performance, achieving the highest area under the curve (AUC) on the training set (AUC = 0.987) and maintaining robust accuracy in external validation (AUCs of 0.826 and 0.721, respectively). Through a two-stage integrated design, sparse feature selection based on generalized linear models boosting (glmBoost) is combined with robust nonlinear gradient boosting modeling (GBM), structurally achieving complementary dimension reduction and complex relationship modeling. This design effectively mitigates overfitting risks in high-dimensional data while leveraging GBM’s capability to capture complex interactions among key features, thereby enhancing the model’s generalization performance and predictive accuracy. The model’s rigor is doubly validated: through its robust performance maintained on an independent validation cohort, and through subsequent SHAP-based interpretability analysis that thoroughly elucidates its decision logic. Four hub genes were involved in glmBoost + GBM model construction: ADORA3, CTSK, AGER, and FLT3, The SHAP analysis scatter plot is shown in Figure 5B. SHAP interpretable analysis revealed distinct functional contributions: ADORA3 (SHAP value = 0.135), CTSK (SHAP value = 0.116), AGER (SHAP value = 0.091), FLT3 (SHAP value = 0.074) emerged as the most influential predictors. Among these, CTSK exhibited positive regulatory effects while ADORA3, AGER, FLT3 exhibited negative regulatory effects, in a representative sample, features were ranked according to their contribution (SHAP value). In a representative experiment, the analysis showed that ADORA3 (SHAP = −0.278) and CTSK (SHAP = −0.198) were the strongest positive drivers, significantly pushing the model’s predicted value from the baseline (E[f(x)] = 0.473) to a higher risk level (Figure 5C). Figures 5D,E shows the contribution (Figure 5D) and direction (Figure 5E) of the SHAP analysis. To further assess the viability of these core genes as diagnostic biomarkers, we constructed ROC curves. Figure 5F illustrates that all four hub genes exhibited strong diagnostic discrimination capability, with the area under the curve above 0.70. CTSK and ADORA3 exhibit the highest diagnostic efficacy, with AUC values of 0.819 and 0.814, respectively, followed by FLT3, and AGER indicating significant promise for clinical application.

Based on the SHAP analysis, the four most prominent potential target genes (CTSK, ADORA3, AGER, and FLT3) were prioritized for further investigation. To elucidate the molecular interactions between amiodarone and these targets, systematic molecular docking was conducted, followed by molecular dynamics simulations. The docking results demonstrated that AMD exhibited strong spontaneous binding affinity to all tested proteins (Vina Score ≤ −6.0 kcal/mol), with the strongest binding observed for AGER (−7.6 kcal/mol) and ADORA3 (−7.5 kcal/mol). Detailed interaction analysis revealed that amiodarone stabilizes within the active pockets of each protein through multiple non-covalent interactions, including hydrogen bonds, extensive hydrophobic contacts, π-π stacking, and van der Waals forces, providing preliminary insight into its potential multi-target mechanism. The results of molecular docking are shown in Table 3 and Figure 6.

To further investigate the dynamic binding characteristics of the complexes, we performed 50 ns molecular dynamics simulations and binding free energy calculations. Analyses of RMSD (Brüschweiler, 2003; He et al., 2023), RMSF (Baskin, 1968), Rg (Davoudmanesh and Mosaabadi, 2018), and SASA (Richmond, 1984; Castrosanto et al., 2023) consistently indicated that the amiodarone-protein complexes maintained good structural stability throughout the simulation, with the CTSK-amiodarone complex displaying the highest conformational rigidity. Free energy landscape (FEL) analysis thermodynamically confirmed that the CTSK system possessed the most concentrated low-energy conformational basin. More importantly, MM/PBSA calculations quantitatively verified that the binding of amiodarone to ADORA3, CTSK, and FLT3 was highly thermodynamically favorable (ΔG_bind = −47.61, −39.45, and −36.01 kcal/mol, respectively), with Van Der Waals interactions identified as the major energetic contributor. Together, these results provide computational evidence at both the atomic dynamical level and the thermodynamic quantitative level, supporting the stable and preferential binding of AMD to the aforementioned target proteins and offering a solid foundation for its potential multi-target pharmacological effects. Analysis of Gibbs FEL maps indicates that all four protein-ligand complex systems exhibit highly stable binding conformations. The free energy distributions for each system reveal a single, continuous, and narrowly concentrated minimum energy region. The deep blue energy traps are broad and smooth, suggesting minimal conformational fluctuations and high structural rigidity in the bound complexes, with no apparent metastable states or conformational transition pathways. Therefore, AMD forms strong and stable specific interactions with CTSK, ADORA3, AGER, and FLT3, providing crucial evidence of conformational stability for subsequent drug design or mechanism studies. The results of molecular dynamics simulations are shown in Table 4 and Figure 7.

We used a CCK-8 assay to assess AMD’s impact on BEAS-2B cell viability and find the best intervention concentration. In Figure 8A, cell viability was not significantly affected by AMD concentrations from 1 μM to 10 μM. Cell viability peaked at 10 μM AMD concentration. In the 10–25 μM range, cell viability significantly decreased. In following trials, 10 μM AMD was found to be the most efficient intervention concentration for modulating BEAS-2B cells without causing severe cytotoxicity. The wound healing assay demonstrated a similar pro-migratory effect induced by AMD. While wound closure rates were comparable at the initial stages, the AMD-treated group displayed a significantly accelerated wound closure over 24 h. The relative wound healing area was substantially smaller in the AMD group compared to the control, further confirming that AMD stimulation enhances the lateral migration of BEAS-2B cells in a two-dimensional model (Figure 8B). To functionally assess the impact of AMD on the motility and invasive potential of BEAS-2B cells, we performed transwell migration and invasion assays alongside a wound healing assay. At the 0-h time point, no discernible difference in cell migration or invasion was observed between the control and AMD-treated groups. However, after a 24-h incubation period, AMD-treated cells exhibited a significantly enhanced migratory and invasive capacity compared to the control. The significant increase in the number of cells migrating through Matrigel-coated membranes suggests that AMD exposure enhances cell motility and confers a pronounced pro-invasive phenotype (Figure 8C). These functional assays demonstrate that AMD intervention significantly enhances the migratory, invasive, and wound healing capabilities of pulmonary epithelial cells, thereby revealing the toxicological mechanisms underlying AIPF.

Furthermore, we analyzed the mRNA expression levels of four hub genes (CTSK, ADORA3, AGER, and FLT3) in AMD-induced BEAS-2B cells using qRT-PCR technology. Quantitative analysis indicated that, in comparison to the control group, CTSK mRNA expression shown significant upregulation (p < 0.05), whereas both ADORA3 and AGER displayed substantial downregulation (p < 0.05, p < 0.01). Nevertheless, under identical treatment settings, FLT3 expression levels exhibited no statistically significant alterations. The results indicate that the overexpression of CTSK and the downregulation of ADORA3 and AGER may be pivotal in the pulmonary epithelial cell response to AMD and the probable advancement of PF (Figure 9A). Western blot analysis revealed differential expression profiles of target proteins between control group and AMD group.​ Compared to the control group, the expression of CTSK was markedly upregulated in the AMD group (p < 0.05). In contrast, the level of ADORA3 and AGER showed a distinct downregulation (p < 0.05, p < 0.05). The expression of FLT3 appeared comparable between the two groups, indicating no substantial change (Figure 9B). Immunofluorescence analysis further validated the protein-level alterations of these hub genes and delineated their distinct subcellular localizations in BEAS-2B cells. Consistent with the mRNA data, CTSK protein expression was significantly upregulated and exhibited a characteristic granular pattern co-localized with lysosomal markers within the cytoplasm. Conversely, ADORA3 displayed a marked downregulation and was predominantly localized to the plasma membrane, exhibiting a perimenbranous ring-like staining. AGER protein expression was also substantially reduced, showing a combined localization pattern with signal present both at the plasma membrane and diffusely throughout the cytoplasm. In contrast, FLT3 protein levels showed no significant change following AMD treatment, and it maintained its expected specific localization to the cell membrane (Figure 9C). These results collectively confirm the differential expression of CTSK, ADORA3, and AGER at the protein level and define their precise cellular compartments, providing crucial morphological evidence for their potential roles in AMD-induced epithelial cell injury and the subsequent pathogenesis of pulmonary fibrosis.

In summary, these data demonstrate that AMD exposure reprograms pulmonary epithelial cells through upregulation of CTSK and downregulation of ADORA3 and AGER, leading to an enhanced migratory and invasive phenotype, highlighting a pivotal role for CTSK-driven epithelial dysfunction in the progression of PF.