A Network Pharmacology and Molecular Dynamics Simulation-Based Study of Qing Run Hua Jie Decoction in Interstitial Pneumonia Treatment

Introduction

As a group of diffuse parenchymal diseases characterized by inflammation and fibrosis of the pulmonary parenchyma,^1,2 interstitial pneumonia (IP) has over 200 different pathologies with known and unknown causes.³ Idiopathic pulmonary fibrosis is the most common and severe symptom in IP.³ Identifying the cause is crucial to the eradication of IP, however, only symptomatic treatment is often applied due to unclear causes. Current regular treatments for IP include immunosuppressive drugs, anti-inflammatory and anti-fibrosis, which effectively slows down the hyperimmunity, inflammation and fibrotic process to a certain extent, but can hardly cure IP.^4,5 Therefore, it would be meaningful to develop new treatments for IP patients to improve their quality of life.

The herbs used in traditional Chinese medicine (TCM) exert a holistic regulatory and synergistic effect through a multi-target and multi-pathway intervention strategy. In this regard, the employing of network pharmacology for illustrating the mechanism of drug action gets increasingly popular. As an integrated discipline, network pharmacology includes pharmacology, network physiology, bioinformatics and cheminformatics.⁶ Network pharmacology offers new ideas for TCM and novel drug study, since its research strategy echoes with the holistic view of TCM theory.⁷ Network pharmacology has been commonly used in medical research currently. Chen et al⁸ constructed a network of herbs and components by screening potential targets of TCM for cough variant asthma treatment through network pharmacology and determined the efficacy by seeking relevant pathways. Tang et al⁹ found that 3 substances in Tripterygium wilfordii, Sedum acre and Hypericum patulum can bind to the zinc-binding site of matrix metalloproteinase-9, which exerts a crucial effect on acute skin inflammation. In this way, the side effects of tetracycline can be avoided, while the efficacy can be improved through a high-throughput virtual screening method of network pharmacology, which can be further developed as a new therapy for acute skin inflammation.⁹ These studies confirm that network pharmacology offers new insight into disease treatment and helps to illustrate TCM theory in terms of biological networks.

As a TCM preparation, Qing Run Hua Jie (QRHJ) decoction is mainly composed of 10 herbs including Radix Adenophorae, Polygonati Rhizoma, Scutellariae Radix, Forsythiae Fructus, Centella Asiatica, Radix Paeoniae Rubra, Aurantii Fructus, Fritillariae Thunbergii Bulbus, Glycyrrhiza uralensis, and Houttuyniae Herba. Several studies have reported the therapeutic effect of the herbs in QRHJ decoction on IP treatment; however, the potential molecular mechanism remains unclear.¹⁰ Our study developed a drug-active ingredient–gene interaction network using network pharmacology and random walk with restart (RWR) analysis. We also revealed the causal mechanism of the effect of QRHJ decoction on IP treatment further by molecular dynamics (MD) simulation, which offers a theoretical basis for developing novel drugs for IP treatment.

Materials and Methods

Acquisition of Active Ingredients of QRHJ Decoction and Corresponding Targets

We searched the 10 herbs in QRHJ decoction in the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, https://tcmsp-e.com/). The active ingredients of 10 herbs were selected with oral bioavailability ≥30% and drug likeness ≥0.18. TCMSP was also used to find the information of the corresponding targets for the screened compounds.

Collection of IP Disease-Related Genes

