The Potential Pharmacological Mechanism of Medicinal Herbs for the Treatment of Cough and Cold via Network Pharmacology, Molecular Docking and Simulation

Cough and cold are among the most common upper respiratory tract disorders, primarily of viral origin, and remain significant causes of morbidity across all age groups. While symptoms are often self-limiting and seasonal, they can impair quality of life and lead to complications in vulnerable populations. Traditional medicinal systems have long used botanicals for respiratory care because of their multi-target therapeutic potential. This study investigated five medicinal herbs, Curcuma longa, Glycyrrhiza glabra, Terminalia chebula, Adhatoda vasica, and Ocimum sanctum using network pharmacology and molecular docking to elucidate their mechanisms against cough and cold. Bioactive compounds were retrieved from phytochemical databases (IMPPAT, Dr. Duke’s), screened for drug-likeness via SwissADME, and mapped to potential human protein targets using SwissTarget Prediction and BindingDB. Disease-related targets were identified from DisGeNET and GeneCards, and overlapping targets were analyzed using cytoscape-based bioactive-target and protein–protein interaction (PPI) networks. GO and KEGG enrichment analyses were performed, followed by molecular docking and simulations with key hub proteins. Ninety-three compounds met the drug-likeness criteria, yielding 56 overlapping targets including AKT1, EGFR, SRC, PTGS2, MMP9, MMP2, and KDR. Enrichment analysis indicated the involvement of PI3K-Akt, EGFR, and focal adhesion pathways, which are central to inflammation, immune regulation, and epithelial repair. Molecular docking revealed strong (−5.2 to −9.66 kcal/mol) for quercetin, luteolin, eugalic acid, and kaempferol with hub proteins, suggesting potent multi-target interactions. MMP2/9 protein with quercetin, luteolin, and kaempferol complexes showed good stability during molecular dynamics (MD) simulations. Estimation of the binding free energy using the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) method validated the inhibitor potential of the identified compounds. These findings provide molecular evidence for the synergistic multi-pathway effects of the studied herbs, supporting their use as safe, affordable, and evidence-based complementary therapies for cough and cold.

Introduction

Cough and cold are among the most common and prevalent upper respiratory tract disorders, often occurring together as part of a viral infection. Cough is caused by various respiratory viruses and its symptoms result from an innate immune response in the respiratory epithelium. Influencing factors include age, immunity, stress, and seasonal variation [1]. The cough acts as a vital protective reflex, clearing the airways of mucus, pathogens, and irritants. Although transient, it can become chronic, lasting over eight weeks in adults, due to nerve hypersensitivity and abnormal neural signaling in both central and peripheral pathways. Persistent coughing can lead to inflammation and impairment in Respiratory function [2]. The common cold is defined by the National Institute for Health and Care Excellence (NICE) as a self-limiting upper respiratory infection. It is characterized by symptoms, such as nasal congestion, sneezing, and sore throat. Despite its typically mild and self-resolving nature, the common cold remains the most widespread human illness, reflecting the complex interplay between respiratory pathogens and host immune responses [3].

Historically, medicinal herbs have played a crucial role in conventional healthcare systems, serving as the cornerstone of therapeutic approaches in many cultures. For centuries, traditional medicinal systems, such as Ayurveda, Western herbalism, Unani, and Traditional Chinese Medicine (TCM), have used botanical remedies to treat a variety of illnesses, especially respiratory tract infections such as colds and coughs. Although these illnesses are frequently mild and self-limiting, they can seriously hinder day-to-day functioning, disrupt sleep, and potentially result in more severe complications in susceptible groups, including children, the elderly, and those with compromised immune systems [4]. Medicinal herbs provide a poly pharmacological approach, in contrast to traditional pharmaceutical drugs, which usually concentrate on isolated active compounds. This indicates that several components of a single herb may function in concert to provide medicinal benefits. Cough and cold, which involve intricate pathophysiological processes such as viral infection, inflammation, mucus overproduction, bronchoconstriction, and immune dysregulation, benefit greatly from such multi-targeted mechanisms [5]. In addition to treating symptoms, these pharmacological actions may help the body’s natural capacity to heal infections more quickly [6].

The consumer desire for safer, natural alternatives to synthetic drugs has led to a resurgence of interest in herbal medicines or supplements in recent decades. Furthermore, the increasing antibiotic resistance and side effects of traditional treatments have prompted more research into botanicals in integrative medicine. The effectiveness of several medicinal herbs in treating upper respiratory tract infections is currently supported by numerous clinical and preclinical studies [7]–[9]. The current study aimed to explore the level of contribution of well-known medicinal herbs, such as Curcuma longa (turmeric), Glycyrrhiza glabra (licorice), Terminalia chebula (haritaki), Adhatoda vasica (vasaka), and Ocimum sanctum (tulsi), extracted and formulated into respiratory wellness tablets (Koff-N-Cold, Greenspace). These plants represent a scientifically and traditionally validated group of medicinal herbs with multi-target effects on respiratory health. They modulate central pathways (NF-κB, MAPKs, PI3K/AKT, TGF-β/Smad, and Nrf2) [10]–[13] and gene networks involved in inflammation, mucus hypersecretion, oxidative stress, immune dysregulation, fibrosis, and even viral infections. Unlike single-target drugs, they offer a holistic and synergistic approach to diseases, such as cough, cold, bronchitis, asthma, and Chronic Obstructive Pulmonary Disease (COPD). Their diverse bioactive profiles make them promising candidates for use in novel anti-inflammatory, antitussive, antiviral, and antioxidant therapies. Studying these plants can uncover new mechanisms, identify bioactive molecules, and validate safe, affordable alternatives or complements to conventional treatments, addressing both unmet clinical needs and growing interest in evidence-based traditional medicine by focusing on their phytochemical and mechanisms of action through network pharmacology and molecular docking simulation.

Materials and Methods

Collection of Bioactive Compounds

The bioactive constituents of Curcuma longa (turmeric), Glycyrrhiza glabra (licorice), Terminalia chebula (haritaki), Adhatoda vasica (vasaka), and Ocimum sanctum (tulsi) herbs were identified using freely accessible phytochemical repositories, including the Indian Medicinal Plants, Phytochemistry and Therapeutics (IMPPAT) database [14], [15] and Dr. Duke’s Phytochemical and Ethnobotanical Database [14]. These resources provide comprehensive information on the phytochemicals present in medicinal plants. The retrieved compounds were subsequently evaluated for their pharmacokinetic properties, focusing on drug-likeness and lead-likeness, using SwissADME. All phytoconstituents were evaluated against Lipinski’s Rule of Five (RO5) to assess drug-likeness. The criteria included a molecular weight below 500 Da, log P value under 5, topological polar surface area (TPSA) less than 140 Ų, and a maximum of 10 hydrogen bond acceptors and 5 hydrogen bond donors. In addition, compounds, including lead-likeness and oral bioavailability, were further assessed for pharmacokinetic suitability. Compounds that fulfilled these criteria were prioritized for subsequent computational analyses.

Screened Bioactive Targets Prediction

Bioactive targets were identified using SwissTargetPrediction [15] and BindingDB [16] databases. Potential protein targets for the selected bioactive compounds were identified by querying the canonical SMILES and 2D-SDF structures retrieved from the PubChem database. Target prediction was performed using SwissTargetPrediction and BindingDB, with the species restricted to Homo sapiens. From SwissTargetPrediction, only targets with a probability score ≥ 0.50 were retained, while in BindingDB, targets were selected based on a chemical similarity threshold of ≥ 0.50. The resulting lists of predicted targets were compiled and curated to remove redundancies, ensuring a nonduplicative set of high-confidence target proteins for further analysis. All de-duplicated targets were standardized and unified using UniProtKB databases and saved as potential targets for further analysis [17], [18]

Collection of Potential Cough and Cold Related Protein Targets

To identify key molecular targets involved in the treatment of cough and cold, protein targets were systematically retrieved from publicly available databases including DisGeNET and GeneCards [18], using specific keywords such as “bronchitis,” “respiratory allergy,” “tussive,” and “cough and cold.” Simultaneously, the predicted targets of the bioactive compounds present in the five medicinal herbs were collected from phytochemical databases. The overlap between disease-associated and compound-related targets was visualized using the BioinfoGP platform through a Venn diagram analysis. To ensure data accuracy and standardization, all target entries were curated using the UniProt database. This integrative approach facilitated the identification of shared molecular targets, highlighting the potential mechanisms by which medicinal herbs may exert therapeutic effects against cough and cold symptoms.

Bioactive-Target Network Construction

