Network Pharmacology and Experimental Verification to Explore Cinnamomi Cortex against Steroid-induced Osteonecrosis of the Femoral Head

Summary

This study aims to explore the potential mechanisms of Cinnamomi Cortex in the treatment of osteonecrosis of the femoral head by integrating network pharmacology, molecular dynamics simulations, and animal experiments.

Abstract

Steroid-induced osteonecrosis of the femoral head (SONFH) causes severe pain and limited mobility, significantly impairing patients' quality of life. Cinnamomi Cortex (CC) has been shown to effectively alleviate this condition, yet its mechanism of action remains unclear. This study aims to identify the active compounds of CC and explore their mechanisms in SONFH. Active constituents were screened using the HERB 2.0, PubChem, and SwissADME databases, and their corresponding targets were predicted using the Swiss Target Prediction database. Targets for SONFH were identified by intersecting targets from GEO, DisGeNET, GeneCards, and OMIM databases with the compound-related targets. A protein-protein interaction (PPI) network was constructed using the STRING database, and GO and KEGG enrichment analyses were conducted via the DAVID database. The most promising compound-target interactions were validated through molecular docking (MD) and molecular dynamics simulations. The researchers identified 563 potential targets, including 61 SONFH-related targets, with AKT1, HIF-1α, and STAT3 serving as central nodes. KEGG enrichment analysis highlighted the HIF-1α signaling pathway as a key mechanism. Furthermore, animal experiments demonstrated that the active fraction of CC effectively mitigated femoral head structural damage in a mouse model with SONFH. The findings suggest that CC may improve SONFH by coordinating hypoxia adaptation and regulating angiogenesis and osteogenesis.

Introduction

Steroid-induced osteonecrosis of the femoral head (SONFH) is a critical osteoarticular disorder, which is caused by the disruption or reduction in blood supply to the femoral head from various¹. The death of bone cells and bone marrow components leads to structural damage within the femoral head and functional impairment of the hip joint². Clinically, patients typically present with hip pain and limited mobility³. Without effective treatment, 80% of SONFH patients will progress to femoral head collapse, requiring total hip replacement⁴. It places a severe psychological pressure and a heavy economic burden on patients⁵. Clinical practice recommends using a combination of anticoagulants, fibrinolytic medicines, vasodilators, and lipid-lowering medicines to treat SONFH, which have demonstrated some clinical potential but generally limited therapeutic efficacy⁶.

Traditional Chinese medicine (TCM) has increasingly demonstrated its positive impact on treating SONFH⁷. TCM shows potential for symptom relief, disease control, and enhanced joint mobility and patients’ quality of life⁸. Clinically, CC and its components have been extensively utilized in improving SONFH, showing notable effectiveness⁹. Numerous bioactive substances found in CC have a variety of biological actions, including anti-inflammatory, antioxidant, angiogenesis-promoting, and microcirculation-improving properties¹⁰. Previous research in bone tissue engineering and metabolism has established that cinnamaldehyde, a principal active constituent of CC, can modulate signaling pathways related to bone remodeling¹¹. The initial study indicated that Yougui Pills (YGPs) have therapeutic potential for SONFH by promoting angiogenesis and enhancing inflammatory responses¹². As the sovereign drug in YGPs, CC can warm yang and benefit qi and unblock meridians, producing beneficial therapeutic effects on SONFH¹³. Nevertheless, the precise mechanisms underlying its efficacy remain unclear. Clarifying the multi-target regulatory network that CC uses to ameliorate SONFH would not only help to understand its pharmacological nature but also encourage its sensible clinical use and the creation of related amelioration. Integrating bioinformatics with network pharmacology provides an effective approach to elucidating mechanisms of action¹⁴.

Network pharmacology can identify bioactive components in herbs and predict the associations between these drug components and gene targets¹⁵. MD is employed to validate the binding interactions between candidate active compounds and key therapeutic targets. Molecular dynamics uses Newtonian mechanics to evaluate ligand-receptor binding stability and flexibility through movement simulation¹⁶. The authors comprehensively investigated the potential effects of CC on SONFH using network pharmacology, MD, and molecular dynamics simulations. The findings provide a reference for future in-depth research into the pharmacodynamic material basis and mechanism of CC in improving SONFH. The study flowchart is shown in Figure 1.

Protocol

All experimental protocols were approved by the Animal Experimental Ethics Committee of Zhejiang Chinese Medicine University (IACUC-20240708-22) and complied with the Guide for the Care and Use of Laboratory Animals issued by the National Institutes of Health. In this experiment, twenty 10-week-old female C57BL/6J mice weighing 20 to 22 g were used. These mice were obtained from the animal center of Zhejiang Chinese Medical University. See the Table of Materials for a list of all the reagents, equipment, and software used in this protocol.

Screening of active compounds and targets of CC
Chemical compounds of CC were identified using the Herb 2.0¹⁷ with the keyword “CC.” Chemicals without a PubChem ID or with the same ID were filtered out. Then merge and remove duplicates to obtain the target sites corresponding to CC. Perform an initial screening using the PubChem database¹⁸ based on Lipinski's rule of five (Mw ≤ 500, miLogP ≤ 5, nOHNH ≤ 5, nOH ≤ 10)¹⁹. The authors use the PubChem database to determine the SMILES representation for each chemical compound. Using the SwissADME database, the authors select compounds with the restriction of “High” GI absorption and ≥2 “Yes” values in Druglikeness²⁰. Utilize the Swiss-Target Prediction database²¹ to extract target proteins that have a probability higher than zero. Then merge and remove duplicates to obtain the target sites corresponding to CC.

Compilation of gene targets associated with SONFH
SONFH-related targets were obtained from databases: DisGeNET database²²; GeneCards database²³; OMIM database²⁴, and GEO database²⁵. The authors downloaded and deduplicated the SONFH genes from the OMIM and DisGeNET databases. A total of 233 disease-related targets were obtained from the GeneCards database after removing duplicates. All retrieved targets had relevance scores greater than 0 and were included for subsequent analyses. The minimum relevance score among the retrieved genes was 6.48. Differentially expressed genes (DEGs) related to SONFH were obtained from the GSE123568 series on the GPL15207 platform in the GEO database. Batch effects were corrected using the limma package in R, and differentially expressed genes (DEGs) were identified with criteria of |logFC| > 1 and P < 0.05²⁶. The ggplot2 package was utilized to create a volcano plot for visualizing the distribution of DEGs, and a heatmap was produced to present the findings. The SONFH disease target library was established by removing duplicate targets using the Venn package in R.

