Network Pharmacology-Based Strategy Combined with Molecular Docking and in vitro Validation Study to Explore the Underlying Mechanism of Huo Luo Xiao Ling Dan in Treating Atherosclerosis

Introduction

Atherosclerosis (AS), a lipid-driven chronic inflammatory vascular disease characterized by progressive thickening of the blood vessel wall, is one of the main causes of death from cardiovascular disease. The pathogenesis of AS involved a variety of possible mechanisms, including lipid infiltrations, injury responses, inflammatory responses, immunological responses, and hemodynamic effects.^1–3 AS can be induced by long-term elevated blood lipids and cause a variety of cardiovascular-related complications, such as thrombosis, ischemia, and angina pectoris.⁴ Therefore, effective prevention and treatment against AS are of great significance in the prevention of certain cardiovascular diseases.⁵ Most individuals with AS need to be treated with lipid-lowering drugs (such as statins).⁴ However, statins have many inevitable side effects, for example, liver/kidney dysfunction and rhabdomyolysis.⁶ Moreover, new options now exist for the prevention of AS that are not optimized for statin therapy.⁷

Traditional Chinese Medicine (TCM) is a unique health resource in China and has been applied to human beings for more than two thousand years. Recent studies show that TCM has the potential to prevent and treat a variety of diseases, such as COVID-19,⁸ cancer,⁹ ischaemic stroke,¹⁰ primary dysmenorrhea,¹¹ and mild active ulcerative colitis.¹² In the treatment of AS, TCM shows unique advantages, including the overall regulatory function of the body and the little toxic and side effects.^13–15

Huo Luo Xiao Ling Dan (HLXLD) is a classic formula in the medical records of Yi-Xue-Zhong-Zhong-Can-Xi-Lu by Zhang Xichun,¹⁶ consisting of Angelica sinensis, Salvia miltiorrhiza, Commiphora myrrha, and Boswellia sacra (the detail information about the formula of HLXLD was listed in Table 1).¹⁷ In recent years, increasing evidence suggests that HLXLD may be effective in the treatment of AS. HLXLD up-regulates the protective factor TGF-β to inhibit the immune-inflammatory response in ApoE knock-out mice with early AS.¹⁸ HLXLD treatment can decrease the values of TC, TG, HDL, TC/HDL, and TG/HDL in serum, indicating the potential to treat AS.¹⁹ Moreover, studies have shown that the active ingredients derived from HLXLD, such as quercetin,²⁰ luteolin,²¹ dihydrotanshinone I,²² β-sitosterol²³ and tanshinone IIA,²⁴ exhibit potent pharmacological activities against AS based on individual single compounds of HLXLD. Like other TCM, HLXLD has a complex mechanism involving many complicated relationships between components and targets, targets and diseases, targets and pathways in the treatment of certain diseases.

Table 1 Detailed Information About the Formula of HLXLD

A modern research strategy combining network science, computer science, mathematics, and network pharmacology has been widely used in TCM research. This research strategy can help to clarify the action mechanisms and pathways of TCM from the aspects of molecular level or molecular network regulation, and screen out active ingredients or disease targets from a large number of data.^25,26 Additionally, molecular docking technology is a theoretical simulation method for studying the intermolecular interactions (such as small-molecule ligands and biomacromolecules), and predicting their binding modes and affinity.^27,28 The combination of network pharmacology prediction and molecular docking-based strategy has been applied in a large number of studies to discover the potential pharmacological mechanism of TCM. For instance, Liu et al identified 95 active compounds and 265 predicted targets in Bu Yang Huan Wu decoction by network pharmacology analysis and revealed that the mechanism of Bu Yang Huan Wu decoction in treating AS is related to inflammation and apoptosis pathways. In the early stage, we have conducted a network pharmacology approach to identify the active ingredients of Huai Hua San and explore the mechanisms underlying the treatment of ulcerative colitis.²⁹

In the present study, as shown in Figure 1, we adopted the method of network pharmacology to find the active compounds and anti-AS targets of HLXLD by constructing a multi-dimensional network diagram. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis combined with molecular docking and cell verification experiments were carried out to investigate the complex effects of botanical drugs on the prescriptions. These findings may provide a scientific basis and ideas for follow-up research.

Figure 1 Schematic representation of the proposed mechanism in HLXLD against AS. Network pharmacology was performed to identify the active compounds and key targets of HLXLD against AS. Molecular docking and molecular dynamics simulation were conducted to predict the interaction and connection stability between active components and candidate targets. Cell experiments confirmed that lipid accumulation, inflammation and MAPK signaling pathway were alleviated by candidate compound dihydrotanshinone I in ox-LDL-induced RAW264.7 cells.

Materials and Methods

Data Preparation

Active Ingredients of HLXLD

The compounds of HLXLD were searched from four public databases including Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, https://tcmsp-e.com/), The Encyclopedia of Traditional Chinese Medicine (ETCM, http://www.tcmip.cn/ETCM/index.php/Home/Index/index.html), Bioinformation Analysis Tool for Molecular mechanism of Traditional Chinese Medicine (BATMAN-TCM, http://bionet.ncpsb.org.cn/batman-tcm/) and Traditional Chinese Medicine Database@Taiwan (TCM, Database@Taiwan, http://tcm.cmu.edu.tw/zh-tw/review.php?menuid=3). The OB (oral bioavailability) and DL (drug-likeness) values of compounds were obtained from TCMSP, and only compounds meeting the screening criteria of OB >30% and DL ≥ 0.18 were retained for further research.