A bioactive-target (BA-TAR) interaction network was constructed using Cytoscape v3.9 [19] to investigate the multi-component and multi-target relationships between the phytoconstituents and their potential therapeutic targets associated with cough and cold. This network facilitates the visualization of complex interactions and enables the identification of key nodes within the system. To evaluate the topological importance of each node, the network was subjected to centrality analysis using the cytoNCA plugin, applying the “degree value” parameter as a primary metric. Bioactive compounds exhibiting high degree values were recognized as core components, suggesting their pivotal role in modulating the key targets involved in the pathophysiology of cough and cold.

Protein–Protein Interactions (PPI) Network Construction and Hub Gene Identification

To explore the potential interactions among the identified anti-cough and cold targets, the intersecting genes were analyzed using the STRING v11.5 database [20], with the organism set to Homo sapiens and a minimum confidence score threshold of 0.4 to ensure biologically relevant associations. The resulting protein–protein interaction (PPI) data were transferred into Cytoscape v3.9.1 for network visualization and topological assessment. Key network parameters, including degree, betweenness, and closeness centralities, were calculated to evaluate the importance of each node within the PPI network. Hub genes were identified by selecting nodes with centrality values that exceeded their respective median thresholds. Additionally, the MCODE plugin [21] was employed to perform module detection, allowing the identification of densely connected subnetworks that may represent functional clusters relevant to the therapeutic mechanisms of medicinal herbs.

Gene Ontology (GO) and KEGG Pathway Enrichment Analyses

Functional enrichment analysis was performed on the identified core protein targets using DAVID online tools to elucidate the potential mechanisms by which bioactive compounds exert therapeutic effects against cough and cold. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted by inputting the official gene symbols of intersecting targets, with Homo sapiens selected as the reference organism. The enrichment results were filtered based on statistical significance criteria, including the FDR-adjusted p-value Q (Q ≤ 0.05), and sorted by gene count in descending order. The top 10 terms from each GO category, including biological process (BP), cellular component (CC), and molecular function (MF), as well as the top 10 KEGG pathways, were selected for interpretation [22]. To visually represent the functional relevance of these targets, an SR Plot was employed to generate bubble plots for GO terms and Sankey diagrams for KEGG pathways, facilitating a clearer understanding of the biological pathways and molecular functions implicated in cough and cold treatments.

Molecular Docking

Molecular docking was performed to evaluate the binding interactions between key phytoconstituents and core protein targets identified through protein–protein interaction (PPI) network analysis. The crystal structures of eight hub proteins, AKT1 (1UNQ), EGFR (1XKK), SRC (2H8H), PTGS2 (5F19), MMP9 (1GKC), MMP2 (1CK7), ERBB2 (3PP0), and KDR (3VHE), were obtained from the Protein Data Bank (PDB), prioritizing human-derived X-ray structures with high resolution, co-crystallized ligands, and structural completeness. Phytoconstituents including quercetin, ellagic acid, kaempferol, apigenin, and luteolin were sourced in 2D format from the PubChem database and converted to 3D structures using OpenBabel. Protein preparation for docking was carried out using AutoDock 4.2.1, involving removal of non-protein entities and water molecules, extraction of native ligands, addition of polar hydrogens, and correction of missing atoms. Docking simulations were executed by defining grid boxes around the active sites of each protein, and the interaction results were visualized and analyzed using BIOVIA Discovery Studio Visualizer to determine the binding conformations and affinities of the phyto ligands. After docking, conformations with the best docking scores and RMSD (root-mean-square deviation) values ≤ 1 were selected for further studies [23].

Molecular Simulation

Molecular dynamics simulations of the protein–ligand complexes were performed using the Desmond simulation package. All computations were executed on a workstation configured with an Intel i7 processor and 32 GB of RAM, operating on Ubuntu 18.0. Prior to simulation, the protein–ligand systems were prepared using the Protein Preparation Wizard in Maestro, which involved the addition of hydrogen atoms, optimization of bond orders, and correction of structural inconsistencies. The solvated systems were generated within an orthorhombic box, maintaining a 10 Å buffer distance around complexes, and the explicit water environment was modeled using the TIP3P water model. To ensure electrostatic neutrality, counter ions (Na+ and Cl−) were introduced. The force field parameters were assigned using OPLS-2005.

The prepared systems were subjected to MD simulations for a total duration of 100 ns under NPT ensemble conditions, maintaining a constant temperature of 300 K and a pressure of 1.01 bar using a Nose Hoover chain thermostat and Martyna–Tobias–Klein barostat algorithms, respectively. Post-simulation analyses were performed using the Simulation Interaction Diagram (SID) module, which provided detailed trajectories of structural and dynamic parameters, including ligand–protein interaction, torsional flexibility, root mean square deviation (RMSD), root mean square fluctuation (RMSF), secondary structure histogram (SSH), and secondary structure evolution (SSE) [24].

Calculation of Binding Free Energy

The Molecular Mechanics Generalized Born Surface Area (MM-GBSA) binding free energy of the simulation trajectories was calculated using the following formula [25]:

$$\begin{array}{r}\mathrm{\Delta }{\mathrm{G}}_{\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{d}}=\mathrm{\Delta }{\mathrm{G}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{x}}-\mathrm{\Delta }{\mathrm{G}}_{\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}\mathrm{o}\mathrm{r}}-\mathrm{\Delta }{\mathrm{G}}_{\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}}\end{array}$$

Here, ΔG_bind indicates the binding free energy, whereas ΔG_complex, ΔG_receptor, and ΔG_ligand refer to the free energy of the complex, target proteins, and ligands, respectively.

Results

Screened Bioactive Compounds

The active phytoconstituents of the five medicinal herbs showed an initial pool of 524 unique compounds: 168 from Curcuma longa, 157 from Glycyrrhiza glabra (licorice), 32 from Terminalia chebula (haritaki), 51 from Adhatoda vasica (vasaka), and 116 from Ocimum sanctum (tulsi). These compounds were subsequently evaluated for drug-likeness and lead-likeness using SwissADME. Applying filters such as Lipinski’s Rule of Five and allowing a maximum of one violation for lead-likeness, 93 compounds were found to meet the essential pharmacokinetic and physicochemical criteria for oral bioavailability. Among these, 31 compounds from C. longa, 42 from G. glabra, 2 from T. chebula, 15 from A. vasica and 9 from O. sanctum were shortlisted for further pharmacological analysis. A detailed overview of the screened phytochemicals is provided in Table I.

Table I. All the Compounds Which Has Passed Drug Likeness and Lead Likeness