Establishment of the PPI network
Based on the anticipated targets of the active components of CC and SONFH-related targets, a Venn diagram was created by Venny 2.1.0²⁷. With the aim of acquiring PPI data, the authors imported the targets into the STRING database²⁸. The limitation of the organism’s screening criterion was “Homo sapiens” with the confidence index ≥0.4. The authors used Cytoscape 3.7.2 to establish the PPI network key targets. Degree centrality (DC) was used to examine PPI network key targets. The parameters filter is more than twice the median value. With the purpose of specifying the key targets further, the add-on to Cytoscape was employed.

Enrichment analyses were conducted using gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) databases
The common targets of CC and SONFH were analyzed for GO and KEGG enrichment analyses through the DAVID database²⁹. In order to visualize the top 10 GO terms and the top 20 KEGG pathways, Wei Sheng Xin³⁰ was used to draw the enrichment dot bubble chart.

Network construction
The molecular mechanisms of CC in ameliorating SONFH were elucidated through the herb-compound-target (H-C-T) and compound-target-pathway (C-T-P) networks. Cytoscape 3.10.3 software illustrated the networks. The H-C-T network was developed using the active compounds from CC and their shared targets. Then use the Analyze Network tool. To enhance understanding of the relationships among pathways, compounds, and targets, the top 20 pathways, along with their related targets and compounds, were organized into a C-T-P network by Network Tools.

MD verification
The core targets: HIF-1α, STAT3, ESR1, AKT1, SRC, ERBB2, CASP3, and EGFR were obtained from the PPI analysis. The target was searched by the Uniport database³¹ with the limitation of Human and reviewed by the RCSB PDB database³². Human protein structures with relatively high resolution were selected. The active ingredients were identified from the C-T-P network. Their corresponding 3D structures were downloaded in mol2 format from the PubChem database. Each file was opened in Chem3D and energy minimized. Then the authors used the AutoDockTools³³ software to perform the dehydration, hydrogenation, and calculate Gasteiger charges of the receptor protein. Both ligands and receptors were kept in the PDBQT format. Based on the spatial location and binding capability, the authors evaluate the feasibility and the stability of docking. Via the AutoDock Vina, the docking grid box and macromolecule were predicted. In this research, all the rotatable bonds of the ligand were allowed to rotate freely, while the receptors were set as rigid. The docking score indicates the binding affinity between the receptor and the ligand; the lower the score, the higher the binding affinity. In accordance with this principle, the authors select the conformations with the most favorable binding energy to investigate binding interactions of ligands with proteins. The visual results were exhibited by PyMOL.

Molecular dynamics simulation
Molecular dynamics simulations were performed using GROMACS³⁴ with the CHARMM36 force field. The protein-ligand complex was centered in a cubic periodic simulation box with a minimum solute-box distance of 1.0 nm, and the system was solvated using the SPC216 water model. The system was neutralized by adding appropriate Na⁺ and Cl⁻ counterions. Energy minimization was performed sequentially using the steepest descent algorithm followed by the conjugate gradient method to remove unfavorable atomic contacts. Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method with a cutoff distance of 10 Å. All covalent bonds involving hydrogen atoms were constrained using the LINCS algorithm. Temperature was maintained at 300 K using the V-rescale thermostat, and pressure was controlled at 1 bar using the Berendsen barostat with isotropic coupling. Position restraints were applied to the ligand with a force constant of 1000 kJ·mol⁻¹·nm⁻² during the equilibration stage. After energy minimization, the system was equilibrated with 2 ns NVT and 2 ns NPT simulations, followed by a 100 ns production molecular dynamics simulation with a 2-fs time step. Atomic coordinates were recorded every 10 ps for subsequent trajectory analysis.

Experimental validation in a SONFH mouse model
Animals and the SONFH model establishment
Twenty female C57BL/6J mice, aged 10 weeks and weighing 20 to 22 g, were obtained from the animal center of Zhejiang Chinese Medical University. The mice underwent a 7-day acclimation period in a controlled environment and were given unrestricted access to food and water. All mice were randomly allocated to the groups (n = 5 per group): the control group, the SONFH group, the SONFH + CC low-dose (CC-Low, 7.5 g/kg/d) group, and the SONFH + CC high-dose (CC-High, 15 g/kg/d) group. The SONFH model was established as previously described. Briefly, mice in the SONFH, SONFH + CC-Low, and SONFH + CC-High groups received two intravenous injections of lipopolysaccharide (LPS; 20 µg/kg) on day 0. Subsequently, three intramuscular injections of methylprednisolone (MPS; 40 mg/kg) were administered at 24-hour intervals, starting 24 h after the LPS injections. Mice in the treatment groups were administered CC extract via oral gavage daily for 6 weeks, starting from the day of the last MPS injection. The control group mice received equivalent volumes of saline at corresponding time points.

Sample collection and preparation
At the end of the 6-week treatment period, mice were euthanized. Bilateral femoral heads were carefully dissected and harvested. For each animal, the left femoral head was fixed for 48 h at 4 °C for subsequent decalcification and paraffin embedding. The right femoral head was fixed and then directly used for Micro-CT scanning without decalcification.

Histological analysis (ABH Staining)
After decalcification and paraffin embedding, femoral head sections were stained using an Alcian Blue/Hematoxylin (ABH) protocol to assess osteonecrosis changes. Briefly, sections were stained in 1% Alcian Blue (pH 2.5) for 30 min, rinsed, and then counterstained with Harris hematoxylin. After dehydration and mounting, slides were imaged under a light microscope. Quantification of osteonecrosis was performed by two blinded observers using ImageJ software. The ratio of empty lacunae was calculated as (number of empty lacunae/total number of lacunae) × 100% in three randomly selected high-power fields per sample within the subchondral region. The ratio of pyknotic nuclei was determined similarly.