Targets Prediction of HLXLD and as

The study was performed in accordance with good clinical practice standards and the tenets of the Declaration of Helsinki 2002, and approved by the clinical research ethics committee in the Second Xiangya Hospital of Central South University (No. 2021–248). The approval certificate was provided in the supporting information. Targets prediction of HLXLD obtained from three public databases, including TCMSP, BATMAN-TCM (score cutoff was set to 20, and adjusted P-value was set to 0.05) and the webtool of SwissTargetPrediction (SwissTargetPrediction, http://www.swisstargetprediction.ch/, Prob value >0). We selected “Homo sapiens” as a species. The chemical information was obtained from the PubChem database (PubChem, https://pubchem.ncbi.nlm.nih.gov/). The protein target names were converted to their corresponding official gene symbols by the UniProt database (https://www.uniprot.org/). Besides, targets of AS were explored by the keyword of “atherosclerosis” from five public databases, including GeneCards: The human gene database (https://www.genecards.org/, Ver. 5.7, Relevance score ≥1), Online Mendelian Inheritance in Man (OMIM, https://omim.org/, Updated September 4th, 2021), PharmGKB (https://www.pharmgkb.org/), Therapeutic Target Database (TTD, http://db.idrblab.net/ttd/, Updated September 29th, 2021), and DrugBank (https://go.drugbank.com/). Only “Homo sapiens” proteins linked to atherosclerosis were selected. Venn diagram of disease targets in five databases were drawn using R software. The Venn diagram of common targets of HLXLD and AS were drawn on the website of http://bioinformatics.psb.ugent.be/webtools/Venn/.

Bioinformatics Analysis

Protein–Protein Interaction (PPI) Network Construction

PPI network of common targets of HLXLD and AS was constructed using the functional protein association networks database (STRING, https://cn.string-db.org/), with species limited to “Homo sapiens”. The minimum required interaction score was set to 0.98, and the independent target protein nodes were hidden. To obtain the potential core target network, the PPI network was visualized by Cytoscape 3.8.2.

Core targets were screened out by calculating six topological parameters of nodes, including betweenness centrality (BC), closeness centrality (CC), degree centrality (DC), eigenvector centrality (EC), network centrality (NC), and local average connectivity (LAC), through CytoNCA plugin in Cytoscape software, which can be utilized to quantitatively assess the calculation results, evaluate the accuracy by statistical measures and integrate biological data with topological data to detect specific nodes. After calculating six parameters using the CytoNCA plugin, only nodes with all six parameters higher than the corresponding median values were retained to further construct the core PPI network.

Enrichment Analysis

After converting the common official gene symbols of HLXLD and AS to Entrez ID by R 4.1.1 software, GO and KEGG pathway enrichment analyses were carried out using R 4.1.1 software with related R packages (colorspace, stringi, DOSE, clusterProfiler, ggplot2, enrichplot, pathview, BiocManager, and org.Hs.eg.db). Q values <0.05 were considered statistically significant and retained.

Network Construction

To illustrate the relationship between botanical drugs and compounds, diseases and targets, as well as targets and signal pathways, the four botanical drugs, candidate active compounds, common target genes, disease name and top 20 pathways were introduced into Cytoscape 3.8.2 to construct a network of Compounds-Disease-Targets (C-D-T) and Compounds-Pathway-Targets (C-P-T) network.

Molecular Docking

The compounds with the highest degree in the C-D-T network and target proteins with the highest degree in the core network were chosen to conduct molecular docking with AutoDock 1.5.6 Vina. The process of the molecular docking experiment was as follows:

First, receptor and ligand preparation of 3D structure: The 3D structural receptors were downloaded from Protein Data Bank (PDB, https://www.rcsb.org/) meeting the standard as far as possible: Organism(s) from “Homo sapiens”, no mutation(s), with higher resolution and containing small-molecule ligand(s). The 2D structural ligands were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), and then imported into Chem3D software to obtain the 3D structure of ligands. The receptor was saved in pdb format, and the ligand was saved in mol2 format.

Second, receptor and ligand preparation of pdbqt format: The macromolecule receptor was opened in Pymol software to remove water molecules, co-crystallized ligands, and ions. Subsequently, the macromolecule receptor was inputted into AutoDockTools 1.5.6 to add hydrogens and Kollman partial charges, and saved in pdbqt format. Then, the ligand was opened in AutoDockTools 1.5.6 and saved in pdbqt format.

Third, prediction of mating pockets of macromolecule receptor: The pdbqt format of the receptor was imported to AutoDockTools 1.5.6 to construct mating pockets of docking. Using the self-carrying molecules surrounding the protein crystal structure as the geometric center, a box with the largest range was formed. Spacing (angstrom) was set to 1, energy_range was set to 5, num_modes was set to 20, the numbers of points in the mating box (centre_x, centre_y, centre_z), and the size of the mating box (size_x, size_y, size_z) were set on the macromolecule and set to keep the protein completely covered by mating box, respectively. These parameters were saved in a configuration (config) file of each protein mating box for further docking by Vina.

Docking and visualization: Molecular docking and calculation of docking affinity were operated with Vina software. During the process, pdbqt files of protein and ligand were also used in addition to the config file containing the information of the mating box. The program performed a conformational search for the ligand in the box range and finally scored according to their conformation, orientation, position and energy of the ligand by Vina. After Vina finished running, the binding affinity score between the macromolecule receptor and small-molecule ligand was output as log.txt files. The best combination mode of receptor and ligand was presented at the first line in log.txt with the lowest docking binding free energy. The interaction between receptor and ligand was visualized and hydrogen bonds were displayed by Pymol (3D) and LigPlus (2D).

Molecular Dynamics Simulation

After docking, the molecular dynamics (MD) simulation strategy was carried out to simulate the binding stability of the receptor and ligand. The entire system added AmberFF99SB and gaff force field parameters to proteins and molecules to set the force fields of carbon atoms and oxygen atoms with different bond angles. Subsequently, with the protein as the center, a 1-nm cubic water box and Na⁺ were added to make the system stable. In the end, MD simulations were performed.

The MD simulation process was as follows: 1) Two-step energy minimization. We firstly limited the protein, minimized the energy of water molecules, and then released the protein, thereby minimizing the energy of the entire system. In the first energy minimization, a total of 5000 cycles were used, and the steepest descent method was used for 1500 cycles. The second energy was minimized for a total of 5000 cycles, and the steepest descent method was firstly used to make 2000 cycles. 2) System balance. The system is heated from 0°C to 100°C, and the Langevin temperature control method was used to balance 100 ps. Then, we proceeded to the boost balance process, during which a total of 100 ps were balanced, and the isotropic Berendsen voltage control method was used. 3) Dynamic simulation. In this dynamic process, a constant temperature of 100°C was used to perform a molecular dynamics simulation of the overall system. The temperature and voltage control methods were the same as in the previous stage. The cutoff distance between van der Waals energy and short-range electrostatic energy was 10 Å, and the Particle-Mesh-Ewald (PME) method was used to calculate the long-range electrostatic energy.