Compound name	PubChemID	Bioactive SMILES	Molecular weight (g/mol)
Ellagic acid	5281855	C1=C2C3=C(C(=C1O)O)OC(=O)C4=CC(=C(C(=C43)OC2=O)O)O	302.19
Kaempferol	5280863	C1=CC(=CC=C1C2=C(C(=O)C3=C(C=C(C=C3O2)O)O)O)O	286.24
Quercetin	5280343	C1=CC(=C(C=C1C2=C(C(=O)C3=C(C=C(C=C3O2)O)O)O)O)O	302.23
Luteolin	5280445	C1=CC(=C(C=C1C2=CC(=O)C3=C(C=C(C=C3O2)O)O)O)O	286.24
Apigenin	5280443	C1=CC(=CC=C1C2=CC(=O)C3=C(C=C(C=C3O2)O)O)O	270.24
Cirsilineol	162464	COC1=C(C=CC(=C1)C2=CC(=O)C3=C(C(=C(C=C3O2)OC)OC)O)O	344.3
Eupalitin	5748611	COC1=C(C(=C2C(=C1)OC(=C(C2=O)O)C3=CC=C(C=C3)O)O)OC	330.29
Cianidanol	9064	C1[C@@H]([C@H](OC2=CC(=CC(=C21)O)O)C3=CC(=C(C=C3)O)O)O	290.27
(-)-Epicatechin	72276	C1[C@H]([C@H](OC2=CC(=CC(=C21)O)O)C3=CC(=C(C=C3)O)O)O	290.27
Isoliquiritigenin	638278	C1=CC(=CC=C1/C=C/C(=O)C2=C(C=C(C=C2)O)O)O	256.25
Pinocembrin	68071	C1[C@H](OC2=CC(=CC(=C2C1=O)O)O)C3=CC=CC=C3	256.25
Formononetin	68071	C1[C@H](OC2=CC(=CC(=C2C1=O)O)O)C3=CC=CC=C3	256.25
Naringetol	439246	C1[C@H](OC2=CC(=CC(=C2C1=O)O)O)C3=CC=C(C=C3)O	272.25
Glabroiso-flavanone A	11221431	Oc1ccc(c(c1)O)C1COc2c(C1=O)ccc1c2C=CC(O1)(C)C	338.4
Liquiritigenin	114829	C1[C@H](OC2=C(C1=O)C=CC(=C2)O)C3=CC=C(C=C3)O	256.25
7-Acetoxy-2-methylisoflavone	268208	CC1=C(C(=O)C2=C(O1)C=C(C=C2)OC(=O)C)C3=CC=CC=C3	294.3
7-Methoxy-2-methyl-3-phenyl-4H-chromen-4-one	354368	CC1=C(C(=O)C2=C(O1)C=C(C=C2)OC)C3=CC=CC=C3	266.29
[(7S,8S)-7-hydroxy-5,6,7,8-tetrahydro-3H-pyrrolizin-1-yl]methyl 2-hydroxy-2-[(1S)-1-hydroxyethyl]-3-methylbutanoate	10197	C[C@@H](C(C(C)C)(C(=O)OCC1=CCN2[C@@H]1[C@H](CC2)O)O)O	299.36
Kumatakenin	5318869	COC1=CC(=C2C(=C1)OC(=C(C2=O)OC)C3=CC=C(C=C3)O)O	314.29
Texasin	5281812	COC1=CC=C(C=C1)C2=COC3=CC(=C(C=C3C2=O)O)O	284.26
Prunetin	5281804	COC1=CC(=C2C(=C1)OC=C(C2=O)C3=CC=C(C=C3)O)O	284.26
Glyzarin	44257206	CC1=C(C(=O)C2=C(O1)C(=C(C=C2)O)C(=O)C)C3=CC=CC=C3	294.3
Glabrone	5317652	CC1(C=CC2=C(O1)C=CC(=C2O)C3=COC4=C(C3=O)C=CC(=C4)O)C	336.3
Glyzaglabrin	5317777	C1OC2=C(O1)C(=C(C=C2)C3=COC4=C(C3=O)C=CC(=C4)O)O	298.25
7-Hydroxy-2-methyl-3-phenyl-4H-chromen-4-one	5380976	CC1=C(C(=O)C2=C(O1)C=C(C=C2)O)C3=CC=CC=C3	252.26
(1E,4E)-1-(4-Hydroxy-3-methoxyphenyl)-5-(4-hydroxyphenyl)-1,4-pentadiene-3-one	10469828	COC1=C(C=CC(=C1)/C=C/C(=O)/C=C/C2=CC=C(C=C2)O)O	296.3
4,10-Epizedoarondiol	24834047	CC(=C1C[C@H]2[C@@H](CC[C@@]2(C)O)[C@@](CC1=O)(C)O)C	252.35
Bisacurone	14287397	C[C@@H](CC(=O)C=C(C)C)[C@H]1C[C@@H]([C@@](C=C1)(C)O)O	252.35
Bisdemethoxycurcumin	5315472	C1=CC(=CC=C1/C=C/C(=O)CC(=O)/C=C/C2=CC=C(C=C2)O)O	308.33
Demethoxycurcumin	5469424	COC1=C(C=CC(=C1)/C=C/C(=O)CC(=O)/C=C/C2=CC=C(C=C2)O)O	338.35
Piperine	638024	C1CCN(CC1)C(=O)/C=C/C=C/C2=CC3=C(C=C2)OCO3	285.34
Procurcumadiol	14633012	CC1=CC(=O)C(=C(C)C)C[C@]2([C@H]1CC[C@]2(C)O)O	250.33
Anisotine	442884	CNC1=C(C=C(C=C1)C2CCN3C2=NC4=CC=CC=C4C3=O)C(=O)OC	349.38
Vasicolinone	627712	CN(C)C1=CC=CC=C1C2CCN3C2=NC4=CC=CC=C4C3=O	305.4
Vasicoline	626005	CN(C)C1=CC=CC=C1C2CCN3C2=NC4=CC=CC=C4C3	291.4
Methyl 2-(methylamino)-5-(1,2,3,9-tetrahydropyrrolo[2,1-b]quinazolin-3-yl)benzoate	5316460	CNC1=C(C=C(C=C1)C2CCN3C2=NC4=CC=CC=C4C3)C(=O)OC	335.4
Desmethoxyaniflorine	101995827	CN(c1ccccc1[C@]1(O)CCn2c1nc1ccccc1c2=O)C	321.38
Vasnetine	9997066	COC(=O)c1ccccc1NC1CCn2c1nc1ccccc1c2=O	335.36
Anthraquinone	6780	C1=CC=C2C(=C1)C(=O)C3=CC=CC=C3C2=O	208.21
Rosmarinic acid	5281792	C1=CC(=C(C=C1C[C@H](C(=O)O)OC(=O)/C=C/C2=CC(=C(C=C2)O)O)O)O	360.3
Isoflavone	72304	C1=CC=C(C=C1)C2=COC3=CC=CC=C3C2=O	222.24
Methylenedioxybenzoyl ethyl PABA	785868	CCOC(=O)C1=CC=C(C=C1)NC(=O)C2=CC3=C(C=C2)OCO3	313.3
Licoriphenone	21591149	CC(=CCC1=C(C=C(C(=C1OC)CC(=O)C2=C(C=C(C=C2)O)O)O)OC)C	372.4
Glycyrrhisoflavanone	5317762	CC1(C=CC2=C(O1)C(=CC(=C2)C3COC4=C(C3=O)C(=CC(=C4)O)OC)O)C	368.4
Semilicoisoflavone B	5481948	CC1(C=CC2=C(O1)C(=CC(=C2)C3=COC4=CC(=CC(=C4C3=O)O)O)O)C	352.3
Glabridin	124052	CC1(C=CC2=C(O1)C=CC3=C2OC[C@H](C3)C4=C(C=C(C=C4)O)O)C	324.4
Glabranin	124049	CC(=CCC1=C2C(=C(C=C1O)O)C(=O)C[C@H](O2)C3=CC=CC=C3)C	324.4
7,4’-Dihydroxyflavan	158280	C1CC2=C(C=C(C=C2)O)O[C@@H]1C3=CC=C(C=C3)O	242.27
Kanzonol R	131753027	CC(=CCC1=C(C=CC(=C1O)C2CC3=C(C=C(C=C3OC)O)OC2)OC)C	370.4
Glabrene	480774	CC1(C=CC2=C(C=CC(=C2O1)C3=CC4=C(C=C(C=C4)O)OC3)O)C	322.4
Glycyrin	480787	CC(=CCC1=C(C=C2C(=C1OC)C=C(C(=O)O2)C3=C(C=C(C=C3)O)O)OC)C	382.4
Shinpterocarpin	10336244	CC1(C=CC2=C(O1)C=CC3=C2OC[C@@H]4[C@H]3OC5=C4C=CC(=C5)O)C	322.4
Glabroisoflavanone B	11405466	COc1ccc(c(c1)O)C1COc2c(C1=O)ccc1c2C=CC(O1)(C)C	352.38
Glabrocoumarin	11427657	CC1(C=CC2=C(C=CC(=C2O1)C3=CC4=C(C=C(C=C4)O)OC3=O)O)C	336.3
Liqcoumarin	11378967	CC1=CC(=O)OC2=C1C(=C(C=C2)C(=O)C)O	218.2
Lauric acid	3893	CCCCCCCCCCCC(=O)O	200.32
Pratol	5320693	COC1=CC=C(C=C1)C2=CC(=O)C3=C(O2)C=C(C=C3)O	268.26
Licoricone	5319013	CC(=CCC1=C(C=C(C(=C1OC)C2=COC3=C(C2=O)C=CC(=C3)O)O)OC)C	382.4
Licoisoflavone A	5281789	CC(=CCC1=C(C=CC(=C1O)C2=COC3=CC(=CC(=C3C2=O)O)O)O)C	354.4
Glycyrrhisoflavone	5317764	CC(=CCC1=C(C(=CC(=C1)C2=COC3=CC(=CC(=C3C2=O)O)O)O)O)C	354.4
4-Hydroxychalcon	5282362	C1=CC=C(C=C1)/C=C/C(=O)C2=CC=C(C=C2)O	224.25
Licoisoflavone B	5481234	CC1(C=CC2=C(O1)C=CC(=C2O)C3=COC4=CC(=CC(=C4C3=O)O)O)C	352.3
Licoflavonol	5481964	CC(=CCC1=C(C2=C(C=C1O)OC(=C(C2=O)O)C3=CC=C(C=C3)O)O)C	354.4
(+)-Curcumenol	167812	C[C@H]1CC[C@@H]2[C@]13CC(=C(C)C)[C@](O3)(C=C2C)O	234.33
(1E)-1,7-bis(4-hydroxy-3-methoxyphenyl)hept-1-ene-3,5-dione	10429233	COC1=C(C=CC(=C1)CCC(=O)CC(=O)/C=C/C2=CC(=C(C=C2)O)OC)O	370.4
(4S,5S)-Germacrone-4,5-epoxide	91753231	CC1=CCC[C@@]2([C@@H](O2)CC(=C(C)C)C(=O)C1)C	234.33
(6Z)-6,10-dimethyl-3-propan-2-ylidenecyclodec-6-ene-1,4-dione	14191392	CC1CC/C=C(\CC(=O)C(=C(C)C)CC1=O)/C	234.33
1,7-Bis(4-hydroxyphenyl)-1,4,6-heptatrien-3-one	10447050	C1=CC(=CC=C1/C=C/C=C/C(=O)/C=C/C2=CC=C(C=C2)O)O	292.3
2-(Hydroxymethyl)anthraquinone	87014	C1=CC=C2C(=C1)C(=O)C3=C(C2=O)C=C(C=C3)CO	238.24
2-Hepten-4-one, 6-(2-hydroxy-4-methylphenyl)-2-methyl-	10955433	CC1=CC(=C(C=C1)C(C)CC(=O)C=C(C)C)O	232.32
Bicyclo(4.1.0)heptan-3-one, 1-methyl-4-(1-methylethylidene)-7-(3-oxobutyl)-, (1S,6R,7R)-	153845	CC(=C1C[C@@H]2[C@H]([C@@]2(CC1=O)C)CCC(=O)C)C	234.33
Cinnamyl cinnamate	1550890	C1=CC=C(C=C1)/C=C/COC(=O)/C=C/C2=CC=CC=C2	264.3
Curcumin	969516	COC1=C(C=CC(=C1)/C=C/C(=O)CC(=O)/C=C/C2=CC(=C(C=C2)O)OC)O	368.4
Curcumol	14240392	C[C@H]1CC[C@@H]2[C@]13C[C@H]([C@](O3)(CC2=C)O)C(C)C	236.35
Curzerenone	3081930	CC1=COC2=C1C(=O)[C@@H]([C@@](C2)(C)C=C)C(=C)C	230.3
Cyclocurcumin	69879809	COC1=C(C=CC(=C1)/C=C/C2=CC(=O)CC(O2)C3=CC(=C(C=C3)O)OC)O	368.4
Elemicin	10248	COC1=CC(=CC(=C1OC)OC)CC=C	208.25
Epiprocurcumenol	10263440	CC1=CC(=O)C(=C(C)C)C[C@H]2[C@@H]1CC[C@]2(C)O	234.33
Ethene;2-methoxyphenol;phenol;propanedial	122130015	COC1=CC=CC=C1O.C=C.C=C.C1=CC=C(C=C1)O.C(C=O)C=O	346.4
Furanodienon	6442374	C/C/1=C/C(=O)C2=C(C/C(=C\CC1)/C)OC=C2C	230.3
Geranyl butyrate	5355856	CCCC(=O)OC/C=C(\C)/CCC=C(C)C	224.34
Germacr-1(10)-ene-5,8-dione	6441391	C[C@H]1CC/C=C(/CC(=O)[C@@H](CC1=O)C(C)C)\C	236.35
Isoprocurcumenol	14543197	CC(=C1CC2C(CCC2(C)O)C(=C)CC1=O)C	234.33
Linalyl isobutyrate	6532	CC(C)C(=O)OC(C)(CCC=C(C)C)C=C	224.34
Linalyl propionate	61098	CCC(=O)OC(C)(CCC=C(C)C)C=C	210.31
Procurcumenol	189061	CC1=CC(=O)C(=C(C)C)C[C@H]2[C@H]1CC[C@]2(C)O	234.33
Turmeronol A	11117927	CC1=C(C=C(C=C1)[C@@H](C)CC(=O)C=C(C)C)O	232.32
Vasicinol	442934	C1CN2CC3=C(C=CC(=C3)O)N=C2[C@H]1O	204.22
Vasicol	92470596	O=C1[C@H](O)CCN1Cc1ccccc1N	206.24
(3S)-3-hydroxy-2,3-dihydro-1H-pyrrolo[2,1-b]quinazolin-9-one	442935	C1CN2C(=NC3=CC=CC=C3C2=O)[C@H]1O	202.21
7-Methoxyvasicinone	101995828	COc1ccc2c(c1)c(=O)n1c(n2)[C@@H](CC1)O	232.24
Flavone	10680	C1=CC=C(C=C1)C2=CC(=O)C3=CC=CC=C3O2	222.24
Heptacosane	11636	CCCCCCCCCCCCCCCCCCCCCCCCCCC	380.7