Micro-computed tomography (Micro-CT) analysis
The three-dimensional bone microstructure of the femoral head was analyzed using a high-resolution micro-CT scanner. Fixed samples were scanned at a resolution of 10 µm (70 kV, 114 µA). A standardized spherical volume of interest (VOI) encompassing the primary weight-bearing region was reconstructed and analyzed using CTAn software. The following morphometric parameters were quantified: Bone Volume/Total Volume (BV/TV), Trabecular Thickness (Tb.Th), and Trabecular Separation (Tb.Sp).

Immunofluorescence (IF) staining
To evaluate the bone microenvironment, immunofluorescence staining was performed on paraffin sections. After antigen retrieval and blocking, sections were incubated overnight at 4 °C with the following primary antibodies: rabbit anti-HIF-1α (1:200), rabbit anti-ALP (1:300), and rabbit anti-VEGF (1:150). After washing, sections were incubated with a mixture of fluorophore-conjugated secondary antibodies: Alexa Fluor 488-labeled goat anti-rabbit IgG, Alexa Fluor 555-labeled goat anti-rabbit IgG. Nuclei were counterstained with DAPI. Images were captured using a fluorescence microscope under consistent exposure settings. Relative fluorescence intensity for each marker was quantified using ImageJ software across three fields per sample.

Statistical analysis
All quantitative data are presented as mean ± standard deviation (SD). One-way analysis of variance (ANOVA) followed by Tukey’s post hoc test was performed using GraphPad Prism software (version 9.0) to determine statistical significance among groups. A P-value of less than 0.05 was considered statistically significant.

Results

Identification of active components and predicted targets of CC and development of the “H-C-T” network.
To explore the mechanisms of key components and targets of CC in SONFH amelioration, the authors used the HERB 2.0 database to identify all CC components screened. 85 probable bioactive compounds were found out of the 209 candidate active constituents of CC (Supplementary Table 1). Using the Swiss Target Prediction platform and eliminating duplicate and invalid entries, 563 putative targets corresponding to these CC components were obtained (Supplementary Table 2). The authors have identified the active components and targets of CC.

To investigate the molecular mechanisms involved in the onset and progression of SONFH and to identify potential biomarker targets. Analysis of the GSE123568 dataset in the GEO database identified DEGs associated with SONFH, including 207 upregulated and 216 downregulated genes. The volcano plot of the 425 DEGs is shown in Figure 2A, whereas the expression patterns of the top 60 DEGs ranked by expression level are presented as a heatmap in Figure 2B, with color intensity reflecting the log-transformed expression values. The databases DisGeNET, OMIM, and GeneCards provided 711 targets linked to SONFH. Combining the 425 DEGs from the GEO database, a total of 1136 targets are linked to SONFH. Following the elimination of duplicates, 1115 SONFH-related targets were confirmed, as illustrated in Figure 2C. The authors have preliminarily established a disease target library for SONFH.

To identify targets and primary components, an intersection analysis between the 563 targets of CC and the 1115 SONFH-related targets identified 61common targets (Figure 3A and Supplementary Table 3), which were brought into Cytoscape 3.10.3 to build the “H-C-T” network (Figure 3B). This network comprised 147 nodes and 413 edges. Topological analysis indicated that (s)-4-nonanolide (degree:16), Isohomogenol (degree:15), Neryl acetate (degree:14), Melilotocarpan A (degree:14), 3-methoxycinnamaldehyde (degree:13) showed the highest connectivity with protein targets.

Development and topological evaluation of the PPI network.
To further explore the targets of CC in improving SONFH. Cytoscape 3.10.3 was used to view the PPI network after 61 identical targets were put into the STRING database. There were 463 edges and 60 nodes in the network. In Figure 4A, nodes transition in color and size from lighter and smaller to darker and larger, indicating an increase in degree from low to high. The MCODE plug-in facilitated the clustering of these targets into two functional modules (Figure 4A), among which cluster 1 contained 23 nodes and 216 edges, with a score of 19.636. To further characterize the network architecture, DC, BC, and CC were used as topological metrics. With the median values of these parameters set as cut-off points, the authors identified 8 active targets (Figure 4B). The CytoHubba plugin utilized the MCC algorithm to rank nodes, identifying the top 10 genes. (Figure 4C).

GO enrichment analysis
To further explore the mechanisms of CC in ONFH amelioration, a GO enrichment analysis was conducted on 61 potential targets, utilizing BP, CC, and MF. A total of 435 GO terms were enriched, including 293 BP terms, 41 CC terms, and 101 MF terms. These findings are shown in Supplementary Tables 4, 5, and 6. Figure 5A visually presents the top 10 enriched BP, CC, and MF terms using a bubble chart. The findings indicated that BP terms were mainly associated with signal transduction, positive regulation of transcription by RNA polymerase II, negative regulation of apoptotic process, positive regulation of DNA-templated transcription, and negative regulation of transcription by RNA polymerase II. The analysis revealed that the MF terms were mainly involved in protein binding, metal ion binding, and identical protein binding. The CC results demonstrated that the majority of targets were chiefly located in the cytoplasm, membrane, and plasma membrane.

KEGG analysis
To shed information on the potential biochemical mechanisms by which CC may improve SONFH, KEGG pathway enrichment analysis was conducted to elucidate the biochemical mechanisms through which CC may ameliorate SONFH. Using a screening criterion of P < 0.05 and FDR < 0.05, 48 KEGG enrichment entries were identified (Supplementary Table 7). The top 20 significantly enriched pathways, determined by Fold Enrichment and count, were visualized using bubble and bar plots (Table 1, Figure 5B and 5C). Using Cytoscape 3.10.3, a C-T-P network with 167 nodes and 589 edges was constructed (Figure 5D). CC's effects on SONFH are mainly linked to Kaposi sarcoma-associated herpesvirus infection, HIF-1 signaling pathway, and Lipid and atherosclerosis, as indicated by GO functional and KEGG pathway enrichment analyses.

MD Verification
To evaluate the potential interactions between active compounds and key targets, the 10 compounds from the C-T-P network, including anethole, melilotocarpan A, neryl acetate, caprylic acid, isohomogenol, 3-methoxycinnamaldehyde, (s)-4-nonanolide, myristicin, borneol, and cinnamyl acetate, were docked with the key targets AKT1, HIF-1α, STAT3, ESR1, CASP3, SRC, and EGFR. The complete docking results for all compound–target pairs are summarized in Supplementary Table 8, and detailed information on the compounds, targets, and docking parameters is provided in Supplementary Table 9. Figure 6 presents a heatmap of the binding energy distribution.

