Integrating network pharmacology, molecular docking, and molecular dynamics simulations to explore potential compounds and mechanisms of Coptis chinensis in treating streptococcal infections

Screening of effective active ingredients and target prediction in c. chinensis

In this study, a multi-database integrative strategy was employed to systematically identify the bioactive compounds of C. chinensis and predict their potential targets, ensuring both data comprehensiveness and accuracy. Initially, the TCMSP database (https://www.tcmsp-e.com/tcmsp.php) was used to screen the active compounds of C. chinensis, with oral bioavailability (OB) ≥ 20% and drug-likeness (DL) ≥ 0.1 set as selection criteria to ensure favorable pharmacokinetic properties. To address potential limitations of a single database and to increase compound coverage, data were further supplemented with information from the TCMID database and relevant literature.

Target prediction for the selected bioactive compounds was conducted using five databases: tCMSP (old.tcmsp-e.com), DrugBank (go.drugbank.com), SwissTargetPrediction (swisstargetprediction.ch), BATMAN-TCM (ncpsb.org.cn), and TargetNet (targetnet.scbdd.com). Specific thresholds were applied where applicable, including Probability > 0 for SwissTargetPrediction and Score cutoff > 20% with P < 0.05 for BATMAN-TCM. These databases incorporate diverse methodologies—including experimental data, computational prediction, machine learning, and network pharmacology—complementing each other to enhance both the breadth and reliability of target identification. To standardize data formats, target names retrieved from TCMSP and TargetNet were converted into official gene symbols using the UniProt database. All predicted targets were subsequently integrated and duplicates removed, resulting in a comprehensive target library for C. chinensis.

Target acquisition and key target intersection in streptococcal disease

To comprehensively collect streptococcal disease–associated targets, four databases were queried: GeneCards (http://www.genecards.org), OMIM (https://omim.org/), TTD (https://db.idrblab.net/ttd/), and DisGeNET (http://www.disgenet.org). These databases encompass diverse layers of information, including genetic, clinical, and pharmacological associations, thereby providing a multidimensional overview of disease-related genes. After integrating the results and removing duplicate entries, a relatively comprehensive dataset of streptococcal disease–related targets was constructed. The predicted targets of C. chinensis bioactive compounds were then intersected with the disease-related targets to identify potential therapeutic targets.

PPI network construction and hub gene identification

The intersection of predicted drug targets and disease-related targets was input into the STRING database (https://cn.string-db.org) with the species set to Homo sapiens and the minimum required interaction score set to high confidence (>0.7); other parameters were kept at default. After removing disconnected nodes (PPI) network data were exported in TSV format and imported into Cytoscape 3.9.1 for network construction and visualization.

To identify potential key targets, topological features of the network nodes were further analyzed using the CytoNCA plugin in Cytoscape, focusing on three main metrics: Degree, Betweenness Centrality, and Closeness Centrality. Based on a comprehensive evaluation of these metrics, the top 10 hub genes were selected as key candidate targets for subsequent analyses.

GO and KEGG pathway enrichment analysis

To elucidate the functional characteristics and associated signaling pathways of the potential targets, GO and KEGG enrichment analyses were performed on 180 intersecting genes using the DAVID database (https://davidbioinformatics.nih.gov/). The analyses were conducted through the “Functional Annotation” module, with the species set to Homo sapiens and the background gene list corresponding to the human reference genome. GO analysis included annotations in three categories: biological process (BP), cellular component (CC), and molecular function (MF), while KEGG analysis was used to identify relevant signaling and disease pathways. Given the moderate size of the gene set, a significance threshold of p < 0.05 without multiple testing correction was applied to retain more biologically meaningful entries for subsequent analyses. GO results were ranked by gene count and the top ten terms of each category (BP, CC, MF) were selected for visualization. All KEGG pathways meeting the criteria were included. Enrichment results were visualized using the Bioinformatics online platform (https://www.bioinformatics.com.cn) through GO bubble plots and KEGG classification bar charts, illustrating the functional distribution and pathway characteristics of the targets.

Construction of drug ingredient-target-pathway network

The top 10 pathways with the most targets identified in the KEGG pathway enrichment analysis, along with target association information and the corresponding numbers of active ingredients from C. chinensis, were imported into Cytoscape 3.9.1 software to construct a drug ingredient-target-pathway network diagram. This network diagram can provide insights into the mechanisms by which C. chinensis exerts its therapeutic effects.

Molecular docking validation

To evaluate the binding potential between candidate active compounds and key target proteins, this study selected 10 critical protein targets, whose three-dimensional structures were primarily obtained from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB; https://www.rcsb.org/) database. Priority was given to crystal structures with co-crystallized ligands or resolutions better than 2.5 Å. Active sites were defined based on the co-crystallized ligand binding pockets and relevant literature reports, and further refined with grid predictions from AutoDock Tools to ensure accuracy and biological relevance of the docking regions. Molecular docking of 15 potential active compounds was performed using AutoDock 4.2, resulting in a total of 28 protein–ligand complexes. Binding free energy was used as the primary evaluation metric to preliminarily screen the binding affinities of all complexes.

Molecular dynamics simulation and Molecular Mechanics Poisson-Boltzmann Surface Area binding free energy calculation

Protein–ligand complexes with favorable binding free energies from molecular docking screening were further subjected to molecular dynamics (MD) simulations using GROMACS 2024.2 to evaluate their stability and binding behavior under dynamic physiological conditions. The protein was parameterized with the CHARMM36-jul2021 force field, while ligand parameters were generated using CGenFF. TIP3P water molecules were used as the solvent model. Following energy minimization and equilibration under NVT and NPT ensembles, production MD simulations were conducted for 200 ns at 310 K and 1 atm. Each simulation was independently repeated three times to ensure data reliability. Stability and flexibility of the complexes were assessed by analyzing RMSD, RMSF, hydrogen bond counts, and radius of gyration (Rg) from the simulation trajectories (Mejia-Gutierrez et al., 2021). Detailed binding mode analyses were based on conformations from the final stage of MD simulations, focusing on hydrogen bonds, hydrophobic interactions, and π-π stacking between ligands and key residues within the protein binding pockets to elucidate the molecular basis of stable binding. Additionally, binding free energies were calculated using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method over the 191–200 ns interval, providing a quantitative theoretical basis for binding affinity evaluation.

Active ingredients of c. chinensis and potential targets for streptococcal disease

From the TCMSP database, 48 active ingredients of Coptis chinensis were retrieved. After applying the screening criteria of OB > 20% and DL > 0.1, 24 active ingredients were selected for further analysis. The target genes for these ingredients were identified using five databases: 202 targets from TCMSP, 47 from DrugBank, 621 from SwissTargetPrediction, 101 from BATMAN-TCM, and 260 from TargetNet. After removing duplicates and converting protein names to gene symbols, a total of 871 unique targets were identified.

For streptococcal disease, 1,198 target genes were gathered from DisGeNET, supplemented by a literature review. Additional relevant targets were retrieved from GeneCards, OMIM, and TTD databases. After eliminating duplicates and converting protein names to gene symbols, a total of 1,210 target genes related to streptococcal disease were identified.

When intersecting the 871 Coptis chinensis targets with the 1,210 streptococcal disease targets, 180 common targets were identified, which could potentially be involved in Coptis chinensis’s therapeutic effect on streptococcal disease (Fig. 2, Table S1).

Construction of protein–protein interaction (PPI) network and hub gene screening

A protein–protein interaction (PPI) network was constructed using the STRING database, comprising 179 nodes and 1,546 edges, with an average node degree of 17.2 and an average local clustering coefficient of 0.52 (Fig. 3A). In this network, nodes represent protein targets, and edges indicate potential functional or physical interactions between targets. The size and color intensity of each node are proportional to its degree value, with higher-degree nodes positioned closer to the network core.

Topological properties of the network nodes were systematically analyzed using the CytoNCA plugin. Node importance was comprehensively assessed based on Degree, Betweenness, and Closeness centrality metrics, leading to the identification of 10 key hub genes: TP53, ACTB, AKT1, TNF, HIF1A, IL6, MYC, EGFR, IL1B, and ALB. Among these, TP53 exhibited the highest degree value (125) and was located at the network core, suggesting its potential pivotal regulatory role in the anti-streptococcal mechanism of C. chinensis (Fig. 3B, Table S2).

GO and KEGG pathway enrichment analysis

To further investigate the functional characteristics and biological pathways involving the 180 screened common target genes, GO and KEGG enrichment analyses were performed (Fig. 4, Table S3). GO analysis revealed that these targets were mainly involved in key biological processes (BP) such as immune response, signal transduction, and cellular regulation, including “inflammatory response”, “cellular response to lipopolysaccharide”, “protein phosphorylation”, and “positive regulation of transcription”, suggesting their potential roles in modulating host immunity and pathogen-related signaling pathways. In terms of cellular components (CC), targets were enriched in “extracellular space”, “plasma membrane”, “extracellular exosome”, and “cytoplasm”, indicating their important roles in signal sensing and transduction. At the molecular function (MF) level, significant enrichment in “protein kinase activity”, “enzyme binding”, “ATP binding”, and “protein binding” was observed, implying that these targets may regulate infection-related processes through kinase-mediated signaling.

KEGG pathway classification analysis revealed that the common targets participate in multiple biological systems, including metabolism, environmental information processing, cellular processes, organismal systems, and human diseases. Among these, environmental information processing was dominated by “signal transduction” pathways (19 pathways); organismal systems showed significant enrichment in “immune system” (20 pathways), endocrine system (15 pathways), and nervous system (five pathways); cellular processes encompassed pathways related to cell community behaviors, cell growth and death, and substance transport. In the disease category, enriched pathways were mainly associated with cancer (26 pathways), viral and bacterial infections (22 pathways), and parasitic infections (six pathways), indicating strong enrichment trends. These findings are highly consistent with the GO functional analysis, highlighting the core roles of common targets in immune regulation and signal transduction. Furthermore, many of these pathways include hub targets identified in the protein–protein interaction network analysis (AKT1, MAPK3, TNF), suggesting that these proteins may act as key nodes mediating signal transduction and constitute crucial pathway hubs for C. chinensis intervention against streptococcal infection.

Drug component-target-pathway network analysis

To better visualize the relationships between C. chinensis components, target genes, and Streptococcus-related signaling pathways, a drug component-target-disease pathway network was constructed, comprising 132 nodes and 807 edges (Fig. 5). In this network, blue square nodes represent common targets, green oval nodes indicate Streptococcus-associated signaling pathways, and red hexagonal nodes denote the potential active components of Coptis chinensis. The network analysis revealed that several active compounds, including quercetin, coptisine (DPEC), epiberberine, and palmatine, were associated with the highest number of targets. These active components may serve as the key bioactive compounds responsible for the therapeutic effects of C. chinensis against Streptococcus infections. Through multiple mechanisms, they collectively contribute to the pharmacological activity of C. chinensis, highlighting its potential clinical applications.

Molecular docking and binding energy analysis

The docking results indicated that all complexes exhibited binding free energies lower than −6 kcal/mol, suggesting strong interactions between these active compounds and their target proteins (Table 1). Notably, seven complexes showed binding free energies below −9 kcal/mol, indicating exceptionally high binding affinities. Among them, the Palmidin A-SRC complex demonstrated the strongest binding affinity, with a binding free energy of −10.7 kcal/mol. Additionally, quercetin, epiberberine, and tetrandrine were identified as key ligand molecules, while SRC, MMP9, and AKT1 emerged as critical target proteins for these active compounds.

Molecular dynamics simulation and MM-PBSA binding free energy calculation

Palmidin A interacts with the binding pocket of SRC protein but exhibits weak binding affinity and moderate overall stability. As shown in Fig. 6A, the small molecule remains relatively fixed within the binding pocket, while the overall protein conformation undergoes significant changes. Centroid distance analysis indicates that the distance between the small molecule and the active site ranges from 1.75 to 2.75 nm, suggesting a rigid binding mode, whereas the protein backbone RMSD fluctuates considerably (Fig. 6C). Residue RMSF trajectories reveal high flexibility in the binding pocket region, particularly with pronounced fluctuations at the C-terminal (Fig. 6D). Additionally, the hydrogen bond count and radius of gyration remain stable, indicating no conformational collapse in the system (Figs. 6E, 6F). Binding mode analysis (Fig. 6G) suggests that Palmidin A primarily interacts with key residues through hydrogen bonding, hydrophobic interactions, and π–π stacking, maintaining a certain degree of binding at the end of the simulation. MM-PBSA calculations (191–200 ns) reveal a binding free energy (ΔG) of approximately −24.08 kJ/mol for Palmidin A with SRC protein, corresponding to a binding affinity (Ki) of approximately 6.06 µM, indicating weak binding. Among the interaction components, van der Waals interactions (−135.24 kJ/mol) contribute the most, whereas electrostatic interactions (−27.02 kJ/mol) are relatively minor (Table 2).

In contrast, the MMP9-quercetin complex demonstrates higher stability over the 200 ns simulation. Conformational superposition analysis (Fig. 7A) shows minimal changes in the small molecule within the active pocket, with quercetin dynamically anchored through a hydrophobic core (Phe110, Leu188) and a hydrogen bond network (Tyr219). Its benzene ring forms π–π stacking with His226 (Fig. 7G). The protein backbone RMSD stabilizes at 0.3–0.4 nm after 30 ns (Fig. 7C), while the RMSD of quercetin remains below 0.1 nm, indicating rapid conformational convergence. Residues at the binding site (Leu181–Tyr223) exhibit low flexibility (Fig. 7D), and the distance between the small molecule and the active site remains stable at 0.3–0.5 nm (Fig. 7B). The hydrogen bond count and radius of gyration show minimal fluctuations (Figs. 7E, 7F), further supporting the tight binding of the complex. MM-PBSA calculations (191–200 ns) indicate that the binding free energy (ΔG) of quercetin with MMP9 is −107.9 kJ/mol, with a corresponding binding affinity of 5.77 × 10⁻⁹ nM, suggesting strong binding. The binding energy is mainly derived from van der Waals forces (−179.6 kJ/mol) and electrostatic interactions (−48.9 kJ/mol) (Table 2).

To assess the stability of the simulations, three independent simulations were conducted for each complex system (Figs. 8 and 9). Multi-trajectory analysis indicates that the stability of the Palmidin A-SRC complex is independent of initial conditions, with its binding mode achieved through the coordination of dynamic hydrogen bond networks and protein conformational adjustments (Fig. 8A). Meanwhile, the three RMSD curves of the MMP9-quercetin complex nearly overlap (Fig. 9A), with a maximum centroid distance deviation of only 0.3 nm (Fig. 9C). The hydrogen bond network (Fig. 9D) and radius of gyration (Fig. 9E) exhibit consistent dynamic patterns, further confirming the stability of the complex. These findings demonstrate good reproducibility of the simulation data, aligning with the technical standards of molecular dynamics research and ensuring the reliability of the study conclusions (Shirts et al., 2017).