Potential Targets in the Treatments of Cough and Cold

Among the 93 shortlisted bioactive compounds, 26 were found to have associated target predictions, whereas 62 did not yield any target predictions. Only SwissTargetPrediction targets with a probability score of ≥ 0.5, and BindingDB targets with a chemical similarity threshold of ≥ 0.5, were retained, and targets below these thresholds were excluded from the target assignment. Additionally, five compounds lacked target data entirely in both databases, indicating either insufficient structural similarity to known ligands or lack of experimental evidence in the referenced datasets.

Cough- and cold-related target proteins were retrieved from the DisGeNET and GeneCards databases, resulting in the identification of 5,696 associated genes. To determine the potential therapeutic relevance of all five medicinal herbs, a target overlap analysis was performed using Venny 2.1.0 by comparing the predicted compound-associated targets with the disease-related targets shown in Fig. 1. This intersection reveals 112 common targets. Subsequent manual curation and literature-based evaluation were conducted to assess the biological relevance of these targets in cough and cold pathophysiology. As a result, 56 unrelated targets were excluded, leaving 56 validated targets for further investigation. These targets were linked to 13 bioactive compounds, which were selected for downstream network construction and pharmacological analysis.

Fig. 1. Venn diagram showing the overlap between cough and cold-related disease targets and Koff-N-Cold tablet targets.

BA-TAR Network Construction

To explore the interaction between the bioactive compounds and overlapping targets associated with cough and cold, a bioactive–target (BA-TAR) network was constructed using Cytoscape v3.9.1. The network comprised 56 validated target proteins and 13 bioactive compounds that were integrated to visualize the multi-target pharmacological landscape. A network Analyzer was employed to quantify the network’s structural properties, revealing 74 nodes and 142 edges, with an average of 3.942 neighbors per node. Degree centrality analysis was conducted on the 26 bioactive compounds initially considered, identifying Quercetin, Kaempferol, Ellagic acid, Luteolin, and Apigenin as the top five candidates linked to 41, 20, 19, 18, and 17 targets, respectively (Table II). These compounds demonstrated significant connectivity within the network, suggesting their central role in modulating key biological targets related to cough and cold stress. Consequently, they were prioritized for subsequent molecular docking analysis. The resulting network, depicted in Fig. 2, underscores the multi-component and multi-target therapeutic potential of the herbs.

Table II. The Degree Analysis of Screened Bioactive Compounds

PubChemID	Bioactive compound	Degree
5280343	Quercetin	41
5280863	Kaempferol	20
5281855	Ellagic acid	19
5280445	Luteolin	18
5280443	Apigenin	17
10680	Flavone	7
439246	Naringenin	4
969516	Curcumin	3
5469424	Demethoxycurcumin	3
354368	7-Methoxy-2-methyl-3-phenyl-4H-chromen-4-one	3
5281792	Rosmarinic acid	3
5280378	Formononetin	2
68071	Pinocembrin	2

Fig. 2. BA-TAR network with 74 nodes and 142 edges linking screened bioactive compounds with cough and cold targets.

PPI Network Analysis

The 56 intersecting targets shared between cough and cold pathology and the bioactive compounds of the five herbs were mapped to the STRING database, with the species restricted to Homo sapiens, to construct a protein–protein interaction (PPI) network. The resulting network consisted of 56 nodes and 337 edges, with an average local clustering coefficient of 0.48 and an average node degree of 12, as depicted in Fig. 3A. To enable more detailed visualization and topological analysis, the STRING output was imported into Cytoscape v3.9.1. After removing two disconnected nodes, the refined PPI network included 54 nodes and 337 edges (Fig. 3B). To identify the key regulatory elements within the network, hub gene analysis was performed using degree centrality as the selection metric (Table III). Genes with degree values greater than the median cut-off (≥ 20) were classified as core targets. Based on this criterion, AKT1, EGFR, SRC, PTGS2, MMP9, MMP2, ERBB2, KDR, and IL2 were identified as hub genes and were subsequently selected for molecular docking studies.

Fig. 3. PPI network of potential anti-cough and cold targets. (A) STRING database PPI network (confidence score ≥ 0.40) and (B) Cytoscape v3.9.1-mapped PPI network. Nodes depicting the target and edges represent the interaction between protein targets. The darker nodes correspond to the higher degree and greater therapeutic importance.

Table III. The General Description of the PPI Network Hub Genes

Hub genes	Degree	Betweenness centrality	Closeness centrality
AKT1	37	0.160702	0.742857
EGFR	32	0.110372	0.693333
SRC	30	0.132473	0.675325
PTGS2	30	0.117895	0.675325
MMP9	30	0.081828	0.675325
MMP2	24	0.018792	0.604651
ERBB2	24	0.021742	0.604651
KDR	23	0.011017	0.597701
IL2	21	0.031676	0.577778
PDGFRB	19	0.005752	0.565217