Overall, several compound-target pairs exhibited relatively favorable binding affinities, with lower binding energy values indicating more stable predicted interactions. Notably, the docking results showed clear heterogeneity in binding affinities among both compounds and targets. Certain compounds, such as melilotocarpan A, consistently exhibited relatively strong binding affinities across multiple targets, whereas others displayed more moderate or weak interactions, suggesting that different compounds may contribute unequally to the predicted pharmacological effects.

In addition, variability was also observed across targets. For example, certain targets such as STAT3 and EGFR demonstrated relatively moderate or weak binding affinities with multiple compounds, with some values approaching -5.0 kcal/mol, compared to targets such as AKT1 or SRC. This pattern indicates that not all core targets necessarily function as direct high-affinity binding partners of the identified compounds and may instead play indirect regulatory roles within the interaction network. Such differences may be related to variations in structural compatibility between ligands and protein binding sites, as well as the inherent characteristics of the targets.

Among the evaluated pairs, AKT1 exhibited the lowest binding energy with melilotocarpan A (-9.6 kcal/mol), suggesting a potentially favorable interaction. Figure 7 illustrates the predicted binding mode of this complex. Specifically, SER205 of AKT1 forms hydrogen bonds with melilotocarpan A, anethole, and myristicin. In addition, melilotocarpan A was predicted to form hydrogen bonds with GLY309 of HIF-1α, LEU438 and THR440 of STAT3, SER433 and ARG412 of ESR1, and ARG500 and GLU510 of SRC.

To identify representative candidate compounds, the lowest binding energy for each target was used, and a binding energy threshold of ≤ -5.0 kcal/mol was set to indicate relatively stable interactions. Based on these criteria, melilotocarpan A, anethole, and myristicin were identified as potential key compounds. Overall, the majority of active components of CC showed potential interactions with the selected therapeutic targets, with some compounds demonstrating relatively stronger binding trends across multiple targets. However, it should be noted that MD is a simplified computational approach that provides preliminary predictions of potential interactions and may not fully account for protein flexibility and complex biological environments. Therefore, these results should be interpreted with caution, and the observed interactions do not constitute definitive evidence of direct binding.

Molecular dynamics simulation
To further assess the complexes' stability, the authors selected the AKT1-melilotocarpan A, HIF-1α-melilotocarpan A, and STAT3-melilotocarpan A complexes for molecular dynamics simulations. Among the predicted targets, STAT3 and HIF-1α were selected for MD simulations based on their central roles in the PPI network and their biological relevance to SONFH pathology. STAT3 showed one of the highest connectivity degrees in the PPI network, indicating a potential regulatory role in multiple signaling pathways. HIF-1α was selected because KEGG enrichment analysis identified the HIF-1 signaling pathway as one of the most relevant pathways associated with angiogenesis and hypoxia adaptation in osteonecrosis.

At the same time, the authors focused on AKT1 as a representative target because of their strong docking scores and direct relevance to the pathological features of SONFH, particularly hypoxia-induced angiogenesis and bone regeneration. Due to computational resource limitations, three representative complexes were selected for detailed MD simulations, which is a commonly adopted strategy in network pharmacology-based molecular simulation studies. Root mean square deviation (RMSD) effectively evaluates the conformational stability of protein-ligand complexes, where lower values indicate greater structural stability. As shown in Figure 8A, the AKT1-melilotocarpan A complex completed relaxation within the first 12 ns and subsequently reached a stable plateau, with only a transient fluctuation observed around 50-55 ns before rapidly returning to equilibrium. The overall RMSD throughout the simulation was about 0.369 nm. Figure 8B illustrates that the HIF1A-melilotocarpan A complex underwent relaxation within the first 16 ns, followed by the maintenance of a stable plateau. A brief fluctuation occurred between 38 and 42 ns, after which rapid re-stabilization was observed. The overall RMSD for this complex was 0.25 nm. In Figure 8C, STAT3-melilotocarpan A complex underwent relaxation within 0.3 ns, and then reached a stable plateau. The brief fluctuation occurred between 5 ns and 80 ns. The overall simulation was about 0.21 nm. As shown in The AKT1-melilotocarpan A, HIF-1α-melilotocarpan A, and STAT3-melilotocarpan complexes maintained stability without notable alterations, indicating a relatively stable combination.

The flexibility of amino acid residues in the proteins was assessed using root mean square fluctuation (RMSF). Figure 9A, the mean global residue fluctuation was about 0.17 nm, suggesting that stability was mainly maintained by the structural core and binding pocket. The higher flexibility was largely confined to intrinsically mobile regions such as terminal segments and loop regions rather than the pocket core. It is well known that these regions are flexible in nature. Figure 9B shows that approximately 94% of residues exhibited RMSF values below 0.20 nm, and only about 6% exceeded 0.30 nm, indicating that the overall protein backbone and internal residues fluctuated within a stable range. Figure 9C shows that the mean global residue was about 0.17 nm. No fluctuations exceeding 0.15 nm were identified in the data, indicating the binding pocket maintained a stable conformation throughout the simulation. Both the AKT1-melilotocarpan A and HIF-1α-melilotocarpan A complexes have been demonstrated to interact strongly and stably.

Hydrogen bonds are critical for promoting protein-ligand interactions. Figure 10A shows that the AKT1-melilotocarpan A complexes typically formed one hydrogen bond, with the number of hydrogen bonds varying between 0 and 2. Between melilotocarpan A and HIF-1α are separated by less than 0.35 nm. AKT1 frequently interacts in proximity with melilotocarpan A. The HIF-1α-melilotocarpan A complexes form 0 to 2 hydrogen bonds, primarily one hydrogen bond, as shown in Figure 10B. In Figure 10 C, the number of hydrogen bonds of STAT3-melilotocarpan A complex is zero to four. Most of the time, A atoms between melilotocarpan A and HIF-1α are separated by less than 0.35 nm. HIF-1α-melilotocarpan A complexes establish 1-3 pairs of close contacts and sometimes even 7-9 pairs. This suggests effective hydrogen-bonding interactions between the small molecule and target proteins. In summary, both the AKT1-melilotocarpan A and HIF-1α-melilotocarpan A complexes demonstrate stability and strength.