In the process of molecular dynamics, the human body temperature was used for the system. The molecular dynamics time was controlled at 100 ns, and the RMSD value of the complex system was monitored in real time until the project goal was reached. Subsequently, the average structure of the protein-molecule meeting the expectations was taken out for subsequent binding mode analysis.

In vitro Validation

Cell Culture

The mouse mononuclear macrophage leukemia cell line (RAW264.7) was purchased from Beijing Beina Chuanglian Institute of Biotechnology (Beijing, China) and cultured using Dulbecco’s Minimum Essential Medium (DMEM) medium containing 10% fetal bovine serum (FBS), 100μg/mL streptomycin, and 100U/mL penicillin. Cells were incubated at 37°C in a humidified atmosphere of 5% CO₂ and 95% air.

Oil Red O Staining

RAW264.7 cells were cultured in 6-well culture plates with sterile coverslips to assess foam cell formation. The cells were pretreated with dihydrotanshinone I with a concentration of 10, 100, 1000 nM for one hour, then treated with 50 μg/mL ox-LDL for 24 h. After treatment, the cells were washed three times with PBS, fixed with 50% isopropanol for one minute, and stained with oil red O staining solution for 10 min, followed by discarding the dye solution and washing it three times. Foam cells were photographed under a microscope and analyzed using an image analysis system.

Cytokine Enzyme-Linked Immunosorbent Assay (ELISA)

RAW264.7 cells were cultured in 6-well culture plates and pretreated with dihydrotanshinone I with a concentration of 10, 100, 1000 nM for one hour, then treated with 50 μg/mL ox-LDL for 24 h. The levels of TNF-α, IL-1β, and MCP-1 in cell supernatants were determined by an ELISA assay kit meeting the manufacturer’s protocol. Quantitative determinations were performed in three independent experiments.

Western Blot Analysis

