To gather information on the active ingredients of BKN, we conducted a search using keywords such as “Qingdai,” “Pingbeimu,” and “Baiguoren.” This search was performed on the TCM system pharmacology technology platform (TCMSP, http://tcmspw.com/tcmsp.php) and the China national knowledge infrastructure (CNKI) reviews (Ru et al., 2014). To filter out irrelevant components, we applied specific criteria based on oral bioavailability (OB) ≥ 30% and drug-likeness (DL) ≥ 0.18. Once the active components were identified through screening in the TCMSP database, their SDF structure formats were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). We used PharmMapper database (http://www.lilab-ecust.cn/pharmmapper/) (Liu et al., 2010) to predict the potential targets of BKN through SDF structure. Further analysis was carried out by selecting BKN component targets with a Norm Fit value of ≥0.5. For standardization purposes, we selected species “Homo sapiens” in uniport (https://www.uniprot.org/) (Uszkoreit et al., 2021) to convert the Uniprot ID into the gene symbol.

We obtained microarray data for influenza from two sources, GSE34205 and GSE42026, which were downloaded from the GEO database. The microarray platforms used were GPL6947 (Illumina HumanHT-12 V3.0 expression beadchip) and GPL570 (Affymetrix Human Gene Expression Array). These datasets included samples from pediatric influenza patients as well as healthy individuals without any other diseases. Both groups had a sample size of more than 10. We obtained gene expression and clinical data, and performed gene symbol annotation and data correction using Perl code. GSE34205 (Ioannidis et al., 2012) was used as the analysis set. This study excluded children suspected or confirmed to have multiple bacterial infections, potential chronic diseases (i.e., congenital heart disease or renal insufficiency), immunodeficiency, or receiving systemic steroids or other immunomodulatory therapies. Among them, all influenza virus infected individuals were confirmed and diagnosed through microbiology, and all blood samples were collected within 42–72 h after hospitalization. A total of 28 influenza virus patients were included in the study, with a median time of 5.5 (1.4–21) months. Blood samples from 12 healthy children with a median age of 18.5 (10.5–26) months who underwent elective surgery or outpatient visits were used as control samples (Supplementary Table S1). And GSE42026 (Herberg et al., 2013) was used as the validation set, which includes 19 H1N1/09 samples and 33 healthy controls. The LIMMA software package (Ritchie et al., 2015) was used to identify differentially expressed genes (DEGs) between the pediatric influenza group and the control group as disease targets in network pharmacology. A heatmap was generated based on these DEGs using R software’s pheatmap package.

To establish a scientifically sound relationship between compounds in herbs and their targets, we utilized Cytoscape 3.9.1 to construct an herb–component–target network that visually represented this relationship. We conducted a topological analysis to identify key components within this network. Immune related genes were collected through the ImmPort database (https://www.immport.org). The overlapping genes between immune related genes and BKN-DEGs are defined as BKN-ImmuneDEGs. Selection and visualization of BKN-ImmuneDEGs by using R packages “ggplot2” and “pheatmap”. Analyzing the gene network of BKN-ImmuneDEGs through the GeneMANIA database (https://genemania.org/).

To explore the core targets of BKN’s anti influenza effect, a protein-protein interaction (PPI) network was built using the STRING (https://cn.string-db.org/; Version 12.0) (Feng et al., 2019) database to examine the functional interaction between proteins, with setting the species to H. sapiens and network confidence score ≥0.4. Furthermore, the Cytohubba plug-ins were utilized to identify the core targets of BKN for combating pediatric influenza.

For the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis, we imported the gene symbols into the Metascape platform (http://metascape. org/), the organism used was H. sapiens (Zhou et al., 2019). The GO and KEGG pathway enrichment analysis findings were saved, and bioinformatics software (http://www.bioinformatics.com.cn/) and R version 4.3.1 finished the visualization.

Using R packages (such as “limma,” “pheatmap,” and “ggpubr”), we retrieved the expression levels of BKN-related genes from the influenza group and the normal group and performed differential expression analysis. Box plots and a heat map were utilized to present the results. BKN-related core genes (BKN-RCG) were classified as genes with a p-value less than 0.05. The BKN core genes were found on the chromosomes using perl code, and we then utilized the R package “Rcircos” to visualize them as circle plots. Additionally, using the “cor” command, correlation coefficients for each BKN-RCG were computed and displayed.

We used the R package “ConsensusClusterPlus” to cluster influenza samples according to the expression of BKN-RCG with a k-means clustering method, Euclidean distance type, and maximum of nine clusters (Wilkerson and Hayes, 2010). Heat maps were used to compare the expression levels of the generated clusters. Lastly, PCA was used to illustrate the differences between the subtypes.

To investigate the expression of gene sets, we utilized the weighted gene co-expression network analysis (WGCNA) approach. The R package “WGCNA” was employed to construct a co-expression network for DEGs. We selected the top 25% of genes with the highest variation to ensure biologically relevant gene modules and determined the soft-thresholding power based on the standard scale-free topology criterion (scale-free R ² ≥ 0.9), which is a widely accepted approach (Langfelder and Horvath, 2008). By determining the optimal soft power, we generated a weighted adjacency matrix and transformed it into a topological overlap matrix (TOM). TOM is a matrix used to describe the degree of topological overlap between nodes in complex networks, which is widely used to analyze the characteristics such as the stability, robustness and community structure of the network. Larger TOM values indicate stronger stability and robustness of the network, while smaller TOM values indicate weaker stability and robustness of the network (Shuai et al., 2021). Using a hierarchical clustering tree technique with a minimum module size set at 100, we identified modules based on TOM difference measure (1-TOM). From these co-expression modules, we could extract key modules that exhibited high correlation with specific phenotypes or diseases. Hub genes were then identified based on internal connectivity within the key module and their correlation with feature vectors of the key module (Langfelder and Horvath, 2008).

We applied the “VennDiagram” program to identify intersected genes between DEGs, diseases, and clusters. Subsequently, we used three machine-learning algorithms to screen important biomarkers for pediatric influenza: RF (Ishwaran and Kogalur, 2010; Wang H. et al., 2016; Perera et al., 2020), LASSO logistic regression (Cheung-Lee and Link, 2019a; Fernández-Delgado et al., 2019), and SVM-RF (Huang et al., 2018a). The SVM algorithm identifies the optimal variables by deleting SVM-generated eigenvectors (Liu et al., 2023). LASSO utilizes regularization to reduce prediction errors (Yang et al., 2018). RF can handle high-dimensional data, establish predictive models, and predict the importance of each variable (Blanchet et al., 2020). Finally, the genes obtained by all three algorithms were ultimately crossed to screen for the characteristic genes. To validate the usefulness of these hub genes comprehensively, the dataset from GSE42026 served as our validation set. The predictive capability of these algorithms was evaluated using ROC curves, and area under curve (AUC)was calculated accordingly.

We obtained the 3D structures of BKN-core components from Pubchem. And the 3D structure of feature genes was obtained through the AlphaFold Protein Structure Database (https://alphafold.com/). Initial conformations for each component were established using energy minimization techniques in Chem3D 14.0 software. AutoDockTools was utilized for molecular docking, with the best four combinations of docking selected and visualized using Pymol and Discovery Studio 2019.

MD simulation was performed on the receptor-ligand with the lowest binding energy from molecular docking results. GROMACS 2023.3 software was used to carry out MD simulations, utilizing amber14sb force field and general Amber force field (GAFF) to generate parameter and topology data for proteins and ligands respectively. The simulation box size was optimized based on protein atom distance greater than 1 nm, followed by water molecule filling at a density of one to achieve electrical neutrality through Cl⁻ and Na⁺ ion replacement methods. Using the steepest descent approach, 5.0 × 10⁴ steps of energy optimization were carried out to minimize the overall system’s energy consumption and, finally to lower the overall system’s unreasonable contact or atom overlap. Then, the system’s temperature was stabilized through first-phase equilibration using the NVT ensemble at 300 K for 100 ps. The NPT ensemble was used to mimic second-phase equilibration at 100 ps and 1 bar. To thoroughly pre-equilibrate the simulation system, the main goal of the simulation is to optimize the interaction between the target protein and the solvent and ions. MD simulation ran for 50 ns at a temperature of 310.15 K and 1 atm of pressure in an isothermal and isostatic ensemble. The V-Rescale and C-Rescale techniques were used to regulate the temperature and pressure, respectively. The temperature and pressure coupling constants were 0.1 ps. The Van der Waals force was computed using the Leonard-Jones function, and the nonbond truncation distance was established at 1.4 nm. The LINCS algorithm limited the bond length of every atom. Using a 0.16 nm Fourier spacing and the Particle Mesh-Ewald method, the long-range electrostatic interaction was computed. Finally, we used gmx_mmpbsa to calculate the binding free energy of the compound (http://jerkwin.github.io/gmxtool). Notably, to further demonstrate the stability of the binding between the small molecule and the protein, all MD simulations were performed twice. In the data presentation, the RMSD, RMSF, Rg, H-bond, and SASA metrics represent the mean values from the two simulations, while the MMPBSA analysis was conducted on the segments of the trajectories that exhibited greater stability.

The BKN active components were screened using the TCMSP database. Twenty prospective BKN compounds, eight from Pingbeimu, eight from Qingdai, and four from Baiguoren, were ultimately gathered (Supplementary Table S2). PharmMapper and UniProt databases were used to predict potential targets. Finally, 200 targets in Pingbeimu, 250 targets in Qingdai, and 229 targets in Baiguoren were collected.

By computing the difference in gene expression between 22 normal samples and 28 influenza infection samples, we were able to select DEGs of pediatric influenza. Ultimately, the limma package in R software was used to identify 3,819 DEGs from the GSE34205 dataset. Figure 2A shows that 1747 upregulated and 2072 downregulated genes were found. Figure 2B displays the top 10 genes with the greatest up- and downregulation.

58 potential targets of BKN in combating pediatric influenza were shown in Supplementary Table S3. We utilized Cytoscape and its plugins to construct “herb-component-target” network, with 79 nodes and 497 edges (Figure 3A). β-sitosterol was the most prolific, with 42 potential targets, secondly isovitexin (38), and pingbeimine C (34), bisindigotin, indican, beta-sitosterol, sitosterol, coniferin, Peimisine, and sevcoridinine as illustrated in Figure 3B. The course of influenza disease is closely related to immune activity. Therefore, we further investigated the correlation between immune related genes and selected genes that overlap with BKN related DEGs as BKN-ImmuneDEGs. A total of 17 BKN ImmuneDEGs were obtained, including nine upregulated genes (CTSG, CHIT1, RNASE3, PPARG, PLAU, TYMP, IGF1, NR1I3, and RARA) and eight downregulated genes (PIK3CG, SYK, FGFR1, IGF1R, NR3C2, PPARA, TGFBR1, and LCK), as illustrated in Figures 3C,D. Further import these BKN-ImmuneDEGs into GeneMANIA to construct a gene network (Figure 3E). The gene mainly involved in response to phosphatidylinositol 3-kinase signaling, regulation of MAP kinase activity, regulation of interleukin-4 production, positive regulation of T cell activation, negative regulation of MAPK cascade, and intracellular receptor signaling pathway.

Using the intersection targets from Cytoscape and STRING, we built a PPI network with 144 edges and 47 nodes, with an average node degree of 4.97 (Figure 4A). A p-value of less than 1.0e-16 indicated significant clustering in the PPI network. Twenty-nine key genes were initially filtered out of the data based on degree values larger than or equal to the standard degree value (Figure 4B). Then, we used eight random algorithms from the CytoHubba plugin to further determine the top 10 central genes, and identified the core targets through intersection, as illustrated in Figure 4C. We use R to draw an Upset graph to display the results of each algorithm. The core targets include PPARG, MMP2, GSK3B, PARP1, CCNA2, and IGF1, as illustrated in Figure 4D.

Based on biological process (BP), cellular component (CC), and molecular function (MF), GO enrichment analysis was examined. To further investigate the various mechanisms by which BKN contributes to pediatric influenza treatment, we ran a GO enrichment analysis of the 58 common targets. Three groups containing a total of 521 items were found: 445 BP, 48 MF, and 28 CC (Supplementary Table S4). As bubble charts and lollipop charts, the top 10 enriched BP terms, MF terms, and CC terms are displayed (Figures 5A,B). Result showed that the BP terms were mainly involved the cellular response to lipid, response to wounding, response to alcohol, response to hormone, and regulation of leukocyte cell-cell adhesion. Highly enriched MF terms were carboxylic acid binding, nuclear receptor activity, protein kinase binding, endopeptidase activity, and oxidoreductase activity. In addition, the CC terms included secretory granule lumen, specific granule, ficolin-1-rich granule lumen, and protein kinase complex. We evaluated KEGG pathway enrichment analysis of 58 common targets (with P < 0.05 as the significance level). 52 pathways had targets that were substantially enriched (Supplementary Table S5). Gene counts were used to filter the top 17 significantly enriched pathways (Figures 5C,D). KEGG enrichment results showed that nucleotide metabolism, AMPK signaling pathway, PI3K-Akt signaling pathway, ovarian steroidogenesis, and pathways in cancer.

We obtained 29 BKN core genes through network pharmacology analysis (Supplementary Table S6). To further analyze the network pharmacology results, we obtained 29 BKN-RCG by intersecting the BKN target with DEGs. Among these genes, a total of 18 genes are highly expressed in pediatric influenza, and a total of 11 genes are highly expressed in the normal group (Figures 6A,B). Figure 6C shows the specific chromosomal locations of the BKN-RCG. As shown in Figures 6D,E, correlation analysis between every two BKN-RCG in pediatric influenza samples revealed a substantial correlation between the BKN-RCG, with the correlation being predominantly positive.

To identify BKN-related molecular clusters in influenza. We used all 58 intersection genes related to influenza in BKN. Studies have indicated that the most stable and credible clustering findings (Figure 7A) and cumulative distribution function (CDF) curves fluctuating within a minimum range (Figure 7B) were obtained when k = 2, and the consistency score was greater than 0.9 (Figure 7C). A significant difference between the two clusters was found using principal component analysis (PCA) (Figure 7D). The expression differences of BKN-related intersection genes between Cluster 1, Cluster 2 as illustrated in Figure 7E.

We constructed the co-expression network through WGCNA package for identifying significant module genes of pediatric influenza. Grouping genes with comparable expression patterns into the same gene module through average linkage clustering. We screened the top 25% of genes with the highest fluctuations for WGCNA analysis. The results showed that when the soft threshold was 3 and the scale-free R2 was 0.9, we identified co-expressed gene modules (Figure 8A). A total of four distinct gene co-expression modules were generated via WGCNA analysis while showing the heatmap of the TOM (Figures 8B–D). The module with the smallest p-value is the one most related to the disease. The WGCNA result shows that the blue module was highly related to pediatric influenza (Figure 8E). Meanwhile, the blue module was positively correlated with module-associated genes (Figure 8F).

Additionally, the WGCNA algorithm was used to generate gene modules and co-expression networks for three BKN-related molecular clusters to identify the critical gene modules linked to pediatric influenza. The results revealed the co-expression gene modules when the scale-free R2 was 0.9 and the soft threshold was 3 (Figure 9A). While displaying the TOM heatmap, the dynamic cutting algorithm produced three co-expression modules in various hues (Figures 9B–D). When the genes within the three modules were compared to ascertain the degree of gene co-expression similarity between influenza and controls, the 5119 genes in the turquoise module were most strongly linked to the BKN-related molecular cluster (Figure 9E). Moreover, the turquoise module showed a high correlation with genes related to modules (Figure 9F).

Taking the intersection of DEGs, WGCNA, and BKN-related cluster WGCNA, we identified 168 key genes (Figure 10A) and performed further exploration of them. Lasso regression, RF, and SVM-RFE algorithms were adopted to screen candidate hub genes for ROC evaluation. According to the SVM-RFE algorithms display N = 2, Root Mean Squared Error (RMSE) (0.24) is the smallest (Figure 10B). This indicated that the model performance was optimal under this combination of features. Therefore, the two optimal candidate genes were screened out from SVM-RFE. The relationship between the error rate, the number of classification trees, and the two genes in descending order of relative relevance was found using RF in conjunction with feature selection (Figures 10C,D). Interferon alpha-inducible protein 27 (IFI27) and Otoferlin (OTOF) rank first and second respectively, suggesting their significance in classification tasks. The optimal λ value is selected by using 10-fold cross-validation to achieve feature selection. And 16 predicted genes were chosen using LASSO regression analysis from the set of statistically significant univariate variables. Figure 10E illustrates the Lasso regression results. The three algorithms identified IFI27 and OTOF as hub genes with overlap (Figure 10F).

We constructed a column chart using IFI27 and OTOF to obtain individual score tables (Figure 11A). The total of these distinctive genes’ expression scores can be used to calculate treatment sensitivity and forecast the risk rate of these genes in the development of pediatric influenza. Two indicators of the accuracy of this prediction are the proximity of the dashed and solid lines in the calibration curve (Figure 11B) and the distance between the red and gray lines in the decision curve (Figure 11C). The ROC curves for IFI27 and OTOF (Figures 12A,B) demonstrated their likelihood as important biomarkers with AUCs of 0.965 and 0.935, respectively. This suggests that the biological markers had a good predictive value accuracy. Within the GSE42026 validation set, there was a significant difference (P < 0.01) in the expression of IFI27 and OTOF between the pediatric sepsis and control groups (Figures 12C,D). The ROC curves for IFI27 and OTOF, with AUCs of 0.971 and 0.912, respectively, suggested that their likelihood as valuable biomarkers in the GSE42026 validation set (Figures 12E,F).

We conducted molecular docking to investigate whether the proteins encoded by the pediatric influenza hub genes we obtained can bind and function with the active components of BKN. The protein structures of IFI27 (AF-P40305-F1-v4) and OTOF (AF-Q9HC10-F1-v4) were downloaded from the AlphaFold protein structure database. The results indicate that the hub genes and active components could form stable structures, and the binding energy of most docking combinations is lower than −5.0 kcal/mol (Figure 13A). This implies that most of the hub genes and active components could form stable structures. The best docking combinations (Peimisine- IFI27 and Peimisine-OTOF) are shown in Figures 13B,C.

In view of the results of molecular docking, the following MD studies mainly focus on the relationship between Peimisine-IFI27 and Peimisine-OTOF, respectively. MD reveals the stability of protein-ligand complexes in physiological aqueous solutions at 37°C (310.15 K) and in silico under standard atmospheric pressure (a pressure of 1 atm). To assess their performance, we used Peimisine-OTOF and Peimisine-IFI27 for 50 ns dynamics simulations. We ran the molecular dynamics simulation twice to ensure reproducibility. And the standard deviation was used to represent the accuracy of repetitive measurements. The results showed that the value of the standard deviation can prove that the results in this study had excellent reproducibility. The detailed results of the two simulations can be found in Supplementary Tables S7 and S8. The initial and final structures of every trajectory in PDB format can be found in Supplementary Material.

It can be observed that all two complexes showed good stability, the average RMSD values for Peimisine- IFI27 and Peimisine-OTOF are as follows: 0.737 nm and 0.729 nm. After 25 ns of MD simulation, the RMSD-time curves became steady (Figure 14A), and the RMSD value of Peimisine-IFI27 was stable at 0.8 nm, while the Peimisine-OTOF was stable at 0.7 nm. In general, the system is said to be stable once the RMSD value stabilizes, and the fluctuation range is smaller than 0.2 nm (Saranya et al., 2023). The majority of the atomic RMSF values varied between 0.2 and 1.2 nm during the MD simulation, the average RMSF values for Peimisine-IFI27 and Peimisine-OTOF are as follows: 0.304 nm and 0.386 nm. Additionally, a nearly horizontal line was created by the Rg values, indicating a relatively stable complex system, the average Rg values for Peimisine-IFI27 and Peimisine-OTOF are as follows: 1.512 nm and 4.555 nm. (Figures 14B,C). Moreover, the Rg values of Peimisine- IFI27 were smaller than Peimisine-OTOF. In conclusion, these findings demonstrate the superior conformational stability of those two complexes and provide compelling evidence for the validity of docking. In addition, to further evaluate the binding interactions in protein-ligand complexes, we quantified the solvent-accessible surface area (SASA) to indicate changes in molecular conformation or exposure during the simulation process. The average SASA values for complexes Peimisine- IFI27 and Peimisine-OTOF during MD simulation are as follows: 73.792 nm² and 1018.103 nm². The SASA values exhibited minor fluctuations throughout the course of the simulation, demonstrating their stable interaction (Figure 14D).

To determine the strength and stability of protein-ligand interactions, we computed the difference between the free energies of bound and unbound states utilizing Molecular Mechanics–Poisson Boltzmann Surface Area (MMPBSA). We selected 10 ns trajectories to choose segments with steady RMSD values for MMPBSA energy analysis. The sum of the entropy (–TΔS) and enthalpy (ΔH) changes was the total binding energy (Wilson Alphonse and Kannan, 2023). We utilized the gmx_mmpbsa.bsh script to evaluate the entropy change in each complex at intervals of 1000 ps. The total binding energy was dissected into four distinct components as shown in Table 1. Table 1 presents the outcomes of protein-ligand binding energies. In the system comprising protein-ligand complexes, Peimisine-OTOF and Peimisine-IFI27 proteins exhibited negative binding free energies with small molecules (−126.691 kJ/mol and −125.930 kJ/mol respectively), indicating a more stable association between Peimisine and IFI27. Notably, Van der Waals interaction played a significant role in this main energy interaction.

Finally, we determined the change in the alteration in hydrogen bond count within the two complexes throughout the 50 ns simulation process. Hydrogen bonds served as indicators of the affinity between the ligand and protein. The docked complex exhibited a stable pattern of hydrogen bonds (Figure 15) (Cong et al., 2023). Our findings indicated that there is fluctuation in the number of hydrogen bonds within both complexes, ranging from 0 to 2. The average count of hydrogen bonds observed during MD simulation for Peimisine-IFI27 and Peimisine-OTOF are as follows: 0.041 and 0.346 respectively.

To demonstrate the stability of the binding between the small molecule Peimisine and two proteins during the MD process, the python package Prolif and MDAnalysis was employed to assess the bonding interactions between the small molecules and the proteins throughout the MD simulations (Figures 16A,B). The results indicated that in the IFI27 protein, the residues Val14, Val37, Ala43, Met44, Leu102, Thr103, Ile106, and Ile110 maintained hydrophobic and van der Waals interactions with the small molecule for most of the simulation time. In contrast, in the OTOF protein, the residues Trp1325, Tyr1474, Gln1482, Ile1484, and Lys1587 exhibited similar stable hydrophobic and van der Waals interactions with the small molecule throughout most of the MD simulation process. The stable interactions of these residues further corroborate that Peimisine maintains a stable binding with either the IFI27 or OTOF protein throughout the MD process.