CC inhibits the inflammatory reaction and apoptosis of SONFH mice
Histological analysis using ABH staining revealed pronounced osteonecrotic changes in the SONFH group, characterized by increased empty lacunae, pyknotic nuclei, and disrupted trabecular architecture, whereas these pathological features were noticeably alleviated following CC treatment (Figure 11A). Quantitative analysis confirmed that the ratios of empty lacunae and pyknotic nuclei were significantly reduced in the CC-treated groups compared with the SONFH group (Figure 11B–C). Micro-CT analysis demonstrated severe trabecular bone loss in the SONFH group, which was partially reversed by CC administration (Figure 11D). Consistently, CC treatment significantly increased BV/TV and Tb. While decreasing Tb.Sp relative to the SONFH group (Figure 11E–G). Immunofluorescence staining showed elevated HIF-1α expression and reduced ALP levels in the SONFH group, accompanied by altered VEGF expression. CC treatment modulated HIF-1α and VEGF expression and enhanced ALP signals, indicating improved bone microenvironment and osteogenic activity (Figure 11H).

DATA AVAILABILITY:
The datasets used or analyzed during the current study are available at the link: https://zenodo.org/records/19730455.

Figure 1: Flowchart of the network pharmacology study strategy for CC amelioration of SONFH. Figure 1 shows that the study identifies potential targets of active ingredients from CC from the database, with SONFH as the core disease. After intersecting these with disease-related targets, a "component–target–pathway" network and GO, KEGG terms were constructed. Molecular docking and Molecular dynamics simulation are then employed to evaluate the binding affinity between key components and core targets, with final confirmation conducted through experimental verification. Please click here to view a larger version of this figure.

Figure 2: Screening of intersecting targets associated with SONFH. (A) Volcano plot depicting the distribution of differentially expressed genes in disease samples. Red points denote upregulated genes, blue points indicate downregulated genes, and grey points represent genes lacking significant differential expression. (B) Heatmap showing the expression patterns of the 60 differentially expressed genes, with columns corresponding to samples and rows corresponding to genes. (C) Venn diagram displaying the overlap of disease-related targets obtained from different databases. Please click here to view a larger version of this figure.

Figure 3: Screening of intersecting targets between CC and SONFH. (A) Venn diagram illustrating the distribution of 61 common targets shared between the predicted targets of active compounds in CC (yellow) and SONFH-related disease targets (purple). (B) The Herb-Compound-Target (H-C-T) network illustrates the interactions between compounds and their corresponding targets. The blue square node represents the disease, green square nodes denote active compounds, and orange square nodes indicate common targets. Edges denote the interactions between compounds and their targets. Please click here to view a larger version of this figure.

Figure 4: Candidate target identification through PPI analysis. (A) PPI network clustered with the MCODE plugin. (B) Schematic workflow of topological screening within the PPI network. (C) Key genes extracted from the PPI network by the CytoHubba plugin. Please click here to view a larger version of this figure.

Figure 5: GO enrichment analysis results and KEGG pathway enrichment analysis for the 61 common targets. (A) Bubble chart displaying the top 10 GO enrichment analysis terms for BP, CC, and MF. (B) Bubble plot depicting the top 20 significantly enriched KEGG pathways. (C) Distribution of the top 20 enriched pathways based on KEGG functional classification. (D) Illustrative C-T-P network depicting the potential mechanisms by which CC may improve osteonecrosis of the femoral head. Please click here to view a larger version of this figure.

Figure 6: Binding energy heatmap of the interactions between the active compounds of CC and the key targets (kcal/mol). Please click here to view a larger version of this figure.

Figure 7: Binding modes of key targets with specific active compounds. AKT1-melilotocarpan A (A1), HIF-1α-melilotocarpan A (B1), STAT3-melilotocarpan A (C1), ESR1-melilotocarpan A (D1), CASP3-melilotocarpan A (E1), SRC-melilotocarpan A (F1), EGFR-melilotocarpan A (G1), AKT1-anethole (H1), AKT1-myristicin (I1). (A2), (B2), (C2), (D2), (E2), (F2), (G2) and (H2) respectively illustrate their 2D binding modes. Please click here to view a larger version of this figure.

Figure 8: RMSD of MD. (A) RMSD values for AKT1-Melilotocarpan A complexes. (B) The RMSD values of HIF-1α-Melilotocarpan A complexes. (C) The RMSD values of STAT3-Melilotocarpan A complexes. Please click here to view a larger version of this figure.

Figure 9: RMSF of MD. (A) The RMSF values of AKT1-Melilotocarpan A complexes. (B) The RMSF values of HIF-1α-Melilotocarpan A complexes. (C) The RMSF values of STAT3-Melilotocarpan A complexes. Please click here to view a larger version of this figure.

Figure 10: H-bonds of MD. (A) The H-bond values of AKT1-Melilotocarpan A complexes. (B) The H-bond values of HIF1A-Melilotocarpan A complexes. (C) The H-bond values of AKT1-Melilotocarpan A complexes. Please click here to view a larger version of this figure.

Figure 11: CC inhibits inflammatory response and apoptosis in the femoral head, thereby improving SONFH. (A) Representative femoral head sections stained with ABH. Empty lacunae are indicated by black arrowheads, and pyknotic nuclei are indicated by arrows. (B–C) Quantitative analysis of the ratio of empty lacunae (B) and the number of pyknotic nuclei (C). n = 5. (D) Representative three-dimensional micro-CT reconstructions of femoral heads. (E–G) Quantification of micro-CT parameters, including bone volume fraction (BV/TV), trabecular separation(Tb.Sp), and trabecular thickness (Tb.Th). (H) Representative immunofluorescence staining of HIF-1α, VEGF, and ALP in femoral head sections; nuclei were counterstained with DAPI. Data are presented as mean ± SD (n = 5). *P < 0.05, **P < 0.01, ***P < 0.001. Please click here to view a larger version of this figure.