Clusters Network Analysis

To further validate the findings from KEGG enrichment and protein–protein interaction (PPI) network analyses, a cluster network analysis was performed on the intersecting targets using the MCODE plugin in Cytoscape v3.9.1. This analysis revealed two significant modules: Cluster 1, comprising 15 nodes and 95 edges, and Cluster 2, comprising 9 nodes and 16 edges. Notably, the top four hub genes identified in Section 3.4 were distributed within these two functional clusters, reinforcing their central roles in the network. Specifically, AKT1 was located within Cluster 1, highlighting its importance in highly interconnected regulatory pathways, whereas PTGS2 emerged as the core target in Cluster 2. The organization of these targets within distinct, yet interconnected modules, as illustrated in Fig. 4, supports their relevance to the molecular mechanisms underlying the therapeutic action of the herbs.

Fig. 4. Cluster analysis of the intersecting protein targets within the PPI network.

GO and KEGG Enrichment Analysis

The 56 potential target genes associated with cough and cold were submitted to the David database for Gene Ontology (GO) annotation and KEGG pathway enrichment analysis, with the significance threshold set at P ≤ 0.05. As shown in Fig. 5, Cellular Component (CC) domain enrichment analysis revealed that a significant proportion of genes were localized to membrane-associated and intracellular structures. The most prominently enriched terms included plasma membrane (43 genes), cytosol (39 genes), nucleus (37 genes), cytoplasm (36 genes), and membrane (35 genes), indicating the involvement of these genes in both the structural and functional compartments that are critical for cellular signaling and transport. In the Molecular Function (MF) category, the enriched terms were predominantly associated with binding and enzymatic activities. The top terms were protein binding (45 genes), ATP binding (25 genes), and identical protein binding (20 genes). Additionally, several kinase-related functions, such as protein kinase activity, histone H2AXY142 kinase activity, and protein tyrosine kinase activity, were significantly enriched, reflecting the role of these genes in post-translational modifications and intracellular signaling cascades. Biological Process (BP) analysis identified enrichment in a wide range of regulatory pathways. Notably, signal transduction (30 genes), negative regulation of apoptotic process (26 genes), ephrin receptor signaling pathway (20 genes), and EGFR signaling pathway (19 genes) were among the most significant. Other enriched pathways included insulin receptor signaling, PI3K/AKT pathway activation, and positive regulation of cell proliferation, suggesting that these genes participate in the essential mechanisms underlying cellular growth, communication, and survival. These findings underscore the functional significance of the gene set, particularly in membrane dynamics, molecular interactions, and signal transduction pathways, which may be central to the underlying biological processes of the studied conditions.

Fig. 5. GO annotation chart.

To further define the signaling mechanisms associated with candidate targets, KEGG pathway enrichment analysis was performed. Of the 56 putative targets implicated in cough and cold pathology, 54 genes (96.4%) were significantly enriched in multiple biological pathways (p ≤ 0.05), underscoring their functional relevance. The top-ranked pathways, as illustrated in Fig. 6 and Table IV, revealed substantial gene enrichment in pathways related to cancer (14 genes), the PI3K-Akt signaling pathway (13 genes), focal adhesion (11 genes), proteoglycans in cancer (11 genes), and EGFR tyrosine kinase inhibitor resistance (10 genes). Other significantly enriched pathways included lipid and atherosclerosis, the Rap1 signaling pathway, prostate cancer, endocrine resistance, and platelet activation, each reflecting diverse yet interconnected biological processes. Several key genes were observed to be involved in several of these pathways, including EGFR, PDGFRB, AKT1, PIK3R1, ERBB2, SRC, KDR, PTK2, IGF1R, MET, MMP9, and GSK3B. The convergence of these genes across multiple signaling routes suggests their integral role in mediating inflammatory, immune, and proliferative responses relevant to respiratory tract disorders. These findings support the hypothesis that herbs exert their therapeutic effects via a multi-target, multi-pathway mechanism, primarily through the modulation of signaling cascades associated with inflammation, immune regulation, cellular adhesion, and tissue remodeling.

Fig. 6. KEGG Pathway Analysis.

Table IV. Top 10 Highly Enriched KEGG Pathways

Pathway ID	Pathway name	p-Value	Gene count	Enriched gene IDs
hsa05200	Pathways in cancer	4.82E-07	16	PDGFRB, CAMK2B, GSK3B, MMP1, MMP2, PIK3R1, F2, PTGS2, MMP9, EGFR, IL2, PTK2, IGF1R, ERBB2, AKT1, MET
hsa04151	PI3K-Akt signaling pathway	2.09E-07	14	PDGFRB, GSK3B, SYK, PIK3R1, EGFR, IL2, PIK3CG, PTK2, IGF1R, ERBB2, KDR, AKT1, TEK, MET
hsa04510	Focal adhesion	2.69E-09	13	PDGFRB, GSK3B, SRC, PIK3R1, EGFR, PTK2, IGF1R, MYLK, ERBB2, KDR, AKT1, FYN, MET
hsa05205	Proteoglycans in cancer	3.52E-08	12	CAMK2B, SRC, MMP2, ERBB2, KDR, AKT1, PIK3R1, MET, MMP9, PTK2, EGFR, IGF1R
hsa01521	EGFR tyrosine kinase inhibitor resistance	4.14E-11	11	PDGFRB, GSK3B, SRC, AXL, ERBB2, KDR, AKT1, PIK3R1, MET, EGFR, IGF1R
hsa05417	Lipid and atherosclerosis	6.60E-07	11	CAMK2B, LYN, GSK3B, MMP1, SRC, MMP3, AKT1, PIK3R1, MMP9, PTK2, HSPA1A
hsa04015	Rap1 signaling pathway	5.22E-06	10	PDGFRB, ADORA2A, SRC, KDR, AKT1, TEK, PIK3R1, MET, EGFR, IGF1R
hsa05215	Prostate cancer	1.31E-07	9	PDGFRB, GSK3B, ERBB2, MMP3, AKT1, PIK3R1, MMP9, EGFR, IGF1R
hsa01522	Endocrine resistance	1.42E-07	9	SRC, MMP2, ERBB2, AKT1, PIK3R1, MMP9, PTK2, EGFR, IGF1R
hsa04611	Platelet activation	9.18E-07	9	LYN, SYK, SRC, AKT1, FYN, PIK3R1, F2, PIK3CG, MYLK

Molecular Docking

Based on the top ten hub gene targets and their involvement in key signaling pathways, eight protein targets, AKT1, EGFR, SRC, MMP9, MMP2, PTGS2, KDR, and STAT3 were selected for molecular docking simulations against five phytoconstituents from the herbs: quercetin, ellagic acid, kaempferol, apigenin, and luteolin. The docking results revealed binding energy values ranging from −3.6 to −11.6 kcal/mol, indicating varying degrees of interaction stability between the ligands and target proteins. In general, a lower binding energy reflects greater molecular interaction strength and enhanced stability of the ligand–target complex. A binding energy below −5.0 kcal/mol is typically considered indicative of strong binding affinity and potential therapeutic relevance. All studied complexes exhibited binding energies below this threshold, suggesting favorable interactions. Notably, luteolin, quercetin, and kaempferol displayed the most stable interactions with MMP9 (−9.52, −9.32, and −9.43 kcal/mol, respectively), while luteolin also showed strong binding to MMP2 (−9.66 kcal/mol), as detailed in Table V and Fig. 7. These results support the hypothesis that these compounds significantly contribute to the therapeutic efficacy of herbs through multi-target molecular interactions.

Table V. The Binding Energy (kcal/mol) of Each Ligand Against the Selected Targets

Targets	Compounds	Binding energy (kcal/mol)	X, Y, Z Coordinates
AKT1	Eugalic acid	−5.22	16.03, 24.0, 13.63
AKT1	Quercetin	−5.47	16.03, 24.0, 13.63
EGFR	Eugalic acid	−6.79	17.18, 33.93, 38.42
EGFR	Quercetin	−7.1	17.18, 33.93, 38.42
SRC	Eugalic acid	−7.13	21.03, 20.19, 58.54
SRC	Quercetin	−8.05	21.03, 20.19, 58.54
PTGS2	Apigenin	−8.39	22.60, 40.99, 39.56
MMP9	Quercetin	−9.32	43.97, 15.07, 140.35
MMP9	Kaempferol	−9.43	37.23, 15.77, 145.08
MMP9	Luteolin	−9.66	40.33, 15.77, 145.08
MMP2	Quercetin	−8.37	58.16, 82.00, 146.27
MMP2	Kaempferol	−8.88	51.12, 91.55, 143.23
MMP2	Luteolin	−9.52	58.16, 82.00, 146.27
ERBB2	Eugalic acid	−7.97	17.79, 18.32, 26.32
KDR	Eugalic acid	−6.13	–23.75, –0.93, –15.07
KDR	Quercetin	−8.72	23.75, –0.93, –15.07