IP disease-related genes (Category = Protein Coding, Relevance score >6.0) were obtained from the GeneCards database (https://www.genecards.org/).

Building of Protein–Protein Interaction (PPI) Network

PPI network analysis was performed on target genes and IP disease-related genes using the STRING database (https://string-db.org/), and interactions with confidence score ≥0.7 were selected to construct PPI network.

RWR Analysis

The PPI network was subjected to RWR¹¹ using the R package dnet¹¹ with the common gene seeds of drug targets and disease-related genes as parameters. The restart probability of 0.85 was adopted, and the Laplacian method was applied to normalize the adjacency matrix to obtain the affinity score of each gene with the seed. We selected the top 50 node genes in the affinity coefficient for subsequent functional analysis and used Cytoscape software to build the drug-active ingredient–gene interaction network.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analyses

GO functional enrichment and KEGG pathway enrichment were performed on the top 50 node genes in affinity coefficient using the R package clusterProfiler¹² (P<0.05).

Drug Similarity Analysis and ADMET Evaluation

The structures of the core compounds were screened to predict whether these compounds obey Lipinski’s rule of five, the rules of Pfizer and Golden triangle.^13,14 The physicochemical and pharmacokinetic properties of core compounds were evaluated using ADMETlab 2.0 (https://admetmesh.scbdd.com/), which systematically assesses the absorption, distribution, metabolism, excretion, and toxicity (ADMET) of compounds.^13,15,16

Molecular Docking

Network topology characteristics of PPI networks were analyzed using Network Analyzer provided by Cytoscape software. Node genes with degree greater than 20 were selected as hub genes, and subsequent analysis was carried out. From the Protein Data Bank (PDB, https://www.rcsb.org/), we selected the appropriate structure and downloaded the pdb files. Pymol was applied to remove redundant structures, ligands, and surrounding water molecules, etc. We used AutoDock Tools software to add hydrogen atoms and convert the files to pdbqt files, applied ProteinsPlus (https://proteins.plus/) to predict the structure pocket, and set the box size and center coordinates. We downloaded the sdf file of the small molecule from the PubChem (https://pubchem.ncbi.nlm.nih.gov/) and converted it to a pdbqt file using AutoDock Tools software. Finally, docking simulations of small molecules were performed on protein pockets using AutoDock vina. 2D plots of the interactions of small molecules with surrounding residues were displayed using Ligplot software.

Molecular Dynamic (MD) Simulation

The MM/PBSA method was used to calculate the binding free energy of receptor-ligand. Based on the scoring of the docking free energy, we selected the protein-small molecule complex with the highest binding free energy. MD simulations were performed using Gromacs (version 2021.3) to study the intrinsic kinetic properties of the small molecule-bound protein. The topology files of proteins were constructed using the CHARMM36 force field and TIP3P water model, while the topology and parameter files for the small molecules were produced using the CGenFF server (https://cgenff.umaryland.edu/). We constructed the topology file of the complex and placed the complex inside a dodecahedral box with a boundary distance of 1.5. Next, we neutralized the charge in the complex and added the corresponding number of ions. After the energy minimization using the 50,000-step STEEP method, we performed two equilibration steps, first the NVT equilibration at 100ps and then the NPT equilibration at 100ps. The temperature of the system was kept at 300 K by using a Berendsen thermostat with the pressure kept at 1 standard atmosphere. The MD simulation was performed for 50 ns, and the system coordinates were saved every 1 ps.

Trajectory Analysis

We used Gromacs software to analyze the trajectories obtained from MD simulation. Root mean square displacement (RMSD) of the structure was calculated using gmx rms to confirm the convergence of the simulation and the stability of the protein. We measured the gyration radius of the structure using gmx gyrate and calculated the root mean square fluctuation (RMSF) of the residue positions using gmx rmsf. We used gmx hbond to calculate the number of hydrogen bonds generated by the binding of small molecules and proteins during the simulation. Principal component analysis (PCA) and intercorrelation analysis were performed on the system using the Bio3D platform in R language.

Cell Culture

Human embryo lung fibroblast (HELF) MRC-5 cells were cultured in DMEM (Gibco, USA) containing 10% fetal bovine serum (FBS) (Sigma, USA) and 1% double antibody (Hyclone, double antibody), and placed in a sterile incubator at 37°C and 5% CO₂. MRC-5 cells were purchased from BeNa Culture Collection (China). TGF-β1 (5 ng/mL) was added to induce MRC-5 cells for 48 h to promote cell fibrosis.

Cell Viability Assay

The effects of luteolin and diosgenin on the viability of MRC-5 cells were determined by CCK-8. Cells were seeded in 96-well plates and treated with different concentrations of luteolin (0, 5, 10, 20, 40 μM) and diosgenin (0, 5, 10, 20 μM) for 24 h. After that, 10 μL CCK-8 solution (Solarbio, China) was added to each well and incubated at 37°C for 2 h. After that, the optical density (OD) value at 450 nm was detected by microplate reader. Luteolin and diosgenin were purchased from Merck (USA).

Cellular Thermal Shift Assay (CETSA)

MRC-5 cells were seeded in a 10 cm culture dish. After 24 h, the cells were collected and placed in a RIPA lysate containing a protease inhibitor. After centrifugation at 12,000 rpm for 15 min, the supernatant was collected. The cells were divided into three groups and treated with DMSO, 10 μM luteolin and 10 μM diosgenin for 1 h at room temperature. The cells were then heated at 50, 55, 60, 65, 70, 75, 80 and 85°C for 3 min and cooled at room temperature for 3 min. The control group was placed at room temperature without heating treatment. After that, the heated lysate was centrifuged at 15,000 rpm for 20 min at 4°C, and the precipitate and soluble fraction were separated in an ice bath. Finally, the expression of EGFR and HIF1A was analyzed by Western blot.

Western Blot

The cells were divided into three groups: cells without TGF-β1 (Control group), cells with TGF-β1 induction and solvent addition (TGF-β1 group), and cells with TGF-β1 induction and different drugs that had no effect on cell viability. Cells were collected after 24 h of culture. RIPA lysis buffer (Beyotime, China) was added, and cell concentration was detected by BCA method (Thermo Fisher, USA). The protein (including the protein obtained by CETSA method) was added to SDS-PAGE gel for separation and converted to PVDF membrane. Next, the membrane was blocked with 5% skim milk and incubated with the primary antibody overnight at 4°C. After that, the membrane was incubated with secondary antibody goat anti-rabbit at room temperature for 1 h and developed for exposure. Antibodies included α-SMA (ab7817, Abcam), Collagen I (ab260043, Abcam), EGFR (ab52894, Abcam), HIF1A (ab179483, Abcam), GAPDH (ab9485, Abcam), Goat Anti-Rabbit IgG H&L (HRP) (ab6721, Abcam). GAPDH was used as an internal reference.

Data Processing

The data were expressed as mean ± SD and analyzed using GraphPad Prism 8.0 software. Paired t-test was used to compare the data between the two groups. P<0.05 was considered significantly different.

Results

Screening of Active Ingredients and Targets of QRHJ Decoction

QRHJ decoction has the efficacy of clearing heat, moistening the lung, resolving dampness and detoxifying, which is effective in relieving the clinical symptoms of IP, promoting the absorption of lung inflammation and improving lung function. A total of 194 active ingredients of the 10 herbs in the QRHJ decoction were collected from the TCMSP database (Table S1). The targets of the 194 active ingredients were then collected, and 153 active ingredients and 127 targets were obtained after removing target genes repeated values (Table S2).

PPI Network Construction and RWR Analysis

A total of 362 IP-related disease genes were acquired from the GeneCards (Table S3). A total of 127 targets of QRHJ decoction and 362 IP disease-related disease genes were mapped to the PPI network through STRING (Figure 1A), producing a total of 438 nodes and 9070 interactions.

Figure 1 PPI network analysis of QRHJ decoction in IP treatment. (A) PPI network of drug target genes and IP disease-related genes. (B) Venn diagram of drug target genes and IP disease-related genes. (C) Interaction networks of active ingredients and genes of 10 herbs in QRHJ decoction. The red hexagonal nodes represent remedies, the light cyan diamond nodes represent herbs, the blue triangular nodes represent active ingredients in herbs, and the green circular nodes represent genes that may respond to drugs.

The IP disease-related genes and the targets of QRHJ decoction active ingredients were intersected, and 20 target genes for drug targeting IP were obtained (Figure 1B). Taking the 20 genes as seeds, RWR was carried out on the PPI network to confirm the affinity coefficients of each node gene with the seeds (Table S4). The top 50 genes in affinity coefficient and their corresponding active ingredients were determined to build the drug-active ingredient–gene interaction network (Figure 1C).

GO and KEGG Enrichment Analyses

GO and KEGG enrichment analyses were performed on the top 50 genes in affinity coefficient. GO analysis showed that these genes were mainly enriched in protease binding, signaling receptor activator activity, integrin binding, receptor ligand activity, RNA polymerase II-specific DNA-binding transcription factor binding, factor binding DNA-binding transcription, p53 binding, ubiquitin protein ligase binding, ubiquitin-like protein ligase binding, growth factor activity, cytokine receptor binding, disordered domain-specific binding, and protein phosphatase binding and other biological functions (Figure 2A). KEGG analysis revealed that these genes were mainly enriched in the signaling pathways including AGE-RAGE signaling pathway in diabetic complications, lipid and atherosclerosis, proteoglycans in cancer, fluid shear stress and atherosclerosis, human cytomegalovirus infection, prostate cancer, HIF-1 signaling pathway, PI3K-Akt signaling pathway, alcoholic liver disease, TNF signaling pathway, non-alcoholic fatty liver disease, and hepatitis C (Figure 2B). In the results of KEGG analysis, the HIF-1 signaling pathway was associated with viral infection and inflammatory vesicles, while the TNF signaling pathway was associated with hypoxia, and these two signaling pathways may play an important role in the mechanism of QRHJ treatment of IP.

Figure 2 GO and KEGG enrichment analyses of top 50 genes in affinity coefficient in drug-active ingredient-target network. (A) Bubble chart of GO enrichment analysis. (B) Bubble chart of KEGG enrichment analysis.

Molecular Docking Validation

Statistical analysis of the topological properties of the PPI network revealed that the genes with the highest core degree in the interaction network of top 50 genes were interleukin-6 (IL-6), epidermal growth factor receptor (EGFR) and hypoxia-inducible factor 1 subunit alpha (HIF1A), with core degrees of 24, 23 and 21, respectively, indicating that these 3 target genes were in the core position (Figure 3). In the drug–gene interaction network (Figure 1C), there were 5 active ingredients (wogonin, oroxylin a, luteolin, quercetin, paeoniflorin) targeting IL-6, 2 active ingredients (luteolin, quercetin) targeting EGFR and 3 active ingredients (diosgenin, baicalein, quercetin) targeting HIF1A. Molecular docking was performed on the target proteins IL-6, EGFR, HIF1A and their corresponding active pharmaceutical ingredients. For IL-6, we chose the structure of PDB ID: 4CNI, C chain. For EGFR, we chose PDB ID: 5HG8, A chain. For HIF1A, we chose the structure predicted by AlphaFold. Pocket prediction was performed for these 3 structures, and the results are demonstrated in Table 1. From the predicted results, appropriate pockets were chosen for molecular docking, and the docking results are demonstrated in Table 2.

Table 1 Pocket Prediction Table for IL6, EGFR and HIF1A Structures

Table 2 Molecular Docking Results of IL-6, EGFR and HIF1A with Active Small Molecules

Figure 3 Statistics of core degree in the interaction network graph, with the horizontal coordinate indicating degrees and the vertical coordinate indicating gene names.

Drug-Likeness Analysis and ADMET Evaluation of the Selected Compounds

Drug-likeness assessment is a key factor in screening compounds with drug-like properties. As shown in Table 3, wogonin, oroxylin a, luteolin, quercetin, paeoniflorin, diosgenin and baicalein conformed to the filtration standards of the Lipinski’s rule of five, Pfizer’s rules and the Golden triangle rule, while diosgenin conforms to the filtration standards of the Lipinski’s rule of five and the Golden triangle rule.

Table 3 Physicochemical Properties of the Core Compounds

The ADMET evaluation of the pharmacokinetic properties of these selected 7 core compounds is described in Table 4. Wogonin, oroxylin a and diosgenin have high Caco-2 permeability than the others. Wogonin, luteolin, quercetin, paeoniflorin and baicalein have high oral bioavailability indicated by F30% than the others. Intestinal absorption analysis (HIA) revealed good absorption of paeoniflorin. With the exception of diosgenin, almost all selected compounds could penetrate the blood–brain barrier (BBB). Wogonin, oroxylin a, luteolin, quercetin, and baicalein showed inhibitory activities on CYP1A2 and CYP2C9 enzymes. In addition, oroxylin a and paeoniflorin showed lower clearance (CL) than other compounds. Luteolin, quercetin and baicalein showed longer half-life (T1/2) than other compounds. All 7 compounds were not hERG blockers and human hepatotoxicity (H-HT). Only paeoniflorin showed possible carcinogenicity toxicity. Except for paeoniflorin, other drugs such as luteolin, quercetin and diosgenin were predicted to have little respiratory toxicity.

Table 4 ADMET Properties of the Core Compound

Results of MD Simulation

Scoring was based on the affinity of molecular docking, with lower negative values representing tighter binding of proteins and small molecules. As shown in Table 2, binding of HIF1A and diosgenin scored the lowest affinity (−9.2 kcal/mol), followed by the binding of EGFR and luteolin (−8.5 kcal/mol), thus these two complexes were chosen for MD simulation. The binding effect of diosgenin and HIF1A pocket is shown in Figure 4A and B. Diosgenin interacts hydrophobic with residues Asp97, Leu175, Thr176, Ser177, Gly179, Pro226, Pro228, Asn232, Ile233, Ile235, Tyr340, Cys856 of HIF1A. The binding effect of luteolin and EGFR pocket is shown in Figure 4C and D. Hydrogen bonds were formed during the interaction between luteolin and residues Ser720, Met793 of EGFR, while hydrophobic interactions occurred between luteolin and residues Leu718, Gly719, Phe723, Val726, Ala743, Cys775, Met790, Gln791, Leu792, Leu844, Cys856.

Figure 4 Molecular docking simulation analysis results. (A) 3D plot of the molecular model of diosgenin in the key binding site of HIF1A pocket. (B) 2D diagram of the molecular model of molecular docking of HIF1A binding to diosgenin. (C) 3D plot of the molecular model of luteolin in the key binding site of EGFR pocket. (D) 2D diagram of the molecular model of molecular docking of EGFR binding to luteolin.

RMSD was used to measure the stability of the protein, the convergence of the simulation, the deviation of the atomic position from the starting position and the conformational stability of protein.^17,18 As shown in Figure 5A and B, the RMSD analysis of protein atoms suggested that HIF1A-diosgenin complex and luteolin-EGFR complex fluctuated in the early stage and then gradually stabilized, indicating that the protein conformation did not change markedly after the binding of luteolin and EGFR, and their binding was relatively stable. Since EGFR was composed of 1186 amino acids, the RMSD curve showed greater fluctuations than the HIF1A complex.

Figure 5 Results of MD simulation analysis. (A) RMSD diagram of protein HIF1A backbone atoms during MD simulations. (B) RMSD diagram of protein EGFR backbone atoms during MD simulations. (C) Gyration radius diagram of the protein HIF1A during MD simulation. (D) Gyration radius diagram of the protein EGFR during MD simulation. (E) RMSF diagram of protein HIF1A residues during MD simulation. (F) RMSF diagram of protein EGFR residues during MD simulation. (G) Number of hydrogen bonds in small molecule diosgenin binding HIF1A protein during MD simulation. (H) Number of hydrogen bonds in small molecule luteolin binding protein EGFR during MD simulation.

Gyration radius indicated the volume and shape of protein molecules and the compactness of the overall protein structure. Larger gyration radius of the same system represented a bulkier system.¹⁷ As Figure 5C and D showed, the system of HIF1A-diosgenin complex and luteolin-EGFR complex was still in a less stable state at the beginning. However, with the progress of the simulation process, the gyration radius gradually decreased, which suggested that the protein conformation gradually tightened and stabilized. This result was consistent with RMSD, indicating that the binding of small molecules did not affect the stability of protein.

RMSF was used to characterize the change in protein amino acid residues compared to the reference conformation during the simulation, indicating the flexibility (degree of freedom) of the atoms.^17,18 As shown in Figure 5E, RMSF analysis of protein residues showed that some regions of HIF1A, such as residues 517–526, 612–619, had greater residue flexibility. As shown in Figure 5F, the RMSF analysis of protein residues demonstrated that the intermediate regions of EGFR, such as regions of residues 746–758, 867–876, and 917–924, had greater residue flexibility. Structurally, these areas are loop and therefore were more flexible. Figure 5G and H show the number of hydrogen bonds of small-molecule binding proteins in the simulation process of 50 ns, whose results showed that the number of hydrogen bonds of HIF1A-diosgenin was 1 for the most part of the time, and the number of EGFR-luteolin hydrogen bonds was 2 or 3 for most of the time.¹⁸

These results suggested that diosgenin could stably bind HIF1A, and luteolin could stably bind to targeted protein EGFR, which might be the potential mechanism of QRHJ decoction Qingreshen Decoction to exert its pharmacological effects.

Trajectory Analysis

PCA can examine the relationships among different sampled conformations in the trajectories of MD simulation.¹⁹ PCA deeply revealed the conformational differences, and the resulting principal components (orthogonal eigenvectors) indicated the variance of the structural distribution. The percentage of the total mean square displacement (or variance) of the atomic position fluctuations obtained in each dimension is featured by their eigenvalues. Three to five dimensions were usually sufficient to capture more than 70% of the variance in a given MD trajectory. Thus, a few principal components were enough to offer a valid description, while retaining most of the variance in the original distribution. As shown in Figure 6A, HIF1A-diosgenin complex has two different conformation groups along PC1 (42.93%). As shown in Figure 6B, the EGFR-luteolin complex was grouped in different conformations along PC1 (33.62%), one centered at –30 and the other at +10. The PC1-3 motion trajectories of HIF1A-diosgenin are presented in Figure 6C, and the PC1-3 motion path of EGFR-luteolin is shown in Figure 6D, in which the loop was mainly moved. The dynamic correlation matrix represents the correlation between specific atoms of each amino acid in the protein. Figure 7A demonstrates the cross-correlation between the HIF1A protein residues during the simulation process, and Figure 7B shows the cross-correlation between the EGFR protein residues during the simulation process, with positive values (blue) indicating that the residues move in the same direction during the simulation process and negative values (red) indicating that they move in the opposite direction. A darker blue color indicated a stronger correlation.

Figure 6 PCA diagram based on trajectories. (A) Trajectory based principal component analysis diagram of HIF1A-diosgenin complex. (B) Trajectory based principal component analysis diagram of EGFR-luteolin complex. (C) Motion trajectory diagram of HIF1A-diosgenin complex ① PC1, ② PC2, ③ PC3. (D) Motion trajectory diagram of EGFR-luteolin complex ① PC1, ② PC2, ③ PC3.

Figure 7 Residue intercorrelation diagram. (A) HIF1A residue cross-correlation diagram. (B) EGFR residue cross-correlation diagram. The blue part indicates synergy, and the red part indicates relative motion. The horizontal and vertical coordinates represent the amino acid sequence label.

The Effects of QRHJ Active Ingredients (Luteolin and Diosgenin) on HELF Cell Fibrosis

Based on the above studies, we used the CETSA method to verify the stability of the target protein induced by the compound. The results showed that the protein expression of EGFR and HIF1A decreased gradually with the increase in temperature. Under the intervention of luteolin and diosgenin, the thermal stability of EGFR and HIF1A proteins was higher than that of the control group (Figure 8A). These studies suggested that EGFR and HIF1A may be the direct targets of luteolin and diosgenin. Pathologically, IP is a common representative of idiopathic pulmonary fibrosis, and the final feature of IP is pulmonary interstitial fibrosis.²⁰ Therefore, referring to the previous research methods,²¹ the inhibitory effect of drugs on TGF-β1-induced HELF cell fibrosis was explored through HELF cells, which could reflect the mechanism of drugs in IP from a specific perspective. Firstly, we analyzed the effect of changes in luteolin and diosgenin concentration on cell viability by CCK-8 method. The results showed that when the concentration of luteolin and diosgenin reached 20 μM, the activity of MCR-5 cells was affected (Figure 8B). Therefore, we chose the drug concentration (10 μM) that was inactive to the cells to treat the cells and detected the effect on the fibrosis of MCR-5 cells. After TGF-β1 treatment, the expression of α-SMA and Collagen I increased, indicating that TGF-β1 successfully induced fibrosis in MCR-5 cells. After treatment with luteolin and diosgenin, the protein expression decreased. It indicated that luteolin and diosgenin could inhibit TGF-β1-induced MCR-5 cell fibrosis (Figure 8C). Therefore, we speculated that one of the mechanisms by which QRHJ exerted therapeutic effects in IP may be to inhibit the progression of IP by inhibiting cell fibrosis.

Figure 8 The effect of QRHJ effective active ingredients on HELF cell fibrosis. (A) CETSA was performed on MCR-5 cells treated with luteolin, diosgenin or DMSO; (B) CCK-8 was used to detect the effect of different concentrations of luteolin and diosgenin on the activity of HELF cells; (C) Western blot was used to determine the expression of fibrosis-related proteins induced by TGF-β1 after luteolin and diosgenin treatment. *P<0.05.

Discussion

QRHJ decoction moistens the lung and relieves dryness through Radix Adenophorae, Polygonati Rhizoma and Fritillariae Thunbergii Bulbus, and also clears away heat and detoxifies by Scutellariae Radix, Centella Asiatica, Forsythiae Fructus and Houttuyniae Herba. In the decoction, Aurantii Fructus can boost Qi, and Houttuyniae Herba activates blood, which can smooth Qi and relieve damp toxins. Scutellariae Radix is bitter and dry, and Houttuynia cordata is clear and beneficial, which can facilitate the elimination of damp toxins.¹⁰ QRHJ decoction is effective for IP with symptoms of fever with dampness; however, its mechanism remains unclear. Our study used network pharmacology and MD methods to explore the potential mechanism of QRHJ decoction improving clinical symptoms of IP and to offer a theoretical basis for the research of novel drugs for IP treatment.

We found that IL-6, EGFR and HIF1A, the active ingredients in QRHJ decoction, were at the core position in the top 50 genes interaction network through network pharmacology analysis. IL-6 is a cytokine that exerts effects on inflammation, B cell maturation and inflammation-associated cancers, and induces the expression of several proteins correlated to acute inflammation. IL-6 can stimulate B cells to produce antibody and has various biological functions in metabolism, tissue regeneration and immunity.^22–24 EGFR is an ErbB of receptor tyrosine kinases that converts extracellular signals (eg, growth factors and cytokines) into cellular responses and involves in non-small cell lung cancer progression.²⁵ HIF1A is a major regulator of cellular oxygen homeostasis and may affect pulmonary fibrosis and inflammation.^26,27 Therefore, IL-6, EGFR, and HIF1A, which are directly targeted by the active ingredients in QRHJ decoction, may be the key targets for IP treatment.

The binding relationships between the target proteins IL-6, EGFR, HIF1A and the corresponding drug active ingredients were validated through molecular docking in our study. Cell experiments further showed that HIF1A and diosgenin, EGFR and luteolin were stably bound. Diosgenin and luteolin could inhibit the promoting effect of TGF-β1 on HELF cell fibrosis. Then, HIF1A-diosgenin complex and EGFR-luteolin complex were selected for MD simulation for their lowest binding affinity score. The results of MD simulation demonstrated that residues of HIF1A and diosgenin formed hydrophobic interactions, and the residues of luteolin and EGFR formed hydrogen bond interactions and hydrophobic interactions. After HIF1A was bound to diosgenin and luteolin to EGFR, the conformation of the protein did not change notably, and became gradually compact and stable. Molecular dynamics simulation is mainly used to analyze the conformational changes and stability of ligand-protein complexes after docking. RMSD represents the relatively stable binding of proteins and compounds, indicating the position difference of the entire structure over time.²⁸ RMSF is an indicator of the flexibility of intra-protein residues.²⁹ The radius of gyration can obtain the flexibility and scalability of the protein. The structural stability of the protein is related to hydrogen bond integrity and stability.³⁰ From the overall MD simulation (RMSD, RMSF, radius of gyration, number of hydrogen bonds), we concluded that diosgenin and luteolin were relatively stable complexes, showing good binding affinity. The combination of drugs and proteins did not affect the protein conformation. Based on MD molecular simulation, trajectory analysis was used to further analyze the specific movement mode within the molecule and the correlation between amino acid residues, contributing to the understanding and development of allosteric small-molecule inhibitors.

Diosgenin is a steroid saponin with neuroprotective,³¹ anticancer,³² atherosclerosis,³³ hypolipemia³⁴ and treatment of diabetes.³⁵ Du et al³⁶ found that diosgenin alleviates lung inflammation and fibrosis in silicosis by promoting autophagy in alveolar macrophages, which is similar to our findings that diosgenin has a potential role in treating lung inflammation and fibrotic damage. Wu et al³⁷ found that diosgenin has the potential to regulate aerobic glycolysis and M1 polarization of macrophages by inhibiting mammalian target rapamycin complex 1 (mTORC1)/HIF1A. Luteolin is also named as lignans, and EGFR is one of the 9 core targets of luteolin in anti-inflammation.³⁸ Studies have reported that luteolin exerts a crucial effect on the inhibition of mutant EGFR and is a potential agent in cancer therapy.³⁹ The anti-inflammatory effects of luteolin are also significant. Wang et al⁴⁰ found that luteolin can hinder the expression of M1 type surface markers on macrophages and proinflammatory mediators, promote M2 type polarization, and influence the expression of p-STAT3 and p-STAT6, thereby exerting an anti-inflammatory effect. Luteolin can also enhance the autophagy in mice with LPS-mediated myocardial injury.⁴¹ Lee et al⁴² discovered that luteolin can reduce the expression of EGFR and repress EGF-mediated EGFR autophosphorylation in ER-negative breast cancer cells, thus inhibiting the EGFR autophosphorylation-induced MAPK and PI3K/Akt signaling pathways, which can promote the malignant progression of cancer. This is similar to our results that luteolin could inhibit LPS-induced inflammatory injury in the lungs. These indicated that QRHJ decoction could inhibit the expression of inflammatory cytokines in LPS-induced inflammatory damage cells by targeting HIF1A via the active ingredient diosgenin, and targeting EGFR via the active ingredient luteolin, which may be the potential mechanisms of QRHJ decoction treating IP. However, this still needs to be explored in subsequent experiments.

We revealed the key compounds, target genes and related pathways of the active ingredients in the 10 herbs of QRHJ decoction in IP treatment through network pharmacology. We also verified the relationship between key compounds and target genes through molecular docking, MD simulation and trajectory analysis, revealing the potential therapeutic mechanism of QRHJ decoction in IP treatment. However, our study demonstrated the mechanism of QRHJ decoction for IP treatment only at the molecular level, which had some limitations and needed further clinical experiments for validation.