Term		Fold Enrichment	P-Value	Count	User Ids
Bladder cancer		28.44	1.23E-07	7	CREBBP,CXCL8,NOS2,MMP2,STAT3,F2,PTGS2, HIF1A,ESR1,MMP9,EGFR,MTOR,VEGFA,CASP3, ERBB2,EP300,PPARG,NFE2L2,BCL2L1
HIF-1 signaling pathway		16.66	6.56E-10	11	SRC,MMP2,STAT3,HIF1A,ESR1,MMP9,EGFR, MTOR,VEGFA,CASP3,ERBB2,KDR,PDCD4,PTPN6
EGFR tyrosine kinase inhibitor resistance		16.66	3.84E-07	8	CCR1,CREBBP,CXCL8,SRC,CASP3,STAT3 ,EP300,TYK2,PTGS2,HIF1A,MTOR,VEGFA
Adherens junction		12.54	1.64E-05	7	CREBBP,ABCB1,CASP3,ERBB2,STAT3,PDCD4, EP300,PTGS2,MMP9,EGFR,MTOR,VEGFA
Endocrine resistance		11.78	2.35E-05	7	CREBBP,NOS2,NOS3,ERBB2,STAT3, SERPINE1,EP300,HIF1A,EGFR,MTOR,VEGFA
AGE-RAGE signaling pathway in diabetic complications		11.55	2.63E-05	7	CCR1,CXCL8,SRC,CASP3,STAT3,CXCR2,PTGS2,EGFR,MTOR,VEGFA
Proteoglycans in cancer		11.43	1.19E-10	14	CXCL8,SRC,NOS3,CASP3,STAT3,PPARG,MMP9,BCL2L1,NFE2L2
Relaxin signaling pathway		10.25	1.02E-05	8	OXTR,NOS2,NOS3,ERBB2,PTAFR,KDR,NOS1,EGFR,VEGFA
Kaposi sarcoma-associated herpesvirus infection		10.20	1.45E-08	12	SRC,ERBB2,STAT3,KDR,EGFR,MTOR,BCL2L1,VEGFA
Thyroid hormone signaling pathway		9.56	7.64E-05	7	NOS2,SRC,NOS3,MMP2,NOS1,MMP9,EGFR,VEGFA
Fluid shear stress and atherosclerosis		8.21	0.000177	7	CREBBP,CXCL8,SRC,CASP3,STAT3,EP300,TYK2,MMP9
Hepatitis B		8.18	4.41E-05	8	CREBBP,STAT3,EP300,PTPN6,TYK2,EGFR,MTOR,BCL2L1
JAK-STAT signaling pathway		7.93	5.35E-05	8	SRC,VDR,STAT3,CYP3A4,ESR1,EGFR,MTOR,VEGFA
Human cytomegalovirus infection		7.34	6.2E-06	10	CREBBP,CASP3,EP300,TYK2,PTGS2,EGFR,MTOR,VEGFA
Lipid and atherosclerosis		6.94	0.000035	9	CXCL8,SRC,MMP2,ERBB2,MMP9,EGFR,VEGFA
MicroRNAs in cancer		6.25	2.06E-06	12	CREBBP,SRC,ERBB2,EP300,PTPN6,PTPRF,EGFR
Chemical carcinogenesis - receptor activation		6.14	0.000265	8	SRC,MMP2,ERBB2,ESR1,MMP9,EGFR,MTOR
Pathways in cancer		5.94	6.64E-10	19	CXCL8,NOS3,CASP3,MMP2,STAT3,SERPINE1,VEGFA
Calcium signaling pathway		5.90	0.00011	9	KAT2B,CREBBP,SRC,EP300,HIF1A,ESR1,MTOR
Human papillomavirus infection		4.00	0.00324	8	SRC,NOS3,MMP2,KDR,MMP9,NFE2L2,VEGFA

Table 1: The KEGG Enrichment Results of the Top 20 Enriched Pathways.

Supplementary Table 1: Basic information of active compounds in CC.Please click here to download this file.

Supplementary Table 2: The targets of 85 active components in CC.Please click here to download this file.

Supplementary Table 3: Information on 61 CC-SONFH common targets.Please click here to download this file.

Supplementary Table 4: Results of the biological process category terms from GO enrichment analysis.Please click here to download this file.

Supplementary Table 5: Results of the cellular component categorization terms from GO enrichment analysis.Please click here to download this file.

Supplementary Table 6: Results of the molecular function category terms from GO enrichment analysis.Please click here to download this file.

Supplementary Table 7: Results of the pathways from KEGG enrichment analysis.Please click here to download this file.

Supplementary Table 8: Molecular docking binding energies (kcal/mol) of compound-target pairs.Please click here to download this file.

Supplementary Table 9: Details of targets and components for molecular docking.Please click here to download this file.

Discussion

In this study, the authors obtained 85 bioactive compounds from CC through several databases and performed screening, among which the main chemical components include anethole, melilotocarpan A, and myristicin. To identify targets, the authors intersected 563 drug CC targets with 1116 SONFH-related targets, ultimately getting 61 targets. The authors imported 61 common targets into the STRING and DAVID databases to build a PPI network and investigate possible pharmacological pathways. The results revealed 10 key targets in the PPI network, including HIF-1α and STAT3. Ischemia and hypoxia are two of the main pathogenic characteristics of SONFH, according to earlier research³⁵. CC may exert positive effects in SONFH through multiple biological processes, including angiogenesis and transcriptional regulation. This could be accomplished through the HIF-1α signaling pathway, fluid shear stress and atherosclerosis, lipids and atherosclerosis, the AGE-RAGE signaling pathway in diabetic complications, and other signaling pathways, according to KEGG enrichment analysis. This suggests that CC may act on SONFH through these key targets and the associated pathways. To further clarify whether the components of CC are tightly coupled with the targets, based on the above network pharmacology research results, namely the screening results of the active ingredients of the compound prescription and the screening results of the targets in the PPI network, the authors conducted molecular dynamics simulations of the melilotocarpan A-AKT1, melilotocarpan A-HIF-1α, and STAT3-melilotocarpan complexes. This shows melilotocarpan A binds to AKT1, HIF-1α, and STAT3 steadily, so CC may indeed exert positive effects by acting on these key targets. All in all, through a multi-component, multi-target, multi-pathway mechanism, CC may exert protective effects on the quality of life in SONFH patients.