Fig. 7. Docking conformations of quercetin, kaempferol, luteolin, apigenin, and eugalic acid with key regulatory proteins (AKT1, EGFR, SRC, PTGS2, MMP9, MMP2, ERBB2, and 3VHE) targets. For each target, the three-dimensional binding pose is shown alongside the two-dimensional interaction profile, highlighting hydrogen bonds, hydrophobic contacts, and other stabilizing interactions that contribute to ligand affinity.

Molecular Simulations

1CK7–Kaempferol Complex

Kaempferol forms multiple persistent interactions with 1CK7. 2D contact analysis showed hydrogen bonding to Ala390 (98% occupancy) and Thr397 (80% occupancy), with additional water-mediated contacts of Pro388, Arg403, and Ile395. Torsion profiles indicated that most rotatable bonds (especially hydroxyl groups) are conformationally restricted, and kaempferol remains in a stable bound conformation. The protein Cα RMSD stabilizes around 4.0–6.5 Å after equilibration (~20 ns), while the ligand RMSD stays between ~3.0–5.0 Å, with no large drift or dissociation as observed in Fig. 8A. RMSF analysis showed that most protein residues fluctuated by ~1–3 Å, with higher peaks (~5–7 Å) in flexible loop regions (~residues 50, 250, and 400). Crucially, binding-site residues (Ala390, Thr397, and Arg403) exhibited lower fluctuations, reflecting local stabilization by kaempferol (Fig. 8B).

Fig. 8. RMSD and RMSF analysis of the 1CK7–Kaempferol complex, showing backbone stability and residue flexibility during 100 ns of simulation.

1CK7–Luteolin Complex

Luteolin engages in hydrogen bonds (via its hydroxyls) with residues SER367, ALA390, LEU391, and THR399 of 1CK7. These interactions occur with high probability (~ 35–75%) over the trajectory. Torsional analysis showed that the rotatable bonds of luteolin adopt largely fixed conformations once bound, with only minor populations in alternate states, indicating a limited conformational drift. The protein Cα RMSD stabilizes around ~5.0–6.0 Å after ~20 ns, while luteolin’s fitted RMSD remains below ~4.0 Å, indicating the ligand remains well-bound during protein fluctuations (Fig. 9A). RMSF analysis revealed low fluctuations (<3 Å) for key binding-site residues (SER367, ALA390, LEU391, and THR399), while more distal loop regions showed higher mobility. This pattern reflects a rigidified binding pocket around luteolin and flexible loops elsewhere (Fig. 9B).

Fig. 9. RMSD and RMSF analysis of the 1CK7–Luteolin complex, highlighting stable binding and localized fluctuations in loop regions.

1CK7–Quercetin Complex

Quercetin established a robust interaction network with 1CK7. The 2D contact map revealed strong hydrogen bonds between ALA390 (95% occupancy) and ARG403 (77%), along with additional contacts with SER405, LEU370, and LEU479. Torsion profiles indicated that quercetin’s rotatable bonds mainly occupied narrow, preferred conformations (sharp histogram peaks), with some flexibility in solvent-exposed angles. This suggests that the ligand adopts its conformation without a large strain. In RMSD analysis, the protein backbone gradually deviated (blue curve rising from ~2 to ~10–12 Å over 100 ns), reflecting conformational adjustment during the run. Quercetin’s fitted RMSD stabilizes around 5–8 Å, indicating that the ligand stays bound as the protein moves, as demonstrated in Fig. 10A. RMSF showed that most residues fluctuated ~2–5 Å, but certain loop regions (residues ~400–500) exceeded 8 Å, highlighting flexible domains. In contrast, regions with persistent secondary structures remained low in fluctuation, as shown in Fig. 10B.

Fig. 10. RMSD and RMSF analysis of the 1CK7–Quercetin complex, indicating consistent conformational stability with minimal deviations.

1GKC–Kaempferol Complex

Kaempferol in 1GKC maintained multiple stable contacts for over 100 ns. Ligand-protein contact analysis showed hydrogen bonds, hydrophobic interactions, ionic contacts, and water bridges, with several residues interacting >30% of the time. Torsion analysis indicated that the aromatic core of kaempferol was rigid, whereas solvent-exposed hydroxyl substituents showed more flexibility. Overall, the torsion density plots revealed minimal strain, which is consistent with a stable protein-bound conformation. The protein backbone RMSD stabilizes at ~1–3 Å, which is consistent with a folded protein under physiological conditions. The fitted ligand RMSD remained similarly low, indicating that kaempferol remained firmly docked without drifting away (Fig. 11A). RMSF profiles showed the highest flexibility at the protein N- and C-termini and in the loop regions, whereas α-helices and β-strands exhibited low fluctuations. In the ligand, the hydroxyl-bearing peripheral atoms fluctuated more, whereas the aromatic rings remained rigid. This pattern indicates that kaempferol’s core interactions are stable and its flexibility is localized, as depicted in Fig. 11B.

Fig. 11. RMSD and RMSF analysis of the 1GKC–Kaempferol complex, revealing moderate residue fluctuations with overall structural stability.

1GKC–Luteolin Complex

Luteolin forms persistent interactions with 1GKC. Throughout the 100 ns, the ligand engages in hydrogen bonds, hydrophobic contacts, ionic interactions, and water-mediated bridges; multiple residues contact luteolin for >30% of the trajectory. Torsional analysis shows that the aromatic core of luteolin remains constrained, while exposed hydroxyl groups sample more conformations. Overall, the ligand backbone remained at low strain, allowing luteolin to adapt within the pocket. The protein RMSD stabilized at approximately 1–3 Å, indicating no large structural drift during binding. Luteolin’s fitted RMSD also remained low, showing that the ligand remained bound without significant displacement, as shown in Fig. 12A. RMSF revealed that the flexible regions of the protein were mainly at the termini and loops, whereas structured helices/strands were rigid. In the ligand, the hydroxyl groups and peripheral atoms fluctuate locally, while the aromatic scaffold is rigid, suggesting strong core interactions with some peripheral flexibility, as observed in Fig. 12B.

Fig. 12. RMSD and RMSF analysis of the 1GKC–Luteolin complex, demonstrating sustained stability of the protein–ligand interaction throughout the trajectory.

1GKC–Quercetin Complex

Quercetin established a stable binding network in 1GKC. It formed hydrogen bonds, hydrophobic contacts, ionic interactions, and water bridges with multiple residues for >30% of the simulation. Torsional analysis indicated that the core scaffold of quercetin was rigid, while solvent-exposed rotatable bonds showed varying flexibility. The torsion density plots suggest that the ligand experiences no significant torsional strain and adopts a stable-bound conformation. The protein RMSD levels off at ~1–3 Å, which is typical of an equilibrated globular protein, confirming no major conformational change during binding. Quercetin’s fitted RMSD aligns closely with the protein backbone, indicating that the ligand remains securely docked, with no major diffusion observed in Fig. 13A. RMSF profiles showed the greatest flexibility at the protein termini and loops, while α-helices/β-strands remained rigid. Ligand RMSF highlights the mobility of the hydroxyl groups, whereas the aromatic ring system remains rigid, reinforcing its role in hydrophobic interactions, as shown in Fig. 13B.

Fig. 13. RMSD and RMSF analysis of the 1GKC–Quercetin complex, showing stable complex formation with flexibility restricted to terminal and loop regions.

Molecular Mechanics Generalized Born Surface Area (MM-GBSA)

The binding free energy of the simulated trajectories was estimated using the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) approach. A more negative ΔG_bind value reflects stronger interactions and greater stability of the protein–ligand complex [24]. The MM-GBSA binding energy of 1CK7–Kaempferol, 1CK7–Luteolin, 1CK7–Quercetin, 1GKC–Kaempferol, 1GKC–Luteolin, and 1GKC–Quercetin was −105.23, -90.45, -91.50, -84.68, −81.49, and −87.53 kcal/mol, respectively. Energy decomposition showed that Van Der Waals and Lipophilic energies contributed the most to the negative binding energy (Table VI).

Table VI. MM-GBSA Energies for Protein-Ligand Complexes