Total cellular protein was collected by using a high-efficiency radio immunoprecipitation assay (RIPA) buffer with 0.1% PMSF (BOSTER Biological Technology; Wuhan, China) and the total protein concentration of each sample was determined using a BCA kit. The protein samples were separated by SDS-PAGE using 8–10% gradient gel. The separated protein on the gel was transferred onto 0.45 μm PVDF membranes. The expression of the target protein was normalized to GAPDH protein. Antibodies against p-STAT3 (#9145, Cell Signaling Technology, Danvers, MA, USA), STAT3 (#4904), p38 (#8690), and p-p38 (#4511) were purchased from Cell Signaling Technology (Danvers, MA, USA), and antibodies against ERK (ab17942), p-ERK (ab278538), JNK (ab199380) and p-JNK (ab124956) were purchased from Abcam (Cambridge, UK).

Statistical Analysis

The results were presented as means ± S.E.M and data comparison of multiple groups was adopted using one-way ANOVA and Student-Newman-Keuls test. The results were considered statistically significant when the P value was less than 0.05.

Results

Potential Active Components in HLXLD

HLXLD is composed of Angelica sinensis, Salvia miltiorrhiza, Commiphora myrrha, and Boswellia sacra. After ADME screening (OB >30% and DL ≥0.18), a total of 108 potentially active compounds were retrieved from four databases (Table 2). Among them, six are from Angelica sinensis, 66 from Salvia miltiorrhiza, 29 from Commiphora myrrha, four from Boswellia sacra, and three (beta-sitosterol, stigmasterol, isoimperatorin) from two or more botanical drugs.

Table 2 108 Bioactive Compounds of HLXLD

Common Targets of HLXLD and AS

After the initial search from TCSMP, BATMAN-TCM, and SwissTargetPrediction, 6014 targets were predicted from 108 bioactive compounds in HLXLD. Subsequently, we removed duplicate values and converted protein names into gene symbols, thereby obtaining 961 targets (Supplement Table 1). Besides, 1463 AS-related targets were searched from five databases, including 1387 in GeneCards (Relevance score ≥1), 25 in PharmGkb, 24 in TTD, 26 in DrugBank, and one in OMIM. After we removed duplicate values, 1436 targets were obtained (Figure 2A, Supplement Table 1). Then, the 321 common targets of HLXLD and AS were achieved by R software for subsequent study (Figure 2B, Supplement Table 1).

Figure 2 Venn diagrams. (A) AS-related targets came from five databases. (B) The intersection genes of identified AS-related targets and targets of HLXLD.

C-D-T Network Construction and Topological Network Analysis

To obtain a further understanding of the relationship between the 108 compounds and 321 intersection target genes on a system level, the C-D-T network with 429 nodes and 1754 edges was constructed by Cytoscape 3.8.2 (Figure 3). The turmeric V nodes represent four botanical drugs from HLXLD. The circular nodes of red, blue, green, and yellow stand for compounds from Angelica sinensis, Salvia miltiorrhiza, Commiphora myrrha and Boswellia sacra, respectively. The circular purple nodes stand for AS target genes. Larger node sizes indicate higher degree values. The top six compound nodes with the largest size and highest degree were quercetin (MOL000098, degree = 116), dihydrotanshinone I (MOL007101, degree = 76), pelargonidin (MOL007101, degree = 61), luteolin (MOL000006, degree = 55), guggulsterone (MOL001175, degree = 55) and β-sitosterol (MOL000358, degree = 53) (Table 2, Figure 3). The top eight gene nodes were PTGS2 (degree = 74), PTGS1 (degree = 43), RXRA (degree = 41), MCOA2 (degree = 38), SCM5A (degree = 37), PTPN1 (degree = 37), HSP90AA1 (degree = 36) and ESR1 (degree = 34) (Figure 3). These top compounds and genes may be the critical nodes in the network and possess an important anti-AS effect.

Figure 3 C-D-T network of HLXLD in the treatment of AS. It contains 429 nodes and 1754 edges. The turmeric V nodes represent four botanical drugs from HLXLD. The circular nodes of red, blue, green and yellow stand for compounds from Angelica sinensis, Salvia miltiorrhiza, Commiphora myrrha and Boswellia sacra, respectively. The circular purple nodes stand for AS-related target genes. Larger node sizes indicate higher degree values.

PPI Network Construction and Analysis

To obtain the PPI network from the STRING database, “Homo sapiens” was selected in the organism column, and the minimum required interaction score was set to 0.98 (Figure 4A). It contains 207 nodes and 405 edges, in which nodes represent protein names and edges represent protein interaction. To obtain the core PPI network, the initial PPI network was imported in Cytoscape 3.8.2 for visualization (Figure 4B). Then, the CytoNCA plugin was used to centrally analyze and evaluate the PPI network. In the first screening in Cytoscape 3.8.2, the selection criteria based on corresponding median values were set as follows: BC >8.27, CC >0.0236, DC >2.00, EC >0.01313, NC >0, and LAC >0. Only gene targets that met the criteria were retained and colored yellow in Figure 4B. After the first screening, the new network containing 66 nodes and 207 edges was retrieved (Figure 4C). Likewise, in the network in Figure 4C, only targets meeting the selection criteria (BC >467.16, CC >0.02395, DC >7.00, EC >0.07230, NC >0, and LAC >0) were retained and colored yellow in Figure 4C. After the second screening, the core PPI network with 16 nodes and 107 edges was obtained (Figure 2D). In the core PPI network, larger node sizes indicate higher degree values. The top eight target nodes were STAT3 (degree=25), HSP90AA1 (degree=21), TP53 (degree=19), AKT1 (degree=17), SRC (degree=17), RELA (degree=15), JUN (degree=15), and ESR1 (degree=15) (Figure 5), which may play a crucial role in the anti-AS effect of HLXLD.

Figure 4 PPI network of HLXLD-AS. (A) The interactive PPI network obtained from STRING database with species limited to “Homo sapiens”, the minimum required interaction score set to 0.98, and the independent target protein nodes hidden. (B) Original PPI network from STRING database imported to Cytoscape 3.8.2 to obtain a new network. It contains 207 nodes and 405 edges. (C) PPI network screened from (B) in Cytoscape 3.8.2. It contains 66 nodes and 207 edges. (D) Core PPI network screened from (C) in Cytoscape 3.8.2. It contains 16 nodes and 107 edges. Larger node sizes indicate higher degree values.

Figure 5 The node degree distribution diagram of Figure 4D. The x-axis shows gene names and y-axis displays the degree value of related targets.

GO and KEGG Pathway Enrichment Analysis

A total of 321 common targets were input into R software to carry out GO and KEGG enrichment analysis. The GO enrichment analysis results were reflected in biological process (BP), cellular component (CC), and molecular function (MF). Under the condition of P<0.05 and Q<0.05, 3420, 108, and 293 items were obtained from these three aspects, respectively. The top ten items of BP, CC, and MF were taken to make a visual bubble diagram with R software, as shown in Figure 6. The size of the dot indicates the number of genes contained under this entry. The color of the point represents the degree of enrichment, and the color from red to blue corresponds to the Q value from small to large. The top four entries of BP were enriched in response to molecule of bacterial origin (GO: 0002237), response to lipopolysaccharide (GO: 0032496), epithelial cell proliferation (GO: 0050673), and positive regulation of anion transport (GO: 1903793). The top four entries of CC were enriched in membrane raft (GO: 0045121), membrane microdomain (GO: 0098857), external side of plasma membrane (GO: 0009897), and focal adhesion (GO: 0005925). The top four entries of MF were enriched in signaling receptor activator activity (GO: 0030546), receptor ligand activity (GO: 0048018), DNA-binding transcription factor binding (GO: 0140297), and RNA polymerase II-specific DNA-binding transcription factor binding (GO: 0061629). The results showed that targets of HLXLD in the treatment of AS play a crucial role in BP, CC, and MF.

Figure 6 The bubble diagram of GO enrichment analysis of HLXLD-AS genes. The top ten GO entries of BP, CC and MF were shown, respectively. The abscissa shows gene ratio and the ordinates displays enriched GO entries. The size of the dot is the number of genes contained under this entry, the color of the point represents the degree of enrichment and color from red to blue corresponds to Q value from small to large.

For 321 common targets of HLXLD in the treatment of AS, KEGG enrichment analysis was carried out with R software, and a total of 187 items were obtained (the results met P<0.05 and Q<0.05). The results ranked in the top 30 according to Q value were visualized (Figure 7). The size of the dot represents the number of genes contained under this entry, the color of the point represents the degree of enrichment, and the color from red to blue corresponds to the Q value from small to large. The top four pathways with the highest gene counts were lipid and atherosclerosis, fluid shear stress and atherosclerosis, MAPK signaling pathway, and PI3K-Akt signaling pathway. The detailed pathway map of the MAPK signaling pathway is drawn in Figure 8. To elucidate the interrelationships of 108 compounds, the C-P-T network containing 453 nodes and 2191 edges was constructed with the top 20 pathways and 321 target genes on a system (Figure 9) by Cytoscape 3.8.2, in which lipid and atherosclerosis pathway (degree = 62), fluid shear stress and atherosclerosis (degree = 53), MAPK signaling pathway (degree = 52) and PI3K-Akt signaling pathway (degree = 48) had higher gene counts. These findings indicates that HLXLD exerts a protective role in the treatment of AS through multi-pathways and multi-targets.

Figure 7 The bubble diagram of KEGG enrichment pathway of HLXLD-AS genes. The top 30 pathways were shown. The abscissa shows gene ratio and the ordinates displays enriched GO entries. The size of the dot is the number of genes contained under this entry, the color of the point represents the degree of enrichment and the color from red to blue corresponds to Q value from small to large.

Figure 8 The map of MAPK signaling pathway. There are 31 target genes (red color) enriched in MAPK signaling pathway.

Figure 9 C-P-T network of HLXLD in the treatment of AS. It contains 453 nodes and 2191 edges. The rose diamond nodes represent four botanical drugs from HLXLD. The circular nodes of red, blue, green and yellow stand for compounds from Angelica sinensis, Salvia miltiorrhiza, Commiphora myrrha and Boswellia sacra, respectively. The circular purple nodes stand for AS-related target genes. The light blue V nodes stand for the top 20 pathways in Figure 7. Larger node sizes indicate higher degree values.

Molecular Docking Analysis

A molecular docking strategy was used to study the interaction between ligand and receptor, and to predict the binding mode and affinity.²⁸ In the present study, candidate target proteins were selected for molecular docking analysis based on being from the human species and owning higher degree values in the core PPI network (Figure 4D). Candidate active compounds were selected for docking based on owning higher values of degree in the C-D-T network (Figure 3). Thus, the top eight candidate target proteins with the highest degree values (STAT3, PDB ID: 5AX3; HSP90AA1, PDB ID: 1UYD; TP53, PDB ID:2VUK; AKT1, PDB ID: 3O96; SRC, PDB ID:2H8H; RELA, PDB ID: 6YOW; ESR1, PDB ID: 1QKT; JUN, PDB ID: 1JUN) and six active compounds with degree >53 (Table 2, Figure 10) were chosen to conduct molecular docking (Table 3, Supplement Figure 1). According to the minimum binding energy, as shown in Table 3, STAT3-dihydrotanshinone I (−8.5 kcal/mol), HSP90AA1-dihydrotanshinone I (−10.1 kcal/mol), TP53-luteolin (−7.9 kcal/mol), AKT1-luteolin (−9.8 kcal/mol) were chosen to further analyze their binding mode, binding affinity and critical interaction by LigPlus 2.2 (2D) and Pymol 2.5 (3D) (Figure 11). According to Figure 11A-D, dihydrotanshinone I interacts with STAT3 by hydrogen bonds to Leu438(A) and Lys370(A) residues, hydrophobic bonds to Asp369(A), His437(A), Leu436(A), Val490(A), Asp371(A), Thr440(A) and His457(B) residues; dihydrotanshinone I interacts with HSP90AA1 by hydrophobic bonds to gly135(A), Leu107(A), Tyr139(A), Trp162(A), Phe138(A), Val150, Met98(A), Asn51(A) residues; luteolin interacts with TP53 by hydrogen bonds to Ser99(A), Lys101(A), Gln100(A), Arg196(B), Ser166(A), Ala138(B) residues, hydrophobic bonds to Glu198(B), Met237(B), Asn235(B), Asp186(B), Gly199(B) and Leu201(B) residues; luteolin interacts with AKT1 by hydrogen bonds to Ser205(A), Thr211(A), Ile290(A), Gln79(A), Tyr272(A) residues, hydrophobic bonds to Trp80(A), Leu264(A), Val270(A), Asn54(A), Val271(A), Leu210(A) residues, respectively.

Table 3 Docking Binding Energy Results of Key Targets and Main Active Components (Kcal/Mol)

Figure 10 The chemical structure and molecular formula of quercetin, dihydrotanshinone I, pelargonidin, luteolin, guggulsterone and β-sitosterol.

Figure 11 The picture of molecule docking model. (A-D) Molecular docking of dihydrotanshinone I with STAT3, dihydrotanshinone I with HSP90AA1, luteolin with TP53 and luteolin with AKT1, respectively. The hydrogen bonds were represented by yellow dotted lines, and the length was marked around the lines.

Molecular Dynamics Simulation

Molecular dynamics simulation was used for studying the binding stability between the protein and the molecule. In the present study, STAT3-dihydrotanshinone I, HSP90AA1-dihydrotanshinone I, TP53-luteolin, and AKT1-luteolin were chosen to further analyze their binding stability. It can be seen from the RMSD graph (Figure 12): (1) The system of molecules and proteins is constantly fluctuating in the range of 0–10 ns during the kinetics process, indicating that the system is unstable and the protein is combined with their corresponding protein architecture when the phase changes constantly; (2) In the range of 10–100 ns, the RMSD value has been in a stable state, indicating that the binding between the protein and the molecule is extremely stable, and there is a good binding state between the two systems.

Figure 12 The picture of RMSD value of molecular complex system in 100ns time interval. (A) The dynamic stability of dihydrotanshinone I with STAT3 (blue line) and with HSP90AA1 (red line). (B) The dynamic stability luteolin with TP53 (blue line) and with AKT1 (red line). DHT: dihydrotanshinone I.

Dihydrotanshinone I Suppressed the Accumulation of Lipid and Inflammation Induced by Ox-LDL in RAW264.7 Cells

The activation of macrophages, as indicated by lipid accumulation and increased inflammatory cytokines secretory, plays a key role in the pathological process of AS. In the present study, oil red O staining was employed to determine the intracellular lipid accumulation in RAW264.7 macrophage-derived foam cells. As shown in Figure 13A, there was a large number of lipids accumulated in the ox-LDL group compared with the control (untreated) and dihydrotanshinone I groups, which means that the foam cells model was successfully induced by 50 μg/mL ox-LDL. Furthermore, the accumulation of lipids in ox-LDL-induced foam cells was decreased by the addition of dihydrotanshinone I in a dose-dependent manner. To investigate the anti-inflammatory effect of dihydrotanshinone I, we evaluated the secretory of several inflammatory cytokines in RAW264.7 cells using ELISA analyses. As shown in Figure 13B-D, the increased secretory of TNF-α, IL-1β, and MCP-1 induced by ox-LDL in RAW264.7 cells was significantly reduced by pretreating with dihydrotanshinone I.

Figure 13 The effect of dihydrotanshinone I on lipid accumulation and the secretory levels of inflammatory factor in ox-LDL-induced RAW264.7 cells. (A) Oil red O staining of RAW264.7 cells; (B-D) the effect of dihydrotanshinone I on the secretory levels of TNF-α, IL-1β, MCP-1, respectively. DHT: dihydrotanshinone I. Data are mean±S.E.M. n=3. **P < 0.01 vs Control; #P < 0.05, ##P < 0.01 vs ox-LDL.

Dihydrotanshinone I Inhibited the STAT3 Phosphorylation and MAPK Signaling Pathway Activation in Ox-LDL-Treated RAW264.7 Cells

The results of network pharmacological analysis indicated that STAT3 is the key target and the MAPK signaling pathway is one of the important pathways of HLXLD against AS. To further verify the potential mechanism of dihydrotanshinone I in the treatment of AS, the phosphorylation of STAT3 and MAPK signaling key proteins (ERK, JNK, and p38) were determined. As shown in Figure 14, the phosphorylation expression of STAT3, ERK, JNK, and p38 were markedly increased in ox-LDL-treated RAW 264.7, whereas pretreating with dihydrotanshinone I can significantly inhibit the phosphorylation of STAT3, ERK, JNK, and p38. These results suggested that dihydrotanshinone I attenuated AS by inhibiting the activation of STAT3 and MAPK signaling pathways.

Figure 14 The effect of dihydrotanshinone I on the protein secretion of p-STAT3 and MAPK signal pathway in ox-LDL-induced RAW264.7 cells. (A) Representative immunoblotting images of p-STAT3, STAT3 and GAPDH; (C) Representative immunoblotting images of p-ERK, ERK, p-JNK, JNK, p-p38, p38 and GAPDH; (B, D-F) Gray value statistics of corresponding proteins. DHT: dihydrotanshinone I. Data are mean±S.E.M. n=3. **P < 0.01 vs Control; #P < 0.05, ##P < 0.01 vs ox-LDL.

Discussion

TCM has been widely applied in treating multiple diseases in clinical practice in China. However, the underlying mechanism has not been easily elucidated clearly from the perspective of modern medicine for most of the TCM. The reason may lie in that the overall effect of TCM originates from the interaction between different components but is not limited to the superposition of single components.³⁰ In recent years, network pharmacology provides a new method that appropriately represents the overall diagnosis and treatment concepts for investigating the compatibility mechanism of TCM,²⁶ for example, the mechanism study of Spica Prunellae against thyroid-associated ophthalmopathy,³¹ Bu Yang Huan Wu Decoction in the treatment of ischemic stroke,³² and Liu Wei Di Huang Pill in treating type 2 diabetes mellitus.³³ Given the complex pathological mechanism of AS, more and more TCM with the characteristics of multiple targets are developed to treat AS. More importantly, the mechanism of several TCM against AS has been elucidated using network pharmacology. Zhang et al found that the Danshen-Shanzha herb-pair can target 41 proteins and 16 signal pathways related to inflammation, lipid metabolism, and endothelial protection in the treatment of AS.¹⁵ Yeong et al screened 18 compounds in Gyejibokryeong-hwan and 213 AS-related targets, suggesting that the TNF, HIF-1, FoxO, and PI3K-Akt signal pathways mediate the beneficial effects of Gyejibokryeong-hwan on AS.³⁴

HLXLD, a well-known traditional Chinese herbal prescription, possesses the ability to prevent AS, but its pharmacological mechanism has not been elucidated.^18,19 Those successful examples mentioned above in the application of TCM for treating AS inspired us to explore the mechanisms of HLXLD in AS by using network pharmacology analysis. The results of the present study showed that 108 active ingredients of HLXLD were found to act on 321 different targets associated with AS. According to the C-D-T network, quercetin, dihydrotanshinone I, pelargonidin, luteolin, guggulsterone, and β-sitosterol are highly connected with targets, which can be defined as vital compounds in HLXLD. Previous studies have determined that these compounds exhibit protective effects in the cardiovascular system,^23,35–39 indicating the beneficial effects of HLXLD in treating AS can be attributed to the combination of these core compounds.

As far as each core active compound is concerned, quercetin, dihydrotanshinone I, pelargonidin, luteolin, guggulsterone, and β-sitosterol have been reported to possess excellent anti-AS effects. The mechanism of quercetin against AS involves gut metabolic regulation, enhancing autophagy and delaying senescence through the MST1 pathway.²⁰ Moreover, quercetin can restore the dysfunction of endothelial and lipid metabolism by regulating the ABCA1, LXR-α and PCSK9 expressions, and attenuating NADPH oxidase-dependent superoxide (O^2.-) production.⁴⁰ Dihydrotanshinone I, whose mechanisms against AS involved in inhibiting LOX-1 mediated by NOX4/NF-κB signaling pathways, stabilizes vulnerable plaque by suppressing RIP3-mediated necroptosis of macrophage.^22,36 Besides quercetin and dihydrotanshinone I, pelargonidin exhibits potential preventive effects toward AS via inhibiting HASMC proliferation and migration.³⁷ Luteolin can treat AS through a variety of mechanisms including lipid metabolism renovation, inactivation of STAT3 and NLRP3 inflammasomes, macrophage polarization, as well as inhibition of AMPK-SIRT1 signaling pathway.^21,41 TMA/TMAO/FMO3 pathway mediates the protective effect of guggulsterone against AS.³⁹ Moreover, its antioxidant and anti-inflammatory activities may be another important mechanism.^23,42

The main targets in the core PPI network, STAT3, HSP90AA1, TP53, and AKT1, play key roles in the occurrence and development of AS. STAT3, the protein node with the highest degree value in the core PPI network, is the signal transducer and activator of transcription 3 and plays an essential role in the process of AS including endothelial cell dysfunction, macrophage polarization, inflammation and immunity.⁴³ HSP90AA1, a heat shock protein, plays a protective role in the homeostasis of the vessel wall but has an impact on immunoinflammatory processes in pathological conditions involved in the development of AS.⁴⁴ P53 is a tumor suppressor protein, and the expression and transcriptional activity of P53 can be tightly regulated to treat AS.^45,46 AKT1, as a critical node in PI3K-PKB/Akt and pPI3K/Akt/mTOR signaling pathway, is involved in the progress and development of AS.^47,48

Taken together, the above results showed that the core components and key targets of HLXLD in the C-D-T network play an important role in the treatment of AS. To integrate functions of 321 protein target genes, the GO enrichment analysis and KEGG pathway enrichment analysis were further conducted. GO results suggested that 321 target genes were enriched in BP, CC and MF, and mainly in BP. KEGG results suggested that 296 of the 321 target genes were enriched in 187 signal pathways, indicating a therapeutic effect of HLXLD on AS through multiple pathways. Among them, 62 of the 321 target genes enriched in lipid and atherosclerosis pathways, indicating lipid metabolism plays a pivotal role in the process of AS. Furthermore, research showed that the regulation of low-density lipoprotein (LDL) and cholesterol (LDL-C) is still the primary approach for the management and prevention of AS.⁷ These findings are congruent with our view that lipid and atherosclerosis act as an important role in the anti-AS mechanism of HLXLD. In addition, MAPK and PI3K-AKT signaling pathways may be excellent therapeutic targets for HLXLD against AS. As one of the main information transmission pathways in cells, MAPK signaling pathway can regulate a variety of cell functions, such as proliferation, differentiation, transformation and death, and a series of inflammatory reactions in AS through activation.⁴⁹ For example, lncRNA VINAS regulates AS by modulating NF-κB and MAPK signaling.⁵⁰ The inhibition of microRNA‐103 attenuates inflammation and endoplasmic reticulum stress in AS by disrupting the PTEN‐mediated MAPK signaling.⁵¹ Other studies have shown that TCM or TCM extract can regulate AS through the MAPK pathway, for instance, Tianxiangdan Granule alleviates atherosclerosis via p38 MAPK signaling pathways. Icariin attenuates high-cholesterol diet-induced AS in rats by preventing p38 MAPK signaling pathway.⁵² PI3K-AKT pathway can participate in the regulation of AS by affecting macrophage polarization, lipid metabolism, autophagy and other functions,⁵³ which is a potential target for the treatment of AS. Wang et al suggested that nepsilon-carboxymethyl-lysine induces foam cell apoptosis in diabetic atherosclerosis by inhibiting the PI3K/AKT pathway.⁵⁴ Danhong Injection, a TCM preparation, attenuates high-fat-induced AS and macrophage lipid accumulation by regulating the PI3K/AKT Insulin Pathway.⁵³

Numerous studies demonstrated that a TCM network pharmacology approach can clarify the interaction of multi-component to multi-target and multi-pathway. In this study, eight key target genes (STAT3, HSP90AA1, TP53, AKT1, SRC, RELA, ESR1and JUN) and four key pathways (lipid and atherosclerosis, fluid shear stress and atherosclerosis, MAPK signaling pathway, and PI3K-Akt signaling pathway) were regulated by three or more compounds in HLXLD. However, these findings were predicted by a network pharmacological approach and require further determination. Therefore, computer‑aided prediction (molecular docking strategy) and experimental research verification (in vitro validation experiment) were conducted to further validate the prediction by network pharmacology analysis. Molecular docking results suggested that quercetin, dihydrotanshinone I, pelargonidin, luteolin, guggulsterone and β-sitosterol could be stably combined with their corresponding representative targets. In vitro experiments showed that the high phosphorylation levels of STAT3, ERK, JNK and p38 were reversed by dihydrotanshinone I, further validating the results of C-D-T network and KEGG. Previous research reports showed that mechanisms of dihydrotanshinone I against AS are involved in the inhibition of LOX-1 mediated by NOX4/NF-κB signaling pathways, and stabilizing vulnerable plaque by suppressing RIP3-mediated necroptosis of macrophage.^22,36 The present study determined that dihydrotanshinone I suppressed macrophage activation via inhibiting STAT3 and MAPK signaling pathways. Our results verify the new targets and pathways of dihydrotanshinone I in the treatment of AS.

However, the present study has several limitations based on the network pharmacology evaluation method guidance⁵⁵ and published literature.⁵⁶ Firstly, limited database usage. There are more TCM databases in addition to the several ones that we used to search for bioactive ingredients and potential targets, such as Traditional Chinese Medicine on Immuno-Oncology (TCMIO, http://tcmio.xielab.net/), SymMap v2 (http://www.symmap.org/) and Traditional Chinese Medicine Integrated Database (TCMID, https://ngdc.cncb.ac.cn/databasecommons/database/id/437). Secondly, the content of the main bioactive ingredients is lacking. In the collection of chemical compounds of HLXLD, we conducted the study based on public databases, which lack relevant quantitative data. Probably, some of them may have only trace content and contribute little to the therapeutic effect of HLXLD. Therefore, future research on content determination should be carried out. Thirdly, the specifically binding relationship between active compounds and their targets was not verified. Molecular docking is a strategy to predict the combination of active compounds with core targets, rather than validate the network data. Therefore, more methods, such as Surface Plasmon Resonance (SPR) technology, can be used to validate whether active compounds bind specifically to their targets. Lastly, the validation experiment was insufficient. Dihydrotanshinone I, though determined as one of the most important bioactive ingredients of HLXLD against AS, could not completely stand for HLXLD. Thus, animal experiments and clinical research verification need to be further carried out to explore the potential molecular mechanism of HLXLD in treating AS.

In conclusion, we firstly utilized network pharmacology and molecular docking strategy to explore the underlying mechanism of HLXLD in treating AS. The main active components are quercetin, dihydrotanshinone I, pelargonidin, luteolin, guggulsterone and β-sitosterol, as well as key targets are STAT3, HSP90AA1, TP53, AKT1, SRC, RELA, ESR1 and JUN. The interaction between active components and targets is dominated by hydrophobic interaction. Additionally, the protective effect and underlying mechanism of dihydrotanshinone I on AS were validated by cell experiments. The above results indicate that HLXLD is effective in the AS treatment through multi-component, multi-target and multi-pathway, which provides a scientific basis for explaining the detailed mechanism of HLXLD against AS.