The research group has focused on elucidating the mechanisms by which YGPs exert therapeutic effects in bone-related diseases. The previous research indicated that YGPs enhance bone formation and improve trabecular microarchitecture in the femoral head by activating β-catenin to inhibit osteoclastogenesis and promote osteogenesis during rabbit SONFH³⁶. Furthermore, the researchers have found that YGPs exert therapeutic effects on SONFH primarily by alleviating inflammation and promoting angiogenesis¹². YGPs can not only ameliorate SONFH but also improve osteoporosis. By inhibiting the IL-17/NF-κB signaling pathway and decreasing Th17 immune responses, YGPs successfully stop ovariectomy-induced bone loss³⁷. Through animal experiments, the researchers confirmed that Cornus officinalis, one of the components of YGPs, exerts effects on SONFH by inhibiting the secretion of inflammatory mediators and bone cell apoptosis³⁸. CC, which can warm yang and benefit qi and unblock meridians, is the principal herb in YGPs, playing a crucial role in improving SONFH and relieving patient suffering. In summary, these studies reflect the long-standing focus on elucidating the mechanisms of YGPs in osteonecrosis and on pinpointing CC as a principal herb worth mechanistic investigation in SONFH.

Among the predicted active ingredients, anethole, myristicin, and other ingredients have exhibited biological activity against bone-related diseases. Anethole has been demonstrated to inhibit the elevation of bone resorption markers. Anethole ultimately reduces osteoclast differentiation and resorptive function by inhibiting downstream pathways and key regulatory factors³⁹. Research indicates that anethole exhibits anti-inflammatory properties and inhibits adipogenic differentiation in human bone marrow mesenchymal stem cells (hBMSCs)⁴⁰. Myristicin inhibits inflammation and protects vascular smooth muscle cells by suppressing the PI3K/Akt and NF-κB pathways⁴¹. Ferroptosis is one of the key mechanisms underlying SONFH pathogenesis. Ferroptosis is a type of programmed cell death that is caused by lipid peroxidation dependent on iron. Myristicin prevents ferroptosis, thereby protecting the redox homeostasis of osteoblast cell membranes⁴². The other ingredient of borneol has been confirmed to inhibit the formation of actin rings, which is a feature of resorbing osteoclasts that reflects cell polarization, within 30 min⁴³. These findings indicate that these components may be very important for the protective effect of CC in SONFH, and it is worthy of further exploration. Within all these ingredients, Melilotocarpan A was selected as the representative ligand because it exhibited the most stable binding affinities with multiple core targets in the MD analysis. Among the top-ranked compounds, Melilotocarpan A showed consistently strong binding interactions with several key proteins in the PPI network, suggesting that it may represent a major bioactive component of CC.

In addition, the results of the PPI network revealed that, especially HIF-1α and STAT3, may be the core targets. HIF-1α is the central mediator of the cellular response to hypoxia. In SONFH experiments, activating HIF-1α enhances angiogenesis and bone repair, thereby alleviating femoral head necrosis⁴⁴. Additionally, research has found that HIF-1α can exert a protective effect on cartilage tissue under hypoxic conditions⁴⁵. Upregulating HIF-1α suppresses ferroptosis and activates STAT3 to accelerate the differentiation of bone marrow monocytes into osteoclasts⁴⁶. STAT3 is a signal transduction and transcription activator. HBMSCs promote the conversion of STAT3 to p-STAT3, thereby enhancing chondrocyte proliferation, migration, and anti-apoptotic effects under hypoxic conditions⁴⁷. STAT3 enhances the osteogenic differentiation and anti-apoptotic effects of hBMSCs while increasing their vascular endothelial growth factor (VEGF) secretion to promote bone microvascular regeneration⁴⁸. AKT1 is a protein kinase that acts as a key mediator of angiogenesis signaling. AKT1 is involved in bone angiogenesis-osteogenesis coupling and trabecular bone formation⁴⁹. All in all, these hub targets converge on angiogenesis, hypoxia adaptation, inflammation, and bone remodeling, suggesting that CC may exert protective effects against SONFH by modulating these targets and their associated signaling pathways.

To explore the mechanism of CC in the amelioration of SONFH, the researchers conducted GO analysis and KEGG enrichment analysis. GO results show that the target genes are mainly enriched in biological functions such as response to hypoxia, cellular response to lipopolysaccharide, and positive regulation of angiogenesis. KEGG enrichment analysis suggested that the pharmacological effects of CC in SONFH are predominantly linked to the AGE-RAGE signaling pathway in diabetic complications, HIF-1α signaling pathway, Lipid and atherosclerosis, and other signaling pathways. SONFH occurs when blood flow to the femoral head is interrupted, leading to a lack of nutritional support for the femoral head and triggering cell apoptosis⁵⁰. HIF-1 consists of HIF-1α and HIF-1β, with HIF-1α playing a key role in the transcriptional response to hypoxic and ischemic environments⁴⁷. Previous studies have demonstrated that, through the VEGF/AKT/mTOR signaling cascade, HIF-1α coordinates osteogenesis and angiogenesis, thereby critically enhancing the osteogenic differentiation of adipose-derived stem cells. After ischemia occurs, tissue hypoxia develops, leading to a marked elevation of HIF-1α within cells. which boosts VEGF expression levels and consequently promotes vascular repair and regeneration⁵¹. Hypoxia disrupts mitochondrial energy production and increases the formation of reactive oxygen species (ROS). These changes promote osteoclastogenesis, inhibit osteoblast activity, and ultimately induce osteocyte apoptosis. High HIF-1α expression shifts cellular energy metabolism from oxidative phosphorylation to glycolysis, reduces ROS production, and thereby attenuates apoptosis in osteoblasts and osteocytes^52,53. Cancellous bone supports the hematopoietic function of bone marrow and stores minerals such as calcium. Increased activity of HIF-1α can directly increase cancellous bone formation⁵⁴. Additionally, HIF-1α maintains chondrocyte survival and homeostasis under hypoxic conditions⁴⁵. Consequently, the HIF-1α signaling pathway has become a therapeutic target for ONFH, with the potential therapeutic application for ONFH.