Complexes	ΔG_Bind (kcal/mol)	Coulomb (kcal/mol)	Covalent (kcal/mol)	H bond (kcal/mol)	Lipo (kcal/mol)	Packing (kcal/mol)	Solv_GB (kcal/mol)	vdW (kcal/mol)
1CK7–Kaempferol	−105.23	−26.97	3.07	−1.00	−59.47	0	32.36	−53.20
1CK7–Luteolin	−90.45	19.31	3.73	−0.71	−58.07	0	37.28	−53.35
1CK7–Quercetin	−91.50	−35.81	3.95	−1.00	−59.03	0	42.80	−42.42
1GKC–Kaempferol	−84.68	−17.96	3.25	−1.34	−47.03	0	21.92	−43.51
1GKC–Luteolin	−81.49	−12.56	1.37	−1.79	−50.28	0	30.13	−48.25
1GKC–Quercetin	−87.52	−20.54	1.49	−1.28	−44.37	0	21.84	−44.65

Discussion

The underlying pathogenesis of cough and cold involves mucosal irritation, excessive mucus secretion, nasal congestion, and activation of sensory nerves, which contribute to coughing, sneezing, and sore throat [26]. If untreated or recurrent, these conditions can lead to complications, such as sinusitis, bronchitis, or secondary bacterial infections. Multiple targets and pathways are known in treatment algorithms for cough and cold. Scientific evidence and clinical studies of medicinal plant extracts and their phytochemical constituents have been well established.

Cough and cold involve airway inflammation, oxidative stress, and remodeling; EGFR is upregulated in asthmatic airways and drives goblet cell hyperplasia and airway hyperresponsiveness; and inhibition of EGFR (or upstream Src kinases) in allergic models suppresses airway inflammation, mucus metaplasia, fibrosis, and bronchoconstriction [27]. SRC family kinases transactivate EGFR to trigger downstream PI3K/Akt and NF-κB signaling; therefore, targeting Src/EGFR blunts inflammatory cell influx and hyperreactivity. Similarly, AKT1 (a PI3K effector) mediates cell survival and cytokine production, promoting inflammation and hyperresponsiveness. PTGS2 encodes COX-2, a key enzyme in airway inflammation. Its induction by cytokines drives prostaglandin production and inflammatory signaling in asthma and viral lung injury [28]. COX-2 products (like PGE₂) can also modulate bronchial tone; therefore, PTGS2 is central to cough/cold inflammation.

MMP9/2 degrades the extracellular matrix, enabling inflammatory cell migration and airway remodeling. MMP-9 is consistently elevated in asthma and cough-variant asthma and correlates with disease severity [29]. Inhibition of MMP-9 reduces airway inflammation and hyper-responsiveness. MMP-2 promotes remodeling by inducing MMP-2 release from fibroblasts, driving the fibroblast-to-myofibroblast transition, and general tissue injury elevates MMP-2 to facilitate smooth muscle migration [30]. Thus, MMP9/2 is a logical target for limiting airway structural changes. ErbB2 (HER2) is an EGF-family receptor that is crucial for epithelial repair. In healthy airways, ErbB2 signaling promotes wound closure, and the loss of ErbB2 delays epithelial repair in asthma, facilitating pathogen invasion [31]. Targeting ERBB2 can enhance barrier repair and reduce inflammation. KDR is a VEGF receptor that is expressed in the airway endothelium and epithelium. Paradoxically, VEGF/KDR signaling helps maintain normal airways: disruption of KDR leads to goblet cell metaplasia and mucus overproduction [32]. Conversely, excess VEGF acting through KDR causes angiogenesis and airway remodeling, and transgenic lung VEGF induces Th2 inflammation, vascular and smooth muscle remodeling, and hyperresponsiveness [33]. Thus, KDR is implicated in both inflammation-induced remodeling and goblet cell differentiation.

Curcuma longa modulates multiple genes and signaling pathways, including AKT1, to deliver anti-inflammatory, antitussive, bronchodilator, and anti-remodeling benefits under respiratory conditions. It inhibits the PI3K/AKT pathway (suppressing AKT1 activity) while enhancing PTEN [34], thereby reducing NF-κB activation and oxidative stress in airway cells, which is crucial for controlling inflammation in COPD and pollution-induced injury models [35]–[37]. This cascade effect leads to decreased transcription of genes such as NFKB1/p65, IL-6, IL-8/CXCL8, TNF-α, COX-2 (PTGS2), MCP-1, and MIP-2, effectively lowering mucus secretion, neutrophil/eosinophil recruitment, cough reflex sensitivity, and bronchoconstriction [38]–[40]. It promotes IκBα stabilization via IKKβ inhibition, further restricting NF-κB nuclear translocation [38], and accelerates the degradation of transcription factors such as SP1, downregulating TLR2, which helps diminish innate immune hyperactivation [39]. Additionally, AKT1 suppression helps limit TGF-β-driven fibrosis and airway remodeling [40].

Glycyrrhiza glabra blocks major inflammatory pathways, such as NF-κB and MAPKs (ERK, p38, JNK), suppressing the expression of pro-inflammatory genes, including TNF-α, IL-1β, IL-6, IL-8/CXCL8, COX-2, iNOS, MCP-1, and chemokines CXCL1/2in lung epithelial cells and immune cells [41], [42]. Licorice also modulates the PI3K/AKT/mTOR pathway, reducing p-Akt signaling to alleviate airway inflammation, hyper-responsiveness, oxidative stress, and fibrosis, which are key drivers of bronchitis and chronic lung injury. It downregulates the IL-17/STAT3 axis and TGF-β/Smad2/3 signaling, thereby reducing eosinophil-driven inflammation, IgE production, and airway remodeling in asthma and bronchitis models [43], [44]. Its flavonoid groups acts as a peripheral and central antitussive, effectively suppressing citric acid-induced cough comparable to codeine [45]. Glycyrrhizin also demonstrates antiviral properties, blocking viral replication and inflammatory injury in the respiratory epithelium, including SARS-CoV-2, by interfering with viral entry and dampening HMGB1–TLR4 signaling [46].

Terminalia chebula exerts its respiratory benefits at the gene target level by modulating key inflammatory, immune, and antiviral pathways. Its extracts suppress the pro-inflammatory enzyme genes, such as iNOS and COX-2, in macrophages, thereby reducing nitric oxide and prostaglandin-mediated bronchial inflammation [47]. It inhibits the NF-κB pathway by preventing IκBα phosphorylation and degradation and blocking p65 nuclear translocation, thereby downregulating downstream genes, such as IL-8, MCP-1, and TNF-α in epithelial and immune cells [48]. Terminalia chebula also constrains MAPK signaling (ERK, p38, JNK), further suppressing the transcription of inflammatory genes, such as COX-2, iNOS, and TNF-α, and reducing reactive oxygen species. Additionally, its bioactive compounds inhibit TLR4/MyD88 signaling, attenuating the activation of the NLRP3 inflammasome, which helps control innate immunity in viral or bacterial infections [49]. On the antiviral front, its constituents block viral protease activity, such as SARS-CoV-2 M, and diminish NF-κB/MAPK-driven inflammatory responses during respiratory viral infection [50].

Adhatoda vasica (Vasaka) downregulates iNOS and COX-2 thereby reducing nitric oxide and prostaglandin-mediated bronchial inflammation [51]. Vasaka inhibits TRPA1 expression and the associated endoplasmic reticulum stress, preventing the overexpression of MUC5AC mucin and protecting epithelial integrity [52]. It also suppresses the upstream stress-related genes AGR2 and TRPV3, further curbing mucin production and cellular damage. Vasaka demonstrates broad transcriptomic modulation, attenuating pro-inflammatory signals such as TGF-β, IL-6, IFN-γ, and hypoxia-inducible HIF-1α, which are implicated in inflammation, fibrosis, and acute respiratory distress [53].

Ocimum sanctum suppresses pro-inflammatory cytokine genes, including TNF, IL-6, MIP-1α, and MCP-1, by inhibiting NF-κB nuclear translocation in monocytes/macrophages, thereby reducing inflammation in bronchial tissues [54]. Tulsi also decreased TLR2 and TLR4 gene expression, further attenuating the innate immune NF-κB axis [55]. In alveolar epithelial cells exposed to bacterial infection, it downregulates HIF-1α and NFκB genes while simultaneously upregulating IFN-γ, thereby enhancing pathogen clearance and cell survival [56]. Additionally, Tulsi activates the Nrf2 antioxidant defense pathway, increasing the transcription of genes like HO-1, SOD, and catalase to protect lung tissues from oxidative stress-related damage [57].

GO enrichment analysis of 54 intersecting targets demonstrated that herbs might display anti- immunomodulatory, antiallergic, antitussive, bronchodilator, and virus-associated respiratory inflammation effects. The bioactive compounds in the herb target genes were predominantly localized to the plasma membrane, cytosol, and extracellular space, as indicated by the GO enrichment analysis. These genes are associated with critical molecular functions such as protein and ATP binding, kinase activity, and receptor signaling. Biologically, they regulate signal transduction, apoptosis, and metabolic signaling.