In the AGE-RAGE signaling pathway, Advanced glycation end-products (AGEs) result from the non-enzymatic covalent cross-linking of carbohydrates with proteins, fats, or other biological macromolecules⁵⁵. AGEs can damage the cells and the tissues via inflammatory and oxidant damage. Through four pathways, signals are transduced by the AGE-RAGE interaction. They are 1) JAK-2-STAT1, (2) PI3K-AKT, (3) MAPK-ERK, and (4) NADPH oxidase-ROS⁵⁶. In the end, the phosphorylated NF-κB enters the nucleus to transcribe the expression of proinflammatory cytokines, growth factors, profibrotic cytokines, and oxidative stress. The study has shown that AGE-RAGE promotes osteoblast apoptosis via the MAPK signaling pathway and the activation of oxidative stress, and inhibits osteogenic differentiation by suppressing the levels of endoplasmic reticulum stress sensors, as well as through DNA methylation/the Wnt pathway⁵⁷. In addition, studies have indicated that AGEs and RAGEs take part in vascular calcification⁵⁸.

In the lipid and atherosclerosis signaling pathway, oxidized lipids activate PPARγ. It is a master regulator of adipogenesis that inhibits osteogenic differentiation⁵⁹. In addition, Several studies have shown that lipids affect the differentiation and maturation of osteoclasts and osteoblasts, thereby disrupting bone homeostasis⁶⁰. Multiple foundational studies have demonstrated the significant involvement of apoptosis in the physiological and pathological mechanisms of SONFH. Subsequent in vivo animal experiments further confirmed that the active fraction of CC effectively suppresses apoptosis and inflammation in SONFH, thereby mitigating disease progression.

Consistent with these findings, network pharmacology and experimental analyses suggest that CC may improve SONFH by coordinating hypoxia adaptation, regulating angiogenesis, and modulating metabolic processes. Quantitative micro-CT analysis revealed that CC intervention significantly restored bone volume fraction (BV/TV), reduced trabecular separation (Tb.Sp, Tb.Th), and promoted the recovery of bone microstructure. KEGG enrichment analysis further highlighted the HIF-1α signaling pathway as a key mechanism underlying these effects. Collectively, these animal experiments indicate that CC may alleviate SONFH by regulating hypoxia-related and inflammation-related biological processes and promoting angiogenesis and bone remodeling, potentially through key pathways such as the HIF-1α signaling pathway.

Although some important preliminary findings were obtained, this study still has some limitations. The authors relied solely on modern bioinformatics methods, including network pharmacology and MD, as well as some animal experiments, to explore the role of CC in SONFH. Therefore, the reliability and accuracy of the predictions need to be further verified in vivo and in vitro experiments. This study provides preliminary clues into the mechanisms by which CC ameliorates SONFH, which require further experimental validation. By integrating network pharmacology with MD, the authors systematically mapped the candidate bioactive ingredients, key targets, and enriched pathways of CC in SONFH. In conclusion, this study investigated the potential protective effects of CC in SONFH through integrative network pharmacology and in vivo validation. The results indicate that the mechanisms were primarily associated with modulation of the HIF-1 signaling pathway and promotion of osteogenic regeneration. As illustrated in the corresponding figure, CC amelioration significantly modulated the expression of key proteins, including HIF-1α, VEGF, and ALP, in a dose-dependent manner, further supporting its role in coordinating hypoxia adaptation, angiogenesis, and osteogenic activity. These findings were corroborated by micro-CT analysis, which demonstrated that CC intervention effectively restored bone volume fraction and improved trabecular microstructure.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
4% PFA	biosharp	BL539A
Absolute ethenol	Sinopharm Group Co.Ltd	10009218
Alcian Blue Stain Kit	Solarbio	G1563
ALP	ariggo	ARG57422
AutoDock Vina	Scripps Research		Version 1.2.7
AutoDockTools	Scripps Research		Version 1.5.7
Bioinformatics online platform			https://www.bioinformatics.com.cn
CHARMM36 force field			Used for MD parameterization
CTD			https://ctdbase.org
Cytoscape	Cytoscape Consortium		Version 3.10.3; plugins: CytoHubba, CytoNCA
DAPI Staining Solution	Beyotime Biotech Inc	C1006-50mL
DAVID	NCI		https://david.ncifcrf.gov
Discovery Studio Visualizer	BIOVIA		BIOVIA Discovery Studio 2020; docking visualization
DisGeNET			https://disgenet.com
EDTA Decalcification	Beyotime Biotech Inc	C0167-3L
Endogenous Peroxidase Blocking Solution	Beyotime Biotech Inc	P0100A
GeneCards			https://www.genecards.org
GEO (GSE123568)	NCBI		https://www.ncbi.nlm.nih.gov/geo/
Goat anti-rabbit IgG(H+L)(Alexa Fluor 488)	CST	4409
Goat Serum	Beyotime Biotech Inc	C0265
GraphPad Prism	GraphPad (Dotmatics)		Version 10; statistical analysis and graphing
GROMACS	GROMACS		Version 2022; CHARMM36 force field
Herb 2.0			http://herb.ac.cn/v2/
HIF-1α	Hangzhou Huaan	HA721997
high resolution micro-CT equipment	Bruker	SkyScan
Lipinski Rule Filter			Mw ≤ 500; miLogP ≤ 5; HBD ≤ 5; HBA ≤ 10
lipopolysaccharide	Sigma-Aldrich	L4516
methylprednisolone	Sinopharm Group Co.Ltd	CATOCCAD302504100MG
Neutral balsam	biosharp	BL704A
OMIM			https://www.omim.org
Paraffin	Sinopharm Group Co.Ltd	C416770020
Particle Mesh Ewald (PME)			Long-range electrostatics
PubChem	NIH		https://pubchem.ncbi.nlm.nih.gov
PyMOL	Schrödinger		Protein visualization and preparation
R	CRAN		Version 4.4.3
RCSB PDB			https://rcsb.org
STRING			https://string-db.org
SwissADME	Swiss Institute		http://www.swissadme.ch
SwissTargetPrediction	Swiss Institute		http://www.swisstargetprediction.ch
TIP3P water model			Used for solvation
UniProt			https://www.uniprot.org
VEGF	Hangzhou Huaan	ET1604-28
Weishengxin			https://www.bioinformatics.com.cn
Windows 11	Microsoft		Used for data processing, molecular docking, and MD simulations
Xylene	Sinopharm Group Co.Ltd	10023418