Complementary KEGG pathway analysis revealed that these targets were significantly enriched in the PI3K-Akt signaling pathway, which plays a central role in airway inflammation, mucus hypersecretion, and epithelial repair. Experimental studies in asthma and COPD models have shown that PI3K inhibitors can reduce inflammatory cell influx, cytokine release, and mucus production, supporting their therapeutic potential in respiratory disorders [58], function in the lung to modulate leukocyte migration, bind and regulate inflammatory mediators, and activate innate immune signaling via Toll-like receptors, thus fine-tuning the pulmonary inflammatory response to infection [59]. Other pathways, such as focal adhesion [60], RAPL signaling [61], [62], and EGFR [63], [64] related networks, which are often mentioned in non-respiratory contexts, directly contribute to preserving epithelial barrier integrity and promoting tissue repair in the airway mucosa. Additionally, platelet activation can amplify airway inflammation through mediator release and interactions with immune cells [65]. Collectively, these interconnected pathways underscore regulatory networks that converge on inflammation resolution, immune modulation, and epithelial regeneration, thereby linking cancer-related signaling to respiratory biology. Their enrichment in our analysis suggests the potential mechanisms by which herbal compounds exert anti-cough, anti-inflammatory, and airway restorative effects. Notably, EGFR resistance, Rap1 signaling, and focal adhesion pathways are closely associated with cell survival, inflammatory responses, and immune regulation. Core genes, including EGFR, AKT1, PDGFRB, ERBB2, PIK3R1, MMP2, and MMP9, function as central signaling hubs across multiple pathways, emphasizing their integrative roles. Together, these findings indicate that medicinal herbs may act through the coordinated modulation of membrane-associated and intracellular signaling, ultimately supporting immune homeostasis and epithelial repair in the context of cough and cold management.

Among the eight docked proteins, MMP-9 and MMP-2 were prioritized for molecular dynamics simulations because they demonstrated the most favorable binding energies with luteolin, quercetin, and kaempferol (all below −9 kcal/mol) and are critically involved in bronchial inflammation, airway remodeling, and pulmonary pathologies related to cough, cold, and respiratory infections. Elevated levels of MMP-9 have been observed in the sputum and bronchoalveolar lavage of asthmatic patients and are positively correlated with the severity of airway inflammation, while MMP-9-deficient animal models display significantly reduced inflammatory responses [66], [67]. Similarly, MMP-2 is implicated in allergic airway inflammation and tissue remodeling, and acts as a key mediator in bronchial extracellular matrix degradation and repair mechanisms [68]. These findings provide a compelling rationale for advancing molecular dynamics to examine the stability, flexibility, and structural implications of quercetin, luteolin, and kaempferol binding under dynamic, physiologically relevant conditions.

The MD data indicated that kaempferol, luteolin, and quercetin all form stable, specific complexes with both 1CK7 (MMP-2) and 1GKC (MMP-9) without disrupting the protein structure. The persistent hydrogen-bond networks we observed were consistent with previously reported flavonoid–protease interactions. For example, kaempferol’s high-occupancy bonds with Ala390 and Thr397 mirror known anchoring modes in metalloproteinase sites, and quercetin’s extensive polar contacts reflect its extra hydroxylation. These binding patterns agree with the literature showing that quercetin, luteolin, and kaempferol can inhibit MMPs, and kaempferol has been experimentally shown to reduce MMP-2/MMP-9 activity and expression during cough-related inflammation [69]. The fact that our simulations yielded no ligand dissociation or stable complexes supports the hypothesis that these compounds are effective binders of MMP-2/9. Moreover, the structure–activity relationship for these flavonoids is evident: the number and arrangement of hydroxyl groups likely modulate the binding affinity. Previous studies have noted that the inhibitory effects of flavonoids depend on the hydroxyl pattern [70]. In our results, the extra OH of quercetin afforded additional interactions (e.g., with Arg403 and Ser405), which could enhance the binding enthalpy, although it also correlated with slightly greater protein flexibility in 1CK7. This trade-off is characteristic: more hydroxyls can increase polar contacts but also impose entropic costs.

Stability metrics highlight key features of molecular recognition. Lower ligand RMSF and low protein fluctuations at the binding site are hallmarks of tight binding. Indeed, studies have found that reduced residue mobility correlates with stronger ligand interactions. In our 1CK7 complexes, binding-site residues showed decreased RMSF relative to the loops, signifying localized stabilization. The 1GKC complexes were more rigid, with a protein RMSD of approximately 1–3 Å and compact radii of gyration (Rg). This greater compactness is associated with stronger binding; as noted in similar MD analyses, a smaller Rg typically implies a tightly packed complex with favorable interactions. Thus, the compact and low-RMSF profiles of the 1GKC systems suggest particularly stable binding. In contrast, the higher RMSD and RMSF in some 1CK7 regions (especially with quercetin) may reflect the inherent domain movements of the full-length gelatinase, but these remained within the normal dynamic range without unfolding.

From an entropic perspective, the restricted torsion of quercetin, luteolin, and kaempferol likely enhances affinity. Each ligand maintained a single or narrow range of conformers when bound, implying that they incur less conformational entropy loss than a more flexible molecule would. This behavior is advantageous for the binding free energy. In fact, studies of protein–ligand thermodynamics emphasize that limiting ligand flexibility upon binding can improve affinity (by reducing the entropy penalty) [71]. Apart from the minor rotation of peripheral hydroxyl bonds, the ring systems remained fixed for quercetin, luteolin, and kaempferol. Consequently, they likely paid a smaller entropic cost and instead gained the enthalpy of the hydrogen bonding and van der Waals contacts.

The preserved secondary structure further confirmed that binding was non-disruptive. None of the ligands induced unfolding, and the helices and strands persisted unbroken. This indicates that these compounds act as classical inhibitors (stabilizing the enzyme–inhibitor complex) rather than as denaturing agents. In enzyme inhibition contexts, maintaining the protein fold is necessary for specific inhibition, and our results reflect this. For instance, the ligand-stabilized binding site remains consistently structured, suggesting that binding does not allosterically destabilize the enzyme. This finding is in line with stable RMSD/RMSF behavior, as other MD studies have shown that stable secondary structure evolution correlates with effective ligand binding.

Conclusion

This study comprehensively investigated the molecular mechanisms underlying the therapeutic effects of five herbs (Respiratory wellness formulation – Koff-N-Cold of Greenspace) in managing cough and cold using network pharmacology and molecular docking. Of the 524 phytochemicals screened, 93 passed drug-likeness filters and 13 bioactive compounds were identified to target 56 disease-associated proteins. Key hub genes, such as AKT1, PTGS2, MMP9, MMP2, EGFR, and SRC, were mapped via protein–protein interaction analysis. GO and KEGG enrichment analyses revealed that these targets are involved in critical biological processes and pathways, such as PI3K-Akt signaling, cytokine-mediated inflammation, and airway remodeling. Molecular docking demonstrated strong binding affinities (−5.22 to −9.66 kcal/mol) for quercetin, kaempferol, luteolin, apigenin, and ellagic acid against major inflammatory and remodeling targets, validating their therapeutic potential. Molecular dynamics simulations confirmed that the conformational stability of MMP2/9 with quercetin, kaempferol, and luteolin showed consistent hydrogen bonding, hydrophobic interactions, and preservation of secondary structural elements over extended trajectories. Collectively, these findings underscore the multitarget therapeutic profile of the polyherbal formulation, suggesting that its phytoconstituents act synergistically to regulate pro-inflammatory cytokines, inhibit matrix metalloproteinases, and stabilize bronchial remodeling enzymes. By addressing both symptomatic relief and the molecular basis of airway inflammation, the formulation demonstrated strong potential as an anti-bronchitis, anti-allergic, and antitussive intervention, thereby offering a promising avenue for respiratory health management and warranting further translational validation in clinical settings.

Limitations and Future Work

Our study presents a systematic computational strategy for the identification of phytoconstituents with therapeutic potential in the management of cough and cold. By integrating network pharmacology, molecular docking, and molecular dynamics simulations, we proposed possible mechanisms of action through which the active compounds may modulate inflammatory mediators, airway remodeling proteins, and immune-related pathways. Drug-target networks were constructed using publicly available databases, which may not comprehensively capture all biologically relevant interactions, leading to possible omission of key targets. Furthermore, the predictive accuracy of network-based analyses would benefit from integration with disease-specific transcriptomic and proteomic data derived from respiratory tissues, which remains limited in current repositories. Preclinical studies employing cellular assays and animal models are necessary to evaluate the efficacy, pharmacokinetic behavior, and safety profiles of the identified compounds.

Acknowledgment

The authors gratefully acknowledge Hirehal Greenspace Herbs for their support in facilitating this study.

Conflict of Interest

The authors declare no conflict of interest with respect to the research, authorship, and/or publication of this article.